用有限元法等数值方法分析大型结构的瞬态动力响应, 一般的做法是引入人工边界, 将结构-无限地基系统划分为有限域和无限域. 无限域被截去, 有限域用有限元法等数值方法模拟. 为了消除波在人工边界上的虚假反射, 需要在有限域的人工边界上施加人工边界条件(artificial boundary condition, ABC)来描述波动在无限域中的传播^[1]. ABC是基于无限域的控制方程推导得到的. 一般假定无限域是线性的, 因为通常认为无限域中只存在向无穷远传播的波, 不存在向有限域传播的波.

刚性基岩上卧成层土是一种普遍存在的地基形式, 又称其为半开放的波导模型. 与开放的半空间模型中只存在传向无穷远的行波相比, 波导模型由于截止频率的存在, 介质中同时存在截止频率以下的近场快衰波和截止频率以上的传向无穷远的行波^[2-3]. 对于半空间模型, 行波的振幅随着传播距离的增大快速减弱, 所以即使不考虑介质材料阻尼的吸收效应, 其无限域也可以被假定为小应变, 即可以被假定为线性的. 而对于波导模型, 由于能量被困在介质中, 行波的振幅衰减得较为缓慢, 如果仅考虑几何散射不考虑介质吸收效应, 将无限域假定为小应变是有问题的, 所以更合理的做法应该是将波导模型考虑为有材料阻尼的线弹性介质.

本文的研究对象是具有确定尺寸参数的多层波导结构. Wang等^[4]研究表明, 当多层波导内部结构尺寸存在随机缺陷时, 弹性波会产生局部化现象. 该现象使得在原本可以传播的频域范围内弹性波被禁止传播, 进而多层波导结构中会产生不均匀的应力分布. 鉴于问题的复杂性以及篇幅的限制, 本文的研究范围不包括上述多层波导中随机失谐引起的弹性波局部化问题.

如有限域含有非线性等许多情况下, 需要ABC是时域方法. 相较于半空间模型的时域ABC已有大量的研究工作^[5-13], 波导模型的时域ABC研究工作较少. 由于波导模型的ABC需要能够同时模拟快衰波和行波的能量传递, 直接在时域下建立其ABC是有挑战的. 根据是否精确地处理波导模型无限域定解问题, 以下将已有的波导模型的ABC分为精确的ABC和近似的ABC进行阐述.

精确的ABC一般是首先基于解析法或半解析半离散法获得精确的频域动力刚度, 再通过时域化方法建立时域ABC. 解析法一般对单层介质标量波^[14-16]适用, 多层介质的解析解难以获得. 半解析半离散法只对无限域的人工边界进行离散, 对于单层介质和多层介质都适用, 如边界元法(BEM)^[17]、薄层法(TLM)^[18-21]和比例边界有限元法(SBFEM)^{[3, 22-26]}. BEM需要基本解, 而多层介质格林函数的基本解是复杂的, 基本解的存在与否限制了BEM在工程实际中的应用. 相比于BEM, TLM和SBFEM不需要基本解, 能够自动满足无穷远辐射边界条件. 将频域动力刚度变换到时域, 一种方法是时间卷积计算^{[3, 25]}, 该方法不仅耗时, 而且不适用于非线性分析. 学者们通过时间局部化方法建立了精确稳定的时间局部的ABC. Zhao等^[15]基于解析法获得2D单层波导标量波的动力刚度, 再用有理函数逼近方法建立精确的时域ABC. Prempramote等^[22]将解析法获得的动力刚度用双渐近连分式展开, 建立了2D单层标量波的时域ABC. Li等^[16]基于解析法获得2D和3D均匀层标量波的动力刚度, 用动力刚度连分式展开的时间局部化方法建立了精确的时域ABC. Liu等^[26]基于半解析半离散法, 利用算子分离法获得2D单层标量波的时域ABC. Wang等^[23]基于SBFEM和动力刚度的双渐近连分式展开^[22]建立了2D库水的时域ABC, 并将方法嵌入到通用有限元软件ABAQUS中用来计算坝-库水的动力相互作用. 但是上述精确的时域ABC都是针对单层波导标量波获得. 根据作者的调查, 目前几乎没有关于多层波导稳定的精确时域ABC的相关报导. 因为对于多层波导, 其精确的时域ABC易发生失稳, 目前还没有较好地抑制其失稳的方法.

