可液化地基地下结构地震反应特征简化有效应力分析

赵 凯, 王秋哲, 王彦臻, 庄海洋, 陈国兴

(南京工业大学 岩土工程研究所， 南京 210009)

大量震害现象表明地下结构在强地震作用下可能会出现严重的震害以及次生灾害[1-8]。由于地下工程结构受周围土体约束，土-结构动力相互作用关系十分复杂。采用动力有效应力法分析可液化地基地下结构的地震反应特性，可以更科学地解释强地震动作用下土-地下结构体系的运动模式及破坏机理。

目前国内外学者对可液化地基地下结构动力有效应力分析方法进行了系统研究，归纳起来主要有两大类：① 不排水解耦有效应力法：假定在动力作用期间土体中的孔隙水封闭在土骨架中，不向外排出。在总应力分析法的基础上，增加残余孔隙水压力或残余变形的计算。陈国兴等[9-13]基于ABAQUS软件平台，实现了修正Davidenkov本构模型与修正Byrne孔压增量模型的有效应力算法，通过地下结构非线性精细化建模和高效的有限元数值模拟平台，系统地研究了输入地震动特性、场地条对复杂地铁地下车站结构体系非线性地震反应特征、结构损伤与空间效应特性的影响。② 基于Biot动力固结理论的双相介质耦合场统一算法：即采用一定的数值积分格式直接对Biot动力方程进行逐步积分计算的方法。Prevost等[14-15]提出了一种将Biot方程应用于饱和多孔介质动力分析的有限元计算方法，并应用于可液化地基中地下结构的地震反应分析。Khoshnoudian等[16]采用考虑运动硬化的动弹塑性模型描述土的特性，基于Biot部分动力固结模型，数值模拟了可液化地基中内衬隧道的地震反应特性。黄雨等[17]采用了Armstrong-Frederick 型非线性运动硬化准则和非关联流动准则来描述动荷载作用下砂土的循环活动性以及液化强度，基于Biot动力固结理论对液化场地桩基础的地震反应进行了三维有效应力分析。王刚等[18]基于能够模拟饱和砂土液化后大应变响应的弹塑性循环本构模型，采用DIANA-SWANDYNE II对阪神地震中破坏的大开车站进行了有效应力分析，说明了考虑液化大变形的土与地下结构动力相互作用分析方法及其有效性。邹炎等[19]采用Yang等[20]提出的砂土液化大变形本构模型，基于OpenSees对饱和砂土场地的地震反应进行了非线性动力有限元分析，阐述了地基的液化机理。上述研究中均采用复杂的黏弹塑性本构模型(如边界面模型和嵌套面模型)来描述土的动应力-应变关系。这将会导致Biot动力固结方程求解难度大，数值计算收敛困难等难题，且对于使用者的要求较高[21-22]。

大量试验研究表明，Masing型非线性动本构模型能较为理想地描述各类土体的非线性动力特性，模型参数少且物理意义明确，从而使得该类模型成为应用最为广泛的黏弹性动本构模型。然而，该类动本构模型仅提供了一维非线性剪应力-剪应变关系，且不能够描述循环剪切和压缩作用引起的塑性体应变[23]。为此，本文给出了适宜于将一维应力-应变关系向三维空间扩展的等效剪应变算法及加卸载判据，并基于Masing法则构造了三维应力空间中的Davidenkov本构模型；同时补充了笔者提出的剪应力和正应力差耦合作用下剪切-体积应变耦合模型作为Biot动力固结方程中计算超静孔隙水压力的源项；随后基于FLAC3D软件平台，建立了适宜于分析可液化地基中土-地下结构相互作用问题的有效应力算法。最终利用该算法初步研究了可液化地基中地下结构周围场地的地震反应规律。

1 可液化场地动力有效应力分析方法

1.1 控制方程

结合Biot动力固结控制方程，考虑地震荷载作用下饱和土体积变化的相容条件，建立土骨架-孔隙流体两相介质动力有效应力分析方法。假设均质各向同性多孔介质，孔隙水的运动符合Darcy定律，土骨架和孔隙水是可压缩的，控制方程如下

(1)

(2)

式中：∇为拉普拉斯算子；上标“·”为对时间求导；p为超孔压；ns为孔隙率；k为渗透系数；u为土单元位移矢量；ε为土单元应变张量；G为动剪切模量；v为泊松比；γw为孔隙流体重度；Kf为流体压缩系数。

(3)

(4)

双向耦合剪切条件下剪切-体积应变耦合模型表达式如下

(5)

λ=RES-RESt

(6)

1.2 循环本构方程

采用Masing法则构造三维应力空间中的Davidenkov本构模型描述土骨架的非线性和滞回性[28-29]，如图1所示。初始骨架曲线的表达式为

