基于ALE方法的彈性圓柱殼入水時的流固耦合模擬

施紅輝,周 棟,溫俊生,賈會霞

(浙江理工大學(xué) 機(jī)械與自動控制學(xué)院,杭州 310018)

海上作戰(zhàn)能力是一個國家軍事強(qiáng)弱的重要體現(xiàn),是維護(hù)國家主權(quán)的重要保障。因此,超空泡技術(shù)應(yīng)運(yùn)而生,國際上涌現(xiàn)出很多超空泡應(yīng)用技術(shù),比如俄羅斯超空泡魚雷“疾風(fēng)”,美國海軍作戰(zhàn)中心采用超空泡作戰(zhàn)技術(shù)設(shè)計的“適應(yīng)高速水下彈”的射彈,這些超空泡減阻技術(shù)得到了廣泛的軍事應(yīng)用[1]。此外,物體撞擊水面時會遭受到巨大的沖擊載荷,嚴(yán)重影響其彈道穩(wěn)定性甚至引起結(jié)構(gòu)破壞。因此,物體入水問題的研究對于解決魚雷、反潛導(dǎo)彈、破障炮彈的入水運(yùn)動問題具有至關(guān)重要的意義[2]。

1929年,Von Karman[3]最早提出了平板撞水問題的漸進(jìn)理論,該理論基于線性化自由表面和體邊界條件。1932年,Wagner[4]對Von Karman的公式在考慮“水堆”效應(yīng)的前提下進(jìn)行了修改。到了2004年,Aquelet等[5]通過有限元方法預(yù)測了楔形體撞擊水面時的局部脈沖載荷。Khazraiyan等[6]基于ALE方法,模擬了低速射彈的入水過程,得到了射彈入水過程中壓力的變化以及射彈速度與加速度隨著入水時間的變化。孫琦等[7]采用ALE方法,對三維彈性彈丸的撞水過程進(jìn)行了流固耦合模擬。康德等[8]采用ALE方法,模擬得到了三維長方體破片在水介質(zhì)中高速運(yùn)動時的速度衰減規(guī)律、墩粗變形規(guī)律以及沖擊波傳播過程。段宇等[9]通過具體的實(shí)驗(yàn)圖像分析發(fā)現(xiàn)低速圓柱殼入水之后受到的阻力具有明顯的不同,大致規(guī)律是:平頭圓柱殼受到的阻力最大,其次是圓錐圓柱殼,而圓球頭圓柱殼受到的阻力最小。這一規(guī)律對相關(guān)研究具有很好的借鑒性。路中磊等[10-11]先是做了開放腔體圓柱殼的入水實(shí)驗(yàn),而后又進(jìn)行相關(guān)數(shù)值計算,其進(jìn)行的對比研究發(fā)現(xiàn),入水速度對空泡的流動規(guī)律具有關(guān)鍵的影響,開放腔體圓柱殼的入水空泡隨著速度的增加從波動流動轉(zhuǎn)變?yōu)樵苹鲃印?/p>

朱棒棒、施紅輝等[12]曾經(jīng)研究了剛性圓柱體入水過程中空泡面積、空氣攜帶量以及空泡面閉合時間的變化。溫俊生、施紅輝等[13]研究了頭型對回轉(zhuǎn)體垂直入水時空泡面閉合的影響。這些研究的應(yīng)用背景主要是射彈或炸彈的入水。在研究用導(dǎo)彈入水攻擊潛艇時,必須考慮兩端封閉圓柱殼體的入水問題。因此,本文在之前研究的基礎(chǔ)上,選擇了一組不同頭型彈性薄殼圓柱體,利用ANSYS/LS-DYNA有限元軟件,采用ALE方法,對入水過程進(jìn)行了三維流固耦合模擬,主要對比了不同頭型彈性圓柱殼的變形、最大壓力、入水運(yùn)動規(guī)律,對比了平頭剛性圓柱殼與不同頭型彈性圓柱殼的的入水運(yùn)動過程。

1 數(shù)值模型

1.1 控制方程

本文假設(shè)流體為不可壓縮流體,不可壓縮牛頓流體的質(zhì)量和動量方程的ALE形式可以表征為如下形式[14]:

·v=0

(1)

(2)

式中:t為時間;v為流體速度;w為網(wǎng)格速度;f為流體體力項(xiàng);流體應(yīng)力張量表示為如下形式:

(3)

式中:I為單位矩陣,p為壓力,Re為雷諾數(shù)。

1.2 材料模型

圓柱殼模型材料采用線彈性材料,水采用Null材料模型,圓柱殼模型屈服函數(shù)數(shù)學(xué)描述如下:

σ=Eε

(4)

式中:σ為應(yīng)力分量,E為彈性模量,ε為應(yīng)變分量。狀態(tài)方程采用Gruneisen狀態(tài)方程:

(5)

式中:c為介質(zhì)聲速;S1=1.921,S2=-0.096,S3=0;γ0=0.35,為Gruneisen常數(shù);a為γ0和介質(zhì)壓縮比μ1的一階體積修正量;U1為單位體積內(nèi)能。μ1形式如下:

(6)

式中:ρ0為常溫狀態(tài)下水的初始密度,ρ為水的當(dāng)前密度。

1.3 模型、計算區(qū)域及邊界條件

模型的二維截面如圖1所示。由左至右分別為平頭彈性圓柱殼、90°錐角彈性圓柱殼、120°錐角彈性圓柱殼、平頭圓柱實(shí)體。長度L均為80 mm,殼體厚度δ均為2 mm。材料的彈性模量E=70 GPa,泊松比ν=0.33,密度ρs=2 700 kg/m3。

圖1 模型示意圖

計算區(qū)域示意圖如圖2所示,其中,圖2(a)為二維截面示意圖,圖2(b)為三維正交示意圖。為提高計算效率,建模時,在ANSYS/LS-DYNA中采用cm-g-μs單位制。由于本文的計算模型為均質(zhì)回轉(zhuǎn)體模型,因此本文只建立了四分之一模型進(jìn)行計算。XOZ與YOZ平面設(shè)置為約束邊界條件,流場其余邊界均設(shè)置為無反射邊界條件(排除壁面反射波對計算的干擾)。圓柱殼體直徑Dn=40 mm,初始位置位于其頭部距自由界面10 mm的空氣域中。計算區(qū)域設(shè)置為圓柱域,包括空氣域和水域兩部分,其中,圓柱域半徑為6Dn,空氣域長度為5Dn,水域長度為8Dn。圓柱殼的入水初速度為50 m/s

圖2 計算區(qū)域示意圖

1.4 網(wǎng)格劃分

本文的計算方法采用ALE方法,選取SOLID 164單元對水、空氣及圓柱殼體進(jìn)行定義,采用Lagrange網(wǎng)格對圓柱殼體進(jìn)行離散化,采用Euler網(wǎng)格對流體域進(jìn)行離散化,通過用戶自定義關(guān)鍵字在流體和固體之間施加耦合算法。固體與流體的網(wǎng)格最小尺寸均為2 mm。本文只在自由面及圓柱殼附近加密網(wǎng)格,這樣可大大提高計算效率。網(wǎng)格單元數(shù)如表1所示。

表1 網(wǎng)格單元數(shù)

2 數(shù)值驗(yàn)證

考慮到實(shí)驗(yàn)會存在測量與計算誤差,模擬會存在精度問題,導(dǎo)致實(shí)驗(yàn)與模擬之間會存在差異,筆者通過對文獻(xiàn)[15]中的圓柱殼模型以1.98 m/s的初速度入水過程進(jìn)行計算來驗(yàn)證本文計算方法的準(zhǔn)確性。在圓柱殼下設(shè)置了5個壓力監(jiān)測點(diǎn),監(jiān)測點(diǎn)上的沖擊壓力峰值ps對比如圖3所示,圖中,橫坐標(biāo)為壓力檢測點(diǎn)到圓心的距離d。數(shù)值模擬結(jié)果與文獻(xiàn)[15]中實(shí)驗(yàn)結(jié)果的最大誤差為8.7%,該誤差小于10.0%,本文認(rèn)為這個范圍可以接受,故采用該方法進(jìn)行數(shù)值模擬。

圖3 數(shù)值模擬結(jié)果與文獻(xiàn)[15]實(shí)驗(yàn)結(jié)果的對比

3 計算結(jié)果及分析

3.1 不同頭型彈性圓柱殼入水過程水相圖