多层波导精确的时域ABC难以获得, 学者们通过对无限域定解问题作适量的简化, 建立了多层波导近似的时域ABC. 近似的低阶时空局部的ABC由于其形式简单, 在工程中被广泛使用, 如黏性边界^[5,27]、多次透射边界^[7,11]以及黏弹性边界^[9-10]. 这类低阶时空局部ABC是基于半空间模型建立的, 用于波导模型时精度较低. Hagstrom等^[28]建立了多层波导标量波的高阶时空局部ABC. 完美匹配层法^[29-30]是近年的研究热点, 它是通过在人工边界外附加一层具有耗能作用的缓冲区来吸收透射波, 但吸收层具有一定厚度, 层厚度和耗能参数选择不当时, 会导致精度低或发生失稳. Zhao等^[32]基于半解析半离散法和一种全频域收敛的动力刚度连分式展开构建了多层波导标量波的时域ABC, 并对其稳定性进行了证明. 吴利华等^[33]将多层波导考虑为达朗贝尔黏弹性介质, 利用和文献[32]类似的思路建立了考虑阻尼的多层波导标量波的时域ABC. 高毅超等^[34]将弹性动力学矢量方程强行解耦为两个标量波方程, 基于SBFEM和双渐近连分式^[22]构建了一种模拟单层波导矢量波的近似ABC. Liu等^[35]将文献[34]的方法嵌入到开源有限元软件OpenSees中, 分析了地下车站的地震动响应. Liu等^[31]基于半解析半离散法和算子分离法建立了多层波导矢量波的近似ABC, 但该方法有时会出现失稳现象. 就作者所知, 目前几乎没有专门针对多层波导矢量波的稳定时域ABC的相关研究.

本文的目标是建立多层波导矢量波稳定的时域ABC. 借鉴文献[34]的做法将多层波导矢量波动方程强行解耦, 得到$x$方向和$y$方向的两个标量波动方程. 将文献[33]构建达朗贝尔黏弹性多层波导标量波的时域ABC的思路用于解耦的$x$方向和$y$方向的两个标量波动方程, 构建含有瑞利阻尼的线弹性多层波导矢量波的稳定的时域ABC.

用有限元法分析如图1所示的刚性基岩上卧成层地基上结构在动力载荷作用下平面内的时程响应. 由于有限元法只能分析有限区域, 需要引入人工边界将左右两个半无限域截去. 为了避免波在人工边界处引起虚假反射, 需要在有限域的人工边界处施加一个人工边界条件(ABC)来模拟波在无限域中的传播.

有限域可以是非线性的, 其有限元方程为

其中, 下标b和i分别表示人工边界部分和有限域除了人工边界的其他部分, $ M$是质量矩阵, $ C$是阻尼矩阵, $ K$是刚度矩阵, $ u$是平面内位移向量, $ F_{\rm i} ( u_{\rm i} ,{\dot{{ u}}}_{\rm i} )$是非线性内力, $ b$是体力, ${ f}_{\rm i}$是施加在有限域的载荷, ${ f}_{\rm b} $是施加在有限域人工边界上的力, 即本文要建立的ABC, 用来表示被截掉的无限域和有限域的相互作用.

一般认为无限域是线性的, 若无限域包含非线性, 需要对其进行等效线性化处理, 使其满足线性条件. 无限域的波动方程为

其中, 变量上方的点表示对时间求导, ($x$, $y)$ 表示笛卡尔坐标, 上标T表示转置; ${U}=\{u_x \ \ u_y\}^{\rm T}$, $u_x $, $u_y $分别为$x$方向和$y$方向的位移, $u_x$和$u_y $在不同土层的分界面以及人工边界处是连续的; $\rho $是质量密度, $\lambda $是拉梅常数, $G$是剪切模量, 每一层的材料参数是常量, 不同层的材料参数可不同. 另外, 无限域的上表面物理边界条件是应力为0, 下表面物理边界条件是位移为0. 土层的初始状态为静止.

将无限域作为研究对象, 建立时域吸收边界条件. 文献[33]建立了模拟刚性基岩上卧多层介质中标量波传播的ABC, 本文将该方法扩展到矢量波. 方法的思路是: (1) 将矢量波动方程简化为$x$方向和$y$方向解耦的两个标量波动方程; (2) 基于比例边界有限元法, 得到两个标量波动方程模态空间下半离散的频域动力刚度方程; (3) 通过文献[33]中提出的连分式来表示频域动力刚度; (4) 引入辅助变量, 将连分式时域化, 得到时域吸收边界条件; (5) 所建立的ABC可以和有限域的有限元方程无缝耦合, 通过数值积分算法求解.

高毅超等^[34]通过将矢量波动方程强行解耦, 基于双渐近连分式建立了单层波导模型的ABC. 本文借鉴其做法, 忽略无限域波动方程(2)中$x$和$y$的耦合项, 简化后得到关于$x$方向和$y$方向的两个标量方程分别为

其中, $x$和$y$解耦简化处理后的应力-应变关系为

$\begin{eqnarray*} &&\sigma _x =(\lambda +2G)\varepsilon _x =(\lambda+2G)\frac{\partial u_x }{\partial x} \\&&\sigma _y =(\lambda +2G)\varepsilon _y =(\lambda +2G)\frac{\partial u_y }{\partial y} \\&&\tau _{xy} =G\gamma _{xy} =G\frac{\partial u_y }{\partial x} \\&&\tau _{yx} =G\gamma _{yx} =G\frac{\partial u_x }{\partial y} \end{eqnarray*}$

经过解耦简化处理后的上下表面物理边界条件分别为

其中, $H$为半无限域的总层高.

根据比例边界有限元法(SBFEM)^[23]建立模态空间下的频域动力刚度方程.

沿着土层$y$方向半离散简化处理后的无限域波动方程(3), 同时考虑简化后的物理边界条件式(4)及介质交界面的连续条件, 得到无限域简化弹性波动方程的比例边界有限元(SBFE)方程

