段旭鵬，孫衛(wèi)平，魏猛，楊永,*

1.西北工業(yè)大學(xué) 航空學(xué)院，西安 710072 2.中航通用飛機(jī)有限責(zé)任公司，珠海 519040

水陸兩棲飛機(jī)的水面高速滑行是力學(xué)現(xiàn)象十分復(fù)雜的過程，涉及到固、液、氣三者的高速劇烈作用。飛機(jī)水上起飛較陸上起飛有很大不同，飛行員要在劇烈變化的水動力和高強(qiáng)度滑流作用下保持對飛機(jī)姿態(tài)和速度的有效控制，這對滑行安全性提出了更為苛刻的要求。

飛機(jī)水上起飛是一個動態(tài)過程，隨著速度的增加，水線和水動力不斷變化，氣動力從無到有逐漸建立，要了解這一動態(tài)過程需要借助空氣動力學(xué)、水動力學(xué)以及飛行力學(xué)等學(xué)科的相關(guān)理論知識。早期水上飛機(jī)除了要進(jìn)行風(fēng)洞試驗(yàn)以外，船體的設(shè)計(jì)還需要通過大量的水動試驗(yàn)進(jìn)行評估。例如20世紀(jì)90年代問世的CL-215[1]水陸兩棲滅火飛機(jī)，其船體尺寸和外形是通過一系列水池試驗(yàn)確定的。隨著計(jì)算流體力學(xué)(CFD)技術(shù)的發(fā)展，研究者不斷提高數(shù)值求解的精度，開始逐漸采用數(shù)值模擬方法，而減少試驗(yàn)的工作量。在過去的幾年里，國際水池拖曳會議(ITTC)開展了很多船舶標(biāo)模的阻力預(yù)測計(jì)算研討會[2]。美國航空航天學(xué)會資助的阻力預(yù)測會議(DPW)致力于提高解算器的求解能力，使得CFD對跨聲速氣動問題的阻力計(jì)算精度不斷得到改善[3]?，F(xiàn)在，CFD技術(shù)已經(jīng)具備了研究水上飛機(jī)相關(guān)問題的模擬能力。

兩相流的數(shù)值研究雖然比較多，但是針對水陸兩棲飛機(jī)的研究還不多見。上海交通大學(xué)的Qiu和Song[4]利用解耦式算法，基于商業(yè)軟件開展了水陸兩棲飛機(jī)起飛過程的動態(tài)模擬。該方法將氣動力和水動力分開計(jì)算，有較高的計(jì)算效率，但是尚未考慮螺旋槳動力影響問題。飛機(jī)水上迫降與此類似，北京航空航天大學(xué)的Qu等[5]計(jì)算了支線客機(jī)的水上迫降問題，采用VOF(Volume Of Fluid)、6DOF(6 Degree Of Freedom)以及全局動網(wǎng)格技術(shù)研究了上單翼、高平尾飛機(jī)的水上迫降性能。

隨著航空和航海技術(shù)的發(fā)展，對螺旋槳空氣滑流和高速船體的水動力特性研究開展的都比較多，但是將這兩方面內(nèi)容耦合研究的工作也不多見。目前中國正在研制大型水陸兩棲飛機(jī)，對于此類以渦槳為動力的大型飛機(jī)來說，如圖1所示，水面滑行將產(chǎn)生大量的噴濺，隨著速度的增加，噴濺基線不斷后移，后機(jī)身周圍產(chǎn)生涌起和凹陷，水動力矩隨之變化，加上螺旋槳滑流的非定常作用，使得水上飛機(jī)的受力較陸基飛機(jī)更為復(fù)雜，設(shè)計(jì)和駕駛難度更大。要了解這一復(fù)雜流固耦合問題，需要對空氣中的螺旋槳滑流和氣水交界面同時研究。

數(shù)值研究上述問題，需要VOF方法、動網(wǎng)格技術(shù)、六自由度方程求解和螺旋槳滑流模擬等多項(xiàng)技術(shù)的耦合。由于兩相流計(jì)算本身較為復(fù)雜，如果再直接模擬螺旋槳的非定常旋轉(zhuǎn)，那么計(jì)算將付出高昂的代價。國外對螺旋槳理論的研究起步較早，早期的螺旋槳理論主要為動量理論和葉素理論[6-9]，隨后基于上述理論方法國外學(xué)者提出了“激勵盤”模型[10-11]。該模型具有方法簡單、計(jì)算量少的特點(diǎn)，在一定的精度范圍內(nèi)可以替代3D實(shí)槳產(chǎn)生的滑流效應(yīng)，同時可以大大簡化用于模擬螺旋槳動態(tài)旋轉(zhuǎn)的計(jì)算量。因此用激勵盤模型來代替螺旋槳非定常計(jì)算的策略對于簡化水陸兩棲飛機(jī)的水面滑行問題來說是一個十分理想的選擇。

圖1 CL-215水陸兩棲飛機(jī)水面高速滑行時刻照片F(xiàn)ig.1 Photograph of CL-215 amphibious aircraft taxiing at high speed on water

本文以O(shè)penFOAM為平臺，在兩相流動態(tài)解算器內(nèi)部加入激勵盤模型，對兩相流計(jì)算模塊和激勵盤方法進(jìn)行了驗(yàn)證和確認(rèn)，最后采用建立的方法研究了水陸兩棲飛機(jī)水面高速滑行的非定常運(yùn)動學(xué)與動力學(xué)特性。

