李述濤 劉晶波 寶 鑫,2)

?(清華大學(xué)土木工程系，北京 100084)

?(軍事科學(xué)院國防工程研究院，北京 100036)

引言

無限域中波動問題的求解是眾多物理領(lǐng)域所面對的重要研究課題.針對此問題最初的分析方法一般是將截?cái)噙吔缛〉阶銐蜻h(yuǎn)的位置，在截?cái)噙吔绠a(chǎn)生的反射波返回感興趣區(qū)域前完成計(jì)算過程[1].這一方法具有較高的精度，但計(jì)算成本過高.人工邊界技術(shù)是目前常用的無限域介質(zhì)數(shù)值模擬方法，將人工邊界引入無限域空間的近場模型中，可吸收截?cái)噙吔缣幍耐庑胁▌覽2-3].在眾多的人工邊界技術(shù)中，黏彈性人工邊界基于全空間均勻介質(zhì)中的柱面波和球面波理論推導(dǎo)得到，可等效為在每個邊界節(jié)點(diǎn)上施加的遠(yuǎn)端固定的并聯(lián)彈簧--阻尼器系統(tǒng)，既具有較高的計(jì)算精度，又可在多種有限元軟件中進(jìn)行集成，獲得了較為廣泛的工程應(yīng)用[4-5].但由于離散彈簧--阻尼器的施加需依賴程序語言對模型進(jìn)行前處理，實(shí)用性受到一定影響.為解決這一問題，劉晶波等[6]和谷音等[7]發(fā)展了黏彈性人工邊界單元，其思路是利用有限單元的矩陣等效構(gòu)造等效人工邊界單元.該方法僅需指定模型最外層單元的密度、彈性模量、泊松比和阻尼系數(shù)，即可模擬集中黏彈性人工邊界，且計(jì)算精度不變.

隨著黏彈性人工邊界單元(以下簡稱邊界單元)在實(shí)際工程中的廣泛應(yīng)用[8-12]，其數(shù)值穩(wěn)定性問題隨之而來.研究人員在引入邊界單元進(jìn)行波動分析時，一般采用隱式無條件穩(wěn)定的時域逐步積分算法，較少采用顯式算法[13-14].這是因?yàn)檫吔鐔卧哂信c內(nèi)部介質(zhì)不同的質(zhì)量密度、剛度和阻尼，導(dǎo)致在邊界單元區(qū)，顯式算法的穩(wěn)定條件不同于內(nèi)部計(jì)算區(qū).由于邊界單元區(qū)與相鄰的內(nèi)部計(jì)算區(qū)的力學(xué)參數(shù)完全不同，使得顯式算法的穩(wěn)定性分析遇到技術(shù)上的困難，長久以來未能給出與內(nèi)部計(jì)算區(qū)類似的解析形式的數(shù)值積分穩(wěn)定性條件，這在一定程度上影響了采用邊界單元時對顯式算法穩(wěn)定性的判斷，而數(shù)值算法穩(wěn)定條件的缺失也限制了邊界單元在顯式動力分析中的應(yīng)用.

顯式算法的解耦特性回避了組集整體剛度陣和求解聯(lián)立方程組等問題[15-18]，在大規(guī)模動力問題數(shù)值分析中的計(jì)算效率顯著高于隱式算法，特別是對于大范圍場地模型的并行計(jì)算問題，因此有必要對引入黏彈性人工邊界單元時顯式算法的穩(wěn)定性開展研究.

數(shù)值穩(wěn)定性的分析和改善研究一直伴隨著人工邊界理論的發(fā)展[19-22].針對透射邊界，可以改變與邊界物理量相關(guān)的計(jì)算格式，通過濾波提高邊界穩(wěn)定性[23].對于多次透射邊界，可以在邊界區(qū)設(shè)置黏性阻尼，以消除多次透射邊界的高頻振蕩失穩(wěn)現(xiàn)象[24].黏彈性人工邊界單元本身為物理上穩(wěn)定的系統(tǒng)，將其引入計(jì)算系統(tǒng)中，并不影響整體的物理穩(wěn)定條件.其失穩(wěn)機(jī)理是顯式時域逐步積分算法所導(dǎo)致的，而與邊界單元理論無關(guān).為改善數(shù)值穩(wěn)定條件，需要全面考慮逐步積分過程中邊界節(jié)點(diǎn)與內(nèi)部單元節(jié)點(diǎn)之間的耦合效應(yīng)，以此為依據(jù)調(diào)整滿足顯式時域積分?jǐn)?shù)值穩(wěn)定條件的臨界時間步長，從而提高人工邊界數(shù)值計(jì)算的穩(wěn)定性.

