基于SPH-FEM方法的地下結(jié)構(gòu)侵徹爆炸數(shù)值模擬研究

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

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

0 引言

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

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

1 SPH-FEM耦合方法

1.1 SPH算法與FEM算法

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

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

(1)

(2)

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

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

FEM算法即傳統(tǒng)有限元方法,在處理小變形問題時具有較高精度和計算速度。但侵徹爆炸近區(qū)具有大變形和非連續(xù)特征,采用FEM方法,容易發(fā)生網(wǎng)格畸變,導(dǎo)致結(jié)果失真或計算終止等問題。

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

1.2 SPH-FEM耦合方式

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

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

2 計算工況及本構(gòu)模型

2.1 計算工況

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

2.2 模型建立

運用傳統(tǒng)有限元方法完成建模,模型埋深10 m,縱深10 m,在侵徹爆炸近區(qū)4 m×4 m×7 m立方體區(qū)域?qū)EM網(wǎng)格置換為SPH粒子,為避免影響計算精度,對SPH粒子附近及拱結(jié)構(gòu)網(wǎng)格適當加密;為實現(xiàn)地下半無限域模擬,模型四周設(shè)置無反射邊界條件。具體布置如圖1所示。

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

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

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

2.3 本構(gòu)模型及參數(shù)

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

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

(3)

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

表2 炸藥的材料參數(shù)Table 2 Material parameters of explosives

鋼筋與混凝土、錨桿與巖石通過關(guān)鍵字*CONSTRAINED_LAGRANGE_IN_SOLID實現(xiàn)相互協(xié)調(diào)變形、共同受力。

3 模型驗證

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

(4)

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

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

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

(5)

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

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

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

以上驗證說明,本構(gòu)模型選擇和參數(shù)確定合理,數(shù)值模擬結(jié)果可信。

4 結(jié)果與分析

4.1 侵徹爆炸聯(lián)合作用過程分析

Mises強度理論認為材料內(nèi)一點的Mises等效應(yīng)力達到某材料極限強度時破壞,其取值為:

(6)

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

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

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

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

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

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

4.3 直墻拱破壞模式分析

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

(7)

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

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

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

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

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

圖9 結(jié)構(gòu)關(guān)鍵點等效應(yīng)力時程曲線Fig.9 Time history curve of effective stress of key points of structure

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

5 結(jié)論

應(yīng)用SPH-FEM耦合方法實現(xiàn)了鉆地彈打擊地下結(jié)構(gòu)侵徹和爆炸連續(xù)全過程數(shù)值模擬,得到了以下結(jié)論:

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

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

3)侵徹階段結(jié)構(gòu)等效應(yīng)力峰值不足爆炸階段的50%,體現(xiàn)了僅侵徹作用對結(jié)構(gòu)毀傷效果有限,但侵徹造成的初始損傷在爆炸聯(lián)合作用下,結(jié)構(gòu)出現(xiàn)多次應(yīng)力峰值,毀傷效果加劇。

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