地震作用下可液化场地管道上浮动力响应及影响因素分析

顾世祥，张继虎，梅 伟，杨 帆，单夺权，童保林，杜文琪

（1.云南省水利水电勘测设计研究院，云南 昆明 650021；2.武汉大学水资源与水电工程科学国家重点实验室，工程风险与防灾研究所，湖北 武汉 430072）

引言

输水管道因具有安全可靠、施工便捷以及土地资源占用较少等优点，已广泛应用于长距离输水工程。近年来，地震作用下场地液化造成的管道上浮破坏现象频繁发生，在2010 年智利Maule 地震、2011 年日本Great East Japan 大地震和2011 年新西兰Christchurch 地震中，场地液化均导致了管道的大面积上浮破坏现象［1］，其中前二者均生成了大量的长持时地震波。目前我国水利建设正逐步面向西南高地震烈度区，如滇中引水二期工程玉溪段地震烈度高且断裂带发育，近场地震效应显著。因此有必要进一步开展地震作用下管道上浮动力反应的相关研究，从而提升长距离输水管道系统的抗震性能。

可液化场地管道上浮破坏的研究方法主要可以概括为三类：（1）物理试验，如邹德高［2］通过振动台试验指出砂土液化后土体抗剪强度减小是导致管道上浮发生的直接原因，CASTIGLIA 等［3］利用振动台试验分析了管道自重和地震幅值、时程对上浮反应的影响等；（2）静力有限元数值模拟，多以折减土弹簧刚度的方式来模拟液化，并施加荷载边界模拟管道上浮。如林均歧等［4］采用非线性增量有限单元法，探讨了液化区长度和轴向力等对管道上浮的影响；陈艳华等［5］基于土弹簧模型分析了管道上浮的影响因素。静力分析方法虽然简单易用，但不能考虑地震特性及孔压增长、扩散、消散等性质，且所得结果偏大［6］；（3）动力时程数值模拟，可以反映土体孔压和管道上浮的时程变化，但对土体本构模型要求较高。如SAEEDZADEH 等［7］运用PLAXIS 软件开展了二维数值模拟，分析了土体相对密度、地下水位和管道埋深等对管道上浮位移的影响，但其选用的液化土体本构模型较为简单。地震作用下可液化场地管道的上浮动力反应，与管土材料和地震动特性密切相关，对于数值模拟计算而言，还涉及到土体本构是否可以正确描述液化后土体性质这一问题。因此，本文注意到现有研究中存在的两个问题：第一，对于液化土体变形和孔压变化的描述不够合理，如静力有限元分析方法甚至无法描述土体孔压变化，所得结果准确度略有不足；第二，地震动持时和近场脉冲特性对于管道上浮以及管道破坏模式的影响尚未得到充分研究。

针对以上两个问题，本文基于美国太平洋地震工程研究中心（Pacific Earthquake Engineering Research Center，PEER）主导开发的OpenSees开源有限元软件［8］，采用PDMY液化土体本构模型，建立可液化场地埋地管道二维分析模型。基于目标反应谱、谱匹配等方法从NGA-West2地震动数据库选取地震波。首先验证了数值模型方法的有效性，然后系统分析了管土特性、长短持时与脉冲非脉冲地震动特性对管道上浮反应的影响规律。

1 地震动选取

1.1 基于目标反应谱选取地震动

为研究管土特性对管道上浮的影响，基于云南省滇中引水二期工程玉溪段某场地设计反应谱，选取具有一致谱形的地震动记录。该场地抗震设防烈度为8度，设计基本地震加速度值0.3 g，Ⅳ类场地。根据反应谱在NGA-West2数据库选择和调幅一定数量的地震动记录，为充分考虑地震动不确定性，共选取15条地震动，所选地震动加速度反应谱及其震级-断层距分布如图1所示。

图1 选取地震动的加速度反应谱与震级-断层距分布Fig.1 Acceleration response spectra and distribution of moment magnitude versus rupture distance

1.2 长短持时地震动的选取

持时作为地震动主要特性之一，其对框架结构、岩土边坡、重力坝等结构响应和破坏的影响得到了研究人员较为广泛的关注［9－12］，但其对管道上浮破坏的影响研究较少。本文采用谱匹配方法，选取长短持时地震动，统一调幅至0.3 g，分析地震动持时对管道上浮反应的影响。对于长短持时界限本文以重要持时Ds5-75=25 s为阀值［9］。基于Arias强度的Ds5-75重要持时，以地震动Arias强度达到75%对应时刻与5%对应时刻之差作为地震动持时指标。Ds5－75重要持时不随地震动缩放变化，且能够较好地反映地震动能量特性，在工程中应用较为广泛［13］。

