半平面問題遠域輻射阻尼的時域邊界元法模擬研究

張 梅， 唐 波,2， 李宏軍*， 郝文秀， 龔海天

(1.河北農(nóng)業(yè)大學 城鄉(xiāng)建設學院，保定 071001； 2.上海強勁地基工程股份有限公司，上海 200233；3.新加坡清水建設株式會社，新加坡 339509)

1 引 言

當無限及半無限域受到動力荷載作用時，應力波向遠域介質傳播，能量發(fā)生逸散，形成輻射阻尼。輻射阻尼對受力體的動力響應具有重要影響，甚至較之受力體本身阻尼更為重要[1]，但半無限域問題中遠域輻射條件較難滿足，致使遠域輻射阻尼的模擬尚沒有較為完善的方案。實際工程中諸如巖石基坑爆破開挖、邊坡穩(wěn)定、結構抗震以及結構-地基動力相互作用[2]等問題均會涉及到遠域輻射阻尼。

目前處理遠域輻射阻尼的方法主要有有限元法與邊界元法。由于實際問題應力波的傳播距離是有限的，無論采用哪種數(shù)值方法，都可以通過截取足夠大區(qū)域(完全包含應力波的影響范圍)進行計算，從而避免輻射阻尼的處理，其操作方法簡單，但是離散單元多，計算量較大，對于不同時程需求的問題往往需要多次改變模型，計算繁復。對于有限元法，若截取區(qū)域較小(未完全包含應力波的影響范圍)，需要設置人工邊界[3]或直接采用比例邊界有限元法[4,5]，以考慮遠域輻射阻尼的影響。人工邊界的引入可以很好地模擬各種無限域的輻射阻尼，但也存在缺點,即低階邊界精度低和高階邊界穩(wěn)定性較差[6]；比例邊界元法兼具有限元法和邊界元法的優(yōu)點，是處理遠域輻射阻尼的一種較為理想的方法，但由于需要確定固定的比例中心，在求解具有偏心域或平行側邊界的問題時較為困難。對于無限域問題，采用全空間基本解的邊界元法，能夠自動滿足無窮遠處的輻射條件，目前國內(nèi)外對于該問題的研究已經(jīng)較為成熟[7-10]。對于半無限域問題，國內(nèi)外學者通過采用基于半空間基本解的邊界元法來滿足無限長邊界上的輻射條件，并進行了較為深入的研究[11-13]，然而動力學問題的半空間基本解較為復雜，且邊界不規(guī)則問題的基本解不易得到，甚至不存在[14,15]，使得該方法的通用性不強。在處理半無限域動力學問題時，有較多學者采用基于較為成熟的全空間基本解的時域邊界元法[15,16]，但由于半無限邊界的存在，通常只能截取部分有限邊界進行動力分析，不能完全考慮輻射阻尼影響，有些學者在邊界離散的基礎上采用無窮邊界元模擬半無限邊界[15,17]，對半空間進行頻域動力分析，既保證了結果的準確性，又提高了求解效率，但該方法不能求解非線性問題，且需要在頻域和時域間多次轉換，計算繁復。綜上，目前國內(nèi)外對于半無限域輻射阻尼的模擬尚沒有較為完善的方案。

本文以全空間基本解的TD -BEM[7-10]為基礎，結合應力波的傳播特性，在時域內(nèi)提出了一種新的單元,即自適應半無限邊界單元，專門用于離散遠域半無限邊界，以考慮遠域輻射阻尼，近域邊界及無限邊界仍采用常規(guī)TD -BEM處理。該單元實質上是一個可隨時間-空間參數(shù)自動調整單元大小的有限單元，因此，可以采用目前已經(jīng)成熟的數(shù)值方法[7,8]對該單元進行處理，理論上不具有困難。

2 半平面問題的時域邊界積分方程

