基于饱和多孔介质的复杂自由场地震响应分析

董 云，楼梦麟

（同济大学 土木工程防灾国家重点实验室，上海 200092）

场地地震液化反应分析及其有效预测对确定地基动力稳定性和进行工程结构抗震安全性评价具有重要意义.考虑到对于实际岩土材料这种常见的多孔多相介质，用传统的单相介质理论已经不能满足实际地震工程的要求，因此，研究多孔介质的动力瞬态响应问题十分重要［1-3］.

对于饱和土体，Biot［4］首先讨论了流体饱和多孔介质中的波传播问题，建立了考虑固相和液相相互作用的动力控制方程，其理论和结果被广泛引用并成为之后众多模型的参考标准.Bowen［5］基于连续介质力学的混合物理论发展和完善了土力学耦合方程体系.Zienkiewicz等［6-7］基于 Terzaghi有效应力的概念，从土壤整体的平衡方程和流体的运动方程出发，建立了饱和土的广义Biot理论模型，同时导出了适用于不同动力环境下的各种形式的控制方程.近年来，Kraaijeveld和邸元等用单位分解有限元法和有限元-有限体积混合方法对饱和多孔介质的动力问题也进行了深入研究［8-9］.

本文采用不可压缩饱和多孔介质的广义Biot理论模型，运用流体饱和两相多孔介质振动问题的有限元方程，并以港珠澳工程沉管隧道段自由场为例，研究了基于多孔介质理论的复杂实际场地的地震瞬态响应问题，并进行了大量的数值分析，得出了一些有益的结论.

1 饱和多孔介质的动力方程

根据混合物理论的空间平均方法和Terzaghi有效应力原理可以得到以增量形式表示的饱和多孔介质的本构关系［10］：

土体运动方程为

流体运动方程为

对于地震荷载这样非高频的情况，可以采用以固相的位移和孔隙水压力为基本变量的所谓u-pf形式［11］，则有

用有限元的离散化方法以形函数表示基本变量，采用Galerkin的关于权函数取为形函数的概念，最终得到以矩阵表示的饱和多孔介质的u-pf形式的动力方程［12-13］.在时域范围离散后，动力耦合方程即可采用Newmark法或 Wilson-θ法求解，如式（6），（7）所示.

式中：M，C，K分别为体系的质量矩阵、阻尼矩阵和刚度矩阵；Q为固液两相耦合矩阵；S为饱和介质Biot系数矩阵；H为孔压矩阵和分别为固相和液相的荷载向量列阵.

2 复杂自由场数值模型

本文以港珠澳工程沉管隧道段实际自由场为例，研究饱和多孔介质在地震荷载作用下的地震响应及液化影响.一般认为土体的抗液化强度和振动孔隙水压力与振动前的初始应力状态有关，为此首先要计算初始应力场，得到各个单元的应力水平.图1为实际自由场详勘纵断面图，其勘探深度范围内覆盖层厚约120m.主要为全新世及晚更新世海相、陆相沉积物，下伏基岩主要为燕山期花岗岩和震旦系片麻状混合花岗岩、混合花岗岩，自上而下可主要分为淤泥质土、粉质黏土、砂土和花岗岩.从隧道场地揭示的水文地质条件分布看，花岗岩裂隙较发育，裂隙水与上部覆盖层主要含水层内的地下水沟通良好，在地震过程中不考虑液化的影响，在分析中以单相介质模拟，不考虑孔隙水压力，其他土层用本文所述饱和多孔介质模拟.图1中位置A和位置B分别表示场地土层厚度最大处和最小处，位置C表示地形相对平缓处.在地震荷载输入过程中，由于持续时间很短，孔隙水来不及排出，故在数值模拟中地表采用不排水边界.

图1 自由场详勘纵断面图Fig.1 Longitudinal profile of free field

