流固耦合地震波动问题的显式谱元模拟方法1)

孔曦骏 邢浩洁 李鸿晶,2)

* (南京工业大学工程力学研究所,南京 211816)

† (中国地震局地球物理研究所,北京 100081)

引言

流固耦合地震波动是地震工程领域中近年来备受重视的一个问题.处理流固耦合问题的一般性做法是认为流体和固体两种介质的耦合作用仅发生在交界面上,因而只需引入界面协调条件以实现波动能量的交换.这种传统分析方法要求流体和固体中的波动以不同的方程进行描述,建立的数值格式不统一,界面耦合处理亦很复杂,应用上十分不方便.陈少林等[1-2]利用两相介质理论发展了流固耦合问题的统一计算框架,将流体和固体分别视为饱和多孔介质退化后的特殊形态,从而通过统一的数值格式实现对流固耦合波动的模拟.流固耦合问题统一计算框架避免了不同模块之间的交互操作,易于实现并行计算,为构建高效率的时空解耦显式算法以实现大规模的流固耦合波动模拟提供了可能.

流固耦合地震波动模拟方法主要有解析法、半解析法、边界元法、时域及频域有限元法、有限差分法等.杜修力等[3]、赵成刚等[4]和王进廷等[5-6]基于集中质量有限元和外推人工边界发展出一套适用于重力坝/拱坝-库水-饱和场地-基岩系统地震反应分析的数值模拟方法,全面探讨了流固耦合界面处理、显式时间积分格式、地震动输入、大型复杂系统建模等问题.李伟华[7]和陈少林等[8-9]对饱和两相介质场地中的波传播与散射理论、显式有限元方法以及人工边界条件进行了进一步研究.李亮等[10-11]、宋佳等[12-13]基于u-p格式的饱和波动方程研究了时域全显式有限元模拟方法以及饱和波动的应力型人工边界条件.梁建文等[14]、巴振宁等[15]、刘中宪等[16]基于多极子快速间接边界元法研究了饱和层状场地对地震波的二维和三维散射.Wang 等[17-19]、Zhao等[20-21]采用半解析法、频域有限元法研究了海水-饱和海床-基岩系统的地震反应以及单桩海上风机、防波堤结构的动力反应.Chen 等[22-23]使用频域有限元法研究了考虑土-结构-海水相互作用时单桩基础以及海堤的频域动力响应.Sun 等[24]通过曲线网格有限差分法研究了考虑海水作用时三维海底场地的地震反应特征.Liu 等[25]、Bao 等[26]、宝鑫等[27-28]基于可压缩流体-弹性固体(考虑非线性)的流固耦合模式,将自主开发的固体介质黏弹性人工边界、流体介质动力人工边界、边界子结构/内部子结构的地震动输入方法、成层介质侧边界地震动输入的一维化时域方法等技术与商业有限元软件相结合,实现了岛礁-海水系统地震反应的有限元分析.陈国兴等[29-30]对典型的实际海域岛礁场地以及海峡场地的二维剖面进行了非线性地震反应分析,探讨了海域场地条件对地震动参数的影响规律.

现有方法中,解析法/半解析法通常只适用于比较规则或由特定函数描述的几何形状,同时由于计算过程中高阶级数项数值很小容易导致溢出计算机精度范围,使得这类方法能够模拟的频率上限不高.有限差分法目前常用高阶差分格式来提高模拟波场的精度并且能够模拟不均匀介质,但是其处理复杂地表地形及介质交界面时往往精度不高并且容易引起稳定性问题,另外它在模拟面波以及界面波方面的效果较差.边界元法能够模拟很宽的频带范围,且近年来发展出多极子技术来提高计算速度,已被成功地应用于大型河谷饱和场地以及城市尺度沉积盆地的地震动模拟,但是其计算模型主要为单一介质的不规则区域或横观各向同性的水平成层及圆弧成层场地.有限元法最大优势在于通过灵活的单元划分来实现对各种复杂几何形状的建模,以及从建立单元到集成为整体计算模型过程的规整性,为开发大型计算程序带来了便利,在应用于场地地震波动模拟时又开发出了集中质量有限元技术以提高计算效率,不过,该方法目前主要使用基于一阶插值函数的四边形单元(二维)或六面体单元(三维),这使得它模拟波场的精度比较有限,高频段模拟结果往往并不准确.

