两个月前,我准备把一个二维斜裂纹板的循环拉伸算例迁到自编程序里跑,结果被网格重划分折磨得够呛:每扩展一个增量步就要重新生成网格,裂尖附近还要层层加密,算出来的扩展路径又对网格取向特别敏感。后来我干脆把目光转向了近场动力学,基于近场动力学搭了一个二维疲劳裂纹扩展模型,把数值模拟的程序实现从头到尾走了一遍。这篇是这个系列的第一篇,先把引言、基础理论和代码框架说清楚,给后面的实现文章定个调。
如果你正准备入门近场动力学,想自己写一套数值模拟程序,或者只是对“疲劳裂纹除了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$ 时,键发生不可逆断裂,之后力密度恒为零。可以把这个模型理解成无数根微弹簧互相拉扯,每根弹簧只承受轴向力,拉断后就永久失效。用弹簧类比的好处是程序实现特别简单,坏处是键基模型里只有一个弹性常数,导致经典泊松比被限制,二维情况下取 $\nu=1/3$ 附近才能和经典弹性力学自洽。
2.3 二维键刚度和临界伸长率的推导思路
键刚度 $c$ 不是随便取的值,它决定了整个离散化之后的宏观弹性模量。确定方法是让近场动力学在均匀变形下的应变能密度和经典弹性力学一致。
以二维平面应力问题为例,取单位厚度 $h=1$,材料点体积 $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)} $$
由于键基近场动力学限制泊松比 $\nu=1/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_{n+1} = 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=键强度退化, 2=S-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$。如果取 $h=1$,体积累计就是 $\Delta x^2$;键刚度公式里的 $h$ 不要去掉。写完之后用一个均匀拉伸小块做自检:10×10节点小板的等效弹性模量与理论值误差在百分之几到十几的范围内算合理,差几倍就是单位或厚度出了问题。
5.2 m=δ/Δx取多大,m=3的由来
近场范围 $\delta$ 与离散间距 $\Delta x$ 的比值 $m=\delta/\Delta x$ 是近场动力学里最重要的数值参数之一。二维情况下平均邻居数约为:
$$ N_{neighbor} \approx \pi m^2 - 1 $$
取 $m=3$ 时,每个点大约有27个邻居。这就是文献里常说的“经验取值”:既能体现非局部效应,又不至于让计算量爆炸。
$m$ 取太小,比如 $m=1$ 或 $m=2$,裂纹路径会明显偏向网格方向,近场动力学的优势荡然无存;$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曲线时会非常方便。至于裂纹的连通区域分析和尖端识别,我会在后续专门写一篇后处理文章,但第一版代码就应当考虑到这个需求,别把断裂循环数这个信息丢掉。
如果让我重新走一遍这个系列,我一定会先把表面修正、连通性分析和参数管理这三件事的框架打好,再往里面塞疲劳模型。因为它们比加一个损伤模型难改得多。近场动力学的门槛其实不在那个积分方程本身,而在离散化的细节和断裂结果的后处理上。第一篇把这些基础理顺,后面写求解器、写疲劳累积、写裂纹路径提取,就都是沿着一条已经铺好的路往前走了。下一篇就从准静态求解器的组装开始写。