本文針對人工邊界典型位置處的局部子系統(tǒng)建立運(yùn)動方程，利用傳遞矩陣譜半徑分析方法，推導(dǎo)得到解析形式的局部子系統(tǒng)穩(wěn)定性條件表達(dá)式.通過研究解析解中各物理參數(shù)對穩(wěn)定性條件的影響，給出通過增加人工邊界單元質(zhì)量密度以實(shí)現(xiàn)改善顯式算法穩(wěn)定性條件的方法.本文研究對于提高黏彈性人工邊界單元在動力顯式計(jì)算中的適用性以及提高大范圍工程場地波動分析的計(jì)算效率均具有重要意義.

1 黏彈性人工邊界單元等效物理參數(shù)設(shè)置

黏彈性人工邊界單元是在無限域半空間模型的截?cái)噙吔缟舷蛲夥ň€方向延伸出的一層厚度為h的單元，通過設(shè)置等效物理參數(shù)來模擬黏彈性人工邊界. 二維黏彈性人工邊界單元的等效密度、等效彈性模量、等效泊松比和等效阻尼系數(shù)分別為[6]

其中，G為內(nèi)部介質(zhì)的剪切模量，ρ 為內(nèi)部介質(zhì)的質(zhì)量密度，CS和CP分別為內(nèi)部介質(zhì)的S 波和P 波波速，R為散射波源至邊界節(jié)點(diǎn)的距離，αT與αN分別為切向與法向黏彈性人工邊界參數(shù)，對于二維黏彈性人工邊界，其推薦值分別為0.5 和1.

黏彈性人工邊界單元厚度h的魯棒性良好，可靈活取值而對精度影響很小[6].由于顯式算法的穩(wěn)定條件受模型中的最小單元尺寸制約，在實(shí)際計(jì)算分析中，建議將邊界單元的厚度設(shè)定為與內(nèi)部單元尺寸一致，如圖1 所示，以排除單元尺寸對整體穩(wěn)定性的影響.

圖1 二維黏彈性人工邊界單元示意圖Fig.1 Diagrams of 2D viscoelastic artificial boundary element

2 邊界局部子系統(tǒng)穩(wěn)定性分析

為改善數(shù)值積分穩(wěn)定性，首先應(yīng)對失穩(wěn)機(jī)理開展研究.若不考慮人工邊界的影響，可使用馮諾依曼方法對計(jì)算區(qū)域進(jìn)行穩(wěn)定性分析[25]；引入邊界單元并進(jìn)行顯式時域逐步積分計(jì)算時，則需要考慮邊界單元對整體計(jì)算系統(tǒng)穩(wěn)定性的影響.此時可以建立包含人工邊界在內(nèi)的整體系統(tǒng)的運(yùn)動方程，采用傳遞矩陣譜半徑分析方法，獲得滿足穩(wěn)定條件的臨界時間步長并進(jìn)行失穩(wěn)判斷.但在實(shí)際應(yīng)用中，整體模型傳遞矩陣的建立和分析工作量巨大，通常難以施行.

針對透射人工邊界條件，關(guān)慧敏等[21]和Kamel等[26]提出了一種通過分析局部子系統(tǒng)傳遞矩陣特征值來研究透射人工邊界穩(wěn)定性的方法.由于該方法處理的子系統(tǒng)包含節(jié)點(diǎn)數(shù)量較多，只能給出系統(tǒng)穩(wěn)定條件的數(shù)值解，難以給出解析解，不適合對穩(wěn)定性條件的影響因素進(jìn)行定量化的參數(shù)分析，制約了針對穩(wěn)定性改善方法的研究.