谱匹配方法是根据选定的长持时地震动记录，对一系列短持时地震动进行缩放，选取谱形与长持时地震动谱形几何均值最小的短持时地震动。首先将选定的长持时地震动反应谱等时间间隔离散为L1，L2，L3，…，Ln，均值为。然后将从NGA-West2数据库中选出的第j条持时小于25 s的短持时地震动反应谱，对应离散为S1j，S2j，S3j，…，Snj，均值为-Sj，缩放系数k=，两谱形方差SSE计算为：

式中：n为所考虑的反应谱周期数，缩放系数限制为k＜5，以避免低强度地震动被过度缩放［14］。取SSE方差最小的第j个短持时地震动为匹配地震动。所选长持时地震动为2008 年汶川邛崃地震记录，而短持时地震动为2002年Baja地震记录（NGA-West2数据库中RSN2005）。图2为谱匹配后长短持时地震动及其反应谱，可以看出两者谱形基本吻合。

图2 谱匹配后长短持时地震动时程及其反应谱Fig.2 Acceleration time histories of a long-duration and short-duration record pair and response spectra

1.3 脉冲地震动选取

近断层脉冲地震动表现为强方向性，速度时程中含有大幅度脉冲，其对于地面建筑以及土石坝的安全具有严重威胁［15-16］。本文基于SHAHI［17］提出的基于连续小波变换的脉冲地震动识别算法来选取脉冲波，该方法以4阶Daubechies（Db4）为母小波函数，小波函数可表示为：

式中：Φs,l和φ分别小波基函数和母小波；s为伸缩因子；l为平移因子。通过对实际地震动做Φs,l的连续小波变换，得到不同位置和尺度对应的小波系数：

式中：c(s,l)为小波系数；f(t)为地震信号。基于BAKER［18］对于脉冲地震动的判定方法，即可选取一组合适近场脉冲波。图3（a）和图3（b）分别为脉冲波与非脉冲波，其中脉冲波选自1999 年台湾Chi-Chi 地震（RSN1513），非脉冲波选自1995年日本Kobe地震（RSN1120）。两地震动震级（＞6）和断层距（＜20 km）虽然都在相同范围，但脉冲地震动表现出明显的强速度脉冲效应。

2 数值模型的建立与验证

2.1 土体本构模型

本文选取OpenSees 中nDMaterial 材料库的Pressure Depend Multi Yield（PDMY）和Pressure Independ Multi Yield（PIMY）两种材料分别模拟可液化砂土及粉质粘土的动力学特性。PDMY 模型是一种多屈服面弹塑性模型，可以模拟砂土非线性滞回及剪胀特性，能够有效模拟砂土的液化行为。该模型由YANG 等［19］提出并基于OpenSees 实现了数值开发，得到了较为广泛的应用。PIMY 模型是一种弹塑性本构模型，服从Von Mises屈服准则，其塑性行为只表现在偏应力-应变反应上，体积应力-应变为线弹性且独立于偏分量，可有效模拟粘土等材料力学特性。

2.2 有限元模型

OpenSees 没有前后处理界面，一般采用Matlab 或GID 用于前处理和计算结果的可视化。本文通过GID作为前处理软件进行网格划分，然后将单元和节点信息导入OpenSees建立有限元模型。土体单元采用四节点u-p 格式流固耦合单元，每个单元节点具有三个自由度，包括两个位移自由度和一个孔压自由度。u-p 格式动力流固耦合方程虽然忽略了孔隙流体和固体之间的加速度差，但能够满足一般砂土的地震分析需求［20］。管道单元采用四节点实体单元，每个单元节点包括两个位移自由度。管道模型采用弹性本构，材质取球墨铸铁，密度7 200 kg/m3，弹性模量156 GPa，泊松比0.3。图4 所示管道埋置在相对密度（Dr）为30%的松砂层中，上覆1 m 厚粉质粘土层，底层为10 m 密砂层（Dr=70%）。土体参数结合云南省滇中引水二期工程玉溪段某场地工程地质勘察报告和模型推荐值综合确定，具体参数信息见表1。

图4 计算模型（单位：m）Fig.4 Computational model

表1 土体材料参数Table 1 Soil parameters

2.3 边界条件与动力求解

鉴于管道为薄壁结构，为保证计算收敛性，绑定土体与管道节点的位移自由度，模型左右两侧同高度节点绑定位移自由度形成周期性边界，用以模拟无限地基条件［21］。地下水平面设与地表平齐，地表为自由排水边界，其余边界均为不透水边界。

