COMSOL复现液氮致岩石热损伤：三场耦合建模与调试全记录

上个星期接了个复现委托，客户发来一篇关于液氮致岩石热损伤的文献，要求把里面的热-力学-损伤耦合模型在COMSOL里重新搭起来。刚看到需求时觉得不难，无非就是温度场带着应力场跑，再加一个损伤变量进去。真正动手才发现，这个“简单”背后全是细节：文献里三步并两步跳过的求解设置、没给全的材料参数、还有那个几乎一定会把初学者坑死的参考温度。这篇文章就把我从接单到交付的全过程写出来，既包含COMSOL建模的完整思路，也把反复调试中积累的教训一并放里面，希望让后来者少走几趟弯路。

标题里的“复现”两个字，实际上是很微妙的工作。它不像做原创研究那样可以自由发挥，又不像抄作业那样直接把截图抄过来，而是要根据文献的结果图反推模型设置，再用自己的方式把物理过程重新落地。尤其是液氮致裂岩石这种强非线性、多物理场耦合问题，任何一处边界条件没设对，最后图片方向都对但数值全部偏离。这篇博客适合三类人看：一是正在做岩石热损伤方向的研究生，二是接COMSOL仿真代做/协作项目的工程师，三是想系统理解热-力-损伤三场耦合建模逻辑的软件爱好者。

1. 先想清楚液氮对岩石到底做了什么，再去碰软件

1.1 液氮致裂的物理过程：热冲击与“外壳收缩”

液氮温度大约是77 K，也就是零下196℃左右。把它注入井筒后，岩石壁面被瞬间冷却，表面那一层温度骤降，开始强烈收缩。麻烦的是岩体内部还保持初始地温，于是表面在“拼命往里缩”，内部却在“顶住不让缩”，这种变形不协调在壁面附近产生巨大的拉应力。

岩石这种材料有个让人头疼的特点：抗压能力极强，但抗拉能力很弱。花岗岩的抗压强度可能上百兆帕，抗拉强度往往只有几个兆帕到十几兆帕。热冲击产生的拉应力一旦超过这个阈值，壁面附近就会萌生微裂纹。这个画面其实特别好想象——就像把刚烧热的陶瓷碗突然丢进冰水，外壁先收缩、内壁还胀着，温差足够大时碗壁就咔啦一下炸出裂纹。液氮致裂岩体，内里的逻辑完全一样。

理解了这个机理，建模型时就知道核心变量不是“岩石最终有多冷”，而是“近壁面瞬时温度梯度有多大”。这意味着求解必须做瞬态分析，稳态温度场在这里毫无意义。COMSOL里把传热模块设为“固体传热，瞬态”，力学模块算准静态，这样配置是符合物理本质的。

1.2 温度梯度决定损伤位置，而不是“温度有多低”

很多人第一次做这类模型时，容易下意识觉得损伤应该出现在温度最低的区域。但实际算出来会发现，损伤区并不是均匀分布在整体低温区，而是集中在井筒壁面附近一个很薄的环带里。

原因在于热扩散深度。岩石的导热系数通常不高，花岗岩大概在2.5到3.5 W/(m·K)，热扩散率大约1.4×10⁻⁶ m²/s。用热扩散深度公式δ ≈ sqrt(4·a·t)估一下，100秒时冷锋只推进1.2厘米左右，600秒时约2.9厘米。换句话说，绝大部分岩体在计算时间内根本没感受到低温，只有贴壁的薄层在“扛伤害”。这个薄层试图收缩又被外部大块岩体约束，于是拉应力全部集中在表层，损伤自然从这里向外逐层推进。

所以建模的时候不要一开始就画一个巨大的岩石模型，尺寸要根据损伤扩展深度确定。合理的几何模型应保证外围边界距离井壁足够远，至少是预计损伤深度的10到20倍以上，否则边界约束会人为抬高应力值，算出来的损伤区完全是假象。

1.3 损伤演化的正反馈逻辑：刚度退化和应力重分布

热损伤不是一次性事件。微裂纹萌生后，局部岩石的弹性模量会下降，即所谓刚度退化。模量一降，这一小块岩体就“变软”，周围应力会自动重新分配，原本由它承担的载荷转移到邻近区域。如果邻近区域也已经接近抗拉强度极限，就会跟着进入损伤状态。这样一层接一层向外推动，形成损伤带的时间演化过程。

在数学模型里，这种正反馈可以用简单的关系写出来：

