COMSOL双重介质注浆模拟：裂隙岩体建模全流程解析

做注浆模拟的朋友应该都有过这个疑惑：现场返回的注浆量、扩散范围，和数值模拟结果怎么总是对不上？

我前几年做坝基帷幕灌浆项目时也趟过这个坑。当时用的还是“等效渗透率”思路，把整个岩体当成均质连续介质算，结果预测的浆液扩散半径比现场实测小了一大截。反复排查后发现，问题不在参数，而在建模假设——真实岩体里浆液的流动路径，同时受两套系统控制：一套是张开度大的裂隙网络，另一套是密布孔隙的岩块基质。这两套系统的渗透特性可以差三四个数量级，硬要折算成一个等效渗透率，等于把主通道和盲道混在一起，物理图像直接糊了。

后来我改了思路，用COMSOL搭建双重介质注浆模型，把裂隙和基质的流动分开建模，再通过质量和压力耦合起来。调整之后的压力场、扩散趋势和现场反馈基本能对上。这篇文章就把这套建模方法的完整思路拆开讲，从双重介质理论、参数计算，到COMSOL里的具体操作、网格处理、求解器设置，再到我踩过的几个典型坑，一次性梳理清楚。

内容比较适合正在做裂隙岩体注浆、压水试验、地下水渗流模拟的同行，也适合刚接触多孔介质和裂隙流动的COMSOL新手。

1. 双重介质模型：为什么裂隙和基质必须分开算

1.1 裂隙体系和基质孔隙体系，根本不是一个世界

先理解“双重介质”在说什么。岩体不是一坨均匀的东西，它同时存在两套储渗系统：

• 裂隙系统：节理、层理、断层破碎带这类张开结构，渗透率通常很高，但孔隙度很低，储存流体能力弱。它的作用主要是“通道”，浆液沿着裂隙跑得飞快。
• 基质系统：完整岩块内部的微孔隙和微裂缝，渗透率极低，但孔隙度大，储存流体的能力强。它的作用是“储层”，或者说“迟滞区”，浆液进入基质后速度会明显减慢。

用一句话概括：裂隙决定流动在哪里发生，基质决定流动有多快结束。注浆过程中，压力首先驱动浆液沿裂隙网络快速推进，同时裂隙壁面上的压力会把一部分浆液压入基质孔隙，形成渗滤和扩散。如果只把两套系统合成一个平均渗透率，裂隙的“快速通道”效应和基质的“吸水缓冲”效应都会被抹平，结果是扩散路径不对、注浆压力时程不真、浆液损失量估算误差大到离谱。

1.2 Warren-Root模型和系统间的交换机制

双重介质概念最早由Barenblatt等人在1960年提出，后来Warren和Root给出了更便于工程应用的双重连续介质模型。它的核心思想很简单：在同一个空间点上，存在两个相互叠加的“连续体”——裂隙连续体和基质连续体。两个连续体各自的压力、速度、饱和度各算各的，两者之间通过一个**窜流项（exchange term）**进行质量交换。

窜流项的本质是描述基质和裂隙之间的“压力差驱动流动”。公式可以写成：

q_exchange = α_s · (p_f - p_m)

其中 α_s 是形状因子（shape factor），和裂隙间距、裂隙几何形状有关；p_f 和 p_m 分别是裂隙和基质中的压力。这个式子告诉我们，基质向裂隙补给还是裂隙向基质泄流，完全取决于二者之间谁的压力更高。

在COMSOL里，当我们使用“Darcy定律”接口并勾选裂隙流动选项时，这个窜流交换其实是被自动处理的——它不是靠手动加一个源项，而是通过基质域和裂隙边界之间的连续质量守恒耦合起来的。但理解这个交换机制仍然很重要，因为它直接影响你对结果的判断。比如你发现基质区域压力迟迟不涨，多半是窜流交换太弱，原因可能是裂隙壁面积太小，或者基质渗透率被设得极低。

1.3 COMSOL实现双重介质的路线选择

用COMSOL建双重介质模型，主要有两条路线：

路线A：单物理场 + 裂隙流动边界条件（推荐）