为了考虑多孔介质的非线性，在静力初始场的分析中采用Duncan Chang模型.而在动力分析的过程中，国内外动力试验已经证明，土体的动应力和动应变具有非线性性质，其剪切模量G和阻尼系数随着剪应变幅度而变，为此在采用Hardin Drnevich模型时考虑这一非线性关系.具体的土层详勘数据、剪切模量和阻尼比随着剪切应变变化规律曲线参见参考文献［14］.

3 饱和多孔介质的动力场分析

对于复杂场地的地震作用分析，直接从基岩面输入地震动，能更准确地反映场地的地形、地貌对地面运动的影响.目前已有的土层特性勘探技术己经具备了直接进行土层地震反应分析的条件.图2为基岩输入加速度时程曲线，该时程曲线为工程场地地震安全评价［15］专项研究中所得出的120年超越概率为10%的土层下卧基岩自由场人工地震波.

图2 加速度时程曲线和Fourier谱Fig.2 Acceleration time-history curves and their Fourier spectra

图3 地表加速度对比结果Fig.3 Comparison results of PGA

3.1 地表加速度的对比分析

从基岩面输入加速度分别分析了场地土层单相介质自由场和固液两相介质自由场的地震响应.图3是两种情况下不同位置处自由场沿深度的加速度对比结果.从图3可以看出，不同位置处，固液两相介质自由场地表加速度峰值沿深度均小于单相介质自由场的计算结果，位置A处相差最大，达到13.5%；位置B和位置C处分别为3.75%和1.61%.主要原因可能是由于考虑了振动孔隙水压力的影响，计算体系更具柔性而引起的.地表以下，位置A处多相结果基本上均小于单相结果，而位置B和C结果基本接近.从图3还可以看出，实际地形对计算结果有较大影响，譬如位置A处砂土层厚度迅速增大，其加速度结果明显较其他位置处大，单相场和两相场的差别也较明显.可见对于饱和多孔介质复杂自由场来说，地表加速度峰值由于考虑了振动孔隙水压力的影响会小于传统单相场的计算结果，而两者相差的幅度可能主要取决于场地的地形条件.

3.2 两相场孔压比的分析

场地地震研究中，孔压比是评价地震液化的一个重要指标.图4为根据公式（7）计算的实际自由场不同位置处孔压计算结果，限于篇幅，仅取基岩以上3层土中点位置处的孔压曲线.由图可知，在基岩地震波输入下，不同位置处自上而下峰值孔压比均逐渐减小，说明上部土层更容易发生地震液化，这与实际震害调查结果相似.主要原因在于：一般情况下，土层中有效应力随埋深的增大而增大，土层的地震响应又随埋深的增大而减小，而根据Terzaghi有效应力原理，孔压比与土层中有效应力成反比，故此自上而下土层中峰值孔压比随埋深的增大而减小.

图4中淤泥质土的孔压比峰值在0.40左右，粉质黏土和砂土的更小，分别在0.20和0.05左右，根据一般认识，该自由场没有液化，这也与文献［12］中的勘察结果一致.但是不同位置处孔压比峰值的差别还是比较明显的，如位置A处淤泥质土和粉质黏土的孔压比明显大于位置B的孔压比，而位置C处的孔压比则介于两者之间.为此进一步研究场地地形对其计算结果的影响程度，分别以图1中A，B，C断面处的土层参数模拟3个水平成层自由场进行分析并和实际场地的计算结果进行对比，表1给出不同位置处不同类型土中点处的孔压比对比结果，两相表示多孔饱和介质场地孔压比，单相表示传统单相介质水平自由场的计算结果，以百分比来表示传统水平自由场和实际场地孔压比的差异程度.

图4 孔压比时程曲线的对比结果Fig.4 Comparison results of time-history curves of pore water pressure

