引言

在无限介质或半无限介质的有限元波动数值模拟中，需人为引入一种人工边界，以获得有限计算模型。这种人为引入的边界会在边界上产生外行波的假反射。因此，需在人工边界上建立人工边界条件，以消除人工边界产生的反射。自20世纪60年代末以来，国内外学者针对全局与局部人工边界已进行了大量研究，并获得了一系列研究成果（廖振鹏，1997）。在已建立的人工边界条件中，基于对单侧波一般运动学特征一维描述的、直接模拟的多次透射公式（MTF），物理概念简单，便于在计算机上实现时空解藕的高精度波动有限元或有限差分数值模拟（Liao等，1984；廖振鹏等，1984；Liao，1996；廖振鹏，1996）。

将集中质量动力时域有限元方法与MTF结合，可实现时空解藕的波动数值模拟，且精度可控。与其它局部人工边界条件一样，稳定实现MTF是波动有限元或有限差分数值模拟的研究重点。MTF结合集中质量动力时域有限元方法模拟失稳现象，包括高频振荡失稳和低频飘移失稳，周正华等（2001）给出了消除多次透射公式飘移失稳的措施。本文将对MTF在波动有限元数值模拟中引起的高频振荡失稳问题进行讨论，并提出相应的消除高频振荡失稳的措施，以稳定实现MTF；在此基础上，通过三维波动数值模拟，检验了提出的消除高频振荡失稳措施的有效性。

1.   多次透射公式高频失稳机理及稳定措施

有关MTF的基本理论及基本公式的导出已有详细介绍（Liao，1996），文中不再赘述。Liao等（1992）和廖振鹏等（1992）基于有限元离散网格对波动数值模拟影响的分析，导出MTF反射系数的一般表达公式，并由此解释了MTF产生高频振荡失稳的主要特征及其机理。此外，谢志南等（2008）通过分析一维波动模型的离散模型，进一步讨论了MTF高频振荡失稳机理。研究结果表明：在对波动有限元数值模拟或有限差分数值模拟有意义的频段内，MTF一般不会产生高频振荡失稳，而对数值模拟无意义的高频段则会产生高频振荡失稳；MTF高频振荡失稳的实质是数值模拟中无意义的高频波动在人为引入的人工边界上发生反射放大，并在数值模拟有限区域内，由人工边界引起多次反射导致在人工边界上发生不断的反射放大，数值模拟值越来越大，最终产生高频振荡失稳现象。在动力有限元或有限差分数值模拟中，MTF随着计算时间的增长有可能会引起数值模拟结果的高频振荡失稳，且振荡的频率远远超过了波动数值模拟有意义的频段，这种现象一般由动力有限元或有限差分数值模拟引起。引起MTF高频振荡失稳的高频波动分量，可能源于入射波中包含此高频成分，或因时域逐步数值积分舍入误差而随机产生的高频成份。若在数值模拟的波动中不含数值模拟无意义的高频波动分量，或这样的波动高频分量甚小，加之与计算随机舍入误差产生的波动高频分量在形成之初其值甚小，则在波动有限元或有限差分数值模拟中，所有这些极其微小的高频波动分量起初并不显现，只有在计算区内经过多次反射放大后，才能在有限元或有限差分数值模拟的结果中显示出来（图 1）。因此，区别于动力有限元稳定条件不满足引起的失稳，由MTF引起的高频振荡失稳并非在数值模拟之初就开始，而是经历一段时间逐步放大直至失稳。基于对MTF高频振荡失稳的上述认识，认为可通过引入高频滤波措施消除MTF产生的高频振荡失稳。廖振鹏（1996）提出了1种消除MTF高频振荡失稳的平滑措施，即通过平滑方法来消除波动的高频分量。由于平滑滤波措施是在MTF波动高频成份形成之后再作滤波，其滤波措施效果不总有效，并不能完全消除MTF在数值模拟中引起的高频振荡失稳。鉴于平滑滤波措施不够理想，关慧敏等（1997）提出了加边界阻尼层的高频滤波措施，但这一措施对于抑制MTF波动高频成份仍不理想。

图 1  MTF数值模拟高频振荡失稳

