基于分段线性化模型的一维流变固结分析

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

(郑州大学 土木工程学院，河南 郑州 450001)

在基础设施建设中，软黏土地基的工后沉降问题一直都是工程专家和学者们热议的话题之一.这主要是因为软黏土具有高压缩性、低渗透性以及黏滞性等特点，使得其变形是一个非常复杂的过程.在工程实践及试验研究过程中，学者们[1-4]发现，经典的太沙基固结理论并不总能很好地与试验结果或观测数据相吻合，而常常出现不同程度的偏离现象.因此不少学者开始从不同的角度入手对该理论进行修正.一方面，该理论假定地基中的渗流符合Darcy定律，但由于软黏土的低渗透性，水在软黏土地基中的渗流过程并不是总能用Darcy定律描述.为此，有学者[5-12]先后考察了以Hansbo渗流及其简化形式、指数渗流等形式为代表的非Darcy渗流对一维固结的影响，并对经典太沙基理论所不能解释的实测固结度与荷载相关的现象给出了较为合理的解释.另一方面，经典太沙基理论常常假定材料是线弹性的，而实际的饱和软黏土很少如此，所以本构关系的差异性也是导致上述偏离现象的一个因素.为此，一些学者[13-16]引入孔隙比和有效应力的半对数或双曲线经验关系来描述软黏土固结过程中变形的非线性特征.同时相关研究发现，对于软黏土这种黏滞性材料，为反映其变形的弹黏塑性或黏弹性，与时间相关的黏性部分也应被充分考虑到其应力应变关系中[17-18].

回顾现有的黏弹性或弹黏塑性本构模型，大致有两种形式：一种是基于元件理论而提出的本构关系，另一种是基于弹黏塑理论建立的弹黏塑本构模型.前者为更好地模拟试验结果(特别是加载初期)，通常会将弹簧、牛顿黏壶、圣维南体等基本元件或基本模型进行多次组合，这就使得此类本构模型中常常有较多参数，从而为求解过程增加了许多不确定的因素[19].较之前者，后者的本构模型得到了不同程度的优化，参数的数量有所减少，而且还往往具有明确的物理意义.这方面较早的工作可归功于Bjerrum[20]，他在1967年提出了一维的等时间线模型，并将土的变形过程分为可恢复的瞬时压缩和不可恢复的延时压缩两部分.而Yin[21-22]认为土体的黏塑变形部分是不可分开的，提出了“等效时间”的概念，据此建立了一个与时间相关的弹黏塑性(EVP)本构模型.后来，姚仰平[23]在Bjerrum瞬时压缩线的基础上，提出了相对瞬时压缩线的概念，并结合超固结统一硬化模型(UH模型)，提出了可以考虑时间效应的UH模型.这些模型很多已成功应用于饱和软黏土的一维固结分析中，其中基于EVP和考虑时间效应的UH模型等弹黏塑性模型的一维固结理论，已成功模拟了加载初期由软黏土的黏滞性而引起的孔压增高现象[24-25].

实际上，经典太沙基理论及其上述修正理论仅适用于小变形的情况，而相当多的软黏土变形过程应属于大变形的范畴，所以不断有学者提出基于大变形的修正固结理论，并给出某些特殊条件下的解析解或数值解[26-30].为了将大变形理论和小变形理论纳入统一体系,Fox等人进行了很有意义的尝试，并于1997年基于分段线性化方法提出了CS2固结模型[31].和以往的固结模型(理论)不同的是，CS2固结模型并没有去构造复杂的方程式，而是以1个土单元作为研究对象，以单元压缩量等于孔隙水的净流出作为基本条件而进行分析.而今，CS2模型已经被广泛拓展，例如，考虑竖向排水过程的径向固结模型[32]，离心加载下的固结模型[33]，多层土固结模型[34]，电渗固结模型[35-36]等.同时，CS2模型还扩展到其他领域，例如，污染底泥覆盖修复的计算模型[37-38]，以及不同条件下土体固结与污染物迁移的耦合模型[39-43]等.这些推广应用充分说明了CS2模型的准确性和灵活性.

但是，现有的CS2固结模型并没有考虑饱和软黏土的流变特性.因此笔者拟引入考虑时间效应的UH本构关系来描述软黏土的弹黏塑性，从而改进CS2固结模型，并分析UH模型参数对一维固结过程的影响.

