李述濤，劉晶波，寶 鑫，陶西貴，陳一村，肖 蘭，賈藝凡

(1. 清華大學(xué)土木工程系，北京 100084；2. 軍事科學(xué)院國防工程研究院，北京 100036)

采用數(shù)值方法計算土-結(jié)構(gòu)地震反應(yīng)或近場波動問題時，需要從半無限介質(zhì)中切取出有限的近場區(qū)域進(jìn)行計算，同時需在截斷邊界處設(shè)置人工邊界來模擬外行散射波向無窮遠(yuǎn)的輻射效應(yīng)[1 ? 4]。較為常見的粘彈性人工邊界[5]計算精度較高但前處理操作較為復(fù)雜，需遍歷所有邊界節(jié)點(diǎn)施加彈簧-阻尼器系統(tǒng)，在此基礎(chǔ)上發(fā)展了粘彈性人工邊界單元技術(shù)[6 ? 7]，它是沿邊界法線向外延伸的一層單元，通過賦予該層單元等效物理參數(shù)即可模擬粘彈性人工邊界對于外行散射波的吸收作用。由于該方法具有良好的計算精度和魯棒性，且前處理過程簡單，在實(shí)際工程計算中得到較多的應(yīng)用[8 ? 12]。

隨著人工邊界技術(shù)的不斷發(fā)展，針對人工邊界條件的穩(wěn)定性研究逐漸得到重視[13 ? 14]，目前針對此類問題的主要研究方法包括：Trefethen 等[15]提出的GKS 定理—基于偏微分方程初邊值問題有限差分格式的穩(wěn)定性理論，給出初邊值問題多層線性差分格式穩(wěn)定的充要條件；Liao 等[16]研究了離散模型中穩(wěn)態(tài)波動的完備解，給出了多次透射邊界的穩(wěn)定性條件；Kamel 等[17]和關(guān)慧敏等[18]通過傳遞矩陣譜半徑分析積分格式的穩(wěn)定性。理論研究表明，如果不全面考慮逐步積分過程中人工邊界節(jié)點(diǎn)與內(nèi)部節(jié)點(diǎn)運(yùn)動方程的耦合效應(yīng)，穩(wěn)定性準(zhǔn)則可能會失效。

由于粘彈性人工邊界或粘彈性人工邊界單元本身均為物理上穩(wěn)定的系統(tǒng)，將其引入計算系統(tǒng)中，并不影響整體的物理穩(wěn)定條件。但采用顯式時域逐步積分方法進(jìn)行計算時，受粘彈性人工邊界單元粘性阻尼、剛度及幾何尺寸的影響，整體計算模型的數(shù)值積分穩(wěn)定性條件將發(fā)生改變，目前尚未給出明確而實(shí)用的含粘彈性人工邊界條件影響的顯式算法穩(wěn)定性準(zhǔn)則，影響了數(shù)值分析時穩(wěn)定時間步長的判斷和選取，進(jìn)一步限制了粘彈性人工邊界單元在顯式動力分析中的應(yīng)用。因此目前在引入粘彈性人工邊界單元進(jìn)行波動問題分析時，多采用隱式的無條件穩(wěn)定的時域逐步積分算法，以避免顯式算法帶來的數(shù)值穩(wěn)定性問題。隱式算法需要求解聯(lián)立方程組，計算效率不高，并不適合大規(guī)模波動問題計算。隨著顯式時域逐步積分算法在工程結(jié)構(gòu)和大范圍場地分析問題中的廣泛應(yīng)用[19 ? 22]，有必要對含有粘彈性人工邊界單元的系統(tǒng)進(jìn)行顯式時域逐步積分算法穩(wěn)定性研究。

