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

图1刚性基岩上卧达朗贝尔黏弹性成层介质

Fig. 1D'Alembert viscoelastic multilayered media overlying rigid bedrock

(1)$ \begin{eqnarray} \label{eq1} G\frac{\partial ^2u}{\partial x^2}+G\frac{\partial ^2u}{\partial y^2}=\rho \ddot{{u}}+\eta \dot{{u}} \end{eqnarray}$

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

(2)$ \begin{eqnarray} \label{eq2} \tau _y =0, \ \ u=0 \end{eqnarray}$

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

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

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

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

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

(3)$ \begin{eqnarray} \label{eq3} E_1 \frac{\partial ^2{u}_{\rm b} }{\partial x^2}-E_2 {u}_{\rm b} -E_3 {\ddot{{u}}}_{\rm b} -E_4 {\dot{{u}}}_{\rm b} ={\bf 0} \end{eqnarray}$

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

(4)$ \begin{eqnarray} \label{eq4} {f}_{\rm b} =\int_0^H {m{N}^{\rm T}\tau _x} {\rm d}y=m{E}_1 \frac{\partial {u}_{\rm b} }{\partial x} \end{eqnarray}$

其中, ${E}_1 =\int_0^H {{N}^{\rm T}G{N}{\rm d}y} $, ${E}_2=\int_0^H {\frac{\partial {N}^{\rm T}}{\partial y}G\frac{\partial {N}}{\partial y}{\rm d}y} $, ${E}_3 =\int_0^H {{N}^{\rm T}\rho {N}{\rm d}y}$和${E}_4 =\int_0^H {{N}^{\rm T}\eta {N}{\rm d}y} $, $u_{\rm b}$ 表示人工边界上的结点位移向量, $H$表示无限域的总高度, $N$ 是线性形函数, T表示矩阵转置, $m=1$ 和$m=-1$分别表示被截掉的左、右无限域的外法线方向. 值得注意的是, 对于每个单元, 式(3)中的阻尼矩阵$E_{4}$是质量矩阵$E_{3}$的$\alpha$倍. 所以达朗贝尔黏弹性介质相当于与质量成正比的瑞利阻尼形式^[32,43].

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

(5)$ \begin{eqnarray} \label{eq5} {\tilde{{f}}}_{\rm b} ={S}^\infty {\tilde{{u}}}_{\rm b} \end{eqnarray}$

顶标$\sim $表示频域下的值, ${S}^\infty$表示物理空间下半无限域的频域动力刚度.

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

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

(6)$ \begin{eqnarray} \label{eq6} {E}_2 { \varPhi }={E}_3 { \varPhi \varLambda }^2 \end{eqnarray}$

其中, ${ \varPhi }$是特征向量矩阵, ${ \varLambda }^2$是对角形式的特征值矩阵, 其对角线上的元素为$\omega _n^2 $. 考虑半无限域模态缩减, 有: ${ \varPhi }_n^{\rm T} {E}_3 { \varPhi }_n ={I}$和${ \varPhi }_n^{\rm T} {E}_2 { \varPhi}_n ={ \varLambda }_n^2 $. ${ \varPhi }_n$表示包含前$n$阶特征向量的矩阵, ${ \varLambda }_n^2$表示由前$n$阶特征值组成的对角阵, $I$表示单位阵.

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

(7)$ \begin{eqnarray} {u}_{\rm b} ={ \varPhi }_n {q} \end{eqnarray}$

(8)$ \begin{eqnarray} {f}_{\rm b} ={E}_3 { \varPhi }_n {p} \end{eqnarray}$

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

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

(9)$ \begin{eqnarray} \label{eq9} {p}=m{e}_1 \frac{\partial {q}}{\partial x} \end{eqnarray}$

(10)$ \begin{eqnarray} \label{eq10} {p}=m{e}_1 \frac{\partial {q}}{\partial x} \end{eqnarray}$

其中${e}_1 ={ \varPhi }_n^{\rm T} {E}_1 { \varPhi }_n$, ${e}_4 ={ \varPhi }_n^{\rm T} {E}_4 { \varPhi }_n$.

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

(11)$ \begin{eqnarray} \label{eq11} {\tilde{{p}}}={S\tilde{{q}}} \end{eqnarray}$

其中${S}={ \varPhi }_n^{\rm T} {S}^\infty { \varPhi }_n$为模态空间下半无限域的频域动力刚度.

当半无限域是单层介质时, ${E}_1 =c^2{E}_3 $, $c=\sqrt {G/\rho } $; ${E}_4 =\alpha {E}_3 $; 所以$e_1 =c^2{I}$, ${e}_4 =\alpha I$. 合式(9) $\sim\!$式(11),