1 考虑黏滞效应的改进CS2模型

假定饱和软黏土地基在自重和地面超载q0作用下已固结完成，现于地面施加均布超载q，地基将沿竖向发生渗流和变形.对于初始时刻(图1a)，地基初始厚度为H0，静水面高度为Hw且假定后者在固结过程中保持不变.以地基底面(即固定的基准面)为坐标零点建立z坐标系，规定垂直向上为坐标正方向.将地基在竖向上均匀划分为N个单元，且每个单元都具有相同的横截面面积和初始厚度L0；取每个单元的中心点为该单元的节点，且以节点的初始高度z0,j代表该单元的初始位置.经过时间t后(图1b)，地基厚度变为Ht，且每个单元的厚度及位置高度均随着地基固结而变化.在固结过程中，仍假定每个单元的节点始终处于相应单元的中心点处.

a 初始状态b 固结过程中

1.1 单元中点总应力

参考Fox等[31]的研究，单元j中点的总应力包括顶部原外荷载q0、新施加的外荷载q所产生的附加应力、可压缩地基自重产生的自重应力以及孔隙水压力，即

(1)

γt,j=(Gs+et,j)γw/(1+et,j)

(2)

式(1)、(2)中:Lt,j是单元j在时刻t的厚度；γt,j为单元j在时刻t的饱和重度;γw为孔隙水重度；et,j为单元j在时刻t的孔隙比；Gs是土粒的相对密度，当不计自重对固结进程影响时，可取Gs=1[31].

1.2 本构关系

假定固结过程中地基土的渗流符合Darcy定律，且渗透系数k为常数.为考虑软黏土的黏滞性对固结进程的影响，引入Yao等[23]考虑时间效应的UH模型描述土体的应力应变关系，相应的一维应变增量dεv及塑性应变增量dεv,p可分别表示为

(3)

(4)

式中：dεv,e和dεv,sp分别为由竖向有效应力增量dσ′产生的弹性应变增量和塑性应变增量；dεv,tp为延时应变增量，即由时间作用产生的塑性应变增量.相应的表达式为

(5)

(6)

(7)

式中:Cs、Cc、Cα分别为回弹指数、压缩指数、次固结系数；e0为初始孔隙比；t、ta分别为真实时间和老化时间；t0为参考时间，可取相应的单位时间；M为临界状态应力比；Mf为潜在破坏应力比，即

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

(8)

(9)

χ=M2/(12(3-M))

(10)

(11)

(ta+t0)/t0=R-α

(12)

其中，α=(Cc-Cs)/Cα.

(13)

(14)

式中:Rt,j为单元j在时刻t的超固结参数.

1.3 渗流速度及沉降

以单元j和单元j+1为例(图2)，相邻单元之间的相对渗流速度为

图2 相邻单元渗流示意图Fig.2 Fluid flow between contiguous elements

(15)

式中:在t时刻相邻单元节点间的水力梯度it,j和等效渗透系数kt,s,j分别为

it,j=(ht,j-1-ht,j)/(zt,j-1-zt,j)

(16)

(17)

kt,j表示单元j在时刻t的渗透系数；ht,j、zt,j分别表示在t时刻单元j节点处的总水头和高度，且ht,j=zt,j+ut,j/γw；ut,j表示单元j节点处在时刻t的孔隙水压力.

对于上、下边界，当其不透水时，其边界流速为0；当其透水时，其渗流速度可分别表示为

vt,rf,1=-kt,1(Hw-ht,1)/(Ht-zt,1)

(18)

vt,rf,N=-kt,N(ht,N-Hw)/zt,N

(19)

如图2所示，当经历时间增量Δt后，单元j的高度变为

Lt+Δt,j=Lt,j-(vt,rf,j-1-vt,rf,j)Δt

(20)

则在t+Δt时刻单元j的孔隙比、地基总厚度、地基沉降量可分别表示为

et+Δt,j=[Lt+Δt,j(1+e0,j)]/L0-1

(21)

(22)

St+Δt=H0-Ht+Δt

(23)

1.4 孔隙压力、有效应力及固结度

(24)

为考察整个地基中的孔压变化情况，引入以超孔压表示的平均固结度，即

