AN INPUT METHOD FOR EXPLOSION PROBLEMS BASED ON THE ACOUSTIC SUBSTRUCTURE OF EXPLOSION SOURCE UNDER UNDERWATER EXPLOSIVE LOADING
-
摘要: 无限域吸收边界和爆炸波输入是水下爆炸载荷数值模拟的两个关键问题. 本文借鉴基于内部子结构的地震波动输入方法和爆源子结构多尺度分析方法, 考虑水下爆炸载荷与地震载荷同属于波动问题, 提出一种水下爆炸作用下爆源子结构的爆炸波输入方法, 该方法首先通过对爆源区域进行水下爆炸波自由场的波场分解, 并利用地震波的内部子结构输入方法, 将该自由波场运动转换为等效爆炸波输入载荷, 从而实现了水下爆炸问题中冲击波的输入; 本文采用圆形爆源子结构区域, 并采用AUTODYN软件中一维自由场模型计算水下爆炸作用下该区域的自由波场压力时程. 进一步, 基于连分式近似方法提出了一种模拟无限域水体辐射效应的高精度时域人工边界条件, 该吸收边界可设置于离结构和爆源子结构较近处. 本文提出的分析方法通过圆形爆源子结构进行水下爆炸载荷的转化, 并采用高精度吸收边界大大减少计算区域, 既保证了计算精度, 又降低了单元数量, 具有较高的计算效率和较强的实用性. 通过数值算例分析, 验证了本文模型与方法的准确性, 模拟了水下爆炸作用下的冲击波和气泡脉动阶段测点的压力时程曲线, 并研究了圆形结构对水下近场爆炸波散射效应的影响规律.Abstract: The infinite domain absorption boundary and explosion wave input method are two keys of numerical simulation research on underwater explosion. This paper draws on the seismic wave input method based on the internal substructure and the multiscale analysis method for explosion problems based on the substructure of explosion source and proposes a shock wave input method for underwater explosion based on the substructure of explosion source considering that the underwater explosion load and the seismic load belong to the same fluctuation problem. The method first decomposes the shock wave field in the explosion source region and the free wave field motions were transformed into the equivalent explosive loads, which enables the input of shock waves in the underwater explosion problem. In this paper, a circular substructure of explosion source is used, and a one-dimensional model in AUTODYN software is used to calculate the free shock wave pressure in this region under the action of the underwater explosion. Further, based on the continuous fractional approximation method, a high-precision time-domain artificial boundary condition is proposed to simulate the radiation effects of infinite domain, which can be placed close to the structure and explosive source substructure is proposed. The method proposed in this paper transforms the underwater explosion load through the substructure of explosion source and adopts the high accuracy absorption boundary to greatly reduce the calculation area, which not only ensures the calculation accuracy, but also reduces the number of elements, with high calculation efficiency and practicability. Finally the numerical example is analyzed to verify the accuracy of the model and method of this paper, The pressure time history curves of measuring points during shock wave and bubble pulsation stage under underwater explosion are simulated, and the effect of circular structure on the scattering effect of underwater near-field explosion wave is studied.
-
引 言
桥梁、码头、大坝等水工结构是影响国计民生和国防军事的重要工程, 常是恐怖分子的重点打击对象, 因此水下爆炸问题一直是研究热点[1-4]. 现阶段水下爆炸的分析方法主要有物理试验、经验公式、数值模拟三种方式. 目前的经验公式主要适用于远场爆炸[5], 基于经验公式计算水下结构的爆炸响应对爆炸距离和结构形状有限制[6-7]; 试验研究常开展小药量炸药试验且试件小, 也未得到有效使用的设计载荷公式[8-10]; 近几十年, 随着计算方法和计算机技术发展进步, 水下爆炸数值模拟推动了水中爆炸研究的快速发展.
数值模拟方法可以对水下爆炸的各个阶段及各自涉及的问题进行系统分析, 研究人员多采用一些商业软件如ABAQUS, LS-DYNA, AUTODYN等进行水下爆炸的数值模拟研究 [11-14]. 尹训强等[15]利用ANSYS/LS-DYNA模拟了近场水下爆炸载荷作用下的近岸结构的动力响应及冲击波在近岸场地的传播特性和规律; Barra等[16]采用ALE算法模拟了水下爆炸气泡脉动过程的水动力响应; 肖秋平等[17] 利用 AUTODYN程序模拟水下爆炸冲击波,分析了模型网格的大小和状态方程的选取对数值模拟结果精度的影响; 闫秋实等[18]通过AUTODYN模拟了近场情况下水下爆炸对高桩码头钢筋混凝土桩毁伤情况, 并分析了单桩在一定距离和不同爆炸深度时的损伤模式; 胡延东等[19]建立了结构在静水压力作用下内爆数值模型, 模拟了裂缝处固−水/气流体的弱耦合, 分析了压力脉冲的传播规律; 张阿漫等[20]采用SPH和RKPM耦合的方法, 对近场水下爆炸作用下的结构损伤进行了数值模拟,并针对水下爆炸载荷引起的大变形, 提出了一种基于RKPM壳体公式的接触检测算法. 但数值模拟方法需要建立一定范围区域的模型, 且区域网格尺寸需要足够小以免过滤掉冲击波峰值, 导致计算量大, 时间长,效率低. 本文方法采用高精度时域人工边界条件吸收散射波, 并且人工边界可以设置的离结构近处, 保证了精度,提高了计算效率.
另外, 基于声固耦合的数值方法也常用于水下爆炸载荷的计算[21-24]. 姚熊亮等[25-26]通过ABAQUS声固耦合法模拟了水下爆炸载荷作用下舰船结构的动力响应, 分析了不同工况对结构动力响应的影响并给出模拟处理方法, 计算结果与实测数据基本吻合. 金乾坤等[27]采用ABAQUS声−固耦合作用对舱段动态响应进行了数值模拟,对比了加速度和速度响应的模拟和试验结果, 二者符合得较好. 声固耦合方法能够提高计算效率, 保证冲击波传播的计算精度, 但是在水下爆炸方面的研究主要集中在远场爆炸, 对近场水下爆炸问题研究的较少.
爆炸波是由内部向四周传播的球状冲击波, 与地震波同属于波动问题, 因此针对爆炸载荷的输入问题可以借鉴地震波输入方法的思想[28]. 刘晶波等[29]在人工边界处输入地震动等效载荷, 提出了一种适用于黏弹性人工边界的地震波动输入方法, 即波动法. 基于波动法的原理, 刘晶波等[30-32]又提出基于人工边界子结构的地震波动输入方法, 也称为人工边界子结构法. 但是以上两种方法中人工边界均参与了等效输入载荷的计算, 计算复杂, 由此刘晶波等[33]又提出了基于内部子结构的地震波动输入方法, 此方法在模型内部取子结构对子结构模型的单元节点施加自由场位移并进行动力分析得出等效输入地震载荷, 使地震波的外源输入问题转化为内源波动问题. 有效提高了地震波输入的效率,同时摆脱了输入方法对于人工边界条件的限制,具有较高的计算精度. 本文基于声学波动方程和内部子结构地震波输入方法, 提出了水下爆炸作用下的爆源子结构输入方法, 模拟远场爆炸同时也可以有效模拟近场水下爆炸问题.
1. 基于声学的水下爆炸有限元方法
1.1 水体有限元方程
水体的声学方程为
$$ \frac{{{\partial ^2}p}}{{\partial {x^2}}} + \frac{{{\partial ^2}p}}{{\partial {y^2}}} = \frac{1}{{{c^2}}}\frac{{{\partial ^2}p}}{{\partial {t^2}}} $$ (1) 式中, c表示水中声速.
基于有限元方法, 采用四边形线性单元对计算水域进行有限元离散[34], 则方程(1)的有限元形式为
$$ {{{\boldsymbol{M}}\ddot {\boldsymbol{p}}}} + {\boldsymbol{Kp}} = {\boldsymbol{F}} $$ (2) 其中
$$ {\boldsymbol{M}} = \iint {{\boldsymbol{N}}_{}^{{\rm{T}}} {\boldsymbol{N}}{\text{d}}V} $$ (3) $$ {\boldsymbol{K}} = \frac{1}{{{c^2}}}\iint {\left( {\frac{{\partial {\boldsymbol{N}}_{}^{{\rm{T}}} }}{{\partial x}}\frac{{\partial {\boldsymbol{N}}}}{{\partial x}} + \frac{{\partial {\boldsymbol{N}}_{}^{{\rm{T}}} }}{{\partial y}}\frac{{\partial {\boldsymbol{N}}}}{{\partial y}}} \right)}{\text{d}}V $$ (4) 式中, N为形函数矩阵, 上标T表示向量转置, M和K分别为水体质量和刚度矩阵, F表示力向量.
1.2 水下爆炸爆源子结构输入方法
根据地震波动内部子结构输入理论[32], 水下爆炸以爆炸源为中心向四周传播的冲击波在进行载荷输入时只与自由场有关, 而与水中结构无关. 参考李述涛等[28]针对土体中爆炸提出的爆源子结构理论, 可首先建立水下爆炸载荷下水体−结构模型(图1)相对应的爆源−自由场模型, 如图2所示, 该模型为图1中的黄色区域. 需要指出的是, 爆源−自由场模型与水体−结构数值模型中对应位置处的网格划分应保持一致. 另外, 根据爆炸波属于球面波将该二维自由场模型设置为圆形区域.
如图2所示, 数值模型中的单元节点按位置分为五类, 内部节点(I)、外部节点(E)和三层子结构节点(由内至外记为A, B, C). 水下爆炸载荷作用下, 爆源−自由场有限元模型的动力方程可表示为
$$ \begin{split} & \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{M}}_{{\text{EE}}}}}&{{{\boldsymbol{M}}_{{\text{AE}}}}}&0&0&0 \\ {{{\boldsymbol{M}}_{{\text{EA}}}}}&{{{\boldsymbol{M}}_{{\text{AA}}}}}&{{{\boldsymbol{M}}_{{\text{BA}}}}}&0&0 \\ 0&{{{\boldsymbol{M}}_{{\text{AB}}}}}&{{{\boldsymbol{M}}_{{\text{BB}}}}}&{{{\boldsymbol{M}}_{{\text{CB}}}}}&0 \\ 0&0&{{{\boldsymbol{M}}_{{\text{BC}}}}}&{{{\boldsymbol{M}}_{{\text{CC}}}}}&{{{\boldsymbol{M}}_{{\text{IC}}}}} \\ 0&0&0&{{{\boldsymbol{M}}_{{\text{CI}}}}}&{{{\boldsymbol{M}}_{{\text{II}}}}} \end{array}} \right]\left\{ {\begin{array}{*{20}{c}} {{{{\boldsymbol{\ddot p}}}_{\text{E}}}} \\ {{{{\boldsymbol{\ddot p}}}_{\text{A}}}} \\ {{{{\boldsymbol{\ddot p}}}_{\text{B}}}} \\ {{{{\boldsymbol{\ddot p}}}_{\text{C}}}} \\ {{{{\boldsymbol{\ddot p}}}_{\text{I}}}} \end{array}} \right\}+ \\ &\quad \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{K}}_{{\text{EE}}}}}&{{{\boldsymbol{K}}_{{\text{EA}}}}}&0&0&0 \\ {{{\boldsymbol{K}}_{{\text{AE}}}}}&{{{\boldsymbol{K}}_{{\text{AA}}}}}&{{{\boldsymbol{K}}_{{\text{AB}}}}}&0&0 \\ 0&{{{\boldsymbol{K}}_{{\text{BA}}}}}&{{{\boldsymbol{K}}_{{\text{BB}}}}}&{{{\boldsymbol{K}}_{{\text{BC}}}}}&0 \\ 0&0&{{{\boldsymbol{K}}_{{\text{CB}}}}}&{{{\boldsymbol{K}}_{{\text{CC}}}}}&{{{\boldsymbol{K}}_{{\text{CI}}}}} \\ 0&0&0&{{{\boldsymbol{K}}_{{\text{IC}}}}}&{{{\boldsymbol{K}}_{{\text{II}}}}} \end{array}} \right]\left\{ {\begin{array}{*{20}{c}} {{{\boldsymbol{p}}_{\text{E}}}} \\ {{{\boldsymbol{p}}_{\text{A}}}} \\ {{{\boldsymbol{p}}_{\text{B}}}} \\ {{{\boldsymbol{p}}_{\text{C}}}} \\ {{{\boldsymbol{p}}_{\text{I}}}} \end{array}} \right\} = \left\{ {\begin{array}{*{20}{c}} {{{\boldsymbol{F}}_{\text{E}}}} \\ {{{\boldsymbol{F}}_{\text{A}}}} \\ {{{\boldsymbol{F}}_{\text{B}}}} \\ {{{\boldsymbol{F}}_{\text{C}}}} \\ {{{\boldsymbol{F}}_{\text{I}}}} \end{array}} \right\} \end{split} $$ (5) 以节点集A为分界线, 将自由场模型分为内、外两部分. 对于外部区域, 采用子结构思想, 通过在节点集B上施加等效爆炸载荷FB, 即可保证节点集B外的区域运动与自由波长相一致, 即
$$ {{\boldsymbol{p}}_{\text{B}}} = {\boldsymbol{p}}_{\text{B}}^0{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} ,{\kern 1pt} {{\boldsymbol{p}}_{\text{C}}} = {\boldsymbol{p}}_{\text{C}}^0{\kern 1pt} ,{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {{\boldsymbol{p}}_{\text{E}}} = {\boldsymbol{p}}_{\text{E}}^0,{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {{\boldsymbol{F}}_{\text{C}}} = {\boldsymbol{0}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} ,{{\boldsymbol{F}}_{\text{E}}} = {\boldsymbol{0}} $$ (6) 式中, 上标0表示爆炸波自由场运动.
对于内部区域节点, 不考虑爆炸波在该区域的输入过程, 因而节点A和I的计算过程保持静止, 即
$$ \begin{array}{*{20}{c}} {{{\boldsymbol{p}}_{\text{A}}} = {\boldsymbol{0}}},&{{{\boldsymbol{p}}_{\text{I}}} = {\boldsymbol{0}}} \end{array} $$ (7) 将式(5)和式(6)代入式(4)整理得到
$$ \left. \begin{gathered} {{\boldsymbol{F}}_{\text{A}}} = {{\boldsymbol{M}}_{{\text{AB}}}} \cdot {\boldsymbol{\ddot p}}_{\text{B}}^0 + {{\boldsymbol{K}}_{{\text{AB}}}} \cdot {\boldsymbol{p}}_{\text{B}}^0 \\ {{\boldsymbol{F}}_{\text{B}}} = {{\boldsymbol{M}}_{{\text{BB}}}} \cdot {\boldsymbol{\ddot p}}_{\text{B}}^0 + {{\boldsymbol{M}}_{{\text{BC}}}} \cdot {\boldsymbol{\ddot p}}_{\text{C}}^0 + {{\boldsymbol{K}}_{{\text{BB}}}} \cdot {\boldsymbol{p}}_{\text{B}}^0 + {{\boldsymbol{K}}_{{\text{BC}}}} \cdot {\boldsymbol{p}}_{\text{C}}^0 \\ \end{gathered} \right\} $$ (8) 由以上推导可以看出, 如图2所示的自由场模型, 只需对节点A和B施加由式(8)所计算的等效载荷FA和FB, 即可完成水下爆炸波在子结构区域外的输入. 因此, 根据A, B, C节点上的质量和刚度矩阵以及B, C节点上自由场的压力及其压力对时间的两阶导数, 即可得到等效载荷FA和FB. 需要指出的是, 本文模型中B和C节点上自由场的压力是通过有限元软件ANSYS-AUTODYN一维水下爆炸模型计算得到.
将以上求得的等效载荷FA和FB施加在如图1所示的水体−结构相互作用数值模型的对应节点处, 即可将水下爆炸自由场运动通过等效载荷的形式输入至子结构外部的区域.
需要指出的是, 当爆炸波自由波场到达结构时会产生散射波场, 为防止散射波场到达计算水域边界时产生反射, 需要在水域截断边界层设置人工边界条件进行吸波, 从而使散射波场透过人工边界无反射传出. 传统的黏弹性人工边界[35]易于施加, 但其精度较低; 为保证计算精度, 往往将黏弹性人工边界设置在距离结构较远的位置. 本文将采用一种高精度的圆形时域人工边界条件模拟无限水域, 该人工边界条件可设置于离结构近处, 具体推导过程见第2节.
2. 高精度时域人工边界条件
2.1 人工边界频域解析解
极坐标系下, 式(1)在频域下控制方程为
$$ \frac{{{\partial ^2}P}}{{\partial {r^2}}} + \frac{1}{r}\frac{{\partial P}}{{\partial \theta }} + \frac{1}{{{r^2}}}\frac{{{\partial ^2}P}}{{\partial {\theta ^2}}} + \frac{{{\omega ^2}}}{{{c^2}}}P = 0 $$ (9) 式中,
$ \omega $ 为圆频率,$ P $ 为p的频域值.基于分离变量方法和无穷远辐射边界条件, 动水压力P的解为
$$ P = \sum\limits_{n = 1}^\infty {{d_n}{\rm{H}}_\beta ^{(2)}\left( {kr} \right){\phi _n}} $$ (10) $$ {\phi }_{n} = \left\{\begin{array}{lc}\mathrm{cos}(\beta \theta ),& n\;{\rm{is}}\;{\rm{odd}}\\ \mathrm{sin}(\beta \theta ),& n\;{\rm{is}}\;{\rm{even}}\end{array}\right. $$ (11) 式中, dn为待定系数, k = ω/c,
$ \beta $ 为n/2的整数部分,${\rm{H}}_\beta ^{(2)}( \cdot )$ 表示$ \beta $ 阶第二类汉开尔函数.式(10)两边同时对坐标r求偏导数, 并将结果与式(10)联立消去待定系数, 整理后可得到人工边界r = R处的频域动力刚度关系为
$$ F = \sum\limits_{n = 1}^N {{F_n}{\phi _n}} $$ (12) $$ {F_n} = - {S_n}{P_n} $$ (13) $$ {S_n} = - {S_0}\frac{{\bar \omega {\rm{H}}{{_\beta ^{(2)}}^\prime }\left( {\bar \omega } \right)}}{{{\rm{H}}_\beta ^{(2)}\left( {\bar \omega } \right)}} $$ (14) $$ {P_n} = \int_0^{2\text{π} } {P{\phi _n}{\text{d}}\theta } $$ (15) 式中, N为选取模态阶数; n = 1时
$ \delta = 2 $ , n > 1时$ \delta = 1 $ ;${\rm{H}}_{\beta} ^{\left( 2 \right)}{'}\left( x \right) = {\text{d}}{\rm{H}}_{\beta} ^{\left( 2 \right)}\left( x \right)/{\text{d}}x$ ;${S_0} = 1/(\delta \text{π})$ ;$ \bar \omega = kR $ ; Sn表示人工边界动力刚度系数;$F = r\dfrac{{\partial P}}{{\partial r}}$ .2.2 频域动力刚度的连分式
根据Li等[36]提出的连分式近似方法, 动力刚度系数式(14)的连分式形式为
$$ \frac{{{S_n}}}{{{S_0}}} = \left(\frac{1}{2} + {\text{i}}\bar \omega \right) - \dfrac{{{{(1/2)}^2} - {\beta ^2}}}{{2(1 + {\text{i}}\bar \omega ) - \dfrac{{{{(3/2)}^2} - {\beta ^2}}}{{2(2 + {\text{i}}\bar \omega ) - \dfrac{{{{(5/2)}^2} - {\beta ^2}}}{{2(3 + {\text{i}}\bar \omega ) - \ddots }}}}}} $$ (16) 通过引入辅助变量
$ \hat S_n^{} $ , 式(16)的显式表达递推关系为[36]$$ \frac{{{S_n}}}{{{S_0}}} = \frac{1}{2} + {{\rm{i}}} \bar \omega - \left[ {{{\left( {1/2} \right)}^2} - {\beta ^2}} \right]\hat S_{n,1}^{ - 1} $$ (17) $$ {\widehat{S}}_{n\text{, }j} = 2(j + \mathrm{i}\overline{\omega })-\left[{\left(j + 1/2\right)}^{2}-{\beta }^{2}\right]{\widehat{S}}_{b,j + 1}^{-1} $$ (18) $$ \hat S_{n,J{\text{ + 1}}}^{ - 1} = 0 $$ (19) 式中,
$ j = 1, 2\cdots , J $ ; J表示连分式的阶数.2.3 连分式的时域实现
将式(17) ~ 式(19)代入式(13), 并引入整理如下辅助变量方程
$$ {\hat{S}}_{n\text{, }j}{P}_{n,j} = {P}_{n,j-1} $$ (20) $$ {P_{n,0}} = {P_n} $$ (21) 进一步模态的动力刚度关系可表示为如下形式
$$ \frac{1}{2}{P_n} + {\rm{i}}{\bar \omega _m}{P_n} - \left[ {{{(1/2)}^2} - {\beta ^2}} \right]{P_{n,1}} = {F_n} $$ (22) $$ - {P_{n,j - 1}} + 2(1 + {\rm{i}}\bar \omega ){P_{n,j}} - \left[ {{{\left( {j + 1/2} \right)}^2} - {\beta ^2}} \right]{P_{n,j{\text{ + }}1}} = 0 $$ (23) $$ {P_{n,J{\text{ + }}1}} = 0 $$ (24) 将式(22) ~ 式(24)傅里叶逆变换到时域可得
$$ \frac{1}{2}{p_n} + {\dot p_n} - \left[ {{{(1/2)}^2} - {\beta ^2}} \right]{p_{n,1}} = {f_n} $$ (25) $$ - {p_{n,j - 1}} + 2{p_{n,j}} + 2{\dot p_{n,j}} - \left[ {{{\left( {j + 1/2} \right)}^2} - {\beta ^2}} \right]{p_{n,j{\text{ + }}1}} = 0 $$ (26) $$ {p_{n,J{\text{ + }}1}} = 0 $$ (27) 写成矩阵形式为
$$ {{\boldsymbol{K}}_n}{{\boldsymbol{p}}_n} + {{\boldsymbol{C}}_n}{{\boldsymbol{\dot p}}_n} = {{\boldsymbol{f}}_n} $$ (28) $$ {{\boldsymbol{u}}_n} = {\left\{ {\begin{array}{*{20}{c}} {{p_n}}&\left| {} \right. & {\begin{array}{*{20}{c}} {{p_{n,1}}}&{{p_{n,2}}}& \cdots &{{p_{n,J}}} \end{array}} \end{array}} \right\}^{\rm{T}}} $$ (29) $$ {{\boldsymbol{f}}_n} = {\left\{ {\begin{array}{*{20}{c}} {{f_n}}&\left| {} \right. & {\begin{array}{*{20}{c}} 0& \cdots &0 \end{array}} \end{array}} \right\}^{\rm{T}}} $$ (30) $$\begin{split} &{{\boldsymbol{K}}_n}=\left[\begin{array}{c:c} {K_{n 0}} & {\boldsymbol{K}}_{n1}\\ \hdashline {\boldsymbol{K}}_{n2} & {\boldsymbol{K}}_{n3} \end{array}\right]=\\ &\qquad {S_0}\left[\begin{array}{c:ccccc} 1/2 & k_{1} & 0 & \cdots & 0 & 0 \\ \hdashline -1 & 2 & k_{2} & \cdots & 0 & 0 \\ 0 & -1 & 2 & \cdots & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 & \cdots & 2 & {k_{j}} \\ 0 & 0 & 0 & \cdots & -1 & 2 \end{array}\right]_{(J+1) \times(J+1)} \end{split}$$ (31) $$ \begin{split} & \boldsymbol{C}_n=\left[\begin{array}{c:c} C_{n 0} & {\boldsymbol{C}}_{n 1} \\ \hdashline {\boldsymbol{C}}_{n 2} & {\boldsymbol{C}}_{n 3} \end{array}\right]= \\ &\qquad S_0\left[\begin{array}{c:ccccc} 1 & 0 & 0 & \cdots & 0 & 0 \\ \hdashline 0 & 2 & 0 & \cdots & 0 & 0 \\ 0 & 0 & 2 & \cdots & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 & \cdots & 2 & 0 \\ 0 & 0 & 0 & \cdots & 0 & 2 \end{array}\right]_{(J+1) \times(J+1)} \end{split} $$ (32) 式中,
$ {k_j} = - [{(j - 1/2)^2} - {\beta ^2}] $ .2.4 时域人工边界与内域的耦合
有限元方程(2)可以重新表示为
$$ \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{M}}_{ii}}}&{{{\boldsymbol{M}}_{ib}}} \\ {{{\boldsymbol{M}}_{bi}}}&{{{\boldsymbol{M}}_{bb}}} \end{array}} \right]\left\{ {\begin{array}{*{20}{c}} {{{{\boldsymbol{\ddot p}}}_i}} \\ {{{{\boldsymbol{\ddot p}}}_b}} \end{array}} \right\} + \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{K}}_{ii}}}&{{{\boldsymbol{K}}_{ib}}} \\ {{{\boldsymbol{K}}_{bi}}}&{{{\boldsymbol{K}}_{bb}}} \end{array}} \right]\left\{ {\begin{array}{*{20}{c}} {{{\boldsymbol{p}}_i}} \\ {{{\boldsymbol{p}}_b}} \end{array}} \right\} = \left\{ {\begin{array}{*{20}{c}} {{{\boldsymbol{F}}_i}} \\ {{{\boldsymbol{F}}_b}} \end{array}} \right\} $$ (33) 式中, 下标b和i分别表示边界节点和除边界外其他节点;
$ {{\boldsymbol{F}}_b} $ 的表达式为$$ {{\boldsymbol{F}}_b} = - \int_0^{2\text{π} } {{\boldsymbol{N}}_b^{\text{T}}{f_{}}{\text{d}}\theta } $$ (34) 式中f表示F的时域形式, Nb表示与人工边界结点相关的全局形函数行向量.
边界压力和模态函数可以分别离散为
$$\qquad\qquad\qquad p = {\boldsymbol{N}}{{\boldsymbol{p}}_b} $$ (35) $$\qquad\qquad\qquad {{\boldsymbol{\varphi }}_n} = {{\boldsymbol{N}}_b}{{\boldsymbol{T}}_n} $$ (36) 式中,
$ {{\boldsymbol{\varphi }}_n} $ 为边界结点处模态列向量.将式(35)和式(36)代入式(15)的时域形式, 之后将结果代入式(32)的时域形式和式(28), 并将得到的结果代入式(33)和式(34), 整理可得到人工边界条件与内域耦合的有限元方程为
$$ \begin{split} & \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{M}}_i}}&{\boldsymbol{0}}&{\boldsymbol{0}} \\ {\boldsymbol{0}}&{{{\boldsymbol{M}}_i}}&{\boldsymbol{0}} \\ {\boldsymbol{0}}&{\boldsymbol{0}}&{\boldsymbol{0}} \end{array}} \right]\left\{ {\begin{array}{*{20}{c}} {{{{\boldsymbol{\ddot p}}}_i}} \\ {{{{\boldsymbol{\ddot p}}}_b}} \\ {{{{\boldsymbol{\ddot p}}}_c}} \end{array}} \right\} + \left[ {\begin{array}{*{20}{c}} {\boldsymbol{0}}&{\boldsymbol{0}}&{\boldsymbol{0}} \\ {\boldsymbol{0}}&{{\boldsymbol{C}}_b^\infty }&{{\boldsymbol{C}}_{bc}^\infty } \\ {\boldsymbol{0}}&{{\boldsymbol{C}}_{cb}^\infty }&{{\boldsymbol{C}}_c^\infty } \end{array}} \right]\left\{ {\begin{array}{*{20}{c}} {{{{\boldsymbol{\dot p}}}_i}} \\ {{{{\boldsymbol{\dot p}}}_b}} \\ {{{{\boldsymbol{\dot p}}}_c}} \end{array}} \right\}+ \\ &\qquad \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{K}}_i}}&{{{\boldsymbol{K}}_{ib}}}&{\boldsymbol{0}} \\ {{{\boldsymbol{K}}_{bi}}}&{{{\boldsymbol{K}}_b} + {\boldsymbol{K}}_b^\infty }&{{\boldsymbol{K}}_{bc}^\infty } \\ {\boldsymbol{0}}&{{\boldsymbol{K}}_{cb}^\infty }&{{\boldsymbol{K}}_c^\infty } \end{array}} \right]\left\{ {\begin{array}{*{20}{c}} {{{\boldsymbol{p}}_i}} \\ {{{\boldsymbol{p}}_b}} \\ {{{\boldsymbol{p}}_c}} \end{array}} \right\} = \left\{ {\begin{array}{*{20}{c}} {{{\boldsymbol{F}}_i}} \\ {\boldsymbol{0}} \\ {\boldsymbol{0}} \end{array}} \right\} \end{split} $$ (37) $$ {\boldsymbol{p}}_c^{\rm{T}} = \left\{ {\begin{array}{*{20}{c}} {{\boldsymbol{p}}_1^{\infty {\rm{T}}}}&{{\boldsymbol{p}}_2^{\infty {\rm{T}}}}& \cdots &{{\boldsymbol{p}}_N^{\infty {\rm{T}}}} \end{array}} \right\} $$ (38) $$ {\boldsymbol{p}}_n^{\infty {\rm{T}}} = \left\{ {\begin{array}{*{20}{c}} {{p_{n,1}}}&{{p_{n,2}}}& \cdots &{{p_{n,J}}} \end{array}} \right\} $$ (39) $$ \begin{split} & {\left[\begin{array}{c:c} \boldsymbol{K}_b^{\infty} & \boldsymbol{K}_{b c}^{\infty} \\ \hdashline \boldsymbol{K}_{c b}^{\infty} & \boldsymbol{K}_c^{\infty} \end{array}\right]=} \\ &\qquad {\left[\begin{array}{c:ccccc} \boldsymbol{K}_b^{\infty} & \boldsymbol{K}_{b 1}^{\infty} & \boldsymbol{K}_{b 2}^{\infty} & \cdots & \boldsymbol{K}_{b(N-1)}^{\infty} & \boldsymbol{K}_{b N}^{\infty} \\ \hdashline \boldsymbol{K}_{1 b}^{\infty} & \boldsymbol{K}_1^{\infty} & \boldsymbol{0} & \cdots & \boldsymbol{0} & \boldsymbol{0} \\ \boldsymbol{K}_{2 b}^{\infty} & \boldsymbol{0} & \boldsymbol{K}_2^{\infty} & \cdots & \boldsymbol{0} & \boldsymbol{0} \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ \boldsymbol{K}_{(N-1) b}^{\infty} & \boldsymbol{0} & \boldsymbol{0} & \cdots & \boldsymbol{K}_{N-1}^{\infty} & \boldsymbol{0} \\ \boldsymbol{K}_{N b}^{\infty} & \boldsymbol{0} & \boldsymbol{0} & \cdots & \boldsymbol{0} & \boldsymbol{K}_N^{\infty} \end{array}\right]} \end{split} $$ (40) $$ {\boldsymbol{K}}_b^\infty = {\boldsymbol{W}}\sum\limits_{n = 1}^N {\left( {{K_{n0}}{{\boldsymbol{\varphi }}_n}{\boldsymbol{\varphi }}_n^{{\rm{T}}} } \right)} {{\boldsymbol{W}}_b} $$ (41) $$ {\boldsymbol{K}}_{bn}^\infty = {K_{n1}}{{\boldsymbol{W}}_b}{{\boldsymbol{I}}_n} $$ (42) $$ {\boldsymbol{K}}_{nb}^\infty = {K_{n2}}{{\boldsymbol{W}}_b}{{\boldsymbol{I}}_n}^{\rm{T}} $$ (43) $$ {\boldsymbol{K}}_n^\infty = {{\boldsymbol{K}}_{n3}} $$ (44) $$ \begin{split} & {\left[\begin{array}{c:c} \boldsymbol{C}_b^{\infty} & \boldsymbol{C}_{b c}^{\infty} \\ \hdashline \boldsymbol{C}_{c b}^{\infty} & \boldsymbol{C}_c^{\infty} \end{array}\right]=} \\ &\qquad {\left[\begin{array}{c:ccccc} \boldsymbol{C}_b^{\infty} & \boldsymbol{C}_{b 1}^{\infty} & \boldsymbol{C}_{b 2}^{\infty} & \cdots & \boldsymbol{C}_{b(N-1)}^{\infty} & \boldsymbol{C}_{b N}^{\infty} \\ \hdashline \boldsymbol{C}_{1 b}^{\infty} & \boldsymbol{C}_1^{\infty} & \boldsymbol{0} & \cdots & \boldsymbol{0} & \boldsymbol{0} \\ \boldsymbol{C}_{2 b}^{\infty} & \boldsymbol{0} & \boldsymbol{C}_2^{\infty} & \cdots & \boldsymbol{0} & \boldsymbol{0} \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ \boldsymbol{C}_{(N-1) {b}}^{\infty} & \boldsymbol{0} & \boldsymbol{0} & \cdots & \boldsymbol{C}_{N-1}^{\infty} & \boldsymbol{0} \\ \boldsymbol{C}_{N b}^{\infty} & \boldsymbol{0} & \boldsymbol{0} & \cdots & \boldsymbol{0} & \boldsymbol{C}_N^{\infty} \end{array}\right]}\end{split} $$ (45) $$ {\boldsymbol{C}}_b^\infty = {{\boldsymbol{W}}_b}\sum\limits_{n = 1}^N {\left( {{C_{n0}}{{\boldsymbol{\varphi }}_n}{\boldsymbol{\varphi }}_n^{{\rm{T}}} } \right)} {{\boldsymbol{W}}_b} $$ (46) $$ {\boldsymbol{C}}_{bn}^\infty = {C_{n1}}{{\boldsymbol{W}}_b}{{\boldsymbol{I}}_n} $$ (47) $$ {\boldsymbol{C}}_{nb}^\infty = {C_{n2}}{{\boldsymbol{W}}_b}{{\boldsymbol{I}}_n} $$ (48) $$ {\boldsymbol{C}}_n^\infty = {{\boldsymbol{C}}_{n3}} $$ (49) $$ {{\boldsymbol{W}}_b} = \int_0^{2\text{π} } {{\boldsymbol{N}}_b^{\text{T}}{{\boldsymbol{N}}_b}{\text{d}}\theta } $$ (50) 式中, In表示NB × J阶的矩阵, 第一列
$ {{\boldsymbol{\varphi }}_n} $ 而其余元素为零; NB表示边界节点数目;$ {K_{n1}} $ 和$ {K_{n2}} $ 分别表示$ {{\boldsymbol{K}}_{n1}} $ 和$ {{\boldsymbol{K}}_{n2}} $ 的第一个元素;$ {C_{n1}} $ 和$ {C_{n2}} $ 分别表示$ {{\boldsymbol{C}}_{n1}} $ 和$ {{\boldsymbol{C}}_{n2}} $ 的第一个元素.式(37)可通过数值积分方法Newmark-β法进行求解.
3. 数值模型算例
3.1 方法验证
本文的计算模型采用二维圆形区域和圆形结构. 如图3所示, 爆炸点位于结构左侧外侧, a表示结构半径; 结构上三个测点位置A, B, C; (xc, yc)表示结构中心坐标; (xo, yo)表示爆炸点坐标; (xm, ym)表示结构表面测点坐标; W表示炸药装药质量; R表示爆炸距离.
为验证本文提出的水下爆炸作用下爆源子结构输入方法的准确性, 选取如图4所示的自由场模型进行计算分析对比验证. 在ABAQUS中建立自由场模型, 半径r = 0.6 m, 网格为0.01 m,建立三圈子结构ABC, 在爆源位置D处加一个如图5所示的脉冲载荷, 从而计算得出整个自由场区域内的响应, 该响应作为参考解以验证本文方法的有效性, 本文选取爆源外一点E为观测点. 进一步利用本文计算方法计算观测E点的压力时程, 图6和图7为本文方法与参考解的对比. 由图中可以看出, 本文方法的计算结果与参考解吻合很好.
如图8所示, r表示水域半径的大小, 通过选取两个不同大小水域的模型, 比较结构上一点的压力时程, 由图可以看出, 本文方法不受模型水域大小的影响, 能提高计算效率, 同时保证计算精度.
在水下爆炸的数值模拟中网格大小对模拟结果有很大的影响, 网格越小计算精度越高, 但对计算机硬件要求较大, 计算时间增长, 因此需要找到一个适宜的网格尺寸, 既能保证计算效率又能保证计算精度. 图9表示, 在同一工况下, 结构上一点测点压力时程峰值随网格尺寸的变化, 可以看出随着网格尺寸的变小压力峰值减少, 网格尺寸为0.008 m与0.01 m时压力峰值变化幅度小于3%, 由于网格尺寸小于0.01 m时计算量较大, 0.01 m网格也具有较高的精度, 因此采用0.01 m的网格大小进行计算.
3.2 方法应用
水下爆炸的过程分为炸药的炮轰、冲击波的传播、气泡脉动三个阶段, 结构表面测点A, B, C的全过程总压力时程曲线如图10所示, 其中, W = 12 g, a = 0.2 m, (xc, yc) = (0.2 m, 0), (xo, yo) = (−0.25 m, 0). 计算全过程压力时程需要计算的时间尺度较长, 由图10可以看出, 冲击波阶段和气泡脉动阶段测点处压力峰值从结构左侧传播到结构右侧的过程中压力值均逐渐衰减. 因此后文只计算分析冲击波阶段的压力时程.
结构表面不同测点的压力时程曲线如图11所示, 其中, W = 12 kg, a = 0.2 m, (xc, yc) = (0.2 m, 0), (xo, yo) = (−0.25 m, 0). 可以看出, 爆炸波从结构左侧传播到结构右侧的过程中压力值逐渐衰减.
本节主要研究的是爆炸波传播过程中由结构引起的散射效应. 图12为结构表面三个测点处自由场模型和有结构模型的压力时程曲线, 自由场模型测点处压力时程表示入射波的压力, 有结构模型测点处的压力表示总压力. 其中a = 0.2 m, W = 12 kg, (xc, yc) = (0.2 m, 0), (xo, yo) = (−0.25 m, 0). 由图可以看出, 三个测点中θ = 0处的入射压力大于总压力, θ = π处的入射压力小于总压, θ = 0.5π处的入射压力接近总压. 这主要是因为结构左侧散射压力方向与入射压力方向相同, 而结构右侧散射压力方向与入射压力方向相反. 图13表示随时间结构一周入射压力和总压力在x方向上的压力(Fx)积分时程,
${F_x} = - \displaystyle\int_0^{2\text{π} } {pa\cos \theta {\rm{d}}\theta }$ ,${F_x}_i = - \displaystyle\int_0^{2\text{π} } {{p_i}a\cos \theta {\rm{d}}\theta }$ 其中Fx 和Fxi分别表示总压力和入射压力在x方向上的积分, p和pi分别表示总压力和入射压力. 由图可以看出, 结构一周总压力在x方向上的积分大于入射压力在x方向上的积分, 这说明在水下爆炸载荷的作用下, 散射波压力会增强结构表面的爆炸波载荷.为了进一步评估结构对水下爆炸波的散射效应, 引入了无量纲参数
$ \delta = {{({F_x} - {F_{xi}})} \mathord{\left/ {\vphantom {{({F_x} - {F_{xi}})} {{F_x}}}} \right. } {{F_x}}} $ , 参数δ随圆形结构半径的变化如图14所示; δ随炸药量的变化如图15所示. 其中W = 12 kg, R = 0.25 m. 可以看出, 结构对爆炸波的散射效应随半径的增大而增大, 随炸药装药质量的增大而减小.4. 结论
本文基于内部子结构的地震波动输入方法和声学波动方程, 提出了一种模拟水下爆炸载荷作用下水−结构相互作用的爆源子结构输入方法. 该方法采用AUTODYN软件中一维自由场模型计算水下爆炸作用下爆源子结构区域自由波场压力时程, 并采用该压力时程进行水下爆炸波的等效载荷输入; 另外, 本模型基于连分式近似方法, 建立了一种高精度的时域人工边界条件, 该边界条件可设置于结构近处, 提高了计算效率. 最后, 基于本文方法研究了水下爆炸载荷作用下圆形结构周围散射波压力的变化规律, 结果表明: 在水下爆炸载荷作用下, 散射波压力会增大结构的爆炸载荷, 并且结构对爆炸波的散射效应随半径的增大而增大, 随炸药装药质量的增大而减小.
-
-
[1] 胡春红, 冯新, 李昕等. 水下爆炸作用下结构响应的数值计算研究综述. 工程爆破2007, 13(1): 28-34 Hu Chunhong, Feng Xin, Li Xin, Zhou Jing. Revier of numerical simulation of structural. Engineering Blasting, 2007, 13(1): 28-34(in Chinese))
[2] Li GL, Shi DY, Wang LF, et al. Measurement technology of underwater explosion load: A review. Ocean Engineering, 2022(254): 111383
[3] Ghoshal R, Mitra N. Underwater oblique shock wave reflection from submerged hydraulic structures. Ocean Engineering, 2020(209): 107324
[4] Zhang ZF, Wang LK, Silberschmidt VV. Damage response of steel plate to underwater explosion: Effect of shaped charge liner, International Journal of Impact Engineering, 2017(103): 38-49.
[5] 闫秋实, 张志杰, 王丕光等. 水下爆炸荷载作用下圆柱结构反射压力解析计算方法研究. 工程力学, 2022(39): 247-256 (Yan Qiushi, Zhang Zhijie, Wang Piguang, et al. Reach on Analytical Method of Circular Cylindrical. Engineering Mechanics, 2022(39): 247-256 (in Chinese) [6] Cole RH. Underwater Explosion. New Jersey: Princeton University Press, 1948: 118-127.
[7] Bobrovskii SV, Gogolev VM, Zamyshlyaev BV, et al. Effect of thermal decomposition on the scabbing velocity due to strong shock waves in solid media. Soviet Mining Science, 1976, 12(3): 273-279 doi: 10.1007/BF02594869
[8] Zhuang TS, Wang MY, Wu J, et al. Experimental investigation on dynamic response and damage models of circular RC columns subjected to underwater explosions. Defence Technology, 2020, 16(4): 856-875 doi: 10.1016/j.dt.2019.10.015
[9] 庄铁栓, 伍俊, 许文轩等. 水中爆炸冲击波在高桩圆柱结构上的分布规律试验研究. 中国测试, 2018, 44(10): 60-66 (Zhuang Tieshuan, Wu Jun, Xu Wenxuan, et al. Experimental investigation of distributing rule on high column stake for shock wave in underwater explosion. China Measurement &Test, 2018, 44(10): 60-66 (in Chinese) [10] Yang GD, Wang GH, Lu WB, et al. Experimental and numerical study of damage characteristics of RC slabs subjected to air and underwater contact explosions. Marine Structures, 2019(66): 242-257
[11] Geers TL, Doubly asymptotic approximations for transient motions of submerged structures. Journal of the Acoustical Society of America, 1978, 64(5): 1500-1508.
[12] Hai L, Ren XD. Computational investigation on damage of reinforced concrete slab subjected to underwater explosion. Ocean Engineering, 2020(195): 106671
[13] 宫翔飞, 刘文韬, 张树道等. 水下爆炸近场峰值压力的数值模拟. 爆炸与冲击, 2019, 39(4): 8 Gong Xiangfei, Liu Wentao, Zhang Shudao, Numerical simulation of near-field peak pressure of underwater explosion. Explosion and Shock Waves, 2019, 39(4): 8 (in Chinese)
[14] Zhang ZF, Wang C, Zhang AM, et al. SPH-BEM simulation of underwater explosion and bubble dynamics near rigid wall, 2019, 62(7): 12. [15] 尹训强, 吕硕硕, 王桂营. 水下爆炸作用下近岸场地动态响应数值模拟. 科学技术与工程, 2022, 22(5): 2009-2015 (Yin Xunqiang, Lü Shuoshuo, Wang Guixuan. Numerical simulation on dynamic response nearshore site under underwater explosion. Science Technology and Engineering, 2022, 22(5): 2009-2015 (in Chinese) [16] Barras G, Souli M, Aquelet N, et al. Numerical simulation of underwater explosions using an ALE method. The pulsating bubble phenomena. Ocean Engineering, 2012(41): 53-66
[17] 肖秋平, 陈网桦, 贾宪振等. 基于 AUTODYN的水下爆炸冲击波模拟研究. 舰船科学技术, 2009, 31(2): 6 (Xiao Qiuping, Chen Wangye, Jia Xianzhen, et al. Numerical study of underwater explosion shockwave based on AUTODYN. Ships Science and Technology, 2009, 31(2): 6 (in Chinese) [18] 闫秋实, 宁素瑜, 杜修力等. 水中近场爆炸作用下钢筋混凝土桩毁伤效应研究. 北京工业大学学报, 2019, 45(2): 153-159 Yan Qiushi, Ning Suyu, Du Xiuli, et al. Damage effect for a typical reinforced concrete pile under the nearfield explosion in water. Journal of Beijing University of Technology. 2019, 45(2): 153-159 (in Chinese))
[19] Hu YD, Zhao YF, Zhao M, et al. Three phases fluid-structure interactive simulations of the deepsea ceramic sphere’s failure and underwater implosion. Ocean Engineering, 2022(245): 110494
[20] Peng YX, Zhang AM, Ming FR. Numerical simulation of structural damage subjected to the near-field underwater explosion based on SPH and RKPM. Ocean Engineering, 2021, 222: 108576
[21] 张浩, 赵梓斌, 李上明. 基于总场公式的水下冲击结构响应分析方法. 工程力学, 2021, 38(11): 220-228, 247 (Zhang Hao, Zhao Zibin, Li Shangming. Analysis method of the structural response to underwater shock based on total field formulation. Engineering Mechanics, 2021, 38(11): 220-228, 247 (in Chinese) [22] 李程辉. 基于声固耦合法的水下爆炸冲击波在复杂边界附近空化特性研究[硕士论文]. 哈尔滨工程大学, 2018 Li Chenghui. Investigation on cavitation of underwater explosion shockwave near compled CEL method [Master Thesis]. Harbin Engineering University, 2018(in Chinese))
[23] Liu L, Zhang JQ, Song CM, et al. An automatic approach for the acoustic analysis of three-dimensional bounded and unbounded domains by scaled boundary finite element method. International Journal of Mechanical Sciences, 2019(151): 563-581
[24] Wang PG, Lu RR, Yan QS, et al. The internal substructure method for shock wave input in 2D fluid-structure interaction analysis with unbounded domain using doubly-asymptotic ABC. Mechanics of Advanced Materials and Structures, 2022
[25] 姚熊亮, 张阿漫, 许维军. 声固耦合方法在舰船水下爆炸中的应用. 哈尔滨工程大学学报, 2005, 6: 707-712 (Yao Xiongliang, Zhang Aman, Xu Weijun. Application of coupled acoustic_structural analysis to warship underwater explosion. Journal of Harbin Engineering University, 2005, 6: 707-712 (in Chinese) [26] 姚熊亮, 张阿漫, 许维军等. 基于ABAQUS软件的舰船水下爆炸研究. 哈尔滨工程大学学报, 2006(1): 37-41 (Yao Xiongliang, Zhang Aman, Xu Weijun. Research on warship underwater explosion with ABAQUS software. Journal of Harbin Engineering University, 2006(1): 37-41 (in Chinese) [27] 金乾坤, 丁刚毅. 舰船舱段水下冲击响应有限元分析//第十二届现代数学和力学会议论文集, 2010 Jin Qiankun, Ding Gangyi, A finite element analysis of ship sections subjected to underwater shock//Proceedings of the 12th Conference on Modern Mathematics and Mechanics, 2010(in Chinese))
[28] 刘晶波, 李述涛等. 基于爆源子结构的爆炸问题多尺度分析方法. 振动与冲击, 2021, 40(20): 63-72.3 (Liu Jingbo, Li Shutao, et. al A multiscale analysis method for explosion problems based on the substructure of explosion source. Journal of Vibration and Shock, 2021, 40(20): 63-72.3 (in Chinese) [29] 刘晶波, 吕彦东 结构–地基动力相互作用问题分析的一种直接方法. 土木工程学报, 1998, 31(3): 55-64 Liu Jingbo, Lü Yandong. A direct method for analysis of dynamic soil-structure interaction. China Civil Engineering Journal. 1998, 31(3): 55-64 (in Chinese)
[30] 刘晶波, 谭辉, 宝鑫等. 土–结构动力相互作用分析中基于人工边界子结构的地震波动输入方法. 力学学报, 2018, 50(1): 32-43 (Liu Jingbo, Tan Hui, Bao Xin, et al. The seismic wave input method for soil-structure dynamic interaction analysis based on the substructure of artificial boundaries. Chinese Journal of Theoretical and Applied Mechanics, 2018, 50(1): 32-43 (in Chinese) [31] 李述涛, 刘晶波, 宝鑫等. 人工边界子结构地震动输入方法在 ABAQUS 中的实现. 自然灾害学报, 2020, 49(4): 133-141 Li Shutao, Liu Jingbo, Bao Xin, et al. Implementation for seismic wave input method based on the artificial boundary substructure in ABAQUS. Journal of Natural Disasters. 2020, 49(4): 133-141(in Chinese))
[32] 刘晶波, 宝鑫, 谭辉等. 土-结构动力相互作用分析中基于内部子结构的地震波动输入方法. 土木工程学报, 2020, 53(8): 87-96 (Liu Jingbo, Bao Xin, Tan Hui, et al. Seismic wave input method for soil-structure dynamic interaction analysis based on internal substructure. China Civil Engineering Journal, 2020, 53(8): 87-96 (in Chinese) [33] Liu JB, Xin B, Dong Y, et al. The internal substructure method for seismic wave input in 3D dynamic soil-structure interaction analysis. Soil Dynamics and Earthquake Engineering, 2019(127): 261-267
[34] 韩昌瑞, 孙伟, 张玉华. 有限元分析与方法. 北京: 科学出版社, 2022 Han Changrui, Sun Wei, Zhang Yuhua. Finite Element Analysis and Methods. Beijing: Science Press, 2022 (in Chinese))
[35] 刘晶波, 王振宇, 杜修力等. 波动问题中的三维时域黏弹性人工边界. 工程力学, 2005, 22(6): 46-51 Liu Jingbo, Wang Zhenyu, Du Xiuli, et al. Three-dimensional visco-elastic artificial boundaries in ime domain for wave motion problems. Engineering Mechanics, 2005, 22(6): 46-51 (in Chinese)