干地热这行绕不开的一个问题就是热储层里到底长什么样。尤其是增强地热系统EGS地面往地下打井注入冷水靠岩石里的天然裂缝和人工压裂裂缝形成循环换热通道再把高温流体抽上来。听起来简单可一旦要定量预测注采压差、热突破时间、裂缝重新闭合或者扩张问题就全堆到数值模拟上了。我这两年一直在做COMSOL里裂缝地层的THM耦合仿真从最初的等效连续介质模型一路做到离散裂缝模型再到随机复杂裂缝网络中间踩了无数坑。今天这篇就当是给自己作个阶段总结也希望能给正在用COMSOL做地热、干热岩、增强地热系统模拟的朋友省点时间。这个内容适合地质工程、岩石力学、地热开发方向的研究生和工程师尤其是已经开始用COMSOL但发现“纯连续介质热—流—固耦合”根本描述不了裂缝主导流动的人。1. 为什么增强地热系统研究必须做THM耦合模拟1.1 EGS的地下水—热—力耦合本质先理清一个概念。增强地热系统跟传统热水型地热最大的区别在于EGS的储层渗透率太低了必须通过水力压裂或者化学刺激把地下的天然裂缝激活、扩展形成一个人工的裂缝换热网络。这也就是说整个热储的本质就是“裂缝系统”。注入冷水之后裂缝里的流体压力升高岩石骨架被撑开裂缝开度变大冷水把岩石冷却岩石产生热收缩应变应力场重新分布应力一变化裂缝开度又跟着变渗透率也变进而影响流体流动路径和热交换效率。温度场、流场、应力场三者互相咬合这就是THM耦合的全部内涵。如果只算流体流动和温度忽略力学效应最直接的后果就是热突破时间预测不准。裂缝在高应力区可能慢慢闭合注水压力上不去在低应力区可能被拉开注入水直接“短路”到生产井导致产水温度快速下降。这些现象本质上都是力学行为在主导纯流动—传热仿真完全看不到。具体展开THM三场的相互作用我习惯画这样一张“闭环”温度场影响流场水的粘度和密度都随温度变化温度梯度还会产生热自然对流流场影响温度场流体对流是热量输运的主要机制尤其裂缝中流速快对流换热强度远高于基质流场影响应力场孔隙水压力上升导致有效应力下降相当于岩石骨架“变得更软”温度场影响应力场温度降低引起热收缩产生拉应力反之升温产生压应力应力场影响流场和温度场应力调整改变裂缝开度进而改变渗透率渗透率一动流动路径和换热面积全变应力场影响温度场的另一个途径是损伤演化比如岩石破裂这通常涉及非线性力学。所以在地热领域做THM耦合难点不在于怎么把一个单独的达西流动或者传热算出来而在于怎么把这三个物理场在一个包含离散裂缝的几何模型里稳定地迭代求解出来。1.2 三场耦合的数学描述与关键假设在COMSOL里面做THM耦合本质上就是求解一组偏微分方程组。我不打算把教科书公式抄一遍但关键的控制关系和方程选择值得用“人话”快速过一遍这样后面设置模块和边界条件才能有数。热传递部分能量守恒方程骨架的导热 流体的对流换热 内能变化率如果你用COMSOL的“多孔介质传热”模块默认就是局域热平衡假设即流体和骨架在同一位置温度相同。这个假设对低渗透基质和裂缝而言基本成立但在流速极高的大裂缝或者注入井近井附近可能需要考虑局部非热平衡也就是流体温度和岩石温度分别求解。我在实际中一般先按局域热平衡算只有在研究注入流体温度瞬态冲击时才切到双温度模型。流体流动部分基质采用达西定律v -K/μ(H) ∇(pρgz)裂缝采用“立方定律”修正的切向达西流q_f -b³/(12μ) f(开度修正) ∇(p_f)需要注意岩石渗透率不是常数。在力学变形影响下裂缝开度和基质孔隙度都在变。我常用的裂缝开度随有效应力变化关系是负指数形式b b_ref exp(-α σ_eff)σ_eff σ_n - p这里σ_n是作用在裂缝面上的法向正应力p是缝内流体压力α是经验系数。这个公式的好处是物理意义清晰、参数少坏处是它没有考虑剪切滑移引起的开度增加。如果做的是压裂后滑移显著的地块建议用更复杂的BB模型但那种模型参数标定起来很麻烦初期不建议碰。力学部分就是标准的固体力学平衡方程加本构关系∇·σ ρ_rock b 0σ D:(ε-ε_T-ε_p)其中ε_T是热应变ε_p是孔压引起的应变。COMSOL的固体力学模块自带热膨胀和孔隙压力选项直接勾上就行不用自己手动添加。这套方程组的难点在于三场之间的非线性耦合非常强。温度场和压力场变化慢但应力场对应变变化极其敏感而裂缝开度一个微小的变化会让裂缝渗透率产生立方级的变化——开度变大1%渗透率变大3%。这种跨数量级的反馈让数值求解很容易发散。我在后面实操部分会具体说怎么对付。2. 离散裂缝模型为什么等效介质方法在地热研究中不够用2.1 连续介质方法的局限等效渗透率的先天不足做地热模拟最古老也最简单的方法是等效连续介质模型也就是把裂缝和基质放在一个格子里平均出等效渗透率张量、等效孔隙度、等效比表面积然后用标准的达西方程和传热方程求解。这个方法在煤储层、油气藏开发中长期在用处理得当确实高效但对EGS来讲有暗藏的问题第一裂缝是高度各向异性的。一组NE走向的裂缝渗透率在走向上可能是横向上高34个数量级。用等效张量当然可以描述各向异性但真实裂缝网络往往有成百上千条不同产状的缝等效张量根本没法描述这类“沿着某一条大缝快速突破”的流动通道。在EGS里有时候恰恰是那一条高渗透的大裂缝决定了热突破时间等效介质平均完把关键通道抹没了。第二热交换面积。EGS换热效率的核心指标之一是“岩石—流体接触面积”。等效介质模型用的是等效比表面积而裂缝型储层的换热面积完全由裂缝网络决定比如裂缝间距、开度、密度和走向分布。搞不准裂缝几何就永远算不准取热功率。第三力学响应。裂缝的变形破坏行为与基质岩石完全不同裂缝面在剪应力下可以滑移在拉应力下可以张开等效介质用基质弹模统一描述的话根本模拟不了裂缝的开度演化就更谈不上模拟渗透率随应力的变化了。2.2 COMSOL中DFN模型的两类实现策略在COMSOL中做离散裂缝网络Discrete Fracture NetworkDFN核心思路很直接把裂缝当成嵌入岩石基质中的低维实体——三维模型里是二维面二维模型里是线。也就是说研究区域是2D或3D岩石域裂缝是镶嵌在其中的边界fracture interfaces或者嵌入的实体。我就从这个角度说说COMSOL里实现DFN的两条路线。路线一是“边界单元法”把裂缝当作岩石内部的不连续边界用COMSOL的“吸附面”或者更准确说是“接触界面/内部边界”来定义。流体的切向流动沿边界进行法向泄漏leak-off则通过界面条件计算。这种方法的优势是裂缝厚度不需要在网格里真实体现数值稳定性好但对裂缝开度沿走向变化、裂缝交叉处的流量分配模拟不够直观。路线二是“嵌入单元法”利用COMSOL的“嵌入式实体”Embedded Discrete Element子域将裂缝作为单独的2D域或1D域嵌入岩石域中裂缝域拥有自己独立的材料属性、本构方程和未知量。这样裂缝开度可以作为裂缝域的一个额外参数随应力、温度、压力动态更新。这个方法更接近物理本质在THM耦合中推荐使用。在实际建模中我通常用路线二。选定裂缝为2D域在三维模型里是2D面控制方程为裂缝切向达西流动 切向传热 面内变形或者简化成弹簧界面。然后在岩石域里使用标准的达西流动、固体力学和传热方程。两类域在边界上通过通量守恒和温度、压力连续条件耦合。如果你用的是较新版本的COMSOL比如6.0以上达西定律接口里面就自带了“裂缝”特征传热接口里面也有“薄膜/裂缝”特征设置起来比我早期做的时候方便太多。老版本则需要用PDE模块手动定义裂缝面上的控制方程工作量大而且容易出错。2.3 二维模型起步三维验证的路线我做项目时有个明确的策略先在2D剖面模型上搞定物理逻辑再扩展到3D。原因很简单2D模型里裂缝就是线段几何处理简单网格数量少反馈快THM三场耦合的难点在于耦合逻辑而不是几何维度2D剖面可以很直观地查看温度前锋、应力云图和裂缝开度变化调试方便迁移到3D时只需要把裂缝单元从线段换成面片、把边界条件从曲线换成曲面即可。当然2D模型有局限真实裂缝倾角、三向应力状态无法在2D中完全还原。所以如果是学术论文级别的定量研究还是要上3D。但前提是2D已经把物理机理搞清楚、参数灵敏度摸透了。从成本和效率来看先2D后3D是性价比最高的路径。3. 随机复杂裂缝网络从统计分布到几何建模3.1 裂缝发育的统计规律与分布函数既然是“随机复杂裂缝”第一步是把地下裂缝的几何参数用统计分布描述出来。做EGS的人常用的裂缝统计体系大致包括四个参数位置、长度、走向/倾向、开度。位置分布通常假设为均匀随机分布也就是说裂缝中心在研究采集窗口内任意位置出现的概率一样。但对于一些多期构造控制的地区裂缝位置倾向于聚簇此时可以用泊松聚簇过程来模拟。长度分布是整个体系里最关键也最难标定的一个参数。实测资料和很多文献都表明裂缝长度常遵循幂律分布少数长裂缝主导渗流网络。我经常用的采样公式L L_min × (1 - U) ^ (-1/(D_f - 1))式中U是(0,1)均匀随机数D_f是分形维度一般在1.5~2.5之间。D_f越大小裂缝占比越大D_f越接近1长裂缝占比越高。对于EGS模拟建议D_f取1.82.0既有足够数量的短裂缝又能保留十几条贯穿性长裂缝。走向分布用Fisher分布或双正态分布来描述。Fisher分布的集中度参数κ控制裂缝走向的离散程度κ0表示完全随机各向同性κ5~10表示相对定向排列。这个参数直接决定渗透率张量的各向异性程度。我在设置时会结合区域应力方向来定主走向——次生裂缝一般平行于最大主应力方向所以在力学各向异性明显的地区倾向分布的平均方位角和σ可以选得更定向一些。裂缝开度b多假设为对数正态分布。典型范围0.1mm到1mm均值0.5mm对数值标准差0.5左右。但要注意开度与长度往往存在相关性长裂缝扩张机会更多开度更大。我会在开度采样时叠加一个与长度L相关的项ln(b) ln(b0) β ln(L/L_ref) εβ通常取0.3~0.5ε是零均值高斯噪声。3.2 MATLAB生成随机裂缝网络的完整流程我平时都是MATLAB COMSOL组合因为生成和预处理随机裂缝用MATLAB写脚本最灵活。整个流程分四步参数设定、裂缝采样、几何输出、COMSOL导入。参数设定块把研究区范围、裂缝数量、长度分布参数、走向分布参数、开度分布参数全部定义成变量。比如一个1000m×500m的2D研究区Lx 1000; Ly 500; % 研究区尺寸m N_f 80; % 裂缝条数 L_min 30; L_max 250; % 裂缝长度上下限m Df 1.8; % 长度幂律分布分形维 theta0 60; kappa 8; % 走向均值与Fisher集中度deg b_mean 0.5; b_std 0.4; % 开度对数正态均值与标准差mm裂缝采样块生成每条裂缝的端点坐标和开度crack struct(); for i 1:N_f xc rand * Lx; yc rand * Ly; % 长度采样 U rand; L L_min * (1 - U)^(-1/(Df-1)); L min(L, L_max); % 走向采样用高斯近似Fisher theta theta0 randn * (sqrt(2)/kappa * 180/pi); theta_rad deg2rad(theta); % 形成裂缝线段端点 half L/2; dx half*cos(theta_rad); dy half*sin(theta_rad); crack(i).p1 [xc-dx yc-dy]; crack(i).p2 [xcdx ycdy]; % 开度采样对数正态 crack(i).b exp(log(b_mean) b_std*randn) * 1e-3; % 转成m end几何输出块就是把这些线段写成一个文本文件每行记录一条裂缝的两个端点和开度值方便COMSOL读取。也可以直接在LiveLink for MATLAB环境下调用COMSOL的建模指令把每一段线加入几何序列。导入COMSOL时我踩过的坑是直接用多边形导入可能产生大量微小碎屑尤其是两条裂缝交叉点不在整数坐标上时。解决办法是先在MATLAB里做一次“截断处理”把超出研究区边界的线段裁剪到边界内再把交叉点附近的微小线段合并。这块自动处理逻辑我之后单独开篇讲这里先给结论。3.3 裂缝率与网络连通性检验生成裂缝网络后强烈建议做一步“直视”检验。随机裂缝网络最糟糕的情况是裂缝密度太低、没有形成贯穿连通路径这样注入水和生产井之间形不成有效循环通道模拟结果就全是“注不进去、采不出来”。我习惯用一个简单的指标来快速判断裂缝面积密度I——单位面积内裂缝总长度之和以及连通组分析。连通组分析说白了就是把每一根裂缝当作图节点交叉或接触表示连边然后看最大连通子图覆盖了多少条裂缝。如果连通组基本覆盖全网络、并且连接了注入井和生产井的位置才算合格。如果连通性差可以增加裂缝条数或者增大长度分布的D_f让更多长裂缝把孤立的短裂缝拼接起来。这个步骤非常重要很多文章的模拟结果看起来“很漂亮”但细看裂缝网络根本没连通结果的可信度要大打折扣。天然裂缝网络是簇聚且连通的不会是一堆互不相连的散线。4. 实操COMSOL中建立THM-DFN耦合模型全流程4.1 模块选择、材料参数与边界条件设定打开COMSOL时第一步是选择模块。做THM三场耦合需要的模块组合是固体力学Solid Mechanics→ 对应力学场达西定律Darcys Law→ 对应流动场多孔介质传热Heat Transfer in Porous Media→ 对应温度场这三个模块在COMSOL的“多物理场耦合”里可以直接组合成“THERMAL FLOW SOLID”的组合名称不确定也没关系关键是勾选好交叉耦合项热膨胀thermal expansion、孔隙弹性孔压、热导率依赖于孔隙压力的修正有时候可以省略。材料参数是我另一个想特别提醒的地方。别直接拿文献里的参数一套就开跑一定要结合你研究区块的地质特征。下面是我常用的一套花岗岩型EGS储层参数供参考参数值单位岩石密度2700kg/m³岩石弹性模量40GPa岩石泊松比0.25-岩石热导率3.0W/(m·K)岩石比热900J/(kg·K)岩石初始渗透率1e-16m²岩石孔隙度0.02-基质热膨胀系数8e-61/K水密度998kg/m³水粘度按温度相关函数Pa·s水的比热4180J/(kg·K)温度相关的粘度函数我建议不要用常数因为注入冷水和储层热水温差常常超过150℃粘度可以差一个数量级流动阻力差别巨大。最简单的办法是给水的动力粘度设置一个温度表函数μ(T)从COMSOL材料库里调用或者自己插几组数据点。边界条件方面二维模型为例力学边界模型四周设置为法向滚支roller support可以模拟深部岩体不受远场构造应力干扰的条件。如果需要模拟地应力上表面地表方向施加垂直主应力的prestress左右边施加水平主应力。流动边界注入井位置设置为恒定注入压力或者恒定注水流量生产井位置设置为定压开采压力低于注入井5001500Pa模拟生产压降。温度边界岩石域四周设置为绝热或恒定初温注入井处入口温度恒定比如30℃模型初始温度按地温梯度设置比如4000m深处对应180~220℃。4.2 在COMSOL中为裂缝赋予独立特性这是整个建模最核心的一步把MATLAB生成的裂缝线段转成COMSOL中的2D域并赋予其独立于岩石基质的物理属性。我在2D模型中会把所有裂缝线段合并成一条多段线polyline然后在“几何”里将多段线转换成一个或多个“内部边”对象最后用“分割”命令把研究区模型划分为岩石域和裂缝域。新版本的COMSOL里你可以在达西定律接口内添加“裂缝”特征具体操作是在“达西定律”节点下添加“裂缝”子节点选择由裂缝线段组成的边界作为“裂缝”的依附面在“裂缝属性”中填写裂缝开度b、渗透率kf和储水系数S。kf可以用立方定律估算k_f b² / 12开度0.5mm的裂缝k_f大约是2.08×10⁻⁸ m²比基质渗透率高8个数量级流体会优先走裂缝——这正是我们想要的物理。传热部分用“多孔介质传热”接口下的“薄膜/裂缝”特征设置裂缝的热导率和厚度并让裂缝中的热传递同时考虑径向传导与流体沿裂缝方向的对流。这样温度场就能精确展现流体沿裂缝快速输运热水的过程。力学部分如果做全三维弹性力学耦合需要把裂缝面定义为“接触界面”并用界面刚度替代裂缝开度对压力的依赖。但在2D中一个很实用的替代方案是把裂缝域当做一个线弹性各向同性材料其弹性模量显著低于基质泊松比设为0.35这样它的法向压缩响应实际上就模拟了裂缝闭合。这个方法虽然有点“数值技巧化”但收敛性好物理上又可解释很适合工程尺度模拟。4.3 多物理场耦合配置开度—渗透率动态联动光把裂缝画出来、写上参数还不行THM耦合里最灵魂的一步是让裂缝开度随应力动态变化并反馈到渗透率。我在COMSOL里是这样实现的在“变量”中定义一个全局变量或边界相关变量b_eff表示当前有效裂缝开度从固体力学接口提取裂缝处的法向应力sigma_n从达西定律接口提取裂缝处孔隙压力p_f计算有效正应力 σ_eff σ_n - p_f用负指数公式更新开度 b_eff b_ref × exp(-α σ_eff)再把变量b_eff传递给裂缝域的渗透率定义处kf b_eff^2 / 12。这个耦合看着简单实际运算中非常容易出现数值振荡。因为开度变化会影响渗透率渗透率变化会影响压力场压力场反作用应力场应力场再影响开度——这个循环里任一步有滞后或者过度修正都会导致求解器来回摆动。解决振荡问题我总结出三件法宝第一给开度更新加一个“延迟”。准确说是在COMSOL中用“离散化层”把这个变量设置为“前几步值”而不是隐式同步耦合。通常在非线性求解器设置里把“耦合”改成“分离”步骤让温度和压力先行求解然后力学弹性再更新开度完成一个交错迭代。这样虽然单步精度略有损失但整体鲁棒性大幅提升。第二开度变化范围设置上下限。b_min取初始开度的0.1倍b_max取初始开度的5倍防止压力场异常波动导致渗透率跨几个数量级爆表。第三打开COMSOL的全耦合求解器时把非线性残差容差从默认的0.001放宽到0.005或者0.01。地热模型在场变量空间梯度极大的情况下默认容差经常导致假发散人工放宽一个量级后反而更容易得到物理合理的解。4.4 网格划分策略与求解器选择网格划分是THM模型里最影响成败的因素之一。裂缝宽度在亚毫米级岩石体积在公里级跨11个数量级全用细网格直接算是不可能的。我的做法是局部加密多重尺度基质域自由三角网格最大单元尺寸50m最小单元尺寸5m裂缝附近使用边界层网格boundary layer mesh第一层厚度0.05m增长率1.2层数610条注入井和生产井点附近球形或矩形细化区单元尺寸控制在3~5m裂缝交叉点手动设置点网格加细尺寸1~2m。这样的网格策略下一个1000m×500m的二维模型大约产生3~8万个自由度跑起来非常快。三维模型如果照猫画虎自由度大概率会破百万请做好服务器算力的准备。求解器选择上瞬态问题我一般用“分离”求解器。具体配置步骤1达西定律指定阻尼因子0.7步骤2多孔介质传热阻尼因子0.9步骤3固体力学阻尼因子0.5步骤4耦合更新阻尼因子0.3。分离求解器看起来迂回实际比全耦合求解器稳太多。你用全耦合可能会在第一步就告诉你“找不到解”而分离求解器即便提示某一步收敛慢也能逐步把场推过去。时间步长方面建议初始步长0.1天最大步长不超过总模拟时长的1/100如果使用“自适应时间步”务必设置最大步长限制否则求解器会自动跳到几个月一步直接把热突破细节抹平了。4.5 后处理怎么提取热突破曲线与裂缝流量后处理这部分是我在实际工程交付中最看重的。模拟算完不能只导出一张温度云图交差关键要提取那些能量化、能对比的曲线和指标。热突破曲线是EGS的核心。具体操作在“派生值”里创建生产井处的“点平均值”表达式提取产出温度T_prod随时间的变化数据导出成CSV然后画成T_prod-t曲线。热突破时间通常定义为产出温度下降10%~20%的时刻。除了热突破裂缝流量分配也值得查看。在达西定律接口下创建一个边界积分表达式对每一条主干裂缝边界积分Q_fracture ∫ q_f × t ds这里q_f是裂缝中的切向通量t是边界单位切线方向。这样可以得到每一条裂缝的注采贡献占比看清是不是某一条裂缝“霸占”了整个流量。这在优化布井位置时特别有用。应力后处理也不能忽视。在固体力学接口下绘制σ1最大主应力或σn法向应力沿裂缝的分布图。如果发现某条裂缝局部应力突变就要警惕这里可能出现裂缝面塌陷或张拉破坏对应结果也需要修正。5. 常见问题与排查技巧实录5.1 收敛困难从“不收敛”到“假收敛”的三级排查做THM耦合仿真遇到不收敛是常态关键在于区分是哪种原因导致的。第一如果连稳态计算都不收敛优先检查几何和材料参数。最常见的坑是单位不统一。COMSOL默认单位制下渗透率用m²压力用Pa如果你从文献里直接抄了个渗透率的mD数值忘了换算成m²量级差了约9个数量级整个方程组的病态程度会直接拉满。换算关系1mD ≈ 9.87×10⁻¹⁶m²。第二瞬态计算中某个时间步突然不收敛一般是局部物理状态跳变了。比如局部压力低于水的蒸气压或者裂缝开度更新超过上限。此时最常用的排查办法是把该时间步的中间场导出来看看哪里出现异数就去那里检查边界条件和材料参数。第三最难处理的“假收敛”就是求解器报收敛但结果在物理上明显不合理。比如产出温度一直恒定、压力场完全对称但裂缝开度已经变成负值了。这种假收敛多出现在耦合变量更新没有真正反馈到渗透率定义里。检查方法很简单在裂缝域画一个渗漏率或开度的等值线图看它是否与压力、应力场变化一致。如果不一致基本可以确定耦合定义漏了某一条反馈链路。我维护的排查速查表长这样症状可能原因解决办法稳态就不收敛单位制错乱/参数量级偏差用COMSOL内置单位检查工具核实某时间步发散局部压力/温度越界缩小时间步或给场设置上下限压力场锯齿状振荡网格尺寸在裂缝附近不连续加密裂缝边界层网格开度为负应力计算未排除孔压影响检查有效应力计算表达式温度前锋过缓裂缝域传热未激活对流项确认薄膜特征选择包含轴向对流流量守恒误差大边界条件重复定义检查是否在井点处同时施加了压力和流量5.2 裂缝交叉处的网格与流量分配难题裂缝网络到处是交叉点而在交叉点处网格处理不好直接导致局部流场震荡甚至发散。关键在于如何在交叉点保持通量守恒。我的经验是在COMSOL里两条裂缝交叉等于两条边的公共端点。此时如果你用标准的有限元离散两条裂缝之间的流量分配是靠节点通量守恒自动满足的。但问题在于交叉点附近网格需要z字形过渡质量不好造成数值扩散。解决办法是在几何建模时你就把每条裂缝在交叉点处打断保证每条裂缝的端点严格重合然后在网格划分时给所有交叉点位置添加“点网格尺寸”控制通常设为该区域最大网格尺寸的1/10。这样就可以避免由于几何求交时的微小间隙导致通量泄漏。另外提一点如果裂缝数量多到上百条交叉点数量爆炸网格尺寸自动控制容易失效。这时候需要关闭“自动细化网格”功能手动设定全局尺寸比例和局部细化区域把网格设计权抓在自己手里。5.3 从二维到三维迁移三维DFN的特别注意事项三维模型里裂缝从线变面几何变成平面多边形。生成方式与2D类似但需要仔细处理裂缝产状倾向和倾角以及面与面的相交线。三维裂缝网络的网格量是二维的数十倍如果还连带求解瞬态THM硬件门槛非常高。我建议三维模型要做好两个“缩减”模型尺寸缩减先做400m×400m×400m的块体或者实际储层区域的一半裂缝数量缩减只保留长度排名前20%~30%的主裂缝短裂缝用等效介质填充。这样既保持了主裂缝控制的优势流路径又将计算量降至可控范围。还有一点三维模型力学边界建议使用无限单元或增加缓冲区。否则因为压缩体积受模型大小影响应变计算会偏大导致裂缝开度更新偏激进。我在早期三维模型里吃过这个亏热突破时间比单井示踪试验结果提前了很多后来加了两层缓冲区才基本对得上。结尾一些实操后的小体会做过几次完整的THM-DFN地热模拟之后我最大的体会是“宁可参数糙别让反馈断掉”。刚上手时总想把每一个参数都弄得精确无比结果模型耦合链里某一步忘记更新渗透率整个过程全白算。反而是先构建一条“粗糙但完整”的温度—压力—应力—渗透率反馈链再逐步精修局部参数更容易得到可靠结果。另外模型的物理合理性比数值精度更重要——如果产水温度曲线和现场试井数据差得离谱大概率不是求解器精度不够而是裂缝连通网络或者渗透率反馈逻辑出了问题。此时不要去调网格回去审视裂缝几何。最后分享一个小技巧COMSOL的“结果导出”不要只用默认表格建议把裂缝开度、有效主应力、温度、孔隙压力几个关键变量绑定到同一组导出探针时间同步导出。这样你做参数敏感性分析时数据对齐会方便很多后期画图、写论文都能省出大量时间。地热储层数值模拟的路很长希望这篇总结能帮你少踩几个我踩过的坑。