基于ABAQUS的三维粘弹性边界单元及地震动输入方法研究

苑举卫 杜成斌 陈灯红

（河海大学工程力学系,南京 210098）

在结构-地基动力相互作用分析中为模拟近场能量向无限域的散射,通常需要在近场有限离散模型的人工截断处设置人工边界,以正确反映结构-地基系统的动力特性.其中局部人工边界由于其时空解耦特性而受到广泛重视,目前较常用的有粘性边界[1]、粘弹性边界[2-4],以及廖振鹏[5]提出的透射边界等.粘性边界只考虑了边界的耗能作用,而未考虑介质的弹性恢复作用,在应用中容易发生结构整体飘移,精度较低；透射边界虽具有较高精度,但在实际应用中一般仅限于二阶精度以内,并且编程较复杂、计算中可能引起高频失稳；粘弹性边界能同时模拟波的散射效应和半无限地基的弹性恢复能力,且能克服粘性边界引起的低频漂移问题,稳定性好.

在考虑地基辐射阻尼的情况下,将地震动引入计算模型存在振动输入和波动输入两种方式,前者适用于结构刚度相对较小、输入波频率低、结构尺寸远较其需考虑的最短波长为小的情况,目前一般工业和民用建筑物和构筑物多采用这种方式；而对于重要的大坝工程,由于坝体尺寸和重量都很大,地基对结构体系动态特性的影响较大,这种情况一般都作为开放系统的波动问题求解,即将地震动输入作为外源问题引入计算区域,波场分解法和等效边界力法基于不同的思想均有效的建立了适用于粘弹性边界的外源输入方法[6-7].

ABAQUS作为功能最强大的有限元软件之一,可以分析各种固体力学、结构力学系统,尤其是能够处理非常复杂的问题和模拟高度非线性问题.ABAQUS还为用户提供了方便灵活的二次开发平台,包括若干用户子程序（user subroutines）以及实用程序（utility routines）,为解决实际工程的复杂问题提供了方便.

为方便地实现粘弹性人工边界的设置及地震动的有效波动输入,本文利用有限元软件ABAQUS提供的自定义单元子程序UEL结合隐式积分算法[8],通过定义边界单元对整个系统的雅克比矩阵贡献来实现粘弹性人工边界.基于波场分离方法,将人工边界上的总波场分解为无局部场地效应影响的自由场与局部场地效应引起的散射场,散射场由粘弹性人工边界吸收,则作用于人工边界的无局部场地效应影响的自由场可由设计地震动位移时程和速度时程统一表达,利用子程序DSLOAD、UTRACLOAD描述边界的加载历程.将子程序 UEL、DSLOAD、UT RACLOAD以及设计地震动等数据汇集成 SDAB.FOR,结合在ABAQUS中建立的有限元模型,可方便有效地进行结构-地基相互作用的分析.

1 粘弹性人工边界的有限元实现

粘弹性人工边界描述的是由近场向远域散射的外行波在人工截断边界处满足的应力条件,可表示如下

式中,σ表示边界点的应力向量；u、˙u分别为边界点的位移和速度向量；KB、CB分别为以粘弹性人工边界弹簧系数、粘性系数为元素的对角矩阵.粘弹性人工边界的实质是在人工边界上施加切向和法向方向上的弹簧和阻尼器,其系数按式（2）计算

式中,KBN、KBT分别表示法向和切向弹簧系数；CBN、CBT分别为表示法向和切向阻尼系数；rb为散射波源至人工边界的距离；cp、cs分别为P波和S波波速；G为介质剪切模量；λ为拉梅常数；ρ为介质质量密度；α、β为无量纲参数,分别取0.8、1.1.

粘弹性人工边界作为一种应力连续分布的人工边界条件,可采用与普通单元类似的形函数进行离散,即一致粘弹性人工边界,经离散后,粘弹性人工边界单元的刚度矩阵[9]如式（3）,阻尼矩阵具有相同的形式,只需将刚度系数K用相应的阻尼系数C替换即可.