針對黏彈性人工邊界單元，李述濤等[27]對一維和二維有限元系統(tǒng)最高階模態(tài)進(jìn)行了分析，證明了可以采用極小的子系統(tǒng)推導(dǎo)得出整體系統(tǒng)的穩(wěn)定性條件.該子系統(tǒng)中只包含一個運(yùn)動節(jié)點(diǎn)，所得的穩(wěn)定性條件是整體系統(tǒng)數(shù)值穩(wěn)定的充分條件.此外，該方法給出了系統(tǒng)穩(wěn)定條件的解析解，可以直觀且定量地分析各物理參數(shù)對數(shù)值穩(wěn)定性的影響，為數(shù)值積分穩(wěn)定性的改善創(chuàng)造了條件.

引入黏彈性人工邊界單元后，可以選擇對整體模型邊界處具有典型特征的節(jié)點(diǎn)子系統(tǒng)開展研究.一旦確定了子系統(tǒng)形式，即可根據(jù)顯式逐步積分算法格式，將子系統(tǒng)的運(yùn)動方程寫成如下形式

其中向量U由子系統(tǒng)各節(jié)點(diǎn)的加速度、速度和位移組成；A是傳遞矩陣，上標(biāo)p是自然數(shù)，t=p?t；?t為時間步長.如果滿足下列兩條件，則積分格式(2)是穩(wěn)定的[28].條件1：ρ(A)1，ρ(A)是傳遞矩陣A的譜半徑；條件2：如果A具有多重特征值，則該特征值的模小于1.由此可將子系統(tǒng)的穩(wěn)定性分析歸結(jié)為傳遞矩陣A的譜半徑計(jì)算.

圖2 所示是均勻離散的二維有限元整體模型中的兩處典型節(jié)點(diǎn)子系統(tǒng).其一位于模型的側(cè)邊或底邊，由2 個內(nèi)部介質(zhì)單元和2 個邊界單元組成.另一子系統(tǒng)位于模型角點(diǎn)處，由1 個內(nèi)部介質(zhì)單元和3 個邊界單元組成.

圖2 二維模型中的兩種節(jié)點(diǎn)子系統(tǒng)Fig.2 Two kinds of nodal subsystems in 2D model

文獻(xiàn)[27]采用上述研究思路，將側(cè)邊子系統(tǒng)與角點(diǎn)子系統(tǒng)從整體模型中分離出來，在四周施加固定約束.利用傳遞矩陣的譜半徑分析時域逐步積分算法的穩(wěn)定性條件，最終推導(dǎo)出整體系統(tǒng)的穩(wěn)定性條件.分析結(jié)果表明，采用邊界單元時，整體模型數(shù)值積分的穩(wěn)定性條件要比內(nèi)部系統(tǒng)的穩(wěn)定性條件更為嚴(yán)格，這在一定程度上限制了顯式數(shù)值積分方法的應(yīng)用.

現(xiàn)有數(shù)值積分方法的穩(wěn)定性研究表明，當(dāng)介質(zhì)的質(zhì)量密度增大時(相當(dāng)于介質(zhì)波速降低或結(jié)構(gòu)自振周期變長)，系統(tǒng)的數(shù)值積分穩(wěn)定性條件將得到改善.在原始的黏彈性人工邊界單元中，等效質(zhì)量密度?ρ=0，但由于邊界單元具有較良好的魯棒性，適當(dāng)調(diào)整參數(shù)對人工邊界的模擬精度影響不大.因此可以通過適當(dāng)增加邊界單元的質(zhì)量密度以達(dá)到改善整體模型數(shù)值積分穩(wěn)定性的目的.

下面首先給出考慮邊界單元質(zhì)量密度時，角點(diǎn)子系統(tǒng)和側(cè)邊子系統(tǒng)的穩(wěn)定性條件.

2.1 角點(diǎn)子系統(tǒng)穩(wěn)定性條件

角點(diǎn)子系統(tǒng)如圖3 所示.設(shè)內(nèi)部介質(zhì)單元的彈性模量為E，剪切模量為G，密度為ρ，泊松比為μ，且無阻尼；黏彈性人工邊界單元的密度為，彈性模量為，阻尼系數(shù)為，泊松比為0，單元邊長均為L.

