汪 振,吳茂林,戴文留

(1.海軍工程大學(xué) 兵器工程學(xué)院,湖北 武漢 430033;2.江南工業(yè)集團(tuán)有限公司,湖南 湘潭 411207)

高速入水問(wèn)題[1-3]廣泛存在于空投魚(yú)雷、水雷等各類軍事工程應(yīng)用中。隨著海軍的日益發(fā)展,水中兵器,特別是新型反潛武器,在未來(lái)海戰(zhàn)中將發(fā)揮更為重要的作用。此類武器相比投放魚(yú)雷、水雷等問(wèn)題,射程更遠(yuǎn),入水速度相對(duì)較高;彈體在高速撞水時(shí)會(huì)受到巨大的沖擊,這種沖擊載荷會(huì)造成彈體結(jié)構(gòu)及其內(nèi)部器件的破壞。此外,結(jié)構(gòu)體高速入水時(shí)的流固耦合作用變得更加復(fù)雜,極大地增加了計(jì)算結(jié)構(gòu)體入水載荷的難度。

在入水理論研究方面,Von Karman[4]提出附加質(zhì)量概念并推導(dǎo)了入水載荷的理論計(jì)算公式。Wagner[5]在文獻(xiàn)[4]理論的基礎(chǔ)上提出了斜升角模型的平板理論,該理論成為現(xiàn)今理論研究此類問(wèn)題的基礎(chǔ)。魏卓慧等[6]采用附加質(zhì)量法針對(duì)截錐形彈體入水問(wèn)題建立了理論計(jì)算模型,分析了各種參數(shù)對(duì)入水載荷的影響。

試驗(yàn)研究方面,陳誠(chéng)等[7]研究了超空泡航行器斜入水時(shí)所受沖擊載荷。陸麗睿等[8]研究了不同射彈斜入水空泡及彈道特性。目前,鮮有從理論或試驗(yàn)角度研究平頭倒圓角彈體高速入水時(shí)的載荷特性。

隨著計(jì)算機(jī)以及有限元方法的發(fā)展,用數(shù)值方法研究彈體高速入水問(wèn)題,是目前新型的研究方法。本文利用有限元仿真軟件LS-DYNA研究該問(wèn)題。該軟件能夠有效地處理流體和流固耦合問(wèn)題,經(jīng)過(guò)眾多工程應(yīng)用的驗(yàn)證,具有較高的可靠性。黃志剛等[9]基于該軟件計(jì)算了平頭錐形回轉(zhuǎn)體以100 m/s入水時(shí)的結(jié)構(gòu)強(qiáng)度;李上明等[10]利用該方法研究了楔形體入水問(wèn)題。

本文主要以未來(lái)海上平臺(tái)發(fā)射大口徑反潛武器為背景,應(yīng)用LS-DYNA中多介質(zhì)ALE方法對(duì)大口徑深彈入水進(jìn)行數(shù)值模擬,計(jì)算其入水載荷特性,探究了大口徑平頭彈體以不同角度、速度以及攻角入水的載荷特性。

1 有限元方法及模型的建立

1.1 流固耦合方法

LS-DYNA[11-12]程序中實(shí)現(xiàn)流固耦合的方法有3種:共節(jié)點(diǎn)算法、接觸算法、歐拉-拉格朗日耦合算法(ALE算法)。其中,ALE算法兼具拉格朗日法和歐拉法的優(yōu)點(diǎn),既能有效地描述結(jié)構(gòu)邊界運(yùn)動(dòng),又可使內(nèi)部網(wǎng)格單元獨(dú)立于物質(zhì)實(shí)體而存在。網(wǎng)格可以在求解的過(guò)程中適當(dāng)調(diào)整,避免出現(xiàn)網(wǎng)格的畸變。

本文所有計(jì)算模型網(wǎng)格采用八節(jié)點(diǎn)六面體單元。對(duì)空氣域和水域采用多介質(zhì)ALE算法,允許網(wǎng)格中存在多種ALE介質(zhì),彈體網(wǎng)格采用Lagrange網(wǎng)格。在進(jìn)行建模與網(wǎng)格劃分時(shí),彈體的幾何模型以及網(wǎng)格可以與氣體域重合。采用罰函數(shù)的約束方式實(shí)現(xiàn)液體和固體各個(gè)力學(xué)參量的傳遞。

罰函數(shù)約束方式是通過(guò)追蹤彈體拉格朗日節(jié)點(diǎn)和歐拉或ALE流體物質(zhì)位置間的相對(duì)位移d。系統(tǒng)運(yùn)算時(shí)檢查節(jié)點(diǎn)對(duì)流體單元表面的貫穿量,流體所受力F的大小與貫穿量成正比,即

