崔高航 王兆亮 劉守花

(東北林業(yè)大學(xué)土木工程學(xué)院，黑龍江哈爾濱 150040)

0 引言

人工邊界條件又稱為透射或吸收邊界條件，理論上對連續(xù)介質(zhì)的模擬更加的精確;目的是保證能量在通過人工邊界的時候，能與連續(xù)介質(zhì)一樣能夠完全的透射或者吸收，而不會由于人工邊界的設(shè)置而發(fā)生反射。

ABAQUS［1］作為通用的有限元軟件，分析計算功能和模擬性能在國際上處于領(lǐng)先的地位，它的分析過程、單元模型和材料模型等的種類比較齊全。ABAQUS在分析彈性問題或者非線性的組合問題的時候，獲得的分析結(jié)果相對于其他的有限元軟件都是比較準(zhǔn)確的。本文基于ABAQUS軟件平臺和透射邊界理論，利用Microsoft Visual Studio 2008平臺的 Intel(R)Fortran Compiler 11.0.061編寫位移迭代程序與ABAQUS軟件的DISP邊界條件子程序接口鏈接，實(shí)現(xiàn)透射邊界條件的施加;同時，基于ABAQUS軟件中動力無限元理論，利用四節(jié)點(diǎn)無限單元引入阻尼(阻尼矩陣)給出無限元邊界模型;最后，考慮足夠遠(yuǎn)的尺寸，建立符合進(jìn)度要求遠(yuǎn)置邊界模型，并將其分析結(jié)果參考為精確解用于對比分析;最后建立相同尺寸的固定邊界模型。算例分析，驗(yàn)證了本文在成層土體中施加透射邊界方法的合理性和有效性。此外，通過與其他邊界模型數(shù)據(jù)結(jié)果對比分析，結(jié)果表明，通過本文方法所施加的透射邊界在處理環(huán)境振動問題中的精度要優(yōu)于其他邊界條件。

1 透射邊界

1.1 透射邊界定義

透射邊界是通過給出單向波動的一般表達(dá)式和人工視波速ca(通常取剪切波速，即ca=cs)，建立了透射理論和相應(yīng)的多次透射公式［2］(MTF)(如圖1所示)。理論依據(jù)是單側(cè)行波解的一般表達(dá)式，通過對波在穿越邊界的過程分析，獲得可行的位移邊界條件(方式:內(nèi)節(jié)點(diǎn)位移代替邊界節(jié)點(diǎn)位移)。透射邊界概念與理論是我國學(xué)者廖振鵬等于20世紀(jì)80年代初提出的，具有物理概念明確、透射效果好、計算簡單，且很容易和有限元等數(shù)值解法相結(jié)合等優(yōu)點(diǎn)，近年來得到廣泛的應(yīng)用。

圖1 多次透射（MTF）公式的幾何關(guān)系

N階多次透射公式(MTF)為:

1.2 透射邊界的穩(wěn)定性

任何人工邊界都存在著穩(wěn)定性問題。透射邊界的穩(wěn)定性問題也是保證近場波動數(shù)值模擬精度的關(guān)鍵問題。其中，高頻振蕩和低頻飄移是透射邊界模型精度主要的數(shù)值失穩(wěn)現(xiàn)象。

參考相關(guān)文獻(xiàn)采取處理措施如下:

1)消除高頻振蕩失穩(wěn)的措施。

［3］～［5］在有限元計算中增加與應(yīng)變速度相關(guān)的阻尼，通常采用瑞利阻尼(C=αM+βK)形式。而通用有限元軟件ABAQUS具有瑞利阻尼模式，在建立分析模型時，可按屬性參數(shù)直接設(shè)置。并要求在保證數(shù)值積分穩(wěn)定的情況下，β值盡量取的小一些。研究表明:施加相關(guān)阻尼的措施，可以較好的限制高頻振蕩失穩(wěn)，而且對低頻波動影響可以降低到一定的范圍，從而獲得較好的數(shù)值穩(wěn)定性。

2)消除零頻飄移失穩(wěn)的措施。