圖3 角點(diǎn)子系統(tǒng)Fig.3 Corner subsystem

二維平面單元的集中質(zhì)量矩陣、剛度矩陣、阻尼矩陣分別參考文獻(xiàn)[6,29-30]，按節(jié)點(diǎn)編號組裝矩陣后，得到角點(diǎn)子系統(tǒng)中5 號節(jié)點(diǎn)的運(yùn)動方程如下

根據(jù)式(4)將式(3)展開，并寫成傳遞矩陣形式

其中

根據(jù)穩(wěn)定性判別條件ρ(A1)1，得出角點(diǎn)子系統(tǒng)穩(wěn)定性條件的解析表達(dá)式為

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

其中，γ1為角點(diǎn)子系統(tǒng)穩(wěn)定性參數(shù)，整理后如式(11)所示.

由式(10)和式(11)可以發(fā)現(xiàn)，當(dāng)內(nèi)部介質(zhì)的P 波波速CP和單元尺寸L給定后，角點(diǎn)子系統(tǒng)數(shù)值穩(wěn)定條件僅與質(zhì)量比/ρ、內(nèi)部介質(zhì)的泊松比μ和單元尺寸與散射波源至邊界距離之比L/R有關(guān).由于內(nèi)部介質(zhì)的泊松比μ ∈[0,0.5]，由式(11)可以清楚的看出，增大人工邊界單元的質(zhì)量密度將有助于改善數(shù)值穩(wěn)定性條件.

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

側(cè)邊子系統(tǒng)如圖4 所示，其模型參數(shù)與角點(diǎn)子系統(tǒng)相同.按節(jié)點(diǎn)編號組裝矩陣后，得到側(cè)邊子系統(tǒng)中5 號節(jié)點(diǎn)的運(yùn)動方程如下

參照上一節(jié)方法，將式(12)寫成傳遞矩陣的形式

圖4 側(cè)邊子系統(tǒng)Fig.4 Edge subsystem

其中

計(jì)算傳遞矩陣A2的特征值，并根據(jù)穩(wěn)定性判別條件ρ(A2)1，得出側(cè)邊子系統(tǒng)穩(wěn)定性條件的解析表達(dá)式為

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

其中γ2為側(cè)邊子系統(tǒng)穩(wěn)定性參數(shù)，整理后如式(18)所示.

由式(18)可知，側(cè)邊子系統(tǒng)穩(wěn)定性系數(shù)γ2同樣為一個無量綱的解析表達(dá)式，可以用來定量討論質(zhì)量密度比/ρ 的變化對數(shù)值穩(wěn)定性條件的影響.

3 數(shù)值計(jì)算穩(wěn)定性條件的改善方法

以上通過理論推導(dǎo)，分別給出了角點(diǎn)和側(cè)邊子系統(tǒng)數(shù)值積分穩(wěn)定性條件，其數(shù)值穩(wěn)定性系數(shù)γ1和γ2均為質(zhì)量密度比/ρ、內(nèi)部介質(zhì)泊松比μ、單元尺寸與散射波源至邊界距離之比L/R的函數(shù).當(dāng)無量綱參數(shù)?tCP/L小于或等于穩(wěn)定性系數(shù)時，相應(yīng)的子系統(tǒng)滿足數(shù)值積分的穩(wěn)定性條件.此外，整體系統(tǒng)的數(shù)值穩(wěn)定性除需滿足人工邊界區(qū)的角點(diǎn)子系統(tǒng)和側(cè)邊子系統(tǒng)的穩(wěn)定條件外，還需滿足內(nèi)部系統(tǒng)的穩(wěn)定性條件.顯式時域逐步積分算法中內(nèi)部單元的穩(wěn)定性條件為[25]

其中，γ0=1.

3.1 人工邊界單元質(zhì)量密度對穩(wěn)定性條件的影響