Figure 1.  High-frequency oscillation instability by MTF numerical simulation

下载: 全尺寸图片 幻灯片

基于平滑滤波与加边界阻尼层措施的启发，尝试通过引入措施对MTF波动高频成份在形成之初就进行压制，避免因人工边界对数值模拟中无意义的高频波动成份的放大而影响数值模拟结果，并通过在整个计算区内引入阻尼效应以消除高频振荡失稳。大量的数值模拟结果显示，在动力有限元数值模拟中，若考虑分析模型介质的阻尼性质，则能有效地抑制MTF引起的高频振荡失稳（杨宇等，2014），同时亦可通过具有高频能耗特性的积分格式消除MTF引起的高频振荡失稳（李小军等，2007；唐晖等，2010；杨宇等，2014）。因此，在分析模型中引入阻尼效应或积分格式的能耗特性，可将无意义的波动高频成份在形成之初消除。

已有的阻尼模型研究结果表明，粘性阻尼与应变速率成正比，且随着频率变高而变大，尤其对高频波动具有较好的抑制效能，但对波动的低频成份影响却很小。因此，本研究将在整个模型计算区内施加粘性阻尼，滤掉数值模拟中不需要考虑的波动高频成份，来实现MTF在动力有限元或有限差分数值模拟中的应用。若在数值模拟分析模型中已具有此类阻尼，且介质阻尼够大，MTF高频振荡失稳现象会自动得到抑制，无需再施加粘性阻尼；若分析模型中无此类阻尼，则可适当施加，其实施办法是在内节点运动方程中，加上用单元阻尼阵[C]所表示的阻尼力，[C]与有限单元刚度矩阵[k]成正比，即：

式中，ω^*为参考圆频率（rad/s），取值通常为大于波动模拟有意义的截止频率；β为无量纲阻尼系数，取值为小正数，决定了引入阻尼的大小；[C]为单元阻尼就矩阵；[k]为单元刚度矩阵。

针对引入阻尼效应消除MTF高频振荡失稳的措施，可通过模态分析对集中质量时域动力有限元方法内节点计算精度的影响进行简单说明。由模态分析可知，动力有限元内节点的运动可通过叠加各振型运动而求得。对于给定的β和ω^*值，第j振型的阻尼比d_j可表示为：

式中，ω_j为第j振型的固有圆频率，小于ω^*。当β取值很小，且ω^*取值比ω_j高时，阻尼比d_j很小，则对动力有限元数值模拟精度的影响很小，可忽略不计。因此，引入与应变速率成正比的粘性阻尼这一措施，只对波动数值模拟中无意义的高频分量具有较强的抑制作用，而不影响动力有限元数值模拟中有意义的较低频段内的计算精度。

2.   数值试验

本文采用时域波动有限元数值模拟的时空解耦方法（廖振鹏等，1984；Liao，1998），即MTF与集中质量有限元相结合的显式时域逐步积分方法，通过典型算例的稳定实现，验证所提措施对MTF引起的高频振荡失稳的抑制效果，同时证明该措施对数值模拟精度影响很小，工程意义上可忽略不计。

分析研究基于三维波源问题的典型算例，考虑在均匀、各向同性线弹性半无限介质自由表面，作用一竖向近似脉冲力产生的波动。采用直角坐标系oxyz，其中oxy平面位于弹性半空间自由表面，z轴垂直向下，且坐标系的原点与竖向近似脉冲力的作用点相同。设定竖向近似脉冲力的作用方向与z轴一致，幅值为10⁴kN，归一化竖向近似脉冲力的时间函数为p（t），脉冲宽度为0.4s，如图 2所示。通过傅里叶分析得到竖向近似脉冲力的截止频率ω_c约为10Hz。

图 2  归一化竖向近似脉冲力时程

Figure 2.  Time-history of the normalized vertical approximate pulse force

下载: 全尺寸图片 幻灯片