參考文獻(xiàn)［5］～［7］認(rèn)為透射邊界數(shù)值計算中的零頻和零波數(shù)分量進(jìn)入計算區(qū)，是導(dǎo)致人工邊界飄移失穩(wěn)的原因。在MTF計算中，通過引入一個遠(yuǎn)小于1的正數(shù)算子γ對多次透射公式(式(1))進(jìn)行修正。修正后的多次透射公式為:

該算子一般取γ≤0.05的正數(shù)，但不宜過小，主要結(jié)合經(jīng)驗(yàn)和試算確定。從實(shí)用的觀點(diǎn)來看，MTF具有普適性、可控性、時空解耦性，在保證模擬精度的情況下，能夠在計算機(jī)上穩(wěn)定實(shí)現(xiàn)。

1.3 透射邊界的模型實(shí)現(xiàn)

一般情況下，MTF的計算點(diǎn)x=－jcaΔt與離散節(jié)點(diǎn)x=－nΔx不一致，故需通過線性內(nèi)插方法獲得計算點(diǎn)位移。MTF的數(shù)值實(shí)現(xiàn)方法有兩種:空間內(nèi)插和時間內(nèi)插。本文采用第一種空間插值形式(見圖1)。

在工程應(yīng)用中，MTF一般采取一階、二階透射公式。經(jīng)過實(shí)踐證明，MTF的精度是可控的，二階透射公式已經(jīng)具有較高的工程可用精度。故本文在下面的算例中采用空間三點(diǎn)兩插值形式，取人工邊界區(qū)前n(透射公式階數(shù):n=2)時刻的位移計算人工邊界點(diǎn)位移。

圖2 二維土體模型人工邊界計算區(qū)

參考圖2，二階透射公式如下:

經(jīng)過插值后的二階透射公式變?yōu)?

本文通過考慮波源問題進(jìn)行建模(見圖3)，外行波場不通過波場分離獲得，全波場［2］即為穿越人工邊界的外行波場，MTF可直接由全位移表出。透射邊界在ABAQUS中的實(shí)現(xiàn)步驟如下:

1)劃分離散網(wǎng)格，獲得具有時步積分格式的內(nèi)節(jié)點(diǎn)運(yùn)動方程。要求空間步距為Δx，且應(yīng)滿足在有意義的波長內(nèi)含有Δx的數(shù)量不應(yīng)少于6個～8個，時間步距Δt應(yīng)滿足穩(wěn)定條件即Δt≤Δx/cmax，cmax為單元的最大波速值［2］。

2)利用動力響應(yīng)分析中的有限差分法［8］，獲得內(nèi)節(jié)點(diǎn)運(yùn)動方程:

并對式(6)求解。其中，利用MATLAB數(shù)學(xué)軟件獲得單元的剛度矩陣K;單元質(zhì)量矩陣M由集中質(zhì)量模型通過MATLAB編程獲得;F(t)通過集中力的等效荷載列陣獲得。利用ABAQUS軟件對土體模型進(jìn)行模態(tài)分析，提取模型前2階固有頻率f1，f2，選取合適的阻尼比ε，進(jìn)而確定α，β參數(shù)，最終獲得瑞利阻尼矩陣C=αM+ βK。

3)在計算區(qū)域內(nèi)選取計算點(diǎn)和離散點(diǎn)，通過空間內(nèi)插方法，并且用計算點(diǎn)數(shù)據(jù)來表示邊界節(jié)點(diǎn)的運(yùn)動方程，MTF由該運(yùn)動方程表示。

4)基于ABAQUS軟件平臺的DISP［9］子程序接口，通過FORTRAN語言來編寫位移迭代的程序，編寫的程序代碼如下:

SUBROUTINE DISP(U，KSTEP，KINC，MODE，NOEL，JDOF，COORDS)

INCLUDE‘TOUSHE_BIANJIE.INC’

DIMENSION U，TIME，COORDS

定義位移邊界條件U(根據(jù)邊界逐個節(jié)點(diǎn)位移時程編輯迭代程序)

RETURN

END

程序運(yùn)行合格后，通過Job模塊General選項(xiàng)卡中的User subroutine file路徑，完成模型透射邊界條件的施加。

2 無限元邊界

2.1 無限元

