采用黏弹性人工边界单元时显式算法稳定性的改善研究1)

删除或更新信息，请邮件至freekaoyan#163.com(#换成@)

本站小编 Free考研考试/2022-01-01

李述涛^*^,^?, 刘晶波^*, 宝鑫^,^*^,²⁾

^*清华大学土木工程系,北京 100084

^?军事科学院国防工程研究院,北京 100036

IMPROVEMENT OF EXPLICIT ALGORITHMS STABILITY WITH VISCO-ELASTIC ARTIFICIAL BOUNDARY ELEMENTS ¹⁾

Li Shutao^*^,^?, Liu Jingbo^*, Bao Xin^,^*^,²⁾

^*Department of Civil Engineering,Tsinghua University,Beijing 100084,China

^?Institute of Defense Engineering,AMS,PLA,Beijing 100036,China

通讯作者: 2) 宝鑫,助理研究员,主要研究方向：地下结构抗爆抗震与海洋工程抗震. E-mail:baox@tsinghua.edu.cn

收稿日期:2020-06-2接受日期:2020-09-2网络出版日期:2020-11-18

基金资助:

1) 国家自然科学基金.51878384
国家自然科学基金.U1839201
国家重点研发计划.2018YFC1504305
博士后创新人才支持计划.BX20200192
清华大学“水木****”计划.2020SM005

Received:2020-06-2Accepted:2020-09-2Online:2020-11-18
作者简介 About authors

摘要
黏弹性人工边界单元是目前常用的处理半无限空间波动问题的数值模拟方法,可有效吸收计算区域内产生的外行波动.黏弹性人工边界单元具有与内部介质不同的质量密度、刚度和阻尼,受其影响,对整体模型进行显式时域逐步积分时,在边界区域易发生失稳现象,影响整体系统显式积分的计算效率. 针对该问题目前尚无行之有效的解决方法.本文针对二维黏弹性人工边界单元,建立可代表整体系统典型特征的侧边子系统和角点子系统,利用传递矩阵谱半径分析方法,基于传统中心差分格式,推导得到局部子系统稳定性条件的解析解.在此基础上通过研究解析解中各物理参数对稳定性条件的影响,给出通过增加人工边界单元的质量密度,以改善采用黏弹性人工边界单元时显式算法稳定性的方法.均匀和成层半空间波动问题算例分析表明,将内部单元质量密度设置为人工边界单元质量密度的上限,可以在保证黏弹性人工边界计算精度的前提下,有效改善整体系统显式时域逐步积分的数值稳定性,大幅提高计算效率.
关键词： 显式时域逐步积分;黏弹性人工边界单元;数值稳定性改善;传递矩阵;谱半径

Abstract
Viscoelastic artificial boundary elements are commonly applied in the analysis of semi-infinite wave propagation problems, which can accurately absorb the scattered waves generated in the calculation domain. However, the mass density, the stiffness and the damping coefficient of the viscoelastic artificial boundary element are different from those of the internal domain. Therefore, the instability often occurs in the boundary region when the explicit time-domain stepwise integration is performed in the overall model, so, the calculation efficiency of the explicit integral for the overall system is affected. Currently, there is no effective solution to this problem which remains to be settled to conduct efficient large-scale wave propagation simulation. For the two dimensional viscoelastic artificial boundary element, we establish the edge subsystem and corner subsystem which can represent the typical characteristics of the overall system, by the analysis method of transfer matrix spectral radius based on the central difference format commonly used, we derive the analytical solutions of the stability conditions of the edge subsystem and corner subsystem. After that, we analyze the influence of various physical parameters of the two dimensional viscoelastic artificial boundary element on the stability conditions, and obtain the method for improving the stability condition of the explicit algorithm by increasing the mass density of the viscoelastic artificial boundary element. The homogeneous and layered half-space examples show that, set the mass density of the internal element as the upper limit of the mass density of the viscoelastic artificial boundary element, the method proposed in this paper can effectively improve the numerical stability of the explicit time-domain integration when using the viscoelastic artificial boundary elements, without affecting the calculation accuracy, and the calculation efficiency can be significantly improved in the explicit dynamic analysis. The model size (distance from the scattered wave source to the artificial boundary) has a little effect on the stability of the explicit integral which can be ignored.
Keywords：explicit time domain integration;viscoelastic artificial boundary element;improvement of numerical stability;transfer matrix;spectral radius

PDF (12594KB)元数据多维度评价相关文章导出EndNote|Ris|Bibtex 收藏本文
本文引用格式
李述涛, 刘晶波, 宝鑫. 采用黏弹性人工边界单元时显式算法稳定性的改善研究¹⁾. 力学学报[J], 2020, 52(6): 1838-1849 DOI:10.6052/0459-1879-20-224
Li Shutao, Liu Jingbo, Bao Xin. IMPROVEMENT OF EXPLICIT ALGORITHMS STABILITY WITH VISCO-ELASTIC ARTIFICIAL BOUNDARY ELEMENTS ¹⁾. Chinese Journal of Theoretical and Applied Mechanics[J], 2020, 52(6): 1838-1849 DOI:10.6052/0459-1879-20-224

引言

