跨介質航行器彈性舵翼空化流固耦合仿真分析

杜特專，黃仁芳，王 暢

(1.中國科學院力學研究所流固耦合系統(tǒng)力學重點實驗室，北京 100190;2.中國科學院大學工程科學學院，北京 100049)

0 引言

跨介質航行器高速入水過程中，在與水接觸的短時間內航行器頭部產(chǎn)生短時高幅值沖擊波，沖擊波在航行器內傳播，誘導航行器頭部塑性變形。入水之后航行器與水介質相互作用，誘導與航行器側面接觸的水介質空化演化與脫落，航行器在頭部流體阻力及水介質空化演化的共同作用下航行。與此同時，航行器在時空分布的水動力載荷作用下，結構發(fā)生變形，并產(chǎn)生長時的振動效應。航行器入水過程中速度較高，翼面、舵面等關鍵結構表面能夠產(chǎn)生空泡。空泡的非定常演化及潰滅會在結構表面形成復雜的壓力脈動，促使結構發(fā)生顯著的動態(tài)響應。同時，航行器在穿越自由表面過程、浸水特性的變化會顯著影響結構自身的動態(tài)特性，與通常的機翼氣動彈性問題具有很大的區(qū)別，給耦合分析帶來很大的難度。

隨著跨介質高速航行器輕量化和高速的需求，空化流固耦合問題逐漸成為重要問題。針對包含空化的流固耦合問題，已有研究主要集中于柔性水翼和復合材料螺旋槳等對象。Young[1-2]利用邊界元方法模擬空化流動，結合有限元分析結構的變形，模擬了水輪機葉片在空化和無空化條件下的耦合效應，并利用實驗進行了驗證，分析了復合材料螺旋槳的性能。Amromin等[3]利用梁模型求解水翼振動，分析了附著空化條件下水翼彈性振動的特征。Wu等[4-6]利用實驗和數(shù)值模擬研究了柔性NACA66水翼的空化流激振動幅值與不同空化類型的關系。Ducoin等[7]運用商業(yè)流體力學計算軟件CFX 建立了2自由度的動網(wǎng)格模型來模擬柔性水翼的變形，分析了柔性水翼動力學的規(guī)律。任懷勛[8]基于商業(yè)軟件ANSYS及腳本控制對水翼空化問題進行了流固耦合計算，分析了三維水翼在空化流動作用下的流體動力和結構載荷特性。Campbell等[9]采用開源軟件OpenFOAM對大變形的柔性葉片進行數(shù)值模擬，與實驗結果吻合良好。已有研究工作為空化流固耦合提供了有效的研究方法，加深了研究人員對于空化流動作用下結構振動特征的認識，為進一步深入研究空化流固耦合效應奠定了良好的基礎。

當前對于超空泡流動條件下結構與水汽界面的相互作用的研究還較少，翼/舵與水氣界面的撞擊產(chǎn)生的沖擊力可能激起較大的結構振動，給結構安全性帶來了很大的挑戰(zhàn)。本文基于徑向基函數(shù)方法建立了流固耦合計算方法，對航行器舵翼在超空化流動下的動態(tài)響應問題進行數(shù)值仿真，分析姿態(tài)角(攻角和側滑角)對空泡形態(tài)以及舵翼的振動響應。

1 計算模型

流固耦合仿真根據(jù)求解方法可以分為3種：一是以固體求解器為主，流動過程采用簡化模型處理；二是以流體求解器為主，結構振動采用簡化的模型進行處理；三是流體和結構分別采用不同的求解器，利用數(shù)據(jù)接口達成耦合計算。當前對于空化流固耦合問題的分析主要以后兩種方法為主。本文采用第2類方法，以流體求解器為主建立松耦合的計算模型。帶空化的流固耦合計算模型主要包括流體計算模型、結構模型和流固耦合接口，下面將分別介紹這3個模型(接口)及其控制方程。

1.1 流體模型

流動包含了多相界面、水汽之間的相變以及湍流流動等問題，具體控制方程如下

(1)

(2)

(3)

(4)

ρ=αvρv+(1-αv)ρl

(5)