1 數(shù)理模型和計(jì)算方法

1.1 兩相流控制方程

VOF方法是一種自由面捕捉方法，其思想較為簡單：每一個計(jì)算網(wǎng)格單元中都額外包含一項(xiàng)代表水的體積分?jǐn)?shù)α，由初始?xì)馑唤缑鎺缀挝恢么_定初值。位于交界面的網(wǎng)格單元既含有氣體又有水，因此0<α<1；位于交界面之上的網(wǎng)格單元都是氣體，因此α=0；位于交界面以下的網(wǎng)格單元都是水，因此α=1。在OpenFOAM中，牛頓流體不可壓守恒形式的雷諾平均Navier-Stokes方程組為

(1)

(2)

式中：U為速度矢量；ρ為控制體內(nèi)的流體密度；p為壓力；g為重力加速度；μ為動力黏性系數(shù)；μt為湍流動力黏性系數(shù)。式(1)為質(zhì)量守恒方程，式(2)為動量守恒方程。

VOF方法是在Navier-Stokes方程的基礎(chǔ)上求解體積分?jǐn)?shù)的，α隨著速度場U運(yùn)動的輸運(yùn)方程為

(3)

式中：α=Vwater/Vtotal，數(shù)值介于0和1之間，其中Vwater為水在單元內(nèi)所占體積，Vtotal為單元總體積，因此流體空間內(nèi)的密度和黏性可表示為

ρ=ρa(bǔ)ir(1-α)+ρwaterα

(4)

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

(5)

式中：ρa(bǔ)ir和ρwater分別為空氣和水的密度；μair和μwater分別為空氣和水的動力黏性系數(shù)。

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

(6)

兩相流的湍流模型方程形式與單項(xiàng)流完全相同，唯一需要注意的就是密度和黏性需要用式(4)和式(5)表示。

1.2 數(shù)值求解方法

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

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

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

1) 根據(jù)科朗數(shù)計(jì)算時間步長。

2) 求解網(wǎng)格運(yùn)動方程，對網(wǎng)格進(jìn)行變形。

3) 更新兩相流運(yùn)動參數(shù)，包括密度和黏性系數(shù)。

4) 采用MULES方法求解式(3)，得到體積分?jǐn)?shù)。

5) 采用PISO算法求解速度和壓力。

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

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

(7)

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

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

1.3 動網(wǎng)格方法

(8)

式中：γ為常數(shù)；u為網(wǎng)格變形速度。令xnew為新時刻網(wǎng)格坐標(biāo)，xold為上一時刻網(wǎng)格坐標(biāo)，Δt為時間步長，則有

xnew=xold+Δtu

(9)

1.4 6DOF方程

6DOF是指剛體六自由度運(yùn)動，即任意給定剛體可以在三維空間沿3個轉(zhuǎn)軸平移和滾轉(zhuǎn)自由運(yùn)動，稱之為六自由度運(yùn)動。圖2顯示了剛體6個自由度的軸系和運(yùn)動形式，其中X、Y、Z表示在笛卡兒坐標(biāo)系中沿3個坐標(biāo)軸的平動，θX、θY、θZ表示繞3個坐標(biāo)軸的轉(zhuǎn)動。求解動網(wǎng)格方程之前，需要通過求解1.1節(jié)中的控制方程得到流場結(jié)果，隨后利用流場結(jié)果對剛體各物面進(jìn)行積分可得到剛體的受力以及對質(zhì)心的力矩，其具體表達(dá)式為

圖2 6DOF軸系和運(yùn)動形式Fig.2 Axis and movement form of 6DOF

(10)

(11)

式中：Fflow和Mflow為流場對物面的作用力和力矩；Fext和Mext為重力和其他外力及其作用力矩，對于本文的研究對象來說，F(xiàn)ext包括螺旋槳拉力和升降舵偏轉(zhuǎn)帶來的平尾附加壓力。在兩相流動態(tài)求解過程中，利用松耦合形式引入人工虛擬質(zhì)量項(xiàng)來對Navier-Stokes方程與六自由度方程進(jìn)行聯(lián)合求解。

1.5 激勵盤的添加

對于以渦槳為動力的水陸兩棲飛機(jī)來說，飛機(jī)起飛時螺旋槳滑流對飛機(jī)動力學(xué)影響不可忽視?，F(xiàn)有商業(yè)軟件、開源軟件中的VOF方法，無論是顯式、半顯式半隱式或者全隱式的時間推進(jìn)格式，其時間步長都很小，計(jì)算效率都不高，而且計(jì)算量比較大。如果再對螺旋槳進(jìn)行非定常模擬，計(jì)算耗時幾乎無法接受。