得到对角形式的动力刚度矩阵$S$对角线上的元素为

(12)$ \begin{eqnarray} \label{eq12} S=c\sqrt{\omega_n^2 -\omega ^2+\alpha {\rm i}\omega } \end{eqnarray}$

其中, i是虚数单位, $\omega $表示圆频率.

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

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

3.1 标量连分式

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

(13)$ \begin{eqnarray} \label{eq13} S=c\omega _n \sqrt{1+z^2+\beta z} \end{eqnarray}$

其中, $z={\rm i}\omega /\omega _n $, $\beta =\alpha /\omega _n $.

函数$\sqrt {1+z^2+\beta z} $有如下恒等式

\begin{eqnarray*} &&\sqrt {1+z^2+\beta z} -(1+z)=\\&&\qquad\frac{(\beta-2)z}{2(1+z)+\left[ {\sqrt {1+z^2+\beta z} -(1+z)}\right]} \end{eqnarray*}

可以将其写成连分式

(14)$ \begin{eqnarray} \label{eq14} &&\sqrt {1+z^2+\beta z} =(1+z)+\\&&\qquad\dfrac{(\beta -2)z}{2(1+z)+\dfrac{(\beta -2)z}{2(1+z)+\dfrac{(\beta -2)z}{2(1+z)+\ddots }}} \end{eqnarray}$

显然, $\sqrt {1+z^2+\beta z} $可以等价为若干种连分式,但不是所有的连分式都能保证在整个频率域都收敛^[19].本文提出的连分式(14)在整个频率域都是收敛的.经推导得到式(14)的非收敛域为$\left| z \right|=1$, Re$(z)<-0.5\beta $,Re($z)$表示$z$的实部.图2中左复平面的红线为连分式的非收敛域.可以看出在整个频率域$\omega =[0, \infty )$, $z$坐落在收敛域的虚轴上,所以连分式(14)对整个频率域都是收敛的.

图2

图2连分式的收敛域

Fig. 2The convergence domain of the continued fraction

为了直观看出连分式的收敛率, 定义总绝对误差

\begin{eqnarray*} \varepsilon =\frac{\sum\limits_{l=0}^N {\left(\left| {CF-\sqrt {1+z^2+\beta z} } \right|\right)} }{N+1} \end{eqnarray*}

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

图3

图3连分式的收敛率

Fig. 3The convergence rate of the continued fraction

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

(15)$ \begin{eqnarray} \left.\begin{array}{ll} \begin{array}{l} S=g_0 +{\rm i}\omega h_0 -{\rm i}\omega S_1^{-1} \\ S_j =g_j +{\rm i}\omega h_j -{\rm i}\omega S_{j+1}^{-1} \\ \end{array},& \ \ j=1, 2, {\cdots}, J \end{array}\ \right\} \end{eqnarray}$

其中, $S_{J+1}^{-1} =0$, $g_0 =c\omega _n$, $h_0 =c$.

当$j$是奇数时

\begin{eqnarray*} g_j =\frac{1}{h_0 -0.5\alpha /\omega _n h_0}, \ \ h_j =\frac{1}{g_0 -0.5\alpha /\omega _n g_0} \end{eqnarray*}

当$j$是偶数时

\begin{eqnarray*} g_j =2g_0, \ \ h_j =2h_0 \end{eqnarray*}

3.2 矩阵连分式

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

(16)$ \begin{eqnarray} {S}={g}_0 +{\rm i}\omega {h}_0 -{\rm i}\omega {S}_1^{-1} \end{eqnarray}$

(17)$ \begin{eqnarray} {S}_j ={g}_j +{\rm i}\omega {h}_j -{\rm i}\omega {S}_{j+1}^{-1},\ \ j=1, 2, {\cdots}, J \end{eqnarray}$

(18)$ \begin{eqnarray} {S}_{J+1}^{-1} ={\bf 0} \end{eqnarray}$

(19)$ \begin{eqnarray} {g}_0 =c{ \varLambda }_n ,\ \ {h}_0 =cI \end{eqnarray}$

当$j$是奇数时

(20)$ \begin{eqnarray} \label{eq18} \left.\begin{array}{l} {g}_j =({h}_0 -0.5{e}_4 { \varLambda }_n^{-1} {h}_0 )^{-1} \\ {h}_j =({g}_0 -0.5{e}_4 { \varLambda }_n^{-1} {g}_0 )^{-1} \end{array}\ \right\} \end{eqnarray}$

当$j$是偶数时

