以写促学| comsol大坝边坡稳定性分析

文摘 2024-09-20 13:22 辽宁

这个模拟花了很久学习，来来回回改了很多次

过程中感受到了严谨认真的重要性

一个小小的错误（没点相关域，少了一个小数点，多复制了一个空格）

就会导致结果出现很多错误，而且很难检查出来

科研挺好的😂，能帮人培养起来严谨的习惯

感谢小李同志的支持😁 (ง •_•)ง (ง •_•)ง (ง •_•)ง

▼

一、模型参数及分析原理介绍

边坡稳定性分析是预测各种荷载和环境条件下土壤的沉降、变形和滑动的重要技术。

本案例来源于COMSOL中岩土力学模块土壤部分的堤坝边坡稳定性分析。

具体模型规格如下所示：大坝的高度L4为12m，水位的高度Hw为10m，渗流高度为4m。大坝的全部宽度L1+L2+L3为12m+5m+12m=29m。实际上，在建立模型时，我们会把堤坝下方再加一个土壤与以避免边界效应。

Fig.1 大坝简图

本例采用平面应变近似方法处理模型，通俗地说，平面应变近似方法就是“以面代体，默认将一个面的情况无限复制就得到整体的情况”“‘平面应力问题’即有一个方向应力为零，‘平面应变问题’即有一个方向应变为零。”

在这个案例中，重力和静水压力的影响也要考虑在内。除此之外，土壤的孔隙压力用达西定律来建模，莫尔库伦准则用于弹塑性分析。采用的方法为强度降低法(strength reduction method)，莫尔库伦材料参数为安全系数（FOS）的函数，安全系数不断增加，材料参数和土体的抗剪强度不断降低最终导致土壤不稳定。在给定荷载的情况下，安全系数不断增大的过程中，路堤就会发生坍塌。关于莫尔库伦参数与安全系数的函数如下所示：

莫尔库伦屈服函数和相关的塑性势为：

其中

参数化的内聚力和内摩擦角分别为：

其中，c是材料内聚力，其他两个u和s分别为非饱和土和饱和土的内摩擦角，p是达西流给出的孔隙压力。

二、建模及研究过程

为了考虑孔隙压力和重力荷载对斜坡稳定性的影响，我们进行了三项静态研究，分三步走。

第一、仅计算达西定律来求解孔隙压力曲线。

第二、根据静水压力、孔隙压力和重力荷载计算堤坝的初始应力。

第三、通过添加初始应力和应变节点来考虑这些初始应力。

饱和土和非饱和土的内摩擦角不同，因此使用了取决于孔隙压力的表达式。在非饱和土壤区域不施加外部应力，因为孔隙被认为是相互连接的，并且处于恒定的大气压力下。

为了为保持求解的网格客观性，在材料模型中引入了长度标度模型，将应变局部化限制在预定宽度内。这里使用的是隐式梯度非局部和长度标度 lint = 0.1 米的隐含梯度非局部塑性模型。

2.1建立二维模型和输入参数

建立二维模型，选择达西定律和固体力学进行稳态研究。

在全局定义中输入各种参数，然后定义插值函数并更改函数名称为cond，在下方输入单位，变元t的单位是m，函数cond的单位是m/s。

在有限元分析中，插值函数可用于估算单元内部的变量，如应力和位移。这些内部值通过插值函数在元素节点值之间进行估算，提高了计算的精度。

Fig.2 输入参数

创建二维几何体，在几何中添加多边形，输入数据（如果要复制已有数据，需要多少选多少即可，选中空白栏也会导致数据错误，出错过一次）。这个时候需要添加圆角消除几何不连续性，在几何当中选择倒圆角，选择需要研究的角，即角4，将倒角半径设置为5，点击构建所有对象。

Fig.3 构建多边形

在水库水位和可能的渗流水位处增加点，以分隔大坝两侧，即定义两个点，分别是水位高度在大坝上的那个点和渗流高度在大坝上的那个点，点击构建选定对象即可看到图上的两个点。

Fig.4 增加水位点和渗流点

定义以下若干变量，这个变量在上文中有描述，解释了安全系数和莫尔库伦材料函数之间的关系。

在主菜单定义中的非局部耦合当中选择最大值，选择域1。

Fig.5 定义变量

2.2定义土壤和水的材料参数

下面，我们要添加两个空白材料，一个用于土壤，一个用于水，然后相应地重命名它们。对于水材料，域选择保持为空(如下图)。（注意土体的选择不要为空，此处出错造成多次返工，且检查未发现）

在全局定义当中添加空材料，分别命名为Soil和Water。

在组件材料中选择更多材料当中的多孔材料，首先在均质材料中选择多孔材料中的土壤，然后右键单击多孔材料，选择流体，在材料栏选Water。右键单击多孔材料，选择固体，选择Soil，在固体中的体积分数一栏中改成1-psi。

Fig.6 定义材料

设置达西流，在达西流的多孔基体当中的渗透率模型中选择水力传导率，并将水里传导率改成K。

Fig.7 更改水力传导率

在全局定义中输入材料的值。土壤的密度为rho_soil，杨氏模量为E_soil，泊松比为nu_soil。水的密度为rho_wat。然后在大坝的下游和上游的水下部分增加压力水头，并且加上重力，将重力改为高度。

Fig.8 增加压力水头和重力

2.3划分网格

接下来，划分网格，单击网格1，在序列类型中选择用户控制网格。单击大小，在预定义栏中选择合适的细化程度。右键单击网格1，单击大小表达式，并且把大小一栏的表达式拖到大小的下面，在大小表达式中输入：

