采用黏弹性人工边界时显式算法稳定性条件

刘晶波; 宝鑫; 李述涛; 王菲

doi:10.11883/bzycj-2021-0196

采用黏弹性人工边界时显式算法稳定性条件

doi: 10.11883/bzycj-2021-0196

刘晶波^1,,
宝鑫^1, ,,
李述涛^{1, 2},
王菲^{1, 3}

1.
清华大学土木工程系，北京 100084
2.
军事科学院国防工程研究院，北京 100036
3.
陆军工程大学国防工程学院，江苏南京 210007

基金项目: 国家自然科学基金（U1839201，51878384）；国家重点研发计划（2018YFC1504305）；中国博士后科学基金（2020M680575）；博士后创新人才支持计划（BX20200192）；清华大学“水木学者”计划（2020SM005）

详细信息

作者简介:
刘晶波（1956－　），男，博士，教授，liujb@tsinghua.edu.cn

通讯作者:
宝　鑫（1992－　），男，博士，助理研究员，baox@tsinghua.edu.cn

中图分类号: O345; TU311
计量
- 文章访问数: 455
- HTML全文浏览量: 253
- PDF下载量: 42
- 被引次数: 16
出版历程
- 收稿日期: 2021-05-17
- 修回日期: 2021-07-17
- 网络出版日期: 2022-03-02
- 刊出日期: 2022-04-07

Stability conditions of explicit algorithms when using viscoelastic artificial boundaries

LIU Jingbo^1
,,
BAO Xin^{1
, ,},
LI Shutao^{1, 2},
WANG Fei^{1, 3}

1.
Department of Civil Engineering, Tsinghua University, Beijing 100084, China
2.
Institute of Defense Engineering, Academy of Military Science, Beijing 100036, China
3.
College of Defense Engineering, Army Engineering University of PLA, Nanjing 210007, Jiangsu, China

摘要

摘要: 黏弹性人工边界是处理无限域波动问题常用的数值模拟方法。采用显式时域逐步积分算法进行计算时，受黏弹性人工边界的阻尼、刚度等影响，人工边界区的稳定性比内部计算域的更严格，尚无明确、实用的稳定性判别准则，这限制了黏弹性人工边界在显式动力分析中的应用。针对二维黏弹性人工边界，利用基于局部子系统的稳定性分析方法和基于传递矩阵谱半径的稳定性判别准则，给出了可代表整体模型局部特征的不同边界子系统的稳定性条件解析解。通过对比分析不同计算区域的稳定性条件及其影响因素，证明了整体模型的稳定性由角点子系统控制。在此基础上，获得了含黏弹性人工边界的整体模型在显示动力计算中的统一稳定性判别准则和简化实用计算方法。在实际应用中，令积分时间步长满足稳定性条件，即可顺利完成整体模型的动力计算。以上研究可为将黏弹性人工边界应用于显式动力计算时积分时间步长的合理选取提供参考。
- 显式时域逐步积分 /
- 黏弹性人工边界 /
- 数值稳定性 /
- 局部子系统 /
- 传递矩阵
Abstract: Viscoelastic artificial boundary is a commonly used numerical simulation method to deal with the wave propagation problems in an infinite domain. When the explicit time-domain stepwise integration algorithm is adopted for such numerical analysis, the stability conditions of the artificial boundary area are more stringent than those of the internal domain due to the influence of the damping and stiffness of the viscoelastic artificial boundary. However, there is currently no clear and practical stability criterion for this problem, which affects the reasonable selection of the integral time step when using the viscoelastic artificial boundaries, and further restricts the application of viscoelastic artificial boundary in the explicit dynamic analysis. Aiming at the two-dimensional (2D) viscoelastic artificial boundary, two typical types of boundary subsystem that can represent the localized characteristics of the overall numerical model, namely the edge boundary subsystem and the corner boundary subsystem, were established and their motion equations as well as the transfer matrixes were obtained according to the stability analysis method based on the local subsystem. Then through the stability criteria based on the spectral radius of the transfer matrix, the analytical solutions of the stability conditions of different local subsystems were derived. Through the comparative analysis of the stability conditions of different calculation areas and their influencing factors, it is found that the stability of the overall model is controlled by the corner boundary subsystem. On that basis, a uniform stability criterion and a simplified practical calculation method of the stability condition for the overall model with 2D viscoelastic artificial boundary in explicit dynamic calculations were proposed. In practical applications, the dynamic calculation of the overall model can be successfully completed once the integral time step meets the proposed stability condition of the numerical system. This study provides theoretical guidance for the reasonable selection of the integral time step when applying 2D viscoelastic artificial boundaries in explicit dynamic calculations.
- explicit time domain integration /
- viscoelastic artificial boundary /
- numerical stability /
- local subsystem /
- transfer matrix

HTML全文

无限介质的波动辐射效应是波动数值模拟中需要考虑的关键问题，人工边界技术是处理此类问题最常用的技术手段。该方法在无限或半无限介质中截取有限的近场计算域，并在截断边界处施加人工边界条件，以吸收计算域内产生的外行波动。由于采用不同的数学、物理或力学原理，常用的人工边界技术可分为透射边界^[1-2]、黏性边界^[3]、黏弹性边界^[4-6]、边界元^[7-8]和完美匹配层^[9-10]等。其中，在黏弹性人工边界技术中，将介质中单侧波动的偏微分方程转化为施加于截断边界上的应力边界条件，并且等效为空间解耦的力学系统，物理意义清晰、实用性强，且具有较好的模拟精度和良好的鲁棒性。近年来，已被研究人员集成于Marc^[11]、LS-DYNA^[12]、ANSYS^[13]、ADINA^[14]、Nastran^[15]、ABAQUS^[16]等通用有限元软件，并应用于大坝^[17]、桥梁^[18]、核工程^[19]、隧道^[20]和地铁车站^[21]等建构筑物与基础的动力相互作用分析及工程场地^[22]的动力响应计算，取得了合理准确的计算结果。

采用显式时域逐步积分算法对含有黏弹性人工边界的整体模型进行计算分析时，受黏弹性人工边界阻尼、刚度等的影响，人工边界区的稳定性需比内部计算域的稳定性条件更严格。目前，尚无明确、实用的考虑黏弹性人工边界影响时显式时域逐步积分算法稳定性的判别准则，难以确定合理的数值积分时间步长，这一定程度限制了黏弹性人工边界在显式动力分析中的应用。与隐式算法相比，显式动力计算方法不需要求解耦联方程组，计算工作量较小、计算效率较高，在大规模复杂系统的动力分析尤其爆炸问题数值模拟中应用广泛。鉴于此，有必要对使用黏弹性人工边界时显式时域逐步积分算法的稳定性开展研究工作，以促进黏弹性人工边界在大规模显式动力分析中的应用。

本文中，针对含黏弹性人工边界的数值模型在显式算法中的稳定性问题，利用基于局部子系统的数值稳定性分析方法，给出可代表整体模型局部特征的不同边界子系统的稳定性条件解析解，比较分析不同计算区域的稳定性条件及其影响因素，给出整体模型在显式动力计算中的统一稳定性判别准则和简化实用计算方法。

1. 黏弹性人工边界模型与参数设置

黏弹性人工边界由在截断边界处设置的空间解耦的弹簧-阻尼器构成，如图1所示。

图 1 黏弹性人工边界

Figure 1. Schematic diagram of viscoelastic artificial boundaries

下载: 全尺寸图片幻灯片

二维情况下，黏弹性人工边界等效物理系统的弹簧刚度K_B和阻尼系数C_B分别为^[23]：

$\begin{gathered} {K_{{\text{BT}}}} = {\alpha _{\text{T}}}\frac{G}{R}\sum {L_i},\qquad {C_{{\text{BT}}}} = \rho {c_{\text{S}}}\sum {L_i} \hfill \\ {K_{{\text{BN}}}} = {\alpha _{\text{N}}}\frac{G}{R}\sum {L_i},\qquad {C_{{\text{BN}}}} = \rho {c_{\text{P}}}\sum {L_i} \hfill \\ \end{gathered}$

(1)