谱元法[31-34]可以看作一种特殊的高阶有限元法,其优点在于使用高阶单元(如四阶以上的单元)进行波动模拟并且从理论上严格地导出集中质量矩阵,成功地解决了目前基于低阶单元(常用一阶单元)的集中质量有限元法模拟波场的精度不高以及现有基于经验的集中质量方法缺乏严格数学基础的问题.Duczek 等[35]证明了谱元法基于节点积分直接导出的集中质量矩阵,与传统基于行和集中、对角元素放大等经验方法得出的集中质量矩阵结果一致.近二三十年来的研究和应用已较为清楚地表明,谱元法是有限元法用于各类波动问题数值模拟时一种最为高效(既能高精度地模拟波场又能灵活地对复杂几何形状进行建模且大规模问题的计算量可接受)的形式.

本文以饱和两相介质波动的Biot 方程及其向理想流体(孔隙率为1 情形)、弹性固体(孔隙率为0 情形)的两种退化形式方程[1-2,32,36]为基础,采用高精度谱元法与多次透射公式人工边界条件[37-40]相结合,发展流固耦合地震波动问题的高效数值模拟方法.通过与传递函数法[41]以及集中质量有限元法的结果进行比较,验证本文方法的高效性和地震波动模拟结果的宽频带特性.

1 流固耦合问题的控制方程

考虑如图1 所示流固耦合问题,该问题模型由理想流体、饱和多孔介质和弹性固体三种介质组成,亦可以退化为只含有两种介质之间耦合的情形.

图1 流固耦合地震波动问题示意图Fig.1 Schematic diagram of seismic wave motion problem in a fluidporoelastic-solid model

本文流固耦合地震波动问题的控制方程由Biot 所提出的饱和多孔介质波动方程及其向理想流体、弹性固体的两种退化形式共同组成.这样组合的优点是饱和波动的求解模式可以直接应用于理想流体和弹性固体,并且饱和多孔介质之间的界面连续条件亦可以转化为任意两种介质界面情形,接下来分别进行介绍.

考虑二维情形,当场变量为固相骨架和液相流体的位移时,Biot 流体饱和多孔介质波动方程为

式中,u,U分别为固相骨架位移、液相流体位移的列向量,u=(ux,uy)T,U=(Ux,Uy)T,变量上方一个、二个圆点分别表示对时间求一阶、二阶导数,上标“T”表示矩阵或向量转置;Ls和Lw为微分矩阵,D1,D2和D3为介质的力学参数,Lw=[∂/∂x,∂/∂y],D2=Q,D3=R

这里β为孔隙率,ρs为固相骨架密度,ρw为液相流体密度,参数b=β2μ0/k0,μ0为孔隙流体的动黏度系数,k0为达西渗透系数;A,G,Q,R为Biot 常数,可通过如下表达式计算

其中λ和μ为固相骨架的拉梅常数,和M为表征饱和多孔介质压缩性的常数.在实际应用中,相关常数可以通过试验来测出.

Biot 方程(1)当孔隙率为1 时退化为理想流体方程.此时孔隙率β=1,以及常数α=1,且由于不存在固相骨架,动黏度系数为零,导致参数b=0,固相模量亦置为0.于是Biot 常数为:A=0,G=0,Q=0,R=M.方程(1)改写为α

Biot 方程(1)当孔隙率为0 时退化为弹性固体的波动方程.此时孔隙率β=0,以及常数α=0,且由于不存在固相骨架,动黏度系数为零,导致参数b=0.于是Biot 常数为A=λ,G=μ,Q=0,R=0.方程(1)改写为

类似上述由Biot 方程导出理想流体和弹性固体方程的做法,不同介质交界面连续条件同样可以由两种不同力学参数饱和多孔介质的交界面这种一般情形,分别退化到理想流体–饱和多孔介质、饱和多孔介质-弹性固体、理想流体–弹性固体三种交界面情形,如图2 所示.

图2 流固耦合问题的不同介质交界面Fig.2 Different interfaces in the fluid-poroelastic-solid problem

两种不同力学参数饱和多孔介质交界面的连续条件为:

(a) 交界面两侧介质的法向总应力相等;

(b) 交界面两侧固相骨架的切向应力相等;

(c) 交界面两侧液相流体的压强相等;

(d) 交界面两侧固相骨架的法向、切向位移相等;

(e) 交界面两侧的法向运动连续.

