考虑孔压积累-消散耦合的等效有限元方法

钱建固 ，林志果 ，马 霄

（1.同济大学 地下建筑与工程系，上海 200092；2.同济大学 岩土及地下工程教育部重点试验室，上海 2 000923）

1 引 言

交通荷载往复作用下饱和软土路基会发生显著的工后沉降，已提出的预测路基长期沉降的方法主要有统计分析法、弹塑性数值方法和基于试验建立的经验显式模型分析法。统计分析法基于统计和概率分析，难以推广到一般的工程案例中，也难以对非常规的工程问题或不同的工况做出准确的分析。弹塑性数值法理论较为精确，但参数取值的困难，在大数目循环荷载下的计算效率依然是其在工程实践中推广使用的阻碍。显式模型具有物理概念明确、参数少、计算方便等特点，利于在工程界的推广应用。

Monismith 等[1]建立了循环累积应变与加载次数之间关系的经验公式。Li 等[2]总结了之前众多的室内试验改进了上述模型，改进模型中考虑了动应力水平qd的影响。Chai 等[3]在以上工作的基础上，考虑了土体初始静偏应力对首次循环加载应变的影响。黄茂松等[4]在上海地区第④层饱和软黏土等向、偏压固结不排水循环加载试验基础上建立了同时反映等向、偏压固结情况下轴向循环累积塑性应变和孔压发展规律的模型。

显式模型一般与分层总和法结合，但分层总和法无法满足路基沉降的位移协调问题。马霄等[5]提出了等效有限元法，克服了上述缺陷，但他采用的瞬时积累逐步消散的孔压模型与实际中孔压积累与消散同时进行并不相符。针对以上不足的改进是本文所做的工作。

2 经验显式模型与等效有限元

2.1 循环累积塑性应变及孔压显式模型

黄茂松等[4]在上海地区第④层饱和软黏土等向、偏压固结不排水循环加载试验基础上建立了同时反映等向、偏压固结情况下轴向循环累积塑性应变和孔压发展规律的模型。

累积轴向应变显式模型：

式中：η=qd/pc，qd为动偏应力；pc为围压；N为循环次数；a1、m、c、b为试验确定的常数；pa为大气压。

累积孔压显式模型：

式中：u*为孔压；η=qd/pc；au1、mu、cu、bu为试验确定的常数，其他符号意义同前。

2.2 单元应变的计算

交通荷载循环作用下土体单元应变包括两部分，如图1 所示[5]。εr为轴向循环累积塑性应变，由于这部分应变对应不排水过程，因此径向应变为-εr。对于二维平面问题，通过坐标变换由轴向和径向应变得到全局坐标系下x、y 方向的应变，公式为

式中：α为轴向应变方向与Y 轴的夹角。

第二部分应变为孔压消散引起的体应变εv。由e-lnp 曲线可得平均主应力增量引起的体应变增量表达式。

式中：p0为前期固结压力；e0为初始孔隙比。以u=0时εv=0为边界条件，对式（4）积分可得

式中：u为累积孔压；U为固结度。在平面应变状态下对应的x、y 方向应变为

图1 单元应变Fig.1 Strain of element

2.3 等效有限元方法

所谓等效有限元法[5]，指借鉴初应变方法的计算思路来实现显式模型的有限元计算。显式模型与一般的本构模型同样都是描述材料在外力作用下的变形特性，不同的是显式模型中出于简化计算的需要没有直接建立应力-应变的关系，而是建立了加载次数与应变的关系，所以不能直接用于有限元计算。

由显示模型可以得到单元在循环荷载下的总应变表达式，以此作为初应变，假设一个弹性矩阵，可以得到初应变{ε0} 对应的初应力{σ0} 与初应变引起的等效节点荷载向量{ Fε0}：

由叠加原理知在原平衡条件上增加的初应力引起的节点位移向量可由式（8）解出。

由于弹性矩阵的选取对总刚度矩阵与等效节点荷载向量的影响完全同步，可以选取任意弹性矩阵而不会影响节点位移向量。

以上过程相当于用初应变和弹性矩阵算出一个等效节点荷载向量，让这个等效节点荷载在有限元模型中“挤”和“拉”出交通荷载作用引起的变形，通过最小位能原理保证了这种变形是受力平衡且几何连续的，从而实现了以显式模型为本构模型的等效有限元计算。

3 孔压消散模型

3.1 瞬时积累逐步消散孔压模型

土体单元体应变计算公式（5）中孔压消散值的计算需要选取孔压消散模型。以往都是采用瞬时积累逐步消散的模式，如图2 所示。OA为基于不排水循环加载试验建立的孔压显式模型，即式（2）计算值。比如求 t0时刻孔压消散值，瞬时积累逐步消散的孔压模型假定全部孔压在 t0/2时刻（图中的C点）全部积累，这意味着用孔压荷载DEF 来代替孔压荷载DOC。显然DOC 大于DEF，且DOC 到 t0的时间长。因此，瞬时积累逐步消散的孔压模型将会低估了孔压消散值。

