王玉凱，趙建新，韓國(guó)柱

(中國(guó)人民解放軍陸軍工程大學(xué)石家莊校區(qū)火炮工程系，河北 石家莊 050003)

目前，針對(duì)各種侵徹問(wèn)題一般采用有限元法進(jìn)行解決。單元和網(wǎng)格作為有限元方法的基礎(chǔ)，對(duì)數(shù)值計(jì)算的規(guī)模、精度等有不可忽視的影響，為了得到客觀、理想的仿真計(jì)算結(jié)果，研究人員需要分析其變化對(duì)仿真結(jié)果的影響程度。門(mén)建兵等[1]通過(guò)對(duì)比不同尺寸網(wǎng)格的數(shù)值計(jì)算結(jié)果，分析了單元尺寸對(duì)混凝土侵徹?cái)?shù)值模擬的影響，獲得了較為合理的網(wǎng)格尺寸。鄧記松[2]通過(guò)模擬研究了壓力容器的應(yīng)力分析設(shè)計(jì)，明確了單元類(lèi)型、單元技術(shù)和網(wǎng)格密度對(duì)有限元分析結(jié)果有直接影響。林華令等[3]采用拉格朗日網(wǎng)格描述法對(duì)混凝土侵徹問(wèn)題進(jìn)行了數(shù)值模擬分析，考察了拉伸失效模式、網(wǎng)格尺寸和銷(xiāo)蝕應(yīng)變對(duì)模擬結(jié)果的影響，得到的結(jié)論是網(wǎng)格尺寸為5.0 mm時(shí)較為合適。

為保證仿真模擬結(jié)果的客觀性，排除人為設(shè)定因素的影響，本文利用ABAQUS軟件對(duì)槍彈侵徹松木靶板過(guò)程進(jìn)行仿真模擬，探究單元類(lèi)型和網(wǎng)格尺寸對(duì)仿真模擬結(jié)果的影響，明確了適用于槍彈侵徹過(guò)程仿真模擬的單元類(lèi)型和網(wǎng)格尺寸。

1 數(shù)值模擬

1.1 模型構(gòu)建

彈、靶幾何結(jié)構(gòu)及其所受載荷具有明顯的軸對(duì)稱性，建模時(shí)可只構(gòu)建1/4彈、靶的仿真模型，以便于仿真模型的建立和減少計(jì)算量。為了模擬真實(shí)情況，由彈頭殼、鉛套和彈芯三部分組成的彈頭模型參照7.62 mm槍彈彈頭建成，如圖1所示。同時(shí)，構(gòu)建尺寸為250 mm×100 mm×25 mm的松木靶板1/4板模型，其材料參數(shù)可通過(guò)查閱含水量為12%左右的俄羅斯紅松木的資料獲得[4]。根據(jù)彈、靶材料力學(xué)性能及特點(diǎn)，分別采用Johnson-Cook模型和改進(jìn)的Hashin失效模型來(lái)描述彈丸和靶板的動(dòng)態(tài)力學(xué)響應(yīng)。為更準(zhǔn)確地描述槍彈侵徹靶板的損傷破壞過(guò)程，選擇面面侵蝕接觸來(lái)處理槍彈與靶板的接觸問(wèn)題。

圖1 槍彈和靶板幾何形狀

1.2 單元類(lèi)型

ABAQUS/EXPLICIT中有豐富的實(shí)體單元庫(kù)，針對(duì)不同的問(wèn)題可以選擇不同的單元類(lèi)型進(jìn)行解決。對(duì)于三維動(dòng)力學(xué)問(wèn)題，六面體單元和四面體單元是兩種應(yīng)用較為廣泛的單元類(lèi)型。由于槍彈侵徹松木靶板屬于動(dòng)態(tài)力學(xué)大變形問(wèn)題，彈、靶在高速撞擊侵徹過(guò)程中均發(fā)生損傷變形，而松木靶板的幾何構(gòu)型為規(guī)則的長(zhǎng)方體，因此為了便于對(duì)靶板進(jìn)行網(wǎng)格劃分，節(jié)約計(jì)算成本以及保證計(jì)算準(zhǔn)確性，分別選擇三維線性六面體單元C3D8R和三維修正四面體單元C3D10M作為靶板的單元類(lèi)型[5]。