在“Darcy定律”物理场中，基质区域按常规多孔介质域设置，裂隙则用“内部边界”上的“裂隙流动”条件来表示。在2D模型里裂隙是一条线，在3D模型里裂隙是一个内部曲面。COMSOL会在裂隙边界上自动附加沿切向的达西流动方程，同时考虑边界两侧的法向压力连续性。

这条路线最大的优势是：不用手动处理基质-裂隙的窜流交换，因为两个区域共享同一套计算网格和同一个压力场，质量守恒天然满足。适合大多数工程尺度的注浆和渗流模拟。

路线B：双物理场耦合（高阶灵活）

分别建立两个Darcy物理场，一个算基质，一个算裂隙，中间通过边界条件手动耦合。这条路灵活度更高，可以模拟非平衡窜流、不同温度场，但实现难度大，收敛性也更娇气。除非你在研究多机制耦合，否则我不建议新手一开始就走这条路。

实际建模时，我先用路线A跑通整个框架，确认压力分布和扩散趋势没问题后，再根据研究需求往里面叠加化学反应、温度场等额外物理场。

2. 建模前的方案规划：从地质条件到COMSOL的映射

2.1 几何简化与裂隙网络取舍

拿到一个实际工程，第一步不是急着打开COMSOL，而是先做地质抽象。现实中裂隙可能密密麻麻，全部建模既不现实也没必要。我一般的取舍原则是：

• 只保留张开度大于0.1mm、延展长度超过模型尺寸1/10的主控裂隙。
• 网状细支裂隙通过降低基质渗透率来等效吸收。
• 裂隙倾角、间距尽量统计平均，避免过度细节导致网格质量崩坏。
• 在二维模型里，先把裂隙表示成折线/曲线，确认流动主方向清晰后，再考虑扩展到三维面裂隙。

举个例子，我那次坝基灌浆模拟，野外节理统计显示有两组优势裂隙，一组倾角55°，一组倾角80°。我就只保留这两组，把它们简化成交叉的网络，产物非常清晰，计算量也小。建模的目标永远是“够用”，不是“真实世界的完美映射”。

2.2 物理场选型：Darcy定律为何够用

注浆流渗流速度通常很慢，雷诺数远小于1，惯性力可以忽略，这种情况下达西定律是适用的。COMSOL的“Darcy定律”接口就是为这个场景设计的，它计算基于多孔介质中流体的压力梯度与渗流速度的关系：

u = -(K / μ) · (∇p + ρg∇D)

K是渗透率，μ是动力粘度，p是压力，D是高程。

如果你的场景里流速较快（比如裂隙特别宽、注浆压力特别高），或者需要模拟裂隙内的非线性流（福希海默流动），那就得考虑Brinkman方程或Navier-Stokes与Darcy的耦合。但从工程注浆角度，绝大多数情况下Darcy定律就足够了。用不到方程越复杂，参数越难标定，这是我坚持先用简单模型的原因。

2.3 单位制与量纲统一

COMSOL虽然支持自定义单位，但我在所有模型里一律使用SI基础单位：

一个特别容易出错的地方是渗透率单位。地质工程里经常碰到“达西（D）”和“毫达西（mD）”，但COMSOL默认界面里如果选择SI单位制，渗透率就是m²。换算关系是：1 D ≈ 1e-12 m²，1 mD ≈ 1e-15 m²。我在最初建模时吃过这个亏，把10 mD的渗透率直接填成10，结果压力场完全不正常，后来检查单位才查出来。

另一个易错的是压力单位。如果入区压力用MPa，出口压力用kPa，边界条件里必须统一转换，否则计算出来的速度场会差出三个数量级。最笨也最安全的办法：所有物理量都先换算到SI，不偷懒。

3. 参数设置与计算逻辑

3.1 基质渗透率与孔隙率的取值逻辑

基质渗透率通常来自现场压水试验或室内岩芯渗透试验。如果手头数据有限，可以参考经验值：完整花岗岩的渗透率大约在1e-18 ~ 1e-16 m²，砂岩在1e-15 ~ 1e-13 m²，泥岩则更低。

孔隙率与渗透率之间不是简单的线性关系。在COMSOL里如果愿意，可以直接用Kozeny-Carman公式估算：

K = ε³ / (C · S² · (1-ε)²)

其中ε是孔隙率，S是比表面积，C是经验常数。但说实话，对于工程模拟，我更推荐直接输入现场实测渗透率，因为Kozeny-Carman的估算偏差经常超过一个数量级。