在计算区域边界外侧拓展一层六面体单元,单元外侧节点施加固定约束,一致粘弹性人工边界单元的刚度矩阵、阻尼矩阵通过ABAQUS的自定义单元子程序UEL,结合动力分析所采用的 Hilber-Hughes-T aylor隐式积分方法,以边界单元对整个系统的矩阵贡献的形式写入.单元的厚度只有象征意义,没有限制,对结果没有影响.

其中,Sb为三维一致粘弹性人工边界的面积[10].

2 地震动输入方式

2.1 粘弹性边界外源波动输入

基于波场分离,将人工边界上的总波场分解为无局部场地效应影响的自由场与局部场地效应引起的散射场,局部场地效应引起的散射场由粘弹性人工边界条件可自动满足,将外源地震波引入计算区域,只需在人工边界处施加自由场荷载,把近域地基作为半无限空间自由场地基的子结构,由其在入射地震波作用下的边界相互作用力[11]给出:

式中,n为边界外法线方向余弦向量；KB、CB均与式（1）含义相同分别为自由场位移矢量 、速度矢量和应力矢量,当在近域地基底边边界输入的地震波波阵面为平面时,自由场应力矢量可由速度场表示.

假定地震波为垂直入射的平面P波和平面S波,则根据波动传播规律及波场应力状态边界相互作用力（4）可由设计地震动位移时程和速度时程统一表达,以分布荷载的形式施加于人工边界,其中边界法向荷载采用子程序DSLOAD描述,切向荷载由子程序UTRACLOAD描述.

2.2 粘弹性边界内源地震动输入

在计算区域人工截断处设置粘弹性人工边界来反映无限区域的情况下,地震荷载的引入还有另一种输入方式,即考虑粘弹性边界的内源地震动输入方法.在均匀输入情况下,或者空间不均匀地震动较小的情况下,地震动产生的荷载只与结构质量有关,而与地基质量无关,但是系统的响应包含地基质量,这也是与常规无质量地基不同的地方[12].

因此,本文将考虑3种地震动输入方式,即考虑粘弹性边界条件的外源波动输入（输入方式1）和内源输入（输入方式2）,及常规的无质量地基输入（输入方式3）,比较不同输入方式对结构-地基系统动力响应的影响.

3 外源波动算例验证

考察均匀三维弹性半空间,从底部人工边界分别竖直入射平面S波、P波的动力响应.入射位移脉冲波方程为

从均匀三维弹性半空间截取长和宽均为400m、高度为600m的有限区域作为计算区域,用有限单元法进行空间离散,3个方向网格尺寸均为20 m,共离散为13671个结点,12000个六面体单元.在计算区域外侧及底侧拓展一层六面体单元,构造一致粘弹性边界单元.自由场荷载作用于人工截断处各表面.材料参数取为:弹性模量E=4.88GPa,密度ρ=2000 kg/m3,泊松比ν=0.22.P波波速cp=1669.05m/s,S波波速cs=1000m/s,时间步长取为Δt=0.01 s.

该问题的自由地表位移解析解为考虑行波延迟后放大2倍的入射位移时程.图1、图2分别为各观测点在平面S波入射时的水平向位移响应和P波入射时的竖向位移响应.可见入射波到达地表后与反射波叠加,地表的响应为入射波的2倍；底面观测点在经历入射波和经地表反射的波后响应为 0,波传向远域,均与解析解一致.数值模拟结果符合波的传播规律,且有良好的模拟精度,说明本文的三维粘弹性边界单元及外源波动输入子程序SDAB.FOR是准确的.

图1 S波入射时底面点和地表点的位移时程

图2 P波入射时底面点和地表点的位移时程

4 工程应用

某水电站位于我国澜沧江上游河段,坝体为抛物线变厚度双曲拱坝,坝底高程953 m,坝顶高程1245 m,最大坝高292m,坝顶弧长 935m,坝顶厚 12m,坝底厚73m,坝体混凝土动态弹性模量Ed=20GPa,密度 ρd=2400kg/m3,泊松比 νd=0.189,地基弹性模量Ed=27.3 GPa,密度 ρd=2400 kg/m3,泊松比 νd=0.25.

根据场地地震安全性评价成果,工程区设计烈度为Ⅸ度,水平向设计地震动峰值加速度为0.308g,竖向设计地震动峰值加速度取水平向的2/3,即0.205 g,根据规范谱人工合成地震波如图3所示,计算总时间为10 s,时间步长为0.005s.

