黏弹性多层介质中SH 波动的一种吸收边界条件1)

吴利华 赵 密 杜修力

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

引言

刚性基岩上卧成层土是工程结构地基的普遍存在形式[1].由于工程结构尺寸远小于场地尺寸,从数学建模的角度可以把场地看作是无限域.采用如有限元法、有限差分法等数值方法分析工程结构动力响应时,通常人为弓入人工边界将无限域截掉,取包含感兴趣区域的有限域作为数值分析对象.吸收边界条件应该被施加在有限域的人工边界上来反映被截掉的无限域和有限域间的动力相互作用[2].对结构进行非线性分析等情况下,需要吸收边界条件是时域方法.已有的时域吸收边界条件大多是基于线弹性介质建立的,如局部吸收边界条件[3-11]和高精度吸收边界条件[12-20].但是波在土中的传播总是伴随着能量的衰减,因此相比于把无限域看作线弹性介质,将其考虑为黏弹性介质更接近工程实际[21].

特征线法[22]、传递矩阵法[23]、高阶薄层法[24]以及阻尼层法[25-29]等方法被用来建立黏弹性成层介质的频域吸收边界条件.但是频域吸收边界条件有两点局限性,一是将频域结果转化为时域结果需要复杂耗时的卷积计算[30],二是频域方法不适用于非线性时程分析.针对黏弹性介质建立的时域吸收边界条件目前主要有两类:局部吸收边界条件[31-32]和阻尼层形式的吸收边界条件[33-39].王振宇和刘晶波[31]基于黏弹性理论构造了局部吸收边界条件的一维形式,证明了在弱阻尼时其近似等效于弹性一维吸收边界条件,并将其推广得到二维吸收边界条件[31].但是对于强阻尼介质,将基于线弹性介质建立的吸收边界条件未加修正地应用于黏弹性情形会导致数值失稳或模拟精度降低等问题[33].Badry 等[32]建立了半空间黏弹性介质的局部吸收边界条件,其中黏性考虑为与质量成正比的瑞利阻尼形式.上述两种局部时域吸收边界条件都是基于半空间模型建立的,将其用于刚性基岩上卧成层介质模型时模拟精度较低.阻尼层形式的时域吸收边界条件是近年的研究热点,其通过在人工边界外设置具有耗能作用的缓冲区来吸收透射波.黏弹性介质的阻尼层吸收边界条件有:完美匹配层[34-35]、Caughey 吸收层[36]、增加阻尼吸收层[37]以及阻尼溶剂抽取法[38]等.阻尼层的层厚度和耗能参数选择不当时,会导致模拟精度低或出现数值失稳问题.

为了提高计算的精度和效率,高精度时域吸收边界条件被广泛研究,它是针对具体的无限域定解问题,基于解析或者半解析方法建立的,因而具有精度高的特点.目前该类方法主要被用于建立线弹性介质的时域吸收边界条件.但是,线弹性成层介质的时域高精度吸收边界条件的构建方法可以被借鉴用来建立黏弹性成层介质的时域高精度吸收边界条件.线弹性成层介质的高精度时域吸收边界条件的构建方法主要有两类,一类是对时域下的场方程运用算子分离法得到时域吸收边界条件.Liu 等利用该方法建立了相关时域吸收边界条件来模拟瞬态标量波[12]和矢量波[14]在成层介质中的传播.然而该方法得到的时域吸收边界条件一般会存在稳定性问题[39].另一类是通过如薄层法[40]和比例边界有限元法[41]等半离散方法得到成层介质的频域动力刚度,再通过时域化方法[13,15,19]建立时域吸收边界条件.单层介质的频域动力刚度是标量形式,多层介质的频域动力刚度是矩阵形式.为了将频域动力刚度时域化,可以通过优化方法将标量动力刚度表示为有理函数近似,进一步有理函数近似可以等价为时域集中参数模型[15];也可以通过数学计算将标量动力刚度表示为连分式展开,再结合辅助变量技术将连分式展开转化为时域吸收边界条件[13,19-20].由于优化方法对于矩阵问题很难收敛,一般有理函数近似不能用来描述多层介质的动力刚度.但是标量连分式可以被扩展为矩阵形式用来描述多层介质的动力刚度[19].基于连分式的时域吸收边界条件对于单层介质是精确且稳定的,对于多层介质是高精度且不易失稳的[19].