• 应力分布由温度场和损伤后的刚度共同决定；
• 当前应力超过强度阈值，损伤变量增加；
• 损伤变量增大导致弹性模量下降，又反过来改变应力状态。

这个回路由温度场、应力场、损伤变量三部分组成，COMSOL里正好可以用三个物理场接口拼接实现。关键的建模选择在于，损伤变量不能简单地作为后处理结果来看，而必须当作真正的状态变量参与求解。后面第三章会详细讲具体操作。

1.4 建模前先把问题分解成“可求解的准静态热力链”

热冲击导致岩石损伤这个问题，虽然看起来复杂，时间尺度却很清晰：传热是秒到分钟量级，而应力波传播是微秒量级。两者的响应速度差了五六个数量级，所以力学问题完全可以按准静态处理，忽略惯性项。这就把原本的全瞬态动力学简化成了“瞬态传热+准静态力学”的弱耦合链。

实际建模顺序可以分三步走：

1. 先单独求解温度场，确认冷锋推进速度、壁面温度变化是否合理；
2. 再把温度场作为热载荷加载到力学模型中，先不考虑损伤，看看最大拉应力出现在哪、量级多大；
3. 最后加入损伤演化，观察损伤区的萌生和扩展。

这个“先分后合”的思路，既是理解问题的好方法，也是后面处理收敛问题时的救命稻草。很多初学COMSOL的朋友一上来就把所有物理场全部耦合在一起，结果第一次求解直接发散，根本分不清是哪个环节出了问题。先解耦再耦合，定位Bug会容易得多。

2. 模型选型与COMSOL接口：三场耦合怎么搭骨架

2.1 几何与轴对称：井筒模型简化原则

液氮注入井筒的场景，几何上具有天然的轴对称性：一个圆筒，液氮在筒内，岩体向外延伸。只要不研究井底拐角处三维效应，完全可以建成二维轴对称模型，也就是在COMSOL里选择“二维轴对称”空间维度，画一个矩形代表岩体剖面，左边边界是井筒壁，右边是远场。

用轴对称模型能把计算量降一个量级，并且云图非常直观，温度和损伤都能用径向切面展示。以我这次做的模型为例，井筒半径取0.05米，岩体径向从0.05米延伸到2米，高度取2米。这个尺寸足够保证边界效应不影响近壁区。如果客户要求三维效果，可以在二维基础上旋转生成，但求解成本会翻数十倍，不建议在调试阶段使用。

边界条件的物理含义要提前想清楚：

• 左边界（井壁）：液氮载荷，体现为固定低温或对流换热；
• 上边界和下边界：在真实地层中是对称面，法向位移给零；
• 右边界（远场）：距井壁足够远，温度保持初始地温，位移可以自由或给法向约束。

这种边界设置与“无限大地层中一个冷源”的物理图像是对应的。切忌在远场边界上加全固定约束，否则热收缩完全被限制，应力状态会严重失真。

2.2 固体传热接口：热载荷不是“边界温度”就是“热流”

COMSOL中“固体传热”接口的核心方程是：

ρCp·∂T/∂t = ∇·(k∇T)

对于纯导热问题，这个接口设置很简单：初始温度设为地温，左右边界设置不同的热边界条件。岩体的热物性参数中，导热系数k和比热容Cp对结果影响最直接，尤其是k决定了冷锋推进速度，必须认真确定。

在复现文献的时候，液氮与壁面的换热方式有两种常见的处理：

• 一是直接指定壁面温度77 K，适合文献使用了“液氮充分浸润、壁面完全达到液氮温度”的假设；
• 二是给一个对流换热系数和流体温度，更适合模拟工程中液氮与岩壁之间存在气膜、换热并非完全充分的情况。

我在实际项目中，第一轮校准时先用固定壁温，因为这样最容易与文献云图对比。等到整体趋势对上了，再做对流换热系数的敏感性分析，看看结果对边界条件的敏感程度。这比一上来就追求“真实边界”更能快速锁定模型。

2.3 固体力学接口：热膨胀与“参考温度”的坑

固体力学接口里，热应变的计算依赖公式：

ε_thermal = α·(T - T_ref)

这里最关键的是T_ref，即热应变为零的参考温度。在岩体热损伤问题里，T_ref应该设置为岩石尚未受扰动的初始地温，而不是0℃或室温。如果这个参数设错，哪怕温差只有几十度，算出来的初始应力就是错的，整个模型从一开始就报废。

