基于SPH-FEM方法的地下結構侵徹爆炸數值模擬研究

張文堂，孫惠香，袁英杰，孫慧穎，康 婷

(1 空軍工程大學航空工程學院，陜西 西安 710038；2 95445部隊，云南 大理 672100；3 中天西北建設集團有限公司，陜西 西安 710077)

0 引言

鉆地武器打擊地下掩體一般分為侵徹和爆炸兩個階段,是高烈度、大規模、高時效的復雜毀傷過程。侵爆試驗受成本高、難度大、危險性大和測量困難限制,難以實時開展并滿足研究需求。數值模擬在侵徹、爆炸領域得到越來越廣泛的應用。地下結構在單一侵徹或爆炸作用下的響應分析已有大量研究[1-8]。隨著Swegle等[9]將SPH方法應用到水下爆炸;Libersky等[10]將SPH方法應用到材料沖擊試驗模擬開始,SPH方法解決了許多固體侵徹、爆炸問題。強洪夫等[11]將SPH方法應用到低速碰撞花崗巖領域并進行了仿真分析,通過對比發現SPH方法具有良好的適應性和準確性;王海兵等[12]在試驗基礎上采用SPH算法較好的模擬了不同彈速撞擊花崗巖的試驗結果;楊廣棟等[13]采用SPH算法先模擬彈體侵徹,再采用Lagrange算法模擬裝藥爆炸作用下混凝土靶的毀傷效應,提出初始侵徹毀傷對混凝土內爆的重要影響,但侵徹成坑后填充炸藥的方法與實際侵爆連續過程仍有差距;胡英國等[14]應用SPH-FEM方法解決了巖體爆破近區精細化數值模擬的問題,與試驗監測結果擬合度較好。現有研究只針對單一侵徹或爆炸作用下結構動力響應,或局限于巖石、混凝土靶單一材料的毀傷分析,對大型地下洞室多材料結構在遭受精確制導武器侵徹爆炸連續作用下結構動力響應研究較少。

在SPH算法中,由于質點之間不存在網格關系,因此可避免大變形時網格扭曲造成的精度破壞問題,為解決FEM算法侵爆近區高烈度、大變形以及連續性模擬難題,文中基于SPH-FEM方法建立了鉆地彈侵爆地下直墻拱結構的模擬模型,模擬了鉆地彈打擊地下結構的侵徹爆炸全過程,通過爆腔應力對比了侵徹爆炸與卸荷爆破的差異,分析了直墻拱結構內部應力情況和鉆地武器侵徹爆炸聯合作用下毀傷模式。

1 SPH-FEM耦合方法

1.1 SPH算法與FEM算法

SPH 算法的基本原理是將連續的流體或固體用相互作用的質點組來描述,各個質點上承載各種物理量,通過求解質點組的動力學方程和跟蹤每個質點的運動軌跡,求得整個系統的力學行為。

SPH算法是模擬流體流動的一種拉格朗日型粒子方法,通過任意分布的粒子來求解具有各種邊界條件的積分方程或偏微分方程。每個粒子代表一個具有獨自物理特性的插值點,用規則的內插函數計算出其力學特性。

(1)

(2)

式中:<>為近似符號;i,j為粒子編號;m,ρ分別為粒子的質量、密度;N為計算區域內粒子總數;Wij=W(|xi-xj|,h)為核函數,h為光滑長度,通常選用三次樣條函數。

SPH算法的無網格自適應特性避免了網格畸變計算終止問題,但由于算法耗時巨大,難以應用于大規模工況的模擬。

FEM算法即傳統有限元方法,在處理小變形問題時具有較高精度和計算速度。但侵徹爆炸近區具有大變形和非連續特征,采用FEM方法,容易發生網格畸變,導致結果失真或計算終止等問題。

兩種算法的優缺點,運用SPH-FEM耦合算法,在侵爆近區使用SPH粒子模擬鉆地彈和巖體的變形與破碎,解決網格畸變問題,在侵爆遠區采用常規FEM方法保證計算精度、縮短計算時間。

1.2 SPH-FEM耦合方式

鉆地彈打擊過程分為侵徹和爆炸兩個階段,為實現侵徹和爆炸沖擊波產生剝落物飛濺及拋射過程,將彈體和侵徹爆炸近區巖體設置成SPH粒子,以解決FEM因變形過大導致網格畸變運算終止的問題。鉆地彈打擊地下工程模擬工況體量龐大,考慮計算時間成本和計算精度要求,在地下工程除侵徹爆炸近區外采用有限元網格,提高計算效率。