基于简化后的本构关系得到无限域的人工边界上$x$和$y$方向的等效结点力分别为

其中, ${ u}_{\rm b}^x $, ${ u}_{\rm b}^y $分别为结点位移向量; $ N$是线性形函数; $m=1$ 和$m=-1$分别表示被截掉的左、右无限域的外法线方向; 以及

$\begin{eqnarray*} &&{ E}_1^x =\int_{\rm 0}^H {{ N}^{\rm T} (\lambda +2} G){ N}{\rm d}y, \ \ { E}_{\rm 1}^y =\int_{\rm 0}^H {{ N}^{\rm T} G} \, { N}{\rm d}y \\&&{ E}_{\rm 2}^x =\int_{\rm 0}^H {\frac{\partial \, { N}^{\rm T}}{\partial y}} G\frac{\partial \, { N}}{\partial y}{\rm d}y, \ \ { E}_{\rm 2}^y =\int_{\rm 0}^H {\frac{\partial \, { N}^{\rm T}}{\partial y}} (\lambda +2G)\frac{\partial \, { N}}{\partial y}{\rm d}y \\&&{ E}_{\rm 3}^x =\int_{\rm 0}^H {{ N}^{\rm T}} \rho \, { N}{\rm d}y, \ \ \, { E}_{\rm 3}^y ={ E}_{\rm 3}^x \end{eqnarray*}$

引言的第二段阐述了波导模型应该考虑阻尼的吸收效应. 由于瑞利阻尼是目前应用最广泛的阻尼矩阵的数值模型, 本文假定无限域土层含有瑞利阻尼. 瑞利阻尼定义阻尼矩阵是质量阵和刚度阵的线性组合^[36], 即${C}_{\rm Rayleigh} =a_0 \, { M}+a_1 \, { K}$, 其中$\left\{{{\begin{array}{c} {a_0 } \\ {a_1 } \\ \end{array} }} \right\}=\dfrac{2\xi }{\omega _{i} +\omega _j }\left\{{{\begin{array}{c} {\omega _{i} \omega _j } \\ 1 \\ \end{array} }} \right\}$, $\xi $为阵型阻尼比, $\omega _{i} $, $\omega _j$为选择用于确定阻尼系数$a_{0}$和$a_{1}$的频率点. 瑞利阻尼的引入体现在式(5a)和式(5b)的第四项, 需要说明的是, 这里无限域瑞利阻尼的质量阵和刚度阵都是通过半离散$x$和$y$解耦的标量控制方程得到的, 注意其与二维有限元法中的瑞利阻尼相区别.

为了使问题简单化, 利用模态分解法^[36], 将物理空间下的SBFE方程变换到模态空间. 计算如下广义特征值问题

$\begin{eqnarray*} &&{ E}_2^x {\varPhi}_n^x ={ E}_3^x {\varPhi }_n^x [{\varLambda}_n^2 ]^x, \ \ \, { E}_2^y {\varPhi }_n^y ={ E}_3^y {\varPhi }_n^y [{\varLambda}_n^2 ]^y \\&&{\varPhi }_n^{\,x{\rm T}}{ E}_3^x {\varPhi }_n^x = I, \ \ {\varPhi }_n^{\,y{\rm T}}{ E}_3^y {\varPhi }_n^y = I \\&&{\varPhi }_n^{\,x{\rm T}}{ E}_{\rm 2}^x {\varPhi }_n^x =[{\varLambda }_n^2 ]^x, \ \ {\varPhi }_n^{\,y{\rm T}}{ E}_{\rm 2}^y {\varPhi }_n^y =[{\varLambda }_n^2 ]^y \end{eqnarray*}$

其中, ${\varPhi }_n^x $, ${\varPhi }_n^y$分别表示$x$和$y$方向包含前$n$阶特征向量的矩阵; $[{\varLambda }_n^2]^x$, $[{\varLambda }_n^2]^y$分别表示$x$和$y$方向由前$n$阶特征值组成的对角阵, 其对角线上的元素为$\omega _n^2$; $ I$表示单位阵.

人工边界上$x$方向和$y$方向的结点力、位移在模态空间下和物理空间下的关系为

其中, ${ q}_{\rm b}^x $, ${ q}_{\rm b}^y $和${ p}_{\rm b}^x $, ${ p}_{\rm b}^y $分别表示$x$和$y$方向模态空间下的位移向量和力向量.

将式(7)代入式(5)和式(6), 得到$x$和$y$方向上模态空间下的SBFE方程

其中, ${ e}_{\rm 1}^x ={\varPhi}_n^{\,x{\rm T}} \, { E}_1^x {\varPhi}_n^x $, ${ e}_{\rm 1}^y ={\varPhi }_n^{\,y{\rm T}} \, { E}_1^y {\varPhi }_n^y $.

假定模态空间下无限域人工边界处$x$方向和$y$方向的结点力-位移关系为

其中, 定义${ s}_x^\infty $和${ s}_y^\infty $分别为模态空间下半无限域的$x$和$y$方向的频域动力刚度; 上标$\sim $表示频域下的值.

