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

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

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全文浏览量: 252
- PDF下载量: 42
- 被引次数: 0
出版历程
- 收稿日期: 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全文

当冲击作用于不同密度或不同可压缩性2种物质的扰动界面时产生Richtmyer-Meshkov不稳定性(RMI)。这种不稳定性理论上由Richtmyer发现并描述^[1], 由Meshkov从实验中证实^[2]。该界面不稳定性问题在许多自然现象及科学和工程领域中起着重要作用^[3-8], 如超新星爆炸、磁化等离子体、磁约束、太阳磁化层、地下盐矿、火山岛及外壳与内部流体混合导致中子收益降低的激光驱动惯性约束聚变和冲击波与火焰相互作用导致的爆燃转爆轰等。此外, RMI也可能从受冲击的金属表面产生喷射物。

RMI的演化通常经历由不稳定模式的振幅h_k和波长λ=2π/k描绘的若干阶段。对于kh_k≪1, 扰动随kUt呈正比例增长, U为激波作用后的界面运动速度。当kh_k达到某一值, 非线性使增长率降低, 驱动模式耦合, 且增长率随着扰动谱宽的增大而减小。然后, 由于尖钉下落(重流体进入轻流体)比气泡上升(轻流体进入重流体)快, 界面变得不对称。对于宽的不稳定谱, 非线性最终导致产生湍流混合层。RMI的脉冲性质令问题复杂, 使得RMI定性上与常见的Rayleigh-Taylor不稳定性(RTI)不同。由于冲击的可压缩性、复杂的物质特性以及后期的非线性运动直至湍流混合, RMI演化的计算是困难的。当然, 随着计算机技术的迅猛发展, 这可以采用多维高分辨率流体力学模拟来进行, 但它们计算强度大, 无法用于工程设计优化研究。因此, 目前实际应用中通常采用捕捉较低分辨率时不稳定流动主要特征的简化“混合模型”^[9]。杨玟等^[10-11]对此进行了尝试, 将传统的k-ε模型应用于界面不稳定性引起的混合, 取得了令人满意的结果。

但是, 由于与RMI相关的其它物理过程非常复杂, 较复杂的混合模型(如k-ε)也难以直接应用到工程设计中。目前, 很多实际应用中对RMI诱发混合现象的处理都非常简单, 假设混合层宽度以指数形式t^θ_i增长。而大量实验研究表明该比例关系仍不确定^[3-6], 因为考虑压缩性的计算是困难的, 它们与实验不符。即使指数律粗略满足, 但不同工况下θ_i的差别也很大, 它显著依赖于初始扰动谱。由此可见, 工程设计中对RMI诱发混合现象的处理过于粗糙。

本文中, 在简单介绍描述作用于混合层中产生的气泡和尖钉的浮阻力模型基础上, 采用该模型对激波管低压缩情况和激光加载高压缩情况下的RMI诱发混合层宽度(气泡与尖钉宽度之和)进行计算, 验证模型和选取参数的有效性。

1. 模型介绍

目前, 典型的浮阻力模型可写为如下形式^[12]:

$\left.\left(\rho_{i}+C_{\mathrm{a}} \rho_{j}\right) \frac{\mathrm{d} v_{i}}{\mathrm{~d} t}-\beta\left(\rho_{i}-\rho_{j}\right) a(t)\right] V_{i}=-C_{\mathrm{d}} \rho_{i} v_{i}\left|v_{i}\right| A_{i} \quad i, j=1, 2 ; i \neq j$

(1)

式中：下角标i, j表示2种不同的流体，下角标为1时表示重流体(尖钉)，为2时表示轻流体(气泡)；ρ_i为重流体/轻流体的密度；v_i是尖钉/气泡的渗透速度，且v_i=dh_i/dt，h_i表示尖钉/气泡的瞬时宽度；C_a是附加质量力系数；β是浮力产生的模型常数；C_d是阻力系数；a(t)为激波脉冲加速度；V_i为尖钉/气泡的体积，A_i为尖钉/气泡的截面积。方程左端第一项为惯性力，第二项为浮力，右端为阻力。关于模型的详细论述可参考文献[13-14]，这里不再重复。对于Richtmyer-Meshkov不稳定性，通常认为冲击简单地给予界面上的气泡和尖钉一个脉冲，则它们随后的运动可以由惯性力和阻力相等来得到(加速度为零)。因此脉冲加速度情况是有启发性的，可以用来研究不稳定性的惯性特性。

本文所求解的模型方程是一组二阶常微分方程, 将它们简化为一阶微分方程:dh_i/dt=v_i; dv_i/dt= -f_iC_dv_i|v_i|/h_i。采用四阶Runge-Kutta方法进行求解。

2. 结果分析与讨论

采用上述模型和数值方法, 对关注的激波管低压缩情况和激光加载高压缩情况下模型的性能进行了考察。这2种工况下RMI产生的机理不同:对于弱冲击, 主要贡献来自于压力梯度和密度梯度不重合引起的旋涡沉积; 对于强冲击, 存在激波在经折射后产生了显著的反射, 这产生增长率的振荡, 但它们最终衰减。

2.1 激波管低压缩的情况下的模拟