(21)$ \begin{eqnarray} \label{eq21} {g}_j =2{g}_0 , \ \ {h}_j =2{h}_0 \end{eqnarray}$

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

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

根据式(9) $\sim\!$式(11), 得到模态动力刚度方程

(22)$ \begin{eqnarray} \label{eq22} {Se}_1^{-1} {S}-{ \varLambda }_n^2 +\omega ^2I-{\rm i}\omega {e}_4 ={\bf 0} \end{eqnarray}$

将式(16)代入式(22), 并令i$\omega $的零次项和二次项分别为0, 得到关于$g_{0}$ 和$h_{0}$的如下方程

(23)$ \begin{eqnarray} \label{eq23} {g}_0 {e}_1^{-1} {g}_0 ={ \varLambda }_n^2, \ \ {h}_0 {e}_1^{-1} {h}_0 =I \end{eqnarray}$

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

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

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

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

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

将式(16) $\sim\!$式(21)代入式(11), 同时引进辅助变${S}_j {\tilde{{q}}}_j ={\tilde{{q}}}_{j-1} $, $j=1$, 2, ${\cdots}$, $J$, 得到

(24)$ \begin{eqnarray} \label{eq24} {\tilde{{p}}}={g}_0 {\tilde{{q}}}+{\rm i}\omega {h}_0 {\tilde{{q}}}-{\rm i}\omega {\tilde{{q}}}_1 \end{eqnarray}$

(25)$ \begin{eqnarray} \left.\begin{array}{l} {\tilde{{q}}}_{j-1} ={g}_j {\tilde{{q}}}_j +{\rm i}\omega {h}_j {\tilde{{q}}}_j -{\rm i}\omega {\tilde{{q}}}_{j+1},\ \ j=1, 2, {\cdots}, J\\ {\tilde{{q}}}_0 ={\tilde{{q}}}, \ \ {\tilde{{q}}}_{J+1} ={\bf 0} \end{array}\right\} \end{eqnarray}$

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

(26)$ \begin{eqnarray} \label{eq26} \left[ {{\begin{array}{*{20}c} {{c}_{\rm bb}^\infty } & {{c}_{\rm bc}^\infty } \\ {{c}_{\rm cb}^\infty } & {{c}_{\rm cc}^\infty } \\ \end{array} }} \right]\left\{ {{\begin{array}{*{20}c} {{\dot{{q}}}} \\ {{\dot{{q}}}_{\rm c}} \\ \end{array} }} \right\}+\left[ {{\begin{array}{*{20}c} {{k}_{\rm bb}^\infty } & {{k}_{\rm bc}^\infty } \\ {{k}_{\rm cb}^\infty } & {{k}_{\rm cc}^\infty } \\ \end{array} }} \right]\left\{ {{\begin{array}{*{20}c} {q} \\ {{q}_{\rm c}} \\ \end{array} }} \right\}=\left\{ {{\begin{array}{*{20}c} {p} \\ {{\bf 0}} \\ \end{array} }} \right\} \end{eqnarray}$

其中,

\begin{eqnarray*} &&\left[ {{\begin{array}{cc} {{c}_{\rm bb}^\infty } & {{c}_{\rm bc}^\infty } \\ {{c}_{\rm cb}^\infty } & {{c}_{\rm cc}^\infty } \\ \end{array} }} \right]=\left[ {{\begin{array}{c;{2pt/2pt}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c} {\ \ {h}_0\ \ } & {\ \ -I\ \ } & {{\bf 0}} & \cdots & {{\bf 0}} & {{\bf 0}} \\ \hdashline[2pt/2pt] {{\bf 0}} & {{h}_1 } & {-I} & \cdots & {{\bf 0}} & {{\bf 0}} \\ {{\bf 0}} & {{\bf 0}} & {{h}_2 } & \cdots & {{\bf 0}} & {{\bf 0}} \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ {{\bf 0}} & {{\bf 0}} & {{\bf 0}} & {{\bf 0}} & {{h}_{J-1} } & {-I} \\ {{\bf 0}} & {{\bf 0}} & {{\bf 0}} & {{\bf 0}} & {{\bf 0}} & {{h}_J } \\ \end{array} }} \right] \\&& \left[ {{\begin{array}{cc} {{k}_{\rm bb}^\infty } & {{k}_{\rm bc}^\infty } \\ {{k}_{\rm cb}^\infty } & {{k}_{\rm cc}^\infty } \\ \end{array} }} \right]=\left[ {{\begin{array}{c;{2pt/2pt}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c} {\ \ {g}_0 \ \ } & {\ \ {\bf 0}\ \ } & {{\bf 0}} & \cdots & {{\bf 0}} & {{\bf 0}} \\ \hdashline[2pt/2pt] {-I} & {{g}_1 } & {{\bf 0}} & \cdots & {{\bf 0}} & {{\bf 0}} \\ {{\bf 0}} & {-I} & {{g}_2 } & \cdots & {{\bf 0}} & {{\bf 0}} \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ {{\bf 0}} & {{\bf 0}} & {{\bf 0}} & {-I} & {{g}_{J-1} } & {{\bf 0}} \\ {{\bf 0}} & {{\bf 0}} & {{\bf 0}} & {{\bf 0}} & {-I} & {{g}_J } \\ \end{array} }} \right] \\&& {q}_{\rm c} =\left\{ {{\begin{array}{*{20}c} {{q}_1^{\rm T} } & {{q}_2^{\rm T} } & \ldots & {{q}_{J-1}^{\rm T} } & {{q}_J^{\rm T} } \\ \end{array} }} \right\}^{\rm T} \end{eqnarray*}

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

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

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