从均匀、各向同性的线弹性半无限介质中获得一方盒形有限计算区域，其上边界为半空间自由表面，底边界与侧边界为人工边界，几何尺寸为20m×20m×20m。有限计算区域内介质的力学参数分别为：剪切波速200m/s、密度1.85×10³kg/m³、泊松比0.25。用空间步距为Δx的均匀立方体有限单元对有限计算区进行离散，根据波动有限元数值模拟的精度要求，Δx定为2m，并基于集中质量有限元方法，形成有限计算区内节点的常微分方程组；采用时域显式数值积分方法完成方程组的数值积分（廖振鹏，1996），按方程组数值积分的稳定性准则确定时间步距Δt为0.002s。人工边界节点的离散方程为2阶MTF。由于MTF计算点的位置$x = - j{c_a}\Delta t$（其中j取1或2；${c_a}$为人工波速，取200m/s）与有限元节点的位置通常不一致，因此，需采用具有2阶精度的三点插值方法内插得到计算点的位移。竖向近似脉冲力作用于有限计算区自由表面中点所引起的该点位移反应如图 3所示，图示结果表明，MTF引起了波动有限元数值模拟结果的高频振荡失稳，且随着计算时间变长，失稳振荡幅值越来越大。在数值模拟中采用本文提出的消除MTF高频振荡失稳的措施，即在数值模拟中引入与应变速率成正比的粘性阻尼，且参考圆频率ω^*取为本文分析截止频率对应圆频率的3倍，若振型的固有圆频率等于参考圆频率ω^*且其阻尼比为3%，则由式（2）可确定β为0.06；若该振型的固有圆频率小于截止频率对应的圆频率，其阻尼比小于1%，数值模拟结果见图 4。由图 4可以看出，本文所提的措施能有效地抑制MTF引起的高频振荡失稳。

图 3  竖向近似脉冲力作用下作用点处的位移反应时程

Figure 3.  The displacement reaction of the point under the normalized vertical approximate pulse force

下载: 全尺寸图片 幻灯片

图 4  引入措施后竖向近似脉冲力作用下作用点处的位移反应时程

Figure 4.  The displacement reaction of the point under the normalized vertical approximate pulse force which can eliminate high-frequency oscillation instability

下载: 全尺寸图片 幻灯片

为检验消除MTF高频振荡失稳措施对数值模拟计算精度的影响，本文通过扩大有限计算区域，获得数值模拟精确解。扩大有限计算区域旨在减小甚至消除人工边界的引入对数值模拟计算精度的影响，几何尺寸为240m×240m×120m，其离散单元大小、模型介质力学参数、数值积分方法、计算时间步距以及MTF的实施均与上述分析相同，且未引入与应变速率成正比的粘性阻尼，作用力仍然在模型自由表面的中点，数值模拟结果如图 5所示。为方便对比分析，图中还给出了采用消除MTF高频振荡失稳措施的计算结果。若以最大体波波速v_p=350m/s计算，则在0.6s内，作用点处的数值模拟结果未受到人工边界的影响。图 5结果表明，2种分析模型的计算结果基本相同；数值精确解的峰值为0.043286m，而采用消除MTF高频振荡失稳措施的峰值为0.043276m，相对误差约为2.3‰，说明本文提出的措施对计算精度的影响很小，可以忽略不计。

图 5  2种模型数值模拟结果对比

Figure 5.  The comparison of results from two numerical simulation models

下载: 全尺寸图片 幻灯片

综上所述，本文提出的消除MTF高频振荡失稳措施行之有效，且对数值模拟分析精度影响甚微。

3.   结语

本文通过对多次透射公式（MTF）高频振荡失稳机理的讨论，提出了在有限元或有限差分波动数值模拟中针对MTF高频振荡失稳的一种消除措施，即在整个计算区域内施加与应变速率成正比的粘性阻尼；并以此为基础，从理论上论证了这一稳定措施的有效性及其对数值模拟精度的影响；最后，通过数值模拟分析进行了检验。结果表明：

（1）引入与应变速率成正比的粘性阻尼，能有效地消除MTF引起的高频振荡失稳，且这一消除措施只对波动数值模拟无意义的高频分量具有抑制作用。

（2）若β的取值很小，且参考圆频率ω^*大于波动数值模拟截止圆频率ω_c时，这一消除高频振荡失稳措施对有意义的较低频段内的波动数值模拟精度影响很小，可忽略不计。