孔隙水壓力及炸藥埋深對堤壩爆炸效應的影響分析

張智超，陳育民，劉漢龍，王維國

（1.河海大學 巖土力學與堤壩工程教育部重點實驗室，南京 210098；2.河海大學 土木與交通學院，南京 210098）

1 引 言

我國水利水電工程眾多，許多重要的水工建筑位于人口密集的地區，三峽水利樞紐工程的攔河大壩高達185 m，最高蓄水水位可達175 m，南水北調工程有著近3000 km長的引水渠道，這些都不可避免地涉及到堤壩的建設和運營。加之2011年“中央一號文件”關于加快水利改革發展的決定，更是將堤壩的興修和維護推向了新的高潮。然而，正是由于水利事業的發展與繁榮，其在給人民群眾帶來巨大的經濟效益的同時，也往往成為戰爭或者恐怖分子的重點打擊對象。一旦此類堤壩工程在爆炸等襲擊下發生破壞乃至潰壩，將會給人民的生命財產造成巨大的損失。從另一方面來說，人為地利用爆炸荷載對堰塞壩進行破壞，往往還能起到事半功倍的效果，汶川地震中形成的多個堰塞壩，就是通過爆破拆除而成功地化險為夷，避免了壩體突然潰決而導致重大洪災的發生，保證了下游百姓的生命和財產安全[1－2]。因此，開展爆炸荷載下堤壩破壞形態的研究，具有重要的實用意義。

巖土體在爆炸荷載下的破壞形態最直觀地表現為爆炸彈坑。以往的爆炸彈坑研究大多集中于水平場地，而專門針對堤壩爆炸彈坑的研究尚不多見。Schmidt和 Holsapple[3]通過理論分析和試驗驗證認為，在離心機中模擬爆炸是采用室內小模型模擬原型的最為有效的方法，論證了模型與原型在爆炸條件下的相似關系，并就爆炸彈坑問題進行了一系列研究；穆朝民等[4－5]進行了自由場地爆炸彈坑的試驗與數值分析；Sauseville[6]、Zimmie 等[7－8]對堤壩爆炸彈坑及爆炸后的穩定性等進行了離心機試驗模擬；劉軍等[9]基于 LS-DYNA軟件模擬了土石壩在爆炸荷載作用下的動力響應，但由于其施加的爆炸荷載為TM5-585-1手冊確定的經驗時程曲線，對觸地爆炸下空氣沖擊波引起的感生地沖擊未作考慮，因此，與實際情況有所偏差，特別是在近地表區域，直接地沖擊和感生地沖擊將發生復雜的疊加和混合，對分析模擬結果的影響更不宜忽略[10]。

影響堤壩爆炸破壞效應的因素有很多[11]，本文首先選取炸藥埋深進行分析；同時，筑壩材料本身的性質也不容忽略，例如，濕土較之干土在爆炸荷載下的動力特性就有很大區別，動載作用下含水土體的孔隙水壓力發展會對土體結構的動力響應產生巨大影響，這種孔壓效應在以往的爆炸彈坑研究中尚不多見。因此，本文在LS-DYNA軟件框架內，建立空氣、炸藥與均質堤壩3種物質的模型，綜合考慮了爆炸的直接地沖擊與感生地沖擊，利用多物質ALE法，分別就影響爆炸彈坑形成的兩個重要因素——炸藥埋深和土體孔隙水壓力進行模擬分析，通過觀察后處理中的物質流動變形，得到堤壩在不同爆炸工況下的彈坑破壞形態，以期能夠對水利工程中的堤壩安全設計以及爆破拆除堰塞壩等工程應用提供參考。

2 算法和材料模型

2.1 多物質ALE算法

ALE方法是在材料域和空間域外引入參考域，在參考網格上進行控制方程的求解。它綜合了拉格朗日方法與歐拉方法的優點，既解決了拉格朗日方法下材料的嚴重變形，又克服了歐拉方法下移動邊界引起的復雜性問題，為爆炸沖擊波與結構的相互作用問題提供了較好的解決方案。

