对称地下结构抗震分析的边界强制反应位移法

2021-05-25韩润波许成顺许紫刚蒋家卫杜修力

工程力学 2021年5期

韩润波，许成顺，许紫刚，蒋家卫，杜修力

(北京工业大学城市与工程安全减灾教育部重点实验室，北京 100124)

长久以来，人们普遍认为由于周边土体约束效应和地震动沿深度的衰减特征，地下结构具有良好的抗震性能，因此对其抗震安全性缺乏足够的重视，研究工作很少。1995 年日本阪神地震中出现了大规模地下结构震害现象，特别是Daikai车站发生了塌毁震害[1]，地下工程抗震问题才引起工程界和学术界的广泛关注[2]。

物理模型试验是研究地下结构抗震问题中有效且直观的手段，地下结构抗震试验常用的方法有普通振动台试验[3-4]、离心机振动台试验[5-7]及拟静力试验[8-9]。在构件及结构的拟静力试验中，地震荷载是通过对构件或结构施加强制水平位移实现的，通常无法反映土-结构相互作用。土-地下结构体系拟静力试验中，地震荷载是通过对土-结构体系侧边界施加特定分布的强制水平位移实现的，此类试验几何缩尺比较大且土体应力水平可控，可近似还原真实土体应力状态及可反映土-结构相互作用，可用于地下结构的抗震性能研究，特别是对结构细部构造要求较高的地下结构抗震试验较为适用。综上，大型土-地下结构体系拟静力试验是地下结构抗震物理模型试验一种发展趋势。JSCE committee (1992)[9]开展了首例土-地下结构体系的拟静力试验，试验将单层双跨的矩形钢筋混凝土地下结构缩尺模型置于盛有干砂的土箱中，推动两侧刚性墙对土-地下结构体系进行静力往复剪切试验以观察地下结构地震作用响应。图1为该试验方案示意图，此试验中采用的刚性墙推覆方式无法准确体现土体在地震作用下的剪切变形模式。徐琨鹏[10]开展了自由场土体拟静力推覆试验，试验中采用了实际中操作性较强的刚性板的推覆方式，此种推覆方式仅能实现倒三角形侧边界位移分布形式。

图 1 JSCE committee 开展的土-结构体系拟静力试验方案示意图[9]Fig. 1 Illustration of quasi-static test of soil-structure system carried out by JSCE committee[9]

土-地下结构体系拟静力推覆试验可较准确地得到地下结构地震破坏模式并揭示地下结构震害机理，为克服上述拟静力试验的局限性，作者拟开展可较好模拟实际地震中土体剪切变形模式且侧边界位移分布形式可调的土-地下结构体系拟静力推覆试验。拟开展试验中，模型几何缩尺比为1/5，采用层状剪切箱作为土-结构体系的容器，通过竖向加载系统在土-结构体系表面施加恒定的竖向荷载以提高土体应力水平，通过水平加载系统对不同高度处的层状剪切箱侧面施加一定分布的强制水平位移来模拟水平地震作用，并逐步增大边界强制水平位移至结构破坏，以期获得结构构件弹性-开裂-屈服-弹塑性-承载力下降的全过程能力曲线并评价地下结构的抗震性能。上述拟静力试验装置已发表相关专利[11]，试验方案示意图如图2 所示。

图 2 土-结构体系拟静力推覆试验方案示意图Fig. 2 Illustration of quasi-static pushover test of soil-structure system

然而受限于当前试验技术手段，土-地下结构体系拟静力推覆试验中的水平地震作用只能考虑在土-结构体系侧边界上施加强制水平位移的方式实现，无法体现实际中的地震惯性力作用[12]，需要通过数值模拟方法定量研究两种不同荷载形式下地下结构受力反应的差异。因此，本文提出一种采用侧边界强制施加水平位移来模拟水平地震作用的简化分析方法，即边界强制反应位移法，并详细介绍了本文方法的实施步骤、基本特点。以1995 年日本阪神地震中遭受严重破坏的Daikai车站为算例，将本文方法、传统的强制反应位移法及动力时程方法计算结果进行了对比分析，验证了本简化分析方法的有效性。本文方法同样采用侧边界强制水平位移的方式施加水平地震作用，因此，可为土-地下结构体系拟静力推覆试验提供一定的参考。