我见过不少初学COMSOL的朋友，把材料属性面板的热膨胀系数填好之后，忘了给“参考温度”赋值。软件默认用293.15 K，但在深部地热场景中初始温度可能高达350 K以上，一上来模型内部就有几十兆帕的虚假热应力，损伤变量初始值直接飘红。这类问题表面上看起来是“报错”，本质上全是物理概念不清晰。

力学边界条件方面，井壁是自由表面，没有额外机械力；远场边界约束要够“软”。通常用“法向约束”模拟远端岩体对内部区域的包围作用，而不是用“固定约束”把所有自由度全部锁死。

2.4 损伤变量的实现：域ODE+瞬时损伤驱动量

这是整个模型中最“COMSOL味儿”的部分。损伤变量D是一个0到1之间的标量，表示材料刚度退化的程度。D=0表示完好，D=1表示完全失去承载力。为了让它以状态变量的身份参与瞬态求解，最干净的做法是使用COMSOL的“域常微分方程（ODE）”接口。

定义一个驱动量D_star，表示当前应力状态下材料“应该”达到的损伤程度：

D_star = max(0, 1 - ft / max(σ1, ft))

其中ft是岩石抗拉强度，σ1是最大主应力（注意符号约定：拉为正）。当σ1远大于ft时，D_star趋近1；当σ1不超过ft时，D_star为0。

然后损伤变量D的演化方程写成：

d(D,t)/dt = (max(D_star, D) - D) / τ

这个式子的妙处在于：

• 当D_star大于当前D时，D会向D_star缓慢靠拢，也就是说损伤随时间逐渐累积；
• 当D_star小于当前D时，分子为0，D保持不变，实现损伤不可逆；
• τ是特征松弛时间，控制损伤发展的速度，数值越小损伤越“瞬时”，但数值越难收敛。

弹性模量随损伤退化的关系用：

E_eff = E0·(1 - 0.99·D) + E0·1e-6

后面的微小项是为了防止完全退化导致刚度矩阵奇异。把这个E_eff填入材料的“弹性矩阵”中，就实现了损伤对力学行为的反向作用。

2.5 耦合回路怎么串起来：温度→应力→损伤→刚度→应力

三场耦合的逻辑链可以总结成一句话：

温度场引起热应变，热应变产生应力，应力超过强度阈值后驱动损伤变量增长，损伤变量反过来降低弹性模量，进而改变下一步的应力分布。

COMSOL里实现这个回路不需要非常复杂的耦合设置。传热接口和固体力学接口之间通过“热膨胀”子节点自动耦合；力学接口和损伤ODE之间通过变量引用耦合，也就是应力变量出现在ODE方程里，E_eff变量出现在材料属性里。求解时用“全耦合”模式让所有因变量同步更新，或者用“分离步”手动控制求解顺序。

如果只是做“单向耦合”，即温度→应力→后处理算损伤，那COMSOL是杀鸡用牛刀，任何有限元软件都能做。真正的难点在于“损伤→刚度→应力”这一条反向路径，它让问题从线性变成强非线性，也让求解收敛成为一场战斗。

3. 复现文献时参数和边界条件怎么“考古”

3.1 材料参数从哪来：文献表格、岩性优选、敏感性验证

复现文献的最大难题，往往不是模型搭建，而是参数考古。很多文献写材料属性时很任性，只给“某型号花岗岩”或者“参数见表”，表里却只有密度和弹性模量，比热容、导热系数、热膨胀系数一概不写。这种时候只能靠岩性知识补全。

以中等强度花岗岩为基准，一套典型的参数组合如下：

参数补齐后，不要急着直接跑最终模型。先做一轮敏感性测试：把每个参数上下浮动20%，观察温度场和损伤范围的变化。哪几个参数影响大、哪几个影响小，一目了然。我这次做下来，热膨胀系数和抗拉强度是对损伤区范围影响最大的两个参数，而比热容的影响相对微弱。这些信息在向客户解释“为什么复现结果和文献存在偏差”时非常有用。

3.2 液氮壁面边界：阶跃温度 vs 对流换热 vs 斜坡加载

壁面热载荷怎么加，是复现成败的分水岭。有三种典型做法：

1. 阶跃温度边界：一开始就把井壁设成77 K。这个设置最暴力，物理上与液氮瞬间完全浸润等效，但数值上极不友好。边界温度从300 K瞬间跳到77 K，初始应变极大，存在严重的数值波动，几乎必然导致第一步或前几步时间推进失败。

