劉晶波，寶 鑫，李述濤,2，王 菲,3

（1. 清華大學土木工程系，北京 100084；2. 軍事科學院國防工程研究院，北京 100036；3. 陸軍工程大學國防工程學院，江蘇 南京 210007）

無限介質(zhì)的波動輻射效應是波動數(shù)值模擬中需要考慮的關鍵問題，人工邊界技術是處理此類問題最常用的技術手段。該方法在無限或半無限介質(zhì)中截取有限的近場計算域，并在截斷邊界處施加人工邊界條件，以吸收計算域內(nèi)產(chǎn)生的外行波動。由于采用不同的數(shù)學、物理或力學原理，常用的人工邊界技術可分為透射邊界、黏性邊界、黏彈性邊界、邊界元和完美匹配層等。其中，在黏彈性人工邊界技術中，將介質(zhì)中單側(cè)波動的偏微分方程轉(zhuǎn)化為施加于截斷邊界上的應力邊界條件，并且等效為空間解耦的力學系統(tǒng)，物理意義清晰、實用性強，且具有較好的模擬精度和良好的魯棒性。近年來，已被研究人員集成于Marc、LS-DYNA、ANSYS、ADINA、Nastran、ABAQUS等通用有限元軟件，并應用于大壩、橋梁、核工程、隧道和地鐵車站等建構筑物與基礎的動力相互作用分析及工程場地的動力響應計算，取得了合理準確的計算結(jié)果。

采用顯式時域逐步積分算法對含有黏彈性人工邊界的整體模型進行計算分析時，受黏彈性人工邊界阻尼、剛度等的影響，人工邊界區(qū)的穩(wěn)定性需比內(nèi)部計算域的穩(wěn)定性條件更嚴格。目前，尚無明確、實用的考慮黏彈性人工邊界影響時顯式時域逐步積分算法穩(wěn)定性的判別準則，難以確定合理的數(shù)值積分時間步長，這一定程度限制了黏彈性人工邊界在顯式動力分析中的應用。與隱式算法相比，顯式動力計算方法不需要求解耦聯(lián)方程組，計算工作量較小、計算效率較高，在大規(guī)模復雜系統(tǒng)的動力分析尤其爆炸問題數(shù)值模擬中應用廣泛。鑒于此，有必要對使用黏彈性人工邊界時顯式時域逐步積分算法的穩(wěn)定性開展研究工作，以促進黏彈性人工邊界在大規(guī)模顯式動力分析中的應用。

本文中，針對含黏彈性人工邊界的數(shù)值模型在顯式算法中的穩(wěn)定性問題，利用基于局部子系統(tǒng)的數(shù)值穩(wěn)定性分析方法，給出可代表整體模型局部特征的不同邊界子系統(tǒng)的穩(wěn)定性條件解析解，比較分析不同計算區(qū)域的穩(wěn)定性條件及其影響因素，給出整體模型在顯式動力計算中的統(tǒng)一穩(wěn)定性判別準則和簡化實用計算方法。

1 黏彈性人工邊界模型與參數(shù)設置

黏彈性人工邊界由在截斷邊界處設置的空間解耦的彈簧-阻尼器構成，如圖1 所示。

圖1 黏彈性人工邊界Fig. 1 Schematic diagram of viscoelastic artificial boundaries

二維情況下，黏彈性人工邊界等效物理系統(tǒng)的彈簧剛度和阻尼系數(shù)分別為：

表1 二維黏彈性人工邊界參數(shù)的數(shù)據(jù)[23]Table 1 The values of two-dimensional viscoelastic artificial boundary coefficients[23]

2 采用黏彈性人工邊界時顯式算法的穩(wěn)定性分析方法

2.1 基于局部子系統(tǒng)的數(shù)值積分穩(wěn)定性分析方法

數(shù)值穩(wěn)定性是制約顯式動力計算時間步長的主要因素，以往針對含人工邊界條件的計算系統(tǒng)的穩(wěn)定性研究，多基于整體模型或由若干排（列）節(jié)點構成的復雜節(jié)點系統(tǒng)，通過數(shù)值方法判斷其穩(wěn)定性。因未能給出解析形式的人工邊界穩(wěn)定性判別準則，此類方法的實用性仍有所欠缺。李述濤等總結(jié)和歸納了現(xiàn)有離散模型數(shù)值積分穩(wěn)定性研究工作，得出以下結(jié)論：數(shù)值積分算法的穩(wěn)定性由計算模型的最高階頻率即系統(tǒng)的截止頻率控制；截止頻率對應的振型一般呈現(xiàn)局部節(jié)點系相鄰節(jié)點交錯振動的模態(tài)，二維平面應變情況下，其振動形式如圖2 所示。

