土—岩界面优势流潜蚀数值模拟 [1]模型简介：使用数值模拟软件COMSOL，复现EI论文(豆红强,谢森华,简文彬,等.球状风化花岗岩类土质边坡土-岩界面优势流潜蚀特性研究[J].岩土力学,2024,45(04):950-960.) [2]案例内容：完整数值模型一个（包括模型边界条件设置、云图结果、后处理数据等），理论方法及具体操作讲解全过程 [3]模型特色：（1）基于连续介质力学，推导出因渗流侵蚀导致的孔隙率演变方程、可侵蚀细颗粒含量方程、液化细颗粒浓度方程；（2）采用VG模型描述土体非饱和特性，采用BC模型描述裂隙优先流特征，并采用双重介质降雨边界/空气单元边界；（3）使用动态渗透系数方程描述侵蚀后的渗透系数因孔隙率的变化；（4）采用Cividini and Gioda（2004）土壤细颗粒侵蚀方程描述细颗粒的侵蚀作用 COMSOL为6.2版本

在岩土工程领域，土—岩界面优势流潜蚀现象的研究对于理解边坡稳定性等问题至关重要。今天咱们就借助COMSOL 6.2软件，来复现一篇EI论文中关于球状风化花岗岩类土质边坡土 - 岩界面优势流潜蚀特性的研究。

模型简介

本次模拟依据豆红强、谢森华等学者发表在《岩土力学》2024年第45卷第04期的论文 “球状风化花岗岩类土质边坡土 - 岩界面优势流潜蚀特性研究”。COMSOL作为强大的多物理场仿真软件，能很好地助力我们实现模型的复现。

模型特色

1. 关键方程推导：基于连续介质力学，这模型推导出了几个关键方程。
   - 孔隙率演变方程：它描述了因渗流侵蚀，土体孔隙率是如何变化的。假设我们用 $\phi$ 表示孔隙率，随着渗流侵蚀的进行，孔隙率的变化可能涉及到流速、颗粒冲刷等因素。例如简化的形式可能是 $\frac{d\phi}{dt} = f(v, C)$，其中 $v$ 是渗流速度，$C$ 是与颗粒冲刷相关的参数，$t$ 为时间。这个方程反映了孔隙率随时间的动态演变，而在COMSOL中，我们可以通过偏微分方程模块来定义类似这样的方程。
   - 可侵蚀细颗粒含量方程：用来追踪土体中可被侵蚀的细颗粒数量变化。设可侵蚀细颗粒含量为 $m$，方程可能像 $ \frac{dm}{dt} = -k \cdot v \cdot m$，$k$ 是一个与侵蚀能力相关的系数，负号表示随着渗流，细颗粒含量在减少。在COMSOL里，通过定义源项或通量来体现这个方程对模型的影响。
   - 液化细颗粒浓度方程：关乎液化后细颗粒在土体中的浓度分布。用 $C{l}$ 表示液化细颗粒浓度，方程形式或许是 $\nabla \cdot (D \nabla C{l}) - v \cdot \nabla C_{l} = S$，$D$ 是扩散系数，$S$ 是源项，在COMSOL中，类似这种扩散 - 对流方程可以在相关的物理场接口中设置。
2. 特性描述模型：
   - VG模型：用来描述土体非饱和特性。在COMSOL中，通过在多孔介质流物理场接口中，选择合适的参数化方式来体现VG模型。例如，对于土 - 水特征曲线，我们可以根据VG模型的参数来设置相对渗透率和饱和度之间的关系。
   - BC模型：用于刻画裂隙优先流特征。在模型构建时，对于代表裂隙的区域，通过设置特殊的渗透率张量等参数来模拟BC模型所描述的优先流现象。同时采用双重介质降雨边界/空气单元边界，在COMSOL里，我们可以通过多物理场耦合，分别定义降雨边界条件和空气单元边界条件来实现。比如在降雨边界，设置降雨强度随时间的变化，在空气单元边界设置气压等相关条件。