首先采用上述模型对4种不同激波脉冲加速度情况下气泡和尖钉宽度进行了计算。图 1给出了所采用的4种加速度曲线，g为重力加速度。脉冲加速度a约为150g，持续时间t₀约为10 ms。这些曲线为LANL的Dimonte等LEM(Linear Electric Motor)实验的测量曲线^[15]。实验中流体和脉冲加速度的性质参数见表 1，其中R为密度比，R=(1+A)(1-A)，A为Atwood数，A=(ρ₂-ρ₁)(ρ₂+ρ₁)，We为韦伯数，Re为雷诺数。对于每一种情况，通过调整阻力系数C_d和初始振幅h_i0来使随时间变化的解与实验数据相符。但是，数值实验发现：在大多数情况下h_i0对结果的影响远小于C_d的影响。

图 1 计算采用的4种不同脉冲加速度曲线

Figure 1. Four kinds of impulsive accelerations used in the calculation

下载: 全尺寸图片幻灯片

表 1 实验中采用的流体和脉冲加速度性质参数

Table 1. Fluid combinations and characteristics for impusive accerleration experiments

No.	流体1	流体2	ρ₁/(g·cm^-3)	ρ₂/(g·cm^-3)	R	A	We	Re
1	H₂O	CCl₂F₂	1.000	1.57	1.57	0.22	4 000	2 600
2	SF₆	C₄H₁₀	0.067	0.81	12.10	0.85	1 100	8 000
3	SF₆	CCl₂F₂	0.067	1.57	23.40	0.92	11 000	23 000
4	SF₆	CCl₂F₂	0.032	1.57	49.10	0.96	6 000	25 000

下载: 导出CSV

| 显示表格

给出了4种加速度驱动下气泡和尖钉宽度随位移Z的变化, $Z=\iint a \mathrm{~d} t^{\prime} \mathrm{d} t$ , 激波作用时Z≈Ut。由图可见, 4种加速度情况下计算的气泡和尖钉宽度与实验基本吻合。计算中阻力系数C_d的取值为3.67±0.73, 与文献[16]中分析得到的C_d的不确定度1.2接近。从图中还可看出:气泡和尖钉的不对称性随着密度比R的增大而增大。此外, 本文中还对实验结果按指数律h_i=h_i0t^θ_i进行了拟合, 其中h_i0的取值范围为0.5~1.0 cm。R=49.1时, θ₁≈0.85, θ₂≈0.33; R=23.4时, 指数迅速下降, θ₁≈0.45, θ₂≈0.24; R=1.57时, θ₁≈0.28, θ₂≈0.22。由此可见, 指数θ_i随密度比变化而变化, 但具体变化规律还未从数值模拟和实验中最终确定, 这主要是由于θ_i对实验初始条件敏感, 需要计算和实验之间更直接的比较。

图 2 气泡和尖钉宽度随位移的变化

Figure 2. The width of bubble and spike with displacement

下载: 全尺寸图片幻灯片

2.2 激光加载高压缩情况下的模拟

为了考察模型在高压缩情况下的性能，我们进一步对Nova激光器上马赫数Ma>10的实验进行了模拟。实验采用一靶丸装置在Nova激光器上进行^[17]。流体1由厚度为125 μm、初始密度为1.7 g/cm³的铍烧蚀层组成。流体2是未压缩密度为0.12 g/cm³的泡沫。波速为46 km/s的入射冲击与界面相互作用产生反射稀疏波和速度为3 km/s的透射激波。界面经加速后速度为56 km/s，物质被压缩后，ρ₁=2 g/cm³，ρ₂=0.5 g/cm³，A=-0.6。这些参数通过对比热比γ₁=1.8和γ₂=1.45的流体求解理想的黎曼问题得到。

图 3给出了Nova实验中计算的加速度曲线。由图可见，激光驱动在4 ns后停止，这导致泡沫减压，由于A < 0而产生Rayleigh-Taylor(RT)分量，因此冲击压缩后流动是亚音速的，本文模型是适用的。图 4给出了混合区总宽度H随位移Z的变化(由于实验不能分辨气泡和尖钉，因此给出了总振幅H)。从图中可看出：混合区总宽度的计算值与实验值吻合，而且C_d=2.0和C_d=5.36的曲线之间包括了全部的实验数据。但是，阻力系数C_d的不确定度约为3.36，明显大于低压缩情况的值(约为1.46)。此外，拟合得到总的混合宽度以指数为0.5的指数律增长，这超过了激波管低压缩时得到的指数，推测其原因可能是：(1)激光驱动随时间减小，使得压力降低、界面减速，这导致扰动膨胀，并引入RT分量(因为Aa>0)。这些影响可能显著增加推测的指数；(2) A=0.6时Nova上的初始扰动比激波管上的更对称，如果指数对初始条件敏感，这可能导致不同的指数。

图 3 Nova实验中的加速度曲线

Figure 3. Acceleration history for Nova experiment

下载: 全尺寸图片幻灯片

图 4 混合区宽度随位移的变化

Figure 4. Variaion of total width with displacement

下载: 全尺寸图片幻灯片

3. 结论

采用浮阻力模型对激波管低压缩和激光加载高压缩情况下Richtmyer-Meshkov不稳定性诱发的物质渗透边界的演化过程进行了计算, 计算结果与实验吻合得较好。这表明本研究中模型参数的选取、方程中现象学比例因子的添加和模型假设是合适的。但是由于实验测量的局限性, 模型中的一些问题仍然是突出的, 包括阻力项的大小和形式、压缩的影响、“附加质量”的描述等。为了更好地评估模型, 需要一些实验上的完善。首先, 气泡和尖钉必须单独分辨, 因为它们的表现相当不同, 尤其在A较大的情况。其次, 实验持续时间应当延长至足以揭示模型的差别为止。尽管如此, 本文模型仍明显优于当前实际应用中所采用的经验公式(本研究也显示指数θ_i随工况的不同而显著变化)。

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