承压水降压引起上覆层状土体二维黏弹性固结分析

Two-Dimensional Viscoelastic Consolidation Analysis of Overlying Layered Soil Caused by Pressure Drop of Confined Water

摘要: 针对承压层的水位下降引起的上覆多层结构土体沉降问题，结合二维固结方程、Merchant流变模型与水位变化边界条件，提出了考虑土骨架黏弹性的多层弱透水层二维渗流固结数学模型。针对该模型，由Fourier-Laplace变换推导了Fourier-Laplace变换域内问题的解析解，利用Stehfest数值反演方法得到了问题在在时域内的真实解。在此基础上，对本文给出的多层弱透水层固结度、沉降量解答与已有一维解析解、二维有限单元数值解进行了对比，验证了本文方法的正确性。通过进一步的算例分析，研究了土体黏滞性、刚度比、以及承压层水位下降幅度对固结过程的影响规律。

Abstract: For the settlement problem of the overlying multi-layer soil caused by the water level drawdown in the confined layer, a two-dimensional seepage consolidation mathematical model of the multi-layer aquitard considering the viscoelasticity of the soil skeleton is proposed in combination with the two-dimensional consolidation equation, the Merchant rheological model and the boundary condition of water level change. For this model, the analytical solution of the problem in the Fourier Laplace transform domain is derived from the Fourier Laplace transform, and the real solution of the problem in the time domain is obtained using the Stehfest numerical inversion method. On this basis, the consolidation degree and settlement of multi-layer aquitard given in this paper are compared with the existing one-dimensional analytical solution and two-dimensional finite element numerical solution, and the correctness of this method is verified. The influence of viscosity, stiffness ratio and amplitude of the water level drawdown on consolidation process is studied through further analysis of several examples.

1. 引言

地面沉降发展时间长、影响范围广、治理难度高，造成区域性社会经济损失难以逆转 [1] [2] 。大量研究表明，引发区域性地面沉降的重要因素之一是地下水资源的无节制开采 [3] [4] [5] 。因此，研究地下水位变化作用下土层的固结变形具有重要的理论和现实意义。

传统外荷载作用下土体发生的固结沉降问题已得到较为详细深入的研究 [6] [7] [8] [9] [10] ，而针对水位变化引发的土体固结沉降问题的相关研究相对较少 [11] 。Tseng等 [12] 提出了含水层水位下降恒定值时，含水层一维弹性固结的解析解。Tao和Xie [13] 针对含水层中的瞬时抽水情况，提出了类似的一维弹性固结的计算方案。Liu等 [14] 考虑了弱透水层的黏弹性，提出了水位变化下的一维单层固结的半解析解。Yang等 [15] 用Galerkin有限单元法研究水位变化作用下的土层沉降问题，并通过移动网格法来提升所得到数值解的精度。Hall和Fox [16] 开发了一个一维大应变数值模型，用于地下水位下降引起的固结问题。Lo等 [17] 推导了具有流体通量边界条件的非饱和土体的一维弹性解。

然而，上述解析研究主要针对一维单层弹性固结问题，为了更真实地反映地面沉降过程，研究多维、多层土体的黏弹性固结问题具有一定的现实意义。为此，本文建立承压层水位变化下的多层弱透水层的二维渗流固结数学模型，利用Laplace-Fourier变换方法推导了该问题的在变换空间的解析解，利用Stechfest反演算法，得出地表沉降和固结度在时间域内的半解析解答。对比已有一维解析解来验证该解答的正确性，结合具体算例分析土体力学参数及水文地质参数对固结过程的影响规律。

本文结构安排如下：第二节介绍含水层系统模型，第三节进行数学求解，第四节通过算例对比验证方法可靠性，第五节进行地面沉降影响因素分析，第六节为本文总结。

2. 问题描述

2.1. 层状弱透水层系统

假设某开采承压含水层，上方覆盖有多层的弱透水层，如图1。上覆多层弱透水层的总厚度设为Ms，其中M1、M2、∙∙∙、Mn分别为自上而下各弱透水层的厚度，Mc为承压层厚度，含水层系统内初始水头均为H， $\Delta h\left(x,t\right)$ 为t时刻承压含水层内的水头变化值。

Figure 1. Layered aquitard model under water level change

图1. 水位变化下层状弱透水层系统示意图

2.2. 控制方程

假设土颗粒、孔隙水无压缩性，地下水的运动符合Darcy定律，则根据流量守恒原理，弱透水层的固结控制方程如下：

$\frac{k_m}{{\gamma }_w}\left(\frac{{\partial }^2{h}_m}{\partial z^2}+\frac{{\partial }^2{h}_m}{\partial x^2}\right)=-\frac{\partial {\epsilon }_v{}_m}{\partial t}$ (1)