3. 动态渗透系数方程：为了描述侵蚀后渗透系数因孔隙率的变化，使用动态渗透系数方程。假设渗透系数 $k$ 与孔隙率 $\phi$ 的关系为 $k = k0 \cdot \phi^n$，$k0$ 是初始渗透系数，$n$ 是一个经验指数。在COMSOL中，我们可以在材料属性设置中，将渗透系数定义为孔隙率的函数，随着孔隙率在模拟过程中的变化，渗透系数也相应改变。
4. 细颗粒侵蚀方程：采用Cividini and Gioda（2004）土壤细颗粒侵蚀方程描述细颗粒的侵蚀作用。在COMSOL里，将这个方程以源项或通量的形式添加到相关的物理场方程中，以准确模拟细颗粒在渗流作用下的侵蚀过程。

案例内容 - 完整数值模型搭建

1. 模型边界条件设置：在COMSOL中创建几何模型后，就需要设置边界条件。比如在边坡的底部，设置为固定边界条件，防止模型在重力方向上的移动。代码实现可能类似这样：

model.geom(1).feature('bnd1').set('name', 'Bottom Fixed', 'geom_id', [1 2 3]); % 假设底部边界几何ID为1、2、3
model.physics('solid').bnd(1).set('ux', 0, 'uy', 0, 'uz', 0); % 在结构力学物理场中设置位移为0

这段代码首先给底部边界命名为 “Bottom Fixed”，然后在结构力学物理场中，将该边界在x、y、z方向的位移都设为0，实现固定边界条件。

土—岩界面优势流潜蚀数值模拟 [1]模型简介：使用数值模拟软件COMSOL，复现EI论文(豆红强,谢森华,简文彬,等.球状风化花岗岩类土质边坡土-岩界面优势流潜蚀特性研究[J].岩土力学,2024,45(04):950-960.) [2]案例内容：完整数值模型一个（包括模型边界条件设置、云图结果、后处理数据等），理论方法及具体操作讲解全过程 [3]模型特色：（1）基于连续介质力学，推导出因渗流侵蚀导致的孔隙率演变方程、可侵蚀细颗粒含量方程、液化细颗粒浓度方程；（2）采用VG模型描述土体非饱和特性，采用BC模型描述裂隙优先流特征，并采用双重介质降雨边界/空气单元边界；（3）使用动态渗透系数方程描述侵蚀后的渗透系数因孔隙率的变化；（4）采用Cividini and Gioda（2004）土壤细颗粒侵蚀方程描述细颗粒的侵蚀作用 COMSOL为6.2版本

对于降雨边界，在多孔介质流物理场中设置：

model.physics('pfd').bnd(4).set('type', 'rainfall', 'rainfall', 0.001); % 假设降雨边界ID为4，降雨强度设为0.001 m/s

这里将边界类型设为 “rainfall”，并设置降雨强度为0.001 m/s 。

1. 云图结果：模拟完成后，我们可以生成各种云图来直观展示结果。例如，孔隙率云图可以帮助我们看到整个模型中孔隙率的分布情况。在COMSOL的后处理模块中，选择 “派生值” - “表面绘图组”，然后选择要绘制云图的变量为 “孔隙率 $\phi$”。通过调整颜色映射和等高线设置，得到清晰的孔隙率云图，如下代码可实现部分后处理操作：

model.result.create('surf1', 'Surface');
model.result('surf1').feature('pg1').set('expr', 'phi'); % 假设孔隙率变量名为phi
model.result('surf1').run;

上述代码创建了一个表面绘图组 “surf1”，并将其表达式设为孔隙率变量 “phi”，最后运行该绘图组得到孔隙率云图。

1. 后处理数据：除了云图，还可以提取具体的数据进行分析。比如提取某一点的孔隙水压力随时间的变化。在COMSOL中，通过 “派生值” - “点求值”，选择要分析的点以及物理场变量 “孔隙水压力 $p_w$”，就可以得到该点孔隙水压力随时间的变化数据。同样可以用代码实现自动化提取：

model.result.create('point1', 'Point');
model.result('point1').feature('ev1').set('expr', 'p_w', 'pos', [0.5 0.5 0.5]); % 假设在坐标[0.5, 0.5, 0.5]处提取孔隙水压力
data = model.result('point1').evaluate;

这段代码创建了一个点求值对象 “point1”，在坐标 [0.5, 0.5, 0.5] 处提取孔隙水压力变量 “p_w” 的值，并将结果存储在 “data” 中以便后续分析。

通过以上在COMSOL 6.2中的操作，我们就能完整复现土 - 岩界面优势流潜蚀数值模型，从理论方法到具体操作，深入探究这一复杂的岩土工程现象。