早在1977年的時候，無限元的概念已經(jīng)被Zienkiewicz提出，無限元在幾何上是無限大的但是由有限個單元組成的［10］。無限元法在解決半無限和無限域等空間問題時經(jīng)常與有限元方法聯(lián)合使用，常被看作是對有限元法進(jìn)一步的補(bǔ)充，故有“半有限元”的稱號。本文在ABAQUS有限元軟件的荷載分析步模塊，建立隨時間變化的脈沖荷載，荷載作用點(diǎn)的位置按模型尺寸考慮。由于涉及動力分析問題，故采用動力分析中的無限元建立無限元模型并進(jìn)行分析(見圖4)，要求模擬單元采用四節(jié)點(diǎn)無限單元(見圖5)，并且確保節(jié)點(diǎn)應(yīng)按逆時針的規(guī)則進(jìn)行編號和單元的第一個面應(yīng)為有限元和無限元的交接面［11］。

圖3 透射邊界模型示意圖

圖4 無限元模型示意圖

圖5 四節(jié)點(diǎn)無限單元

2.2 無限元的原理

邊界面反射并引起分析網(wǎng)格能量疊加，這是波在傳播過程中通過人工截斷邊界時，所面臨的最大問題。在ABAQUS動力分析中，無限元對這種波的傳播問題有比較好的處理方法。通過對線性土體中一維波的傳導(dǎo)進(jìn)行分析，獲得波傳導(dǎo)的動力平衡方程為:

其中，ρ為密度;E為彈性模量;x為軸向坐標(biāo)。

則當(dāng)波傳導(dǎo)到邊界截面上時，波的形式為:

反射波的形式為:

為了達(dá)到不出現(xiàn)反射的情況，可以在邊界上設(shè)置一個阻尼邊界條件，使波動能量能通過邊界，向無窮遠(yuǎn)處傳播。

邊界上的應(yīng)力:

在邊界上設(shè)定一個阻尼邊界條件:

為了消除反射波的影響，假定不存在反射波，并聯(lián)立式(9)和式(10)則可得到:

也就是說，只要邊界阻尼參數(shù)選擇合適，就可模擬無反射的情況。ABAQUS軟件在動力分析無限元中考慮了阻尼的設(shè)置，從而可以保證波在通過截斷邊界時，波動能量被無限域吸收，從而對分析區(qū)域的影響可控制在允許范圍之內(nèi)。

3 算例分析

3.1 參數(shù)設(shè)置

本文將成層土體分兩層，上層厚度5 m，下層厚度15 m。通過ABAQUS通用有限元軟件建立遠(yuǎn)置邊界模型，同時將計算結(jié)果作為精確解用于對比，模型尺寸取1 000 m×200 m，單元尺寸取Δx=Δy=1 m;其他邊界模型尺寸取60 m×20 m，單元尺寸取Δx=Δy=1 m;其中，無限元邊界取 Δx=30 m，Δy=1 m。荷載為作用于模型頂面中心(即(0，20)或(0，200))沿負(fù)Y方向的脈沖荷載，最大值Fmax=1.5×104N，荷載幅值曲線如下(見圖6)。

圖6 脈沖荷載幅值曲線

其中，分析步采用ABAQUS/Explicit動力顯式分析步，計算分析步總長度step=6 s。表1列出所需要的模型參數(shù)。

表1 材料參數(shù)表

3.2 結(jié)果分析

在地面，取距震源5 m點(diǎn)處;在地下，取距震源深度為3 m，5 m和10 m點(diǎn)處，作為觀察點(diǎn)提取Y方向的位移。圖7～圖10給出了在不同觀察點(diǎn)處采用不同邊界條件的位移反應(yīng)時程曲線。

圖7 地表距震源5 m點(diǎn)計算結(jié)果對比圖

圖8 地下距震源3 m點(diǎn)計算結(jié)果對比圖

圖9 地下距震源5 m點(diǎn)計算結(jié)果對比圖

1)從圖7～圖10可以看出，垂向3 m，5 m和10 m點(diǎn)處的位移幅值，隨著與震源距離的增大，出現(xiàn)明顯的衰減趨勢;分析圖7和圖9，地面5 m點(diǎn)處振動幅值要小于地下5 m點(diǎn)的振動幅值，說明振動在傳播過程中，在不同方向與震源相同距離處，振動衰減程度是不同的。對比圖7和圖9中，波動能量在水平方向的衰減幅度要高于豎向，這與波在土中的傳播規(guī)律相吻合。