无限域中波动问题的求解是众多物理领域所面对的重要研究课题. 针对此问题最初的分析方法一般是将截断边界取到足够远的位置,在截断边界产生的反射波返回感兴趣区域前完成计算过程^[1]. 这一方法具有较高的精度,但计算成本过高. 人工边界技术是目前常用的无限域介质数值模拟方法,将人工边界引入无限域空间的近场模型中,可吸收截断边界处的外行波动^[2-3]. 在众多的人工边界技术中,黏弹性人工边界基于全空间均匀介质中的柱面波和球面波理论推导得到,可等效为在每个边界节点上施加的远端固定的并联弹簧-阻尼器系统,既具有较高的计算精度,又可在多种有限元软件中进行集成,获得了较为广泛的工程应用^[4-5]. 但由于离散弹簧-阻尼器的施加需依赖程序语言对模型进行前处理,实用性受到一定影响. 为解决这一问题,刘晶波等^[6]和谷音等^[7]发展了黏弹性人工边界单元,其思路是利用有限单元的矩阵等效构造等效人工边界单元. 该方法仅需指定模型最外层单元的密度、弹性模量、泊松比和阻尼系数,即可模拟集中黏弹性人工边界,且计算精度不变.

随着黏弹性人工边界单元 (以下简称边界单元) 在实际工程中的广泛应用^[8-12],其数值稳定性问题随之而来. 研究人员在引入边界单元进行波动分析时,一般采用隐式无条件稳定的时域逐步积分算法,较少采用显式算法^[13-14]. 这是因为边界单元具有与内部介质不同的质量密度、刚度和阻尼,导致在边界单元区,显式算法的稳定条件不同于内部计算区. 由于边界单元区与相邻的内部计算区的力学参数完全不同,使得显式算法的稳定性分析遇到技术上的困难,长久以来未能给出与内部计算区类似的解析形式的数值积分稳定性条件,这在一定程度上影响了采用边界单元时对显式算法稳定性的判断,而数值算法稳定条件的缺失也限制了边界单元在显式动力分析中的应用.

显式算法的解耦特性回避了组集整体刚度阵和求解联立方程组等问题^[15-18],在大规模动力问题数值分析中的计算效率显著高于隐式算法,特别是对于大范围场地模型的并行计算问题,因此有必要对引入黏弹性人工边界单元时显式算法的稳定性开展研究.

数值稳定性的分析和改善研究一直伴随着人工边界理论的发展^[19-22]. 针对透射边界,可以改变与边界物理量相关的计算格式,通过滤波提高边界稳定性^[23]. 对于多次透射边界,可以在边界区设置黏性阻尼,以消除多次透射边界的高频振荡失稳现象^[24]. 黏弹性人工边界单元本身为物理上稳定的系统,将其引入计算系统中,并不影响整体的物理稳定条件. 其失稳机理是显式时域逐步积分算法所导致的,而与边界单元理论无关. 为改善数值稳定条件,需要全面考虑逐步积分过程中边界节点与内部单元节点之间的耦合效应,以此为依据调整满足显式时域积分数值稳定条件的临界时间步长,从而提高人工边界数值计算的稳定性.

本文针对人工边界典型位置处的局部子系统建立运动方程,利用传递矩阵谱半径分析方法,推导得到解析形式的局部子系统稳定性条件表达式.通过研究解析解中各物理参数对稳定性条件的影响,给出通过增加人工边界单元质量密度以实现改善显式算法稳定性条件的方法.本文研究对于提高黏弹性人工边界单元在动力显式计算中的适用性以及提高大范围工程场地波动分析的计算效率均具有重要意义.

1 黏弹性人工边界单元等效物理参数设置

黏弹性人工边界单元是在无限域半空间模型的截断边界上向外法线方向延伸出的一层厚度为$h$的单元,通过设置等效物理参数来模拟黏弹性人工边界. 二维黏弹性人工边界单元的等效密度 $\tilde{\rho }$、等效弹性模量 $\tilde {E}$、等效泊松比 $\tilde {\mu}$ 和等效阻尼系数 $\tilde {\eta }$ 分别为^[6]

(1)$ \left.\begin{array}{l} \tilde {\rho } = 0 , \ \ \ \tilde {E} = \alpha _{N} h\dfrac{G}{R} \\ \tilde {\mu } = 0 , \ \ \ \tilde {\eta } = \dfrac{\rho R}{2G}\Big (\dfrac{C_{S} }{\alpha_{T} } + \dfrac{C_{P} }{\alpha _{N} }\Big ) \end{array}\!\! \right \}$

其中,$G$ 为内部介质的剪切模量,$\rho $ 为内部介质的质量密度,$C_{S}$ 和 $C_{P}$ 分别为内部介质的 S 波和 P 波波速,$R$为散射波源至边界节点的距离,$\alpha _{T}$ 与 $\alpha _{N}$ 分别为切向与法向黏弹性人工边界参数,对于二维黏弹性人工边界,其推荐值分别为 0.5 和 1.

黏弹性人工边界单元厚度 $h$ 的鲁棒性良好,可灵活取值而对精度影响很小^[6]. 由于显式算法的稳定条件受模型中的最小单元尺寸制约,在实际计算分析中,建议将边界单元的厚度设定为与内部单元尺寸一致,如图1 所示,以排除单元尺寸对整体稳定性的影响.

