水平非均匀场综合管廊地震响应及离散元模拟研究*

黄德龙 宗钟凌 刘强 岑航 汤爱平

1.江苏海洋大学土木与港海工程学院 连云港222005

2.哈尔滨工业大学土木工程学院 150090

引言

城市地下综合管廊的建设可以有效解决直埋管线所带来的“马路拉链”、“空中蜘蛛网”等“城市病”［1］。截止到2020 年，中国已有25 个试点城市开始推进管廊的建设并逐步投入使用。这些城市，有50%位于Ⅶ度设防区，有20%位于Ⅷ度设防区，因此管廊的抗震研究势在必行。

在1995 年日本阪神大地震时，地铁车站等地下结构遭遇严重毁坏，学者们才开始关注地下结构的抗震问题［2］。管廊等地下结构体系的破坏主要以外壁混凝土开裂以及内容物破坏为主。其破坏的直接原因主要是土-结相互作用（SSI），与地面建筑物所产生的惯性效应不同［3，4］。因此，场地振动、地震引起土体失效以及土-结相互作用是现阶段研究地下结构抗震的热点问题。近些年，国内外学者们对管廊等地下结构的抗震性能进行了大量动力试验和数值模拟，部分学者还进行了现场震害调查，提出了一系列理论和技术方法。

对于试验研究，振动台试验是主要采取的测试手段［5-11］，得出了均匀土体中管廊的地震响应规律以及管廊-土-地面结构（TSF）体系的响应规律，并对一致和非一致激励下土-结相互作用机理进行了深入研究。也有一些学者利用离心机试验进行地下结构抗震的研究［12-16］，探究了管廊和地震波之间的相互影响机制、土体致密化响应特点，并同样对土-结相互作用进行了深入探索。对于数值模拟研究，学者们同样主要关注均匀土或竖向成层土内的土-结相互作用、管廊及节点的地震响应、管廊抗震设计简化计算方法以及管廊-土-地面建筑之间的动力相互影响［4，17-25］。对于震害现场调查研究［1］，主要是针对地震作用下管廊的破坏形式进行总结和分析。

综合管廊可跨河床或不同地质单元，并且只要不是膨胀土，就可选用挖出土回填，造成管廊周围土体的水平非一致性。上述试验、数值模拟及现场震害调查中，鲜有学者考虑到穿越水平不均匀土体的管廊地震响应机理问题，并且多是运用有限元理论，无法精确模拟土颗粒与结构的接触，更无法模拟土颗粒迁移规律。本文基于振动台试验以及改进离散元对上述不足进行了探究。

1 试验与模拟方法

本文首先利用振动台试验得到管廊的动力响应，并将试验与改进离散元方法的结果进行对比，以验证改进离散元方法的适用性。

1.1 振动台试验

本次振动台试验在哈尔滨工业大学结构与抗震实验中心进行。试验首先预制缩尺管廊，见图1，然后采用单向振动台，并利用叠层剪切箱，将砂土与黏土分两侧放入剪切箱，管廊埋置深度为0.4m。针对El-Centro 波、Kobe 波和人工波，分别经过调幅（0.1g／0.2g／0.4g／0.8g／1.0g）进行加载。试验中分别监测上覆土体的竖向位移、管廊外壁土压力以及管廊的应变，其中在砂土和黏土内各敷设7 个土压力计，为E-1～E-14，用来测量管廊侧壁（E-1～E-3 处于砂土，E-4～E-6 处于黏土）受到的土压力增量以及管廊上（E-10～E-14）、下（E-8～E-10）的土压力增量。具体的试验设置方案参照黄德龙等［26］关于穿越水平非均匀场地管廊振动台试验。

图1 缩尺管廊试验前准备Fig.1 Preparation work before scale utility tunnel test

1.2 离散元模拟

1.改进离散元方法基本原理

