水上飛機起飛歷程的水動力特性研究

段旭鵬，孫衛(wèi)平，魏猛，左仔濱，楊永

(1.西北工業(yè)大學(xué) 航空學(xué)院，西安 710072)(2.中航通飛研究院有限公司 總體部，珠海 519040)(3.中國特種飛行器研究所 水動力研究中心，荊門 448035)

0 引 言

隨著海洋經(jīng)濟的發(fā)展，人類對海洋資源的開發(fā)越來越深入，快速實施對遠海島嶼的補給和海上事故的救援顯得越來越重要，水上飛機提供了一種高效便捷的解決方案。為了滿足在水面滑行、起飛和降落的需要，水上飛機必須具備許多不同于陸基飛機的設(shè)計特點。最為顯著的就是機身底部與陸基飛機截然不同，它背負著整架飛機在水面上高速滑行，需要具備優(yōu)良的水面滑行特性。

水上飛機的一個重要研制工作就是船體設(shè)計。因為飛機要從水上起飛，船體的設(shè)計必須具備較高的滑行效率，否則會影響飛機的起飛性能。事實上，水阻力在起飛過程的大部分時間里都比氣動阻力大很多。在發(fā)動機的拉力選擇上，水上飛機的最小拉力通常取決于水上起飛阻力的需求。因此，為了使飛機的總體設(shè)計更為合理，起飛過程水動力特性的準(zhǔn)確分析至關(guān)重要，否則飛機會背負來自發(fā)動機的額外不必要的重量，大大削弱飛機巡航性能。

在低速滑行階段飛機受到的氣動力很小，對飛機的影響可以忽略不計，大部分的重量由水面承載，因此飛機的阻力主要來自水動阻力。隨著飛行器的加速，水動阻力逐漸增大，在某一時刻達到最大值，稱為阻力峰(hump drag)。這一般發(fā)生在起飛速度的25%～50%[1]。一旦通過這個區(qū)域，氣動力開始主導(dǎo)飛機的運行，對水動升力的要求逐漸降低，水動阻力開始逐漸減小，直到起飛離水。

為了使水上飛機能夠順利起飛，船體需要將水動阻力控制得盡可能小，尤其是在高速滑行階段。為此水上飛機與滑行艇類似，采用了一個特別的設(shè)計，叫做斷階[2]。水面滑行(planing)的方式可以幫助提供更多的水動升力，但是因為滑行船體趨向于降低自身的縱傾角(即迎角)[3]，導(dǎo)致增大了浸濕面積進而加大了粘性阻力。為了解決這一問題，船海領(lǐng)域的研究者巧妙地設(shè)計了“斷階”，通過在幾何外形上的階躍使水在船體繞流過程產(chǎn)生了流動分離，造成了浸濕面積的顯著減小從而解決了這一問題。但是，由于在高速滑行前斷階處的流動變得更為復(fù)雜，水上飛機設(shè)計者們必須清楚地了解斷階處的滑行性能。水上飛機船體滑行性能的評估以往通常是由水池拖曳試驗得到的。

兩相流的數(shù)值研究雖然比較多，但是針對水上飛機的研究還不多見。Qiu L J等[4]利用解耦式算法，基于商業(yè)軟件開展了水陸兩棲飛機起飛過程的動態(tài)模擬。該方法將氣動力和水動力分開計算，有較高的計算效率，但是尚未開展流動現(xiàn)象的研究。Shen Z等[5]對OpenFOAM添加了重疊網(wǎng)格模塊，并對KCS標(biāo)模船體的低速自主航行進行了計算模擬，計算采用VOF(volume of fluid)、動態(tài)重疊網(wǎng)格以及6DOF(6 degree of freedom)等技術(shù)模擬了船體在螺旋槳旋轉(zhuǎn)、舵面偏轉(zhuǎn)條件下非定常動態(tài)氣水兩相流問題，得到了較好的結(jié)果，然而基于OpenFOAM平臺的改進方法對航空級高速問題的計算效率和精度如何尚無結(jié)論。飛機水上迫降與此類似，Qu Q L等[6]計算了支線客機的水上迫降問題，采用VOF、6DOF以及全局動網(wǎng)格技術(shù)研究了上單翼、高平尾飛機的水上迫降性能。