由表1可以看出，不同位置处不同土类之间孔压比的差异程度明显.在场地动力分析中，通常假定各类土是水平成层的，忽略了地形起伏的影响，这种简化对大范围的自由场研究是可取的，譬如对某一地区的地震区划等，但是对于具体的工程而言，地形的变化却不能忽略，尤其对于地下长大工程来说.如本文所示，沿着土层深度增加实际场地和水平成层场地的差异增大，这是因为中上部土层近乎水平，而砂土层起伏较为明显，地形的变化对砂土计算结果的影响较为显著.由图1可以看出位置AB之间距离670m，位置AC之间距离1400m，位置C更远离位置A.由表1可见，与位置A的距离不同，孔压比的大小也不同.位置C相对于单相水平成层场的差异程度要小于位置B，可见地形的影响有一定的范围，离地形变化复杂处越远，受其影响越小.

图4中，随着基岩地震荷载的输入，孔压迅速增长，到18s左右达到峰值并趋于水平，这一变化规律是因为数值模拟中只针对孔压动力增长过程，采用的是每经过一个应力循环进行孔压增量计算的办法，所以无法模拟峰值后因孔隙水排出和土层流动引起的孔压消散过程.

3.3 渗透系数的影响

在饱和多孔介质的运动方程中，流体的渗透系数是一个关键参数，为此本文研究在不同渗透系数下场地的地震响应.在不改变其他参数的情况下，仅改变粉质黏土的渗透系数，将其从原来的3×10-8m·s-1变为3×10-6m·s-1，对这两种情况进行研究.输入地震荷载之后，发现地表加速度几乎没有变化.图5为两种情况下位置A处各层土中点位置处的孔压比曲线对比.

从图5可以看出，在不同渗透系数下，不同土类的孔压比曲线也发生了变化.粉质黏土的孔压比峰值从0.27减小到0.21，淤泥质土稍微有所变小，而砂土几乎没有变化.而且还可以明显看出粉质黏土达到孔压比峰值的时间提前了，淤泥质土也有同样的变化.可见渗透系数增大，孔压比可能会变小，这主要是因为土体的透水性增大，相对于渗透系数较小的情况来说，每一应力循环步孔压增量变小，导致峰值孔压比减小，且孔压比曲线迅速趋于水平.

3.4 不同震级下的地震响应

震害资料显示，大震时土体液化情况增多，本文所输入地震荷载峰值约为0.15g，按国内情况来讲属于“中震”，为此将本文地震荷载按比例放大，使其峰值达到0.40g，研究在这两种情况下饱和两相介质自由场的地震液化情况.图6为两种情况下（0.15g，0.40g）位置 A处各层土中点位置处的孔压比曲线。

表1 实际场地与水平成层场孔压比对比结果Tab.1 Comparison results of pore pressure with horizontal layered sites

图5 不同渗透系数下的孔压比对比结果Fig.5 Comparison results of pore water pressure with different permeability coefficients

图6 不同震级下的孔压比对比结果Fig.6 Comparison results of pore water pressure under different earthquake magnitudes

图6显示在“大震”作用下，自由场上部淤泥质土孔压比达0.8，可以认为已接近液化状态，而中部的粉质黏土的孔压比峰值大于0.5，且超过“中震”情况下上部淤泥质土的孔压比，也存在潜在的液化风险.可见对于同一场地，不同震级下的孔压比存在较大差异，需要根据场地当地具体情况选择合适的荷载级别进行地基动力稳定性和工程结构抗震安全性评价.

4 结论

本文采用不可压缩饱和多孔介质的广义Biot理论模型，导出流体饱和两相多孔介质振动问题的有限元方程，直接在每一个应力循环步中进行孔压增量计算，并以港珠澳工程沉管隧道段自由场为例，分析并研究了不同情况下实际自由场的地震响应，得到以下结论：

（1）在基岩地震输入下，饱和两相多孔介质自由场地表加速度峰值相对于传统的单相介质自由场来说有所降低，且不同地形下两者的差异程度不同，这主要是由于考虑了振动孔隙水压力的影响，计算体系更具柔性而引起的.

（2）对于实际工程场地而言，多数情况下并不是理想化的水平成层体系，本文的计算结果显示，地形地貌等参数对饱和两相多孔介质自由场的地震响应影响显著.

（3）多孔介质的渗透系数增大，峰值孔压比可能会减小，且达到峰值的时间提前了，这主要是因为透水性增大，每一应力循环步孔压增量变小，导致最终峰值孔压比减小.

