0.   引言

可液化地基地下结构的地震反应受周围土层复杂的力学过程控制，主要表现为土的滞回非线性和残余孔压增长引起土的性状变化（动剪切模量降低，发生循环软化或者液化）。因此，科学地解析强地震动作用下可液化地基-地下结构体系的运动模式及破坏机理，适宜采用动力有效应力分析方法。

目前国内外学者对可液化场地有效应力分析方法进行了系统研究，归纳起来主要有两大类。

（1）强耦合非线性有效应力分析法

基于流-固耦合框架，采用先进的弹塑性本构模型（如边界面模型和嵌套面模型）计算孔压增长和应力–应变反应，如庄海洋等^[1]改进了Elgamal等^[2]和Yang等^[3]提出的嵌套面硬化规则和弹塑性模量确定方法，基于ABAQUS软件平台对地铁地下车站结构体系精细化建模，系统地研究了地铁地下车站结构体系的非线性地震反应特性。基于两相介质的Biot动力固结理论，Khoshnoudian等^[4]采用考虑各向同性运动硬化的动弹塑性模型描述土的特性，模拟了可液化地基中内衬隧道的地震反应特性。Bao等^[5]通过FE-FD DBLEAVES编码程序对大型地铁地下隧道进行了二维地震反应分析。一些学者探讨砂土液化后的土体性态，建立了可描述砂土液化大变形行为的本构模型。王刚等^[6]基于DIANA-SWANDYNE Ⅱ软件平台，采用能够模拟饱和砂土液化后大应变响应的弹塑性循环本构模型，建立二维模型对阪神地震中破坏的大开车站进行了有效应力分析。王睿等^[7]基于OpenSees有限元计算平台，将液化大变形本构映射在三维空间，发展了一种动力时程有效应力分析法，并将其用于分析可液化地基中桩基础三维模型的震动响应，对比离心机振动台试验，验证了此方法的有效性。基于临界状态框架建立的更为先进的土体本构模型适宜模拟液化土的行为特性，然而存在模型参数较多，求解Biot动力方程难度大，数值计算收敛困难等难题。

（2）弱耦合非线性有效应力分析法

采用半经验的基于应力或应变的孔压模型与滞回应力–应变模型相耦合方法预测土的强度变化。Tropeano等^[8]提出了一种弱耦合形式的有效应力分析方法，通过SCOSSA平台设计程序进行一维地震反应分析；Olson等^[9]基于Deepsoil V6.1软件平台建立一维分析模型，使用GQ/H土骨架曲线耦合Dobry孔压模型，对水平地表液化地震反应进行了广泛验证。Finn等^[10]基于其提出的一种弱耦合形式有效应力分析方法设计计算软件TARA3，可实现场地模型的二维有效应力分析^[11]。陈国兴等^[12]、王彦臻等^[13]实现了Davidenkov本构模型与修正Byrne孔压增量模型耦合的有效应力算法，通过地下结构精细化建模，应用高效的有限元数值模拟平台ABAQUS^[14-15]，系统地研究了输入地震动特性、场地条件等对复杂地铁地下车站结构体系非线性地震反应特征、结构损伤与空间效应特性的影响。目前，有效应力分析法研究多基于一维剪切模型或二维平面假定，难以充分考虑地震发生时，土–地下结构共同运动产生的空间效应^[16]。建立精细化三维数值模型，简洁高效、深入全面地揭示地下结构关键部位在地震反应中的受力和变形特点十分必要。

陈国兴等^[12]、王彦臻等^[13]将土的非线性滞回模型与剪切-体积应变耦合的孔压增量模型相结合，考虑循环加载过程中循环刚度退化与超静孔隙水压力增长的耦合关系，提出了新的一维场地地震反应的弱耦合非线性有效应力分析方法，通过1995年阪神地震Port Island液化场地井下台阵强震记录的反演验证了该有效应力分析法的有效性。本文采用适宜于将一维应力–应变关系和孔压增量模型向三维空间扩展的等效剪应变算法及加卸载判据，将该弱耦合非线性有效应力分析方法拓展至三维土–地下结构动力相互作用分析。针对已完成可液化场地三层三跨地铁车站结构振动台试验资料^[17]，开展了对应模型试验的三维数值模拟。通过三维数值模拟和振动试验结果的对比，既能验证数值分析结果的可靠性，同时也能验证模型试验结果的合理性。