由于國外對水上飛機的需求不強，最近幾十年一直沒有新的機型問世，相應(yīng)的技術(shù)攻關(guān)停滯已久。特別是水上飛機的水動力特性方法研究還停留在第二次世界大戰(zhàn)期間的水池拖曳試驗方法上。

在20世紀(jì)興起的新興學(xué)科——CFD領(lǐng)域幾乎沒有相關(guān)問題的研究。隨著算法理論的成熟和計算機硬件的發(fā)展，數(shù)值研究飛行器水面高速航行問題已成為可能。本文采用水池拖曳試驗和數(shù)值計算兩種研究手段，研究飛機滑行中的平衡態(tài)阻力、升沉和俯仰三者的變化歷程，重點分析斷階滑行階段的氣-水兩相流力學(xué)特征，證明現(xiàn)代數(shù)值計算方法應(yīng)用于水上飛機研制的可行性，以期為航空級兩相流高速CFD數(shù)值計算方法體系的建立提供技術(shù)支撐。

1 試驗和計算模型

設(shè)計加工一個某型水上飛機的縮比模型，在高速水動力航空科技重點實驗室中進行實驗，如圖1所示。模型外形如圖2所示，主要尺寸如表1所示。水上飛機為適應(yīng)水上起降的需求，采用高置上單翼，4發(fā)渦槳發(fā)動機，采用T尾布局以避免水面噴濺對升降舵的沖擊。翼下兩側(cè)安裝有浮筒，襟翼偏至起飛位置。船體被斷階分成前體和后體兩部分，前體周圍安裝有抑波槽和壓浪板，用以抑制噴濺對高速旋轉(zhuǎn)螺旋槳的沖擊。試驗和計算均為無動力模擬，因此3D數(shù)模中未設(shè)計螺旋槳。

圖1 水池拖曳試驗中模型靜止時刻的照片

圖2 水上飛機3D模型

參 數(shù)數(shù)值參 數(shù)數(shù)值全機長/m3.75機翼翼展/m3.9船體寬/m0.35機翼面積/m21.7排水量/kg49

2 自由拖曳水池試驗方法

水上飛機水池拖曳試驗的方法與文獻[7]中描述的方法類似。不同之處在于，一般的水池拖曳試驗只使用船體模型，而本次試驗使用的是帶氣動布局的全機模型，盡管如此，試驗中仍要進行升力和力矩的卸載。其原因是在水池拖曳試驗過程中只能保證水的動力學(xué)相似，而空氣的動力學(xué)相似是無法做到的。Jr J.D.Anderson[8]指出兩種流動要達到動力學(xué)相似必須滿足①流線的幾何相似；②無量綱的速度、密度、溫度等流場特征物理量在全流域分布相同；③力系數(shù)相同。要滿足上述要求，需要做到：①模型幾何相似；②相似參數(shù)相同。理想的模型試驗要使得主導(dǎo)水動力相似的Fr和主導(dǎo)氣動力相似的Re、Ma對于模型和實機運動而言滿足同時相等，然而這在水池中是無法做到的，因為水池拖曳試驗僅能保證Fr相等。

三個相似參數(shù)定義如下：

(1)

(2)

(3)

式中：ρ為密度;μ為動力學(xué)粘性系數(shù);U為來流速度;L為參考長度。

為了解決這一問題，對于某一拖曳速度，在水池拖曳試驗開始前先將模型提離水面，利用拖車將模型在水面上方拖行，可以測得機翼和平尾提供的模型氣動升力和氣動力矩，在此基礎(chǔ)上利用風(fēng)洞試驗和數(shù)值計算得到的實機飛行氣動力系數(shù)CL和CM，進行相減分別得到升力和力矩差量，以此作為卸載即可補齊由于空氣流動不相似帶來的氣動力差異。從而達到了“氣動力相似”的效果，測得的無量綱水動阻力、縱傾角度和無量綱升沉即可直接應(yīng)用到真實飛機上。

3 數(shù)理模型和計算方法

采用OpenFOAM平臺作為計算流體力學(xué)模擬工具對相同尺度的模型進行數(shù)值計算。因為本文是對數(shù)值方法進行驗證，因此計算狀態(tài)采用與水池拖曳試驗完全相同的條件，基于URANS的兩相流解算器開展計算研究。