图1

新窗口打开|下载原图ZIP|生成PPT
图1二维黏弹性人工边界单元示意图

Fig. 1Diagrams of 2D viscoelastic artificial boundary element

2 边界局部子系统稳定性分析

为改善数值积分稳定性,首先应对失稳机理开展研究. 若不考虑人工边界的影响,可使用冯诺依曼方法对计算区域进行稳定性分析^[25];引入边界单元并进行显式时域逐步积分计算时,则需要考虑边界单元对整体计算系统稳定性的影响.此时可以建立包含人工边界在内的整体系统的运动方程,采用传递矩阵谱半径分析方法,获得满足稳定条件的临界时间步长并进行失稳判断. 但在实际应用中,整体模型传递矩阵的建立和分析工作量巨大,通常难以施行.

针对透射人工边界条件,关慧敏等^[21]和 Kamel等^[26]提出了一种通过分析局部子系统传递矩阵特征值来研究透射人工边界稳定性的方法.由于该方法处理的子系统包含节点数量较多,只能给出系统稳定条件的数值解,难以给出解析解,不适合对稳定性条件的影响因素进行定量化的参数分析,制约了针对稳定性改善方法的研究.

针对黏弹性人工边界单元,李述涛等^[27]对一维和二维有限元系统最高阶模态进行了分析,证明了可以采用极小的子系统推导得出整体系统的稳定性条件.该子系统中只包含一个运动节点,所得的稳定性条件是整体系统数值稳定的充分条件.此外,该方法给出了系统稳定条件的解析解,可以直观且定量地分析各物理参数对数值稳定性的影响,为数值积分稳定性的改善创造了条件.

引入黏弹性人工边界单元后,可以选择对整体模型边界处具有典型特征的节点子系统开展研究. 一旦确定了子系统形式,即可根据显式逐步积分算法格式,将子系统的运动方程写成如下形式

(2)${\pmb U}^{p + 1} = {\pmb A \pmb U}^p$

其中向量 ${\pmb U}$ 由子系统各节点的加速度、速度和位移组成;${\pmb A}$ 是传递矩阵,上标 $p$ 是自然数,$t = p\Delta t$;$\Delta t$ 为时间步长. 如果满足下列两条件,则积分格式 (2) 是稳定的^[28].条件 1：$\rho ({\pmb A}) \leqslant 1$,$\rho ({\pmb A})$ 是传递矩阵 ${\pmb A}$ 的谱半径;条件 2：如果 ${\pmb A}$ 具有多重特征值,则该特征值的模小于 1.由此可将子系统的稳定性分析归结为传递矩阵 ${\pmb A}$ 的谱半径计算.

图2所示是均匀离散的二维有限元整体模型中的两处典型节点子系统. 其一位于模型的侧边或底边,由 2 个内部介质单元和 2 个边界单元组成. 另一子系统位于模型角点处,由 1 个内部介质单元和 3 个边界单元组成.

图2

新窗口打开|下载原图ZIP|生成PPT
图2二维模型中的两种节点子系统

Fig. 2Two kinds of nodal subsystems in 2D model

文献^[27] 采用上述研究思路,将侧边子系统与角点子系统从整体模型中分离出来,在四周施加固定约束. 利用传递矩阵的谱半径分析时域逐步积分算法的稳定性条件,最终推导出整体系统的稳定性条件. 分析结果表明,采用边界单元时,整体模型数值积分的稳定性条件要比内部系统的稳定性条件更为严格,这在一定程度上限制了显式数值积分方法的应用.

现有数值积分方法的稳定性研究表明,当介质的质量密度增大时 (相当于介质波速降低或结构自振周期变长),系统的数值积分稳定性条件将得到改善. 在原始的黏弹性人工边界单元中,等效质量密度 $\tilde {\rho } = 0$,但由于边界单元具有较良好的鲁棒性,适当调整参数对人工边界的模拟精度影响不大. 因此可以通过适当增加边界单元的质量密度以达到改善整体模型数值积分稳定性的目的.

下面首先给出考虑边界单元质量密度时,角点子系统和侧边子系统的稳定性条件.

2.1 角点子系统稳定性条件

角点子系统如图3 所示. 设内部介质单元的弹性模量为 $E$,剪切模量为 $G$,密度为 $\rho $,泊松比为 $\mu $,且无阻尼;黏弹性人工边界单元的密度为 $\tilde {\rho }$,弹性模量为 $\tilde {E}$,阻尼系数为 $\tilde {\eta }$,泊松比为 0,单元边长均为 $L$.

图3

新窗口打开|下载原图ZIP|生成PPT
图3角点子系统

Fig. 3Corner subsystem

二维平面单元的集中质量矩阵、刚度矩阵、阻尼矩阵分别参考文献 ^[6,29-30],按节点编号组装矩阵后,得到角点子系统中 5 号节点的运动方程如下