τ=G·γ=Gmax·γ·[1-H(γ)]

(7)

其中

(8)

图1 修正Davidenkov应力-应变关系示意图Fig.1 Schematic illustration of hypsteretic stress-straincurve of modified Denvidenkov model

由Masing法则构造的滞回曲线表达式为

(9)

式中：τ、γ分别为剪应力和剪应变；Gmax为最大动剪切模量；A、B和γc为土的试验参数；γc为加卸载转折点处的应变。

Masing型非线性动本构模型的核心是时变剪切模量的计算，其准确性与等效剪应变的选取密切相关，即等效剪应变算法的选取是本构关系由一维空间向三维空间扩展的关键问题。以应变偏量的第二不变量描述的等效剪应变能综合反映三维空间中土体的应力-应变关系，而直接采用应变偏量的第二不变量作为等效剪应变时，其数值只能描述一维应力-应变关系曲线中应变大于零的部分，需将等效剪应变改由增量形式通过加卸载判断后再次计算[30]

(10)

(11)

(12)

(13)

1.3 土骨架模型与孔压增量模型的耦合方式

动力有效应力分析方法可以合理地考虑土骨架与孔隙水两相介质相互作用以及场地液化过程中土体力学性质的变化。如果在动荷载施加过程中孔隙流体的排水通道受阻，对于同一土骨架模型，可以采用上一状态的骨架曲线来预测当前状态的骨架曲线，即初始剪切模量随结构孔压(亦即随有效应力变化)的累积而减小。

(14)

2 模型验证

Martin等基于干砂的常幅循环直剪试验，建立了剪应变幅值相关的体应变增量模型。Byrne[31]将Martin等采用多项式给出的体应变增量表达式改为指数函数形式，该模型仅有两个参数且物理意义也更为明确，广泛用于地下结构地震反应分析[32-33]。

(15)

式中，参数C1和C2可由修正标准贯入锤击数Nl或相对密度Dr(以百分数表示)确定

C1=7 600(Dr)-2.5或C1=8.7(N1)-1.25

(16)

C2=0.4/C1

(17)

图2 不同分析模型预测值与实测值的对比Fig.2 Comparison between the measured and predicted excess pore pressure time-histories based on different models

表1 模型计算参数Tab.1 Mechanical parameters of the soil element

图2为本文模型、修正Byrne模型以及室内试验获得的超孔压增长时程曲线。由图可知，采用修正Byrne模型能够较好地预测一维剪应力路径下超孔压的发展，而在正应力差和剪应力耦合作用下(即椭圆形和圆形应力路径)的预测结果远小于试验结果。笔者建立的模型预测结果与试验数据均吻合较好，能够反映应力路径变化对超孔压发展的影响规律。这是由于修正Byrne模型由直剪试验的数据建立，难以反映正应力差和剪应力的耦合变化引起土单元的超孔压发展规律。饱和砂土的孔压增长规律与应力路径密切相关，对于RCS相同而a/b不同的双向耦合剪切试验，对试样所施加应力路径不一致，因此受到的剪切效果也不相同。值得一提的是，应力孔压(超孔压波动值)随着循环周数增加逐步放大，这是由于结构孔压的增长引起土骨架模量衰减导致的。

图3以Chen等给出的应变控制不排水循环三轴加载试验验证本文模型的适用性。可以看出，本文模型预测的结果对于不同相对密度，不同加载类型的试验结果均较为理想。

图3 应变控制不排水循环三轴试验预测值与实测值对比

3 算 例

3.1 计算模型建立

选取某砂质场地矩形截面隧道为算例，土层总厚度为60 m，底面为不透水基岩，隧道为三跨结构，总宽L= 40 m，边跨宽16 m，中跨宽8 m，高为10 m。顶板、底板及侧墙厚度均取1 m，隧道埋深为1 m，结构采用C45混凝土。模型如图4所示，采用四节点等应变有限差分单元，网格大小约为0.25 m × 0.25 m，土-结相互作用接触面为法向硬接触，接触面切向采用摩擦接触。土-结接触面法向刚度kn= 350 GPa，剪切刚度ks= 350 GPa，摩擦系数μ=0.7。场地的剪切波速Vs= 203 m/s，渗透系数K= 1×10-5cm/s，饱和重度γ= 19.8 kN/m3，临界阻尼比D取值为5%。模型计算参数见表1。

图4 模型网格示意图Fig.4 Grid layout of computation domain in detail