運用ANSYS/LS-DYNA求解器實現無網格Lagrange求解器和Lagrange處理器的耦合。首先運用LS-Prepost處理器將已經建立好的數值模型中彈體和侵徹爆炸近區巖體的網格單元轉換為SPH粒子,經過必要的調試和運行[15],SPH粒子離散間距設置為網格單元尺寸的一半。然后通過設置關鍵字*CONTACT_TIED_NODES_TO_SURFACE_CONSTRAINED_OFFSET將巖體損傷近區與遠區耦合,實現應力波從近區向遠區的有效傳遞,通過比對相鄰粒子與單元的應力時程獲取最佳參數,實現SPH-FEM耦合計算。

2 計算工況及本構模型

2.1 計算工況

某地下洞庫為直墻拱結構,埋深10 m,周圍為三級花崗巖,縱深10 m,等間距布置有2 m長加固錨桿,直徑2.2 cm;被覆拱為鋼筋混凝土結構,凈跨為14 m,內拱矢高4.5 m,內墻高2 m,鋼筋直徑為1.8 cm。

2.2 模型建立

運用傳統有限元方法完成建模,模型埋深10 m,縱深10 m,在侵徹爆炸近區4 m×4 m×7 m立方體區域將FEM網格置換為SPH粒子,為避免影響計算精度,對SPH粒子附近及拱結構網格適當加密;為實現地下半無限域模擬,模型四周設置無反射邊界條件。具體布置如圖1所示。

圖1 地下洞庫口部截面圖(單位:cm)Fig.1 Section of underground cavern mouth(unit:cm)

選擇GBU-28鉆地彈作為侵徹爆炸彈[16],彈丸總長5.84 m,圓柱部直徑37 cm,彈體壁厚12.77 mm,導彈總重2 300 kg,其中裝藥406 kg,CRH=R/D=3.0。建立1∶1數值模型如圖2所示,將模型置換為SPH粒子。

圖2 鉆地彈模型Fig.2 Model of ground penetrating projectile

2.3 本構模型及參數

鉆地彈彈體材料為93鎢合金彈,選擇Johnson-Cook本構方程描述[17],參數選擇見表1;巖體選擇HJC本構方程描述,SPH方法中各參數見文獻[11];混凝土材料選擇C&K材料描述[18],鋼筋和錨桿選用率相關塑性隨動強化模型(*MAT_PLASTIC_KINEMATIC),參數見文獻[19]。TNT炸藥采用*MAT_HIGH_EXPLOSIVE_BURN材料模型結合*EOS_JWL狀態方程描述[20]:

表1 鎢合金彈殼的材料參數Table 1 Material parameters of tungsten alloy cartridge case

(3)

式中:P為爆轟壓力;V為相對體積;E0表示初始體積內能;A、B、R1、R2和ω為炸藥常數。炸藥參數取值見表2。

表2 炸藥的材料參數Table 2 Material parameters of explosives

鋼筋與混凝土、錨桿與巖石通過關鍵字*CONSTRAINED_LAGRANGE_IN_SOLID實現相互協調變形、共同受力。

3 模型驗證

由于缺乏同等當量試驗實測資料,采用侵徹深度、爆炸峰值壓力計算公式驗證,以確保模型建立和參數確定的正確性。文獻[21]中通過對比多組國內外巖體侵徹數據發現,Bernard公式(C)對于不同RQD值的巖體均具有較高計算精度。Bernard公式(C)的無量綱形式表示為:

(4)

圖3為彈頭侵徹深度時程曲線和速度時程曲線,由圖可見:t=15 ms時起爆炸藥;彈頭在侵入巖體 時速度值降為0。模擬中彈殼壁薄,剛度比實際剛度小,裝藥強度小,另外從彈頭開始裝藥,與實際裝藥有差異,裝藥強度比彈體強度低,在侵徹過程中,彈體在鉆地后會發生擠壓變形,比實際的剛性彈侵徹下的阻力變大,因此模擬侵徹深度與理論計算值相比較小,誤差為1.9%,模擬精度較高。

圖3 彈頭速度、深度時程曲線Fig.3 Time history curve of warhead velocity and depth

計算爆炸遠區巖體峰值壓力判斷模型有效性其壓力峰值σp可表示為[22]:

(5)

式中:R為爆炸傳播距離;c為介質中沖擊波波速;f為耦合系數;ρ為巖石密度;W為裝藥質量;n為介質衰減系數。

選擇爆炸遠區巖體單元作為研究對象,式(5)計算的單元峰值壓力值與模擬值對比見表3。由表3可知,模擬單元壓力值均大于計算峰值,因為侵徹過程中已經形成初始應力,最大誤差為7.9%。

