考虑非Darcy渗流的砂井地基弹黏塑性固结分析

2021-03-22刘忠玉徐传永朱新牧张家超夏洋洋

哈尔滨工业大学学报 2021年5期

刘忠玉，徐传永，朱新牧，张家超，夏洋洋

(郑州大学土木工程学院,郑州 450001)

沿海地区由于饱和软黏土地基的高压缩性，在荷载作用下建筑容易产生较大沉降，因此,需对地基进行预先处理，以加快地基的固结速率，其中设计理论大多基于砂井固结理论.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]中，体积应变增量的表达式为

(1)

(2)

(3)

(4)

M=6sinφ/(3-sinφ),

(5)

(6)

χ=M2/12(3-M),

(7)

(8)

(ta+t0)/t0=R-α,

(9)

α=(Cc-Cs)/Cα.

(10)

1.2 控制方程推导

(11)

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

logK=logK0-(e0-e)/Ck.

(12)

式中：K0为初始孔隙比e0对应的初始渗透系数；Ck为渗透指数，反映渗透系数的变化程度.

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

K=K0exp[-(1+e0)εvln 10/Ck].

(13)

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

(14)

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

(15)

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

(16)

根据有效应力原理，可得

(17)

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

(18)

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

(19)

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

(20)

(21)

(22)

2 方程的有限差分解法

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

(23)

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

(24)

(25)

(26)

(27)

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

(28)

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

3 解法验证

图1 本文数值解与耿雪玉[37]解析解对比(Z=1，X=2.5)

4 参数分析

4.1 次固结系数Cα的影响

图2给出次固结系数Cα对某些特殊位置处孔压的影响曲线.可以看出，在流变固结初期(t=1 000 min)，Cα>0时，远离排水边界处的孔压都超过了初始值p0=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倍.很明显，土体的黏滞效应延缓了地基的固结过程.

图2 Cα对孔压的影响

图3给出不同次固结系数Cα时的砂井地基平均固结度Up随时间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]对一维流变固结分析时所得出的结论一致.

图3 Cα对平均固结度的影响

4.2 初始超固结参数R0的影响

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

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

图4 R0对孔压的影响

图5 R0对平均固结度的影响

4.3 压缩指数Cc的影响

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

图6 Cc对孔压的影响

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

图7 Cc对平均固结度的影响

4.4 回弹指数Cs的影响

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

图8 Cs对孔压的影响

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

图9 Cs对平均固结度的影响

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

4.5 渗透指数Ck的影响

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

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

图10 Ck对孔压的影响

图11 Ck对平均固结度的影响

4.6 Hansbo渗流参数m和i1的影响

为分析Hansbo渗流模型对固结的影响，图12，13分别给出了不同m或i1时的孔压分布曲线.二图表明，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或i1值越大，同一位置处的孔压越大，即孔压消散越慢.例如，图13(a)中，z=5 m，r=0.9 m，在Darcy渗流模式下(i1=0时)算得的孔压值为33 kPa，当i1取5，10和15时，对应Hansbo渗流模式下计算的结果分别为47，58和67 kPa，大约分别是Darcy渗流下孔压的1.42，1.76和2.03倍，可见，Hansbo渗流延缓了地基中孔压的消散.究其原因是随着固结时间的延长，较大范围土体内的孔压有所消散，水力梯度值也变得较小，当小于直线渗流起始水力梯度i1时，渗流是很缓慢的，且随着Hansbo渗流参数的增大，渗流速度更小，孔隙水更不易排出，所以，孔压消散也更为缓慢.

图12 m值对孔压的影响

图13 i1值对孔压的影响

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

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