基于爆源子結構的爆炸問題多尺度分析方法

李述濤，寶 鑫，劉晶波，陳葉青

(1.清華大學 土木工程系，北京 100084；2.軍事科學院 國防工程研究院，北京 100036)

爆炸荷載具有強沖擊、短持時、大變形等特點，顯式有限元方法是目前進行不同介質中炸藥爆炸及爆轟波傳播數值模擬的主要研究方法。在此類分析中，有限元模型的網格尺寸對計算結果的影響較為顯著[1-5]，通常需要設置細密的網格以捕捉爆炸荷載中的高頻成分。顯式時域逐步積分算法中，滿足穩定性條件的最小時間步長與模型中的最小單元尺寸密切相關[6-8]，網格尺寸與計算效率的平衡長久以來是爆炸問題數值模擬研究中的重要挑戰之一。

在近爆源問題中，若炸藥距離目標較近，一般考慮建立炸藥-介質-目標的整體計算模型，在可以接受的計算成本范圍內，盡可能地提高單元密度，以獲得較為準確的計算結果。而當爆源距離目標較遠時，例如水下爆炸對大壩的影響、巖土中爆炸對深埋隧道的作用等問題，即便通過人工邊界技術截取有限的計算區域，整體模型的空間尺度仍然較大。盡管隨著地沖擊波傳播距離增加，峰值持時會逐漸增大，高頻成分逐漸衰減，目標及其周邊介質的網格尺寸也可以適當增大[9-10]。但由此帶來的單元尺寸過渡問題將給建模帶來一定困難，特別是三維模型，同時，炸藥與周邊介質的局部高密度網格仍將嚴重制約整體模型的計算效率。此時采用整體模型進行計算分析恐難以施行。

解決上述問題的思路之一是將爆源與目標分離，根據有限元理論構建一種能夠轉化爆炸荷載的方法，通過荷載和網格尺度的轉換，間接完成整體模型的計算。近期的一些地震動輸入問題的研究成果給予上述思路以支持和啟發：Liu等[11]根據波場分解理論，提出了基于內部子結構的地震波動輸入方法(internal substructure method,ISM)和基于人工邊界結構的地震波動輸入方法(boundary substructure method,BSM)[12-13]。兩種方法均是利用子結構將外源地震波動轉化為內源形式下的等效地震荷載，將地震動輸入至子結構所包圍的內部區域。Bao等[14]針對深埋地下結構的情況，對內部子結構地震波動輸入方法進行了改進，使用封閉式的內部子結構實現地震波場由外至內的輸入，并針對典型的土-結構相互作用問題進行了計算，驗證了內部子結構方法的實用性[15]。

與地震反應分析不同的是，爆炸過程的數值模擬本身屬于內源波動問題，爆炸波由模型內部產生，并向外界傳播。針對此類問題，可以借鑒內部子結構地震波動輸入方法的思想，將近爆源波場轉化為爆炸等效荷載，并構建適宜的網格尺度過渡的方法，實現大范圍空間的爆炸波傳播模擬。Li等[16]提出一種場地地震反應多尺度分析方法，利用人工邊界子結構將粗網格的自由波場運動轉化為作用在細密網格模型中的等效地震荷載，實現了大范圍工程場地的地震反應計算。

基于以上思路，本文發展了一種基于爆源子結構的爆炸問題多尺度分析方法，該方法將爆炸問題中的整體模型拆分為近爆源小尺度模型和大尺度模型，利用黏彈性人工邊界[17-18]良好的吸波性能，將近爆源小尺度模型中的自由波場運動轉化為大尺度模型中的等效爆炸荷載，從而完成針對大范圍爆炸問題的多尺度計算分析。通過與整體模型直接計算的結果對比，驗證了本文方法的計算精度。本文方法克服了整體模型難以解決的網格過度與計算效率問題，在滿足計算精度的前提下，大幅提高了計算效率，具有較高的實用性。

1 爆源子結構理論

以典型深埋地下隧道結構為例，建立爆炸荷載作用下的二維爆源-介質-目標模型，如圖1所示。在截斷邊界處設置人工邊界單元以吸收外行散射波[19]。為保證計算精度，人工邊界一般設置在距離結構較遠的位置。