圖5 和圖6 分別給了出角點(diǎn)子系統(tǒng)和側(cè)邊子系統(tǒng)數(shù)值積分穩(wěn)定性系數(shù)γ1和γ2隨單元質(zhì)量密度比的變化曲線.其中單元尺寸與散射波源至邊界距離之比L/R分別為1 和1/100，內(nèi)部介質(zhì)泊松比μ分別為0，0.25，0.3，0.4 和0.5.由圖5 和圖6 可以看出，增大質(zhì)量密度比可以提高穩(wěn)定性系數(shù).對于角點(diǎn)和側(cè)邊兩種子系統(tǒng)，內(nèi)部介質(zhì)泊松比μ的取值越大，穩(wěn)定性越寬松.

圖5 角點(diǎn)子系統(tǒng)穩(wěn)定性系數(shù)隨密度比變化情況Fig.5 The variation of corner subsystem stability coefficient with density ratio

圖6 側(cè)邊子系統(tǒng)穩(wěn)定性系數(shù)隨密度比變化情況Fig.6 The variation of edge subsystem stability coefficient with density ratio

整體系統(tǒng)的數(shù)值穩(wěn)定條件應(yīng)同時滿足人工邊界子系統(tǒng)和內(nèi)部單元的穩(wěn)定性條件.圖7 給出了當(dāng)L/R=1 和μ=0.25 時，整體系統(tǒng)的數(shù)值積分穩(wěn)定條件.圖中陰影區(qū)域即為滿足整體系統(tǒng)數(shù)值積分穩(wěn)定性的區(qū)域.

圖7 整體系統(tǒng)數(shù)值穩(wěn)定區(qū)域Fig.7 The region of numerical stability in the overall system

以上分析可以得出以下兩點(diǎn)結(jié)論：(1)人工邊界區(qū)的數(shù)值積分穩(wěn)定性條件由角點(diǎn)子系統(tǒng)控制；(2)采用黏彈性人工邊界單元時，增大邊界單元的質(zhì)量密度可以有效改善顯式數(shù)值積分的穩(wěn)定性條件.

3.2 距離R 對穩(wěn)定性條件的影響

影響人工邊界子系統(tǒng)穩(wěn)定性的另外一個因素是散射波源至人工邊界的距離R.由于穩(wěn)定性系數(shù)計(jì)算式(11)和式(18)的分子和分母中均含有反映距離R影響的無量綱系數(shù)L/R，因此，難以像質(zhì)量密度一樣，直接給出對穩(wěn)定性影響效果的明確判斷.圖8 和圖9 給出了泊松比等于0.25，質(zhì)量密度比分別為0 和1 時，3 種系統(tǒng)穩(wěn)定性系數(shù)γ0，γ1和γ2隨R/L的變化曲線.

圖8 穩(wěn)定性系數(shù)隨R/L的變化情況Fig.8 The variation of stability coefficient with R/L

從圖8 中可以看出，當(dāng)散射波源至人工邊界的距離與單元尺寸之比R/L增大時，角點(diǎn)和側(cè)邊兩種子系統(tǒng)的穩(wěn)定性參數(shù)只在初始階段略為增大，當(dāng)R大于約5 倍的L后，穩(wěn)定性系數(shù)基本保持不變.因此，增大模型尺寸對整體系統(tǒng)數(shù)值穩(wěn)定性的改善幾乎不產(chǎn)生影響.另外，當(dāng)泊松比μ取為其他值時，可以得到完全相同的結(jié)論.

3.3 改進(jìn)的黏彈性人工邊界單元

基于以上分析，可確定改善顯式數(shù)值算法穩(wěn)定性具體措施(方法)和原則：(1)通過增大黏彈性人工邊界單元的質(zhì)量密度達(dá)到改善算法穩(wěn)定性的目的；(2)保證人工邊界區(qū)子系統(tǒng)的穩(wěn)定性條件達(dá)到或略寬于內(nèi)部單元系統(tǒng)穩(wěn)定性條件，使得整體系統(tǒng)的穩(wěn)定性由內(nèi)部系統(tǒng)控制和判斷.

