Biot固结问题的有限元求解

学术科技 2024-10-24 22:22 山东

点击上方「蓝字」关注我们

近年来，多孔介质流固耦合问题越来越受到人们的重视。它在很多工程领域都有广泛的应用，如软土地基固结沉降，地下流体开采所导致的地面沉降，煤层气的耦合渗流和突出，边坡和坝基的稳定性问题，城市垃圾填埋及核废料处理，生物体软组织变形研究，等等。多孔介质流固耦合理论源自 Biot 的饱和土固结理论。本文以二维饱和土体的 Biot 固结问题为例，介绍有限元方法和 FEtch 系统在多孔介质流固耦合问题求解中的应用。

控制方程

微分形式

二维直角坐标系下，不考虑固体骨架和孔隙流体的压缩性，忽略惯性力和体力的作用，准静态的 Biot 方程可写成如下较为简单的形式。

平衡方程

\frac{\partial \sigma_{x x}}{\partial x}+\frac{\partial \sigma_{x y}}{\partial y}-\frac{\partial p}{\partial x} = 0 \ \left(\text{in } \Omega\right) \tag{1}

\frac{\partial \sigma_{x y}}{\partial x}+\frac{\partial \sigma_{y y}}{\partial y}-\frac{\partial p}{\partial y} = 0\ \left(\text{in } \Omega\right) \tag{2}

\frac{\partial \epsilon}{\partial t}-\frac{\partial}{\partial x}\left(\frac{k}{\gamma_w} \frac{\partial p}{\partial x}\right)-\frac{\partial}{\partial y}\left(\frac{k}{\gamma_w} \frac{\partial p}{\partial y}\right)=0 \ \left(\text{in } \Omega\right) \tag{3}

几何方程

\begin{aligned} \varepsilon_{xx}&=\frac{\partial u}{\partial x},\\ \varepsilon_{yy}&=\frac{\partial v}{\partial y},\\ \varepsilon_{x y}&=\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x},\\ \epsilon&=\varepsilon_{xx}+\varepsilon_{yy} \end{aligned}

本构方程（平面应变）

\left(\begin{array}{c}\sigma_{x x} \\ \sigma_{y y} \\ \sigma_{x y} \end{array}\right)=\frac{E}{(1+\nu)(1-2\nu)}\left(\begin{array}{ccc} 1-\nu & \nu & 0 \\ \nu & 1-\nu & 0 \\ 0 & 0 & 0.5-\nu \end{array}\right)\left(\begin{array}{l} \varepsilon_{x x} \\ \varepsilon_{y y} \\ \varepsilon_{x y} \end{array}\right)

以上式子可以简记为

\nabla\cdot(\boldsymbol{\sigma}-p\boldsymbol{I})=\boldsymbol{0}\ \left(\text{in } \Omega\right)

\frac{\partial \epsilon}{\partial t}-\nabla\cdot \frac{k}{\gamma_w}\nabla p=\boldsymbol{0} \ \left(\text{in } \Omega\right)

\boldsymbol{\varepsilon}=\left(\nabla \boldsymbol{u}+\nabla^{T} \boldsymbol{u}\right)/2,\quad\epsilon =\boldsymbol{\varepsilon}:\boldsymbol{I}

\boldsymbol{\sigma}=\boldsymbol{D} \boldsymbol{\varepsilon}

其中，、为位移的分量，为孔压（以压为正）。、、为应变的分量，为体应变，、、为有效应力 的分量。为梯度算子，为刚度矩阵，为单位矩阵。参数为弹性模量，为泊松比，为渗透系数，为水的容重。、为直角坐标系下的坐标分量。

边界条件

考虑两类边界条件:

第一类边界条件，给定位移和孔压

u=u_0,\ v=v_0,\ p=p_0\ \left(\text{on } \Gamma_1\right)

第二类边界条件，给定外力和流量

\begin{aligned} T_{x}=\left(\sigma_{x x}-p,\sigma_{x y}\right)\cdot \mathbf{n}=t_x,\\ T_{y}=\left(\sigma_{x y},\sigma_{y y}-p\right)\cdot \mathbf{n}=t_y,\\ \frac{k}{\gamma_w}\frac{\partial u}{\partial \mathbf{n}} \cdot \mathbf{n}=q_0\ \left(\text{on } \Gamma_2\right) \end{aligned}

其中，、为边界上的位移，为边界上的孔压；、分别为边界上的外力在和方向的分量，为边界上的流量（以流入为正）。为边界上的单位外法向量。

弱形式

有限元模型的建立依赖于待解偏微分方程的弱形式。根据变分原理，对力学平衡方程 (1) 和 (2) 两边分别乘以位移的变分，并在求解域内积分，得到其等效积分形式

\begin{aligned} &\int_{\Omega}\left(\frac{\partial \sigma_{x x}}{\partial x}+\frac{\partial \sigma_{x y}}{\partial y}-\frac{\partial p}{\partial x}\right) \delta u d\Omega\\ +&\int_{\Omega}\left(\frac{\partial \sigma_{x y}}{\partial x}+\frac{\partial \sigma_{y y}}{\partial y}-\frac{\partial p}{\partial y}\right) \delta v d\Omega=0 \end{aligned}

对上式进行分部积分，并代入边界条件，得到弱形式