圖2 基于局部子系統(tǒng)的穩(wěn)定性分析Fig. 2 Stability analysis based on local subsystems

在上述基礎上，李述濤等提出了一種利用局部子系統(tǒng)估算整體有限元模型數(shù)值穩(wěn)定性的方法。根據(jù)該方法，對于任意整體有限元模型，如能較準確地判斷其最高階振型，則可利用振型節(jié)點的空間分布規(guī)律，截取由相鄰振型節(jié)點包圍的最小局部子系統(tǒng)（見圖2 中的子系統(tǒng)a）。對該子系統(tǒng)的邊界節(jié)點施加與整體模型最高階振型一致的約束條件，并進行穩(wěn)定性分析，獲得的局部子系統(tǒng)穩(wěn)定性條件即為整體模型的穩(wěn)定性條件。而對于不能準確判斷整體模型最高階振型的情況，可選取由相鄰單元構成的最小子系統(tǒng)（見圖2 中的子系統(tǒng)b）。對該子系統(tǒng)的全部邊界節(jié)點施加固定約束，分析所有b 型局部子系統(tǒng)的穩(wěn)定性條件，其下限值為整體有限元模型穩(wěn)定性條件的上限逼近。

基于局部子系統(tǒng)的穩(wěn)定性分析方法，提供了一種從理論上估算含人工邊界的復雜離散系統(tǒng)數(shù)值穩(wěn)定性的技術手段。對于整體模型的內(nèi)部計算域，可以采用子系統(tǒng)a 直接獲穩(wěn)定性條件，而對于邊界區(qū)域，由于人工邊界的影響，難以準確判斷最高階振型，無法從整體計算模型中分割出子系統(tǒng)a，此時可以采用b 型子系統(tǒng)進行穩(wěn)定性分析。

2.2 基于傳遞矩陣譜半徑的穩(wěn)定性判別方法

穩(wěn)定性條件因所采用的數(shù)值積分方法而異，為解決實際工程應用中黏彈性人工邊界的穩(wěn)定性問題，以在通用有限元軟件中應用較廣泛的蛙跳格式中心差分時域逐步積分算法為基礎，分析含黏彈性人工邊界的整體計算模型的數(shù)值穩(wěn)定性。對于其他顯式時域逐步積分算法，采用相似的步驟即可得到對應的穩(wěn)定性條件。

數(shù)值積分方法的穩(wěn)定性僅與積分格式、時空離散步長和計算系統(tǒng)的力學參數(shù)有關，而與外力向量無關。如滿足以下條件，則積分格式是穩(wěn)定的：(1) ρ() ≤ 1，其中ρ()為傳遞矩陣的譜半徑，即ρ() =max|λ|，λ為傳遞矩陣的第個特征值，對于本文分析的二維情況，=1～4；(2) 如具有多重特征值，則該特征值的模小于1。

結(jié)合以上理論與方法，可將黏彈性人工邊界在顯式時域逐步積分算法中的穩(wěn)定性分析分解為3 個步驟：（1）選取能體現(xiàn)整體有限元模型不同局部特征的子系統(tǒng)；（2）建立各子系統(tǒng)的運動方程，結(jié)合時域逐步積分算法推導傳遞矩陣，并求解譜半徑；（3）根據(jù)譜半徑的模小于等于1 的限制條件求解各子系統(tǒng)的穩(wěn)定性條件，比較后選取最嚴格的穩(wěn)定性條件作為整體模型最大穩(wěn)定積分時間步長的上限估計。

另外，在滿足了離散模型動力計算的穩(wěn)定性條件下，為進一步滿足精度要求，離散化網(wǎng)格的尺寸Δ需滿足以下條件：

式中：λ為離散網(wǎng)格中波動傳播的最短波長，為介質(zhì)中的最小波速，為波動問題數(shù)值模擬的截止頻率。