1.   一维弱耦合有效应力分析法回顾

陈国兴等^[12]提出的弱耦合场地反应有效应力分析法中，土的初始骨架曲线由Davidenkov模型描述：

式中，

骨架曲线上的时变剪切模量为

由广义non-Masing法则^[12]构造的应力–应变滞回曲线表达式为

应力–应变滞回曲线上的时变剪切模量为

其中，

式中τ，$\gamma$为剪应力和剪应变；G_max为初始剪切模量；$\gamma_r$为参考剪应变；A，B为土性相关的无量纲常数；n_c为滞回圈比例系数；${\tau _{\text{c}}}$和${\gamma _{\text{c}}}$为当前加载—卸载反转点的剪应力和剪应变；${\tau _{{\text{ex}}}}$和${\gamma _{{\text{ex}}}}$为历史最大剪应力和剪应变；不规则加载—卸载期间，滞回圈向上弯曲时"±"取正、向下弯曲时取负，如图 1所示。

图  1  Davidenkov模型骨架曲线图

Figure  1.  Davidenkov model skeleton curve

下载: 全尺寸图片 幻灯片

孔压的增长则是通过耦合由循环剪切作用引起的体应变增量模型及回弹模量公式来刻画^[18]，

式中，C₁，C₂，C₅，m，n为与土性有关的拟合参数，${\gamma _{{\text{tv}}}}$为体积门槛剪应变，E_r为土的回弹模量。

2.   弱耦合有效应力分析法的三维拓展

上述一维场地地震反应有效应力分析法的核心是土骨架时变剪切模量的计算和孔压增长的预测，而两者均是剪应变的函数。因此，将该有效应力分析法由一维空间向三维空间扩展的关键问题是选取能够反映三维应力应变状态的等效剪应变。以应变偏量的第二不变量描述的等效剪应变能综合反映三维空间中土体的应力，应变关系，而直接采用应变偏量的第二不变量作为等效剪应变时，只能描述一维应力–应变关系曲线中应变大于零的部分，需将等效剪应变改由增量形式通过加卸载判断后再次计算：

式中$\gamma _{{\text{eq}}}^t$，$\gamma _{{\text{eq}}}^{t + \Delta t}$为t和$t + \Delta t$时刻的等效剪应变。加载时，sign = 1；卸载时，sig = -1。$\Delta \gamma _{{\text{incre}}}^{t + \Delta t}(e_{ij}^o)$为等效剪应变增量，$\gamma _{{\text{gen}}}^t(e_{ij}^o)$为$t$时刻的广义剪应变，由t时刻应变偏量的第二不变量$J_{2\varepsilon }^t(e_{ij}^o)$确定。

由式（8）可知，更新剪切模量时只需得到当前时刻应变与转折点处剪应变差值$(\gamma - {\gamma _{\text{c}}})$，因此定义式中应变偏张量$e_{ij}^o$：以加卸载转折点处的应变张量${e_{ij, c}}$为起点，转折点过后计算得到的应变张量为${e_{ij}}$终点的张量，

在等效剪应变的算法中，当$\Delta \gamma _{{\text{incre}}}^{t + \Delta t}(e_{ij}^o) < 0$时可直接判断应力发生转向，对于三维应力空间中加卸载判断简单有效。

应力–应变滞回曲线段t+Δt时刻的时变体积模量${K^t}$可由下式确定，

式中，$\nu$为泊松比。

由此，对t时刻的三维空间应力张量$\sigma _{ij}^t$进行更新，

