做comsol仿真模拟水力压裂。 岩石变形方程、流体渗流方程、应变方程以及相场扩散方程的求解分别采用4 个指定模块。
打开COMSOL Multiphysics新建模型时，总有人被满屏的物理场接口搞懵。水力压裂这玩意儿，说穿了就是固体变形和流体撕逼的过程。今天咱们拆解下四个核心方程对应的模块怎么玩。
岩石变形交给固体力学模块
先看岩石怎么被压到变形。固体力学模块处理的是应力-应变关系，核心方程长这样：
rho * d2u/dt2 = ∇·σ + F_body σ = C : ε # 本构关系 ε = 0.5*(∇u + (∇u)^T) # 应变张量这里的C其实是弹性矩阵，对应着材料属性里的杨氏模量和泊松比。设置时最容易翻车的是边界条件——记得把模型底部设为固定约束，否则整个地层会像果冻一样乱晃。
流体渗流用达西定律模块
压裂液在裂缝里的流动用达西定律描述：
% 达西速度方程 v = -(k/mu) * (∇p + rho_fluid*g)这里k是渗透率，实际模拟时要用变量重定义技巧——随着裂缝张开，k值会指数级暴增。建议用阶跃函数过渡，避免求解器抽风：
k = k0 + (k_fracture - k0) * flc2hs(phi, 0.5)flc2hs这个函数能平滑处理相场变量phi的突变，比简单if判断稳定得多。
应变方程藏在损伤力学里
很多人找不到应变演化方程，其实它寄生在损伤变量里。用非线性结构材料模块时，勾选"损伤"选项会自动激活：
// 损伤演化伪代码 ddamage/dt = Gc*(1 - damage) - damage*strain_energy注意损伤阈值要配合岩石的抗拉强度设置。遇到过不收敛的情况？把初始损伤值设为1e-4而不是0，问题秒解。
相场扩散用PDE模块魔改
COMSOL自带的相场模块适合标准裂纹，水力压裂得自己写控制方程：
d(phi)/dt = -Gc*div(grad(phi)) + Gc/epsilon*(1 - phi) + lambda*(strain_energy - critical_energy)epsilon这个参数控制裂缝宽度，建议从0.1倍单元尺寸开始试。求解时打开自适应网格，看到裂缝路径开始蛇皮走位的时候，赶紧把单元尺寸减半。
最后说个骚操作：在求解器配置里把固体力学和相场方程分两组求解，用分离式步进比全耦合快三倍不止。但记得在固体力学步里冻结相场变量，反过来也同理，不然会解出个四不像。
仿真跑到70%突然发散？别慌，把瞬态求解器的BDF方法换成广义alpha，时间步长放大五倍继续跑——这招救过我的毕业设计。