(25)

2 模型验证

2.1 与文献[25]的对比

胡晶[25]已将UH模型引入到一维固结分析中，并用有限差分法进行了分析.这里取其算例进行算法验证，即地基厚度H0=1 m，Hw=H0，渗透系数k=3.63×10-7m·min-1，沿深度方向均匀分布的初始应力σ0=10 kPa，外荷载q=90 kPa，其他参数见表1.取N=50，Δt=0.1 min，不同次固结系数Cα对应的平均固结度Up与时间t的变化曲线如图3所示.很明显，本文计算结果与文献[25]的一致.

表1 分析所用参数[25]Tab.1 Parameters in analysis[25]

图3 平均固结度随时间的变化曲线Fig.3 Degree of consolidation versus time

2.2 UH模型的验证

李西斌等[44-45]曾利用GDS固结试验系统对萧山黏土进行流变固结试验，其中4组试样在800～1 600 kPa荷载级的试样高度改变量S随时间t的变化曲线示于图4.按本文方法模拟时，取N=50，Δt=0.1 s，土体参数见表2.其中，试样高度H0、压缩指数Cc、回弹指数Cs和次固结系数Cα根据文献[44-45]中的试验结果选用，初始超固结参数R0和渗透系数k根据试算选用，黏土的有效内摩擦角参考文献[25]取为25°.模拟结果示于图4.很明显，UH模型能够较好地描述萧山黏土的固结特性.

表2 模拟试验过程所用参数Tab.2 Parameters in simulation test

图4 试样高度改变量随时间的变化曲线Fig.4 Height of sample versus time

3 参数分析

以下分析中，取H0=Hw=5.0 m，向上单面排水，初始孔隙比e0=1，沿深度方向均匀分布的初始有效应力σ0=100 kPa，外荷载q=100 kPa，渗透系数k=1.0×10-7m·min-1，临界状态应力比M=1.2，N=50，Δt=0.1 min.

3.1 初始超固结参数的影响

不考虑时间效应时，初始超固结参数R0的倒数即为超固结比(OCR)，所以它反映了地基土体的初始超固结程度.为考察R0对地基土固结过程的影响，这里在压缩指数Cc=0.4，回弹指数Cs=0.04和次固结系数Cα=0.008条件下，分别取R0=0.1，0.5和0.9进行了计算.图5给出了不同初始超固结参数R0时底部不排水面附近超孔压ue随时间t的变化曲线.该图表明，与传统太沙基一维固结理论或线性黏弹性固结理论不同的是，在加载的初期，底部不排水面附近的超孔压出现了升高现象，并且R0越大，超孔压峰值就越大，达到峰值的时间也越长.

图6给出了R0=0.9时地基内部不同位置处的超孔压ue随时间的变化曲线.很明显，在加载初期，地基下半部分的超孔压均出现了不同程度的升高现象，只不过在远离不排水面(z较大)处的超孔压峰值相对较小，达到峰值需要的时间也较短，即较高位置处的超孔压更快地进入了消散状态.这种现象并不能用Mandel-Cryer效应进行解释，Yin等[21-22,24]将其归为软黏土的黏性(蠕变)效应，认为它是由于排水不畅而产生的应力松弛现象而导致的.胡晶等[25]同意这一观点，并从理论上对此进行了解释.

图5 初始超固结参数R0对底部超孔压的影响Fig.5 Influence of R0 on excess pore pressure at bottom

图6 地基内部超孔压与时间的关系(R0=0.9)Fig.6 Excess pore pressure in the ground versus time(R0=0.9)

图7给出了不同初始超固结参数R0时以超孔压计算的平均固结度Up随时间t的变化曲线.从图7可以看出，R0越小(即超固结程度越高)，Up就越小，即地基整体孔压消散越快，这也符合李西斌[44-45]试验所揭示的规律.另外，图7表明，当R0=0.9时，加载初期出现了负的平均固结度，这是因为该阶段尽管上部靠近排水面处地基中的超孔压处于消散状态，但远离排水面的下部地基则处于更为显著的超孔压累积过程中.

初始超固结参数R0对地基沉降S的影响示于图8.图8表明，R0越小(即超固结程度越高)，地基沉降就越小，也就是说，前期对地基进行较大荷载的预压可有效地减小地基的沉降量.这符合人们对堆载预压法的定性认识.

