沿海地区由于饱和软黏土地基的高压缩性，在荷载作用下建筑容易产生较大沉降，因此, 需对地基进行预先处理，以加快地基的固结速率，其中设计理论大多基于砂井固结理论.Barron^[1]首先对砂井地基给出了等应变和自由应变两种假定下径向固结度的计算公式，并指出当固结时间较长时，两者假定下的计算差异很小.在此基础上，Hansbo等^[2-3]分别研究了井阻作用对砂井固结的影响.Tang等^[4-6]在等应变条件下分别考察了变荷载、影响区形状等因素对砂井地基固结机制的影响，进一步完善了砂井固结理论.同时，为得到更为精确的固结解答，Berry等^[7-9]考虑土体非线性压缩关系，采用不同关系式修正砂井径向固结方程，并给出解析解答.然而，上述研究均忽略了软黏土的流变特性和渗流的非达西特性，因而地基沉降的理论值与实际观测值常存在较大差异.

固结试验和现场观测^[10-12]表明，软黏土具有显著的流变性.为了更好地反映土体流变固结机制，Taylor等^[13]首先采用Kelvin模型分析土体的主次固结耦合过程，随后有学者引入了Merchant模型^[14]、广义Voigt模型^[15]等.但由于元件模型通过基本流变元件的组合而成，本构方程有时含有较多参数，且形式比较复杂.因此，刘忠玉等^[16-17]引入了分数阶导数模型.同时，有学者基于经典的弹塑性理论，提出了弹黏塑性本构模型.其中，Yin等^[18]引入等效时间的概念，给出了一个非线性的EVP本构模型，较好地描述了土体的应力-应变-时间关系.姚仰平等^[19-20]则在修正剑桥模型的基础上，采用统一硬化参数，得到了考虑时间效应的统一硬化(UH)模型.上述理论多应用于一维固结情况，同时也开始应用于砂井地基等多维固结分析.例如，刘兴旺等^[21-23]先后引入Merchant模型和四元件模型描述砂井地基黏弹性固结问题，并分析了黏弹性参数对固结过程的影响.但目前还少有文献对砂井地基的弹黏塑性变形特性进行深入分析.

另外，已有试验^[24-26]表明，某些黏性土中的渗流速度与水力梯度之间存在非线性关系，有别于Darcy渗流，称之为非Darcy渗流.为深入探讨非Darcy渗流下土体的固结机制，学者们^[27-29]分别引入指数渗流、起始水力梯度渗流、Hansbo渗流等模式来描述砂井地基的非Darcy渗流，并通过与Darcy渗流下的结果对比，分析相应参数对地基整体孔压消散速率的影响.最近，刘忠玉等^[30-31]分别在传统Terzaghi固结理论和经典Biot固结理论的基础上，引入Hansbo渗流给予修正，分析了Hansbo渗流参数对多维弹性饱和地基固结的变化规律.不过上述文献并没有考虑土体应力-应变的非线性关系，因而属于线弹性理论范畴.同时需要注意的是，在土体固结时，随着孔隙比的减小，渗透系数也应当是变化的.Taylor^[32]首先提出了渗透系数与孔隙比的经验关系式，之后不断有学者^[33-34]通过试验研究影响渗透系数变化的主要因素.同时，也不断有学者^[35-36]讨论渗透系数变化对固结过程的影响.但是，这些固结分析仅限于渗流或变形的非线性方面，并未考虑土体的流变特性，并不能准确地分析多维饱和软黏土地基的固结过程.

总之，上述研究成果丰富和发展了对固结机理的理解.但相对而言，在多维固结分析中同时考虑渗流非Darcy特性、渗透系数变化及土体流变特性等方面影响的报道还较少.为此，在自由应变条件下，拟引入考虑时间效应的UH本构模型和Hansbo渗流模型，同时考虑变渗透系数的影响，重新推导理想砂井地基固结方程，并初步探讨相关参数对软黏土地基流变固结特性的影响机制.

1 分析模型及推导控制方程 1.1 应力-应变关系

在考虑时间效应的UH模型^[19]中，体积应变增量的表达式为