结合式(8)$\sim\!$式(10), 得到$x$和$y$方向模态空间下的动力刚度方程

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

为了将动力刚度时域化, 用连分式展开来表示式(11)中的${ s}_x^\infty $和${ s}_y^\infty $.

半无限域是单层介质时, $e_{\rm 1}^x =c_{\rm p}^2 \, { I}$, $e_{\rm 1}^y =c_{\rm s}^2{ I}$, 其中, $c_{\rm p} =\sqrt {(\lambda +2G)/\rho } $, $c_{\rm s} =\sqrt {G/\rho }$. 此时, 动力刚度${ s}_x^\infty $和${ s}_y^\infty $矩阵是对角形式的, 其对角线上的元素统写为

其中, $z={\rm i}\omega /\omega _n $, $\beta =a_0 /\omega _n +a_1 \omega _n $, i是虚数单位, $\omega $表示圆频率. $x$和$y$方向分别有$c=c_{\rm p} $和$c=c_{\rm s} $.

式(12)中的$\sqrt {1+z^2+\beta z} $可以被表示为若干种递归形式, 作者在文献[32,33]中提出了一种如下形式的递归式, 并且证明了其在频率域$[0, \ \infty ]$是指数收敛的. 所以后文基于该连分式建立的时域ABC也应该是稳定的.

由于$x$方向和$y$方向的SBFE方程形式相同, 以下仅给出$x$方向的推导过程, $y$方向的推导同$x$方向. 以下推导思路类似文献[33], 本文仅给出与文献[33]不同的地方以及相关重要公式, 详细推导过程可参见文献[33]. 根据式(12)和式(13), 无限域$x$方向模态空间的动力刚度可表示为如下的连分式

其中

当$j$为奇数时, 有

当$j$为偶数时, 有

和文献[33]一样, 连分式(14)也可用来近似地表示多层介质的动力刚度. 对于多层介质, 连分式中的${ g}_0^x $ 和${ h}_0^x$通过求解如下黎卡提方程得到

式中, ${ g}_j^x$, ${ h}_j^x $的关系式同单层介质.

同文献[33]一样, 通过引入辅助变量的方式将连分式时域化. 将式(14)代入式(10a), 同时引入辅助变量${ s}_j {\tilde{{ q}}}_j^x ={\tilde{{ q}}}_{j-1}^x$, 得到

其中, ${\tilde{{ q}}}_{J+1}^x ={ 0}$, ${\tilde{{ q}}}_0^x ={\tilde{{ q}}}_{\rm b}^x $.

将式(18)傅里叶变换到时域, 得到$x$方向模态空间下的时域ABC

其中

$ \begin{array}{ q}_{\rm c}^x =\left\{ {{\begin{array}{c@{\quad }c@{\quad }c@{\quad }c@{\quad }c} {{ q}_1^{\,x{\rm T}} } & {{ q}_2^{\,x{\rm T}} \;} & \ldots & {{ q}_{J-1}^{\,x{\rm T}} } & {{ q}_J^{\,x{\rm T}} } \\ \end{array} }} \right\}^{\rm T} \\ { c}_{\rm bc}^{x\infty} =\left[ {{\begin{array}{c@{\quad }c@{\quad }c@{\quad }c@{\quad }c} {-{ I}} & { 0} & \cdots & { 0} & { 0} \\ \end{array} }} \right] \\{ k}_{\rm cb}^{x\infty} =\left[ {{\begin{array}{c@{\quad }c@{\quad }c@{\quad }c@{\quad }c} {-{ I}} & { 0} & \cdots & { 0} & { 0} \\ \end{array} }} \right]^{\rm T} \\{ c}_{\rm cc}^{x\infty} =\left[ {{\begin{array}{c@{\quad }c@{\quad }c@{\quad }c@{\quad }c} {{ h}_1^x } & {-{ I}} & \cdots & { 0} & { 0} \\ { 0} & {{ h}_2^x } & \cdots & { 0} & { 0} \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ { 0} & { 0} & { 0} & {{ h}_{J-1}^x } & {-{ I}} \\ { 0} & { 0} & { 0} & { 0} & {{ h}_J^x } \\ \end{array} }} \right] \\ { k}_{\rm cc}^{x\infty} =\left[ {{\begin{array}{c@{\quad }c@{\quad }c@{\quad }c@{\quad }c} {{ g}_1^x } & { 0} & \cdots & { 0} & { 0} \\ {-{ I}} & {{ g}_2^x } & \cdots & { 0} & { 0} \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ { 0} & { 0} & {-{ I}} & {{ g}_{J-1}^x } & { 0} \\ { 0} & { 0} & { 0} & {-{ I}} & {{ g}_J^x } \\ \end{array} }} \right] \end{array} $

同理可求出$y$方向模态空间下的时域ABC

将式(7b)代入式(19)和式(20), 得到物理空间下的ABC

式(20)中, 人工边界上单个结点的水平自由度和竖向自由度是空间解耦的, 但是人工边界结点间的水平自由度和结点间的竖向自由度都是空间耦联的. 其中