本文考慮人工邊界和內(nèi)部單元節(jié)點(diǎn)的運(yùn)動耦合效應(yīng)，利用傳遞矩陣的譜半徑分析時域逐步積分格式(中心差分)的穩(wěn)定性，給出不同位置處局部節(jié)點(diǎn)系統(tǒng)顯式時域逐步積分算法穩(wěn)定條件的解析解，通過各子系統(tǒng)穩(wěn)定性條件及內(nèi)部介質(zhì)穩(wěn)定性條件的對比分析，給出考慮粘彈性人工邊界條件影響的整體耦合系統(tǒng)顯式時域逐步積分算法的穩(wěn)定性條件。

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

粘彈性人工邊界單元是在模型截斷邊界處沿法線向外延伸的一層單元(見圖1)，人工邊界單元的厚度為h，質(zhì)量密度為0，單元最外層節(jié)點(diǎn)固定。通過賦予單元等效物理參數(shù)來模擬粘彈性人工邊界。二維粘彈性人工邊界單元等效彈性模量、等效泊松比和等效阻尼系數(shù)分別為[2]：

圖1 人工邊界單元尺寸示意圖Fig. 1 Artificial boundary elements size diagram

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

盡管粘彈性人工邊界單元等效剛度矩陣的推導(dǎo)過程引入了邊界單元厚度h遠(yuǎn)小于寬度L的假設(shè)，如圖1 所示，劉晶波等[6]和谷音等[7]通過進(jìn)一步的研究表明，邊界單元厚度的魯棒性良好，可靈活取值而對精度影響很小。由于顯式算法的穩(wěn)定條件受模型中的最小單元尺寸制約，在實(shí)際計算分析中，建議將邊界單元的寬度設(shè)定為與內(nèi)部單元尺寸一致，以排除邊界單元尺寸對整體穩(wěn)定性的影響，如圖2 所示。此時粘彈性人工邊界單元等效彈性模量、等效剪切模量可以寫為：

圖2 適用于顯式算法的人工邊界單元Fig. 2 Artificial boundary elements adapted to explicit algorithms

使用通用有限元軟件對粘彈性人工邊界單元賦予材料參數(shù)時，由于是各向同性介質(zhì)，輸入等效彈性模量即可。另外還包括式(1)計算得到的等效阻尼系數(shù)，以及泊松比(=0)和質(zhì)量密度(=0)。由于大多數(shù)軟件不支持將密度設(shè)為0，可采用近似0 的小數(shù)代替。

2 顯式時域逐步積分算法穩(wěn)定性分析方法

時域逐步積分算法穩(wěn)定性分析的目的是獲得時域逐步積分計算時滿足穩(wěn)定性要求的時間離散步長 ?t。最大的穩(wěn)定時間步長 ?t與單元的尺寸和系統(tǒng)的物理性質(zhì)有關(guān)。為此在進(jìn)行穩(wěn)定性分析時，通常采用均勻介質(zhì)和均勻離散化網(wǎng)格模型，如圖3 所示。

圖3 均勻離散化網(wǎng)格模型Fig. 3 Uniform discretization mesh model

當(dāng)實(shí)際力學(xué)模型介質(zhì)非均勻，單元尺寸不相同時，可以根據(jù)最小尺寸的單元或綜合考慮單元尺寸及其介質(zhì)性質(zhì)來確定整體計算模型的穩(wěn)定離散時間步長 ?t。算法的穩(wěn)定性通?？刹捎靡韵聨追N分析方法：

① 不考慮自由邊界和人工邊界的影響，取無限計算區(qū)域模型，采用馮諾依曼方法進(jìn)行穩(wěn)定性分析。這一方法可以獲得內(nèi)部系統(tǒng)的數(shù)值計算穩(wěn)定條件，但無法考慮含人工邊界條件時系統(tǒng)的穩(wěn)定性。

② 取實(shí)際的計算模型，考慮人工邊界的影響，通過對整體系統(tǒng)時域逐步積分方法的傳遞矩陣進(jìn)行特征值分析，以獲得考慮耦合人工邊界條件影響的整體系統(tǒng)的數(shù)值計算穩(wěn)定性條件。但由于涉及到整體模型傳遞矩陣的建立和分析，這一方法在實(shí)際工作中通常是不可行的，特別是對大規(guī)模的計算問題。