式中，$C_{ijkl}^t{\text{ = }}$

3.   可液化地铁车站数值模型

3.1   三维数值模型

传感器布置如图 2所示，相应传感器分别由字母A、LD及W标识，分别标识地基土的加速度反应、位移反应及孔压发展。模型结构具体尺寸及主、次观测面布置如图 3所示，考虑微粒混凝土的承受能力及结构整体的受力性能，模型结构附加配重为480 kg，具体试验设计见文献[17]。

图  2  传感器布置

Figure  2.  Layout of sensors

下载: 全尺寸图片 幻灯片

图  3  地基土–车站结构的三维有限元分析模型

Figure  3.  3D finite element model for soil-subway station structure

下载: 全尺寸图片 幻灯片

建立可液化场地三层三跨地铁车站结构的三维有限元模型，如图 3所示。ABAQUS/Explicit求解器中，条件稳定增量步与单元尺寸直接相关，各部件网格划分需以均匀为主。根据不同构件的几何形状，模型土的网格尺寸采用0.05 m划分，可配合后续对比分析时能选到模型地基中传感器的对应位置；模型结构墙、板的网格尺寸为0.025 m；模型结构中柱的竖向网格尺寸为0.01 m，横截面网格为0.003 m。整个结构体系的单元数为101200个，节点个数121300，自由度共363900。模型土和模型结构的单元类型均采用八结点线性六面体减缩单元（C3D8R单元），增强型沙漏控制。

3.2   材料模型及参数

（1）土体

采用第1节中论述有效应力法分析振动台试验中的黏土与饱和南京细砂动力特性，Davidenkov模型与剪切–体应变耦合增量模型相关参数列于表 1中饱和南京细砂的孔压增量模型的参数已由文献[18]给出。试验给出的震前的平均剪切波速为92.1 m/s，模拟时模型场地剪切波速随深度的修正如图 3所示，最大剪切模量由G =ρV_s²计算。

表  1  黏土与南京细砂的模型参数

Table  1.  Model parameters of clay and Nanjing fine sand

土层	密度/(g·cm-3)	动泊松比	Davidenkov模型						孔压增量模型
			A	B	a1	a2	a3		C1	C2	C3	γtv/%	m	n
黏土	1.75	0.49	1.2	0.35	2.5×10-4	0	0		—	—	—	—	—	—
南京细砂	1.83	0.49	1.02	0.43	4.1×10-4	0.5	0.45		1.051	0.143	1.25	0.02	0.345	6.689

下载: 导出CSV 

| 显示表格

（2）混凝土

动力荷载作用下的混凝土材料可由混凝土塑性损伤模型（CDP模型）描述，采用非关联多向硬化塑性和各向同性损伤弹性相结合的方法^[19]来描述混凝土破碎产生的不可恢复的损伤。Chen等^[14]采用CDP模型来模拟微粒混凝土的动力特性，模拟的应变反应规律与试验结果基本一致。模型结构的总体积为0.106 m³，考虑附加配重后的惯性效应，微粒混凝土的等效密度为6.486 g/cm³。假设混凝土受荷由拉伸转向压缩，裂缝闭合则压缩刚度恢复；受荷由压缩转向拉伸时，出现压碎裂纹则拉伸刚度不可恢复，微粒混凝土CDP模型的参数列于表 2中。

表  2  微粒混凝土CDP模型的参数

Table  2.  Model parameters of micro-concrete

弹性模量Eo /GPa	拉伸变量ωt	压缩变量ωc	泊松比	ψ/(°)	ε	fbo/fc	Kc	u
7.5	0	1	0.18	36.31	0.1	1.16	2/3	0.005

下载: 导出CSV 

| 显示表格

（3）土–结构接触面模拟与边界处理