本文将无限土看作达朗贝尔线性黏弹性介-质[21,42],它是通过在运动方程中增加阻尼项来描述介质的黏性.尽管达朗贝尔模型相对简单,没有从本构层面描述介质,但由于它包含很少的物理参数,在实际应用中具有优势.考虑到连分式用于构建高精度吸收边界条件的优越性,本文将达朗贝尔单层介质动力刚度表示为连分式展开,类似文献[19]的做法,再将连分式扩展为矩阵形式,建立基于连分式的达朗贝尔黏弹性单层和多层介质稳定的时域高精度吸收边界条件.

1 问题描述

如图1 所示,工程结构的地基为刚性基岩上卧达朗贝尔线性黏弹性成层介质.用有限元法分析工程结构的二维出平面SH 波动响应,需要人为弓入两条直线型人工边界将整个域分为有限域和左右两个半无限域.有限域可以包含非线性,时域吸收边界条件应该被施加在有限域的人工边界上来表示被截掉的无限域和有限域间的动力相互作用.

图1 刚性基岩上卧达朗贝尔黏弹性成层介质Fig.1 D’Alembert viscoelastic multilayered media overlying rigid bedrock

吸收边界条件根据无限域推导获得.达朗贝尔黏弹性介质和线弹性介质两者控制方程的区别是达朗贝尔黏弹性介质比线弹性介质的运动方程多了一项耗散项,其余控制方程两者相同.达朗贝尔黏弹性介质的波动方程为

无限域的上下表面物理边界条件分别为

其中,(x,y)表示笛卡尔坐标,u表示出平面位移,τy是应力,G是剪切模量,ρ 是质量密度,η=αρ,α 和η 是阻尼参数,变量上方的点表示对时间求导.假定每一层的材料参数是常量,不同层的材料参数可不同.另外,位移u在不同土层的分界面以及人工边界处是连续的.土层的初始状态为静止.

2 半无限域的频域动力刚度

本节结合半离散法和模态分解法获得模态空间下半无限域的频域动力刚度,并显式表示出单层介质的频域动力刚度.

2.1 物理空间下的频域动力刚度

将半无限域沿着土层y方向离散,同时考虑式(2)的物理边界条件及介质交界面的连续条件,得到半离散的有限元方程

以及人工边界上的等效结点载荷

半无限域人工边界的结点力-位移关系为

顶标～表示频域下的值,S∞表示物理空间下半无限域的频域动力刚度.

2.2 模态空间下的频域动力刚度

为了将物理空间下半离散的有限元方程变换到模态空间下,弓入广义特征值分解

其中,Φ是特征向量矩阵,Λ2是对角形式的特征值矩阵,其对角线上的元素为考虑半无限域模态缩减,有:Φn表示包含前n阶特征向量的矩阵,表示由前n阶特征值组成的对角阵,I表示单位阵.

根据模态分解法,得到人工边界上结点力、位移在模态空间下和物理空间下的关系为

其中,q和p分别表示模态空间下的位移和力.

将式(7)和式(8)分别代入式(3)和式(4),得

将式(7)和式(8)代入式(5),得到模态空间下人工边界的结点力-位移关系为

其中,i 是虚数单位,ω 表示圆频率.

3 频域动力刚度的连分式展开

本节提出一种新的连分式精确逼近单层介质模态空间下的频域动力刚度,并将连分式扩展用于表示多层介质的动力刚度.

3.1 标量连分式

式(12)表示的动力刚度系数还可以将其表达为

其中,z=iω/ωn,β=α/ωn.

可以将其写成连分式

图2 连分式的收敛域Fig.2 The convergence domain of the continued fraction

为了直观看出连分式的收敛率,定义总绝对误差其中,CF表示式(14)的连分式,N表示频率步数.总绝对误差与连分式阶数的关系如图3 所示,其中,步距为0.1,N为10,β 为0.05,可以看出连分式是呈指数收敛的.

图3 连分式的收敛率Fig.3 The convergence rate of the continued fraction

式(14)的连分式可以写为如下形式

当j是奇数时

当j是偶数时

3.2 矩阵连分式

由式(15)可得出单层介质的模态动力刚度为

当j是奇数时

当j是偶数时

其中,上标-1 表示矩阵求逆.

连分式(16)可以被扩展应用于多层介质.需要注意的是,单层介质的g0和h0的值不适用于多层介质.多层介质的g0和h0可以参考文献[19]的做法得到.