式中：G和ρ分别为内部介质的剪切模量和质量密度，c_S和c_P分别为内部介质的剪切波与压缩波波速， $\sum {L_i}$ 为与黏弹性人工边界相连的有限元节点所代表的边界长度，L为有限单元的长度（对于侧边节点， $\sum {L_i}$ =L，对于角点， $\sum {L_i}$ =L/2）；R为波源至人工边界节点的距离，由于土-结构动力相互作用问题中的波源通常非点波源，因此R一般取为平均意义上的常值；下标T和N分别代表切向和法向；α_T和α_N分别为切向和法向的黏弹性人工边界参数，文献[23]中推荐的数据见表1。

表 1 二维黏弹性人工边界参数的数据^[23]

Table 1. The values of two-dimensional viscoelastic artificial boundary coefficients^[23]

参数	范围	建议
α_T	0.35～0.65	0.5
α_N	0.80～1.20	1.0

下载: 导出CSV

| 显示表格

2. 采用黏弹性人工边界时显式算法的稳定性分析方法

2.1 基于局部子系统的数值积分稳定性分析方法

数值稳定性是制约显式动力计算时间步长的主要因素，以往针对含人工边界条件的计算系统的稳定性研究，多基于整体模型或由若干排（列）节点构成的复杂节点系统^[24-25]，通过数值方法判断其稳定性。因未能给出解析形式的人工边界稳定性判别准则，此类方法的实用性仍有所欠缺。李述涛等^[26] 总结和归纳了现有离散模型数值积分稳定性研究工作，得出以下结论：数值积分算法的稳定性由计算模型的最高阶频率即系统的截止频率控制；截止频率对应的振型一般呈现局部节点系相邻节点交错振动的模态，二维平面应变情况下，其振动形式如图2所示。

图 2 基于局部子系统的稳定性分析

Figure 2. Stability analysis based on local subsystems

下载: 全尺寸图片幻灯片

在上述基础上，李述涛等^[26]提出了一种利用局部子系统估算整体有限元模型数值稳定性的方法。根据该方法，对于任意整体有限元模型，如能较准确地判断其最高阶振型，则可利用振型节点的空间分布规律，截取由相邻振型节点包围的最小局部子系统（见图2中的子系统a）。对该子系统的边界节点施加与整体模型最高阶振型一致的约束条件，并进行稳定性分析，获得的局部子系统稳定性条件即为整体模型的稳定性条件。而对于不能准确判断整体模型最高阶振型的情况，可选取由相邻单元构成的最小子系统（见图2中的子系统b）。对该子系统的全部边界节点施加固定约束，分析所有b型局部子系统的稳定性条件，其下限值为整体有限元模型稳定性条件的上限逼近。

基于局部子系统的稳定性分析方法，提供了一种从理论上估算含人工边界的复杂离散系统数值稳定性的技术手段。对于整体模型的内部计算域，可以采用子系统a直接获稳定性条件，而对于边界区域，由于人工边界的影响，难以准确判断最高阶振型，无法从整体计算模型中分割出子系统a，此时可以采用b型子系统进行稳定性分析。

2.2 基于传递矩阵谱半径的稳定性判别方法

稳定性条件因所采用的数值积分方法而异，为解决实际工程应用中黏弹性人工边界的稳定性问题，以在通用有限元软件中应用较广泛的蛙跳格式中心差分时域逐步积分算法为基础^[27]，分析含黏弹性人工边界的整体计算模型的数值稳定性。对于其他显式时域逐步积分算法，采用相似的步骤即可得到对应的稳定性条件。

根据所采用的时域逐步积分算法，二维有限元模型节点运动量的时域递推公式为：

$\left[ {\begin{array}{*{20}{c}} {{{\dot u}_x}\left( {t + {{\Delta t}/ 2}} \right)} \\ {{{\dot u}_y}\left( {t + {{\Delta t} / 2}} \right)} \end{array}} \right]= \left[ {\begin{array}{*{20}{c}} {{{\dot u}_x}\left( {t - {{\Delta t} / 2}} \right)} \\ {{{\dot u}_y}\left( {t - {{\Delta t} / 2}} \right)} \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} {{{\ddot u}_x}\left( t \right)} \\ {{{\ddot u}_y}\left( t \right)} \end{array}} \right]\Delta t,\qquad \left[ {\begin{array}{*{20}{c}} {{u_x}\left( {t + \Delta t} \right)} \\ {{u_y}\left( {t + \Delta t} \right)} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{u_x}\left( t \right)} \\ {{u_y}\left( t \right)} \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} {{{\dot u}_x}\left( {t + {{\Delta t} / 2}} \right)} \\ {{{\dot u}_y}\left( {t + {{\Delta t} /2}} \right)} \end{array}} \right]\Delta t$

(2)

式中：u、 $\dot u$ 和 $\ddot u$ 分别为位移、速度和加速度，t为时间， ${{\Delta }}t$ 为积分时间步长，x和y的方向见图2。

对离散的动力系统进行稳定性分析时，可根据显式时域逐步积分格式，运动方程为：

${\boldsymbol{U}}\left( {t + \Delta t} \right) = {\boldsymbol{A}} \cdot {\boldsymbol{U}}\left( t \right) + {\boldsymbol{P}}\left( t \right)$

(3)

式中：A为传递矩阵，P为外力向量，向量U由系统各节点的位移、速度和加速度组成。本文中，采用显式时域逐步积分算法：

${\boldsymbol{U}}\left( t \right) = \left[ {\begin{array}{*{20}{c}} {{u_x}\left( t \right)} \\ {{u_y}\left( t \right)} \\ {{{\dot u}_x}\left( {t - {{\Delta t} / 2}} \right)\Delta t} \\ {{{\dot u}_y}\left( {t - {{\Delta t}/ 2}} \right)\Delta t} \end{array}} \right]$

(4)

数值积分方法的稳定性仅与积分格式、时空离散步长和计算系统的力学参数有关，而与外力向量P无关。如满足以下条件，则积分格式是稳定的^[28]：(1) ρ(A) ≤ 1，其中ρ(A)为传递矩阵A的谱半径，即ρ(A) =max|λ_i|，λ_i为传递矩阵A的第i个特征值，对于本文分析的二维情况，i=1～4；(2) 如A具有多重特征值，则该特征值的模小于1。

结合以上理论与方法，可将黏弹性人工边界在显式时域逐步积分算法中的稳定性分析分解为3个步骤：（1）选取能体现整体有限元模型不同局部特征的子系统；（2）建立各子系统的运动方程，结合时域逐步积分算法推导传递矩阵，并求解谱半径；（3）根据谱半径的模小于等于1的限制条件求解各子系统的稳定性条件，比较后选取最严格的稳定性条件作为整体模型最大稳定积分时间步长的上限估计。

另外，在满足了离散模型动力计算的稳定性条件下，为进一步满足精度要求，离散化网格的尺寸 $\Delta x$ 需满足以下条件^[29]：

$\Delta x \text{≤} \left( {\frac{1}{6}\sim \frac{1}{8}} \right){\lambda _{\min }},\qquad {\lambda _{\min }} = \frac{{{c_{\min }}}}{{{f_{\max }}}}$

(5)

式中： $\lambda_{\min }$ 为离散网格中波动传播的最短波长，c_min为介质中的最小波速，f_max为波动问题数值模拟的截止频率。

3. 采用黏弹性人工边界时显式算法的稳定性分析

下面，针对二维情况，分析采用黏弹性人工边界时显式时域逐步积分算法的稳定性。由于稳定性条件与计算系统的单元尺寸和物理特性密切相关，为不失一般性，在数值算法的稳定性分析中通常采用规则几何模型、均一介质和均匀离散网格进行建模分析。对于不规则网格或非均匀介质，可根据有限元网格的最小尺寸及不同区域的介质材料参数，采用本文方法分别计算稳定性条件，并选取其下限值作为最大积分时间步长的选取依据。

二维半无限空间近场有限元模型如图3(a)所示，不考虑介质阻尼，内部均匀介质在中心差分算法中的稳定性条件为^[26]：

图 3 二维半无限空间近场有限元模型及二维平面应变单元

Figure 3. Two-dimensional semi-infinite near-field finite element model and two-dimensional plane strain element