（4）震害资料和计算结果显示，不同震级下自由场液化的风险不同，建议场地地震液化评估时根据实际情况选取不同输入荷载.

（5）本文孔压动力增长过程采用的是每经过一个应力循环进行孔压增量计算的办法，无法模拟峰值后因孔隙水排出和土层流动引起的孔压消散过程，所以孔压比曲线达到峰值后趋于水平.

［1］Vessia G，Venisti N.Liquefaction damage potential for seismic hazard evaluation in urbanized areas［J］.Soil Dynamics and Earthquake Engineering，2013，31（8）：1094.

［2］Susumu T，Hirofumi T，Tadashi H，et al.Liquefaction at a housing site reclaimed using well-graded gravelly soil ［J］.Japanese Geotechnical Journal，2010，5（2）：377.

［3］Sarah D，Soebowo E.Liquefaction due to the 2006 Yogyakarta Earthquake：field occurrence and geotechnical analysis ［J］.Procedia Earth and Planetary Science，2013，6：383.

［4］Biot M A.Theory of propagation of elastic waves in a fluidsaturated porous solid I low-frequency range［J］.The Journal of the Acoustical Society of America，1956，28（2）：168.

［5］Bowen R M.Incompressible porous media models by use of the theory of mixtures［J］.International Journal of Engineering Science，1980，18（9）：1129.

［6］Zienkiewicz O C，Shiomi T.Dynamic behaviour of saturated porous media：the generalized Biot formulation and its numerical solution［J］.International Journal for Numerical and Analytical Methods in Geomechanics，1984，8（1）：71.

［7］Edelman I，Wilmanski K.Asymptotic analysis of surface waves at vacuum／porous medium and liquid／porous medium interfaces［J］.Continuum Mechanics and Thermodynamics，2002，14（1）：25.

［8］Kraaijeveld，Huyghe，Remmers，et al.Two-dimensional mode I crack propagation in saturated ionized porous media using partition of unity finite elements［J］.Applied Mechnics，2013，80（2）：35.

［9］邸元，唐小微.饱和多孔介质大变形分析的一种有限元-有限体积混合方法［J］.计算力学学报，2008，25（4）：483.DI Yuan，TANG Xiaowei.A finite element-finite volume hybrid method for large deformation analysis of porous media［J］.Chinese Journal of Computational Mechanics，2008，25（4）：483.

［10］Zienkiewicz O C，Taylor R L.The finite element method for solid and structural mechanics［M］.Burlington：Butterworth-Heinemann，2005.

［11］Johari A， Khodaparast A R. Modelling of probability liquefaction based on standard penetration tests using the jointly distributed random variables method［J］.Engineering Geology，2013，158：1.

［12］Takashi K，Junichi K，Takeshi S.Relationship between limitingshear strain and reduction of shear moduli due to liquefaction in large strain torsional shear tests［J］.Soil Dynamics and Earthquake Engineering，2013，49：122.

［13］楼梦麟，董云.大型直流电压发生器塔架结构地震反应分析［J］.防灾减灾工程学报，2013，1（2）：35.LOU Menglin，DONG Yun.Seismic response analysis of largescale DC voltage generator tower structure［J］.Journal of Disaster Prevention and Mitigation Engineer，2013，1（2）：35.

［14］中国地震局地壳应力研究所.港珠澳大桥工程场地设计地震动参数研究报告［R］.上海：中国地震局地壳应力研究所，2009.Institute of Crustal Stress of China Seismological Bureau.Research report for design parameters of ground motion of engineering site of Hong Kong-Zhuhai-Macao Bridge［R］.Shanghai：Institute of Crustal Stress of China Seismological Bureau，2009.

［15］中国地震局.GB17741—2005工程场地地震安全性评价［S］.北京：中国建筑工业出版社，2005.China Seismological Bureau.Evaluation of seismic safety for engineering sites［S］.Beijing：China Architecture ＆Building Press，2005.