采用粘弹性人工边界单元时显式算法稳定性分析

李述涛，刘晶波，宝 鑫，陶西贵，陈一村，肖 兰，贾艺凡

(1. 清华大学土木工程系，北京 100084；2. 军事科学院国防工程研究院，北京 100036)

采用数值方法计算土-结构地震反应或近场波动问题时，需要从半无限介质中切取出有限的近场区域进行计算，同时需在截断边界处设置人工边界来模拟外行散射波向无穷远的辐射效应[1 − 4]。较为常见的粘弹性人工边界[5]计算精度较高但前处理操作较为复杂，需遍历所有边界节点施加弹簧-阻尼器系统，在此基础上发展了粘弹性人工边界单元技术[6 − 7]，它是沿边界法线向外延伸的一层单元，通过赋予该层单元等效物理参数即可模拟粘弹性人工边界对于外行散射波的吸收作用。由于该方法具有良好的计算精度和鲁棒性，且前处理过程简单，在实际工程计算中得到较多的应用[8 − 12]。

随着人工边界技术的不断发展，针对人工边界条件的稳定性研究逐渐得到重视[13 − 14]，目前针对此类问题的主要研究方法包括：Trefethen 等[15]提出的GKS 定理—基于偏微分方程初边值问题有限差分格式的稳定性理论，给出初边值问题多层线性差分格式稳定的充要条件；Liao 等[16]研究了离散模型中稳态波动的完备解，给出了多次透射边界的稳定性条件；Kamel 等[17]和关慧敏等[18]通过传递矩阵谱半径分析积分格式的稳定性。理论研究表明，如果不全面考虑逐步积分过程中人工边界节点与内部节点运动方程的耦合效应，稳定性准则可能会失效。

由于粘弹性人工边界或粘弹性人工边界单元本身均为物理上稳定的系统，将其引入计算系统中，并不影响整体的物理稳定条件。但采用显式时域逐步积分方法进行计算时，受粘弹性人工边界单元粘性阻尼、刚度及几何尺寸的影响，整体计算模型的数值积分稳定性条件将发生改变，目前尚未给出明确而实用的含粘弹性人工边界条件影响的显式算法稳定性准则，影响了数值分析时稳定时间步长的判断和选取，进一步限制了粘弹性人工边界单元在显式动力分析中的应用。因此目前在引入粘弹性人工边界单元进行波动问题分析时，多采用隐式的无条件稳定的时域逐步积分算法，以避免显式算法带来的数值稳定性问题。隐式算法需要求解联立方程组，计算效率不高，并不适合大规模波动问题计算。随着显式时域逐步积分算法在工程结构和大范围场地分析问题中的广泛应用[19 − 22]，有必要对含有粘弹性人工边界单元的系统进行显式时域逐步积分算法稳定性研究。

本文考虑人工边界和内部单元节点的运动耦合效应，利用传递矩阵的谱半径分析时域逐步积分格式(中心差分)的稳定性，给出不同位置处局部节点系统显式时域逐步积分算法稳定条件的解析解，通过各子系统稳定性条件及内部介质稳定性条件的对比分析，给出考虑粘弹性人工边界条件影响的整体耦合系统显式时域逐步积分算法的稳定性条件。

1 粘弹性人工边界单元等效物理参数及尺寸设置

粘弹性人工边界单元是在模型截断边界处沿法线向外延伸的一层单元(见图1)，人工边界单元的厚度为h，质量密度为0，单元最外层节点固定。通过赋予单元等效物理参数来模拟粘弹性人工边界。二维粘弹性人工边界单元等效弹性模量、等效泊松比和等效阻尼系数分别为[2]：

图1 人工边界单元尺寸示意图Fig. 1 Artificial boundary elements size diagram

其中：设G为内部介质剪切模量； ρ为介质质量密度；CS和CP分别为介质S 波和P 波波速；R为散射波源至边界节点的距离； α=αN/αT， αT与 αN分别为切向与法向粘弹性人工边界参数，对于二维粘弹性人工边界，其推荐值分别为0.5 和1，根据式(1)，此时为0。

尽管粘弹性人工边界单元等效刚度矩阵的推导过程引入了边界单元厚度h远小于宽度L的假设，如图1 所示，刘晶波等[6]和谷音等[7]通过进一步的研究表明，边界单元厚度的鲁棒性良好，可灵活取值而对精度影响很小。由于显式算法的稳定条件受模型中的最小单元尺寸制约，在实际计算分析中，建议将边界单元的宽度设定为与内部单元尺寸一致，以排除边界单元尺寸对整体稳定性的影响，如图2 所示。此时粘弹性人工边界单元等效弹性模量、等效剪切模量可以写为：

