本人985研究生，comsol裂隙岩体热-流-固耦合数值模拟建模 二维模型，地热开采，注超临界二氧化碳采热，注水采热

裂隙岩体里的地热开采就像在石头缝里捞火锅底料——既要保持岩体结构不崩，又要让热量跟着流体跑得欢。今天咱们用COMSOL搞个二维模型，试试注水和注超临界二氧化碳这两把刷子到底哪家强。

先摆弄几何模型。裂隙这玩意儿得用离散裂隙网络（DFN）方法建模，直接在Geometry里用线段画几条随机分布的裂缝。这里有个骚操作：用MATLAB生成裂隙坐标后，直接在COMSOL里调用LiveLink对接。比如生成主裂隙的代码长这样：

theta = rand(10,1)*pi; %随机角度
L = 0.1 + 0.3*rand(10,1); %裂隙长度0.1-0.4m
comsol.model.geom('geom1').feature.create('frac1', 'Line');
comsol.model.geom('geom1').feature('frac1').set('p1', [0.5-L(1)/2*cos(theta(1)), 0.5-L(1)/2*sin(theta(1))]);

这代码在模型中心区域生成10条随机走向的裂隙，注意长度控制在模型尺寸的1/5以内，防止边界效应捣乱。

物理场设置才是重头戏。得同时开固体力学、达西流和非等温流动三个模块。裂隙处理要用到裂缝流接口，关键参数是渗透率各向异性——法向渗透率得比切向低两个数量级。有个坑要注意：在材料属性里得用if语句区分基质和裂隙区域：

if (dom==2) //裂隙域
   k_xx = 1e-12;
   k_yy = 1e-14;
else //岩块域
   k_xx = 1e-16;
   k_yy = 1e-16;
}

超临界CO₂的物性参数得用插值函数处理。推荐调用NIST数据库的数据，做成.txt文件导入。特别是粘度随温度变化曲线，在80-120℃区间会出现反常拐点，这个不处理好计算结果绝对跑偏。

耦合方程方面，岩体变形会影响渗透率，这里用指数型经验公式：

本人985研究生，comsol裂隙岩体热-流-固耦合数值模拟建模 二维模型，地热开采，注超临界二氧化碳采热，注水采热

k = k0 exp(α (σ/σ0 - 1))

在COMSOL里用变量耦合实现，把固体力学模块的应力分量σ作为达西流的输入参数。注意设置双向耦合时步长别太大，建议初始步长控制在1e4秒以内。

求解器设置建议分两步走：先稳态计算原始地应力场，再用瞬态分析开采过程。遇到不收敛时，把固体力学和流动模块解耦计算几次再重新耦合，相当于给模型做热身运动。

最后对比注水和注CO₂的效果，重点看出口温度随时间变化曲线。超临界CO₂由于携热能力强，前三个月采热效率能高出20%，但后期容易引发裂隙闭合导致产能骤降。这时候就需要调整注入压力参数，在模型里试几个方案比选最优值。

模拟完记得验证网格独立性。在裂隙尖端局部加密到0.1mm后发现流速变化小于3%，说明初始网格设置够用了。要是老板问起计算时间，就说工作站跑了8小时——其实偷偷用了简化模型先试算，这个秘籍一般人不告诉。