基质孔隙率对压力扩散的影响相对温和，但对浆液渗滤量影响很大。在注浆模拟中，基质孔隙率除了质量守恒外，还决定了浆液进入基质后的“储留空间”，取值需要尽量贴近岩性实测。

3.2 裂隙渗透率与立方定律

裂隙的渗透率计算通常用平行板模型的立方定律：

k_f = b_f² / 12

b_f是水力开度（m）。注意这里算出来的是渗透率，单位m²，不是导水系数。比如一条张开度0.5mm的裂隙：

k_f = (5e-4)² / 12 ≈ 2.08e-8 m²

这个值和基质渗透率1e-16 m²相比，高出8个数量级。所以浆液在裂隙里流动，本质上是“高速公路”。

在COMSOL的“裂隙流动”边界条件里，需要分别填裂隙孔隙率ε_f和裂隙渗透率k_f。裂隙孔隙率通常取0.5~1.0之间，它的意义是裂隙内能容纳流体的体积比例。如果裂隙内部还有充填物，就取低一点。

3.3 浆液本构模型与时变粘度

注浆浆液不是水，也不是恒粘度流体。水泥基浆液在泵送和渗透过程中，粘度会随着水化反应逐步升高，最终失去流动性。化学浆液的行为更复杂，可能是宾汉姆流体，有屈服应力，也可能呈幂律特征。

对于第一版双重介质注浆模型，我建议用“等效时变粘度”来近似浆液的凝结过程。先把浆液当成牛顿流体，但让动力粘度随时间从初始值增加到终凝值：

η(t) = η₀ · (1 + (η_final/η₀ - 1) · S(t))

其中S(t)是一个从0平滑上升到1的函数，可以用COMSOL自带的平滑阶跃函数flc1hs或flc2hs实现。这里有一个容易忽略的重点：不要用普通Heaviside阶跃，必须用平滑过渡。否则粘度从0时刻直接跳变到终值，数值求解器会因时间导数剧烈振荡而崩掉。过渡时间dt_gel我一般取凝胶时间的5%~10%，比如凝胶时间600s，过渡时间就设60s。

浆液的密度也会影响重力项，但在高压注浆场景下重力影响相对小，密度取一个定值1200~1500 kg/m³就够了。要精确的话可以定义ρ = ρ_water + (ρ_slurry - ρ_water) · S(t)。

3.4 边界条件与初始条件

边界条件要尽量还原现场注浆孔的工艺参数：

• 注浆孔：最常见的是定压力边界，设p_in = 2e6 Pa（2MPa），对应现场泵压；或者定流量边界，对应泵送的体积流量。
• 模型远端/排水边界：通常设为定压力边界，p_out = 0，模拟影响半径之外的压力平衡。
• 对称边界：如果模型是对称的，可以切一半，对称面设置零通量。
• 初始条件：初始压力场一般设为0（相对大气）或静水压力分布。

我自己的习惯是，先用稳态求解器跑一个固定粘度的“启动工况”，把压力场稳下来，再把这个稳态解作为瞬态求解的初始值。这样瞬态求解的收敛性会好很多，能省掉不少调试时间。

4. COMSOL实操流程：从新建模型到出图

4.1 几何构建与裂隙内边界处理

打开COMSOL后，先建一个2D模型。几何上画一个矩形域代表岩体，比如50m × 50m；然后画裂隙线。注意裂隙线必须完全落在矩形域内，端点不能悬空或超出域边界太多，否则后续边界条件会出错。

这里有个小技巧：裂隙交叉处尽量不要出现“尖角交点”，最好做一点倒角或直接把交叉点当成普通交点处理。裂隙线越多，交叉角度越小，网格质量越差，求解越容易出问题。我一般在CAD里把裂隙线先整理成“折线公差不超过1cm”的干净线，再导入。

4.2 全局参数与变量定义

在“全局定义”里新建参数表，把前面计算的参数都写进去。一个典型的参数表如下：

这里建议把k_f直接写成表达式，而不是填一串数字，方便后面改裂隙开度时自动更新。

4.3 物理场接口与裂隙流动设置

在添加物理场时，选择“Darcy定律”。然后：