图2 适用于显式算法的人工边界单元Fig. 2 Artificial boundary elements adapted to explicit algorithms

使用通用有限元软件对粘弹性人工边界单元赋予材料参数时，由于是各向同性介质，输入等效弹性模量即可。另外还包括式(1)计算得到的等效阻尼系数，以及泊松比(=0)和质量密度(=0)。由于大多数软件不支持将密度设为0，可采用近似0 的小数代替。

2 显式时域逐步积分算法稳定性分析方法

时域逐步积分算法稳定性分析的目的是获得时域逐步积分计算时满足稳定性要求的时间离散步长 ∆t。最大的稳定时间步长 ∆t与单元的尺寸和系统的物理性质有关。为此在进行稳定性分析时，通常采用均匀介质和均匀离散化网格模型，如图3 所示。

图3 均匀离散化网格模型Fig. 3 Uniform discretization mesh model

当实际力学模型介质非均匀，单元尺寸不相同时，可以根据最小尺寸的单元或综合考虑单元尺寸及其介质性质来确定整体计算模型的稳定离散时间步长 ∆t。算法的稳定性通常可采用以下几种分析方法：

① 不考虑自由边界和人工边界的影响，取无限计算区域模型，采用冯诺依曼方法进行稳定性分析。这一方法可以获得内部系统的数值计算稳定条件，但无法考虑含人工边界条件时系统的稳定性。

② 取实际的计算模型，考虑人工边界的影响，通过对整体系统时域逐步积分方法的传递矩阵进行特征值分析，以获得考虑耦合人工边界条件影响的整体系统的数值计算稳定性条件。但由于涉及到整体模型传递矩阵的建立和分析，这一方法在实际工作中通常是不可行的，特别是对大规模的计算问题。

③ 为对透射人工边界的数值计算稳定性进行有效分析，Kamel 等[17]和关慧敏等[18]提出了一种对含透射人工边界条件子系统的传递矩阵进行特征值分析，以获得波动问题中透射人工边界稳定性的分析方法。由于分析中采用的子系统为沿人工边界切向(宽度方向)上若干排节点和沿法向(深度方向)上若干排节点组成的节点系，子系统自由度较多，仅能靠数值方法对人工边界的稳定性进行分析和判断，未给出具有解析形式的更易于使用的人工边界稳定性准则。

从以上有限元离散模型数值积分稳定性准则的研究工作中可以发现以下两个特点：1)算法的稳定性与模型的截止频率(系统最高阶自振频率)有关。2)截止频率相应的振型呈现局部节点系相邻节点交错振动的形态。由以上两个特点可以判断，整体系统的截止频率可通过对局部子系统模型的分析得到，进一步可通过局部子系统模型的数值稳定性分析，获得整体系统的数值稳定性条件。

首先以一维剪切梁问题为例证明以上判断。图4 为一维剪切梁的离散模型及最高阶振型示意图，梁的剪切模量为G，密度为 ρ，横截面面积为A，单元长度为L。取两个子系统进行分析。子系统a：取振型两节点之间的子系统，两端施加约束，由于约束是施加于振型节点(振幅为零)之上，因此并不改变剪切梁有限元模型的自振频率，如图5(a)所示；子系统b：取两个单元，在两端施加约束，此时子系统的自振频率与原系统的截止频率不相同，如图5(b)所示。

图4 一维剪切梁振动示意图Fig. 4 One-dimensional shear beam vibration diagram

图5 一维局部子系统Fig. 5 One-dimensional local subsystem

根据子系统的刚度和质量，可以建立子系统a的单自由度运动方程，当采用中心差分法进行时域逐步积分计算时，一维子系统a 的数值稳定条件为：

其中，Tn为子系统的自振周期，将式(3)中第三式代入式(4)，得到子系统数值积分稳定性条件为：

式(5)即为一维连续介质波动有限元分析时中心差分算法的稳定性条件，可见采用子系统a 获得的局部稳定性条件即为连续介质整体模型的稳定性条件。

一维子系统b 的刚度kb、质量mb、自振频率ωb和中心差分稳定性条件分别为：

对比式(5)和式(6)，可以发现一维子系统b的稳定性条件比实际整体模型的稳定条件更为宽松，但仍然可以给出稳定时间积分步长 ∆t的一个上限估计。

