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

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

1 引 言

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

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

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

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

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

(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 遠(yuǎn)域邊界的處理方法

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

(4)

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

圖2 自適應(yīng)半無限邊界單元

4 數(shù)值處理

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

4.1 數(shù)值離散

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

(5)

式中ξ為自然坐標(biāo)，取值范圍為[-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表達(dá)式為

(8)

式中a,b=1或2。

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

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

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

(9)

(10)

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

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

(11)

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

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

4.3 組裝

(12)

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

5 驗 證

5.1 算例描述

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

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

5.2 對比驗證方案

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

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

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

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

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

5.3 計算結(jié)果及討論

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

圖5 近場加載r =10 m處計算結(jié)果

圖6 近場加載r =20 m處計算結(jié)果

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

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

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

6 結(jié) 論

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

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