羅文莉 陳皓宇

(上海飛機(jī)設(shè)計研究院飛機(jī)結(jié)構(gòu)強(qiáng)度工程技術(shù)所,上海 201210)

水上迫降是指陸基飛行器在水面上的可控緊急降落,屬于典型的入水沖擊問題,世界各國民用運輸機(jī)適航條例均要求跨水域飛行必須通過水上迫降適航審定。早在1929 年,Karman[1]就對水上飛機(jī)入水沖擊的水載荷進(jìn)行過理論研究,早期眾多理論方法都集中于研究二維物體垂直入水問題,但對于復(fù)雜三維民用飛機(jī)的水上迫降問題,理論分析無法得到應(yīng)用。1940 年代開始,國外相繼開展了一系列動力相似模型試驗,用以研究民用飛機(jī)水上迫降特性[2-3]。然而動態(tài)模型試驗具有較強(qiáng)的隨機(jī)性,且無法準(zhǔn)確模擬氣動力的影響,同時也無法較好地解釋迫降過程中的作用機(jī)理。20 世紀(jì)后期,隨著計算理論和計算機(jī)硬件技術(shù)的發(fā)展,成本低廉且展示度高的數(shù)值仿真方法逐步取代了模型試驗。越來越多的研究通過數(shù)值仿真方法模擬飛機(jī)水上迫降,在此過程中不斷總結(jié)提煉[4],逐漸完善仿真方法,使之成為未來研究水上迫降問題的重要手段。

絕大部分的數(shù)值仿真均基于以下幾類方法：面元法(panel method, PM),有限元法(finite element method, FEM),光滑粒子水動力學(xué)法(smooth particle hydrodynamics, SPH) 和有限體積法(finite volume method, FVM)。2001 年,Vladimir 等[5]基于動量守恒原理采用PM 對福克614 飛機(jī)的水上迫降壓力、變形和漂浮特性開展了研究。2008 年,Shoji等[6]基于LY-DYNA 軟件使用FEM 對縮比機(jī)身等直框段垂直入水進(jìn)行了模擬,壓力峰值數(shù)據(jù)和試驗吻合較好但時間歷程相差較大。2010 年,賀謙等[7]研究了水平速度、下沉率、俯仰角對迫降時飛機(jī)機(jī)身載荷的影響,但仿真方法未經(jīng)驗證。2013 年,徐文岷等[8]采用FEM 考慮了水、空氣和飛機(jī)結(jié)構(gòu)之間的耦合作用,對民用支線飛機(jī)水上迫降進(jìn)行了模擬,但得到的入水初期運動姿態(tài)與試驗相差較大。SPH 方法由于其不依賴網(wǎng)格的特點優(yōu)勢在水上迫降仿真中得到了應(yīng)用。2006 年,Climent 等[9]針對CN-235 飛機(jī)的研究表明在水上迫降這一類以水平速度為主導(dǎo)的運動中,以本構(gòu)方程為控制方程的SPH 法無法模擬出氣穴效應(yīng)和負(fù)壓力,從而影響對運動姿態(tài)和沖擊力的模擬。之后,文獻(xiàn)[10-12] 開發(fā)了基于弱可壓縮連續(xù)方程為控制方程的SPH,模擬負(fù)壓的精度得到了提升。

PM 計算的原理一般是基于動量守恒或勢流理論,忽略黏性、可壓縮性、重力和氣墊效應(yīng)等等,多用于簡單外形的物體入水問題,在精確描述民用飛機(jī)水上迫降動力學(xué)問題時存在一定難度。FEM 絕大部分基于商用有限元分析軟件,對于氣動力的模擬存在困難,另外,計算結(jié)果依賴于流固耦合算法,常見的問題是容易引起流體滲漏,導(dǎo)致流體和固體載荷不統(tǒng)一。SPH 避免了由于網(wǎng)格劃分產(chǎn)生的問題,但在流體與固體的邊界處理上存在較大難度,在模擬精度上存在劣勢,且建模復(fù)雜度較高。實際水上迫降過程中,飛行速度必須維持在可以保持飛機(jī)平穩(wěn)下降的水平,氣動力的影響無法忽略。此外,由空氣和水面壓差導(dǎo)致的機(jī)身尾部的吸力也是影響迫降運動狀態(tài)和全機(jī)載荷的重要因素,因此FVM 具有較大的優(yōu)勢。