图6 二维节点系统最高阶振型图Fig. 6 Highest-order modal pattern of two-dimensional nodes system

同样取两个子系统进行分析，二维问题局部子系统由4 个单元构成，如图7(a)和图7(b)所示。其中二维子系统a 的4 个角点为振型节点，位移为0，4 个边点的位移条件可由振型和4 边形等参元的性质确定，子系统a 对应的边界条件为：

图7 二维局部子系统Fig. 7 Two-dimensional local subsystem

二维子系统a 自振频率和中心差分稳定条件为：

式(9)即为二维波动问题有限元分析时中心差分算法的稳定性条件，可见同样可以由子系统a的局部稳定性分析获得二维有限元模型数值算法的整体稳定性条件。

二维子系统b 由4 个单元构成，在周边节点施加固定约束，子系统b 运动方程为：

同样，二维子系统b 给出了整体有限元模型稳定时间积分步长的上限估计。

一维和二维子系统分析的结果表明：1)若能较为准确地判断系统截止频率对应的振动模态(振型)，则可以利用最高振型的特点和振型节点的分布规律，从整体系统中分离出由阵型节点所包围的局部子系统，对振型节点施加约束条件，然后对该子系统进行稳定性分析，获得的局部子系统稳定性条件即为整体模型的数值稳定条件；2)若不能准确判断截止频率对应振动模态和相应振型节点的位置，也可选取一个能体现整体有限元模型特征的最小的子系统，对该子系统的边界节点施加约束，通过分析获得子系统的数值稳定性条件，从而获得整体系统稳定性判据中各物理量的关系式，且该条件是整体系统稳定性条件的一个逼近。再通过扩展数值算例进行修正，确定合适的系数，即可得到整体系统的数值积分稳定性条件。

当采用粘弹性人工边界单元时，内部系统的数值稳定条件已知，但难以确定人工边界区截止频率对应的振型，因而仅能采用类似于b 型的子系统进行数值计算的稳定性分析。

对非均质子系统进行稳定性分析时，可根据显式逐步积分算法格式，将系统运动方程写成如下形式：

积分格式的稳定性问题与外力向量Qp无关。如果满足下列两条件，则积分格式(12)是稳定[23]：条件1： ρ(A)≤1 ， ρ(A)是传递矩阵A的谱半径，即ρ(A)=max|λi|；条件2：如果A具有多重特征值，则该特征值的模小于1。因而可将子系统的稳定性分析归结为传递矩阵A的形成及谱半径计算，按上述条件即可建立稳定性准则。

3 人工边界子系统稳定性分析

采用粘弹性人工边界单元时，人工边界区的子系统有两种形式，一种在侧面或底面切取，如图8 所示，另一种在角点处切取，切取边界可设置固定边界(见图9)或自由边界(见图10)。

图8 侧边固定边界子系统Fig. 8 Side fixed boundary subsystem

图9 角点固定边界子系统Fig. 9 Corner fixed boundary subsystem

图10 角点自由边界子系统Fig. 10 Corner fixed boundary subsystem

3.1 侧边固定边界子系统稳定性分析

侧边固定边界子系统由两个粘弹性人工边界单元和两个内部介质单元组成，如图8 所示。设内部介质单元的弹性模量为E，剪切模量为G，密度为 ρ，泊松比为μ，且无阻尼；粘弹性人工边界单元弹性模量为，阻尼系数为，泊松比为0，密度为0，单元边长均为L。四节点正方形平面单元的集中质量矩阵、刚度矩阵、阻尼矩阵分别见文献[6]和文献[24]，按节点编号进行矩阵组装后得到子系统中5 号节点运动方程如下：