不同頭型彈性圓柱殼入水過程水相圖如圖4所示,每幅圖之間的時間間隔為1 ms,t=0 ms對應(yīng)的時間為入水0 ms時刻(即圓柱殼頭部剛接觸水面時刻),依此類推。從圖中可以清晰地觀察到,3種不同頭型彈性圓柱殼入水超空泡的直徑從大到小的順序?yàn)?平頭彈性圓柱殼,120°錐角彈性圓柱殼,90°錐角彈性圓柱殼。

圖4 不同頭型入水過程水相圖

3.2 不同頭型彈性圓柱殼入水過程中的變形情況

當(dāng)彈性圓柱殼從空氣中進(jìn)入水中時,其下表面與上表面均會發(fā)生一定程度的變形。不同頭型彈性圓柱殼入水過程中殼體的變形情況如圖5所示。

圖5 不同頭型彈性圓柱殼變形情況

圖5中每幅圖之間的時間間隔為0.2 ms,t=0.1 ms對應(yīng)的時間為入水0.1 ms時刻,依此類推。對比發(fā)現(xiàn),平頭彈性圓柱殼上、下表面均發(fā)生明顯變形,而錐角彈性圓柱殼(120°與90°錐角)只有上表面變形顯著。

為了更加直觀地呈現(xiàn)不同頭型圓柱殼的變形情況,本文從3種不同頭型的彈性圓柱殼的圓心位置出發(fā),先是對比了平頭彈性圓柱殼上、下表面的變形,然后進(jìn)一步對比了3種不同頭型彈性圓柱殼上表面的變形情況。

圖6對比了平頭彈性圓柱殼上、下表面變形ld。ld>0,表示變形與運(yùn)動方向相反;ld<0,表示變形與運(yùn)動方向相同。下表面的變形存在3個波峰,波峰出現(xiàn)的時間依次為入水0.1 ms時刻、入水0.6 ms時刻、入水1.1 ms時刻。上表面的變形存在2個波峰,波峰出現(xiàn)的時間依次為入水0.1 ms時刻、入水0.7 ms時刻。下表面的振動頻率約為2 000 Hz,上表面的振動頻率約為1 667 Hz。平頭彈性圓柱殼下表面變形大于上表面,上、下表面均表現(xiàn)為凹陷狀態(tài)。

圖6 平頭彈性圓柱殼上、下表面中心點(diǎn)變形對比

圖7對比了不同頭型彈性圓柱殼上表面的變形ld。依然規(guī)定:ld>0,表示變形與運(yùn)動方向相反;ld<0,表示變形與運(yùn)動方向相同。120°錐角和90°錐角彈性圓柱殼上表面的變形與平頭彈性圓柱殼上表面的變形有所不同,它們的變形只存在一個波峰。觀察到平頭彈性圓柱殼上表面的變形程度最大,120°錐角彈性圓柱殼上表面的變形次之,而90°錐角彈性圓柱殼體上表面的變形最小。平頭彈性圓柱殼第1個波峰出現(xiàn)的時間最早,為入水0.1 ms時刻(第2個波峰為0.7 ms時刻);120°錐角彈性圓柱殼波峰出現(xiàn)時間次之,為入水0.3 ms時刻;波峰出現(xiàn)最晚的是90°錐角彈性圓柱殼,其在入水0.5 ms時刻出現(xiàn)波峰。平頭、120°錐角、90°錐角這3種頭型彈性圓柱殼上表面的變形均處于向內(nèi)凹陷的狀態(tài),并隨時間逐漸趨于一致。

圖7 不同頭型彈性圓柱殼上表面中心點(diǎn)的變形對比

3.3 平頭圓柱殼入水過程的壓力變化

為了研究彈性圓柱殼產(chǎn)生變形的原因,本文在圓柱殼下表面沿徑向方向設(shè)置了5個壓力監(jiān)測點(diǎn)分別為r1=40 mm,r2=8 mm,r3=12 mm,r4=16 mm,r5=50 mm,如圖8所示。

圖8 壓力檢測點(diǎn)示意圖

圖9(a)、9(b)分別為平頭彈性圓柱殼下表面和平頭彈性圓柱實(shí)體下表面在每個壓力監(jiān)測點(diǎn)下的壓力(ps)變化圖。

圖9 壓力變化圖

