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

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

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

1 引 言

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

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

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

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

在分析半平面問題時，可將其邊界劃分為有限域邊界Γ1、無窮遠邊界Γ∞及一端延伸到無窮遠處的半無限邊界Γ2和Γ3，如圖1所示。由于Γ∞上的場點都位于無限遠處，有限時間段內任意瞬時發出的應力波波前只能到達有限遠處，對Γ∞上的場點無任何影響，因此在計算時無需考慮Γ∞。在初始條件為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 遠域邊界的處理方法

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

(4)

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

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

4 數值處理

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

4.1 數值離散

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

(5)

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

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

(6,7)

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

(8)

式中a,b=1或2。

4.2 影響系數的求解

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

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

(9)

(10)

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

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

(11)

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

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

4.3 組裝

(12)

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

5 驗 證

5.1 算例描述

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

圖4 動力荷載函數圖像

5.2 對比驗證方案

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

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

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

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

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

5.3 計算結果及討論

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

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

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

圖7 遠場加載節點N計算結果

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

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

6 結 論

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

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