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

崔高航 王兆亮 劉守花

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

0 引言

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

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

1 透射邊界

1.1 透射邊界定義

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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é)點位移時程編輯迭代程序)

RETURN

END

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

2 無限元邊界

2.1 無限元

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

圖3 透射邊界模型示意圖

圖4 無限元模型示意圖

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

2.2 無限元的原理

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

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

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

反射波的形式為:

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

邊界上的應力:

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

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

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

3 算例分析

3.1 參數(shù)設置

本文將成層土體分兩層，上層厚度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列出所需要的模型參數(shù)。

表1 材料參數(shù)表

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

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

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

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

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

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

4 結語

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

參考文獻:

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

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

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

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

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

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

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

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

［9］李 帆，巴振寧，肖成志，等.波動數(shù)值模擬中人工透射邊界的實現(xiàn)技術［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.