图2 瞬时积累逐步消散的孔压模式示意图Fig.2 Sketch of accumulating instantly and dissipating gradually pore water pressure model

3.2 孔压积累与消散相互耦合的孔压模型

张建民等[6]将大沙基固结理论与不排水条件下震动孔压增长模式相耦合，提出了描述自然排水条件下震动孔压长消的微分方程式：

式中：u*为震动孔压增长模式；Cvx、Cvy分别为水和竖直方向的固结系数。

张建民等[6]通过 Sturm-Liouvill 系统以及Duhamel 积分与多重Fourier 展开给出了式（9）水平自由界面的饱和土层的解析解。求解的边界条件为底部与两端不排水，顶部自由排水：

式中：h、l 分别为求解土层的厚度与长度，求得的解析解表达式为

式中：

式中：x0、y0、t、u*分别为求解的水平、竖直坐标、时间与震动孔压增长模式（即文中孔压显示模型）。

文中显式模型 u*建立了孔压与加载次数的关系，而式（12）要求 u*为与时间的关系。只需规定加载次数随时间的增长模式N=f(t)，即可将文中显示模型（见式（2））变为与时间的关系 u*(t)。因此式（5）中uU 这一项即可按下式计算：

3.3 两种孔压模式计算结果对比

这里计算图3 中A 点的孔压长消及固结度。土层只有上方为排水边界，土体各参数取表1、2 中的淤泥质黏土，孔压显示模型中的η 取0.15，pc取150 kPa。将参数代入式（2）、（11）中可得A 点孔压变化，如图4 所示。图中，曲线OB 所反映的孔压变化规律与试验[7]得到的规律一致。两种孔压模型所计算的固结度的差异，如图5 所示。

图3 计算区域Fig.3 Computational region

图4 孔压变化规律Fig.4 Rules of pore pressure change

图5 固结度随时间变化规律Fig.5 Rules of degree of consolidation change over time

4 算例分析

上海浦东国际机场第一跑道于1999年10月建成投入运营，长4 000 m，宽60 m。岩土工程勘察资料表明，机场地基主要由黏性土、粉性土和砂土组成，土层分布如图6 所示，各层参数见表1。由表可知，粉土层的固结系数远大于黏土，可将其看做排水面。第①层土在地下水位线以上，不考虑其孔压累积。第④层土上层是砂质粉土，近似认为是自由排水面，下层是黏土，近似认为不排水。第⑤层土排水条件很差，而且动应力比很小，累积的孔压也很小，计算中忽略其孔压值。

表1 土层物理力学参数Table 1 Physico-mechanical parameters of every stratum

图6 浦东国际机场地质断面图(单位:m)Fig.6 Engineering geological profile of Pudong International Airport(unit:m)

4.1 有限元模型

浦东机场第一跑道长4 000 m，宽60 m，可以用平面应变有限元模型来分析。建模选取了一个典型的跑道截面，考虑对称性，对跑道中轴线一侧建模。模型中跑道1/2 的宽度为30 m，结构层总厚度为100 cm，其中混凝土道面厚46 cm，其下为3 层18 cm 厚的二灰碎石层。模型中路基土体厚40 m，宽100 m。实际土层分布第2 层中上面2 个亚层经常缺失，模型中简化为1 层。

4.2 动应力场计算

本文荷载大小及作用位置参照空客A300 大飞机。空客A300 最大起飞重量171.4 t，降落时进场速度为249 km/h。A300 采用三点式起落架，起落架有1 组鼻轮和2 组主轮组成，如图7 所示。考虑到主轮承担了大部分荷载，而且鼻轮距离主轮较远，计算时忽略鼻轮承受的荷载。飞机在起降过程中，运行速率较大，引起的土体动荷载大于静载。为了简单起见，这里通过经验关系式（14）[8]将动载等效为静载：

图7 飞机轮载分布Fig.7 Wheel load distribution of aircraft

将以上等效荷载加到3D 有限元模型中提取轮载下方跑道截面的应力场，再通过式（15）计算整个截面内各单元的最大偏应力：

4.3 显示模型参数选取

循环累积轴向应变及累积孔压显式模型的参数见表2。表中，黏土和粉土的显式模型参数根据钱建固等[9]不排水心形应力路径循环扭剪试验结果拟合得到。

表2 循环累积轴向应变及累积孔压公式参数Table 2 Parameters of cylic accumulated axial strain and cumulative pore pressure formulas

4.4 计算结果

图8为飞机动载作用下跑道路基竖向沉降云图。