1. 在“多孔介质域”设置里指定基质域的孔隙率eps_rock和渗透率K_rock。
2. 右键物理场，添加“裂隙流动”边界特征，并选中所有裂隙线。
3. 在裂隙流动设置里填入裂隙孔隙率eps_f和裂隙渗透率k_f。

这一步是双重介质建模的核心，千万别漏。很多新手在几何里画了裂隙线，但物理场设置里忘了勾选“裂隙流动”，结果裂隙完全不起作用，算出来的结果跟均质介质一模一样。

流体属性里，密度填rho_f，动力粘度填eta0 * (1 + (eta_final/eta0 - 1) * flc1hs(t - t_gel, dt_gel))。这样就实现了粘度随时间上升的时变效应。如果后面要追踪浆液前沿，这里再叠加一个浓度相关函数，详见第5.4节。

边界条件设置：在注浆点/注浆孔边界加压p_in，在出口加压p_out。如果是点源，可以添加“点”特征，选择“源/汇”，输入质量流量。

4.4 网格划分策略

网格是COMSOL模拟最容易拖垮收敛的一环，尤其带裂隙的时候。裂隙线虽然几何维度低，但在物理上它却是一条关键通道，网格必须足够细才能解析裂隙内的速度变化。

我的网格策略是：

• 全局“自由三角形网格”最大单元尺寸设为2m。
• 所有裂隙线设置“边网格”最大单元尺寸0.1m。
• 注浆孔附近再局部加密，最大单元尺寸0.05m。

如果裂隙是曲曲折折的，第二步的细节就非常关键。裂隙周围没有细网格，哪怕物理场设置正确，误差也会把浆液扩散路径完全抹掉。还有一种经验性验算：裂隙网格最大尺寸最好不超过裂隙开度的100倍。对于0.5mm开度的裂隙，0.1m的网格其实已经很粗了，但在2D模型里裂隙是“线”，不是“带”，所以网格实际控制的是沿裂隙方向的分辨率，0.1m通常够用。

4.5 求解器与瞬态时间控制

稳态工况用默认求解器就行。瞬态工况要注意时间步的设置：

• 时间范围：0 ~ 1800s（对应模拟30分钟注浆）。
• 时间步长：如果粘度是平滑过渡的，可以直接用默认的自适应步长；但保险起见，我会在“瞬态求解器”设置里把初始步长设为1s，最大步长不超过60s。
• 容差：默认“物理场控制容差”通常可用，但如果遇到收敛困难，我会把相对容差调严到1e-4。

求解器选择上，PARDISO对多物理场问题比较皮实。多核机器上记得开启PARDISO多线程，速度提升很可观。

5. 常见问题与排查技巧实录

5.1 裂隙流动没生效，怎么识别和修复

症状：算出来的压力云图是均匀圆环扩散，完全看不出裂隙的“指状”优势通道。

原因八成是裂隙边界没有激活“裂隙流动”条件。检查方法是在结果里画速度场流线，如果流线在裂隙线附近依然横穿，而不是沿着裂隙展开，那就是裂隙没起作用。

修复办法：回到“Darcy定律”物理场，右键添加“裂隙流动”特征，选中几何里的裂隙边，填好渗透率和孔隙率。另外，如果裂隙线画在了域外面，或者没有完全嵌入域内部，COMSOL不会把它识别为内部边界，这时候需要回到几何处理。

5.2 压力、流量对不上，先怀疑单位换算

症状：注浆压力2MPa，计算出来的流速只有1e-8 m/s，小得离谱。

解决思路：检查渗透率单位。若把现场报告的“10 mD”直接填成10，而模型单位是m²，实际就只有10 m²，那是完全错误的，应该填 10e-15 m²。这种量纲错误经常导致压力场分布看起来正常，但速度场小几个数量级。常规自查：岩体渗透率在1e-18~1e-13 m²之间，如果你填的渗透率大于1e-10 m²，先停一下，确认单位。

5.3 瞬态不收敛，最大的坑是粘度突变

症状：瞬态求到某一步后报错“求解器在时间步长上发散”，或干脆一开始就退出。

最常见的元凶是粘度从初始到终凝的阶跃变化。我前面反复强调用flc1hs平滑过渡，就是为了避免这个坑。当你发现初始步长无论怎么调小都不救不回来时，先看粘度函数是否连续。另外，初始压力场设为“零场”也会让第一个时间步内的压力建立过程非常暴力，推荐先用稳态求解器求一个初始压力分布，再把它作为瞬态初始值。

