Comsol模拟液氮水力压裂。 热-流-固-损伤耦合模型。 展示结果为损伤分布和热场分布，可以模拟非均匀的热损伤。
液氮压裂这玩意儿最近在油气开采领域挺火的。传统水力压裂用水会产生地层污染，改用液氮不仅环保还能利用低温特性搞事情。不过低温引发的热应力场和岩石损伤耦合起来贼复杂，不多物理场耦合模型还真玩不转。
咱在COMSOL里搭模型时得把四个模块揉在一起：传热、达西流、固体力学，还有自定义的损伤演化方程。关键是要让温度场影响材料属性，损伤变量反过来又改变渗流通道。举个栗子，岩石导热系数参数得写成温度函数：
k_rock = 2.5*(1 - 0.003*(T - 293)) + 0.5*d; //d是损伤变量这段代码实现了导热系数随温度降低而下降，同时损伤导致孔隙率增大又让导热更差。注意温度单位是开尔文，损伤量d在0到1之间变化。
损伤演化方程最让人头秃。我们用了修正的Drucker-Prager准则，把温度应变考虑进去：
double effective_stress = sqrt(1.5)*deviatoricStress; double thermal_strain = alpha*(T_ref - T); damage_rate = (effective_stress > 35e6 + 0.7*thermal_strain) ? 0.02 : 0;这里thermal_strain是热膨胀系数引发的附加应变，当等效应力超过温度修正的阈值时损伤开始累积。注意条件运算符后面跟的损伤速率是经验值，得根据岩石类型调整。
网格划分有讲究。压裂液注入区域要用边界层网格，像洋葱皮那样套个七八层。特别是近井筒位置，网格尺寸得压到5mm以下才能捕捉到温度骤降导致的应力集中。遇到过不收敛的情况？试着把瞬态求解器的初始步长改成1e-5秒，再打开自动牛顿迭代就好了。
仿真结果里能看到明显的蝴蝶状损伤区。有意思的是低温让损伤带宽度比常规压裂窄了40%，但纵向延伸更深。温度场分布更骚气——离井筒3米处温度从20℃骤降到-196℃，但受岩石非均质性影响，温度等值线像狗啃的一样参差不齐。这说明预设的材料参数场必须包含随机分布缺陷，用个随机函数给杨氏模量加点料：
E = 20e9 * (1 + 0.1*random(-0.3,0.3)); //随机波动30%这么一搞，损伤路径就开始魔性分叉了，跟实际岩芯压裂的实验结果对得上号。
后处理时建议同时显示温度梯度和损伤云图。用交叉滤色器功能能让两者叠加显示，明显看到高温度梯度区对应着损伤前沿。导出数据时记得选高斯平滑，不然云图边缘会有马赛克。
