趙偉鵬，張英浩，高 婷，梅 凱

（1.沈陽航空航天大學(xué)，遼寧 沈陽 110136；2.武漢第二船舶設(shè)計研究所，湖北 武漢 430205）

我國是一個海洋資源極其豐富的國家，據(jù)統(tǒng)計，截至2015年，全國海岸線總長為18 604.3 km，其中人工岸線總長為12 747.8 km，自然岸線總長5 856.5 km[1]。綿長的海岸線為屬地經(jīng)濟的發(fā)展和海洋資源的開發(fā)提供了更多機會，但同時海洋活動的越發(fā)頻繁和國家之間對海洋權(quán)益競爭的日趨激烈，都對海上偵查巡邏和應(yīng)急救援工作提出了更高的難度和挑戰(zhàn)。隨著飛機設(shè)計與建造技術(shù)的不斷進步，飛機的應(yīng)用需求越來越明確化，飛機的功能分類越來越多樣化，對設(shè)計和建造的要求也越來越精細。因此出現(xiàn)了電動飛機、地效飛機、水陸兩棲飛機和艦載機等一批專門針對水氣兩相界面應(yīng)用需求的特種飛機。這些飛機憑借著高效和快速的優(yōu)勢在海上巡航飛行、海上應(yīng)急救援和森林防火等方面扮演著相當(dāng)重要的角色。除了滿足國家需求外，隨著人們物質(zhì)生活水平的提高，飛機也作為一項重要工具而被廣泛應(yīng)用于體育運動、旅游觀光和海上運輸?shù)戎T多民用領(lǐng)域。但不管是面對日益繁雜的海上巡防任務(wù)，還是確保有關(guān)民用航空上下游產(chǎn)業(yè)的有序發(fā)展，都對飛機近海面飛行的安全提出了更高的要求，因此迫切需要探究飛機近水面飛行氣動特性。飛機近海面高速飛行過程，既受水面效應(yīng)影響[2-3]，又同時包含了氣體、液體、剛性固體三者的高速劇烈作用，因此總是處于一個非線性強、耦合特性強的復(fù)雜環(huán)境中。在近海面飛行過程中飛機狀態(tài)稍有偏差，強大的水動載荷極有可能對浮筒連接桿、機翼和發(fā)動機等重要構(gòu)件造成結(jié)構(gòu)性破壞，從而使飛機的漂浮特性進一步變差，且由于海上救援難度偏大，所以近海面飛行的事故率和事故嚴重程度都要遠遠高于常規(guī)飛行階段。作為近海面飛行的代表——艦載機，美國海軍發(fā)生的艦載機事故中有80%是發(fā)生在近海面的著艦過程中，BROOKS C J等[4-5]對1979—2006年發(fā)生在加拿大和1981—2011年發(fā)生在美國的直升機水面墜機事故進行了詳細的調(diào)查和分析后發(fā)現(xiàn)其死亡率分別高達23%和26%，因此很多學(xué)者對飛機的近海面飛行展開研究。米百剛等[6]基于非結(jié)構(gòu)動態(tài)網(wǎng)格重構(gòu)技術(shù)，對NACA0012二維翼型自由空間、近地效應(yīng)，以及近水面效應(yīng)的非定常運動流場進行了數(shù)值計算分析。韓龍等[7]使用滑動網(wǎng)格法對三維等直翼NACA0015在波浪條件下的水面效應(yīng)進行了數(shù)值模擬與分析。段旭鵬等[8]模擬過程全面考慮了飛機高速滑行時的水動力、氣動力、水面效應(yīng)和螺旋槳的動力影響。SCHWEIKHARD Y S等[9]采用離散渦模型計算了二維機翼的動態(tài)地面效應(yīng)；BARBER T J[10]提出水面上發(fā)生的任何變形很可能是由翼尖渦旋引起的，而不是機翼下方壓力分布增加引起的；MORYOSSEF Y等[11]在假設(shè)無粘和有粘流體的情況下，對地面效應(yīng)中的振蕩翼型進行了數(shù)值模擬。IM Y H等[12]使用基于LU因子分解算法和高階逆風(fēng)方案的歐拉代碼，對波浪表面上的二維翼型進行了數(shù)值模擬。

