基于流固耦合的降載空化器高速入水數值研究

李 鴻，李光磊，劉志遠，鄭可心，高曉初，李 濤

(1.哈爾濱工程大學 航天與建筑工程學院，哈爾濱 150001;2.哈爾濱汽輪機廠有限責任公司，哈爾濱 150046)

結構入水是典型的高度非線性瞬態流固耦合過程。由于理論研究對彈性體入水過程求解具有一定局限性，實施等比試驗具有一定的困難，目前國內外學者大多采用數值模擬方法對入水問題進行研究。

Karman[1]基于動量守恒定律，將飛機模型簡化為二維楔形體，建立邊界元模型進行計算，并引入附加質量的概念，得到了入水沖擊載荷的計算方法，為今后的研究奠定了基礎。潘光等[2]使用Msc.Dytran軟件對魚雷入水過程進行了數值仿真，得到了不同工況下的沖擊壓力峰值和壓力分布，對比分析了不同入水速度和不同入水角度對魚雷殼體所受沖擊載荷的影響。Shi等[3]基于任意拉格朗日-歐拉(arbitrary Lagrangian-Eulerian,ALE)算法對水下航行器高速入水過程進行了數值模擬，基于罰函數法計算了耦合面的沖擊載荷，結果表明入水沖擊載荷隨入水速度的增加而增加，且其峰值與速度的平方成正比。Derakhshanian等[4]分別采用有限體積歐拉法、有限差分歐拉法及ALE法對彈體模型斜入水沖擊問題進行了對比研究，同時結合試驗分析了不同慣性矩及中心對動力學特性的影響。孫玉松等[5]基于多介質ALE方法對彈體高速入水沖擊載荷進行了分析，得到了頭型對入水載荷峰值及其持續時間的影響規律，分析了入水角度對入水沖擊載荷的影響。Gao等[6]基于雷諾平均數方程(Reynolds-averaged Navier-Stokes,RANS)和六自由度剛體運動方程對彈體高速斜入水問題進行了研究，結果表明：速度衰減規律主要受彈體長度及頭部形狀的影響，入水速度及角度對其影響不大。Chen等[7]開展了不同頭型的空化器20°斜入水試驗研究，同時結合ANSYS Fluent仿真分析了入水速度與質心的加速度響應及軸向抨擊載荷的耦合規律、阻力系數與入水位移的關系等。盧丙舉等[8]采用混合模型并結合動網格技術對超空泡航行器高速入水過程進行研究，分析了不同入水速度及航行器頭部是否含有通氣口對入水沖擊載荷的影響規律，發現頭部通氣可以有效降低入水瞬間的軸向過載。Panciroli等[9]基于光滑粒子流體動力學(smoothed particle hydrodynamics,SPH)方法并結合試驗研究了楔形體以自由落體方式入水的沖擊載荷，對比分析了剛體模型和彈性體模型對楔形體入水沖擊力的影響。Yan等[10]基于SPH方法對水下航行器高速入水瞬間的沖擊載荷進行研究，結合試驗分析了入水速度、入水角度及攻角等對入水軸向和徑向沖擊加速度的影響。Ahmadzadeh等[11]基于耦合的歐拉-拉格朗日(coupled Eulerian-Lagrangian,CEL)方法對小球自由落體沖擊入水過程進行研究，考慮了層流延伸的黏性效應，忽略湍流效應對結構沖擊載荷的影響，分析了不同密度球體入水的位移及空泡形狀的變化規律。

結構入水瞬間遭遇較高的抨擊載荷會造成結構破壞，過載過高會引起內部零器件的失靈。為實現結構入水緩沖降載，本文采用泡沫鋁材料對空化器進行了入水降載設計；基于MSC.Dytran的一般耦合方法，對比研究了鋁合金空化器與降載空化器的入水響應特性，探討了泡沫鋁的降載效能。

1 理論基礎

1.1 耦合機理

基于一般耦合方法的流固耦合計算通過一層固定在結構上的封閉耦合面傳遞流場壓力及固場位移，不斷交替求解；其中，固場采用基于拉格朗日體系的有限元法，流場采用基于歐拉體系的有限體積法，耦合面兩側的質量及動量輸運滿足守恒關系。

1.2 關鍵技術

對于結構入水這類瞬態動力響應、高度非線性問題的研究，時間域的離散采用顯示時間積分法。計算時間步為n內顯示求解方法的運動微分方程有如下形式

Man=Fe,n-Cvn-Kdn=Fr

(1)

式中：Fe,n為外載荷矢量；Fr為慣性力矢量；M、C與K分別為結構質量陣、阻尼陣與剛度陣;an、vn與dn分別為節點的加速度、速度與位移；加速度可以通過質量陣求逆再與慣性力矢量求積得出，單元質量分布在節點上，各自由度下節點的加速度求解歸結為式(2)的一元一次方程

