做多物理场耦合仿真的朋友应该都有这种感觉：单场分析永远好算，一旦把热、流、固串在一起，计算就开始闹脾气。如果还要把损伤加进去，那基本就是数值模拟里的“地狱模式”。但我这两年做下来，恰恰是这种“地狱模式”最出活——THMD模型（热-流-固-损伤耦合模型）在页岩气开采、地热开发、核废料深地质处置、深部巷道稳定性分析这些场景里，已经成了绕不开的分析工具。这篇博文就从Comsol实现角度出发，把我搭THMD模型时踩过的坑、试过的路子、验证过的参数配置，以及为什么这么配置的底层逻辑，一次性写清楚。适合已经会Comsol基础操作、想往多场耦合方向进阶的工程师和研究生参考。
1. 项目概述：THMD模型到底在算什么样的物理过程
1.1 什么是THMD耦合：一场“牵一发动全身”的多场联动
THMD是Thermo-Hydro-Mechanical-Damage的缩写，对应的是温度场、渗流场、应力场、损伤场四个物理过程的交叉耦合。在自然界和工程场景里，这四个过程从来不是独立的：温度变化引起岩体膨胀或收缩，膨胀收缩改变应力状态，应力状态变化导致微裂隙萌生和扩展，裂隙扩展又改变了渗透率和渗流传热特性——这是一个完整的闭环因果链。传统做法是一次只算其中两三个场，热-流-固三场耦合已经被认为“够复杂了”，为什么还要把损伤加进来？
因为损伤往往是整个耦合系统的“放大器”。没有损伤时，岩石的渗透率基本是常数或者仅随应力小幅变化；一旦损伤变量开始增长，渗透率可能成百上千倍地跳升，热量迁移路径彻底改变，应力场也会因为刚度的局部退化而重新分布。少了损伤这一环，等于把系统里最关键的非线性反馈源头给阉割了。
THMD的本质就是把这四个场的控制方程联立起来，在同一套几何模型、同一套时间框架下同步求解。这种模型能揭示一个纯粹的热-流-固模型看不到的现象：比如注冷水诱发热开裂、损伤区形成优势渗流通道、应力集中区随时间推移发生渐进破坏等。
1.2 为什么必须四场耦合，只做热-流-固行不行
先说结论：如果你的研究对象是完整致密岩体，应力远低于破坏强度，温度变化也不剧烈，那么热-流-固三场耦合确实够用，因为损伤根本不会启动，加上损伤项无非是让计算多了个永远为0的变量，徒增成本。
但实际的深部工程几乎不存在“永远安全”的情况。地热井注冷水时井筒附近温差可以达到几十甚至上百摄氏度，热应力叠加在原岩应力和孔隙压力上，局部拉应力很容易超过岩体抗拉强度；页岩储层压裂改造时水力裂缝尖端应力强度因子飙升，损伤过程和流体流动直接耦合；核废料处置库周围膨润土和岩体在热-水-力联合作用下长期演化，损伤累积决定屏障的长期有效性。这些场景如果不把损伤加进去，计算得到的温度场、渗流场和应力场都会过度乐观，实际工程中可能已经发生损伤的区域在模型里却显示“完好”。
另一个更技术性的理由是：损伤会反过来改变热-流-固三场的基本参数。从本构层面看，损伤后的弹性模量、导热系数、渗透率都应该是状态相关变量。只有把损伤变量纳入本构方程，参数退化才有物理依据，后续结果分析才有说服力。
1.3 这套模型的核心应用场景
我实际接触过的THMD模型应用主要分布在几个方向。地热储层开发是重头戏，尤其是增强型地热系统，冷水注入热储层后岩石温度骤降，热开裂形成的裂隙网络直接决定采热效率；页岩气和水合物开采紧随其后，压裂过程本身就是损伤驱动的渗流-力学耦合过程；深部矿山和隧道工程里开挖卸荷引起围岩损伤演化，损伤区渗透率改变后可能带来突水隐患；核废料深地质处置库的长期安全性评估也需要THMD框架，它要回答的不是“某年某月有没有问题”，而是“数百年演化尺度上屏障层会不会失效”。
这四个应用场景的共同特征是：多场强耦合、非线性强、时间尺度跨度大、安全容错空间小。也正是因为这些特征，纯解析解或经验公式根本顶不住，数值模拟成了唯一可行的手段，而Comsol这种多物理场平台比较适合在统一框架下实现耦合逻辑。
2. 模型搭建：从控制方程到Comsol模块选型
2.1 四场控制方程与它们的“互相咬合”关系
作为建立THMD模型的起点，列出四个物理场的基本控制方程并理清场间耦合关系是必要步骤。我在建模时习惯先把每个场的方程写在纸上，再把耦合项单独圈出来，因为这一步直接决定后面Comsol里接口要怎么选、变量要怎么定义。
热传导方程采用含对流项的形式：
ρC_p ∂T/∂t + ρ_f C_f u·∇T = ∇·(λ_eff ∇T) + Q_T这里的核心是热传导与流体对流的耦合。λ_eff是等效导热系数，我一般写成λ_eff = λ_matrix(1-φ) + λ_f φ，其中φ是孔隙率。损伤发生后λ_eff还会随损伤变量d退化，这部分我在4.1节会展开说明。
渗流场用达西定律的扩展形式：
S_s ∂p/∂t + ∇·(-k/μ (∇p + ρ_f g∇D)) = Q_f关键参数是渗透率张量k，它必须写成应力状态和损伤变量的函数。这正是热-流-固-损伤耦合的核心枢纽：压力场p影响有效应力，损伤重新分布应力，应力改变裂隙开度，裂隙开度直接改渗透率。
固体力学用有效应力形式的平衡方程：
∇·σ' + ρ g + F_s = 0 σ' = σ_total - α_B p I其中α_B为比奥系数，I为单位张量。热应变项加到本构方程里：
σ' = D(d) : (ε - ε_th - ε_p) ε_th = α_T(T - T_ref)ID(d)是随损伤退化的刚度矩阵，α_T为热膨胀系数。这一步的物理含义是温度变化产生热应变，而热应变中能与实际变形协调的部分才转化为热应力。
损伤演化方程最灵活，可以用Comsol的系数型偏微分方程接口施加。采用一个简化但物理合理的指数型演化律：
∂d/∂t = A * max(0, (χ - χ_0)/χ_0)^m * (1 - d)其中χ是等效应变或等效应力比，χ_0是损伤阈值，A和m是材料参数。max(0,·)保证损伤不反向愈合，这个“单向性”对THMD稳定性很重要。
耦合关系可以归纳成一张总览表：
| 耦合方向 | 控制变量 | 力学本质 | 在Comsol中的实现方式 |
|---|---|---|---|
| 热→力 | 温度变化ΔT | 热应变量 | 固体力学中添加热膨胀节点 |
| 力→热 | 应力导致的变形 | 热固耦合（焦耳效应弱） | 双向耦合，但常忽略力学产热 |
| 流→力 | 孔隙压力p | 有效应力原理 | 固体力学中添加孔隙压力载荷 |
| 力→流 | 体积应变ε_v | 孔隙度-渗透率变化 | 自定义渗透率表达式 |
| 热→流 | 温度影响μ/ρ | 热浮力驱动对流 | 通过材料属性表达式实现 |
| 损伤→力 | 刚度退化D(d) | 损伤力学本构 | 弹性模量E=(1-d)E_0 |
| 损伤→流 | 渗透率倍增 | 裂隙连通效应 | k=k_0 exp(βd)等表达式 |
| 损伤→热 | 导热路径改变 | 非均质传热 | λ_eff随d的函数 |
2.2 Comsol物理场接口怎么选：我的配置清单
Comsol里实现THMD模型有多条路，不存在唯一正确答案。我在不同项目里试过几种组合，最终沉淀下来的配置方案是：固体力学接口 + 达西定律接口 + 固体传热接口 + 系数型PDE接口。
固体力学接口选择“固体力学(solid)”，研究类型选瞬态或稳态视需求而定。如果是长期蠕变或热冲击瞬态问题，用瞬态；如果是初步静态分析，用稳态。这里要提醒的是，如果选择稳态但你的损伤演化是时间的函数，那么d变量无法收敛，因为稳态求解器找不到稳定点。瞬态能合理反映损伤演化路径，所以做准确定量分析必须用瞬态。
达西定律接口比“自由和多孔介质流动”接口要轻量得多。后者适合孔隙通道中速度场比较复杂的情形，但真让THMD模型里用的不是局部流动细节，而是压力场的空间分布和渗流速度的宏观趋势，所以达西定律正好够用且更容易收敛。
固体传热接口选“固体传热(heat)”即可，注意把热物性写成孔隙率或损伤变量的函数。如果流体流速较高、对流占主导，固体传热接口的“对流”速度场需要手动填到达西速度计算出来的值。
系数型PDE接口在Comsol里写作“系数形式PDE(c)”，我只用它定义一个标量场d。系数设置中时间依赖项系数da设为1，质量系数ma设为0，守恒通量项系数c设为0（损伤扩散可忽略），源项f写成上文的损伤演化表达式。在这里面还可以顺手加几个“特征值”逻辑让损伤变量总是单调递增。
2.3 损伤变量的实现：系数型PDE是最灵活的路子
很多初学者第一反应是想用内建的“损伤”模块，但Comsol新版本虽然有一些断裂力学界面，真正做连续损伤演化时还是系数型PDE最灵活。原因有二：一是损伤演化方程的形式千变万化，内置接口很难覆盖所有本构模型；二是损伤场往往需要和多个场交叉引用，自定义PDE变量在全局变量和表达式引用上没有任何限制。
举例来说，我定义损伤驱动量χ时，可以写成一个基于最大拉应变和等效压应力的组合表达式。这个表达式可以在Comsol的“变量”节点里写成：
chi = sqrt(⟨ε_1⟩² + ⟨ε_2⟩² + ⟨ε_3⟩²)这里的⟨·⟩表示麦考利括号，即只取正值的部分。最大主拉应变代表了拉开型微裂纹的主导驱动力。如果你想让模型更精细，还可以把剪应变和孔隙压力项加进去，例如：
chi = sqrt(⟨ε_1⟩² + ⟨ε_2⟩² + ⟨ε_3⟩²) + a_p * (p / p_ref)后面那个a_p可以理解成一个经验耦合系数。这个系数的标定如果手头没有实验数据，就先用数量级估算，然后在参数敏感性阶段分析它到底影响多大。
损伤阈值χ_0也很关键。在玄武岩、花岗岩、砂岩这些不同岩性上，χ_0数量级从1e-5到2e-4不等。我建议在初期试算阶段不要追求精确标定，先取一个偏小的阈值让损伤“能发生”，把整个耦合逻辑跑通后，再通过参数扫描逐步逼近实验观测值。
注意：损伤变量在边界上不能随意设为零。如果边界固定为d=0，而边界附近应力又很大，那么边界单元会形成人为“损伤屏障”，不真实的应力重分布会随之而来。我一般让所有边界默认零通量Neumann条件，同时在下边界加滚轴约束限制刚体位移即可。
3. 几何、网格与边界条件：数值实现的关键细节
3.1 几何简化原则：能二维不三维，能对称不全场
THMD模型的计算开销往往很大，几何上的第一原则是“能简化就简化”。
以地热注入井模型为例，如果没有强烈的地层各向异性或断層走向差别，完全可以采用二维轴对称模型代替三维全模型。Comsol的“二维轴对称”选项可以让计算量直接降到三维的十分之一以下，而井筒周围径向温度、孔压和应力变化在轴对称假设下物理上是合理的。
第二个简化原则是切分1/2或1/4对称模型。前提是载荷、边界条件和高渗透带分布都对称。这里要小心，如果损伤演化后形成了非对称的裂缝网络，对称假设就不成立了——这种情况建议还是算全模型。
第三点是在远离扰动区的区域可以只保留一个较大的“远场圆”或“远场矩形”，不需要模拟到地层边界。远场尺寸取关注区域特征尺寸的3~5倍即可，然后再在远场边界上施加等效原岩应力、原位温度与给定孔压。
我吃过亏的地方在于把几何切得太小，结果损伤区扩展到了模型边界，后续所有结果都要作废。建议先做一次小的“试跑”，看看损伤带尺寸占模型尺寸的比例，低于30%比较安全。超过这个比例就要扩大几何范围。
3.2 网格策略：损伤区加密与网格敏感性
多场耦合模型里网格疏密的选择直接影响收敛性和计算量。纯弹性材料问题网格要求相对低，但引入损伤后，局部软化会导致强烈的非线性变形，损伤区模拟结果对网格尺寸异常敏感。这不是Comsol特有的问题，而是连续损伤模型的通病。损伤带宽度取决于网格特征长度，网格越粗，损伤带越宽，局部应变越被低估。
我的做法是先用均匀粗网格快速扫描物理场规律，再在损伤预计发生的区域做局部加密。判断损伤预计发生区的经验依据是应力集中部位、温度梯度最大的地方、渗流高压力梯度带。
具体网格参数可以参考这组我常用的设置：损伤区最大单元尺寸设为(1/5~1/10)的特征裂纹带估计宽度，其它区域尺寸放大到3~5倍，并从损伤区边界向外用渐变拉伸过渡。以二维平面应变问题为例，损伤区如果预计宽度是0.2 m，网格尺寸设为0.02~0.04 m就能基本稳定收敛。
一个非常重要的工程细节：在塑性或者损伤软化区，如果网格严重不加密，计算出的损伤变量可能表现出“锯齿状”或“条带状”分布，这通常是网格敏感性的信号。先用同一套几何跑5倍加密和10倍加密，如果损伤区范围和极端值偏差小于10%，说明网格够了。不满足的话就需要引入正则化手段，这一点放到第6节细说。
3.3 边界条件与初始场的设置：别让模型“起跑就歪”
边界条件设置是新手翻车最多的地方。尤其是多孔介质里的力学边界，常常有人忘记区分总应力边界和有效应力边界。
正确做法是：如果模型边界在地层内部，力学边界用法向应力或位移约束；如果边界是自由表面（比如地表或巷道壁面），则作用的是总应力边界σ_total = σ_sur (总外力)，内部本构方程使用有效应力σ' = σ_total - α_B p。在Comsol里直接写总应力边界就能正确考虑孔压的贡献。
温度边界相对简单。要么给固定温度边界，比如井壁恒温T_inj；要么给对流通量边界，比如地表对流散热。这里有一个我反复验证过的设置：注入井壁面做热对流边界比固定温度更现实，对流换热系数取10~100 W/(m²·K)可以模拟出井壁与流体的换热温差。
初始场要防止不协调。初始温度、初始孔压和初始地应力之间不能互相矛盾，否则第一时步会产生很大的“初始扰动”。我的习惯是先用稳态求解“裸算”一次无注入的平衡态，把它作为瞬态问题的初值，这样温度场、孔压场和应力场自动满足控制方程，能有效避免开场震荡。
4. 求解设置与稳定性控制
4.1 求解器选择与耦合策略
THMD模型是刚性方程组的典型代表，因为热扩散、渗流、应力调整的时间常数相差好几个数量级。岩石导热时间尺度以小时到天计，渗流压力扩散可能在几分钟内完成，力学平衡几乎瞬时。直接显式时间推进大概率会让时间步长小到无法承受，所以必须用隐式求解器。
Comsol里有全耦合和分离两种解算策略。全耦合法则把四个场联立为一个大方程组，直接求解雅可比矩阵。好处是同步收敛、物理上严格一致；坏处是矩阵规模巨大、内存占用高、非线性迭代极易发散，尤其在损伤快速软化阶段。分离法则是按顺序分别求解每个场，然后用迭代修正逐次逼近耦合解。对复杂THMD问题，分离法稳定性更好，调整空间更大。
我实践下来的建议是：前期试算和参数扫描阶段使用分离法，并严格限制每个场的内迭代次数在10次以内；对最终模型、确认参数無误后，再用全耦合方法打磨精度。全耦合求解器设置中内存选“直接”求解器（如PARDISO），并在求解器配置里推荐容差设为1e-4，不要一开始就用1e-6，否则计算量暴涨但精度提升有限。
4.2 时间步与非线性控制参数
时间步控制对损伤演化型问题至关重要。损伤演化在应力超过阈值的那一瞬间往往发生局部突变，时间步过大时求解器直接越过这个突变点，要么错过损伤启动时间，要么非线性迭代完全不收敛。
推荐策略是使用自适应步长，最大步长限制在特征演化时间的1/5以下。以地热注入热冲击问题为例，温度骤变特征时间可能只有几小时，初始时间步长取0.01小时，然后让求解器自动调整到1~5小时，全程不超过最大步长。
非线性控制参数方面，我习惯把阻尼因子从默认1降到0.7左右。这会让求解器更早进入线搜索，减少发散概率。同时限制每时步的Newton迭代最大次数为8~10次，如果在最大迭代次数内不收敛就自动缩小步长重算。
4.3 加速收敛的小技巧
几个实操中至关重要的手段值得单独列出。
第一，引入“连续性过渡”技巧。在损伤阈值附近设计一个平滑的过渡函数，让损伤演化从0平滑启动，而不是硬性突跳。比如用tanh函数写一个门槛函数smooth(χ-χ_0)，这样可以大幅减少非线性迭代次数。这个过渡带设置得过宽会改变损伤启动位置和时间，所以宽度不能超过阈值的5%。
第二，使用扫描辅助参数法。如果有多组参数需要组合分析，先固定一组基参数，用“辅助扫描”逐渐逼近目标参数组合。例如损伤阈值从大往小扫描，让模型从无损状态逐步逼近损伤状态，每步起始解都靠近上一步解，这种“延续法/扫描法”反而比直接设最终参数更稳。
第三，回流到网格和几何检查。如果某时步一直不收敛，先不要怀疑求解器参数，检查是不是损伤区部分单元发生严重扭曲变形。在固体力学里设置“变形几何”或使用自适应网格重划分也能挽救部分收敛失败的情况——但自适应网格在损伤PDE上要非常小心，它可能导致损伤场映射到新网格时人为扩散或损失，所以我在高精度定量分析中宁可手动加细网格也不用自动重划分。
5. 结果分析与参数敏感性
5.1 应该提取哪些关键结果
模型跑通之后的重点从“算出来”变成了“怎么看”。我的习惯是直接从以下结果类别抓取工程可量化指标：
温度场关注的是井壁/边界附近的温度骤降区深度和热影响半径，绘制不同时间点的温度径向剖面能直观反映热锋面的推进速度。这个数据直接关联热采效率和油管热应力评价。孔压场则关注超孔隙压力扩展范围，把p/Δp_max做归一化云图能快速看到哪些区域可能出现张性破坏。力学输出重点看有效应力、偏应力和最大主应力的时空演化，尤其是井壁环向应力和径向应力的组合，它能决定井壁稳定性。损伤场要统计损伤区体积分数、最大损伤变量峰值和损伤带的走向，再叠加渗透率放大系数看到优势和危险渗流通道。
另外一个容易忽略的输出是能量守恒诊断。把注入热量、产出热量、热储存变化、力学功和渗流摩擦耗能全部积分出来，如果有超过5%的能量不守恒，说明热边界条件或者流体属性处理上有问题。
5.2 参数敏感性分析：哪些参数最“要命”
THMD模型动辄十几二十个参数，每个参数都做敏感性分析不现实。结合多个项目的试算结果，我按影响程度把参数排了个优先级：
| 参数 | 影响程度 | 敏感性表现 | 建议标定优先级 |
|---|---|---|---|
| 损伤阈值χ_0 | 极高 | 决定损伤起始时间和空间位置 | 必须实验标定 |
| 渗透率损伤放大系数β | 极高 | 决定损伤区是否成为超级渗流通道 | 必须实验标定 |
| 热膨胀系数α_T | 高 | 直接决定热应力幅值 | 首轮即定 |
| 弹性模量E_0 | 中高 | 影响应力分布模式 | 首轮即定 |
| 初始渗透率k_0 | 高 | 决定压力扩散和热量传递路径 | 首轮即定 |
| 比奥系数α_B | 中 | 影响孔压-应力耦合强度 | 次轮标定 |
| 损伤演化系数A/m | 中 | 控制损伤增长速率 | 次轮标定 |
| 导热系数λ | 中 | 影响热锋面推进速度 | 次轮标定 |
| 热对流换热系数h | 中低 | 影响边界温度动态 | 边界实测 |
实际标定过程中最常出现的情况是：χ_0和β配合不当造成损伤区“蔓延成灾”或“完全不发生”。如果你发现模型里损伤区像雪崩一样扩展到了远场，先检查β是否过大或χ_0是否过小。反过来如果怎么提高温度都迟迟不出现损伤，看看χ表达式里的麦考利括号和温度权重是否写错了。
参数敏感性还能和“逆问题”结合。如果现场监测数据中有温度与位移时程，可以通过Comsol的参数估计模块反演损伤阈值。这一做法对地下工程支护优化挺有价值，但反演模型的数量级误差容易失控，还是要以室内试验标定和局部监测数据联合校核为准。
6. 常见问题与排查技巧实录
6.1 损伤不演化的常见原因
这个现象几乎每个THMD新手都会遇到。明明应力很大，但损伤变量始终保持零。原因排查顺序建议按以下逻辑：先确认有效应力场计算结果足够大，再看损伤驱动量的表达式与阈值对比是否正确，最后检查损伤演化源项里的时间依赖系数有没有赋值。
最常见的原因出现在麦考利括号上。很多人写损伤驱动量时忘了只取拉伸/微裂相关分量，而是把压缩应力也乘以权重加了进了χ。如果阈值卡得很高，压缩分量又贡献有限，自然就很难触发。要检查χ云图，看它是否在预计的应力集中区出现了峰值。如果该有峰值的区域没有峰值，那问题出在χ定义或子模型变量链接不上。
损伤演化表达式中如果用了分段函数或if语句，可能在梯度变化处造成非物理转向，求解器会把它当成不连续点而自动跳过。建议用平滑化函数替代硬阈值判断。
6.2 计算不收敛怎么定位
首先要区分“瞬态首步不收敛”和“后程某一时步不收敛”。首步不收敛大概率是初始场不协调或边界条件冲突；后程不收敛则集中在损伤演化或材料退化导致刚度矩阵发生非正定。
定位方法我总结为“分场剥离法”：先固定损伤场d=0，只做热-流-固三场求解，如果能收敛说明问题出在损伤与其它场的耦合上；接下来启用损伤但把β设为零，看看是否因为渗透率倍增引发压力崩溃；最后再逐步恢复参数。
同时可以在求解器日志里观察残差曲线。如果残差在降到1e-3量级后反复震荡但不下降，多半是损伤变量在局部单元出现振荡，需要回到网格加密和正则化策略。
6.3 网格依赖与正则化
连续损伤模型的最大软肋是网格依赖性。同一套参数，网格加密后损伤带变窄、损伤阈值表现可能完全不同。更不可接受的是，损伤区体积分数随网格加密可能显著变化。
工程上接受的几种主流处理方式：一是引入特征长度l_ch，把损伤演化本构写成依赖单元尺寸的表达式（梯度损伤或非局部积分模型）；二是裂纹带模型，认为损伤带宽度是材料常数，当网格单元尺寸变化时相应调整损伤参数使得耗散能守恒；三是在损伤PDE中加入微尺度扩散项，让损伤场存在一个非局部扩散长度，其实等效于在系数型PDE里设一个很小的扩散系数，比如1e-4到1e-3 m²/s，这样损伤带宽度与网格解耦。
我必须提醒一点：扩散项加得过大，损伤会无中生有地扩展到未损伤区，这也是很多人犯了却不知道的错误。建议先把初始扩散系数设为零跑基准案例，再按需要逐步增大，并通过云图检查它有没有改变损伤区的基本形态。
7. 后续扩展方向与个人实践体会
最后聊一点我在多次项目中的实际体会。THMD模型最大的价值不在于“把四个场算出来了”，而在于它能让你看清一个复杂系统里“哪些环节在主导演化”。我做过一个模拟项目X，当时一开始只做热-流-固耦合，预测的渗透率变化才两倍；后来把损伤加进去，渗透率短期提升了近两个数量级，这才和现场监测的趋势对上了。从那以后，凡是涉及开裂和长期演化的题目，我都默认必须考虑损伤。
后续扩展路径上，我认为最值得尝试的方向是把THMD与离散裂隙网络结合。连续损伤能描述微裂纹的均匀演化，但无法准确描述一条宏观裂缝的起裂、转向和剪切滑移。可以考虑在损伤变量超过某个临界值后的单元节点嵌入不连续面，把连续损伤模型升级为“连续-离散耦合”模型。另外把THMD嵌入随机场框架、引入岩石非均质性参数分布，也能得到更可信的概率型预测而非单一确定性结果。
脚本、参数和模型文件如果整理得干净，后续复用价值极高。我建议每个项目都把参数表单独做成一个Excel文件，用Comsol的“全局参数定义”统一引用，并且把每个关键参数的中文备注和参考依据写在注释里。当你三个月后回头看自己的模型时，会感谢当初多写了这几行注释。
如果这篇内容能帮你把THMD模型的架构搭建起来，我的目的就达到了。余下的是工程细节的打磨——那部分要交给时间、实验和现场数据，一步一步来。