改進(jìn)的黏彈性人工邊界單元參數(shù)設(shè)置如下：設(shè)內(nèi)部介質(zhì)剪切模量為G，質(zhì)量密度為ρ，剪切波速和壓縮波速分別為CS和CP；內(nèi)部介質(zhì)單元厚度和邊界單元厚度均為L，邊界節(jié)點(diǎn)至散射波源的距離為R，邊界單元質(zhì)量密度為；αT與αN分別為切向與法向黏彈性人工邊界參數(shù)，二維模型建議αT取0.5，αN取1.0.穩(wěn)定性條件改善后，二維黏彈性人工邊界單元質(zhì)量密度、等效剛度、泊松比和等效阻尼分別為

穩(wěn)定性改善后，整體系統(tǒng)穩(wěn)定性條件基本受內(nèi)部介質(zhì)區(qū)域控制，可采用最小單元尺寸與內(nèi)部介質(zhì)壓縮波速的比值估算滿足穩(wěn)定性條件的臨界時間步長.

4 穩(wěn)定性數(shù)值算例及精度分析

4.1 均勻半空間模型

建立如圖9 所示的均勻半空間模型.模型尺寸為320 m×160 m，內(nèi)部介質(zhì)的密度為2500 kg/m3，剪切波速為500 m/s，泊松比為0.3，有限元網(wǎng)格尺寸為2 m×2 m，模型側(cè)邊和底邊設(shè)置黏彈性人工邊界單元，選取圖中A,B,C,D四點(diǎn)作為觀測點(diǎn)，考察邊界單元參數(shù)變化后，模型節(jié)點(diǎn)位移的計(jì)算精度.采用持時為0.2 s 的脈沖作為動力載荷，沿y軸正向施加于模型中點(diǎn)處，載荷時程曲線如圖10 所示.

圖9 均勻半空間模型Fig.9 Homogeneous half space model

圖10 動力載荷時程曲線Fig.10 Time history of the dynamic load

表1 給出了均勻半空間模型穩(wěn)定性改善前后臨界時間步長比較情況.穩(wěn)定性條件改善前，邊界單元的質(zhì)量密度為0，整體系統(tǒng)穩(wěn)定性受角點(diǎn)子系統(tǒng)控制，按式(9)計(jì)算得到角點(diǎn)子系統(tǒng)滿足穩(wěn)定性條件的臨界時間步長為0.62 ms.按照此時間步長可順利完成模型動力計(jì)算.而如果將時間步長增大，例如取0.85 ms，計(jì)算發(fā)生失穩(wěn).

表1 均勻半空間模型穩(wěn)定性條件改善前后臨界時間步長比較Table 1 Comparison of critical time steps before and after the improvement of the stability in the homogeneous half space model

采用改進(jìn)的黏彈性人工邊界單元，即將人工邊界單元密度設(shè)為與內(nèi)部介質(zhì)一致，其值為2500 kg/m3，整體模型穩(wěn)定性變?yōu)槭軆?nèi)部區(qū)域控制，采用內(nèi)部介質(zhì)穩(wěn)定性條件(?t=2.14 ms)可順利完成對整體模型的顯式時域逐步積分計(jì)算，不會發(fā)生失穩(wěn).從圖11 可以看出，穩(wěn)定性條件改善后，觀測點(diǎn)A,B,C,D四點(diǎn)的位移時程曲線與改善前相比吻合良好，說明邊界單元自身質(zhì)量密度的適當(dāng)增大對波動問題計(jì)算精度的影響很小.另一方面，由于數(shù)值積分的時間步長比改善前增大了243.3%，總計(jì)算時間相較于原時長縮短了70.9%左右，大幅提高了顯式時域逐步積分的計(jì)算效率.

圖11 均勻介質(zhì)模型中4 個觀測點(diǎn)y 方向位移比較Fig.11 Comparison of displacements in y direction of four nodes in the homogeneous model

4.2 成層半空間模型

