企业官网建设流程全解析

如果你已经用 PFC5.0 颗粒流离散元做过几次岩土模拟，大概率会在同一个路口迷路：试样生成完了，墙也建好了，一加载，应力应变曲线却完全不像土。我整理这套岩土单轴、双轴及三轴测试算例代码集，就是为了让“跑通一个测试”这件事不再靠运气。这篇文章会把代码集的目录逻辑、关键写法和每组测试背后的力学思路拆开讲。适合两类人：一类是刚开始做离散元、被建模流程劝退的新手；另一类是已经能跑通简单模型，但总在伺服控制、参数标定和后处理之间反复折腾的岩土工程师。代码集是骨架，你能读懂为什么这样写，才能把它改成自己的东西。

1. 为什么算例代码集要把三种测试拆开维护

1.1 单轴、双轴、三轴分别对应哪种岩土受力状态

很多人以为单轴、双轴、三轴只是“多几面墙”的关系，其实三者对应的是完全不同的岩土受力状态，代码写法和结果解读逻辑也不一样。

单轴压缩通常模拟无侧限约束条件下的强度，试样只受轴向加载，侧向自由变形。它最常用于确定无侧限抗压强度（UCS），对黏性土、岩石类脆性材料特别有意义。在 PFC 里，这个测试的难点不是“不会建墙”，而是怎么让侧向约束彻底消失，以及怎么控制端部摩擦。

双轴测试在离散元里更多指的是二维平面应变模型。矩形试样在 x 方向施加恒定围压，y 方向做轴向加载，另一个方向要么被约束、要么忽略。它与单轴的最大区别在于侧向应力不为零，土体在围压下的剪切特性、剪胀行为才能表现出来。边坡、地基承载力这类平面应变问题，用双轴测试比用单轴更有说服力。

三轴测试则是三维圆柱试样的经典室内试验，σ2 = σ3，由围压控制，σ1 由轴向加载控制。它最接近真实土体的受力状态，也是离散元参数标定里最常用的一组数据。

把这三种测试放在一张表里看会更清楚：

看清楚这一点后，你就会明白为什么代码不能混着写。单轴的边界条件是“无侧限”，双轴的边界条件是“恒围压”，三轴的边界条件是“轴对称恒围压”。三套边界逻辑不同，强行塞进一个脚本，后面改参数时一定会互相干扰。

1.2 代码集按“生成-加载-后处理”三层组织

我整理代码集时用的目录结构很简洁：

pfc5_test_cases/ ├── 01_common/ │ ├── sample_generation.p2dat │ ├── contact_parameters.csv │ └── measure_setup.p2dat ├── 02_ucs/ │ ├── ucs_main.p2dat │ └── ucs_postprocess.fis ├── 03_biaxial/ │ ├── biax_main.p2dat │ ├── servo_control.fis │ └── biax_postprocess.fis ├── 04_triaxial/ │ ├── triax_main.p2dat │ ├── triax_confine.p2dat │ └── triax_postprocess.fis └── 05_postprocess/ ├── plot_stress_strain.py └── read_history.py

这样拆分的好处是：三种测试共用同一套试样生成和接触参数，差异只体现在边界加载上。如果你把单轴、双轴、三轴分别做成三个完整的独立脚本，后续改一次孔隙率或粒径范围，就要同步改三份文件，很容易漏。

实际使用中，我先单独跑通01_common里的试样生成，确认孔隙率、初始接触状态正常，再分别把三套加载逻辑挂上去。这个顺序能减少大量调 bug 时间。

2. 所有测试共用的模型生成骨架

2.1 颗粒级配与目标孔隙率：先算体积再生成

不管是单轴、双轴还是三轴，第一步都是生成一个初始稳定的颗粒体。这里最容易犯的错是：还没想清楚目标孔隙率就随手生成一堆颗粒，接着发现试样体积不对、空隙太多或太少。

