做注浆仿真的人一定见过这种让人崩溃的画面明明设置的是水泥浆结果算出来的云图跟清水差不多哗啦一下扩散到整个裂隙而现场注浆记录却明明白白告诉你浆液最远只跑了不到一米。问题多半出在材料本构——你把宾汉姆流体当成牛顿流体处理了。这期想聊聊我在 Comsol 里模拟宾汉姆流体浆液扩散的一些经验包括从最基础的模型选型、参数设置到后处理怎么把扩散半径这个工程指标挖出来以及我踩过的几个让人头秃的坑。文章里用到的思路适合做注浆加固、裂隙防渗、水泥基材料研究的同学也适合那些刚开始学非牛顿流体仿真、想找一套稳妥入门路线的朋友。1. 先搞懂宾汉姆流体浆液为什么会半流半固1.1 牙膏和水泥浆的共同点不施加足够压力它就不动理解宾汉姆流体你只需要回想一下早上挤牙膏的场景牙膏立在牙刷上不会自己流下来但你用力一挤它立刻乖乖出来。它的启动需要克服一个门槛这个门槛在流变学里就叫屈服应力yield stress通常记作 τ₀。宾汉姆流体的本构关系可以写成当 τ ≤ τ₀ 时流体不流动剪切应变率为 0 当 τ τ₀ 时剪切应力与剪切速率满足 τ τ₀ μₚ·γ̇。这里 μₚ 是塑性黏度γ̇ 是剪切速率。如果用表观黏度来理解宾汉姆流体的等效黏度并不恒定而是η_eff μₚ τ₀ / γ̇也就是说剪切速率越低表观黏度越大流体越稠剪切速率趋于 0 时表观黏度趋近无穷大所以浆液能像固体一样停在原地。这和水的行为完全不同。水的黏度恒定压力一撤就静止但不具备保持形状的能力水泥浆、黏土浆、泡沫、泥浆这类宾汉姆流体却天然自带一个刹车装置。你做注浆扩散模拟时如果忽略这个特性扩散半径、压力场、填充效果都会严重失真。1.2 屈服应力在注浆工程里到底决定了什么工程上关心的核心问题就一句话浆液能扩散多远、需要多大的注浆压力。而这两个问题都和 τ₀ 强相关。浆液在裂隙或孔隙中流动压力梯度就是驱动力。如果沿径向的压力梯度不足以克服屈服应力浆液就走不动了。计算结果往往是存在一个极限扩散半径 R_max无论注浆时间怎么延长半径都不再增长。这很好理解——压力沿程衰减越往前压力梯度越小当梯度降到临界值以下前沿就冻住了。反过来如果你想扩大扩散半径有两条路加大注浆压力或者降低屈服应力比如加高效减水剂。做仿真分析时这正好对应两个可以自由调整的参数入口压力 P_in 和屈服应力 τ₀。通过参数扫描很容易就能画出一张工程实用曲线。我常用的水泥基注浆浆液参数范围如下具体以你的浆液配合比为准参数典型取值范围说明浆液密度 ρ1400~1800 kg/m³水灰比越低密度越大塑性黏度 μₚ10~100 mPa·s水灰比 0.8~1.0 时常取 20~50 mPa·s屈服应力 τ₀5~30 Pa常用的低屈服注浆材料可能更小注浆压力 P_in0.1~2 MPa视深度和地层而定太大容易劈裂地层做仿真时建议先用这些典型值跑通流程再根据你的实测流变数据替换。2. 建模前先选路线三种模拟方案该怎么挑很多初学者打开 Comsol 的第一反应是找注浆或者浆液扩散的现成案例但 Comsol 更像一个工具箱你需要自己决定用哪一个物理接口来组装这个问题。我做过几次之后总结出三条比较实用的路线。2.1 方案 A层流 非牛顿流体本构适合裂隙、孔洞通道如果你的目标对象是岩体裂隙、结构缝、管道这类几何边界清晰的通道层流Laminar Flow接口非常合适。裂隙中的浆液流速通常很低雷诺数很小流动处于层流甚至蠕动流区间直接用层流接口就很稳。关键是在层流物理场的流体属性节点里本构关系可以选择非牛顿流体下的 Bingham 模型。Comsol 自带了这个模型你只需要输入塑性黏度和屈服应力。这个方案的好处是物理意义直接边界条件直观后处理能直接给出速度场、压力场、剪切速率分布。它在模拟边界清晰的通道内流动时很顺手但如果浆液在推进过程中存在复杂的前沿形态单靠层流接口算出来的还是一个流动区域内的稳态/瞬态流场并不自动追踪浆液边界。2.2 方案 BDarcy 渗流 等效黏度适合多孔介质渗透注浆如果你的对象是砂层、土体这类孔隙性介质Darcy 定律接口更符合注浆机理。多孔介质渗透注浆的经典解析公式比如 Maag 公式就是在达西定律基础上推导的。问题在于 Darcy 接口默认材料是牛顿流体。要让宾汉姆浆液生效需要把动力黏度写成跟剪切速率相关的等效黏度也就是 η_eff μₚ τ₀ / |γ̇|。但这个表达式在剪切速率为 0 时会发散所以工程上经常加一个极小量来避免除零例如eta_eff mu_p tau0 / (1e-3 shear_rate)其中 shear_rate 在 Comsol 里可以使用速度场变量的组合来估算也可以直接用自定义变量近似。这个方案的好处是计算量小、稳定适合大尺度区域缺点是剪切速率的物理定义在达西框架下不算严格结果适合做趋势研究不适合追求精确界面形态。2.3 方案 C移动网格或相场法追踪浆液前缘进阶玩法前两种方案看着是满场都是浆液实际上边界默认处处都有浆液流动。如果你想直观看到浆液从孔口慢慢往外推进、前沿像章鱼触手一样延伸就需要显式追踪相界面。两种主流做法移动网格/变形几何把浆液区域单独划出来区域边界随流动推进。优点是界面锐利缺点是大变形时网格容易畸变计算容易崩。相场或水平集接口在固定网格上用体积分数场描述浆液与空气/水的界面界面预测更稳定但需要额外设置表面张力、两相物性场计算成本明显增加。我的建议很直接如果只是在讨论模型原则上先用方案 A 跑通如果你后续要出漂亮的扩散前锋图再上方案 C。不要一上来就搞移动网格否则大概率被网格畸变折磨到怀疑人生。3. 拿层流方案完整跑通一次几何、参数、求解器下面以方案 A 为例把一条最稳妥的路径走一遍。我用的是 Comsol 6.2 版本如果你用的是 6.4 或 5.6界面可能略有差异但物理场节点逻辑基本一致。3.1 几何创建与边界条件设计第一步是在模型向导里选择二维空间维度添加层流spf物理场。几何可以很朴素一段矩形区域表示裂隙的剖面。比如长 5 m、开度 0.02 m 的裂隙注浆孔在左端孔直径按 0.1 m 处理即可。如果真的要精细模拟裂隙面粗糙度可以加一些微小起伏但初期建议先做光滑平直裂隙这是为了排除几何因素干扰先验证材料模型本身。边界条件设计边界位置条件类型取值建议注浆孔口压力入口先给 0.5 MPa换算成相对压力裂隙远端压力出口0 Pa相对压力裂隙上下壁面壁面无滑移对称轴如有对称可减少一半计算量请特别留意一点注浆压力是瞬态驱动还是稳态维持如果你只关心浆液最终扩散范围稳态求解就够了如果你关心注入一定时间后的浆液分布就必须做瞬态求解压力入口和出口保持不变把时间拉长观察。3.2 在 COMSOL 里把宾汉姆流体写进去这一步是整个模型的核心。在物理场研究树中展开层流找到流体属性节点点开后本构关系下拉框里选Bingham plastic不同版本叫法可能略有不同有的叫 Yield stress fluid 或 Bingham-Papanastasiou但思路一回事。有三个参数需要填密度 ρkg/m³塑性黏度 μₚPa·s屈服应力 τ₀Pa如果你用的模块版本里找不到 Bingham 内置模型可以改成用户定义黏度然后在黏度表达式里写mu_p tau0 / (1e-3 spf.sr)这里的spf.sr就是层流接口计算出的剪切速率1e-3是为了防止除零。使用内置模型和自定义表达式的差别主要在于内置模型可能有更精细的非线性平滑处理但自定义表达式的可调控空间更大。别忽视量纲问题。工程数据里塑性黏度经常用 mPa·sComsol 方程体系默认单位是 Pa·s转换系数是 1 mPa·s 0.001 Pa·s。填入表观黏度前先做一次单位换算否则结果可能差了三个数量级而毫不自知。3.3 网格、求解器与时间步设置网格建议分两步走第一步用默认物理场控制网格单元大小选细化或较细化快速验证模型能否收敛。 第二步再把裂隙壁面附近加密因为宾汉姆流体在壁面附近剪切速率变化剧烈那里是黏度变化最敏感的区域。求解器方面稳态问题直接用默认的稳态研究即可。瞬态问题要注意时间步的设置经验如果注浆压力是恒定压力流体启动很快初期时间步要小比如从 0.01 s 起步后期流动趋稳可以让求解器自动放长时间步直接在瞬态求解器的时间步进里选择 BDF 方法容差保持默认但把初始步长手动设为 1E-3 s 到 1E-2 s 之间比较保险。我遇到过很多次一开始就红屏的情况十有八九是初始步长太大。压力越高、屈服应力越接近临界值初始步长就要越小。3.4 先检查流场合理性再谈扩散计算完成后别急着看扩散半径先在通用后处理里切三个图验证物理合理性速度场云图浆液应该从入口向出口延伸呈现一个楔形速度剖面壁面附近速度低、中心速度高。压力场云图高压区在注浆孔口附近沿流动方向平滑下降。剪切速率分布入口和壁面附近剪切速率高远离入口处剪切速率低。如果速度场出现大面积发散或震荡条纹返回到网格和时间步那一节调参数不要继续往下做后处理否则结论必然是垃圾进垃圾出。4. 从结果里挖出扩散规律剪切速率、扩散半径和压力当模型能稳定出结果以后仿真才真正开始有意义。下面讲几个我常用的后处理思路它们决定了你能从模型里带走哪些工程结论。4.1 用剪切速率画活性和死区宾汉姆流体最直观的特征就是有的地方在流有的地方已经冻住。在后续处理里定义一个变量活性指数最简单的方式是直接用剪切速率阈值Active if(spf.sr 0.001, 1, 0)阈值取多少要看剪切速率的量级。可以把剪切速率云图先画出来看分布范围再选一个比最大剪切速率小三四个数量级的值作为阈值。然后用表面绘图里的表达式过滤功能把 Active0 的区域显示成灰色这就是浆液死区。如果把死区区域叠加到裂隙几何上看你会发现一个很有趣的现象死区往往出现在两条条带里——一条是壁面附近的极薄层因为黏附壁面所以几乎不动另一条是扩散前锋远端那里压力梯度已经衰减到不足以克服屈服应力。这种做法比直接看速度场更能表现宾汉姆流体的特殊性。4.2 扩散半径-时间曲线怎么提取工程上最关心的扩散半径 R(t)本质上是一个随时间变化的界面位置。在方案 A 的框架下最方便的做法是定义一组探针或直接使用派生值 体积积分来追踪。更简单的思路在裂隙中心线上画一条截线某一时刻沿截线提取速度大小把速度下降到入口速度 1% 的位置定义为前锋位置。从多个时间点提取后就能画出一条 R-t 曲线。实测下来的曲线通常会呈现两个阶段初期快速扩展入口附近压力梯度大浆液像挤牙膏一样快速前冲后期渐进停滞压力梯度衰减扩散半径增长越来越慢最后趋近某个极限值 R_max。这个 R_max 就是注浆方案设计的核心参数。你可以把数值模拟得到的 R_max 与经典解析公式对比如果相差在 10% 以内说明模型参数标定基本可靠如果差得远先去检查屈服应力和压力边界条件的换算。4.3 参数化扫描屈服应力与注浆压力如何匹配Comsol 的参数化扫描功能非常适合做这类分析。把 τ₀ 设为参数取一组值比如 5、10、15、20、25 Pa扫描后叠加画出各条 R-t 曲线马上能看到屈服应力对扩散半径的压制效果。同样把入口压力 P_in 设为参数从 0.2 MPa 扫到 1 MPa可以得到注浆压力—极限扩散半径曲线。这两组扫描做完你基本就能回答在这个地层参数下想达到 3 米扩散半径需要多大压力、多少注浆量这类问题。有一组很有意思的临界规律当入口压力产生的压力梯度刚超过屈服应力阈值时浆液可能纹丝不动。不要误以为这是数值错误这恰恰是宾汉姆流体的物理特性。估算方法很简单对开度为 b 的裂隙临界压力梯度大致是 2τ₀/b。所以先手算一个临界压力再去设置参数扫描的范围才能避开无意义的全零结果。5. 避坑实录收敛性、网格畸变、参数魔改5.1 找不到一致的初始值或直接 NaN问题出在哪这是评论区出现频率最高的问题。遇到这种报错先按下面顺序排查确认所有材料参数都已填写且单位正确。缺失密度或黏度是最常见的低级错误。检查初始值设置。把层流的初始压力设为入口压力的 50%初始速度设为 0往往能大幅提高收敛率。检查自定义黏度表达式里是否有除零风险。1e-3这个保护项看起来小实际上在低速区作用很大。如果用了瞬态求解把初始步长再缩小两个数量级。收敛性跟初始步长密切相关很多 NaN 都是初期步长过大导致的。如果上述都试过还是报错就退回到层流 牛顿流体先算一遍确认模型框架没问题再把本构切换成 Bingham。这种从简到繁的排查思路能帮你快速定位是哪一步引入的问题。5.2 移动网格一时爽网格畸变火葬场老实说用移动网格追踪浆液前锋时我最开始被网格自交和负雅可比逼疯过。后来总结出的经验是移动网格适合中小变形和边界位移可控的场景比如盾尾注浆中浆液沿环隙推进不适合大范围自由扩散因为浆液前沿越高网格就被拉得越厉害扭曲程度往往超过网格重划分的合理范围。如果你真的需要做前沿追踪建议优先考虑网格平滑方法选超弹性它对大变形的容忍度比 Laplacian 平滑高很多把移动边界设置成指定网格位移而不是指定网格速度控制更直接定期在求解设置里加网格重构事件让 Comsol 自动重新划分网格如果模型域很大、变形很剧烈就改用相场或水平集接口不要在移动网格上死磕。做仿真要记得目标是什么。如果你要的是扩散前锋形态和半径相场法虽然计算量更大但稳定性和物理一致性都远好于硬凹移动网格。5.3 参数扫描时的小技巧先粗后细先稳态后瞬态参数扫描最忌讳一次撒大网。我现在的习惯是先用粗网格、稳态求解跑一遍全部参数组合筛出有物理意义的大致范围锁定感兴趣的区域后加密网格并改瞬态逐个参数精细计算在同一研究中复用结果在参数扫描里设置扫描表格时小心未求解的组合不要依赖先前解。另外如果想批量提取 R-t 曲线做后续绘图可以用 Comsol 的 App 开发器录制一段批处理也可以在 Linux 上用命令行启动 comsol batch 跑批量任务。2024 年之后的版本对 Java 客户端支持得不错配合脚本做参数寻优已经是很常规的玩法。5.4 说点个人体会我自己最初做这个仿真时总觉得把模型调得越高级越有成就感后来发现真正有用的是把物理规律搞对。一个用层流接口快速跑通的二维单裂隙模型比一个半途崩掉的带移动网格三维模型更有说服力。宾汉姆流体扩散这个问题里屈服应力是关键中的关键任何让计算结果偏离屈服应力影响的细节都是次要细节。浆液在裂隙里走走停停的过程看起来也像极了一个性格固执的人在缓慢前进压力给够了它往前挪一步压力一放松它就立刻原地驻留你以为它彻底停了其实它只是因为没攒够压力而已。这个想法也是我把这篇标题叫奇妙之旅的原因。最后再分享一个小经验如果你手头有现场注浆数据一定拿它来校核仿真。把给料流量、注浆时间和最终扩散半径三个数据对上你才会真正信任这个模型。后续想继续扩展的话还可以把水泥浆的时变性黏度随时间增长加进宾汉姆模型或者在裂隙网络里做随机裂隙分布下的扩散规律研究那是另一个有意思的话题了。