3.1 兩相流控制方程

VOF方法是一種自由面捕捉方法，其思想較為簡單：每一個計算網(wǎng)格單元中都額外包含一項代表水的體積分數(shù)α，由初始氣水交界面幾何位置確定初值。位于交界面的網(wǎng)格單元既含有氣體又含有水，因此0<α<1；位于交界面之上的網(wǎng)格單元都是氣體，因此α=0；位于交界面以下的網(wǎng)格單元都是水，因此α=1。在OpenFOAM中，牛頓流體不可壓守恒形式的雷諾平均N-S方程組如式(4)和式(5)所示，其中式(4)是質(zhì)量守恒方程，式(5)是動量方程。

(4)

(5)

式中：U為速度矢量；ρ為控制體內(nèi)的流體密度；p為壓力；μ為動力粘性系數(shù)；μt為湍流動力粘性系數(shù)。

VOF方法是在NS方程的基礎(chǔ)上求解體積分數(shù)的，α隨著速度場U運動的輸運方程如方程(6)所示。

(6)

其中α定義為Vwater/Vtotal，數(shù)值介于0和1之間，因此流體空間內(nèi)的密度和粘性可由式(7)和式(8)表達。

ρ=ρair(1-α)+ρwaterα

(7)

μ=μair(1-α)+μwaterα

(8)

α函數(shù)在一個無限薄的界面上從1變成0，因而難以對α進行梯度的估計，進而導(dǎo)致氣水交界面的數(shù)值耗散較大。為了改善這個問題，OpenFOAM引入一個虛擬壓縮項，即加入一個虛擬速度場W，并使得W垂直于界面，保證最終結(jié)果與原始方程保持一致。因此α輸運方程的最終形式為

(9)

兩相流的湍流模型方程形式與單項流完全相同，密度和粘性需要用式(7)和式(8)表示。

3.2 剛體外力模塊

3.3 數(shù)值求解方法

采用不可壓等溫不混溶的非定常兩相流解算器，附帶可選擇的網(wǎng)格運動和網(wǎng)格拓撲重構(gòu)模塊并耦合6DOF方程的求解。對于本文研究的問題不對網(wǎng)格細化重構(gòu)，只執(zhí)行網(wǎng)格運動變形。解算器采用VOF方法計算每個控制體單元內(nèi)的氣、水所占組分，用以捕捉氣水交界面即自由面。

用半隱式MULES(Multidimensional Universal Limiter with Explicit Solution)方法[9]求解含有體積分數(shù)α的方程(6)。在非定常計算中盡量控制時間步長，因為時間步長越短，壓力與速度的耦合越強，這是PISO(Pressure-Implicit with Splitting of Operators)算法[10]的典型特征。該算法具有較高的效率和精度，但是對網(wǎng)格質(zhì)量較為敏感。壓力與速度的耦合是通過傳統(tǒng)的多重網(wǎng)格處理的，同時大量采用Rhie-Chow[11-12]插值方法，提高了PISO算法對各類網(wǎng)格的適應(yīng)性。

時間導(dǎo)數(shù)項采用歐拉離散格式進行離散，傳導(dǎo)項和耗散項的體積分通過高斯公式轉(zhuǎn)換成面積分數(shù)進行計算。在求解非定常自由面時，把動量方程中的傳導(dǎo)項線性離散并加入限制器，體積分數(shù)方程的對流項采用vanLeer格式進行計算。

具體計算步驟如下：

步驟一：根據(jù)科朗數(shù)計算時間步長；

步驟二：求解網(wǎng)格運動方程，對網(wǎng)格進行變形；

步驟三：更新兩相流運動參數(shù)，包括密度和粘性系數(shù)；

步驟四：采用MULSE方法求解方程(6)，得到體積分數(shù)；

秦立紅等[47]通過大鼠腎上腺嗜鉻細胞瘤克隆化細胞株(PC12)實驗，發(fā)現(xiàn)從仙草中分離出的咖啡酸、3-(4-乙氧基-3-羥基-苯基) 烯丙酸、咖啡酸乙酯、山奈酚、高山黃芩素、2-十六烷基-十八烷酸、熊果酸和豆甾醇等8個單體化合物有較高的抗缺氧活性。