圖9(a)與圖6變形圖中波峰出現(xiàn)的時間一致,且均有3個波峰。對比圖9(a)與圖9(b)發(fā)現(xiàn),殼體下表面受到的壓力大概是實(shí)體的3.5倍?？梢岳斫?在流體沖擊力的作用下,物體高速入水時實(shí)體下發(fā)生的變形要小于殼體,因此實(shí)體受到的壓力小于殼體受到的壓力。由圖9可以觀察到,實(shí)體壓力變化曲線只存在2個波峰,與殼體前2個波峰出現(xiàn)的時間相同。不同頭型彈性圓柱殼最大峰值壓力pmax的對比情況見表2。圓柱殼頭部受到的最大壓力隨著錐角的增大而增大(平頭圓柱殼等效為180°錐角),這說明壓力受幾何形狀與彈性形變的變形影響較大,這與文獻(xiàn)[15]中的結(jié)論是一致的。

表2 不同頭型彈性圓柱殼最大壓力變化

3.4 不同頭型圓柱殼的入水運(yùn)動規(guī)律

不同頭型圓柱殼入水過程中加速度as的對比如圖10所示。圓柱殼撞擊水面時,會瞬間產(chǎn)生巨大的流體沖擊力,其加速度會迅速上升。隨后,圓柱殼的動能傳遞給液面,克服表面張力使液面產(chǎn)生變形,圓柱殼周圍的空氣隨著液體向圓柱殼頭部周圍流動被攜帶進(jìn)入水中,液體發(fā)生汽化,在空氣和水蒸氣的共同作用下,圓柱殼的周圍形成空泡,所以圓柱殼的加速度又會減小。最后,不同頭型彈性圓柱殼的加速度趨向于一個相同的值。平頭彈性圓柱殼加速度在入水0.1 ms時刻達(dá)到峰值,而且峰值最大;其次是平頭剛體圓柱殼,略小于平頭彈性圓柱殼;120°錐角彈性圓柱殼在入水0.2 ms時刻加速度達(dá)到峰值,該峰值約為平頭彈性圓柱殼加速度峰值的一半;90°錐角彈性圓柱殼在入水0.3 ms時刻加速度達(dá)到峰值,其峰值最小。

圖10 不同頭型圓柱殼入水加速度對比

入水過程中速度vs的對比如圖11所示。由圖可見,速度變化規(guī)律與加速度的變化規(guī)律一致。不同頭型圓柱殼入水后速度vs從大到小的順序?yàn)?90°錐角彈性圓柱殼,120°錐角彈性圓柱殼,平頭剛體圓柱殼,平頭彈性圓柱殼。在入水初期,速度均下降明顯,下降趨勢隨時間逐漸變緩。

圖11 不同頭型圓柱殼入水速度對比

在3.3節(jié)中,本文發(fā)現(xiàn),壓力是造成圓柱殼下表面變形的主要原因。本文進(jìn)一步分析發(fā)現(xiàn),由于物體穿過介質(zhì)界面時速度迅速衰減,所以其上表面的變形主要是由慣性作用造成的。

4 結(jié)論

通過對不同頭型彈性圓柱殼入水過程的數(shù)值模擬并與文獻(xiàn)中的結(jié)果進(jìn)行對比,本文得到了如下結(jié)論:

①圓柱殼上表面在慣性作用下發(fā)生變形,下表面在流體沖擊力作用下發(fā)生變形。圓柱殼的上、下表面均處于凹陷狀態(tài)。

②平頭圓柱殼上表面變形曲線存在2個波峰,其振動頻率約為1 667 Hz;下表面變形曲線存在3個波峰,其振動頻率為2 000 Hz。錐角越大,最大變形出現(xiàn)的時間越早,上表面的變形越大。

③平頭彈性圓柱殼變形波峰和壓力波峰出現(xiàn)的時間相一致。平頭圓柱殼下表面壓力曲線存在3個波峰,平頭圓柱實(shí)體下表面壓力曲線存在2個波峰,兩者下表面受到的最大壓力均隨其頭部錐角的增大而增大。

④錐角越大,圓柱殼體的加速度越大,加速度峰值出現(xiàn)的時間越早。剛性圓柱殼入水后的速度略大于彈性圓柱殼。本文的計算結(jié)果表明,圓柱殼的錐角越小,入水時受到的沖擊壓力越小,入水后物體的速度衰減越慢,這也符合一般的常識。