ALE方法分為單物質ALE方法和多物質ALE方法。多物質ALE方法容許在一個網格中包含多種物質材料，通過跟蹤每種材料的邊界，在相應的單元中進行物質交換和輸送[12]。爆炸荷載下的堤壩動力響應問題涉及到土體、空氣和炸藥多種物質的大變形問題，需要同時處理界面大變形和網格大變形，而多物質 ALE方法的以下幾個特點就很好地解決了這些問題：

（1）允許一個網格內含有多種物質，物質界面可以穿過網格，即炸藥爆炸產物、土體和空氣能夠相互侵入原來彼此所占據的空間，形成新的物質分布形態。因此，可以通過物質的流動狀態來判斷彈坑的形成過程及最終形態，而不需要通過網格變形或者單元刪除的方法。

（2）用重分重映處理網格的大變形，可以避免爆炸作用下單元網格畸變可能導致的計算中斷。

（3）用level set方法、VOF方法等界面追蹤方法處理界面的大變形，克服了歐拉法難以精確描述移動邊界的缺陷。

因此，利用多物質ALE法能夠較為合理地對堤壩爆炸彈坑進行模擬分析。

2.2 土體材料模型

數值計算的核心是本構模型，其對計算結果的準確性有重要影響。LS-DYNA中的土體本構模型[13]MAT_FHWA_SOIL可以模擬諸多材料特性，如爆炸沖擊加載下的應變率效應、應變硬（軟）化、考慮孔隙效應的彈性本構、修正Mohr-Coulomb屈服面以及孔隙水壓力效應等，是較為合適的計算土體爆炸響應的模型。它的優越性已在文獻[14]中充分驗證，因此，本文選取其來描述均質土壩筑壩材料在爆炸荷載下的力學行為。

為了確保數值計算的有效性和穩定性，當剪應力較小時，該模型將標準Mohr-Coulomb屈服面修正為一個光滑的曲面，并且垂直于壓力軸，其表達式為

式中：P為壓力；φ為內摩擦角；J2為偏應力張量第二不變量；K(θ)為張量平面角的函數；c為黏聚力；AHYP為決定修正后的 Mohr-Coulomb屈服面和標準的Mohr-Coulomb屈服面相似程度的參數。

當AHYP=0時，式（1）表示的是標準Mohr-Coulomb屈服面，當AHYP采用較大值時，修正后的屈服面明顯偏離標準Mohr-Coulomb屈服面，對于數值模擬來講，AHYP的取值應該小于ccotφ，一般按照式（2）進行選?。?/p>

圖1為標準Mohr-Coulomb屈服面與修正后的屈服面對比。除了在低應力區外，二者幾乎一致。

圖1 標準與修正后的Mohr-Coulomb屈服面對比Fig.1 Standard and modified Mohr-Coulomb yield surfaces

同時，為了令屈服面在低圍壓狀態下的形狀呈現為更符合實際的三角形，將標準庫倫函數K()θ改進為[15－16]

含水土體的孔壓發展對計算結果具有重要影響?？讐旱脑鲩L會降低有效應力，并由此導致抗剪強度的下降。對于非飽和土，當氣體體積在加載過程中被壓縮到0時，孔隙水壓力將開始增加，使得土骨架的有效應力降低。模型采用以下關系式來計算孔隙水壓力u：

式中：εv為體積應變；Ksk為孔壓系數[14]；ncur為當前孔隙率；D2為氣體孔隙坍塌前控制孔隙水壓力的材料常數。

參數Ksk影響空氣孔隙坍塌后的孔壓-體積應變曲線曲率。當D2相對于Ksk非常大時，在體積應變達到空氣孔隙體積之前，孔壓幾乎為0；但當D2降低時，孔壓開始增加。D2通過Skempton法中計算孔隙水壓力的參數B得到：

式中：K為體積模量；n為孔隙率；s為飽和度。對于飽和土體，D2=0[17]，于是式（4）簡化為

因此，參數Ksk直接影響土體孔壓的發展，進而影響到土體有效應力的變化。本文將通過調整參數Ksk的大小來描述孔壓上升以及有效應力下降的程度，以此來研究孔壓效應對爆炸彈坑的影響。筑壩材料參數值將在各個算例中具體給出。

2.3 炸藥材料模型及參數

采用高能炸藥燃燒模型和 JWL狀態方程來模擬所使用的TNT炸藥[13]。JWL狀態方程能夠精確描述爆炸過程中爆轟產物的壓力、體積、能量特性，其關系式為