1 相关地下结构抗震简化分析方法

地下结构抗震常用简化分析方法中水平地震作用大多采用力或加速度方式加载[13-19]，采用位移方式加载的简化分析方法相对较少，仅有强制反应位移法[20]、整体强制反应位移法[21]等。强制反应位移法首先计算地震作用下自由场位移，随后将上一步计算的自由场位移施加在土-结构拟静力计算模型侧边界处，认为此时拟静力模型中结构反应即为其地震响应，该方法示意如图3(a)所示。但该方法计算结果误差较大，是由于其施加的自由场变形只能在拟静力模型侧边界及附近土体范围内被准确传递，到土体中部位置时，变形衰减较大，Tateishi[20]研究表明该方法在边界处施加的位移传递到模型中部时，应变已衰减20%～50%，并提出修正的强制反应位移法，该方法示意图如图3(b)所示，将自由场一维地震反应计算得到的地层应变换算成等效节点荷载，随后将该等效荷载施加到土-结构整体模型上，提高了强制反应位移法的精度，但该方法中水平地震作用采用力的方式加载，本质上与反应加速度法一致。陈之毅等[21]提出了整体强制反应位移法，该方法示意图如图3(c)所示，首先进行自由场一维土层地震反应分析，读取自由场地下结构顶、底板位置处土体达到最大水平相对位移时刻时土层水平位移沿深度的分布，将该位移分布按深度强制施加到土-结构拟静力模型上的全部土体上，计算结构在地震作用下的响应。该方法拟静力模型中整体土层水平位移分布均与自由场地震反应分析得到的水平位移分布一致，提高了强制反应位移法的精度及可靠性，但结构周围土体也同样被强制施加自由场土体变形，因此对土-结构相互作用考虑不充分。随后陈之毅等[22]提出了修正的整体强制反应位移法，通过结构剪切变形修正系数近似考虑土-结构相互作用。整体强制反应位移法中荷载施加方式是对同高度土体施加同一水平位移，该方法并不适用于大型土-地下结构体系拟静力推覆试验。

图 3 相关抗震设计方法模型示意图[20-21]Fig. 3 Illustration of models of relative seismic design methods[20-21]

因此，结合大型土-地下结构体系拟静力推覆试验的需求，本文借鉴传统的强制反应位移法中水平地震作用的施加方式，发展了一种边界强制反应位移法，可为土-地下结构体系拟静力推覆试验提供一定的参考。

2 边界强制反应位移法

2.1 计算模型

由柔度系数法可知，在分析已知场地条件下的地下结构地震反应时，土与结构刚度比已经确定，则自由场动力计算时地下结构顶、底板外表面对应的土体最大水平位移差(以下简称为RD 值)很大程度地影响了地下结构地震反应。当采用在土-地下结构体系侧边界处施加强制水平位移以模拟水平地震作用时，由于侧边界强制位移向土体中部传递过程中存在衰减性，地层变形传递到结构位置处附近时，结构所受的地层变形作用已与边界处施加的地层变形作用不一致，这也是传统的强制反应位移法误差较大的原因。本文中边界强制位移指有限元模型中侧边界上施加的强制水平位移。为修正边界强制位移传递中的衰减，本文方法建立两个拟静力有限元模型，模型1 是自由场土体模型(以下称为附加全自由场模型)，如图4(a)所示，模型2 是土-结构相互作用模型，如图4(b)所示，模型1 及模型2 尺寸一致。对于形状对称的地下结构而言，拟静力模型中结构中线处土体变形反应较能代表结构所受到的地震荷载作用。