(1)

式中:ki為基于流體、彈體節(jié)點(diǎn)質(zhì)量模型特性的剛度系數(shù);K和V分別為被穿透單元的體積模量和體積;A為連接節(jié)點(diǎn)的結(jié)構(gòu)單元平均面積;P為罰函數(shù)因子,為設(shè)置的PFAC數(shù)值。

ALE算法下的控制方程包括質(zhì)量、動(dòng)量和能量守恒方程:

(2)

(3)

(4)

式中:ρ為流體密度;t為時(shí)間;vi為流體速度;wi為對(duì)流速度;xi為歐拉坐標(biāo);σij為應(yīng)力張量;bi為作用在流體上的體積力,體積力相對(duì)較小時(shí)可以忽略;e為比總能;下標(biāo)i和j為網(wǎng)格坐標(biāo)。

1.2 模型的建立

水和空氣均選擇材料模型*MAT_NULL來(lái)描述[13-14],狀態(tài)方程均采用Gruneisen狀態(tài)方程來(lái)描述,即

(5)

(6)

當(dāng)材料膨脹時(shí),壓力為

p=ρ0C2μ+(γ0+αμ)E

(7)

參數(shù)數(shù)據(jù)如表1所示。

表1 水與空氣的狀態(tài)方程參數(shù)

彈體的單元類型SOLID164,算法上采用Lagrange算法,選擇模型*MAT_PLASTIC _KINEMATIC來(lái)定義材料。材料為45#鋼,密度為7 850 kg/m3,彈性模量為210 GPa,泊松比為0.3,屈服應(yīng)力為355 MPa,塑性失效應(yīng)變?yōu)?.35。

為節(jié)約計(jì)算資源,建立1/2模型,并對(duì)對(duì)稱面進(jìn)行約束。為保證其模擬無(wú)限水域的情況,空氣及水域非對(duì)稱面采取無(wú)反射邊界條件設(shè)置。模型如圖1所示。

圖1 彈體入水模型及網(wǎng)格劃分模型

其中彈體尺寸:最大截面半徑r=0.15 m,總長(zhǎng)L=1.2 m,頭部為球體的回轉(zhuǎn)體;為避免截?cái)嗟挠?jì)算域邊界對(duì)彈體入水產(chǎn)生影響,水域截?cái)噙吔鐬閺楏w半徑的6倍。下面對(duì)相關(guān)重要參數(shù)以及結(jié)果進(jìn)行驗(yàn)證。

1.3 數(shù)值模擬參數(shù)驗(yàn)證

1.3.1 網(wǎng)格驗(yàn)證

針對(duì)網(wǎng)格尺寸對(duì)計(jì)算精度的影響進(jìn)行研究。為保證計(jì)算過(guò)程的穩(wěn)定性,對(duì)彈體、水域以及空氣域網(wǎng)格采用1∶1的網(wǎng)格比例來(lái)劃分,通過(guò)改變網(wǎng)格單元的大小來(lái)確定較為合適的網(wǎng)格。建立了以下3種精度的網(wǎng)格模型,如表2所示。表中,N為網(wǎng)格數(shù)目,l為單位網(wǎng)格尺寸。

不同網(wǎng)格模型對(duì)稱面局部放大如圖2所示。

表2 網(wǎng)格種類表

圖2 不同模型局部網(wǎng)格圖

將不同模型進(jìn)行數(shù)值計(jì)算,并將沖擊載荷數(shù)據(jù)濾波處理,得到入水沖擊載荷曲線,如圖3所示,圖中,an為載荷。

圖3 不同網(wǎng)格密度載荷圖

從圖3可以得出,模型2和模型3的網(wǎng)格密度達(dá)到一定程度時(shí),數(shù)值計(jì)算結(jié)果因網(wǎng)格引起的誤差在1%的范圍內(nèi),可以認(rèn)為此時(shí)的結(jié)果已經(jīng)收斂。對(duì)于網(wǎng)格更小的模型3,數(shù)值計(jì)算時(shí)間急劇增加。綜合考慮,既能保證計(jì)算結(jié)果的準(zhǔn)確性又節(jié)約計(jì)算機(jī)資源減少計(jì)算時(shí)間,本文模型的網(wǎng)格大小采用模型2的網(wǎng)格大小。

1.3.2 PFAC數(shù)值驗(yàn)證

