COMSOL裂缝注浆模拟：非饱和多孔介质Richards方程建模实践

干工程的老伙计们都有种感觉：注浆泵一开，浆液从管口出去之后就彻底“放飞自我”了。你以为它乖乖沿着裂缝走，结果它在地层里绕了个大圈；你以为终于堵住了，一降压，它又从另一条通道溜回来。这背后的物理过程，说白了就是非饱和多孔介质里的流体运动。我用COMSOL搭了一个裂缝注浆的二维模型，把混凝土或岩体简化成孔隙-裂隙双介质，核心用Richards方程，目的只有一个：让模型自己学会“呼吸”——也就是真实反映饱和度随时间上升、回落、再向深处扩散的动态过程。今天就把这套模型的建模思路、参数设置和踩坑记录完整拆开讲一遍，适合正在做裂缝灌浆模拟、地质灾害注浆加固模拟，以及刚接触COMSOL地下水流方向的朋友参考。

1. 裂缝注浆模拟到底难在哪

1.1 浆液跑的路径：非饱和多孔介质里的“可见”与“不可见”

很多人一说注浆，脑子里浮现的是“把浆液压进缝里”，好像裂缝是根水管。但真实的混凝土或者岩体，在注浆前往往处于非饱和状态：孔隙里既有水，又有空气。浆液进入后，它不是在一条光滑管道里跑，而是在毛细力、重力、粘性力三者拉扯下，沿着不同尺寸的孔隙网络慢慢推挤空气、驱替水。这个过程比饱和流复杂得多，因为饱和度一变，渗透率也跟着变，而且是强烈非线性地变。

我一开始习惯性地用达西定律去做这个模型，结果发现完全对不上。饱和达西流动假设介质中的孔隙始终被水充满，渗透率是常数，但注浆伊始，裂缝周围很可能只有残留水膜，大部分孔隙是被空气占据的。浆液要想往前推进，先得克服毛细压力把空气挤走。这时候模型必须能给出一个“负压吸水”的状态，也就是压力水头为负、饱和度小于1的状态，这正是Richards方程擅长描述的场景。你可以把非饱和多孔介质想象成一块半干的毛巾，水倒上去以后，毛巾先局部湿润，然后水分在毛细作用下慢慢向四周扩散，而不是像水管一样瞬间流满。这个过程，就是模型的“呼吸感”来源。

1.2 裂缝的加入：让本来均匀的问题变成“正面刚”

如果说非饱和基质是“难缠”，那加上裂缝就是“更难缠”。裂缝在几何上实在太薄了，动辄0.1到几毫米，如果真按实体建模，网格尺寸会被迫压到毫米级以下，整个计算量直接爆炸。更麻烦的是，裂缝和基质在渗透性上差了不止一个数量级：基质渗透率可能是10^-18甚至更低，裂缝等效渗透率却能到10^-11，这种极端差异对数值求解器是个严峻考验。

做这个模型之前，我做了一次纯基质的均质模型作为对照，结果非常“平滑”——饱和度从注浆点开始，像滴墨水一样均匀扩散开。但现实工程里，注浆压力一到，沿着裂缝方向总是有股浆液“抢跑”的势头，裂缝变成了输浆的主通道，基质反而成了储存库。这种“缝优先、孔隙滞后”的行为，必须在模型里显式体现。我的做法是把裂缝设置为内部边界，让流动方程在裂缝边界上以“切向流+法向交换”的形式表达，既避免了极细网格，又保留了裂缝的高导流性和渗漏特征。

2. 建模思路：物理场选错，后面全是白干

2.1 为什么用 Richards 方程而不是达西定律

COMSOL里地下水流相关的物理场接口有不少：达西定律、Brinkman方程、Richards方程，还有两相达西流动。裂缝注浆这个场景，我最终选了Richards方程作为基质流场的核心方程，原因主要有三个。

