条形荷载下饱和软土地基二维非线性固结分析

Two-dimensional nonlinear consolidation analysis of saturated soft clay foundation under strip load

摘要

研究了均质黏性土地基在条形荷载作用下的二维非线性固结变形规律. 基于Terzaghi固结理论，考虑土体渗透性和压缩性的非线性关系，建立了土体二维非线性固结控制方程. 利用有限差分法求解，并通过解析解和数值解验证求解的合理性，分析了压缩指数与渗透指数比值C_c/C_k、初始荷载q₀、渗透系数及其各向异性和地基深度H对固结行为的影响，探讨了不同位置的沉降变化. 结果表明，各因素均显著影响固结特性. C_c/C_k越大，非线性固结速率越慢，沉降速度越慢. 初始荷载q_{\mathrm{0}}对非线性固结速率影响与C_c/C_k有关. 渗透系数越大，固结速率越快，沉降速度也越快，竖向渗透系数影响比水平渗透系数大. 地基深度H越大，非线性固结速率越慢，沉降量越大. 地基沉降分布呈现中轴线最大、向两侧递减的规律.

Abstract

This paper investigates the two-dimensional nonlinear consolidation deformation law of a homogeneous clay foundation under strip loading. Based on Terzaghi’s consolidation theory and considering the nonlinear relationships between soil permeability and compressibility， the governing equations for two-dimensional nonlinear soil consolidation are established. The finite difference method is employed to solve the equations， and the rationality of the solution is verified through analytical and numerical solutions. The effects of the compression index to permeability index ratio C_c/C_k， initial load q₀， permeability coefficient and its anisotropy， as well as foundation depth H on consolidation behavior are analyzed， and the settlement variations at different locations are examined. The results indicate that all these factors significantly influence the consolidation characteristics. A higher C_c/C_k ratio leads to slower nonlinear consolidation and settlement rates. The impact of initial load on the nonlinear consolidation rate is related to C_c/C_k. A higher permeability coefficient results in faster consolidation and settlement rates， and the vertical permeability coefficient has a greater influence than the horizontal one. A greater foundation depth H corresponds to slower nonlinear consolidation rates and larger settlement amounts. The settlement distribution of the foundation follows a pattern where it is maximum at the central axis and it decreases toward both sides.

Graphical abstract

关键词

固结 / 条形荷载 / 二维非线性 / 平均固结度 / 地表沉降

Key words

consolidation / strip load / two-dimensional nonlinearity / average degree of consolidation / surface settlement

固结现象在土体工程领域极为普遍，例如路基、坝体、建筑物基础等工程都会涉及地基固结问题，地基固结会对工程的安全性和稳定性产生重要影响，因此，研究地基固结变形规律具有重要的工程意义. 传统的Terzaghi固结理论^［1］假设土体的渗透性和压缩性为线性关系，然而，实际上土体的性质通常表现出非线性特征，使得传统理论在某些情况下难以准确描述固结过程. 近年来，非线性固结理论得到了广泛关注，并逐渐应用于工程实践中.

众多学者在非线性达西渗流问题方面进行了深入研究，提出了多种近似解析解和数值方法. Kim等^［2］提出了饱和多层土体在时变荷载下一维非线性固结的解析解. Chen等^［3］采用差分求积法来分析具有部分排水边界和任意加载条件的多层土体的一维非线性固结问题. Oliveira等^［4］提出了一种考虑蠕变的一维饱和黏土固结模型，采用非线性黏性定律，能够准确模拟试验中的时间次生效应和长期沉降. 李传勋等^［5］提出了一种压缩指数与渗透指数比值不等于1的软土一维非线性固结近似解析解. 江留慧等^［6］推导了变荷载下考虑土体固结系数变化的双层地基一维非线性固结近似解析解. 夏长青等^［7］研究了连续排水边界下成层地基的一维非线性固结问题，推导了平均固结度和沉降的解析解. 以上研究仅考虑土体竖直方向的固结，忽略了水平方向的固结效应，而许多工程问题中土体的固结是具有二维甚至三维特性的.

因此目前众多学者对二维固结也进行了许多研究. 张添等^［8］研究了分段循环荷载作用下二维非饱和土的固结特性，分析了渗透系数和荷载特征参数对沉降的影响. 胡安峰等^［9］研究了隧道周围饱和软土的非线性固结问题，建立了考虑土体压缩性和渗透性非线性变化的二维固结控制方程，并利用交替隐式差分法进行求解. 刘奔奔^［10］研究了堆载作用下黏土砂层交替分布地基的二维固结理论，建立了考虑待固结土层相互影响的单元和整体固结计算模型. 黄传志等^［11］基于Biot固结理论，针对二维固结问题，给出了有限厚地基和半无限地基在不同荷载作用下的孔隙水压力、应力和位移的解析解. 艾智勇等^［12］研究了竖向简谐荷载作用下二维层状饱和地基的动力响应，并分析了荷载圆频率、荷载作用深度和地基成层性对地基竖向位移的影响. 关聪^［13］研究了考虑Hansbo渗流的二维Biot固结理论，并通过数值分析探讨了非Darcy 渗流参数和地基厚度对土体固结和Mandel效应的影响. Wang等^［14-16］通过傅里叶正弦级数和拉普拉斯变换，提出了针对不同初始和边界条件的非饱和土体二维平面应变固结的半解析解. Shen等^［17］提出了针对非饱和土体在条形分布渗透边界条件下的二维平面应变固结半解析解. Sun等^［18］提出了一种基于整体分析的二维模型来研究层状地基中各黏土层的相互作用对土壤固结的影响. Wang等^［19］通过半解析方法研究了随时间和深度变化的应力扩散对非饱和土体二维平面应变固结特性的影响. 这些研究涵盖了不同类型的土体和荷载条件下的固结行为，更接近实际工程，不过在非线性方面的考虑不够充分.

在此基础上，本文基于Terzaghi-Rendulic二维固结理论、Butterfield^［20］和谢康和等^［21］提出的双对数压缩、渗透模型，构建了非线性固结控制方程. 通过这一研究，深入探讨地基固结与沉降问题，并分析土体的压缩指数、渗透指数、初始荷载、渗透系数及其各向异性和地基深度等因素对土体非线性固结特性的影响规律，以期望为更精确地预测地基变形和稳定状态提供理论支持.

1 固结控制方程

1.1 固结控制方程推导