$ \mathrm{d} \varepsilon_{\mathrm{v}}=\left\{\begin{array}{l} \mathrm{d} \varepsilon_{\mathrm{v}}^{\mathrm{e}}+\mathrm{d} \varepsilon_{\mathrm{v}}^{\mathrm{sp}}+\mathrm{d} \varepsilon_{\mathrm{v}}^{\mathrm{tp}}, \mathrm{d} \sigma^{\prime} \geqslant 0 ,\\ \mathrm{~d} \varepsilon_{\mathrm{v}}^{\mathrm{e}}+\mathrm{d} \varepsilon_{\mathrm{v}}^{\mathrm{tp}}, \mathrm{d} \sigma^{\prime}<0, \end{array}\right. $

(1)

$ \mathrm{d} \varepsilon_{\mathrm{v}}^{\mathrm{e}}=C_{\mathrm{s}} \mathrm{d} \sigma^{\prime} /\left[\left(1+e_{0}\right) \sigma^{\prime} \ln 10\right], $

(2)

$ \mathrm{d} \varepsilon_{\mathrm{v}}^{\mathrm{sp}}=\frac{C_{\mathrm{c}}-C_{\mathrm{s}}}{\left(1+e_{0}\right) \ln 10} \frac{M^{4}}{M_{\mathrm{f}}^{4}} \frac{\mathrm{d} \sigma^{\prime}}{\sigma^{\prime}}, $

(3)

$ \mathrm{d} \varepsilon_{\mathrm{v}}^{\mathrm{tp}}=C_{\mathrm{\alpha}} \mathrm{d} t /\left[\left(1+e_{0}\right)\left(t_{\mathrm{a}}+t_{0}\right) \ln 10\right]. $

(4)

式中：dε_v^e和dε_v^sp分别表示有效应力增量dσ′作用下产生的弹性体应变增量和塑性体应变增量，dε_v^tp表示时间作用产生的塑性体应变增量，t_a和t₀分别表示老化时间和单位时间，C_c为压缩指数，C_s为回弹指数，C_α为次固结系数，M和M_f分别为临界状态应力比和潜在破坏应力比，且

$ M =6 \sin \varphi /(3-\sin \varphi), $

(5)

$ M_{\mathrm{f}}=6\left[\sqrt{\frac{\chi}{R}\left(1+\frac{\chi}{R}\right)}-\frac{\chi}{R}\right], $

(6)

$ \chi=M^{2} / 12(3-M) $

(7)

$ R=\frac{\sigma^{\prime}}{{\sigma^{\prime}}_{\mathrm{c} 0}} \cdot \exp \left.[-\left(1+e_{0}\right) \varepsilon_{\mathrm{v}}^{\mathrm{p}} \ln 10 /\left(C_{\mathrm{c}}-C_{\mathrm{s}}\right)\right], $

(8)

$ \left(t_{\mathrm{a}}+t_{0}\right) / t_{0}=R^{-\alpha}, $

(9)

$ \alpha =\left(C_{\mathrm{c}}-C_{\mathrm{s}}\right) / C_{\mathrm{\alpha}}. $

(10)

其中：φ为土的内摩擦角；R表示超固结参数；σ′_c0类似先期固结压力，若不考虑时间效应即为先期固结压力，此时超固结参数的初值R₀为超固结比的倒数；e₀为初始孔隙比；ε_v^p=ε_v^sp+ε_v^tp为塑性体应变.

1.2 控制方程推导

设某均质饱和软黏土层厚度为H，底部不透水，顶部透水，并且在自重等荷载作用下已经固结完成.其中均匀设置打穿软土层的砂井，其半径为r_w，影响区半径为r_e.现于砂井地基表面瞬时施加均布荷载p₀，并假定在固结过程中土体不发生侧向变形^[5].为便于计算，设初始有效应力σ′₀和与之对应的初始孔隙比e₀沿土层深度不变^[37].在半径r深度z处，时刻t的有效应力和超孔隙水压力(以下简称孔压)分别为σ′(r, z, t)和u(r, z, t)，渗流速度为$\overrightarrow v = \left( {{v_r}, {v_z}} \right)$，水力坡降为$\overrightarrow i = \left( {{i_r}, {i_z}} \right)$.假定渗流的非Darcy特性可用Hansbo渗流模型来描述，即

$ \vec{v}=\left\{\begin{array}{ll} K|\vec{i}|^{m-1} /\left(m i_{1}^{m-1}\right) \cdot \vec{i}, & |\vec{i}| \leqslant i_{1}, \\ K\left[1-i_{1}(m-1) /(m|\vec{i}|)\right] \cdot \vec{i}, & |\vec{i}|>i_{1}. \end{array}\right. $

(11)

式中：m为由试验确定的常数；i₁为直线段起始水力梯度；K为渗透系数，这里采用Taylor经验关系式来描述其与孔隙比e的关系，即

$ \log K=\log K_{0}-\left(e_{0}-e\right) / C_{\mathrm{k}}. $

(12)

式中：K₀为初始孔隙比e₀对应的初始渗透系数；C_k为渗透指数，反映渗透系数的变化程度.

考虑体积应变ε_v与孔隙比e的关系，可得

$ K=K_{0} \exp \left[-\left(1+e_{0}\right) \varepsilon_{\mathrm{v}} \ln 10 / C_{\mathrm{k}}\right]. $

(13)

为便于后面计算，将式(11)改写成形如Darcy定律的表达式

$ \vec{v}=K^{\prime} \vec{i}. $

(14)

式中K′称为视渗透系数，即

$ \begin{array}{c} K^{\prime}=\left\{\begin{array}{ll} K|\vec{i}|^{m-1} /\left(m i_{1}^{m-1}\right), & |\vec{i}| \leqslant i_{1}, \\ K\left[1-i_{1}(m-1) /(m|\vec{i}|)\right], & |\vec{i}|>i_{1} . \end{array}\right. \end{array} $

(15)

式中：$|\vec i| = \sqrt {i_r^2 + i_z^2} ;{i_r} = - \frac{1}{{{\gamma _w}}}\frac{{\partial u}}{{\partial r}};{i_z} = - \frac{1}{{{\gamma _w}}}\frac{{\partial u}}{{\partial z}}$.

考虑土体渗流连续性条件，可得

$ \frac{\partial \varepsilon_{\mathrm{v}}}{\partial t}=\frac{1}{r} \frac{\partial}{\partial r}\left(r v_{r}\right)+\frac{\partial v_{z}}{\partial z}. $

(16)

根据有效应力原理，可得

$ \frac{\partial u}{\partial r}=-\frac{\partial \sigma^{\prime}}{\partial r}, \frac{\partial u}{\partial z}=-\frac{\partial \sigma^{\prime}}{\partial z}. $

(17)

由式(1)、(11)、(16)和(17)联立可得

$ \frac{1}{r} \frac{\partial}{\partial r}\left(\psi r \frac{\partial \sigma^{\prime}}{\partial r}\right)+\frac{\partial}{\partial z}\left(\psi \frac{\partial \sigma^{\prime}}{\partial z}\right)=\beta \frac{\partial \sigma^{\prime}}{\partial t}+\lambda . $

(18)

式中：ψ=K′/γ_w；λ=C_α/[(1+e₀)(t_a+t₀)ln 10]，由于该参数同时含有次固结系数和老化时间，该参数反映了土体流变性和超固结特性；β可表示为

$ \beta=\left\{\begin{array}{ll} \frac{C_{\mathrm{s}}+\left(C_{\mathrm{c}}-C_{\mathrm{s}}\right) M^{4} / M_{\mathrm{f}}^{4}}{\left(1+e_{0}\right) \sigma^{\prime} \ln 10}, & \mathrm{~d} \sigma^{\prime} \geqslant 0, \\ \frac{C_{\mathrm{s}}}{\left(1+e_{0}\right) \sigma^{\prime} \ln 10}, & \mathrm{~d} \sigma^{\prime}<0. \end{array}\right. $

(19)

本课题的初始条件及边界条件如下:

$ \sigma^{\prime}(r, z, 0)={\sigma^{\prime}}_{0}, r_{w} \leqslant r \leqslant r_{e}, 0 \leqslant z \leqslant H , $

(20)

$ \left\{\begin{array}{l} \sigma^{\prime}\left(r_{\rm{w}}, z, t\right)=\sigma_{0}^{\prime}+p_{0}, 0 \leqslant z \leqslant H, t>0, \\ \sigma^{\prime}(r, 0, t)=\sigma_{0}^{\prime}+p_{0}, r_{\rm{w}} \leqslant r \leqslant r_{\rm{e}}, t>0, \end{array}\right. $

(21)

$ \left\{\begin{array}{l} \frac{\partial \sigma^{\prime}\left(r_{\rm{e}}, z, t\right)}{\partial r}=0,0 \leqslant z \leqslant H, t>0, \\ \frac{\partial \sigma^{\prime}(r, H, t)}{\partial z}=0, r_{\rm{w}} \leqslant r \leqslant r_{\rm{e}}, t>0. \end{array}\right. $

(22)

2 方程的有限差分解法

由于求解上述非线性微分方程解析解比较困难，采用隐式有限差分法进行求解.首先，在r_w≤r≤r_e内，将土层沿径向以Δr为步长由内向外离散，节点从0到N编号；在0≤z≤H内，土层沿竖向以Δz为步长从上往下离散，节点编号从0到Z.同时以Δt为时间步长离散时间.这样，方程(18)离散为

$ \begin{array}{c} {\sigma^{\prime}}_{i, j}^{\prime k+1}={\sigma^{\prime}}_{i, j}^{\prime k}+\varphi_{i^{\prime}, j}^{k+1}\left(\psi_{i^{\prime}+1 / 2, j}^{k+1} r_{i^{\prime}+1 / 2} A_{i^{\prime}, j}^{k+1}-\psi_{i^{\prime}-1 / 2, j}^{k+1}\right. \cdot \\ \left.r_{i^{\prime}-1 / 2} B_{i^{\prime}, j}^{k+1}\right)+\eta_{i^{\prime}, j}^{k+1}\left(\psi_{i^{\prime}, j+1 / 2}^{k+1} C_{i^{\prime}, j}^{k+1}-\psi_{i^{\prime}, j-1 / 2}^{k+1} D_{i^{\prime}, j}^{k+1}\right)+E_{i^{\prime}, j}^{k+1}. \end{array} $

(23)

式中：x_{i^′, j}^k表示变量x在kΔt时刻，点(r_w+i^′Δr，jΔz)处的值，且${A_{{i^\prime },j}^{k + 1} = \sigma _{{i^\prime } + 1,j}^{\prime k + 1} - \sigma _{{i^\prime },j}^{\prime k + 1};B_{{i^\prime },j}^{k + 1} = \sigma _{{i^\prime },j}^{k + 1} - \sigma _{{i^\prime } - 1,j}^{\prime k + 1};\varphi _{{i^\prime },j}^{k + 1} = {a_1}/\left( {\beta _{{i^\prime },j}^{k + 1}{r_{{i^\prime }}}} \right),{a_1} = \Delta t/{{(\Delta r)}^2};\eta _{{i^\prime },j}^{k + 1} \\= {a_2}/\beta _{{i^\prime },j}^{k + 1},{a_2} = \Delta t/{{(\Delta z)}^2};\psi _{{i^\prime } \pm 1/2,j}^{k + 1} = K_{{i^\prime } \pm 1/2,j}^{k + 1}/{\gamma _w};{r_{{i^\prime }}} = {r_{\rm{w}}} + {i^\prime }\Delta r;{r_{{i^\prime } \pm 1/2}} = {r_{\rm{w}}} + \left( {{i^\prime } \pm 1/2} \right)\Delta r;E_{i,j}^{k + 1}\\ = - \lambda _{{i^\prime },j}^{k + 1}\Delta t/\beta _{{i^\prime },j}^{k + 1};C_{{i^\prime },j}^{k + 1} = \sigma _{{i^\prime },j + 1}^{\prime k + 1} - \sigma _{{i^\prime },j}^{\prime k + 1};D_{{i^\prime },j}^{k + 1} = \sigma _{{i^\prime },j}^{\prime k + 1} - \sigma _{{i^\prime },j - 1}^{k + 1}.}$

同时，初始条件(20)，边界条件(21)和(22)的离散结果如下：

$ {\sigma^{\prime}}_{i, j}^{0}={\sigma^{\prime}}_{0}, i^{\prime}=0,1,2, \cdots, N ; j=0,1,2, \cdots, Z. $

(24)

$ \left\{\begin{array}{l} {\sigma^{\prime}}_{i, 0}^{\prime k}={\sigma^{\prime}}_{0}+p_{0}, i^{\prime}=0,1, \cdots, N ; k=1,2, \cdots, \\ {\sigma^{\prime}}_{0, j}^{ k}={\sigma^{\prime}}_{0}+p_{0}, j=0,1, \cdots, Z ; k=1,2, \cdots. \end{array}\right. $

(25)

$ \begin{array}{l} {\sigma^{\prime}}_{N,j}^{k + 1} = {\sigma^{\prime}}_{N,j}^k - \varphi _{N,j}^{k + 1}\left( {\psi _{N + 1/2,j}^{k + 1}{r_{N + 1/2}} + \psi _{N - 1/2,j}^{k + 1}} \right. \cdot \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \left. {{r_{N - 1/2}}} \right)B_{N,j}^{k + 1} + E_{N,j}^{k + 1} + \eta _{N,j}^{k + 1}\left( {\psi _{N,j + 1/2}^{k + 1}} \right. \cdot \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \left. {C_{N,j}^{k + 1} - \psi _{N,j - 1/2}^{k + 1}D_{N,j}^{k + 1}} \right). \end{array} $

(26)

$ \begin{array}{l} {\sigma^{\prime}}_{{i^\prime },Z}^{k + 1} = {\sigma^{\prime}}_{{i^\prime },Z}^k + \varphi _{{i^\prime },Z}^{k + 1}\left( {\psi _{{i^\prime } + 1/2,Z}^{k + 1}{r_{{i^\prime } + 1/2}}A_{{i^\prime },Z}^{k + 1} - } \right.\\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \left. {\psi _{{i^\prime } - 1/2,j}^{k + 1}{r_{{i^\prime } - 1/2}}B_{{i^\prime },Z}^{k + 1}} \right) + E_{{i^\prime },Z}^{k + 1} - \eta _{{i^\prime },Z}^{k + 1} \cdot \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \left( {\psi _{{i^\prime },Z + 1/2}^{k + 1} + \psi _{{i^\prime },Z - 1/2}^{k + 1}} \right)D_{{i^\prime },Z}^{k + 1}. \end{array} $

(27)

这样，由式(23)~(27)构成封闭的方程组，用迭代法可解得各节点的有效应力，然后利用有效应力原理可得其孔压值.最后，引入按孔压定义的砂井地基的平均固结度式(28)，来定量地描述地基中孔压的整体消散情况.

$ {U_{\rm{P}}} = 1 - \frac{{\int_0^H {{\rm{d}}} z\int_{{r_{\rm{w}}}}^{{r_{\rm{e}}}} 2 {\rm{ \mathsf{ π} }}ru{\rm{d}}r}}{{\int_0^H {{\rm{d}}} z\int_{{r_{\rm{w}}}}^{{r_{\rm{e}}}} 2 {\rm{ \mathsf{ π} }}r{p_0}{\rm{d}}r}} = 1 - \frac{{2\Delta z\sum\limits_{j = 0}^{Z - 1} {\sum\limits_{{i^\prime } = 0}^{N - 1} {S_{{i^\prime },j}^k} } }}{{{p_0}H\left( {r_{\rm{e}}^2 - r_{\rm{w}}^2} \right)}}. $

(28)

式中：${S_{{i^\prime },j}^k = \frac{{u_{{i^\prime } + 1,j}^k - u_{{i^\prime },j}^k}}{{3\Delta r}}\left( {r_{{i^\prime } + 1}^3 - r_{{i^\prime }}^3} \right) + \left\{ {\frac{1}{4}\left( {u_{{i^\prime },j + 1}^k - u_{{i^\prime },j}^k} \right) + } \right.\left. {\frac{1}{2}\left[ {u_{{i^\prime },j}^k - \left( {u_{{i^\prime } + 1,j}^k - u_{{i^\prime },j}^k} \right){r_{{i^\prime }}}/\Delta r} \right]} \right\}\left( {r_{{i^\prime } + 1}^2 - r_{{i^\prime }}^2} \right).}$

最后，将上述算法编成了Fortran程序.

3 解法验证

令m=1或i₁=0，M=M_f，C_α=0，C_c=C_k，本课题将退化成耿雪玉^[37]在自由应变假定下曾讨论过的符合Darcy定律的砂井地基非线性固结问题.为便于对比，引入如下无量纲参数：Z=z/H，X=r/r_w，U=u/p₀，T_v=[(μ₁²+λ₁²L)C_h0t]/r_w²，N_σ=(σ′₀+p₀)/σ′₀，C_h0=K₀(1+e₀)σ′₀ln 10/(γ_wC_c)，同时取C_k=0.6，C_c=0.6，L=0.01，e₀=1.5，K₀=6.0×10^-7m/min，σ′₀=50 kPa，λ₁=π/2，μ₁=0.28，r_w=0.2 m，r_e=1.0 m，Δr=0.01 m，H=2 m，Δz=0.01 m，Δt=0.1 min，迭代精度为10^-7.对应N_σ=1，2，4和8，p₀分别取0，50，150和350 kPa，按本文方法得到Z=1，X=2.5处的孔压U随时间因素T_v的变化曲线示于图 1，同时图 1中也附上了耿雪玉^[37]的解析解.可见两者基本一致，这说明本文算法是有效的.

4 参数分析

为探讨地基土的流变效应、非达西特性及渗透系数变化等因素对砂井地基固结进程的影响，以下采用单因素法进行分析，其中参数分别取σ′₀=50 kPa，p₀=150 kPa，K₀=6.0×10^-7 m/min，γ_w=10 kN/m³, M=0.772，C_c=0.6，C_α=0.04，C_s=0.12，C_k=0.8，R₀=0.75，e₀=1.5，m=1.5，i₁=5，r_w=0.05 m，r_e=1.0 m，N=50，H=10 m，Z=30，Δt=0.1 min，迭代误差为10^-7.

4.1 次固结系数C_α的影响

图 2给出次固结系数C_α对某些特殊位置处孔压的影响曲线.可以看出，在流变固结初期(t=1 000 min)，C_α>0时，远离排水边界处的孔压都超过了初始值p₀ =150 kPa，这说明该阶段这些位置处的孔压不是消散而是升高.并且随着次固结系数C_α的增大，这种现象变得更加明显.如图 2(a)所示，当t=1 000 min，且C_α分别等于0.02，0.04和0.06时，z=5 m，r=0.9 m处的孔压值分别为153，164，173 kPa，即大约是初始孔压的1.02，1.09，1.15倍.这种孔压升高现象不同于二维或三维Biot固结课题中的Mandel-Cryer效应，因为本课题仍属于Terzaghi固结理论范畴，刘忠玉等^[38]称之为“类曼德尔效应”.实际上一维固结初期也会出现该现象^[18].Yin等^[18]将其归结为土的黏滞效应或主次固结耦合效应，认为在土体蠕变过程中，由于靠近不排水面处，土体排水不畅，孔隙水出现拥挤情况，导致远离不排水面处的孔压升高，有效应力减小.这也可用来解释本文砂井地基固结初期的孔压升高现象.

同时由图 2(a)可知，在固结中后期，孔压的消散速率随着C_α的增大而变得明显缓慢.例如，当不考虑黏滞性时，在z=5 m，r=0.9 m处的孔压值为12 kPa；当考虑黏滞性时，对应C_α等于0.02，0.04和0.06，同一位置处的孔压值分别为27，47和76 kPa，分别为不考虑黏滞性时的2.25，3.92和6.33倍.很明显，土体的黏滞效应延缓了地基的固结过程.

图 3给出不同次固结系数C_α时的砂井地基平均固结度U_p随时间t的变化曲线.可以看出，在流变固结初期，考虑黏滞效应(C_α >0)的曲线都出现了上升段(固结度甚至为负).如t=1 000 min时，与C_α=0，0.02，0.04，0.06对应的平均固结度值分别为0.080，0.059，-0.016，-0.077.这表明随着土体次固结系数的增大，地基中较大范围内都出现了孔压升高的现象.同时，考虑黏滞效应的平均固结度变化曲线都在不考虑土体黏滞效应(C_α=0)计算曲线的上方，也说明了土体黏滞性对砂井地基中整体孔压消散具有滞后作用.这与胡晶等^{[20, 38]}对一维流变固结分析时所得出的结论一致.

4.2 初始超固结参数R₀的影响

图 4给出了t=1 000和100 000 min两个时刻下不同初始超固结参数R₀时特殊位置处的孔压分布曲线.可以看出，在流变固结过程中，R₀越大(即初始超固结程度越低)，孔压越大，且类曼德尔效应越明显.例如图 4(b)中，t=1 000 min时，与R₀取0.50，0.75和1.00对应的z=1 m，r=0.5 m处的孔压值分别等于145，160和170 kPa.这与胡晶等^[20]在一维流变固结下得出的结论一致，他们认为土体所受的应力历史不同，或者说前期固结应力不同，固结速度也会有所不同，即增大超固结比(或减小初始超固结参数R₀)，会导致主固结与次固结的变形同时有所减小，从而加快孔压的消散速率.

图 5给出初始超固结参数R₀对地基平均固结度U_p的影响曲线.与图 3类似，图 5表明，平均固结度随时间并非都逐渐增大，而是在固结前期，固结度可有所减小，到某个时刻达最小值，且随着R₀的增大，固结度减小得更为明显.例如，当t=1 000 min，R₀=0.50，0.75，1.00时，对应的平均固结度值分别为0.072，-0.016，-0.070.从图 5整体来看，R₀取值越小，U_p-t曲线就越靠左下方，可见初始超固结参数R₀越小(即初始超固结程度越强)，整体孔压消散越快.这里也与胡晶等^{[20, 38]}对一维固结分析得出的结论相同.

4.3 压缩指数C_c的影响

为分析压缩指数C_c对固结过程的影响，图 6给出了次固结系数C_α=0和C_α=0.04两种情况下，C_c分别取0.6，0.8和1.0时，两个特殊时刻孔压u的分布曲线.可以看出，当不考虑土体黏滞性时，在整个固结过程中，压缩指数C_c越小(即土体的压缩性越低)，孔压消散就越快.而当考虑黏滞性(C_α=0.04)时，在固结的前期(t=1 000 min)，远离排水面的大部分范围内，压缩指数C_c越小(即土体的压缩性越低)，孔压值反而越大，即类曼德尔效应越明显.例如，从图 6(a)中C_α=0.04，t=1 000 min时的孔压分布曲线可以得出，在z=5 m，r=0.9 m处，与C_c=0.6，0.8和1.0对应的孔压值分别为164，159和155 kPa，均明显大于初始孔压值，而C_α=0时，与不同压缩指数对应的同一时刻同一位置处的孔压值均约为150 kPa.这与胡晶等^[20]分析压缩指数对土体一维固结影响时的结果类似.

图 7反映了压缩指数C_c对地基平均固结度U_p随时间变化的影响.可以看出，当不考虑土体黏滞性(C_α=0)时，每一时刻的平均固结度都随着压缩指数C_c的增大而减小，即压缩性越低，地基中孔压的整体消散越快.而当考虑黏滞效应(C_α=0.04)时，几条曲线是相互交叉的.即在固结前期(对应图 7中交点的左边)，固结度可能为负值，且压缩指数越大，地基的平均固结度越大.结合图 6可以看出，当土体压缩性较低时，更容易在地基中较大范围内出现类曼德尔效应；而在固结中后期(对应图 7中交点的右边)，平均固结度随压缩指数变化的趋势与C_α=0时相同.同时，从图 7也可以看出相同压缩指数时，考虑土体黏滞性后达到同一固结度需要的时间更长.

4.4 回弹指数C_s的影响

为考察回弹指数C_s对饱和软黏土地基流变固结过程的影响，图 8给出了z=5 m或r=0.5 m处不同回弹指数C_s时的孔压分布曲线.可以看出，虽然回弹指数C_s的变化对初期孔压升高有一定的影响，但与次固结系数C_α、初始超固结参数R₀和压缩指数C_c相比并不显著.例如，图 8(b)中，t=1 000 min，取C_s =0.06，0.12，0.30时，对应z=8 m，r=0.5 m处的孔压值分别为160，161，164 kPa，其差别并不大.

图 9给出了不同回弹指数C_s时的平均固结度随时间的变化曲线.很明显，其变化趋势与图 5类似，即在流变固结初期，随着回弹指数C_s的增大，平均固结度减小的现象也比较明显，同时在固结中后期，达到相同固结度所需时间随着回弹指数的增大也会相应增加.这与胡晶等^{[20, 38]}在一维固结分析时得出的结论类似，即回弹指数越大，土体可恢复的弹性变形部分也越大，这往往会使得地基固结时间延长.

对比图 3，5，7和9可以看出，考虑时间效应的UH本构模型中各参数对固结过程的影响程度有所不同.很明显，在本文讨论的参数范围内，相比次固结系数C_α，初始超固结参数R₀和压缩指数C_c，回弹指数C_s对固结度的影响相对较弱.

4.5 渗透指数C_k的影响

为探讨砂井地基流变固结过程中渗透指数C_k的影响，图 10给出了不同C_k时孔压的分布曲线.可以看出，t=1 000 min时，对应C_k=0.48，0.60，0.80时的几条曲线差别不大，特别在图 10(b)中几乎重叠在一起，在z=1 m，r=0.5 m处，对应不同C_k的孔压值分别为160，161，161 kPa.可见，孔压值均超过了初始孔压，但与前面讨论的UH模型参数C_α、R₀和C_c相比，渗透指数C_k的变化对固结初期类曼德尔效应的影响相对较弱.

图 11反映了渗透指数C_k对砂井地基平均固结度的影响.可以看出，在固结中后期，随着渗透指数C_k的减小，同一时刻的平均固结度U_p也减小.实际上，渗透指数C_k反映了渗透系数K对固结过程中孔隙比变化的敏感程度.在固结过程中，随着有效应力增大，土体变得密实，孔隙比减小，渗透系数也减小.而且，若减小渗透指数C_k，根据式(12)可知，渗透系数减小得更多，从而导致整体孔压消散更为缓慢.例如，当C_k取0.48，0.60，0.80时，对应平均固结度达到90%所需时间分别为189 700，292 600，460 400 min.这一点与李刚等^[36]得出的结论一致，他认为随着固结的进行，整个土层有效应力增加，导致土体渗透系数降低，渗流量减少，引起整体孔压消散滞后，地基达到相同固结度的时间也增加.可见，在分析饱和软黏土的固结进程时，应当考虑渗透系数随固结降低所产生的延滞影响.

4.6 Hansbo渗流参数m和i₁的影响

为分析Hansbo渗流模型对固结的影响，图 12，13分别给出了不同m或i₁时的孔压分布曲线.二图表明，Hansbo渗流对固结初期的类曼德尔效应有影响，但并不明显.例如，图 12(b)中，t=1 000 min时，z=1 m，r=0.5 m处，与m=1.0，1.2，1.5和1.8对应的孔压值分别为160，160，161和161 kPa.刘忠玉等^[30]曾对考虑Hansbo渗流但不计黏滞性的砂井地基固结进程做过分析，并没有发现类曼德尔效应.因此，结合前述分析主、次固结耦合变形机制对地基固结过程的影响规律可知，类曼德尔效应主要是由土体的黏滞效应或主次固结耦合机制引起的.

图 12，13表明，在固结后期(t=100 000 min)时，两种渗流模式计算的孔压值差异比较明显，且m或i₁值越大，同一位置处的孔压越大，即孔压消散越慢.例如，图 13(a)中，z=5 m，r=0.9 m，在Darcy渗流模式下(i₁=0时)算得的孔压值为33 kPa，当i₁取5，10和15时，对应Hansbo渗流模式下计算的结果分别为47，58和67 kPa，大约分别是Darcy渗流下孔压的1.42，1.76和2.03倍，可见，Hansbo渗流延缓了地基中孔压的消散.究其原因是随着固结时间的延长，较大范围土体内的孔压有所消散，水力梯度值也变得较小，当小于直线渗流起始水力梯度i₁时，渗流是很缓慢的，且随着Hansbo渗流参数的增大，渗流速度更小，孔隙水更不易排出，所以，孔压消散也更为缓慢.

同时，与前述UH模型参数有关分析图类似的是，图 12，13中孔压沿竖向分布曲线表明，在超过一定深度后，孔压沿竖向基本不再变化，且Hansbo渗流参数越大，孔压不变化的范围也越大.这说明超过某个深度后，水力梯度的竖向分量很小，渗流速度的竖向分量也很小，因此，对于厚度较大的砂井地基，考虑土体黏滞性和非Darcy渗流时径向渗流为主要的排水途径.

图 14给出m或i₁取不同数值时，平均固结度随时间的变化曲线.可以看出，在流变固结中后期，按Hansbo渗流模式计算的平均固结度曲线均在按Darcy渗流模式计算的上方，即地基的固结比基于Darcy定律计算的要慢，并且这一现象随着Hansbo渗流参数的增大而变得更为明显.这也表明在实际工程的固结分析中考虑土体渗流非Darcy特性的必要性，否则将高估砂井地基固结的程度.这与文献[27-30]的结论类似.

5 结论

1) 固结初期出现的孔压异常升高现象(类曼德尔效应)主要由土体的黏滞效应或主次固结耦合机制引起.增大次固结系数C_α和初始超固结参数R₀、或减小压缩指数C_c都会使该现象变得更加明显，但其他参数对此现象的影响较弱.

2) UH模型参数的变化对砂井地基固结过程均有不同程度的影响.随着次固结指数C_α、初始超固结参数R₀和压缩指数C_c等参数的增大，砂井地基的固结变慢.与这些模型参数相比，回弹指数C_s变化引起的固结度变化较小.

3) 渗透指数C_k对砂井地基固结进程的影响主要体现在中后期，且该参数越小，孔压消散就越慢，从而砂井地基固结也就越慢.

4) 与Darcy渗流相比，Hansbo渗流使得砂井地基固结变慢.同时考虑Hansbo渗流和土体黏滞性时，砂井地基中的排水更是以径向渗流为主.因此，应该重视渗流非达西特性对砂井地基固结进程的影响.