上述界面连续条件可用表达式分别表示为

式中,上标“A”和“B”指示交界面两侧的不同介质,下标“N”和“T”分别指示变量的法向分量和切向分量;σ为固相骨架应力,τ为平均孔压,τ=-βP,β为孔隙率,P为孔隙流体压力.

式(6)所给出的饱和多孔介质A-饱和多孔介质B 这种一般情形下的界面连续条件向理想流体-饱和多孔介质、饱和多孔介质-弹性固体、理想流体-弹性固体三种具体情形下交界面连续条件的退化,可按照由方程(1)退化到方程(4)或方程(5)的基本规律以及由图2 所给出的参数取值方式进行.不过,在退化过程中式(6)的一些项需要随着固相或液相成分的消失而舍去,相关细节可参见文献[1-2,5,36],本文不再赘述.

综上,内域波动方程式(1)、式(4b)和(5a)和界面连续条件方程(6)组成了本文流固耦合地震波动问题的控制方程.

2 流固耦合地震波动的谱元分析方法

为提高流固耦合地震波动问题的求解精度和计算效率,不同于现有的集中质量有限元法,本文采用能够天然地导出集中质量高阶单元的谱元法来求解方程式(1)、式(4b)、式(5a)和式(6)组成的系统在地震作用下的动力反应.谱元法可以看作一种特殊的高阶有限单元,因此其基本步骤与有限元法一致,区别仅在于它采用的是基于谱插值方式构建的高阶单元.对上述流固耦合方程系统的空间离散,由于退化关系的存在,使得只需详细探讨对于饱和两项介质波动Biot 方程(1)和一般情形界面连续条件方程(6)的谱元离散,然后根据退化规则即可退化到其他情形.本节首先论述内域波动方程的谱元离散及其时域逐步积分格式,然后在此基础上介绍界面连续条件的离散化表示,最后探讨人工边界条件与地震动输入,这样形成一套流固耦合地震波动的谱元分析方法.

2.1 谱元离散方程

首先考虑饱和波动方程(1)的谱元离散.与有限元法类似,谱元离散亦需在原方程的等效积分弱形式上进行.对方程(1)的每一项乘以任意的测试函数v或V,对空间变量在全域上进行积分;对于方程右边关于空间导数的各项进行分部积分,利用格林公式将域积分转化为边界积分,引入自由地表条件使得边界积分项等于零而从方程中消失.得到等效积分弱形式方程如下

在弱形式方程(7)中,将计算模型Ω划分为各个谱单元Ωe的组合,则全域积分转换为各个谱单元积分的和.将各个谱单元的积分、微分运算转化到[-1,1]×[-1,1]的标准参考单元Λ上进行.

根据Galerkin 离散化法则,对以单元形式表示的弱形式方程进行数值离散,单元内的各个场变量和测试函数均采用统一的谱插值模式来构建单元插值函数.谱单元的插值函数在标准参考坐标Λ(ξ,η)∈[-1,1]×[-1,1]内进行定义,其中单元场变量的插值函数表达式为

式中,ux,i,uy,i,Ux,i,Uy,i(i=1,2,···,ng)为待求解的离散网格节点上的场函数值,Ni(ξ,η) 为单元内第i个节点的形函数,ng为谱单元的节点数目.

形函数表达式为

式中,ngllξ和ngllη分别表示ξ,η方向的单元节点数目.lm(ξ),ln(η)均为Lagrange 插值基函数,ξj,ηj为Gauss-Lobatto-Legendre (简称GLL)高精度数值积分节点[35,39].本文采用常用的四阶谱单元,其GLL 节点为[-1,-0.654 653 670 7,0,0.654 653 670 7,1].

根据上述谱单元场函数的定义,可以完成以单元形式表示的弱形式方程的数值离散,然后将各个单元离散方程集成为总体离散方程.于是,得到饱和两项介质波动Biot 方程(1)采用谱元法进行数值离散的运动方程为

式中,u和U为单元节点位移的列向量,Ms为固相质量矩阵,Mw为液相质量矩阵,C为流固耦合阻尼矩阵(用于计算固液两相之间的黏性阻力),Kss为固相刚度矩阵,Kww为液相刚度矩阵,Ksw和Kws为代表固液两相之间本构力的耦合矩阵;Fs为作用在固相上的外荷载矩阵,Fw为作用在液相上的外荷载矩阵,FRs和FRw分别为单元之间作用在固相和液相上的相互作用力,主要体现在不同介质的分界面上,若单元在同一介质内时,该力可被忽略.上述各个整体矩阵由相应的单元矩阵集成得到