現(xiàn)有研究多以地效飛機或大型水陸兩棲飛機為研究對象，對輕型電動飛機的近海面飛行性能研究鮮少。然而，受設(shè)計布局、飛行性質(zhì)等因素的影響，飛機受到的自由液面影響也會有所不同。電動飛機具有體積小，載重量少，操作靈活、起飛和降落對水域要求較低等特點，相較于地效飛機和水陸兩棲飛機，它的機翼面積較小，體量較輕，因此難以測量其在近水面飛行時受到的水面效應(yīng)影響程度，且本身穩(wěn)定裕度小，所以升力和阻力更易受到外界條件的干擾。本文研究的電動飛機具有自升式滑翔飛行能力，滿足的機動能力包括：正常飛行的機動；緩“8”字，急上升，轉(zhuǎn)彎；坡度不大于60°的急轉(zhuǎn)彎和失速。目前以該型飛機為基礎(chǔ)已進一步設(shè)計研發(fā)了水上飛機[13]。因此，考慮該飛機的結(jié)構(gòu)特點，研究飛機的近海面飛行特性勢在必行，既可以為飛行員的訓(xùn)練提供必要的理論支撐，也可以為其他相關(guān)類型飛機的設(shè)計提供參考。

此外，現(xiàn)有研究視角多置于二維、三維機翼。飛機迎角不同、距水面高度不同和海況不同時，其受到的自由液面影響也有很大的不同，但從前人研究成果來看，還缺乏對這些方面的綜合對比。本文基于STAR-CCM＋軟件，采用歐拉多相流模型和VOF法，該方法是一種建立在歐拉網(wǎng)格下的界面捕捉方法。對某型電動飛機貼近海面飛行時的非定常動態(tài)運動流場進行數(shù)值模擬研究，對表征水面效應(yīng)的指標——升力系數(shù)、阻力系數(shù)進行可靠計算，分析不同迎角、不同離水高度和不同海況下水面效應(yīng)的影響，以期能對水上飛機設(shè)計和飛機近海面飛行性能分析等提供參考。

1 數(shù)值方法

1.1 控制方程和湍流模型

以雷諾平均方程作為求解電動飛機近海面飛行的基本方程，其具體形式如下。

式中，ρ為流體密度；μ為流體粘度；p為靜壓；fi為單位質(zhì)量的質(zhì)量力；ui和uj為速度分量。

湍流模型選取可實現(xiàn)的K-Epsilon兩層模型[14]，它是低雷諾數(shù)方法的替代方法，可以將K-Epsilon模型應(yīng)用于粘性影響層（包括粘性子層和緩沖層）。在壁面旁邊的層中，湍流耗散率ε和湍流粘度μ均被指定為壁面距離的函數(shù)。系統(tǒng)將對整個流體域的湍動能方程進行求解，下面給出該模型的數(shù)學(xué)表達式[15-17]。

式中，v代表速度；μ代表動態(tài)粘度；σk、σε、Ct、Cε1和Cε2都代表模型系數(shù)；Pk和pε分別代表兩種不同的K-Epsilon湍流模型變體式；Sk、Sε代表用戶自定義源項；ε0是源項中用于抵消湍流衰減的環(huán)境湍流值；f2代表相應(yīng)變體式的阻尼函數(shù)。

1.2 VOF多相模型及VOF波建模

本文選用兩相流液面模擬方法是界面捕捉法——VOF法，通過定義函數(shù)來區(qū)分不同的流體介質(zhì)，進而精準捕捉空氣—水交界面。界面捕捉法相較于界面追蹤法的優(yōu)勢在于易于實現(xiàn)、二次編程難度低和結(jié)果準確度高。VOF多相模型[18-19]是一種簡單的多相模型，適用于模擬數(shù)值網(wǎng)格（能夠求解混合物的相之間的交界面）上的多個不混溶流體流，此模型捕捉自由液面的思想如下：每一個計算網(wǎng)格單元中的水的體積分數(shù)為α。位于水、氣交界面處的網(wǎng)格單元既含有氣體又含有水，因此0＜α＜1；液面以上的網(wǎng)格單元都是氣體，因此α=0；位于交界面以下的網(wǎng)格都是水，因此α=1。