ani=Fr,ni/Mi

(2)

加速度在各時間步內為常數，采用中心差分法進行時間推進，進一步得到節點的速度及位移響應的表達式為

v(n+1)/2=v(n-1)/2+an(Δt(n+1)/2+Δt(n-1)/2)/2

(3)

dn+1=dn+v(n+1)/2Δt(n+1)/2

(4)

結構高速入水問題中，流體內部會產生十分復雜的沖擊波。在數值求解中，網格的離散造成沖擊波速的不連續，會造成壓力、密度及節點加速度和能量的跳躍。基于顯式中心差分方法對動力學方程進行時間積分，由于沒有算法阻尼的概念會引起高頻的數值噪聲。因此，壓力項中人工體積黏性Q的引入使得強間斷的沖擊波連續化，對抑制數值噪聲至關重要，如

(5)

式中：d為單元特征尺寸；c為材料聲速。本文取CQ=0，CL=0.06。

2 仿真與驗證

2.1 數值模型簡介

降載空化器包括空化盤、泡沫鋁及氣缸3個部件，其初始位置距水面20 mm；具體地，空化盤頭部厚度6 mm，直徑26.6 mm，空化盤長度90 mm，泡沫鋁內外徑分別為3 mm與13.3 mm，氣缸壁厚11.25 mm，其余詳細尺寸見圖1；空化盤與氣缸采用鋁合金材料，空化盤采用DMAT24模型描述其應力應變關系，并通過Cowper-Symonds模型考慮其高速入水抨擊下的應變率響應式(6)，不考慮材料的失效；氣缸采用剛體模型MATRIG，具體參數見表1。

表1 鋁合金材料屬性Tab.1 Material properties of Aluminum Alloy

圖1 數值模型(mm)Fig.1 Numerical Model(mm)

(6)

研究表明，當泡沫鋁應變率低于10 s-1時對其動態應力影響不大，本文采用FOAM1材料模型描述其準靜態壓縮時的應力-應變行為[12]，如圖2所示。該模型的三維空間屈服面為球面，用主應力表示為式(7)，其中半徑Rs=f(Re)，Re滿足式(8)，函數f即泡沫鋁的本構方程。

圖2 泡沫鋁本構模型Fig.2 Constitutive model of foam Aluminum

(7)

(8)

流場模型的描述利用了Graded-Meshes技術(如圖1)，空化器附近及其運動軌跡上的網格進行了加密；全局歐拉域尺寸為0.52 m×0.52 m×1.34 m，水深0.76 m；盒型加密區尺寸為0.20 m×0.20 m×1.16 m，水深0.60 m。結構入水問題涉及跨介質航行，為考慮初始大氣壓強的影響，本文考慮了空氣的作用?？諝鈮毫εc密度及比內能關系采用γ率狀態方程描述(式(9))，γ=1.4，通過給定比內能e=2.068×105J/kg以及密度ρ=1.225 kg/m3進行大氣壓力的初始化；水的壓力p與密度ρ及比內能e的關系采用如下(式(10))多項式狀態方程描述

p=(γ-1)ρe

(9)

p=a1μ+a2μ2+a3μ3+(b0+b1μ+b2μ2+b3μ3)ρ0e

(10)

式中:μ=ρ/ρ0-1；ρ0為初始密度，詳細參數見表2。

表2 水的材料參數Tab.2 Material properties of water

本文考慮流體材料的黏性，其應力張量tij定義為式(11)；sij為應力偏張量(式(12))，ed,ij為應變偏張量，μ為流體黏性系數；p為壓力，其與密度關系滿足式(13)，體積模量K=a1。

tij=-p·δij+sij

(11)

sij=2μ·ded,ij/dt

(12)

dp/dt=K/ρ(dρ/dt)

(13)

2.2 研究內容

相較于斜入水，結構物垂直入水將遭受更高的抨擊載荷。本文以空化器垂直入水為研究背景，分別對鋁合金空化器與降載空化器以不同初速度的入水過程進行了數值模擬，二者構型一致，差別在于鋁合金空化器整體均為鋁合金材質，質量為1.77 kg，降載空化器質量為1.70 kg。表3為具體工況描述，字母A與F分別表示鋁合金空化器與降載空化器，數字表示入水速度。

表3 計算工況Tab.3 Study Case

2.3 數值方法驗證

本文以Chen等的試驗研究為基礎，對入水速度v=70 m/s及90 m/s的工況進行了數值模擬，驗證了一般耦合方法的可靠性，圖3為驗證模型。阻力系數指標定義如式(14)，其中空化器彈頭半徑R=10.55 mm，表4結果顯示誤差在1.2%內。另外，70 m/s情況下的軸向力的對比也基本吻合，如圖4所示。

圖3 驗證模型Fig.3 Verification model