在分析半平面問題時，可將其邊界劃分為有限域邊界Γ1、無窮遠邊界?！藜耙欢搜由斓綗o窮遠處的半無限邊界Γ2和Γ3，如圖1所示。由于?！奚系膱鳇c都位于無限遠處，有限時間段內(nèi)任意瞬時發(fā)出的應力波波前只能到達有限遠處，對Γ∞上的場點無任何影響，因此在計算時無需考慮?！蕖Ｔ诔跏紬l件為0和忽略體力的條件下，彈性動力學半無限平面問題的時域邊界積分方程為

(1)

(cs/cd)(Fi k/Ld+Ji kLdNd)Hd]

(2)

(3)

圖1 半無限平面邊界

式中Ei k=δi k,Fi k=δi k/r2,Ji k=-r,ir,k/r2

Hw=H[cw(t-τ)-r]

3 遠域邊界的處理方法

有限時間段內(nèi)的動力學問題，應力波傳播距離亦為有限遠，對于超過Γ1范圍的應力波，應考慮無限域的輻射阻尼，反之，始終處于Γ1范圍的應力波不需考慮此效應。實際工程中，涉及無限域的輻射阻尼問題較為常見，為了較好解決這一效應的模擬問題，本文將半無限邊界Γ2和Γ3分別采用一個自適應半無限邊界單元進行數(shù)值離散，半無限邊界與有限邊界公共點A和應力波的波前位置B分別作為自適應半無限邊界單元的內(nèi)外側節(jié)點，構建自適應半無限邊界單元，如圖2所示，P點為任意應力波源點，cd(t-τ)為τ時刻作用在源點P處的應力波在t時刻傳播的最遠距離。根據(jù)幾何關系，節(jié)點B的坐標通式為

(4)

本文提出一種新的單元，即自適應半無限邊界單元，專門用于離散遠域半無限邊界，該單元外側節(jié)點B是一個始終處于應力波波前位置的動態(tài)節(jié)點，保證計算區(qū)域在任何情況下都恰好包含應力波的影響范圍，從理論上為遠域輻射阻尼的模擬提供保障。

圖2 自適應半無限邊界單元

4 數(shù)值處理

首先將邊界進行數(shù)值離散，再求出各單元的位移和面力影響系數(shù)，最后進行組裝形成代數(shù)方程組并求解邊界點位移。由于邊界Γ1為有限邊界，其數(shù)值處理方法與常規(guī)TD -BEM的處理方法完全相同，不再贅述，本節(jié)著重介紹半無限邊界單元的數(shù)值處理。

4.1 數(shù)值離散

空間上，根據(jù)應力波在彈性介質中的傳播特性，將兩側的半無限邊界轉化為兩個自適應半無限邊界單元。假定面力和位移在半無限邊界單元上線性變化，其插值函數(shù)表達式為

(5)

式中ξ為自然坐標，取值范圍為[-1,1]。

時間上，將時間域[0,t]離散為M個步長為Δt=t/M的時間單元。且令時間節(jié)點tm=mΔt，其中m取值范圍為0,1,…,M。假定在每一時間步上面力為常量，位移線性變化，對于任意時間單元，面力和位移的插值函數(shù)分別為

(6,7)

離散后，單元面力影響系數(shù)gi k與位移影響系數(shù)hi k表達式為

(8)

式中a,b=1或2。

4.2 影響系數(shù)的求解

通過數(shù)值離散，將邊界積分方程的半無限邊界部分轉化為自適應半無限邊界單元積分，積分中奇異性情況及求解方法如下。

(1) 無奇異性。采用高斯數(shù)值積分法計算。

(9)

(10)

式中 表示Hadamard主值積分的計算符號。

(3) 空間奇異性。此時r→0，奇異性的表現(xiàn)形式為1/rn→∞。這類奇異性可通過考慮P波及S波的相互作用，將兩種波對應的數(shù)學表達式進行代數(shù)求和，采用分子有理化的數(shù)學方法對奇異積分解析求解，可表示為

(11)

(4) 雙重奇異性。即同時具有波前和空間奇異性。這類奇異性需要在波前及空間奇異性處理的基礎上，再通過在空間上影響系數(shù)組裝時，考慮變形協(xié)調條件聯(lián)系奇異點兩側單元影響，將兩側單元相應表達式進行代數(shù)求和，即可解析互消。

