0 引言

在二维或三维固结初期，土体内部的孔压可能不降反升，该现象首先由Mandel在1953年发现，后来Cryer^[1]在研究土球受均布压力径向向外排水时也发现此现象，故而称为Mandel-Cryer效应.这一效应并不能用Terzaghi一维或拟多维固结理论解释，而Biot从较严格的固结机理出发推导出的多维固结理论却能很好地解释这一效应，因此，学者们对此展开了大量研究.文献[2]通过试验再次证实了Mandel-Cryer效应，文献[3-6]通过多种方法求解了圆球土样的Biot固结方程，并对Mandel-Cryer效应的影响因素进行了分析，研究结果表明，减小泊松比、渗透系数、杨氏模量等参数或远离排水边界都会延迟Mandel-Cryer效应.但是，这些学者对Biot固结理论的研究都基于Darcy渗流，鲜有学者在研究Biot固结理论时考虑渗流的非线性.早在1925年，Terzaghi就指出大塑性黏性土中的渗流可能偏离Darcy定律.文献[7-12]的渗透试验证实了这一论断，并分别提出不同形式的非Darcy渗流数学拟合公式，其中应用较广的当属Hansbo渗流^[7].

目前考虑Hansbo渗流对固结性状影响的研究多集中于一维问题，例如文献[13-21]先后引入Hansbo渗流或其简化形式、指数渗流等修正了Terzaghi一维固结理论，成功解释了后者所不能解释的一些试验现象.随着研究的深入，已有学者开始讨论非Darcy渗流对砂井地基固结的影响^[22-23]，但在Biot固结理论框架内讨论非Darcy渗流影响的报道还很少.因此，为深入探讨非Darcy渗流对固结过程的影响，引入Hansbo渗流模型代替Darcy定律，基于Biot固结理论重新推导了圆球土样的固结方程，并探讨了Hansbo渗流参数对固结进程的影响.

1 圆球土样Biot固结方程的推导及求解 1.1 控制方程

设半径为a的饱和圆球土样为均质且各向同性的线弹性体，其弹性模量为E，泊松比为μ.假定土颗粒和孔隙水不可压缩，土球内部没有初始孔压，且土样外表面完全透水.现于其外表面沿径向作用均布压力q，这样土样内的渗流和变形都是球对称的.在外荷载作用下，t时刻距球心r处的径向位移为u，体应变为θ，孔压为p.并且规定径向位移u沿半径方向向外为正，反之为负；体应变θ以收缩为正，反之为负；压力q以压为正，以拉为负.则其平衡方程、物理方程和几何方程分别为

$ \frac{{\partial {\sigma _{\rm{r}}}}}{{\partial r}} + \frac{2}{r}\left( {{\sigma _{\rm{r}}} - {\delta _{\rm{ \mathit{ τ} }}}} \right) + \frac{{\partial p}}{{\partial r}} = 0, $

(1)

