SH 波入射凹陷半球形谷地的數值分析

瞿宗新 蔣 格

(南京航空航天大學土木工程系，江蘇南京 210016)

給定輸入波動，研究散射體及其附近介質中的動力反應稱為散射問題。對這類問題的分析可以采用解析法和數值法求解。解析法主要有波函數展開法，數值法主要有邊界積分法、有限元法和有限差分法等［1］。

局部地形對彈性波散射問題的精確解極為有限，其中有半圓形、半橢圓凹陷地形對平面SH波散射的級數解和半圓形凸起地形對平面SH波散射的解答。原則上講，數值法可應用于各種不規則地形［2］。

波動的時域數值模擬具有廣闊的應用前景，完善這一數值模擬技術的主要問題是在保證結果可靠的基礎上提高計算效率，而就算法而言，提高計算效率的重要途徑是發展顯示方法。采用集中質量方法，并通過時域顯示積分構成的一種顯示集中質量有限元技術，這一技術能夠大大的提高計算效率［3］。

本論文考慮地震波斜入射情形，通過理論分析、數值模擬，來實現SH波入射三維凹陷半球形谷地的數值分析方法。本文最后給出了分析算例和數值結果，并對其進行了討論。

1 基本原理

1.1 控制方程

地震波入射到土體介質的一個速度異常區，除了有入射的平面波外，還會產生從異常區各個方向傳播的散射波，包括反射、折射波等，這個異常區也稱之為散射體。在xy平面內，SH波所表示的位移u(x，y，t)與xy平面垂直。在一般常態下，位移函數u(x，y，t)滿足運動方程:

應力和應變方程為:

1.2 質點運動方程

將半空間無限均勻地基介質內的一有限土體通過有限元軟件離散成八結點實體單元。并劃分為兩個區域，即內部計算區和人工邊界區(如圖1所示)，單元結點分為內結點、人工邊界點、與基礎相聯的點這三類［5］。根據結構動力學原理，采用集中質量有限元法，內結點的振動方程如下:

其中，uI為內結點I的位移向量;uI為速度向量;uI為加速度向量;N為與I結點相鄰的其他結點總數;MI為內結點I的質量矩陣;KIL，CIL分別為結點I與相鄰結點L之間的剛度矩陣和阻尼矩陣;FI為結點I所承受的外荷載向量。

圖1 模型示意圖

由顯式差分格式:

則式(4)變成:

1.3 SH波散射數值模擬的基本步驟

1)前處理過程:進行土體模型的建模和網格劃分，將單元和結點坐標的數據導出，并提取出人工邊界結點;

2)自由場位移的計算:利用傳遞矩陣法，計算出人工邊界區的自由場響應，將其作為計算土層動力反應時該邊界的輸入;

圖2 輸入脈沖的波形圖

圖3 輸入脈沖的幅值頻譜圖

3)質點總位移的計算:將自由場響應作為輸入，并利用集中質量有限元方法，通過顯示數值積分方法來計算區域內結點的總位移場。對于側邊界區的質點單元，通過總位移場減去自由場得到質點的散射場，根據透射公式可以得到邊界區下一時步的散射場，再加上自由場位移可得到邊界結點下一時步的總位移［6］。而底邊界可以依據線性波動傳播直接給出入射場，進而得到底邊界區的散射場，再利用透射公式進行求解。

2 算例

為了驗證此數值方法的有效性，我們取一個模型，模型形式為三維凹陷半球體的長方體均勻介質土體模型，如圖1所示。并將所得的結果與解析解進行對比分析驗證。

2.1 模型參數

圖4 圖形結構示意圖

選取人工計算區域為80 m×80 m×40 m的長方體均勻介質，以上表面中心為球形的半球體部分為凹陷土體，半徑為10 m。將距邊界10 m范圍以內的介質取為人工邊界區。介質的剪切波速和質量密度分別為 cvs=400 m/s，ρv=3 000 kg/m3，彈性模量E=1.248 ×109Pa，泊松比 v=0.25，阻尼比 ξ=0.02。

通過ANSYS軟件進行建模，模型采用八結點單元進行網格劃分。在模型左下角處輸入脈沖位移，脈沖寬度為0.2 s，時間步距為0.000 5 s，計算時步為8 192步，有限單元的大小取為1 m×1 m×1 m，有限單元模型的總結點數為267 311個，模型單元總數為254 040個(見圖1)。

本算例主要給出圖1中A，F兩點所對應位置的位移響應時程圖和頻譜圖。E為半球形球心，坐標為(40，40，0)。其余各點的坐標依次為 A(0，40，0)，F(50，40，0)。

輸入脈沖位移的波形圖和輸入脈沖位移的幅值譜圖分別見圖2，圖3。

2.2 圖形結果

圖形結果見圖4。

2.3 結果分析

由波動理論可知，當SH波入射自由表面時，不論入射角和入射波頻率的差異，自由表面的響應(位移，速度，加速度)總是入射波的兩倍，也就是自由表面的動力放大因子為2［7，8］。由圖4可知，波峰的峰值都接近于2，與理論解相符合。

通過比較30°和60°下A，F兩點的位移時程曲線和不同頻率下表面各點的放大系數圖，可以得出如下結論:

1)SH波以不同角度入射時，這點的位移反映波形和達到波峰的時刻均不同。SH波以一定角度入射時，不同點的位移波形也不同，波形起始時刻與質點位置有關。

2)位于散射體附近點的放大系數明顯高于周圍的一些點。隨著無量綱頻率的增加，譜放大系數的變化越來越復雜。

3 結語

本文實現了一種對SH波散射三維凹陷半球形谷地問題的數值解法，并通過計算結果驗證了方法的可行性。算例中給出了SH波分別以30°角和60°角斜入射三維凹陷半球形谷地時表面一些點的位移時程圖和放大系數圖。該數值方法具有簡便高效的優點，比傳統的解析方法更具有優越性。該數值計算方法在土體—結構動力相互作用分析領域具有廣闊的應用前景。

［1］梁建文，張郁山，顧曉魯.圓弧形層狀凹陷地形對平面SH波的散射［J］.振動工程學報，2003，16(2):158-165.

［2］袁曉銘，廖振鵬.圓弧形凹陷地形對平面SH波散射問題的級數解答［J］.地震工程與工程振動，1993，13(2):1-11.

［3］廖振鵬，劉 恒，謝志南.波動數值模擬的一種顯示——一維波動［J］.力學學報，2009，41(3):350-360.

［4］林 宏，劉殿魁.半無限空間中圓形孔洞周圍SH波的散射［J］.地震工程與工程振動，2002，22(2):9-16.

［5］陳少林，唐 敢，劉啟芳，等.三維土—結構動力相互作用的一種時域直接方法［J］.地震工程與工程振動，2010，30(2):24-31.

［6］黃祖光，劉 云.吸收邊界條件下平面SH波散射的數值模擬［J］.山西建筑，2005，31(2):25-26.

［7］梁建文，巴振寧.彈性層狀半空間中沉積谷地對入射平面SH波的放大作用［J］.地震工程與工程振動，2007，27(3):1-9.

［8］李山有，王學良，周正華.地震波斜入射情形下水平成層半空間自由場的時域計算［J］.吉林大學學報，2003，33(3):372-376.