基于UPFs的珊瑚岛礁场地地震反应分析方法研究

朱圣华, 尹训强, 王桂萱(大连大学 建筑工程学院, 辽宁 大连 116622)

0 引言

珊瑚岛礁相关理论技术的研究对推进我国南海开发建设具有重大意义,与此同时,由于面临着复杂海域地震地质环境,珊瑚岛礁场地地震稳定性评价是近年来地震工程研究的热点问题[1-3]。其中,如何合理描述能够反映珊瑚砂动力特性的地震反应计算模型是进行珊瑚岛礁抗震安全分析的关键技术问题[4-5]。

针对陆域工程场地非线性地震反应分析国内外学者开展了全面系统的研究,取得了许多有益的成果[6-8]。相对而言,对于珊瑚岛礁场地非线性反应分析方法的报道较少,但近几年也逐渐成为国内外学者的研究热点。胡进军等[9-10]对典型岛礁场地地震效应开展了研究,在考虑海域特殊场地环境及礁岩特殊构造的基础上,通过建立一维的土层分析模型,讨论了脉冲型地震动作用下地震反应特征;李天男[11]结合地质勘测资料,建立一维场地土层剖面模型,利用DEEPSOIL研究不同地震类型输入条件下礁坪区等四类岛礁场地的地震动特性参数。徐长琦[12]选用OpenSees中开发的塑性本构模型模拟其通过物理实验测得的饱和钙质砂特性,并实现了海水和孔隙水压力的模拟。通过对岛礁资料的整理,建立简化的二维分析模型,求索了礁灰岩倾角与上覆珊瑚砂厚度的不同对场地反应的影响规律。陈国兴等[4]从网格划分、人工边界条件以及材料本身特性出发建立二维模型,主要分析了地表加速度谱形、持时以及地表峰值加速度放大效应的发展规律。张巍等[5]在文献[4]的基础上,利用FLAC建立了海水—岛礁场地的耦合计算模型,对比分析了海水与地震动综合作用下的地震反应规律。然而,鉴于珊瑚岛礁海洋工程地质体及钙质砂动力力学特性的特殊性,如何合理且高效地模拟珊瑚砂非线性动力特性及远场无限地基辐射阻尼的影响,以及提高计算模型的可操作性,仍值得进一步研究。

为此,本文以通用有限元软件ANSYS为研究工具,基于UPFs二次开发灵活性的特点,建立了适用于珊瑚砂地基条件下的岛礁场地地震反应分析时域计算模型,并通过自由场动态分析验证其良好合理有效性。进而,结合国内某实际珊瑚岛礁场地,探究珊瑚岛礁场地非线性地震反应规律。

1 珊瑚岛礁场地地震反应分析模型

如图1所示,珊瑚岛礁场地地震反应分析模型由珊瑚岛礁、近场地基以及半无限地基所组成,其中,岛礁四周还要受海水的作用。如前所述,构成岛礁场地的珊瑚砂的非线性动力特性、远场无限地基辐射阻尼以及海域岛礁动水压力的影响是需要考虑的关键技术问题。

图1 珊瑚岛礁场地地震反应分析模型示意图Fig.1 Schematic diagram of seismic response analysis model for coral island site

1.1 珊瑚砂非线性动力特性

珊瑚砂作为海洋生物死后残骸经过长期的生物及理化作用产生的特殊岩土类介质,具有易破碎、棱角高、多孔隙等特性,与陆域工程场地有较大差异。陈国兴等[4]通过对某南海岛礁珊瑚砂试样的动三轴试验,结合理论基础给出了修正之后Matasovic本构模型的参考系数,同时验证了珊瑚砂与其较好的相适性,如图2所示。另外,该本构模型可通过等价线性法[13]进行分析,该方法概念明确且优势突出,仍是场地地震反应分析不可替代的主流方法。

图2 珊瑚砂典型的G/Gmax-γ和D-γ曲线Fig.2 G/Gmax-γ and D-γ curves of coral sand

