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

施紅輝,周 棟,溫俊生,賈會(huì)霞

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

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

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

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

1 數(shù)值模型

1.1 控制方程

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

·v=0

(1)

(2)

式中:t為時(shí)間;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 模型、計(jì)算區(qū)域及邊界條件

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

圖1 模型示意圖

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

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

1.4 網(wǎng)格劃分

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

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

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

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

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

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

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

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

圖4 不同頭型入水過(guò)程水相圖

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

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

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

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

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

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

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

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

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

3.3 平頭圓柱殼入水過(guò)程的壓力變化

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

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

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

圖9 壓力變化圖

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

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

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

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

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

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

圖11 不同頭型圓柱殼入水速度對(duì)比

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

4 結(jié)論

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

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

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

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

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