圖1 爆源-介質-目標模型示意圖Fig.1 Source-medium-target model

炸藥起爆后，一部分能量導致土體發生相變(固態到液態的變化)和產生塑性變形，一部分耗散于“成坑效應”，除此之外，大部分爆炸能量瞬時轉化為在介質中傳播的應力波，迅速向四周擴散，形成地沖擊效應[20-21]。部分應力波經過復雜介質，傳播到達地下結構處，造成結構損傷破壞，而到達截斷邊界處的地沖擊應力波和經地下結構散射及經自由地表反射的波動將被人工邊界有效吸收。圖1中近爆源區域包含爆炸空腔、塑性破碎區和彈塑性區等，集中了炸藥瞬時爆炸產生的大量爆轟產物和迅速向四周膨脹的高頻沖擊波，其化學反應和動力響應的過程極為復雜[22]。針對該區域的有限元數值模擬需利用流固耦合算法，劃分細密網格，充分考慮介質大變形和材料失效等非線性問題[24-26]。

由于介質阻尼的存在，當地沖擊應力波遠離爆源區域時，應力波的持時逐漸增大，其中的高頻成分逐漸衰減，在有限元計算時，遠離爆源區域的單元尺寸可以適當增大。此外，研究人員所關心的爆炸波與深埋地下結構的動力相互作用一般也發生在這一區域。鑒于此，可基于波場分解理論，將整體模型分離為含有爆源的小尺度內源模型和不含爆源的介質-目標相互作用大尺度模型，并在兩者的交界區域將爆炸波轉化為爆炸荷載，從而實現爆炸過程的多尺度分析。

根據地震波動輸入法理論[27]，入射地震波在轉化為等效輸入荷載時，僅與自由波場有關。同樣，以爆源為中心向四周擴散的地應力波在進行波動-荷載轉化時只與自由場有關，而與地下結構無關。鑒于此，可首先建立與實際爆源-介質-結構模型對應的爆源-自由場模型，該模型為圖1中黃色區域(近爆源區域)，在近爆源區域四周設置二維黏彈性人工邊界單元，如圖2所示。需要注意的是，為保證等效爆炸荷載空間位置的對應，爆源-自由場模型與爆源-介質-結構整體計算模型中對應位置處的網格劃分和材料屬性應保持一致。

圖2 二維爆源-自由場有限元模型示意圖Fig.2 The 2D source-free field finite element model

(1)

(2)

式中，上標0表示自由波場。

對于自由場模型，當應力波通過子結構區域后，可以認為再無反射波進入內部區域，此時假設內部區域節點處于靜止狀態，即

uC=0,uI=0

(3)

將式(3)代入式(2)，得到

(4)

圖3 二維爆源子結構模型Fig.3 The 2D substructure model of explosion source

爆源子結構的運動方程可以寫成如下形式

(5)

式中，上標S表示子結構模型。根據有限元理論，若子結構模型的網格尺寸和材料屬性與自由場模型相應位置處完全一致，那么以下關系成立

(6)

(7)

以上推導得到的爆源子結構方法可利用較小的計算成本獲得近爆源場地的等效爆炸荷載。需要說明的是該方法是基于有限元理論推導得到的，其計算精度由有限元理論支撐。對于爆炸荷載等強沖擊、短持時的荷載形式，有限元顯式時域逐步積分理論中有明確的關于穩定性條件要求，無阻尼和有阻尼情況下，穩定極限分別如下：

無阻尼穩定極限

(8)

有阻尼穩定極限

(9)

式中：ωmax為系統最高頻率；ξ為最高頻率模態的臨界阻尼部分。在實際研究過程中，穩定性條件也可使用以下定義[28]

(10)

式中：Le為模型最小單元尺寸；CP為材料壓縮波速。計算前可利用上式估算滿足穩定性條件的臨界時間步長。

2 基于爆源子結構的多尺度分析方法