2. 斜坡加载：在0.5到2秒内让壁面温度线性从初始地温降到77 K。这个策略兼顾物理合理性和数值稳定性。液氮注入后到完全浸润本来就需要一小段时间；即使真实工况中换热极快，这1秒的过渡相对几百秒的总计算时长而言可以忽略，但对数值求解的帮助是巨大的。

3. 对流换热边界：给一个换热系数h和流体温度77 K。这种处理最接近真实物理，但需要知道h的量级。文献如果没给，就要自己按经验设置或做参数扫描。

我在实操中通常先用斜坡加载，等模型完全跑通后再尝试阶跃或对流边界，评估结果差异。先把一个稳定的版本掌握在手里，再去测试更复杂的边界条件，这是调试的基本素养。

3.3 力学边界：为什么“固定约束”最容易把模型算死

如果把远场边界全部固定，岩石受冷收缩时无处变形，热应力会被人为放大到不可收拾的程度。轻则损伤区异常扩大，重则求解器直接发散。不少初学者遇到“中途不收敛”第一反应是调时间步长，实际上根子在边界条件。

正确的做法是：

• 远场边界（径向最外侧）：采用法向约束，让径向位移为零但允许切向滑动；
• 上下两个对称面：法向位移为零；
• 井壁：完全自由，不施加任何位移约束。

这样设置的意义在于模拟无限大地层对目标区域的弹性约束，而不是把岩体“焊”在边界上。我统计过自己处理的多个液氮模型，凡是一跑就发散的，十有八九是边界约束过强；把固定约束改成法向约束后，求解器立刻变得温顺。

3.4 网格策略：边界层只解决“看得准”，不解决“伤得起”

岩石热损伤问题中，温度梯度和应力梯度都集中在近壁区，所以这里必须细化网格。经验做法是：在井壁处布置8到12层边界层网格，第一层厚度控制在0.1毫米量级，增长因子1.2左右。离开壁面后网格逐渐变稀，远场区域甚至可以用较粗的网格。

但要注意一个常见的误区：网格细只是让结果更精确，并不能挽救不合适的边界条件或错误的物理模型。如果远场固定约束产生了虚假高应力，网格再密也只是把一个错误的结果算得更精美而已。网格问题必须在“物理设置正确”的前提下讨论。

建议至少做三组不同网格密度的试算，比较损伤区宽度和最大应力的变化。如果细网格和粗网格的结果差异超过5%，说明网格还不够细，需要继续加密。这一步叫网格无关性验证，是复现工作中能让客户信任的硬指标。

3.5 关于“复现”的诚实工作流：先标定、再对比、后交付

复现文献模型听起来像是“照猫画虎”，但实际操作中，文献往往省略了大量关键细节。参数缺失、边界模糊、图例不清，几乎是常态。一个可靠的复现工作流是：

1. 把文献结果图中可读的特征提取出来：损伤区厚度、冷锋推进距离、损伤萌生时间；
2. 用合理的默认参数搭一个初步模型；
3. 跑通后，调整材料参数和损伤演化速率，让计算结果与文献特征对齐；
4. 记录每次调整对结果的影响，形成参数敏感性报告；
5. 把最终模型和报告一并交付给客户，说明哪些参数是“标定”出来的，而不是文献直接给定的。

这个过程其实很考验专业判断：哪些差异属于合理范围，哪些差异说明模型结构本身有误。完全一致不现实，但“量级正确、趋势吻合、特征位置对应”是可以做到的。

4. 接单调试实录：一个不收敛模型从报错到跑通的全过程

4.1 第一次运行直接发散：时间步长和初始温度跳跃

我给客户搭的第一版模型，第一次求解只算了几步就报错停掉，提示“找不到一致初值”或“时步可能过小”。这是一类非常典型的失败模式，原因几乎可以锁定在两个方面：

一是边界温度阶跃加载。壁面从300 K瞬间跳到77 K，边界上的温度梯度在初始时刻趋向无穷大，力学计算也产生无穷大的热应变率，没有哪个数值求解器能平稳处理这种跳变。

二是时间步长设置不合理。COMSOL自动时间步进在遇到强非线性时，会不断缩小步长试图收敛。如果模型本身存在刚性源项，步长会缩到10⁻⁶秒量级仍不收敛，计算永远推不进去。

解决办法当然是组合拳：用斜坡函数做温度过渡，给初始时间步长一个合理值（我常用1e-3秒），并开启“自动步进”的同时限制最小步长，防止求解器在无效区间空洞循环。

4.2 把全耦合拆开跑：传热先行、力学跟进、损伤最后