2007 年,Streckwall 等[13]使用基于FVM 和流體體積(volume of fluid, VOF)法的COMET 軟件對飛機(jī)水上迫降過程進(jìn)行了模擬,對比試驗表明計算結(jié)果較為準(zhǔn)確,并直觀體現(xiàn)了飛機(jī)入水過程中的水面運動和飛機(jī)表面壓強(qiáng)分布；2012 年,Zhang 等[14]模擬了一種支線客機(jī)的水上迫降過程,表明吸力對于飛機(jī)水上迫降的姿態(tài)影響較大,且空氣模型是吸力產(chǎn)生的必要條件。同一時期,屈秋林等[15]和Guo等[16]詳細(xì)地研究了尾吊高平尾、翼吊低平尾和翼身融合三種不同布局飛機(jī)的水上迫降特性,認(rèn)為正常布局飛機(jī)的低平尾具有抑制飛機(jī)過度上仰的作用,從而有利于減緩過載和局部壓力。2019 年,吳宗成等[17]基于FVM 采用滑移動網(wǎng)格研究了波浪對水上迫降特性的影響。由于水上迫降過程中的水平位移較大,目前基于FVM 的方法大多采用動網(wǎng)格實現(xiàn)對運動過程的模擬,通過網(wǎng)格變形和重構(gòu)模擬飛機(jī)和計算域的運動,需消耗龐大的計算資源。本文在已有FVM 方法的基礎(chǔ)上采用全流場運動網(wǎng)格的方法,添加約束保持水面高度不變,使流場相對水面自由運動,大大縮減了計算量。首先通過典型算例對該方法進(jìn)行驗證,然后模擬分析了某型飛機(jī)的水上迫降運動過程。

1 仿真方法

本文采用ANSYS FLUENT 軟件,基于非定常RANS 方程,使用加強(qiáng)壁面處理的可實現(xiàn)k-ε湍流模型,壓力速度耦合采用SIMPLEC 算法。各控制方程中擴(kuò)散項采用二階中心差分格式離散,動量方程的對流項采用三階守恒律的單調(diào)迎風(fēng)格式,湍動能、湍流耗散率和能量方程對流項采用二階迎風(fēng)格式離散,非定常項采用一階隱式格式離散,水氣體積分?jǐn)?shù)離散采用Modified HRIC 格式以精確捕捉水氣界面。采用全流場運動網(wǎng)格模擬模型和周圍的流場運動。在空氣和水面兩相流的界面定義中運用了VOF模型,即同一個網(wǎng)格內(nèi)空氣和水的體積分?jǐn)?shù)之和為1,對于第1 相流空氣的體積分?jǐn)?shù)v1,如果一個網(wǎng)格中的v1= 0 表示不存在空氣,v1= 1 表示該網(wǎng)格中只有空氣,0< v1<1 則該網(wǎng)格中同時存在空氣和水。使用用戶自定義函數(shù)定義水面位置和流場初始壓力分布。同時使用六自由度模型求解平動和轉(zhuǎn)動方程確定模型的位置和姿態(tài)。當(dāng)模型為對稱運動時,可將側(cè)滑和滾轉(zhuǎn)方向的運動約束住,從而簡化為三自由度模型,方程為

式中,Vi表示x或y方向的速度,Fi表示x或y方向的合力,t表示時刻,m表示質(zhì)量,M表示俯仰力矩,I表示俯仰慣量,θ表示俯仰角,n和n ?1 分別表示當(dāng)前和上一個時間步長。由此可獲得模型當(dāng)前的位置和姿態(tài)為

2 二維垂直入水算例

2.1 算例參數(shù)

首先對兩個典型的簡單二維模型垂直入水進(jìn)行模擬,模型參數(shù)和試驗數(shù)據(jù)均取自試驗[18-19]。算例中,圓柱直徑200 mm,高度200 mm,重12.5 kg,入水速度為0.989 m/s,在圓心角0?,7.5?,15?和30?的位置布置了壓力傳感器；楔形體長333 mm,截面寬度80 mm,高度45 mm,楔角20?,重1.639 kg,入水速度為1.28 m/s,采用高速攝像設(shè)備獲取運動過程中的位移和速度參數(shù)。圖1 所示為試驗相關(guān)設(shè)備及模型參數(shù)。

圖1 二維算例試驗示意圖(續(xù))Fig.1 Experiment diagram of the two-dimensional example(continued)

圖1 二維算例試驗示意圖Fig.1 Experiment diagram of the two-dimensional example