离散元（PFC）的计算思想是使每个元素均满足牛顿第二定律，用中心差分的方法求解各元素的运动方程，得到研究对象的整体运动形态，而地震波对结构的作用是一个能量不断积聚的过程，和PFC的思想具有一定的差异性。因此本研究中改进的离散元方法是基于能量的思想，利用微分正交法（DQM）和Newmark-β 方法来改进离散元，提高计算精度和计算效率。

综合管廊外壁的总势能V是管廊壁的应变能U、动能K和外部土体做的功W之和，应变能可以写成：

动能可表示为：

由于土体的压力而产生的外部功可以通过下式获得：

使用Hamilton原理推导出结构的控制方程，原理如下：

为了解决工程领域中出现的初始值或者边界值问题，有很多数值方法，其中一个较好的数值方法是DQM法。与其他数值方法相比，DQM 法有以下几个优点：（1）是一种解决非线性微分方程导数近似的精确方法；（2）可以满足各种边界条件，并且需要更少的方程和编程工作；（3）准确性和收敛性都很高。基于上述优点，近年来，DQM法在结构和动力学问题分析的数值求解中越来越受欢迎。基于此方法，结构控制方程的矩阵形式［27］可以写成式（5）：

式中：KL、KNL、C、M、db和dd分别代表线刚度矩阵、非线性刚度矩阵、阻尼矩阵、质量矩阵、边界点和域点。

利用Newmark-β方法来获得结构的时程响应序列，公式（5）可以改写如下［28］：

基于DQM 及Newmark-β 的计算原理，对PFC2D程序进行二次开发，以求解式（6）非线性离散元运动方程。PFC2D程序中开发自定义求解原理，需先对模型设置特有的编号，避免与标准模型冲突。在C++程序中，依次修改源文件中CM_vep：：CM_vep（）、CM_vep：：Name（）、CM_vep：：PropNames（）、CM_vep：：ReturnProp（）、CM_vep：：AcceptProp（）、CM_vep：：SaveRestore（）等成员函数，完成源文件的修改。最后，将源文件编译成动态链接库文件（DLL 文件），可供PFC2D程序调用。

2.离散元方法的设置

图2 为管廊的PFC2D离散元模型，此模型融入管廊及其土体的Newmark动力求解方法，基于黄景崎等［29］的研究，在模型外侧施加等效荷载来实现试验中叠层剪切箱的作用效果。对于模型的构建，Potyondy等［30］指出：由于离散元对土体粒子数目无法按照真实情况来模拟，虽然土体颗粒的尺度因子对离散元是有影响的，但是当L'／d＞32，即剪切箱的最大边长与颗粒粒径之比大于32时，此时宏观响应与颗粒尺寸几乎无关，即此时颗粒的大小对管廊的响应几乎没有影响。所以模拟时为了减少粒子数目和计算所需要的时间，采用比较大的球形颗粒对土粒进行建模。模型分为两个区，接近管廊的区域颗粒小，其最大的L'／d的值为180；离管廊较远的区域为颗粒较大区，其最大的L'／d为125，两者都大于32，所以此模型可用来模拟管廊及周围土体的地震响应。本研究粒子数为14733。

图2 管廊的二维离散元模型Fig.2 2D discrete element simulation model of utility tunnel

整个布置区域与试验振动台的竖向截面尺寸一致，其中管廊和土箱分别采用wall来模拟。管廊周围相近的土体孔隙比都采用0.15，砂土采用颗分的形式定义，阻尼比为0.3，黏性土的阻尼比为0.45。对于构建砂土模型来说，颗粒与颗粒之间采用线弹性模型，如图3a。此模型需要土体的刚度系数，以土体的基床系数为依据［31］，基于本试验的土体物理性质，法向刚度和切向刚度取为3.3 ×104N／m，滑动摩擦系数取为0.4；构建黏土模型颗粒与颗粒之间同样采用线弹性模型，其中法向刚度取为2.8 ×

图3 离散元粒子与管廊接触模型Fig.3 The contact models between discrete element particles and utility tunnel