求解过程分为重力分析、弹塑性分析和动力分析三部分。重力分析获得初始应力场和孔压场，然后更新材料至塑性进行弹塑性分析，最后进行动力分析。动力计算采用基底一致激励方式从模型底部施加水平加速度时程；采用Newmark 隐式求解法，取γ=0.5，β=0.25；为反映塑性土体的能量耗散，消除高频噪声，施加5%的瑞利阻尼。为兼顾计算收敛性及计算精度，动力计算时间步长取为0.005 s。

2.4 模型有效性验证

综上所述，管道上浮动力反应的准确描述很大程度上依赖于液化土体本构模型。为验证有限元模型的准确性，对CHIAN等［22］所做圆形隧道上浮离心机试验进行了数值模拟。试验原型为埋置在相对密度为45%的Hostun 砂中的圆形隧道，隧道直径5 m，埋深7.5 m，底部施加幅值为0.1 g、频率为0.75 Hz 的27 s 正弦波。具体试验流程及模型信息等可参考文献［22－23］，限于篇幅本文不再说明。图5 显示了本文数值模拟与离心机试验以及CHIAN 等［22］基于Flac2D 软件进行数值模拟的结果对比，可以看出本文有限元模型能够较好反映地下结构的上浮反应和可液化土体的孔压变化特性。

图5 数值模拟与离心机试验验证Fig.5 Comparison of numerical results and the centrifuge test

3 土体和管道特性的影响

3.1 土体相对密度

土体相对密度对于液化触发具有重要影响，一般来说，相对密度越大，地震作用下土体颗粒剪缩量越少，有效应力减小量也相应较小，液化越不易触发。为考虑土体相对密度的影响，将中间土层相对密度分别设置为30%（松砂）、50%（中密砂）和70%（密砂），依次输入所选15条地震波进行非线性时程分析。不同相对密度条件下的土性参数信息见表1。图6 给出了NGA-West2 数据库编号RSN6911 地震波作用下埋深4 m和10 m处的孔压时程图，可以看出：虽然在埋深4m处区域液化均得到了触发，但相对密度越大，土体液化触发越迟、液化持续时间越短并且液化范围越小。由图7（a）可以看出：相对密度越大，管道上浮量及破坏风险越小。图7（b）中所示散点为15 条地震动作用下管道上浮位移均值，可见：Dr=30%上浮位移均值几乎为Dr=70%上浮位移均值的4倍。根据拟合公式，发现管道上浮位移与土体相对密度呈幂函数关系。有鉴于此，在管线铺设过程中，应保证回填土层及上覆土层具有良好的压实度，以降低管道上浮破坏风险。

图6 相对密度对孔压影响Fig.6 Effect of density ratio of soil on pore pressure

图7 土体相对密度对管道上浮的影响Fig.7 Effect of density ratio of soil on pipe uplift

3.2 管径尺寸

在考虑场地完全液化的拟静力分析方法以及地下结构抗浮稳定分析方法中，往往假设单位长度管道受到的上浮力满足阿基米德定律［4］：

式中：F为单位长度管道所受到的上浮力；γsat为可液化土体饱和重度；D为管道直径。虽然完全液化假设过于简化，但是仍然可以看出管径对于上浮反应的重要影响。本文考虑水利水电工程中常用的四种管径：1.0 m、1.3 m、1.6 m 和2.0 m，以埋深2 m，土体相对密度Dr=30%为例来分析管径的影响。由图8（a）可以看出：管径越大，管道发生上浮破坏风险越大。图8（b）中所示散点为各管径下上浮位移均值，可以看出：2 m 管径的上浮位移均值几乎是1 m 管径上浮位移均值的4 倍。根据拟合发现管道上浮位移和管径大体呈二次函数关系，这与CHIAN等［24］试验结果认为圆形结构上浮位移与直径呈平方关系近似。考虑到输水管道管径一般与水力计算及工程成本相关，为减轻上浮破坏风险，在液化场地管道设计中应尽可能兼顾管径尺寸对于工程效益和管道上浮破坏风险的影响。

图8 管径对管道上浮的影响Fig.8 Effect of pipe diameter on pipe uplift

3.3 埋深

液化触发与土层深度密切相关，随深度增加，地震液化风险一般也随之降低。KOSEKI等［25］对日本1993年Kushiro-Oki 和1994 年Hokkaido-Toho-Oki 地震中管线上浮破坏现象进行了调查，发现大部分由上浮导致破坏的管线埋深在1m 左右。为了研究管道埋深对管道上浮反应的影响，以管道管径为2 m，中间层砂土相对密度为50%为例，分别研究1.0 m、2.0 m、3.0 m 和4.0 m 四种埋深条件对管道上浮反应的影响。图9（a）显示了不同埋深下上浮位移累积概率分布，可以看出随着埋深的增加，管道发生上浮破坏的风险显著降低。这是因为随着管道上覆土层厚度的加大，管道周围土体液化程度会随之降低（图6中已有表明），因此上浮破坏现象得到了有效抑制。图9（b）显示了管道上浮位移与埋深大体上呈线性函数关系，这与CHIAN 等［24］试验观察结果是类似的。