③ 為對透射人工邊界的數(shù)值計算穩(wěn)定性進(jìn)行有效分析，Kamel 等[17]和關(guān)慧敏等[18]提出了一種對含透射人工邊界條件子系統(tǒng)的傳遞矩陣進(jìn)行特征值分析，以獲得波動問題中透射人工邊界穩(wěn)定性的分析方法。由于分析中采用的子系統(tǒng)為沿人工邊界切向(寬度方向)上若干排節(jié)點(diǎn)和沿法向(深度方向)上若干排節(jié)點(diǎn)組成的節(jié)點(diǎn)系，子系統(tǒng)自由度較多，僅能靠數(shù)值方法對人工邊界的穩(wěn)定性進(jìn)行分析和判斷，未給出具有解析形式的更易于使用的人工邊界穩(wěn)定性準(zhǔn)則。

從以上有限元離散模型數(shù)值積分穩(wěn)定性準(zhǔn)則的研究工作中可以發(fā)現(xiàn)以下兩個特點(diǎn)：1)算法的穩(wěn)定性與模型的截止頻率(系統(tǒng)最高階自振頻率)有關(guān)。2)截止頻率相應(yīng)的振型呈現(xiàn)局部節(jié)點(diǎn)系相鄰節(jié)點(diǎn)交錯振動的形態(tài)。由以上兩個特點(diǎn)可以判斷，整體系統(tǒng)的截止頻率可通過對局部子系統(tǒng)模型的分析得到，進(jìn)一步可通過局部子系統(tǒng)模型的數(shù)值穩(wěn)定性分析，獲得整體系統(tǒng)的數(shù)值穩(wěn)定性條件。

首先以一維剪切梁問題為例證明以上判斷。圖4 為一維剪切梁的離散模型及最高階振型示意圖，梁的剪切模量為G，密度為 ρ，橫截面面積為A，單元長度為L。取兩個子系統(tǒng)進(jìn)行分析。子系統(tǒng)a：取振型兩節(jié)點(diǎn)之間的子系統(tǒng)，兩端施加約束，由于約束是施加于振型節(jié)點(diǎn)(振幅為零)之上，因此并不改變剪切梁有限元模型的自振頻率，如圖5(a)所示；子系統(tǒng)b：取兩個單元，在兩端施加約束，此時子系統(tǒng)的自振頻率與原系統(tǒng)的截止頻率不相同，如圖5(b)所示。

圖4 一維剪切梁振動示意圖Fig. 4 One-dimensional shear beam vibration diagram

圖5 一維局部子系統(tǒng)Fig. 5 One-dimensional local subsystem

根據(jù)子系統(tǒng)的剛度和質(zhì)量，可以建立子系統(tǒng)a的單自由度運(yùn)動方程，當(dāng)采用中心差分法進(jìn)行時域逐步積分計算時，一維子系統(tǒng)a 的數(shù)值穩(wěn)定條件為：

其中，Tn為子系統(tǒng)的自振周期，將式(3)中第三式代入式(4)，得到子系統(tǒng)數(shù)值積分穩(wěn)定性條件為：

式(5)即為一維連續(xù)介質(zhì)波動有限元分析時中心差分算法的穩(wěn)定性條件，可見采用子系統(tǒng)a 獲得的局部穩(wěn)定性條件即為連續(xù)介質(zhì)整體模型的穩(wěn)定性條件。

一維子系統(tǒng)b 的剛度kb、質(zhì)量mb、自振頻率ωb和中心差分穩(wěn)定性條件分別為：

對比式(5)和式(6)，可以發(fā)現(xiàn)一維子系統(tǒng)b的穩(wěn)定性條件比實(shí)際整體模型的穩(wěn)定條件更為寬松，但仍然可以給出穩(wěn)定時間積分步長 ?t的一個上限估計。

圖6 二維節(jié)點(diǎn)系統(tǒng)最高階振型圖Fig. 6 Highest-order modal pattern of two-dimensional nodes system