本文選用的激勵盤模型是在螺旋槳槳盤區(qū)域定義一個具有一定厚度的扁平圓盤，在圓盤邊界包絡(luò)的計(jì)算網(wǎng)格內(nèi)添加體積力源項(xiàng)，通過修改解算器中的控制方程以達(dá)到向流場注入動量的目的。Svenning[19]對較早版本的OpenFOAM單相流解算器進(jìn)行了激勵盤的添加，本文將在其工作的基礎(chǔ)上，將激勵盤模型加入到OpenFOAM 2.3.1中的動態(tài)兩相流解算器。主要改進(jìn)包括：激勵盤的空間位置從一次性計(jì)算改為隨飛機(jī)剛性實(shí)時運(yùn)動；體積力從固定不變改為實(shí)時更新，以適應(yīng)動態(tài)計(jì)算問題。動量的添加是通過添加體積力源項(xiàng)實(shí)現(xiàn)的，體積力的分布按照復(fù)雜程度可以分為均勻、非均勻以及基于螺旋槳性能分析的多種形式。為平衡精度與算法復(fù)雜度，本文采取Goldstein[20]的非均勻分布：

(12)

(13)

(14)

式中：fbX和fbθ分別為軸向和切向體積力；r為激勵盤中某個確定點(diǎn)的徑向坐標(biāo)；RP為槳盤外徑；RH為槳盤內(nèi)徑；AX和Aθ為常數(shù)。令Δ為激勵盤的厚度，則AX和Aθ的表達(dá)式為

(15)

(16)

(17)

其中:T為拉力；Q為扭矩。根據(jù)式(12)～式(16)，在已知螺旋槳旋轉(zhuǎn)產(chǎn)生的拉力T和扭矩Q的情況下，就可以得到槳盤沿軸向和切向的體積力fbX和fbθ的分布。此時，修改動量方程式(2)，以源項(xiàng)的形式將體積力fbX和fbθ添加到方程右端，式(2)就變成式(17)所示的新的動量方程。此時按照1.1節(jié)和1.2節(jié)介紹的方程和方法即可求得滑流作用下的流場結(jié)果。

2 算法的驗(yàn)證與確認(rèn)

在進(jìn)行水陸兩棲飛機(jī)水面高速滑行的數(shù)值模擬之前，應(yīng)首先對算法的氣動力和水動力計(jì)算模塊進(jìn)行耦合驗(yàn)證。受到設(shè)備條件的限制，水池拖曳試驗(yàn)尚無法對氣動力進(jìn)行準(zhǔn)確測量。因此本文將采用氣/水動力分別考核的思路來對算法進(jìn)行驗(yàn)證。具體來說，是采用Wigley水池拖曳試驗(yàn)和螺旋槳單槳風(fēng)洞試驗(yàn)的數(shù)據(jù)結(jié)果對新開發(fā)的動態(tài)激勵盤兩相流方法分別進(jìn)行水動力和氣動力方面的驗(yàn)證和確認(rèn)。

2.1 Wigley水池拖曳標(biāo)模計(jì)算

在船舶領(lǐng)域，Wigley標(biāo)模外形簡單有豐富的試驗(yàn)數(shù)據(jù)[21]，是考核水面航行體計(jì)算有效性的典型算例。該標(biāo)模有3種船體截面形狀：拋物線型、直線型、矩形，船體長寬比也有不同參數(shù)。本文選取長寬比為10的拋物線型船體為研究對象，如圖3所示。模型方程如式(18)所示，參數(shù)見表1。

(18)

采用笛卡兒網(wǎng)格對計(jì)算域離散，并對水面加密。中等網(wǎng)格的船體面網(wǎng)格最小尺度為40 mm，即1%L，首層附面層厚度為3.3 mm，對應(yīng)的y+值約為200，膨脹比率為1.3，如圖3所示。在此基礎(chǔ)上，為了進(jìn)行網(wǎng)格收斂性研究，還生成了另外兩套密度不同的網(wǎng)格，3套網(wǎng)格的體單元數(shù)目分別為：80萬、120萬、360萬。計(jì)算使用的弗洛德數(shù)表達(dá)式為

(19)

式中：V為航行體運(yùn)動速度;g重力加速度的大小。

計(jì)算時Fr=0.267，雷諾數(shù)Re=6.66×106，湍流模型為SST(Shear Stress Transport)。對3套不同密度、不同y+的網(wǎng)格進(jìn)行了計(jì)算，得到的結(jié)果如表2所示。試驗(yàn)得到的CD=0.004 16，雖然y+=60的加密網(wǎng)格計(jì)算精度最高，但是當(dāng)y+達(dá)到200，壁面網(wǎng)格增長率低于1.3以后，結(jié)果精度隨網(wǎng)格量的增加將變得不明顯。3套網(wǎng)格對自由面的捕捉精度較好，相比于試驗(yàn)結(jié)果，差異主要在于船艏處的波高，如圖4所示?；谏鲜鼋Y(jié)果可以得出，y+越小結(jié)果越準(zhǔn)確，然而與單相流氣動力計(jì)算有所不同，水動力計(jì)算的y+達(dá)到200后，就能得到較為理想的結(jié)果。

圖3 Wigley計(jì)算網(wǎng)格和結(jié)果波形Fig.3 Computational grid and wave result of Wigley hull

表1 Wigley船體模型參數(shù)Table 1 Parameters of Wigley hull model

表2 Wigley計(jì)算網(wǎng)格收斂性Table 2 Computational grid convergence of Wigley hull

圖4 Wigley船身波形(Fr=0.267，Re=6.66×106)Fig.4 Wave profile on Wigley hull (Fr=0.267，Re=6.66×106)

2.2 兩相流激勵盤模型驗(yàn)證

