干工程的老伙计们都有种感觉注浆泵一开浆液从管口出去之后就彻底“放飞自我”了。你以为它乖乖沿着裂缝走结果它在地层里绕了个大圈你以为终于堵住了一降压它又从另一条通道溜回来。这背后的物理过程说白了就是非饱和多孔介质里的流体运动。我用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 /Pan 1.6l 0.5。如果是砂质岩体θs可能到0.35α也可能大一个数量级n到2.0。我见过不少新手把COMSOL自带的默认参数直接灌进模型结果算出来饱和度云图要么是一团糊要么是单薄的同心圆完全看不出裂缝的影响。问题往往出在相对渗透率曲线上如果α太大基质在很轻微的负压下就开始大量导水浆液会绕过裂缝直接均匀扩散模拟结果就变成了“均质土球吸水”这和实际不符。所以每次建模我都会先算一下裂缝和基质各自在目标饱和度区间内的相对渗透率确保特征长度相差明显模型才会表现出“先缝后孔”的推进层次。3. 实操搭建从空白模型到第一个“呼吸”结果3.1 参数全局表先建模再谈感情我在搭建这类模型时的习惯是先把全局参数表写满哪怕某些参数后面不用也把名字占好位。这就像施工前先画好脚手架后面调参时才知道哪里动了。参数表里至少要有以下几类参数符号数值示例单位说明模型高度H1.0m二维断面高度模型宽度L2.0m二维断面宽度裂缝开度b0.0005m裂缝法向厚度裂缝渗透率折减f_f0.31粗糙度折减系数基质渗透率K_m1e-17m²混凝土或岩体基质渗透率基质孔隙率eps_p0.11有效孔隙率浆液密度rho_f1800kg/m³水泥浆密度浆液动力粘度mu_f0.05Pa·s等效牛顿粘度注浆压力p_in0.3MPa注浆孔处表压初始压力水头h0-3m对应非饱和状态分析时长t_end3600s模拟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 注浆压力总对不上单位换算和远场位置注浆压力换算这个坑值得单独说一次。工程上注浆压力表通常用MPaCOMSOL里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模型本来是稳定可靠的每加一个物理场非线性程度都会上一大截网格和求解器要求也会水涨船高。最好先把纯流体模型的结果与现场实测或经验公式校准完再逐步加耦合项不然一旦发散你根本不知道是哪一层物理场把数值解带崩的。最后再分享一个我个人的使用习惯拿到现场压水实验数据后别急着做全尺寸三维模型先在二维骨架模型里做三组参数敏感性分析——基质渗透率、裂缝开度、注浆压力。这三组扫描结果能帮你快速判断该工程是“基质控制型”还是“裂缝控制型”。如果是基质控制型现场采取细颗粒浆液、长时间低压注浆更有效如果是裂缝控制型就得抓重点裂缝、调整浆液流变参数和注浆节奏。有了模型做基础现场调整方案的时候心态会稳很多因为你大概知道每一步操作在理论层面会产生什么后果。