表3 數值模擬與計算公式結果對比Table 3 Comparison of results of numerical simulation and empirical formula

以上驗證說明,本構模型選擇和參數確定合理,數值模擬結果可信。

4 結果與分析

4.1 侵徹爆炸聯合作用過程分析

Mises強度理論認為材料內一點的Mises等效應力達到某材料極限強度時破壞,其取值為:

(6)

式中:σ1,σ2,σ3分別指第一、第二、第三主應力。

圖4為計算模型的Mises應力云圖,由于采用cm-g-ms單位制計算,圖中所示單位為105MPa。由圖4可見,在侵徹爆炸過程中,不斷有巖體粒子向上向外拋出,侵爆產生的高壓氣團急速壓縮周圍巖體,并在巖體中形成空腔。

圖4 侵徹爆炸不同時刻Mises應力云圖Fig.4 Mises stress nephogram at different moments of penetration and explosion

圖4(a)初始t=8 ms純侵徹時,彈體應力較大,巖體中僅形成細條狀彈坑,此時彈頭處最大等效應力為1.436 GPa,超過屈服強度1.2 GPa,彈頭與炸藥銜接處沿徑向張開,這是因為炸藥與彈殼材料強度數值上相差較大,強沖擊作用下兩種材料界面處出現應力集中,包裹炸藥處合金截面較小使得應力較大而最先屈服,并逐漸發生變形;圖4(b)為t=15 ms時炸藥爆炸,彈尾出現較大應力變化,細條狀彈坑伴隨彈腹擴張向兩側擴展并呈現紡錘狀,這與巖體內柱形裝藥爆炸初期爆腔形狀一致,彈頭繼續向深部侵徹并發生較大變形,這是因為受爆炸壓力和彈體反作用力疊加作用,大部分彈體單元等效應力超過屈服強度;圖4(c)為t=29.5 ms時,鉆地彈主體基本解體,零星破片隨大量巖體向外拋出,應力波繼續向深部巖體傳播,紡錘狀應力云圖擴大,原因是近地表巖體由于受到侵徹初始擠壓影響,密度變大,應力波傳播較快,在橫向傳播距離更長,出現上寬下窄的類紡錘形;圖4(d)為t=46 ms時,拋擲物變少,被覆結構直墻拱底應力發生變化,表明前驅應力波已穿過地下結構。SPH粒子與六面體網格界面處應力云圖連續,進一步證明了巖體損傷近區SPH算法與FEM算法耦合有效,參數設置合理。

4.2 封堵作用對侵徹爆炸的影響

通常FEM算法模擬侵徹爆炸過程是分段進行的,即首先模擬彈體侵徹巖體過程,侵徹結束后再在侵蝕算法形成彈坑處安置等當量炸藥,模擬爆炸過程,該方法存在以下問題:一是對于侵徹效應的模擬采用侵蝕算法,即單元應力或應變達到破壞準則即被刪除,如圖5所示,網格的刪除造成了質量不守恒,整個過程只能大致符合侵徹爆炸的物理過程;二是由于侵徹通道單元的刪除,形成了侵徹孔,對爆炸應力波有卸載作用,類似于卸壓爆破,忽略了彈殼及侵徹中坍落巖體對彈坑的封堵作用,分段模擬的方法形成的Mises應力云圖如圖6所示,形成等效峰值應力為246.8 MPa。實際上,堅硬彈體的高速沖擊對巖體進行了快速而劇烈的擠壓,彈頭近區巖石密度升高,對應單元即使擠壓破碎失效也同樣會對彈體產生握裹阻力。SPH-FEM耦合算法模擬侵徹爆炸過程是連續的,可以考慮侵徹后彈體尾部失效散落巖體對爆炸應力波的封堵作用和彈頭近區巖石的握裹力。考慮封堵作用下爆心等效峰值應力為437 MPa,比卸荷爆破增加了77%。對比圖4(c)和圖6(a)可以發現,侵徹爆炸連續作用下爆腔形狀較卸壓爆破模擬的橢圓爆腔形狀反映了彈頭形狀的擠壓變形和近地端半徑更大,更符合彈體侵徹爆炸的實際情況。

4.3 直墻拱破壞模式分析

拱體混凝土材料采用K&C模型,定義比例損傷因子δ來表示其損傷情況,其表示為:

(7)

式中:λ為材料損傷函數值;λm為材料初始損傷量。LS-Prepost中用標記為“plastic strain”的變量云圖代替比例損傷因子表示的損傷云圖,如圖7所示。