式中,N为形函数矩阵,J=J(x,y;ξ,η)为雅可比矩阵,|J|为雅可比行列式,矩阵Bs=LsN,矩阵Bw=LwN,为由方向导数组成的矩阵,n为沿外边界外法线的方向矢量,相关表达式为

根据第1 节所揭示的波动微分方程退化关系,按照对应的参数取值规则,可以方便地由这里的饱和波动谱元离散运动方程(10)分别退化得到理想流体、弹性固体的离散运动方程.因此在计算程序中,只需适当地调整饱和多孔介质的单元特性矩阵及总体特性矩阵中的相关参数取值,就可以得到理想流体、弹性固体的单元特性矩阵及总体特性矩阵.这种统一的处理方式首先由陈少林等[1-2,36]提出,使得流固耦合问题的离散化建模过程大为简化.

2.2 时域逐步积分

传统有限元法中的总体质量矩阵Ms,Mw和流固耦合阻尼矩阵C均为非对角矩阵,显式的时步积分格式需要对质量矩阵(Ms,Mw)进行人为的行和集中等措施才能实现,但流固耦合阻尼矩阵C通常不做处理,依然为非对角矩阵.因此,基于集中质量有限元的流固耦合地震波动问题一般选取具有一阶计算精度的中心差分-单边差分法或者形式繁琐但具有二阶计算精度的中心差分-Newmark 常平均加速度法.

谱元法基于GLL 积分直接导出集中质量矩阵.文献[35]从理论角度严格证明了谱元质量模型与传统基于行和集中、对角元素放大等经验方法得出的集中质量矩阵在本质上是一致的,但具有坚实的数学基础.然而很少有研究关注到Biot 方程中的流固耦合阻尼矩阵C在应用谱元法进行GLL 积分后同样为对角矩阵,这一点可以从式11(a)、式11(b)与式11(c)的对比中直接看出: 质量矩阵和阻尼矩阵的积分公式只有系数上的区别,进行GLL 数值积分后皆为对角矩阵.

如何理解流固耦合阻尼矩阵C为对角矩阵的物理本质对波场问题的数值计算具有重要意义.有限元的行和集中质量矩阵被认为是忽略了单元内加速度的变化,即假定单元内惯性力为常量.基于GLL 积分的流固耦合阻尼矩阵C为对角矩阵,则可以类似地理解为忽略了单元内固、液两相速度之差的变化,即假定单元内由两相速度差导致的黏性阻力为常量.

考虑到本文中总体质量矩阵Ms,Mw和流固耦合阻尼矩阵C均为对角矩阵,可以采用计算精度为二阶且形式简洁的中心差分法来处理谱元离散方程(10)中的时间导数

其中,上标p表示t=pΔt时刻.将式(13)代入谱元离散方程(10),得到时域逐步积分格式

其中,u和U表示被计算节点的固、液相位移;ue和Ue表示与该计算节点相关单元内节点的固、液相位移.时域逐步积分格式(14)是关于整个谱元离散模型的时间步进格式.式中含有模型的整体特性矩阵,实际编程时这样做很不经济,且对于大规模问题往往难以实现.考虑到这里总体质量矩阵Ms,Mw和流固耦合阻尼矩阵C均为对角矩阵,所以式(14)可以写成关于各个节点运动ux,uy,Ux,Uy的完全显式的时间积分格式.以编号为k的节点固相x向位移ux为例,其显式逐步积分格式可表示为

其中,nc为与节点k相关的节点的数目.

2.3 界面连续条件处理

图3 为两种不同力学参数饱和多孔介质的交界面在谱元离散情况下的相互作用力示意图.在图3中,介质交界面一条单元边上有ngll个不等间距的计算节点(图中圆点ngll=5);FRs和FRw分别为单元之间作用在固相和液相上的相互作用力;下标中的N和T分别表示交界面节点处作用力的法向分量和切向分量;A和B及对应的上标分别表示图中上下介质交界面任意一个计算节点;e1,e2,e3,e4 为单元编号.

图3 界面连续条件的谱元离散Fig.3 Spectral-element discretization of the continuous condition on the interface

结合式(6)交界面节点的应力和位移连续条件,可得到对应的相互作用力与位移连续条件,即

