在描述多孔介质渗流问题时，多数采用Darcy定律。该定律主要用于描述层流问题，当雷诺数较大时，不再适用。黏性土试验结果表明，当水力梯度较小时，土体中的渗流存在偏离Darcy定律的现象，称之为非Darcy渗流^[1]。TEH等将非Darcy渗流方程引入到固结方程中，分析了非Darcy渗流对固结性状的影响^[2-3]。李传勋等^[4-10]基于Hansbo非Darcy渗流模型，考虑变荷载与应力历史的影响，建立了以超孔压为变量的软土固结模型，分析了Hansbo模型参数对地基固结速率的影响。崔莉红等采用渗流-固结-溶质运移耦合模型研究了华北平原弱透水层非Darcy渗流固结过程，并对咸水下移的影响进行了分析^[11]，指出非Darcy渗流在增加黏性土固结时间的同时，也减缓了咸水下移的速度。李菲菲等基于竖井地基固结理论，假设孔隙水服从指数渗流模式，建立了新的竖井地基固结模型，分析了渗流模型对地基固结的影响^[12]。LIU等结合Hansbo非Darcy渗流与Terzaghi固结理论，考虑施工变荷载的影响，用有限体积法重新推导了Terzaghi固结方程，分析了非Darcy渗流与施工速度对黏土固结过程的影响^[13]。然而，土体在固结过程中，体积不断减小，孔隙比随之降低，由此导致渗透系数降低，土体排水能力减弱，进而引起渗流速度下降，对土体固结过程产生影响。因此，在固结分析时，应考虑渗透系数随固结变化所产生的影响。本文考虑Hansbo非Darcy渗流与渗透系数变化的耦合影响，建立了二维固结控制方程。通过数值计算分析了非线性渗流对土体沉降与超孔压的影响。

1 渗透系数变化与非线性渗流下的固结控制方程 1.1 考虑渗透系数变化的非线性渗流问题描述

Hansbo^[14-15]基于试验给出渗流速度v与水力梯度i之间的关系式为

$ v = \left\{ {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {c{i^m}}\\ {k\left( {i - {i_0}} \right)} \end{array}}&{\begin{array}{*{20}{c}} {i < {i_L}}\\ {i \ge {i_L}} \end{array}} \end{array}} \right. $

(1)

式中：c、k分别为曲线段与直线段的渗透系数，cm/s；m为试验常数，取值范围为1.2~1.5；i₀为计算起始水力梯度；i_L为线性渗流起始水力梯度，取值范围为4~10。

正常固结土的孔隙比e与平均有效应力p′的关系可表示为

$ e = {e_0} - {C_c}\lg \left( {\frac{{p'}}{{{{p'}_0}}}} \right) $

(2)

式中：e₀为初始孔隙比；C_c为压缩指数；p′₀为初始平均有效应力，kPa。

Taylor^[16]基于试验结果，给出渗透系数k与孔隙比e的关系为

$ \lg k = \lg {k_0} - \frac{{{e_0} - e}}{{{C_k}}} $

(3)

式中：k₀为初始渗透系数，cm/s；e₀为初始孔隙比；C_k为渗透系数变化指数。将式(2) 代入式(3) 得

$ k = {k_0}{\left[ {\frac{{p'}}{{{{p'}_0}}}} \right]^{ - {\alpha _c}}} $

(4)

式中：α_c=C_c/C_k，取值范围为0.5~1.5。

利用式(1) 中v及其导数在分界点i_L处的连续条件，可求得参数c=k/(mi_L^m-1)，i₀=(m-1)i_L/m。

将式(4) 代入式(1) 得到考虑渗透系数变化的非Darcy渗流方程为

$ \left\{ \begin{array}{l} v = {k_0}{\left( {\frac{{p'}}{{{{p'}_0}}}} \right)^{ - {\alpha _c}}}{i^m}/\left( {mi_L^{m - 1}} \right),\;\;\;\;\;\;\;\left| i \right| < {i_L}\\ v = {k_0}{\left( {\frac{{p'}}{{{{p'}_0}}}} \right)^{ - {\alpha _c}}}\left( {i - \frac{{m - 1}}{m}{i_L}} \right),\;\;\;\;\;\left| i \right| \ge {i_L} \end{array} \right. $