第一，Richards方程本身就把饱和-非饱和状态统一在一个框架里。方程的压力变量是压力水头h，当h大于等于0时，介质视为饱和，饱和度等于1；当h小于0时，介质处于非饱和状态，饱和度和相对渗透率都由特征曲线决定。这意味着同一个模型可以覆盖“初始非饱和—浆液到达后局部饱和—停浆后饱和度回落”的完整呼吸周期。

第二，相对渗透率随饱和度变化的规律，在COMSOL里可以直接选van Genuchten模型，配套的参数包括残余水饱和度θr、饱和水饱和度θs、经验参数α和n等。这些参数在岩土工程和混凝土渗透性研究里已经被大量标定过，即使是工程人员也能从经验表里凑出一套合理的初始值，不需要从头做实验。第三，Richards方程对于注浆这种“单相液体宏观推进”的问题足够用了。我不需要显式追踪气液界面，因为空气在开放边界始终可以排出，气相压力近似视为大气压，这样计算负担小很多，收敛性也比完整两相流好。

当然，如果后面要研究浆液的流变性、颗粒沉积或气液界面形态，那就得升级到两相达西流动甚至相场方法。但作为第一步，“看趋势、定参数、估扩散范围”这套需求，Richards方程是最划算的选择。

2.2 裂缝怎么进模型：内部边界与裂缝开度

裂缝建模是本项目最需要动脑的一步。我的方案是把裂缝在几何里画成一条线段，然后在物理场里把它定义成“内部边界”，给这条边界赋予裂缝属性：法向厚度（开度b）、等效渗透率K_f、孔隙率ε_f。计算时，裂缝内的流动沿切向发生，而裂缝与基质之间的质量交换以法向通量形式自动耦合。

等效渗透率我采用经典的平行板立方定律估算：K_f = b² / 12。这个公式把裂缝理想化为两片平行板之间的光滑缝隙，对大多数混凝土冷缝、岩体张开裂隙来说，作为一阶估计是足够的。比如开度0.5mm，算下来等效渗透率大约是2.08×10^-8 m²，比混凝土基质的10^-17量级高了9个数量级，足以体现出“浆液抢跑”的效应。

但要注意，立方定律没有考虑裂缝壁面的粗糙度和填充物，所以它会给出偏高的导流能力。实际模型里我通常会对计算出的K_f乘一个0.1到0.5的折减系数，具体数值根据现场压水试验或者经验判断。这个折减系数可以做成全局参数，后续做参数扫描时非常方便，我后面还会再提。

2.3 材料参数怎么给：VG曲线不是文档模板

很多人一看到van Genuchten那套参数就打怵，其实它的物理含义很直白：α大致对应“进气值”的倒数，决定介质在多大负压下开始失饱和；n则控制饱和度-压力曲线的陡峭程度，n越大，曲线越“刚”，稍微吸点水就迅速饱和。我常用的一组初值（以混凝土类材料为例）大概是：θs = 0.12，θr = 0.01，α = 1e-4 /Pa，n = 1.6，l = 0.5。如果是砂质岩体，θs可能到0.35，α也可能大一个数量级，n到2.0。

我见过不少新手把COMSOL自带的默认参数直接灌进模型，结果算出来饱和度云图要么是一团糊，要么是单薄的同心圆，完全看不出裂缝的影响。问题往往出在相对渗透率曲线上：如果α太大，基质在很轻微的负压下就开始大量导水，浆液会绕过裂缝直接均匀扩散，模拟结果就变成了“均质土球吸水”，这和实际不符。所以每次建模，我都会先算一下裂缝和基质各自在目标饱和度区间内的相对渗透率，确保特征长度相差明显，模型才会表现出“先缝后孔”的推进层次。

3. 实操搭建：从空白模型到第一个“呼吸”结果

3.1 参数全局表：先建模，再谈感情