104N／m，切向刚度取为2.0 ×104N／m，滑动摩擦系数取为0.4；对于砂土与管廊之间的接触（即ball 与wall 接触）采用线性平行粘结（linearpbond）模型，如图3b。基于式（3）中的土压力，可以计算接触法向刚度和切向刚度都取为2 × 106N／m，摩擦系数取0.4，阻尼比取为0.2；对于黏土与管廊之间的接触，也同样采用线性平行粘结模型，接触法向刚度和切向刚度都取为1.0 ×106N／m，摩擦系数取0.4，阻尼比均取为0.4。

3.离散元与试验结果对比

基于DQM 计算方法，加速度、位移和应变的响应结果与试验对比如图4 所示。可以发现，加速度和位移的计算结果吻合较好。

图4 离散元模拟和试验结果对比（El-Centro，PGA ＝1.0g）Fig.4 Comparison of numerical simulation and experimental results（El-Centro，PGA ＝1.0g）

在模拟中，由于是二维模型并且土体颗粒尺寸与真实颗粒相差较大，导致颗粒之间的刚度和阻尼无法按照真实情况模拟。这种情况下，粒子的宏观运动特性可以准确实现，微应变等微观变形响应的时程曲线将很难与试验准确对应，但是其幅值均在可接受范围内。

2 试验与模拟结果讨论

2.1 竖直位移

图5 为输入PGA 从0.1g 调幅至1.0g（增幅0.2g）的三种地震波后，利用激光位移计所测得的砂土和黏土竖向位移。可以看出砂土数值增大明显，砂土区发生了剪沉现象（沉降），而黏土区数值呈减小趋势，发生了剪胀现象。随着加速度的增大，剪沉和剪胀增加量逐渐减小，即剪沉和剪胀虽然逐渐累积，但累积速度在放缓。到最后加载PGA ＝1.0g的人工波后，发现砂土依然在剪沉，黏土同样也发生了轻微剪沉现象，但是其量值很小。说明随着加载的进行，砂土的密度逐渐增大，而黏土的密度逐渐减小到一定值后还有增加的趋势。对于三种波在PGA ＝0.8g时出现了与其他地震动强度不同的情况，即砂土出现了明显的变形上浮，主要是因为三种波在PGA ＝0.4g均发生了明显的剪沉。

图5 不同土体表面的竖向位移Fig.5 Vertical displacement of different soil surfaces

土体变形将会导致土压力发生变化，而管廊所受的土压力又是管廊发生变形的直接因素。图6 为试验中砂土发生了剪沉现象，可以发现砂土流向与剪切箱接触边的中心，并逐渐向下流动。

图6 试验观测砂土沉降Fig.6 Sand settlement observation

图7 为土体发生竖向位移的改进离散元二维截面模拟情况。此程序无法实现内部管廊的转动作用，但是从图7a 中可以发现左侧砂土区发生了剪沉。土体与管廊之间出现了明显的空隙，之后砂土区左侧土体补充此空隙，但是有明显的滞后效应。这与图6 砂土区剪切箱边砂土发生的流动变化相一致。并且图7b 土颗粒的运动也可以说明砂土主要是向下沉降，而黏土是向侧壁扩展，进而导致剪胀现象的发生。

图7 土体竖向移动模拟Fig.7 Simulation of vertical movement of the soil

2.2 管廊上部土压力增量