正确顺序是先定模型尺寸和粒径范围，再算目标孔隙率对应的颗粒体积。以二维双轴试样为例，假设试样宽 50 mm、高 100 mm，目标孔隙率 0.18，半径范围 0.5~1.0 mm。颗粒骨架体积约等于试样体积乘以 (1 - n)，也就是：

颗粒总体积 = 宽 × 高 × (1 - 0.18)

PFC5.0 里可以直接用ball distribute按孔隙率生成，但要注意这个孔隙率是生成时的几何孔隙率，不是最终平衡后的孔隙率。生成后必须用一个测量区域（measure region）去实测，必要时通过删除局部颗粒或调整半径来校准。

这段代码是示意核心逻辑：

def generate_sample ball distribute radius 0.5e-3,1.0e-3 ... porosity 0.18 ... box 0,0.05 0,0.1 cycle 2000 end

我建议在生成后先查看孤立颗粒：接触数小于 3 的颗粒通常意味着试样内部有悬浮状态，这些颗粒会在加载初期突然移动，造成应力曲线前段异常。把它们删掉后重新平衡，再进入加载阶段。

2.2 接触模型怎么选：线性模型、平行粘结模型与滞后接触

接触模型选错，后面所有结果都会失真。PFC5.0 中常用的是线性接触linear，它只有法向刚度、切向刚度和摩擦系数，适合模拟砂土、碎石这类无粘结材料。对于岩石或胶结性土，需要改用平行粘结linearpbond，它能在接触点传递拉力和弯矩，模型才可能出现脆性破坏。

新手最容易犯的错是上来就挂高级模型，比如流体耦合、滞后接触、完整三维颗粒形状，结果一个参数变化引发连锁反应，最后根本定位不到问题。我的建议是先用最简单的模型把边界、伺服和加载逻辑跑通，再逐步替换接触本构。

如果目标是做砂土双轴或三轴，linear配合合适的摩擦系数和刚度足够应对大多数标定需求。做岩石单轴时再换成linearpbond，并单独标定粘结强度、粘结刚度。

2.3 参数标定为什么必须放在测试代码之前

离散元的微观参数和室内试验的宏观参数不是一一对应关系，这是不少刚入门朋友最难接受的一点。你没法直接说“土的弹模是 50 MPa，所以接触刚度就是 50 MPa”，必须通过数值试算去建立映射关系。

参数标定不能等三轴测试跑完再反推。一个合理的流程是：先用单轴压缩定弹性模量和破坏强度，再用双轴或三轴的不同围压水平校核摩擦角和剪胀特性。如果顺序反过来，你会发现改摩擦系数会同时影响弹模和峰值强度，改刚度又会影响破坏模式，最后陷入无意义的循环调参。

代码集里单独放了一张参数表，而不是把参数散落在各个脚本里：

contact_parameters.csv kn,ks,fric,pb_kn,pb_ks,pb_ten,pb_coh

每次跑一个算例前，脚本自动读取这张表，保证三组测试用的是同一套参数，否则标定无从谈起。

3. 单轴压缩测试：最容易被低估的边界条件

3.1 加载墙的生成与端部约束处理

单轴测试看起来最简单：上下各一面加载墙，把试样放在中间，让上面的墙以恒定速率向下移动。但这里有一个非常隐蔽的坑：试样生成时通常会有一个临时侧向边界墙，加载步骤开始后如果忘记删除或移走它，就等于偷偷给试样加了侧向约束，单轴变成了侧限压缩。

正确做法是，在平衡阶段保留临时边界来稳定试样，加载阶段开始前把侧向边界删除，只保留上下加载墙。代码逻辑类似：

def setup_ucs wall delete walls range group 'lateral' wall property fric 0.05 range group 'top' wall property fric 0.05 range group 'bottom' end