观察式(13)可知，x和y方向运动方程是解耦的且形式相同，因此可只对其中一个方向的运动方程进行稳定性分析。显式时域逐步积分格式(中心差分[25]如下：

式(19)即为侧边局部人工边界子系统稳定性条件的解析表达形式，观察可知该子系统的稳定性条件与内部介质压缩波速、泊松比、单元尺寸和模型大小有关。

3.2 角点固定边界子系统稳定性分析

角点处固定边界子系统中由3 个粘弹性人工边界单元和1 个内部介质单元组成，如图9 所示。将四节点正方形平面单元质量、阻尼和刚度矩阵按节点编号组装后得到5 号节点的运动方程如下：

由于该运动方程的坐标耦联，需对整体方程进行稳定性分析。将式(1)代入式(20)，然后根据式(14)将式(20)展开，写成如下传递矩阵的形式：

其中：

式(24)即为角点处局部系统稳定性条件的解析表达形式，观察可知影响角点局部人工边界子系统稳定性的参数与侧边相同，但各自贡献不同。

3.3 角点自由边界子系统稳定性分析

为进一步研究角点处自由边界子系统的稳定性，如图10 所示。研究对象为编号为1、2、4、5 的四个节点。子系统中既有单元内部节点，也包含边界处节点。将四节点正方形平面单元质量、阻尼和刚度矩阵按节点编号组装后，可得到子系统的刚度、质量和阻尼矩阵，均为8×8 阶对称矩阵。将整体运动方程按照式(14)中心差分格式展开，得到16×16 阶传递矩阵。该传递矩阵无法得出解析形式的稳定性条件，可代入参数计算数值解。考虑稳定性条件受系统截止频率影响，自由边界子系统截止频率要低于固定边界子系统，由此判断后者稳定性条件应比前者严格，下一节将进行验证。

4 稳定性条件比较

为比较3 种人工边界子系统的数值计算稳定条件，采用两组参数进行分析，具体参数见表1。

表1 两组参数值(无量纲)Table 1 Two sets of parameter values (dimensionless)

侧边固定和角点固定子系统的稳定条件可由式(19)和式(24)获得，角点自由边界子系统的稳定条件可采用数值方法获得。3 种子系统谱半径的计算结果见图11。当谱半径小于等于1 时，数值计算满足稳定性条件，由图11 可以直观地比较两组参数条件下3 种子系统的稳定条件。

表2 中对3 种局部子系统和内部系统满足稳定性条件的临界时间步长进行了定量比较。内部系统的稳定性条件( ∆t=L/CP)未考虑粘弹性人工边界单元对数值算法的影响，最为宽松；侧边和角点子系统考虑了粘弹性人工边界单元质量、等效刚度和阻尼对稳定性的影响，稳定性条件要比无人工边界单元时更为严格。四周固定的角点处边界节点只共享1/4 内部单元质量，其子系统截止频率最高，因此稳定性条件最为严格。以上3 种子系统数值积分稳定条件的对比表明，采用粘弹性人工边界单元时，系统数值积分的稳定条件由角点区控制。

图11 三种子系统的谱半径对比Fig. 11 Comparison of spectral radius of three subsystems

表2 稳定性条件比较(无量纲)Table 2 Comparison of stability conditions (dimensionless)

5 算例验证

5.1 均匀半空间模型

为验证以上稳定性分析的准确性，按第一组参数建立有限元近场模型，模型尺寸为320×160，内部介质密度为2500，剪切波速为500，泊松比为0.3，网格尺寸为2×2，模型侧边和底边最外层单元是粘弹性人工边界单元，如图12 所示。采用持时为0.2 的脉冲作为动力荷载，施加于模型中点处，时程曲线如图13 所示。分别按照内部介质数值稳定条件(∆t=0.0021)、侧边子系统稳定条件(∆t=0.00163)、固定边界角点子系统稳定条件(∆t=0.00623)，采用固定时间步长的显示时域逐步积分算法对整体模型进行计算。

图12 均匀半空间算例模型图Fig. 12 Homogeneous half apace example model diagram

图13 动力荷载时程曲线Fig. 13 Dynamic load time history curve

按照内部介质数值稳定条件(∆t=0.0021)计算时，由于不满足侧边(底边)子系统稳定条件，因此P 波传播到距离波源最近的底边节点时发生失稳，如图14 所示；按照侧边子系统稳定条件(∆t=0.00163)计算时，由于不满足角点处子系统稳定条件，波动传播到模型角点处时发生失稳，如图15所示，其中U为介质中的位移，以上两种稳定条件均无法完成整体有限元模型的计算。

图14 按内部稳定条件计算时底边失稳状态(0.11 时刻)Fig. 14 The unstable state of the bottom edge when calculated according to internal stability conditions (time 0.11)

按照角点处子系统稳定条件(∆t=0.000623)计算时，可顺利完成整体模型的动力显式计算，粘弹性人工边界单元很好模拟了外行波向无穷远的辐射，结果如图16 所示。此外，通过进一步的计算分析发现，当整体稳定条件取值略大于角点子系统稳定条件时(例如∆t=0.00085)，模型角点处也发生失稳，失稳状态与图15 相同。

图15 按侧边子系统稳定条件计算时角点失稳状态(0.19 时刻)Fig. 15 The unstable state of the bottom edge when calculated according to the stability conditions of side subsystem (time 0.19)

图16 按角点子系统稳定条件计算结果Fig. 16 Results calculated by the stability conditions of the corner subsystem

以上计算分析验证了采用粘弹性人工边界单元时，整体模型显式数值积分算法的稳定性由角点区域控制，式(24)给出的角点处子系统稳定条件是整体模型数值稳定的充分条件。

5.2 成层半空间模型

为满足实际场地计算需要，对成层半空间算例进行验证。成层半空间有限元模型尺寸为320×160，上半部分内部介质密度2000，剪切波速300，泊松比0.27，下半部分内部介质密度为2500，剪切波速为500，泊松比为0.3，整体模型网格尺寸为2×2，模型侧边和底边最外层单元是粘弹性人工边界单元，如图17 所示。脉冲荷载施加于模型中心，如图13 所示。

表3 给出了成层半空间局部子系统和内部系统满足稳定性条件的临界时间步长。模型中的上层介质只有侧边子系统，无角点子系统，下层介质两者均存在。

图17 成层半空间算例模型图Fig. 17 Layered half space example model diagram

表3 成层半空间稳定性条件比较(无量纲)Table 3 Comparison of stability conditions in layered half space (dimensionless)

经比较发现，上层介质中满足内部系统稳定条件的时间步长(∆t=0.0037)大于侧边子系统稳定时间步长(∆t=0.0028)。而二者均比下层介质内部系统稳定的时间步长 (∆t=0.0021)要宽松。因此，当采用上层介质稳定性条件对整体模型计算时，首先发生失稳的是下层介质的内部系统，如图18所示。

图18 按上层介质稳定条件计算时下层介质内部系统失稳状态(0.1 时刻)Fig. 18 The unstable state of the internal system in lower media when calculated according to the stability conditions in upper media (time 0.1)

采用下层介质内部系统稳定条件(∆t=0.0021)计算时，整体系统失稳状态与图14 相同；采用下层介质侧边子系统稳定条件计算时(∆t=0.00163)，失稳状态和图15 相同。

按照下层角点子系统稳定条件(∆t=0.000623)计算时，可顺利完成整体模型的动力显式计算，分层介质中，粘弹性人工边界单元也可很好地模拟外行波向无穷远的辐射，结果如图19 所示。成层半空间算例进一步验证了采用粘弹性人工边界单元时，整体模型显式数值积分算法的稳定性仍然由角点区域控制。

图19 按下层介质角点子系统稳定条件计算结果Fig. 19 Results calculated by the stability conditions of the corner subsystem in lower media

6 结论与展望

本文将整体模型数值稳定性问题合理转移到若干局部子系统中，充分考虑粘弹性人工边界单元和内部单元节点间的运动耦合效应，通过传递矩阵谱半径分析方法推导出各局部子系统显式时域逐步积分(中心差分)数值稳定条件的解析解和数值解。通过计算软件验证理论分析的可靠性。具体结论如下：

(1) 对于大规模数值计算问题，可选取局部的子系统并对其进行显式时域逐步积分算法的稳定性分析，该稳定性条件与整体系统稳定性条件相同或近似。

(2) 采用粘弹性人工边界单元时，受人工边界单元质量、刚度和阻尼影响，整体系统的稳定性条件与内部介质的稳定性条件不同，前者的稳定性条件更为严格，需使用更小的积分步长以满足整体系统的数值稳定。

(3) 采用粘弹性人工边界单元时，整体模型显式数值积分算法的稳定性由角点区域控制，本文给出了角点处子系统数值积分的稳定性条件，该稳定性条件是整体模型数值积分稳定的充分条件。此外，本文给出的稳定性条件是以正方形平面单元为对象推导的，同样适用于矩形平面单元，由于系统稳定条件受最小单元尺寸影响，可使用矩形单元的最小边长作为参数计算稳定性条件。

下一步展望：

(1) 对于三维粘弹性人工边界单元，同样可以利用本文提出的传递矩阵谱半径分析方法对显式时域逐步积分算法的稳定性条件进行分析，三维问题涉及的局部子结构更为特殊，传递矩阵的生成和特征值求解更加复杂，需进一步开展研究。

(2) 相比隐式算法，显式算法的解耦特性对于求解大范围复杂工程场地问题更有优势。本文的研究成果为在显式算法中合理使用粘弹性人工边界单元提供了理论依据。可在此基础上进一步开展分析和研究，以改善使用粘弹性人工边界单元时显式算法的稳定性，提高大范围工程场地问题的计算效率。