PFAC數(shù)值設(shè)置大小是流固耦合計(jì)算載荷的一個(gè)重要參數(shù),合理的數(shù)值大小才能使得計(jì)算的結(jié)果更加合理真實(shí)。對(duì)PFAC數(shù)值f分別取0.05,0.1,0.15,0.2,0.4進(jìn)行測(cè)試,計(jì)算的載荷曲線如圖4所示。

從圖4中可以看出,峰值載荷并不是隨著PFAC數(shù)值的增加而增加而是呈現(xiàn)出先增加后減小的趨勢(shì),但是載荷數(shù)值波動(dòng)在一個(gè)較小的范圍內(nèi)。同時(shí)隨著PFAC數(shù)值的增加,載荷峰值到來(lái)的時(shí)間也在向前推移。當(dāng)PFAC數(shù)值f=0.1,0.15,0.2時(shí),載荷峰值大小的變化在1%的范圍內(nèi),特別是PFAC數(shù)值f=0.15和f=0.2時(shí),載荷曲線基本重合,在f=0.1時(shí)取值最大,但是不能簡(jiǎn)單地取極大值,要結(jié)合能量變化來(lái)確定較為合理的取值。下面通過(guò)接觸滑移能來(lái)進(jìn)一步判斷,在無(wú)摩擦的情況下,滑移能應(yīng)該穩(wěn)定在0附近,不能出現(xiàn)數(shù)值較大的正滑移能和負(fù)滑移能,滑移能EH的變化如圖5所示。

圖4 不同PFAC數(shù)值載荷曲線

圖5 不同PFAC數(shù)值滑移能曲線

從圖5中可以看出,當(dāng)PFAC數(shù)值過(guò)小,為0.05時(shí),滑移能為較大的正值,此時(shí)計(jì)算的耦合力不足,并且產(chǎn)生滲透;PFAC數(shù)值過(guò)大,為0.4時(shí),滑移能數(shù)值上為較大的負(fù)值,同時(shí)容易產(chǎn)生負(fù)體積,導(dǎo)致計(jì)算失敗。PFAC數(shù)值f=0.1,0.15,0.2時(shí),滑移能變化趨勢(shì)基本一致,數(shù)值上穩(wěn)定在0附近的一個(gè)較小負(fù)值,此時(shí)取值已經(jīng)較為合理。觀察載荷數(shù)據(jù),f=0.1,0.15時(shí),載荷基本完全一致。本文取PFAC數(shù)值f=0.15。彈體與水域的耦合如圖6所示,圖6(a)為f=0.05部分液體發(fā)生滲透的情況,圖6(b)為f=0.15時(shí)無(wú)滲透的現(xiàn)象。

圖6 彈體與水域耦合

2 理想模型的理論驗(yàn)證

2.1 載荷驗(yàn)證

根據(jù)前文建立的彈體垂直入水模型,彈體半徑為0.15 m,總長(zhǎng)1.2 m,圓柱段h=1.05 m,頭部為半球的回轉(zhuǎn)體。對(duì)其進(jìn)行數(shù)值計(jì)算,得到其入水載荷曲線,并與相關(guān)理論進(jìn)行對(duì)比,以驗(yàn)證模型的正確性及精度。

Lee等[15]對(duì)入水過(guò)程進(jìn)行研究時(shí)發(fā)現(xiàn),彈體撞水后會(huì)進(jìn)入穩(wěn)定航行階段,該階段入水結(jié)構(gòu)體耦合界面上的壓力相對(duì)較小,入水阻力大小穩(wěn)定。在建立的入水空泡求解模型中,認(rèn)為此時(shí)結(jié)構(gòu)體入水的動(dòng)能損失量完全用于排開(kāi)結(jié)構(gòu)體周圍流體并形成空泡,并且得到了實(shí)驗(yàn)與理論的驗(yàn)證。考慮一個(gè)初速為v的拋射體垂直進(jìn)入流體,由牛頓第二定律得到?jīng)_擊力:

(8)

式中:m為拋射體的質(zhì)量;z為侵徹方向;t為時(shí)間;ρ為流體密度;Amax為拋射體最大截面積;CD=CD(v)是與速度有關(guān)的阻力系數(shù)。計(jì)算低速?zèng)_擊時(shí)的阻力常用恒定的阻力系數(shù)。速度越高,阻力系數(shù)越大,阻力系數(shù)取決于馬赫數(shù)。

對(duì)于球體入水,當(dāng)馬赫數(shù)低于0.5時(shí),穩(wěn)定航行階段阻力系數(shù)CD=0.384。在穩(wěn)定航行階段,半球頭彈體垂直入水時(shí)與球體入水工況基本一致。在200 m/s工況下,球頭彈體入水阻力系數(shù)取CD=0.384。在高速入水不計(jì)重力的情況下,其入水阻力表達(dá)式為