大型水陸兩棲飛機(jī)在國內(nèi)FL12風(fēng)洞進(jìn)行了螺旋槳單槳測力試驗(yàn)，本節(jié)采用試驗(yàn)數(shù)據(jù)結(jié)合非定常計(jì)算的數(shù)值模擬結(jié)果對激勵盤模型進(jìn)行驗(yàn)證。如圖5所示，F(xiàn)L12風(fēng)洞是4 m×3 m閉口回流式低速風(fēng)洞，試驗(yàn)?zāi)Ｐ蜑?∶15的六葉單槳和短艙的組合構(gòu)型，動力模塊采用電機(jī)驅(qū)動，使用機(jī)械天平測定拉力系數(shù)。試驗(yàn)的來流風(fēng)速為40 m/s，模型迎角和側(cè)滑角均為0°，螺旋槳槳葉角為31°，轉(zhuǎn)速為10 000 r/min。

圖5 單槳測力低速風(fēng)洞試驗(yàn)現(xiàn)場照片F(xiàn)ig.5 Photograph of single propeller force measurement in low speed wind tunnel

首先通過非定常單槳計(jì)算對試驗(yàn)狀態(tài)進(jìn)行復(fù)現(xiàn)，目的是得到足夠豐富的流場數(shù)據(jù)用以驗(yàn)證激勵盤的計(jì)算能力。這里采用滑移動網(wǎng)格技術(shù)模擬槳盤轉(zhuǎn)動，湍流模型為SST，非定常時間步長采用槳盤旋轉(zhuǎn)3.6°所用時間dt=6×10-5s。非定常單槳計(jì)算得到的拉力系數(shù)Tc為0.059 5，比試驗(yàn)結(jié)果小10%。

體積力的添加不影響計(jì)算的復(fù)雜度，在耦合激勵盤的兩相流非定常計(jì)算過程中，時間推進(jìn)步長僅僅受到VOF方法限制，對于水上飛機(jī)的高速滑行問題來說，通常是在10-3s量級。相對而言，螺旋槳旋轉(zhuǎn)的非定常模擬需要采用更為復(fù)雜的動網(wǎng)格方法，比如重疊網(wǎng)格、滑移網(wǎng)格，即便采用非慣性系簡化方法，對時間步長的選擇都有嚴(yán)格的限制。對于本文研究的問題來說，如果計(jì)算3D實(shí)槳的旋轉(zhuǎn)，時間步長一般在10-5s量級。因此激勵盤的選用從時間成本上來看，可以使計(jì)算保持較為經(jīng)濟(jì)的運(yùn)行方式，具有較強(qiáng)的工程實(shí)用價值。

激勵盤的驗(yàn)證計(jì)算采用了3套不同密度的網(wǎng)格，網(wǎng)格節(jié)點(diǎn)數(shù)分別為稀疏網(wǎng)格64萬，中等網(wǎng)格200萬，細(xì)化網(wǎng)格620萬。以非定常計(jì)算得到的拉力和扭矩作為輸入進(jìn)行計(jì)算，結(jié)果如圖6所示：截取流場中心剖面觀察馬赫數(shù)(Ma)云圖，3D實(shí)槳與激勵盤二者的計(jì)算結(jié)果總體上是相似的，流場經(jīng)過槳盤加速流動，形成包裹發(fā)房的一束高速射流，即“滑流”。細(xì)微的差別在于，前者在繞過發(fā)房的頭部邊界層內(nèi)有小范圍速度減緩，隨著氣流的行進(jìn)滑流邊界開始變得模糊，滑流區(qū)域開始向中心收縮；與此相對，后者滑流強(qiáng)度較為均勻，邊界更為清晰平直。此外為進(jìn)行定量對比，在圖6(a)的坐標(biāo)系下，在槳盤中心后方3倍半徑處，從發(fā)房表面沿氣流垂向向上分別截取線段分析速度型,如圖6(c)所示，其中Z/R是以槳盤半徑R無量綱化得到的流場內(nèi)的位置高度。3套 網(wǎng)格的結(jié)果差異較小，說明激勵盤計(jì)算對空間離散精度要求不高。激勵盤滑流束狀直徑略小，其馬赫數(shù)峰值較非定常結(jié)果小0.01，這對復(fù)雜的全機(jī)兩相流非定常計(jì)算來說是可以接受的。

圖6 激勵盤與3D實(shí)槳非定常模擬結(jié)果對比Fig.6 Comparison between simulation results of actuator disk and unsteady computation of 3D propeller

3 水上飛機(jī)兩相流數(shù)值模擬結(jié)果及分析

國內(nèi)外學(xué)者對螺旋槳滑流現(xiàn)象以及滑流和地效的耦合問題已經(jīng)做了較為深入的研究，對水面高速航行體的研究也取得了較為豐富的成果[22-24]。本節(jié)將基于對以上成果的研究認(rèn)識，嘗試進(jìn)一步探索在滑流、地效、自由面三者的共同作用下，飛機(jī)高速滑行的非定常運(yùn)動過程。利用經(jīng)過驗(yàn)證的動態(tài)激勵盤兩相流算法進(jìn)行高速滑行數(shù)值模擬，本節(jié)首先給出非定常計(jì)算的時間歷程，然后分析該階段滑流、拉力、平尾控制力對全機(jī)運(yùn)動學(xué)和動力學(xué)特性的影響，最后給出自由面受到擾動之后的典型特征。