(27)$ \begin{eqnarray} \label{eq27} &&\left[ {{\begin{array}{*{20}c} {{C}_{\rm bb}^\infty } & {{C}_{\rm bc}^\infty } \\ {{C}_{\rm cb}^\infty } & {{C}_{\rm cc}^\infty } \\ \end{array} }} \right]\left\{ {{\begin{array}{*{20}c} {{\dot{{u}}}_{\rm b} } \\ {{\dot{{q}}}_{\rm c} } \\ \end{array} }} \right\}+\left[ {{\begin{array}{*{20}c} {{K}_{\rm bb}^\infty } & {{K}_{\rm bc}^\infty } \\ {{K}_{\rm cb}^\infty } & {{K}_{\rm cc}^\infty } \\ \end{array} }} \right]\left\{ {{\begin{array}{*{20}c} {{u}_{\rm b} } \\ {{q}_{\rm c}} \\ \end{array} }} \right\}=\left\{ {{\begin{array}{*{20}c} {{f}_{\rm b} } \\ {{\bf 0}} \\ \end{array} }} \right\}\end{eqnarray}$

$ \begin{eqnarray}\left[ {{\begin{array}{*{20}c} {{C}_{\rm bb}^\infty } & {{C}_{\rm bc}^\infty } \\ {{C}_{\rm cb}^\infty } & {{C}_{\rm cc}^\infty } \\\end{array} }} \right]=\left[ {{\begin{array}{c;{2pt/2pt}c@{\ }c@{\ }c@{\ }c@{\ }c} {\ \ {E}_3 { \varPhi }_n {h}_0 { \varPhi }_n^{\rm T} {E}_3\ \ } &{\ \ -{E}_3 { \varPhi }_n\ \ } & {{\bf 0}} & \cdots & {{\bf 0}} & {{\bf 0}} \\\hdashline[2pt/2pt] {{\bf 0}} & {{h}_1 } & {-I} & \cdots & {{\bf 0}} & {{\bf 0}} \\ {{\bf 0}} & {{\bf 0}} & {{h}_2 } & \cdots &{{\bf 0}} & {{\bf 0}} \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ {{\bf 0}} & {{\bf 0}} & {{\bf 0}} & {{\bf 0}} & {{h}_{J-1} } & {-I} \\ {{\bf 0}} & {{\bf 0}} & {{\bf 0}} & {{\bf 0}} & {{\bf 0}} & {{h}_J } \\\end{array} }} \right] \end{eqnarray} $

$ \begin{eqnarray} &&\left[ {{\begin{array}{*{20}c} {{K}_{\rm bb}^\infty } & {{K}_{\rm bc}^\infty } \\ {{K}_{\rm cb}^\infty } & {{K}_{\rm cc}^\infty } \\\end{array} }} \right]=\left[ {{\begin{array}{c;{2pt/2pt}c@{\ }c@{\ }c@{\ }c@{\ }c} {\ \ {E}_3 { \varPhi }_n {g}_0 { \varPhi }_n^{\rm T} {E}_3\ \ } &{\ \ {\bf 0}\ \ } & {{\bf 0}} & \cdots & {{\bf 0}} & {{\bf 0}} \\\hdashline[2pt/2pt] {-{ \varPhi }_n^{\rm T} {E}_3 } & {{g}_1 } & {{\bf 0}}& \cdots & {{\bf 0}} & {{\bf 0}} \\ {{\bf 0}} & {-I} & {{g}_2 } & \cdots & {{\bf 0}} & {{\bf 0}} \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ {{\bf 0}} & {{\bf 0}} & {{\bf 0}} & {-I} & {{g}_{J-1} } & {{\bf 0}} \\ {{\bf 0}} & {{\bf 0}} & {{\bf 0}} & {{\bf 0}} & {-I} & {{g}_J } \\\end{array} }} \right]\end{eqnarray} $

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

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