土–结构的动力相互作用主要由结构和地基的相对刚度及两者接触界面变形的传递能力决定。受刚度较大的模型土箱的侧向约束作用，地基与地下结构之间的变形传递能力较强，土与结构的相互作用方式采用tie接触形式，并定义刚度较大的模型结构为主面，刚度较小的模型地基为从面；结构梁、板与柱子的接触对中，网格较大的梁、板为主面，网格较小的柱子为从面。

考虑模型地基–地下结构体系的惯性作用以及叠层土箱的剪切效应，动力分析步中模型土体四周及底部均设置黏弹性人工边界，并在模型土体底部施加x方向的地震动。

模拟土–地下结构初始应力场。首先模拟地下结构应力分布：约束土体四周及底部的法向位移，使用传递分析向地下结构模型中传递应力结果以平衡结构模型中的应力集中；其次模拟土体应力分布：同样约束条件下，使用传递分析向土体模型中传递应力结果以平衡土体模型的位移。平衡后土体表层最大位移7.024×10^-6 m。本次计算时间增量步设置为t=2×10^-6。

（4）加载工况

试验分别选取2008年汶川地震发生时松潘台站和什邡八角地震台站的加速度记录，经1︰2时间压缩及幅值调整后，共计4个工况，如表 3所示，未调幅的振动台面输入地震动的加速度时程及傅里叶谱如图 4所示。

表  3  振动台试验加载工况

Table  3.  Schemes of shaking table tests

地震动	峰值加速度/g	工况	持时/s
松潘波	0.1	SP-1	100
	0.5	SP-2	100
什邡八角波	0.1	SF-1	100
	0.5	SF-2	100

下载: 导出CSV 

| 显示表格

图  4  输入地震动

Figure  4.  Input motions

下载: 全尺寸图片 幻灯片

4.   结果对比分析

4.1   地基土孔压反应对比

图 5对比了PGA为0.1g，0.5g什邡波与松潘波作用下，模型结构侧边及结构底部（孔压计W9、W10、W6和W7）测得的孔压比时程与数值模拟得到的孔压比时程。

图  5  数值模拟与试验给出的孔压比时程曲线的对比

Figure  5.  Comparison between predicted and measured results for time histories of pore pressure ratios

下载: 全尺寸图片 幻灯片

其中，图中灰色图层为地基土孔压比从快速发展至完全稳定的时间区间，为数值模拟结果。可以看出，对于不同频谱特性及强度的输入地震动，数值模拟与试验得到的孔压比峰值均能基本吻合。PGA = 0.1g的什邡波和松潘波作用下，模型地基土的孔压比基本小于0.4，砂土尚未液化并拥有较高的强度；PGA = 0.5g什邡波和松潘波作用下，各测点的孔压比均达到1或接近1，砂土的液化程度较高。对比各测点的孔压比时程曲线可知，数值模拟和试验得到的孔压比时程均有明显的门槛现象且出现时间也十分接近，当输入地震动强度增强时该门槛值变小，这与孔压增量模型中引入的门槛剪应变的意义一致。其中，PGA = 0.1g松潘波作用下，试验给出的孔压比在快速发展后趋于平缓，随后再次进入快速发展阶段，数值模拟结果也较好地反映了该工况下孔压比的发展规律。PGA = 0.5g松潘波作用下，孔压计W6处试验得到的孔压比时程存在明显的振荡反应，使得试验得到的孔压比的发展时间较数值模拟得到的孔压比的发展时间更长。由于孔压在瞬时降低后迅速恢复，而数值计算时采用的孔压增量模型主要描述了孔压比的累积过程，但此阶段数值模拟得到的孔压比均无明显增长。此外，由图 5数值模拟得到的不同位置的孔压比时程可知（图中灰色图层），随着输入地震动峰值加速度的增大，地基土孔压比突破门槛的时间更快，土体完全液化所需的时间也更短。

4.2   地基土加速度反应对比

图 6对比了数值模拟与试验给出的A1测点（模型结构中线的地表处）的加速度时程及其傅氏谱。