图7 初始超固结参数R0对平均固结度的影响Fig.7 Influence of R0 on average degree of consolidation

图8 初始超固结参数R0对地基沉降的影响Fig.8 Influence of R0 on ground settlement

3.2 回弹指数的影响

为考察回弹指数对固结进程的影响，在本节取压缩指数Cc=0.4，次固结系数Cα=0.008，超固结参数R0=0.8，回弹指数依次取为0.04、0.08和0.20.图9给出了不同回弹指数时底部不排水面附近的超孔压ue随时间t的变化曲线.从该图可以发现，回弹指数越大，底部不排水面附近的超孔压越容易出现升高的现象，同时还发现该值的变化对达到孔压峰值所用的时间影响很小.

图10给出了回弹指数对以超孔压计算的平均固结度Up的影响曲线.从图10不难看出，回弹指数越大，地基超孔压的整体消散就越慢，进而使得固结进程变慢；同时还发现回弹指数对平均固结度的影响并不是全程的，后期影响几乎是微弱的.

回弹指数对地面沉降的影响示于图11.从图11可以看出，在固结的后期，地基沉降的差异性越来越明显，随着该值增大，同一时刻的地基总沉降会更大.

图9 回弹指数对底部超孔压的影响Fig.9 Influence of Cs on excess pore pressure at bottom

图10 回弹指数对平均固结度的影响Fig.10 Influence of Cs on average degree of consolidation

图11 回弹指数对地基沉降的影响Fig.11 Influence of Cs on ground settlement

3.3 次固结系数对固结进程的影响

为考察次固结系数对固结进程的影响，这里取压缩指数Cc=0.4，回弹指数Cs=0.04，超固结参数R0=0.8，而次固结系数Cα分别取0.004、0.008和0.016进行了计算，结果分别示于图12至图14.

图12表明，在加载初期，地基底部不透水面附近的超孔压升高现象随次固结系数的增大而更加显著，但次固结系数影响的主要是超孔压的峰值，而对于达到该峰值需要的时间则影响较小.

图13表明，次固结系数Cα越大，同一时刻的平均固结度就越小，即地基超孔压的整体消散就越慢.所以，软黏土的黏滞性延缓了地基超孔压的整体消散过程.

图12 次固结系数Cα对底部超孔压的影响Fig.12 Influence of Cα on excess pore pressure at bottom

图13 次固结系数Cα对平均固结度的影响Fig.13 Influence of Cα on average degree of consolidation

图14表明，次固结系数对地基沉降的影响是显著的.随着次固结系数的增大，同一时刻的地基沉降量增大，且地基沉降进入稳定状态所需的时间也明显变长.根据定义，次固结系数是指孔隙比与时间曲线尾端的斜率.由孔隙比与变形的转化关系，则沉降与时间曲线尾端的斜率也能反映次固结系数的相对大小.从图14可以发现，次固结系数增大，尾端斜率会相应增大，从而相同时间间隔所产生的变形就会增大；同时还可以发现，在尾端近似直线段出现之前，次固结系数对沉降的影响也是不容忽视的，这说明次固结系数对主固结阶段是有影响的，这也印证了按主次固结耦合进行固结分析的必要性[46].

图14 次固结系数Cα对地基沉降的影响Fig.14 Influence of Cα on ground settlement

4 结论

通过引入考虑时间效应的UH本构关系，修正了CS2固结模型，据此分析了UH模型参数对饱和软黏土地基固结进程的影响.参数分析表明：

(1)UH模型中初始超固结参数R0、回弹指数Cs和次固结系数Cα均会影响固结初期出现的超孔压升高现象.

(2)地基中超孔压的整体消散随初始超固结参数的增大而变慢，而地基沉降量会随该值的增大而增大.

(3)回弹指数对以超孔压计算的平均固结度的影响主要表现在固结初期，并且后者会随前者的增大而减小；回弹指数对地基沉降的影响主要表现在固结后期，该值的增大会使地基沉降量变大.

(4)增大次固结系数会使地基超孔压的整体消散变得缓慢，但同时会使得地基沉降量变大，而且地基沉降进入稳定状态所用的时间明显变长.