等价线性法是将珊瑚砂的非线性动力问题近似地描述为多次迭代的线性分析方法,珊瑚砂力学性质将在迭代步之间进行校正,以迭代的方式保证与应力应变相关联的珊瑚砂的力学特征。也就是,通过等效剪应变γeff取0.65倍的最大剪应变γmax参照图2所示的G/Gmax-γ和D-γ曲线重新得出修正剪切模量和阻尼比,直到收敛后的应变与珊瑚砂性质相一致。

另外,具体计算中,为反映随岛礁场地深度的变化,在计算时可以每个单元赋值不同的初始计算参数。依据工程经验,绝大数情况下5～6步以内即可完成收敛。

1.2 黏性人工边界

等价线性法的线性化处理将珊瑚砂场地在每次动力计算中都可视为线性材料,因此,可采用工程场地中常用的黏性人工边界来考虑无限地基辐射阻尼的影响。

黏性人工边界[14-16]场地模型的目的是模拟外源地震动的输入以及消耗散射波的波动能量。如图1所示,在二维计算模型中,需在截断外边界节点处施加的等效荷载如式(1)。

(1)

对层状分布较显著的珊瑚岛礁场地,需在进行动力时程分析之前,通过自由场分析获取截断外边界节点处的响应,进而,利用式(1)计算相应节点处的等效荷载力,从而作为珊瑚岛礁场地的地震动输入条件。

1.3 Westergaard动水压力附加质量

考虑到海域岛礁位置的特殊性,需要注意的是周边海水对岛礁的影响。其中最为重要的是海水与岛礁接触的耦合作用,海水荷载的作用会造成土体的形变或运动,而土体的形变又会反过来影响海水荷载的分布和大小。目前对于此类耦合作用问题,大多工程考虑其仅发生在海水与岛礁的接触面,通过引入接触面的平衡协调关系来解决问题。

Westergaard是公认的研究土木水利工程领域研究水坝耦合作用的先驱,其对附加动水压力的简化因为简单实用在各国的抗震规范中得到了广泛的运用。但是其限制性假设条件过多,对于结果也存在一定的影响,Clough教授对其进行了改良使之适应于任意形状河床场地,能够不需拘泥于水平向的地震加速度。本文所用具体表达式如下:

(2)

式中:a为地震动加速度;h是总的水的高度;y为该节点水的深度;ρ是海水的密度大小。

2 基于UPFs的岛礁场地地震反应分析模型开发

2.1 珊瑚砂等价线性单元

等价线性法的计算特点是动力时程分析过程的多次迭代,与传统有限元软件的计算荷载步内的迭代有较大的差异性,为解决此问题,基于通用有限元软件ANSYS,并运用其UPFs的二次开发特性,创建了一种新的珊瑚砂等价线性单元描述岛礁近场场地的非线性动力特性,该单元形式可实现修正的Matasovic本构模型的嵌入,是已有等参单元材料模型的修改。

所创建单元可通过修改FORTRAN接口子程序UserElem.F编译完成,具体修改编译过程如图3所示。其中,通过子程序不仅可以进行弹性模量E、泊松比v以及密度ρ等单元参数的获取,还能同时创建多个单元并由KEYOPT关键选项进行转换。值得注意的是,珊瑚砂等价线性单元与标准单元的不同关键在于不仅通过读取外部文件即可获取计算参数,且能够将等价线性迭代过程中阻尼比D和剪切模量G的变化进行传递,并通过自主设置的误差阈值(一般取0.05)作为迭代是否终止的收敛依据。

图3 珊瑚砂等价线性单元接口子程序的编译过程Fig.3 The compiling process of equivalent linear element interface subroutine of coral sand

2.2 黏性人工边界单元

实现黏性人工边界场地模型在ANSYS软件的嵌入,需解决外源地震动的输入和散射波的波动能量耗散两个问题。传统的方法是利用弹簧阻尼器单元进行模拟,然而,由于等效荷载的计算繁琐以及每个外边界节点的等效物理元件参数的不同,导致可操作性较差。为此,基于黏性边界的相关理论,推导并开发了黏性人工边界单元用于考虑截断外边界的能量耗散及外源地震动的输入,该单元类型依附于已划分场地网格的外侧,易于操作,是对已有节点单元矩阵的继承与修正。