3.1 水上飛機(jī)的計(jì)算網(wǎng)格

在不影響主要研究目標(biāo)的情況下，對飛機(jī)外形進(jìn)行了簡化，去掉了浮筒、發(fā)房和起落架艙等部件。計(jì)算網(wǎng)格仍然采用對水面適應(yīng)性較好的笛卡兒網(wǎng)格。雖然結(jié)構(gòu)網(wǎng)格同樣具有流向延展性的特點(diǎn)，但是受物面網(wǎng)格拓?fù)涞挠绊?，在機(jī)身近水面附近的網(wǎng)格延伸方向往往會被迫改變，導(dǎo)致自由面計(jì)算的誤差增加。幸運(yùn)的是，笛卡兒網(wǎng)格一方面可以嚴(yán)格按照水平方向進(jìn)行各向異性的網(wǎng)格加密，另一方面在水面與物面相交處，它具有非結(jié)構(gòu)網(wǎng)格的特點(diǎn)，可以靈活對物面進(jìn)行貼合，而不改變自由面的網(wǎng)格拓?fù)浞较颉?/p>

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

圖7 水陸兩棲飛機(jī)計(jì)算網(wǎng)格Fig.7 Computational grid of amphibious aircraft

3.2 飛機(jī)水面滑行的非定常計(jì)算歷程

飛機(jī)水上起飛過程分為兩個階段，排水運(yùn)動階段和斷階滑行階段。飛機(jī)啟動后，受螺旋槳拉力作用開始低速航行，船體的斷階繞流造成底部壓力降低和流動分離，使得飛機(jī)受到向下的“吸”力，飛機(jī)下沉水線上升。隨著速度的增加，水下流動分離范圍逐漸增加，負(fù)壓達(dá)到極值，飛機(jī)下沉至最低點(diǎn)，阻力達(dá)到峰值即“阻力峰”。利用螺旋槳的剩余拉力飛機(jī)通過阻力峰繼續(xù)加速，氣動力開始逐漸起作用，飛機(jī)受氣動升力和水動升力共同作用克服重力不斷抬升，斷階后部開始進(jìn)入空氣，飛機(jī)從雙斷階滑水逐漸過渡為單斷階滑行，當(dāng)氣動升力增加到與重力相等時飛機(jī)離水完成起飛。

本文聚焦起飛過程中的斷階滑行階段。飛機(jī)在斷階滑行過程中物理機(jī)制涉及3個領(lǐng)域，包括空氣動力學(xué)、水動力學(xué)和飛行力學(xué)。對于斷階滑行階段，飛機(jī)的動態(tài)穩(wěn)定性至關(guān)重要，受穩(wěn)定性邊界影響，必須對飛機(jī)加以控制才能避免進(jìn)入“海豚跳”[25]的危險狀況。對于這種復(fù)雜過程的非定常數(shù)值計(jì)算來說，則需要考慮平尾的控制力，才能避免大幅振蕩從而得到穩(wěn)定的計(jì)算結(jié)果。如圖8所示，除了螺旋槳滑流以外，飛機(jī)的滑行還需要在螺旋槳拉力和升降舵附加力的共同作用下達(dá)到平衡。這兩種外力的添加是通過修改OpenFOAM六自由度動態(tài)計(jì)算函數(shù)，添加剛體外力模塊實(shí)現(xiàn)的。相對于模擬真實(shí)的升降舵偏轉(zhuǎn)，該方法進(jìn)一步提高了計(jì)算效率。因此，本節(jié)的非定常計(jì)算過程是首先利用Navier-Stokes方程耦合VOF和激勵盤方法求解兩相流問題獲得飛機(jī)的氣動力和水動力，然后通過添加的剛體外力模塊施加螺旋槳拉力和升降舵附加力，最后將物面受到的流體作用力和剛體外力一起代入到飛行力學(xué)6DOF方程進(jìn)行計(jì)算，即可得到飛機(jī)滑行的運(yùn)動特性。

非定常計(jì)算歷程如圖9所示，圖中Zcoord為自由液面的Z坐標(biāo)即自由面高度，Zcoord=0表面液面未受擾動，t為時間，AOA為迎角。計(jì)算速度設(shè)定為40 m/s，給定4°初始迎角，在螺旋槳滑流、拉力和平尾附加力的共同作用下，經(jīng)過15 s物理時間飛機(jī)穩(wěn)定在5.5°～5.8°之間，重心高度較靜浮狀態(tài)升高約1.28 m。全機(jī)阻力建立得較快，5 s時間即達(dá)到穩(wěn)定，阻力系數(shù)約為0.45，注意這是全機(jī)阻力，包括氣動阻力和量值較大的水動阻力?？傋枇εc全機(jī)重力之比約為0.16。根據(jù)迎角變化歷程曲線，選擇3個俯仰最為劇烈的典型時刻觀察飛機(jī)姿態(tài)和流場情況。首先可以明顯看到流線通過激勵盤以后的收縮和旋轉(zhuǎn)，其次斷階排開水形成噴濺自由面，自由面高度隨迎角增加而增加。t1時刻迎角為7.5°時的液面最大高度可達(dá)0.45 m，t2時刻迎角為4.5°時的液面最大高度為0.33 m。隨著計(jì)算時間的推移，飛機(jī)縱搖幅值逐漸減小趨于收斂，最終在t3時刻達(dá)到穩(wěn)定(俯仰振蕩幅值不超過2°[25])。此時的平尾附加下壓力為8 000 N，等效偏角約為-18°。

