遠場水下爆炸沖擊作用下重力壩動力穩定性的簡化時程分析方法

黃宇峰，盧文波，王高輝，陳 明，嚴 鵬

(1.武漢大學水資源與水電工程科學國家重點實驗室，湖北，武漢 430072；2.武漢大學水工巖石力學教育部重點實驗室，湖北，武漢 430072)

高壩大庫由于其顯著的政治經濟效益，極易成為局部戰爭[1-2]以及恐怖襲擊[3-4]的首要目標，高壩大庫一旦失事，上游庫區大量洪水下泄，將會給下游的人民帶來不可估量的災難。重力壩作為一種被廣泛采用的壩型，其抗爆性能以及安全評價越來越受到重視。

數值模擬和模型試驗方法在結構抗爆的研究中得到了廣泛的應用[5-10]，國內外學者采用這兩種方法，對近場和接觸水下爆炸作用下重力壩的動態響應[11-13]、毀傷特性[14-15]以及破壞模式[16-19]進行了大量的研究。胡晶等[12]進行了擋水結構水下爆炸作用的離心模型試驗，分別分析了水下爆炸沖擊波及氣泡脈動對結構的作用機理；王高輝等[15]對水下爆炸沖擊作用下重力壩的毀傷機理進行了研究，并據此得到了混凝土重力壩的毀傷模式；LI 等[19]研究了水下接觸爆炸下重力壩的破壞模式，探討了靜水壓力以及起爆位置對重力壩的損傷特性的影響。水下近場以及接觸爆炸工況下沖擊波峰值壓力大，但衰減得很快，因此重力壩受到的破壞主要為局部毀傷破壞，主要表現為迎爆面及附近的爆炸成坑破壞及沖切破壞、下游壩面的拉伸破壞等，不會對重力壩的整體穩定性造成較大的威脅；而遠場起爆工況下，雖然壩面受到的峰值壓力較小，但由于其衰減較慢，作用時間更長，影響范圍更廣，對重力壩壩的整體穩定性造成的威脅更大。因此，需重點關注重力壩在遠場水下爆炸沖擊下的穩定安全評價問題。

近年來有關爆炸沖擊等荷載作用下的結構動力穩定性能的研究取得了較大的進展[20-22]，在重力壩的穩定計算方面，目前的研究主要集中于復雜地質條件下的深層抗滑穩定計算以及地震作用下重力壩的動力穩定性計算，其主要計算方法有剛體極限平衡法以及有限元法。剛體極限平衡法[23]將大壩視為剛體，以滑移體上的抗滑力與滑動力之間的關系來評價重力壩的穩定性，其力學概念明確，計算方法簡便，是重力壩穩定分析的主要方法。近年來有學者提出了一種動力有限元時程法，以有限元方法計算區域內單元的應力，積分得到滑移面上的滑動力與抗滑力，再結合極限平衡理論，計算得到穩定安全系數，該方法在重力壩抗震動力穩定分析中應用比較廣泛[24-25]。

綜上，目前尚缺乏對遠場水下爆炸沖擊作用下重力壩穩定性的研究。本文正是聚焦于此問題，綜合考慮了遠場水下爆炸沖擊波的傳播及其與重力壩的相互作用過程，基于波動理論以及平面波假設，建立了遠場水下爆炸沖擊荷載作用下重力壩動應力場的簡化理論計算模型，并結合時程分析方法，建議了一種遠場水下爆炸沖擊作用下重力壩沿建基面的穩定安全系數的計算方法，為重力壩的抗爆安全評價提供理論參考。

1 遠場水下爆炸沖擊波的傳播及其與重力壩的相互作用過程