(9)

沖擊加速度表達(dá)式為

(10)

式中:r為最大截面半徑。

下面利用該理論進(jìn)行計(jì)算:彈體為45#鋼,ρs=7 850 kg/m3,海水ρ=1 025 kg/m3,半徑r=0.15 m,總長(zhǎng)L=1.2 m,圓柱段長(zhǎng)度h=1.05 m,入水速度v=200 m/s。

體積為

質(zhì)量為

m=V1ρs=637.8 kg

加速度為

針對(duì)上述球頭彈體垂直入水的有限元模型,將數(shù)值計(jì)算的載荷數(shù)據(jù)進(jìn)行平滑濾波處理,結(jié)果如圖7所示。

圖7 球頭彈體入水載荷曲線

球頭彈體垂直入水達(dá)到2 000 μs時(shí),已經(jīng)完全進(jìn)入穩(wěn)定航行的開(kāi)空泡階段。對(duì)2 000～6 000 μs的載荷數(shù)據(jù)取平均,得到的加速度約為880 m/s2,與理論計(jì)算值基本一致,說(shuō)明該方法數(shù)值計(jì)算準(zhǔn)確度較高。

2.2 空泡驗(yàn)證

另有學(xué)者[16]對(duì)魚(yú)雷垂直入水初期的空泡外形進(jìn)行了研究,發(fā)現(xiàn)除了魚(yú)雷頭部沾水以外,入水初期的空泡可以精確地用拋物面表示。其空泡外形為

(11)

式中:db為彈體或魚(yú)雷直徑;CD為穩(wěn)定開(kāi)空泡階段入水阻力系數(shù)。

(12)

由于魚(yú)雷入水頭部為球頭,可以利用上述理論與球頭彈體入水進(jìn)行空泡外形理論計(jì)算。計(jì)算過(guò)程如下:海水密度ρ=1 025 kg/m3;入水速度v=200 m/s;彈體直徑db=0.3 m;最大截面Amax=0.070 65 m3;彈體體積V=0.081 2 m3;穩(wěn)定航行加速度a=880 m/s2;45#鋼密度ρs=7 850 kg/m3;穩(wěn)定航行阻力:FD=ρsVa=5.5×105N。

由式(12)得到穩(wěn)定開(kāi)空泡階段阻力系數(shù)CD=0.383 6。由式(11)得出理論計(jì)算空泡外形曲線為y=9.24x2,單位為m。

數(shù)值計(jì)算得,6 600 μs時(shí)彈體尾部與水面平齊,此時(shí)空泡外形如圖8所示。以彈頭頂點(diǎn)為左邊原點(diǎn)建立橫軸為x、縱軸為y的坐標(biāo)系,將理論計(jì)算得到的空泡外形以及數(shù)值模擬得到的空泡外形在該坐標(biāo)系中進(jìn)行對(duì)比,結(jié)果如圖9所示。兩者曲線基本一致,說(shuō)明了數(shù)值模擬計(jì)算精度較高。

圖8 彈體入水空泡

圖9 空泡對(duì)比

3 大口徑彈體高速入水仿真

采用上述方法以及調(diào)試驗(yàn)證的數(shù)據(jù)及參數(shù)對(duì)某大口徑彈體進(jìn)行數(shù)值計(jì)算,探究彈體的入水速度及角度對(duì)入水載荷的影響。

建立大口徑彈體幾何模型,并對(duì)其進(jìn)行網(wǎng)格劃分,如圖10所示。

圖10 彈體模型及網(wǎng)格劃分

彈體入水有限元模型如圖11所示。

圖11 彈體入水有限元模型

對(duì)于該彈體,在不考慮攻角的情況下,研究入水角度以及入水速度對(duì)彈體入水時(shí)所受載荷的影響。以入水角度β為45°,50°,55°,60°,入水速度在150～190 m/s之間每隔10 m/s取一個(gè)速度值進(jìn)行數(shù)值計(jì)算。

圖12(a)～圖12(e)分別為彈體在不同角度,在同一速度下入水時(shí)所受到的徑向載荷和軸向載荷曲線。圖中,在an=0上方的4條曲線為軸向載荷,方向與入水速度方向相反;下方4條曲線為徑向載荷,方向垂直于彈軸,在圖11中使彈體產(chǎn)生順時(shí)針力矩的為負(fù),使彈體產(chǎn)生逆時(shí)針力矩的為正。

圖12 不同速度彈體入水載荷曲線

從圖12中可以得出:

①在同一入水速度的情況下,其軸向載荷隨入水角度的增加而增加;徑向載荷峰值隨著入水角度的改變有一定的波動(dòng)。