式中：h_m是第m个弱透水层水头变化值(1 ≤ m ≤ n)，z为该位置距离弱透水层顶面的垂直距离，x为该位置距离水位变化中心的水平距离，t为含水层系统内固结进行的时间，k_m为第m个弱透水层的渗透系数，ε_vm为第m个弱透水层内单元体的体积应变，H为承压含水层初始水头，γ_w是水的重度，取10 kN/m³。

弱透水层土体的流变性可以Merchant模型描述(见图2)，如果模型中胡克弹簧的弹性模量分别为E_0m和E_1m (单位为MPa)，牛顿黏壶的黏滞系数为η_1m (单位为MPa∙d)，则该模型应力应变关系如下式：

$\epsilon =\frac{{\sigma }_0}{E_0}+{\displaystyle {\int }_0^t\frac{{\sigma }_0}{{\eta }_1}{\text{e}}^{-\frac{E_1}{{\eta }_1}\left(t-\tau \right)}}\text{d}\tau$ (2)

将该式代入(1)，可得出基于Merchant流变模型的弱透水层控制方程如下：

$\frac{k_m}{{\gamma }_w}\left(\frac{{\partial }^2{h}_m}{\partial z^2}+\frac{{\partial }^2{h}_m}{\partial x^2}\right)=\frac{\partial h_m}{\partial t}{G}_m\left(0\right)+{\displaystyle {\int }_0^t\frac{\text{d}{G}_m\left(t-\tau \right)}{\text{d}\left(t-\tau \right)}}\frac{\partial h_m}{\partial \tau }\text{d}\tau$ (3)

Figure 2. The mth rheological model of aquitard layer

图2. 第m个弱透水层流变模型示意

2.3. 定解条件

假设弱透水层系统内初始无水头变化，顶部边界自由透水，下部边界受承压含水层水位变化影响，该多层系统的初始条件及上、下边界条件如下：

$h_m|{}_{t=0}=0$ (4)

$h_1|{}_{\text{z}=0}=0$ (5)

$h_n|{}_{z=M_s}=\Delta h\left(x,t\right)$ (6)

在第 $m-\text{1}$ 与m个土层的交界处(2 ≤ m ≤ n)，还应该满足如下层间连续条件：

$h_{\left(m-1\right)}|{}_{z={\displaystyle \underset{l=1}{\overset{m-1}{\sum }}{M}_l}}=h_m|{}_{z={\displaystyle \underset{l=1}{\overset{m-1}{\sum }}{M}_l}}$ (7)

$k_{\left(m-1\right)}\left(\frac{\partial h_{(m-1)}}{\partial z}+\frac{\partial h_{(m-1)}}{\partial x}\right)|{}_{z={\displaystyle \underset{l=1}{\overset{m-1}{\sum }}{M}_l}}=k_m\left(\frac{\partial h_m}{\partial z}+\frac{\partial h_m}{\partial x}\right)|{}_{z={\displaystyle \underset{l=1}{\overset{m-1}{\sum }}{M}_l}}$ (8)

3. 问题解答

下面将采用Laplace-Fourier变换对多层弱透水层黏弹性固结模型的控制方程和相应定解条件进行变换和求解，并通过相应的逆变换得到固结模型在真实物理域的解答。

对式(3)~(8)进行Laplace-Fourier变换可得：

$\frac{{\text{d}}^{\text{2}}{\hat{h}}_{\text{m}}}{{\text{dz}}^{\text{2}}}=I_m{\hat{h}}_{\text{m}}$ (9)

${\hat{h}}_{\text{m}}|{}_{t=0}=0$ (10)

${\hat{h}}_{\text{1}}|{}_{z=0}=0$ (11)

${\hat{h}}_n|{}_{z=M_{\text{s}}}=\alpha \left(w,s\right)$ (12)

${\hat{h}}_{\left(m-1\right)}|{}_{z={\displaystyle \underset{l=1}{\overset{m-1}{\sum }}{M}_l}}={\hat{h}}_m|{}_{z={\displaystyle \underset{l=1}{\overset{m-1}{\sum }}{M}_l}}$ (13)

$k_{(m-1)}\left(\frac{\partial {\hat{h}}_{(m-1)}}{\partial z}+iw{\hat{h}}_{(m-1)}\right)|{}_{z={\displaystyle \underset{l=1}{\overset{m-1}{\sum }}{M}_l}}=k_m\left(\frac{\partial {\hat{h}}_m}{\partial z}+iw{\hat{h}}_m\right)|{}_{z={\displaystyle \underset{l=1}{\overset{m-1}{\sum }}{M}_l}}$ (14)