圖8 水陸兩棲飛機(jī)外力示意圖Fig.8 Sketch of external forces of amphibious aircraft

圖10顯示了穩(wěn)定滑行時的滑流速度剖面，來流馬赫數(shù)為0.12，經(jīng)過加速后最大可達(dá)0.28。在槳盤后等距切出若干剖面，觀察流過機(jī)翼和襟翼的滑流形狀，可以看到在激勵盤所在區(qū)域流場得到加速，加速區(qū)域?yàn)閳A環(huán)形狀，這與真實(shí)槳盤的效果相同。順航向在內(nèi)、外發(fā)房中心位置分別截取兩個剖面。觀察外側(cè)剖面的流線，可以更加明顯地看到，由于槳盤的加速作用，流管出現(xiàn)明顯的收縮，隨后受機(jī)翼和增升裝置影響，滑流在機(jī)翼前后緣產(chǎn)生明顯的上洗和下洗，最后沿自由面流向后方。

圖9 水陸兩棲飛機(jī)滑行計(jì)算歷程Fig.9 Computation history of amphibious aircraft during taxiing

圖10 滑流剖面的速度分布和流線形狀Fig.10 Velocity field distribution and slipstream at a certain slice

3.3 滑流與動力影響分析

目前除了俄羅斯的Be200以外，大部分水上飛機(jī)都是以渦槳為動力的，這是因?yàn)樗巷w機(jī)更加注重低速飛行性能，更重要的是利用渦槳的低速大拉力特點(diǎn)可以更快地起飛離水以提高安全性。螺旋槳滑流對水上起飛性能的影響較大，在設(shè)計(jì)中必須加以考慮。本節(jié)通過3個狀態(tài)有無滑流和動力的對比研究，嘗試揭示螺旋槳對水上起飛動態(tài)過程的影響規(guī)律。

3個狀態(tài)分別為：CC(Clean Case)—無動力無滑流、SO(Slipstream Only)—無動力有滑流、SF(Slipstream & Forces)—有動力有滑流并附加平尾操縱力。來流速度為40 m/s，采用與3.1節(jié)相同的網(wǎng)格，計(jì)算得到的運(yùn)動學(xué)和動力學(xué)參數(shù)時間歷程如圖11、圖12所示。

圖11中右側(cè)縱軸為飛機(jī)重心升沉變化(CG(Center of Gravity) rise)。從圖中可知，與無動力無滑流狀態(tài)CC相比，滑流產(chǎn)生一個較大的低頭力矩，使得狀態(tài)SO的迎角降低約2.5°。分析圖12曲線，左側(cè)縱軸為氣動力升力系數(shù)(AerodynamicCL)，右側(cè)縱軸為氣動阻力系數(shù)(AerodynamicCD)和總阻力系數(shù)(TotalCD)，通過能量的注入，升力系數(shù)大幅增加,從1增加到1.8，但是由于沒有平尾的附加控制，結(jié)合圖11，飛機(jī)的低頭抵消了升力系數(shù)的增加，使得飛機(jī)的重心沒有明顯的抬升。觀察圖12的阻力系數(shù)曲線，以右側(cè)縱軸為刻度，可以看到，飛機(jī)總阻力明顯高于氣動阻力，說明盡管飛機(jī)處于斷階滑水狀態(tài)，浸潤面積已經(jīng)很小，但是由于水的密度和黏性遠(yuǎn)遠(yuǎn)大于空氣，水阻力仍然占據(jù)相當(dāng)高的比重。對比氣動阻力發(fā)現(xiàn)，高速滑流流過機(jī)翼時，大幅提高了壓差阻力使得氣動阻力系數(shù)增加。然而氣水動總阻力沒有顯著變化，是因?yàn)轱w機(jī)低頭降低了水面滑行的水阻力，氣動力阻力的增加與水動力阻力的減小使得總阻力基本不變。

圖11 滑流影響下的飛機(jī)俯仰和重心升沉變化歷程Fig.11 History of aircraft trim angle and rise of center of gravity under effect of slipstream

圖12 滑流影響下的飛機(jī)氣動力和總阻力變化歷程Fig.12 History of aircraft aerodynamic forces and total drag under effect of slipstream

為盡量模擬真實(shí)起飛狀態(tài)，添加了螺旋槳的軸向拉力，由于水陸兩棲飛機(jī)采用位置較高的上單翼布局，拉力線高于重心高度，因此不可避免地會進(jìn)一步增加飛機(jī)的低頭力矩。為此，必須對飛機(jī)進(jìn)行控制才能保證高速滑行過程中飛機(jī)能夠保持一個合理的姿態(tài)。通過添加剛體隨體力模擬升降舵上偏18°的效果，此即狀態(tài)SF。觀察圖11的藍(lán)色曲線，飛機(jī)迎角回到理想的5.8°左右，縱搖幅值變小，變得更加穩(wěn)定。圖12顯示，迎角的增加使得對滑流的利用更為有效，藍(lán)色的氣動升力系數(shù)顯著增加到2左右，如此高的升力系數(shù)得益于滑流和地效的共同作用，這比無滑流和地效狀態(tài)提高了接近一倍。飛機(jī)重心平穩(wěn)提升，浸潤面積減小，水動阻力降低，總阻力系數(shù)隨之減小到0.5左右。正是由于迎角的增加，氣動阻力與單獨(dú)滑流作用的狀態(tài)相比是有所增加的。