2.2 求解過程

2.2.1 流場設(shè)置

對計算域進(jìn)行網(wǎng)格劃分,由于模型均左右對稱,因此計算時采用半模,對模型附近尤其是底部與水面交界處的網(wǎng)格進(jìn)行適當(dāng)加密。對稱面采用對稱邊界條件,模型為無滑移壁面條件,其余邊界均為壓力入口邊界條件。計算域的網(wǎng)格劃分結(jié)果如圖2 所示。

圖2 二維算例流場網(wǎng)格Fig.2 Mesh of the two-dimensional example

2.2.2 時間步長設(shè)置

采用第1 節(jié)中的方法進(jìn)行求解參數(shù)設(shè)置。由于入水過程為動態(tài)變化過程,時間步長的設(shè)置需要兼顧對運動特征的捕捉、計算穩(wěn)定性以及計算效率。為了保證計算穩(wěn)定性,需在每一個時間步長內(nèi)迭代計算達(dá)到收斂,然后開始下一個時間步長。因此時間步長應(yīng)選定一個合理的區(qū)間,太大或太小都可能會導(dǎo)致無法收斂。在確定數(shù)值時一般會參考特征長度除以特征速度得到的特征時間,在此基礎(chǔ)上低兩個數(shù)量級。另外由于入水沖擊過程中速度變化較為劇烈,因此可將時間步長設(shè)置為可變的,設(shè)定最大和最小時間步長以及步長變化系數(shù),計算過程中根據(jù)收斂情況自動更改時間步長。以二維圓形算例為對象,特征長度為0.2 m,特征速度為0.989 m/s,計算特征時間約為0.2 s,將初始時間步長設(shè)置為10?4s,分別對最大時間步長10?2s,10?3s,和10?4s 進(jìn)行影響研究。圖3 所示為不同時間步長計算得到的圓弧角0?(即最低點) 壓力系數(shù)Cp與試驗結(jié)果的對比,可以看出,不同時間步長下,對于壓力峰值和變化趨勢的捕捉結(jié)果基本一致,在入水后期由于速度大幅降低,壓力數(shù)值較小,計算結(jié)果出現(xiàn)小幅波動,但此時并不是沖擊問題的關(guān)注點,影響可以忽略。綜上,將時間步長設(shè)置為可變的可以降低計算結(jié)果對時間步長的依賴性,同時保證計算精度?？紤]到計算效率,可以適當(dāng)增加最大時間步長,因此在后續(xù)的計算中,采用初始時間步長10?4s,時間步長范圍設(shè)定為10?3s～10?6s。

圖3 不同時間步長壓力計算結(jié)果對比Fig.3 Pressure results of different time step

2.3 仿真結(jié)果

圖4顯示了圓形入水過程中圓弧角0?(即最低點) 與7.5?位置壓力數(shù)據(jù)的計算結(jié)果與試驗結(jié)果對比,圖5 顯示了楔形入水過程中位移、速度的模擬結(jié)果和試驗結(jié)果對比,可以看出該仿真方法對于壓力峰值、壓力變化趨勢、垂向位移和垂向速度均能夠達(dá)到較高的準(zhǔn)確度。

圖4 圓形入水壓力結(jié)果對比Fig.4 Pressure comparison of circle example

圖5 楔形入水結(jié)果對比Fig.5 Result comparison of wedge example

在工程實際應(yīng)用中,不僅關(guān)心飛機(jī)在水上迫降過程中運動和受力情況,同時也需關(guān)注水面變形情況,為確定水線位置以及后續(xù)的漂浮特性提供依據(jù)。因此能否準(zhǔn)確模擬水面位置和動態(tài)變化趨勢也十分重要。圖6 所示為仿真與試驗中模型入水過程中水面變化情況的對比,水面的變形過程得到了很好的復(fù)現(xiàn)。雖然無法捕捉到液體的破碎和噴濺,但這些現(xiàn)象對于宏觀運動趨勢影響較小,在更關(guān)注運動姿態(tài)和受力形式的飛機(jī)水上迫降問題研究中可以忽略。

圖6 水面變形對比Fig.6 Water surface deformation comparison

3 三維水上迫降算例

3.1 算例參數(shù)