式中 ${\hat{h}}_{\text{m}}\left(z\text{,}w\text{,}s\right)={\displaystyle {\int }_{\text{-}\infty }^{\infty }{\displaystyle {\int }_{\text{0}}^{\infty }{h}_{\text{m}}\left(z,x,t\right){\text{e}}^{-\text{i}wx}{\text{e}}^{-st}\text{d}t\text{d}x}}$ ， $\alpha \left(w,s\right)$ 是水位变换 $\Delta h\left(x,t\right)$ 的Fourier-Laplace变换，其中w、s分别为Fourier-Laplace域内的变换参数， $I_m=\frac{s{\gamma }_w}{k_m}\left(\frac{1}{s{\eta }_{1m}+E_{1m}}+\frac{1}{E_{0m}}\right)+w^2$ 。

常微分方程(9)的通解为：

${\hat{h}}_{\text{m}}=A_{\text{m}}{\text{e}}^{\sqrt{I_{\text{m}}}z}+B_{\text{m}}{\text{e}}^{-\sqrt{I_{\text{m}}}z}$ (15)

式中：A_m、B_m为待定系数，由各定解条件确定。

将变换后的边界条件式(10)~(14)代入式(15)，并整理成矩阵形式可得：

$\left[\begin{array}{ccccccc}V& & & & & & \\ L_1& F_1& & & & & \\ & \ddots & \ddots & & & & \\ & & L_{\text{m}}& F_m& & & \\ & & & \ddots & \ddots & & \\ & & & & L_{\text{n-1}}& F_{\text{n-1}}& \\ & & & & & & P\end{array}\right]\left[\begin{array}{c}{A}_1\\ B_1\\ ⋮\\ A_2\\ B_2\\ ⋮\\ A_{\text{n}}\\ B_{\text{n}}\end{array}\right]=\left[\begin{array}{c}0\\ ⋮\\ \alpha \left(w,s\right)\end{array}\right]$ (16)

式中：

$V=\left[\begin{array}{cc}1& 1\end{array}\right]$ ,

$L_m=\left[\begin{array}{cc}{\text{e}}^{\sqrt{I_m}{\displaystyle \underset{l=1}{\overset{m}{\sum }}{M}_l}}& {\text{e}}^{-\sqrt{I_m}{\displaystyle \underset{l=1}{\overset{m}{\sum }}{M}_l}}\\ k_m\sqrt{I_m}{\text{e}}^{\sqrt{I_m}{\displaystyle \underset{l=1}{\overset{m}{\sum }}{M}_l}}& -k_m\sqrt{I_m}{\text{e}}^{-\sqrt{I_m}{\displaystyle \underset{l=1}{\overset{m}{\sum }}{M}_l}}\end{array}\right]$ ,

$F_m=\left[\begin{array}{cc}-{\text{e}}^{\sqrt{I_{\left(m\text{+}1\right)}}{\displaystyle \underset{l=1}{\overset{m}{\sum }}{M}_l}}& -{\text{e}}^{\sqrt{I_{\left(m+1\right)}}{\displaystyle \underset{l=1}{\overset{m}{\sum }}{M}_l}}\\ -k_{\left(m\text{+}1\right)}\sqrt{I_{\left(m\text{+}1\right)}}{\text{e}}^{\sqrt{I_{\left(m\text{+}1\right)}}{\displaystyle \underset{l=1}{\overset{m}{\sum }}{M}_l}}& k_{(m+1)}\sqrt{I_{\left(m\text{+}1\right)}}{\text{e}}^{-\sqrt{I_{\left(m\text{+}1\right)}}{\displaystyle \underset{l=1}{\overset{m}{\sum }}{M}_l}}\end{array}\right]$ ,

$P=\left[\begin{array}{cc}{\text{e}}^{\sqrt{I_n}{M}_{\text{s}}}& {\text{e}}^{-\sqrt{I_n}{M}_{\text{s}}}\end{array}\right]$ .

通过式(16)求解待定系数，将其代入式(15)可以得到各弱透水层在Fourier-Laplace域内的解，再对其进行相应逆变换，可以得到多层土体固结模型在真实物理域域内的解答：

$h_m\left(z,x,t\right)=\frac{1}{4{\text{π}}^{2}i}{\displaystyle {\int }_{-\infty }^{\infty }{\displaystyle {\int }_{\alpha -i\infty }^{\alpha +i\infty }{\hat{h}}_m\left(z,w,s\right){\text{e}}^{iwx}{\text{e}}^{st}\text{d}s\text{d}w}}$ (17)

弱透水层系统顶部计算原点处的沉降可用下式计算：

$S_c={\displaystyle {\int }_0^{M_{\text{s}}}{\epsilon }_{\text{v}m}\left(z,0,t\right)\text{d}z}$ (18)