(5)

1.2 平衡方程

在静荷载下(惯性力为零)，平面问题的平衡方程、有效应力原理、本构方程及大变形几何方程可以表示如下

$ \left\{ \begin{array}{l} {\mathit{\boldsymbol{L}}^{\rm{T}}}\mathit{\boldsymbol{S}} + {f_0} = 0\\ \mathit{\boldsymbol{S}} = \mathit{\boldsymbol{S' + T}}{\mathit{\boldsymbol{p}}_{ow}}\\ \mathit{\boldsymbol{S' = DE}}\\ \mathit{\boldsymbol{E}} = {\mathit{\boldsymbol{E}}_L} + {\mathit{\boldsymbol{E}}_N} \end{array} \right. $

(6)

式中：L为微分算子矩阵，S为Kirchhoff应力矩阵，f₀为体力荷载矩阵，S′为有效应力矩阵，T为选项矢量，p_ow为孔隙水压力，D为本构矩阵，E为Green应变矩阵，E_L为几何线性应变矩阵，E_N为几何非线性应变矩阵，各矩阵的具体表达形式见文献[17]。

1.3 连续方程

在饱和土体中，单元体内水量的变化率在数值上等于土体积的变化率^[18]，即

$ \frac{{\partial {v_x}}}{{\partial x}} + \frac{{\partial {v_y}}}{{\partial y}} = \frac{\partial }{{\partial t}}\left( {\frac{{\partial u}}{{\partial x}} + \frac{{\partial v}}{{\partial y}}} \right) $

(7)

式中：v_x、v_y分别为x与y方向的渗流速度，u、v分别为x与y方向的位移。式(5) 可表示为

$ \left\{ \begin{array}{l} {v_\psi } = {S_\psi }{k_{\psi 0}}{\left( {\frac{{p'}}{{{{p'}_0}}}} \right)^{ - {\alpha _c}}}{\left| {{i_\psi }} \right|^m}/\left( {mi_L^{m - 1}} \right),\;\;\;\;\;\;\left| {{i_\psi }} \right| < {i_L}\\ {v_\psi } = {S_\psi }{k_{\psi 0}}{\left( {\frac{{p'}}{{{{p'}_0}}}} \right)^{ - {\alpha _c}}}\left( {\left| {{i_\psi }} \right| - \frac{{m - 1}}{m}{i_L}} \right),\;\;\;\;\left| {{i_\psi }} \right| \ge {i_L} \end{array} \right. $

(8)

式中：下标ψ可取x或y；当水力梯度方向与微元体正方向相同时，S_ψ为1；方向相反时，S_ψ为-1。将式(8) 对ψ求导得