本文研究中的遠場起爆工況為起爆距離大于50 倍的藥包半徑的水下起爆工況。炸藥在上游庫區起爆后，會以球面波的形式向四周傳播。當沖擊波傳播至上游壩面后，會與壩面相互作用，一部分沖擊波會透射進入壩體內部，以應力波的形式繼續向下游壩面傳播； 當應力波傳播至下游壩面時，會在下游壩面發生反射，形成反射縱波與反射橫波；反射波繼續向壩上游側傳播，當其傳播至上游壩面或是建基面等結構面時，又會發生透反射，如此循環往復，直至波的能量完全衰減。遠場水下爆炸的沖擊波的傳播及其與壩體的相互作用過程十分復雜，為了便于研究，本文對其進行了簡化：假設沖擊波傳播至上游壩面時，以平面波的形式作用于壩體[26]；由于波在傳播過程中不斷的衰減，特別是在壩體混凝土與其他介質分界面處發生透反射后，波的能量將大幅削弱，因此在研究時只考慮沖擊波在壩面的透反射以及隨后應力波在下游壩面的反射對壩體穩定性的影響，忽略反射波在壩基截面等其他分界面的透反射以及之后壩內波的傳播等因素。

簡化之后的遠場水下爆炸沖擊波的傳播及其與重力壩的相互作用過程如圖1 所示。

圖1 遠場水下爆炸沖擊波的傳播及其與重力壩的相互作用過程示意圖Fig.1 Schematic diagram of far-field underwater explosion shock wave propagation and its interaction with gravity dam

1.1 水下爆炸沖擊波與壩面的相互作用

自由場水下爆炸沖擊波峰值壓力與時程曲線一般采用Cole 經驗公式計算[27]：

傳播至上游壩面的沖擊波會在壩面發生透反射，反射壓縮波與入射壓縮波疊加，導致壩面壓力陡增，如圖2 所示。根據平面彈性波在不同介質間的透反射公式[28]，大壩表面實際受到的壓力為：

圖2 爆炸沖擊波與壩體的相互作用Fig.2 Interaction between explosion shock wave and gravity dam

式中：pm、分別為大壩表面實際峰值壓力以及入射峰值壓力；p(t)為大壩表面實際壓力時程曲線；ρ1、ρ2分別為水和壩體混凝土的密度；CP1、CP2分別為波在水中和混凝土中的波速。波在水中的波速可取1500 m/s，在混凝土中的波速CP根據式(6)計算得到，E、ρ、μ分別為壩體混凝土的彈性模量、密度和泊松比。

1.2 爆炸應力波在下游壩面的反射

透射進入壩體內部的應力波為縱波，當其傳播至下游壩面處時，透射波作為新的入射縱波IP，在下游壩面處發生反射，形成反射縱波RP和反射橫波RS，如圖2 所示。

由于波傳播過程中能量的耗散，傳播至下游壩面的應力波可視為彈性平面波。為了得到應力波在下游壩面處的反射系數，基于彈性平面波的假設，本文對入射縱波IP、反射縱波RP和反射橫波RS在下游壩面處引起的單元體應力狀態進行了分析，如圖3 所示。圖中σx1、σx2、σx3分別為入射縱波、反射縱波以及反射橫波引起的正應力，τxy1、τxy2、τxy3分別為入射縱波、反射縱波以及反射橫波引起的剪應力。

圖3 不同應力波作用下交界面處單元體應力狀態Fig.3 Stress state of unit body at the interface under the action of different stress waves

遠場起爆條件下爆炸荷載周期較長，一般可達十幾甚至幾十毫秒，其波長與重力壩的截面寬度較為接近，其引起的結構微元的慣性效應較小，為了簡化計算，本文在分析不同應力波作用下交界面處單元體的應力狀態時，忽略了慣性效應的影響。

入射縱波所引起的三角形單元應力狀態如圖3(a)所示，根據靜力平衡條件[29]可得：

式中：α 為入射角；AC=ABcosα，BC=ABsinα。

將AC、BC代入式(7)，解得：

對單元體ADE，根據平衡方程可得反射縱波引起的正應力與剪應力分別為：

對單元體AFG，根據平衡方程可得反射橫波引起的正應力與剪應力分別為：

式中：α′為縱波反射角；β 為橫波反射角，根據波傳播的Snell 定律，入射角與反射角之間滿足以下關系：