t时刻沉降量与最终沉降量可以定义随时间变化的弱透水层平均固结度：

$U\left(t\right)=\frac{S_c\left(t\right)}{S_c\left(\infty \right)}$ (19)

式中： $U\left(t\right)$ 为t时刻平均固结度。

4. 算例对比验证

式(17)给出了多层弱透水层的固结控制方程在时域内的解答。基于Stechfest逆变换，利用MATLAB编写了计算程序实现了该解答在真实物理域的计算，为了验证上述方法的可行性，通过具体算例与一维现有半解析解进行对比。

假定承压含水层中为均匀水位下降，此时可将本文解退化为一维问题的解。假设承压含水层中水位瞬时下降10 m，初始水头14 m，其上覆有两层弱透水层土体，弱透水层参数见表1。将本文解与已有一维固结解进行计算对比，结果见图3。可以看出，本文解与杨伟涛 [18] 的一维固结半解析解吻合良好。

Table 1. Calculation parameters for aquitard layer

表1. 弱透水层计算参数

Figure 3. Comparison with one-dimensional solution

图3. 本文解与一维解的对比

5. 影响因素分析

5.1. 土体黏滞性对固结发展的影响

本文使用Merchant流变模型来综合考虑土体在有效应力变化时的瞬时变形和蠕变变形。弱透水层土体黏滞系数取不同的值，取值情况见表2，其余参数同第三节验证算例。使用本文方法对弱透水层在不同黏滞系数下的固结度进行计算，结果见图4。可以看出，在土体其余因素不发生变化时，随着土体黏滞性的增加，其固结完成的时间分别增长了300 d和500 d，说明了在一定水位变化状态下，土体黏滞性的增大会大幅增加固结变形的完成时间，且这种效应在固结的中后期比较明显，表现出明显的变形滞后性。

Table 2. Parameter values for Viscosity Coefficient of aquitard Layers

表2. 弱透水层黏滞系数取值

Figure 4. Comparison of consolidation degree under the change of viscosity coefficient

图4. 粘滞系数变化下的固结度对比

5.2. 刚度比对固结发展的影响

无量纲参数刚度比 $b=E_{0m}/E_{1m}$ ，刚度比的大小影响了土体在应力变化时滞后变形所占的比值。对弱透水层土体的刚度比取不同的值，取值情况见表3，其余参数同第四节验证算例。使用本文方法对弱透水层固结时其顶部x = 1 m处不同刚度比下固结度进行计算，结果见图5。可以看出，随着刚度比的增大，水位变化时弱透水层的沉降发展将更为滞后，相应的固结完成时间就有所增长。

Figure 5. Comparison of consolidation degree under the change of stiffness ratio

图5. 刚度比变化下的固结度对比

Table 3. Parameter values for stiffness ratio

表3. 刚度比取值

5.3. 承压水位变化幅度对固结发展的影响

为了研究水位下降幅度对弱透水层固结沉降的影响，对验证算例中承压层降深取不同的值，取值情况见表4，其余参数同第三节验证算例。使用本文方法对不同排水量下的沉降量进行计算，结果见图6。可以看出，当水位降深由8 m增大至14 m时，弱透水层的累计沉降量增大了0.24 m，固结完成时间也相应增加，说明承压层水位变化幅度的增大会增加上覆弱透水层的最终沉降量及固结最终完成的时间。

Table 4. Variant values for confined groundwater drawdown

表4. 承压水降深幅度取值

Figure 6. Subsidence development under different drawdown

图6. 不同降深下的沉降发展

6. 结论

1) 将多层弱透水层固结模型与黏弹性模型相结合，考虑土体流变性，建立了水位变化作用下多层弱透水层二维的渗流固结模型，使用积分变换法得到了固结模型在变换域内的解析解，在此基础上使用了相应的数值逆变换法求得了真实物理域内的半解析解。

2) 使用MATLAB软件对求得的多层土体黏弹性固结解进行了编程计算，通过具体的算例将本文解与已有一维单层解进行了对比，验证了本文方法的正确性。

3) 结合具体算例进一步讨论了弱透水层土体黏滞性、刚度比和承压含水层的水位变化幅度对固结发展的影响，计算结果表明：弱透水层土体的黏滞性越大，在固结发展的中后期其固结速率会减小，相应的固结完成时间增加；而土体刚度比的增加会使土体固结发展更加滞后，同时，承压层水位下降幅度的增加会使弱透水层土体的应力变化增加，进而导致土体最终沉降量以及固结发展完成的时间有所增长。

4) 值得注意的是，实际情况中地下水位降深往往是一个关于时间t的函数，下一步可以在本文基础上进行地下水非瞬时降深引发的地基土固结沉降问题研究。