3 采用黏彈性人工邊界時顯式算法的穩(wěn)定性分析

下面，針對二維情況，分析采用黏彈性人工邊界時顯式時域逐步積分算法的穩(wěn)定性。由于穩(wěn)定性條件與計算系統(tǒng)的單元尺寸和物理特性密切相關，為不失一般性，在數(shù)值算法的穩(wěn)定性分析中通常采用規(guī)則幾何模型、均一介質(zhì)和均勻離散網(wǎng)格進行建模分析。對于不規(guī)則網(wǎng)格或非均勻介質(zhì)，可根據(jù)有限元網(wǎng)格的最小尺寸及不同區(qū)域的介質(zhì)材料參數(shù)，采用本文方法分別計算穩(wěn)定性條件，并選取其下限值作為最大積分時間步長的選取依據(jù)。

二維半無限空間近場有限元模型如圖3(a)所示，不考慮介質(zhì)阻尼，內(nèi)部均勻介質(zhì)在中心差分算法中的穩(wěn)定性條件為：

當采用黏彈性人工邊界時，因難以確定人工邊界區(qū)的最高階振型，應采用b 型子系統(tǒng)進行穩(wěn)定性分析。此時，根據(jù)截取位置的不同，人工邊界區(qū)的b 型局部子系統(tǒng)可分為側(cè)邊（底邊）子系統(tǒng)和角點子系統(tǒng)，如圖3(a)所示。

為建立局部子系統(tǒng)的運動方程，先給出規(guī)則二維平面應變等參數(shù)單元（見圖3(b)）的剛度矩陣：

圖3 二維半無限空間近場有限元模型及二維平面應變單元Fig. 3 Two-dimensional semi-infinite near-field finite element model and two-dimensional plane strain element

式中：為8×8 型的單位矩陣。

3.1 側(cè)邊子系統(tǒng)的穩(wěn)定性條件

在圖3(a)所示的側(cè)邊子系統(tǒng)中，僅節(jié)點1 為自由節(jié)點，其余節(jié)點均施加了固定約束。利用黏彈性人工邊界物理參數(shù)（式(1)）及平面應變單元的剛度和質(zhì)量矩陣（式(7)～(8)），進行矩陣組裝，得到側(cè)邊子系統(tǒng)中節(jié)點1 的運動方程：

3.2 角點子系統(tǒng)的穩(wěn)定性條件

3.3 穩(wěn)定性條件的比較

由于整體模型的穩(wěn)定性受不同局部區(qū)域中最嚴格的穩(wěn)定性條件控制，根據(jù)式(6)、(21) 和(26)，可將采用黏彈性人工邊界時顯式算法的穩(wěn)定性條件統(tǒng)一改寫為：

側(cè)邊和角點子系統(tǒng)的穩(wěn)定性系數(shù)γ和γ均為無量綱系數(shù)，僅與內(nèi)部介質(zhì)的泊松比、波源距與有限單元長度比/有關。不同區(qū)域的穩(wěn)定性系數(shù)γ 隨和/的變化情況，如圖4 所示；當/=1,5,+∞時穩(wěn)定性系數(shù)γ 隨泊松比的變化，當泊松比=0.2,0.3,0.4 時穩(wěn)定性系數(shù)γ 隨/的變化，如圖5～6 所示。

圖4 穩(wěn)定性條件的比較Fig. 4 Comparison of stability conditions

圖5 不同R/L 時穩(wěn)定性系數(shù)γ 隨泊松比μ的變化Fig. 5 Variations of stability coefficient γ with Poisson’s ratio μ under different R/L

圖6 不同泊松比μ時穩(wěn)定性系數(shù)γ 隨R/L 的變化Fig. 6 Variations of stability coefficient γ with R/L under different Poisson’s ratio μ

由圖4～6 可見，側(cè)邊和角點子系統(tǒng)的數(shù)值積分穩(wěn)定性均隨泊松比的增大變得寬松，且隨著波源距與單元尺寸比/的增大先略為放松，當/＞5 時基本保持不變。此外，在不同情況下，角點子系統(tǒng)的穩(wěn)定性系數(shù)均最小，對整體計算模型的穩(wěn)定性起控制作用，其穩(wěn)定性條件即為采用黏彈性人工邊界時顯式時域逐步積分算法的統(tǒng)一穩(wěn)定性條件?？杀硎緸椋?/p>