同樣取兩個子系統(tǒng)進(jìn)行分析，二維問題局部子系統(tǒng)由4 個單元構(gòu)成，如圖7(a)和圖7(b)所示。其中二維子系統(tǒng)a 的4 個角點(diǎn)為振型節(jié)點(diǎn)，位移為0，4 個邊點(diǎn)的位移條件可由振型和4 邊形等參元的性質(zhì)確定，子系統(tǒng)a 對應(yīng)的邊界條件為：

圖7 二維局部子系統(tǒng)Fig. 7 Two-dimensional local subsystem

二維子系統(tǒng)a 自振頻率和中心差分穩(wěn)定條件為：

式(9)即為二維波動問題有限元分析時中心差分算法的穩(wěn)定性條件，可見同樣可以由子系統(tǒng)a的局部穩(wěn)定性分析獲得二維有限元模型數(shù)值算法的整體穩(wěn)定性條件。

二維子系統(tǒng)b 由4 個單元構(gòu)成，在周邊節(jié)點(diǎn)施加固定約束，子系統(tǒng)b 運(yùn)動方程為：

同樣，二維子系統(tǒng)b 給出了整體有限元模型穩(wěn)定時間積分步長的上限估計。

一維和二維子系統(tǒng)分析的結(jié)果表明：1)若能較為準(zhǔn)確地判斷系統(tǒng)截止頻率對應(yīng)的振動模態(tài)(振型)，則可以利用最高振型的特點(diǎn)和振型節(jié)點(diǎn)的分布規(guī)律，從整體系統(tǒng)中分離出由陣型節(jié)點(diǎn)所包圍的局部子系統(tǒng)，對振型節(jié)點(diǎn)施加約束條件，然后對該子系統(tǒng)進(jìn)行穩(wěn)定性分析，獲得的局部子系統(tǒng)穩(wěn)定性條件即為整體模型的數(shù)值穩(wěn)定條件；2)若不能準(zhǔn)確判斷截止頻率對應(yīng)振動模態(tài)和相應(yīng)振型節(jié)點(diǎn)的位置，也可選取一個能體現(xiàn)整體有限元模型特征的最小的子系統(tǒng)，對該子系統(tǒng)的邊界節(jié)點(diǎn)施加約束，通過分析獲得子系統(tǒng)的數(shù)值穩(wěn)定性條件，從而獲得整體系統(tǒng)穩(wěn)定性判據(jù)中各物理量的關(guān)系式，且該條件是整體系統(tǒng)穩(wěn)定性條件的一個逼近。再通過擴(kuò)展數(shù)值算例進(jìn)行修正，確定合適的系數(shù)，即可得到整體系統(tǒng)的數(shù)值積分穩(wěn)定性條件。

當(dāng)采用粘彈性人工邊界單元時，內(nèi)部系統(tǒng)的數(shù)值穩(wěn)定條件已知，但難以確定人工邊界區(qū)截止頻率對應(yīng)的振型，因而僅能采用類似于b 型的子系統(tǒng)進(jìn)行數(shù)值計算的穩(wěn)定性分析。

對非均質(zhì)子系統(tǒng)進(jìn)行穩(wěn)定性分析時，可根據(jù)顯式逐步積分算法格式，將系統(tǒng)運(yùn)動方程寫成如下形式：

積分格式的穩(wěn)定性問題與外力向量Qp無關(guān)。如果滿足下列兩條件，則積分格式(12)是穩(wěn)定[23]：條件1： ρ(A)≤1 ， ρ(A)是傳遞矩陣A的譜半徑，即ρ(A)=max|λi|；條件2：如果A具有多重特征值，則該特征值的模小于1。因而可將子系統(tǒng)的穩(wěn)定性分析歸結(jié)為傳遞矩陣A的形成及譜半徑計算，按上述條件即可建立穩(wěn)定性準(zhǔn)則。

3 人工邊界子系統(tǒng)穩(wěn)定性分析