端部摩擦也值得注意。加载墙的摩擦系数如果很大，试样端部几乎不能横向滑移，破坏时会形成明显的锥形或双锥形破坏带，和真实无侧限压缩试验中的端部效应一致。若想模拟理想光滑压板，就把加载墙摩擦设小；若想匹配室内试验，需要根据压板材质单独标定。

3.2 加载速率与准静态判定

PFC 是显式动力求解，加载速率一旦太快，惯性力就会参与受力平衡，最典型的表现是峰值强度虚高、峰后曲线锯齿状明显。判断是否处于准静态，常用的指标是动能与总应变能的比值，建议控制在 5% 以下。

经验上，单轴加载速率可以先用“轴向应变率”来设置，比如让压板速度等于 0.001 乘以试样高度每秒，也就是每秒产生 0.1% 的应变。这是一个保守且稳定的起点。跑完后查看系统动能曲线，如果动能一直偏高，就把速度降一半再跑。

代码里可以用这样的伪代码判断：

def check_quasistatic kinetic_energy = mech.kinetic strain_energy = mech.strain ratio = kinetic_energy / (strain_energy + 1e-12) if ratio > 0.05 local msg = '通过结构' endif end

加一个自动判断比人眼盯曲线更可靠，尤其当你批量跑多组围压时，这个检查能帮你过滤掉无效结果。

3.3 应力应变曲线计算与破坏模式识别

计算轴向应力时，很多人直接用“轴向力除以初始截面积”，这在弹性小变形阶段问题不大，但在大变形阶段会有误差，因为试样侧向膨胀后真实截面积已经变了。更合理的做法是用当前截面积计算。

轴向应变由压板位移换算：

axial_strain = (初始高度 - 当前高度) / 初始高度

应力则取加载墙的接触力合力，除以试样当前截面积。在 PFC 里获取墙的接触力并不复杂，关键是确定“当前截面面积”。对于圆柱试样，可以按当前径向平均位置估算面积；对于二维模型，直接用当前宽度乘以单位厚度。

破坏模式的识别不能只看应力峰值。单轴低围压下，脆性模型容易出现贯穿剪切带或劈裂；如果接触参数给的粘结强度过高，可能出现整体“炸开”的异常破坏，这时候要检查是不是加载速率过大或粘结刚度与颗粒刚度不匹配。

4. 双轴测试里的平面应变与伺服控制

4.1 双轴和单轴的本质差异

双轴测试不是单轴测试加一面墙，而是让试样在一个方向承受恒定围压，另一个方向做轴向压缩。这里的核心是“恒定围压”，它需要墙根据实时应力调整位置，而不是固定不动。固定墙在加载过程中会因为试样膨胀产生不可控的应力增长，那就不是双轴试验了。

在二维双轴模型里，通常 x 方向是围压方向，y 方向是轴向加载方向。试样受轴向压缩时会在 x 方向膨胀，如果左右两面墙不做调整，围压会急剧上升。伺服控制的作用就是让墙始终把围压稳定在目标值附近。

4.2 围压伺服控制的核心逻辑

伺服控制的逻辑并不复杂：墙上测得当前应力，和目标围压比较，得到误差，再用误差乘一个增益系数得到墙的移动速度。写成代码就是：

def servo_wall while running measured = stress_of_wall('side') error = target_confine - measured if abs(error) > tolerance wall_vel = gain * error else wall_vel = 0.0 endif cycle 1 endloop end

这里的gain是整个伺服控制里最敏感的参数。gain太小，墙来不及响应，围压可能长期偏低；gain太大，墙来回振荡，围压曲线像心电图。一个合理的估算方法是根据当前墙壁接触刚度和接触数量，算出控制墙体位移到目标应力所需的参考速度数量级，再乘以一个折减系数让系统稳定。

实际调参时，我一般先设一个很小的gain，观察围压是否缓慢收敛；收敛但太慢，就逐步放大；出现振荡，立即减半。手工调一轮之后，把稳定的参数写回代码，后续批量跑围压就很少出问题。