图9 管道埋深对管道上浮的影响Fig.9 Effect of burial depth on pipe uplift

4 地震动特性的影响

4.1 地震动持时的影响

将1.2 节选出的长短持时地震动分别输入模型进行时程分析。液化分析中超静孔隙水压力比ru常被用来进行液化判断：

式中：Δu为超静孔隙水压力；σ'v为初始竖向有效应力。本文以ru达到1作为液化状态，监测与管底等高程的A点孔压时程。

图10（a）和图10（b）分别为短持时和长持时作用下A 点超孔压比时程和管道顶部结点上浮位移时程。可以看出长持时地震动作用下管线上浮位移明显大于短持时地震动作用下管线上浮位移，峰值位移可达短持时对应位移的2 倍左右。主要原因在于较长持时意味着地震动有效作用时间增加，超静孔压长时间维持在较高水平。图10 中长持时地震动作用下液化时间几乎是短持时作用下液化时间的2.5 倍，从而导致管道上浮位移大大增加。因此在液化场地的管道抗震设计中，应考虑地震动持时特性对其上浮反应的重要影响。

图10 孔压比和上浮位移时程曲线Fig.10 Evolutions of ru for the monitor node A and uplift displacement of the pipe versus time

4.2 速度脉冲的影响

为消除1.3节选取的脉冲波与非脉冲波之间地震动幅值（PGA）差异，本文将其统一调幅至0.3 g，以分析地震动脉冲特性对液化场地管道动力响应的影响。由于地震动的速度脉冲效应往往能够使结构发生较大的方向性变形，因此同时监测了管道顶部结点的上浮和横向位移。

图11（a）和图11（b）分别为脉冲波和非脉冲波作用下管道上浮、横向位移时程和监测点A的超孔压比时程，可以看出脉冲波作用下管道上浮破坏风险相对较小，而横向破坏风险却大大增加。其原因在于脉冲波一般持时较短、周期较长，如前所述地震动持时对液化持续时间具有重要影响，因此脉冲波孔压消散较快，导致上浮位移较小，文献［26］也指出了脉冲波作用下的自由场地孔压消散特性。而脉冲波速度峰值较大，表现为较大的瞬时能量，导致管道产生较大的横向位移。因此脉冲地震动作用下，液化场地管道具有上浮破坏和横向破坏两种破坏模式，而且横向破坏风险高于上浮破坏风险。应当指出：由于管道上浮属于多力耦合问题，除孔压之外，由于“棘轮效应”所导致的显著横向位移也会影响到管道的上浮反应程度［27－28］，两者之间相互作用仍有待于进一步研究。

图11 孔压比、管道上浮位移和横向位移时程Fig.11 Evolutions of ru for the monitor node A and uplift and lateral displacement of the pipe versus time

5 结论

本文基于OpenSees有限元软件，采用多屈服面砂土液化PDMY 模型，研究了土体相对密度、管径和管道埋深等管土特性，地震动持时和近断层脉冲等地震特性对于液化场地管道上浮破坏的影响，主要结论如下：

（1）土体相对密度越小，管径越大，管道埋深越浅，管道上浮破坏的风险越大，在此基础上，分别提出了土体相对密度、管径和管道埋深对管道上浮位移的经验拟合公式。因此在管径和埋深受限于水力计算和经济环境等限制而难以得到有效调整时，须保证管道回填土体的压实度以减小上浮风险。

（2）长持时地震动作用下，超静孔压消散缓慢，地震动作用时间较长，导致管道上浮位移大于短持时地震动作用。在地震动幅值相同（PGA=0.3 g）和反应谱形状基本无差别的情况下，长持时地震动作用下管道上浮位移是短持时地震动作用下管道上浮位移的2 倍左右，液化持续时间是短持时地震动作用下液化持续时间的2.5倍左右。

（3）近断层脉冲地震动作用下管道上浮破坏和横向破坏两种破坏模式同时存在。由于脉冲波持时较短以及峰值速度较大，液化场地孔压消散较快，导致管线上浮破坏风险小于非脉冲地震动；而由于强速度脉冲效应，其对管道横向破坏的风险大于非脉冲地震动。

（4）鉴于地震动本身的复杂性与不确定性，地震动特性对管道上浮破坏的影响仍有待进一步深入研究，包括地震动持时对管道上浮影响的机理与具体量化指标，以及脉冲地震动作用下管道上浮破坏和横向破坏的耦合作用等。