所开发二维黏性人工边界单元如图4所示,通过节点1和2、节点2和3、节点3和4所建立的线单元实现人工边界的施加。以节点2和3组成的单元为例:首先,计算该单元的控制高度h2,即1′和2′所围成高度,即此单元的控制区域,然后,通过阻尼器元件系数计算边界节点在该单元控制高度内对阻尼阵的贡献并将其团聚到相应节点,并将其累计到整体阻尼矩阵中。根据如图4所示的坐标系,二维黏性边界单元阻尼矩阵的表达式如下:

(3)

式中:CBN和CBT分别为阻尼器元件的切向系数和法向系数;hi(i=1～4)为边界节点控制区域的高度(三维情况为控制面积A)。

图4 二维黏性边界单元示意图Fig.4 Sketch of two-dimensional viscous boundary element

另外,在黏性边界单元接口子程序中,可获取该单元内节点的坐标信息及每一荷载步的响应,从而可方便施加式(1)所示节点等效荷载力。

2.3 附加质量法在ANSYS中的实现

由于ANSYS中没有对动水压力进行直接施加的模块,对于式(2)需要运用语言编程实现用单元来模拟附加质量从而达到对动水压力施加在节点上的效果进行替代。其中主要利用如图5所示的结构点质量单元mass21,此单元具有三向的位移及扭转共六个自由度,特别是在静态解中不需考虑任何效应。

图5 mass21单元示意图Fig.5 Schematic diagram of mass21 unit

ANSYS中利用命令流将有限元常规分析进行自动化或者此类分析模型的建立通过编程脚本来完成。此外,还有着参数的交互输入以及页面友好、操作简单等优势。最主要的是其可以避免重复性的工作,尤其是对于本文进行的多种工况运算的模型响应分析有着得天独厚的优点。

2.4 时域模型的建立

珊瑚岛礁场地地震反应时域模型在ANSYS的建立,关键在于珊瑚砂等价线性单元及黏性人工边界单元用户子程序的编制,并通过子程序的编译连接实现新建单元的嵌入。进而,可结合ANSYS成熟的前后处理及强大的求解能力,完成所开发时域模型的建立,详细操作过程列于表1。

表1 所开发时域模型的建立流程详解

3 自由场分析

3.1 材料参数

依据现场原位测试结果[17]及室内动三轴实验数据[4],岛礁场地的分层较为清晰,本文所建模型各分层所用的材料参数列于表2。

表2 场地土层材料参数

3.2 有限元计算模型验证

基于ANSYS的简化验证自由场有限元模型如图6所示,底部为半无限地基,四周为自由边界,选取截断尺寸长宽为100 m×54 m,各土层自上而下的深度分别取为3 m、5 m、5 m、5 m、36 m。同时,建立基于SuperFLUSH/2D的二维有限元模型输入文件,利用该模型计算时,两个方向的响应通过动力计算获得。为了满足简谐波在频率影响内对传播网格尺寸的要求,在验证分析中,将珊瑚礁灰岩处的网格高度划分为2.0 m,上覆各软弱砂层部分的划分为1.0 m。

图6 自由场分析有限元模型Fig.6 Finite element model of free field analysis

3.3 输入地震动的选取

考虑到日本国土与岛礁海域位处版块交界处,结构断裂活动频繁,且海域地震动特性与陆地存在差异。因此,以近20年日本公开的K-NET(Kyoshin network)强震台网中KNG201-206深海台站的强震记录为基础,综合考虑震中距离、震级、波形频率等因素,选择其中3条地震动记录KNG2010102251405、KNG2051801150312、KNG2060710010221作为地震动输入(下面均简称KNG201、KNG205、KNG206),具体的加速度时程及傅氏谱曲线如图7。其中,用KNG201作为验证的输入地震波,水平向与竖直向加速度幅值分别为2.28 m/s2和0.76 m/s2,时间步长取为0.01 s,总持时40 s。