全耦合求解虽然看起来“一步到位”，但一旦发散，很难判断是哪条耦合路径出了问题。我的调试策略是先拆开：

第一步，只算纯传热。关闭固体力学和损伤ODE，算温度场。这一步最容易，几乎不可能失败，目的是验证网格、时间和热边界设置是否合理，同时得到完整的温度场演化数据。

第二步，激活固体力学接口，但把损伤演化方程关掉，只算温度和应力场。观察最大拉应力的量级和出现位置，确认是否与文献描述一致。这时候如果应力场异常，问题出在力学边界或热膨胀设置，和损伤没有关系。

第三步，再激活损伤ODE，做完整三场耦合。经过前两步验证，发散原因就只可能集中在损伤演化和刚度退化引起的非线性上，针对性强得多。

这套“由简到繁”的调试思路，适用于所有多物理场耦合问题。看起来多花一点时间，实际上一轮跑通节省的调试时间远大于此。

4.3 损伤变量“突跳”或振荡：正则化与下限保护

三场耦合打开之后，可能遇到新的麻烦：D值在某些区域突然从0跳到接近1，云图上出现一块“爆炸式”损伤区，随后整个求解步长被压得极小，计算卡死。

问题出在损伤驱动的正反馈上。一旦某点损伤发生，E_eff下降，应力重新分配，可能导致邻近区域应力反而上升，进一步加剧损伤，形成雪崩。这在物理上对应宏观裂纹贯通前的不稳定扩展，但数值上必须用正则化手段抑制病态行为。

我的做法有两个：

• 一是给损伤演化加松弛时间τ，不让D瞬间跳到D_star，而是用一阶滞后让损伤有时间发展。τ取0.1到1秒，既能保持物理趋势，又能显著提升稳定性。
• 二是给E_eff加下限保护，不允许模量退化到零，否则刚度矩阵会出现零主元，求解器直接崩掉。前面提到的E_eff = E0·(1 - 0.99·D) + E0·1e-6就是干这个用的。

这两个手段在文献中通常被归入“粘性正则化”或“非局部损伤”的范畴。严格意义上它们改变了原始方程的数学性质，但在薄损伤带的模拟中，这是工程上普遍接受的做法。

4.4 速度对比与结果稳定性：三种求解顺序的实测差异

调试过程中，我对比了三种求解策略的性能：

最终我采用的是第二套策略：开启“分离”求解器，把传热和力学/损伤分为两个组，在每个时间步内先更新温度场，再更新力学和损伤。这样既保留了双向耦合效应，又把非线性控制在可收敛的范围内。

需要注意的是，这种“不完全全耦合”的做法会引入微小的时间滞后误差。验证方法是把时间步长减半，看结果是否基本不变。如果减半后损伤区宽度和温度场曲线变化不大，说明滞后误差可以接受。

4.5 收敛不代表可信：网格无关性和时间步收敛性检查

模型收敛之后，不能立刻交给客户。要先通过两关检查：

第一关是网格无关性验证。把边界层数量和整体网格密度同时提高一档，重新计算，比较损伤区宽度、最大拉应力、损伤变量云图的差异。如果差异在5%以内，说明网格精度足够；如果差异显著，需要继续加密，直到结果稳定。这个检查和前面4.1里的初始调试不是一回事——前者是确保算得动，这里是确保算得准。

第二关是时间步收敛性检查。把最大时间步长缩小一半再算一遍，观察关键变量随时间的变化是否一致。如果温度曲线出现明显锯齿或者损伤变量不光滑，说明时间离散精度不够，需要加密。

这两关都过了，模型才算“可信”。否则只凭“能算完、云图好看”，根本不敢拿出去跟文献对比。

5. 后处理、验证与交付：怎样才算“复现成功”

5.1 客户最关心的三张图：温度云图、第一主应力、损伤因子

仿真项目的交付，最直观的载体永远是图。对液氮致岩石热损伤这个问题，有三张图是客户一定会盯着的：

第一张是温度场云图。核心看点在于冷锋的推进形态：壁面附近应该是紧密排列的等温线，向外逐渐稀疏。这是温度梯度的直接体现，也是判断热边界条件是否合理的重要证据。如果冷锋推进太快，说明导热系数设大了或者几何尺寸太小；如果太慢，则可能是比热容过大。

第二张是最大主应力云图。重点看拉应力集中的位置和量级。理想状态下，拉应力最大值出现在井壁附近，随着温度梯度发展而向内迁移。如果某个角点出现孤立的高应力点，基本可以断定那里存在几何尖角或约束过强，需要回头修几何或边界条件。

