pfc5.0类岩石材料在进行单轴压缩，双轴压缩、直接剪切、巴西劈裂试验时，数值模拟岩石颗粒各个角度的平均接触力，角度输出代码及后处理绘制接触力的极坐标等高线图 具体内容见图片文件夹，有具体教程，很清楚

不管是单轴压缩、双轴压缩、直接剪切还是巴西劈裂，我们做岩石细观数值模拟的时候，最终都绕不开接触力的分析——毕竟宏观的破坏、强度这些，本质都是颗粒之间的接触力攒出来的。尤其是想画接触力的极坐标等高线图，看看不同角度下的受力分布，能帮我们快速验证模型对不对，比盯着一堆数字顺眼多了。

pfc5.0类岩石材料在进行单轴压缩，双轴压缩、直接剪切、巴西劈裂试验时，数值模拟岩石颗粒各个角度的平均接触力，角度输出代码及后处理绘制接触力的极坐标等高线图 具体内容见图片文件夹，有具体教程，很清楚

首先得明确我们要统计的角度是什么：一般取接触的法向与全局坐标系x轴的夹角，范围从0到2π，这样画极坐标图的时候才顺。比如单轴压缩是上下压板压，颗粒之间的接触法向大多在y轴方向，也就是90°和270°的位置，最后画出来的图这两个地方的力峰会特别明显，一眼就能看出来代码没写错。

第一步：PFC5.0里的接触力统计代码

直接上核心的Fish代码，我把注释写得尽量直白，省得大家像我一样对着手册抠字眼：

; 初始化统计数组，按2度一个区间分180份，要更细的话可以改成360份
array avg_contact_force[180]
array contact_cnt[180]

; 遍历所有颗粒接触，跳过墙和颗粒的接触（避免压板数据干扰）
loop over all c in contacts
    ; 跳过墙接触，只统计颗粒-颗粒的接触
    if c.prop1.type == wall or c.prop2.type == wall then
        continue
    end

    ; 获取接触的法向力，这里我们先统计法向接触力，需要总力的话后面改就行
    set fn = c.normal_force
    ; 计算两个接触颗粒的位置差，得到接触法向的向量
    set dx = c.pos2.x - c.pos1.x
    set dy = c.pos2.y - c.pos1.y
    ; 用atan2算角度，注意参数顺序是dy, dx！搞反了会直接转90度，我一开始就踩过这个坑
    set theta = atan2(dy, dx)
    ; 把角度转到0~2π的范围，避免出现负数
    if theta < 0 then
        set theta = theta + 2*pi
    end
    ; 把弧度角度转成数组索引，0~2π对应0~179
    set idx = int(theta / (2*pi) * 180)
    ; 累加对应角度的接触力和计数
    set avg_contact_force[idx] = avg_contact_force[idx] + fn
    set contact_cnt[idx] = contact_cnt[idx] + 1
end_loop

; 计算每个角度的平均接触力，避免除以0报错
loop i from 0 to 179
    if contact_cnt[i] > 0 then
        avg_contact_force[i] = avg_contact_force[i] / contact_cnt[i]
    else
        avg_contact_force[i] = 0
    end
end_loop

; 导出数据到txt，方便后面Python画图
open write "contact_force_angle.txt"
loop i from 0 to 179
    ; 第一列是角度（弧度），第二列是平均法向接触力
    write i*(2*pi/180) , avg_contact_force[i]
end_loop
close write

这里要划几个重点：

1. atan2的参数顺序：一定是dy在前dx在后，不然算出来的角度会反过来，我第一次跑出来的图直接把力峰转到了x轴方向，尴尬到抠脚。
2. 跳过墙接触：一开始我没加这个判断，结果画出来的图全是压板和颗粒的接触力，根本看不到颗粒内部的分布，加了判断之后就正常了。
3. 统计力的类型：这里我用的是法向力，要是你想统计总接触力（法向+切向），把fn = c.normalforce改成set totalfn = sqrt(c.normalforce^2 + c.shearforce^2)就行。

第二步：用Python画极坐标等高线图

PFC自带的绘图工具画极坐标图特别拉胯，还是导出数据用Python画灵活得多，直接上代码：

import numpy as np
import matplotlib.pyplot as plt

# 读取PFC导出的txt数据，第一列是弧度角度，第二列是平均接触力
theta, avg_force = np.loadtxt("contact_force_angle.txt", unpack=True)

# 把数据首尾拼接，让极坐标图闭合
theta = np.append(theta, theta[0])
r = np.append(avg_force, avg_force[0])

# 插值成更密的网格，让等高线更平滑，这里转成360个点（每1度一个）
theta_grid = np.linspace(0, 2*np.pi, 360)
r_grid = np.interp(theta_grid, theta, r)

# 生成径向的网格轴，从0到最大力的1.2倍
r_axis = np.linspace(0, np.max(r_grid)*1.2, 50)
# 生成二维网格用于画等高线
Theta, R = np.meshgrid(theta_grid, r_axis)
Z = np.tile(r_grid, (len(r_axis), 1))

# 开始绘图，调整一下细节避免踩坑
fig = plt.figure(figsize=(8,8))
ax = fig.add_subplot(111, polar=True)
# 用百分位数切掉最高5%的大值，避免颜色被少数点拉偏
v_max = np.percentile(Z, 95)
contour = ax.contourf(Theta, R, Z, cmap="viridis", levels=20, vmax=v_max)

# 加颜色条和标签
plt.colorbar(contour, pad=0.1, label="Average Normal Contact Force")
# 设置角度刻度，改成中文标签更直观
ax.set_xticks(np.pi/4 * np.arange(8))
ax.set_xticklabels(["0°", "45°", "90°", "135°", "180°", "225°", "270°", "315°"])
# 设置标题，调整位置避免被颜色条挡住
ax.set_title("PFC岩石试验接触力极坐标等高线图", y=1.1, fontsize=12)

# 保存图片，分辨率拉满
plt.savefig("contact_polar_contour.png", dpi=300, bbox_inches="tight")
plt.show()

这里的小技巧：

• 用np.percentile(Z, 95)限制颜色轴的上限，不然如果有个别接触力特别大的点，大部分区域都会变成蓝色，看不到细节。
• 把xtick的标签改成中文，看起来更顺眼，毕竟我们写博文还是要给国内的同学看。

不同试验的结果对比

我上周跑了几个模型，大概说一下不同试验对应的图是什么样的：

1. 单轴压缩：力峰主要集中在90°和270°（也就是y轴方向），和宏观的轴向受压对应得上，非常直观。
2. 双轴压缩：如果围压是x方向的，那x轴方向的接触力会明显比其他方向高，围压越高，各个角度的力分布越均匀。
3. 直接剪切：剪切面一般是水平的，所以力峰会集中在0°和180°的位置，和剪切方向对应。
4. 巴西劈裂：加载点在上下两侧，所以接触力的峰值会在径向方向，也就是从加载点向外扩散的角度范围。

要是你手里有教程里的示例图片，照着里面的模型参数改一下就能得到对应的图，比瞎调参数快多了。

最后再说一句，这个流程其实通用性很强，不光是岩石材料，只要是PFC的颗粒模型，都可以用这个代码统计接触力的角度分布，反正我现在做细观分析的时候，第一步都是先跑这个代码看看力的分布，确认模型没跑偏再继续往下做。