图7 输入的加速度时程曲线Fig.7 Time history curve of input acceleration

3.4 计算结果分析

通过对不同高程处的加速度反应谱及加速度幅值随深度的变化规律来考察所开发模型,其中反应谱选取高程为3 m,-15 m,-51 m 3个代表位置阻尼比为5%的结果对比。

图8为不同方向沿高程的加速度幅值分布情况,可以看出,本文所开发计算模型的自由场2个方向加速度幅值的数值模拟结果与商用软件SuperFLUSH/2D的计算结果整体趋势相同,曲线基本吻合,从数值上来看差别也较小。

图9为2种计算模型在不同高程(3 m、-15 m和-51 m)处2个方向的加速度反应谱对比,新开发模型的计算结果与SuperFLUSH/2D 的计算结果总体趋势及数值基本一致,在2个方向高程的反应谱峰值相差最大分别为0.058 m/s2和0.039 m/s2,进一步验证了本文所开发模型的合理性及有效性。

图8 加速度峰值分布图Fig.8 The distribution of peak acceleration

图9 高程变化的加速度反应谱对比Fig.9 Comparison of acceleration response spectra with different elevation

4 某珊瑚岛礁场地地震反应分析

4.1 工程概况及计算模型

综合崔永圣等[17]地质勘测资料及现场原位试验结果以及陈国兴等[4]给出的岛礁场地地形图,建立了适用于珊瑚砂地基条件下的岛礁场地地震反应分析时域几何模型(图10)。主体部分宽度约为2 500 m,根据工程经验将底部礁灰岩(基岩)部分取至5 000 m,将整体模型进行网格剖分得到图11所示有限元模型,共有85 936个节点,94 867个单元。从地表水平峰值加速度和反应谱两方面求索不同地震动输入条件下珊瑚岛礁场地地震反应特性。

图10 珊瑚岛礁几何模型Fig.10 The geometric model of coral island

4.2 场地加速度放大特性研究

本文给出了在不考虑海水影响时,在输入Denali波的地表PGA放大结果以及输入KNG201、KNG205、KNG206地震波的地表PGA的放大结果(图12)。通过对比可以看出,在无水条件时,本文模型的总体趋势与其相接近,再一次论证了本文开发珊瑚岛礁时域模型的科学性与可靠性。

图11 珊瑚岛礁有限元模型Fig.11 The finite element model of coral island

图12 输入不同地震波条件下主体PGA放大系数曲线Fig.12 The PGA amplification factor curves of main body under different input ground motions

地震波在传播过程中会发生多次不同角度不同方向的折射以及反射,使得岛礁不同地形处的PGA放大系数也会存在明显的差别。从图中可以看出灰砂岛场地的顶角节点相对于其右边凹陷处节点的放大系数变化尤为显著,呈现逐渐递减趋势,两者最大相差接近1倍,原因在于顶部节点不仅处于地表凸起位置还会受到航道的临空效应影响,而中部节点受盆地地形效应的因素影响。此外,从航道底部与两侧岛礁的接触节点到航道中间点处的PGA系数表现出幅度较小的逐渐变大,可能是由于其处于珊瑚礁基岩处受地形以及上覆海水综合因素的影响。

从场地主体PGA变化曲线来看,输入中低频成分较为丰富的KNG201波时反应最小,地表峰值加速度的最大放大系数仅为1.796。然而随着中高频的增多,场地反应变化较大,当输入KNG206波时,地表峰值加速度的最大放大系数可达到4.732。而输入波形成分与KNG206波相近的KNG205波时,表现出的放大效应规律以及数值都较为接近。

将灰砂岛、航道、焦坪(图中简称HS、HD、JP)3个典型场地构造在考虑航道动水压力时输入同一地震动不同幅值(0.10g、0.20g、0.30g)情况下加速度变化系数随高程的变化特征在图13中表示出来。不同高程的上覆珊瑚砂土层对PGA变化系数具有显著影响,输入地震动在珊瑚砂土层向上传播过程中,放大系数变化曲线的数值随高程值变大而呈现出增大的趋势,体现出明显的“砂层放大”效应,但是随着地震动幅值的变大,砂层的放大效应会逐渐减弱。主要是存在地震动与土体弱化之间的相互作用,地震动会带来土体弱化,而土体弱化又反作用于地震动,具体体现在地表剪应变的增大带来土体弱化加剧。可见输入不同的地震波、岛礁场地的地形地貌及航道海水共同影响着PGA系数放大效应。