5.4 浆液前沿怎么追踪

等效时变粘度只能描述“整个区域粘度都随时间增加”，无法捕捉“浆液走到哪里了”。要真正追踪浆液背影，需要再加一个物理场：多孔介质中的稀物质传递（Transport of Diluted Species in Porous Media）。

设置方法：定义浓度变量c，初始值c=0，注浆入口处固定浓度c=c0。粘度改为依赖浓度：

η = η0 · (1 + (η_final/η0 - 1) · c/c0)

这样浆液推进到的地方粘度高，未到达的地方粘度低，“前沿”就自然产生了。扩散系数在浆液渗透过程中通常很小，甚至可以不考虑扩散，只考虑对流。

在使用这个方案时，需要注意数值振荡。对流占优时，浓度场容易发生“数值过冲”，浓度超过1或者变负。解决方法是加密网格，或者在“稀物质传递”界面里选择“对流稳定”选项（COMSOL默认会加，但有时强度不够，需要手动调）。

5.5 与Fluent多孔介质设置的对比速查

有朋友问过：“这个模型用Fluent能不能做？”当然能，但两者的工作流差别很大。这里列个对比表：

一句话总结：如果你的核心关注点是“浆液在多孔介质和裂隙中的渗流过程 + 多物理场耦合”，COMSOL更顺手；如果还要考虑泵送管内的湍流、弯管阻力、高速射流等细节，Fluent那边更强势。具体选哪个，取决于注浆系统的哪个环节是你的主战场。

6. 结果解读与工程扩展思路

6.1 从压力场和流线看注浆主通道

模型跑通后，第一个要看的图是压力分布云图。一个典型的双重介质注浆模型，压力云图会呈现明显的“各向异性”：沿裂隙方向压力衰减慢，等势线被拉长；垂直裂隙方向压力衰减快，等势线密集。如果你看到的压力云图是完美同心圆，说明裂隙没有真正参与流动，需要回头检查裂隙流动设置。

接着看速度流线。流线会明显集中在裂隙网络内部，基质区域的箭头短而稀疏。这说明浆液的主体流量走的是裂隙通道，基质主要承担压力和渗滤扩散。

6.2 扩散半径与注浆量评估

在瞬态结果中，可以用浓度场/粘度场识别浆液前沿。工程上通常参考“注浆量达到设计值”或“压力突增”来判断结束时机。模型中，你可以提取不同时刻的浆液扩散面积，换算成等效扩散半径，再和现场考察结果比对。

如果你做的是定流量注浆，还能直接从出口压力随时间的上升趋势里看到“粘度上升导致压力抬升”的曲线。这个曲线对工艺控制挺有参考价值——压力出现拐点的时刻，对应浆液可泵送性的临界时间，现场施工可以在接近这个时刻前控制好注浆速率。

6.3 模型升级方向

第一版双重介质注浆模型跑通后，后续可以扩展的方向很多：

• 渗透率时变封堵：把基质渗透率定义成浓度/时间的函数，模拟浆液渗滤导致的孔隙堵塞。
• 温度场耦合：如果浆液水化放热明显，可以加入固体传热，研究温度对粘度演化和岩体热应力的影响。
• 两相流/多相流：把空气、水和浆液一起考虑，注浆初始阶段的驱动响应会更真实。
• 结构应力耦合：浆液压力导致裂隙开度变化，反过来又影响渗流场，这就是流固耦合的范畴，COMSOL的“固体力学 + Darcy定律”接口可以做。

老实说，第一次建双重介质模型时，我花了差不多两三个晚上才把压力分布跑到与现场趋势一致。最花时间的地方不是COMSOL操作，而是对“裂隙和基质到底怎么分配渗透率”的理解。一旦在脑子里面把“通道”和“储层”的分界线画清楚了，建模就像搭积木一样顺。建议新手从最简单的单裂隙二维模型入手，先把网格、物理场、求解整套流程跑通，再逐步加裂隙数量和物理过程。宁可模型简单、结果真实，也别一开始就搞一个“数字孪生”，最后被困在网格和收敛性里出不来。