图8 累积沉降云图(单位:m)Fig.8 Nephogram of cumulated strains(unit:m)

从图8 中计算结果显示，路基浅层竖向沉降较大，最大沉降出现在两侧轮载作用之间；在路基浅层，沉降主要分布在轮载作用范围附近；而在路基深层，由于应力扩散，沉降的分布范围有所扩大。随着时间的增加，沉降的量值与分布范围也有所增加。

4.5 计算结果对比分析

王广德等[10]监测了该机场第一跑道自建成以来至2009年的沉降数据，其中一项为P265 断面上的差异沉降，可用于与本文结果作对比分析，对比结果如图9 所示。

图9 差异沉降对比Fig.9 Comparison of differential settlement

实测数据与等效有限元计算得到的差异沉降发展规律基本相同，但数值大小存在一定的差异，造成这种差异的原因可能是：（1）交通荷载之外其它的因素也对差异沉降的发展有贡献，如飞机降落时产生的冲击荷载，以及机场跑道自重荷载的不均匀或排水条件的不均匀造成的差异沉降。（2）建模与实际的差异。本文计算模型对土层做了简化，实际的土层分层在表层附近差异较大，常有亚层缺失，将第2 层土简化为一层并按粉土参数取值可能使计算得到的差异沉降偏低。本文用拟静力方法计算了等效动应力，也可能带来一定的误差。

5 结 论

（1）等效有限元计算得到的差异沉降发展规律与实测数据基本相同，在跑道投入运营初期差异沉降发展较快，运营一段时间以后增速放缓，趋于稳定。

（2）考虑孔压积累-消散相互耦合的孔压模型所计算的沉降值比瞬时积累逐步消散的孔压模型所计算的沉降值大。这种差异在运营初期更加明显，随着时间的推移二者逐步接近。

（3）本案例中，孔压消散引起的沉降占总沉降约10%的比例。随着时间的推移，固结沉降所占的比例还会有所增长。

[1]MONISMITH C L,OGAWA N,FREENE C R.Permanent deformation characteristics of subsoil due to repeated loading[J].Transportation Research Record,1975,537:1－17.

[2]LI D,SELIG E T.Cumulative plastic deformation for fine-grained subgrade soils[J].Journal of Geotechnical Engineering,1996,122(12):1006－1013.

[3]CHAI J C,MIURA N.Traffic-load-induced permanent deformation of road on soft subsoil[J].Journal of Geotechnical and Geoenvironmental Engineering,2002,128(11):907－916.

[4]黄茂松,姚兆明.循环荷载下饱和软黏土的累积变形显式模型[J].岩土工程学报,2011,33(3):325－331.HUANG Mao-song,YAO Zhao-ming.Explicit model for cumulative strain of saturated clay subjected to cyclicloading[J].Chinese Journal of Geotechnical Engineering,2011,33(3):325－331.

[5]马霄,钱建固,韩黎明,等.交通动载下路基长期运营沉降等效有限元分析[J].岩土工程学报,2013,35(增刊2):910－913.MA Xiao,QIAN Jian-gu,HAN Li-ming,et al.Equivalent finite element method for long-term settlement of subgrade induced by traffic load[J].Chinese Journal of Geotechnical Engineering,2013,35(Supp.2):910－913.

[6]张建民,谢定义,刘家沛.饱和砂层震动孔隙水压力长消的解析算法[J].水利学报,1992,47(12):70－80.ZHANG Jian-min,XIE Ding-yi,LIU Jia-pei.An analytical algorithm for calculating pore water pressure of sandy soil under condition of vibration[J].Journal of Hydraulic Engineering,1992,47(12):70－80.

[7]HYODO M,YASUHARA K,HIRAO K.Prediction of clay behavior in undarained and partially drained cyclic triaxial tests[J].Journal of Soils and Foundations,1992,32(4):117－127.

[8]马霄.交通动载下软土路基长期运营沉降等效有限元分析[D].上海:同济大学,2014.

[9]钱建固,王永刚,张甲峰,等.交通动载下饱和软黏土累计变形的不排水循环扭剪试验[J].岩土工程学报,2013,35(10):1790－1798.QIAN Jian-gu,WANG Yong-gang,ZHANG Jia-feng,et al.Undrained cyclic torsion shear tests on permanent deformation responses of soft saturated clay to traffic loadings[J].Chinese Journal of Geotechnical Engineering,2013,35(10):1790－1798.

[10]王广德,韩黎明,柴震林,等.上海浦东机场一跑道地基沉降规律[J].工程地质学报,2012,20(1):131－137.WANG Guang-de,HAN Li-ming,CHAI Zhen-lin,et al.Subgrade settlement rules of first runway of Pudong Airport in Shanghai[J].Journal of Engineering Geology,2012,20(1):131－137.