采用粘彈性人工邊界單元時，人工邊界區(qū)的子系統(tǒng)有兩種形式，一種在側(cè)面或底面切取，如圖8 所示，另一種在角點(diǎn)處切取，切取邊界可設(shè)置固定邊界(見圖9)或自由邊界(見圖10)。

圖8 側(cè)邊固定邊界子系統(tǒng)Fig. 8 Side fixed boundary subsystem

圖9 角點(diǎn)固定邊界子系統(tǒng)Fig. 9 Corner fixed boundary subsystem

圖10 角點(diǎn)自由邊界子系統(tǒng)Fig. 10 Corner fixed boundary subsystem

3.1 側(cè)邊固定邊界子系統(tǒng)穩(wěn)定性分析

側(cè)邊固定邊界子系統(tǒng)由兩個粘彈性人工邊界單元和兩個內(nèi)部介質(zhì)單元組成，如圖8 所示。設(shè)內(nèi)部介質(zhì)單元的彈性模量為E，剪切模量為G，密度為 ρ，泊松比為μ，且無阻尼；粘彈性人工邊界單元彈性模量為，阻尼系數(shù)為，泊松比為0，密度為0，單元邊長均為L。四節(jié)點(diǎn)正方形平面單元的集中質(zhì)量矩陣、剛度矩陣、阻尼矩陣分別見文獻(xiàn)[6]和文獻(xiàn)[24]，按節(jié)點(diǎn)編號進(jìn)行矩陣組裝后得到子系統(tǒng)中5 號節(jié)點(diǎn)運(yùn)動方程如下：

觀察式(13)可知，x和y方向運(yùn)動方程是解耦的且形式相同，因此可只對其中一個方向的運(yùn)動方程進(jìn)行穩(wěn)定性分析。顯式時域逐步積分格式(中心差分[25]如下：

式(19)即為側(cè)邊局部人工邊界子系統(tǒng)穩(wěn)定性條件的解析表達(dá)形式，觀察可知該子系統(tǒng)的穩(wěn)定性條件與內(nèi)部介質(zhì)壓縮波速、泊松比、單元尺寸和模型大小有關(guān)。

3.2 角點(diǎn)固定邊界子系統(tǒng)穩(wěn)定性分析

角點(diǎn)處固定邊界子系統(tǒng)中由3 個粘彈性人工邊界單元和1 個內(nèi)部介質(zhì)單元組成，如圖9 所示。將四節(jié)點(diǎn)正方形平面單元質(zhì)量、阻尼和剛度矩陣按節(jié)點(diǎn)編號組裝后得到5 號節(jié)點(diǎn)的運(yùn)動方程如下：

由于該運(yùn)動方程的坐標(biāo)耦聯(lián)，需對整體方程進(jìn)行穩(wěn)定性分析。將式(1)代入式(20)，然后根據(jù)式(14)將式(20)展開，寫成如下傳遞矩陣的形式：

其中：

式(24)即為角點(diǎn)處局部系統(tǒng)穩(wěn)定性條件的解析表達(dá)形式，觀察可知影響角點(diǎn)局部人工邊界子系統(tǒng)穩(wěn)定性的參數(shù)與側(cè)邊相同，但各自貢獻(xiàn)不同。

3.3 角點(diǎn)自由邊界子系統(tǒng)穩(wěn)定性分析

為進(jìn)一步研究角點(diǎn)處自由邊界子系統(tǒng)的穩(wěn)定性，如圖10 所示。研究對象為編號為1、2、4、5 的四個節(jié)點(diǎn)。子系統(tǒng)中既有單元內(nèi)部節(jié)點(diǎn)，也包含邊界處節(jié)點(diǎn)。將四節(jié)點(diǎn)正方形平面單元質(zhì)量、阻尼和剛度矩陣按節(jié)點(diǎn)編號組裝后，可得到子系統(tǒng)的剛度、質(zhì)量和阻尼矩陣，均為8×8 階對稱矩陣。將整體運(yùn)動方程按照式(14)中心差分格式展開，得到16×16 階傳遞矩陣。該傳遞矩陣無法得出解析形式的穩(wěn)定性條件，可代入?yún)?shù)計算數(shù)值解?？紤]穩(wěn)定性條件受系統(tǒng)截止頻率影響，自由邊界子系統(tǒng)截止頻率要低于固定邊界子系統(tǒng)，由此判斷后者穩(wěn)定性條件應(yīng)比前者嚴(yán)格，下一節(jié)將進(jìn)行驗證。

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

