说起来上周还在为怎么把PFC5.0里颗粒接触力按角度统计出来头疼，翻了好几篇教程终于摸清楚门道，今天把整个流程捋一遍，顺便把踩过的坑都标出来

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这里要划几个重点atan2的参数顺序一定是dy在前dx在后不然算出来的角度会反过来我第一次跑出来的图直接把力峰转到了x轴方向尴尬到抠脚。跳过墙接触一开始我没加这个判断结果画出来的图全是压板和颗粒的接触力根本看不到颗粒内部的分布加了判断之后就正常了。统计力的类型这里我用的是法向力要是你想统计总接触力法向切向把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, unpackTrue) # 把数据首尾拼接让极坐标图闭合 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, polarTrue) # 用百分位数切掉最高5%的大值避免颜色被少数点拉偏 v_max np.percentile(Z, 95) contour ax.contourf(Theta, R, Z, cmapviridis, levels20, vmaxv_max) # 加颜色条和标签 plt.colorbar(contour, pad0.1, labelAverage 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岩石试验接触力极坐标等高线图, y1.1, fontsize12) # 保存图片分辨率拉满 plt.savefig(contact_polar_contour.png, dpi300, bbox_inchestight) plt.show()这里的小技巧用np.percentile(Z, 95)限制颜色轴的上限不然如果有个别接触力特别大的点大部分区域都会变成蓝色看不到细节。把xtick的标签改成中文看起来更顺眼毕竟我们写博文还是要给国内的同学看。不同试验的结果对比我上周跑了几个模型大概说一下不同试验对应的图是什么样的单轴压缩力峰主要集中在90°和270°也就是y轴方向和宏观的轴向受压对应得上非常直观。双轴压缩如果围压是x方向的那x轴方向的接触力会明显比其他方向高围压越高各个角度的力分布越均匀。直接剪切剪切面一般是水平的所以力峰会集中在0°和180°的位置和剪切方向对应。巴西劈裂加载点在上下两侧所以接触力的峰值会在径向方向也就是从加载点向外扩散的角度范围。要是你手里有教程里的示例图片照着里面的模型参数改一下就能得到对应的图比瞎调参数快多了。最后再说一句这个流程其实通用性很强不光是岩石材料只要是PFC的颗粒模型都可以用这个代码统计接触力的角度分布反正我现在做细观分析的时候第一步都是先跑这个代码看看力的分布确认模型没跑偏再继续往下做。