随着我国城市化进程的快速推进及一带一路和西部大开发战略的实施，涌现出大量长大地下结构工程．综合管廊作为长大线状地下结构，承担着收纳和保护各类管道及缆线的作用，属城市生命线工程，在灾害来临时保证其安全、完整对城市功能运转与灾后修复具有重大意义．数值模拟作为一种有效分析手段常被用于长大地下结构地震反应特征和破坏机理的研究．文献[1]提出了一种使隧道接头在张拉状态得到释放的改进设计，建立了地下隧道土-结构相互作用模型，模拟了隧道在地震作用下的地震动响应．文献[2]对浅层隧道开展有限元模拟，用全新的土体和隧道衬砌本构模型，更真实地再现其在动力作用下的力学特性．文献[3]针对上海市供水盾构隧道进行真三维数值分析，对隧道的抗震性能进行了深入研究，给出了隧道整体地震安全性评价．文献[4]基于ANSYS有限元软件，通过全面的参数分析探讨紧邻多孔交叠盾构隧道体系的地震响应规律．文献[5-6]采用文献[7-9]提出的修正Davidenkov本构模型，以武汉三阳路和苏埃海底盾构隧道为工程背景，分析了盾构隧道纵向管节张开量及内力响应．文献[10]采用总应力法，通过调整沿隧道纵向管环间的连接螺栓数量，分析了连接螺栓数量沿隧道纵向管环的非一致分布对纵向环间张开量的影响及抗震机理．综上可知，虽然国内外学者采用数值分析方法已对长大地下结构地震反应特性进行了一些有益的探讨．然而，介于本构模型的局限性及二次开发的难度较大，目前尚未有采用有效应力本构开展针对可液化场地长大地下结构地震反应特性的研究．基于文献[11-12]提出的剪切-体应变耦合的孔压增量模型，在ABAQUS有限元软件中独立进行了二次开发，以南京江北新区综合管廊为实际工程背景，基于广义反应位移法研究可液化场地综合管廊非线性地震反应特性．1 数值分析模型1.1　工程概况选取南京江北新区地下综合管廊(浦乌路段)为工程实例，工程场址位于长江沿岸．综合管廊采用预制单舱的形式，外轮廓尺寸为5.3 m(高)×4.15 m(宽)，埋深约为15 m，图1为管廊横断面尺寸及实物图．管廊的纵向长度为1 500 m，每个预制节段长度为1.5 m，共计1 000个预制管廊节段．图2为管廊所处工程场地的土层分布示意图，图中：D为场地距离；H为高程．由图2可知：该管廊结构穿越土体复杂，包括粉质黏土、碎石粉质黏土、粉细砂和粉砂．10.13245/j.hust.220811.F001图1管廊横断面尺寸及实物图10.13245/j.hust.220811.F002图2管廊所处工程场地的土层分布示意图1.2　地震动选取《工程场地地震安全性评价》(GB 17741—2005)[13]中规定：针对工程场地的地震安全性评价中考虑的潜在震源区域范围应不小于工程场地外围150 km，近场区范围应不小于工程场地及外延25 km．按照该项规定，本研究选取的潜在震源区域范围及历史地震分布如图3所示．由图3可知：潜在震源区域内有三条主要断裂带，分别为淮阴-响水断裂带、郯城-庐江断裂带和六安断裂带．根据历史地震统计资料，自17世纪以来，区域内发生过三次震级在6～6.25级之间的地震，分别是1624年扬州地震、1831年淮南地震和1979年溧阳地震．此外，1970年1月～2013年12月区域内共发生了306次震级超过2级的地震，其中超过70%的地震震源深度集中在20～30 km之间．10.13245/j.hust.220811.F003图3场地潜在震源区范围及历史地震分布由《建筑抗震设计规范》(GB 50011—2010)[14]可知：工程场地所在的南京市浦口区抗震设防烈度为7度，设计基本加速度为0.1g，g为重力加速度，设计地震分组为第一组(近震)．此外，有限元数值模型中输入基岩面距离地表约80 m．综合考虑上述因素，确定了如下选波原则：a. 矩震级在6~7级之间；b. 震源深度在20~30 km之间；c. 震中距小于150 km；d. 使用地下基岩记录，且基岩埋深在50~150 m之间；e. 地震动水平方向峰值加速度在0.1~0.2g之间．根据上述筛选条件，利用地下基岩地震动选波软件(SELECTWAVE)成功得到了一组满足上述所有条件的天然地震动记录，根据地震波记录地点将其命名为静冈波．图4为输入地震动图谱，图中：T为时间；f为频率；a为地震动加速度；A为傅里叶幅值．静冈波的震级、震源深度和震源距分别为6.5级、23 km和24.705 km，地震动的水平向和竖向峰值加速度(PGA)分别为1.2410.13245/j.hust.220811.F004图4输入地震动图谱m∙s-2和0.57 m∙s-2．1.3　材料模型与参数由于目前场地非线性地震反应数值模拟中通常采用的总应力方法无法准确描述强震作用下由孔压增加引起的土体非线性行为，砂土中的细粒含量也会显著影响孔压积累所导致的土体液化特性[15]，因此，砂性土采用文献[12]提出的弱耦合有效应力分析方法，该方法结合了超静孔隙水压力增长模型[11]和扩展的non-Masing滞回本构模型[10]，建立了不规则循环荷载下循环刚度退化与超静孔隙水压力生成之间的耦合关系．黏性土及基岩采用修正Davidenkov本构模型(总应力方法)，应力-应变加卸载曲线示意图如图5所示，图中：γex和τex分别为历史最大剪应变和剪应力；γc和τc分别为当前加-卸载反转位置的剪应变和剪应力．该模型共包含13个参数，其中：Gmax，μ和β分别为剪切模量、泊松比和阻尼系数；A，2B，γr，a，C1，C2，C3，m和n为拟合参数；σv0为土体初始有效围压，根据实际场地地应力分析计算得出；γth为门槛剪应变，取0.005%．土体本构参数取值如表1所示．10.13245/j.hust.220811.F005图5修正的Davidenkov模型的应力-应变关系曲线示意图10.13245/j.hust.220811.T001表1土体本构参数材料Gmax/MPaμA2Bγr / %βaC1C2C3γth/%mn粉质黏土39.260.41.090.820.0370.003———————粉砂68.540.41.050.840.0430.0030.50.9210.1631.250.0050.3452 244.2粉细砂174.650.41.040.80.0640.0030.50.9210.1631.250.0050.3452 244.2碎石粉黏土249.530.410.70.0770.003———————花岗岩12800.410.760.10.003———————1.4　高效数值分析方法建立基于广义反应位移法模型对综合管廊结构部的连接方式进行改进，采用Timoshenko梁单元(B31)模拟预应力筋的预紧力，如图6所示．为反应场地效应，二维自由场模型根据图2真实场地条件建立，模型主要采用4节点四边形单元(CPE4R)进行网格划分，在土层交界面处用3节点三角形单元(CPE3)进行网格过渡，进而保证网格的整体质量，整个二维模型中CPE4R单元为2.192 52×105个，CPE3单元为246个．网格最大尺寸hmax取截止频率对应波长的1/8~1/10[16]．人工边界采用文献[17]提出的黏弹性人工边界等效荷载的改进算法，并选择水平方向作为地震动输入方向．管廊结构采用C50混凝土预制，弹性模量Ep为33.5 GPa，抗压强度fc为23.1 MPa，抗拉强度ft为1.89 MPa．同样采用梁单元(B31)模型两段相邻管廊间的预应力钢筋，抗拉刚度fptk为1 860 MPa，弹性模量Ep为200 GPa，钢绞线截面面积为973 mm²，预紧力值为480 MPa．10.13245/j.hust.220811.F006图6计算模型示意图将二维自由场模型进行动力时程计算，提取综合管廊位置处系列位移时程，再将其作用于管廊模型中土弹簧的自由端．由于管廊周围的土体属性差别较大，因此在地基床系数选取时上下左右各自取值．计算地基弹簧刚度采用现行国家标准《城市轨道交通结构抗震设计规范》(GB 50909—2014)[18]．2 结果与分析2.1　改进广义反应位移法的验证建立局部三维管廊精细化有限元模型验证上述方法的有效性，取图2中250 m~500 m段作为研究对象，进行三维实体动力时程分析．模型中土体和综合管廊结构均采用实体单位，三维模型尺寸为250 m(长)×约80 m(深)×37 m(宽)，模型主要采用八节点线性六面体剪缩积分单元(C3D8R)建立，共7.740 92×105个单元．管廊与土体之间采用Tie约束，相邻管段采用实体单位施加预应力进行模拟，局部三维管廊精细化有限元模型如图7所示．将静冈波加速度峰值调整为0.3g(g为重力加速度)，分别采用广义反应位移法和三维时程法对比相邻管段张开量．图8为加速度峰值为0.3g时广义反应位移法和三维时程法管廊张开量极值对比，图中：O为管节间张开量极值；L为管节所在位置．由图8可知：两种方法得到的张开量极值走势相似，位于粉质黏土、粉砂和粉细砂交界位置(约300 m处)土层交汇、刚度变化较大，出现了张开量突变，其他位置处张开量较小；同时，三维时程分析得到的张开量极值明显小于广义反应位移法，广义反应位移法计算结果更保守，这可能是由于反应位移法中结构为三向地震动输入，而三维时程分析法仅输入水平和竖向地震动，因此采用上述改进后的广义反应位移法结果可靠，偏保守的计算结果更适用于工程设计及论证．10.13245/j.hust.220811.F007图7局部三维管廊精细化有限元模型10.13245/j.hust.220811.F008图8广义反应位移法和三维时程法管廊张开量极值对比2.2　总应力法与有效应力法对比根据文献[11-12]提出的剪切-体应变耦合的孔压增量模型，基于ABAQUS软件的二次开发平台，采用Fortran编译了上述孔压增量模型本构．同样将静冈波加速度峰值调整为0.3g，图9为采用有效应力法计算的二维自由场地孔压比云图，图中：横向比例尺为1:1×104；纵向比例尺为1:5×103．由图9可知：当加速度峰值为0.3g时，仅有少量的区域发生了液化，主要集中在300 m位置附近(粉砂和粉细砂为主)．10.13245/j.hust.220811.F009图9二维自由场地孔压比云图(色标单位：m)图10为当加速度峰值为0.3g时总应力法与有10.13245/j.hust.220811.F010图10采用总应力法与有效应力法的张开量极值对比效应力法下综合管廊张开量极值对比．由图10可知：采用总应力法与有效应力法计算结果相近，张开量极值突变位置相同，均发生在土层交界处或透镜体位置．当采用有效应力法计算时，在出现液化的300 m位置附近明显张开量较采用总应力法时偏大，且当加速度峰值为0.3g时张开量已经超过了15 mm，因此须要在容易液化区域进行抗减震设计以消除安全隐患．2.3　消能节点减震措施由管廊纵向地震反应分析结果可知：在局部位置(土层软硬交界面或液化区域)相邻管节张开量偏大，必须在该位置采取相应抗震或减震措施以满足张开量防水要求．基于ABAQUS软件，自行编译形状记忆合金(SMA)简化本构模型，在综合管廊张开量最大位置处(约300 m处)布设形状记忆合金阻尼器柔性减震节点[19]，SMA简化本构模型如图11所示，图中：σ为应力；γ为应变；EA为完全奥氏体相的SMA弹性模量；α为弹性模型折减系数；εa为马氏体开始时的应变；εb为进入马氏体相变，即进入第二加载阶段加载结束时的应变；εc为按奥氏体弹性模量卸载的终止应变；εd为再次按奥氏体弹性模量卸载的开始应变．10.13245/j.hust.220811.F011图11SMA简化本构模型在上节计算得到的张开量最大位置处设置消能节点，图12为采用有效应力法，且当加速度峰值为0.3g时设置消能节点前后管廊张开量极值对比．可以看出：在设置消能节点的位置张开量极值放大，但正是由于该节点处吸收了大量的能量，明显改善了其他区域管廊相邻管节间的张开量，因此在局部位置设施耗能减震装置能有效地提高长大地下结构的抗震性能．10.13245/j.hust.220811.F012图12设置消能节点管廊张开量极值对比3 结论以南京江北新区地下综合管廊(浦乌路段)为工程背景，针对预制综合管廊管段的连接方式对广义反应位移法接头形式进行了改进．建立精细化三维实体模型，验证了改进反应位移法的有效性．基于孔压增量模型，对本构进行了二次开发，分析了强震作用下总应力法与有效应力法综合管廊纵向地震反应．针对管段张开量较大位置处设置消能减震节点，分析其减震性能，得到如下结论．a. 相较于精细化三维时程分析，改进的反应位移法计算结果偏于保守，但整体响应趋势一致，相互佐证了方法的准确性．b. 在强震作用下，二维自由场地在砂类土区域发生了液化，且综合管廊在液化区的地震响应较其他区域明显强烈．较总应力方法计算得到的结果，采用有效应力法除在液化区有明显的增大外，其他区域并不明显．c. 在设置消能节点位置有明显的管节张开量放大，其他位置处的张开量明显降低，表明消能节点的设置可有效减小接头张开量值，提高管廊地震安全储备．