下载: 全尺寸图片幻灯片

$\Delta t \text{≤}\frac{L}{{{c_{\text{P}}}}}$

(6)

当采用黏弹性人工边界时，因难以确定人工边界区的最高阶振型，应采用b型子系统进行稳定性分析。此时，根据截取位置的不同，人工边界区的b型局部子系统可分为侧边（底边）子系统和角点子系统，如图3(a)所示。

为建立局部子系统的运动方程，先给出规则二维平面应变等参数单元（见图3(b)）的刚度矩阵^[30]：

${{\boldsymbol{K}}^{\text{e}}} = \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{K}}_{11}}}&{{{\boldsymbol{K}}_{12}}}&{{{\boldsymbol{K}}_{13}}}&{{{\boldsymbol{K}}_{14}}} \\ {{{\boldsymbol{K}}_{21}}}&{{{\boldsymbol{K}}_{22}}}&{{{\boldsymbol{K}}_{23}}}&{{{\boldsymbol{K}}_{24}}} \\ {{{\boldsymbol{K}}_{31}}}&{{{\boldsymbol{K}}_{32}}}&{{{\boldsymbol{K}}_{33}}}&{{{\boldsymbol{K}}_{34}}} \\ {{{\boldsymbol{K}}_{41}}}&{{{\boldsymbol{K}}_{42}}}&{{{\boldsymbol{K}}_{43}}}&{{{\boldsymbol{K}}_{44}}} \end{array}} \right] ,\qquad {{\boldsymbol{K}}_{ij}} = \frac{{E\left( {1 - \mu } \right)}}{{4\left( {1 + \mu } \right)\left( {1 - 2\mu } \right)}}\left[ {\begin{array}{*{20}{c}} {{K_1}}&{{K_3}} \\ {{K_4}}&{{K_2}} \end{array}} \right]{\text{ }}\qquad {i,j = 1,2,3,4}$

(7)

$\begin{array}{c} \displaystyle{K_1} = {\xi _i}{\xi _j}\left( {1 + \frac{1}{3}{\eta _i}{\eta _j}} \right) + \frac{{1 - 2\mu }}{{2\left( {1 - \mu } \right)}}{\eta _i}{\eta _j}\left( {1 + \frac{1}{3}{\xi _i}{\xi _j}} \right),\qquad {K_2} = {\eta _i}{\eta _j}\left( {1 + \frac{1}{3}{\xi _i}{\xi _j}} \right) + \frac{{1 - 2\mu }}{{2\left( {1 - \mu } \right)}}{\xi _i}{\xi _j}\left( {1 + \frac{1}{3}{\eta _i}{\eta _j}} \right), \\ \displaystyle {{K_3} = \frac{\mu }{{1 - \mu }}{\xi _i}{\eta _j} + \frac{{1 - 2\mu }}{{2\left( {1 - \mu } \right)}}{\eta _i}{\xi _j},\qquad {K_4} = \frac{\mu }{{1 - \mu }}{\eta _i}{\xi _j} + \frac{{1 - 2\mu }}{{2\left( {1 - \mu } \right)}}{\xi _i}{\eta _j}} \end{array}$

式中：E为杨氏模量，μ为泊松比，ξ和η为局部坐标，ξ_i,j=±1，η_i,j=±1，下标i和j为节点。二维平面应变单元的质量阵的集中质量模型为^[30]：

${{\boldsymbol{M}}^{\text{e}}}{\text{ = }}\frac{1}{4}\rho {L^2}{\boldsymbol{E}}$

(8)

式中：E为8×8型的单位矩阵。

3.1 侧边子系统的稳定性条件

在图3(a)所示的侧边子系统中，仅节点1为自由节点，其余节点均施加了固定约束。利用黏弹性人工边界物理参数（式(1)）及平面应变单元的刚度和质量矩阵（式(7)～(8)），进行矩阵组装，得到侧边子系统中节点1的运动方程：

$\left[ {\begin{array}{*{20}{c}} {{m_1}}&0 \\ 0&{{m_1}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{{\ddot u}_x}} \\ {{{\ddot u}_y}} \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} {{c_1}}&0 \\ 0&{{c_2}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{{\dot u}_x}} \\ {{{\dot u}_y}} \end{array}} \right]{\text{ + }}\left[ {\begin{array}{*{20}{c}} {{k_1}}&0 \\ 0&{{k_2}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{u_x}} \\ {{u_y}} \end{array}} \right] = \left[{\begin{array}{*{20}{c}} 0 \\ 0 \end{array}} \right]$

(9)

${m}_{1}=\frac{\rho {L}^{2}}{2},\text{ }{c}_{1}=\rho {c}_{\text{P}}L,\text{ }{c}_{2}=\rho {c}_{\text{S}}L,\text{ }{k}_{1}=\rho {c}_{\text{S}}^{2}\left[\frac{2(3-4\mu)}{3(1-2\mu )}\text+\frac{{\alpha }_{\text{N}}L}{R}\right],\text{ }{k}_{2}=\rho {c}_{\text{S}}^{2}\left[\frac{2(3-4\mu)}{3(1-2\mu )}\text+\frac{{\alpha }_{\text{T}}L}{R}\right]$

(10)

联立式(2)、(9)，得到侧边子系统的传递矩阵：

$\left[ {\begin{array}{*{20}{c}} {{u_x}\left( {t + \Delta t} \right)} \\ {{u_y}\left( {t + \Delta t} \right)} \\ {{{\dot u}_x}\left( {t + {{\Delta t} / 2}} \right)\Delta t} \\ {{{\dot u}_y}\left( {t + {{\Delta t} / 2}} \right)\Delta t} \end{array}} \right] = {\boldsymbol{A}}\left[ {\begin{array}{*{20}{c}} {{u_x}\left( t \right)} \\ {{u_y}\left( t \right)} \\ {{{\dot u}_x}\left( {t - {{\Delta t} /2}} \right)\Delta t} \\ {{{\dot u}_y}\left( {t - {{\Delta t}/ 2}} \right)\Delta t} \end{array}} \right],\quad {\boldsymbol{A}} = \frac{1}{{{m_1}}}\left[ {\begin{array}{*{20}{c}} {{m_1} - {k_1}{{(\Delta t)}^2}}&0&{{m_1} - {c_1}\Delta t}&0 \\ 0&{{m_1} - {k_1}{{(\Delta t)}^2}}&0&{{m_1} - {c_1}\Delta t} \\ { - {k_2}{{(\Delta t)}^2}}&0&{{m_1} - {c_2}\Delta t}&0 \\ 0&{ - {k_2}{{(\Delta t)}^2}}&0&{{m_1} - {c_2}\Delta t} \end{array}} \right]$

(11)

求解传递矩阵A的特征值，得到：

$\begin{aligned} & {\lambda _{1,2}} = 1 - \frac{{{k_1}{{(\Delta t)}^2}}}{{2{m_1}}} - \frac{{{c_1}\Delta t}}{{2{m_1}}} \pm \frac{{\sqrt {k_1^2{{(\Delta t)}^4} - 4{k_1}{m_1}{{(\Delta t)}^2} + 2{k_1}{c_1}{{(\Delta t)}^3} + {c_1^2}{{(\Delta t)}^2}} }}{{2{m_1}}} \\ & {\lambda _{3,4}} = 1 - \frac{{{k_2}{{(\Delta t)}^2}}}{{2{m_1}}} - \frac{{{c_2}\Delta t}}{{2{m_1}}} \pm \frac{{\sqrt {k_2^2{{(\Delta t)}^4} - 4{k_2}{m_1}{{(\Delta t)}^2} + 2{k_2}{c_2}{{(\Delta t)}^3} + {c_2^2}{{(\Delta t)}^2}} }}{{2{m_1}}} \end{aligned}$

(12)

稳定性判别条件 $\rho \left( {\boldsymbol{A}} \right){\text{=}}\max \left| {{\lambda _i}} \right| \text{≤} 1$ 可等效为传递矩阵A的任意特征值的模均小于等于1，即：

$\left\{\Delta t|\left|{\lambda }_{i}\left(\Delta t\right)\right|\text{≤} 1\qquad i=1～4,\text{ }\Delta t＞0\right\}$