(28)$ \begin{eqnarray} \label{eq28} &&\left[ {{\begin{array}{*{20}c} {{M}_{\rm ii} } & {{M}_{\rm ib} } \\ {{M}_{\rm bi} } & {{M}_{\rm bb} } \\ \end{array} }} \right]\left\{ {{\begin{array}{*{20}c} {{\ddot{{u}}}_{\rm i} } \\ {{\ddot{{u}}}_{\rm b} } \\ \end{array} }} \right\}+\left[ {{\begin{array}{*{20}c} {{C}_{\rm ii} } & {{C}_{\rm ib} } \\ {{C}_{\rm bi} } & {{C}_{\rm bb} } \\ \end{array} }} \right]\left\{ {{\begin{array}{*{20}c} {{\dot{{u}}}_{\rm i} } \\ {{\dot{{u}}}_{\rm b} } \\ \end{array} }} \right\}+\\&&\qquad \left[ {{\begin{array}{*{20}c} {{K}_{\rm ii} } & {{K}_{\rm ib} } \\ {{K}_{\rm bi} } & {{K}_{\rm bb} } \\ \end{array} }} \right]\left\{ {{\begin{array}{*{20}c} {{u}_{\rm i} } \\ {{u}_{\rm b} } \\ \end{array} }} \right\}=\left\{ {{\begin{array}{*{20}c} {{f}_{\rm i} } \\ {-{f}_{\rm b} } \\ \end{array} }} \right\} \end{eqnarray}$

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

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

(29)$ \begin{eqnarray} \label{eq29} && \left[ {{\begin{array}{*{20}c} {{M}_{\rm ii} } & {{M}_{\rm ib} } & {{\bf 0}} \\ {{M}_{\rm bi} } & {{M}_{\rm bb} } & {{\bf 0}} \\ {{\bf 0}} & {{\bf 0}} & {{\bf 0}} \\ \end{array} }} \right]\left\{ {{\begin{array}{*{20}c} {{\ddot{{u}}}_{\rm i} } \\ {{\ddot{{u}}}_{\rm b} } \\ {{\ddot{{q}}}_{\rm c}} \\ \end{array} }} \right\}+\left[ {{\begin{array}{*{20}c} {{C}_{\rm ii} } & {{C}_{\rm ib} } & {{\bf 0}} \\ {{C}_{\rm bi} } & {{C}_{\rm bb} +{C}_{\rm bb}^\infty } & {{C}_{\rm bc}^\infty } \\ {{\bf 0}} & {{C}_{\rm cb}^\infty } & {{C}_{\rm cc}^\infty } \\ \end{array} }} \right]\left\{ {{\begin{array}{*{20}c} {{\dot{{u}}}_{\rm i} } \\ {{\dot{{u}}}_{\rm b} } \\ {{\dot{{q}}}_{\rm c}} \\ \end{array} }} \right\}+ \\ \\&&\qquad\left[ {{\begin{array}{*{20}c} {{K}_{\rm ii} } & {{K}_{\rm ib} } & {{\bf 0}} \\ {{K}_{\rm bi} } & {{K}_{\rm bb} +{K}_{\rm bb}^\infty } & {{K}_{\rm bc}^\infty } \\ {{\bf 0}} & {{K}_{\rm cb}^\infty } & {{K}_{\rm cc}^\infty } \\ \end{array} }} \right]\left\{ {{\begin{array}{*{20}c} {{u}_{\rm i} } \\ {{u}_{\rm b} } \\ {{q}_{\rm c}} \\ \end{array} }} \right\}=\left\{ {{\begin{array}{*{20}c} {{f}_{\rm i} } \\ {{\bf 0}} \\ {{\bf 0}} \\ \end{array} }} \right\} \end{eqnarray}$

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

5 数值算例

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

5.1 单层介质

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

(30)$ \begin{eqnarray} &&f(t)=16[Z(t/0.4)-4Z(t/0.4-0.25)+\\&&\qquad 6Z(t/0.4-0.5)-4Z(t/0.4-0.75)+\\&&\qquad Z(t/0.4-1)] \end{eqnarray}$

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

图4

图4单层数值算例

Fig. 4The numerical example of the single-layered medium

图5

图5狄拉克脉冲

Fig. 5Dirac pulse