$ \left\{ \begin{array}{l} {\varepsilon _{\rm{r}}} = \frac{1}{{2G}}\left[ {{\sigma _{\rm{r}}} - \frac{\mu }{{1 + \mu }}\left( {{\sigma _{\rm{r}}} + 2{\sigma _{\rm{ \mathit{ τ} }}}} \right)} \right], \\ {\varepsilon _{\rm{ \mathit{ τ} }}} = \frac{1}{{2G}}\left[ {{\sigma _{\rm{ \mathit{ τ} }}} - \frac{\mu }{{1 + \mu }}\left( {{\sigma _{\rm{r}}} + 2{\sigma _{\rm{ \mathit{ τ} }}}} \right)} \right], \end{array} \right. $

(2)

$ {\varepsilon _{\rm{r}}} = \frac{{\partial u}}{{\partial r}};{\varepsilon _{\rm{ \mathit{ τ} }}} = \frac{u}{r}, $

(3)

式中：G为剪切模量；ε_r、ε_τ分别为径向和环向应变；σ_r、σ_τ分别为径向和环向应力.

联立方程组(1)~(3)，可得用位移和孔压表示的平衡微分方程为

$ \frac{{{\partial ^2}u}}{{\partial {r^2}}} + \frac{1}{r}\frac{{\partial u}}{{\partial r}} - \frac{2}{{{r^2}}}u - \frac{{\left( {1 + \mu } \right)\left( {1 - 2\mu } \right)}}{{E\left( {1 - \mu } \right)}}\frac{{\partial p}}{{\partial r}} = 0. $

(4)

体应变θ和径向位移u的关系为

$ \theta = - \left( {\frac{{\partial u}}{{\partial r}} + \frac{2}{r}u} \right), $

(5)

$ \frac{{\partial \theta }}{{\partial t}} = - \frac{1}{{{r^2}}}\frac{\partial }{{\partial r}}\left( {{r^2}\frac{{\partial u}}{{\partial t}}} \right). $

(6)

设t时刻距球心r处的径向水力梯度为i，径向渗流速度为v，且二者之间的关系可以采用Hansbo渗流公式描述为

$ v = \left\{ {\begin{array}{*{20}{l}} {c{i^m}, }&{i \le {i_1}, }\\ {k\left( {i - {i_0}} \right), }&{i > {i_1}, } \end{array}} \right. $

(7)

式中：m为试验确定的常数；c、k分别为指数和线性关系式中的渗透系数；i₀和i₁分别为直线渗流的起始水力梯度和计算起始水力梯度，且i₀=i₁(m-1)/m，c=k/(mi₁^m-1).如令m=1，i₀=0，则式(7)可退化为Darcy定律表达式.为便于应用，仿照Darcy定律，将式(7)改写为

$ v = {k_1}i = - {k_1}\frac{{\partial p}}{{{\gamma _{\rm{w}}}\partial r}} = - k'\frac{{\partial p}}{{\partial r}}, $

(8)

式中：$ i = - \frac{1}{{{\gamma _{\rm{w}}}}}\frac{{\partial p}}{{\partial r}}$，其中γ_w为水的重度；k₁=k'γ_w，可称为视渗透系数，它是水力梯度i的函数，即

$ {k_1} = \left\{ {\begin{array}{*{20}{l}} {c{{\left| i \right|}^{m - 1}}, }&{\left| i \right| \le {i_1}, }\\ {k\left( {1 - {i_0}/\left| i \right|} \right), }&{\left| i \right| > {i_1}.} \end{array}} \right. $

(9)

考虑渗流连续性，即土骨架的体积变化=孔隙体积变化=流入流出单元体的水量差，则

$ \frac{{\partial \theta }}{{\partial t}} = - \frac{1}{{{r^2}}}\frac{\partial }{{\partial r}}\left( {{r^2}v} \right). $

(10)

将式(8)代入式(10)，可得

$ \frac{{\partial \theta }}{{\partial t}} = - \frac{1}{{{r^2}}}\frac{\partial }{{\partial r}}\left( {k'{r^2}\frac{{\partial p}}{{\partial r}}} \right). $

(11)

联立式(6)、(11)，可得

$ \frac{{\partial u}}{{\partial t}} = k'\frac{{\partial p}}{{\partial r}}. $

(12)

将式(12)代入式(4)，可得

$ \frac{{\partial u}}{{\partial t}} = {c_{\rm{v}}}k'\frac{{{\gamma _{\rm{w}}}}}{k}\left( {\frac{{{\partial ^2}u}}{{\partial {r^2}}} + \frac{2}{r}\frac{{\partial u}}{{\partial r}} - \frac{2}{{{r^2}}}u} \right). $

(13)

式中：C_v=kE(1-μ)/[γ_w(1+μ)(1-2μ)].方程组(12)、(13)的边值条件及初始条件分别为

$ - \frac{E}{{\left( {1 + \mu } \right)\left( {1 - 2\mu } \right)}}\left[ {\left( {1 - \mu } \right)\frac{{\partial u}}{{\partial r}} + 2\mu \frac{u}{r}} \right]{|_{r = a}} = q, $

(14)

$ u\left( {0, t} \right) = 0;u\left( {r, 0} \right) = 0, $

(15)

$ p\left( {0, t} \right) < \infty ;p\left( {a, t} \right) = 0. $

(16)

为了便于讨论引入下列无量纲化参数：

$ R = \frac{r}{a};U = \frac{u}{a};T = \frac{{{C_{\rm{v}}}t}}{{{a^2}}};I = \frac{{i{\gamma _{\rm{w}}}a}}{q};P = \frac{p}{q};{I_1} = \frac{{{i_1}{\gamma _{\rm{w}}}a}}{q}. $

(17)

则式(12)~(16)可以化为

$ \frac{{\partial P}}{{\partial R}} = \frac{{E\left( {1 - \mu } \right)}}{{\psi q\left( {1 + \mu } \right)\left( {1 - 2\mu } \right)}}\frac{{\partial U}}{{\partial T}}, $

(18)

$ {R^2}\frac{{\partial U}}{{\partial T}} = \psi \left( {{R^2}\frac{{{\partial ^2}U}}{{\partial {R^2}}} + 2R\frac{{\partial U}}{{\partial R}} - 2U} \right), $

(19)

$ - \frac{{E/q}}{{\left( {1 + \mu } \right)\left( {1 - 2\mu } \right)}}\left[ {\left( {1 - \mu } \right)\frac{{\partial U}}{{\partial R}} + 2\mu \frac{U}{R}} \right]{|_{R = 1}} = 1, $

(20)

$ U\left( {0, T} \right) = 0;U\left( {R, 0} \right) = 0, $

(21)

$ P\left( {0, T} \right) < \infty ;P\left( {1, T} \right) = 0, $

(22)

式中: $ \psi = \left\{ {\begin{array}{*{20}{l}} {{I^{m - 1}}/(mI_1^{m - 1}), }&{|I| \le {I_1}}\\ {1 - {I_1}\left( {m - 1} \right)/\left( {|I|m} \right), }&{|I| > {I_1}, } \end{array}} \right.I = - \frac{{\partial P}}{{\partial R}}.$

1.2 有限差分法求解

上述微分方程组为变系数抛物型方程，取得其解析解较为困难.因此，拟采用Crank-Nicolson有限差分法求其数值解.即在0≤R≤1范围内以ΔR为步长，将土球由球心沿半径方向离散为N层，而对时间则按步长ΔT进行离散，则式(19)~(21)可离散为

$ {A_b}U_{b - 1}^{j + 1} + {B_b}U_b^{j + 1} + {C_b}U_{b + 1}^{j + 1} = - {A_b}U_{b - 1}^j + \left( {2 - {B_b}} \right)U_b^j - {C_b} \cdot U_{b + 1}^j, $

(23)

$ \left( {{A_N} + {C_N}} \right)U_{N - 1}^{j + 1} + \left( {{B_N} - \eta {C_N}} \right)U_N^{j + 1} = \left( { - C - {A_N}} \right)\\U_{N - 1}^j + \left( {2 - {B_N} + \eta {C_N}} \right)U_N^j + {C_N}Q, $

(24)

$ U_0^j = 0;U_b^0 = 0, $

(25)

式中: $ {A_b} = \frac{{\left( {1 - b} \right)\psi \Delta T}}{{2b\;\Delta {R^2}}};{C_b} = - \frac{{\left( {1 + b} \right)\psi \Delta T}}{{2b\;\Delta {R^2}}};$$ {B_b} = \frac{{1 + (1 + {b^2})\psi \Delta T}}{{{{\left( {b\Delta R} \right)}^2}}};\eta = \frac{{4\mu \cdot \Delta R}}{{1 - \mu }};$$ Q = \frac{{4\Delta R\left( {1 + \mu } \right)\left( {1 - 2\mu } \right)}}{{E\left( {1 - \mu } \right)}}$，下标和上标分别表示空间和时间节点.

这样由式(23)~(25)可求得位移，然后再利用式(18)和式(22)求孔压，可离散得

$ P_b^j = P_{b + 1}^j - \frac{{E\left( {1 - \mu } \right)}}{{\psi q\left( {1 + \mu } \right)\left( {1 - 2\mu } \right)}}\frac{{U_b^{j + 1} - U_b^j}}{{\Delta T}} \cdot \Delta R, $

(26)

$ P_0^j < \infty ;P_N^j = 0. $

(27)

为分析圆球土样整体孔压消散情况，引入按孔压定义的平均固结度为

$ {U_t} = 1 - \frac{{4\pi \int_0^a {p{r^2}{\rm{d}}r} }}{{{p_0} \cdot 4\pi {a^3}/3}} = 1 - 3\int_0^1 {P{R^2}{\rm{d}}R} = 1 - 3\\ \sum\limits_{b = 0}^{N - 1} {\int_{{R_b}}^{{R_{b + 1}}} {P{R^2}{\rm{d}}R} } . $

(28)

假定孔压P在区间[R_b, R_b+1]内为线性分布，则在j时刻的平均固结度为

$ {U_t} = 1 - \sum\limits_{b = 0}^{N - 1} {\left[ {\frac{{3\left( {P_{b + 1}^j - P_b^j} \right)}}{{4\Delta R}}\left( {R_{b + 1}^4 - R_b^4} \right) + \frac{{P_b^j{R_{b + 1}} - P_{b + 1}^j{R_b}}}{{\Delta R}}\left( {R_{b + 1}^3 - R_b^3} \right)} \right].} $

(29)

2 算法验证

用Matlab软件编制程序，为了验证本文算法及所编程序的正确性，对m=1.0或I₁=0即Darcy渗流时的圆球土样固结过程进行了数值模拟.模拟时，取ΔR=0.01，ΔT=0.000 01，得到的土球外边界处的位移随时间的变化曲线示于图 1.Cryer^[1]曾给出了其解析解，为便于对比，将Darcy渗流时圆球土样的边界位移也示于图 1中.可以看出，本文有限差分解与Cryer的解析解基本一致，表明本文算法及所编程序是有效的.

3 参数分析

文献[3-6]讨论了弹性模量和泊松比对土体固结过程的影响，本文仅讨论Hansbo渗流参数的影响，分析中取E=5 MPa，μ=0.35，q=1 MPa.

3.1 Hansbo渗流对孔压的影响

为考察Hansbo渗流对Mandel-Cryer效应的影响，图 2给出了m不同时孔压P随时间T的变化曲线.很明显，在加载的初期，孔压都是升高的，达到峰值后再开始消散，即出现了Mandel-Cryer效应.在Darcy渗流条件下，在T=0.017时孔压P达到其峰值1.12.而在Hansbo渗流条件下，当m=1.2、1.5和1.8时，孔压峰值分别为1.21、1.31和1.38，对应的时间分别为0.022、0.027和0.031.因此，Hansbo渗流增强了Mandel-Cryer效应，m越大，该效应越明显，即孔压峰值就越大，且达到峰值所需要的时间也越长.

图 3给出了I₁不同时孔压P随时间T的变化曲线.I₁对Mandel-Cryer效应的影响规律和m值相似，即I₁越大，孔压峰值就越大，且达到峰值所需要的时间也越长.

另外，图 2和图 3还表明，在孔压达到峰值后转入消散的固结中后期，Hansbo渗流将使圆球土样内的孔压消散滞后，并且其参数越大，这种滞后越明显.

3.2 Hansbo渗流对固结度的影响

平均固结度U_t反映了圆球土样的固结程度，也就是圆球土样在相应时刻的孔压整体消散程度.为考察Hansbo渗流对平均固结度的影响，图 4和图 5分别给出了不同m或I₁时平均固结度U_t随时间T的变化曲线.很明显，基于Hansbo渗流计算的固结度曲线都在基于Darcy渗流计算的曲线下方，并且Hansbo渗流参数I₁或m越大，相应的曲线偏离基于Darcy渗流计算的曲线越远，即固结度U_t随着Hansbo渗流参数m和I₁的增大而减小，也就是说圆球土样内孔压的整体消散随Hansbo渗流参数m和I₁的增大而减缓，这与文献[16]中Hansbo渗流对一维固结影响的结论相似.

3.3 Hansbo渗流对边界位移的影响

记U(1, ∞)为圆球土样边界处的最终位移，则U(1, T)/U(1, ∞)为边界处T时刻的无量纲位移，类似于一维固结时按变形定义的固结度U_s.按照文献[16, 18, 21]的分析，对于一维弹性饱和地基来说，不论Darcy渗流还是Hansbo渗流，按变形定义的固结度U_s和按孔压定义的固结度U_t是一致的，但对于圆球土样却不是这样的.图 6和图 7分别给出了不同Hansbo渗流参数下U(1, T)/U(1, ∞)随时间T的变化曲线.很明显，这些曲线与图 4和图 5所示的U_t随T的变化曲线尽管趋势一样，但数值相差却很大；尽管Hansbo渗流参数越大，同一时刻的边界位移会越小，但从图 6和图 7的几条曲线几乎重合来看，其数值相差无几，所以Hansbo渗流对边界位移的影响很小.

4 结论

基于Hansbo渗流模型推导了圆球土样Biot固结方程，给出了方程的Crank-Nicolson有限差分格式，并用Matlab软件编写了有关程序.计算结果表明：和Darcy渗流相比，在固结初期，Hansbo渗流将增强Mandel-Cryer效应，增大孔隙水压力的峰值，并延长孔隙水压力达到峰值的时间；在固结中后期，Hansbo渗流将使圆球土样内孔压消散明显滞后.随着Hansbo渗流参数的增大，土样平均固结度减小，即土样固结速率减慢，但Hansbo渗流对边界位移的影响很小.