在近爆源區域，由于炸藥釋放能量的時間極短，爆轟產物和周邊介質的動力相互作用極為復雜，即便計算步長滿足穩定性條件，計算結果也會隨網格尺寸發生較大變化。一般認為近爆源區域的網格密度越大，計算精度越高。為獲得更加精確的計算結果，需要在爆源子結構的基礎上，對近爆源區域的網格進行加密處理，加密后的爆源子結構局部見圖4。

圖4 二維爆源自由場單元加密模型局部圖Fig.4 Local diagram of the 2D source-free field model with high-density mesh

對網格加密后的近爆源區域進行動力計算，可獲得較為精確的爆炸自由場運動，將爆源子結構所對應的節點位移數據從小尺度的近爆源區域自由場模型中提取出來，再施加到大尺度爆源子結構模型上，通過動力分析獲得等效爆炸荷載。

以上即是多尺度分析的實現方法。需要注意的是，近爆源區域的動力計算是在小尺度的網格上進行的，而爆源子結構是大尺度網格模型，因此在近爆源自由場提取節點位移時，不需遍歷子結構區域對應的所有節點，而是在節點群中“跳躍”式地間隔提取，此操作大幅降低了需要處理的節點數據。隨后再將提取的位移數據一一施加于爆源子結構的對應位置處，通過動力計算獲得大尺度場地-目標模型中的等效爆炸荷載。這種處理方法既體現了多尺度分析的思想，又具有兩點優勢：一是通過“跳躍”式地采集和加載位移數據實現了網格尺度過渡，避免了采用不規則形狀單元的傳統網格過渡方式，既降低了建模難度，又提高了計算精度。二是在有限元理論中，排除網格尺度效應影響。節點位移的提取和加載與單元尺寸無關，而節點力則與單元尺寸相關。本文正是利用了節點位移這一不受網格尺寸影響的物理量，通過爆源子結構將其轉化為節點力，從而保證大尺度場地-目標模型的計算精度。

另外，本文網格尺寸過渡的對象是近爆源區域的拉格朗日網格。利用流固耦合算法對近爆源區域進行數值計算時，需要在炸藥和釋放爆轟產物區域劃分歐拉網格。為保證計算精度，歐拉網格的尺寸一般要小于或等于拉格朗日網格。在近爆源區域，拉格朗日和歐拉網格可能互相重疊。由于歐拉網格節點沒有位移運動，在提取自由場運動數據時需要注意區分。

3 具體操作流程

爆炸問題的多尺度分析方法的核心思想是利用爆源子結構理論，將近爆源區域的場地運動轉化為等效爆炸荷載，同時實現由小尺度網格至大尺度網格的網格尺寸過渡。該方法的操作流程如圖5所示，具體步驟如下：

圖5 二維爆源子結構多尺度方法流程圖Fig.5 The implementation steps of the 2D multiscale method based on the explosion source substructure

步驟1建立近爆源區域的爆源-自由場小尺度模型，模型四周均設置人工邊界條件。根據計算精度要求和穩定性條件進行網格劃分，開展動力計算，提取場地運動的時程數據，即節點A(■)和節點B(◆)的位移時程數據(速度時程或加速度時程亦可)。

4 算例驗證

4.1 均勻半空間自由場算例

首先利用二維自由場算例對本文方法的有效性進行驗證，二維模型如圖6所示，長寬均為500 m，炸藥尺寸為0.4 m×0.4 m，炸藥的幾何中心位于埋深100 m處，模型兩側和底部設置二維黏彈性人工邊界單元。自由場介質參照常用土壤進行參數設置，密度取為2 000 kg/m3，剪切波速為300 m/s，泊松比為0.3。選取A、B、C和D點觀察計算結果。

圖6 二維均勻半空間自由場模型Fig.6 The 2D homogeneous half-space free field model

利用基于爆源子結構的多尺度方法計算整體模型的動力反應時，首先要考慮如何截取近爆源區域。近爆源模型爆炸反應計算需要模擬爆炸波在全無限空間中的輻射，一旦人工邊界的計算精度不夠，外行爆炸波場到達截斷邊界時會產生反射波，直接導致自由波場運動出現誤差，對后續計算分析影響較大。黏彈性人工邊界單元是利用單元矩陣等效原理構造等效單元來模擬黏彈性人工邊界，而黏彈性人工邊界是基于柱面波和球面波理論推導的，在彈性波場中具有較高的計算精度。因此，應盡可能將人工邊界設置在彈塑性區和彈性區的交界處，在降低計算成本的同時，又可保證獲得較為準確的計算結果。本算例截取的近爆源區域的尺寸為60 m×60 m，炸藥位于近爆源區域的幾何中心。近爆源區域模型如圖7所示。