如图4(b)所示, 引入人工边界, 截取有限域用有限元法计算, 网格尺寸为2 m $\times$ 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

图6单层介质的吸收边界条件的精度

Fig. 6Accuracy 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

图7单层介质的吸收边界条件的稳定性

Fig. 7Stability of proposed absorbing boundary condition for the single-layered medium

5.2 多层介质

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

图8

图8五层数值算例

Fig. 8The numerical example of the five-layered media

表1五层介质的材料参数

Table 1 The material parameters of the five-layered media

新窗口打开|下载CSV

分析方法与5.1节单层介质算例类似. 因为表示多层介质动力刚度的连分式是通过简单扩展单层介质的连分式获得的,必然会有误差引进. 为了使吸收边界条件对于多层介质具有高精度的特性,人工边界的位置到感兴趣区域的距离$L$需要被讨论, 分别取$L=0.5$ $H$和$L=0.75$ $H$.从上到下每层竖向网格的尺寸为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.5$ $H$, $L=0.75$ $H$的结果如图9所示(注:$L$表示人工边界到感兴趣区域的距离). 可以看出两种边界位置的结果几乎完全重合,说明$L=0.5$ $H$时所提出的吸收边界条件已具有很高的精度.图9还展示了$L=0.5$ $H$边界位置处黏弹性边界的结果, 可以看出黏弹性边界精度较低.

图9

图9五层介质的吸收边界条件的精度

Fig. 9Accuracy of proposed absorbing boundary condition for the five-layered media

(2)方法的稳定性

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

图10

图10五层介质的吸收边界条件的稳定性

Fig. 10Stability of proposed absorbing boundary condition for the five-layered media

6 结论

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

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

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

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

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

Kennett

Seismic Wave Propagation in Stratified Media. ANU E Press , 2009

[2]

杜修力

. 工程波动理论与方法. 北京 : 科学出版社 , 2009

( Du

Xiuli

. Theories and Methods of Wave Motion for Engineering. Beijing : Science Press , 2009 (in Chinese))

[3]

Lysmer

, Kuhlemeyer

Finite dynamic model for infinite media

Journal of the Engineering Mechanics Division, 1969 , 95 ( 4 ): 859 - 878

[4]

Higdon

Absorbing boundary conditions for acoustic and elastic waves in stratified media

Journal of Computational Physics, 1992 , 101 ( 2 ): 386 - 418

[5]

杜修力 , 赵密 , 王进廷 .

近场波动模拟的人工应力边界条件

力学学报, 2006 , 38 ( 1 ): 49 - 56

( Du

Xiuli

, Zhao

, Wang

Jintin

A stress artificial boundary in FEA for near-field wave problem

Chinese Journal of Theoretical and Applied Mechanics , 2006 , 38 ( 1 ): 49 - 56 (in Chinese))

[6]

Skelton

, Adams

SDM

, Craster

Guided elastic waves and perfectly matched layers

Wave Motion, 2007 , 44 ( 7-8 ): 573 - 592

[7]

Hagstrom

, Mar-Or

, Givoli

High-order local absorbing conditions for the wave equation: Extensions and improvements

Journal of Computational Physics, 2008 , 227 ( 6 ): 3322 - 3357

[8]

, Zhao

A local time-domain transmitting boundary for simulating cylindrical elastic wave propagation in infinite media

Soil Dynamics & Earthquake Engineering, 2010 , 30 ( 10 ): 937 - 946

[9]

章旭斌 , 谢志南 , 廖振鹏 .

SH波动模拟中透射边界反射放大失稳研究

哈尔滨工程大学学报, 2019 ( 6 ): 1 - 6

( Zhang

Xubin

, Xie

Zhinan

, Liao

Zhenpeng

Study on reflection amplification instability of transmitting boundary in SH wave simulation

Journal of Harbin Engineering University , 2019 ( 6 ): 1 - 6 (in Chinese))

[10]

Hagstrom

, Kim

Complete radiation boundary conditions for the Helmholtz equation I: Waveguides

Numerische Mathematik, 2019 , 141 ( 4 ): 917 - 966

[11]

Medvinsky

, Tsynkov

, Turkel

Direct implementation of high order BGT artificial boundary conditions

Journal of Computational Physics, 2019 , 376 : 98 - 128

[12]

Liu

, Xu

Discrete artificial boundary conditions for transient scalar wave propagation in a 2D unbounded layered media

Computer Methods in Applied Mechanics and Engineering, 2002 , 191 ( 27-28 ): 3055 - 3071

[13]

Prempramote

, Song

, Tin-Loi

, et al .

High-order doubly asymptotic open boundaries for scalar wave equation

