趙立杰,田孟偉,李景奎,王明陽,劉達
1. 沈陽航空航天大學 航空宇航學院,沈陽 110136 2. 遼寧通用航空研究院,沈陽 110136 3. 沈陽航空航天大學 民用航空學院,沈陽 110136
隨著航空事業持續發展對環境帶來的負面影響,各國營造綠色航空的呼聲越來越高。新能源電動飛機的出現為航空的綠色化提供了一條光明的技術途徑[1]。
水上飛機不僅可以充分利用中國豐富的水域資源以滿足消防滅火等方面的特殊需求,對環境的影響也可以降到最低。而低空領域的逐漸開放推動了通用航空產業的迅速發展。水上電動飛機的出現是通航領域的一次大膽嘗試,以現有雙座電動飛機為基礎進行改裝設計,可以充分利用現有電推進技術,極大地推進綠色低碳航空產業的發展,填補相應的市場空白[2]。水上電動飛機的主要優勢體現在兩個方面:水面起降方式減少了建設機場占用的土地面積;電推進系統零排放避免了排放廢氣造成的空氣質量下降以及殘留油料對水體環境的污染。
浮筒式水上飛機相比船身式飛機結構更加簡便,易于改裝制造,且具有良好的橫向穩定性。相比較陸基飛機,水上飛機主要變化體現在結構和起降方式兩個方面:① 首先增加的浮筒結構使飛行時巡航阻力更高;其次由于大濕潤表面積,在滑行時會受較大水動阻力;這就在設計時對浮筒的外形有更高的要求。② 水面起降方式更加復雜,是一個動態的液氣耦合過程,而在滑行階段主要阻力來源于水動阻力,浮筒的設計需要具有良好的水動性能[3]。
過去水上飛機的研制多停留在水池拖拽實驗階段,但CFD(Computational Fluid Dynamics)技術的快速發展為水上飛行器的數值模擬計算提供了技術支持[4-9]。李新穎等[10]研究了水陸兩棲飛機高性能復合船型在規則波中航行時的耐波性能,得到了其在規則波中航行時的受力及運動響應。段旭鵬等[11]基于OpenFOAM研究了水陸兩棲飛機水面單斷階滑行過程中的氣動力和水動力特性,并給出了滑流和動力的影響規律。趙蕓可等[12]利用開發的整體網格方法對某型水上飛機水面降落的全過程進行了數值模擬,得到了較好的模擬結果。兩相流的數值研究使水上飛機的研制得到發展,目前主要應用在大型水陸兩棲飛機,針對雙浮筒式飛機的應用還比較少。
本文首先以設計的基準浮筒為基礎,利用解析函數線性疊加法對現有浮筒縱向截面進行參數化描述,結合近似模型和遺傳算法,以浮筒最大升阻比為目標進行優化設計。然后對加裝新浮筒的某型水上電動飛機利用非穩態多相流(VOF)模型方法進行起飛滑行數值模擬計算,重點探究水上電動飛機起飛階段阻力變化特性、流場噴濺特性以及浮筒壓力分布情況。最后對“阻力峰”下所需電推進系統功率進行驗證計算,總結浮筒式水上飛機起飛的一般性規律,以期為水上電動飛機的設計優化提供一定的借鑒。
浮筒基本外形按照基本剖面組合沿龍骨布置,既需要保持良好的流線外形以減少氣動阻力的影響,又需要要保證良好的水動力性能。浮筒結構設計中最為關鍵的就是斷階設計和縱向截面形狀的設計,要保證前體艏部曲翹程度足夠大,防止滑行時被水浪完全埋沒,前體底部縱剖線彎曲過渡要適中,保持良好的滑行性能,既能起到阻礙噴濺的作用,又能提供良好的浮性,在低速滑行時應該隨吃水深度的細微變化提供足夠的靜浮力和流動升力。
在起飛滑跑階段,當飛機在水中滑行時,水會沿著底部的凸面上移動,因此沿水面的法線方向產生吸力。這意味著飛機在水中加速行駛時,向下吸力會隨著速度的增加而增加。起飛性能將因此受到極大削弱,因此在浮筒的前后體設置斷開的臺階,隨著飛機速度的提高,會在后體處產生“氣袋”,斷階的存在能保證水上飛機在水中滑行時浮筒底部有足夠的空氣通量,使飛機擺脫水對它的吸附力而離水。隨著更多的空氣被吸入并開始向尾部方向延伸,最終使浮筒尾部與水完全分離。臺階的位置可以變化,但是通常放置在浮力中心后約浮筒長度的1/16處,臺階的高度為浮筒最大寬度的5%~9%[13]。
浮筒縱向截面形狀直接決定浸水部分的底面形狀。底部越平坦,飛機離水速度就越小。但是在波浪的撞擊作用下,較平坦的底部將導致嚴重的撞擊,這可能會使水上飛機的起飛離水速度比具有尖銳外形的起飛速度大,后者會穿透波浪而不是越過波浪[14]。這里基準浮筒選取雙凹面外形,它既有一個鋒利的邊緣,可以進入水中并切穿波浪,而且仍然具有相對平坦的曲面。雙凹面的另一個好處是,它可使噴濺向下偏轉,遠離浮筒和機身。
浮筒基準方案的確定通過以下步驟完成:
1) 根據整機排水重量計算浮筒體積。
2) 通過引入最大寬度載荷系數,考慮前后體比例系數、噴濺系數推導出最優浮筒結構的最大寬度和前后體長度。
3) 計算合適的斷階高度、底部斜升角和后緣角。
經過反復計算選擇進行建模,基準浮筒外形如圖1所示,具體參數如表1所示。