仿真方法經(jīng)過簡單外形的驗證之后,為了考慮更接近真實情況的復(fù)雜外形和三維效應(yīng),對水上迫降典型算例NACA TN 2929 模型試驗[3]進(jìn)行了仿真模擬。該經(jīng)典試驗中研究了多種機(jī)身外形對迫降特性的影響。選取其中更接近民用飛機(jī)后機(jī)身外形的Model F 進(jìn)行建模,如圖7 所示,該模型長1.2 m,翼展1.68 m,襟翼偏角60?,重5.67 kg,飛行速度9.144 m/s,初始俯仰角10?。

圖7 三維算例模型Fig.7 Model of three-dimensional example

3.2 求解過程

由于試驗時飛機(jī)為對稱運動,建模時采用半模,機(jī)體表面網(wǎng)格如圖8 所示,對稱面采用對稱邊界條件,機(jī)身表面采用無滑移邊界條件,其余邊界均為壓力入口邊界條件。所有網(wǎng)格均為六面體結(jié)構(gòu)網(wǎng)格,網(wǎng)格總數(shù)約240 萬,計算初始時間步長為10?4s。

圖8 三維算例網(wǎng)格Fig.8 Mesh of the three-dimensional example

3.3 計算結(jié)果

圖9所示為模擬與試驗的對比結(jié)果,其中實線表示試驗數(shù)據(jù),虛線表示仿真結(jié)果。可以看出俯仰角和重心高度的變化趨勢總體上非常吻合。飛機(jī)開始觸水即0～t1時刻,重心高度隨著入水深度的增加逐漸降低,之后由于俯仰角逐漸增大,而重心位置靠近前機(jī)身,因此高度也隨之增加。t2時刻開始,俯仰角逐漸減小即飛機(jī)逐漸低頭,飛機(jī)在氣動升力和水面浮力作用下保持在水面上滑行,入水深度變化較小,高度隨俯仰角的降低而減小。t3時刻開始前機(jī)身逐漸入水,受到?jīng)_擊力俯仰角略有增大,但變化幅度較小,同時飛機(jī)在浮力和重力作用下繞重心上下浮動,因此重心高度逐漸趨于穩(wěn)定。

圖9 三維算例結(jié)果對比Fig.9 Results comparison of the three-dimensional example

4 某型支線飛機(jī)水上迫降模擬

4.1 參數(shù)設(shè)置

經(jīng)過以上幾種算例的驗證,將該方法應(yīng)用于某型支線飛機(jī)的水上迫降過程,進(jìn)行仿真模擬。試驗相關(guān)參數(shù)和數(shù)據(jù)取自文獻(xiàn)[20]。參考縮比模型試驗中的參數(shù)設(shè)置,模型重37.88 kg,飛行速度為17 m/s,初始俯仰角為9?。同樣該飛機(jī)水上迫降過程可認(rèn)為是對稱運動,依然選擇半模進(jìn)行計算,模型采用全封閉的發(fā)動機(jī)短艙,襟縫翼為著陸構(gòu)型。流場網(wǎng)格劃分結(jié)果如圖10 所示,邊界條件設(shè)置與第3 節(jié)相同。

圖10 某型支線飛機(jī)流場網(wǎng)格Fig.10 Mesh of the regional aircraft

4.2 計算結(jié)果對比

圖11所示為計算得到的俯仰角和重心處垂向加速度隨時間的變化曲線與試驗對比結(jié)果。俯仰角變化趨勢與第3 節(jié)中的經(jīng)典模型類似,經(jīng)歷了兩次明顯的上仰過程。重心處垂向加速度在入水后約0.15 s出現(xiàn)峰值,隨后出現(xiàn)波動并逐漸趨于0?。俯仰角和加速度峰值的模擬數(shù)據(jù)與試驗值誤差分別為1.5%和4.5%,表明模擬結(jié)果能夠較好地反映水上迫降運動過程和模型整體受力情況。

圖11 支線飛機(jī)結(jié)果對比(續(xù))Fig.11 Results comparison of the regional aircraft (continued)

圖11 支線飛機(jī)結(jié)果對比Fig.11 Results comparison of the regional aircraft

4.3 運動過程和受力分析

4.3.1 運動過程

結(jié)合圖12 中的水線位置和機(jī)體表面壓力變化過程,可見0.2 s～0.6 s 之間,后機(jī)身和發(fā)動機(jī)短艙底部大部分區(qū)域均有明顯的負(fù)壓,產(chǎn)生較強(qiáng)的后體吸力,這是使得俯仰角增大的主要原因。從圖11(b)加速度曲線可知此時飛機(jī)的加速度到達(dá)峰值,雖然合力仍表現(xiàn)為向上的沖擊力,但由于負(fù)壓區(qū)更靠近尾部,力臂更長,因此合力矩表現(xiàn)為抬頭力矩。