(3)$ \dfrac{(\rho + 3\tilde {\rho })L^2}{4}\left[ \!\!\begin{array}{cc} 1 \!&\! 0 \\ 0 \!&\! 1 \end{array} \!\! \right]\left\{ \!\! \begin{array}{c} {\ddot {u}_x } \\ {\ddot {u}_y } \end{array} \!\!\right\} + \dfrac{\rho R\tilde {E}}{2G}(2C_{S} + C_{P} )\left[\!\! \begin{array}{cc} {\dfrac{3}{2}} \!&\! {\dfrac{1}{8}} \\ {\dfrac{1}{8}} \!&\! {\dfrac{3}{2}} \end{array} \!\!\right]\left\{ \!\!\begin{array}{c} {\dot {u}_x } \\ {\dot {u}_y } \end{array} \!\! \right\} +\\ \left(\!\! \begin{array}{cc} {\dfrac{E(3 - \mu )}{6(1 - \mu ^2)} + \dfrac{3\tilde {E}}{2}} & {\dfrac{ - E}{8(1 - \mu )} + \dfrac{\tilde {E}}{8}} \\ {\dfrac{ - E}{8(1 - \mu )} + \dfrac{\tilde {E}}{8}} & {\dfrac{E(3 - \mu )}{6(1 - \mu ^2)} + \dfrac{3\tilde {E}}{2}} \end{array} \!\! \right)\left\{ \!\! \begin{array}{c} {u_x } \\ {u_y } \end{array} \!\!\right\} = 0$

其中,$\ddot {u}$,$\dot {u}$,$u$ 分别为节点的加速度、速度和位移,下标 $x$,$y$ 表示运动方向. 显式时域逐步积分格式 (中心差分) 如下^[31]

(4)$ \left.\begin{array}{l} \dot {u}\Big (t + \dfrac{\Delta t}{2}\Big) = \dot {u}\Big (t - \dfrac{\Delta t}{2}\Big) + \Delta t\ddot {u}(t) \\ u(t + \Delta t) = u(t) + \Delta t\dot {u}\Big (t + \dfrac{\Delta t}{2}\Big) \end{array} \!\! \right \}$

根据式(4)将式(3)展开,并写成传递矩阵形式

(5)$ \left\{ \begin{array}{c} {u(t + \Delta t)} \\ {v(t + \Delta t)} \\ {\Delta t\dot {u}\Big (t + \dfrac{\Delta t}{2}\Big)} \\ {\Delta t\dot {v}\Big (t + \dfrac{\Delta t}{2}\Big)} \end{array} \right\} = {\pmb A}_1 \left\{ \begin{array}{c} {u(t)} \\ {v(t)} \\ {\Delta t\dot {u}\Big (t - \dfrac{\Delta t}{2}\Big)} \\ {\Delta t\dot {v}\Big (t - \dfrac{\Delta t}{2}\Big)} \end{array} \right\}$

其中

(6)${\pmb A}_1 = \left[ \!\! \begin{array}{cccc} {1 - \dfrac{k_1 }{m_1 }\Delta t^2} & { - \dfrac{k_2 }{m_1 }\Delta t^2} & {1 - \dfrac{3c\Delta t}{2m_1 }} & { - \dfrac{c\Delta t}{8m_1 }} \\ { - \dfrac{k_2 }{m_1 }\Delta t^2} & {1 - \dfrac{k_1 }{m_1 }\Delta t^2} & { - \dfrac{c\Delta t}{8m_1 }} & {1 - \dfrac{3c\Delta t}{2m_1 }} \\ { - \dfrac{k_1 }{m_1 }\Delta t^2} & { - \dfrac{k_2 }{m_1 }\Delta t^2} & {1 - \dfrac{3c\Delta t}{2m_1 }} & { - \dfrac{c\Delta t}{8m_1 }} \\ { - \dfrac{k_2 }{m_1 }\Delta t^2} & { - \dfrac{k_1 }{m_1 }\Delta t^2} & { - \dfrac{c\Delta t}{8m_1 }} & {1 - \dfrac{3c\Delta t}{2m_1 }} \end{array} \!\!\right]$

(7)$\left.\!\!\begin{array}{l} m_1 = \dfrac{(\rho + 3\tilde {\rho })L^2}{4} , \ \ c = \dfrac{\rho L}{2}(2C_{S} + C_{P} ) \\ k_1 = \rho C_{S}^2\left( {\dfrac{(3 - \mu )}{3(1 - \mu )} + \dfrac{3L}{2R}} \right) , \ \ k_{2} = \dfrac{\rho C_{S}^2}{8}\left( {\dfrac{L}{R} - \dfrac{2(1{ + }\mu )}{1 - \mu }} \right) \end{array} \!\!\right\}$

传递矩阵 ${\pmb A}_{1}$ 的谱半径 $\rho ({\pmb A}_1 )$ 的表达式为

(8)$ \rho ({\pmb A}_1 ) = \dfrac{1}{8}\left| {8 - \dfrac{4(k_1 + k_2 )}{m_1 }\Delta t^2 - \dfrac{13c\Delta t}{2m_1 } - } \right. \\ \left. {\sqrt { - 8\left( {8 - \dfrac{13c\Delta t}{m_1 }} \right) + \left[ { - 8 + \dfrac{4(k_1 + k_2 )}{m_1 }\Delta t^2 + \dfrac{13c\Delta t}{2m_1 }} \right]^2} } \right|$

