基于OpenFOAM的水陸兩棲飛機水面高速滑行研究

段旭鵬，孫衛平，魏猛，楊永,*

1.西北工業大學 航空學院，西安 710072 2.中航通用飛機有限責任公司，珠海 519040

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

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

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

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

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

圖1 CL-215水陸兩棲飛機水面高速滑行時刻照片Fig.1 Photograph of CL-215 amphibious aircraft taxiing at high speed on water

本文以OpenFOAM為平臺，在兩相流動態解算器內部加入激勵盤模型，對兩相流計算模塊和激勵盤方法進行了驗證和確認，最后采用建立的方法研究了水陸兩棲飛機水面高速滑行的非定常運動學與動力學特性。

1 數理模型和計算方法

1.1 兩相流控制方程

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

(1)

(2)

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

VOF方法是在Navier-Stokes方程的基礎上求解體積分數的，α隨著速度場U運動的輸運方程為

(3)

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

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

(4)

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

(5)

式中：ρair和ρwater分別為空氣和水的密度；μair和μwater分別為空氣和水的動力黏性系數。

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

(6)

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

1.2 數值求解方法

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

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

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

1) 根據科朗數計算時間步長。

2) 求解網格運動方程，對網格進行變形。

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

4) 采用MULES方法求解式(3)，得到體積分數。

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

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

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

(7)

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

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

1.3 動網格方法

(8)

式中：γ為常數；u為網格變形速度。令xnew為新時刻網格坐標，xold為上一時刻網格坐標，Δt為時間步長，則有

xnew=xold+Δtu

(9)

1.4 6DOF方程

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

圖2 6DOF軸系和運動形式Fig.2 Axis and movement form of 6DOF

(10)

(11)

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

1.5 激勵盤的添加

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

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

(12)

(13)

(14)

式中：fbX和fbθ分別為軸向和切向體積力；r為激勵盤中某個確定點的徑向坐標；RP為槳盤外徑；RH為槳盤內徑；AX和Aθ為常數。令Δ為激勵盤的厚度，則AX和Aθ的表達式為

(15)

(16)

(17)

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

2 算法的驗證與確認

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

2.1 Wigley水池拖曳標模計算

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

(18)

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

(19)

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

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

圖3 Wigley計算網格和結果波形Fig.3 Computational grid and wave result of Wigley hull

表1 Wigley船體模型參數Table 1 Parameters of Wigley hull model

表2 Wigley計算網格收斂性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 兩相流激勵盤模型驗證

大型水陸兩棲飛機在國內FL12風洞進行了螺旋槳單槳測力試驗，本節采用試驗數據結合非定常計算的數值模擬結果對激勵盤模型進行驗證。如圖5所示，FL12風洞是4 m×3 m閉口回流式低速風洞，試驗模型為1∶15的六葉單槳和短艙的組合構型，動力模塊采用電機驅動，使用機械天平測定拉力系數。試驗的來流風速為40 m/s，模型迎角和側滑角均為0°，螺旋槳槳葉角為31°，轉速為10 000 r/min。

圖5 單槳測力低速風洞試驗現場照片Fig.5 Photograph of single propeller force measurement in low speed wind tunnel

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

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

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

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

3 水上飛機兩相流數值模擬結果及分析

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

3.1 水上飛機的計算網格

在不影響主要研究目標的情況下，對飛機外形進行了簡化，去掉了浮筒、發房和起落架艙等部件。計算網格仍然采用對水面適應性較好的笛卡兒網格。雖然結構網格同樣具有流向延展性的特點，但是受物面網格拓撲的影響，在機身近水面附近的網格延伸方向往往會被迫改變，導致自由面計算的誤差增加。幸運的是，笛卡兒網格一方面可以嚴格按照水平方向進行各向異性的網格加密，另一方面在水面與物面相交處，它具有非結構網格的特點，可以靈活對物面進行貼合，而不改變自由面的網格拓撲方向。

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

圖7 水陸兩棲飛機計算網格Fig.7 Computational grid of amphibious aircraft

3.2 飛機水面滑行的非定常計算歷程

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

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

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

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

圖10顯示了穩定滑行時的滑流速度剖面，來流馬赫數為0.12，經過加速后最大可達0.28。在槳盤后等距切出若干剖面，觀察流過機翼和襟翼的滑流形狀，可以看到在激勵盤所在區域流場得到加速，加速區域為圓環形狀，這與真實槳盤的效果相同。順航向在內、外發房中心位置分別截取兩個剖面。觀察外側剖面的流線，可以更加明顯地看到，由于槳盤的加速作用，流管出現明顯的收縮，隨后受機翼和增升裝置影響，滑流在機翼前后緣產生明顯的上洗和下洗，最后沿自由面流向后方。

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

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

3.3 滑流與動力影響分析

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

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

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

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

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

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

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

表3 飛機水面滑行V=40 m/s的航態和氣動力比重對比

3.4 飛機高速滑行時的自由面結構

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

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

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

圖14 穩定滑行時的噴濺結構Fig.14 Spray structure during stable taxiing

3.5 螺旋槳滑流和地面效應對高升力系統的影響

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

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

圖16顯示了剖面與機翼、襟翼交線上的壓力系數(Cp)分布對比,圖中X/c為無量綱弦向位置。關閉發動機的狀態下，地面效應對機翼上表面的影響只在前緣增加了吸力峰，而下表面受到氣流速度阻滯的影響，由伯努利方程可知該區域的壓力就會增大，圖中藍色下半部曲線也證實了這一點?？梢钥吹綁毫Φ脑黾友叵蚁蛞恢背掷m到襟翼縫道底部，使氣流加速進入襟翼縫道，進而使得襟翼上表面的附面層能量得到補充，使分離得到了抑制。在發動機和地面效應同時開啟狀態下，機翼和襟翼上表面吸力沿弦向均勻增加，這是因為水上飛機的發動機安裝位置較高，使得滑流大部分從高升力系統的上方通過，更好地將滑流的流向速度轉化為更多升力以幫助飛機盡快離水。同時觀察黑色曲線的下半部，可以看到在壓力表面，機翼和襟翼前緣正壓力峰值明顯，這得益于機翼下方復合氣流通道的一部分貢獻。觀察速度場，正是在機翼和襟翼的前緣下方存在兩個流速最低的區域，所以再次由伯努利定律可知相應位置的壓力才會顯著增加。除此以外，與文獻[26]的研究結果類似，襟翼前緣這部分正壓力將會更加明顯地對縫道內的氣流加速，氣流繞襟翼頭部的圓弧流動將形成柯恩達效應，從而進一步抑制襟翼上表面氣流分離。因此，可以說水上飛機水面高速滑行時地面效應和螺旋槳滑流的耦合使得機翼高升力系統的升力特性得到了顯著的提升。

圖15 距離對稱面5 m處過內側發動機中心的剖面速度分布和襟翼上表面流線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種狀態的機翼和襟翼壓力系數對比Fig.16 Comparison of pressure coefficient distribution over the wing and flap in 3 cases

4 結 論

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

2) 計算結果顯示，當水陸兩棲飛機達到單斷階高速滑行時，并不是像陸地起飛一樣平穩，而是會伴有小幅的俯仰和升沉振蕩。

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

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

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