International Journal for Numerical Methods in Engineering, 2009 , 79 ( 3 ): 340 - 374

[14]

Liu

, Zhao

Finite element modeling of wave propagation problems in multilayered soils resting on a rigid base

Computers and Geotechnics, 2010 , 37 ( 3 ): 248 - 257

[15]

, Zhao

Stability and identification for rational approximation of frequency response function of unbounded soil

Earthquake Engineering & Structural Dynamics}, 2010 , 39 ( 2 ): 165 - 186

[16]

Zhao

, Du

, Liu

, et al .

Explicit finite element artificial boundary scheme for transient scalar waves in two-dimensional unbounded waveguide

International Journal for Numerical Methods in Engineering, 2011 , 87 ( 11 ): 1074 - 1104

[17]

陈灯红 , 杜成斌 .

结构-地基动力相互作用的时域模型

岩土力学, 2014 ( 4 ): 1164 - 1172

( Chen

Denghong

, Du

Chengbin

A computational model for structure-foundation dynamic interaction in time domain

Rock and Soil Mechanics , 2014 ( 4 ): 1164 - 1172 (in Chinese))

[18]

Zhao

, Li

, Du

, et al .

Time-domain stability of artificial boundary condition coupled with finite element for dynamic and wave problems in unbounded media

International Journal of Computational Methods, 2019 , 16 ( 4 ): 1850099

[19]

Zhao

, Wu

, Du

, et al .

Stable high-order absorbing boundary condition based on new continued fraction for scalar wave propagation in unbounded multilayer media

Computer Methods in Applied Mechanics and Engineering, 2018 , 334 : 111 - 137

[本文引用: 7]

[20]

, Zhao

, Wu

, et al .

High-order absorbing boundary condition based on new continued fraction for transient scalar wave propagation in 2D and 3D unbounded layers

Engineering Computation, 2019 , 36 ( 7 ): 2445 - 2479

[21]

牛滨华 , 孙春岩 . 半无限空间各向同性黏弹性介质与地震波传播. 北京 : 地质出版社 , 2007

( Niu

Binhua

, Sun

Chunyan

. Viscoelastic Medium and Seismic Wave Propagation. Beijing : Geological Publishing House , 2007 (in Chinese))

[22]

Abu-Alshaikh

, Turhan

, Mengi

Two-dimensional transient wave propagation in viscoelastic layered media

Journal of Sound and Vibration, 2001 , 244 ( 5 ): 837 - 858

[23]

Sun

, Luo

Transient wave propagation in multilayered viscoelastic media: theory, numerical computation, and validation

Journal of Applied Mechanics, 2008 , 75 ( 3 ): 031007

[24]

Sun

, Pan

, Gu

High-order thin layer method for viscoelastic wave propagation in stratified media

Computer Methods in Applied Mechanics and Engineering, 2013 , 257 : 65 - 76

[25]

Basu

, Chopra

Perfectly matched layers for time-harmonic elastodynamics of unbounded domains: theory and finite-element implementation

Computer Methods in Applied Mechanics and Engineering, 2003 , 192 ( 11-12 ): 1337 - 1375

[26]

单启铜 , 乐友喜 .

PML 边界条件下二维粘弹性介质波场模拟

石油物探, 2007 , 46 ( 2 ): 126 - 130

( Shan

Qitong

, Yue

Youxi

Wavefield simulation of 2-D viscoelastic medium in perfectly matched layer boundary

Geophysical Prospecting for Petroleum , 2007 , 46 ( 2 ): 126 - 130 (in Chinese))

[27]

Martin

, Komatitsch

An unsplit convolutional perfectly matched layer technique improved at grazing incidence for the viscoelastic wave equation

Geophysical Journal International, 2009 , 179 ( 1 ): 333 - 344

[28]

Josifovski

Analysis of wave propagation and soil-structure interaction using a perfectly matched layer model

Soil Dynamics and Earthquake Engineering, 2016 , 81 : 1 - 13

[29]

Farzanian

, Arbabi

, Pak

PML solution of longitudinal wave propagation in heterogeneous media

Earthquake Engineering and Engineering Vibration, 2016 , 15 ( 2 ): 357 - 368

[30]

韩泽军 , 林皋 , 周小文

等 .

横观各向同性层状地基动应力响应的求解与分析

岩土力学, 2018 ( 6 ): 2287 - 2294

( Han

Zejun

, Lin

Gao

, Zhou

Xiaowen

, et al .

Solution and analysis of dynamic stress response for transversely isotropic multilayered soil

Rock and Soil Mechanics , 2018 ( 6 ): 2287 - 2294 (in Chinese))