VOF波模型用于模擬輕流體（空氣）和重流體（水）之間的交界面上的表面重力波。海浪等級主要由有效波波高和1/10大波波高兩個參數(shù)來界定，參照我國對海浪等級的劃分標準，以及考慮到輕型電動飛機實際飛行情況，本研究主要模擬兩種海況：波高為0 m的平靜海面和波高為0.2 m的小浪海面。其中波浪水面使用斯托克斯波理論[20-21]對波浪建模，此波仿真精確度更高，能更接近海面真實波浪。其中波形和波相速度取決于水深、波高和水流，厄塞爾數(shù)定義如下。

式中，H為波高；λ為波長；d為水深。

2 氣動特性分析方法及其驗證

在對電動飛機近海面高速飛行進行數(shù)值模擬前，應(yīng)首先對本文所用的氣動特性分析方法進行驗證，因飛行過程并不觸水且受到實驗設(shè)備條件的限制，無法通過水池實驗對氣動特性進行準確計算。因此本文以某型電動飛機為例，分別模擬飛機以不同迎角在空氣中平飛時的情況，將升力系數(shù)、阻力系數(shù)和升阻比與風(fēng)洞實驗數(shù)據(jù)對比，驗證本文數(shù)值方法的準確性。

2.1 計算模型

某型電動飛機[22]采用大展弦比上單翼、低平尾、前置螺旋槳、前三點式不可收放起落架布局，幾何模型（半模）如圖1所示，飛機的設(shè)計參數(shù)如表1所示。

圖1 電動飛機模型

表1 某型電動飛機設(shè)計參數(shù)

2.2 計算域及網(wǎng)格劃分

為節(jié)約計算資源，縮短研究周期，且模型具有對稱性特點，本文僅對一半的幾何進行建模（圖2），計算域的范圍為6.5 L×1.5 L×5.0 L，其中L為機身長度，邊界條件設(shè)置如下：機頭前1.5倍機身長度，設(shè)置邊界條件為速度入口；機尾后4倍機身長度，設(shè)置邊界條件為壓力出口，壓力初始值為大氣壓力；計算域的上、下邊界設(shè)置為速度入口，計算域、左右邊界設(shè)置為對稱面；飛機表面定義為無滑移、不可穿透邊界條件，邊界條件設(shè)置為無滑移壁面。

圖2 計算域大小及邊界條件設(shè)置

結(jié)構(gòu)化網(wǎng)格是計算流體力學(xué)數(shù)值模擬中比較常見的網(wǎng)格劃分形式，具有簡單易生成、網(wǎng)格質(zhì)量高和占用計算資源少的優(yōu)點，它的缺點也顯而易見，就是在劃分復(fù)雜幾何模型的時變形較大，不能精確貼合復(fù)雜構(gòu)型，從而降低計算結(jié)果準確性。切割六面體網(wǎng)格[23-24]以兼具結(jié)構(gòu)化網(wǎng)格優(yōu)點同時又能保證對復(fù)雜表面精確貼合的獨特優(yōu)勢，在近幾年得到了廣泛發(fā)展。為既能減少生成的網(wǎng)格數(shù)量又能兼顧計算精度，本文使用切割體、棱柱層網(wǎng)格生成器創(chuàng)建體網(wǎng)格，計算中采用全y＋處理邊界層。飛機表面邊界層厚度定義為0.75 mm以滿足壁面函數(shù)要求，網(wǎng)格增長因子的數(shù)值定義為1.3，設(shè)置棱柱層的層數(shù)為8。切割六面體網(wǎng)格生成前，一方面對機頭、機翼、方向舵和升降舵等重點部位附近進行網(wǎng)格加密，另一方面對水表面附近的網(wǎng)格進行Z軸方向垂向加密，以保證VOF能準確捕捉到自由液面的變化。劃分的網(wǎng)格如圖3所示，圖4為y＋值分布圖。

圖3 飛機機體及加密部件網(wǎng)格展示

圖4 機體表面y+值云圖

2.3 數(shù)值模擬驗證

因本文主要研究內(nèi)容為飛機在水面以上的飛行過程，且受風(fēng)洞水池實驗條件和成本等因素的制約，故本文通過比較數(shù)值仿真結(jié)果和風(fēng)洞試驗測得數(shù)據(jù)，分析模擬方法的可行性。在空氣域中分別模擬-6°～6°共13個工況下飛機的升力系數(shù)、阻力系數(shù)和升阻比，來流風(fēng)速需定義在該型飛機的巡航速度區(qū)間（100～120 km/h）以內(nèi)，模擬中的取值為31 m/s，迎角和側(cè)滑角均設(shè)置為0°。在保證模擬結(jié)果一定精確程度的前提下，對飛機外形結(jié)構(gòu)進行適當(dāng)?shù)暮喕?，刪減掉了起落架以便計算和分析。計算平均值后與風(fēng)洞實驗數(shù)據(jù)對比見圖5。