圖1 浮筒參數示意圖Fig.1 Schematic diagram of float parameters

表1 浮筒關鍵參數Table 1 Key parameters of float
原始浮筒安裝后,斷階滑行時阻力過大,水上飛機滑行時無法保持穩定,縱向穩定性較差,因此對原始浮筒縱向截面形狀進行優化設計。
浮筒水動性能對于曲線的變化十分敏感,側壁曲線彎度對浮筒滑行阻力有較大影響,通過拖船實驗可以得出陡峭的側壁可以減小水流附著的摩擦阻力,加快飛機離水起飛;截面曲線的微小變化也可能對其水動性能產生很大的影響,浮筒流線型設計及優化成為水上飛機起飛性能研究的重要課題。
結合解析函數線性疊加法與神經網絡模型及多島遺傳算法相結合,以某浮筒前體為研究對象,實現了對前體截面形狀的優化設計。而對于完整的迭代優化過程,如果每一步都進行CFD流場數值計算,計算量和所耗時間過大,不符合實際工程要求,因此引入近似模擬技術,先建立神經網絡模型再對模型進行優化計算可以有效縮短計算周期。優化流程如圖2所示,首先基于Isight平臺,將MATLAB、CATIA(Computer Aided Three-dimensional Interactive Application)、ICEM(Integrated Computer Engineering and Manufacturing code)和Fluent組件進行集成,利用拉丁超立方技術對設計變量進行抽樣;然后調用CATIA進行參數化建模;接著進行網格劃分和流場計算,利用樣本和所得水動性能響應構建神經網絡模型;最后利用多島遺傳算法求解。

圖2 截面形狀優化設計流程Fig.2 Section shape optimization design process
優化問題的數學模型為
maxm(ck)=CL/CD
(1)

以解析函數線性疊加法對浮筒剖面形狀進行描述,通過對原始外形坐標點添加擾動對其進行形狀調整,以此獲得水動性能更好的外形形狀[15-16]。
截面曲線的坐標函數為
(2)