(13)

直接求解式(13)较困难，可先将式(13)中的不等式化为等式：

$f\left( {{{\Delta }}t} \right) = \left| {{\lambda _i}\left( {{{\Delta }}t} \right)} \right| - 1 = 0\qquad i = 1\sim 4$

(14)

求其在 ${{\Delta }}t$ 正半轴的零点，得到：

$\Delta t = \left\{ {\begin{array}{*{20}{c}} {0,\;\dfrac{{4{m_1}}}{{\sqrt {4{k_1}{m_1} + c_1^2} + {c_1}}}}&\quad{i = 1,2} \\ {0,\;\dfrac{{4{m_1}}}{{\sqrt {4{k_2}{m_1} + c_2^2} + {c_2}}}}&\quad{i = 3,4} \end{array}} \right.$

(15)

对于i=1,2或i=3,4的情况，式(15)所示的零点分别将 ${{\Delta }}t$ 的正半轴划分为两个区间，检验当 ${{\Delta }}t$ 处于不同区间时函数f ( ${{\Delta }}t$ )的正负，再利用函数的连续性判断满足稳定性条件（式(13)）的 ${{\Delta }}t$ 范围。

令 ${{\Delta }}t$ 为大于0的小量，忽略式(12)中的高阶小量 ${{o(\Delta }}t{{)}}$ ，得到：

${\lambda _{1,2}} = 1 - \frac{{{c_1}\Delta t}}{{2{m_1}}} \pm \frac{{\sqrt {c_1^2 - 4{k_1}{m_1}} \Delta t}}{{2{m_1}}},\qquad{\text{ }}{\lambda _{3,4}} = 1 - \frac{{{c_2}\Delta t}}{{2{m_1}}} \pm \frac{{\sqrt {c_2^2 - 4{k_2}{m_1}} \Delta t}}{{2{m_1}}}$

(16)

将式(10)代入式(16)，得到：

$c_1^2 - 4{m_1}{k_1} = - \left( {\frac{{6 - 10\mu }}{{3 - 6\mu }} + \frac{{2{\alpha _{\text{N}}}L}}{R}} \right){\rho ^2}c_{\text{S}}^2{L^2},\qquad{c_2^2} - 4{m_1}{k_2} = - \left( {\frac{{9 - 10\mu }}{{3 - 6\mu }} + \frac{{2{\alpha _{\text{T}}}L}}{R}} \right){\rho ^2}c_{\text{S}}^2{L^2}$

(17)

由于 $0 \text{≤} \mu \text{≤} 0.5$ ，因此式(17)中的两式均为负值，即 ${\lambda _{1,2}}$ 和 ${\lambda _{3,4}}$ 均为复数，其模为：

$\left| {{\lambda _{1,2}}} \right|{{ = }}\sqrt {1 - \frac{{{c_1}{{\Delta }}t}}{{{m_1}}} + \frac{{{k_1}{{{{(}}\Delta t{{)}}}^2}}}{{{m_1}}}} \approx \sqrt {1 - \frac{{{c_1}{{\Delta }}t}}{{{m_1}}}} \text{＜} 1,\qquad\left| {{\lambda _{3,4}}} \right|{{ = }}\sqrt {1 - \frac{{{c_2}{{\Delta }}t}}{{{m_1}}} + \frac{{{k_2}{{{{(}}\Delta t{{)}}}^2}}}{{{m_1}}}} \approx \sqrt {1 - \frac{{{c_2}{{\Delta }}t}}{{{m_1}}}} \text{＜} 1$

(18)

由函数连续性可知，当 $\Delta t \in \left( {0,\dfrac{{4{m_1}}}{{\sqrt {4{k_1}{m_1} + c_1^2} + {c_1}}}} \right)$ 时， $\left| {{\lambda _{1,2}}\left( {\Delta t} \right)} \right| \text{≤} 1$ ；当 $\Delta t \in \left( {0,\dfrac{{4{m_1}}}{{\sqrt {4{k_2}{m_1} + c_2^2} + {c_2}}}} \right)$ 时， $\left| {{\lambda _{3,4}}\left( {\Delta t} \right)} \right| \text{≤} 1$ 。这两种情况分别满足稳定性条件式(13)。

再令 ${{\Delta }}t\to \infty$ ，此时 ${\lambda _{1,2}}$ 和 ${\lambda _{3,4}}$ 均为实数，且满足：

$\max \left| {{\lambda _{1,2}}} \right| \approx \frac{{{k_1}{{(\Delta t)}^2}}}{{{m_1}}} \to \infty ,\qquad {\rm max} \left| {{\lambda _{3,4}}} \right| \approx \frac{{{k_2}{{(\Delta t)}^2}}}{{{m_1}}} \to \infty$

(19)

由函数连续性可知，当 $\Delta t \in \left( {\dfrac{{4{m_1}}}{{\sqrt {4{k_1}{m_1} + c_1^2} + {c_1}}},{\text{ + }}\infty } \right)$ 时， $\left|{\lambda }_{1,2}\left(\Delta t\right)\right|＞1$ ；当 $\Delta t \in \left( {\dfrac{{4{m_1}}}{{\sqrt {4{k_2}{m_1} + c_2^2} + {c_2}}},{{ + }}\infty } \right)$ 时， $\left|{\lambda }_{3,4}\left(\Delta t\right)\right|＞1$ 。这两种情况均不满足稳定性条件式(13)。

综合上述，可得到满足式(13)的稳定性条件的 $\Delta t$ 范围为：

$\Delta t \text{≤} \left\{ {\begin{array}{*{20}{c}} {\dfrac{{4{m_1}}}{{\sqrt {4{k_1}{m_1} + c_1^2} + {c_1}}}}&\qquad{i = 1,2} \\ {\dfrac{{4{m_1}}}{{\sqrt {4{k_2}{m_1} + c_2^2} + {c_2}}}}&\qquad{i = 3,4} \end{array}} \right.$

(20)

由式(10)可知k₁＞k₂、c₁＞c₂，则式(20)中第1式的稳定性条件更严格，将式(10)代入，整理可得：

$\Delta t \text{≤} {\gamma _1}\dfrac{L}{{{c_{\text{P}}}}},\qquad{\gamma _1}{{ = }}\dfrac{{2\sqrt {\dfrac{{7{{+}}2\tilde \mu }}{3}\dfrac{{{R^2}}}{{{L^2}}}{{+}}{\alpha _{\text{N}}}\tilde \mu \dfrac{R}{L}}-\dfrac{{2R}}{L}}}{{\dfrac{{{{4+}}2\tilde \mu }}{3}\dfrac{R}{L}{{+}}{\alpha _{\text{N}}}\tilde \mu }},\qquad {\tilde\mu} {{= }}\dfrac{{1-2\mu }}{{1-\mu }}$

(21)

3.2 角点子系统的稳定性条件

在图3(a)所示的角点子系统中，除节点1外，其余节点均施加了固定约束。按节点进行矩阵组装，得到角点子系统中节点1的运动方程：

$\left[ {\begin{array}{*{20}{c}} {{m_2}}&0 \\ 0&{{m_2}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{{\ddot u}_x}} \\ {{{\ddot u}_y}} \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} c&0 \\ 0&c \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{{\dot u}_x}} \\ {{{\dot u}_y}} \end{array}} \right]{\text{ + }}\left[ {\begin{array}{*{20}{c}} {{k_3}}&{ - {k_4}} \\ { - {k_4}}&{{k_3}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{u_x}} \\ {{u_y}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} 0 \\ 0 \end{array}} \right]$

(22)

${m}_{2}=\frac{\rho {L}^{2}}{4},\qquad c=\frac{\rho \left({c}_{\text{P}}\text+{c}_{\text{S}}\right)L}{2}\text{,}\qquad {k}_{3}=\rho {c}_{\text{S}}^{2}\left[\frac{3-4\mu }{3(1-2\mu )}\text+\frac{\left({\alpha }_{\text{N}}\text+{\alpha }_{\text{T}}\right)L}{2R}\right],\qquad {k}_{4}=\frac{\rho {c}_{\text{S}}^{2}}{4(1-2\mu )}$

(23)