if(Y>-L4&&X>L1/2&&X<(L1+L2+L3*1.5),0.75,10)

在每个维度的单元数中填入数字50，点击全部构建。

Fig.9 网格划分情况

2.4求解孔隙水压

首先是达西流，在求解1中标明达西流，取消生成默认绘图，然后从研究中禁用固体力学，只求解孔隙压力。点击步骤1，稳态，关闭固体力学的内容。开始计算。

添加二维绘图组，命名为压力水头。右击压力水头，选择等值线。在表达式右上角选择替换表达式，选择图示内容，定义方法改为水平，输入图示内容，然后将等值线类型改为管，半径比例因子勾选并且改为0.1。

Fig.10 堤坝内的压力水头

得到饱和土和非饱和土的分界线：这张图片显示了堤坝内的压力水头，左侧水下墙的压力水头从0到10m不等，渗流面的压力水头为0m。正压头表示正孔隙压力，意思是土壤饱和，非饱和土壤则为零压头。所以零压力水头线表示饱和土和非饱和土的分界线。随后，随着塑性应变的增加和剪切强度的降低，斜坡会坍塌。

2.5根据静水压力、孔隙压力和重力荷载计算堤坝的初始应力

将水压作为边界载荷加在大坝下游侧，在固体力学当中选择边界荷载，如下图所示的两条边界加上压力，并且输入p。选择固定约束，选择最下面的一条边。选择辊支撑，选择两侧的两条边。

Fig.11 添加边界荷载等支撑

在线弹性材料1中的正交设置里勾选减缩积分，加速模拟的计算速度。在物理场的属性中选择外部应力，更改为孔隙压力，更改如图所示的参数（图中未标注Pref为0）。

Fig.12 添加孔隙水压为外部应力

继续添加另一项研究，计算重力和孔隙压力产生的原位应力。

添加研究2，更名为现场应力初始化，然后取消默认绘图。点击计算。在线弹性材料中添加土壤塑性，选择莫尔库伦准则，在塑性势中选择关联，然后再非局部塑性模型中选择隐式梯度，长度尺度选择0.1。

Fig.13 土壤塑性采用莫尔库伦准则

将第二个研究步骤中计算的应力作为初始应力加入，于是在线弹性材料中选择初始应力和应变，用withsol运算符号提取相应数据，如图所示。

Fig.14 提取初始应力和应变

为莫尔-库仑模型添加材料属性。在全局定义的土壤材料中添加内聚力和内摩擦角的参数，如图所示。

Fig.15 为莫尔库伦准则添加材料属性

禁用第二个研究步骤中的土壤塑性以及初始应力和应变节点。为弹塑性分析添加第三个研究步骤。

在第二步的研究中，单击步骤1稳态，在物理场和变量选择中勾选修改步骤的模型配置。禁用达西定律，禁用土壤塑性和初始应力应变。展开因变量值的不求解变量值，输入相关参数即可。

Fig.16 在第二步中禁用三个选项

添加第三个研究。命名为固体力学安全系数，取消生成默认绘图。在这一研究步骤中不求解孔隙压力变量，而是从第一次求解中获取其值，对参数 FOS 进行辅助扫描。取消勾选达西定律。在不求解变量的值中选择同样的内容，在研究扩展中勾选辅助扫描，参数名称为安全系数，参数值列表如图。

Fig.17 在第三步中禁用达西定律

三、结果分析

3.1斜坡坍塌前的滑动面及土颗粒位移

添加二维绘图组，标签改为滑动表面。在标题栏中把标题类型改为标签。要显示滑移面，则需要自定义等值线图。

右键单击滑动表面选择等值线，在替换表达式中选择固体力学、位移中的固体位移。把水平中的定义方法改为水平，在水平栏中输入0 0.1 Inf。在着色和样式中，将等值线改为填充，颜色表选择Linear当中的GrayPrint。然后把颜色表转换改为反转，取消颜色图例的勾选。点击绘制。

Fig.18 土颗粒位移方向及滑移圆

此图显示了土颗粒的位移方向，说明了土壤的滑动现象，右下角的土壤不滑动，因为下边有固定约束。

3.2斜坡坍塌前的等效塑性应变

再添加二维绘图组，命名为等效塑性应变。右键选择表面，在替换表达式中选择高斯点里的solid.epeGp，点击绘制。

Fig.19 坍塌前的等效塑性应变

这张图片显示了坍塌前的等效塑性应变。

3.3最大位移与 FOS 的关系

添加一维绘图组，改为安全系数，选择研究3。右键单击安全系数，选择全局，在y轴数据中输入max……，在x轴数据中改为表达式，最后找到FOS有关的式子选择。在着色和样式中的线标记中选择循环。

Fig.20 添加一维绘图组

单击安全系数，在一维绘图组的图例部分选择左上角，点击绘制。

Fig.21 FOS与最大位移的关系

此图显示的是FOS与最大位移的关系，最大位移在FOS=1.9左右明显增加，这表明斜坡在此时开始坍塌。

3.4堤坝坍塌前的位移量级

创建拉伸数据集，用于三维位移场可视化。在结果中选择数据集，选择二维拉伸。选择研究3，最大值选择L1+L2+L3，Z变量填写Z。在将平面映射到一栏中选择xz平面。

在结果中选择三维绘图组，命名为位移，右键单击选择表面。在替换表达市长选择固体力学中的位移，固体位移，把单位改为mm。单击绘制。

Fig.22 大坝位移的三维视图

这张图借助拉伸数据集显示了位移场的三维可视化。

ULS AGV DMT

一群为了自己幸福和他人幸福而奋斗的朋友和同志，记录我们的生活和思考。