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

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

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

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

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

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

1 數值模型

1.1 控制方程

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

·v=0

(1)

(2)

式中:t為時間;v為流體速度;w為網格速度;f為流體體力項;流體應力張量表示為如下形式:

(3)

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

1.2 材料模型

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

σ=Eε

(4)

式中:σ為應力分量,E為彈性模量,ε為應變分量。狀態方程采用Gruneisen狀態方程:

(5)

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

(6)

式中:ρ0為常溫狀態下水的初始密度,ρ為水的當前密度。

1.3 模型、計算區域及邊界條件

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

圖1 模型示意圖

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

圖2 計算區域示意圖

1.4 網格劃分

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

表1 網格單元數

2 數值驗證

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

圖3 數值模擬結果與文獻[15]實驗結果的對比

3 計算結果及分析

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖8 壓力檢測點示意圖

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

圖9 壓力變化圖

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

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

3.4 不同頭型圓柱殼的入水運動規律

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

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

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

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

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

4 結論

通過對不同頭型彈性圓柱殼入水過程的數值模擬并與文獻中的結果進行對比,本文得到了如下結論:

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

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

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

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