圖4 軸向力對比驗證Fig.4 Comparison and verification of axial force

表4 阻力系數Tab.4 Drag coefficient

CD=F/(0.5ρ0πv2R2)

(14)

3 結果與分析

如下圖5展示了降載空化器180 m/s入水時，從初始位置運動不同距離時的空化效果以及泡沫鋁的壓縮變形。同時給出了泡沫鋁構件特征單元隨降載空化器整體位移的應力響應，如圖6所示，其中單元1靠近空化盤，單元2位于泡沫鋁構件軸向的中間段，單元3靠近氣缸。

(a)x1=0.038 7 m

圖6 泡沫鋁應力響應Fig.6 The stress response of foam Aluminum

結果顯示，當降載空化器位移達到0.038 7 m時，即入水瞬間，靠近空化盤頭部的泡沫鋁(單元1)的應力迅速達到屈服平臺；而后保持平臺應力，當位移達到0.073 9 m時開始發生硬化；隨著入水抨擊載荷的降低，又發生了一定程度的回彈，應力降并伴隨著一定程度的振蕩；隨著遠離空化盤頭部，泡沫鋁的應力響應存在一定滯后；泡沫鋁兩端與鋁合金材料連接處的應力響應幅值高于中間段泡沫鋁。

3.1 加速度響應規律

圖7為空化器不同初速度入水情況下的氣缸加速度與位移耦合規律曲線。從總體上看具有泡沫鋁降載機構的空化器在入水過程中，氣缸的加速度響應會出現數次突變的行為，與泡沫鋁的應力響應趨勢一致。第一次出現在入水瞬間，而后出現加速度平臺并維持一段時間，這是由于泡沫鋁的壓縮變形緩沖所致；當泡沫鋁壓縮硬化到一定程度，氣缸的加速度瞬間發生二次突變；隨著泡沫鋁的回彈，應力程度的降低，氣缸加速度由峰值迅速降低，并伴隨著一定程度的振蕩，經第三次突變達到穩定值。

圖7 加速度與位移關系Fig.7 Acceleration vs displacement

不同初速度入水瞬間產生的第一次加速度突變平臺值基本一致，發生硬化產生二次加速度突變時的位移量相當；入水初速度越大，氣缸的二次加速度峰值越大，且伴隨的微幅振蕩越明顯，波峰衰減到波谷相伴的位移量越大。當速度達到一定程度，會從波谷產生第三次加速度抬升，此后維持穩定值，速度越大，穩定值越大。

鋁合金空化器入水的加速度與位移耦合規律如圖8所示。顯然地，由于鋁合金材質的剛度遠大于泡沫鋁材料，使得其加速度響應振幅過高，經過數個周期快速衰減到穩定值，伴隨的位移量小于降載空化器。

圖8 加速度與位移關系Fig.8 Acceleration vs Displacement

3.2 軸向力響應規律

圖9對比了3組特征速度下(80 m/s，150 m/s，220 m/s)兩種空化器入水的阻力遭遇情況，鋁合金空化器與降載空化器質量相差4%。結果顯示，相同初速度情況下鋁合金空化器入水瞬間遭遇阻力的波峰寬度基本相當，峰值高于具有降載機構的空化器約14%，這是由于泡沫鋁在沖擊瞬間可以起到很好的緩沖作用。入水初速度越大，阻力峰值越大，穩定段阻力的振蕩越明顯。

圖9 軸向力響應Fig.9 Response of axial force

3.3 峰值與平穩值響應

如圖10結果表明，降載空化器與鋁合金空化器的軸向加速度幅值及軸向力幅值與入水速度呈正相關關系；隨著速度的增大，二者軸向力幅值的增長率相當，而加速度幅值的增長率鋁合金空化器大于降載空化器，具體響應統計數據見表5。

圖10 結構響應峰值與入水速度關系Fig.10 Peak value of structure response vs water entry velocity

表5 響應參數統計Tab.5 Response parameter statistics

泡沫鋁的應用對空化器入水降載十分顯著，相較于鋁合金空化器，能夠有效降低軸向峰值加速度至少70%，軸向阻力峰值也略有降低；無論是加速度穩定值還是軸向力穩定值，二者基本相當，差值在±6%。

4 結 論

基于MSC.Dytran一般耦合方法的空化器垂直入水研究，結論如下：

(1)本文空化器構型在80～220 m/s的入水速度內，軸向加速度峰值及軸向力峰值與入水速度呈近似的線性正相關關系。

(2)泡沫鋁的使用能夠有效減緩空化器入水沖擊的軸向加速度及軸向力的振蕩，響應參數緩慢發展到峰值。

(3)泡沫鋁對空化器入水降載十分顯著，軸向加速度幅值降載率至少70%，軸向力幅值降載率可達14%；入水穩定段降載空化器與鋁合金空化器的軸向加速及軸向力相當。