表3顯示了3種狀態(tài)的航態(tài)和氣動力比重的對比數(shù)據(jù)。在螺旋槳滑流的單獨(dú)作用下，飛機(jī)明顯低頭，加上平尾控制力后，迎角增大，氣動升力占比從34%提高到77%，在滑流附加升力的作用下飛機(jī)進(jìn)一步抬升，導(dǎo)致水阻力降低。這說明把迎角控制得盡量高，而不超出穩(wěn)定性上邊界對飛機(jī)的水面高速滑行是十分有利的。顯而易見，滑流產(chǎn)生的附加阻力會使得氣動力阻力占比有一定的增加。

表3 飛機(jī)水面滑行V=40 m/s的航態(tài)和氣動力比重對比

3.4 飛機(jī)高速滑行時的自由面結(jié)構(gòu)

船體設(shè)計(jì)的一個主要目的是，使飛機(jī)有一個較為寬泛的穩(wěn)定滑行迎角范圍，從圖11中可以看到，在3°～8°范圍內(nèi)飛機(jī)均可穩(wěn)定滑行。斷階位置、高度和斷階后的船底后緣角這3個參數(shù)的綜合設(shè)計(jì)影響到單斷階滑行特性。從圖13中可以看到迎角為5.8°時，斷階兩側(cè)可以觀察到明顯的涌起，高度約為0.4 m，涌起形成后向飛機(jī)側(cè)后方流動。斷階脫出的自由面凹陷較深，使得后機(jī)身底部與水面還有較大空間，有利于后體的船底避免來自斷階的水流沖刷而造成阻力的增大。如果仔細(xì)觀察圖9可以看到，自由面形狀受滑行速度和滑行面迎角影響較大，迎角越大自由面涌起和凹陷的高度越大。結(jié)合圖9和圖13可知，在機(jī)身尾部，水面凹陷達(dá)到最深并開始逐漸恢復(fù)。從機(jī)翼后方的自由面來看，螺旋槳滑流盡管經(jīng)過襟翼有強(qiáng)烈下洗，但由于水的密度是空氣密度的1 000倍，除了受網(wǎng)格密度影響而無法模擬的飛沫和細(xì)絲等霧化結(jié)構(gòu)以外，滑流對水面主體形狀的影響是十分微小的。其次船體設(shè)計(jì)還要盡量減小水面沖擊載荷，以降低結(jié)構(gòu)重量。如圖13(b)所示，40 m/s的速度已接近離水速度，此時的水面沖擊壓力可達(dá)174 kPa(計(jì)算參考壓力為標(biāo)準(zhǔn)大氣壓，即100 kPa)。

此外，船體設(shè)計(jì)對噴濺特性的要求也很高，要求在高速滑行、轉(zhuǎn)彎以及波浪水面等各種滑行工況下，噴濺區(qū)域盡量避開螺旋槳、襟翼、尾翼等空氣動力學(xué)控制部件。與光滑粒子流體動力學(xué)(SPH)方法相比，VOF方法的重點(diǎn)在于模擬流體的力學(xué)特性，對自由面細(xì)微結(jié)構(gòu)捕捉能力是不足的。但是，當(dāng)自由液面區(qū)域的空間離散達(dá)到一定規(guī)模后，其自由面的主要結(jié)構(gòu)也是可以在一定程度上體現(xiàn)出來的。從本例的自由液面形狀來看，如圖14所示，基本的噴濺結(jié)構(gòu)已經(jīng)模擬出來，可以清晰看到斷階前受舭彎抑制而改變流向的須狀噴濺和斷階后涌起的主噴濺。

圖13 穩(wěn)定滑行時自由面形狀和船體壓力分布(AOA=5.8°)Fig.13 Free surface and distribution of pressure at hull during stable taxiing (AOA=5.8°)

圖14 穩(wěn)定滑行時的噴濺結(jié)構(gòu)Fig.14 Spray structure during stable taxiing

3.5 螺旋槳滑流和地面效應(yīng)對高升力系統(tǒng)的影響

為了進(jìn)一步分析螺旋槳滑流和地面效應(yīng)二者對高升力系統(tǒng)的影響，還在單項(xiàng)流不可壓解算器simpleFoam中添加了激勵盤模型，并專門進(jìn)行了無滑流無地效、無滑流有地效、滑流地效同時作用的3種狀態(tài)計(jì)算對比。模型迎角統(tǒng)一采用3.3節(jié)真實(shí)高速滑行狀態(tài)的計(jì)算結(jié)果5.8°，來流速度和地面高度同樣與該狀態(tài)相同。網(wǎng)格依然采用3.1節(jié)介紹的笛卡兒網(wǎng)格，與前述計(jì)算不同的是本節(jié)用滑移壁面代替了自由面。