式中：A、B、R1、R2、ω為材料參數；E01為爆轟產物單位體積的內能；V為單位體積裝藥產生的爆轟產物的體積；P為爆炸產生的壓力。炸藥材料參數取值分別為[18]：密度為1.63 g/cm3，爆轟速度為6930 m/s，CJ壓力為21 GPa，A=371 GPa，B=3.23 GPa，R1=4.15，R2=0.95，ω=0.38，E01=8 GPa，V=1。

2.4 空氣材料模型及參數

采用空材料（NULL）模型和線性多項式狀態方程表示空氣的本構關系[13]。將空氣介質簡化為非黏性理想氣體，假設沖擊波的膨脹為絕熱過程，其線性多項式狀態方程可簡化為

式中：ρ為空氣密度；γ為絕熱指數；E為單位初始體積的內能。參數取值分別為[13]ρ=1.25 kg/m3，γ=1.4，E=0.25 MPa。

3 數值模擬方案

基于上述的土體本構模型，利用多物質 ALE法，分別就不同炸藥埋深和不同孔壓上升程度的堤壩彈坑效應進行數值模擬，通過多個模擬結果的比較分析，探討不同爆炸工況和不同筑壩材料特性下的堤壩爆炸成坑效應。

按照文獻[6]中的試驗方案建立堤壩模型，壩高為4 m，壩底寬為20 m，壩頂寬為4 m，左、右坡度比都為 1:2。為了減少計算量，將模型簡化成厚度為一個單元的準二維模型，并根據對稱性，只建立了半模型（見圖2），在后處理中再映射成為整體模型。共劃分為8668個節點，8578個單元。設計藥包尺寸為0.4 m×0.4 m×0.1 m，質量為26 kg。藥包位置在壩體中心線上變化。建模時，首先只劃分空氣和土體單元，然后再根據需要，在計算前將某坐標范圍內原來的單元材料屬性修改為炸藥單元，以此方便地設置不同的炸藥埋深。壩頂和壩坡與空氣接壤，設置為自由邊界；實際情況中壩底與無限大的壩基接壤，因此，設置為無反射邊界，以真實地反映爆炸波在邊界上的透射情況。為了便于不同計算工況的比較，選取了壩體中的參考單元H1830與H1569進行跟蹤，記錄該單元的孔壓時程曲線和土體密度變化時程曲線等。

圖2 有限元模型Fig.2 Finite element model

3.1 考慮炸藥埋深變化的彈坑數值模擬

3.1.1 工況介紹

將藥包中心距離壩頂的距離定義為炸藥埋置深度，在壩頂以上為負，壩頂以下為正。工況分別設置為：①埋深為-0.2 m，即藥包位于壩頂表面，為觸地爆炸工況；②埋深為0.3 m；③埋深為0.6 m；④埋深為0.8 m；⑤埋深為1.3 m。炸藥材料參數如前所述。土體的力學行為通過前述的MAT_FHWA_SOIL材料模型來描述，結合相關經驗[13－14,18]，采用的主要參數取值分別為：密度為1.8 g/cm3，土粒相對密度為2.70，體積模量為68.6 MPa，剪切模量為19.4 MPa，內摩擦角為30°，黏聚力為27.3 kPa，屈服面修正系數為34.5 kPa，含水率為0，孔壓系數Ksk為0，即在炸藥埋深變化的研究中采用干土，不考慮孔隙水壓力的發展，進行總應力分析。

3.1.2 結果分析

以壩中軸線為對稱軸，將原來的半模型映射為完整模型，不同炸藥埋深情況下的堤壩爆炸彈坑形態如圖3所示。圖中深色部分為壩體形態，中軸線上的小方形表示炸藥埋設位置。對模擬結果具體分析如下：

（1）當藥包位于壩頂表面時（圖 3(a)，觸地爆炸），大部分爆炸能量都逸散到空氣中去，作用于土體的能量相對較少，因此，只是由于藥包下方的土體介質被壓實而在壩頂形成一個淺坑。