$\begin{eqnarray*} \begin{array}{l} { u}_{\rm b} =\left\{ {{\begin{array}{*{20}c} {{ u}_{\rm b}^{\,x{\rm T}}} & {{ u}_{\rm b}^{\,y{\rm T}}} \\ \end{array} }} \right\}^{\rm T},\ { q}_{\rm c} =\left\{ {{\begin{array}{*{20}c} {{ q}_{\rm c}^{\,x{\rm T}}} & {{ q}_{\rm c}^{\,y{\rm T}}} \\ \end{array} }} \right\}^{\rm T},\ { f}_{\rm b} =\left\{ {{\begin{array}{*{20}c} {{ f}_{\rm b}^{\,x{\rm T}}} & {{ f}_{\rm b}^{\,y{\rm T}}} \\ \end{array} }} \right\}^{\rm T} \\ { C}_{\rm bb}^\infty =\left[ {{\begin{array}{*{20}c} {{ C}_{\rm bb}^{x\infty} } & { 0} \\ { 0} & {{ C}_{\rm bb}^{y\infty} } \\ \end{array} }} \right], \ \ { K}_{\rm bb}^\infty =\left[ {{\begin{array}{*{20}c} {{ K}_{\rm bb}^{x\infty} } & { 0} \\ { 0} & {{ K}_{\rm bb}^{y\infty} } \\ \end{array} }} \right] \\ { C}_{\rm bc}^\infty =\left[ {{\begin{array}{*{20}c} {{ C}_{\rm bc}^{x\infty} } & { 0} \\ { 0} & {{ C}_{\rm bc}^{y\infty} } \\ \end{array} }} \right], \ \ { K}_{\rm cb}^\infty =\left[ {{\begin{array}{*{20}c} {{ K}_{\rm cb}^{x\infty} } & { 0} \\ { 0} & {{ K}_{\rm cb}^{y\infty} } \\ \end{array} }} \right] \\ { C}_{\rm cc}^\infty =\left[ {{\begin{array}{*{20}c} {{ C}_{\rm cc}^{x\infty} } & { 0} \\ { 0} & {{ C}_{\rm cc}^{y\infty} } \\ \end{array} }} \right], \ \ { K}_{\rm cc}^\infty =\left[ {{\begin{array}{*{20}c} {{ K}_{\rm cc}^{x\infty} } & { 0} \\ { 0} & {{ K}_{\rm cc}^{y\infty} } \\ \end{array} }} \right] \\ { C}_{\rm bb}^{x\infty} ={ E}_3^x {\varPhi }_n^x \, { h}_0^x {\varPhi }_n^{\,x{\rm T}}{ E}_3^x, \ \ \ \ { C}_{\rm bb}^{y\infty} ={ E}_3^y {\varPhi }_n^y \, { h}_0^y {\varPhi }_n^{\,y{\rm T}}{ E}_3^y \\ { K}_{\rm bb}^{x\infty} ={ E}_3^x {\varPhi }_n^x \, { g}_0^x {\varPhi }_n^{\,x{\rm T}}{ E}_3^x, \ \ \ \ { K}_{\rm bb}^{y\infty} ={ E}_3^y {\varPhi }_n^y \, { g}_0^y {\varPhi }_n^{\,y{\rm T}}{ E}_3^y \\ { C}_{\rm bc}^{x\infty} =\left[ {{\begin{array}{*{20}c} {-{ E}_3^x {\varPhi }_n^x }\ \ \ & { 0} & \cdots & { 0} & { 0} \\ \end{array} }} \right] \\ { C}_{\rm bc}^{y\infty} =\left[ {{\begin{array}{*{20}c} {-{ E}_3^y {\varPhi }_n^y }\ \ \ & { 0} & \cdots & { 0} & { 0} \\ \end{array} }} \right] \\ { K}_{\rm cb}^{x\infty} =\left[ {{\begin{array}{*{20}c} {-{ E}_3^{\,x{\rm T}} {\varPhi }_n^x }\ \ \ & { 0} & \cdots & { 0} & { 0} \\ \end{array} }} \right]^{\rm T} \\ { K}_{\rm cb}^{y\infty} =\left[ {{\begin{array}{*{20}c} {-{ E}_3^{\,y{\rm T}} {\varPhi }_n^y }\ \ \ & { 0} & \cdots & { 0} & { 0} \\ \end{array} }} \right]^{\rm T} \\ { C}_{\rm cc}^{x\infty} ={ c}_{\rm cc}^{x\infty}, \ \ { C}_{\rm cc}^{y\infty} ={ c}_{\rm cc}^{y\infty}, \ \ { K}_{\rm cc}^{x\infty} ={ k}_{\rm cc}^{x\infty}, \ \ { K}_{\rm cc}^{y\infty} ={ k}_{\rm cc}^{y\infty} \end{array} \end{eqnarray*}$

本文建立的ABC可与有限元法无缝耦合. 将式(20)代入有限元方程式(1), 得

式(21)可通过标准的时间积分算法如Newmark常平均加速度法等求解.

由于在有限域的人工边界上施加ABC增加的额外自由度为