過內(nèi)側(cè)發(fā)動機(jī)中心做一個平行于對稱面的縱向剖面，圖15顯示了3種狀態(tài)的速度場和襟翼上表面流線，并用來流速度對速度場進(jìn)行無量綱化，即以U/Uinlet表征無量綱速度。當(dāng)關(guān)閉發(fā)動機(jī)時，在地面效應(yīng)的單獨(dú)作用下(圖15(b))機(jī)翼上表面的速度場基本不變；然而機(jī)翼的下表面出現(xiàn)速度阻滯，其流速有一定程度的降低；襟翼上表面流速略微增加，進(jìn)而繞襟翼流動的分離區(qū)相應(yīng)減小。無地效狀態(tài)(圖15(a))之所以會出現(xiàn)襟翼分離區(qū)是因?yàn)轱w機(jī)還沒達(dá)到起飛速度，氣流動壓較低，襟翼后緣的逆壓梯度較大。打開發(fā)動機(jī)并耦合地面效應(yīng)(圖15(c))可以看到，機(jī)翼上表面流速顯著增加，襟翼上表面尾緣區(qū)域形成高速附著流；機(jī)翼下表面形成了一條由低速氣流和高速氣流組成的復(fù)合氣流通道；此時襟翼在高速附著流的作用下，其上表面氣流分離區(qū)完全消失。

圖16顯示了剖面與機(jī)翼、襟翼交線上的壓力系數(shù)(Cp)分布對比,圖中X/c為無量綱弦向位置。關(guān)閉發(fā)動機(jī)的狀態(tài)下，地面效應(yīng)對機(jī)翼上表面的影響只在前緣增加了吸力峰，而下表面受到氣流速度阻滯的影響，由伯努利方程可知該區(qū)域的壓力就會增大，圖中藍(lán)色下半部曲線也證實(shí)了這一點(diǎn)。可以看到壓力的增加沿弦向一直持續(xù)到襟翼縫道底部，使氣流加速進(jìn)入襟翼縫道，進(jìn)而使得襟翼上表面的附面層能量得到補(bǔ)充，使分離得到了抑制。在發(fā)動機(jī)和地面效應(yīng)同時開啟狀態(tài)下，機(jī)翼和襟翼上表面吸力沿弦向均勻增加，這是因?yàn)樗巷w機(jī)的發(fā)動機(jī)安裝位置較高，使得滑流大部分從高升力系統(tǒng)的上方通過，更好地將滑流的流向速度轉(zhuǎn)化為更多升力以幫助飛機(jī)盡快離水。同時觀察黑色曲線的下半部，可以看到在壓力表面，機(jī)翼和襟翼前緣正壓力峰值明顯，這得益于機(jī)翼下方復(fù)合氣流通道的一部分貢獻(xiàn)。觀察速度場，正是在機(jī)翼和襟翼的前緣下方存在兩個流速最低的區(qū)域，所以再次由伯努利定律可知相應(yīng)位置的壓力才會顯著增加。除此以外，與文獻(xiàn)[26]的研究結(jié)果類似，襟翼前緣這部分正壓力將會更加明顯地對縫道內(nèi)的氣流加速，氣流繞襟翼頭部的圓弧流動將形成柯恩達(dá)效應(yīng)，從而進(jìn)一步抑制襟翼上表面氣流分離。因此，可以說水上飛機(jī)水面高速滑行時地面效應(yīng)和螺旋槳滑流的耦合使得機(jī)翼高升力系統(tǒng)的升力特性得到了顯著的提升。

圖15 距離對稱面5 m處過內(nèi)側(cè)發(fā)動機(jī)中心的剖面速度分布和襟翼上表面流線Fig.15 Profile of velocity distribution across center of inner engine at the distance of 5 m to symmetry plane and streamlines over suction surface of flaps

圖16 3種狀態(tài)的機(jī)翼和襟翼壓力系數(shù)對比Fig.16 Comparison of pressure coefficient distribution over the wing and flap in 3 cases

4 結(jié) 論

1) 本文基于開源CFD軟件OpenFOAM發(fā)展的動態(tài)激勵盤兩相流計(jì)算方法可大幅降低計(jì)算復(fù)雜度，減少計(jì)算量，可應(yīng)用于以渦槳為動力的兩棲飛機(jī)水面高速滑行模擬。

2) 計(jì)算結(jié)果顯示，當(dāng)水陸兩棲飛機(jī)達(dá)到單斷階高速滑行時，并不是像陸地起飛一樣平穩(wěn)，而是會伴有小幅的俯仰和升沉振蕩。

3) 滑流在滑行中的作用明顯，顯著增加氣動升力和低頭力矩，因此需要增加升降舵的控制力才能更好地將姿態(tài)控制在較為穩(wěn)定的范圍。在穩(wěn)定性邊界范圍內(nèi)，盡量提高滑行迎角對飛機(jī)的快速離水是十分有利的。

4) 斷階這一水陸兩棲飛機(jī)特有的設(shè)計(jì)特點(diǎn)，不僅可以在高速滑行過程中減小水動阻力還可為飛機(jī)的縱向操縱提供更大的可用空間。

5) 水上飛機(jī)水面高速滑行時地面效應(yīng)和螺旋槳滑流的耦合使得機(jī)翼高升力系統(tǒng)的升力特性得到了顯著的提升。