圖7 二維近爆源區域模型圖Fig.7 The 2D near-source area model

首先利用通用有限元軟件LS-DYNA的流固耦合算法對近爆源區域進行動力反應計算。近爆源區域模型采用小尺度網格離散，炸藥和空氣采用Eulerian單元，網格尺寸4 cm×4 cm，土壤采用Lagrange單元，網格尺寸5 cm×5 cm。土壤介質的本構模型采用帶失效的土壤和可壓縮泡沫塑料模型(soil and foam)，該本構模型可以模擬土體大變形和失效行為。土體材料的力學參數為(單位制cm，g，us)：密度2.0，剪切波速0.03，彈性模量0.004 8，剪切模量0.001 8，泊松比0.3。TNT炸藥采用高爆炸藥材料模擬，該材料利用炸藥的爆轟速度、CJ(chapman-jouguet)面壓力等材料性質，聯合狀態方程確定一定當量炸藥在周圍介質中形成的爆炸沖擊波壓力。主要參數為(單位制cm,g,us)：密度1.64，爆轟速度0.693，爆轟波陣面壓力0.27。炸藥爆轟產物狀態方程采用JWL狀態方程，具體參數為(單位制cm，g，us)：單位體積炸藥內能0.07，A=3.14×10-2，B=3.23×10-2，R1=4.15，R2=0.95，ω=0.3。在空氣模型中，通過塊體單元，采用空材料和線性多項式狀態方程來模擬炸藥爆炸產物流動的空腔，進行多物質流固耦合處理[29-31]。

爆源子結構取為正方形，其各邊至炸藥中心的距離為1 800 cm。子結構及不含爆源的自由場模型采用大尺度網格離散，網格尺寸為90 cm×90 cm。選擇爆源子結構外層節點M和節點N、中間層節點P和節點Q觀察近爆源自由場的運動。

此外，為驗證本文方法的準確性，建立考慮爆源的自由場整體模型。為準確模擬爆炸過程及爆炸沖擊波在近爆源區域的傳播，采用小尺度網格對整體模型進行有限元離散，網格尺寸與多尺度方法中的近爆源區域一致，為5 cm×5 cm。

圖8給出了M點、N點、P點和Q點的豎向位移時程曲線。炸藥起爆后，自由場運動在0.2 s左右到達峰值，此后地沖擊應力波開始回落，由于爆炸空腔和塑性區影響，存在一定程度的殘余位移。

圖8 M、N、P和Q點位移時程曲線Fig.8 The displacement history curves on nodes M,N,P,Q

固定爆源子結構模型內側一周的全部節點，將自由場運動的位移時程施加到子結構最外側和中層節點上，計算內側和中層節點的反作用力。圖9給出了子結構內側節點S和節點T、中間層節點P和節點Q的豎向反作用力時程曲線。

圖9 P、Q、S和T點豎向節點反力時程曲線Fig.9 The history curves of the vertical reaction forces on nodes P,Q,S,T

將爆源子結構內側和中間層節點的反力時程數據施加到不含爆源的自由場大尺度模型的對應位置處，計算自由場模型的動力反應。圖10比較了采用基于爆源子結構的多尺度方法和整體方法計算得到的位移波場快照。

圖10 多尺度方法和整體模型計算結果比較Fig.10 Comparison of the results calculated by the multiscale method and the overall model

從對比結果可以看出，當采用基于爆源子結構的多尺度方法進行計算時，自由波場從子結構所在位置開始向外傳播，到達截斷邊界時，波場能量被人工邊界單元有效吸收。由于多尺度方法沒有考慮爆腔對波場運動的影響，自由表面反射波從子結構區域自上而下直接穿行，但從位移云圖上看，此現象并沒有對計算結果造成較大影響。為進一步比較多尺度方法的計算精度，圖11對自由場模型中A、B、C和D點的位移時程數據進行比較。計算結果表明，對于自由場模型，基于爆源子結構的多尺度方法的計算結果與整體模型的計算結果吻合良好。