其中, $n_{\rm L}$和$n_{\rm R}$分别表示左无限域和右无限域参与计算的模态数. $J$表示连分式的阶数, 通过算例分析统计, 一般可取$J=3$, $J$的取值会在第3节详细阐述. 对于实际地震工程等问题, 根据经验一般取模态数$n\leqslant 3$. 所以额外增加的自由度数一般不会超过36, 与有限域的自由度数相比, 额外增加的自由度数几乎可以忽略不计.

本文建立的时域ABC是对无限域矢量波动的近似模拟, 为了实现高精度, 需要将有限域的人工边界放置在距离兴趣域$L$处, $L$的取值将在第3节通过数值算例详细分析.

本文的ABC是基于全频域收敛稳定的连分式建立的, 所以理论上ABC应该是时域稳定的. 方法的稳定性将在第3节通过数值算例进一步说明.

本节通过数值算例说明方法的精度和稳定性. 本文方法包含3个可变参数: 无限域的模态数$n$, 连分式阶数$J$和人工边界远离兴趣域的距离$L$. 文献[33]表明, 仅用能被载荷激起的无限域的模态数$n$参与计算, 不但能减少计算量, 同时还不影响精度. 下文详细分析$J$和$L$的取值要求.

分析如图2所示的3种模型. Model 1和Model 2都是自由场地, Model 1的总层高$H$为20 m, Model 2的总层高$H$为40 m. 比较Model 1和 Model 2的结果可以考察$L$, $J$取值与$H$的关系. Model 3内含一个尺寸为5$\times $5的地下方洞, 其总层高$H$为20 m. 比较Model 1和Model 3的结果可以考察有无地下结构及地下结构尺寸对$L$, $J$取值的影响. 3种模型都是在地表作用10 m长的水平线载荷, 线载荷为狄拉克脉冲, 其时程及频谱如图3所示. 为了考察狄拉克脉冲宽度是否会影响本文方法的精度和收敛性, 考虑三种狄拉克脉冲: $T=0.4$, 0.2和0.1, 其相应的频谱范围分别为0$\sim$10 Hz, 0$\sim$20 Hz和0$\sim$40 Hz. 图2中每个模型的兴趣域(region of interest, ROI)为蓝虚线框起来的部分, 红实线表示人工边界, 两条红实线框起来的部分是用于计算的有限域. $L$表示人工边界到兴趣域的距离.

为了分析土层材料参数对$L$, $J$取值的影响, 将图2中3种模型分别都考虑为均匀介质和分别都考虑为4层介质. 分别将单层Model 1, 单层Model 2和单层Model 3记为Case-1, Case-2和Case-3; 分别将4层Model 1, 4层Model 2 和4层Model 3记为Case-4, Case-5和Case-6. Case-1$\sim$Case-6土层都被考虑为有瑞利阻尼的线弹性介质, 其土层材料参数见表1. 6种工况的计算时长都为2 s, 输出A点的水平位移时程和水平加速度时程.

用有限元法分析Case-1$\sim$Case-6. 图2所示的Model 1$\sim$Model 3中两条红实线表示的人工边界框起来的部分为有限域, 人工边界到兴趣域的距离被标记为$L$, 在人工边界处施加ABC. 有限域用四边形网格离散, Newmark常平均加速度法用于时间积分计算. 狄拉克脉冲中$T=0.4$时, 有限域的网格尺寸为2.5 m$\times $2.5 m, 时间步长为0.005 s. $T=0.2$和0.1时的网格尺寸和时间步长分别是$T=0.4$时的0.5倍和0.25倍.

用于确定瑞利阻尼系数的频率点的选择方式有较多研究^[37]. 本节的重点是验证本文提出方法的精度和稳定性, 由于瑞利阻尼系数的具体数值不影响本文方法的性能, 所以这里采用简单的频率点的选取方式, 取土层前2阶频率用于确定其瑞利阻尼系数$a_{0}$和$a_{1}$. Case-1$\sim$Case-6土层的阻尼比$\xi $都取为0.05, 表1列出了6种工况的瑞利阻尼系数值. 狄拉克脉冲中$T=0.4$时, Case-1$\sim$Case-6用于计算的无限域模态数$n$分别为2, 3, 2, 1, 3和1.

将足够大有限域的有限元结果作为验证本文方法的参考解, 大有限域尺寸的选取原则是人工边界上的虚假反射波在计算时间内不会污染兴趣域.

用数值算例结果分析$J$和$L$的取值对方法精度的影响. 定义如下时程相对误差

其中, $X_{j}$为基于有限元法和ABC得到的小有限域的时程结果, $Y_{j}$为仅基于有限元法得到的大有限域的时程结果(即参考解), $N$为总时步.