图  6  数值模拟与试验给出的土中加速度计A1，A3的加速度时程与Fourier谱的对比

Figure  6.  Comparison between predicted and measured results for time histories of acceleration and Fourier spectra for points A1 and A3 in model soil

下载: 全尺寸图片 幻灯片

可以看出，数值模拟得到的加速度反应的时程、幅值以及傅氏谱均与试验结果一致。加速计A1处：PGA为0.1g，0.5g什邡波作用下，数值模拟和试验得到的峰值加速度分别相差16.6%，8.1%；PGA = 0.5g松潘波作用下，峰值加速度相差28.7%。加速度计A3处：PGA为0.1g，0.5g什邡波作用下，数值模拟和试验得到的峰值加速度分别相差35.7%，1.8%。对比PGA为0.1g，0.5g什邡波作用下加速度计A1处的Fourier谱可知，随着输入地震动幅值的增大，数值模拟和试验得到的Fourier谱的高频成分不断衰减，低频成分不断放大；PGA = 0.5g的两条波作用下，主频均在1.6 ~ 2.4 Hz。当地基发生液化时，由输入地震动频谱差异引起的地表加速度反应的频域特性差异逐渐减小，傅里叶谱谱形也趋于一致。加速度计A3处，PGA为0.1g，0.5g的什邡波作用下，数值模拟和试验得到的Fourier谱可知，随着输入地震动幅值的增大，Fourier谱除低频发育外，高频成分仍较为丰富，主频仍聚焦在6.8 Hz。

为对比主、次观测面地基土动力反应的差异性，图 7给出了PGA为0.1g，0.5g的什邡波、松潘波作用下数值模拟和试验的对比结果。其中，峰值加速度放大系数随深度的变化关系分别选取主观测面A4，A5，A6，A7加速度计以及两个次观测面对应位置A9，A11，A10，A12加速度计。

图  7  数值模拟与试验给出的土中峰值加速度放大系数随深度的变化关系

Figure  7.  Comparison between predicted and measured results for peak acceleration amplification factors with model soil depth

下载: 全尺寸图片 幻灯片

由图 7中各工况下振动台试验结果可知，对于模型结构埋置深度以下的地基土，主、次观测面的峰值加速度放大系数基本一致，而在模型结构埋置区域至地表处，主、次观测面峰值加速度放大系数的差异性较大。究其原因，可能主要是由振动台模型试验的动力边界模拟不准确造成的。振动台试验采用的框架式叠状剪切箱虽削弱了边界上地震波的反射或散射效应，仍无法消除模型箱对模型土–地下结构体系振动特性的影响^[20-21]，而数值模拟可以较好的模拟半无限场地中地下结构体系的地震反应，因而所得主、次观测面峰值加速度放大系数随深度的变化规律亦一致。隧道深度之上黏土层表面峰值加速度反应，试验与数值结果之间的差异较为显著。观察图 6中模型土表层加速度时程，发现数值计算峰值加速度发生时间位于砂土层孔压比到达峰值之后，较之试验结果提前，推测认为：加速度计A1，A3位于模型土表层，由于表面土层较软，加速度计与土层之间发生较大的相对错动，导致数值模拟结果和试验记录之间存在一定的差距。另外，本文提出的弱耦合有效应力分析方法不能充分模拟表层黏土循环刚度退化与超静孔隙水压力增长的耦合关系也是引起差距的一个方面^[21]。对于表层黏土动力特性的模拟还需进一步的研究。就整体而言，数值模拟和振动台试验得到的加速度时程和的傅氏谱吻合度较高，表明数值模拟的效果较为理想。

4.3   地基土动力反应时频特性

地基土的加速度和孔压比反应，存在较为密切的联系。由Chen等^[15]开展的振动台试验结果显示，地基土的加速度反应达到峰值时的时间几乎与孔压比峰值出现的时刻相近。以PGA = 0.5g什邡波和松潘波两个工况为例，对数值模拟和试验给出的主观测面加速度计A5处的加速度时程进行时-频-能量变换，与临近位置处孔压计W2的孔压比反应进行对比，对比如图 8所示。对于地震反应这一类典型的非平稳信号，HHT变换可自适应的对其进行分解，能更细致的表征动力反应的局部特征^[22]。

