在水平成層土體中波動數值模擬的人工邊界★

崔高航 王兆亮 劉守花

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

0 引言

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

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

1 透射邊界

1.1 透射邊界定義

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

圖1 多次透射（MTF）公式的幾何關系

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

1.2 透射邊界的穩定性

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

參考相關文獻采取處理措施如下:

1)消除高頻振蕩失穩的措施。

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

2)消除零頻飄移失穩的措施。

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

該算子一般取γ≤0.05的正數，但不宜過小，主要結合經驗和試算確定。從實用的觀點來看，MTF具有普適性、可控性、時空解耦性，在保證模擬精度的情況下，能夠在計算機上穩定實現。

1.3 透射邊界的模型實現

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

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

圖2 二維土體模型人工邊界計算區

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

經過插值后的二階透射公式變為:

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

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

2)利用動力響應分析中的有限差分法［8］，獲得內節點運動方程:

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

3)在計算區域內選取計算點和離散點，通過空間內插方法，并且用計算點數據來表示邊界節點的運動方程，MTF由該運動方程表示。

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

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

INCLUDE‘TOUSHE_BIANJIE.INC’

DIMENSION U，TIME，COORDS

定義位移邊界條件U(根據邊界逐個節點位移時程編輯迭代程序)

RETURN

END

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

2 無限元邊界

2.1 無限元

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

圖3 透射邊界模型示意圖

圖4 無限元模型示意圖

圖5 四節點無限單元

2.2 無限元的原理

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

其中，ρ為密度;E為彈性模量;x為軸向坐標。

則當波傳導到邊界截面上時，波的形式為:

反射波的形式為:

為了達到不出現反射的情況，可以在邊界上設置一個阻尼邊界條件，使波動能量能通過邊界，向無窮遠處傳播。

邊界上的應力:

在邊界上設定一個阻尼邊界條件:

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

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

3 算例分析

3.1 參數設置

本文將成層土體分兩層，上層厚度5 m，下層厚度15 m。通過ABAQUS通用有限元軟件建立遠置邊界模型，同時將計算結果作為精確解用于對比，模型尺寸取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))沿負Y方向的脈沖荷載，最大值Fmax=1.5×104N，荷載幅值曲線如下(見圖6)。

圖6 脈沖荷載幅值曲線

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

表1 材料參數表

3.2 結果分析

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

圖7 地表距震源5 m點計算結果對比圖

圖8 地下距震源3 m點計算結果對比圖

圖9 地下距震源5 m點計算結果對比圖

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

圖10 地下距震源10 m點計算結果對比圖

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

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

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

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

4 結語

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

參考文獻:

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

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

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

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

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

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

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

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

［9］李 帆，巴振寧，肖成志，等.波動數值模擬中人工透射邊界的實現技術［J］.自然災害學報，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］費 康，張建偉.ABAQUS在巖土工程中的應用［M］.北京:中國水利水電出版社，2010:44-88.