图4展示了Case-1$\sim$Case-6在连分式阶数$J$分别为1, 2, 3, 4, 5的情况下, $L=H \sim 5H$时对应的$A$点水平加速度时程的相对误差, 这里狄拉克脉冲中$T=0.4$. 为了说明本文方法的优越性, 图4也展示了同样$L$下黏性边界^[5] (VB)、刘晶波等^[9]提出的黏弹性边界(VSB-Liu)和杜修力等^[10]提出的黏弹性边界(VSB-Du)的结果. 黏性边界和黏弹性边界都是经典的ABCs, 因其形式简单且时域稳定, 在工程中被广泛使用. 黏性边界^[5]是基于平面波假定建立的阻尼器边界条件, 其每层的阻尼元件法向参数为$C_N =\rho c_{\rm p}$, 切向参数为$C_{\rm T}=\rho c_{\rm s}$. 黏弹性边界是基于柱面波假定建立的并联弹簧-阻尼器边界条件. 刘晶波等^[9]提出的黏弹性边界每层的弹簧元件法向参数为$K_{\rm N}=2G/r$, 切向参数为$K_{\rm T}=3G/(2r)$, 每层的阻尼元件法向参数为$C_{\rm N} =\rho c_{\rm p}$, 切向参数为$C_{\rm T} =\rho c_{\rm s}$; 杜修力等^[10]提出的黏弹性边界每层的弹簧元件法向参数为$K_{\rm N}=(\lambda +2G)/(3.6r)$, 切向参数为$K_{\rm T}=G/(3.6r)$, 每层的阻尼元件法向参数为$C_{\rm N} =1.1\rho c_{\rm p} $, 切向参数为$C_{\rm T} =1.1\rho c_{\rm s} $. 3个公式中每层的材料参数按实际值取, $r$认为是有限域的竖向对称轴到人工边界的水平距离. 从图中可以看出, 6种工况呈现的结果规律类似. 每种工况下, 黏性边界和2种黏弹性边界结果的相对误差接近, 都远大于本文方法的结果. 本文方法的$J=1$, 2, 3, 4, 5对应的5条曲线都是$J=1$时相对误差较大, $J=2$时相对误差快速变小, $J=3$, 4和5时三者的结果基本重合. 根据大量算例结果(不仅限于图4所列的算例结果)统计, 当$J>3$时, 再增加$J$的值, 基本不再改善精度. 这是因为本文方法是一种近似方法, 连分式只是近似地表示无限域的动力刚度, 即使取很高的阶数, 连分式也不能逼近无限域动力刚度的精确解. 因近似损失的精度需要通过增大人工边界到兴趣域的距离$L$来弥补. 为了减少可变因素, 作者建议在使用本文方法时, 参数$J$可以统一取为3.

进一步详细分析$L$的选取原则. 分别计算Case-1$\sim$Case-6, 当$J=3$, $L$为$H \sim 5H$时$A$点水平位移时程和水平加速度时程各自的相对误差, 这里狄拉克脉冲中$T=0.4$. 将相对误差$\leqslant 5 $时各自对应的$L$值, 以及同样$L$下黏性边界和黏弹性边界结果的相对误差列于表2. 根据表2可以得出以下结论 (表中第2列$u$表示位移, $a$表示加速度; 第3列表头$L$表示人工边界到兴趣域的距离; 第4列表头ABC表示本文方法的结果; 第5列表头VB表示黏性边界^[5]的结果; 第6列表头VSB-Liu表示刘晶波等^[9]提出的黏弹性边界的结果; 第7列表头VSB-Du表示杜修力等^[10]提出的黏弹性边界的结果).

(1)对比表2的最后3列可以看出, 黏性边界和黏弹性边界结果的相对误差接近, 都约为本文方法结果的5$\sim$6倍左右, 除了Case-2和Case-5的加速度时程结果, 其黏性边界和黏弹性边界结果的相对误差是本文方法结果的2$\sim$3倍. 另外, 本文方法与黏性边界和黏弹性边界相比, 仅多出少量的辅助自由度数, 对于Case-1$\sim$Case-6, 其值分别为24, 36, 24, 12, 36, 12. 因此本文方法相比于黏性边界和黏弹性边界, 能够在几乎不增加计算量的情况下使精度提高约2$\sim$6倍. 精度的提高程度与土层的总厚度有关, 原因是本文方法是针对层模型建立的, 其应用在深厚土层和浅土层都具有高精度; 黏性边界和黏弹性边界是针对半空间模型建立的, 由于深厚土层相较于浅土层更接近半空间模型, 因此其在深厚土层中要比在浅土层中的精度高.

(2)对比表2第3列表示的水平位移对应的$L$值和水平加速度对应的$L$值可以看出, 每种工况下, 水平加速度对应的$L$值都是水平位移对应结果的2倍, 除了Case-5是2.5倍. 说明人工边界的位置以加速度为控制指标来确定更为严格.

(3)比较单层的Case-1和Case-3水平加速度对应的$L$值. Case-1和Case-3总层高相同, Case-1为自由场地, Case-3含有5 m$\times $5 m的地下结构. 两者都是$L=3H$时, 水平加速度时程的相对误差能控制在5%以内. 说明$L$的取值基本不受是否含有地下结构或结构尺寸影响. 比较4层的Case-4和Case-6能得出同样的结论.