為保證浮筒截面縱向深度及最大寬度保持不變,以此保證在xl= 0和xl= 0.28處(如圖3紅點)的擾動為0,同時保證xl= 1.00時變形不為0,對型函數進行了修改:
(3)
式中:e為固定常數,e(k)=lg 5/lgxl;m=11/39,而在明確了變形系數對浮筒阻力影響規律的情況下,對設計變量的范圍進行描述是更有效的。c1~c4控制浮筒底部曲線,取值范圍為[-0.002, 0.002],當c> 0時,c越大底部形狀越平坦,反之內凹加大;c5~c10控制側壁的陡峭程度和彎度,取值范圍為[-0.050, 0.050]。

圖3 浮筒截面示意圖Fig.3 Schematic diagram of float section
建立近似模型的前提是進行數據采樣,拉丁超立方設計以其有效的空間填充能力和適用于擬合非線性相應的優點而被廣泛應用。其原理是在設計區間里,將每一維的坐標區間均勻分為數個區間,在每個區間內隨機選取一個點,構成多維空間。在建模中選用LHS(Latin Hypercube Sampling)方法進行數據采集,樣本總數為50個。
在擬合問題中,利用神經網絡在變形系數與目標函數升阻比之間進行映射[17]。利用MATLAB中的神經網絡工具箱進行樣本數據的訓練并使用均方誤差和回歸分析評估其性能。樣本數據中,訓練數據、驗證數據和測試數據按照80%、10%、10%的比例進行分配。主要過程分為兩個部分:① 利用拉丁超立方實驗設計所得的樣本數據進行訓練,其中隱含神經元數量根據Kolmogorov定理確定[18];② 對近似模型進行準確度評價,主要通過規劃問題中的擬合程度進行評價,相關系數R越接近1,模型的擬合程度越好。圖4分別給出了近似模型訓練、驗證、測試和總體的線性擬合評估值,縱坐標為所擬合的模型函數預測值,從結果來看擬合效果良好,證明了模型的準確性和可靠性。


圖4 近似模型評估Fig.4 Approximate model evaluation
為減小計算量,選擇浮筒前體對稱結構的一半(見圖5)進行批量計算,以滿排水、水流速度為20 m/s的流場對模型進行了水動性能分析。利用命令流文件調用ICEM進行網格自動劃分;利用編寫的批量處理文件RunFluent調用Fluent進行求解計算,并將監控阻力和升力系數輸出到DAT文件中。

圖5 流場及網格Fig.5 Flow field and grid
由于在此優化問題中設計點數量較多,各優化目標之間往往存在不可兼容性,因此選用多島遺傳算法。多島遺傳算法一方面只評價設計點,不需要計算函數的梯度;另一方面具有全局性,能夠求解出全局最優解,避免陷入局部最優,從而對各優化目標實現協調優化,盡可能求出相對最優解。
浮筒截面優化前后形狀對比如圖6所示,優化后的外形底部內凹增大,側壁曲線沿縱向分布彎度變小,陡峭趨勢更加明顯,頂部平臺變寬,經過測量截面面積減小了8.9%。
優化后的浮筒性能較原浮筒的水動力性能有明顯改善,表2給出了優化前后浮筒水動力性能相關變化,可見阻力系數減小了1.64%,但主要體現在升力系數的增大,增幅達到19.20%,圖7為原始浮筒和改進后浮筒的壓力云圖,可見頂端峰值壓力減小,整體水阻力減小,水流流速受阻減小,浮筒底部負壓明顯下降。

圖6 原始截面與優化截面外形Fig.6 Original section and optimized section shapes

表2 浮筒前體水動性能對比Table 2 Comparison of hydrodynamic performance of float precursor