式中，CS為反射橫波波速，其計算式為：

不考慮下游庫水，將重力壩下游壩面視為自由面，在自由面上有 σx=0、τxy=0，即：

將式(8)～式(10)代入式(13)，即可得到入射縱波在下游壩面的縱波反射系數k1以及橫波的反射系數k2：

式中，A1、B1、C1、D1各系數的表達式為：

2 壩內動應力計算

應力波在壩體內傳播的過程中會發生能量的衰減，其應力峰值會隨著傳播距離的增大不斷減小，由于巖石和混凝土的物理性質較為接近，可參考爆炸應力波在巖石介質中峰值壓力的衰減公式：

式中：R為起爆距離，即爆源到重力壩上游壩面的距離；pm為壩面壓力峰值；x為應力波在壩內傳播的距離；γ 為壩體混凝土中峰值壓力的衰減系數，在彈性區可取1[30]。

根據前文的假設可知，對于壩體內任意一點M，該處的動應力應為從上游壩面透射進入壩體內部的P 波以及下游壩面反射形成的P 波和SV 波三者的疊加，如圖4 所示。

圖4 壩體動應力計算示意圖Fig.4 Schematic diagram of dam dynamic stress calculation

入射縱波IP在M點引起的動應力為：

式中：LOM為入射縱波傳播路徑OM的長度；t1為入射縱波從壩面傳播至M點所需要的時間，計算式為：

從O1 處入射的縱波，在下游壩面J 處發生反射后，形成的反射縱波RP，經JM 傳播至M 處，在反射縱波的作用下，M 點的動應力為：

從O2處入射的縱波，在下游壩面K處發生反射后，形成的反射橫波RS，經KM傳播至M處，在反射橫波的作用下，M點的動應力為：

在M點的位置給定以后，根據幾何關系，便可確定入射縱波、反射縱波以及反射橫波的傳播路徑，其表達式為：

式中：θ 為混凝土重力壩下游壩坡坡角；ε 為反射縱波與水平面之間的夾角，；φ為反射橫波與水平面之間的夾角，

在分別分析得到入射P 波、反射P 波以及反射SV 波在M點引起的動應力后，將其分別向x方向以及y方向分解，得到爆炸應力波作用下壩體質點的正應力σy與剪應力τyx。

對于入射縱波IP，其表達式為：

對于反射縱波RP，其表達式為：

對于反射橫波RS，其表達式為：

將三組波所引起的應力分別疊加，即可得到應力波作用下M點的正應力與剪應力的表達式：

3 重力壩遠場爆炸動力穩定性簡化時程分析方法

根據第1 節和第2 節推導得出的遠場水下爆炸沖擊作用下重力壩動應力場計算公式，即可求解重力壩建基面上的動應力場，結合時程法，就能得到重力壩在任意時刻沿建基面的抗滑穩定安全系數。

混凝土重力壩遠場爆炸動力穩定性簡化時程分析方法的具體步驟如圖5(a)～圖5(d)所示。

圖5 重力壩遠場爆炸動力穩定性簡化時程法計算流程圖Fig.5 Calculation flow chart of simplified time history method for far-field explosion dynamic stability of gravity dam

1) 首先根據給定的起爆藥量、起爆距離等參數，計算得到水下爆炸沖擊荷載，考慮波在壩面的透反射作用，計算得到壩面沖擊荷載時程曲線。

2) 計算應力波在下游壩面處的反射系數，得到透射波σIP、反射縱波σRP以及反射橫波τRS的應力時程表達式。

3) 分別計算入射縱波、反射縱波和反射橫波在建基面上引起的動應力，然后將三者引起的動應力疊加，得到建基面上任意質點的正應力σy與剪應力τyx。要注意的是，上述的正應力與剪應力是由水下爆炸沖擊荷載引起的，除此之外，建基面上的總應力還應該考慮靜力荷載作用下(重力G、靜水壓力荷載pw以及揚壓力荷載pu)產生的靜應力和。建基面上的靜應力可根據材料力學法確定。這樣，就能得到建基面上總的應力分布。