场地两侧边界采用竖向约束、水平向加弹簧阻尼器的黏弹性边界，场地底部基岩面仅施加竖向约束。假定地震动为基岩面竖直向上传播的SH波，由模型底部边界输入[34]。采用地震工程中常用的Kobe波作为输入地震动，将输入地震波的持续时间截取为30 s。截取后的地震波加速度时程及傅氏谱如图5所示，在输入地震动时将峰值调幅为0.1g，在计算过程中，仅考虑地震波单向垂直入射。

图5 输入地震动原始加速度时程及傅里叶谱Fig.5 Original acceleration time-history and Fourier spectraof the input Kobe wave

3.2 结果分析

图6为场地不同位置处的土单元应力路径((σz-σx)/2 ～τzx)。由图可见，远场P1(距离结构侧墙100 m)在地震作用过程中处于正应力不变条件下规则的往复水平剪切状态，应力路径为直线型。而结构周围场地土单元体经历的应力路径则明显不同，表现为正应力差与剪应力连续变化，P2(距离结构侧墙30 m)的应力路径为斜椭圆型，而P3(距离结构侧墙3 m)为斜直线型。这主要是由于地震波垂直向上传播引起远场土体纯剪应力变化，当地震波传播至在土-结构接触界面时发生反射与透射现象，使地下结构周围场地土的受力状态更为复杂，从而在土单元体中形成比远场土单元体更为复杂的动应力路径，同时产生动剪应力和动正应力，两者不同的耦合作用形成了斜椭圆型和斜直线型应力路径。

图7给出了隧道结构周围场地不同位置处的超孔压增长时程曲线，并将本文模型和修正Byrne模型的结果进行了对比。从图中可以看出，本文模型和修正Byrne模型计算的远场P1的超孔压增长时程曲线基本一致，而在结构周围P2，P3位置处，本文模型较修正Byrne模型计算得出的超孔压发展更快，最终的超孔压值也更大，且超孔压波动性(即应力孔压)也更显著。

图6 隧道周围场地动应力路径Fig.6 Predicted stress paths of the soil elements around the tunnel

图7 隧道周围场地超孔压发展时程Fig.7 Time histories of excess pore pressure around the tunnel

图8描述了隧道结构周围场地的液化进程。由图可见，与修正Byrne模型相比，本文模型计算的土层液化范围明显更大。结构周围出现液化，尤其结构底板标高处的土层液化范围较大，在底板角度位置出现了明显的液化集中区，且顶板和底板最先发生液化，随后液化区域沿着结构周围逐渐扩展形成整片区域。隧道结构的出现改变了自由场的液化特性。距离结构较远区域，液化最先发生在场地表面，并逐渐向下扩展。本文模型和Byrne模型计算得出的最大液化深度达到了隧道结构底板下4 m和2 m。在距离结构100 m处，在地震作用结束后，液化深度相同，均达到1.5 m。

图8 隧道周围场地渐进液化云图Fig.8 Development of liquefaction zone around the tunnel

结合前文室内试验结果可知，剪应力和正应力差的耦合剪切作用显著加大了结构孔压的累积速度，加剧了土骨架模量的衰退，进而导致液化区域的扩大。综合图6，7和8可以看出，基于直剪试验数据建立的修正Byrne模型仅能够较好地反映远场P1的动力响应，而低估了隧道结构周围P2和P3的超孔压发展和液化进程。这与3.1节对于室内试验模拟的规律吻合较好。

4 结 论

本文采用适宜于将一维应力-应变关系向三维空间扩展的等效剪应变算法及加卸载判据，基于Masing法则构造了三维应力空间中的Davidenkov本构模型，补充了剪应力和正应力差耦合作用下剪切-体积应变耦合模型作为Biot动力固结方程中计算超静孔隙水压力的源项，基于FLAC3D软件平台，建立了适宜于分析可液化地基中土-地下结构相互作用问题的有效应力算法，研究了地下结构所处可液化地基的地震反应规律，可以得到以下主要结论：

(1) 饱和砂土液化特性与循环加载模式和应力路径密切相关。对饱和砂土不排水循环扭剪试验模拟结果表明：相较于单一循环直剪试验结果建立的修正Byrne模型，本文方法可以更为合理表征复杂应力路径下饱和砂土的孔压发展规律及液化过程。

(2) 土-地下结构动力相互作用显著改变了地下结构周围场地地震应力以及所形成的应力路径。地震波垂直向上传播引起远场土体保持正应力不变条件下规则的往复水平剪应力，而地下结构周围土体处于往复剪应力和正应力差耦合剪切状态。

(3) 结合对于室内试验模拟的规律可知，采用修正Byrne模型可以较好预测地下结构远场的动力响应，却低估了结构周围场地的超孔压累积速度和液化区域，本文方法能够较好反映土-地下结构动力相互作用对结构周围场地动力响应的影响。