為比較3 種人工邊界子系統(tǒng)的數(shù)值計算穩(wěn)定條件，采用兩組參數(shù)進(jìn)行分析，具體參數(shù)見表1。

表1 兩組參數(shù)值(無量綱)Table 1 Two sets of parameter values (dimensionless)

側(cè)邊固定和角點(diǎn)固定子系統(tǒng)的穩(wěn)定條件可由式(19)和式(24)獲得，角點(diǎn)自由邊界子系統(tǒng)的穩(wěn)定條件可采用數(shù)值方法獲得。3 種子系統(tǒng)譜半徑的計算結(jié)果見圖11。當(dāng)譜半徑小于等于1 時，數(shù)值計算滿足穩(wěn)定性條件，由圖11 可以直觀地比較兩組參數(shù)條件下3 種子系統(tǒng)的穩(wěn)定條件。

表2 中對3 種局部子系統(tǒng)和內(nèi)部系統(tǒng)滿足穩(wěn)定性條件的臨界時間步長進(jìn)行了定量比較。內(nèi)部系統(tǒng)的穩(wěn)定性條件( ?t=L/CP)未考慮粘彈性人工邊界單元對數(shù)值算法的影響，最為寬松；側(cè)邊和角點(diǎn)子系統(tǒng)考慮了粘彈性人工邊界單元質(zhì)量、等效剛度和阻尼對穩(wěn)定性的影響，穩(wěn)定性條件要比無人工邊界單元時更為嚴(yán)格。四周固定的角點(diǎn)處邊界節(jié)點(diǎn)只共享1/4 內(nèi)部單元質(zhì)量，其子系統(tǒng)截止頻率最高，因此穩(wěn)定性條件最為嚴(yán)格。以上3 種子系統(tǒng)數(shù)值積分穩(wěn)定條件的對比表明，采用粘彈性人工邊界單元時，系統(tǒng)數(shù)值積分的穩(wěn)定條件由角點(diǎn)區(qū)控制。

圖11 三種子系統(tǒng)的譜半徑對比Fig. 11 Comparison of spectral radius of three subsystems

表2 穩(wěn)定性條件比較(無量綱)Table 2 Comparison of stability conditions (dimensionless)

5 算例驗證

5.1 均勻半空間模型

為驗證以上穩(wěn)定性分析的準(zhǔn)確性，按第一組參數(shù)建立有限元近場模型，模型尺寸為320×160，內(nèi)部介質(zhì)密度為2500，剪切波速為500，泊松比為0.3，網(wǎng)格尺寸為2×2，模型側(cè)邊和底邊最外層單元是粘彈性人工邊界單元，如圖12 所示。采用持時為0.2 的脈沖作為動力荷載，施加于模型中點(diǎn)處，時程曲線如圖13 所示。分別按照內(nèi)部介質(zhì)數(shù)值穩(wěn)定條件(?t=0.0021)、側(cè)邊子系統(tǒng)穩(wěn)定條件(?t=0.00163)、固定邊界角點(diǎn)子系統(tǒng)穩(wěn)定條件(?t=0.00623)，采用固定時間步長的顯示時域逐步積分算法對整體模型進(jìn)行計算。

圖12 均勻半空間算例模型圖Fig. 12 Homogeneous half apace example model diagram

圖13 動力荷載時程曲線Fig. 13 Dynamic load time history curve