圖7 壓力云圖Fig.7 Pressure cloud diagram
實際安裝后,水上飛機停浮時可見姿態比優化前要小,進行試滑不離水試驗時,姿態比之前低很多,起飛滑跑階段僅出現輕微彈跳,無縱搖發生,驗證了優化結果是有效的。
利用VOF(Volume of Fluid)兩相流模型對水上電動飛機起飛滑行進行了數值仿真,模擬了滑行過程中氣液兩相流動分布狀態以及浮筒的噴濺狀態。對飛機各速度下浮筒相對于水線不同姿態下的阻力系數進行分析,并與試驗計算數據進行對比驗證,論證了數值仿真的可行性,也為雙浮筒式水上飛機的設計優化提供了借鑒。
以某型水上電動飛機為研究對象,對其進行實體1∶1建模,主要考慮飛機在水中滑行的動力特性,對機身的附件部分進行了簡化,專注對浮筒模型的細化。在導入Fluent時用毫米單位進行縮放。整機模型如圖8所示。

圖8 飛機模型Fig.8 Aircraft model
數值計算采用VOF方法。VOF方法通過對一種或多種相互不相容流體間的交界面進行追蹤,對于交界面不同的流體組分共用一套動量方程,計算時在整個流場的每個計算單元內引入水體積分數,通過求解體積分數的連續性方程來確定交界面的位置關系。
湍流模型選用Realizablek-e方程模型,與標準k-e模型相比,湍流黏性由定值改為變量,渦耗散率的輸運方程從精確的方程中推導得到,使方程能夠更加符合湍流的物理特性,可以更好地模擬強逆壓梯度下的邊界層流動以及流動分離。
采用耦合隱式算法對Navier-Stoke方程組進行聯立求解,由于采用隱式格式,因而計算精度與收斂性要優于耦合顯式方法,另一個突出的優點是可以求解全速度范圍,即求解范圍從低速流動到高速流動。采用基于網格中心最小二乘法構建流場梯度,壓力采用PRESTO!格式,動量方程采用二階迎風格式離散;體積分數項選用可壓縮方法,湍流耗散率和湍流動能采用一階逆風格式離散。在計算時,對于涉及表面張力的計算,選擇Body Force Formulation模式,可以保證壓力梯度和動量方程中表面張力的部分平衡,從而提高解的收斂性。
模型控制域依據流場選擇為長方體,飛機模型設置在控制域的中心位置,交界面以上為空氣域,以下為水域,控制域大小具體尺寸如下(L為模型長度):浮筒前部距入口邊界為2L;飛機尾部距離出口邊界為3L;交界面距上下邊界為2L;機身中線距離左右邊界為3L。
網格質量對流場的數值仿真結果十分重要,計算時選用兩套不同尺寸的網格(圖9)行對比,以飛機最小結構(0.4 m)尺寸進行遞減選擇;首先采用粗網格進行快速計算;再用更加精確的網格進行劃分,最小尺寸分別選為0.08 m和0.04 m,網格數量分別為83萬和140萬。
對于控制域采用結構化網格劃分,為更好地模擬自由液面處的浮筒附近水流噴濺狀態以及捕捉尾部流場情況,對浮筒部分及附近區域采用重疊網格進行加密(如圖10中綠色部分所示),利用重疊域保證整體計算域與模型相交部分網格的連續性,總網格數量為360萬。
設定氣液交界面類型為Interior,圖11為初始氣液分布示意圖;水域底部選為移動壁面邊界,壁面函數選擇Scalable Wall Function;飛機前部入口設定為速度入口邊界;出口設定為壓力出口;其余邊界設定為對稱邊界。

圖9 計算網格Fig.9 Computational grid

圖10 局部加密網格Fig.10 Locally encrypted grid