圖12 迫降過程中的水線位置和壓力分布Fig.12 Water surface position and pressure distribution during ditching

0.6 s 左右俯仰角達(dá)到第一個峰值約32?。由于該支線飛機(jī)采用T 型尾翼,高置平尾,因此即使俯仰角到達(dá)峰值時平尾仍然位于水面以上,局部噴濺的水花對平尾的沖擊力較小,無法起到抑制俯仰角繼續(xù)上仰的作用。0.6 s 以后,由于巨大的水阻力使得滑行速度迅速降低、后體吸力減弱,同時中后機(jī)身以及發(fā)動機(jī)短艙著水部位受到較大的沖擊力,產(chǎn)生低頭力矩,從而使得飛機(jī)逐漸低頭。過程中飛機(jī)持續(xù)下沉,中后機(jī)身和發(fā)動機(jī)不斷排開水面,沖擊力出現(xiàn)小幅增大。

0.8 s～1.0 s 之間,隨著飛機(jī)的低頭前機(jī)身逐步入水,在此過程中受到較大的沖擊力產(chǎn)生抬頭力矩,使得1.0 s 以后加速度和俯仰角均出現(xiàn)第二次峰值,但由于此時速度已經(jīng)大幅度降低,第二次峰值相較第一次有明顯降低。此后飛機(jī)逐漸趨于穩(wěn)定,1.2 s時機(jī)身底部的壓力分布已經(jīng)顯著降低且分布較為均勻。

4.3.2 受力分析

圖13(a) 給出了整個機(jī)身以及機(jī)身尾部受到的法向力隨時間的變化曲線,t= 0 s 時刻機(jī)體開始著水,可以看到法向合力的變化趨勢與圖11(b) 的過載曲線一致,并且機(jī)身尾部受力始終為負(fù)且峰值達(dá)到約?1100 N,后體吸力作用明顯。在著水1.2 s 以后受力逐漸趨于穩(wěn)定,表明飛機(jī)逐漸停止上下沉浮運動。圖13(b) 所示為飛機(jī)的俯仰力矩隨時間的變化曲線,使飛機(jī)抬頭為正。入水初期俯仰力矩迅速增加,飛機(jī)出現(xiàn)第一次抬頭,隨后俯仰力矩產(chǎn)生波動并逐漸減小,對應(yīng)圖11(a) 中抬頭趨勢逐漸平緩,0.6 s以后低頭力矩逐漸減小,對應(yīng)圖11(a)中低頭趨勢也逐漸平緩。0.8 s 以后,與運動過程分析一致,俯仰力矩再次增加從而產(chǎn)生第二次抬頭運動,1.2 s 后俯仰力矩趨于0 表明飛機(jī)逐漸停止俯仰運動,達(dá)到靜止?fàn)顟B(tài)。

圖13 迫降過程中的受力隨時間變化Fig.13 Time history of force during ditching

5 結(jié)論

通過對圓形和楔形垂直入水以及三維飛機(jī)入水試驗的模擬,基于全流場運動網(wǎng)格的CFD 水上迫降仿真方法得到了驗證。利用該方法研究了某型支線飛機(jī)水上迫降的運動和受力特性。結(jié)果表明仿真得到的俯仰角和加速度數(shù)據(jù)均與試驗值相近,對后體吸力的準(zhǔn)確模擬證實了其對于迫降運動狀態(tài)和受力特點的重要影響。通過對仿真結(jié)果的進(jìn)一步分析表明,尾吊?高平尾布局的支線飛機(jī)在水上迫降時,后機(jī)身和發(fā)動機(jī)短艙底部的后體吸力會導(dǎo)致飛機(jī)產(chǎn)生明顯的抬頭,由于平尾高置無法抑制抬頭,俯仰角峰值會達(dá)到30?以上,入水后期水阻力使得速度減小,從而減弱吸力,使得飛機(jī)逐漸低頭并在較小的波動后逐漸趨于平穩(wěn)。針對民用飛機(jī)水上迫降特性的研究目的,本文目前僅完成了方法的驗證和典型迫降工況的分析,對于俯仰角、氣動力等不同因素的影響研究還需持續(xù)開展。