按以上方法求解，可以得到各單元的影響系數(shù)。

4.3 組裝

(12)

(e=2～Ne-1)(13)

5 驗 證

5.1 算例描述

圖3 動力荷載作用下的彈性半平面計算模型

圖4 動力荷載函數(shù)圖像

5.2 對比驗證方案

采用3種方案對兩種加載方式的算例進行求解，并將計算結果進行對比驗證。

(1) 采用FEM(ANSYS)計算，取距荷載60 m范圍內(nèi)的區(qū)域作為計算域，共離散了44984個PLANE42單元。

(2) 采用常規(guī)TD -BEM求解，取距原點25 m以內(nèi)區(qū)域進行計算，共離散500個有限邊界單元。

(3) 采用結合自適應半無限邊界單元的時域邊界元法求解(下稱自適應TD -BEM)，與方案(2)相比，僅在有限邊界兩側各添加了一個自適應半無限邊界單元，其他處理保持不變。

此處需要說明，離散系數(shù)β[9](β=cdΔt/lmax，lmax為最大單元長度)是影響TD -BEM計算精度的重要因素之一。對比半平面問題β=0.5～3.0 時的計算結果，發(fā)現(xiàn)β=1.5～2.5時，較為理想。因此，方案(2，3)選取了β=2對有限域邊界進行離散。

5.3 計算結果及討論

近場加載選取地表r=10 m和20 m處節(jié)點的位移響應進行比對，結果如圖5和圖6所示；遠場加載選取溝槽底部N點進行比對，結果如圖7所示(注：圖5～圖7設置為每6個時間步顯示一個結果點，并對橫坐標時間變量進行了無量綱化處理)。

圖5 近場加載r =10 m處計算結果

圖6 近場加載r =20 m處計算結果

圖7 遠場加載節(jié)點N計算結果

計算結果表明,常規(guī)TD -BEM由于僅截取了部分有限邊界進行計算，不能考慮應力波在遠域的輻射條件，應力波出現(xiàn)了回彈，近場入射r=10 m和20 m地表處以及遠場入射的節(jié)點N處，分別在cdt/r0=40，30和30時，回彈應力波剛好到達計算點，節(jié)點在受到回彈應力波的影響后，其位移與FEM結果相比出現(xiàn)較大差異；而自適應TD -BEM得到的結果能夠與FEM結果吻合良好，滿足了應力波在半無限域上的輻射條件，較好地模擬了遠域輻射阻尼效應。

如果想要通過常規(guī)TD -BEM得到更精確的結果，勢必要擴大計算范圍，離散更多的網(wǎng)格。而自適應TD -BEM僅在原有限邊界的基礎上添加兩個自適應半無限邊界單元，即可自動滿足半無限域的輻射條件，在計算時間成本與常規(guī)TD -BEM幾乎相同的情況下，得到了較高精度的計算結果。

6 結 論

本文根據(jù)應力波在彈性介質中的傳播特性，將半無限邊界轉化為自適應半無限邊界單元，該單元可根據(jù)具體的時間-空間參數(shù)自動調節(jié)單元大小，保證計算區(qū)域在任何情況下都恰好包含應力波的影響范圍，滿足了應力波在半無限域上的輻射條件。采用自適應半無限邊界單元的TD -BEM較好地解決了遠域輻射阻尼的模擬問題，且在計算時間成本與常規(guī)TD -BEM幾乎相同的前提下，具有更高的計算精度。

需要說明的是，本文以彈性動力學問題為切入點，重點研究遠域輻射阻尼的模擬方法。對于非線性動力學問題，足夠遠域介質往往處于彈性階段，因此，本文提出的基于TD -BEM的自適應半無限邊界單元模擬遠域輻射阻尼的方法和結論依然是適用的。對于近域非線性問題TD -BEM的詳細介紹可參見文獻[20,21]。