C3D8R單元(圖2(a))是線性六面體減縮積分單元，包含8個(gè)節(jié)點(diǎn)，每個(gè)節(jié)點(diǎn)有6個(gè)自由度，分別為沿x,y,z軸向的平移和繞各個(gè)軸的轉(zhuǎn)動(dòng)。C3D8R單元形狀規(guī)則，具有一定的承受扭曲變形的能力，可以根據(jù)需要設(shè)定扭曲控制的長(zhǎng)度比來(lái)減少單元扭曲變形的影響，劃分出的靶板網(wǎng)格質(zhì)量較高。但是該單元本身存在沙漏數(shù)值問(wèn)題，需要通過(guò)“沙漏控制”設(shè)定“沙漏剛度”來(lái)增加單元的剛度，控制沙漏在模型內(nèi)的擴(kuò)展。在槍彈侵徹松木靶板過(guò)程的仿真模擬中，靶板單元會(huì)發(fā)生嚴(yán)重的扭曲變形，為避免單元扭曲對(duì)計(jì)算結(jié)果的影響，需將扭曲控制的長(zhǎng)度比設(shè)定為0.1。其多項(xiàng)式位移模式為：

u=a1+a2x+a3y+a4z+a5xy+a6yz+a7xz+a8xyz

(1)

v=b1+b2x+b3y+b4z+b5xy+b6yz+b7xz+b8xyz

(2)

w=c1+c2x+c3y+c4z+c5xy+c6yz+c7xz+c8xyz

(3)

式中：u為x方向的節(jié)點(diǎn)位移函數(shù)；v為y方向的節(jié)點(diǎn)位移函數(shù)；w為z方向的節(jié)點(diǎn)位移函數(shù)；ai為x方向節(jié)點(diǎn)位移函數(shù)的插值多項(xiàng)式系數(shù)；bi為y方向節(jié)點(diǎn)位移函數(shù)的插值多項(xiàng)式系數(shù)；ci為z方向節(jié)點(diǎn)位移函數(shù)的插值多項(xiàng)式系數(shù)；x,y,z為整體坐標(biāo)。

自然坐標(biāo)系的形函數(shù)為：

(4)

式中：Ni為形函數(shù)；ξ0=ξiξ,η0=ηiη,ζ0=ζiζ，其中ξ,η,ζ為自然坐標(biāo)，ξi,ηi,ζi為節(jié)點(diǎn)i的自然坐標(biāo)。

C3D10M單元(圖2(b))是修正四面體單元，包含10個(gè)節(jié)點(diǎn)，各節(jié)點(diǎn)的自由度與C3D8R單元相同。該單元同樣需要考慮單元扭曲的影響，設(shè)定合理的扭曲控制參數(shù)，與C3D8R單元不同的是,該單元本身能夠較好地解決沙漏數(shù)值問(wèn)題，不需要通過(guò)設(shè)定沙漏控制參數(shù)來(lái)解決仿真過(guò)程中沙漏的影響。10節(jié)點(diǎn)四面體單元的多項(xiàng)式位移模式為：

u=a1+a2x+a3y+a4z+a5xy+a6yz+a7xz+a8x2+a9y2+a10z2

(5)

v=b1+b2x+b3y+b4z+b5xy+b6yz+b7xz+b8x2+b9y2+b10z2

(6)

w=c1+c2x+c3y+c4z+c5xy+c6yz+c7xz+c8x2+c9y2+c10z2

(7)

自然坐標(biāo)系的形函數(shù)為

Ni=(2Li-1)Li(i=1,2,3,4)

N5=4L1L2

N6=4L2L3

N7=4L1L3

N8=4L4L1

N9=4L4L2

N10=4L4L3

式中：Li為第i節(jié)點(diǎn)單元的面積坐標(biāo)。

圖2 靶板模型單元類(lèi)型

根據(jù)上述兩單元的多項(xiàng)式位移模式和自然坐標(biāo)系的形函數(shù)可以發(fā)現(xiàn)，兩單元的多項(xiàng)式位移插值函數(shù)階次均比較高，由有限元相關(guān)理論可知，基于兩單元的有限元計(jì)算的精度較高，整個(gè)計(jì)算收斂的速度也較快。但是，由于其階次較高，相應(yīng)的計(jì)算量也會(huì)比較大，因此運(yùn)算過(guò)程中計(jì)算精度與計(jì)算量需要平衡好。

1.3 網(wǎng)格劃分

由于Zukas等[6]對(duì)網(wǎng)格劃分的研究結(jié)果表明彈丸半徑方向上的網(wǎng)格數(shù)應(yīng)不少于3個(gè)，因此本文中彈丸半徑上的網(wǎng)格數(shù)劃定為4個(gè)。考慮到靶板實(shí)際的損傷情況以及計(jì)算量，基于Lagrange法應(yīng)用局部網(wǎng)格加密的方法對(duì)彈靶接觸區(qū)域進(jìn)行局部網(wǎng)格細(xì)化，而其他區(qū)域的網(wǎng)格尺寸則相對(duì)較大[7-9]，如圖3所示，其中網(wǎng)格細(xì)化的范圍為20 mm×20 mm。另外，為了便于操作，設(shè)定靶板厚度方向的網(wǎng)格尺寸為1 mm，網(wǎng)格數(shù)量始終為25，其他兩個(gè)方向的網(wǎng)格尺寸相同[10-11]，根據(jù)需要分別設(shè)定為1.00 mm、1.30 mm、4.00 mm、6.67 mm。分別對(duì)這4種尺寸的網(wǎng)格進(jìn)行計(jì)算，綜合分析比較其計(jì)算結(jié)果，找到最合適的網(wǎng)格尺寸。靶板細(xì)化區(qū)域劃分的網(wǎng)格尺寸與單元數(shù)量見(jiàn)表1。