根据稳定性判别条件 $\rho ({\pmb A}_1 ) \leqslant 1$,得出角点子系统稳定性条件的解析表达式为

(9)$ \Delta t \leqslant \dfrac{\sqrt {256m_1 (k_1 + k_2 ) + (13c)^2} - 13c}{8(k_1 + k_2 )}$

将式 (7) 代入式 (9),整理得到

(10)$ \dfrac{\Delta tC_{P} }{L} \leqslant \gamma _1$

其中,$\gamma _1 $ 为角点子系统稳定性参数,整理后如式 (11) 所示.

由式 (10) 和式 (11) 可以发现,当内部介质的 P 波波速 $C_{P}$ 和单元尺寸 $L$ 给定后,角点子系统数值稳定条件仅与质量比 $\tilde {\rho }/\rho $、内部介质的泊松比 $\mu$ 和单元尺寸与散射波源至边界距离之比 $L/R$ 有关. 由于内部介质的泊松比 $\mu \in \left[ {0,0.5}\right]$,由式 (11) 可以清楚的看出,增大人工边界单元的质量密度 $\tilde {\rho}$ 将有助于改善数值稳定性条件.

(11)$\begin{array}{l}\gamma _1 = \Bigg \{ \sqrt 3 (1 - \mu )\Bigg \{3\left( {603 + 338\sqrt 2 \sqrt {1 - 2\mu } \sqrt {1 - \mu } } \right) + 48(1 - 2\mu )\left\{ {13\left( {1 - \mu } \right)\dfrac{L}{R} + \left[ {39(1 - \mu )\dfrac{L}{R} + 2(9 - 7\mu )} \right]\dfrac{\tilde {\rho }}{\rho }} \right\} +\\ \mu \left( { - 4856 - 1014\sqrt 2 \sqrt {1 - 2\mu } \sqrt {1 - \mu } + 2983\mu } \right)\Bigg \}^{\tfrac 12} - 39\left( {1 - \mu } \right)^{\tfrac{3}{2}}\left( {\sqrt 2 \sqrt {1 - 2\mu } + \sqrt {1 - \mu } } \right) \Bigg \} \Bigg / \\ \Bigg [ {2\left( {9 - 7\mu } \right)\left( {1 - 2\mu } \right){ + }39\left( {1 - \mu } \right)\left( {1 - 2\mu } \right)\dfrac{L}{R}} \Bigg ]\end{array}$

2.2 侧边子系统稳定性条件

侧边子系统如图4 所示,其模型参数与角点子系统相同. 按节点编号组装矩阵后,得到侧边子系统中 5 号节点的运动方程如下

(12)$ \dfrac{(\rho + \tilde {\rho })L^2}{2}\left[\!\! \begin{array}{cc} 1 & 0 \\ 0 & 1 \end{array}\!\! \right]\left\{ \!\! \begin{array}{c} {\ddot {u}_x } \\ {\ddot {u}_y } \end{array}\!\! \right\} + \dfrac{\rho R\tilde {E}}{2G}(2C_{S} + C_{P} )\left[\!\! \begin{array}{cc} 1 & 0 \\ 0 & 1 \end{array} \!\! \right]\cdot \\ \left\{ \!\! \begin{array}{c} {\dot {u}_x } \\ {\dot {u}_y } \end{array}\!\! \right\} +\left[ {\dfrac{E(3 - \mu )}{3(1 - \mu ^2)} + \tilde {E}} \right]\left(\!\! \begin{array}{cc} 1 & 0 \\ 0 & 1 \end{array}\!\! \right)\left\{ \!\! \begin{array}{c} {u_x } \\ {u_y } \end{array} \!\! \right\} = 0 $

图4

新窗口打开|下载原图ZIP|生成PPT
图4侧边子系统

Fig. 4Edge subsystem

参照上一节方法,将式 (12) 写成传递矩阵的形式

(13)$ \left\{ \!\! \begin{array}{c} {u(t + \Delta t)} \\ {\Delta t\dot {u} \Big (t + \dfrac{\Delta t}{2} \Big)} \end{array}\!\! \right\} = {\pmb A}_2 \left\{ \!\! \begin{array}{c} {u(t)} \\ {\Delta t\dot {u} \Big (t - \dfrac{\Delta t}{2} \Big)} \end{array}\!\! \right\}$

其中

(14)${\pmb A}_2 = \left[\!\!\begin{array}{cc} {1 - \dfrac{k\Delta t^2}{m_2 }} & {1 - \dfrac{c\Delta t}{m_2 }} \\ { - \dfrac{k\Delta t^2}{m_2 }} & {1 - \dfrac{c\Delta t}{m_2 }} \end{array}\!\!\right] $

(15)$\left. \!\!\begin{array}{l} m_2 = \dfrac{(\rho + \tilde {\rho })L^2}{2} , \ \ c = \dfrac{\rho L}{2}(2C_{S} + C_{P} ) \\ k = \rho C_{S} ^2\left[ {\dfrac{2(3 - \mu )}{3(1 - \mu )} + \dfrac{L}{R}} \right] \end{array}\!\!\right\} $

