侵徹過程數(shù)值模擬單元及網(wǎng)格的試驗研究

王玉凱，趙建新，韓國柱

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

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

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

1 數(shù)值模擬

1.1 模型構(gòu)建

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

圖1 槍彈和靶板幾何形狀

1.2 單元類型

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

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

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é)點位移函數(shù)；v為y方向的節(jié)點位移函數(shù)；w為z方向的節(jié)點位移函數(shù)；ai為x方向節(jié)點位移函數(shù)的插值多項式系數(shù)；bi為y方向節(jié)點位移函數(shù)的插值多項式系數(shù)；ci為z方向節(jié)點位移函數(shù)的插值多項式系數(shù)；x,y,z為整體坐標。

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

(4)

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

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

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)

自然坐標系的形函數(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é)點單元的面積坐標。

圖2 靶板模型單元類型

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

1.3 網(wǎng)格劃分

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

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

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

2 計算結(jié)果

2.1 單元類型的影響

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

圖4 不同單元類型的彈丸速度-時間曲線

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

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

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

圖6 靶板的損傷情況

3 結(jié)論

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

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

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

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