图 4 边界强制反应位移法拟静力计算模型Fig. 4 Calculation model of boundary forced response displacement method

在模型1 侧边界处施加按比例单调递增的自由场地震反应水平位移分布，当模型1 中部(结构中线位置)对应结构顶底板外表面位置处土体水平位移差与自由场土体动力计算得到的RD 值一致时，记录此时施加在模型1 上的边界强制位移U1，将U1 施加在同尺寸的土-结构模型2 侧边界处，此时模型2 中结构反应即为结构在特定地震动作用下的内力及变形响应。

2.2 实施步骤

根据前文所述可知，采用本文方法进行地下结构抗震计算时，应首先求得目标位移差RD 值及侧边界强制水平位移分布形式，随后建立两同尺寸有限元模型进行拟静力分析。因此，本文方法实施步骤如下：

1)求解目标位移差RD 值及侧边界强制水平位移分布。对场地进行自由场地震反应分析，得到自由场在地震作用下的目标位移差RD 值，并取对应地下结构顶、底板外表面位置处土体达到最大水平相对位移时刻的土层水平位移沿深度的分布作为侧边界强制水平位移分布形式。

2)建立拟静力模型并进行地应力平衡。建立附加全自由场及土-结构拟静力模型，模型底边界固定，约束侧边界水平向位移，对模型施加重力荷载进行计算，并进行地应力平衡。

3)重新定义拟静力模型边界条件。为在模型侧边界处施加强制水平位移，取消模型侧边界处水平向约束代之以地应力平衡后的水平支座反力，并约束侧边界各节点竖向位移，模型底边界保持固定。

4)施加侧边界强制位移进行分析。分别在附加全自由场及土-结构拟静力模型侧边界施加第1)步确定的侧边界强制水平位移分布，并按相同比例单调递增加载，直至附加全自由场模型中部(结构中线位置)对应结构顶底板外表面位置处土体水平位移差达到第1)步中求得的目标位移差RD值，此时土-结构模型中结构的变形及内力反应即为特定地震动作用下的结构响应。

3 方法验证

3.1 模型建立

以1995 年日本阪神地震中破坏的Daikai 车站为例进行计算分析。Daikai 车站为两跨中柱地下结构，埋深4.8 m，车站标准横断面如图5 所示。既往研究表明：对于Daikai 车站，中柱及边墙底部是较危险的截面[23]，因而选取此两截面作为控制截面，位置见图5 所示。车站结构周围土层的物理性质如表1 所示[24]，土体本构模型采用基于Davidenkov 骨架曲线的非线性粘弹性模型，各参数详见表2[25]。

图 5 Daikai 车站标准横断面图 /mFig. 5 Standard cross-section of Daikai subway station

采用通用有限元软件ABAQUS 对Daikai 车站标准横断面进行二维地震反应分析，建模时结构及土体均采用平面应变实体单元。因中柱在车站纵向是3.5 m 等间距分布的，需将中柱按一定原则等效成一个柱间距长度的一面纵墙，并和顶板、底板及侧墙一样，取单位长度作为研究对象。顶板、底板和侧墙弹性模量取为30 GPa，密度取2500 kg/m3。为保证等效前后截面的抗弯刚度、抗剪刚度、抗压刚度及截面质量均不改变[26]，等效后的中柱弹性模量取为8.57 GPa，密度取为714 kg/m3。

表 1 土层物理性质表[24]Table 1 Physical properties of soil[24]

表 2 Davidenkov 模型参数A、B 以及γ0 的参考值[25]Table 2 Reference value of A, B and γ0 of Davidenkov model[25]

建立同尺寸的附加全自由场模型及土-结构模型，模型高度为地表至基岩面间距，即39.2 m。模型宽度为结构外宽的3 倍，拟静力模型如图6所示。土-结构模型中土体与结构表面之间的接触相互作用：法向采用硬接触，单元之间互不侵入；切向采用摩擦接触，摩擦系数取0.4。