圖5 數(shù)值模擬結(jié)果與實驗數(shù)據(jù)對比圖

3 近海面飛行氣動特性分析

為了探究飛機距水面高度、迎角和海況不同時飛行氣動特性變化情況，本文以升力系數(shù)、阻力系數(shù)、壓力分布為主要指標展開研究。重點主要有三方面：首先模擬迎角為6°時，飛機從距水面15 m下降到距水面0.7 m的過程，探究不同飛行高度時靜水面對飛機升力系數(shù)、阻力系數(shù)和機體壓力分布的影響；其次在距靜水面高度2 m時分別對迎角為0°、3°、6°和9°的飛機進行平飛模擬，分析迎角改變后液面影響程度的變化；最后在VOF波模型中設(shè)置波高、波長等參數(shù)，模擬飛機同樣以6°迎角在波浪水面上空平飛，比較海況不同造成的影響。

3.1 距水面高度對氣動特性的影響

分析距水面高度對電動飛機近海面飛行氣動特性的影響，設(shè)置飛機平飛速度31 m/s，迎角6°，海浪等級為0。以飛機機頭最前端作為參照，并依據(jù)飛機尾端臨近水面時密集選點的原則選取距水面高度分別為15 m、2.5 m、2.0 m、1.0 m、0.9 m、0.7 m共6個工況，模擬每個工況下飛機的近海面飛行特性，分析升力系數(shù)和阻力系數(shù)的變化情況，得到時歷曲線如圖6所示。飛機機體在各不同高度時受到的剪切阻力和壓差阻力如圖7所示。

圖6 不同離水高度下的飛機升力、阻力系數(shù)時歷曲線

圖7 剪切、壓差阻力匯總圖

從圖6（a）中可看出隨著飛機距水面高度降低，升力系數(shù)呈現(xiàn)出先增大后減小的變化趨勢，從圖6（b）中可看出阻力系數(shù)出現(xiàn)先減小后增大的趨勢。升力系數(shù)增大與阻力系數(shù)減小這兩點符合水面效應(yīng)的典型作用規(guī)律，而后升力系數(shù)的減小與阻力系數(shù)的增大是因為當(dāng)飛機高度下降到0.9 m時，從圖7中的剪切阻力和壓差阻力變化情況，以及圖8飛機觸水前后空氣體積分數(shù)圖來看，飛機雖未觸水（圖8（a）），但水面效應(yīng)對飛機升力的提升效力減弱，升力系數(shù)已不再增大；高度繼續(xù)下降到0.7 m時，機尾部分觸水（圖8（b）），水面效應(yīng)的增升作用徹底消失，故升力系數(shù)減小且飛機受水流影響剪切阻力和壓差阻力同時增加導(dǎo)致阻力系數(shù)增大。結(jié)果說明：飛機在近海面飛行時，水面效應(yīng)區(qū)有效利用的前提是必須嚴格控制飛行高度，離水面過低時輕則增升效應(yīng)減弱，重則觸水后導(dǎo)致阻力急劇增大得不償失，并有水載荷超過飛機材料本身承受能力的危險。圖9所示為飛機觸水前的壓力云圖，可以看出飛機下降過程中在機翼上方出現(xiàn)負壓區(qū)，機翼上方壓力呈現(xiàn)先增大后減小的趨勢，而飛機機頭上的壓力值幾乎一直不變。

圖8 飛機觸水前后空氣體積分數(shù)圖

圖9 不同高度下飛機表面的壓力云圖（靜水）

3.2 迎角變化對氣動特性的影響

減小定義飛機的平飛速度31 m/s，通過模擬不同迎角下飛機近海面的飛行過程，分析升力和阻力系數(shù)的變化情況，研究飛機迎角不同時受到靜水面影響的差異。水上飛機水面正常降落時一般迎角較小，接水前飛行員改平，因此在迎角區(qū)間內(nèi)等差選取模擬迎角為0°、3°、6°、9°，各工況下升力系數(shù)時歷曲線如圖10（a）所示，阻力系數(shù)時歷曲線如圖10（b）所示。