4.3 横向墙位移方向问题与调试顺序

双轴代码里最让我头疼的错误是横向墙位移方向搞反。矩形试样有左右两面墙，一面需要向左移动、一面需要向右移动。如果给两边设成同向速度，试样会被整体平移，围压应力场完全不对称。

调试顺序非常重要。建议先不开轴向加载，只让横向墙伺服，待左右墙应力都稳定在目标围压后，再开始轴向压缩。这样能把“伺服问题”和“加载问题”分开。如果你一上来就同时加载和伺服，曲线一乱，很难判断是墙速方向错了还是增益不对。

另外，每次进入加载阶段前，把墙的速度清零，避免上一轮残留速度干扰新一轮计算。这个操作看着小，但能省下大量重复试错时间。

5. 三轴测试：围压施加方式决定模拟可信度

5.1 常规三轴与真三轴的离散元区别

常规三轴就是 σ2 = σ3，侧向围压用一个值控制，轴向加载产生 σ1。真三轴则需要三个主应力独立控制，在数值模拟中要用三组互相独立的边界系统，难度直线上升。

对大多数土工问题，常规三轴已经够用，代码集里也建议先做常规三轴。真三轴更适合研究应力路径复杂、各向异性明显的问题，如果不是特别需要，不要一开始就上真三轴，否则容易被边界实现细节拖住。

5.2 刚性墙伺服、柔性膜颗粒与周期边界的适用场景

三轴围压的施加方式直接影响破坏模式。PFC 里常见有三种方案，各有利弊：

刚性墙伺服的好处是稳定，三轴代码里可以复用双轴写好的伺服函数，把墙改成圆柱壳或矩形边界即可。缺点是墙体太“硬”，试样侧向变形会被限制，高围压下的鼓胀破坏和真实试验有偏差。

柔性膜颗粒更接近室内试验的橡皮膜，但生成时要在试样外围布置一层特殊颗粒，并设置它们与试样颗粒的接触。膜太刚会限制变形，膜太软又起不到围压作用，需要单独标定。周期边界适合做代表体元，能获得更均匀的宏观应力应变关系，但剪切带会被周期性边界“切断”，不适合研究破坏局部化。

5.3 排水条件在PFC5.0里如何近似处理

PFC5.0 如果不接流体耦合模块，默认模拟的是不包含孔隙水压力的情况，可以近似看作排水条件下颗粒重排列和体积变化的自由发展。换句话说，你看到的三轴曲线是“排水强度”的雏形，不是不排水强度。

要模拟不排水，需要额外引入流体耦合、体积约束或恒定体积剪切。最简单的一种近似方法是约束试样总体积不变，通过调整边界位移来保持体积恒定，但这已经超出基础算例范围。建议先把排水方案跑明白，再考虑不排水扩展。排水条件下，体积应变曲线能清楚反映剪胀剪缩特性，这是标定摩擦系数和颗粒形状的重要依据。

6. 单轴、双轴、三轴结果如何反哺参数标定

6.1 三组曲线分别约束哪类微观参数

不同测试结果对微观参数的敏感度不一样。单轴压缩主要约束弹性模量和破坏强度，双轴测试约束围压依赖性和剪胀，三轴测试校核摩擦角和残余强度。

一个实用的对照关系是：

我在标定时不会一次性把所有参数都丢进优化器，而是分阶段：先用单轴曲线固定刚度，再用不同围压的三轴峰值强度固定摩擦角，最后看体积应变曲线微调孔隙率和颗粒形状。每个阶段只调一个目标涉及的参数，出问题了能快速定位。

6.2 破坏模式的量纲检查与网格尺寸影响

离散元结果受颗粒尺寸影响明显。模型宽度与最大粒径之比至少要大于 5，否则剪切带宽度会被颗粒尺寸扭曲，得到的摩擦角可能失真。颗粒越细，结果越接近连续介质解，但计算量成倍增加，所以需要做一次粒径敏感性分析。