土体的渗透性和压缩性需要分别满足e-\mathrm{l}\mathrm{g}\sigma \mathrm{\text{'}}、 e-\mathrm{l}\mathrm{g}k的非线性关系，即：

e=e_{\mathrm{0}}-C_{\mathrm{c}}\mathrm{l}\mathrm{g}\left(\sigma \mathrm{\text{'}}/{\sigma }_{\mathrm{0}}\mathrm{\text{'}}\right)（1）

e=e_{\mathrm{0}}-C_{\mathrm{k}}\mathrm{l}\mathrm{g}\left(\frac{k_{\mathrm{v}}}{k_{\mathrm{v}\mathrm{0}}}\right)（2）

式中：C_c为压缩指数；C_k为渗透指数；e_{\mathrm{0}}为初始有效应力{\sigma }_{\mathrm{0}}^{\text{'}}对应的初始孔隙比；\sigma \mathrm{\text{'}}为有效应力；k_v为孔隙比e对应的渗透系数；k_v0为初始孔隙比e_{\mathrm{0}}对应的渗透系数.

由式（1）可得：

m_{\mathrm{v}}=\frac{\mathrm{1}}{\mathrm{1}+e_{\mathrm{0}}}\frac{\partial e}{\partial \sigma \mathrm{\text{'}}}=m_{\mathrm{v}\mathrm{0}}\frac{{\sigma }_{\mathrm{0}}\mathrm{\text{'}}}{\sigma \mathrm{\text{'}}}（3）

式中：m_{\mathrm{v}}为体积压缩系数；m_{\mathrm{v}\mathrm{0}}为初始体积压缩系数.

m_{\mathrm{v}\mathrm{0}}=\frac{C_{\mathrm{c}}}{(\mathrm{1}+e_{\mathrm{0}}){\sigma }_{\mathrm{0}}^{\mathrm{\text{'}}}}\frac{\mathrm{1}}{\mathrm{l}\mathrm{n}\mathrm{10}}（4）

式（1）、式（2）联立得任意时刻的土体渗透系数.

k_{\mathrm{v}}=k_{\mathrm{v}\mathrm{0}}{\left(\frac{{\sigma }_{\mathrm{0}}^{\mathrm{\text{'}}}}{\sigma \mathrm{\text{'}}}\right)}^{\frac{C_{\mathrm{c}}}{C_{\mathrm{k}}}}（5）

于是可得土体固结系数为：

C_{\mathrm{v}}=\frac{k_{\mathrm{v}}}{m_{\mathrm{v}}{\gamma }_{\mathrm{w}}}=C_{\mathrm{v}\mathrm{0}}{\left(\frac{{\sigma }_{\mathrm{0}}^{\mathrm{\text{'}}}}{\sigma \mathrm{\text{'}}}\right)}^{\frac{C_{\mathrm{c}}}{C_{\mathrm{k}}}-\mathrm{1}}（6）

式中：{\gamma }_{\mathrm{w}}为孔隙水重度；C_{\mathrm{v}\mathrm{0}}为土体初始固结系数.

C_{\mathrm{v}\mathrm{0}}=\frac{k_{\mathrm{v}\mathrm{0}}}{m_{\mathrm{v}\mathrm{0}}{\gamma }_{\mathrm{w}}}=\frac{\left(\mathrm{1}+e_{\mathrm{0}}\right)k_{\mathrm{v}\mathrm{0}}{\sigma }_{\mathrm{0}}^{\mathrm{\text{'}}}\mathrm{l}\mathrm{n}\mathrm{10}}{C_{\mathrm{c}}{\gamma }_{\mathrm{w}}}

由式（6）可知，当C_{\mathrm{c}}/C_{\mathrm{k}}\ne \mathrm{1}时，土体固结系数会随着有效应力不断变化；当C_{\mathrm{c}}/C_{\mathrm{k}}=\mathrm{1}时，土体固结系数为常数，该比值反映了固结过程中的非线性特征.

由Terzaghi-Rendulic二维固结理论，对于均质饱和软土地基，土体固结控制方程为：

\frac{\mathrm{1}}{{\gamma }_{\mathrm{w}}}\frac{\partial }{\partial x}\left(k_{\mathrm{v}}\frac{\partial u}{\partial x}\right)+\frac{\mathrm{1}}{{\gamma }_{\mathrm{w}}}\frac{\partial }{\partial y}\left(k_{\mathrm{v}}\frac{\partial u}{\partial y}\right)=-\frac{\partial {\epsilon }_{\mathrm{v}}}{\partial t}（7）

式中：k_{\mathrm{v}}为土体渗透系数；{\epsilon }_{\mathrm{v}}为土体体积应变；u为超孔隙水压力.

土体体积应变{\epsilon }_{\mathrm{v}}为：

{\epsilon }_{\mathrm{v}}=\frac{e_{\mathrm{0}}-e}{\mathrm{1}+e_{\mathrm{0}}}（8）

则土体体积应变变化率为：

\frac{\partial {\epsilon }_{\mathrm{v}}}{\partial t}=-\frac{\mathrm{1}}{\mathrm{1}+e_{\mathrm{0}}}\frac{\partial e}{\partial t}（9）

根据有效应力原理：

\sigma \mathrm{\text{'}}=\sigma -u（10）

则：

\frac{\partial \sigma \mathrm{\text{'}}}{\partial t}=\frac{\partial u}{\partial t}（11）

联立式（3）、式（9）、式（11），可得：

\frac{\partial {\epsilon }_{\mathrm{v}}}{\partial t}=-m_{\mathrm{v}}\frac{\partial u}{\partial t}（12）

将式（5）、式（12）代入式（7），即得到本文研究的土体二维非线性固结控制方程：

\frac{\partial }{\partial x}\left(k_{\mathrm{v}\mathrm{0}}{\left(\frac{{\sigma }_{\mathrm{0}}\mathrm{\text{'}}}{\sigma \mathrm{\text{'}}}\right)}^{\frac{C_{\mathrm{c}}}{C_{\mathrm{k}}}}\frac{\partial u}{\partial x}\right)+\frac{\partial }{\partial \mathrm{y}}\left(k_{\mathrm{v}\mathrm{0}}{\left(\frac{{\sigma }_{\mathrm{0}}\mathrm{\text{'}}}{\sigma \mathrm{\text{'}}}\right)}^{\frac{C_{\mathrm{c}}}{C_{\mathrm{k}}}}\frac{\partial u}{\partial y}\right)=

\frac{C_{\mathrm{c}}{\gamma }_{\mathrm{w}}}{\left(\mathrm{1}+e_{\mathrm{0}}\right){\sigma }_{\mathrm{0}}\mathrm{\text{'}}\mathrm{l}\mathrm{n}\mathrm{10}}\frac{\partial u}{\partial t}（13）

1.2 平均固结度的计算

在实际工程中常采用平均固结度来衡量土体中孔隙水压的消散程度，当任意时刻t的孔隙水压被计算后，进一步可得平均固结度为：

U_{\mathrm{t}}=\mathrm{1}-\frac{{\iint }_D^{}u\left(x,y,t\right)\mathrm{d}A}{{\iint }_D^{}u\left(x,y,\mathrm{0}\right)\mathrm{d}A}（14）

式中：A为计算区域D的面积.

1.3 沉降计算

土体内部任意一点土体体积应变为：

{\epsilon }_{\mathrm{v}}={\epsilon }_x+{\epsilon }_y+{\epsilon }_z=\frac{\mathrm{1}-\mathrm{2}\nu }{E}\left({\sigma }_x^{\text{'}}+{\sigma }_y^{\text{'}}+{\sigma }_z^{\text{'}}\right)（15）

式中：{\epsilon }_x\mathrm{、}{\epsilon }_y\mathrm{、}{\epsilon }_z分别为土体在x、y、z方向的应变；{\sigma }_x^{\text{'}}\mathrm{、}{\sigma }_y^{\text{'}}\mathrm{、}{\sigma }_z^{\text{'}}分别为土体在x、y、z方向的有效应力；\nu为土体泊松比，取值0.36；E为土体弹性模量.

根据广义胡克定律及平面应变假设，可得：

{\sigma }_z^{\text{'}}=\nu \left({\sigma }_x^{\text{'}}+{\sigma }_y^{\text{'}}\right)（16）

\left\{\begin{array}{c}{\epsilon }_x=\frac{\mathrm{1}-{\nu }^{\mathrm{2}}}{E}\left({\sigma }_x^{\text{'}}-\frac{\nu }{\mathrm{1}-\nu }{\sigma }_y^{\text{'}}\right)\\ {\epsilon }_y=\frac{\mathrm{1}-{\nu }^{\mathrm{2}}}{E}\left({\sigma }_y^{\text{'}}-\frac{\nu }{\mathrm{1}-\nu }{\sigma }_x^{\text{'}}\right)\end{array}\right.（17）

{\sigma }_x\mathrm{\text{'}}=K_{\mathrm{0}}{\sigma }_y\mathrm{\text{'}}（18）

式中：K_{\mathrm{0}}为土静止土压力系数，规范中正常固结饱和软土的经验值为0.5~0.6^［22］，本文取值0.5.

联立式（15）~式（18）可得：

{\epsilon }_y=\frac{-\nu K_{\mathrm{0}}+\left(\mathrm{1}-\nu \right)}{\left(\mathrm{1}-\mathrm{2}\nu \right)\left(K_{\mathrm{0}}+\mathrm{1}\right)}{\epsilon }_{\mathrm{v}}（19）

同时将式（1）代入式（8），得：

{\epsilon }_{\mathrm{v}}=\frac{C_{\mathrm{c}}}{\left(\mathrm{1}+e_{\mathrm{0}}\right)}\mathrm{l}\mathrm{g}\left(\frac{\sigma \mathrm{\text{'}}}{{\sigma }_{\mathrm{0}}^{\text{'}}}\right)（20）

由此，结合式（20），对式（19）沿深度方向进行积分可得任意时刻轴线位置处的地表固结沉降值：

S_{\mathrm{t}}={\int }_{-H}^{\mathrm{0}}{\epsilon }_y\left(x,y,t\right)\mathrm{d}y（21）

式中：H为地基深度.

1.4 初始条件

采用图1所示的饱和软黏土地基作为计算模型. 地基深度为H=\mathrm{10}\mathrm{m}，地基宽度为\mathrm{2}L=\mathrm{70}\mathrm{m}，本文研究的条形荷载为瞬时均布荷载，覆盖宽度为\mathrm{2}B=\mathrm{20}\mathrm{m}，大小为q_{\mathrm{0}}，模拟实际工程中快速填土或大面积堆载的情况. 假定初始有效应力恒定时，其值为土层总厚度一半处的土体自重应力大小^［23］，即{\sigma }_{\mathrm{0}}^{\text{'}}=\gamma \cdot H/\mathrm{2}=\mathrm{90}\mathrm{k}\mathrm{P}\mathrm{a}. 以对称轴为y轴，水平向右为x轴. 假定地基在自重作用下已完成固结，之后的固结仅在外加荷载作用下发生，以瞬时荷载施加后地基中产生的超静孔隙水压力为初始值，则依据弗洛林理论^［24］，二维土层内超静孔隙水压力u_{\mathrm{0}}的初始分布为：

u_{\mathrm{0}}\left(x,y,\mathrm{0}\right)={\sigma }_{\mathrm{0}}\left(x,y\right)=\frac{\mathrm{1}}{\mathrm{2}}\left({\sigma }_x+{\sigma }_y\right)=

\frac{q_{\mathrm{0}}}{\mathrm{\pi }}\left[\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{t}\mathrm{a}\mathrm{n}\left(\frac{B+x}{y}\right)+\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{t}\mathrm{a}\mathrm{n}\left(\frac{B-x}{y}\right)\right]（22）

式中：σ₀为初始应力；{\sigma }_x和{\sigma }_y分别为x和y方向的应力；q_{\mathrm{0}}为地基表面施加的初始荷载；B为施加瞬时荷载的半宽度.

1.5 边界条件

1） 地基顶面透水\left(y=\mathrm{0}\right)

u\left(x,\mathrm{0},t\right)=\mathrm{0}（23）

2） 地基底面不透水\left(y=-H\right)

\frac{\partial u\left(x,-H,t\right)}{\partial y}=\mathrm{0}（24）

3） 对称轴不透水\left(x=\mathrm{0}\right)

\frac{\partial u\left(\mathrm{0},y,t\right)}{\partial x}=\mathrm{0}（25）

4） 水平向远端透水\left(x=L\right)

u\left(L,y,t\right)=\mathrm{0}（26）

1.6 有限差分法求解

对式（13）进行差分法求解，将计算区域在x和y方向上分别进行离散化，划分成等间距的网格. 设x方向的网格间距为\mathrm{\Delta }x，y方向的网格间距为\mathrm{\Delta }y，时间步长为\mathrm{\Delta }t. 对于节点\left(i,j\right)，其坐标为\left(x_i=i\mathrm{\Delta }x,y_i=j\mathrm{\Delta }y\right)，对应的超静孔隙水压力为u_{i,j}^n，其中n表示时间步数.

对于\frac{\partial u}{\partial x}，采用中心差分格式，在节点\left(i,j\right)处有：

{\frac{\partial u}{\partial x}}_{\left|i,j\right.}\approx \frac{u_{i+\mathrm{1},j}^n-u_{i-\mathrm{1},j}^n}{\mathrm{2}\mathrm{\Delta }x}（27）

对于\frac{\partial }{\partial x}\left(k_{\mathrm{v}\mathrm{0}}{\left(\frac{{\sigma }_{\mathrm{0}}^{\text{'}}}{\sigma \mathrm{\text{'}}}\right)}^{\frac{C_{\mathrm{c}}}{C_{\mathrm{k}}}}\frac{\partial u}{\partial x}\right)，先将k_{\mathrm{v}\mathrm{0}}{\left(\frac{{\sigma }_{\mathrm{0}}^{\text{'}}}{\sigma \mathrm{\text{'}}}\right)}^{\frac{C_{\mathrm{c}}}{C_{\mathrm{k}}}}在节点\left(i,j\right)处的值记为k_{ij}，则有：

\frac{\partial }{\partial x}{\left(k_{\mathrm{v}\mathrm{0}}{\left(\frac{{\sigma }_{\mathrm{0}}^{\text{'}}}{\sigma \mathrm{\text{'}}}\right)}^{\frac{C_{\mathrm{c}}}{C_{\mathrm{k}}}}\frac{\partial u}{\partial x}\right)}_{\left|i,j\right.}\approx

\frac{k_{i+\mathrm{1},j}\frac{u_{i+\mathrm{2},j}^n-u_{i,j}^n}{\mathrm{2}\mathrm{\Delta }x}-k_{i-\mathrm{1},j}\frac{u_{i,j}^n-u_{i-\mathrm{2},j}^n}{\mathrm{2}\mathrm{\Delta }x}}{{\left(\mathrm{\Delta }x\right)}^{\mathrm{2}}}（28）

对于\frac{\partial u}{\partial y}，采用中心差分格式，在节点\left(i,j\right)处有：

{\frac{\partial u}{\partial y}}_{\left|i,j\right.}\approx \frac{u_{i,j+\mathrm{1}}^n-u_{i,j-\mathrm{1}}^n}{\mathrm{2}\mathrm{\Delta }y}（29）

对于\frac{\partial }{\partial y}\left(k_{\mathrm{v}\mathrm{0}}{\left(\frac{{\sigma }_{\mathrm{0}}^{\text{'}}}{\sigma \mathrm{\text{'}}}\right)}^{\frac{C_{\mathrm{c}}}{C_{\mathrm{k}}}}\frac{\partial u}{\partial y}\right)，先将k_{\mathrm{v}\mathrm{0}}{\left(\frac{{\sigma }_{\mathrm{0}}^{\text{'}}}{\sigma \mathrm{\text{'}}}\right)}^{\frac{C_{\mathrm{c}}}{C_{\mathrm{k}}}}在节点\left(i,j\right)处的值记为k_{ij}，则有：

\frac{\partial }{\partial y}{\left(k_{\mathrm{v}\mathrm{0}}{\left(\frac{{\sigma }_{\mathrm{0}}^{\text{'}}}{\sigma \mathrm{\text{'}}}\right)}^{\frac{C_{\mathrm{c}}}{C_{\mathrm{k}}}}\frac{\partial u}{\partial y}\right)}_{\left|i,j\right.}\approx

\frac{k_{i,j+\mathrm{1}}\frac{u_{i,j+\mathrm{2}}^n-u_{i,j}^n}{\mathrm{2}\mathrm{\Delta }y}-k_{i,j-\mathrm{1}}\frac{u_{i,j}^n-u_{i,j-\mathrm{2}}^n}{\mathrm{2}\mathrm{\Delta }y}}{{\left(\mathrm{\Delta }y\right)}^{\mathrm{2}}}（30）

对于\frac{\partial u}{\partial t}，采用向前差分格式，在节点\left(i,j\right)处有：

{\frac{\partial u}{\partial t}}_{}^{\left|i,jn\right.}\approx \frac{u_{i,j}^{n+\mathrm{1}}-u_{i,j}^n}{\mathrm{\Delta }t}（31）

将有限差分近似代入式（13）可得：

\frac{k_{i+\mathrm{1},j}\frac{u_{i+\mathrm{2},j}^n-u_{i,j}^n}{\mathrm{2}\mathrm{\Delta }x}-k_{i-\mathrm{1},j}\frac{u_{i,j}^n-u_{i-\mathrm{2},j}^n}{\mathrm{2}\mathrm{\Delta }x}}{{\left(\mathrm{\Delta }x\right)}^{\mathrm{2}}}+

\frac{k_{i,j+\mathrm{1}}\frac{u_{i,j+\mathrm{2}}^n-u_{i,j}^n}{\mathrm{2}\mathrm{\Delta }y}-k_{i,j-\mathrm{1}}\frac{u_{i,j}^n-u_{i,j-\mathrm{2}}^n}{\mathrm{2}\mathrm{\Delta }y}}{{\left(\mathrm{\Delta }y\right)}^{\mathrm{2}}}=

\frac{C_{\mathrm{c}}{\gamma }_{\mathrm{w}}}{\left(\mathrm{1}+e_{\mathrm{0}}\right){\sigma }_{\mathrm{0}}^{\text{'}}\mathrm{l}\mathrm{n}\mathrm{10}}\frac{u_{i,j}^{n+\mathrm{1}}-u_{i,j}^n}{\mathrm{\Delta }t}（32）

整理后可得到关于u_{i,j}^{n+\mathrm{1}}的方程：

u_{i,j}^{n+\mathrm{1}}=u_{i,j}^n+\frac{\mathrm{\Delta }t\left(\mathrm{1}+e_{\mathrm{0}}\right){\sigma }_{\mathrm{0}}^{\text{'}}\mathrm{l}\mathrm{n}\mathrm{10}}{C_{\mathrm{c}}{\gamma }_{\mathrm{w}}}\times

\left[\frac{k_{i+\mathrm{1},j}\frac{u_{i+\mathrm{2},j}^n-u_{i,j}^n}{\mathrm{2}\mathrm{\Delta }x}-k_{i+\mathrm{1},j}\frac{u_{i-\mathrm{2},j}^n-u_{i,j}^n}{\mathrm{2}\mathrm{\Delta }x}}{{\left(\mathrm{\Delta }x\right)}^{\mathrm{2}}}+\right.

\left.\frac{k_{i,j+\mathrm{1}}\frac{u_{i,j+\mathrm{2}}^n-u_{i,j}^n}{\mathrm{2}\mathrm{\Delta }y}-k_{i,j-\mathrm{1}}\frac{u_{i,j}^n-u_{i,j-\mathrm{2}}^n}{\mathrm{2}\mathrm{\Delta }y}}{{\left(\mathrm{\Delta }y\right)}^{\mathrm{2}}}\right]（33）

在给定初始条件u_{i,j}^{\mathrm{0}}［由式（22）］和边界条件［式（23）~式（26）］后，就可以通过时间迭代求解上述方程，得到不同时间步长下各节点的超静孔隙水压力u_{i,j}^n. 故地基中轴线的地表沉降和按孔压定义的平均固结度可分别由式（34）和式（35）计算得到.

S_{\mathrm{t}}={\int }_{-H}^{\mathrm{0}}{\epsilon }_y\left(x,y,t\right)\mathrm{d}y={\displaystyle \sum _{j=\mathrm{1}}^n}{\epsilon }_y\left(x,y,t\right)\mathrm{\Delta }{y}_j（34）

U_{\mathrm{t}}=\mathrm{1}-\frac{{\iint }_D^{}u\left(x,y,t\right)\mathrm{d}A}{{\iint }_D^{}u\left(x,y,\mathrm{0}\right)\mathrm{d}A}=\mathrm{1}-\frac{{\displaystyle \sum _{i=\mathrm{1}}^m}{\displaystyle \sum _{j=\mathrm{1}}^n}{u}_{i,j}{s}_{i,j}}{{\displaystyle \sum _{i=\mathrm{1}}^m}{\displaystyle \sum _{j=\mathrm{1}}^n}{u}_{\mathrm{0}\left(i,j\right)}{s}_{i,j}}（35）

式中：s_{i,j}=\mathrm{\Delta }{x}_i\mathrm{\Delta }{y}_j为离散点\left(x_i,y_j\right)处空间步长所围面积.

2 模型验证

根据有限差分格式和相关计算流程，编制Matlab计算程序，为验证计算模型和方法的正确性，选取已有解析解和有限元软件COMSOL Multiphysics的数值解与本文差分解进行比较分析.

本文解退化为一般的一维非线性土体固结解. 选取文献［25］算例中的土体基本参数进行计算. 将本文退化解计算结果与文献［25］中土体的非线性固结解析解计算结果进行对比，如图2所示.

将本文二维有限差分解计算结果与COMSOL的计算结果进行对比，基于两种方法计算的平均固结度随时间变化曲线如图2所示. 综合上述两种模型验证，可以看出两种模型具有很好的一致性，表明本文数值计算方法对该模型的计算是合理可靠的.

3 参数分析

按照式（14）和式（21）所定义的平均固结度U_{\mathrm{t}}和沉降S_{\mathrm{t}}进行计算并绘制变化曲线，计算时假定C_{\mathrm{c}}/C_{\mathrm{k}}=1.0，均布初始荷载q_{\mathrm{0}}=\mathrm{100}\mathrm{k}\mathrm{P}\mathrm{a}，初始孔隙比e_{\mathrm{0}}=\mathrm{1.4}，水平渗透系数k_x=\mathrm{4.17}\times {\mathrm{10}}^{-\mathrm{7}}\mathrm{m}/\mathrm{s}，竖直渗透系数k_y=\mathrm{4.17}\times {\mathrm{10}}^{-\mathrm{7}}\mathrm{m}/\mathrm{s}，定义K=k_x/k_y.

3.1 C_{\mathrm{c}}/C_{\mathrm{k}}对固结行为的影响

3.1.1 C_{\mathrm{c}}/C_{\mathrm{k}}对平均固结度的影响

绘制出C_{\mathrm{c}}/C_{\mathrm{k}}分别为0.50、0.75、1.00、1.25和1.50的情况下，平均固结度U_{\mathrm{t}}随时间的变化曲线，如图3所示.从图中可知，不同C_{\mathrm{c}}/C_{\mathrm{k}}下，相同时间时的平均固结度存在差异.在图3（a）中，当C_{\mathrm{c}}/C_{\mathrm{k}}为0.50、0.75、1.00、1.25和1.50时，100 h时的平均固结度分别达到70.92%、56.57%、47.37%、40.77%和35.83%；在图3（b）中，当C_{\mathrm{c}}/C_{\mathrm{k}}为0.50、1.00和1.50时，100 h时的平均固结度分别达到51.11%、47.37%和43.92%. 综合分析可得，C_{\mathrm{c}}/C_{\mathrm{k}}越大，平均固结度越低. 这表明，随着C_{\mathrm{c}}/C_{\mathrm{k}}的增大，意味着压缩性相对增强或渗透性相对减弱，导致孔隙水排出困难，非线性固结速率变慢. 进一步通过对比图3（a）和3（b）发现，当C_{\mathrm{c}}/C_{\mathrm{k}}从0.50提高到1.50时，改变C_{\mathrm{c}}时其平均固结度下降了35.09%，而改变C_{\mathrm{k}}时平均固结度仅下降了7.19%，由此可见在相同比值下，改变C_{\mathrm{c}}对非线性固结速率影响更为显著.

3.1.2 孔隙水压力等值线图

图4呈现了当C_{\mathrm{c}}/C_{\mathrm{k}}=1.00时，不同时间点的部分计算模型的超静孔隙水压分布等值线. 在固结早期阶段，超静孔隙水压的分布主要受荷载产生的附加应力作用和边界条件的制约. 在加载初期，一定深度下内部孔隙水压升高，形成超静孔压峰值. 在x>10 m的区域，受荷载影响较小，超静孔压较小，且在水平方向上，超静孔压随着距离的增大而逐渐减小. 随着时间推移，孔隙水压逐渐消散，其分布和大小也发生变化.例如，t=1 h时，最大孔隙水压值在80 kPa左右，t=10 h时下降至60 kPa左右，t=100 h时仅有40 kPa了.这清晰地展示了超静孔压在固结过程中的动态变化过程.

3.1.3 不同C_c/C_k对特殊点超静孔隙水压变化的影响

为研究C_c/C_k比值变化对地基内超静孔隙水压力的影响，选取位于中心上方的位置（0，-3）、下方的位置（0，-7）和底部的位置（0，-10）探究这些位置超静孔压的变化情况，如图5所示. 3 m处的孔压呈现逐渐消散的趋势，7 m和10 m处的孔压在固结初期会缓慢上升，这是由于荷载中心点下的超静孔压呈上大下小分布，如图5（a）所示，荷载中心点下3 m处的超静孔压约为80 kPa，7 m处的超静孔压约为 60 kPa，10 m处的超静孔压约为50 kPa，固结初期从高孔压向低孔压渗流，向下的渗流作用使得固结初期下部的超静孔压略有上升，达到峰值后开始消散.其中10 m处的孔压峰值出现的时间稍晚于7 m处，因为距离排水边界变远，排水路径变长，不仅孔隙水排出困难，上层的压缩也导致孔隙水向下迁移，使得10 m处孔压持续积累，峰值出现时间延迟.

从图5（a）（b）可知，不同C_c/C_k对初始超静孔压的分布和消散速率有影响，随着C_c/C_k比值增大，孔压消散变慢，且改变C_c时的影响比改变C_k时更为明显. 这进一步表明C_c/C_k对孔隙水压力的消散过程起着关键作用.

3.1.4 C_c/C_k对沉降的影响

绘制C_c/C_k分别为0.50、0.75、1.00、1.25和1.50的情况下，中轴线沉降随时间的变化曲线，如图6所示. 当C_k不变时，从图6（a）可以看出，不同C_c/C_k比值下瞬时沉降存在差异. 当C_c/C_k别为0.50、0.75、1.00、1.25、1.50时，瞬时沉降分别为0.003、0.005、0.006、0.008和0.009 m. 随着C_c/C_k比值增大，沉降速度提高，最终沉降量也增大，分别为0.23、0.35、0.46、0.57和0.68 m，但随着时间增加，不同比值下沉降速度趋于一致，说明在C_k不变的情况下，改变C_c/C_k的比值会显著影响初始瞬时沉降、沉降速度和最终沉降量. 当C_c不变时，从图6（b）可以看到，不同C_c/C_k比值下均有瞬时沉降，大概为0.006 m. 随着C_c/C_k比值的增大，沉降速度降低，但最终沉降值差距并不明显，均接近0.46 m. 这表明在C_c不变的情况下，改变C_c/C_k的比值主要影响沉降速度，对初始瞬时沉降和最终沉降量影响较小.

因此，渗透指数C_k主要影响沉降速度，可加快固结进度，进而使整个地基沉降加快；而压缩指数C_c对初始瞬时沉降、沉降速度和最终沉降量均有显著影响.

3.2 初始荷载q_{\mathrm{0}}对固结行为的影响

3.2.1 初始荷载q_{\mathrm{0}}对平均固结度的影响

图7展示了初始荷载q_{\mathrm{0}}对土体固结度的影响. 在C_c/C_k=0.50和1.50的情况下，分别以初始荷载q_{\mathrm{0}}=50、150和250 kPa进行分析，其对比变化曲线如图7所示. 通过曲线可以发现，当C_c/C_k=0.50时，三种初始荷载下100 h时的平均固结度均为59%左右，非线性固结速率无明显差距；当C_c/C_k=1.50时，三种初始荷载下100 h时的平均固结度分别为53.86%、48.32%和45.25%，差距有所扩大. 综合来看，随着C_c/C_k值增大，三种初始荷载下的固结曲线差距越来越明显，呈现出随着初始荷载q_{\mathrm{0}}增大，固结速率逐渐变小的趋势，说明C_c/C_k与初始荷载之间存在着交互作用，当C_c/C_k>1时，高荷载下土体压缩性显著增加，导致渗透性降低更明显，从而减缓固结速率.

3.2.2 初始荷载q_{\mathrm{0}}对沉降的影响

绘制初始荷载q_{\mathrm{0}}分别为50、150和250 kPa的情况下中轴线处沉降随时间的变化曲线，如图8所示. 从图中可以看出，三种初始荷载下的瞬时沉降分别为0.003、0.007和0.009 m. 随着初始荷载增大，沉降速度提高，最终沉降量也显著增大，分别达到0.199、0.465和0.648 m左右，说明初始荷载q_{\mathrm{0}}会影响初始瞬时沉降、沉降速度和最终沉降量.

3.3 渗透系数及其各向异性对固结行为的影响

3.3.1 水平渗透系数变化对平均固结度的影响

在竖向渗透系数k_y不变的情况下，通过改变水平渗透系数k_x来调整K值，并绘制K值分别为50、10、1.0、0.1和0.01时的平均固结度U_{\mathrm{t}}随时间的变化曲线，如图9所示. 不同K值下，相同时间时的平均固结度存在差异，比如1 000 h时，K为50、10、1.0、0.1和0.01时的平均固结度分别为97.59%、71.64%、53.41%、51.23%和50.87%，随着K从0.01提高到1，平均固结度稍有提高，而随着K从1提高到50，平均固结度提高幅度较大.这表明随着水平渗透系数k_x增大，非线性固结速率会加快，当水平渗透系数大于竖向渗透系数时，随着水平渗透系数的提高，其对非线性固结速率的影响愈发显著. 由此可见，水平渗透系数k_x对平均固结度有影响，且其影响程度和竖向渗透系数密切相关.

3.3.2 水平渗透系数变化对中轴线沉降的影响

在竖向渗透系数k_y不变的情况下，通过改变水平渗透系数k_x来调整K值，并绘制每种情况下中轴线沉降随时间的变化曲线，如图10所示. 从图中可知，五种K值下均有瞬时沉降，分别约为0.003 9、0.004 2、0.004 6、0.004 7和0.004 7 m，随着K值增大，沉降速度有所提高，比如在1 000 h时，五种K值下的沉降分别为0.34、0.30、0.22、0.19和0.19 m.这表明在水平渗透系数大于竖向渗透系数时，沉降速度提升更明显，沉降量也较大，但当水平渗透系数小于竖向渗透系数时，沉降速度提升不明显，沉降量变化也较小. 五种K值下最终沉降值差距不明显，均接近0.35 m.这说明保持k_y不变的情况下，改变K值主要影响沉降速度，对初始瞬时沉降和最终沉降量影响较小，且影响和竖向渗透系数相关.

3.3.3 竖向渗透系数变化对平均固结度的影响

在水平渗透系数k_x不变的情况下，通过改变竖向渗透系数k_y来调整K值，并绘制K分别为10、5.0、1.0、0.5和0.1时的平均固结度U_{\mathrm{t}}随时间的变化曲线，如图11所示. 不同K值下，相同时间时的平均固结度差异明显，比如1 000 h时，K为10、5.0、1.0、0.5和0.1时的平均固结度分别为19.63%、26.40%、53.41%、72.36%和99.66%，随着K减小，平均固结度显著提高.这表明竖向渗透系数k_y增大可有效提高非线性固结速率.

对比图9可知，竖向渗透系数对固结度的影响比水平渗透系数更大，因为前者直接影响主要排水方向，而后者的作用受限于侧向排水条件，故竖向渗透系数对非线性固结速率的影响更为关键.

3.3.4 竖向渗透系数变化对中轴线沉降的影响

在水平渗透系数k_x不变的情况下，通过改变竖向渗透系数k_y来调整K值，并绘制每种情况下中轴线沉降随时间的变化曲线，如图12所示. 首先五种K值下均有瞬时沉降，在0.005 m左右，随着K值减小，沉降速度提高，例如在1 000 h时，五种K值下的沉降分别为0.13、0.15、0.22、0.27和0.35 m，但是最终沉降值差距并不明显，最终沉降量趋势上接近0.35 m.这表明保持k_x不变的情况下，改变k_y主要影响沉降速度，对初始瞬时沉降和最终沉降量影响较小.

综上所述，渗透系数主要通过影响土体排水来加速沉降，竖向渗透系数k_y的提高对加快沉降的作用更为显著，对初始瞬间沉降和最终沉降量的影响相对有限.

3.4 地基深度对固结行为的影响

3.4.1 地基深度对平均固结度的影响

绘制在H分别为10、15、20、25和30 m情况下的平均固结度U_{\mathrm{t}}随时间的变化曲线，如图13所示.由图可知，在不同深度下，相同时间时的平均固结度不同，比如100 h时，H为10、15、20、25和30 m时的平均固结度分别为54.48%、39.81%、33.07%、29.22%和26.82%，随着深度增大，平均固结度逐渐降低，表明随着深度H的增大，非线性固结速率变慢. 这是由于浅层地基接近地表，排水条件较好，固结速率较快，而深度增大后加深的地基排水路径更长，排水条件差，固结速率慢.

3.4.2 地基深度对沉降的影响

绘制不同深度H时沉降随时间的变化曲线，如图14所示. 首先五种深度H下均有瞬时沉降，大概在0.005 m左右，而随着深度H的增大，沉降速度有所提高，比如在1 000 h时，五种深度下的沉降分别为0.35、0.45、0.52、0.55和0.58 m，且最终沉降值也不同，分别0.35、0.46、0.56、0.63和0.70 m. 随着每5 m的深度增加，其沉降值增加越少，故随着深度增加，总沉降量是增加的，但是每米的沉降是减少的，说明地基越深的地方竖向位移越小.

3.5 不同位置沉降量变化

为研究在该均布荷载下，不同位置的竖向沉降差距，分别记录距离中轴线0、6、9和11 m的沉降，并绘制沉降随时间变化曲线，如图15所示. 在200 h时， 四个位置沉降分别达到0.29、0.27、0.22和0.13 m，可以看出，在中轴线处的沉降最大，随着距离中轴线越来越远，沉降也在逐渐降低. 通过数据可知，在11 m处开始存在微量负沉降，即地表隆起，这是由于均 布荷载边缘处会因为土体的泊松比效应产生水平方向的膨胀，边缘处的土体受到的侧向应力较小，更容易向侧面移动或挤出，导致前期边缘处土体向上隆起.

4 结 论

本研究通过对饱和软土地基在条形荷载作用下的二维非线性固结问题的分析，构建了考虑土体渗透性和压缩性非线性关系的控制方程，探讨了多种因素对土体固结行为的影响，得出以下结论：

1）压缩指数与渗透指数的比值C_{\mathrm{c}}/C_{\mathrm{k}}对土体固结行为影响显著，C_{\mathrm{c}}/C_{\mathrm{k}}越大，固结速率和沉降速度越慢，C_{\mathrm{c}}在固结行为的影响作用大于C_{\mathrm{k}}，故遇到高压缩性黏土时，需要延长预压时间或采用排水措施加速固结.

2）初始荷载q_{\mathrm{0}}对土体固结过程的影响与C_{\mathrm{c}}/C_{\mathrm{k}}相关，当C_{\mathrm{c}}/C_{\mathrm{k}}>1，q_{\mathrm{0}}增大，会使非线性固结速率降低，初始瞬时沉降、沉降速度和最终沉降量都增加. 实际工程中需关注土体参数的非线性特性，合理设计加载方案，避免因加载过快或过大导致的不均匀沉降问题.

3）渗透系数及各向异性对固结行为影响明显，水平或竖向渗透系数越大，固结速率和沉降速度越快，竖向渗透系数k_y对固结速率的影响大于水平渗透系数k_x. 工程中可通过竖向排水井显著加速固结.

4）地基深度H越大，固结速率越慢，瞬时沉降不变，最终沉降量变大，但每米的沉降量变小，即深层地基竖向位移较小，故可通过表层加固减少差异沉降，深层地基需要预留充足的工后沉降期.

5）在均布荷载作用下，地基不同位置沉降量不同，中轴线处的沉降最大，需要重点监测，距离中轴线越远，沉降逐渐降低，边缘处可能存在地表隆起，需要关注横向位移.

参考文献

[1] TERZAGHI V K ．Erdbaumechanik：auf bodenphysikalischer grundlage［M］．Leipzig： Franz. Deuticke，1925．

[2] KIM P， KIM H S， PAK C U，et al ．Analytical solution for one-dimensional nonlinear consolidation of saturated multi-layered soil under time-dependent loading［J］. Journal of Ocean Engineering and Science， 2021， 6（1）： 21-29．

[3] CHEN R P， ZHOU W H， WANG H Z，et al. One-dimensional nonlinear consolidation of multi-layered soil by differential quadrature method［J］. Computers and Geotechnics， 2005， 32（5）： 358-369．

[4] OLIVEIRA F S， MARTINS I S M， GUIMARÃES L J N ．A model for one-dimensional consolidation of clayey soils with non-linear viscosity［J］．Computers and Geotechnics，2023，159：105426．

[5] 李传勋， 仇超. 高压缩性软土一维非线性大应变固结解析解［J］. 岩石力学与工程学报，2021，40（11）：2344-2356．

[6] LI C X， QIU C. An analytical solution for one-dimensional nonlinear large-strain consolidation of soft clay with high compressibility［J］. Chinese Journal of Rock Mechanics and Engineering，2021，40（11）：2344-2356．（in Chinese）

[7] 江留慧，李传勋，杨怡青，等 ．变荷载下双层地基一维非线性固结近似解析解［J］．岩土力学，2020，41（5）：1583-1590．

[8] JIANG L H， LI C X， YANG Y Q，et al ．Approximate analytical solutions for one-dimensional nonlinear consolidation of double-layered soil under time-dependent loading［J］．Rock and Soil Mechanics，2020，41（5）：1583-1590．（in Chinese）

[9] 夏长青，胡安峰，崔军，等 ．饱和软土成层地基一维非线性固结解析解［J］．岩土力学，2018，39（8）：2858-2864．

[10] XIA C Q， HU A F， CUI J，et al ．Analytical solutions for one-dimensional nonlinear consolidation of saturated soft layered soils［J］. Rock and Soil Mechanics，2018，39（8）：2858-2864．（in Chinese）

[11] 张添，汪磊，沈思东 ．分段循环荷载作用下二维非饱和土固结特性分析［J］．工程地质学报，2022，30（4）：1010-1018．

[12] ZHANG T， WANG L， SHEN S D ．Analysis of two-dimensional plane strain consolidation characteristics for unsaturated soils under piece-wise cyclic loading［J］．Journal of Engineering Geology，2022，30（4）：1010-1018．（in Chinese）

[13] 胡安峰，龚昭祺，肖志荣，等 ．隧道周围饱和软土二维非线性固结分析［J］．上海交通大学学报，2023，57（12）：1631-1638．

[14] HU A F， GONG Z Q， XIAO Z R，et al ．Two dimensional nonlinear consolidation analysis of saturated soft soil around tunnel［J］. Journal of Shanghai Jiaotong University， 2023，57（12）： 1631-1638．（in Chinese）

[15] 刘奔奔 ．堆载作用下层状地基二维固结理论研究［D］．济南：山东建筑大学，2022．

[16] LIU B B ．Research on two-dimensional consolidation theory of layered foundations under surcharge loading［D］. Jinan： Shandong Jianzhu University， 2022．（in Chinese）

[17] 黄传志，肖原 ．二维固结问题的解析解［J］．岩土工程学报，1996， 18（3）： 47-54．

[18] HUANG C Z， XIAO Y ．Analytic solution of a two dimensional consolidation problem［J］. Chinese Journal of Geotechnical Engineering， 1996， 18（3）： 47-54．（in Chinese）

[19] 艾智勇，慕金晶 ．竖向简谐荷载下二维层状饱和地基的解析层元解［J］．岩土力学，2018，39（7）：2632-2638．

[20] AI Z Y， MU J J ．Analytical layer element solution for two-dimensional multilayered saturated subsoils under a vertical time-harmonic load［J］．Rock and Soil Mechanics， 2018，39（7）： 2632-2638．（in Chinese）

[21] 关聪 ．考虑Hansbo渗流的二维Biot固结分析［D］．郑州：郑州大学，2016．

[22] GUAN C ．Analysis of two-dimensional Biot’s consolidation with hans Hansbo’s flow［D］. Zhengzhou：Zhengzhou University，2016．（in Chinese）

[23] WANG L， XU Y F， XIA X H，et al ．Semi-analytical solutions to two-dimensional plane strain consolidation for unsaturated soil［J］．Computers and Geotechnics，2018，101：100-113．

[24] WANG L， XU Y F， XIA X H，et al ．Semi-analytical solutions to two-dimensional plane strain consolidation for unsaturated soils under time-dependent loading［J］．Computers and Geotechnics，2019，109：144-165．

[25] WANG L， XU Y F， XIA X H，et al ．Semi-analytical solutions to the two-dimensional plane strain consolidation for unsaturated soil with the lateral semi-permeable drainage boundary under time-dependent loading［J］．Computers and Geotechnics，2020，124：103562．

[26] SHEN S D， WANG L， ZHOU A N，et al ．Two-dimensional plane strain consolidation for unsaturated soils with a strip-shaped distributed permeable boundary［J］．Computers and Geotechnics，2021，137：104273．

[27] SUN C S， YU T W， LIU B B，et al. Semianalytic solution for two-dimensional consolidation of interbedded soil［J］. International Journal of Geomechanics， 2024， 24（9） ：04024203．

[28] WANG L， SHEN S D， LI T Y，et al. Two-dimensional plane strain consolidation of unsaturated soils considering the depth-dependent stress［J］. Journal of Rock Mechanics and Geotechnical Engineering，2023，15（6）： 1603-1614．

[29] BUTTERFIELD R. A natural compression law for soils （an advance on e-log p'）［J］. Géotechnique， 1979， 29（4）： 469-480．

[30] 谢康和，齐添， 胡安峰，等 ．基于GDS的黏土非线性渗透特性试验研究［J］．岩土力学，2008，29（2）：420-424．

[31] XIE K H， QI T， HU A F，et al ．Experimental study on nonlinear permeability characteristics of Xiaoshan clay［J］．Rock and Soil Mechanics， 2008， 29（2）： 420-424．（in Chinese）

[32] 建筑地基基础设计规范：GB 50007—2011 ［S］．北京：中国计划出版社，2012．

[33] Code for design of building foundation：GB 50007—2011 ［S］．Beijing： China Planning Press，2012．（in Chinese）

[34] CUI P L， CAO W G， LIU Y Y，et al ．One-dimensional nonlinear creep consolidation of soft soils with time-dependent drainage boundary under construction load［J］．International Journal for Numerical and Analytical Methods in Geomechanics， 2023， 47（9）： 1612-1636．

[35] 弗洛林 . 土力学原理， 第2卷 ： 建筑物地基的变形与稳定 性［M］. 徐志英， 译. 北京： 中国建筑工业出版社，1973.

[36] FLORIN. Principles of soil mechanics， Volume 2： deformation and stability of building foundations ［M］. Translated by XU Z Y. Beijing： China Architecture & Building Press， 1973.（in Chinese）

[37] LI C X， HUANG J S， WU L Z，et al ．Approximate analytical solutions for one-dimensional consolidation of a clay layer with variable compressibility and permeability under a ramp loading［J］．International Journal of Geomechanics，2018，18（11）：06018032．

基金资助

国家自然科学基金资助项目(52078206)

National Natural Science