我在搭建这类模型时的习惯，是先把全局参数表写满，哪怕某些参数后面不用，也把名字占好位。这就像施工前先画好脚手架，后面调参时才知道哪里动了。参数表里至少要有以下几类：

注浆压力p_in不能直接填进Richards方程，需要换算成压力水头：h_in = p_in / (rho_f * g)，0.3MPa算出来大概是17米水头。别看这个转换简单，我见过有人直接把0.3填进“压力水头”里，结果入口只有0.3米水头，浆液推不动，模拟结果全变了。

3.2 几何和网格：该省的地方别使劲，该密的地方别偷懒

几何模型我用一个1×2米的矩形代表混凝土或岩体断面，注浆孔放在左侧中点，一条斜向裂缝从注浆孔延伸到右下侧，裂缝在几何里画成一条线段。这里有个关键操作：裂缝线段必须把矩形内部“劈开”，也就是让线段两边都生成边界，这样内部边界才能在物理场里被识别。

网格这块最考验耐心。裂缝附近的网格必须加密，但又不至于把整个区域都搞成细网格。我的做法是：先给整个域用自由三角形网格，尺寸设为“较细化”，然后在裂缝边界上加一个“边界层网格+附近域细化”的组合。边界层数取3到5层，第一层厚度大约0.002m，这样既能捕捉法向压力梯度，又不会让网格总数超过几万。说句实在话，COMSOL里几何如果简单，网格质量一般不会太差，真正翻车通常发生在裂缝端点：端部网格必须做局部细化，否则饱和度响应在那里会出现奇怪的锯齿状突变。

3.3 边界条件与求解器：注浆压力和注射节奏

边界条件我分三类处理。注浆孔位置设为压力边界条件，压力值随时间可以是阶梯函数或方波，模拟“注浆—保压—停歇”的循环。模型顶部和右侧设为开放边界，允许空气和多余浆液渗出；底部和左侧则设为对称或零通量边界，模拟无限远处的封闭条件。这里的开放边界其实很关键：如果边界封死，注浆开始后内部压力一直累积，饱和度迟早全变1，模型就失去“呼吸”能力了。

求解器方面，Richards方程是非常硬的非线性问题，时间步长建议交给自适应求解器。我通常直接用默认的PARDISO直接求解器，并开启“自动的非常缓慢”的瞬态选项，让COMSOL自己判断什么时候用大步长跳过稳定期，什么时候加密小步长捕捉压力突变。实测下来，这种配置在初值合理的情况下，能稳定算完整个注浆过程。

4. 结果看出花：饱和度、压力和流速怎么组合看

4.1 饱和度云图：看浆液推进到哪块“肺叶”

模型算完的第一步，我总是先看饱和度云图。非饱和介质里的饱和度是0到1之间的数，云图上一片从注浆孔向外扩散的渐变色，刚好对应“浆液湿润前锋”的推进。好的模型表现是：裂缝路径上的饱和度最先冲到接近1，随后沿着裂缝两侧像羽翼一样向外扩散，形成一种“主通道亮，基质渐染”的视觉效果。这个过程很像肺部吸气——气道先打开，肺泡再慢慢充盈，所以我说这个模型会呼吸。

记得我初次算出来时，饱和度云图确实给了我一个惊喜：裂缝方向的前锋已经跑到模型中部，而垂直于裂缝方向只扩散了几厘米。这说明裂缝导流和基质渗吸两种机制被完美分离了。如果你算出来的云图依然像均质模型那样呈现同心圆扩散，那基本可以断定裂缝没有真正参与流动，要么是内部边界没设置成功，要么是裂缝渗透率根本没比基质高出几个量级。

4.2 压力时间曲线：裂缝被浆液卡住的第一信号

在注浆孔和裂缝中部的监测点，我各放了一个“点探针”记录压力水头随时间的变化。非饱和介质里，压力水头从初始负值向正值跃升的那一刻，就是浆液前锋到达的信号。如果你让模型继续算到停浆阶段，还会看到压力回落曲线出现一个平台段，这通常是浆液已经充填了裂缝主通道、正在向基质深层渗吸的表现。