将上述连续条件与节点A和B的固液两相动力平衡方程(参考式10)结合,即可求得图3 中的8 个相互作用力.

需要说明的是,当图3 中的饱和多孔介质向理想流体或弹性固体退化时,式(18)的界面连续条件亦会发生变化[1-2,5,36].将不同介质交界面的连续条件与对应节点的动力平衡方程相结合,即可求出对应的界面相互作用力.

2.4 人工边界条件与地震动输入

本文使用廖振鹏等[42-43]提出的多次透射公式(multi-transmitting formula,MTF)作为局部场地模型的人工边界,该边界具有简单、通用、精度较高的特点.在饱和多孔介质的场地模型中,多次透射边界需要同时应用于边界节点固相位移和液相位移的x,y两个方向分量之上,即

式中,N为透射阶次,为二项式系数,以为例,其表示参考点j的x向固相位移在t=pΔt时刻的波场值,由下式定义

其中,Δt为时间步距,ca为人工波速.这里人工波速指的是MTF 边界所使用的计算波速参数,理论上它可以设定为任意值,在实际计算中通常取为等于或略大于介质物理波速的某个值,本文取为固相纵波波速和液相纵波波速中的大值.

式(20)中参考点的位移通常不在波场的计算节点上,需要通过空间内插来获得,而插值形式又与单元场函数的形式有关.采用本课题组开发的MTF 谱元离散格式实现[38-40],其特点为精度高且插值基函数与单元形函数一致.以为例,对应的MTF 谱元离散格式如下

表1 透射边界(MTF)谱元离散格式中插值点的局部坐标Table 1 Local coordinates of the interpolation points in the spectral-element formulation of multi-transmitting formula

针对外源问题,地震动的输入需要通过人工边界来实现,这也是MTF 的一个优势,应力型边界或PML 边界都不便于实现在人工边界上的地震动输入.MTF 作为位移型人工边界,能方便地对全波场进行入射波场与散射波场的分离.图4 为谱元模型中多次透射公式人工边界条件和地震动输入示意图.其中,等号右边的节点代表之前某时刻求得的全波场解;加号右边节点代表已知的入射波场,本文通过传递矩阵方法(transfer matrix method,TMM)[41-44]求得;全波场解减去入射波场解,则得到之前某时刻的散射波场,即图中加号左边圆形实心节点和矩形实心节点,空心圆点表示MTF 的空间插值参考点,其值由式(21c)插值而得.求得MTF 各参考点的位移后,根据式(21a)可求得p+1 时刻的边界节点位移,再结合波动方程算出的其余内节点位移,则得到p+1时刻完整的全波场解.

图4 多次透射公式人工边界条件与地震动输入Fig.4 Artificial boundary condition and seismic wave input based upon multi-transmitting formula

3 算例验证

选取了二维水平成层、不规则层状界面和任意形状界面的三种场地模型,介质皆为理想流体-饱和两相土体-弹性固体,具体材料参数见表2.二维模型的侧边界和底面采用多次透射人工边界(MTF,本文取为一阶)来模拟无限域,上表面为自由地表,侧边界的入射波场采用传递矩阵方法(TMM)[41,44]得到的自由场,底面的入射波场则为图5 所示的脉冲波.

图5 输入样条脉冲波Fig.5 The input wave of a spline impulse

表2 介质参数表Table 2 Parameters of medium

3.1 水平层状场地

图6 为二维水平成层场地模型,下述所有工况的谱单元尺寸皆满足 ΔxE＜λmin(λmin为三种介质中的最短波长),时间步距为 Δt=0.000 2 s,满足稳定性条件.

图6 二维水平成层场地模型Fig.6 Two-dimensiona model of horizontal layered site

本文方法(S E M) 除了与传递矩阵解析解(TMM)作对比外,还与有限元法(FEM)和统一框架计算理论相结合的结果作对比.其中,有限单元尺寸为 Δx＜λmin/10,对应的时间步距为Δt=0.000 05 s(Δx=1 m)或 Δt=0.000 2 s(Δx≥2 m).测点A,B,C,D位于底部边界和介质交界面的中点处.