(4)比较单层Case-1和Case-2水平加速度对应的$L$值. Case-1和Case-2都是自由场地, Case-2的总层高$H$是Case-1的2倍. 两者都是$L=3H$时, 水平加速度时程的相对误差能控制在5 以内. 同样, 比较4层的 Case-4和Case-5. 为了将水平加速度时程的相对误差控制在5 以内, Case-4要求$L=2H$, Case-5要求$L=2.5H$. 说明$L$的取值大约与土层总层高$H$成正比关系, 关系系数与土层的材料参数有关. 根据大量算例统计, 一般$L=3H$都能将加速度时程相对误差控制在5 以内.

图5展示了本文方法的计算结果满足工程精度要求 (相对误差$\leqslant 5 )$ 时的$A$点水平加速度时程. 其中, $J=3$, Case-1$\sim$Case-6对应的$L$值分别为3$H$, 3$H$, 3$H$, 2$H$, 2.5$H$, 2$H$. 同样$L$下刘晶波等提出的黏弹性边界(VSB-Liu)的结果也被展示在图5中. 从图中可以看出, 本文方法的结果和参考解几乎完全吻合, 而黏弹性边界基本是在0.5 s之后结果的精度较低.

为了研究图3中狄拉克脉冲的宽度对本文方法的计算精度和收敛性的影响, 图6展示了Case-6在$T=0.2$和$T=0.1$两种载荷下, 本文方法($J=1\sim 5$, $L=H\sim 5H)$和VSB-Liu方法的计算结果, 图中纵坐标为式(23)所示的$A$点水平加速度时程的相对误差. $T=0.2$时, 本文方法中模态数$n$取为5; $T=0.1$时, $n$取为10. 可以看出, 图6(a) ($T=0.2$)和图6(b) ($T=0.1$)中结果的规律和图4 ($T=0.4$)中Case-6结果的规律相似, 都是黏弹性边界结果远大于本文方法的结果, 本文方法$J=3$时结果已收敛, 再增大$J$值, 基本不再改善精度; 因近似损失的精度需要通过增大人工边界到兴趣域的距离$L$来弥补, 都是当$L=2H$时, 能满足工程精度要求 (相对误差<5 ), 随着$L$的增大, 结果向0收敛. 因此, 图 3 中狄拉克脉冲的宽度基本不会影响本文方法的计算精度和收敛性.

由于本文的ABC是基于全频域收敛的连分式建立的, 理论上ABC应该是时域稳定的. 通过后验的方法进一步验证本文方法的稳定性. 根据线性系统理论, 若ABC系统的特征值全部分布在复平面左侧, 则说明ABC稳定. 图7展示了Case-1$\sim$Case-6的水平方向上和竖直方向上式(19)所示的模态空间下的ABC系统的特征值在复平面的分布. 总特征值的个数为$n\times (J+1)$. 从图中可以看出ABC特征值的实部都小于0, 且其分布规律为: 特征值分布在一系列半圆上, 且半圆的数目是模态数$n$, 每个半圆上特征值的数目是$J+1$. 由于本算例每一个虚部为0的特征值都是重根, 因此图中呈现出来的是每个半圆上特征值的数目为$J$.

本文将文献[33]中针对标量波提出的方法扩展到矢量波, 建立了一种稳定的近似时域人工边界条件(ABC)来模拟含有瑞利阻尼的线弹性多层波导模型的平面内矢量波动. 相较于文献[33]的方法只能模拟标量波问题, 本文提出的方法不仅适用于矢量波问题, 对标量波问题也同样适用.

提出的方法中影响计算精度和计算效率的参数有3个, 分别为无限域的模态数$n$, 连分式阶数$J$, 和人工边界远离兴趣域的距离$L$. 数值算例表明, 仅用能被载荷激起的无限域的模态数$n$参与计算即可. 由于本文方法是一种近似方法, 一般当$J>3$时, 再增大$J$值, 基本不再改善精度. 作者建议使用本文方法时参数$J$可以统一取为3. 因近似处理损失的精度需要通过增大人工边界到兴趣域的距离$L$来弥补. 确定人工边界位置$L$以加速度为分析指标更为严格. $L$的取值基本与地下结构尺寸无关, 它与土层的总层高$H$约成正比关系, 关系系数与土层的材料参数有关. 一般取$L=3H$都能将加速度时程相对误差控制在5%以内. 本文数值算例结果表明在同样的$L$下, 与工程中被广泛使用的黏性边界和黏弹性边界相比, 本文方法能够在几乎不增加计算量的情况下一般大约能使精度提高2$\sim$6倍, 精度的提高程度与土层的总厚度有关.

本文提出的ABC在理论上是稳定的, 数值算例也进一步验证了其稳定性. 算例表明模态空间下ABC系统的特征值呈一定规律分布在复平面的左侧, 其分布规律为: 特征值分布在一系列半圆上, 半圆的数目是模态数$n$, 每个半圆上特征值的数目是$J+1$.

本文提出的ABC虽然人工边界上单个结点的水平自由度和竖向自由度是空间解耦的, 但是人工边界结点间的水平自由度和结点间的竖向自由度都是空间耦联的. 下一步工作可以考虑将方法进一步优化, 使人工边界上结点间的水平自由度和结点间的竖向自由度都空间解耦, 更方便工程应用. 另外, 可以考虑借鉴本文方法的研究思路建立矢量弹性波精确的时域人工边界条件.