\begin{aligned} \int_{\Omega}\left(\sigma_{x x} \delta \varepsilon_{x x}+\sigma_{y y} \delta \varepsilon_{y y}+\sigma_{x y} \delta \varepsilon_{x y}\right)d\Omega \\ -\int_{\Omega}p\delta \epsilon d\Omega =\int_{\Gamma_2}\left(t_{x} \delta u+t_{y} \delta v\right) d\Gamma \end{aligned}

将本构方程代入上式，即可得到以位移为基本未知量的弱形式。

对渗流方程 (3) 两边分别乘以孔压的变分并在全域上积分，得

\int_{\Omega}\left[\frac{\partial \epsilon}{\partial t}-\frac{\partial}{\partial x}\left(\frac{k}{\gamma_w} \frac{\partial p}{\partial x}\right)-\frac{\partial}{\partial y}\left(\frac{k}{\gamma_w} \frac{\partial p}{\partial y}\right)\right] \delta p d \Omega=0

对于二阶导数的部分，利用分部积分将其中的一阶导数转化到孔压的变分上，并继续代入边界条件，得

\begin{aligned} &\int_{\Omega}\frac{\partial \epsilon}{\partial t} \delta p d \Omega+\int_{\Omega}\left(\frac{k}{\gamma_w} \frac{\partial p}{\partial x} \frac{\partial \delta p}{\partial x}+\frac{k}{\gamma_w} \frac{\partial p}{\partial y} \frac{\partial \delta p}{\partial y}\right) d \Omega \\ =&\int_{\Gamma_2} q_{0} \delta p d \Gamma \end{aligned}

这就是孔压场的最终的弱形式。

离散化

对以上弱形式的力学平衡方程和渗流方程进行空间离散化，可进一步写成矩阵形式

M \dot{U}+SU=F

其中，为位移和孔压的组合向量，为组合向量对时间的导数向量，为质量矩阵，为刚度矩阵，为载荷向量。采用向后差分法进行时域离散，可将上式改写为

M \frac{U^{t+\Delta t}-U^{t}}{\Delta t}+S U^{t+\Delta t}=F^{t+\Delta t}

整理得

(M+S \Delta t) U^{t+\Delta t}=F^{t+\Delta t} \Delta t+M U^{t}

为时间步长，上标代表当前时间步，代表上一时间步。该式即为最终需要求解的线性方程组。

算例1

Tezaghi 固结

Tezaghi 固结问题是土力学的经典算例。考虑厚度为的饱和含水层，水平方向为无限长。底面为不透水边界，顶面为渗流自由面，孔隙压力；含水层材料为各向同性，其弹性模量，泊松比，渗透系数，水的容重为。顶面突然施加的均布载荷。不计重力影响。求 1 小时后土体的变形情况和孔压分布。

前处理

进入 GiD 前处理，建立一个的长方形区域，它的两个对角的角点坐标分别为（0, 0）和（0.2, 1）。网格剖分采用 8 节点等参单元。共使用了 40 个单元，165 个节点。

初始时间步长为 1 s，结束时间为 3600 s，时间步长放大因子为 1.05。

需要注意的是，Biot 问题的数值方法对时间步长的下限有要求，不能太小，否则容易发生数值震荡现象。具体理论请查看相关文献。

计算结果

通过自动改变时间步长，变化范围从 1 s 到 167 s，总共使用了 107 个时间步。对计算结果稍加整理，效果如下。

竖向位移的演化过程

孔压的演化过程

不同深度处位移随时间的变化曲线

土体表面（蓝线）和 0.5 m深处（红线）的位移情况如下图所示，反映了外力作用下土体逐渐压缩的过程。

不同高度处孔压随时间的变化曲线

选取深度为 0.1、0.3、0.5、0.7 和 0.9 m 的位置作为观察点，孔压随时间变化的情况如下图所示。

进一步分析可以发现数值解与解析解符合得很好，进而证明算法和程序的有效性。

算例2

圆柱土体固结

一个半径为的无限长圆柱土体，在其外边界突然施加的均匀径向压缩载荷。圆柱体表面为渗流自由面，孔隙压力。材料参数取值同算例 1。不计重力的影响。求 1 小时后土体的位移和孔压分布。

前处理

进入 GiD 进行前处理。根据圆截面的对称性，取四分之一的区域进行模拟。网格剖分采用 8 节点等参单元。共使用了 300 个单元，961 个节点。

计算结果

对计算结果稍加整理，效果如下。

径向位移的演化过程

径向位移，随时间的变化过程如下图所示。

孔压的演化过程

径向位移随时间的变化曲线

土体表面（红线）和 0.5 m 深处（蓝线）的径向位移情况如下图所示，反映了土体整体逐渐压缩的过程。

孔压随时间的变化曲线

选取径向距离 0.0、0.4 和 0.8 m 的三个位置作为观察点，孔压随时间变化的情况如下图所示。

可以发现，较深位置的孔压有先上升后下降的现象，这就是著名的 Mandel-Cryer 效应。

进一步分析可以发现数值解与解析解符合得很好，进而证明算法和程序的有效性。

往期推荐

一图读懂 FEtch 的工作流程

谈谈有限元程序开发，兼谈 FEtch 入门

通用前后处理软件 GiD 简介

Biot固结问题的有限元求解

控制方程

微分形式

弱形式

离散化

算例1

前处理

计算结果

算例2

前处理

计算结果

推荐阅读