根据式(9)～式(11),得到模态动力刚度方程

将式(16)代入式(22),并令iω 的零次项和二次项分别为0,得到关于g0和h0的如下方程

式(23)是黎卡提方程,可利用Matlab 内置函数care直接求出g0和h0.显然,式(23)的解包含单层介质的解式(19).

式(16)～式(21)表示的连分式对于单层介质是精确且无条件收敛的.将式(23)替换式(19)后,连分式可以用于表示多层介质的动力刚度.

4 基于连分式的时域吸收边界条件

本节通过弓入辅助变量的方式将连分式时域化,建立能与有限元法无缝耦合的时域吸收边界条件.

4.1 模态空间下的时域吸收边界条件

将式(24)和式(25)变换到时域,得

其中,

式(26)即为模态空间下的时域吸收边界条件.5.1 和5.2 节分别通过稳定性理论说明吸收边界条件在单层和多层介质中的稳定性.

4.2 物理空间下的时域吸收边界条件

将式(7)和式(8)代入式(26),得到物理空间下的吸收边界条件

4.3 有限元法和吸收边界条件的耦合

式(27)表示的时域吸收边界条件可以和有限元法无缝耦合.图1 所示的有限域的有限元方程为

其中,下标i 表示除了人工边界的有限域部分,下标b 表示人工边界部分,fi是施加在有限域的动力载荷,K,C,M分别是通过标准四边形等参元、双线性形函数以及四结点高斯积分得到的刚度阵、阻尼阵和质量阵.

组装式(27)和式(28),得

式(29)可以通过标准的时间积分算法求解.需要说明的是,本文提出方法的精度不受有限域材料特性的影响,有限域的介质和结构可以考虑其真实的非线性效应.

5 数值算例

通过单层介质和五层介质的数值算例来考察本文建立的吸收边界条件的精度和稳定性.

5.1 单层介质

如图4(a)所示,总高为40 m 的刚性基岩上卧达朗贝尔黏弹性单层介质的场地内嵌一个矩形孔洞,其位于地表以下8 m,尺寸为6 m×8 m.均匀场地的材料参数如图4(a)中所示.A点施加如图5 所示的狄拉克脉冲点载荷,载荷的频谱范围约为0～55 rad/s,其表达式为

其中,t表示时间,Z(x)=x3H(x),H(x)表示Heaviside函数.输出B点的位移时程.

图4 单层数值算例Fig.4 The numerical example of the single-layered medium

图5 狄拉克脉冲Fig.5 Dirac pulse

如图4(b)所示,弓入人工边界,截取有限域用有限元法计算,网格尺寸为2 m×2 m.在有限域的人工边界上施加本文提出的时域吸收边界条件.由于本文的时域吸收边界条件对于单层介质是精确的,将人工边界紧靠感兴趣区域放置,取人工边界到孔洞的距离L为一个网格尺寸即2 m.将只用有限元法计算足够大的计算域得出的结果作为参考解.采用Newmark常平均加速度法进行时间积分计算,时步为0.01 s.

(1)方法的精度

考察单层介质时本文提出的吸收边界条件的精度.B点的位移时程如图6 所示(注:n表示无限域的模态数,J表示吸收边界条件的阶数).图6(a)是无限域取满模态时,吸收边界条件的阶数分别为5 和10时的结果.可以看出与参考解相比,J=5 时,1.5 s 之后结果的精度降低,当J提高到10 时,在整个分析时间内其结果和参考解一致.以上分析说明了当吸收边条件的阶数达到要求,吸收边界条件具有很高的精度.

图6 单层介质的吸收边界条件的精度Fig.6 Accuracy of proposed absorbing boundary condition for the single-layered medium

无限域为线弹性介质的前5 阶固有圆频率分别为7.16,21.48,35.92,50.55 和65.31 rad/s,为了减少计算量,将无限域模态缩减,取载荷能激起的前4 阶模态参与计算,其结果与参考解的比较如图6(b)所示.可以看出,缩减模态的结果和参考解重合,说明缩减模态,即只取载荷能激起的土层模态数能够不影响精度的同时减少计算量.图6(b)也展示了同样的边界位置处黏弹性边界的计算结果,可以看出,黏弹性边界精度较低.

(2)方法的稳定性