步驟五：采用PISO算法求解速度和壓力。

采用欠松弛迭代方法[13]來增加數(shù)值求解的穩(wěn)定性以提高收斂性。對任意待求解變量，將新迭代得到的數(shù)值與上一步數(shù)值進行加權(quán)平均。其數(shù)學(xué)表達如下，令β為迭代因子，有

φnew=βφnew+(1-β)φold

(10)

式中：φnew為新計算得到的結(jié)果；φold是上一步迭代的結(jié)果。

對于本文選擇的k-ω湍流模型，其中的ω方程可以對整個壁面邊界層進行積分，能更為方便地處理近壁面流場解算問題。解算器提供的壁面函數(shù)可以自動切換低雷諾數(shù)和高雷諾數(shù)求解方程，使得ω的解可以在粘性附面層和對數(shù)近壁面區(qū)域進行靈活適配[14]。

4 計算網(wǎng)格和邊界條件

計算網(wǎng)格采用對水面適應(yīng)性較好的笛卡爾網(wǎng)格。雖然結(jié)構(gòu)網(wǎng)格具有流向延展性的特點，但是受物面網(wǎng)格拓撲的影響，在機身近水面附近的網(wǎng)格延伸方向通常會被迫改變，導(dǎo)致自由面計算的誤差增加。然而，笛卡爾網(wǎng)格一方面可以嚴(yán)格按照水平方向進行各向異性的網(wǎng)格加密，另一方面在水面與物面相交處，它具有非結(jié)構(gòu)網(wǎng)格的特點，可以靈活對物面進行貼合，而不改變自由面的網(wǎng)格拓撲方向。

水上起飛水動力占據(jù)合力的大部分比重，水動力的模擬至關(guān)重要，因此除了機翼前緣等部位還特別對船體表面和周圍近場空間進行加密，如圖3所示。對于本算例，相同首層高度網(wǎng)格對應(yīng)的水動力與氣動力y+恰好相差10倍，水下y+取100是可行的。因此全機物面邊界層首層高度統(tǒng)一設(shè)為0.1 mm，對應(yīng)的氣動力y+為10，水動力y+為100，總網(wǎng)格量為900萬。

圖3 計算網(wǎng)格示意圖

在OpenFOAM中，邊界條件是施加在計算域的patch上的，本文的計算均采用圖4中的邊界條件設(shè)置。藍色長方體代表水存在的區(qū)域，其上為空氣介質(zhì)，整個空間為空氣-水的兩相流控制體。其中入口邊界距離機頭為2倍體長，出口邊界距離機尾5倍體長，計算域?qū)?倍體長，水深3倍體長，空氣流域高度2倍體長。計算過程中來流速度不變，飛機通過動態(tài)網(wǎng)格和6DOF模塊進行自由升沉和俯仰的2自由度運動，與試驗類似地，當(dāng)運動達到穩(wěn)定即認為得到一個自由拖曳的數(shù)值結(jié)果。

圖4 邊界條件示意圖

5 計算結(jié)果分析

飛機水上起飛過程可以細分為四個階段[17]：排水航行階段、過渡階段、滑行階段和起飛階段。飛機啟動后，受螺旋槳拉力作用開始低速航行，飛機重量主要由浮力支撐，與排水型船舶類似，此階段稱之為排水航行階段。在速度達到起飛速度的25%～50%，飛機進入過渡階段，從前體斷階處至后體前部之間區(qū)域開始出現(xiàn)一個“空氣墊”，導(dǎo)致后體滑行面和水動升力減小，低頭力矩隨之減小，縱傾角增加，全機水阻力迅速增加直至達到一個峰值。斷階滑行階段中，速度范圍位于起飛速度的50%～80%，隨著速度的增大，斷階后空氣墊增長，后體滑行面積進一步減小，直至后體離水，形成單斷階滑行，此時在船底上的主要受力為前體水動力，并且氣動力開始起作用，氣動升力逐漸占據(jù)主導(dǎo)，平尾力矩和升降舵控制力開始逐漸發(fā)揮效用。當(dāng)飛機速度大于起飛速度的80%時，飛機進入起飛階段，該階段飛機斷階只有最后一小部分三角形滑行區(qū)位于水面之上，氣動力幾乎完全占據(jù)主導(dǎo)，氣動升力顯著，滑流和地面效應(yīng)共同作用，飛行員可自由調(diào)整升降舵使飛機處于合適的迎角，隨著速度的不斷增加飛機最終離水完成水上起飛。