图8 为El-Centro 波，PGA ＝0.4g 和1.0g，分别对应砂土和黏土区20cm和40cm深度处（管廊上部）的土压力增量ΔσE时程。可以发现深度越深，土压力增量越大；小震时砂土区的土压力较大，而大震时，黏土区的土压力较大；其中部分土压力增量迅速达到一个平衡阶段，然后在该位置震荡直到结束，见图8a PGA ＝0.4g 的曲线，这与Cilingir 和Madabhushi［32］发现的规律一致。Hushmand等［14］给出的解释是土体致密化会导致土压力增量增大，土压力最后可能有残余值。结合图5也可以发现，随着加载的进行，在PGA 增大的同时，砂土和黏土各自发生剪沉和剪胀现象，其残余内力在增大，表现为土压力增量在增大。图9 为利用改进的离散元程序对加载1.0g El-Centro波的土体内部不平衡力变化过程进行模拟，此不平衡力可以反映土压力的大小。由于El-Centro波在3s左右加速度达到最大峰值，因此对模拟的前3s中6 个时刻的土压力进行研究（限于篇幅，图9 中仅给出3 个代表时刻）。可以发现随着加载时刻的进行，加速度逐渐增大，砂土区管廊与土体之间的空隙逐渐增大。由于上部土体的补充，此空隙逐渐减小。可以看出土压力整体在逐渐增大。开始时，由于砂土区下部土体沉降，此区域土体内部不平衡力逐渐增大，即土压力增大。随着加速度的增大，管廊上覆土体的土压力也在增大，达到加速度峰值1.0g时，管廊上部土压力超过下部，在边界以及两种土体分界区土压力最大。

图8 不同土体的土压力增量（El-Centro 波）Fig.8 Earth pressure increments for different soils（El-Centro）

图10 表示利用改进的离散元程序对PGA ＝1.0g时三种波作用下的土压力峰值，可以发现Kobe波的土压力峰值最小。

图10 三种地震波PGA ＝1.0g 时土压力峰值模拟结果（单位：kPa）Fig.10 Simulation of peak earth pressure at 1.0g for three waves（unit：kPa）

2.3 管廊侧壁土压力增量

图11 分别为砂土区和黏土区管廊侧壁沿纵向受到的土压力增量ΔσE时程，可以发现砂土区土体对侧壁土压力增量略大。并且侧壁中间（砂土E-2 和黏土E-5）土压力增量较大，其原因是处于侧壁中间处的土体受到侧壁的约束，残余内力无法自由释放，而导致对侧壁土压力的增大。同样可以发现砂土区E-2 发生了土压力增量的负向残余值，其原因是纵向运动的过程中，砂土与管廊有垂直侧壁方向的脱离趋势。

图11 沿竖直方向分布土压力增量响应（El-Centro，PGA ＝1.0g）Fig.11 Response of earth pressure increments along vertically distributed（El-Centro，PGA ＝1.0g）

柔度比F表示结构周围土体与结构之间的相对刚度［2］。Wang［3］和Hashash等［4］都对矩形地下结构的柔度比进行了深入研究，其计算公式如下：

式中：Gm为在自由场中沿管廊高度土体的平均应变所对应的剪切模量；K 为管廊的刚度系数，可通过结构力学的方法确定，即在限制管廊底部水平运动的同时，在管廊顶部施加单位力所引起的管廊顶部横向位移的倒数；W为管廊截面的宽度；H为管廊截面的高度。

本研究的工况只有砂土和黏土两种，只涉及两个柔度比，即砂土区柔度比F1和黏土区柔度比F2。由于本研究砂土的剪切模量Gm-sand近似为7.7MPa，黏土的剪切模量Gm-clay近似为3.5MPa，而两个工况柔度比所涉及的其他参数均相同，因此F1＞F2。

图12 为砂土和黏土区沿管廊高度分布土压力增量峰值变化曲线。可以发现随着PGA 的增大，土压力增量峰值有增大的趋势，但是由于土体和管廊有分离的趋势，部分土压力增量减小。小震PGA ＝0.1g时，在侧壁中间高度处，土压力增量都要比PGA ＝0.2g时大，这是由于地震波以PGA ＝0.1g开始加载，当突然加载或波形（频率）发生改变时，土压力均将发生显著变化。大震时砂土区土压力增量要比黏土区大，即柔度比在一定范围内越大时，管廊相对于土体越柔，受到的土压力越大。还可以发现随着加速度的增大，土压力增量分布近似由底部（顶部）较大向中间位置较大变化。Hushmand 等［14］指出随着地震作用的增大，土压力增量的分布近似由三角形变为抛物线形，柔度比可以决定侧壁土压力的形状［2］。本文得出的结果与上述结论相符，但是两者不同点在于砂土区管廊下部的土压力增量峰值要大于上部，而黏土区则相反。因此砂土区要着重关注管廊下部，而黏土区上部较危险，这与2.1 节土体竖向变形有密切关系。