图7 为P 波垂直入射情况下各测点的竖向位移.从图7 中不难看出,本文方法与传递矩阵解析解在所有测点的位移都完全重合.有限元法与传递矩阵解析解在测点A和测点B的位移几乎重合(见图7(a)～图7(d)),但在C点和D点的位移则存在一定误差(见图7(e)～图7(h)),尤其是C点的饱和土固相位移和液相位移(见图7(e)和图7(f)).值得注意的是,谱元模型的计算节点总数为2499 个,而有限元模型的计算节点总数为14 883 个,这意味着谱元模型用比有限元模型少约83%的计算节点数,达到了比有限元模型更高的精度.

图7 各测点的竖向位移Fig.7 Vertical displacement at measuring points

通过上述时程曲线观察得出的谱元法模拟结果要优于有限元法模拟结果,其本质原因在于高阶单元比低阶单元能够在更宽的频率范围内精确地模拟波场.数值离散网格对于波传播过程而言相当于低通滤波器,高阶方法的优势在于它能够达到的上限频率更高.

为定量地考察这种频率范围的不同,这里以解析解作为参照,探讨谱元法和有限元法数值解相对于解析解的谱比的差异.谱比计算公式如下

其中,S(f)为谱比,U(f)为频谱幅值.

图8 绘出了在不同网格尺寸条件下谱元法和有限元法计算出的C点固相、液相位移相对于解析解的谱比曲线.若谱比值越接近于1,表示该频率处数值解越精确;反之,若偏离程度越大,则表示该频率处数值解的误差越大.图中结果清楚地表明,谱元法能够在比有限元法宽得多的频带范围内给出精确的模拟结果.本算例中谱元法(四阶单元)能够精确模拟的波动频率上限约为13 Hz,而有限元法(一阶单元)只能达到约5 Hz;谱元法在0～13 Hz 范围内模拟精度极好,而有限元法在0～5 Hz 范围内仍存在少量误差.此外,图中结果还显示出网格尺寸对波动模拟精度的影响规律,比较FEM 网格尺寸分别为1 m,2 m,4 m 的谱比曲线可知,网格尺寸越小,谱比误差越低,即模拟精度越高.但是,缩小网格尺寸对模拟精度的提高远不如提高单元阶次所带来的收益.总之,谱元法是适用于波动数值模拟的一种特殊的高阶有限元,本算例表明高阶单元的最大优势在于能够在宽得多的频带范围内实现对波场的精确模拟.

图8 不同方案情况下C 点的位移谱比Fig.8 Spectral ratio of displacement at point C with different scheme

需要说明的是,这里主要体现的是空间离散方式对波场频散的影响,时步积分方法亦能产生较大影响[45-46].统筹考虑时空离散方法对数值模拟的频散影响具有重要研究意义.鉴于本文的主题是使用显式谱元模拟方法求解大型场地的地震动反应,这里不再深入讨论.

3.2 不规则层状界面场地

图9 为不规则层状界面场地模型,P 波从模型右下角以 θ=10°的倾斜角斜入射至波场中,介质分界面处的斜坡坡角为 α=14°.谱单元尺寸为ΔxE=2～20 m,时间步距为Δt=0.000 1 s.测点A,B,C,D位于介质交界面的角点处.

图9 不规则层状界面模型Fig.9 Irregular layered interface model

图10 为各测点的位移时程图.其中,时程曲线名称对应的符号u和U分别表示固相位移和液相位移,上标b,s,w分别表示基岩、饱和土、理想流体,下标x,y分别表示测点的横向与竖向.从图10中可以看出,由于P 波斜入射以及存在不规则界面的原因,各测点的横向位移出现了不同的波动.基岩和饱和土中的测点(A-D)由于处于坡度较小的坡面上,其横向位移主要来源于入射波的水平分量,幅值较小,见图10(a)～图10(f).然而,在理想流体中,测点C与测点D的液相位移出现了较大的水平分量(见图10(g)和图10(h)).这主要是因为这两个测点位于不规则界面的角点处,且理想流体的厚度相对较小,容易导致散射波的叠加,并在界面之间来回反射.

在图10 中,测点C和测点D的位移曲线没有观察到一致性,主要是因为界面的位移连续性是针对界面的法向和切向的,而不是横向与竖向.在实际编程过程中,界面测点的法向方向为界面上该点所在角的角平分线,切向方向则垂直于法向方向.下面将测点C和测点D的位移沿着法向方向进行分解.根据理想流体和饱和土的位移法向连续条件,有

图10 各测点的位移Fig.10 Displacement at measuring points