计算传递矩阵 ${\pmb A}_{2}$ 的特征值,并根据稳定性判别条件 $\rho ({\pmb A}_2 ) \leqslant 1$,得出侧边子系统稳定性条件的解析表达式为

(16)$\Delta t \leqslant \dfrac{1}{k}(\sqrt {4km_2 + c^2} - c)$

将式 (15) 代入式 (16),整理得到

(17)$\dfrac{\Delta tC_{P} }{L} \leqslant \gamma _2$

其中 $\gamma _2 $ 为侧边子系统稳定性参数,整理后如式 (18) 所示.

由式 (18) 可知,侧边子系统稳定性系数 $\gamma_{2}$ 同样为一个无量纲的解析表达式,可以用来定量讨论质量密度比 $\tilde {\rho }/\rho $ 的变化对数值稳定性条件的影响.

(18)$\gamma _2 = \dfrac{{3}(\mu - 1)^2\left[ {\sqrt {\dfrac{8(3 - \mu ) + 12(1 - \mu )\dfrac{L}{R}}{3(1 - \mu )^2}(1 - 2\mu )\left( {1 + \dfrac{\tilde {\rho }}{\rho }} \right) + \left( {1 + \sqrt {\dfrac{2 - 4\mu }{1 - \mu }} } \right)^2} - \sqrt {\dfrac{2 - 4\mu }{1 - \mu }} - 1} \right]}{2(3 - \mu )(1 - 2\mu ) + 3(1 - \mu )(1 - 2\mu )\dfrac{L}{R}}$

3 数值计算稳定性条件的改善方法

以上通过理论推导,分别给出了角点和侧边子系统数值积分稳定性条件,其数值稳定性系数 $\gamma_{1}$ 和 $\gamma_{2}$ 均为质量密度比 $\tilde {\rho }/\rho $、内部介质泊松比 $\mu $、单元尺寸与散射波源至边界距离之比 $L/R$ 的函数. 当无量纲参数 ${\Delta tC_{P}} /L$ 小于或等于稳定性系数时,相应的子系统满足数值积分的稳定性条件. 此外,整体系统的数值稳定性除需满足人工边界区的角点子系统和侧边子系统的稳定条件外,还需满足内部系统的稳定性条件. 显式时域逐步积分算法中内部单元的稳定性条件为^[25]

(19)$\dfrac{\Delta tC_{P} }{L} \leqslant \gamma _0$

其中,$\gamma _0 = 1$.

本节首先讨论人工边界单元质量密度 $\tilde {\rho }$ 和距离 $R$ 对数值稳定性的影响,并进一步给出改善数值稳定条件的方法.

3.1 人工边界单元质量密度对稳定性条件的影响

图5和图6分别给了出角点子系统和侧边子系统数值积分稳定性系数 $\gamma _{1}$ 和 $\gamma_{2}$ 随单元质量密度比 $\tilde {\rho }/\rho $ 的变化曲线. 其中单元尺寸与散射波源至边界距离之比 $L/R$ 分别为 1 和 1/100,内部介质泊松比 $\mu $ 分别为 0,0.25,0.3,0.4 和 0.5. 由图5 和图6 可以看出,增大质量密度比 $\tilde {\rho }/\rho $ 可以提高稳定性系数. 对于角点和侧边两种子系统,内部介质泊松比 $\mu $ 的取值越大,稳定性越宽松.

图5

新窗口打开|下载原图ZIP|生成PPT
图5角点子系统稳定性系数随密度比变化情况

Fig. 5The variation of corner subsystem stability coefficient with density ratio

图6

新窗口打开|下载原图ZIP|生成PPT
图6侧边子系统稳定性系数随密度比变化情况

Fig. 6The variation of edge subsystem stability coefficient with density ratio

整体系统的数值稳定条件应同时满足人工边界子系统和内部单元的稳定性条件. 图7 给出了当 $L/R =1$ 和 $\mu = 0.25$ 时,整体系统的数值积分稳定条件. 图中阴影区域即为满足整体系统数值积分稳定性的区域.

图7

新窗口打开|下载原图ZIP|生成PPT
图7整体系统数值稳定区域

Fig. 7The region of numerical stability in the overall system

以上分析可以得出以下两点结论：(1) 人工边界区的数值积分稳定性条件由角点子系统控制;(2) 采用黏弹性人工边界单元时,增大边界单元的质量密度可以有效改善显式数值积分的稳定性条件.

3.2 距离 ${\pmb R}$ 对稳定性条件的影响

影响人工边界子系统稳定性的另外一个因素是散射波源至人工边界的距离 $R$. 由于稳定性系数计算式 (11) 和式 (18) 的分子和分母中均含有反映距离 $R$ 影响的无量纲系数 $L/R$,因此,难以像质量密度 $\tilde {\rho}$ 一样,直接给出对稳定性影响效果的明确判断. 图8 和图9 给出了泊松比等于 0.25,质量密度比分别为 0 和 1 时,3 种系统稳定性系数 $\gamma_{0}$,$\gamma _{1}$ 和 $\gamma_{2}$ 随 $R/L$ 的变化曲线.

图8

新窗口打开|下载原图ZIP|生成PPT
图8稳定性系数随$R/L$的变化情况

Fig. 8The variation of stability coefficient with $R/L$

图9