（2）當藥包完全埋置時，由藥包下方土體的壓實和藥包上覆土層的拋擲共同作用產生彈坑。爆炸發生后，壓力波從藥包中心迅速傳出，爆心周圍出現一個近球形的空腔向外擴張，藥包下方土體形成密實的壓縮帶。壓力波到達壩頂自由面后以稀疏波的形式反射回來，其拉伸作用使得藥包上覆土層變得疏松；稀疏波的拉伸以及壓力波和爆生氣體的推動，導致壩頂不斷向上隆起出現鼓包；隨著爆炸產物的進一步發展，爆腔兩側的土體開裂擴展到壩頂，與地表連通，爆炸產物便沿著這些裂隙噴發，并攜帶破碎土塊一齊向上方沖出，于是出現了彈坑的雛形（如圖 3(b)、3(c)、3(d)所示）；此時，由于炸藥埋深的增加，釋放到空氣中的爆炸能量逐漸減少，而直接作用于土體的能量則逐漸增大，因此，彈坑尺寸隨著藥包埋置深度的加深而增大。

（3）當藥包埋深超過某一深度時，釋放到空氣中的能量已經忽略不計，埋深的繼續增加對其不再產生影響；而此時，由于上覆土體壓力較大，主要是由于爆炸的“擴腔”效應使得土體向外擴張，即在爆炸沖擊波作用下，土體質點獲得速度，沿徑向產生位移，炮孔空腔持續擴大；隨著沖擊波的向外傳播，能量迅速衰減，當沖擊波到達壓縮區界面時，對土體的沖擊壓縮過程結束，爆腔擴展至極限值，最終在壩頂產生鼓包，沒有出現明顯的彈坑，而在壩體中則形成一個封閉的地下爆腔（見圖3(e)）。

在實際爆炸彈坑試驗[19]中，部分土體還會回落，形成一個“亂石井”。本文尚未模擬出土塊的回落，但通過多物質ALE法，已經較好地還原了不同埋深條件下爆炸成坑的宏觀現象，反映了數值手段的正確性。因此，在進行堤壩抗爆安全設計或爆破拆除堰塞壩時，能夠通過數值研究手段預先判斷壩體被爆后的形態，為安全設計方案或爆破方案提供一定的指導措施。

在計算中選取了單元H1569（壓密區）與單元H1830（空腔區）進行跟蹤，記錄該單元處的土體密度時程曲線，參考單元的所在位置如圖3(d)所示。圖4為炸藥埋深為0.8 m的工況下，這兩個單元中的土體密度變化的對比，由圖可以看出，藥包下部的單元由于被壓實，土體密度變大，而爆腔區中的單元受到爆炸劇烈的沖擊壓縮作用后，質點向外發生強烈的位移，土體物質被爆炸產物帶走，形成空腔，因此，此單元的土體密度在短暫上升之后迅速降為0。

圖5為不同炸藥埋深工況下堤壩彈坑尺寸的變化關系。隨著埋深的增加，彈坑的直徑和深度都在加大（見圖3(a)～(d)）；但當埋深增加到1.3 m時，無法產生彈坑，而是形成一個封閉的爆腔（見圖3(e)）。由于本文采用的是準二維模型，對裝藥量的影響尚無法進行完全符合實際的定量分析，所以只能夠定性地說明不同炸藥埋深情況下堤壩爆炸彈坑的變化趨勢。

圖4 爆腔與壓密區土體密度變化(炸藥埋深0.8 m)Fig.4 Variations of soil density at explosion cavity and compacting region (burial depth=0.8 m)

圖5 炸藥埋深與彈坑尺寸關系Fig.5 Relationships of burial depth of explosive and crater size

3.2 考慮土體孔壓效應的彈坑數值模擬

3.2.1 工況介紹

幾何模型皆為壩頂觸地爆炸工況中的模型，即藥包中心距離壩頂以上 0.2 m。為了研究爆炸荷載作用下含水土體孔壓的上升、有效應力的下降對堤壩彈坑形態的影響，將 MAT_FHWA_SOIL模型中的孔壓系數分別設置為 Ksk=0（干土）、1、2、5 MPa，其余參數值不變。如前所述，Ksk越大，材料模型中計算的孔隙水壓力就上升得越快，以此來研究土體孔壓發展對爆炸成坑效應的影響。

3.2.2 結果分析

