对于软黏土地基，世界范围内常采用砂井排水固结法进行处理，并已取得了非常不错的效果[1]。针对砂井固结理论的研究，Carrillo将砂井排水固结分为径向渗流和竖向渗流两部分，这样简化了径竖向组合渗流固结方程[2]。Barron首先在Darcy渗流条件下给出了自由竖向应变和等竖向应变两种极端假设下的线弹性砂井径向固结的解析解[3]，对两者计算结果的对比表明在较长时间后差异较小。Yoshikuni等在Barron自由应变假定基础上，引入竖井流量方程并考虑井阻作用得出固结解析解[4]。此后，文献[5-7]分别介绍了考虑井阻及涂抹效应条件下进一步完善的Barron理论解。不过上述砂井固结理论仍属于经典Terzaghi固结理论的范畴，即土体变形和渗流均服从线性规律。但工程实测资料表明：其理论预测值与实际值常常存在较大的差异。许多学者通过软黏土试验发现[8-10]：软黏土渗流存在非达西特性，并提出了不同的渗流模型，较为常见的有Hansbo渗流、指数渗流和三参数渗流等。为此，一些学者开始引入非达西渗流对传统的砂井固结理论进行修正，例如，Hansbo在砂井固结分析时首先将Hansbo渗流引入，给出了仅考虑幂函数渗流段影响的竖井地基固结近似解[11-12]；刘忠玉等同时考虑了低速渗流幂函数曲线段和较高速渗流直线段[13]；Walke等引入的则是指数渗流模式[14]。这些学者的研究均表明土体的非达西特性能够显著延缓砂井地基的固结。需要说明的是：上述研究均假定土体的应力-应变之间是线性关系，即压缩系数为常数。但实际上土的宏观变形主要是因为颗粒间的位置变化，不同的应力条件下，相同应力增量引起的应变增量常表现出非线性。因此很早就有学者开始探讨土体变形非线性对固结性状的影响。文献[15-18]介绍了引入压缩指数描述土骨架变形的非线性关系，并给出了砂井地基固结的解答，其结果表明若不考虑土体的非线性常会高估地基的固结程度。随后，文献[19-20]又分析了涂抹区渗透系数变化、井阻效应及外荷载等因素对地基固结进程的影响，使砂井固结理论更加严密。此外，许多固结试验[21-22]表明：软黏土具有明显的流变特性，它不仅仅是引起次固结沉降的主要原因，而且对主固结进程也会有比较大的影响。为此，国内外的众多学者开始引入相应的本构关系进行流变固结分析。一种是引入元件模型，即将Hooke弹簧、Venant刚塑体和Newton黏壶进行串联或并联构成。比如，赵维炳引入广义Viogt模型，推求了砂井地基流变固结普遍理论解[23]。谢康和团队在引入Merchant模型的基础上，分析了井阻与涂抹效应、边界透水性及变荷载等因素对砂井流变固结性状的影响[24-25]。但是，基于元件模型的微分型本构表达式的阶数有时过高，且参数过多，不便于在工程应用中推广[26]。因此，一些学者引进另外一种方法，即建立黏弹塑性本构模型来分析其流变特性。较为经典的弹黏塑性模型有基于速率效应的Leroueil模型[27]、基于等效时间概念的殷建华模型[28]、基于圧缩曲线的蠕变模型[29]等。其中，殷建华模型考虑了等效时间、参考时间线、瞬时时间线、极限时间线等影响，较为全面地描述了土体的流变特性。后来，姚仰平等引入瞬时压缩线的概念，建立了新的弹黏塑性框架，并在一维蠕变规律的基础上，结合统一硬化(UH)模型提出了考虑时间效应的UH模型[30-31]。该模型相比于修正剑桥模型仅增加了一个次固结系数，且参数均可以通过室内常规试验确定。胡晶等已将其引入到一维固结分析之中，解释了加载初期地基内部出现的孔压上升而应力却松弛的现象[32]。
为进一步探讨软黏土砂井地基的流变固结机理，本文引入考虑时间效应的UH模型描述土体弹黏塑性应力-应变关系，同时考虑Hansbo渗流及变渗透系数的影响对Barron理想径向砂井固结理论进行修正，并初步分析相应参数对固结进程的影响。
1 分析模型
1.1 考虑时间效应的UH模型
本文引入姚仰平等提出的考虑时间效应的UH模型[31]来模拟土体一维非线性竖向变形，即：
(1)
其中
式中：dεv为竖向应变增量；
和
为由有效应力增量dσ′引起的弹性应变增量和塑性应变增量；
为时间作用产生的黏塑性应变增量；Cc为压缩指数；Cs为回弹指数；Cα为次固结系数；ta为老化时间；t0为初始时间，根据实际情况取1 s、1 min、1 h等；M为临界状态应力比；Mf为潜在破坏应力比；φ为土的内摩擦角；e0为初始孔隙比；R为反映超固结程度的参数，在不考虑时间效应时，其初始值R0是超固结比OCR的倒数；
为回弹线与瞬时压缩线交点的竖向应力，类似于先期固结压力；
为非弹性应变。
R与老化时间的关系为：
(2)
其中
1.2 Hansbo渗流模型
引入描述非达西渗流的Hansbo渗流模型[8]代替传统的Darcy渗流模型(图1)，即：
(3)
其中 i0=i1(m-1)/m
式中：m为试验常数；i1为线性渗流的起始水力坡降；i为水力坡降；i0为线性渗流的计算起始水力坡降；c为幂指数形式的渗透系数；k为线性关系式中的渗透系数。
图1 Darcy渗流及Hansbo渗流模型
Fig.1 Darcy’s and Hansbo’s flow models
将i0、c的表达式代入式(3)得：
(4)
1.3 控制方程
设某均质地基为各向同性软黏土，在自重条件下已完全固结。现均匀设置竖向砂井，并打穿土层，然后瞬间在地面上施加均布竖向荷载p0。设砂井直径为2rw，影响区直径为2re，记井径比n=re/rw。
工程实践和理论分析表明：当土层较厚时，径向排水对固结进程起主要作用，竖向渗流的影响可以忽略[13,33]。同时Barron和Richart均认为当时间较大时，自由竖向应变和等竖向应变两种假定下的固结度计算值差异较小[3，34]，有学者认为自由应变假定更接近实际[4-5]。因此本文仅对自由竖向应变假定下的轴对称径向固结问题进行分析，同时假定：
1)土颗粒和水均不可压缩，土体仅发生竖向压密变形。
2)竖向砂井表面及内部的超静孔隙水压u为零。
3)仅考虑土体径向渗流，忽略竖向渗流。
4)不考虑施工过程对竖向砂井周边土体的扰动，即假定无涂抹效应。
5)p0在天然地基中引起的竖向附加应力σ沿深度不变。
6)加载开始时，荷载引起的全部应力由孔隙水承担。
7)固结过程中，到砂井中心距离r处的渗流速度vr可用Hansbo渗流方程(式(4))表示。
引入普遍公认的经验关系式(5)描述固结过程中渗透系数k与孔隙比e的关系：
(5)
式中：e0为初始孔隙比；k0为相应的初始渗透系数；Ck为渗透指数。
在侧限压缩条件下，式(5)可表示为：
(6)
由渗流变形连续性条件可得：
(7)
根据有效应力原理σ=σ′+u，可得：
(8)
对于轴对称问题，根据微元体的径向出水量等同于同一时间内的压缩量，可得控制方程为：
(9)
其中
式中：γw为孔隙水的重度。
式(9)的初始及边界条件为：
rw≤r≤re
(10)
t>0
(11)
t>0
(12)
2 方程的有限差分求解
对控制方程式(9)采用隐式有限差分法求解时，首先以Δr为步长将砂井影响区(rw≤r≤re)土层由内到外均匀离散为N层，且以Δt为步长将时间t均匀离散，并用i′和j分别表示空间结点和时间结点，则：
(13)
其中：
下标i′=1,2,3，…,N-1，上标j=1,2,3，…。
对初始条件式(10)和边界条件式(11)、(12)分别离散化可得：
i′=0,1,2，…,N (14)
j=1,2,…，N (15)
(16)
式(13)～(16)组成封闭方程组，可采用迭代法求解。这里采用Fortran软件编程计算有效应力σ′后，再通过有效应力原理算得孔压u：
(17)
引入按孔压定义的平均固结度描述砂井地基中孔压的消散情况：
(18)
假定j时间点在Δt区间[ri′,ri′+1]内的孔压线性分布，即：
(19)
将式(19)代入式(18)可得时间点j在Δt的固结度为：
(20)
其中
3 算法验证
令Cα=0，M=Mf，m=1或i1=0，本课题即退化为Berry等所分析的课题[15]。这里取其中一算例来验证。为便于与Berry等的课题[15]对比，引用无量纲参数：
X=r/re。空间与时间步长分别取ΔR=0.01，ΔT=10-6，迭代精度为10-7，其余参数取值与Berry相同：
将本文与Berry计算结果示于图2。很明显，两者结果基本重合，这说明本文算法是有效的。
——本文解Cc/Ck=0.5；---本文解Cc/Ck=1.0；……本文解Cc/Ck=1.5；
—Berry解Cc/Ck=0.5；
—Berry解Cc/Ck=1.0；
—Berry解Cc/Ck=1.5。
图2 孔压变化曲线(X=0.57)
Fig.2 Pore pressure variation curves (X=0.57)
4 参数分析
土体的超固结程度、回弹指数、渗透指数、压缩指数等对砂井固结有着不同程度的影响，已有不少学者进行过研究[15-17,35-36]。这里仅着重探讨土体流变特性及渗流非Darcy特性对砂井地基径向固结性状的影响，并且进行下述单因素分析法，参数选取如下：
迭代精度为10-7。
4.1 次固结系数的影响
首先考虑黏滞性对流变固结过程的影响。图3、图4分别给出了Cα=0.00，0.02，0.04，0.06时砂井影响区外边缘处的孔压和径向平均固结度随时间的变化曲线。
——Cα=0.06；---Cα=0.04；……Cα=0.02；-·-·-Cα=0.00。
图3 Cα对砂井影响区边缘处孔压的影响
Fig.3 Influence of Cα on the pore water pressure on the boundary of the influence zone of the sand drain
——Cα=0.06；---Cα=0.04；……Cα=0.02；-·-·-Cα=0.00。
图4 Cα对径向平均固结度的影响
Fig.4 Influence of Cα on radial average degree of consolidation
由图3可观察到：当Cα=0.00(即不考虑土体的流变效应)时，该处的孔压在固结初期变化很小，当时间超过1 860 min后，孔压随时间增长而消散，这与传统砂井固结理论一致。但是在考虑流变效应(Cα=0.02，0.04和0.06)的情况下，却与之有所不同，即在固结初期，该处的孔压不降反升，并在某一时刻达到最大值后再随时间增长而消散；而且次固结系数越大，孔压最大值就越大，达到最大值的时间也越长。例如对应Cα=0.02，0.04，0.06的最大孔压分别为153.51，165.77，175.38 kPa，达到相应的时间分别为1 920，2 360，3 050 min。这一现象与采用比奥固结理论分析饱和黏性土多维固结时的曼德尔效应[37]类似，但产生机理不同，因为本文的理论分析是在太沙基固结理论范畴内进行的。这里暂且称之为“类曼德尔效应”。实际上在考虑流变效应对一维饱和黏性土地基进行固结分析时，也会出现该现象[32,38-41]。文献[38-41]分别将其归结为主、次固结耦合效应或流变效应。这里的数据也支持这一观点。
从图4可以看出：在平均固结度UP达到工程一般要求的90%时，不考虑流变效应(Cα=0.00)时仅需83 350 min，而次固结系数Cα为0.02，0.04，0.06时，对应的时间分别为132 720，210 530，335 440 min，所需时长相较Cα为0.00时依次为1.59，2.53和4.02倍。这表明土体的流变效应对砂井地基的整体孔压消散有延缓作用，且随着Cα的增大延缓作用更为明显。另外，在固结初期(t=100 min)，Cα大于零时固结度甚至为负。例如，Cα=0.06时，固结度UP=-0.1，这说明考虑黏滞性时土体内部较大范围内会出现类曼德尔效应。
4.2 Hansbo渗流模型参数的影响
在考虑流变效应的前提下(Cα=0.04)，图5和图6所示的不同Hansbo渗流参数时砂井影响区外边缘处孔压随时间的变化曲线表明：不论是Darcy渗流还是Hansbo渗流，固结初期都出现了类曼德尔效应，而且随Hansbo渗流参数的增大，孔压最大值略有增大。以图5为例，假定渗流参数i1=5不变，m取1.0，1.2，1.5和1.8时，孔压峰值分别为164.16，164.88，165.77，166.46 kPa，达到相应峰值的时间依次为1 450，1 800，2 360，2 930 min。这说明曼德尔效应与渗流模式没有直接关系，但Hansbo渗流对该效应有一定的增强作用。此外，该处孔压在达到最大值后，就开始随时间增长而消散。但Hansbo渗流参数越大，孔压消散就越慢。如要使孔压消散到同一值u=25 kPa，参数m=1.0，1.2，1.5，1.8对应的时间分别为118 520，135 480，162 060，191 380 min。
——Darcy渗流；---m=1.2；……m=1.5；-·-·-m=1.8。
图5 参数m对砂井影响区边缘处孔压的影响(i1=5)
Fig.5 Influence of parameter m on the pore water pressure on the boundary of the influence zone of the sand drain (i1=5)
——Darcy渗流；---i1=5；……i1=10；-·-·-i1=15。
图6 参数i1对砂井影响区边缘处孔压的影响(m=1.5)
Fig.6 Influence of parameter i1 on the pore water pressure on the boundary of the influence zone of the sand drain (m=1.5)
——Darcy渗流；---m=1.2；……m=1.5；-·-·-m=1.8。
图7 参数m对径向平均固结度的影响(i1=5)
Fig.7 Influence of parameter m on radial mean degree of consolidation (i1=5)
——Darcy渗流；---i1=5；……i1=10；-·-·-i1=15。
图8 参数i1对径向平均固结度的影响(m=1.5)
Fig.8 Influence of parameter i1 on radial mean degree of consolidation (m=1.5)
图7和图8为不同Hansbo渗流参数下径向平均固结度UP随时间t的变化曲线。可见，相较Darcy渗流，按Hansbo渗流计算的径向平均固结度要小，且Hansbo渗流参数越大，径向固结度就越小。如在固结后期，对应同一时刻t=105 min，图7中沿Hansbo渗流参数m增大方向的固结度UP依次为0.79，0.76，0.72，0.70，与刘忠玉等获得的不考虑黏滞性时的结论[13]类似，即Hansbo渗流会延缓土体中孔压的整体消散。
5 结束语
引入考虑时间效应的UH模型描述土体的弹黏塑性变形，同时考虑Hansbo渗流及变渗透系数的影响，基于自由应变假定，重新推导了砂井地基径向固结方程。通过数值分析次固结系数Cα和Hansbo渗流参数m、i1对固结性状的影响，得出如下结论：
1)在固结初期，砂井影响区外边缘附近的孔压出现了明显升高的现象，这与土体的黏滞性或主次固结耦合机制有着密切关系，且次固结系数Cα越大，这种现象越显著。
2)Hansbo渗流增强了固结初期的孔压升高现象，且阻碍了固结中后期的孔压消散。
3)土体的黏滞特性及渗流的非达西特性对砂井地基的整体孔压消散均有着明显的延缓作用。
[1] 殷宗泽.土工原理[M].北京:中国水利水电出版社,2007.
[2] CARRILLO N.Simple Two and Three Dimensional Cases in the Theory of Consolidation of Soils[J].Journal of Mathematics and Physics,1942,21:1-5.
[3] BARRON R A.Consolidation of Fine Grained Soils by Drain Wells[J].Transactions of ASCE,1948,113:718-742.
[4] YOSHIKUNI H,NAKANADO H.Consolidation of Fine-Grained Soils by Drain Well with Filter Permeability[J].Soil and Foundations,1974,14(2):35-46.
[5] ONOUE A.Consolidation by Vertical Drains Taking Well Resistance and Smear into Consideration[J].Soil and Foundations,1988,28(4):165-174.
[6] HANSBO S,JAMIOLKOWSKI M,KOK L.Consolidation by Vertical Drains[J].Geotechique,1981,31(1):45-46.
[7] 谢康和,曾国熙.等应变条件下的砂井地基固结解析理论[J].岩土工程学报,1989,11(2):3-17.
[8] HANSBO S.Consolidation of Clay with Special Reference to Influence of Vertical Drains[C]//Proceeding of Swedish Geotechnical Institute.Stockholm:Swedish Geotechnical Institute,1960:45-50.
[9] SLEPICKA F.Contribution to the Solution of the Filtration Law[C]//International Union of Geodesy and Geophysics,Commission of Subterranean Waters.1960:245-258.
[10] MILLER R J,LOW P E.Threshold Gradient for Water Flow in Clay Systems[J].Soil Society of American Journal,1963,27(6):605-609.
[11] HANSBO S.Aspects of Vertical Drain Design:Darcian or Non-Darcian Flow[J].Geotechique,1997,47(5):983-992.
[12] HANSBO S.Consolidation Equation Valid for Both Darcian or Non-Darcian Flow[J].Geotechique,2001,51(1):51-54.
[13] 刘忠玉,焦阳.基于Hansbo渗流的理想砂井地基固结分析[J].岩土工程学报,2015,37(5):792-801.
[14] WALKER R,INDRARATNA B,RUJIKIATKAMJORN C.Vertical Drain Consolidation with Non-Darcy Flow and Void-Ratio-Dependent Compressibility and Permeability[J].Geotechnique,2012,62(11):985-997.
[15] BERRY P L,WILKNSON W B.The Radial Consolidation of Clay Soils[J].Geotechnique,1969,19(2):253-284.
[16] BASAK P,MADHAV R.Analytical Solutions of Sand Drain Problems[J].Journal of the Geotechnical Engineering Division,ASCE,1978,104 (1):129-135.
[17] LEKHA K R,KRISHNASWAMY N R,BASAK P.Consolidation of Clay by Sand Drain Under Time-Dependent Loading[J].Journal of the Geotechnical and Geoenvironmental Engineering,1998,124 (1):91-94.
[18] 周琦,刘汉龙,陈志波.考虑固结参数变化时砂井地基的非线性径向固结[J].岩土力学,2007,28(增刊1):855-858.
[19] 张海丘,高广运,雷丹.考虑3种非线性关系的径向排水固结解析解[J].工程地质学报,2015,23(4):681-686.
[20] 郭霄,谢康和,卢萌盟,等.直排式真空预压法下竖井地基的非线性固结解析解[J].中南大学学报(自然科学版),2018,49(2):384-392.
[21] 李西斌,谢康和,陈福全.考虑软土流变特性和应力历史的一维固结与渗透试验[J].水利学报,2013,44(1):18-25.
[22] 刘俊新,杨春和,谢强,等.基于流变和固结理论的非饱和红层路堤沉降机制研究[J].岩土力学,2015,36(5):1295-1305.
[23] 赵维炳.广义Voigt模型模拟的饱和土体轴对称固结理论解[J].河海大学学报,1988,16(5):47-56.
[24] 刘兴旺,谢康和,潘秋元,等.竖向排水井地基黏弹性固结解析理论[J].土木工程学报,1998,31(1):10-19.
[25] 王瑞春,谢康和.半透水边界的竖向排水井地基黏弹性固结分析[J].长江科学院院报,2001,18(6):33-36.
[26] 袁静,龚晓南,益德清.岩土流变模型的比较研究[J].岩石力学与工程学报,2001,20(6):772-779.
[27] LEROUEIL S,KABBAJ M,TAVENAS F,et al.Stress-Strain-Strain Rate Relation for the Compressibility of Sensitive Natural Clays[J].Géotechnique,1985,35(2):159-180.
[28] YIN J H,GRAHAM J.Viscous-Elastic-Plastic Modelling of One-Dimensional Time-Dependent Behaviour of Clays[J].Canadian Geotechnical Journal,1989,26(2):199-209.
[29] KUTTER B L,SATHIALINGAM N.Elastic Viscoplastic Modelling of the Rate-Dependent Behaviour of Clays[J].Géotechnique,1992,42(3):427-441.
[30] 姚仰平.土的统一硬化模型及其发展[J].工业建筑,2008,38(8):1-5.
[31] 姚仰平,孔令明,胡晶.考虑时间效应的UH模型[J].中国科学:技术科学,2013,43(3):298-314.
[32] 胡晶,姚仰平.基于考虑时间效应UH模型的一维固结分析[J].北京航空航天大学学报,2015,41(8):1492-1498.
[33] LEO C J.Equal Strain Consolidation by Vertical Drains[J].Journal of Geotechnical and Geoenvironmental Engineering,2004,130(3):316-327.
[34] RICHART JR F E.A Review of the Theories for Sand Drains[C]//Proceeding of ASCE,1957,83(SM3):1-38.
[35] 秦爱芳,李天义,裴杨从琪,等.半渗透边界下非饱和土砂井地基固结特性[J].工程地质学报,2019,27(2):390-397.
[36] 张玉国,万东阳,郑言林,等.考虑径向渗透系数变化的真空预压竖井地基固结解析解[J].岩土力学,2019,40(9):1-9.
[37] MANDEL J.Consolidation Des Sols[J].Geotechnique,1953,9(3):287-299.
[38] 殷建华,Clark J I.土体与时间相关的一维应力-应变性状、弹黏塑性模型和固结分析[J].岩土力学,1994,15(3):65-80.
[39] 殷建华,Clark J I.土体与时间相关的-维应力-应变性状、弹黏塑性模型和团结分析(续)[J].岩土力学,1994,15(4):65-75.
[40] 丁洲祥,袁大军,朱合华.一维大变形主、次固结耦合新模型[J].岩土力学,2010,31(8):2367-2372.
[41] 仇玉良,丁洲祥.一维小变形主、次固结耦合理论模型分析[J].岩土力学,2012,33(7):1957-1964.