在采用黏彈性人工邊界和中心差分法進行動力分析時，根據(jù)黏彈性人工邊界參數(shù)α、α和波源距與單元尺寸之比/及泊松比，即可由式(28)確定穩(wěn)定的時間積分步長。在實際應用中，為了更方便地確定穩(wěn)定時間積分步長 Δ，可采納幾種常見情況的穩(wěn)定性系數(shù)，見表2。

表2 建議的幾種常見情況的穩(wěn)定性系數(shù)Table 2 Recommended stability coefficients for several common cases

在實際工程中，通常滿足材料泊松比 μ≥0.1 、波源距與單元尺寸比/≥5 ，且穩(wěn)定性條件隨泊松比和波源距與單元尺寸比/的增大變得寬松。在實際應用中，為避免復雜的公式計算，可直接采用μ=0.1、/=5對應的穩(wěn)定性條件作為系統(tǒng)最大穩(wěn)定積分時間步長的保守估計：

4 穩(wěn)定性條件驗證

4.1 均勻半空間模型

為驗證以上穩(wěn)定性分析的準確性，采用大型通用有限元軟件ABAQUS 建立均勻半空間模型，如圖7所示。模型尺寸為100 m×50 m，內(nèi)部介質(zhì)的密度為2 000 kg/m，剪切波速為200 m/s，泊松比為0.3。采用四節(jié)點平面應變單元進行有限元離散，網(wǎng)格尺寸為1 m×1 m。在模型的截斷邊界處施加黏彈性人工邊界，并令人工邊界參數(shù)α=0.5、α=1.0。在模型中點（點）處施加持時0.2 s、幅值為1 MN 的脈沖荷載，如圖8 所示。

圖7 均勻半空間模型Fig. 7 The homogeneous half-space model

圖8 脈沖荷載Fig. 8 The impulse load

根據(jù)以上模型參數(shù)，分別采用式(6)、(21)和(26)，給出的穩(wěn)定性條件及式(29)給出的穩(wěn)定性系數(shù)建議值，計算的模型不同區(qū)域的最大穩(wěn)定時間步長見表3。

表3 均勻半空間模型的穩(wěn)定性系數(shù)和最大穩(wěn)定時間步長Table 3 Stability coefficients and maximum stable time steps of the homogeneous model

分別采用不同的固定時間步長進行顯式動力計算，并統(tǒng)計其穩(wěn)定性狀態(tài)，結(jié)果見表4。提取失穩(wěn)時計算模型的位移云圖，如圖9 所示，并比較采用不同固定積分時間步長時模型底部角點（點）的豎向位移，如圖10 所示。

表4 不同固定時間步長時均勻半空間模型的穩(wěn)定性狀態(tài)Table 4 The stability states of the homogeneous model under different fixed time steps

圖9 均勻半空間模型的位移分布Fig. 9 Displacement distributions of the homogeneous half-space model

圖10 均勻半空間模型底部角點的豎向位移Fig. 10 Vertical displacements of the corner point in the homogenous half-space model

由表3～4 和圖9～10，可以得到以下結(jié)論：本文中的采用黏彈性人工邊界時顯式時域逐步積分算法的穩(wěn)定性條件解析解可準確判斷數(shù)值計算中的穩(wěn)定性狀態(tài)；當所積分時間步長不滿足內(nèi)部計算域、側(cè)邊子系統(tǒng)或角點子系統(tǒng)的穩(wěn)定性條件時，相應區(qū)域分別發(fā)生失穩(wěn)，失穩(wěn)狀態(tài)如圖9 所示；當積分時間步長滿足角點子系統(tǒng)的穩(wěn)定性條件時，可順利完成整體模型的數(shù)值計算，即采用黏彈性人工邊界時顯式時域逐步積分算法的穩(wěn)定性條件由角點區(qū)控制；由于穩(wěn)定性系數(shù)的建議值比理論值更保守，采用本文中的穩(wěn)定性系數(shù)建議值確定的積分時間步長可保證顯式動力計算的順利完成。

4.2 成層半空間模型