平台段越平，说明裂缝被浆液“卡”得越好，注浆效果越接近“缝内充填满、基质吸浆量小”的理想状态。平台段越陡，说明浆液还在持续漏失，基质的耗浆量很大，需要考虑改用更浓的浆液或调整注浆压力。这套“压力-时间特征”的判断，在真实注浆施工里也成立，所以模型做出来不仅是为了发论文，是真的能反哺工艺参数设计。

4.3 流速矢量：分清串浆、渗漏和有效扩散

饱和度云图只能告诉我们浆液“到没到”，流速矢量图才能告诉我们浆液“跑得多快、往哪跑”。在COMSOL里，Richards方程的后处理变量自带流速矢量，可以直接叠加在饱和度云图上。你会发现，刚开始注浆时，裂缝内的流速矢量又长又集中，指向裂缝末端；几十秒后，裂缝末端附近开始出现向基质偏离的分叉矢量，说明浆液正在向两侧渗吸。

我特别看重矢量图里是否存在“越过裂缝末端继续高速前进”的迹线。如果裂缝末端压力太高，而基质的储集能力有限，流线会绕过失水区向边界“逃逸”，这在施工上很可能对应“串浆”风险。做注浆方案时，如果模拟显示流线直达远场边界，我会建议现场缩短注浆段长、降低注浆压力，或者在远端加设封堵孔，而不是盲目增加注浆量。

5. 常见问题与排查心得

5.1 不收敛：先别怪求解器，看看初值

我遇到过最多的报错就是“求解器在时间步长x.x秒处不收敛”。新手第一反应往往是调小步长或改求解器，但实际上多半是初始条件给得太离谱。Richards方程的初始压力水头必须是一个能在“非饱和”区间自洽的值——如果你直接给0，模型会默认整个介质已经饱和，饱和度初始就是1，那注浆过程自然没有任何“从非饱和到饱和”的动态，数值解还会因为饱和度突变而振荡。

我现在的习惯是先做稳态求解作为瞬态初值，或者用“辅助扫描”功能计算稳态。具体操作是，先不启用注浆压力，只算重力条件下介质自然达到的水头分布和饱和度分布，得到的结果再作为瞬态求解的初始值。这套流程虽然多花几次求解时间，但收敛稳定性提升非常明显。

5.2 饱和度出现负值或超一：数值震荡的几个“熄火开关”

饱和度越界基本可以归结为两类原因：时间步长太大，或空间离散太粗。前者容易解决，把求解器时间步的“最大步长”设到总分析时长的千分之一，同时开启“事件”检测，让压力变化剧烈时自动加密；后者则需要给裂缝前端、注浆孔附近做网格加密，避免压力梯度跨过太粗的单元。

另外还有个参数层面的隐性原因：θr和θs之间的差值太小。有效饱和度在这个差值下被放大，一点点数值误差就会导致饱和度冒出0或1的区间。我一般会把θs和θr的差值保持在0.05以上，同时给饱和度后处理表达式加一个“clamp”式的截断，不改变求解过程但避免出图时吓人。

5.3 裂缝像不存在：边界网格和特征开关

裂缝不起作用通常有两种情况。第一种是几何上裂缝没有真正“切开”域，导致内部边界没有被识别为内部边界，模型里只有基质流动。检查方法是：在“网格”节点下看裂缝边界处是否生成了边界单元，如果没有，回到几何步骤把线段设为“形成内部边界”。第二种是裂缝渗透率设得不够高，尤其是开了粗糙度折减系数后，K_f被折减到了和基质相近的量级，自然就看不到差异。

