COMSOL B-B本构方程在三维粗糙裂隙岩体建模中的应用及热流固三场耦合研究

comsol B-B本构方程对三维粗糙裂隙岩体建模，并耦合热流固三场 研究其温度场，应力场，渗流场的变化

裂隙岩体里的热流固耦合总是让人头大，尤其是当三维粗糙裂隙遇上非线性本构模型。最近在COMSOL折腾Barcelona Basic本构模型时，发现这玩意儿对岩体塑性变形的描述确实带劲——特别是当温度场掺和进来，整个模型就像被扔进高压锅的土豆，既酥软又充满不确定性。

先说说怎么给裂隙造个"糙脸"。在几何模块里用随机高度函数生成裂隙面，核心代码就三行：

rng(2023); % 固定随机种子 roughness = 0.1*randn(N_points); % 标准差控制粗糙度 surface = mphgeom(model,'rough_surface','script',['createSurface(',roughness,')']);

这相当于给裂隙面喂了把跳跳糖，生成类似天然岩体表面的凹凸结构。网格划分记得用四面体二阶单元，不然粗糙边缘会像狗啃的似的。

B-B本构的关键在弹塑性转换。材料属性设置里有个魔鬼参数——硬化模量H。在固体力学接口的塑性子节点下：

model.material("mat1").propertyGroup("Plasticity").set("H", "200e6+10e3*T");

这里的温度T可不是摆设，热膨胀会让屈服面像气球一样膨胀。实测当温度超过80℃时，等效塑性应变会突然飙车，这时候需要把求解器的自动阻尼因子调到0.5以上才能收敛。

comsol B-B本构方程对三维粗糙裂隙岩体建模，并耦合热流固三场 研究其温度场，应力场，渗流场的变化

渗流场和温度场的私奔关系更精彩。达西定律模块里有个隐藏关卡——渗透率随裂隙开度的立方变化：

k = (b^3)/(12*mu)*exp(-alpha*(sigma_eff/1e6));

sigma_eff是有效应力，这指数项能把应力变化放大十倍。有次忘了温度对流体粘度mu的影响，结果在120℃工况下渗流速度直接飞升三倍，导师盯着结果问是不是核废水泄漏了。

耦合迭代最刺激的是热-力边界条件。在裂隙接触面设置双向耦合：

contact_pair.thermalContact().set("ThermalConductivity", "k_rock*(1-plastic_strain)"); contact_pair.mechanicalContact().set("Penetration", "H*exp(-T/100)");

这相当于给热传导加了塑性应变的反向调节器，当裂隙压密时传热效率反而下降。有一次算例出现温度场倒挂现象，排查半天发现是这公式里的指数项符号搞反了。

最后说个血泪教训：多场耦合的求解顺序千万别头铁用全耦合。先算稳态温度场，再逐步加载机械应力，最后耦合渗流场，这样收敛成功率能从抓阄水平提升到抛硬币水平。当看到三个场的云图像彩虹蛋糕一样层层展开时，那种快感堪比在数值模拟的荒漠里找到绿洲。