观察破坏模式时要有基本判断：单轴无围压情况下出现贯穿剪切带或劈裂是合理的；三轴高围压下出现鼓胀是合理的；如果模型在低围压下就整体飞散，大概率是加载速率过高或接触粘结设置异常。另一点是量纲检查，PFC 中长度、质量、时间单位在导出结果时很容易混，应力单位错误会让标定参数完全不可用。

6.3 代码集的校验流程：先跑通再调参

每套算例在新环境里跑通后，我会强制走一遍校验流程。第一步，载入模型后不加载，让试样自然平衡，检查是否存在过大的初始应力。第二步，做一次小应变加载再卸载，确认弹性阶段可逆。第三步，跑完整破坏过程，观察破坏模式是否符合该围压下的经验规律。

这三个检查都通过之后，才进入参数标定阶段。如果第一或第二步就异常，不要急着调参数，大概率是模型生成阶段有问题，比如孔隙率不对、悬浮颗粒太多、接触模型设置遗漏。这套校验流程不复杂，但能帮你省掉最想骂人的那些调试夜晚。

7. 我整理这套代码集时踩过的五个典型坑

7.1 试样生成完，孔隙率却对不上

ball distribute指定的孔隙率只是生成时刻的几何孔隙率，平衡过程中重叠颗粒被推开，实际孔隙率会变。后来我发现正确做法是生成后用测量区域量测孔隙率，再把悬浮颗粒筛掉，重新平衡，必要时微调颗粒半径。还有一个技巧：先用稍小半径生成，平衡后再放大到目标半径，孔隙率控制更准，但要注意颗粒重叠过大时反而会产生巨大初应力。

7.2 伺服墙振动，围压曲线像心电图

围压曲线振荡通常就是增益系数过大。降低增益是最直接的解决办法，但很多人不知道降到多少合适。我的经验是先通过墙体接触刚度和接触数量估算参考增益数量级，再乘 0.1 起步，逐步增大到收敛速度和稳定性平衡点。每一步调试都要同时记录墙速和围压曲线，而不是只看最终强度。

7.3 加载速率太快，峰值强度虚高

显式求解最怕惯性效应。压板速度设太快，试样内部还没来得及重排列，峰值就被惯性力抬上去了。除了看动能占比，还要看破坏面的形态。惯性效应明显时，破坏往往从加载端开始，而不是在试样最薄弱部位自然发展。降低速度后重新跑一遍，如果峰值强度明显下降，说明之前的速率不可信。

7.4 输出文件太多，后处理脚本找不到数据

我早期吃过亏：每个循环步都输出颗粒坐标，结果一个三轴算例跑出几十 GB 数据，后处理脚本根本读不动。现在代码集里统一规定，墙的接触力和位移按历史曲线输出，颗粒坐标只在特定阶段输出。文件命名也要规范，比如triax_100kPa_top_wall.csv，这样批量处理不同围压时不会混淆。

7.5 参数换了，结果却“无法解释”

如果你只改了一个参数，结果却出现大幅变化，先别急着判断参数敏感，十有八九是随机种子不同。PFC 的颗粒生成带随机性，固定随机种子才能保证同一参数下结果可复现。做敏感性分析时，每个方案至少跑三到五个不同随机种子取平均，才不会把数值噪声当成规律。

最后再分享一个小习惯：每次跑三轴前，我会把上一轮作用于墙体的速度清零，再把所有墙接触力清零后重新计算，而不是直接在上一轮状态上继续加载。这个操作很简单，但能避免伺服计算把上一轮的残余速度当成当前应力误差，省掉一整晚的 debug 时间。整套代码集看似只是把单轴、双轴、三轴三个测试排在一起，真正值钱的其实是那些在边界条件和伺服控制里反复验证过的判断顺序。你把它换成自己的本构参数后，会发现离散元参数标定不再是玄学。