同时，为分析本文方法的计算精度及可靠性，对上述算例另采用传统的强制反应位移法及动力时程方法进行计算，传统的强制反应位移法及动力时程方法有限元模型中结构尺寸、材料参数及土层参数等设定与上述本文方法计算模型一致。传统强制反应位移法模型宽度为结构外宽的3 倍，动力时程方法模型宽度为结构外宽的7 倍，模型高度均为地表至基岩面间距[27-28]。动力时程计算采用振动法输入地震动，计算模型如图7 所示，振动法输入较适合分析存在下卧刚性基岩场地条件的土-结构动力相互作用问题[29-30]，模型侧边界采用捆绑边界，该人工边界将同高度处的侧边界土体节点捆绑在一起作一致的运动，在土-结构动力时程计算中具备较好的计算精度及计算效率[29]。

图 6 拟静力有限元模型Fig. 6 Quasi-static finite element model

图 7 动力时程计算模型[29]Fig. 7 Model of dynamic time history method[29]

将动力时程计算结果作为计算校核的基准，为验证本文方法对于不同特性及不同峰值强度地震动的适用性，计算采用El-Centro 地震动、Kobe地震动、Loma Prieta 地震动分别为0.1g、0.2g、0.3g及0.4g峰值加速度，共计12 种工况，各地震动时程曲线如图8 所示。

3.2 计算过程及结果分析

图 8 输入地震动时程曲线图Fig. 8 Input acceleration time histories

首先计算自由场位移分布，采用动力有限元法进行土层自由场动力反应分析，自由场动力模型与土-结构动力模型采用相同的土体本构、材料参数及边界条件。计算El-Centro 地震动、Kobe 地震动、Loma Prieta 地震动峰值加速度分别为0.1g、0.2g、0.3g及0.4g共12 种工况，每种工况进行自由场动力分析后得到目标位移差RD 值及对应时刻土层水平位移沿深度的分布，各水平位移分布如图9 所示。观察图9 可知，对于同一特性地震动，目标位移差RD 值随着输入地震动峰值加速度的增大而增大，自由场水平位移沿深度分布曲线形式基本一致。

图 9 各地震动不同峰值加速度下自由场水平位移沿深度分布图Fig. 9 Displacement distributions of free-field model under earthquake motions with different peak accelerations

表3 给出了传统的强制反应位移法(M1)及本文方法(M2)计算得到的结构中柱底部截面内力、边墙底部截面内力以及结构层间相对位移，表中计算误差为上述两方法计算值与动力时程方法计算基准值的相对误差，结构内力计算值均为每延米内力计算结果。由表3 可以看出，本算例在不同特性不同峰值强度地震动工况下，传统的强制反应位移法计算误差大部分在20%～50%，而本文方法计算误差基本在15%以内。可知，当模型宽度取值合理时，本文方法计算精度显著高于传统的强制反应位移法计算精度，这是由于传统的强制反应位移法模型侧边界施加的地层变形在向模型中部传递过程中会衰减，最终施加在结构上的土层变形作用已远小于侧边界处土层变形作用，侧边界位移传递衰减规律见3.3 节。本文方法修正了地层变形自侧边界向中部传递中的衰减，从而使结构所受的地层变形作用接近于自由场动力分析得到的最大地层变形作用。

表 3 各分析方法计算结果对比Table 3 Comparison of results calculated by different methods

3.3 侧边界位移传递衰减规律

当在拟静力模型侧边界处施加强制水平位移时，该强制位移向中部传递过程中会衰减[20]。建立3.1 节中所述附加全自由场模型，自侧边界至中线等间距取四个监测面，分别为A1、A2、A3、A4，相邻监测面间距为B/6(其中B为拟静力模型宽度)，如图10 所示。位移衰减程度 γ计算方法见式1，其中 ΔU为相应监测位置处结构顶、底板外表面高度位置的土体水平位移之差。