[31]

王振宇 , 刘晶波 .

成层地基非线性波动问题人工边界与波动输入研究

岩石力学与工程学报, 2004 , 23 ( 7 ): 1169 - 1169

[本文引用: 4]

( Wang

Zhenyu

, Liu

Jingbo

Study on wave motion input and artificial boundary for problem of nonlinear wave motion in layered soil

Chinese Journal of Rock Mechanics and Engineering , 2004 , 23 ( 7 ): 1169 - 1169 (in Chinese))

[本文引用: 4]

[32]

Badry

, Ramancharla

Local absorbing boundary conditions to simulate wave propagation in unbounded viscoelastic domains

Computers & Structures, 2018 , 208 : 1 - 16

[本文引用: 3]

[33]

谢志南 , 郑永路 , 章旭斌

等 .

弱形式时域完美匹配层——滞弹性近场波动数值模拟

地球物理学报, 2019 , 62 ( 8 ): 3140 - 3154

( Xie

Zhinan

, Zheng

Yonglu

, Zhang

Xubin

, et al .

Weak-form time-domain perfectly matched layer for numerical simulation of viscoelastic wave propagation in infinite-domain

Chinese Journal of Geophysics , 2019 , 62 ( 8 ): 3140 - 3154 (in Chinese))

[34]

Ping

, Zhang

, Xu

, et al .

Efficiency of perfectly matched layers for seismic wave modeling in second-order viscoelastic equations

Geophysical Journal International, 2016 , 207 ( 3 ): 1367 - 1386

[35]

Poul

, Zerva

Time-domain PML formulation for modeling viscoelastic waves with Rayleigh-type damping in an unbounded domain: Theory and application in ABAQUS

Finite Elements in Analysis and Design, 2018 , 152 : 1 - 16

[36]

Sembalt

, Lenti

, Ali

A simple multi directional Absorbing Layer method to simulate elastic wave propagation in unbounded domains,

International Journal for Numerical methods in engineering, 2010 , 1 : 1 - 22

[37]

Rajagopal

, Drozdz

, Skelton

, et al .

On the use of absorbing layers to simulate the propagation of elastic waves in unbounded isotropic media using commercially available finite element packages,

NDT & E Int, 2012 , 51 : 30 - 40

[38]

尹训强 , 李建波 , 林皋 .

无限地基相互作用力时域求解的阻尼溶剂逐步抽取法

工程力学, 2015 , 32 ( 9 ): 20 - 26

( Yin

Xunqiang

, Li

Jianbo

, Lin

Gao

Damping solvent stepwise extraction method for evalution of interaction forces of unbounded soil in time domain,

Engineering Mechanics , 2015 , 32 ( 9 ): 20 - 26 (in Chinese))

[39]

Liu

, Li

The control-theory-based artificial boundary conditions for time-dependent wave guide problems in unbounded domain

Communications in Numerical Methods in Engineering, 2005 , 21 ( 12 ): 691 - 700

[40]

Kausel

Thin-layer method: Formulation in the time domain

International Journal for Numerical Methods in Engineering, 1994 , 37 ( 6 ): 927 - 941

[41]

Chen

, Birk

, Song

Transient analysis of wave propagation in layered soil by using the scaled boundary finite element method

Computers and Geotechnics, 2015 , 63 : 1 - 12

[42]

张金波 , 杨顶辉 , 贺茜君

等 .

求解双相和黏弹性介质波传播方程的间断有限元方法及其波场模拟

地球物理学报, 2018 , 61 ( 3 ): 926 - 937

( Zhang

Jinbo

, Yang

Dinghui

, He

Qianjun

, et al .

Discontinuous Galerkin method for solving wave equations in two-phase and viscoelastic media. Chinese Journal of

Geophysics-Chinese Edition , 2018 , 61 ( 3 ): 926 - 937 (in Chinese))

[43]

Chen

, Duan

, Li

Mass proportional damping in nonlinear time-history analysis//3rd International Conference on Material

Mechanical and Manufacturing Engineering (IC3ME 2015), Atlantis Press, 2015

[44]

王笃国 , 赵成刚 .

地震波斜入射下考虑场地非线性、地形效应和土结动力相互作用的大跨连续刚构桥地震响应分析

工程力学, 2017 , 34 ( 4 ): 37 - 46

( Wang

Duguo

, Zhao

Chenggang

Seismic analysis of long-span continuous rigid frame bridge considering site nonlinearity, topography effect and soil-structure dynamic interaction under oblique incidence

Engineering Mechanics , 2017 , 34 ( 4 ): 37 - 46 (in Chinese))