建立如圖12 所示的成層半空間有限元模型.模型尺寸為320 m×160 m，上層介質(zhì)密度為2000 kg/m3,剪切波速為300 m/s，泊松比為0.27，下層介質(zhì)密度為2500 kg/m3，剪切波速為500 m/s，泊松比為0.3，整體模型的有限元網(wǎng)格尺寸為2 m×2 m.脈沖載荷施加于模型的幾何中心，沿y軸正向加載，時程曲線如圖6 所示.通過令邊界單元的質(zhì)量密度與內(nèi)部區(qū)域相等的方法對穩(wěn)定性進(jìn)行改善.同樣選取圖中A,B,C,D四點(diǎn)作為觀測點(diǎn)，考察穩(wěn)定性改善方法的可靠性和改善后半無限域空間波動問題的計(jì)算精度.

圖12 成層半空間模型Fig.12 Layered half space model

表2 對成層介質(zhì)中各子系統(tǒng)和內(nèi)部區(qū)域穩(wěn)定性改善前后的臨界時間步長進(jìn)行了比較.從結(jié)果可以看出，上層介質(zhì)的穩(wěn)定性條件比下層介質(zhì)寬松.由于穩(wěn)定性條件受最不利情況控制，因此成層介質(zhì)模型整體系統(tǒng)的穩(wěn)定性條件由下層介質(zhì)所控制.本算例中的下層介質(zhì)參數(shù)與4.1 節(jié)的均勻半空間模型相同，因此整體系統(tǒng)穩(wěn)定性條件也與上節(jié)算例相同.

表2 成層半空間模型穩(wěn)定性條件改善前后臨界時間步長比較Table 2 Comparison of critical time steps before and after the improvement of the stability in the layered half space model

按照穩(wěn)定性改善前后的臨界時間步長對成層半空間模型進(jìn)行動力分析，結(jié)果如圖13 所示.圖中A,B,C,D四點(diǎn)位移時程曲線吻合良好，再次驗(yàn)證了邊界單元自身質(zhì)量密度的適當(dāng)增大對波動問題計(jì)算精度的影響很小.對于成層半空間模型，顯式時域逐步積分的總計(jì)算時間相較于原時長也縮短了70.9%，計(jì)算效率顯著提高.

圖13 成層介質(zhì)模型中四個觀測點(diǎn)y 方向位移比較Fig.13 Comparison of displacements in y direction of four nodes in the layered model

5 結(jié)論

本文針對采用黏彈性人工邊界單元時顯式時域逐步積分算法穩(wěn)定性改善問題，推導(dǎo)得到了含黏彈性人工邊界單元在內(nèi)的有限元系統(tǒng)在采用顯式時域逐步積分(中心差分)時數(shù)值穩(wěn)定條件的解析解.研究了各物理參數(shù)對穩(wěn)定性條件的影響，提出了改善數(shù)值積分穩(wěn)定性的方法，并通過算例驗(yàn)證了這一改善方法的有效性和對黏彈性人工邊界單元計(jì)算精度的影響.結(jié)論如下：

(1)采用有適當(dāng)質(zhì)量的黏彈性人工邊界單元，可在保證無限域空間波動問題計(jì)算精度的前提下，大幅增大滿足顯式時域逐步積分穩(wěn)定性條件的臨界時間步長，達(dá)到改善顯式數(shù)值積分穩(wěn)定性的目的.

(2) 可將內(nèi)部單元密度設(shè)置為邊界單元密度的上限，此時整體系統(tǒng)的積分穩(wěn)定性條件受內(nèi)部介質(zhì)區(qū)域控制，可通過模型中最小單元尺寸與內(nèi)部介質(zhì)壓縮波速的比值來估算滿足穩(wěn)定性條件的臨界時間步長.

(3)黏彈性人工邊界單元中的物理參數(shù)R，即散射波源至人工邊界的距離，對顯式時域逐步積分算法穩(wěn)定性的影響不大，當(dāng)R很小時，R的增大可以略微改善顯式算法的穩(wěn)定性，當(dāng)R大于約5 倍的單元邊長時，其對顯式算法穩(wěn)定性的影響基本不再變化.

(4)如果顯式時域逐步積分格式發(fā)生變化，依然可以使用本文提出的傳遞矩陣譜半徑分析方法，按照不同積分格式將節(jié)點(diǎn)子系統(tǒng)運(yùn)動方程展開，利用譜半徑不大于1 的判別條件給出穩(wěn)定性結(jié)論，研究分析適用的穩(wěn)定性改善方法.