新窗口打开|下载原图ZIP|生成PPT
图9均匀半空间模型

Fig. 9Homogeneous half space model

从图8 中可以看出,当散射波源至人工边界的距离与单元尺寸之比 $R/L$ 增大时,角点和侧边两种子系统的稳定性参数只在初始阶段略为增大,当 $R$ 大于约 5 倍的 $L$ 后,稳定性系数基本保持不变.因此,增大模型尺寸对整体系统数值稳定性的改善几乎不产生影响. 另外,当泊松比 $\mu$ 取为其他值时,可以得到完全相同的结论.

3.3 改进的黏弹性人工边界单元

基于以上分析,可确定改善显式数值算法稳定性具体措施 (方法) 和原则：(1) 通过增大黏弹性人工边界单元的质量密度达到改善算法稳定性的目的;(2) 保证人工边界区子系统的稳定性条件达到或略宽于内部单元系统稳定性条件,使得整体系统的稳定性由内部系统控制和判断.

改进的黏弹性人工边界单元参数设置如下：设内部介质剪切模量为$G$,质量密度为 $\rho$,剪切波速和压缩波速分别为 $C_{S}$ 和 $C_{P}$;内部介质单元厚度和边界单元厚度均为 $L$,边界节点至散射波源的距离为 $R$,边界单元质量密度为 $\tilde {\rho }$;$\alpha _{T} $ 与 $\alpha _{N} $ 分别为切向与法向黏弹性人工边界参数,二维模型建议 $\alpha _{T} $ 取 0.5,$\alpha _{N} $ 取 1.0. 稳定性条件改善后,二维黏弹性人工边界单元质量密度 $\tilde {\rho }$、等效刚度 $\tilde {E}$、泊松比 $\tilde {\mu }$ 和等效阻尼 $\tilde {\eta }$ 分别为

(20)$ \left.\begin{array}{l} \tilde {\rho } = \rho \\ \tilde {E} = \alpha _{N} h\dfrac{G}{R} \cdot \dfrac{(1 + \tilde {\mu })(1 - 2\tilde {\mu })}{1 - \tilde {\mu }}\\ \tilde {\mu } = 0\\ \tilde {\eta } = \dfrac{\rho R}{2G}\Big (\dfrac{C_{S} }{\alpha _{T} } + \dfrac{C_{P} }{\alpha _{N} } \Big ) \end{array} \!\! \right \}$

稳定性改善后,整体系统稳定性条件基本受内部介质区域控制,可采用最小单元尺寸与内部介质压缩波速的比值估算满足稳定性条件的临界时间步长.

4 稳定性数值算例及精度分析

4.1 均匀半空间模型

建立如图9 所示的均匀半空间模型. 模型尺寸为 320 m $\times$ 160 m,内部介质的密度为 2500 kg/m$^{3}$,剪切波速为 500 m/s,泊松比为 0.3,有限元网格尺寸为 2 m $\times$ 2 m,模型侧边和底边设置黏弹性人工边界单元,选取图中 $A, B, C,D$ 四点作为观测点,考察边界单元参数变化后,模型节点位移的计算精度. 采用持时为 0.2 s 的脉冲作为动力载荷,沿 $y$ 轴正向施加于模型中点处,载荷时程曲线如图10 所示.

图10

新窗口打开|下载原图ZIP|生成PPT
图10动力载荷时程曲线

Fig. 10Time history of the dynamic load

表1给出了均匀半空间模型稳定性改善前后临界时间步长比较情况. 稳定性条件改善前,边界单元的质量密度为 0,整体系统稳定性受角点子系统控制,按式 (9) 计算得到角点子系统满足稳定性条件的临界时间步长为 0.62 ms. 按照此时间步长可顺利完成模型动力计算. 而如果将时间步长增大,例如取 0.85 ms,计算发生失稳.

Table 1
表1
表1均匀半空间模型稳定性条件改善前后临界时间步长比较
Table 1Comparison of critical time steps before and after the improvement of the stability in the homogeneous half space model

新窗口打开|下载CSV

采用改进的黏弹性人工边界单元,即将人工边界单元密度设为与内部介质一致,其值为2500 kg/m$^{3}$,整体模型稳定性变为受内部区域控制,采用内部介质稳定性条件 ($\Delta t= 2.14$ ms) 可顺利完成对整体模型的显式时域逐步积分计算,不会发生失稳.从图11 可以看出,稳定性条件改善后,观测点 $A, B, C, D$ 四点的位移时程曲线与改善前相比吻合良好,说明边界单元自身质量密度的适当增大对波动问题计算精度的影响很小.另一方面,由于数值积分的时间步长比改善前增大了 243.3%,总计算时间相较于原时长缩短了 70.9% 左右,大幅提高了显式时域逐步积分的计算效率.

图11

新窗口打开|下载原图ZIP|生成PPT
图11均匀介质模型中 4 个观测点 $y$ 方向位移比较

Fig. 11Comparison of displacements in $y$ direction of four nodes in the homogeneous model

4.2 成层半空间模型