将图9 中所示各地震动位移分布，输入至图10所示拟静力模型侧边界处，观察A1～A4 各监测面位移衰减程度 γ，如图11 所示。

观察图11 可知，对于不同峰值强度的各地震动，侧边界位移向中部传递过程中衰减程度略有不同，本算例中拟静力模型宽度为3 倍结构外宽时，模型中线位置处 ΔU衰减程度在45%～55%。对于不同峰值强度的同一地震动而言，位移衰减程度随地震动峰值强度的增加而略微增大；越靠近于模型中线位置，位移衰减程度越大，但衰减程度的增幅越来越小，即越靠近中线位移衰减逐渐缓慢。

图 10 拟静力模型监测位置图Fig. 10 Monitoring diagram of quasi-static model

图 11 各边界强制位移分布下位移衰减图Fig. 11 Displacement attenuation under different boundary displacement distributions

分别建立宽度为1.8b、2b、3b、5b及7b(b为结构外宽)的附加全自由场模型，模型参数同3.1节所述，取图9 中峰值强度为0.4g的El-centro 地震动计算得到的位移分布作为侧边界处输入的位移分布，得到不同宽度模型中部位移衰减程度如图12 所示。由图12 可知，随着拟静力模型宽度的增加，模型中部位移衰减程度逐渐增大。采用本文2.2 节所述步骤修正侧边界位移传递过程中的衰减，使各附加全自由场模型中部对应结构顶底板外表面位置处土体水平位移差等于目标位移差RD 值，修正后各模型中部土体水平位移沿深度的分布如图13 所示，图13 中黑实线(自由场)代表峰值强度为0.4g的El-centro 地震动自由场动力分析得到的水平位移分布。观察图13 可知，随着模型宽度的增加，土体中部水平位移分布与自由场动力分析得到的水平位移分布差异逐渐变大。

图 12 不同宽度模型中部位移衰减程度Fig. 12 Displacement attenuation in middle of model with different widths

图 13 修正后各宽度模型中部土体水平位移沿深度分布Fig. 13 Modified displacement distributions in middle of model with different widths

4 拟静力模型宽度影响分析

由上文可知，拟静力模型宽度影响侧边界强制水平位移的衰减程度。为判断本文方法在不同拟静力模型宽度下的适用性，对上述算例分别建立模型宽度为1.8b、2b、3b、5b、7b(b为结构外宽)的拟静力模型进行分析，计算工况为El-Centro地震动峰值加速度分别为0.1g、0.2g、0.3g、0.4g，对比了不同模型宽度时采用本文方法及动力时程方法计算得到的中柱底部弯矩及剪力值、边墙底部弯矩及剪力值、结构层间相对位移，如图14 所示。观察图14 可知，除边墙底部剪力计算结果外，本文方法计算结果误差基本随模型宽度增大而变大，当模型宽度为5b时，最大误差约21%；当模型宽度为7b时，最大误差约29%。这是由于随着模型宽度的增加，采用本文方法修正后的土体中部水平位移分布与自由场动力分析得到的水平位移分布差异逐渐变大，见图13 所示。而对于边墙底部剪力而言，本文方法计算误差则随模型宽度减小而变大，当模型宽度为2b时，最大误差约14.6%；当模型宽度为1.8b时，最大误差约16%。这是由于边界强制位移向模型中部传递中存在衰减，因此作用在边墙处的相应土体位移差 ΔU1大于作用在模型中部处的相应土体位移差 ΔU，而越靠近模型中线，位移衰减越缓慢，因此，模型宽度越大，结构边墙处相应土体位移差 ΔU1越接近于结构中线处相应土体位移差 ΔU，边墙剪力计算误差也就越小。综上，为保证一定的计算精度及稳定性，本文方法中拟静力模型宽度取值范围一般可为结构外宽的2 倍～3 倍。