近场动力学这几年在断裂模拟领域的热度一直不低尤其是做疲劳裂纹扩展的人多多少少都动过用它的念头。传统的有限元处理裂纹要么靠网格重划要么靠扩展有限元里的富集函数去“迁就”裂纹路径一旦遇到多裂纹交汇、分叉或者疲劳引起的缓慢扩展前处理和后处理的精力就明显不够用了。近场动力学的做法完全不同它把材料离散成带有长程相互作用的物质点裂纹只是键断裂后自然涌现的结构不需要预先规定路径。我最早接触这个方向是在做二维疲劳裂纹扩展数值模拟时被XFEM的裂纹追踪问题逼到墙角临时转换思路试了近场动力学。这篇文章是这个系列的第一部分先把“为什么要用近场动力学做二维疲劳裂纹”“理论模型怎么变成程序骨架”“起步阶段最容易在哪几个地方翻车”这几件事讲透。如果你正准备用近场动力学做疲劳断裂仿真或者已经在写了但被各种数值问题卡住这篇文章应该对你有用。1. 为什么选近场动力学来啃疲劳裂纹这块硬骨头1.1 连续介质力学的困境与近场动力学的破局经典连续介质力学里位移场是连续的应力应变关系建立在空间导数的概念上。可裂纹一出现位移场在裂尖处直接不连续应力理论上会趋于无穷大这就是所谓的奇异性。有限元要么在裂尖附近加密网格要么引入特殊的奇异单元来逼近这个状态扩展有限元则通过额外的富集函数把不连续“塞回”单元内部省去网格重划的麻烦。可不管是哪条路扩展方向都需要单独判据多裂纹共同扩展时富集区域相互干扰实现难度直接上一个台阶。近场动力学换了个底层视角。它认为材料点不仅与相邻的点发生作用还会与一定距离范围内的所有点发生非局部相互作用。简单理解就是每根“键”都是一根小弹簧材料变形时弹簧伸长超过某个限度就断开。键断开了裂纹就在这里产生路径完全是计算自组织的不需要预设也没有奇异点需要特殊处理。我用一个比较粗糙的类比帮助理解一块渔网一根线被拉断后它原来承担的力会立刻转移到周围的线上周围线继续说断就断破口就沿着受力最不利的方向扩展。近场动力学就是这种“线断网破”的力学版材料被离散成物质点点与点之间的键负责传力键断则裂纹生。1.2 键基近场动力学的数学基础近场动力学的运动方程是积分形式的这是它和经典偏微分方程最大的区别。对于物质点x运动方程写为ρ(x) ü(x, t) ∫_{H_x} f(η, ξ) dV_{x} b(x, t)这里H_x是以x为中心、半径为δ的邻域δ叫近场范围或horizonξ x - x是x与x之间的初始相对位置η u(x, t) - u(x, t)是两点之间的相对位移f是键力密度b是体积力。后面写键基近场动力学时这根键只有拉伸压缩刚度不区分剪切变形因此对泊松比有天然约束。键力密度f的方向沿当前键的方向大小取决于键的伸长率。伸长率的定义是s (|ξ η| - |ξ|) / |ξ|如果ξ η代表当前两点的相对位置那么s就是“被拉长了百分之多少”。键力密度可以写成f c·s·(ξη)/|ξη|其中c就是微模量相当于单位体积的“弹簧刚度”。当s超过临界值时键断裂此后该键不再传力。程序实现的角度看整个模型的核心数据就是两类一类是材料点的位置、速度、加速度和体积另一类是所有键的初始长度、当前伸长率、受力状态和是否已经断裂。后面所有代码结构都围绕这两类数据展开。2. 核心细节解析从理论到可编程的损伤模型2.1 本构力与微模量的工程换算微模量c是近场动力学里最基础的参数它必须与经典弹性模量E对得上。推导的基本思路是让均匀变形下近场动力学的应变能密度与经典弹性理论一致。二维情况下不同应力状态对应不同表达式其中平面应力状态最常用微模量可以写成c 9E / (π·h·δ³)这里h是二维模型的厚度δ是近场范围。要注意这个结果隐含了泊松比被约束为1/3的假设这是键基近场动力学的固有局限二维平面应力对应固定的泊松比。如果你的材料泊松比离1/3太远键基模型输出的裂纹形态可能偏软或偏硬这时候就要考虑键基模型的修正项或者直接上态型近场动力学。临界伸长率s0也需要换算上断裂能Gc物理意义是“键拉伸到多长会断”。二维平面应力下的常用关系大约为s0 sqrt(π·Gc / (9·E·δ))这个公式里的系数在不同文献里可能略有差异因为推导时如何处理面上键的投影、厚度方向补偿不同作者有不同约定。我的建议是不要直接抄公式。拿到参数后先建一个单边缺口平板的小规模算例用已知的断裂载荷和裂纹起始位置做标定反过来校准s0再上疲劳模型。我一开始图省事直接从文献抄了s0结果静载算例里裂纹提前了快一倍就起裂折腾半天才发现是系数推导口径不一致。2.2 疲劳累积损伤模型给每个键装一个寿命计数器静载断裂只要判断s是否超过临界伸长率但疲劳问题不同。结构在远低于静载强度的循环载荷下也能破坏而且破坏发生在几万甚至几百万个循环之后逐循环模拟完全不现实。工程化的近场动力学疲劳模型思路是给每一根键记录“寿命消耗”。常见的做法是引入损伤变量D初始为0。在每一个“代表性载荷循环”里先通过静力计算得到该键承受的最大伸长率s_max再用某种S-N关系估算该键在这个伸长率水平下能存活的循环数N_f(s_max)。每个代表性循环结束后ΔD 1 / N_f(s_max)D D ΔD当D累计到1键断裂正式退出传力。N_f与s_max的关系通常写成幂律形式例如N_f N0·(s_max/s0)^(-β)N0和β是材料疲劳参数需要实验标定。这本质上相当于把Parish型的Paris裂纹扩展律在键尺度上做了离散化。这里有个很关键的操作细节在这个框架下我们不需要真的算几百万个循环而是每隔若干个循环计算一次最大伸长率然后把ΔD乘以一个跳跃倍数。这个跳跃倍数不能贪大我实测下来过大时键的损伤会在两三次更新之间突然集中增长裂纹呈现出锯齿状跳跃路径看起来极不自然。初学时我习惯先把跳跃倍数设为1跑通模型后逐步放大到10、100观察路径变化稳定后再定最终值。除了损伤累积型模型还有另一派基于强度退化的做法每个循环后按键的损伤量降低其剩余强度当剩余强度被某个静载工况突破时键才断裂。这类模型更适合需要同时考虑疲劳与过载交互的场景。我自己的二维疲劳裂纹程序用的是损伤累积型模型代码结构更简单且对等幅疲劳载荷拟合较好。变幅载荷下需要额外引入载荷交互修正那就是后话了。3. 实操过程与核心环节实现3.1 网格离散化与邻域搜索近场动力学的离散化比有限元简单得多不需要单元直接把计算域切成一排排的均匀网格点每个点带一个体积量二维里就是ΔV Δx²·h。硬件友好的做法是用结构化网格点序号从0到N-1排好相邻点间距Δx。近场范围δ一般取Δx的3倍左右。δ太小非局部效果出不来断裂行为退化成类局部模型路径对网格方向的敏感性急剧上升δ太大邻域内键太多计算开销成倍增长。我用δ 3.0·Δx起步效果比较平衡。邻域搜索是程序里第一个坑。朴素做法是双重循环遍历所有点对距离小于δ就建立键复杂度O(N²)一两万个点时还能忍受到十万点就明显变慢。稍微优化一下就是用格子哈希或者叫空间网格桶把空间按δ大小划分成网格只搜索当前点所在格子及周围格子的点。实现不复杂但能轻松提速几十倍。建立键的数据结构时我强烈建议把每个键的初始长度、两端点编号、当前状态提前存成扁平数组而不是用C里那种松散的vector 。一是减少内存碎片二是OpenMP并行遍历时写冲突少三是cache友好。疲劳模型阶段每个键还需要额外存储历史最大伸长率和累计损伤值。3.2 时间积分与循环加载的准静态化处理近场动力学运动方程含加速度项天然适合显式时间积分我用的最多的是速度Verlet格式。显式积分有一个绕不开的限制时间步要满足类似CFL的条件否则数值解会很快失稳爆炸。步长上限跟网格间距和材料波速有关纵波波速c_wave sqrt(E/ρ)要求大概是Δt 0.8·Δx / c_wave在真实材料单位下这个步长通常非常小。以钢材为例E约200GPa密度约7800kg/m³波速能到5000m/s网格间距如果取1mmΔt上限大约1.6e-7秒。疲劳循环动辄是秒级周期真按时间轴模拟算到天荒地老也跑不完。所以工程处理上都往“准静态”方向靠。实际操作是把一个循环的载荷变化压缩成几步到位每一步施加位移增量后就像做静态松弛一样让系统把动能耗散掉再走下一步。耗散的实现方式可以是滞回阻尼、局部阻尼或者速度重置。我个人用得比较顺的是局部阻尼f_damp -α·|f|·sign(v)其中α是一个0到1之间的常数通常取0.5左右v是点的速度。计算时每一步在合力上叠加上去系统很快就能稳定计算效率比真实时间积分高好几个量级。这个做法在岩土和材料动态模拟里也很常用不算近场动力学专属。3.3 核心代码骨架与疲劳更新实现下面是简化版本的核心力学循环示意形式是类C伪代码省略了网格生成和参数读取部分重点表达键的循环逻辑// 键的更新与受力计算 void computeForces() { for (int b 0; b numBonds; b) { Bond bond bonds[b]; if (bond.broken) continue; Vec dx pos[bond.j] - pos[bond.i]; double dist dx.norm(); double s (dist - bond.length0) / bond.length0; if (s bond.criticalStretch) { bond.broken true; // 静载断裂 bond.damage 1.0; continue; } double f micromodulus * s; Vec dir dx / dist; // 力施加到点上按体积加权 force[bond.i] dir * f * volJ; force[bond.j] - dir * f * volI; } }疲劳部分单独做一个函数循环结束之后调用// 疲劳损伤更新每执行一次代表一个代表性载荷循环 void updateFatigue(double jumpFactor) { for (int b 0; b numBonds; b) { Bond bond bonds[b]; if (bond.broken) continue; double sMax bond.historyMaxStretch; // 本循环内最大伸长率 double Nf N0 * pow(sMax / bond.staticCriticalStretch, -beta); if (Nf 0) { bond.broken true; continue; } bond.damage 1.0 / Nf * jumpFactor; if (bond.damage 1.0) { bond.broken true; bond.damage 1.0; } bond.historyMaxStretch 0.0; // 清理历史记录 } }这两段代码距离生产级还有差距比如surface correction、键的有限体积修正、加载边界层的实现都没有展开但核心思想都在了。拿到骨架后先把静载跑通确认无载荷情况下系统静止、单轴拉伸下应力应变和解析解接近然后再打开疲劳更新逻辑。我习惯把整个求解循环写成这种结构每一载荷步内做一次computeForces→速度Verlet更新→检查能量是否稳定一个载荷循环结束后做一次updateFatigue。这样的模块划分清晰后续想换成态型近场动力学也只需要替换computeForces内部的键力表达式整体框架可以复用。4. 常见问题与排查技巧实录4.1 位移场震荡、能量暴涨怎么办近场动力学显式程序最让人崩溃的问题就是“算着算着点全飞了”位移场像爆炸一样。原因高度集中在三个地方。时间步太大是第一嫌疑。把Δt压小到理论限值的0.5倍甚至0.3倍立刻试一次单轴拉伸算例观察能量曲线是否平滑。第二步看阻尼系数α太小耗散不足系统一直震α太大又会把真实变形路径磨平导致裂纹路径偏移。我的调试顺序永远是先开无裂纹的均匀拉伸把能量震荡调到可以接受再去跑带裂纹的模型。还有一个隐蔽原因初始化时刻系统不是静力平衡的。比如加载边界层上的位移直接一步加到目标值波前会在模型里来回反射很多次才平复。解决办法是分多个子步加载每步加一小段位移让系统有时间松弛。做法简单但非常有效。4.2 裂纹路径“长歪”或者提前分叉二维疲劳裂纹模拟里裂纹应该在指定位置起裂、沿预定对称方向扩展。如果路径长歪先别怀疑算法先检查三件事。第一表面修正做了没有。近场范围在边界处被截断边界上的点邻域不完整等效刚度比内部点低很多裂纹往往就从边界“以为发生了破损”的地方先断。这个误差必须通过修正微模量来处理。最简单的做法是给边界附近的键按邻域体积损失做比例放大复杂一点的做法是先做均匀变形标定计算每个点的修正系数。第二δ和网格尺寸的相对关系。δ/Δx比小于3时裂纹对网格方向很敏感45度方向的键和水平方向键分布差异巨大路径容易沿网格锯齿。把比调大一些会有改善但计算量也随之上来。第三临界伸长率标定是否严谨。疲劳模型中每个键的s0直接决定了起裂时序。我在初期曾经因为把s0调得太低导致整个模型提前进入大面积破损裂纹不是一条主裂纹而是碎成一摊裂纹带。后来用静载单点键破坏分析单独标定了s0问题才消失。4.3 计算效率的优化经验疲劳近场动力学比静载更慢因为每个循环都要重新完成一遍完整的静力松弛过程。优化空间主要有三个。邻域关系必须预计算。最忌讳在力计算里现场搜索邻居数据量一大直接卡死。空间桶的构建成本是一次性的后续所有键的遍历都基于预计算的键数组。循环体强烈建议用OpenMP并行。键与键之间相互独立只有力累加存在写冲突。我通常的做法是每个点维护一个力数组遍历键时只读键信息按点累加时允许原子操作或者先用私有数组归约最后统一合并。实测四核并行能到三倍左右的加速。循环跳跃倍数要按需调节。年轻的时候喜欢追求效率一次性跳几百个循环结果裂纹前缘出现强烈锯齿化路径和实验对不上。后来收敛到跳跃倍率不超过疲劳寿命百分之一的范围精度和效率才算平衡。这一点没有万能公式需要针对自己算例做参数敏感性测试。4.4 疲劳参数标定的一点建议疲劳模型里N0和β这两个参数在多数工程资料里都拿不到现成值。我自己的经验是用宏观点蚀数据反推——如果你手头有实验测得的S-N曲线可以用单键算例把相同应力水平下键的寿命曲线拟出来再做平板试样的整体回代。这套流程做下来虽然费时间但比盲目套文献参数可靠得多。标定完成后一定要保存一份参数记录表。近场动力学的参数之间强耦合往往调了一个就要连带调其他参数。没有记录的情况下跑几个算例后就会陷入“参数不知道是什么导致结果不对”的泥潭。我自己的做法是每个模型文件夹里固定放一个param_log.md记录日期、参数值、算例结果和当时的判断带裂纹的疲劳模型反复迭代时极其有用。5. 最后再分享一点个人体会从接触近场动力学到把一个疲劳裂纹算例稳定跑通中间大概有三分之一的时间都在跟数值稳定性和参数标定较劲。这个领域最大的特点是思路直观但细节极多教科书上的漂亮公式和能出结果的工程代码之间隔着大量看似不起眼的修正和调试工作。我给自己定的路径是先搞静载断裂再上单循环加载最后才开疲劳累积每一步都要保留一个“基线算例”用来回归验证。下一篇我会详细展开二维模型的网格生成、加载边界层的实现和裂纹扩展路径的后处理分析如果时间允许还会把断裂能参数标定的具体算例原原本本整理出来。这里写的很多内容都是踩过坑之后的经验总结希望对正在往这个方向试水的朋友有些帮助。