μ=αvμv+(1-αv)μl

(6)

式中，下標l和v分別代表液態(tài)水和水蒸氣。

本文采用基于有限體積法的SIMPLE方法求解流體方程組，時間采用二階有界中心差分格式，壓力、動量和能量的插值格式分別采用PRESTO、有界中心差分格式、二階迎風格式。

1.2 結構模型

本文關注的舵翼截面為三角型，底邊長3.2 mm，高18 mm，展長40 mm。舵翼的材料為304不銹鋼，彈性模量E=194 GPa，密度為7 930 kg/m3。在結構計算中，將舵翼簡化成懸臂梁，并采用模態(tài)疊加法進行求解。懸臂梁固有頻率ωi和對應的質量歸一化的模態(tài)φi(x)如下

(7)

φi(x)=cosβix-coshβix-?i(sinβix-sinhβix)

(8)

(9)

其中，β1l=1.875,β2l=4.694,β3l=7.855。結構振動方程采用4階龍格-庫塔方法進行求解。

1.3 耦合接口及動網(wǎng)格方法

流場和結構的耦合求解是通過流固界面上的數(shù)據(jù)交換來實現(xiàn)的，但是在流固界面上流場網(wǎng)格和結構網(wǎng)格并不是一一對應的，這就需要通過插值來傳遞位移和受力的信息。本文的插值方法基于徑向基函數(shù)(RBF)[10]，其插值的一般形式為

首先，培養(yǎng)學生勤于思考的習慣。愛因斯坦說“學習知識要善于思考、思考、再思考，我就是靠這個學習方法成為科學家的”，孔子也有名言“學而不思則罔”，所以在學習中思考發(fā)揮著至關重要的作用。思考的過程是對新學知識的吸收、消化、升華的過程。因此，提高閱讀的效果，就需要培養(yǎng)學生勤于思考的閱讀習慣。

(10)

(11)

(12)

式(10)具有N+4個未知系數(shù)，將N個結構節(jié)點的位置和位移代入式(10)可以建立N個方程組，為了使該方程組有解，還需補充4個方程

(13)

通過這N+4個方程可以求出系數(shù)γi和αi，則得到流體節(jié)點的位移為

Δxa=Aasλ

(14)

其中,

(15)

(16)

前面提到，流場網(wǎng)格節(jié)點數(shù)通常較多，對于復雜工程問題，把所有耦合界面上的流體節(jié)點全都取為控制點將會造成計算量過大。對此，可采用貪婪算法減少物面插值基點數(shù)目，從而提高計算效率。該方法在空化流固耦合計算中能夠得到較好的結果[11]。

1.4 幾何模型及網(wǎng)格劃分

本文研究對象的主體是圓臺和圓柱的組合體，在主體尾部對稱分布兩個舵翼(圖1)。

圖1 航行器外型及舵翼幾何外型Fig.1 The geometry of vehicle and the rudder wing

這里不考慮航行器的耦合運動，只分析均勻來流條件下舵翼的流激振動。航行器主體和Z軸負向的舵翼設置為剛體，Z軸正向舵翼設置為彈性體。計算域網(wǎng)格采用分塊結構化網(wǎng)格劃分(圖2、圖3),能夠更好地捕捉水氣界面的形狀和演化。近壁面邊界層網(wǎng)格保證y+≈30。

圖2 航行器網(wǎng)格劃分Fig.2 The grids for the vehicle

圖3 舵翼網(wǎng)格劃分Fig.3 The grids around the rudder wing

2 計算結果

本文模擬了不同來流攻角和側滑角，用以分析不同航行狀態(tài)下空泡與三角型截面翼的相互作用。其中攻角在XOY平面內，側滑角在ZOX平面內，具體計算工況如表1所示。

表1 計算工況Tab.1 Calculation conditions

從計算結果可以看出，在來流攻角為2°～6°時，航行器及舵翼大部分都被包在超空泡內部，僅有錐段末端與水直接接觸(圖4)。隨著攻角的增大，航行器的沾濕面積也不斷增大，但是舵翼仍然處于空泡內。結構振動的模態(tài)位移如圖5所示，由于在給定的參數(shù)范圍內舵翼始終被包裹在空泡內部振動，因此相當于在空氣中振動，受力較小，振動位移也較小，其振動頻率基本上與干模態(tài)基本相同，約為700 Hz。隨著攻角的增大結構振動的平衡位置也不斷增大(圖5)。

