水位下降诱发含水层–弱透水层1维黏弹性固结分析

徐 进，杨伟涛，陈 征，王少伟

(1.烟台大学 土木工程学院，山东 烟台 264005；2.武汉大学 水工岩石力学教育部重点实验室，湖北 武汉 430072)

地面沉降是一种多发的严重环境地质问题，主要诱因是过量开采地下水引起的地下水位降深。如何对地下水位变化诱发的地面沉降灾害进行正确评估与科学防治一直备受学术界关注[1–4]。荷载作用下土体固结问题一直是土力学的传统课题，已得到较深入的研究。Gray[5]较早对瞬时荷载作用下双层地基固结问题进行了研究。栾茂田等[6]针对层状饱和土体，运用分离变量法求解了荷载作用下的Terzaghi 1维固结方程。谢康和[[7]及杨骏[8]等进一步考虑了变荷载作用下的类似问题。冯健雪等[9–10]考虑了连续排水边界条件下自重与线性加载时地基土的1维固结解析解。陈宗基[11]最先将流变模型引入固结分析中，刘加才等[12]基于广义Voigt流变模型，给出了荷载作用下双层黏弹性地基土的1维固结解析解。

与荷载作用不同，含水层中地下水开采引起的水位下降，会导致含水层和弱透水层发生层内渗流及层间越流，伴随这一渗流过程中的孔压和有效应力转化会导致土体发生变形，最终诱发地面沉降[13]。骆冠勇等[13]给出了含水层降水引起的弱透水层1维弹性沉降固结度计算公式。Tseng等[14]考虑了土体重力的影响，给出了含水层中水位下降为定值时弱透水层的1维弹性固结解。Tao等[15]给出含水层水位瞬时下降和线性下降两种模式下的单层土体1维固结弹性解。吴浩等[16]考虑了弱透水层黏性土的结构性，推导了含水层降水所引起的弱透水层1维固结解。最近，Lo等[17]推导了流体通量边界条件下非饱和土的1维弹性解。可以看出，为了反映地面沉降的真实演化特征，现有解析研究不断深入，趋向于考虑水位变化下更复杂的土体变形特性。

地面沉降具有持续时间长的特点。大量现场观测和试验研究表明，控制地下水开采以后，地下水位得到稳定甚至回升，地层变形仍会长期发展。一般认为，产生这种变形滞后现象的重要原因之一在于弱透水层黏性土具有的流变性[18–20]。Liu等[21]考虑弱透水层土体的黏滞性，并给出了水位变化下的1维固结半解析解。Li等[22]在此基础上考虑了弱透水层中流体为非达西流，推导了水位变化下的单层地基土的1维固结黏弹性解。同时，近年来进一步的研究发现，由于颗粒间的滑移、错动以及破碎等细观机理，含水层砂土同样表现出一定的流变性，当层厚较大时，含水层的这种蠕变变形将不可忽略[23–26]。因此，给出水位下降诱发含水层–弱透水层双层系统的1维固结模型，在其中考虑含水层和弱透水层的流变性，对于地面沉降的合理评估具有一定理论意义。

为此，针对地下水位变化引起的土层固结问题，用广义开尔文模型描述饱和弱透水层黏土与饱和含水层砂土的黏滞性，建立固结变形控制方程，给出求解条件，经过Laplace变换等数学推导给出传统矩阵传递法和边界转换法两种计算方法，求得水位下降诱发弱透水层–含水层1维固结的半解析解。采用Stehfest数值方法实现Laplace逆变换计算，验证本文计算方法的正确性。通过算例分析传统矩阵传递法和边界转换法的区别，进一步探讨土体黏滞性、渗透性等力学性质与水文地质性质差异大小以及水位下降速率对土层变形的影响。

1 数学模型

水位变化引起的弱透水层–含水层1维固结模型如图1所示。

图1 水位变化下弱透水–含水层固结示意图Fig. 1 Consolidation schematic diagram of aquitard–aquifer due to groundwater level variation

假定弱透水层顶部与含水层底部总水头相等且为h，弱透水层顶部有足够补给，其水头始终保持不变，在含水层底面t时刻的水位变化为 ∆h(t)。 用参数n表示土层编号（n=1,2 ），每层土厚度依次是H1和H2，且满足总厚度为H=H1+H2；kv1、kv2分别为两土层的渗透系数。

1.1 控制方程

根据达西定律和有效应力原理，可得两层土体的控制方程为[21]：

式中：kvn为土体渗透系数，其中n=1,2 ；hn为 第n层土的超孔隙水头高度；εzn为 竖向应变；z和t分别为土层离地面的距离和时间。

土体的黏滞性用广义开尔文模型描述，结构如图2所示，由一系列虎克弹簧和牛顿黏壶组成。该模型的优点为可以很好地描述流变性相对较弱的含水层砂土的蠕变特性，同时也可以退化成常用的Merchant模型（结构简图见图3），以描述弱透水层黏性土的流变特性。

图2 广义开尔文流变模型Fig. 2 Generalized Kelvin rheological model

图3 Merchant流变模型Fig. 3 Merchant rheological model

基于广义开尔文流变模型，竖向应变εzn可以表示成如下积分形式：

式中，γw为单位体积水的重度， δn(0)为流变模型的柔度系数 δn(t) 在t=0 时 刻的值。δn(t)表达式为：

式中， η1n=E1n/K1n， η2n=E2n/K2n，其中E0n、E1n、E2n、K1n、K2n为广义开尔文体流变模型元件的参数（图2）。

将式（3）代入式（2），再将式（2）代入式（1），则可以得到基于广义开尔文流变模型的固结控制方程为：

1.2 定解条件

在两层土体交界面（即z=H1）处，应满足水头和流量连续条件：

上、下边界条件满足：

初始条件为：

2 问题求解

2.1 矩阵传递法

矩阵传递法是利用Laplace变换，对待求微分方程与定解条件实现转换后，再根据定解条件求解方程得到变量的矩阵表达形式，最后通过矩阵运算法及Laplace逆变换求得解的一种经典解析方法[27]。详细求解过程如下：

首先，进行如下变量代换：

将式（10）代入控制方程（4）及求解条件（5）～（9）得：

根据初始条件（16），对控制方程（11）及边界条件（12）～（15）关于时间t进行Laplace变换得：

对于控制方程（17）进行求解得：

将边界条件（18）～（21）代入通解式（22）并整理成矩阵形式得：

式中：

从而可得超静孔压：

式中，un(z,t) 为t时刻土层各位置的超孔隙水压。

2.2 边界转换法

边界转换法是由Chen等[28]提出的一种新的边界处理方法，其最大的优点在于可将复合边界转化为统一的单一边界条件。本文的渗流连续边界条件（5）及（6）就属于Cauchy边界（即边界上同时作用第1类及第2类边界）[29]。采用边界转换法将渗流连续边界转换为单一边界进行求解。

根据边界转换法，将渗流连续边界转换为第1类边界，并结合土层顶面及底面边界可综合表示为：

基于边界式（26），控制方程式（17）的通解为：

将通解式（27）代入土层顶面边界（20）、流量连续条件（19）及底面边界（21），并整理成矩阵形式：

式中：

2.3 沉降和固结度计算公式

根据式（2）可知，任意时刻双层土的总沉降为：

由式（29）可知，任意时刻下按沉降定义的双层土体平均固结度可表示为：

3 两种解法的验证与对比

为了得到定解问题在真实物理空间的解，需要进行Laplace逆变换，本文选用Stehfest法[30]进行数值逆变换计算。利用MATLAB分别编写了计算程序，为验证本文计算方法和程序的正确性，进行了算例对比验证。

3.1 单层地基下本文解法的验证

Tao等[15]给出了水头瞬时下降时弱透水层变形的弹性解，如下：

式中，S(t)为t时 刻下土体的沉降量，Tv为土体的固结时间因数，h0与hd分别为初始水头高度及抽水后含水层的水头高度，H为弱透水层厚度，N=(2n−1)（n=1,2,3,···），。

将本文解退化到这种单层情形，利用计算程序得出水头瞬时下降时弱透水层变形的弹性解，与文献[15]进行对比，见图4。算例中参数取值如下：H1=H2=5m，kv1=kv2=8.64×10−4m/d，E01=E02=2MPa，E11=E12=5MPa， η11=η12近 似取无穷大（如1 ×1010d−1），初始水头h=12m ，瞬时水位变化 ∆h(t)=10m，hd=2m。由图4可以看出：不考虑黏滞性时，本文两种计算方法所得的解与文献[15]所得的解析解结果均吻合良好，在变形初期（ 25d之前）边界转换法所得结果略优于矩阵传递法。这是由于在数值逆变换中，当时间t较小时会出现计算结果溢出，导致计算精度降低甚至出错（如当t取 到1 ×10−6d），边界转换法在求解中的未知参数少于矩阵传递法，所以在 25d之前边界转换法计算结果优于矩阵传递法。

图4 单层地基土本文解与文献[15]所得弹性解的沉降对比Fig. 4 Comparison of Settlement for single-layer foundation obtained from solutions by the proposed methods and reference [15] method

此外，Liu等[21]基于Merchant流变模型给出了水位变换引起单层土体1维固结的半解析解。为了便于对比，将本文广义开尔文模型退化为Merchant模型，双层解退化成单层解，利用编写的计算程序，分别计算含水层水头瞬时下降，在t=100d 和t=500d时弱透水层变形的黏弹性解，将本文解与Liu等[21]的解对比，如图5所示。算例中 η11=η12= 8 .64×10−4d−1，其余参数保持不变。可以看出，考虑黏滞性时两种情形下本文计算结果都与Liu等[21]的解吻合良好，而且矩阵传递法和边界转换法都表现出良好的一致性。

图5 单层地基土本文解与文献[21]所得解的孔压对比Fig. 5 Comparison of pore pressures for single-layer foundation obtained from solutions by the proposed methods and reference [21] method

3.2 双层地基下本文解法的验证

上述两个退化算例验证了两种特殊情况，原则上，理论计算模型的合理性和适用性仍需要通过现场实测资料或者室内试验结果加以验证。现场的水文地质条件和边界条件复杂，无法完全确定实测结果的控制因素。相对而言，室内模型试验可以对试验条件进行人为控制和简化，去除了很多现场实测中的不确定因素，所得的结果更有助于验证理论解的有效性[31]。

为此，利用本文边界转换法计算程序对村山朔郎的大比例地面沉降模型试验进行了分析[32]，试验模型如图6所示。

图6 地面沉降模型试验示意图Fig. 6 Schematic diagram of model experiments on land subsidence

在试验中，上覆层弱透水层始终保持自由水面，并且水头在含水层顶板以上1.5 m处；下卧含水砂层水头从初始水头h0=1.5m 下 降 ∆h=−1m来模拟水位降深，并在整个试验过程中长期观测地层的沉降量。

图7给出了弹性和黏弹性两种情形下的本文理论解与试验值的对比结果。计算参数如下：上覆弱透水层的渗透系数kv1=0.00006m/d ，弹性模量Es1=0.48MPa ；下层砂土的渗透系数kv2=1.12m/d，弹性模量Es2=20MPa；黏弹性解中，弱透水层的黏滞系数η11=0.01d−1，E0=E1=0.48MPa，含水层砂土参数不变。图7可以看出：前期（含水层抽水开始 40d以内）试验值与弹性解、黏弹性计算值均比较吻合，但是，随着时间增长，弹性解越来越偏离试验值；相比而言，黏弹性解在整体上与 2 00d的试验结果吻合更加良好。图7中同时给出了该双层地基问题的黏弹性半解析数值模拟结果[4]，由图7可以看出本文计算解与该数值结果也保持了较好的一致性，体现了水位下降引起的含水层与弱透水层长期变形特性。

图7 弹性和黏弹性情形下本文解与半解析数值解以及试验值的验证Fig. 7 Verification of the present solutions with the semianalytical numerical solutions and experimental values for elastic and viscoelastic cases

3.3 本文两种计算方法的对比

矩阵传递法是从固体力学规则厚板推广而成的经典解析方法，可操作性强，只要列出矩阵方程，便可通过矩阵运算法则进行求解。但是，随着土层划分数N的增加，求解方程（20）时则有2N个未知数，这会使得计算规模翻倍，计算效率降低。边界转换法通过转换边界，使得层间接触面水头的连续性在解方程之前得到满足，对于N层土而言，只需N+1个参数。因此，相比于矩阵传递法，边界转换法在计算工作量方面具有一定优势，而且在处理混合边界条件时也更加方便。为进一步分析两种方法的计算效率，通过数值算例，分别计算分析两种解法下沉降量的长期发展过程，并将计算结果绘制于图8。

图8 两种解法下的沉降值Fig. 8 Settlement values under two solutions

图8结果表明，两种解法得到的沉降值计算精度相近，但是计算效率相差较大，其计算效率矩阵传递法比边界转换法高74.56%。随着土层数增加，边界转换法的计算优势将更加明显。此外，从图4中看出，在弹性解中，变形初期（ 25d之前）边界转换法的精度略高于矩阵传递法。

4 影响土体沉降因素分析

在上述边界转换法计算程序基础上，通过数值算例进一步探讨了土体的黏滞性、渗透性、层状土性质的差异性以及水头下降速率对水位变化引起土层变形过程的影响。

算例具体描述如下：某双层土地基，土层总厚H=12m，H1= 9m，H2=3m ，初始总水头高度h=14m，t时刻承压层瞬时降水高度 ∆h(t)=12m，土体竖向渗透系数kv1=9×10−4m/d ，kv2= 3 ×10−4m/d。两层土均采用Merchant模型：第1层土模型参数为E01=2000kPa，E11= 5 000kPa ，η11=10−4d−1；第2层土模型参数为E02=4800kPa，E12=4800kPa ，η12= 1 0−3d−1。

4.1 土体黏滞性的影响

黏滞系数η1n是 描述土体流变性的主要参数，其大小反映了土体黏滞性的程度。为了充分探讨黏滞性对土体长期变形的影响，参考文献[21]给出了关于黏滞系数η1n一 般取值范围，本文所采用的黏滞系数η1n取值见表1。算例中其余参数不变。当η1n接近1010倍时，近似认为η1n趋 于∞，此时流变模型等效于弹性模量为E0E1/(E0+E1)的串联弹性模型。用计算程序进行了数值算例，将计算得到的固结度–时间关系绘于图9。

表1 黏滞系数的取值Tab. 1 Values of viscosity coefficients

图9 黏滞系数对固结度的影响Fig. 9 Influence of viscosity cofficients on consolidation process

由图9可以看出：相同时间下，黏滞系数小的固结度低，并且，随着η1的减小，土体的蠕变性越明显；相同水位变化下土层达到最终沉降所需的时间越长，即完成固结所需的时间越长，特别是在后期（t=100d之后）这种黏滞性的影响效应越明显。

4.2 渗透性的影响

渗透系数是描述孔隙水在土体内部流动能力的主要参数，为了探讨其对土体变形的影响，对两层土的渗透系数kvn同时逐级扩大，取值见表2，其余参数保持不变。用计算程序进行了数值算例，将计算得到的渗透系数影响下，固结–时间关系曲线绘于图10。从图10中可以看出：在固结前期（1 00d之前），影响固结的主要因素为渗透系数，渗透系数越大，固结完成的越快；但在后期该影响趋弱，影响的主要因素是土体的黏滞性。

表2 渗透系数取值Tab. 2 Values of permeability coefficients

图10 渗透系数对固结度的影响Fig. 10 Influence of permeability coefficients on consolidation process

4.3 土层性质差异的影响

土体为非均质各向异性体。为了研究土层间力学性质差异的影响，进行了数值算例分析，参数取值为：土层厚度均取6 m，初始总水头高度为h=14m，t时刻承压层瞬时降水高度 ∆h(t)=12m，土体竖向渗透系数kv1=kv2=6×10−4m/d，E01=E02= 2 000kPa，E11=E12=5000kPa，土层的黏滞系数取值如表3所示。

表3 层状土黏滞性的差异Tab. 3 Difference of layered soils in viscosity property

为了探讨土层间水文地质性质差异对固结的影响，在保持其他参数不变的前提下，土层的黏滞系数取η11=η12=6×10−4d−1，土层的竖向渗透系数的取值见表4。土层黏滞系数和渗透性差异对固结度的影响如图11～12所示。

表4 层状土渗透性的差异Tab. 4 Difference of layered soils in permeability

图11 土体黏滞性差异对固结度的影响Fig. 11 Influence of soil viscosity difference on consolidation process

从图11可以看出，土层黏滞性的差异对固结前期（ 400d之前）影响较小，只有在后期这种差异性才逐渐表现出来。区别于黏滞性差异，图12表明：土层的渗透性差异对于固结的影响主要表现在前期（400 d之前）；渗透性差异越大，前期固结度越高，并且在100d之后固结速率开始逐渐变慢；在固结后期，渗透性差异的影响趋弱。

图12 土体渗透性差异对固结度的影响Fig. 12 Influence of soil permeability difference on consolidation process

4.4 水位下降速率的影响

以往研究大多假设抽水作用下含水层的水位下降是瞬时完成的，而实际情形中，地下水水位的降深不可能瞬时发生，应该存在一个渐变的过程[13,16]。Liu等[21]以指数函数∆h(t)=H(1−e−bt)描述抽水作用下水位变化过程，其中，参数b综合反映了井的抽水效率和开采含水层的水文地质条件，且b＞0，单位为d−1。当参数b足 够大（如b＞10000）时，相当于水头下降瞬时完成这种理想情形。

图13给出了参数b取不同值时的固结度变化曲线，其余参数取值同上述算例描述。由图13可以看出：当参数b越小，即抽水水位下降速率越慢时，固结度越小；前期沉降发生时间与参数b有 关，即参数b越小，地面沉降越迟发生，最终沉降量不受水位下降速率的影响。

图13 参数b对固结度的影响Fig. 13 Influence of parameter b on consolidation process

5 结 论

1）利用广义开尔文流变模型和1维固结理论建立了水位变化下弱透水层–含水层越流系统土体变形的控制方程。通过数学推导给出了传统矩阵传递法和边界转换法两种解法，得到了该固结方程的半解析解。编制了计算程序，通过与已有的单层解、双层室内试验数据、数值解进行对比，验证了本文两种解法的正确性。

2）对比两种解法发现，固结初期，同种数值逆变下边界转换法计算精度略高于传统矩阵传递法。边界转换法可以将复杂的混合边界处理为统一的单一边界条件，更适合于双层地基土固结模型中层间连续性条件的处理，当土层数较多时，边界转换法更加适宜。

3）基于边界转换法编写的计算程序进行了数值算例。结果表明：土体的蠕变性越明显，完成固结所需时间越长，土层变形越滞后，而且土层力学性质差异对后期固结影响较大；相反地，水文地质参数及其差异主要影响中前期固结过程，土层渗透性越强，土体固结度越高。值得注意的是，抽水作用下水位下降速率对地面沉降过程也有较大影响，可能会导致土层变形发生的延迟。