联立式(2)、(22)，可得到角点子系统的传递矩阵A及其特征值：

$\begin{split} & {{\boldsymbol{A}} = \dfrac{1}{{{m_2}}}\left[ {\begin{array}{*{20}{c}} {{m_2} - {k_3}{{(\Delta t)}^2}}&{{k_4}{{(\Delta t)}^2}}&{{m_2} - c\Delta t}&0 \\ {{k_4}{{(\Delta t)}^2}}&{{m_2} - {k_3}{{(\Delta t)}^2}}&0&{{m_2} - c\Delta t} \\ { - {k_3}{{(\Delta t)}^2}}&{{k_4}{{(\Delta t)}^2}}&{{m_2} - c\Delta t}&0 \\ {{k_4}{{(\Delta t)}^2}}&{ - {k_3}{{(\Delta t)}^2}}&0&{{m_2} - c\Delta t} \end{array}} \right]} \\ & {{\lambda _{1,2}} = 1 - \dfrac{{{k_3}{{(\Delta t)}^2}}}{{2{m_2}}} - \dfrac{{c\Delta t}}{{2{m_2}}} + \dfrac{{{k_4}{{(\Delta t)}^2}}}{{2{m_2}}} \pm \dfrac{{\sqrt {{{\left( {{k_3}{{(\Delta t)}^2} - {k_4}{{(\Delta t)}^2} + c\Delta t - 2{m_2}} \right)}^2} - 4\left( {m_2^2 - c{m_2}\Delta t} \right)} }}{{2{m_2}}}} \\ & {{\lambda _{3,4}} = 1 - \dfrac{{{k_3}{{(\Delta t)}^2}}}{{2{m_2}}} - \dfrac{{c\Delta t}}{{2{m_2}}} - \dfrac{{{k_4}{{(\Delta t)}^2}}}{{2{m_2}}} \pm \dfrac{{\sqrt {{{\left( {{k_3}{{(\Delta t)}^2} + {k_4}{{(\Delta t)}^2} + c\Delta t - 2{m_2}} \right)}^2} - 4\left( {m_2^2 - c{m_2}\Delta t} \right)} }}{{2{m_2}}}} \\[-100pt] \end{split}$

(24)

采用与第3.1节相同的方法，计算方程 $\left| {{\lambda _i}\left( {{{\Delta }}t} \right)} \right| - 1 = 0\left( {i = 1\sim 4} \right)$ 在 ${{\Delta }}t$ 正半轴的零点，并基于方程的连续性和零点位置，判断满足稳定性条件式(13)的 ${{\Delta }}t$ 所处的区间：

$\Delta t \text{≤} \frac{{4{m_2}}}{{\sqrt {4({k_3} + {k_4}){m_2} + {c^2}} + c}}$

(25)

将式(23)代入式(25)，整理得到：

$\Delta t \text{≤} {\gamma _2}\dfrac{L}{{{c_{\text{P}}}}},\qquad {\gamma _2}{{=}}\dfrac{{4\sqrt {3\left( {10 + 2\tilde \mu + 3\sqrt {2\tilde \mu } } \right)\dfrac{{{R^2}}}{{{L^2}}} + 9\left( {{\alpha_{\text{N}}}{{ + }}{\alpha_{\text{T}}}} \right)\tilde \mu \dfrac{R}{L}-6 \left( {2 +\sqrt {2\tilde \mu } } \right)\dfrac{R}{L}}}}{{\left( {14 + \tilde \mu } \right)\dfrac{R}{L} + 6\left( {{\alpha_{\text{N}}}{{ + }}{\alpha_{\text{T}}}} \right)\tilde \mu }},\qquad \tilde \mu {{=}}\dfrac{{1 - 2\mu }}{{1 - \mu }}$

(26)

3.3 稳定性条件的比较

由于整体模型的稳定性受不同局部区域中最严格的稳定性条件控制，根据式(6)、(21)和(26)，可将采用黏弹性人工边界时显式算法的稳定性条件统一改写为：