圖11 氣液交界示意圖Fig.11 Schematic diagram of gas-liquid junction
計算主要對整機的升阻力系數進行監控,得到水上飛機起飛滑行時阻力系數變化曲線和升阻比曲線(圖12和圖13)。兩組仿真模擬結果與后期實驗計算數據變化趨勢一致,都能夠基本描述飛機水面低速滑行階段的基本動力特征,對比粗網格與細網格模型來看,誤差主要存在兩個地方:前者出現“阻力峰”的速度向后推遲;越過“阻力峰”的后半階段阻力系數存在一定差異,細網格模型在各速度節點的阻力系數更大,總體來看細網格模型比粗網格模型的結果更優,數值吻合程度較好,最大誤差在15%之內,達到了較高的精度,驗證了仿真的可行性和可靠性。
對細網格模型的仿真結果進行了分析總結,主要以流場的噴濺狀態、底部相位分布探究了“阻力峰”出現的機制,并分析了浮筒結構的壓力變化分布,為浮筒內部加強結構的設計提供了依據。

圖12 阻力系數變化曲線Fig.12 Variation curves of drag coefficient

圖13 升阻比Fig.13 Lift-to-drag ratio
運用VOF方法對飛機模型起飛滑跑階段進行了數值模擬計算,對其過程進行了劃分,較好地模擬了自由液面對浮筒的噴濺狀態,對阻力變化特性進行分析總結,探究了“阻力峰”出現的機制,流場精確結果如圖14所示,圖中v為速度。
第1階段飛機受力比較簡單,主要重量由靜浮力支撐;隨著速度增大,縱向姿態角逐漸增大,從阻力曲線可以看出,這一階段的阻力快速增大,大部分來源于水動阻力(包括摩擦阻力、形狀阻力和興波阻力)和一部分空氣阻力。在此過程中浮筒升阻比迅速下降,是因為浮力基本不變,而浮筒受到的“吸力”隨速度的增大而增大。

圖14 流場噴濺示意圖Fig.14 Schematic diagrams of flow field splash
第2階段為過渡階段,從底部壓強分布(如圖15 所示)可以看出,阻力主要來源于浮筒前端受到的沖擊(紅色區域)以及浮筒底部和斷階處受負壓帶來的“吸力”(藍色區域)。隨著速度的提升,氣動升力也慢慢開始發揮主要作用,這是由于水對飛機的支撐逐漸由靜態浮力轉變為水動升力;在受力方面導致水動阻力迅速增大,底部負壓逐漸增大,范圍逐漸向浮筒后部轉移,從流場示意圖可以觀察到水流對浮筒噴濺較為嚴重,大量附著包圍在浮筒兩側;全機阻力在12 m/s左右到達“阻力峰”在越過最大阻力持續加速之后,縱向姿態開始小幅回正,全機總阻力開始減小。

圖15 浮筒底部壓強相位分布Fig.15 Pressure phase distribution at bottom of float
在整個滑水階段浮筒的壓力分布主要受到滑水速度和吃水深度的影響,在第1階段的較低速度下,浮筒開始排水,水流的噴濺會導致浮筒前端受到較大的撞擊,從壓力分布來看產生了局部高壓區域(如圖16所示),底部龍骨附近區域壓力值沿線性分布,逐步減小,整體受力比較均勻。