圖6為不同孔壓上升程度下壩頂觸地爆炸形成的彈坑形態。由圖可以看出，隨著孔壓系數Ksk的增大，爆炸彈坑在尺寸上也呈加大的趨勢，如圖 7所示，當孔壓系數Ksk=0、1、2、5 MPa時，堤壩產生的彈坑直徑分別為1.50、1.60、1.68、3.00 m，彈坑深度分別為0.60、0.70、0.75、0.92 m，說明在爆炸的強烈動荷載作用下，土體的孔壓上升導致了土體抗剪強度的下降，加劇了爆炸過程中的土體變形，因此，在爆炸加載下，含水土體較干土會產生更大的彈坑；孔壓上升得越劇烈，這種增大就越強烈。

為了直觀地顯示孔壓系數Ksk對土體孔壓增長的影響，圖 8給出了不同孔壓系數Ksk下壩體中參考單元 H1569（見圖 3）的超孔壓時程曲線，可以看出，當Ksk=0時，超孔壓為0，而Ksk越大，單元的超孔隙水壓力幅值也越大。

因此，在進行爆破拆除堰塞壩等工程應用時，還必須考慮含水土體的孔壓增長對爆破效果產生的影響，可以通過數值手段預先對爆破拆除效果進行分析，合理地控制爆破參數，以避免堰塞壩的壩基或壩體液化而導致的突然潰壩和引發次生災害[2]。同時可以推斷，在爆炸荷載作用下如果堤壩的開裂深度低于水位線，不僅可能導致堤壩失去使用功能，而且將可能造成洪水傾瀉乃至潰壩的嚴重后果，尤其是對于筑壩材料處于飽和狀態的土體，在爆炸強烈動荷載作用下產生的超孔隙水壓力將進一步加劇爆炸成坑的過程。因此，在遭遇戰爭、爆炸襲擊之前，有必要迅速降低庫水水位，以防庫水倒灌，引起重大洪災；同時還可以根據數值手段大致地判斷彈坑可能的形狀和尺寸，對降水措施提供一定的指導。

圖6 不同孔壓上升程度下的爆炸彈坑形態Fig.6 Blasting craters under different values of Ksk

圖7 不同孔壓系數Ksk下的爆炸彈坑尺寸比較Fig.7 Comparison of blasting crater size under different values of Ksk

圖8 不同孔壓系數Ksk下的土體超孔隙水壓力時程比較Fig.8 Comparison of excess pore water pressure time histories under different values of Ksk

4 結 論

（1）利用多物質ALE法的物質流動變形，能夠較好地模擬出堤壩在爆炸荷載作用下的彈坑形態。

（2）炸藥埋深對堤壩爆炸彈坑形態有重要影響，在埋深較淺時，彈坑的尺寸隨著炸藥埋深的增加而增大；但當埋深超過某一深度，則無法形成彈坑，而是由于爆炸“擴腔”效應產生一個封閉的爆腔。

（3）爆炸荷載作用下含水土體的孔隙水壓力增長對彈坑的形成具有重要影響，孔壓上升程度越強，產生的彈坑尺寸越大。

由于本文采用的是準二維模型，計算結果與實際工程的三維工況存在差別，并且尚未考慮拋擲土體回落并覆蓋至實際彈坑而導致的可見彈坑與實際彈坑的區別，因此，只能定性地說明爆炸成坑的宏觀現象；此外，對堤壩爆炸成坑后的壩體穩定性還未做深入分析。這些都需要進一步開展工作。

[1]梁向前,崔亦昊,魏迎奇. 工程爆破在堰塞湖處理中的應用及實例[J]. 水力發電,2009,35(10): 88－90.LIANG Xiang-qian,CUI Yi-hao,WEI Ying-qi. Treat-ments of landslide dam and its cases by blasting[J].Water Power,2009,35(10): 88－90.

[2]張文煊,吳新霞. 堰塞湖搶險爆破對壩基及保留堰體的液化影響分析[J]. 工程爆破,2008,14(3): 21－23.ZHANG Wen-xuan,WU Xin-xia. Analysis of impact on liquefying of dam foundation and retained weir body in emergency blasting of dammed lake[J]. Engineering Blasting,2008,14(3): 21－23.

[3]SCHMIDT R M,HOLSAPPLE K A. Theory and experiments on centrifuge cratering[J]. Journal of Geophysical Research,1980,85(B1): 235－252.