按照內(nèi)部介質(zhì)數(shù)值穩(wěn)定條件(?t=0.0021)計算時，由于不滿足側(cè)邊(底邊)子系統(tǒng)穩(wěn)定條件，因此P 波傳播到距離波源最近的底邊節(jié)點(diǎn)時發(fā)生失穩(wěn)，如圖14 所示；按照側(cè)邊子系統(tǒng)穩(wěn)定條件(?t=0.00163)計算時，由于不滿足角點(diǎn)處子系統(tǒng)穩(wěn)定條件，波動傳播到模型角點(diǎn)處時發(fā)生失穩(wěn)，如圖15所示，其中U為介質(zhì)中的位移，以上兩種穩(wěn)定條件均無法完成整體有限元模型的計算。

圖14 按內(nèi)部穩(wěn)定條件計算時底邊失穩(wěn)狀態(tài)(0.11 時刻)Fig. 14 The unstable state of the bottom edge when calculated according to internal stability conditions (time 0.11)

按照角點(diǎn)處子系統(tǒng)穩(wěn)定條件(?t=0.000623)計算時，可順利完成整體模型的動力顯式計算，粘彈性人工邊界單元很好模擬了外行波向無窮遠(yuǎn)的輻射，結(jié)果如圖16 所示。此外，通過進(jìn)一步的計算分析發(fā)現(xiàn)，當(dāng)整體穩(wěn)定條件取值略大于角點(diǎn)子系統(tǒng)穩(wěn)定條件時(例如?t=0.00085)，模型角點(diǎn)處也發(fā)生失穩(wěn)，失穩(wěn)狀態(tài)與圖15 相同。

