采用黏彈性人工邊界時顯式算法穩(wěn)定性條件*

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

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

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

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

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

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

黏彈性人工邊界由在截斷邊界處設(shè)置的空間解耦的彈簧-阻尼器構(gòu)成，如圖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é)點構(gòu)成的復(fù)雜節(jié)點系統(tǒng)，通過數(shù)值方法判斷其穩(wěn)定性。因未能給出解析形式的人工邊界穩(wěn)定性判別準(zhǔn)則，此類方法的實用性仍有所欠缺。李述濤等總結(jié)和歸納了現(xiàn)有離散模型數(shù)值積分穩(wěn)定性研究工作，得出以下結(jié)論：數(shù)值積分算法的穩(wěn)定性由計算模型的最高階頻率即系統(tǒng)的截止頻率控制；截止頻率對應(yīng)的振型一般呈現(xiàn)局部節(jié)點系相鄰節(jié)點交錯振動的模態(tài)，二維平面應(yīng)變情況下，其振動形式如圖2 所示。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖3 二維半無限空間近場有限元模型及二維平面應(yīng)變單元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)）及平面應(yīng)變單元的剛度和質(zhì)量矩陣（式(7)～(8)），進(jìn)行矩陣組裝，得到側(cè)邊子系統(tǒng)中節(jié)點1 的運動方程：

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

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

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

側(cè)邊和角點子系統(tǒng)的穩(wěn)定性系數(shù)γ和γ均為無量綱系數(shù)，僅與內(nèi)部介質(zhì)的泊松比、波源距與有限單元長度比/有關(guān)。不同區(qū)域的穩(wěn)定性系數(shù)γ 隨和/的變化情況，如圖4 所示；當(dāng)/=1,5,+∞時穩(wěn)定性系數(shù)γ 隨泊松比的變化，當(dāng)泊松比=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)定性均隨泊松比的增大變得寬松，且隨著波源距與單元尺寸比/的增大先略為放松，當(dāng)/＞5 時基本保持不變。此外，在不同情況下，角點子系統(tǒng)的穩(wěn)定性系數(shù)均最小，對整體計算模型的穩(wěn)定性起控制作用，其穩(wěn)定性條件即為采用黏彈性人工邊界時顯式時域逐步積分算法的統(tǒng)一穩(wěn)定性條件?？杀硎緸椋?/p>

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

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

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

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

4.1 均勻半空間模型

為驗證以上穩(wěn)定性分析的準(zhǔn)確性，采用大型通用有限元軟件ABAQUS 建立均勻半空間模型，如圖7所示。模型尺寸為100 m×50 m，內(nèi)部介質(zhì)的密度為2 000 kg/m，剪切波速為200 m/s，泊松比為0.3。采用四節(jié)點平面應(yīng)變單元進(jìn)行有限元離散，網(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

分別采用不同的固定時間步長進(jìn)行顯式動力計算，并統(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)定性條件解析解可準(zhǔn)確判斷數(shù)值計算中的穩(wěn)定性狀態(tài)；當(dāng)所積分時間步長不滿足內(nèi)部計算域、側(cè)邊子系統(tǒng)或角點子系統(tǒng)的穩(wěn)定性條件時，相應(yīng)區(qū)域分別發(fā)生失穩(wěn)，失穩(wěn)狀態(tài)如圖9 所示；當(dāng)積分時間步長滿足角點子系統(tǒng)的穩(wěn)定性條件時，可順利完成整體模型的數(shù)值計算，即采用黏彈性人工邊界時顯式時域逐步積分算法的穩(wěn)定性條件由角點區(qū)控制；由于穩(wěn)定性系數(shù)的建議值比理論值更保守，采用本文中的穩(wěn)定性系數(shù)建議值確定的積分時間步長可保證顯式動力計算的順利完成。

4.2 成層半空間模型

為進(jìn)一步驗證以上結(jié)論在復(fù)雜場地中的適用性，建立成層半空間模型。如圖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é)點平面應(yīng)變單元對計算模型進(jìn)行有限元離散，網(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)定性應(yīng)由下層介質(zhì)控制。分別采用不同的固定時間步長進(jìn)行顯式動力計算，并統(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)定性條件解析解可準(zhǔn)確判斷采用不同積分時間步長進(jìn)行顯式動力計算時整體模型的穩(wěn)定性狀態(tài)；整體模型中同一層介質(zhì)的穩(wěn)定性由角點區(qū)域控制，在實際工程應(yīng)用中，令積分時間步長滿足式(28)的穩(wěn)定性條件，或?qū)⒎e分時間步長取為由式(29)給出的穩(wěn)定性系數(shù)建議值與內(nèi)部計算域臨界時間步長的乘積，即可順利完成整體模型的顯式動力計算。

5 結(jié) 論

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

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

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

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