圖10 地下距震源10 m點(diǎn)計算結(jié)果對比圖

2)由圖7～圖10中的固定邊界模型計算數(shù)據(jù)曲線，可以發(fā)現(xiàn):在前1 s時間段內(nèi)，波動主要由入射波引起的，在1 s～2 s左右時間段內(nèi)，波動主要由反射波引起的。說明:振動在傳播到固定邊界時，波動能量不能透過邊界，而在半無限土體中形成反射現(xiàn)象，這與設(shè)置邊界條件的初衷不符。

3)采用無限元邊界時，由于模型設(shè)置了阻尼矩陣，相當(dāng)于添加了吸收邊界條件。雖然反射現(xiàn)象依然存在，但相對固定邊界，反射幅值已經(jīng)明顯減弱。

4)通過透射邊界模型結(jié)果與遠(yuǎn)置邊界的模型計算結(jié)果(精確解)進(jìn)行比較，可以發(fā)現(xiàn):透射邊界結(jié)果與精確解更接近，能夠較好的模擬振動波在半無限土體中的傳播。說明，在成層土體中施加透射邊界的方法，可以獲得與精確解總體趨勢符合程度較好的效果。因此，本文施加透射邊界的方法是可行并且有效的。

5)通過對地下5 m點(diǎn)與3 m，10 m點(diǎn)位移時程比較，說明本文基于ABAQUS軟件在不同介質(zhì)分界面處施加透射位移邊界條件也是同樣適用的。

4 結(jié)語

人工邊界的選取將很大程度上決定數(shù)值模型的精度?；诖?，本文在對成層土體建立分析模型時，考慮透射邊界理論和ABAQUS軟件平臺及子程序接口，提出了在成層土體中施加透射邊界方法。該方法通過算例分析，表明是合理的;通過與無限元邊界、固定邊界、遠(yuǎn)置邊界的模型結(jié)果比較，表明精度是可靠的;并且具有較高的實(shí)用價值，通過編程可以擴(kuò)展復(fù)雜的二維問題和三維問題。

參考文獻(xiàn):

［1］莊 茁.ABAQUS非線性有限元分析與實(shí)例［M］.北京:科學(xué)出版社，2005:1-21.

［2］廖振鵬.工程波動導(dǎo)論［M］.第2版.北京:科學(xué)出版社，2002:232-285.

［3］廖振鵬，周正華，張艷紅.波動數(shù)值模擬中透射邊界的穩(wěn)定實(shí)現(xiàn)［J］.地球物理學(xué)報，2002，45(4):533-545.

［4］李小軍，廖振鵬.非線性結(jié)構(gòu)動力方程求解的顯示差分格式的特性分析［J］.工程力學(xué)，1993，10(3):143-148.

［5］李小軍，楊 宇.透射邊界穩(wěn)定性控制措施探討［J］.巖土工程學(xué)報，2012，34(4):641-645.

［6］周正華，廖振鵬.消除多次透射公式飄移失穩(wěn)的措施［J］.力學(xué)學(xué)報，2001，33(4):550-554.

［7］楊 宇.場地地震波動中透射邊界穩(wěn)定性問題研究［D］.北京:中國地震局地球物理研究所，2011.

［8］周昌玉，賀小華.有限元分析的基本方法及工程應(yīng)用［M］.北京:化學(xué)工業(yè)出版社，2006:46-76.

［9］李 帆，巴振寧，肖成志，等.波動數(shù)值模擬中人工透射邊界的實(shí)現(xiàn)技術(shù)［J］.自然災(zāi)害學(xué)報，2010，19(4):49-53.

［10］Zienkiewicz O C.Diffraction and refraction of surface waves using finite and infinite elements［J］.International Journal for Numerical Methods in Engineering，1977(11):1270-1293.

［11］費(fèi) 康，張建偉.ABAQUS在巖土工程中的應(yīng)用［M］.北京:中國水利水電出版社，2010:44-88.