还有一个小坑：Richards方程接口里，裂缝特征可能需要单独勾选参与求解。如果你在物理场设置里只在“域”里看到Richards方程，而在“边界”里没有对应的裂缝流动特征，那这层物理根本没被激活。我一个朋友就是因为漏了这一步，做了两周的模型，裂缝始终是摆设。

5.4 注浆压力总对不上：单位换算和远场位置

注浆压力换算这个坑，值得单独说一次。工程上注浆压力表通常用MPa，COMSOL里Richards方程默认的压力单位是Pa，但压力水头用的是米。很多人在全局参数表里写0.3，然后让入口边界直接等于这个值，结果0.3Pa连重力都克服不了，浆液当然寸步难行。正确做法是设置p_in = 0.3[MPa]，再在压力水头参数里写p_in/(rho_fg[kg/(m^2s^2)])。在COMSOL的变量表达式里，这一步非常直观。

远场边界的位置也值得斟酌。如果模型边界离注浆孔太近，边界效应会提前“吸收”压力波，让饱和度前锋看起来异常保守；离太远则浪费计算资源。一般我取你认为注浆影响半径的2到3倍即可。对一个1米尺度的注浆影响区，2米模型宽度是够用的。如果拿不准，先算一次，然后看边界附近的压力是否还有明显变化，若边界处压力仍在随时间大幅波动，就该把模型拉大。

6. 从单缝到随机裂隙网络：还能怎么玩

6.1 随机裂缝网络：从“一条缝”到“一张网”

单条裂缝模型的意义在于把物理机制吃透，但真实混凝土或岩体里的裂缝往往是一套网络。把单裂缝扩展成裂隙网络，最省事的办法是用随机几何工具生成若干条裂隙线段，赋予不同的开度和方向，再统一添加内部边界特征。这样做的代价是网格量和非线性程度都会显著上升，但回报也很直接：你能看到浆液优先通过哪一组裂隙、在哪个交叉点顶住、哪里因为裂隙网络连通性差而形成“注浆盲区”。

这个过程中有个很有价值的后处理统计：把不同裂隙开度下的注浆吸收量做参数化扫描，就能找到控制注浆效果的主要裂缝参数。比如开度分布的标准差越大，浆液越容易集中到某几条主缝里，基质反而得不到充分加固。这种规律在单一裂缝模型里是看不出来的，必须靠网络模型才能显现。

6.2 多物理场扩展：热、力和化学灌浆

一旦基础流动模型跑通，往上加物理场就顺理成章了。比如化学灌浆，可以耦合稀物质传递接口，让浆液中的固化剂浓度随流场分布，再通过反应源项描述凝胶时间的空间差异。再比如高地温环境下的注浆，可以耦合固体传热，让粘度随温度变化，浆液在高温区变得更稀、更容易流失。做岩体加固时，还能进一步把固体力学加进来，用“注浆压力引起的孔压变化→有效应力改变→裂缝开度变化”这种多场耦合链，模拟注浆压裂或者劈裂注浆的扩展过程。

我个人的建议是，不要为了耦合而耦合。基础Richards模型本来是稳定可靠的，每加一个物理场，非线性程度都会上一大截，网格和求解器要求也会水涨船高。最好先把纯流体模型的结果与现场实测或经验公式校准完，再逐步加耦合项，不然一旦发散，你根本不知道是哪一层物理场把数值解带崩的。

最后再分享一个我个人的使用习惯：拿到现场压水实验数据后，别急着做全尺寸三维模型，先在二维骨架模型里做三组参数敏感性分析——基质渗透率、裂缝开度、注浆压力。这三组扫描结果能帮你快速判断该工程是“基质控制型”还是“裂缝控制型”。如果是基质控制型，现场采取细颗粒浆液、长时间低压注浆更有效；如果是裂缝控制型，就得抓重点裂缝、调整浆液流变参数和注浆节奏。有了模型做基础，现场调整方案的时候心态会稳很多，因为你大概知道每一步操作在理论层面会产生什么后果。