圖10 不同迎角的升力、阻力系數(shù)時歷圖

通過對比分析圖10中不同迎角下飛機的升力、阻力系數(shù)的變化可以看出在不超過失速迎角的前提下，升力、阻力系數(shù)的變化趨勢與迎角的變化趨勢是一致的，這種變化趨勢與在地面飛行時相同[25]。

3.3 海況變化對氣動特性的影響

論文最后分析了海況變化（水面波浪和風(fēng)速）對氣動特性的影響。首先對波浪進行建模，設(shè)置波高為0.2 m，波長為10 m，設(shè)置來流速度31 m/s，飛機迎角6°。同樣選取距水面高度分別為15 m、2.5 m、2.0 m、1.0 m、0.9 m、0.7 m這6個工況進行監(jiān)測。與靜水條件下的結(jié)果對比分析后，發(fā)現(xiàn)在海面有波浪時，飛機的升力、阻力系數(shù)在1.5 s進入穩(wěn)定期后并不像靜水條件下那樣幾乎一成不變，而是像波浪一樣作周期性變化，這說明飛機承受著周期性脈動載荷，但兩者相位并不相同，這種變化與波浪高度有關(guān)。將有波浪時的各工況結(jié)果取穩(wěn)定后4個周期數(shù)據(jù)進行平均，發(fā)現(xiàn)飛機飛行高度為0.9 m時，飛機尾部在每一個波浪周期內(nèi)的某一瞬間觸及波峰，從圖11（a）中發(fā)現(xiàn)飛機不持續(xù)接水時對升力系數(shù)有較小影響，且從圖11（b）中0.9 m時的阻力系數(shù)值依然小于1.0 m時的阻力系數(shù)值這點來看，飛機在0.9 m高度時并不受到周期性波浪力的影響；飛行高度0.7 m的飛機顯然會受到周期性波浪力影響，升力系數(shù)減小，阻力系數(shù)增大，且由于在一個觸水周期內(nèi)觸水深度周期性變化，兩者也出現(xiàn)周期性變化。

圖11 不同高度下升力、阻力系數(shù)時歷圖（波浪）

分別設(shè)置來流風(fēng)速依次為27.7 m/s、30.25 m/s、32.8 m/s、35.35 m/s、37.9 m/s，速度設(shè)置的依據(jù)是以該型飛機巡航速度區(qū)間為上下限并在中間區(qū)段等差取3個不同速度值，高度為距水面15 m，迎角6°，探究不同風(fēng)速對飛機升力系數(shù)和阻力系數(shù)如圖12所示，可以看出隨著風(fēng)速的不斷提高，飛機的升力系數(shù)和阻力系數(shù)都持續(xù)增大，符合空氣動力學(xué)規(guī)律。

圖12 不同風(fēng)速下飛機升力、阻力系數(shù)時歷圖

4 結(jié) 論

本文基于數(shù)值模擬技術(shù)對某型電動飛機模型近海面飛行的氣動特性進行研究，建立了電動飛機近海面飛行的數(shù)值模擬方案，比較分析了距水面高度、迎角和海況不同時氣動特性的變化情況。研究結(jié)果表明本文所用數(shù)值模擬方案可行；電動飛機距水面飛行高度降低時，升力系數(shù)增大，阻力系數(shù)減小，當(dāng)飛機繼續(xù)抵近海面時水面效應(yīng)的增升效果會減弱直至觸水后完全消失，升力系數(shù)要先于阻力系數(shù)發(fā)生變化；飛機近海面飛行迎角變化引起的升力系數(shù)和阻力系數(shù)變化情況與在空氣域中規(guī)律一致；在海面有波浪時，升力系數(shù)、阻力系數(shù)曲線隨波浪作周期性變化，但其平均值與水面無波浪時的結(jié)果基本一致。

本文只是初步分析了電動飛機近海面氣動特性，接下來將考慮安裝浮筒后所帶來的影響和水面效應(yīng)與飛機動力耦合的情況，另外將重點研究著水過程中不同環(huán)境風(fēng)速和浪向?qū)︼w機動穩(wěn)定性的影響。