圖11 多尺度方法和整體模型計算的節點位移比較Fig.11 Comparison of the nodal displacements calculated by the multiscale method and the overall model

4.2 成層半空間自由場算例

進一步驗證本文方法對于成層半空間場地的適用性。算例模型沿豎向分為三層，上層介質的密度為2 000 kg/m3，剪切波速為300 m/s，泊松比為0.3；中層介質的密度為2 500 kg/m3，剪切波速為500 m/s，泊松比為0.27；下層介質的密度為2 700 kg/m3，剪切波速為800 m/s，泊松比為0.26。近爆源區域以及其余模型參數均與均勻半空間算例相同。選取A、B、C和D點觀察計算結果。具體模型和觀測點位置如圖12所示。

圖12 二維成層半空間自由場模型Fig.12 The 2D layered half-space free field model

根據多尺度分析步驟，先建立與上節算例相同的近爆源區域小尺度模型，獲得近爆源自由場運動，隨后建立爆源子結構模型，施加節點位移時程后，通過動力計算獲取節點的反作用力，再施加到對應的成層自由場大尺度模型的相應位置處。

圖13比較了多尺度方法和整體方法計算得到的位移波場快照。自由波場從子結構所在位置開始向外傳播，到達自由表面處發生反射，到達介質分層處同時發生反射和折射，到達截斷邊界位置被人工邊界單元吸收。自由表面和分層界面的反射波均在子結構內部通過，從位移云圖上看，同樣沒有對計算結果造成較大影響。

圖13 多尺度方法和整體模型計算結果比較Fig.13 Comparison of the results calculated by the multiscale method and the overall model

為進一步分析本文方法對于成層介質情況的計算精度，對A、B、C和D點四個觀測點的位移時程數據進行比較。如圖14所示。

圖14 多尺度方法和整體模型計算的節點位移比較Fig.14 Comparison of the nodal displacements calculated by the multiscale method and the overall model

計算結果表明，多尺度方法的計算結果同樣具有較高精度。表1給出多尺度方法和整體模型計算效率的比較。多尺度方法的單元數量比整體模型少95.9%，計算時間減少88.4%，在滿足計算精度的前提下，多尺度方法的計算效率優勢明顯。

表1 多尺度方法和整體模型計算效率比較Tab.1 Comparison of the computational efficiency between the multiscale method and overall model

5 結 論

針對大當量爆炸荷載作用下的大范圍工程場地有限元模型動力計算問題，本文發展了一種基于爆源子結構的多尺度分析方法。該方法在完成爆炸波場轉換的同時，實現網格尺度過渡，從而將大范圍場地整體模型拆分，由近至遠逐步完成動力計算分析。通過二維算例驗證了本文方法的可靠性。本文主要結論如下：

(1)本文提出的基于爆源子結構理論的爆炸荷載多尺度計算方法克服了整體模型單元數量多，計算成本高等問題，不僅具有較高的計算精度，且大幅提高了計算效率。

(2)就多尺度方法本身而言，第一步近爆源自由場地動力計算是耗時最多的。但在炸藥當量給定的情況下，該步計算結果可作為標準荷載，對于不同模型，或是相同模型分析不同位置爆炸時，均可以重復使用，節省大量計算時間。

(3)連續介質波動問題的多尺度分析中，一般采用不規則單元進行網格過渡，這會對計算精度產生一定影響。本文多尺度方法利用爆源子結構實現網格過渡，避免了不規則單元的出現，保證了計算精度。

此外，本文方法在計算過程中涉及大量節點空間信息比對和時程數據的讀寫，手工操作效率低且極易出錯，影響該方法的實用性。本文通過自編程序對限元軟件進行二次開發，利用計算機在較短時間內自動完成模型建模和后處理等工作，極大提高了工作效率，對多尺度方法的推廣應用有重要意義。