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

李述涛 刘晶波 宝 鑫,2)

∗(清华大学土木工程系，北京 100084)

†(军事科学院国防工程研究院，北京 100036)

引言

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

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

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

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

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

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

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

其中，G为内部介质的剪切模量，ρ 为内部介质的质量密度，CS和CP分别为内部介质的S 波和P 波波速，R为散射波源至边界节点的距离，αT与αN分别为切向与法向黏弹性人工边界参数，对于二维黏弹性人工边界，其推荐值分别为0.5 和1.

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

图1 二维黏弹性人工边界单元示意图Fig.1 Diagrams of 2D viscoelastic artificial boundary element

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

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

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

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

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

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

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

图2 二维模型中的两种节点子系统Fig.2 Two kinds of nodal subsystems in 2D model

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

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

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

2.1 角点子系统稳定性条件

角点子系统如图3 所示.设内部介质单元的弹性模量为E，剪切模量为G，密度为ρ，泊松比为µ，且无阻尼；黏弹性人工边界单元的密度为，弹性模量为，阻尼系数为，泊松比为0，单元边长均为L.

图3 角点子系统Fig.3 Corner subsystem

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

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

其中

根据稳定性判别条件ρ(A1)1，得出角点子系统稳定性条件的解析表达式为

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

其中，γ1为角点子系统稳定性参数，整理后如式(11)所示.

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

2.2 侧边子系统稳定性条件

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

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

图4 侧边子系统Fig.4 Edge subsystem

其中

计算传递矩阵A2的特征值，并根据稳定性判别条件ρ(A2)1，得出侧边子系统稳定性条件的解析表达式为

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

其中γ2为侧边子系统稳定性参数，整理后如式(18)所示.

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

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

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

其中，γ0=1.

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

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

图5 角点子系统稳定性系数随密度比变化情况Fig.5 The variation of corner subsystem stability coefficient with density ratio

图6 侧边子系统稳定性系数随密度比变化情况Fig.6 The variation of edge subsystem stability coefficient with density ratio

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

图7 整体系统数值稳定区域Fig.7 The region of numerical stability in the overall system

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

3.2 距离R 对稳定性条件的影响

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

图8 稳定性系数随R/L的变化情况Fig.8 The variation of stability coefficient with R/L

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

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

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

改进的黏弹性人工边界单元参数设置如下：设内部介质剪切模量为G，质量密度为ρ，剪切波速和压缩波速分别为CS和CP；内部介质单元厚度和边界单元厚度均为L，边界节点至散射波源的距离为R，边界单元质量密度为；αT与αN分别为切向与法向黏弹性人工边界参数，二维模型建议αT取0.5，αN取1.0.稳定性条件改善后，二维黏弹性人工边界单元质量密度、等效刚度、泊松比和等效阻尼分别为

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

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

4.1 均匀半空间模型

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

图9 均匀半空间模型Fig.9 Homogeneous half space model

图10 动力载荷时程曲线Fig.10 Time history of the dynamic load

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

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

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

图11 均匀介质模型中4 个观测点y 方向位移比较Fig.11 Comparison of displacements in y direction of four nodes in the homogeneous model

4.2 成层半空间模型

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

图12 成层半空间模型Fig.12 Layered half space model

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

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

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

图13 成层介质模型中四个观测点y 方向位移比较Fig.13 Comparison of displacements in y direction of four nodes in the layered model

5 结论

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

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

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

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

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