采用文献[13]对不同地基范围粘弹性人工边界的精度研究结论,地基部分沿结构各方向延伸一倍坝高就可满足工程精度要求,边界设置粘弹性人工边界.坝体和基岩用三维8结点六面体单元进行有限元离散,共划分31026个单元,35366个结点,有限元网格如图4所示.图5给出了个关键部位的位置示意图,括号中表示对应下游面.

图4 拱坝-地基系统有限元模型 图5 关键部位示意图

在上述3种不同地震动输入方式下,上、下游面关键部位的应力峰值如表1、2所示,其中括号中百分数表示与常规地震动输入方法计算结果的降低幅度,拱冠处上下游的主拉应力时程曲线如图6所示.从表中可以看出,除个别点外,采用无质量地基模型的地震动常规输入方法的响应最大,引入粘弹性人工边界后,坝体的动力响应显著降低,且采用外源波动输入方法计算的响应最小.

表1 上游面各关键部位应力比较 （单位:MPa）

续表1 上游面各关键部位应力比较（单位:MPa）

表2 下游面各关键部位应力比较 （单位:MPa）

外源波动输入方法与常规地震动输入方法相比,无论是主拉应力还是主压应力降幅都较大,最大降幅达53%.相互作用无质量地基与常规地震动输入方法相比,主拉应力降幅较大,其中上游面主拉应力降幅在10%～35%之间,下游面主拉应力降幅在15%～45%之间,上下游面主压应力降幅均25%.这两种输入方法相比在拱冠处的响应较其他部位接近.

5 结 论

本文利用FORT RAN语言编写了适用于粘弹性人工边界的外源波动输入子程序SDAB.FOR,包含粘弹性边界和加载两部分,从而可方便地实现地震动的有效输入.将考虑粘弹性边界条件的外源波动输入和内源地震动输入方法分别应用于某高拱坝的动力响应分析,并与常规无质量地基模型计算结果进行了对比,结果表明,考虑地基辐射阻尼的两种地震动输入方法均会显著的降低结构的动力响应,其中采用外源波动输入方法的降幅更大.

[1]Lysmer J,Kulemeyer R L.Finite Dynamic Model for Infinite Media[J].Journal of Engineering Mechanics Division,ASCE,1969,95（4）:759-877.

[2]Deeks A J,Randolph M F.Axisymmetric Time-domain Transmitting Boundaries[J].Journal of Engineering Mechanics Division,ASCE,1994,120（1）:25-42.

[3]刘晶波,谷 音,杜义欣.一致粘弹性人工边界及粘弹性边界单元[J].岩土工程学报,2006,28（9）:1070-1075.

[4]杜修力,赵 密,王进廷.近场波动模拟的一种应力人工边界[J].力学学报,2006,38（1）:49-56.

[5]廖振鹏.工程波动理论导论[M].北京:科学出版社,2002.

[6]刘晶波,吕彦东.结构-地基动力相互作用问题分析的一种直接方法[J].土木工程学报,1998,31（3）:55-64.

[7]赵建锋,杜修力,韩 强等.外源波动问题数值模拟的一种实现方式[J].工程力学,2007,24（4）:52-58.

[8]Abaqus Theory Manual（Version 6.7）,Dassault Systems Inc,2007.

[9]谷 音,刘晶波,杜义欣.三维一致粘弹性人工边界及等效粘弹性边界单元[J].工程力学,2007,24（12）:31-37.

[10]徐 磊,叶志才,任青文等.基于ABAQUS的粘弹性动力人工边界精确自动施加[J].三峡大学学报:自然科学版,2010,32（1）:20-23.

[11]陈厚群.坝址地震动输入机制探讨[J].水利学报,2006,37（12）:1417-1423.

[12]李 静,陈健云.动力相互作用分析中无质量地基的应用研究[J].世界地震工程,2007,23（2）:58-62.

[13]Zhang Chuhan,Pan Jianwen,Wang Jinting.Influence of seismic input mechanisms and radiation damping on arch dam response[J].Soil Dynamics and Earthquake Engineering,2009,29（9）:1282-1293.