$\dfrac{{c}_{\text{P}}\text{Δ}t}{L} \text{≤} \gamma ,\quad \gamma ={\min}\left\{\begin{array}{*{20}{l}} 1&内部计算域\\ {\gamma }_{1}&侧边子系统\\ {\gamma }_{2}&角点子系统\end{array}\right.$

(27)

侧边和角点子系统的稳定性系数γ₁和γ₂均为无量纲系数，仅与内部介质的泊松比μ、波源距与有限单元长度比R/L有关。不同区域的稳定性系数γ随μ和R/L的变化情况，如图4所示；当R/L=1,5,+∞时稳定性系数γ随泊松比μ的变化，当泊松比μ=0.2,0.3,0.4时稳定性系数γ随R/L的变化，如图5～6所示。

图 4 稳定性条件的比较

Figure 4. Comparison of stability conditions

下载: 全尺寸图片幻灯片

图 5 不同R/L时稳定性系数γ随泊松比μ的变化

Figure 5. Variations of stability coefficient γ with Poisson’s ratio μ under different R/L

下载: 全尺寸图片幻灯片

图 6 不同泊松比μ时稳定性系数γ随R/L的变化

Figure 6. Variations of stability coefficient γ with R/L under different Poisson’s ratio μ

下载: 全尺寸图片幻灯片

由图4～6可见，侧边和角点子系统的数值积分稳定性均随泊松比μ的增大变得宽松，且随着波源距与单元尺寸比R/L的增大先略为放松，当R/L>5时基本保持不变。此外，在不同情况下，角点子系统的稳定性系数均最小，对整体计算模型的稳定性起控制作用，其稳定性条件即为采用黏弹性人工边界时显式时域逐步积分算法的统一稳定性条件。可表示为：

$\Delta t \text{≤} \dfrac{{4\sqrt{3\left( {10 + 2\tilde \mu + 3\sqrt{2\tilde \mu } } \right)\dfrac{{{R^2}}}{{{L^2}}} + 9\left( {{\alpha_{\rm{N}}}{{ + }}{\alpha_{\rm{T}}}} \right)\tilde \mu \dfrac{R}{L}} - 6\left( {2 + \sqrt {2\tilde \mu } } \right)\dfrac{R}{L}}}{{\left( {14 + \tilde \mu } \right)\dfrac{R}{L} + 6\left( {{\alpha_{\rm{N}}}{{ + }}{\alpha_{\rm{T}}}} \right)\tilde \mu }}\dfrac{L}{{{c_{\rm{P}}}}}，\qquad\tilde \mu {{=}}\dfrac{{1 - 2\mu }}{{1 - \mu }}$

(28)

在采用黏弹性人工边界和中心差分法进行动力分析时，根据黏弹性人工边界参数α_N、α_T和波源距与单元尺寸之比R/L及泊松比μ，即可由式(28)确定稳定的时间积分步长。在实际应用中，为了更方便地确定稳定时间积分步长 ${{\Delta }}t$ ，可采纳几种常见情况的稳定性系数，见表2。

表 2 建议的几种常见情况的稳定性系数

Table 2. Recommended stability coefficients for several common cases

R/L				γ
R/L	μ=0.10	μ=0.15	μ=0.20	μ=0.25	μ=0.30	μ=0.35	μ=0.40
1	0.47	0.48	0.49	0.50	0.51	0.53	0.55
5	0.50	0.50	0.51	0.52	0.53	0.55	0.57
10	0.50	0.51	0.51	0.52	0.53	0.55	0.57
20	0.50	0.51	0.52	0.53	0.54	0.55	0.57
50	0.50	0.51	0.52	0.53	0.54	0.55	0.57
+∞	0.51	0.51	0.52	0.53	0.54	0.55	0.57

下载: 导出CSV

| 显示表格

在实际工程中，通常满足材料泊松比 $\mu \text{≥} 0.1$ 、波源距与单元尺寸比 ${R \mathord{\left/ {\vphantom {R L}} \right. } L} \text{≥} 5$ ，且稳定性条件随泊松比μ和波源距与单元尺寸比R/L的增大变得宽松。在实际应用中，为避免复杂的公式计算，可直接采用 $\mu=0.1$ 、 $R / L=5$ 对应的稳定性条件作为系统最大稳定积分时间步长的保守估计：

${{\Delta }}t \text{≤} \frac{1}{2}\frac{L}{{{c_{\text{P}}}}}$

(29)

4. 稳定性条件验证

4.1 均匀半空间模型

为验证以上稳定性分析的准确性，采用大型通用有限元软件ABAQUS建立均匀半空间模型，如图7所示。模型尺寸为100 m×50 m，内部介质的密度为2 000 kg/m³，剪切波速为200 m/s，泊松比为0.3。采用四节点平面应变单元进行有限元离散，网格尺寸为1 m×1 m。在模型的截断边界处施加黏弹性人工边界，并令人工边界参数α_T=0.5、α_N=1.0。在模型中点（点A）处施加持时0.2 s、幅值为1 MN的脉冲荷载，如图8所示。

图 7 均匀半空间模型

Figure 7. The homogeneous half-space model

下载: 全尺寸图片幻灯片

图 8 脉冲荷载

Figure 8. The impulse load

下载: 全尺寸图片幻灯片

根据以上模型参数，分别采用式(6)、(21)和(26)，给出的稳定性条件及式(29)给出的稳定性系数建议值，计算的模型不同区域的最大稳定时间步长见表3。

表 3 均匀半空间模型的稳定性系数和最大稳定时间步长

Table 3. Stability coefficients and maximum stable time steps of the homogeneous model

模型分区	稳定性系数γ		最大稳定时间步长Δt/s
模型分区	计算	建议	计算	建议
内部区域	1.00	0.50	0.002 7	0.001 35
侧边子系统	0.74	0.50	0.002 0	0.001 35
角点子系统	0.52	0.50	0.001 4	0.001 35

下载: 导出CSV

| 显示表格

分别采用不同的固定时间步长进行显式动力计算，并统计其稳定性状态，结果见表4。提取失稳时计算模型的位移云图，如图9所示，并比较采用不同固定积分时间步长时模型底部角点（点B）的竖向位移，如图10所示。

表 4 不同固定时间步长时均匀半空间模型的稳定性状态

Table 4. The stability states of the homogeneous model under different fixed time steps

时间步长 ${{\Delta }}t$ /s	稳定性系数γ	稳定性状态
0.003 0	1.11	内部首先失稳
0.002 7	1.00	侧边首先失稳
0.002 0	0.74	角点首先失稳
0.001 4	0.52	稳定计算
0.001 35	0.50	稳定计算

下载: 导出CSV

| 显示表格

图 9 均匀半空间模型的位移分布

Figure 9. Displacement distributions of the homogeneous half-space model

下载: 全尺寸图片幻灯片

图 10 均匀半空间模型底部角点的竖向位移

Figure 10. Vertical displacements of the corner point in the homogenous half-space model

下载: 全尺寸图片幻灯片

由表3～4和图9～10，可以得到以下结论：本文中的采用黏弹性人工边界时显式时域逐步积分算法的稳定性条件解析解可准确判断数值计算中的稳定性状态；当所积分时间步长不满足内部计算域、侧边子系统或角点子系统的稳定性条件时，相应区域分别发生失稳，失稳状态如图9所示；当积分时间步长满足角点子系统的稳定性条件时，可顺利完成整体模型的数值计算，即采用黏弹性人工边界时显式时域逐步积分算法的稳定性条件由角点区控制；由于稳定性系数的建议值比理论值更保守，采用本文中的稳定性系数建议值确定的积分时间步长可保证显式动力计算的顺利完成。

4.2 成层半空间模型

为进一步验证以上结论在复杂场地中的适用性，建立成层半空间模型。如图11所示，模型尺寸为100 m×100 m，分为上下两层，层厚均为50 m，上层介质的密度、剪切波速和泊松比分别为1 800 kg/m³、150 m/s和0.28，下层介质的密度、剪切波速和泊松比分别为2 200 kg/m³、250 m/s和0.35。采用四节点平面应变单元对计算模型进行有限元离散，网格尺寸为1 m×1 m。在模型截断边界处施加黏弹性人工边界，在模型中点（点A）处施加脉冲荷载（见图8）。

图 11 成层半空间模型

Figure 11. The layered half-space model

下载: 全尺寸图片幻灯片

由以上参数计算得到的成层半空间模型不同区域的稳定性系数和最大稳定时间步长见表5，上层介质及其侧边子系统的稳定性条件比下层介质的稳定性条件更宽松，可以判断该模型的数值积分稳定性应由下层介质控制。分别采用不同的固定时间步长进行显式动力计算，并统计其稳定性状态，结果见表6。提取失稳时计算模型的位移云图，如图12所示，并比较采用不同固定时间步长时模型底部角点（点B）的竖向位移，如图13所示。

表 5 成层半空间模型的稳定性系数与最大稳定时间步长

Table 5. Stability coefficients and maximum stable time steps of the layered model

介质	模型分区	稳定性系数γ		最大稳定时间步长Δt/s
介质	模型分区	计算	建议	计算	建议
上层	内部区域	1.00	0.5	0.003 7	0.000 95
上层	侧边子系统	0.76	0.5	0.002 8	0.000 95
下层	内部区域	1.00	0.5	0.001 9	0.000 95
	侧边子系统	0.79	0.5	0.001 5	0.000 95
	角点子系统	0.59	0.5	0.001 1	0.000 95

下载: 导出CSV

| 显示表格

表 6 采用不同固定时间步长时成层半空间模型的稳定性状态

Table 6. The stability state of the layered model under different fixed time steps

时间步长 ${{\Delta }}t$ /s	稳定性系数γ	稳定性状态
0.002 0	1.05	内部首先失稳
0.001 9	1.00	侧边首先失稳
0.001 5	0.79	角点首先失稳
0.001 0	0.53	稳定计算
0.000 95	0.50	稳定计算

下载: 导出CSV

| 显示表格

图 12 成层半空间模型的位移分布

Figure 12. Displacement distributions of the layered half-space model

下载: 全尺寸图片幻灯片

图 13 成层半空间模型底部角点的竖向位移

Figure 13. Vertical displacements of the corner point in the layered half-space model

下载: 全尺寸图片幻灯片

由表5～6和图12～13，可以得到与第4.1节相似的结论：本文中稳定性条件解析解可准确判断采用不同积分时间步长进行显式动力计算时整体模型的稳定性状态；整体模型中同一层介质的稳定性由角点区域控制，在实际工程应用中，令积分时间步长满足式(28)的稳定性条件，或将积分时间步长取为由式(29)给出的稳定性系数建议值与内部计算域临界时间步长的乘积，即可顺利完成整体模型的显式动力计算。

5. 结　论

以二维黏弹性人工边界为研究对象，利用基于局部子系统的稳定性分析方法研究其在显式时域逐步积分算法中的稳定性，并通过数值算例对稳定性分析的准确性加以验证。具体结论如下。

（1）给出了采用黏弹性人工边界时，模型侧边和角点区域在显式时域逐步积分算法中的稳定性条件解析解，获得了含黏弹性人工边界的整体模型在显示动力计算中的统一稳定性判别准则，在数值模拟时可利用该解析解对整体模型的稳定性进行预判。

（2）考虑黏弹性人工边界影响时，整体模型的数值积分稳定性条件比内部计算域更严格，整体模型的稳定性由角点区域控制。侧边和角点区域的稳定性条件仅与波源距与单元尺寸比R/L、泊松比μ和内部计算域的稳定时间步长L/c_P有关。随着R/L和μ的增大，人工边界区的稳定性变得宽松，当R/L>5时，稳定性条件基本保持不变。

（3）分析稳定性条件的控制因素及其影响规律的，给出了采用黏弹性人工边界时显式时域逐步积分算法稳定性系数的保守建议值。在实际应用中，可直接采用该稳定性系数建议值与内部计算域临界时间步长的乘积，作为显式动力计算最大稳定时间步长的估计。

图 1 黏弹性人工边界

Figure 1. Schematic diagram of viscoelastic artificial boundaries

下载: 全尺寸图片幻灯片

图 2 基于局部子系统的稳定性分析

Figure 2. Stability analysis based on local subsystems

下载: 全尺寸图片幻灯片

图 3 二维半无限空间近场有限元模型及二维平面应变单元

Figure 3. Two-dimensional semi-infinite near-field finite element model and two-dimensional plane strain element

下载: 全尺寸图片幻灯片

图 4 稳定性条件的比较

Figure 4. Comparison of stability conditions

下载: 全尺寸图片幻灯片

图 5 不同R/L时稳定性系数γ随泊松比μ的变化

Figure 5. Variations of stability coefficient γ with Poisson’s ratio μ under different R/L

下载: 全尺寸图片幻灯片

图 6 不同泊松比μ时稳定性系数γ随R/L的变化

Figure 6. Variations of stability coefficient γ with R/L under different Poisson’s ratio μ

下载: 全尺寸图片幻灯片

图 7 均匀半空间模型

Figure 7. The homogeneous half-space model

下载: 全尺寸图片幻灯片

图 8 脉冲荷载

Figure 8. The impulse load

下载: 全尺寸图片幻灯片

图 9 均匀半空间模型的位移分布

Figure 9. Displacement distributions of the homogeneous half-space model

下载: 全尺寸图片幻灯片

图 10 均匀半空间模型底部角点的竖向位移

Figure 10. Vertical displacements of the corner point in the homogenous half-space model

下载: 全尺寸图片幻灯片

图 11 成层半空间模型

Figure 11. The layered half-space model

下载: 全尺寸图片幻灯片

图 12 成层半空间模型的位移分布

Figure 12. Displacement distributions of the layered half-space model

下载: 全尺寸图片幻灯片

图 13 成层半空间模型底部角点的竖向位移

Figure 13. Vertical displacements of the corner point in the layered half-space model

下载: 全尺寸图片幻灯片

表 1 二维黏弹性人工边界参数的数据^[23]

Table 1. The values of two-dimensional viscoelastic artificial boundary coefficients^[23]

参数	范围	建议
α_T	0.35～0.65	0.5
α_N	0.80～1.20	1.0

下载: 导出CSV

表 2 建议的几种常见情况的稳定性系数

Table 2. Recommended stability coefficients for several common cases

R/L				γ
R/L	μ=0.10	μ=0.15	μ=0.20	μ=0.25	μ=0.30	μ=0.35	μ=0.40
1	0.47	0.48	0.49	0.50	0.51	0.53	0.55
5	0.50	0.50	0.51	0.52	0.53	0.55	0.57
10	0.50	0.51	0.51	0.52	0.53	0.55	0.57
20	0.50	0.51	0.52	0.53	0.54	0.55	0.57
50	0.50	0.51	0.52	0.53	0.54	0.55	0.57
+∞	0.51	0.51	0.52	0.53	0.54	0.55	0.57

下载: 导出CSV

表 3 均匀半空间模型的稳定性系数和最大稳定时间步长

Table 3. Stability coefficients and maximum stable time steps of the homogeneous model

模型分区	稳定性系数γ		最大稳定时间步长Δt/s
模型分区	计算	建议	计算	建议
内部区域	1.00	0.50	0.002 7	0.001 35
侧边子系统	0.74	0.50	0.002 0	0.001 35
角点子系统	0.52	0.50	0.001 4	0.001 35

下载: 导出CSV

表 4 不同固定时间步长时均匀半空间模型的稳定性状态

Table 4. The stability states of the homogeneous model under different fixed time steps

时间步长 ${{\Delta }}t$ /s	稳定性系数γ	稳定性状态
0.003 0	1.11	内部首先失稳
0.002 7	1.00	侧边首先失稳
0.002 0	0.74	角点首先失稳
0.001 4	0.52	稳定计算
0.001 35	0.50	稳定计算

下载: 导出CSV

表 5 成层半空间模型的稳定性系数与最大稳定时间步长

Table 5. Stability coefficients and maximum stable time steps of the layered model

介质	模型分区	稳定性系数γ		最大稳定时间步长Δt/s
介质	模型分区	计算	建议	计算	建议
上层	内部区域	1.00	0.5	0.003 7	0.000 95
上层	侧边子系统	0.76	0.5	0.002 8	0.000 95
下层	内部区域	1.00	0.5	0.001 9	0.000 95
	侧边子系统	0.79	0.5	0.001 5	0.000 95
	角点子系统	0.59	0.5	0.001 1	0.000 95

下载: 导出CSV

表 6 采用不同固定时间步长时成层半空间模型的稳定性状态

Table 6. The stability state of the layered model under different fixed time steps

时间步长 ${{\Delta }}t$ /s	稳定性系数γ	稳定性状态
0.002 0	1.05	内部首先失稳
0.001 9	1.00	侧边首先失稳
0.001 5	0.79	角点首先失稳
0.001 0	0.53	稳定计算
0.000 95	0.50	稳定计算

下载: 导出CSV

参考文献(30)

[1]	LIAO Z P, WONG H L. A transmitting boundary for the numerical simulation of elastic wave propagation [J]. International Journal of Soil Dynamics and Earthquake Engineering, 1984, 3(4): 174–183. DOI: 10.1016/0261-7277(84)90033-0.
[2]	廖振鹏, 周正华, 张艳红. 波动数值模拟中透射边界的稳定实现 [J]. 地球物理学报, 2002, 45(4): 533–545. DOI: 10.3321/j.issn:0001-5733.2002.04.011. LIAO Z P, ZHOU Z H, ZHANG Y H. Stable implementation of transmitting boundary in numerical simulation of wave motion [J]. Chinese Journal of Geophysics, 2002, 45(4): 533–545. DOI: 10.3321/j.issn:0001-5733.2002.04.011.
[3]	LYSMER J, KUHLEMEYER R L. Finite dynamic model for infinite media [J]. Journal of the Engineering Mechanics Division, 1969, 95(4): 859–878. DOI: 10.1061/JMCEA3.0001144.
[4]	刘晶波, 王振宇, 杜修力, 等. 波动问题中的三维时域粘弹性人工边界 [J]. 工程力学, 2005, 22(6): 46–51. DOI: 10.3969/j.issn.1000-4750.2005.06.008. LIU J B, WANG Z Y, DU X L, et al. Three-dimensional visco-elastic artificial boundaries in time domain for wave motion problems [J]. Engineering Mechanics, 2005, 22(6): 46–51. DOI: 10.3969/j.issn.1000-4750.2005.06.008.
[5]	刘晶波, 李彬. 三维黏弹性静-动力统一人工边界 [J]. 中国科学: E辑, 2005, 35(9): 966–908. LIU J B, LI B. Three-dimensional viscoelastic static-dynamic unified artificial boundary [J]. Science in China: Series E, 2005, 35(9): 966–908.
[6]	DU X L, ZHAO M. A local time-domain transmitting boundary for simulating cylindrical elastic wave propagation in infinite media [J]. Soil Dynamics and Earthquake Engineering, 2010, 30(10): 937–946. DOI: 10.1016/j.soildyn.2010.04.004.
[7]	BREBBIA C A. The boundary element method for engineers [M]. London: Wiley, 1978.
[8]	金峰, 王光纶, 贾伟伟. 离散元-边界元动力耦合模型在地下结构动力分析中的应用 [J]. 水利学报, 2001(2): 24–28. DOI: 10.3321/j.issn:0559-9350.2001.02.004. JIN F, WANG G L, JIA W W. Application of distinct element-boundary element coupling model in underground structure dynamic analysis [J]. Journal of Hydraulic Engineering, 2001(2): 24–28. DOI: 10.3321/j.issn:0559-9350.2001.02.004.
[9]	BERENGER J P. Perfectly matched layer for the FDTD solution of wave-structure interaction problems [J]. IEEE Transactions on Antennas and Propagation, 1996, 44(1): 110–117. DOI: 10.1109/8.477535.
[10]	KUCUKCOBAN S, KALLIVOKAS L F. A symmetric hybrid formulation for transient wave simulations in PML-truncated heterogeneous media [J]. Wave Motion, 2013, 50(1): 57–79. DOI: 10.1016/j.wavemoti.2012.06.004.
[11]	李彬, 刘晶波. 粘弹性人工边界在Marc中的实现 [C] // 第14届全国结构工程学术会议论文集(第一册). 烟台: 中国力学学会工程力学编辑部, 2005: 303−307.
[12]	胡汛训, 张燎军, 华慧玲. 粘弹性人工边界在LS-DYNA中的实现 [C] // 首届全国水工抗震防灾学术会议论文集. 南京: 中国水力发电工程学会, 2006: 134−139.
[13]	刘晶波, 杜义欣, 闫秋实. 粘弹性人工边界及地震动输入在通用有限元软件中的实现 [C] // 第三届全国防震减灾工程学术研讨会论文集. 南京: 中国土木工程学会, 2007: 43−48.
[14]	张燎军, 张慧星, 王大胜, 等. 黏弹性人工边界在ADINA中的应用 [J]. 世界地震工程, 2008, 24(1): 12–16. ZHANG L J, ZHANG H X, WANG D S, et al. The application of artificial viscous-spring boundary in ADINA [J]. World Earthquake Engineering, 2008, 24(1): 12–16.
[15]	张焜煌, 钱彦岭, 徐慧峰, 等. 基于Nastran的粘性和粘弹性人工边界的模拟与验证 [J]. 兵工自动化, 2009, 28(3): 29–31. DOI: 10.3969/j.issn.1006-1576.2009.03.011. ZHANG K H, QIAN Y L, XU H F, et al. Simulation and verification of viscous and viscous-spring artificial boundary based on Nastran [J]. Ordnance Industry Automation, 2009, 28(3): 29–31. DOI: 10.3969/j.issn.1006-1576.2009.03.011.
[16]	HUANG J Q, DU X L, JIN L, et al. Impact of incident angles of P waves on the dynamic responses of long lined tunnels [J]. Earthquake Engineering and Structural Dynamics, 2016, 45(15): 2435–2454. DOI: 10.1002/eqe.2772.
[17]	WANG J T, ZHANG C H, JIN F. Nonlinear earthquake analysis of high arch dam-water-foundation rock systems [J]. Earthquake Engineering & Structural Dynamics, 2012, 41(7): 1157–1176. DOI: 10.1002/eqe.1178.
[18]	郜新军, 赵成刚, 张延. 多源散射黏弹性叠加人工边界探究及在桥梁工程中的应用 [J]. 土木工程学报, 2010, 43(11): 130–138. DOI: 10.15951/j.tmgcxb.2010.11.005. GAO X J, ZHAO C G, ZHANG Y. A study of viscous-spring superposition artificial boundary for multi-source scattering problems and its application in bridge engineering [J]. China Civil Engineering Journal, 2010, 43(11): 130–138. DOI: 10.15951/j.tmgcxb.2010.11.005.
[19]	李忠诚, 凡红. 基于粘弹性人工边界的核电工程地基动力阻抗分析 [J]. 核动力工程, 2014, 35(2): 67–70. LI Z C, FAN H. Dynamic impedance analysis based on an artificial viscoelastic boundary technology for nuclear power engineering [J]. Nuclear Power Engineering, 2014, 35(2): 67–70.
[20]	王子辉. 饱和两相与单相土互层场地中地铁车站地震反应分析 [D]. 北京: 北京交通大学, 2008: 95−96.
[21]	闫秋实. 典型地铁结构内爆炸流场分布及动力反应研究 [D]. 北京: 清华大学, 2011: 100−102.
[22]	BAO X, LIU J B, LI S T, et al. Seismic response analysis of the reef-seawater system under obliquely incident P and SV waves [J]. Ocean Engineering, 2020, 200: 107021. DOI: 10.1016/j.oceaneng.2020.107021.
[23]	刘晶波, 谷音, 杜义欣. 一致粘弹性人工边界及粘弹性边界单元 [J]. 岩土工程学报, 2006, 28(9): 1070–1075. DOI: 10.3321/j.issn:1000-4548.2006.09.004. LIU J B, GU Y, DU Y X. Consistent viscous-spring artificial boundaries and viscous-spring boundary elements [J]. Chinese Journal of Geotechnical Engineering, 2006, 28(9): 1070–1075. DOI: 10.3321/j.issn:1000-4548.2006.09.004.
[24]	KAMEL A H. A stability checking procedure for finite-difference schemes with boundary conditions in acoustic media [J]. Bulletin of the Seismological Society of America, 1989, 79(5): 1601–1606. DOI: 10.1785/BSSA0790051601.
[25]	关慧敏, 廖振鹏. 局部人工边界稳定性的一种分析方法 [J]. 力学学报, 1996, 28(3): 376–380. DOI: 10.6052/0459-1879-1996-3-1995-344. GUAN H M, LIAO Z P. A method for the stability analysis of local artificial boundaries [J]. Acta Mechanica Sinica, 1996, 28(3): 376–380. DOI: 10.6052/0459-1879-1996-3-1995-344.
[26]	李述涛, 刘晶波, 宝鑫, 等. 采用粘弹性人工边界单元时显式算法稳定性分析 [J]. 工程力学, 2020, 37(11): 1–11,46. DOI: 10.6052/j.issn.1000-4750.2019.12.0755. LI S T, LIU J B, BAO X, et al. Stability analysis of explicit algorithms with visco-elastic artificial boundary elements [J]. Engineering Mechanics, 2020, 37(11): 1–11,46. DOI: 10.6052/j.issn.1000-4750.2019.12.0755.
[27]	Abaqus analysis user’s manual (version 6.14) [Z]. 2013.
[28]	HUGHES T J R. Analysis of transient algorithms with particular reference to stability behavior [C] // BELYTSCHKO T, HUGHES T J R. Computational methods for transient analysis. Amsterdam: Elsevier, 1983: 67−155.
[29]	杜修力. 工程波动理论与方法 [M]. 北京: 科学出版社, 2009: 215−216.
[30]	王勖成, 邵敏. 有限单元法基本原理和数值方法 [M]. 北京: 清华大学出版社, 1997: 66−67.

施引文献

期刊类型引用(12)

1.	常白雪，张元瑞，王少华，彭克锋，虞吉林，郑志军. 梯度多胞材料的动态力学性能分析与设计研究综述. 爆炸与冲击. 2024(08): 5-33 . 本站查看
2.	袁良柱，陈美多，谢雨珊，陆建华，王鹏飞，徐松林. 细观非连续介质的应力波传播研究. 爆炸与冲击. 2024(09): 50-61 . 本站查看
3.	周辉，任辉启，吴祥云，易治，黄魁，穆朝民，王海露. 成层式防护结构中分散层研究综述. 爆炸与冲击. 2022(11): 3-28 . 本站查看
4.	余同希，朱凌，许骏. 结构冲击动力学进展（2010-2020）. 爆炸与冲击. 2021(12): 4-64 . 本站查看
5.	刘冕，王根伟，宋辉，王彬. 负梯度泡沫金属中的局部密实化现象. 高压物理学报. 2020(04): 117-127 . 百度学术
6.	常白雪，郑志军，赵凯，何思渊. 具有恒定冲击载荷的梯度泡沫金属材料设计. 爆炸与冲击. 2019(04): 3-11 . 本站查看
7.	胡俊，任建伟，马巍，刘建华，王爱国. 冲击荷载下含随机缺陷的梯度蜂窝材料的力学性能. 材料导报. 2019(16): 2777-2784 . 百度学术
8.	胡俊，任建伟，王爱国，吴德义. 非线性梯度胞元分布蜂窝材料的冲击力学响应. 材料导报. 2019(24): 4066-4071 . 百度学术
9.	罗伟铭，石少卿，廖瑜，孙建虎. 成层式铝蜂窝夹芯结构冲击响应试验研究. 材料导报. 2018(08): 1328-1332 . 百度学术
10.	常白雪，郑志军，赵凯，虞吉林. 梯度多胞材料耐撞性设计的简化模型和渐近解. 中国科学:物理学力学天文学. 2018(09): 233-241 . 百度学术
11.	Xu-ke Lan，Shun-shan Feng，Qi Huang，Tong Zhou. Blast response of continuous-density graded cellular material based on the 3D Voronoi model. Defence Technology. 2018(05): 433-440 . 必应学术
12.	李志斌，卢芳云. 梯度温度场中多胞材料牺牲层的抗冲击分析. 兵工学报. 2017(12): 2463-2471 . 百度学术