图13 不同幅值KNG201输入条件下各截面PGA放大系数曲线Fig.13 The PGA amplification coefficient curves of different sections under input KNG201 with different amplitudes

4.3 地表加速度反应谱特性研究

为避免输入地震动幅值对岛礁场地地表反应谱特性研究的影响,工程中常常用动力放大系数β对地震要素反应谱进行研究,即加速度反应谱与峰值加速度的比值,其中Sa取阻尼比为0.05时的加速度反应谱进行计算,具体如式(4)所示 :

(4)

图14给出了灰砂岛、航道、焦坪(图中简称HS、HD、JP)3个典型场地考虑航道动水压力时的地表动力放大系数变化曲线图。总的来看,动力放大系数谱的谱形较瘦长,表明珊瑚岛礁场地在传播地震动时对其中的长周期频谱进行了明显的滤波作用,反应在较短周期内的系数谱值较大。输入KNG201波时,反应谱存在多峰特性且峰值向短周期方向变化,在大于0.55 s后总体谱值相近。当输入KNG205波、KNG206波时,反应谱存在单峰特性且都表现为在短周期内的瘦长状,峰值周期与输入地震动相近,在0.2 s左右时谱反应特别明显。建议建设岛礁重大工程的卓越周期应小范围规避0.2 s及0.55 s,避免发生较大地震破坏。

5 结论

本文基于UPFs二次开发灵活性的特点,结合等价线性法及黏性边界相关理论开发新的单元,基于附加质量法,建立了珊瑚岛礁场地地震响应时域分析模型,通过自由场以及已有成果对比验证其可靠性,并将模型应用于国内某实际岛礁场地,对不同输入地震动下场地反应进行数值分析,主要结论如下:

(1) 对于同一输入地震动,有无航道海水两种工况的场地模型整体放大系数变化趋势大体相近,对结果影响较小。

(2) 灰砂岛场地的顶角节点到本区域凹陷处节点的放大系数变化尤为显著,呈现逐渐递减趋势,两者最大相差接近1倍。此外,航道处的放大系数较其他地形处小,从航道底部两侧节点到航道中间节点处的PGA系数表现出幅度较小的逐渐变大。

(3) 不同高程的上覆珊瑚砂土层对PGA变化系数具有显著影响,在高程-15 m土体剪切模量变化处发生放大系数的突变。放大系数变化曲线的数值随高程值变大而呈现出增大的趋势,体现出明显的“砂层放大”效应,但是随着加大输入地震动的强度,砂层的放大效应会逐渐变弱。

图14 不同输入地震动条件下主体地表动力放大系数曲线Fig.14 The dynamic amplification coefficient curves of surface under different input ground motions

(4) 从场地主体PGA变化曲线来看,输入中低频成分较为丰富的KNG201波时反应最小,地表峰值加速度的最大放大系数仅为1.796。然而随着中高频的增多,场地反应变化较大,当输入KNG206波时,地表峰值加速度的最大放大系数可达到4.732。而输入波形成分与KNG206波相近的KNG205波时,表现出的放大效应规律以及数值都较为接近。

(5) 输入KNG201波时,反应谱存在多峰特性且峰值向短周期方向变化,在大于0.55 s后总体谱值相近。当输入KNG205波、KNG206波时,反应谱存在单峰特性且都表现为在短周期内的瘦长状,峰值周期与输入地震动相近,在0.2 s左右时谱反应特别明显。

综上可见,海水压力、不同的输入地震波波形与幅值、岛礁场地的地形地貌及其本身的抗剪强度共同影响着PGA系数放大效应。建议建设岛礁重大工程的卓越周期应小范围规避0.2 s及0.55 s,避免发生较大地震破坏。