做地热模拟的同行应该都遇到过这个场景:项目要求评估一座增强型地热开采(EGS)系统的产能与热突破时间,但注水、换热、采热整个过程在地下空间里牵涉到渗流、传热和井筒流动三个不同的物理过程。刚上手时我在 COMSOL Multiphysics 里其实纠结了很久,直到确定用达西定律接口算储层渗流、多孔介质传热接口算岩石骨架与流体之间的热交换、非等温管道流接口算井筒内沿程流动与散热,三个接口通过速度场和温度场耦合起来做瞬态求解,整个模型才顺起来。这套组合非常适合模拟“低温水注入–热储加热–高温水产出”的完整热循环,能回答采出温度何时下降、热提取功率随年份如何衰减、注采压差是否在设计范围内这些关键问题。下面把建模思路、参数取值范围、网格与求解器设置以及我在调试中踩过的坑完整整理出来,供做地热、地下储能和废弃矿井余热利用的模拟工程师直接参考。
1. 项目整体设计与物理图景
1.1 增强型地热系统的闭环换热逻辑
增强型地热开采和传统水热型地热不一样,源岩体通常是比较致密的花岗岩、砂岩,天然渗透率很低,必须先通过水力压裂在储层里造出人工裂缝网络。模型里整个系统是这么闭环工作的:地表注水泵把低温水通过注入井送下去,水在井筒内沿程流动消耗掉一部分压头,到达井底后进入裂缝网络形成的等效多孔介质区域,沿着高渗透带一边向前推进一边和高温岩石发生对流换热,水被加热后从采出井的井底重新进入井筒,最后以高温流体的形式到达地表。
整个过程看似简单,但温度场、压力场、速度场互相耦合,单纯用某一类接口很难完整刻画。我在建模型时习惯先画一张物理过程表:储层内部对应达西渗流,能量交换对应多孔介质传热,井筒内对应管道流。这张表看起来简单,却能直接决定边界条件加在哪一段,后续调试也能少走很多弯路。
1.2 三个接口的分工与耦合逻辑:达西定律、传热与管道流怎么配合
达西定律接口处理的是热储段的低速渗流。与井筒流不同,储层内的流速通常是毫米到厘米每秒量级,雷诺数很小,惯性效应可以忽略,渗流速度与压力梯度之间呈线性关系,这正好落在达西定律的适用范围里。在模型里我一般用渗透率和孔隙率来表征储层,渗透率决定流量分配,孔隙率控制流体的热容量占比。
多孔介质传热接口处理的是水与岩石骨架之间的热交换,默认采用局域热平衡假设,即同一个位置岩石与水的温度相同。这个假设对裂缝发育的EGS储层是合理的,因为单位体积内的换热面积非常大,流固温度差极小。能量方程里的等效热容和等效导热系数是按孔隙率加权得到的,对流项用的就是达西速度场。
非等温管道流接口处理的则是井筒。这个接口把井筒简化成一维线单元,内部有质量守恒、动量守恒和能量方程,能算出沿程摩擦压降、由于温差带来的密度变化和管壁与周围地层的散热。它和储层的连接点就在井底,数据传递方式是:管道末端给出流入储层的流量和温度,储层边界返回压力和温度给管道。用一个不严谨但容易理解的类比:储层是一整块吸饱热量的大海绵,井筒是插进海绵的吸管,三个接口各管一段,最后在井底交界面上通过流量和压力“握手”。
1.3 为什么必须做瞬态模拟而非稳态求解
EGS方案论证阶段最核心的指标不是初始采出温度,而是热突破时间。所谓热突破,就是注入的冷锋面经过若干年推进之后到达采出井,采出温度开始明显下降的现象。这个时间跟井距、流量、渗透率、热容都有关,可能发生在第10年也可能在第50年,它天然是一个随时间演化的过程,稳态模拟完全回答不了。
从时间尺度上看,注水泵启停造成的是小时级的压力扰动,热锋面在储层里的推进却是年量级的事件,瞬态求解要跨越非常大的时间跨度。好在这类模型中温度场的时间常数远大于压力场,可以在求解器里用分阶段时间步进,前面用较小步长让近井温度场建立起来,后面放大步长直接推进到几十年。
2. 核心细节解析与实操要点
2.1 接口方程与适用前提:局域热平衡假设能省多少事
达西定律接口的核心方程是质量守恒方程和达西速度表达式:速度 u = -(k/μ)(∇p - ρg∇D),其中k是渗透率,μ是动力粘度,D是深度坐标。注意这个速度是体积平均的达西速度,不是裂隙内的真实流速,真实流速还要再除以孔隙率。热储流动数值模拟中,压力边界和流量边界经常换算出错,根源就是速度定义和单位没有统一。
多孔介质传热接口的能量方程可以写成 (ρCp)eff ∂T/∂t + ρf Cp,f u·∇T = ∇·(keff∇T)。右边是热传导,左边第二项是对流项,这一项把达西定律接口算出来的速度场当作已知量,所以两个接口必须在同一个模型里按顺序耦合。如果软件版本提供了“多孔介质流动+多孔介质传热”的预置多物理场耦合,建议直接勾选,避免手动漏掉对流项。
非等温管道流接口则是把井筒沿轴线积分得到的一维模型,方程里包含沿程摩擦压降、动能变化和与管壁环境的热交换。对这个接口,我更关注边界条件设置,不推荐在井筒上做三维网格,因为井筒半径相比储层尺度太小,三维建模会带来巨大的网格量,而且对结果精度几乎没有提升。
2.2 关键参数怎么取才不翻车:从岩石骨架到水的物性
参数是整个模型最容易出问题的地方。先说岩石骨架,致密花岗岩类储层的密度在2700 kg/m³左右,比热容在790 J/(kg·K)上下,导热系数范围2.5到3.5 W/(m·K)。这些参数相对稳定,可以先用常数。
水的物性千万不要直接取常数。水温从20℃升高到200℃,密度下降约15%,粘度变成原来的1/5左右,这两个变化会直接改变达西流速大小和注采压差。我在模型里通常把密度、粘度、比热容、导热系数做成温度相关的插值表,让求解器根据当地温度自动取值;如果担心非线性太强导致不收敛,可以先用常数跑通整个流程,再逐步替换成温度相关函数。
储层物性参数里最敏感的是渗透率和孔隙率。基质渗透率通常取0.01到1 mD,压裂裂缝带高渗区可取100到500 mD,孔隙率在5%到15%之间。这个范围跨度很大,所以建完模型后第一件事一定是做渗透率敏感性扫描,看看热突破时间随渗透率怎么变化。我把常用的初始参数范围整理成一个表,方便直接套用:
| 参数 | 典型取值范围 | 备注 |
|---|---|---|
| 储层深度 | 3000-5000 m | 决定初始地温 |
| 地温梯度 | 0.025-0.04 ℃/m | 乘深度得到初始温度 |
| 初始地温 | 100-200 ℃ | 具体由深度与梯度计算 |
| 注入温度 | 60-90 ℃ | 防止过低造成热冲击 |
| 基质渗透率 | 0.01-1 mD | 影响冷锋推进速度 |
| 裂缝带渗透率 | 100-500 mD | 主换热通道 |
| 孔隙率 | 0.05-0.15 | 控制热容与储集空间 |
| 注入流量 | 10-50 kg/s | 由井泵能力和储层注入能力确定 |
这个表是初始模型的输入,不是最终方案的数值;真正工程上还要根据现场测井、试井数据进行校准。
2.3 边界条件与初始条件:深度符号和温度单位这类细节
注入井条件我习惯用定质量流量而不是定压力,工程上注水是泵排量控制,定质量流量模拟出来的注入压力天然可以作为注水泵选型的参考。采出井则用定压力,例如给定井底流压或者规定采出井压力等于某个值,这样采出流量会随储层温度、粘度变化自动调整,更接近真实生产动态。
初始温度场用“地表温度+地温梯度×深度”写入,但要注意深度符号。我在三维模型里习惯把顶部设为z=0,向下为负;如果计算式里用了正深度值,初始温度会把顶底颠倒,热对流方向直接反掉,这是一个很容易被忽略的隐蔽问题。
热边界上,侧向边界离井足够远时用热绝缘就可以,因为热锋面在模拟年限内基本推不到边界;底部边界建议设置成定温,代表深部热补给。这个差别对20年时间尺度的结果影响不大,对50年尺度则不可忽略。
2.4 重力、粘度突变与容易被忽略的热效应
如果模型在垂直方向跨越上百米,重力项就必须保留。达西定律接口里如果关闭重力项,静水压力梯度不会被正确算出,模型底部压力会显得偏低,影响初始压力场。打开重力项之后,初始压力场要先用静水压力梯度初始化,否则瞬态刚开始的十几步会因为压力不平衡出现虚假振荡。
水的粘度对温度非常敏感,但很多新人在第一次建模时都忽略了这一点。用一个我实测过的对比例子:同一套参数,粘度用常数时注采压差约6 MPa,改用温度相关粘度后压差变成约9 MPa,差别足以影响注水泵选型结论。
焦耳-汤姆逊效应和机械耗散在这个场景里基本可以忽略,低温水注入压力不高、流速也不快,它们带来的温度变化远不如对流换热显著,通常不必在能量方程里额外加入耗散源项。
3. 实操过程与核心环节实现
3.1 先从二维轴对称模型搭起:几何越简单越好
做三维井网之前一定先把二维轴对称模型跑通。以单注单采的简化井对为例:径向长度取2000 m,垂直厚度取500 m,注入井和采出井分别位于两侧,中间用一条30到50 m宽的窄高渗透带代表主裂缝方向。几何建模用长方形加两条管道流线段就可以,不需要在井筒内部画实体网格。
二维轴对称的好处是第一轮调试可以集中在物理场耦合上而不是几何上。我修改得最频繁的其实是注入流量和渗透率,这些参数在二维模型里一秒就能改完,放到三维模型里每次重算可能就是数小时。
3.2 多物理场耦合搭建的稳定顺序:分阶段启动
第一步先加达西定律,跑一个稳态得到贯穿储层的速度场。第二步打开多孔介质传热,先不接管道流,直接给储层一个固定出口温度探头,用小步长瞬态跑一下让温度分布有初步形态。第三步再把非等温管道流接口加进来,在井底边界用耦合条件把管道流量和达西域连接起来。
这个顺序非常关键。直接一步到位同时开三个物理场,非线性很强,求解器往往在第一页时间步里就反复缩短步长甚至报错。分阶段启动相当于给耦合问题一个良好的初值,而不是让求解器从零开始裸跑。
管道流接口的井底边界数据传递要检查流量方向和单位。软件里的物理量单位一般默认给的是国际单位制(m/s、Pa、W),但注水工程习惯用kg/s、MPa,如果外部导入数据,一定要在设置界面确认单位换算系数没有漏掉。
3.3 网格划分的核心依据:对流佩克莱数
传热、流动耦合模型的网格密度不能靠感觉,最好先用一个无量纲数判断:对流佩克莱数 Pe = ρf Cp,f u Δx / keff。Pe 大于1意味着网格单元内对流输运强于传导,容易产生数值伪振荡;小于1则说明当地以热传导为主,网格可以放宽。实际操作中我一般把高渗透带和井筒附近网格控制在2到5 m,外围区域放大到20到50 m,然后做一次网格无关性检验。
网格无关性检验的操作是:把网格数量翻倍再算一遍,对比前20年采出温度曲线,如果温降差异不超过1到2 K,这个网格就够用了。温度场在低渗透区域推进缓慢,这些区域网格稀疏一点不影响整体产能结论,但会影响温度云图的美观度,后处理时如果觉得云图发毛,可以局部加密几层再输出。
3.4 瞬态求解器时间步与容差:跨30年模拟的步长策略
时间步设置上我推荐两级策略。前1到2年用小步长,比如1到7天,目的是把井筒附近的径向温度剖面容错建立起来;之后进入长期生产期,步长放大到30到90天;如果还要看50年后的趋势,最后阶段可以用半年步长。直接全部用自动时间步也可以,但相对容差建议不要小于1e-3,否则计算量会翻好几倍。
求解器的时间步策略还和压力场的响应速度有关。地下渗流压力场是几分钟到几小时量级的响应,这一部分如果跟着温度场的大时间步走,可能会在每一步内出现瞬态压力误差。稳妥做法是先给压力场一个稳态预计算,再把温度场用瞬态推进,压力场的初始场就是稳态解,这样后续压力波动幅度小,非线性迭代次数也少。
3.5 后处理里真正有用的输出指标
不要只看温度云图,工程报告里最常用的是三类指标。第一类是采出温度随时间变化曲线,热突破时间定义为采出温度相比初始值下降5 ℃的时刻,这个时间与井距和流量有近似线性关系,是敏感性分析的重心。第二类是热提取功率,由管口温度和流量直接计算为 P(t)=ṁ c_p (T_prod−T_inj)。第三类是累积产出热能,对功率曲线做时间积分,相当于把地热储层这个“地下银行”的累计收益算清楚。
压力方面,记录注入井井底压力随时间的变化,可以看到由于储层温度降低、水粘度变化导致的注入压力演化。把这些数据导出来和实际监测值比对,如果趋势一致,模型就可以用于方案外推;如果趋势相反,基本可以认为是渗透率或边界条件方向出了问题。
4. 常见问题与排查技巧实录
4.1 残差不收敛、步长反复缩小的应急处理
最典型的场景:模型能启动但每算几十步就报警,时间步从30天被强制缩到0.1天。这种问题九成来自水物性的强非线性。我试过最稳的流程是:先等温算流动、再单独解温、最后再放开完整耦合;同时把水物性从温度相关函数替换成分段常数。等整个模型用常数物性稳定跑完20年,再切换回插值表,收敛概率会高很多。
除此之外,打开求解器设置里的阻尼或增加迭代次数上限也能勉强把残差压下来,但这是治标不治本。根源还是初值给得太差,或者物性参数在某个温度区间变化太陡,导致雅可比矩阵迭代时反复震荡。
4.2 井筒方向反了导致注抽反号
非等温管道流接口是按线段方向计算流量的。如果几何建模时线段坐标方向画反了,或入口出口定义得跟预期相反,耦合边界会出现“该注水的地方在抽水”的现象,表现为储层压力整体下降、采出井温度不升反降。排查办法是在初始静止状态下先给一个很小的流量跑几步,看速度场箭头:正常情况注入井附近箭头应该是从井壁向外发散,采出井附近箭头应该向井内汇聚。
如果方向反了,直接翻转载荷方向或翻转线段,不要靠改边界条件的符号来硬拗,硬拗容易把压力场符号也带乱。这是我自己第一版模型踩过的坑,后处理压力曲线怎么看都别扭,最后发现就是线段方向的问题。
4.3 温度单位混乱和负温度
COMSOL 默认的绝对温度单位是K,但工程上给的数据往往是摄氏度。我见过最多的问题是把“温度”变量当成摄氏度的差直接代入物性插值,结果密度和粘度全部错位。建议全程统一:方程里的绝对温度用T表示,摄氏温度用T-273.15显式转换;输入物性表时表头写“温度(K),密度(kg/m³)”,这样后处理里再转回摄氏度也方便。
初始地温计算时还要注意深度符号。顶面为0、向下为负时,温度表达式是 T0 减去梯度乘以深度,而不是加上。符号错了整个初始温度场会颠倒,模型注定跑不出合理结果。
4.4 模型规模太大导致计算时间爆炸
三维全井网、多裂缝带模型动不动就是上百万自由度,瞬态跑30年可能要好几天。我的降规模办法优先级是:先用二维轴对称验证物理行为,再改用等效连续介质替代离散裂缝,然后把模拟时间缩短到10年先看趋势,最后才考虑并行求解。每一次缩减都保持原来的结论不变,基本就能判断模型复杂度有没有必要。
流场和温度场可以分区域粗化:井筒附近的加密网格保持精细,远场的网格粗化到不影响热锋面前沿跟踪即可。这种做法在工程方案筛选阶段完全够用,精度损失可以控制在可接受范围。
4.5 灵敏度扫描与快速参数反演
参数反演未必需要复杂算法。我常用的是双参数扫描:固定注入流量,扫描渗透率从1 mD到500 mD,同时扫描孔隙率0.05到0.15,输出热突破时间矩阵。这个矩阵对现场工作非常有价值,如果实测产水温度在第8年开始下降,且渗透率扫描结果里只有某个数量级的渗透率能给出相近热突破时间,就能快速反推有效渗透率量级。
具体操作是启用参数化扫描功能,批量计算20到30个组合,跑完后把热突破时间矩阵画成等值线图。这是整套模型里性价比最高的输出之一,很多工程决策就靠这张图。
4.6 遇到问题先查什么:问题速查表
| 现象 | 可能原因 | 排查/对策 |
|---|---|---|
| 注采压差异常偏大 | 水的粘度取了常数且未随温度更新 | 换温度相关粘度,检查插值表范围 |
| 采出温度持续不升高 | 井筒方向定义反了或流量为负 | 检查速度场箭头方向 |
| 温度场出现锯齿状振荡 | 对流佩克莱数过大,网格过粗 | 加密高渗透带网格,或开启稳定化 |
| 瞬态一开始压力剧烈波动 | 初始压力场不是静水压力 | 用稳态预计算压力分布作为初始场 |
| 热突破时间异常短 | 高渗透带宽度或渗透率取太高 | 检查裂缝带几何尺寸与渗透率 |
| 计算中途报浮动点错误 | 物性表外插到超出范围 | 把插值表温度范围覆盖整个模拟区间 |
我实际排查的顺序是:先看方向,再查初值,然后看网格,最后检查物性。这个顺序能把九成的问题圈在一个小时以内。先看方向是因为箭头最直观;初值问题通常表现为早期振荡;网格问题在中后期才会显现;物性问题则整个时间段的数值都会变得“不对劲”。这样做比一头扎进报错日志里效率高得多。
最后提一句我自己的习惯:建EGS模型时我会把每一次参数修改都记录在一个单独的表格里,特别是渗透率和注入流量这两列。因为这套模型到了后期,最难的不是让它算出来,而是保证自己记得清每个结果对应的是哪一套参数。有了这个习惯,返工成本能少一半。