圖3 松木靶板模型網(wǎng)格劃分情況

表1 靶板細(xì)化區(qū)域網(wǎng)格尺寸和數(shù)量

2 計(jì)算結(jié)果

2.1 單元類(lèi)型的影響

圖4所示分別為兩種不同的單元類(lèi)型模擬得到的彈丸速度隨時(shí)間變化的曲線。由圖可以看出，不同單元類(lèi)型模擬得到的彈丸速度衰減過(guò)程以及著彈時(shí)彈丸剩余速度是不同的。C3D8R單元模擬得到的彈丸速度變化是逐步衰減至剩余速度為697.0 m/s，而C3D10M單元模擬得到的彈丸速度變化是突然衰減至剩余速度為675.0 m/s，然后在該值上下波動(dòng)。C3D8R單元模擬得到的彈丸剩余速度值與實(shí)驗(yàn)值698.8 m/s幾乎一致，相差僅為1.8 m/s，而C3D10M單元模擬得到的彈丸剩余速度與實(shí)驗(yàn)值相差23.8 m/s。綜上所述，單元類(lèi)型對(duì)于槍彈侵徹松木靶仿真結(jié)果具有重要影響，C3D8R單元是較為合適的單元類(lèi)型。

圖4 不同單元類(lèi)型的彈丸速度-時(shí)間曲線

2.2 網(wǎng)格尺寸的影響

在分析使用不同網(wǎng)格尺寸模擬彈丸速度-時(shí)間曲線(圖5)時(shí)發(fā)現(xiàn)，仿真得到的4條速度-時(shí)間曲線非常相似，且彈丸著彈時(shí)剩余速度也相差不大，網(wǎng)格尺寸為1.00,1.30,4.00,6.67 mm分別對(duì)應(yīng)的彈丸剩余速度為 697.0,699.6,696.0和695.2 m/s。圖5顯示除了網(wǎng)格尺寸為1.30 mm外，當(dāng)網(wǎng)格尺寸由1.00 mm增加至6.67 mm時(shí)，仿真得到的彈丸剩余速度隨網(wǎng)格尺寸的增加而減小，但衰減幅度卻隨之增大。當(dāng)網(wǎng)格尺寸為1.30 mm時(shí)，仿真得到的彈丸剩余速度與實(shí)驗(yàn)值698.8 m/s更加接近。圖6所示為4種網(wǎng)格尺寸的靶板模型模擬的靶板損傷情況，分析發(fā)現(xiàn)，網(wǎng)格尺寸為1.00 mm和1. 30 mm時(shí)模擬的靶板損傷情況比網(wǎng)格尺寸為4.00 mm和6.67 mm時(shí)更加精細(xì)，彈孔的形狀與實(shí)際也更加接近。

圖5 不同網(wǎng)格尺寸的彈丸速度-時(shí)間曲線

圖6 靶板的損傷情況

3 結(jié)論

本文通過(guò)對(duì)7.62 mm槍彈侵徹松木靶過(guò)程進(jìn)行仿真模擬，分析了單元類(lèi)型和網(wǎng)格尺寸對(duì)仿真結(jié)果的影響，得到如下結(jié)論：

1)運(yùn)用 C3D8R單元?jiǎng)澐志W(wǎng)格的單元數(shù)量明顯少于C3D10M單元，表明模擬仿真時(shí)的計(jì)算量小于C3D10M單元。利用C3D8R單元進(jìn)行網(wǎng)格劃分時(shí)，模擬獲得的彈丸剩余速度與實(shí)驗(yàn)結(jié)果更加接近。

2)單元類(lèi)型一定時(shí)，仿真模擬的彈丸剩余速度和靶板損傷情況隨著網(wǎng)格尺寸的變化而變化。網(wǎng)格尺寸為4.00 mm和6.67 mm時(shí)模擬的結(jié)果與實(shí)際值差別較大，網(wǎng)格尺寸為1.00 mm和 1.30 mm時(shí)模擬的靶板損傷情況與實(shí)際情況很接近，且網(wǎng)格尺寸為1.30 mm時(shí)模擬的彈丸剩余速度值更加接近實(shí)際值。

3)為了得到理想的仿真結(jié)果，較為合理的建模設(shè)定是C3D8R單元和1.30 mm的網(wǎng)格尺寸。