第三张是损伤因子D云图。物理上它应该沿井壁形成一圈连续的“损伤壳”，在几厘米到十几厘米范围内衰减到零。如果损伤区出现零星的“孤岛”或夸张的锯齿状边界，往往是收敛不够或者网格质量不行。

这三张图放在一起，温度场说明“冷到哪”，应力场说明“哪受力”，损伤场说明“哪坏了”。逻辑链条完整，客户一看就懂。

5.2 定量对比的指标：冷锋位置、损伤半径、萌生时间

光看图还不够，复现工作最终要靠数字说话。我建议提取三类定量指标与文献对比：

其中损伤萌生时间是最敏感也最能暴露模型错误的指标。如果文献中说损伤在注液后10秒左右出现，而你的模型60秒才开始损伤，那多半是热膨胀系数或抗拉强度标定有问题，云图形状再相似也不能草率交付。

提取这些指标时，COMSOL里可以定义“探针”，在指定点记录D值随时间的变化曲线。也可以直接在结果导出中沿某条边提取数据，存成文本后用表格工具处理。这些数据文件同时是交付物的一部分，客户拿去做后续对比或二次开发都很方便。

5.3 边界条件敏感性测试：给客户讲清楚模型的“权限”

我通常会在交付前做一组敏感性测试，内容围绕两个问题：如果壁面不是完全达到77 K，损伤区会差多少？如果实际地温比预设高，结果将如何变化？

具体操作方式是在合理范围内扫描液氮换热系数和初始地温，把损伤区宽度、损伤萌生时间随参数的变化整理成表格。这么做的意义很大：一方面能帮客户理解模型的适用范围，明确哪些结果依赖假设条件；另一方面也让“复现偏差”有了合理化的解释空间，比如文献可能在实验过程中未能完全达到液氮温度，模拟中按77 K计算会导致损伤区偏大，这是在物理上说得通的。

给客户交付时，我会明确标注哪些参数来自文献原表、哪些来自岩性经验、哪些是标定结果。这个“参数溯源清单”对你的职业信誉很重要，防止模型被二次转发或改参数时发生误用。

5.4 交付物与二次开发：模型还能迁移到哪些液氮/低温场景

最终交付物通常包括：一个干净的COMSOL模型文件、一份参数说明表、一组后处理导出图、以及关键位置的时程曲线数据。如果客户有脚本能力，还可以附上一份详细的建模步骤说明，方便他自行改参数重算。

这套模型框架的扩展空间很大，比如：

• 把液氮换成一氧化碳或冷水，研究不同介质下的热冲击效果；
• 在边界条件中增加地应力场，考察初始应力对损伤形态的影响；
• 把均匀介质换成多层岩体，研究层理面对损伤扩展路径的影响；
• 进一步引入渗流场，模拟液氮气化后的压力致裂与热损伤的联合作用。

这些都是同一个“热-力-损伤”框架上做加法。

5.5 一个容易忽略的小细节：别让结果动画骗了你

每次接到这类项目，我都会额外提醒客户一件事：结果动画里损伤云图的“突然出现”往往不是真实的物理突变，而是损伤变量到达阈值前后的数值表现。如果动画里看到损伤在几十毫秒内从0跳到满值，不要兴奋，那大概率是正则化没做好或者时间步长切得太粗。

判断标准很简单：把最大时间步长缩小10倍重算一次，看损伤扩展过程是否变得平滑。如果大幅变化，说明原计算的时间分辨率不足，动画只是数值“跳变”而非真实演化。

最后再聊一个接单时常用的沟通技巧。做“复现”类项目，第一次给客户看结果图时，不要直接给最完美的版本，先同时展示一版“文献原始图”、一版“初步复现结果”、一版“参数标定后结果”。让客户直观看到复现流程中每一步带来的差异。这一步能极大降低沟通成本，对方不会一根筋地追问“为什么不是一模一样”，而是会理解仿真与实验、文献之间永远有合理的容差带。

液氮致岩石热损伤这个问题，物理上非常精彩，几十兆帕的拉应力能在几秒内把坚硬的岩石撕开，COMSOL把它变成可以逐步调试的多物理场模型后，反而显得克制可靠。整套模型跑通之后，我看着损伤云图上那圈从井壁慢慢往外爬的“裂纹带”，还是很有成就感的。希望这篇记录能给正在和收敛性搏斗的朋友们一点方向上的帮助。