圖15 按側(cè)邊子系統(tǒng)穩(wěn)定條件計算時角點(diǎn)失穩(wěn)狀態(tài)(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 按角點(diǎn)子系統(tǒng)穩(wěn)定條件計算結(jié)果Fig. 16 Results calculated by the stability conditions of the corner subsystem

以上計算分析驗證了采用粘彈性人工邊界單元時，整體模型顯式數(shù)值積分算法的穩(wěn)定性由角點(diǎn)區(qū)域控制，式(24)給出的角點(diǎn)處子系統(tǒng)穩(wěn)定條件是整體模型數(shù)值穩(wěn)定的充分條件。

5.2 成層半空間模型

為滿足實(shí)際場地計算需要，對成層半空間算例進(jìn)行驗證。成層半空間有限元模型尺寸為320×160，上半部分內(nèi)部介質(zhì)密度2000，剪切波速300，泊松比0.27，下半部分內(nèi)部介質(zhì)密度為2500，剪切波速為500，泊松比為0.3，整體模型網(wǎng)格尺寸為2×2，模型側(cè)邊和底邊最外層單元是粘彈性人工邊界單元，如圖17 所示。脈沖荷載施加于模型中心，如圖13 所示。

表3 給出了成層半空間局部子系統(tǒng)和內(nèi)部系統(tǒng)滿足穩(wěn)定性條件的臨界時間步長。模型中的上層介質(zhì)只有側(cè)邊子系統(tǒng)，無角點(diǎn)子系統(tǒng)，下層介質(zhì)兩者均存在。

圖17 成層半空間算例模型圖Fig. 17 Layered half space example model diagram

表3 成層半空間穩(wěn)定性條件比較(無量綱)Table 3 Comparison of stability conditions in layered half space (dimensionless)

經(jīng)比較發(fā)現(xiàn)，上層介質(zhì)中滿足內(nèi)部系統(tǒng)穩(wěn)定條件的時間步長(?t=0.0037)大于側(cè)邊子系統(tǒng)穩(wěn)定時間步長(?t=0.0028)。而二者均比下層介質(zhì)內(nèi)部系統(tǒng)穩(wěn)定的時間步長 (?t=0.0021)要寬松。因此，當(dāng)采用上層介質(zhì)穩(wěn)定性條件對整體模型計算時，首先發(fā)生失穩(wěn)的是下層介質(zhì)的內(nèi)部系統(tǒng)，如圖18所示。

圖18 按上層介質(zhì)穩(wěn)定條件計算時下層介質(zhì)內(nèi)部系統(tǒng)失穩(wěn)狀態(tài)(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)

采用下層介質(zhì)內(nèi)部系統(tǒng)穩(wěn)定條件(?t=0.0021)計算時，整體系統(tǒng)失穩(wěn)狀態(tài)與圖14 相同；采用下層介質(zhì)側(cè)邊子系統(tǒng)穩(wěn)定條件計算時(?t=0.00163)，失穩(wěn)狀態(tài)和圖15 相同。

按照下層角點(diǎn)子系統(tǒng)穩(wěn)定條件(?t=0.000623)計算時，可順利完成整體模型的動力顯式計算，分層介質(zhì)中，粘彈性人工邊界單元也可很好地模擬外行波向無窮遠(yuǎn)的輻射，結(jié)果如圖19 所示。成層半空間算例進(jìn)一步驗證了采用粘彈性人工邊界單元時，整體模型顯式數(shù)值積分算法的穩(wěn)定性仍然由角點(diǎn)區(qū)域控制。

圖19 按下層介質(zhì)角點(diǎn)子系統(tǒng)穩(wěn)定條件計算結(jié)果Fig. 19 Results calculated by the stability conditions of the corner subsystem in lower media

6 結(jié)論與展望

本文將整體模型數(shù)值穩(wěn)定性問題合理轉(zhuǎn)移到若干局部子系統(tǒng)中，充分考慮粘彈性人工邊界單元和內(nèi)部單元節(jié)點(diǎn)間的運(yùn)動耦合效應(yīng)，通過傳遞矩陣譜半徑分析方法推導(dǎo)出各局部子系統(tǒng)顯式時域逐步積分(中心差分)數(shù)值穩(wěn)定條件的解析解和數(shù)值解。通過計算軟件驗證理論分析的可靠性。具體結(jié)論如下：

(1) 對于大規(guī)模數(shù)值計算問題，可選取局部的子系統(tǒng)并對其進(jìn)行顯式時域逐步積分算法的穩(wěn)定性分析，該穩(wěn)定性條件與整體系統(tǒng)穩(wěn)定性條件相同或近似。

(2) 采用粘彈性人工邊界單元時，受人工邊界單元質(zhì)量、剛度和阻尼影響，整體系統(tǒng)的穩(wěn)定性條件與內(nèi)部介質(zhì)的穩(wěn)定性條件不同，前者的穩(wěn)定性條件更為嚴(yán)格，需使用更小的積分步長以滿足整體系統(tǒng)的數(shù)值穩(wěn)定。

(3) 采用粘彈性人工邊界單元時，整體模型顯式數(shù)值積分算法的穩(wěn)定性由角點(diǎn)區(qū)域控制，本文給出了角點(diǎn)處子系統(tǒng)數(shù)值積分的穩(wěn)定性條件，該穩(wěn)定性條件是整體模型數(shù)值積分穩(wěn)定的充分條件。此外，本文給出的穩(wěn)定性條件是以正方形平面單元為對象推導(dǎo)的，同樣適用于矩形平面單元，由于系統(tǒng)穩(wěn)定條件受最小單元尺寸影響，可使用矩形單元的最小邊長作為參數(shù)計算穩(wěn)定性條件。

下一步展望：

(1) 對于三維粘彈性人工邊界單元，同樣可以利用本文提出的傳遞矩陣譜半徑分析方法對顯式時域逐步積分算法的穩(wěn)定性條件進(jìn)行分析，三維問題涉及的局部子結(jié)構(gòu)更為特殊，傳遞矩陣的生成和特征值求解更加復(fù)雜，需進(jìn)一步開展研究。

(2) 相比隱式算法，顯式算法的解耦特性對于求解大范圍復(fù)雜工程場地問題更有優(yōu)勢。本文的研究成果為在顯式算法中合理使用粘彈性人工邊界單元提供了理論依據(jù)?？稍诖嘶A(chǔ)上進(jìn)一步開展分析和研究，以改善使用粘彈性人工邊界單元時顯式算法的穩(wěn)定性，提高大范圍工程場地問題的計算效率。