4) 將建基面分為n個條塊，每個條塊的長度記為li，取每個條塊中點處的應力代表整個條塊的應力，于是對于每個條塊，其滑動力為抗滑力為，再將每一份塊體的抗滑力與滑動力累加，即可得到沿建基面整體的滑動力與抗滑力，兩者的比值即為重力壩沿壩基面的抗滑穩定安全系數：

式中，f′、c′分別為壩體混凝土與壩基接觸面的抗剪斷摩擦系數及粘聚力。

4 混凝土重力壩動力穩定性分析

4.1 計算模型與參數

本文選取某混凝土重力壩典型壩段作為動力穩定性分析對象，該重力壩高度為80 m，壩基底部寬度為66 m，下游壩面折坡比為1∶0.8，上游正常蓄水位深度為75 m，壩段寬度為10 m。為驗證本文所建議方法的可靠性，建立了混凝土重力壩數值計算模型，如圖6 所示。本文采用聲固耦合方法，聲單元模擬水下爆炸沖擊波在庫水中的傳播，實體單元模擬大壩、地基。壩體混凝土、基巖均假定為線彈性材料，彈性模量為30 GPa，泊松比為0.167，密度為2400 kg/m3，單位類型為C3D8R；庫水選用聲學單元AC3D8R，其體積模量為2.14 GPa，密度為1000 kg/m3。無反射邊界條件是通過acoustic impedance(聲阻抗)的方式施加于水網格的邊界上，因為水體邊界是平面的，所以施加平面無反射邊界。本文荷載施加采用輸入壓力時程曲線的方式，入射波荷載以平面波的形式施加于庫壩耦合面，入射波壓力時程曲線根據式(1)～式(3)計算所得，如圖7 所示。

圖6 混凝土重力壩數值計算模型Fig.6 Numerical calculation model of concrete gravity dam

圖7 上游壩面入射壓力時程曲線Fig.7 Time history curve of incident wave pressure at the upstream dam surface

本文在進行穩定計算時，考慮的荷載組合為：上游壩面靜水壓力、揚壓力、重力壩自重以及水下爆炸沖擊荷載；壩體混凝土與壩基接觸面的抗剪斷摩擦系數f′取0.9、粘聚力c′取0.8 MPa。本文選取的計算工況為起爆藥量W為5000 kgTNT，起爆距離R為100 m。

4.2 方法可靠性驗證

為了驗證本文所建議的簡化時程分析方法的可靠性，基于ABAQUS 數值模擬軟件對遠場水下爆炸作用下重力壩內部的應力波傳播過程進行了模擬，得到了重力壩建基面質點動應力以及沿建基面的穩定安全系數的數值解，并與彈性平面波假設下的解析解進行了對比。由于數值計算中壩體以及基巖材料選取的是線彈性模型，在數值計算中無法模擬出應力波在傳播過程的能量損耗帶來的應力波峰值的衰減，因此選取衰減系數γ 為0 時的解析解與數值解進行對比，如圖8 所示。

圖8 解析解與數值解的對比Fig.8 Comparison between analytical solution and numerical solution

根據圖8(a)、圖8(b)可以發現，解析解與數值解在前期偏差較小，在后期數值解的波形較亂，與解析解之間的偏差較大。根據圖8(c)可以發現，穩定安全系數的數值解與解析解在30 ms之前吻合得較好，兩種方法計算得到的穩定安全系數的變化趨勢一致，在數值上的偏差也較小，前20 ms 兩者的偏差在20%以內，20 ms～30 ms 的偏差在30%以內。在30 ms 后，數值解出現了與解析解不一樣的變化規律，這是因為此時上游壩面反射的應力波傳播至壩基截面，在數值計算中，由于沒有考慮波在傳播過程中的衰減，壩基截面上產生了較大的剪應力，所以導致穩定安全系數顯著下降。但是在實際情況中，波在傳播過程中存在著能量的耗散，上游壩面反射形成的應力波在傳播至壩基截面后壓力很小，對重力壩動力穩定性的影響較小，因此在解析解中忽略了上游壩面反射波的影響。這就導致了在后期數值解與解析解之間的差距較大?？傮w看來，解析解與數值解在變化趨勢上具有良好的一致性，在數值上大部分時間兩者的偏差較小，從側面說明了該解析方法具有較好的可靠性。