图  8  加速度反应的时-频-能量特性及与孔压比时域发展的相关性

Figure  8.  Relationship between time-frequency-energy spectra of accelerations and time histories of EPWP ratio of model soil

下载: 全尺寸图片 幻灯片

由数值和试验给出的加速度反应的时-频-能量谱和孔压比时程关系可得：①数值模拟和试验中，瞬时能量最大值均出现在孔压比的快速发展阶段；②在振动台试验中，循环剪切作用下砂土的瞬时剪胀使得孔压出现振荡现象，并伴随加速度反应瞬时能量及瞬时卓越频率（为某一时刻最大瞬时能量所对应的频率）的增加；③在振动台试验中，长持时大幅度的孔压消散也将引起瞬时卓越频率的增加；④在数值模拟中，PGA = 0.5g松潘波作用下，瞬时卓越频率随孔压的上升呈现明显的衰减特性。

此外，对比数值模拟和试验给出加速度反应的时-频-能量谱可知，PGA = 0.5g什邡波作用下，试验给出的能量最大值出现在约10.5 s，该时刻的瞬时卓越频率约为6.2 Hz；数值模拟得到的能量最大值出现在约9.4 s，该时刻的瞬时卓越频率约5.9 Hz。在PGA = 0.5g松潘波作用下，试验给出的能量最大值出现在约22.0 s，该时刻的瞬时卓越频率约为2.4 Hz；数值模拟得到的能量最大值显现在约22.1 s，该时刻的瞬时卓越频率约为2.0 Hz。综上，循环剪切作用下，虽然砂土瞬时剪胀以及振动过程中的孔压消散对加速度幅值的时域发展特性及孔压消散段的瞬时卓越频率存在影响，但数值模拟得到的能量聚焦时刻以及该时刻对应的瞬时卓越频率与试验结果均接近，考虑残余孔压对动力反应影响的有效应力法能反映振动台试验中地基土动力反应的一些基本特征，可满足实际应用需求。

5.   结论

本文采用适宜于将一维应力–应变关系和孔压增量模型向三维空间扩展的等效剪应变算法及加卸载判据，考虑循环加载过程中土骨架循环刚度退化与超静孔隙水压增长的耦合关系，建立了三维土–地下结构动力相互作用分析的弱耦合有效应力法。针对已完成可液化场地三层三跨地铁车站结构振动台试验资料，开展了对应模型试验的三维数值模拟，有3点主要结论。

（1）土的性状变化显著影响土–地下结构动力相互作用。本文提出的三维弱耦合有效应力法可以合理表征饱和砂土超静孔隙水压力上升和土骨架循环刚度退化对土–地下结构动力相互作用的影响规律。循环剪切作用下，土骨架瞬时剪胀以及振动过程中的超孔压消散对地基土加速度幅值的时域发展特性及孔压消散段的瞬时卓越频率存在一定影响，但数值模拟得到的能量聚焦时刻及该时刻对应的瞬时卓越频率与振动台试验结果均吻合较好。

（2）土骨架有效应力水平显著影响超孔压比的发展规律，同一深度处模型结构底板下部的孔压比峰值小于远离模型结构地基土的孔压比峰值。在模型结构附加质量块后，模型结构下部土体的初始有效围压明显大于原有自由场处土体，造成模型结构下部地基土的孔压比发展较同深度处远离模型结构的地基土孔压比低。振动台试验不完备的密度相似比设计会造成原型与模型结构底部地基土的孔压比分布规律不同。

（3）本文方法所需要的计算参数少且易于通过试验获取、计算效率高，在大型三维可液化场地土–地下结构相互作用分析中有较大的应用潜力。