建立如图12 所示的成层半空间有限元模型. 模型尺寸为 320 m $\times $ 160 m,上层介质密度为 2000 kg/m$^{3}$, 剪切波速为 300 m/s,泊松比为 0.27,下层介质密度为 2500 kg/m$^{3}$,剪切波速为 500 m/s,泊松比为 0.3,整体模型的有限元网格尺寸为 2 m $\times $ 2 m. 脉冲载荷施加于模型的几何中心,沿 $y$ 轴正向加载,时程曲线如图6 所示. 通过令边界单元的质量密度与内部区域相等的方法对稳定性进行改善. 同样选取图中 $A, B, C, D$ 四点作为观测点,考察稳定性改善方法的可靠性和改善后半无限域空间波动问题的计算精度.

图12

新窗口打开|下载原图ZIP|生成PPT
图12成层半空间模型

Fig. 12Layered half space model

表2对成层介质中各子系统和内部区域稳定性改善前后的临界时间步长进行了比较.从结果可以看出,上层介质的稳定性条件比下层介质宽松.由于稳定性条件受最不利情况控制,因此成层介质模型整体系统的稳定性条件由下层介质所控制. 本算例中的下层介质参数与 4.1 节的均匀半空间模型相同,因此整体系统稳定性条件也与上节算例相同.

Table 2
表2
表2成层半空间模型稳定性条件改善前后临界时间步长比较
Table 2Comparison of critical time steps before and after the improvement of the stability in the layered half space model

新窗口打开|下载CSV

按照稳定性改善前后的临界时间步长对成层半空间模型进行动力分析,结果如图13 所示. 图中 $A, B, C,D$ 四点位移时程曲线吻合良好,再次验证了边界单元自身质量密度的适当增大对波动问题计算精度的影响很小.对于成层半空间模型,显式时域逐步积分的总计算时间相较于原时长也缩短了 70.9%,计算效率显著提高.

图13

新窗口打开|下载原图ZIP|生成PPT
图13成层介质模型中四个观测点 $y$ 方向位移比较

Fig. 13Comparison of displacements in $y$ direction of four nodes in the layered model

5 结论

本文针对采用黏弹性人工边界单元时显式时域逐步积分算法稳定性改善问题,推导得到了含黏弹性人工边界单元在内的有限元系统在采用显式时域逐步积分 (中心差分) 时数值稳定条件的解析解.研究了各物理参数对稳定性条件的影响,提出了改善数值积分稳定性的方法,并通过算例验证了这一改善方法的有效性和对黏弹性人工边界单元计算精度的影响. 结论如下：

(1) 采用有适当质量的黏弹性人工边界单元,可在保证无限域空间波动问题计算精度的前提下,大幅增大满足显式时域逐步积分稳定性条件的临界时间步长,达到改善显式数值积分稳定性的目的.

(2) 可将内部单元密度设置为边界单元密度的上限,此时整体系统的积分稳定性条件受内部介质区域控制,可通过模型中最小单元尺寸与内部介质压缩波速的比值来估算满足稳定性条件的临界时间步长.

(3) 黏弹性人工边界单元中的物理参数 $R$,即散射波源至人工边界的距离,对显式时域逐步积分算法稳定性的影响不大,当 $R$ 很小时,$R$ 的增大可以略微改善显式算法的稳定性,当 $R$ 大于约 5 倍的单元边长时,其对显式算法稳定性的影响基本不再变化.

(4) 如果显式时域逐步积分格式发生变化,依然可以使用本文提出的传递矩阵谱半径分析方法,按照不同积分格式将节点子系统运动方程展开,利用谱半径不大于 1 的判别条件给出稳定性结论,研究分析适用的稳定性改善方法.

参考文献原文顺序
文献年度倒序
文中引用次数倒序
被引期刊影响因子

[1]

Alterman

, Aboudi

, Karal

. Pulse propagation in a laterally heterogeneous solid elastic sphere
Geophysical Journal of the Royal Astronomical Society, 1970,21(3):243-260

[本文引用: 1]

[2]

杜修力, 赵密, 王进廷. 近场波动模拟的人工应力边界条件
力学学报, 2006,38(1):49-56

URL [本文引用: 1]

采用平面波和远场散射波混合透射，引入无限介质线弹性本构关系建立了一种应力人工边界条件. 其优点在于边界结点反应与内部有限元结点反应采用相同的积分格式计算，有限元积分方法稳定时不存在人工边界失稳问题. 数值算例表明：边界精度高于现有黏性边界、黏弹性人工边界，以及一、二阶透射人工边界.

( Du

Xiuli

, Zhao

, Wang

Jinting

. 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))

URL [本文引用: 1]

The stability and accuracy of local artificial boundary are two important criterions for evaluating its applied value in engineering. In this paper, a stress artificial boundary condition was developed by mixing the planar waves and the scattered waves transmitting into the linear elastic infinite media through artificial boundary. For the proposed boundary condition, the responses of boundary nodes and finite element nodes are calculated by the same integral method, and its stability condition is also the same as the finite element integral method. Numerical examples demonstrate that the accuracy of the proposed boundary is higher than that of existing viscous boundaries, viscous-spring boundaries, the first-order transmitting boundary and the second-order transmitting boundary.

[3]

邢浩洁, 李鸿晶. 透射边界条件在波动谱元模拟中的实现:一维波动
力学学报, 2017,49(2):367-379