图12 沿管廊侧壁高度土压力增量峰值变化Fig.12 Variations of earth pressure increment peaks along the height of utility tunnel

在地震作用下，地基土体会发生沉降变形。管廊地基土体的沉降变形是破坏管廊稳定性的直接原因。土体一旦出现沉降变形就会使管廊受力发生变化。受力超过允许值，管廊将遭到破坏。由地震作用所产生的管廊侧壁水平向动土压力为［33］：

这表明挡板上的土压力水平分量py随深度增加趋于一个极值，而传统的Coulomb 土压力理论和Rankine土压力理论都认为土压力沿深度线性增加，这说明地震动力作用下的挡板土压力分布与Coulomb土压力理论和Rankine 土压力理论有较大差别。这与图12 管廊侧壁土压力呈现中间大两端小的极值曲线相接近。

分析式（8）可以发现，土压力的极值与侧压力系数K、土体重度γ、土体内摩擦角φ、挡板与土体摩擦角δ、地震加速度系数kv、kh和剪切箱-管廊侧壁净距L有关。

2.4 地表竖向位移与管廊侧壁土压力的关系

图13 为试验中地表竖向位移增量与管廊侧壁土压力增量的关系。可以发现，砂土区随着竖向沉降增量峰值的增大，侧壁土压力增量在开始阶段较大，砂土区侧壁中间E-2 的土压力增量在沉降增量为0.4mm～0.8mm 附近出现极小值。说明在发生沉降初期，侧壁土压力变化明显，侧壁中间土压力增量最大，并且存在极小值。而黏土区的剪胀与土压力增量关系存在极大值，当剪胀量在0.2mm～0.3mm 附近时，侧壁土压力增量出现极大值。

图13 管廊侧壁土压力增量峰值与地表竖向位移之间的关系Fig.13 Relationship between the vertical displacement of the ground surface and peak value of the earth pressure on the side wall of the tunnel

图14 为利用改进的离散元程序对PGA ＝1.0g时El-Centro波作用下砂土和黏土区管廊侧壁土压力峰值的模拟。可以发现砂土区土压力整体要比黏土区土压力大，并且侧壁中部土压力大于侧壁上下两侧的土压力，与图12 曲线相一致。还可以发现管廊侧壁的土压力要比管廊上下壁土压力大，管廊的破坏主要受侧壁土压力控制。

图14 El-Centro 波管廊侧壁土压力峰值模拟（PGA ＝1.0g）（单位：kPa）Fig.14 Simulation of peak earth pressure corresponding to El-Centro wave（PGA ＝1.0g）（unit：kPa）

3 结论

本文基于振动台试验以及改进的离散元数值方法对水平非均匀场地综合管廊地震响应进行探究，通过利用改进的离散元程序对管廊及其周围土体进行模拟，并与试验进行对比，得出如下结论：

1.基于能量思想，并利用微分正交法以及Newmark-β方法得到模型体系的运动控制方程，以此改进离散元程序，建立PFC2D管廊模型。模拟结果与试验结果的误差均在20%以内，即二次开发的离散元数值模拟方法是准确的；

2.砂土区土体发生了剪沉，黏土区发生了剪胀现象，并且伴随着土体密度的变化，管廊发生了轻微转动，砂土区剪沉现象的发生是为了补充管廊下部砂土与管廊的空隙；

3.土体致密化会导致土压力增大，土压力增量可能会出现残余值，土压力的大小与所输入PGA有关，不均匀土体对土压力的影响较大。管廊侧壁的土压力增量呈中间大两端小的极值形，并且在一定范围内，柔度比越大，管廊所受到的土压力越大，管廊侧壁变形也越大；

4.砂土区在发生沉降初期，侧壁土压力变化明显，且之后存在极小值；黏土区发生剪胀，随着剪胀量的减小，土压力增量存在极大值。