其中,下标N表示测点的位移矢量沿着法向的分量,1 表示理想流体,2 表示饱和土,nx和ny表示单位法向矢量在横向和竖向的分量.

同理,对于测点A和测点B,根据弹性固体和饱和土的位移法向连续条件,有

其中,2 表示饱和土,3 表示弹性固体(基岩).

图11 为结合式(23)与式(24)计算得到的各测点法向位移.从图11 中可以看出,弹性固体与饱和土分界面上的测点A和B以及饱和土与理想流体分界面上的测点C和D皆满足位移法向连续条件.

图11 各测点的法向位移Fig.11 Normal displacement at measuring points

图12 给出了土体(基岩与饱和土)表面各计算节点固相位移峰值与输入波幅值的比值,即|max(ut)/AP|.这里 max(ut) 为各节点的位移时程峰值,AP为输入P波的幅值,图中分别表示固相横向位移与固相竖向位移.

图12 土体表面的位移幅值Fig.12 Surface displacement amplitude of soil

从图12 中可以看出,由于土体表面整体上比较平稳,且入射角较小,各节点的固相水平位移幅值较小,而竖向位移峰值较大;基岩表面的竖向位移幅值比饱和土表面的竖向位移幅值要大,这主要是由传播介质的不同以及地形变化所导致的.

3.3 任意形状界面场地

图13 为任意形状界面场地模型,P 波从模型右下角以 θ=5°的倾斜角斜入射至波场中.谱单元尺寸为 ΔxE=3～10 m,时间步距为 Δt=0.000 2 s.测点A,B,C,D位于介质交界面的拐点处.

图13 任意形状界面模型Fig.13 Arbitrary shape interface model

结合式(23)和式(24),图14 给出了各测点的法向位移.从图中可以看出,不同分界面上各测点的法向位移皆满足位移连续条件.

图14 各测点的法向位移Fig.14 Normal displacement at measuring points

图15 给出了饱和土表面各计算节点固相位移峰值与输入波幅值的比值.从图中可以看出,由于饱和土表面起伏较大,部分节点的固相横向位移的成分较大,接近0.5;固相竖向位移峰值随节点的起伏分布有一定变化,但都在1.3 左右.

图15 饱和多孔介质表面的位移幅值Fig.15 Surface displacement amplitude of saturated porous medium

3.2 节与3.3 节两个模型的计算结果证明了本文方法应用于不规则界面复杂介质模型(理想流体-饱和多孔介质-弹性固体)的适用性,可用于大型河流湖泊、海洋区域工程、河谷盆地及大坝工程等场地的地震动数值模拟研究.

4 结论

如何对涉及到理想流体、饱和多孔介质、弹性固体三种介质的流固耦合地震波动问题进行高效的数值模拟是大型河流、湖泊以及海洋范围内各种重要工程的抗震设防迫切需要解决的关键问题.本文基于饱和两相介质波动的Biot 方程及其向理想流体和弹性固体的两种退化形式方程,采用兼具波场模拟的高精度特性(源于谱插值)与复杂区域建模灵活性(源于与有限元相同的计算模式)的谱元法,结合简单高效的廖氏透射边界条件,建立了一种流固耦合地震波动问题的高效数值模拟方法.通过含有水平界面、不规则层状界面和任意界面场地的数值模型算例,证明了本文方法的可行性,并且揭示出与目前大量使用的有限元法相比,该方法能够在宽得多的频带范围内实现对流固耦合场地中地震动的有效模拟.

本文方法除了具备现有流固耦合地震波动有限元模拟方法中完全显式数值计算格式的一般特征之外,还具有从理论上严格地导出集中质量矩阵和高精度谱单元能够很好地平衡地震动模拟精度与大规模数值计算效率这两方面的优势.文中讨论主要针对二维问题进行,可以方便地推广至三维情形.

此项研究对于数值模拟得到的地震动的宽频带特性分析具有理论意义,并且对于大型河流、湖泊以及海洋区域的场地地震反应分析具有工程应用价值.后续工作可以继续研究大尺度谱单元对局部不均匀介质的有效建模、高精度人工边界、地震动输入方式以及时步积分方法等对地震动不同频段特性的影响,或将各向异性介质[15,47-48]纳入统一框架理论,更重要的是考虑饱和两相介质的非线性并使用并行计算技术开发实用的三维计算程序,为解决实际工程问题服务.