5.1 水動力特性和計算驗證

數(shù)值計算和水池拖曳試驗的結(jié)果如圖5～圖7所示。自由俯仰和升沉運動的阻力結(jié)果如圖5所示，其中阻力(D)用模型重量(W)進行無量綱化；縱傾角變化如圖6所示；質(zhì)心升沉隨速度變化如圖7所示，升沉(Rise)用船寬(b)進行無量綱化；速度統(tǒng)一采用速度系數(shù)表達，定義如下：

(11)

式中：g為重力加速度；b為船寬。

模型起飛速度大約在CV=8左右?？傮w而言，計算與試驗吻合較好，特別是在第二、三、四階段，在過渡區(qū)、高速滑行、起飛區(qū)對于復(fù)雜兩相流耦合飛行力學(xué)方程的動態(tài)非定常計算而言，全程阻力計算精度可以控制在10%以內(nèi)。計算反映出了過渡區(qū)復(fù)雜的力學(xué)積分結(jié)果以及起飛前的質(zhì)心離水歷程，可以借助數(shù)值模擬結(jié)果揭示試驗中難以捕捉的流動現(xiàn)象和力學(xué)規(guī)律。

CV在0～2區(qū)間是水上起飛的第一階段，力學(xué)現(xiàn)象較為簡單，阻力與速度基本成線性關(guān)系，升沉俯仰無明顯變化，根據(jù)以往經(jīng)驗，此段無需做過多研究，因此試驗和計算只研究一個速度狀態(tài)。觀察CV在2～4區(qū)間，阻力、俯仰、升沉三條曲線都顯示了該區(qū)域的飛機的動力學(xué)特性變化顯著。阻力先是快速增加隨后斜率放緩最后在CV=3.7達到阻力峰，然后逐漸降低。一旦通過了該階段，阻力開始迅速降低，飛機進入斷階滑行和起飛階段，CV=4.5到CV=6.2阻力變化開始放緩。在此階段隨著速度的增加，氣動力逐漸占據(jù)主導(dǎo)，尤其是氣動升力的增加為水動載荷分擔(dān)了越來越多的重量導(dǎo)致水阻力大幅降低。

圖5 阻力隨速度關(guān)系對比

圖6 縱傾角隨速度關(guān)系對比

圖7 質(zhì)心升沉隨速度關(guān)系對比

5.2 阻力峰產(chǎn)生機制分析

從圖5可以看出：CV在2～4區(qū)間，即當(dāng)速度位于起飛速度的25%～50%飛機進入過渡階段。在此階段，在斷階后方開始出現(xiàn)一個“空氣墊”，導(dǎo)致后體滑行面和水動升力減小，抬頭力矩隨之增大，使得飛機縱傾角顯著增加。隨著速度的增加縱傾角減小，阻力峰出現(xiàn)。

過渡階段滑行過程中的船體浸濕表面發(fā)展歷程如圖8所示。后體的空氣墊(船底藍色區(qū)域)在CV=2.8已在斷階之后大面積形成。CV=3.4時氣墊區(qū)域開始向后延伸，同時與水流進行摻混形成復(fù)雜的水氣交混區(qū)域。如果隱藏掉浮筒，從側(cè)視圖可以看到起落架艙和后機身側(cè)方開始受到噴濺作用。當(dāng)速度達到CV=3.75，雖然后體的浸濕區(qū)域變小，但是觀察側(cè)視圖可知由于前體舭彎處的壓力增加使得水流順抑波槽向后強烈噴出并直接大面積覆蓋船體的側(cè)面突出的起落架艙。

(a1) 底部視圖 (a2) 側(cè)視圖

(a)CV=2.8

(b1) 底部視圖 (b2) 側(cè)視圖

(b)CV=3.4

(c1) 底部視圖 (c2) 側(cè)視圖

(c)CV=3.75

圖8 過渡階段滑行過程中的船體浸濕表面發(fā)展歷程

Fig.8 Development process of hull wetted surface during taxiing in transitional phase

