两个月前我准备把一个二维斜裂纹板的循环拉伸算例迁到自编程序里跑结果被网格重划分折磨得够呛每扩展一个增量步就要重新生成网格裂尖附近还要层层加密算出来的扩展路径又对网格取向特别敏感。后来我干脆把目光转向了近场动力学基于近场动力学搭了一个二维疲劳裂纹扩展模型把数值模拟的程序实现从头到尾走了一遍。这篇是这个系列的第一篇先把引言、基础理论和代码框架说清楚给后面的实现文章定个调。如果你正准备入门近场动力学想自己写一套数值模拟程序或者只是对“疲劳裂纹除了Paris公式之外还有没有其他模拟思路”感兴趣这篇应该对胃口。我会尽量用工程语言讲不绕理论圈子涉及公式的地方也会顺手把物理意义和程序里的对应关系交代清楚。后续文章会围绕这个框架逐步展开所以第一篇的骨架尤其重要。1. 这个坑是怎么开的网格法做疲劳裂纹为什么让人抓狂1.1 一个复现算例引发的折腾事情的起因很简单我想复现一张经典的中心斜裂纹板在循环拉伸下的裂纹扩展路径。听起来不算难但真正动手才发现全是麻烦。如果走传统有限元路线裂纹一旦开始扩展几何边界就变了网格必须跟着重画。重画不是最难受的最难受的是每次网格更新都要把上一步的应力、损伤等状态变量从旧网格映射到新网格映射过程本身就有插值误差。裂纹尖端附近单元质量稍微差一点扩展路径立刻开始飘。后来我换成了扩展有限元裂纹路径不必跟随网格边界看起来省事了很多。但扩展有限元要处理裂尖增强函数、水平集更新、裂尖位置追踪遇到复杂路径或者多裂纹汇合时稳定性问题一个接一个。而且无论哪种方法疲劳问题都要反复提取应力强度因子幅值再把Paris定律代进去每一步都依赖断裂力学里的解析解来“指导”裂纹怎么走。大概折腾到第三周我终于意识到疲劳裂纹扩展的难点不只是“断裂准则选哪个”而是整个模拟框架在结构上就不适合处理裂纹这个动态边界。1.2 疲劳问题比单次断裂更麻烦在什么地方如果只是算一次拉伸断裂哪怕网格重画咬咬牙也能跑完。但疲劳问题有三个维度上的麻烦会同时压过来循环次数数量级太大。一个真实的疲劳算例动辄十万、百万次循环不可能一个循环一个循环地做显式时间积分必须把“循环”当成一个可以跳进的伪时间量这本身就是数值格式上的挑战。应力强度因子幅值提取难。Paris公式里最核心的输入是ΔK而ΔK的计算依赖裂纹尖端应力场的精确描述。裂纹路径一旦不是标准直线K场的解析解就不够用了需要数值积分或者J积分麻烦且容易引入误差。萌生和合并问题不好预判。工业构件里疲劳裂纹往往不只一条萌生位置、多裂纹连通、裂纹之间的屏蔽效应都是网格法最不擅长的场景。我在实践中发现这些麻烦的共同根源是传统方法把“裂纹”当成一个需要预先定义、随时追踪的几何对象。只要这个前提不变网格重划分和裂尖奇异性就永远躲不开。1.3 为什么我会转到近场动力学近场动力学的思路在根上就不一样。它把连续介质力学的基本方程从微分形式改成积分形式物质点之间的相互作用不是通过“相邻单元共享节点”而是通过一个有限半径范围内的一堆“键”来传递。裂纹在这个框架里不是边界条件也不是需要追踪的几何体而是键断裂之后的自然结果。所以近场动力学处理疲劳问题时不需要在每个增量步里问“裂纹尖端在哪”只需要在每个键上判断“你还能承受多少次循环”。Paris公式中对应力强度因子幅值的依赖在这里换成键上循环载荷引起的伸长率幅度损伤在每个键上独立累积累积到一定程度键断掉裂纹自然前进。就是这种“把裂纹降维成键的生死”的思路让我下定决心自己写一套二维近场动力学程序。为了直观对比我把三种思路的核心差异整理成了一个表方法裂纹表征裂尖奇异性处理疲劳实现依赖主要痛点传统有限元几何边界随裂纹扩展重画网格奇异单元/加密重画网格状态映射网格依赖强重画成本高扩展有限元水平集隐式描述裂尖增强函数要先定义裂纹位置再算ΔK多裂纹和三维情况复杂近场动力学一组断键的集合无需特殊处理键损伤累积计算量大边界效应需修正2. 近场动力学改了什么一个积分方程重构断裂问题2.1 从微分方程到积分方程传统连续介质力学的运动方程是局部形式的偏微分方程$$ \rho(\boldsymbol{x}),\ddot{\boldsymbol{u}}(\boldsymbol{x},t) \nabla\cdot\boldsymbol{\sigma} \boldsymbol{b}(\boldsymbol{x},t) $$这个方程成立的前提是位移场足够光滑至少可微。但裂纹出现时位移场在裂纹面两侧直接间断微分关系在裂尖附近失去意义。近场动力学的运动方程换成了积分形式$$ \rho(\boldsymbol{x}),\ddot{\boldsymbol{u}}(\boldsymbol{x},t) \int_{H_x} \boldsymbol{f}(\boldsymbol{\eta},\boldsymbol{\xi}),dV_{x} \boldsymbol{b}(\boldsymbol{x},t) $$其中 $\boldsymbol{\xi}\boldsymbol{x}-\boldsymbol{x}$ 是初始相对位置$\boldsymbol{\eta}\boldsymbol{u}-\boldsymbol{u}$ 是相对位移$H_x$ 是以点 $\boldsymbol{x}$ 为中心、以近场范围 $\delta$ 为半径的区域。在二维问题里$H_x$ 就是一个圆盘。物理图像非常直观每个物质点不是只跟紧挨着的邻居打交道而是跟半径 $\delta$ 范围内的所有物质点通过“键”相互作用。键可以拉长、压缩也可以断裂。键断掉后力传递消失微裂纹自然就出现了。因为整个方程里没有对位移场求导所以在不连续处照样有严格定义。这是近场动力学能绕开裂尖奇异性的根本原因。2.2 键基模型与力-拉伸关系近场动力学里最简单的模型是键基模型。两个物质点之间的键力密度可以写成$$ \boldsymbol{f}(\boldsymbol{\eta},\boldsymbol{\xi}) c,s,\frac{\boldsymbol{\xi}\boldsymbol{\eta}}{|\boldsymbol{\xi}\boldsymbol{\eta}|} $$其中 $s$ 是键的伸长率$$ s \frac{|\boldsymbol{\xi}\boldsymbol{\eta}|-|\boldsymbol{\xi}|}{|\boldsymbol{\xi}|} $$常数 $c$ 是键刚度。当键的伸长率超过临界伸长率 $s_c$ 时键发生不可逆断裂之后力密度恒为零。可以把这个模型理解成无数根微弹簧互相拉扯每根弹簧只承受轴向力拉断后就永久失效。用弹簧类比的好处是程序实现特别简单坏处是键基模型里只有一个弹性常数导致经典泊松比被限制二维情况下取 $\nu1/3$ 附近才能和经典弹性力学自洽。2.3 二维键刚度和临界伸长率的推导思路键刚度 $c$ 不是随便取的值它决定了整个离散化之后的宏观弹性模量。确定方法是让近场动力学在均匀变形下的应变能密度和经典弹性力学一致。以二维平面应力问题为例取单位厚度 $h1$材料点体积 $V_i h,\Delta x^2$。考虑等双轴拉伸应变为 $\varepsilon_0$此时所有键的伸长率都等于 $\varepsilon_0$。近场动力学应变能密度积分可以化简为$$ W_{PD} \frac{1}{2}\int_H c,s^2,\xi,dV \frac{1}{3}\pi c h \delta^3 s^2 $$经典平面应力下的等双轴应变能密度是$$ W_{class} \frac{E s^2}{1-\nu} $$让两者相等得到$$ c \frac{3E}{\pi h \delta^3 (1-\nu)} $$由于键基近场动力学限制泊松比 $\nu1/3$代入后经常写成$$ c \frac{9E}{2\pi h \delta^3} $$这个推导一定要自己走一遍。我见过不少人直接从文献里抄公式结果平面应力、平面应变不分厚度乘没乘也搞混程序结果差出5到10倍多半就出在这一步。临界伸长率 $s_c$ 的确定思路类似让单位面积裂纹完全张开所消耗的能量等于键断裂释放的应变能之和。具体表达式会在后面实现篇里专门推导因为二维和三维不同还与 $\delta/\Delta x$ 有关。引言阶段先把它当成一个可标定的材料参数程序调试时可以先取 $0.01$ 量级跑通流程再精确标定。2.4 裂纹在这里是“结果”不是“边界条件”这是近场动力学最颠覆认知的一点。有限元里你要先有裂纹几何然后网格贴上去扩展有限元里你要用水平集描述裂纹位置但在近场动力学里你只需要给一块完整板料、一组材料参数和一个初始损伤状态裂纹从哪萌生、往哪偏转、怎么分叉全部由键的断裂过程自发生成。代价是后处理变麻烦了。网格法可以直接从几何模型里量裂纹长度近场动力学输出的是一堆散落断键必须自己做连通性分析才能提取“裂纹路径”和“裂纹尖端位置”。这个内容我放在后面专门写但第一篇就要把这种思路转变过来近场动力学的模拟结果不是传统意义上的“裂纹面”而是“损伤场”。3. 疲劳模型把循环计数折算成不可逆损伤3.1 经典Paris定律给我们的参考框架材料疲劳领域最广为人知的是Paris公式$$ \frac{da}{dN} C,(\Delta K)^m $$它描述的是宏观裂纹扩展速率与应力强度因子幅值之间的幂律关系。这个公式简洁、好用但隐含了一个重要前提在计算之前你必须已经知道有一条主导裂纹而且知道它的位置和长度还得能用线弹性断裂力学算出 $\Delta K$。一旦遇到裂纹萌生阶段、多条裂纹汇合、路径在空间自由发展这类问题Paris公式作为“宏观指导”就显得力不从心。近场动力学的疲劳模型可以把Paris定律的幂律思想“下放”到键的层面。宏观裂纹的扩展速率不再直接由 $\Delta K$ 决定而是由每个键的循环伸长率历史决定。键断裂的集合构成宏观裂纹宏观扩展率是大量键失效统计后的涌现结果。3.2 模型A键强度随循环下降这种模型最简单直接每个键在循环载荷作用下峰值伸长率 $s_{\max}$ 超过疲劳门槛 $s_{th}$ 时开始累积损伤。损伤增量写成$$ \Delta D A,(s_{\max}-s_{th})^{\beta},\Delta N $$其中 $\Delta N$ 是这一步跨越的循环数$A$、$\beta$ 是材料常数。键的当前临界伸长率随损伤退化$$ s_{c,N} s_{c0},(1-D) $$当 $s_{\max} \ge s_{c,N}$ 时键断裂。这个模型的物理图像很清晰每根微弹簧的强度随着循环次数增加逐渐下降直到某次循环峰值把它拉断。它的优点是参数少、好实现缺点是没有显式计入平均应力、载荷比和载荷顺序效应。第一版程序用这个模型最容易跑通。3.3 模型B疲劳寿命累积S-N曲线加Miner准则另一种思路更像工程疲劳设计里的做法。假设材料有S-N曲线可以读出任一应力水平对应的寿命 $N_f$那么每个循环造成的损伤就是 $1/N_f$。在近场动力学里把应力幅映射成键的循环峰值伸长率 $s_{\max}$寿命函数写成$$ N_f(s_{\max}) C_1,s_{\max}^{-C_2} $$每个循环步的损伤累积为$$ D_{n1} D_n \frac{\Delta N}{N_f(s_{\max})} $$当累计损伤超过1时键断裂。这就是经典的线性累积损伤准则也叫Miner准则在工程上数据来源充分标定起来比模型A更容易。缺点是线性累积不体现载荷顺序效应高载低载之间的先后顺序对寿命的影响在模型里体现不出来。模型基础数据主要参数优点不足键强度退化裂纹扩展速率试验$A,\beta,s_{th}$形式与Paris律相似易实现不能直接继承S-N数据寿命累积S-N曲线$C_1,C_2$ 或表格工程数据好获得线性累积忽略载荷顺序3.4 需要自己拍板的几个参数无论选哪种模型有几个参数是程序实现里必须早做决定的疲劳门槛 $s_{th}$。它对应传统疲劳理论里的 $\Delta K_{th}$循环峰值低于门槛就不产生损伤。第一版可以先取 $0.4\sim0.5$ 倍的 $s_c$后面再标定。循环跳进步长 $\Delta N$。真实循环数不可能逐个模拟必须一次跨越一批循环。这个参数太大会高估损伤太小则计算量失控具体做法我在第五章展开。拉压不对称性。压缩半循环通常对裂纹扩展贡献很小第一版建议只统计正伸长率峰值即键被拉长的最大量。多轴应力状态。二维多轴时应该取每个键沿键方向的伸长率不要直接用材料点的等效应变否则会丢失方向信息。这点在编程时很容易踩坑。4. 程序框架第一篇先把骨架搭起来4.1 语言与库为什么用C配Eigen疲劳模拟要跑大量循环步每一步又要重新求解位移场对性能有硬性要求。Python做原型很方便但到后面一个算例跑几十万次循环根本扛不住。我最后选了C配合Eigen库处理矩阵和线性代数代码可读性和性能比较均衡。用传统有限元做对比近场动力学的刚度矩阵是非局部的带宽明显更宽但仍然是稀疏的。Eigen的稀疏矩阵和张量操作足够应付二维中等规模问题。如果只是验证小算例Python配Numpy也不是不行但后续要加表面修正、连通性分析、参数扫描时编译语言的迭代优势会越来越明显。4.2 数据结构核心设计程序的第一步是定数据结构。我按“材料—节点—键”三层组织struct Material { double E; // 弹性模量 double nu; // 泊松比键基PD取1/3附近 double delta; // 近场范围 horizon double s_c; // 临界伸长率 double s_th; // 疲劳门槛伸长率 double A; // 疲劳模型系数 double beta; // 疲劳模型指数 int model; // 0脆断, 1键强度退化, 2S-N寿命累积 }; struct Node { Eigen::Vector2d x; // 初始位置 Eigen::Vector2d u; // 当前位移 double V0; // 体积 int fixed; // 是否固定边界点 double damage; // 节点损伤标量用于后处理输出 }; struct Bond { int i, j; // 键连接的两个节点编号 double len0; // 初始键长 double stretch; // 当前伸长率 double damage; // 键的疲劳损伤累积 bool broken; // 是否断裂 };这里有一个关键设计疲劳模型的损伤累积放在Bond上不放在Node上。疲劳裂纹的萌生和扩展是“键级别”的事件节点上的damage只是为了输出云图方便从周围断裂键的比例换算出来不直接参与本构计算。把这个搞混后面更新逻辑会很乱。4.3 邻居搜索与键的存储近场动力学离散化的核心操作是找到每个点半径 $\delta$ 内的所有邻居。最简单的做法是双重循环但二维问题只要网格超过几百个点这个开销就没法接受。工程做法是空间分桶把节点按 $\delta$ 尺寸分到网格桶里搜索时只查相邻格子。伪代码如下// 伪代码邻居搜索和建键 for (auto p : nodes) grid.insert(p); for (auto p : nodes) { for (auto q : grid.query(p.x, delta)) { if (p.id q.id) continue; // 对称键只存一次 double dist (q.x - p.x).norm(); if (dist 1e-8 dist delta) { bonds.push_back({p.id, q.id, dist, 0.0, 0.0, false}); } } }注意键只存一份后续求力时对i和j分别累加对称贡献。这个“只存一半但是访问两次”的模式贯穿整个程序务必统一。4.4 主循环准静态求解加循环跳进疲劳模拟的主循环和常规有限元不同不能用真实时间积分去一个循环一个循环地算。我的框架是int main() { buildGeometry(); buildNeighborListAndBonds(); setBoundaryConditions(); int n 0; while (n maxCycle) { // A. 施加峰值载荷求解准静态力平衡 solveLinearSystem(); // B. 更新所有键的当前伸长率 updateBondStretch(); // C. 用疲劳模型累积损伤标记新断键 updateFatigueDamage(deltaN); // D. 输出VTK和后处理数据 writeVTK(n); n deltaN; } }这里的关键决策是第A步用准静态求解而不是显式时间积分。疲劳加载频率很低惯性效应通常可以忽略所以把一次循环的峰值载荷当成准静态力解一个线性方程组就够了。显式积分在断裂问题里虽然能自然处理不连续但受Courant条件限制时间步极小用它跑一百万次循环没有任何实际可行性。每次有键断裂刚度矩阵都会改变。小算例直接用Eigen的SparseLU重新组装重新分解简单稳妥算例规模变大后再换预条件共轭梯度法等迭代求解器。4.5 边界条件用位移控制更稳定疲劳加载有载荷控制和位移控制两种方式。程序实现里我强烈建议先做位移控制在加载边界点施加固定位移幅值而不是在外边界点施加力。原因是裂纹扩展时试件整体刚度不断下降载荷控制下位移会不断增大数值上容易出现发散位移控制更平滑实际疲劳试验也经常用位移幅值控制。具体施加方法很简单把边界点的自由度在求解方程里约束掉加载端的位移按循环峰值设置即可。5. 二维离散化中那些一不留神就翻车的细节5.1 单位厚度、平面应力与键刚度公式的统一二维问题模型默认取单位厚度但是否在键刚度公式里保留厚度不同文献写法差异很大。程序里最容易犯的错就是把“厚度”重复乘两遍键刚度公式里带一个厚度体积积分里又乘一个厚度结果整体刚度偏大。我的习惯是统一按国际单位制写代码并且把键刚度的最终公式明确写成$$ c \frac{3E}{\pi h \delta^3 (1-\nu)} $$然后所有节点体积一律用 $V_i h,\Delta x^2$。如果取 $h1$体积累计就是 $\Delta x^2$键刚度公式里的 $h$ 不要去掉。写完之后用一个均匀拉伸小块做自检10×10节点小板的等效弹性模量与理论值误差在百分之几到十几的范围内算合理差几倍就是单位或厚度出了问题。5.2 mδ/Δx取多大m3的由来近场范围 $\delta$ 与离散间距 $\Delta x$ 的比值 $m\delta/\Delta x$ 是近场动力学里最重要的数值参数之一。二维情况下平均邻居数约为$$ N_{neighbor} \approx \pi m^2 - 1 $$取 $m3$ 时每个点大约有27个邻居。这就是文献里常说的“经验取值”既能体现非局部效应又不至于让计算量爆炸。$m$ 取太小比如 $m1$ 或 $m2$裂纹路径会明显偏向网格方向近场动力学的优势荡然无存$m$ 取太大边界效应范围变大计算量按平方增长。需要注意的是$\delta$ 本身不只是数值参数它代表材料的一个长度尺度所以做收敛性分析时要关注裂纹扩展路径随 $m$ 的变化而不只是宏观应力。5.3 自由表面和边界的刚度损失近场动力学所有积分都是在有限半径范围内做的靠近自由表面的点周围缺少一部分邻居体元积分被截断导致表面附近整体刚度低于内部。这在断裂模拟里是个很尴尬的问题因为裂纹扩展路径会受这种“表面软化”影响甚至被非物理地吸引到边界附近。处理办法有几个层面。第一版程序可以先不做修正但要把分析区域放到离自由边界足够远的地方精细计算时需要做表面修正因子通过均匀变形测试给每个点乘一个刚度补偿系数。初始裂纹的生成方式也要注意近场动力学里预置裂纹通常不是“切除几何面”而是把跨越裂纹面的键删掉或者把它们的损伤直接设为1。这个操作会让裂纹面附近出现更严重的边界效应所以初始裂纹不要取得太短。5.4 循环跳进的稳定性与ΔN取值循环跳进是疲劳近场动力学程序的核心数值策略也是最容易失控的地方。$\Delta N$ 取得太大可能在一次循环步内让整片键的损伤直接冲到1出现“一夜之间裂纹贯穿板件”的非物理现象取得太小又跑不动几十万次循环。我的经验是做两层限制。第一层是控制单步损伤增量上限比如限制 $\Delta D \le 0.05$如果算出来的增量超过这个值就自动把 $\Delta N$ 减半重算。第二层是实时监控每个键的损伤最大值避免累积出负值或超过1。断键较多时还要做子迭代施加峰值载荷、求解、断键、再求解直到没有新的断键出现。这类似于隐式分析里的“平衡态收敛判断”对疲劳裂纹扩展稳定性非常重要。5.5 输出哪些量才够后处理程序输出不能只存位移云图。疲劳裂纹分析需要看的量包括节点损伤标量、断裂键分布、最大伸长率、累计循环数、裂纹长度。最推荐的数据格式是VTK节点位移和损伤作为点数据断裂键可以作为线数据单独输出。用ParaView打开后能看到损伤云图随时间演变。需要注意一点节点damage直接取该节点所有键的断裂比例裂纹面画出来会很破碎这是正常现象不代表算错了。肉眼判断裂纹路径时反而比只画断键线更直观。6. 第二篇开写之前的三个整理动作6.1 参数接口先统一别让算例参数散在代码里疲劳模型涉及的材料参数非常多弹性模量、泊松比、近场范围、临界伸长率、疲劳门槛、疲劳模型系数、指数、循环跳进步长、输出频率。如果都硬编码在main函数里参数一多很快就乱。我建议一开始就做一个JSON输入文件或一个Config结构体把所有参数集中管理。每个算例一个配置文件程序跑完后把参数自动写进结果头文件。这样做的好处是两周后回来看VTK结果时还能对上“这个算例用的到底是哪组参数”。别看这个小动作不起眼疲劳参数多起来后它能救你很多次。6.2 准备一个最小验收算例近场动力学程序不是写完就能信必须有一个最小验收算例。我准备选一个带初始裂纹的单边缺口板中心或单边预置一条裂纹施加循环拉伸。第一阶段的验收目标是单调拉伸下裂纹能按预期方向张开和扩展第二阶段再切到疲劳模型目标是能得到一条合理的a-N曲线也就是“裂纹长度随循环数增长”的曲线。单调拉伸验证可以用来排查键刚度公式、边界条件和求解器的问题疲劳验证用来排查损伤累积和循环跳进的稳定性。只有这两关都过了再去做更复杂的多裂纹或斜裂纹算例。6.3 提前想好怎么识别和量化裂纹路径近场动力学输出的是大量断键而不是一条显式的裂纹线。要做工程分析必须从断键集合里提取裂纹路径、裂纹尖端位置和扩展长度。这一步的必要性经常被低估等算完拿到一坨散点才反应过来。建议在输出结构里额外保留“每个键断裂时的循环数”后面做裂纹扩展动画和a-N曲线时会非常方便。至于裂纹的连通区域分析和尖端识别我会在后续专门写一篇后处理文章但第一版代码就应当考虑到这个需求别把断裂循环数这个信息丢掉。如果让我重新走一遍这个系列我一定会先把表面修正、连通性分析和参数管理这三件事的框架打好再往里面塞疲劳模型。因为它们比加一个损伤模型难改得多。近场动力学的门槛其实不在那个积分方程本身而在离散化的细节和断裂结果的后处理上。第一篇把这些基础理顺后面写求解器、写疲劳累积、写裂纹路径提取就都是沿着一条已经铺好的路往前走了。下一篇就从准静态求解器的组装开始写。