4.3 大壩動力穩定性分析

圖9 給出了衰減系數γ 取1 時重力壩建基面形心處的動應力時程曲線。從圖9 可以看出，壩基截面形心處質點的動應力時程曲線可以明顯分為以下3 個階段：Ⅰ：透射縱波作用下的受壓階段；Ⅱ：反射縱波作用下的拉剪階段；Ⅲ：反射橫波作用下的壓剪階段。透射進入壩體內部的縱波為壓縮波，在其作用下，質點將處于受壓階段；當應力波傳播至大壩下游面后，會發生反射，形成一組縱波和一組橫波，縱波為拉伸波、橫波為剪切波。反射縱波傳播速度快，因此會先于反射橫波傳播至建基面形心處，在其作用下，該處質點會產生拉應力以及剪應力，該階段質點處于拉剪階段；之后，反射橫波也傳播至該處，導致該處質點產生較大的剪應力與壓應力，該階段質點處于壓剪階段。

圖9 建基面形心處質點動應力時程曲線 (γ=1)Fig.9 Time history curve of particle dynamic stress at the centroid of foundation surface (γ=1)

圖10 給出了衰減系數取1 時建基面上抗滑力、滑動力以及穩定安全系數隨時間變化的規律。根據圖10 可知，重力壩的穩定安全系數變化呈現出先增大后減小再增大并最終保持恒定的變化趨勢。初始時刻重力壩只受到透射壓縮波的作用，透射壓縮波在壩基截面引起壓應力，使得建基面上的抗滑力增大，從而導致初始時刻重力壩的穩定安全系數會略大于靜力穩定安全系數。但反射形成的拉伸波傳播至建基面后，會導致建基面產生拉應力與剪應力，而當剪切波傳至建基面后，又將會引起較大的剪應力以及壓應力，在兩組反射波的共同作用下，建基面上的滑動力顯著增加，同時抗滑力也顯著下降，重力壩的穩定安全系數不斷減小，在22 ms 時達到最小值。之后，隨著應力波的衰減，大壩的穩定安全系數不斷增大，逐漸趨于靜力穩定安全系數。

圖10 穩定安全系數時程曲線圖(γ=1)Fig.10 Time history curve of stability safety factor (γ=1)

5 結論與展望

本文基于平面波假設和波動理論，對遠場水下爆炸沖擊作用下重力壩的穩定安全問題進行了研究，得到以下結論：

(1) 基于遠場水下爆炸沖擊波的傳播及其與壩體的相互作用過程的分析，建立了重力壩爆炸動應力場的簡化理論計算模型，最終建議了一種遠場水下爆炸沖擊作用下重力壩動力穩定性的簡易時程分析方法，并得到了動力有限元數值方法的驗證。該方法計算簡便，能較為快速地計算得到重力壩沿建基面的動力穩定安全系數時程曲線，具有一定的可靠性。

(2) 在遠場水下爆炸沖擊作用下，重力壩沿建基面的抗滑穩定安全系數呈現出先增大后減小再增大并最終保持恒定的變化趨勢。反射橫波與反射縱波的共同作用會導致建基面的滑動力增大，抗滑力下降，使得重力壩的穩定安全系數顯著減小。

本文僅僅只是對水下爆炸打擊下重力壩動力穩定性問題進行了一個初步的研究，在實際的遠場水下爆炸打擊下，重力壩壩基會在沖擊荷載的作用下開裂，基巖的力學性質也會在應力波擾動下發生變化，如何在穩定分析中將這些因素加以考慮，仍有待進一步的研究。