圖4 空泡形態(tài)(2°攻角) Fig.4 The supercavity shape (angle of attack is 2°)

(a)不同攻角下的振動時間曲線

(b)2°攻角下舵翼振動頻譜圖5 舵翼1階振動特征 (工況1～3)Fig.5 The first-order vibration displacement of rudder wing (cases 1～3)

當來流攻角為8°時，由于超空泡向Y軸正向偏移，舵翼的自由端(翼梢)會穿出空泡壁，切割空泡界面，并在尾流區(qū)引起了空泡界面的褶皺(圖6紅色箭頭處)。此時，舵翼翼梢的壓力面為沾濕面，吸力面仍然處于空泡之中，舵翼在液態(tài)水和空泡內水蒸氣的共同作用下振動，其振動特征發(fā)生了較大的變化。從圖7(a)中可以看出，舵翼在切割自由面時受到較大的彎矩，產(chǎn)生較大的模態(tài)位移①。隨后舵翼被空泡完全包圍，受力減小，振動以為也隨之降低，直到第2次切割空泡界面，產(chǎn)生第2個峰值②。經(jīng)過兩個周期后，舵翼一直處于切割空泡界面的狀態(tài)③，在附加質量作用下其振動頻率小于干模態(tài)的700 Hz(圖7(b))。 圖8給出了相同條件下彈性舵翼和剛性舵翼的受力對比，可以看到在t=0.17 s時刻剛性舵翼的受力有1個峰值，而彈性舵翼則峰值較小。這是由于彈性翼振動變形對空化形態(tài)造成影響，導致舵翼沾濕面積較小，從而降低其所受水動力。

圖6 舵翼切割空泡(8°攻角)Fig. 6 The supercavity around rudder wing(angle of attack is 8°)

(a)8°攻角下的振動時間曲線

(b)8°攻角下舵翼振動頻譜圖7 舵翼1階振動特征(工況4)Fig.7 The first-order vibration displacement of rudder wing(case 4)

圖8 8°攻角條件下的剛體和彈性體的受力曲線Fig.8 The lift curves of the rigid wing and the elastic wing under 8° angle of attack

圖9 舵翼一階振動特征(工況5)Fig.9 The first-order vibration displacement of rudder wing(case 5)

對于攻角和側滑角均為2°的航行工況，舵翼的沾濕面積相對于工況4來說更小，僅有翼梢頂端一小部分切割空泡壁。每個振動的峰值對應著舵翼切割空泡壁的過程，而低谷則表明舵翼完全處于泡內。該工況條件下，結構振動始終無法達到一個穩(wěn)定的狀態(tài)(圖9)。這說明在該條件下結構的受力變化最為劇烈，比小攻角條件下的受力大了1個量級，可能對航行的穩(wěn)定性帶來較大隱患。

3 結論

本文在FLUENT的基礎上進行二次開發(fā)，基于徑向基函數(shù)實現(xiàn)了空化流激振動的耦合數(shù)據(jù)傳遞及網(wǎng)格變形，利用梁模型模擬舵翼的振動，建立了空化流固耦合計算方法。利用該方法模擬了包含航行器影響的舵翼流激振動過程，獲得了不同攻角條件下空泡形態(tài)和結構振動特征。計算結果表明，在計算空化數(shù)條件下，當攻角較小時(不大于6°)，舵翼完全處于空泡內部，其振動與空氣中的振動特征基本一致。當攻角為8°時，舵翼的自由端切割空泡界面會導致舵翼一開始的振動幅值比小攻角條件大1個量級，隨后在一個較高的平衡位置上穩(wěn)定振動。對于攻角和偏航角均為2°時，其振動無法達到穩(wěn)定狀態(tài)，這是因為舵翼的自由端始終在切割自由面和完全被氣泡包圍兩種狀態(tài)中切換。