圖16 浮筒壓強分布Fig.16 Float pressure distribution
隨著速度的增大,浮筒底部在水流的作用下由靜態正壓力逐步轉化為負壓,且范圍不斷擴大,當速度達到10 m/s時,浮筒前半段基本覆蓋;而斷階處基本在滑行階段都處于負壓狀態,變化規律為先隨速度增大而增大,在接近駝峰速度時開始出現減小的趨勢,沿橫向來看壓力局部分布呈現為兩端大中間小的狀態,龍骨脊線處的壓強為兩端點壓強的76%。
加速至“駝峰速度”時,由于吃水深度的增加,浮筒后半段也逐漸受到水的吸力,也是此階段浮筒阻力的主要來源;但是由于斷階的存在,水動升力逐漸開始發生作用,斷階后部區域負壓峰值出現減小趨勢,部分區域的壓強出現由負轉正的現象。
通過對浮筒壓力分布進行分析,可以對后續浮筒的外形設計和內部加強結構的布置提供技術支持,例如在浮筒前端、斷階等主要受力點布置加強框。
電動飛機與常規油動飛機在動力性能方面的區別是在飛行中的輸出功率并不會隨速度的改變而改變[19]。因此需要對“阻力峰”狀態下的最大功率進行驗證計算,確認現有動力系統是否可以使飛機順利跨越“阻力峰”而正常起飛。利用仿真結果選取最大阻力狀態作為電動力系統設計的驗證點,起飛滑跑階段所需最大功率為[20]
(4)
式中:M為飛機總質量;g為重力加速度;V為出現“阻力峰”時飛機滑行速度;L/D為飛機升阻比。
電動力系統最大功率計算需要考慮飛機滑行起飛階段電推進系統各部件的效率:
(5)
式中:ηP為滑行階段螺旋槳效率;ηM為滑行狀態下電動機效率;ηC為滑行狀態下控制器效率;ηB為電池效率。
通過計算可知,克服“阻力峰”所需功率約為72 kW,現采用電動機最大功率為90 kW,滿足飛機正常起飛的性能要求。
安裝優化過的浮筒后,對水上飛機進行了實際試飛,試驗結果顯示縱向穩定性范圍變大,可以保持穩定滑行,未出現縱搖狀態,實現了對水上飛機起飛滑跑階段的一般性認識[21],主要考量了水上飛機滑行時縱向姿態角的變化情況,如圖17所示。

圖17 縱向姿態角Fig.17 Longitudinal attitude angle
可見,起飛過程與仿真結果基本吻合,根據速度變化劃分,滑行過程可以分為3個階段。
第1階段速度為0~8 m/s,為排水階段,水上飛機通過浮筒的作用將水推到兩側在水中移動,飛機滑行姿態穩定提升。
第2階段速度為8~20 m/s,為斷階滑水階段(圖18),隨著速度的增加,浮筒會由于斷階繞流產生的負壓而向下偏移,飛機縱向姿態角的不斷增大導致浮筒前體抬起而后體浸濕范圍增加,斷階位置成為主要滑水位置,飛機后的平行尾流越來越明顯,滑行過程速度超過8.3 m/s有困難,需要推桿降姿態以助于增速,接近12.0 m/s時達到最大姿態角度,后續自然增速至16.7 m/s以上,接近20.0 m/s時開始有縱搖趨勢。
第3階段為離水階段(圖19),機翼產生氣動力逐漸使飛機抬升高度,縱向姿態角穩定在8°,隨著在浮筒底部產生的“氣袋”不斷增大,空氣通量隨之增加,浮筒逐漸擺脫水對它的吸附力而離水,空氣不斷被吸入,并開始向尾部方向延伸,最終隨著速度的增大,速度達到22 m/s時浮筒與水面完全分離。

圖18 斷階滑行階段Fig.18 Step-off coasting stage

圖19 離水階段Fig.19 Water separation stage
1) 利用解析函數線性疊加法、拉丁超立方采樣、神經網絡模型和多島遺傳算法的組合策略,實現了對浮筒截面形狀的優化設計,浮筒水動性能有所提升,浮筒升阻比提高21.00%,為相關水下升力體的設計提供了經驗。
2) 基于Fluent中的VOF方法,對加裝浮筒的某型電動水上飛機起飛滑跑階段的力學特征進行了數值模擬計算,著重分析了不同速度下的姿態,阻力特性,全機阻力隨速度增大呈直線遞增,接近12.0 m/s時達到“阻力峰”,這有助于對水上飛機起飛滑跑過程進行一般性認識。
3) 數值模擬與試飛實驗的流場特征均有明顯尾流和噴濺,仿真結果與實驗數據相比,曲線變化趨勢基本吻合,最大誤差控制在15%之內;對“阻力峰”節點下電動力推進系統進行了驗證計算,所需最大功率滿足飛機起飛的性能要求,所得結論可為浮筒式電動水上飛機的研究設計提供借鑒。