通过考察式(26)表示的模态空间下吸收边界条件的特征值分布说明吸收边界条件的稳定性,即特征值全部分布在复平面左侧,说明吸收边界条件是稳定的.图7 是吸收边界条件取模态数n=4 和阶数J=10 时特征值的分布.可以看出,特征值全部分布在复平面左侧,且特征值的分布规律是:特征值分布在一族位于复平面左侧的同心半圆上,半圆的数目是模态数n,每个半圆上特征值的数目是J+1.以上分析说明,吸收边界条件对于单层介质是稳定的.

图7 单层介质的吸收边界条件的稳定性Fig.7 Stability of proposed absorbing boundary condition for the single-layered medium

5.2 多层介质

如图8(a)所示,总高H为40 m 的刚性基岩上卧达朗贝尔黏弹性五层介质的场地内嵌一个矩形孔洞,其位于地表以下8 m,尺寸为6 m×8 m.场地每层的高度如图中所示,每层的材料参数如表1 所示.A点施加如图5 所示的狄拉克脉冲点载荷,输出B点的位移时程.

图8 五层数值算例Fig.8 The numerical example of the five-layered media

表1 五层介质的材料参数Table 1 The material parameters of the five-layered media

分析方法与5.1 节单层介质算例类似.因为表示多层介质动力刚度的连分式是通过简单扩展单层介质的连分式获得的,必然会有误差弓进.为了使吸收边界条件对于多层介质具有高精度的特性,人工边界的位置到感兴趣区域的距离L需要被讨论,分别取L=0.5H和L=0.75H.从上到下每层竖向网格的尺寸为2 m,2 m,3 m,3 m 和4 m,横向网格尺寸为2.5 m,如图8(b)所示.时步为0.004 s.

(1)方法的精度

考察多层介质时本文提出的吸收边界条件的精度.无限域为线弹性五层介质的前3 阶固有圆频率分别为17.21,41.39,67.20 rad/s,取载荷能激起的前2 阶模态参与计算,J取10.两种人工边界位置L=0.5H,L=0.75H的结果如图9 所示(注:L表示人工边界到感兴趣区域的距离).可以看出两种边界位置的结果几乎完全重合,说明L=0.5H时所提出的吸收边界条件已具有很高的精度.图9 还展示了L=0.5H边界位置处黏弹性边界的结果,可以看出黏弹性边界精度较低.

图9 五层介质的吸收边界条件的精度Fig.9 Accuracy of proposed absorbing boundary condition for the five-layered media

(2)方法的稳定性

对于上述五层介质,其吸收边界条件取n=2,J=10 时特征值的分布如图10 所示.特征值都分布在复平面左侧,且特征值的分布规律与单层介质相同,即特征值分布在一族位于复平面左侧的同心半圆上,半圆的数目是模态数n,每个半圆上特征值的数目是J+1.以上分析说明,吸收边界条件对于五层介质是稳定的.

图10 五层介质的吸收边界条件的稳定性Fig.10 Stability of proposed absorbing boundary condition for the five-layered media

6 结论

本文建立了一个高精度的时域吸收边界条件,与有限元法结合用于分析达朗贝尔黏弹性成层介质中出平面SH 波动响应.

结合半离散和模态分解获得半无限域模态空间下的频域动力刚度,并显式地表示出了达朗贝尔黏弹性单层介质模态空间下的频域动力刚度.为了将频域动力刚度时域化,采用一种全频域收敛的连分式精确逼近单层介质模态空间下标量形式的频域动力刚度;将标量连分式扩展为矩阵形式用来表示多层介质的频域动力刚度.通过弓进辅助变量,将模态空间下基于连分式的动力刚度关系转化为时域吸收边界条件,进一步转换到物理空间后得到物理空间下的时域吸收边界条件.

数值算例表明,本文建立的吸收边界条件对于达朗贝尔黏弹性单层介质是精确且稳定的;对于达朗贝尔黏弹性多层介质,为了保证精度,需要将人工边界放置在距离感兴趣区域约为0.5 倍无限域高度的位置处.吸收边界条件仅用能被载荷激起的无限域模态数参与计算,不但能减少计算量,同时还不影响精度.

下一步工作可以通过等效线性化方法[31,44]考虑被截掉的无限域的非线性效应.另外,本文的工作还可以被扩展用于模拟矢量波在无限成层介质中的传播问题.