三個速度下的船體壓力分布如圖9所示。CV=2.8和CV=3.4的壓力分布差異不大，因而全機阻力差異也很小，這與圖5中的壓力曲線是相一致的。不同之處在于后者的起落架艙后部開始受到一定的壓力作用，表明噴濺現(xiàn)象開始出現(xiàn)。當(dāng)速度增加到CV=3.75時，舭彎前部、起落架艙前部的紅色區(qū)域表明壓力增大，前體噴濺基線已經(jīng)從壓浪板內(nèi)側(cè)后移至舭彎處，壓力中心向質(zhì)心靠攏，導(dǎo)致縱傾角減小。前體舭彎處的壓力增加使得水流順抑波槽向后強烈噴出，并沖刷船體側(cè)面和突出的起落架艙，形成紅色的壓力駐點。因此當(dāng)CV=3.75時，盡管后體底面浸濕面積減小，抑波槽水流的噴濺和沖刷導(dǎo)致了阻力峰的最終形成。

(a1) 底部視圖 (a2) 側(cè)視圖

(a)CV=2.8

(b1) 底部視圖 (b2) 側(cè)視圖

(b)CV=3.4

(c1) 底部視圖 (c2) 側(cè)視圖

(c)CV=3.75

圖9 過渡階段滑行過程中的船體壓力分布發(fā)展歷程

Fig.9 Development process of hull pressure distribution during taxiing in transitional phase

5.3 滑行階段中初值對計算結(jié)果的影響

滑行階段初期用不同初值作為輸入計算得到的兩個數(shù)值結(jié)果如圖10所示。

(a) 深吃水作為初值

(b) 淺吃水作為初值

從圖10(a)可以看出：吃水較深噴濺較強不符合物理現(xiàn)象。從圖10(b)可以看出：后體離水清晰，斷階滑行明顯，符合物理規(guī)律。在正常的起飛過程中，當(dāng)飛機通過阻力峰以后，其質(zhì)心應(yīng)該繼續(xù)向上攀升，后體逐漸清晰離水，飛機進入斷階滑行階段。然而在數(shù)值研究中發(fā)現(xiàn)，如果斷階滑行初期如CV=4.5，初始自由面高度設(shè)置過高，如圖10(a)所示，即飛機從吃水較深的狀態(tài)釋放，就可能出現(xiàn)機身被“吸住”無法正常上升的情況，導(dǎo)致出現(xiàn)阻力和航態(tài)與試驗結(jié)果差別較大的情況。然而如果初始將飛機提起一定高度，如圖10(b)所示，使斷階后方有機會順利進氣，那么計算結(jié)果將恢復(fù)與試驗較好的一致性。

其原因正是由于此階段的滑行特點決定的。斷階減阻的原理是通過引入空氣到斷階后方，使液體流動產(chǎn)生分離并逐步抬升后體，導(dǎo)致浸濕面積大幅降低進而達到減阻的目的?？梢娤驍嚯A底部引入空氣是飛機斷階滑行的必要條件。數(shù)值上的靜水滑行可以認為是“鏡面”水面效應(yīng)沒有任何擾動。如果初始吃水較深，飛機斷階被水緊密包裹，那么計算中空氣無法進入到斷階后部，最終便形成了被“吸住”的效果；如果吃水較淺，使得斷階底部一開始便有空氣存在，斷階后的水氣兩相流作用使得數(shù)值上增加了足夠的擾動，那么最終結(jié)果就會符合客觀物理現(xiàn)象(圖例代表水面波浪高度，0代表水平面)。

6 結(jié) 論

(1) 利用CFD技術(shù)可以較為精確地模擬水池拖曳試驗，復(fù)現(xiàn)了水池拖曳試驗的流場形態(tài)和模型動力學(xué)特性，在起飛全過程的各個速度點上阻力誤差均在10%以內(nèi)。

(2) 對于本文研究的對象，縱傾角的極值并未直接導(dǎo)致阻力峰值的出現(xiàn)，機身側(cè)面起落架艙周圍的噴濺作用才是阻力峰出現(xiàn)的根本原因。

(3) 數(shù)值研究斷階滑行階段時需要避免“鏡面”效應(yīng)，給予斷階足夠的擾動才能保證飛機順利進入斷階滑行狀態(tài)。