$ \left\{ \begin{array}{l} \frac{{\partial {v_\psi }}}{{\partial \psi }} = \left( {A + B} \right)\frac{{{k_{{\psi _0}}}}}{{{\gamma _w}}}\frac{{{\partial ^2}{p_{{\rm{ow}}}}}}{{\partial {\psi ^2}}},\;\;\;\;\;\;\;\;\;\;\;\;\left| {{i_\psi }} \right| < {i_L}\\ \frac{{\partial {v_\psi }}}{{\partial \psi }} = \left( {C - D + E} \right)\frac{{{k_{{\psi _0}}}}}{{{\gamma _w}}}\frac{{{\partial ^2}{p_{{\rm{ow}}}}}}{{\partial {\psi ^2}}},\;\;\;\;\;\;\left| {{i_\psi }} \right| \ge {i_L} \end{array} \right. $

(9)

式中：

$ A = \frac{{{\alpha _c}}}{{mi_L^{m - 1}\gamma _w^{m - 1}p'\frac{{{\partial ^2}{p_{{\rm{ow}}}}}}{{\partial {\psi ^2}}}}}{\left( {\frac{{p'}}{{{{p'}_0}}}} \right)^{ - {\alpha _c}}}{\left| {\frac{{\partial {p_{{\rm{ow}}}}}}{{\partial \psi }}} \right|^{m + 1}} $

$ B = \frac{1}{{i_L^{m - 1}\gamma _w^{m - 1}}}{\left( {\frac{{p'}}{{{{p'}_0}}}} \right)^{ - {\alpha _c}}}{\left| {\frac{{\partial {p_{{\rm{ow}}}}}}{{\partial \psi }}} \right|^{m - 1}} $

$ C = \frac{{{\alpha _c}}}{{p'\frac{{{\partial ^2}{p_{{\rm{ow}}}}}}{{\partial {\psi ^2}}}}}{\left( {\frac{{p'}}{{{{p'}_0}}}} \right)^{ - {\alpha _c}}}{\left| {\frac{{\partial {p_{{\rm{ow}}}}}}{{\partial \psi }}} \right|^2} $

$ D = \frac{{\left( {m - 1} \right){i_L}{\gamma _w}{\alpha _c}}}{{mp'\frac{{{\partial ^2}{p_{{\rm{ow}}}}}}{{\partial {\psi ^2}}}}}{\left( {\frac{{p'}}{{{{p'}_0}}}} \right)^{ - {\alpha _c}}}\left| {\frac{{\partial {p_{{\rm{ow}}}}}}{{\partial \psi }}} \right| $

$ E = {\left( {\frac{{p'}}{{{{p'}_0}}}} \right)^{ - {\alpha _c}}} $

式(9) 也可表示为

$ \frac{{\partial {v_\psi }}}{{\partial \psi }} = {H_\psi }\frac{{{k_{{\psi _0}}}}}{{{\gamma _w}}}\frac{{{\partial ^2}{p_{{\rm{ow}}}}}}{{\partial {\psi ^2}}} $

(10)

其中，

$ {H_\psi } = \left\{ \begin{array}{l} A + B,\;\;\;\;\;\;\;\;\;\;\left| {{i_\psi }} \right| < {i_L}\\ C - D + E,\;\;\;\;\left| {{i_\psi }} \right| \ge {i_L} \end{array} \right. $

(11)

将式(10) 代入式(7) 中得到

$ {H_x}\frac{{{k_{{\psi _0}}}}}{{{\gamma _w}}}\frac{{{\partial ^2}{p_{{\rm{ow}}}}}}{{\partial {\psi ^2}}} + {H_y}\frac{{{k_{{y_0}}}}}{{{\gamma _w}}}\frac{{{\partial ^2}{p_{{\rm{ow}}}}}}{{\partial {y^2}}} = \frac{\partial }{{\partial t}}\left( {\frac{{\partial u}}{{\partial x}} + \frac{{\partial v}}{{\partial y}}} \right) $

(12)

写成矩阵形式为

$ {\mathit{\boldsymbol{T}}^{\rm{T}}}\mathit{\boldsymbol{LH}}{\mathit{\boldsymbol{k}}_0}{\mathit{\boldsymbol{L}}^{\rm{T}}}\mathit{\boldsymbol{T}}{p_{ow}} - \frac{\partial }{{\partial t}}{\mathit{\boldsymbol{T}}^{\rm{T}}}\mathit{\boldsymbol{E}} = 0 $

(13)

式中：H为渗流控制矩阵，H=$\left[ {\begin{array}{*{20}{l}} {{H_x}}&0\\ 0&{{H_y}} \end{array}} \right]$；k₀为初始渗透系数矩阵，k₀=$\frac{1}{{{\gamma _w}}}\left[ {\begin{array}{*{20}{l}} {{k_{x0}}}&0\\ 0&{{k_{y0}}} \end{array}} \right]$。

2 渗透系数变化与非线性渗流下的固结有限元方程 2.1 位移与孔压模式选取

采用八结点四边形等参单元进行离散，如图 1所示。其中，位移结点为8个，孔压结点为4个。单元的位移函数与孔压函数为

$ \left\{ \begin{array}{l} \tilde u = \sum\limits_{i = 1}^8 {{N_i}{u_i}} \\ \tilde v = \sum\limits_{i = 1}^8 {{N_i}{v_i}} \\ \tilde p = \sum\limits_{j = 1}^8 {{N_j}{p_j}} \end{array} \right. $

(14)

式中：$\tilde u,\tilde v,\tilde p$为单元内任一点的位移与孔压近似值；u_i、v_i、p_j为结点值；N_i与N_j分别为位移与孔压形函数，可表示为

$ \left\{ \begin{array}{l} {N_{i\left( {1,3,5,7} \right)}} = \frac{1}{4}\left( {1 + {\xi _0}} \right)\left( {1 + {\eta _0}} \right)\left( {{\xi _0} + {\eta _0} - 1} \right)\\ {N_{i\left( {2,6} \right)}} = \frac{1}{2}\left( {1 - {\xi ^2}} \right)\left( {1 + {\eta _0}} \right)\\ {N_{i\left( {4,8} \right)}} = \frac{1}{2}\left( {1 - {\eta ^2}} \right)\left( {1 + {\xi _0}} \right) \end{array} \right. $

(15)

$ {N_{j\left( {1,2,3,4} \right)}} = \frac{1}{4}\left( {1 + {\xi _0}} \right)\left( {1 + {\eta _0}} \right) $

(16)

式中：ξ₀=ξξ_i，η₀=ηη_i。

2.2 平衡方程离散

将有效应力原理、本构方程代入到平衡方程中，结合几何方程得到离散后的平衡方程为

$ {\mathit{\boldsymbol{K}}_{DS}}\Delta {a_e} + {\mathit{\boldsymbol{K}}_C}\Delta p_{ow}^e = - \mathit{\boldsymbol{\bar R}} - {\mathit{\boldsymbol{R}}_{SP}} $

(17)

式中：K_DS为刚度矩阵，K_C为耦合矩阵，R为等效结点荷载矢量，R_SP为等效结点力矢量，Δa_e、Δp_ow^e分别为单元的结点位移与超孔压增量。具体计算式如下

$ {\mathit{\boldsymbol{K}}_{DS}} = {\mathit{\boldsymbol{K}}_D} + {\mathit{\boldsymbol{K}}_S},{\mathit{\boldsymbol{K}}_D} = \int_{{V_0}} {{\mathit{\boldsymbol{B}}^{\rm{T}}}\mathit{\boldsymbol{DB}}{\rm{d}}{V_0}} , $

$ {\mathit{\boldsymbol{K}}_S} = \int_{{V_0}} {{\mathit{\boldsymbol{G}}^{\rm{T}}}\mathit{\boldsymbol{MG}}{\rm{d}}{V_0}} ,{\mathit{\boldsymbol{K}}_C} = \int_{{V_0}} {{\mathit{\boldsymbol{B}}^{\rm{T}}}\mathit{\boldsymbol{T\bar N}}{\rm{d}}{V_0}} , $

$ \mathit{\boldsymbol{\bar R}} = \int_{{V_o}} {{\mathit{\boldsymbol{N}}^{\rm{T}}}{\mathit{\boldsymbol{f}}_0}{\rm{d}}{V_0}} = \int_{{A_o}} {{\mathit{\boldsymbol{N}}^{\rm{T}}}{\mathit{\boldsymbol{q}}_0}{\rm{d}}{A_0}} , $

$ {\mathit{\boldsymbol{P}}_{SP}} = \int_{{V_0}} {{\mathit{\boldsymbol{B}}^{\rm{T}}}\left( {\mathit{\boldsymbol{S' + T}}{p_{ow}}} \right){\rm{d}}{V_0}} . $

2.3 连续方程离散

采用Galerkin加权残值法与差分法分别对连续方程进行空间与时域离散，得到离散后的连续方程为

$ \mathit{\boldsymbol{K}}_C^{\rm{T}}\Delta {a_e} - \theta \Delta t{\mathit{\boldsymbol{K}}_{SP}}\Delta p_{ow}^e = {\mathit{\boldsymbol{R}}_{PQ}} $

(18)

式中：θ为差分系数，取值范围为0.5~1.0，若取θ=0.5，则属中心差分格式；K_SP为渗流矩阵；R_PQ为边界流量等效结点力矢量。具体计算式如下

$ \begin{array}{l} {\mathit{\boldsymbol{K}}_{SP}} = \int_{{V_0}} {\mathit{\boldsymbol{B}}_s^{\rm{T}}\mathit{\boldsymbol{Hk}}{\mathit{\boldsymbol{B}}_s}{\rm{d}}{V_0}} \\ {\mathit{\boldsymbol{R}}_{PQ}} = \Delta t\left( {{\mathit{\boldsymbol{K}}_{SP}}\Delta p_{ow}^e + {\mathit{\boldsymbol{R}}_Q}} \right) \end{array} $

3 非线性渗流的计算方法

由有限元方程的推导可知，在考虑非线性渗流时，渗流矩阵K_SP中增加了一项渗流控制矩阵H，计算时需要根据应力与孔压来更新控制矩阵H。为了计算K_SP，在计算的初始时步，设H为单位阵，在随后的计算中，通过前一时步得到的有效应力与结点孔压来计算H，由此得到K_SP。

为计算水力梯度i_ψ，在局部坐标系下，如图 1所示，设j时刻的结点孔压为p_i^j(i表示结点，j表示时刻)，则j+1时刻的水力梯度i^j+1可按下式计算：

$ \left\{ \begin{array}{l} i_x^{j + 1} = \frac{1}{{{\gamma _w}}}\frac{{\partial p_{ow}^j}}{{\partial x}} = \frac{{\left( {p_{1\# }^j + p_{4\# }^j - p_{2\# }^j - p_{3\# }^j} \right)}}{{4{\gamma _w}a}}\\ i_y^{j + 1} = \frac{1}{{{\gamma _w}}}\frac{{\partial p_{ow}^j}}{{\partial y}} = \frac{{\left( {p_{1\# }^j + p_{2\# }^j - p_{3\# }^j - p_{4\# }^j} \right)}}{{4{\gamma _w}b}} \end{array} \right. $

(19)

当水力梯度i_ψ确定后，$\frac{{\partial {p_{{\rm{ow}}}}}}{{\partial \psi }}$即可求得，而$\frac{{{\partial ^2}{p_{{\rm{ow}}}}}}{{\partial {\psi ^2}}}$可按下式计算：

$ \left\{ \begin{array}{l} \frac{{{\partial ^2}p_{{\rm{ow}}}^{j + 1}}}{{\partial {x^2}}} = \left( {\frac{{\partial p_{{\rm{ow}}}^{j + 1}}}{{\partial x}} - \frac{{\partial p_{{\rm{ow}}}^j}}{{\partial x}}} \right)/\Delta x\\ \frac{{{\partial ^2}p_{{\rm{ow}}}^{j + 1}}}{{\partial {y^2}}} = \left( {\frac{{\partial p_{{\rm{ow}}}^{j + 1}}}{{\partial y}} - \frac{{\partial p_{{\rm{ow}}}^j}}{{\partial y}}} \right)/\Delta y \end{array} \right. $

(20)

4 计算结果与讨论 4.1 数值解与近似解析解比较

针对前几节所建立的有限元方程，编制了二维有限元程序。为验证本文所建立的有限元方程及计算程序的可靠性，参照文献[19]，建立了厚度为20 m(变形模量E_s=0.6 MPa，泊松比v=0.4，初始渗透系数k₀=10^-7 m/s)，顶部排水条件下的一维软黏土计算模型，分别对50 kPa与100 kPa荷载下的固结沉降进行了计算，计算结果如图 2所示。由图可见，在两种荷载作用下，数值解与近似解析解的沉降曲线基本吻合，表明本文所建立的有限元方程及程序较为可信。

4.2 非Darcy渗流对固结的影响分析

为分析非Darcy渗流对固结的影响，参照文献[19-20]中算例，针对长度、厚度分别为20 m与10 m的黏土层，建立了二维有限元模型。其中，重度γ=10 kN/m³，静止土压力系数K₀=0.53，压缩模量E_s=0.6 MPa，泊松比v=0.4，初始渗透系数k₀=1.0×10^-9 m/s。均布荷载作用宽度为5 m，顶部荷载施加部分以外排水，底部及两侧不排水。载荷分10级施加，每级间隔3 d，分别为15、25、35、42.7、50、55、60、65、70、75 kPa，计算总时间为4 000 d。模型采用八结点四边形等参单元进行剖分，共划分为204个单元，671个结点，网格剖分如图 3所示。

为分析非Darcy渗流对固结的影响，令α_c=0，对以下两种工况进行分析：1)m=1.2，变动i_L(4，7，10)；2)i_L=4，变动m(1.2，1.3，1.4)，计算结果如图 4~7所示。由图 4的沉降过程曲线可见，当渗流符合Darcy定律时，土层在荷载作用下所产生的沉降较大，而非Darcy渗流产生的沉降相对较小，且随i_L的增加而略有减小，表明非Darcy渗流对土层固结具有一定的延迟作用，在完成相同固结度时，需要耗时更多。图 5给出500 d时，不同i_L下的超孔压分布等值线。由图可见，在500 d时，非Darcy渗流与Darcy渗流下的超孔压分布规律基本一致，在模型左侧加荷面下的超孔压相对较大。对于非Darcy渗流情况，该处超孔压略大，且随i_L的增大而略有增加，由此可见，非Darcy渗流延缓了超孔压的消散，导致土体固结沉降减小。

当固定参数i_L，变动参数m时，非Darcy渗流对固结过程的影响如图 6与图 7所示。由图 6的沉降过程曲线可见，参数m对固结的影响规律与参数i_L相似，随参数m的增加，沉降减小。与i_L相比，参数m的影响相对较大，其对固结的延迟效应也更为明显。图 7给出500 d时的超孔压分布等值线。由图可见，随m的增加，非Darcy渗流现象越明显，超孔压消散越缓慢，延迟效应也越显著。

4.3 渗透系数变化对固结的影响分析

土体在荷载作用下，因排水而压密，渗透系数随固结发展而减小。因此，固定参数i_L与m(i_L=4.0，m=1.0)，变动参数α_c(0.5，1.0，1.5)，分析渗透系数变化对固结过程的影响，计算结果如图 8与图 9所示。由图 8沉降过程曲线可见，加载初期，渗透系数变化的影响较小；随固结时间的增加，渗透系数变化对沉降的影响逐步显现，且随α_c的增大沉降相应减小。由式(4) 可知，α_c取值较大时，尽管渗透系数略有变化，但沉降减小较为明显，与不考虑渗透系数变化情况相比，达到相同沉降时的历时延长，延迟效应明显。图 9给出500 d时的超孔压分布曲线。由图可见，超孔压分布规律基本相似，随α_c的增加，超孔压消散变缓，主要是由于渗透系数变小所致。

4.4 非Darcy渗流与渗透系数变化对固结的影响

为分析非Darcy渗流与渗透系数随固结变化的耦合影响，针对以下5种工况进行计算：1)i_L=4，m=1.2，α_c=0.5；2)i_L=10，m=1.2，α_c=0.5；3)i_L=4，m=1.2，α_c=1.5；4)i_L=4，m=1.5，α_c=0.5；5)i_L=4，m=1.5，α_c=1.5，将计算结果与Darcy渗流结果进行对比，如图 10与图 11所示。由图 10可见，当考虑非Darcy渗流与渗透系数随固结变化的耦合影响时，在100 d以后，沉降曲线出现较大差异，其中，以工况5下的沉降值为最小。尽管参数m与α_c取值范围变化不大，但两者对沉降的影响较大，而参数i_L的影响则相对较小。可以得出，在黏性土固结分析时，若孔压梯度不大，应考虑非Darcy渗流与渗透系数变化的耦合作用，否则，无法真实地反映黏性土的固结过程。图 11给出500 d时，不同工况下的超孔压分布曲线。如图所示，当同时考虑非Darcy渗流与渗透系数变化的作用时，超孔压分布曲线出现显著差异，在模型左侧超孔压消散明显缓慢，以工况5) 最为明显，由超孔压的分布特点较好地印证了上述的沉降分析。

5 结论

1) 在100 d之后，与Darcy渗流相比，非Darcy渗流下的沉降较小，超孔压消散缓慢，且随固结时间的增加，差异越为明显，非Darcy渗流所引起的延迟效应也越来越显著，其中参数m的影响较大，而参数i_L的影响较小，可忽略。原因在于土体固结过程中，土颗粒越来越密实，孔隙越来越小，导致渗透性降低；而随着超孔压的消散，水力梯度随之减小，导致渗流速度降低，由此导致上述规律。

2) 参数α_c对固结过程具有重要影响，且随α_c的增加，沉降逐渐减小，孔压消散趋于缓慢，延迟效应较为明显。

3) 参数m与α_c对黏土固结过程的影响较大，当m=1.5，α_c=1.5时最为显著。建议在进行黏性土固结分析时，应考虑非Darcy渗流与渗透系数随固结变化的耦合作用。