②入水初期,徑向載荷達(dá)到峰值后慢慢收斂于0,入水角度較大的徑向載荷曲線收斂于0的速度快于入水角度小的曲線。

③入水初期,軸向載荷峰值隨著入水角度增大而增大,但到達(dá)較為穩(wěn)定的開(kāi)空泡階段時(shí),載荷基本一致,說(shuō)明同一速度下入水角度的改變對(duì)峰值載荷影響較大,對(duì)穩(wěn)定載荷影響較小。

圖13(a)～圖13(b)分別為速度相同,彈體在不同角度下入水時(shí)所受到的徑向載荷和軸向載荷曲線。在an=0上方的5條曲線為軸向載荷,與入水速度方向相反,下方5條曲線為徑向載荷,方向垂直于彈軸。在圖11中使彈體產(chǎn)生順時(shí)針旋轉(zhuǎn)的力矩為負(fù),使彈體產(chǎn)生逆時(shí)針旋轉(zhuǎn)的力矩為正。

圖13 不同角度彈體入水載荷曲線

從圖13中可以得出:

①?gòu)楏w在同一入水角度的情況下,徑向載荷收斂于0的速度基本一致,但徑向載荷峰值大小隨著入水角度的增加而增加。

②徑向載荷的峰值出現(xiàn)在撞水瞬間產(chǎn)生的軸向載荷第1次小峰值時(shí)刻。

由于彈體飛行過(guò)程中速度與彈軸的偏差會(huì)產(chǎn)生一定的攻角,設(shè)計(jì)氣動(dòng)特性良好的彈體攻角一般較小。為了進(jìn)一步研究攻角對(duì)彈體入水時(shí)徑向載荷以及軸向載荷的影響,針對(duì)彈體以速度180 m/s,彈軸角45°入水模型,改變速度角度,使速度矢量與彈軸矢量產(chǎn)生正、負(fù)3°,5°,7°的攻角,來(lái)研究攻角δ對(duì)載荷的影響,計(jì)算結(jié)果如圖14所示。圖14(a)為正攻角曲線,圖14(b)為負(fù)攻角曲線。圖中,an=0上方的曲線為軸向載荷曲線,下方的曲線為徑向載荷曲線。

圖14 攻角載荷曲線

從圖14中可以看出:

①對(duì)于具有攻角的彈體入水模型。在較小攻角的情況下,攻角對(duì)彈體的整體入水載荷的影響較小,載荷曲線基本一致。但由于攻角引起的彈體垂向速度分量不一樣,而導(dǎo)致正攻角的入水載荷峰值更大,接近7 000 m/s2,負(fù)攻角的入水載荷峰值只有6 000 m/s2左右。

②攻角對(duì)徑向載荷峰值影響較大,正攻角徑向載荷峰值達(dá)到2 200 m/s2,負(fù)攻角徑向載荷峰值只有1 600 m/s2。

③攻角對(duì)徑向載荷收斂值有一定影響,在入水初期,正攻角會(huì)導(dǎo)致徑向載荷收斂于一個(gè)較小的負(fù)值,即使彈體產(chǎn)生順時(shí)針的力矩,且該力矩隨著攻角的增大而增大。負(fù)攻角會(huì)導(dǎo)致徑向載荷收斂于一個(gè)較小的正值,即使彈體產(chǎn)生逆時(shí)針的力矩,且該力矩也隨著攻角的增大而增大。

4 結(jié)束語(yǔ)

本文通過(guò)對(duì)某大口徑彈體以不同速度與角度以及帶有不同攻角的入水模型進(jìn)行數(shù)值模擬分析,得出以下結(jié)論:

①在同一入水速度的情況下,其軸向載荷隨入水角度的增加而增加;徑向載荷峰值隨著入水角度的改變有一定的波動(dòng)。

②入水初期,徑向載荷會(huì)慢慢收斂于0,且入水角度較大的模型徑向載荷收斂于0的速度較快。

③同一速度下入水角度的改變對(duì)峰值載荷影響較大,對(duì)穩(wěn)定載荷影響較小。

④彈體在不同速度同一入水角度的情況下,徑向載荷曲線收斂于0的速度基本一致。

⑤徑向載荷的峰值出現(xiàn)在撞水瞬間產(chǎn)生的軸向載荷第1次小峰值時(shí)刻。

⑥正攻角會(huì)使彈體產(chǎn)生偏向于空氣中的彎矩;負(fù)攻角會(huì)使彈體產(chǎn)生偏向于水中的彎矩;載荷大小隨著攻角數(shù)值大小的增大而增大。