[4]穆朝民,任輝啟,辛凱,等. 變埋深條件下土中爆炸成坑效應[J]. 解放軍理工大學學報(自然科學版),2010,11(2): 112－116.MU Chao-min,REN Hui-qi,XIN Kai,et al. Effects of crater formed by explosion in soils[J]. Journal of PLA University of Science and Technology (Natural Science Edition ),2010,11(2): 112－116.

[5]穆朝民,任輝啟,李永池,等. 變埋深條件下飽和土爆炸能量耦合系數的試驗研究[J]. 巖土力學,2010,31(5):1574－1578.MU Chao-min,REN Hui-qi,LI Yong-chi,et al.Experiment study of explosion energy coupling coefficient with different burial depths in saturated soils[J].Rock and Soil Mechanics,2010,31(5): 1574－1578.

[6]SAUSEVILLE M J. Explosive cratering using a geotechnical centrifuge on dry,saturated,and partially saturated earth embankments and dams[D]. New York:Rensselaer Polytechnic Institute,2005.

[7]ZIMMIE T F,DE A. Centrifuge modeling to study stability of dams[C]//Proceedings of the Annual Conference of the Association of State Dam Safety Officials. Seattle: Association of State Dam Safety Officials,1996: 630－639.

[8]ZIMMIE T F,SAUSVILLE M J,SIMPSON P T,et al.Blasting studies of dams using the geotechnical centrifuge[C]//Proceedings of the ASDSO Annual Conference on Dam Safety. New Orleans: Association of State Dam Safety Officials,2005: 438－445.

[9]劉軍,劉漢龍,張正珺. 爆炸荷載下土石壩動力響應特征的數值模擬[J]. 防災減災工程學報,2010,30(1): 10－16.LIU Jun,LIU Han-long,ZHANG Zheng-jun. Numerical simulation of dynamic response of an earth and rock-filldam to a blast loading[J]. Journal of Disaster Prevention and Mitigation Engineering,2010,30(1): 10－16.

[10]柳錦春,方秦,還毅,等. 炸藥地面接觸爆炸下土中感生地沖擊的實用計算方法[J]. 解放軍理工大學學報(自然科學版),2010,11(2): 121－124.LIU Jin-chun,FANG Qin,HUAN Yi,et al. Practicable calculating method of indirect ground shock in soil at surface contact explosion[J]. Journal of PLA University of Science and Technology (Natural Science Edition ),2010,11(2): 121－124.

[11]FERRERO V H. Further checks for scaling effects on explosion induced craters[D]. Baltimore: University of Maryland,College Park,MD. 1988.

[12]賈祖朋. 基于MOF界面重構的多物質ALE方法[J]. 計算物理,2010,27(3): 353－360.JIA Zu-peng. A multi-material Arbitrary Lagrangian-Eulerian method based on MOF interface reconstruction[J]. Chinese Journal of Computational Physics,2010,27(3): 353－360.

[13]HALLQUIST J Q. LS-DYNA Keyword User’s Manual(971)[M]. California: Livermore Software Technology Corporation,2007.

[14]WAYNE Y Lee. Numerical modeling of blast-induced liquefaction[D]. Provo: Brigham Young University,2006.

[15]KLISINSKI,M. Degradation and Plastic Deformation of Concrete[D]. Warsaw: Polish Academy of Sciences,1985.

[16]KLISINSKI M,MROZ Z. Description of inelastic deformation and degradation of concrete[J]. International Journal of Solids and Structures,1988,24(4):391－416.

[17]LEWIS B A. Manual for LS-DYNA soil material model 147[R]. Washington,D C: Federal Highway Administration,2004.

[18]劉晶波,杜義欣,閆秋實. 地下箱形結構在爆炸沖擊荷載作用下的動力反應分析[J]. 解放軍理工大學學報(自然科學版),2007,8(5): 520－524.LIU Jing-bo,DU Yi-xin,YAN Qiu-shi. Dynamic response of underground box structures subjected to blast load[J].Journal of PLA University of Science and Technology(Natural Science Edition ),2007,8(5): 520－524.

[19]亨利奇. 爆炸動力學及其應用[M]. 熊建國譯. 北京:科學出版社,1987.