為進一步驗證以上結(jié)論在復雜場地中的適用性，建立成層半空間模型。如圖11 所示，模型尺寸為100 m×100 m，分為上下兩層，層厚均為50 m，上層介質(zhì)的密度、剪切波速和泊松比分別為1 800 kg/m、150 m/s 和0.28，下層介質(zhì)的密度、剪切波速和泊松比分別為2 200 kg/m、250 m/s 和0.35。采用四節(jié)點平面應變單元對計算模型進行有限元離散，網(wǎng)格尺寸為1 m×1 m。在模型截斷邊界處施加黏彈性人工邊界，在模型中點（點）處施加脈沖荷載（見圖8）。

圖11 成層半空間模型Fig. 11 The layered half-space model

由以上參數(shù)計算得到的成層半空間模型不同區(qū)域的穩(wěn)定性系數(shù)和最大穩(wěn)定時間步長見表5，上層介質(zhì)及其側(cè)邊子系統(tǒng)的穩(wěn)定性條件比下層介質(zhì)的穩(wěn)定性條件更寬松，可以判斷該模型的數(shù)值積分穩(wěn)定性應由下層介質(zhì)控制。分別采用不同的固定時間步長進行顯式動力計算，并統(tǒng)計其穩(wěn)定性狀態(tài)，結(jié)果見表6。提取失穩(wěn)時計算模型的位移云圖，如圖12 所示，并比較采用不同固定時間步長時模型底部角點（點）的豎向位移，如圖13 所示。

圖12 成層半空間模型的位移分布Fig. 12 Displacement distributions of the layered half-space model

圖13 成層半空間模型底部角點的豎向位移Fig. 13 Vertical displacements of the corner point in the layered half-space model

表5 成層半空間模型的穩(wěn)定性系數(shù)與最大穩(wěn)定時間步長Table 5 Stability coefficients and maximum stable time steps of the layered model

表6 采用不同固定時間步長時成層半空間模型的穩(wěn)定性狀態(tài)Table 6 The stability state of the layered model under different fixed time steps

由表5～6 和圖12～13，可以得到與第4.1 節(jié)相似的結(jié)論：本文中穩(wěn)定性條件解析解可準確判斷采用不同積分時間步長進行顯式動力計算時整體模型的穩(wěn)定性狀態(tài)；整體模型中同一層介質(zhì)的穩(wěn)定性由角點區(qū)域控制，在實際工程應用中，令積分時間步長滿足式(28)的穩(wěn)定性條件，或?qū)⒎e分時間步長取為由式(29)給出的穩(wěn)定性系數(shù)建議值與內(nèi)部計算域臨界時間步長的乘積，即可順利完成整體模型的顯式動力計算。

5 結(jié) 論

以二維黏彈性人工邊界為研究對象，利用基于局部子系統(tǒng)的穩(wěn)定性分析方法研究其在顯式時域逐步積分算法中的穩(wěn)定性，并通過數(shù)值算例對穩(wěn)定性分析的準確性加以驗證。具體結(jié)論如下。

（1）給出了采用黏彈性人工邊界時，模型側(cè)邊和角點區(qū)域在顯式時域逐步積分算法中的穩(wěn)定性條件解析解，獲得了含黏彈性人工邊界的整體模型在顯示動力計算中的統(tǒng)一穩(wěn)定性判別準則，在數(shù)值模擬時可利用該解析解對整體模型的穩(wěn)定性進行預判。

（2）考慮黏彈性人工邊界影響時，整體模型的數(shù)值積分穩(wěn)定性條件比內(nèi)部計算域更嚴格，整體模型的穩(wěn)定性由角點區(qū)域控制。側(cè)邊和角點區(qū)域的穩(wěn)定性條件僅與波源距與單元尺寸比/、泊松比和內(nèi)部計算域的穩(wěn)定時間步長/有關。隨著/和的增大，人工邊界區(qū)的穩(wěn)定性變得寬松，當/＞5 時，穩(wěn)定性條件基本保持不變。

（3）分析穩(wěn)定性條件的控制因素及其影響規(guī)律的，給出了采用黏彈性人工邊界時顯式時域逐步積分算法穩(wěn)定性系數(shù)的保守建議值。在實際應用中，可直接采用該穩(wěn)定性系數(shù)建議值與內(nèi)部計算域臨界時間步長的乘積，作為顯式動力計算最大穩(wěn)定時間步長的估計。