圖7(a)t=7.5 ms初始侵徹階段時,由于拱體與圍巖無空隙,拱頂外側直接承受上部壓力荷載,損傷明顯;圖7(b)t=13 ms時,損傷由拱頂向下蔓延,拱頂損傷穿透結構,拱肩損傷由外向內延伸,破壞區域占1/2,尚未形成貫穿破壞;圖7(c)t=35 ms時,爆炸產生應力波到達結構,損傷進一步加劇,損傷傳遞至直墻底部,由于底板與結構拱分離澆筑,垂直振動時,結構產生豎向變形,應力波產生的結構振動對底板影響不大,未發生明顯破壞;圖7(d)t=47 ms時,底板與直墻相鄰處發生損傷,中心無損傷變化;圖7(e)t=50 ms時通過結構底部局部放大圖可以看到,底板與圍巖間形成了一定空隙。

為進一步了解底板與圍巖間產生空隙原因,選取圖7(e)單元為研究對象,分別繪制直墻上、中和下部側方對應圍巖單元R1、R2和R3的X向位移時程曲線,如圖8(a);底板1/2、1/4和邊緣處下部對應圍巖單元R4、R5和R6的Y向位移時程曲線,如圖8(b);下述X、Y軸正方向如圖7(e)。由圖8(a)可知,R1、R2和R3單元都有先向X軸正向移動至峰值后逐漸減小,隨后產生負向位移,且隨著單元Y負向深度的增加,X軸向位移峰值逐漸變小且時間延后,這表明在應力波作用下,受成拱效應和振動效應影響,圍巖先受壓向X軸正向位移,而后由于圍巖反向振動發生回彈導致其負向位移;由圖8(b)可知,底板下的圍巖先向Y軸負向移動至峰值后減小,隨后又向Y軸正向移動,且距離直墻腳越近位移越大,這說明上部荷載主要由直墻腳處向下傳遞,較大的壓力使得腳部發生較大位移,由于底板與主體拱結構分離澆筑導致底板下圍巖受動載相對較小,底板中心幾乎沒有位移。綜上所述,圍巖與拱結構的相互作用使得底板受到向內擠壓從而發生翹曲,底板與直墻采用分離式澆筑,塑性鉸彎矩變為擠壓力,可以減輕地板的損傷,直墻腳處雙向受壓損傷較大。

圖8 圍巖位移時程曲線Fig.8 Time history curve of surrounding rock displacement

為從單元應力變化情況具體了解上述損傷破壞傳遞情況,圖9給出了拱頂、拱肩、拱腳和直墻底部單元Von Mises等效應力時程曲線。

圖9 結構關鍵點等效應力時程曲線Fig.9 Time history curve of effective stress of key points of structure

在結構破壞過程中,由于侵徹爆炸聯合作用應力波的疊加,各部位應力出現明顯雙峰或多峰情況。拱頂在爆炸主導階段出現3個峰值,分別為9.01 MPa、13.49 MPa和16.05 MPa,拱肩應力峰值較小,均不大于7.71 MPa,但在破壞過程中多次出現;拱腳在爆炸主導階段出現最大峰值(18.36 MPa);直墻底部應力峰值出現在破壞末期,峰值較小(6.90 MPa)。整個過程荷載作用情況復雜,在炸藥起爆前(t=15.0 ms)各部位應力逐漸增加,但最大值不大于峰值的50%,說明純侵徹作用對結構毀傷效果有限;但后期侵徹爆炸兩者應力波的復雜疊加,使結構內部不僅有高應力產生,而且峰值多次出現,對結構抗性提出更高要求。

5 結論

應用SPH-FEM耦合方法實現了鉆地彈打擊地下結構侵徹和爆炸連續全過程數值模擬,得到了以下結論:

1)實現了侵徹爆炸連續全過程作用的模擬,根據經驗公式驗證了模擬結果的有效性,解決了單一FEM方法侵爆近區高烈度、大變形以及不連續模擬難題。

2)侵徹爆炸毀傷過程由于不刪除網格,與卸壓爆破相比,侵爆近區由于彈體擠壓和破碎巖石升壓更大,這一特點體現了侵徹彈體對巖體內爆的封堵效應,彌補了侵蝕算法對侵徹爆炸過程侵徹孔卸荷缺陷,更準確模擬了侵徹-爆炸作用。

3)侵徹階段結構等效應力峰值不足爆炸階段的50%,體現了僅侵徹作用對結構毀傷效果有限,但侵徹造成的初始損傷在爆炸聯合作用下,結構出現多次應力峰值,毀傷效果加劇。

4)拱頂、拱肩和拱腳成為侵徹-爆炸作用下直墻拱結構的薄弱部位,建議底板與直墻分離澆筑,可以減輕底板的損傷。