大口徑彈體高速入水載荷特性研究

汪 振,吳茂林,戴文留

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

高速入水問題[1-3]廣泛存在于空投魚雷、水雷等各類軍事工程應用中。隨著海軍的日益發展,水中兵器,特別是新型反潛武器,在未來海戰中將發揮更為重要的作用。此類武器相比投放魚雷、水雷等問題,射程更遠,入水速度相對較高;彈體在高速撞水時會受到巨大的沖擊,這種沖擊載荷會造成彈體結構及其內部器件的破壞。此外,結構體高速入水時的流固耦合作用變得更加復雜,極大地增加了計算結構體入水載荷的難度。

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

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

隨著計算機以及有限元方法的發展,用數值方法研究彈體高速入水問題,是目前新型的研究方法。本文利用有限元仿真軟件LS-DYNA研究該問題。該軟件能夠有效地處理流體和流固耦合問題,經過眾多工程應用的驗證,具有較高的可靠性。黃志剛等[9]基于該軟件計算了平頭錐形回轉體以100 m/s入水時的結構強度;李上明等[10]利用該方法研究了楔形體入水問題。

本文主要以未來海上平臺發射大口徑反潛武器為背景,應用LS-DYNA中多介質ALE方法對大口徑深彈入水進行數值模擬,計算其入水載荷特性,探究了大口徑平頭彈體以不同角度、速度以及攻角入水的載荷特性。

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

1.1 流固耦合方法

LS-DYNA[11-12]程序中實現流固耦合的方法有3種:共節點算法、接觸算法、歐拉-拉格朗日耦合算法(ALE算法)。其中,ALE算法兼具拉格朗日法和歐拉法的優點,既能有效地描述結構邊界運動,又可使內部網格單元獨立于物質實體而存在。網格可以在求解的過程中適當調整,避免出現網格的畸變。

本文所有計算模型網格采用八節點六面體單元。對空氣域和水域采用多介質ALE算法,允許網格中存在多種ALE介質,彈體網格采用Lagrange網格。在進行建模與網格劃分時,彈體的幾何模型以及網格可以與氣體域重合。采用罰函數的約束方式實現液體和固體各個力學參量的傳遞。

罰函數約束方式是通過追蹤彈體拉格朗日節點和歐拉或ALE流體物質位置間的相對位移d。系統運算時檢查節點對流體單元表面的貫穿量,流體所受力F的大小與貫穿量成正比,即

(1)

式中:ki為基于流體、彈體節點質量模型特性的剛度系數;K和V分別為被穿透單元的體積模量和體積;A為連接節點的結構單元平均面積;P為罰函數因子,為設置的PFAC數值。

ALE算法下的控制方程包括質量、動量和能量守恒方程:

(2)

(3)

(4)

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

1.2 模型的建立

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

(5)

(6)

當材料膨脹時,壓力為

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

(7)

參數數據如表1所示。

表1 水與空氣的狀態方程參數

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

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

圖1 彈體入水模型及網格劃分模型

其中彈體尺寸:最大截面半徑r=0.15 m,總長L=1.2 m,頭部為球體的回轉體;為避免截斷的計算域邊界對彈體入水產生影響,水域截斷邊界為彈體半徑的6倍。下面對相關重要參數以及結果進行驗證。

1.3 數值模擬參數驗證

1.3.1 網格驗證

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

不同網格模型對稱面局部放大如圖2所示。

表2 網格種類表

圖2 不同模型局部網格圖

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

圖3 不同網格密度載荷圖

從圖3可以得出,模型2和模型3的網格密度達到一定程度時,數值計算結果因網格引起的誤差在1%的范圍內,可以認為此時的結果已經收斂。對于網格更小的模型3,數值計算時間急劇增加。綜合考慮,既能保證計算結果的準確性又節約計算機資源減少計算時間,本文模型的網格大小采用模型2的網格大小。

1.3.2 PFAC數值驗證

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

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

圖4 不同PFAC數值載荷曲線

圖5 不同PFAC數值滑移能曲線

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

圖6 彈體與水域耦合

2 理想模型的理論驗證

2.1 載荷驗證

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

Lee等[15]對入水過程進行研究時發現,彈體撞水后會進入穩定航行階段,該階段入水結構體耦合界面上的壓力相對較小,入水阻力大小穩定。在建立的入水空泡求解模型中,認為此時結構體入水的動能損失量完全用于排開結構體周圍流體并形成空泡,并且得到了實驗與理論的驗證。考慮一個初速為v的拋射體垂直進入流體,由牛頓第二定律得到沖擊力:

(8)

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

對于球體入水,當馬赫數低于0.5時,穩定航行階段阻力系數CD=0.384。在穩定航行階段,半球頭彈體垂直入水時與球體入水工況基本一致。在200 m/s工況下,球頭彈體入水阻力系數取CD=0.384。在高速入水不計重力的情況下,其入水阻力表達式為

(9)

沖擊加速度表達式為

(10)

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

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

體積為

質量為

m=V1ρs=637.8 kg

加速度為

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

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

球頭彈體垂直入水達到2 000 μs時,已經完全進入穩定航行的開空泡階段。對2 000～6 000 μs的載荷數據取平均,得到的加速度約為880 m/s2,與理論計算值基本一致,說明該方法數值計算準確度較高。

2.2 空泡驗證

另有學者[16]對魚雷垂直入水初期的空泡外形進行了研究,發現除了魚雷頭部沾水以外,入水初期的空泡可以精確地用拋物面表示。其空泡外形為

(11)

式中:db為彈體或魚雷直徑;CD為穩定開空泡階段入水阻力系數。

(12)

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

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

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

圖8 彈體入水空泡

圖9 空泡對比

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

采用上述方法以及調試驗證的數據及參數對某大口徑彈體進行數值計算,探究彈體的入水速度及角度對入水載荷的影響。

建立大口徑彈體幾何模型,并對其進行網格劃分,如圖10所示。

圖10 彈體模型及網格劃分

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

圖11 彈體入水有限元模型

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

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

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

從圖12中可以得出:

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

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

③入水初期,軸向載荷峰值隨著入水角度增大而增大,但到達較為穩定的開空泡階段時,載荷基本一致,說明同一速度下入水角度的改變對峰值載荷影響較大,對穩定載荷影響較小。

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

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

從圖13中可以得出:

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

②徑向載荷的峰值出現在撞水瞬間產生的軸向載荷第1次小峰值時刻。

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

圖14 攻角載荷曲線

從圖14中可以看出:

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

②攻角對徑向載荷峰值影響較大,正攻角徑向載荷峰值達到2 200 m/s2,負攻角徑向載荷峰值只有1 600 m/s2。

③攻角對徑向載荷收斂值有一定影響,在入水初期,正攻角會導致徑向載荷收斂于一個較小的負值,即使彈體產生順時針的力矩,且該力矩隨著攻角的增大而增大。負攻角會導致徑向載荷收斂于一個較小的正值,即使彈體產生逆時針的力矩,且該力矩也隨著攻角的增大而增大。

4 結束語

本文通過對某大口徑彈體以不同速度與角度以及帶有不同攻角的入水模型進行數值模擬分析,得出以下結論:

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

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

③同一速度下入水角度的改變對峰值載荷影響較大,對穩定載荷影響較小。

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

⑤徑向載荷的峰值出現在撞水瞬間產生的軸向載荷第1次小峰值時刻。

⑥正攻角會使彈體產生偏向于空氣中的彎矩;負攻角會使彈體產生偏向于水中的彎矩;載荷大小隨著攻角數值大小的增大而增大。