基于Simulink的某型無人直升機艦面懸停特性分析

唐宏清，王華明

（南京航空航天大學航空學院，南京210016）

艦載直升機近艦面跟進懸停是其艦面起降的重要階段，在此期間其平衡特性受艦面流場環境影響顯著。目前常用的研究手段是：首先通過數值仿真確定艦面流場，其次采用葉素理論將直升機的主要氣動面（旋翼、尾槳）離散成若干微元，通過艦面流場與旋翼（尾槳）流場的耦合獲得微元氣動中心的總流場速度，進而確定旋翼（尾槳）氣動力，建立適用于艦面起降的直升機飛行動力學模型[1?5]。將所建模型作為仿真對象，進而可以研究無人直升機的艦面起降仿真及風限圖計算，對于設計自主起降控制率和保障艦面起降安全具有重要的工程價值。

目前Matlab、C、Fortran等程序設計語言仍是直升機飛行動力學建模的主要實現方式，由此建立的程序語言類模型不僅缺少層次感、可重用性，而且很難與后續的直升機飛行控制系統設計與驗證進行對接。Simulink是整合于Matlab中的系統建模和仿真專用平臺，借助其建模有如下特點及優勢：以塊作為基本建模單元，只需依據系統的描述形式（一般為微分方程）用信號線連接不同功能的塊，即可建立塊圖形式的仿真模型；其航空航天塊庫將航空航天領域系統建模與仿真的常用操作、計算程序封裝成一系列可直接使用的塊，顯著地提升了建模效率；其線性分析工具箱克服了牛頓迭代法[6]和優化求解法[7]的缺點，可直觀、高效地對模型進行配平計算。國內學者應用Simulink研究直升機飛行動力學模型的建立、配平和線化雖起步較晚，但仍取得了較大的成果。2006年，于志等借助Simulink平臺搭建了黑鷹（UH?60）直升機的飛行動力學模型，并通過模型配平獲得了該機的平衡特性，最后指出應用Simulink平臺建立直升機飛行動力學模型可以極大地提升模型的封裝性、層次性和可拓展性[8]；隨后，翁智勇也以UH?60直升機為研究對象，應用Simulink建立了該機的飛行動力學仿真模型，以研究直升機懸停和前飛時的操縱響應[9]。從公開發表的文獻來看，目前國內鮮有學者應用Simulink研究無人直升機的艦面起降飛行動力學建模及配平仿真。

本文采用CFD軟件計算SFS2艦在不同風速和風向下尾部飛行甲板上方的流場，在Simulink平臺中建立單旋翼帶尾槳式直升機的艦面起降飛行動力學仿真模型，并使用其專用配平工具——線性分析工具箱，對某型無人直升機懸停于載艦飛行甲板上方時的平衡特性進行仿真計算，以分析艦船尾流場風向、風速對該機近艦面懸停時的影響。

1 艦船流場計算

1.1 計算模型

為了模擬艦面的氣動環境，建立了護衛艦SFS2的外形數模（圖1），采用CFD軟件STAR?CCM+計算了來流風向角在-90°～90°（間隔15°）變化時，艦船飛行甲板上方的流場。SFS2艦長lsh=138.7 m，其他相關數據參考文獻[3]。

圖1 SFS2艦船外形數模Fig.1 Mathematic model of SFS2 ship contour

1.2 邊界條件與網格生成

流場計算域采用長方體（長為10lsh、寬為7lsh、高為5lsh），以簡化后續的網格劃分及入口風向設置。計算域的邊界條件為：速度入口、壓力出口、滑移壁面和無滑移壁面。為了使生成的流體網格大部分為計算性能較好的六面體網格，本文采用Star?CCM+中的切割體網格生成器進行網格劃分。同時，為了能捕捉到飛行甲板和艦船周圍的流場信息，本文使用混合尺寸網格劃分策略，在計算域內部從里到外設置了3個長方體網格加密區，以實現對飛行甲板和艦船周圍的網格密度控制，最終生成的流體網格數量約為343萬個（圖2）。另外，針對艦船尾流場雷諾數高和極不穩定的特點，本文采用K?Epsilon模型模擬尾流場的湍流分量。

圖2 計算域網格劃分Fig.2 Grid division of computing domain

1.3 流場計算結果

圖3展示了來流風向角0°、風速Vw=10 m/s時，SFS2艦船縱向對稱面內的X方向速度云圖。從圖3中可知，由于該艦船高層建筑的阻擋作用，導致了氣流分離的產生，形成了4個大尺度的渦流區（圖3中青色區域）。

圖3 艦船縱向對稱面內X方向速度云圖Fig.3 Velocity cloud chart in X?direction in longitudinal symmetry plane of ship

圖4～5展示了艦船機庫后方15 m處橫截面內X、Z方向的速度分布云圖。由圖4～5可知，X、Z方向的速度基本呈現對稱分布，與實際情況基本一致：機庫后方形成了渦流區（圖4中顏色由青色逐漸變為藍色），在甲板正上方產生了較大的下沖氣流（圖5中藍色區域）。

圖4 機庫后方15 m處橫截面內X方向速度云圖Fig.4 Velocity cloud in X?direction in cross section 15 m behind the hangar

圖5 機庫后方15 m處橫截面內Z方向速度云圖Fig.5 Velocity cloud in Z?direction in cross section 15 m behind the hangar

為了驗證本文計算結果的正確性，選取機庫后方6.858 m處橫截面（25%甲板長度）的X方向速度計算結果與文獻[10]的試驗結果進行對比。圖6為該截面內高10.668 m（35 ft）處的對比結果，從圖6中可以看出，本文的計算結果與試驗結果吻合良好。

圖6 25%甲板長度X方向速度計算結果對比Fig.6 Comparison of calculation results of X?direction ve?locity at 25%deck length

1.4 流場數據處理與應用

由于混合尺寸網格的使用，Star?CCM+導出的流場速度分量為流體域內各個非均勻網格節點的計算結果。為了簡化流場數據使用，精簡數據容量，提高計算效率，本文使用Tecplot軟件的線性插值功能實現流場速度分量從非均勻網格節點至均勻網格節點的映射，最終完成了艦船尾流場數據庫的建立。

為了驗證均勻映射后流場的正確性，選取機庫后方15 m處橫截面內高5 m處的位置，應用Tec?plot軟件對Star?CCM+計算的非均勻網格節點速度與映射后的均勻網格節點速度(VX、VZ)進行了對比，如圖7所示，可以看出兩者十分吻合。

圖7 X方向和Z方向節點速度對比Fig.7 Comparison of node specd in X-and Z-directions

當直升機在艦船附近飛行時，根據實時計算出的氣動中心（如槳葉單元中心、尾槳槳轂中心等）位置矢量進行插值，便可從流場數據庫中獲得氣動中心處的流場速度。將插值獲得的流場速度疊加至氣動中心的原速度上，從而將艦船尾流場耦合至直升機艦面起降飛行動力學模型。

2 直升機艦面起降飛行動力學建模

2.1 坐標系

圖8展示了本文建模使用的部分坐標系及其位置。艦船坐標系（下標sh）的原點與載艦機庫門下方的中點重合，Xsh軸為艦船縱向對稱軸，指向船尾為正，Zsh軸與Xsh軸垂直且正向朝上，Ysh軸依據右手法則確定。機體坐標系（下標B）的原點與機體重心重合，XB軸沿機身縱向對稱軸且正向朝前，ZB軸與XB軸垂直且正向朝下，YB軸由右手法則確定。地面坐標系（下標E）的XE軸指向直升機初始航向，ZE軸指向地心，YE軸由右手定則確定，原點與艦船坐標系的原點重合，因為本文暫不考慮艦船的縱搖、橫搖和沉浮運動，且將艦船的航行速度等價為艦船靜止而氣流沿反方向從船前吹來。飛行動力學建模涉及的其他坐標系參考文獻[11]。

圖8 建模坐標系Fig.8 Coordinate systems of modeling

2.2 Simulink建模

本文利用Simulink平臺基于組件的建模特點，同時結合單旋翼帶尾槳式直升機飛行動力學建模的理論方法[6,12?13]，將整個模型分成旋翼、尾槳、機身、平尾和垂尾5個頂層組件，先用Simulink分別建立各組件的氣動力計算模型并進行驗證，然后組裝成直升機艦面起降飛行動力學模型。

旋翼的氣動力與誘導速度、槳葉揮舞運動和艦面流場密切相關，本文采用動力入流方程[14?16]計算旋翼的誘導速度，采用錐度角、后倒角和側倒角描述槳葉的揮舞運動[17]，利用Matlab Function塊建立旋翼誘導速度模型和槳葉揮舞運動模型。為了計入艦船尾流對旋翼氣動載荷的影響，將槳盤沿周向和徑向離散成若干葉素，先計算各葉素相對地面坐標系的位置矢量，然后從數據庫中提取艦船尾流信息，通過三維插值獲得各葉素處的艦面氣流速度，將其與旋翼流場合成后計算各葉素氣動力，最終獲得整個旋翼的氣動力。為了便于定義循環次數和提取每次循環的計算結果，本文使用For Each Subsystem塊實現單個仿真時間步上的葉素方位角循環計算。

尾槳氣動力計算原理與旋翼類似，可以通過對旋翼組件的重用完成尾槳組件的建模。機身、平尾和垂尾組件通過各自的實時流場信息及氣動力系數插值表實現氣動力的計算。在計入艦船尾流場對尾槳、機身等組件的氣動力影響時，首先根據各組件氣動中心相對地面坐標系的位置矢量進行插值獲得氣動中心處的艦面氣流速度，然后將其與氣動中心的原速度疊加獲得合氣流速度，進而根據合氣流速度計算各組件的氣動力。

將直升機近艦面飛行視為剛體運動，應用Simulink航空航天塊庫中的6?DOF(Euler angles)塊構建直升機運動微分方程（如式（1～2）所示），輸入量為直升機在機體坐標系下的合力及合力矩，輸出量為直升機的位置、姿態、速度和加速度，即有

式中：VB=(u,v,w)T為直升機的運動線速度，ωB=(p,q,r)T為角速度，mh為直升機的質量，Ih為直升機的慣量矩陣，FsumB=(Fx,Fy,Fz)T為機體坐標下直升機各組件氣動力與重力的總和，MsumB=(Mx,My,Mz)T為機體坐標系下各組件氣動力矩的總和。

最終建立的直升機艦面起降飛行動力學Simulink模型如圖9所示。

圖9中旋翼模型（Main_Rotor）的內部結構如圖10所示，圖10中cal_elem_force_moment子系統即為一個For Each Subsystem塊。

圖9 直升機艦面起降飛行動力學模型Simulink塊圖Fig.9 Simulink block diagram of flight dynamics model of helicopter shipboard operations

圖10 旋翼模型Simulink塊圖Fig.10 Simulink block diagram of main rotor model

3 直升機艦面懸停平衡特性仿真

3.1 算例直升機與流場環境

本文以某型單旋翼帶尾槳式的無人直升機為算例，以研究該機的艦面懸停平衡特性。算例無人直升機旋翼構型為蹺蹺板式，總質量為420 kg，旋翼和尾槳主要參數如表1所示。

表1 旋翼和尾槳的主要參數Table 1 Main parameter s of r otor and tail rotor

本文假設近艦面懸停時，算例無人直升機重心在SFS2艦船坐標系中的位置坐標為(15 m,0 m,4 m)，以該點為中心，并結合算例直升機的全機尺寸，從SFS2艦船的流場數據庫中選取一個長方體區域作為配平計算的艦面流場環境。

3.2 基于線性分析工具箱的配平計算

線性分析工具箱是Simulink專門用于配平、線化動力學模型的工具。相比于常用的配平方法，線性分析工具箱的動力學模型配平方法不僅形象簡單、求解速度快，而且對操作點的初值要求極低。因此，本文使用Simulink線性分析工具箱作為直升機飛行動力學模型的配平工具。

應用該工具箱配平主要包含兩個步驟：一是約束操作點，即對模型狀態量、根輸入和根輸出的屬性添加適當約束；二是優化求解，即利用解算函數findop計算出給定約束下的平衡操作點。其配平的關鍵和難點在于根據直升機的飛行狀態對Simulink模型施加適當合理的操作點約束，而且約束不能重復。

3.3 艦面懸停平衡特性分析

將建立的艦船尾流場數據庫作為飛行環境，應用線性分析工具箱對算例直升機展開近艦面懸停配平仿真，獲得了相應的平衡特性如圖11～16所示。圖11展示了算例直升機旋翼總距在不同風向、風速下的平衡計算結果。圖11顯示結果曲線具有對稱性的特點。當甲板上方以正向來流為主時（即較小的風向角），隨著風速的增加，總距逐漸減小；而在大側風下，隨著風速的增大，旋翼總距先減小后增大，其原因是該風況下艦尾飛行甲板上方存在較大的側風分量，直升機的廢阻顯著增大，進而引起了旋翼總距增大。在0°風向角時，旋翼總距基本不變是圖5所示的下沖氣流引起的，即甲板上方的下沖氣流抵消了槳盤平面流量增大造成的總距減小。

圖11 旋翼總距隨風向和風速變化曲線Fig.11 Variation curves of rotor collective pitch with wind direction and wind speed

圖12～13展示了算例直升機橫向平衡特性隨風向、風速的變化。由圖12可知，橫向周期變距和滾轉角隨風向、風速的變化趨勢，呈現出反對稱的特點。當甲板上方以右側風為主時（即風向角小于0°），直升機的橫向周期變距、滾轉角隨著風速的增大均增大（分別對應右壓桿、右滾）；而左側風為主時（風向角大于0°），則相反。上述現象與實際飛行一致，即右側風吹向甲板可以視為大氣靜止而直升機往右側飛行，因而需要右壓桿操縱，產生右滾姿態。圖13顯示甲板側風對算例直升機的滾轉姿態影響顯著，原因是該風況下的機身、垂尾出現了較大的橫向阻力，故需要重力提供反向的橫向力分量以維持橫向平衡，從而形成了較大的滾轉姿態。但由于尾槳拉力（向右）也參與直升機的橫向平衡，所以大風向、風速下的直升機左滾姿態比右滾姿態更大。

圖12 橫向周期變距隨風向和風速變化曲線Fig.12 Variation curves of lateral cyclic pitch with winddirection and wind speed

圖13 滾轉角隨風向和風速變化曲線Fig.13 Variation curves of roll angle with wind direction and wind speed

圖14～15展示了算例直升機縱向平衡特性隨風向、風速的變化。由圖14可知，縱向周期變距、俯仰角隨風向及風速的變化趨勢，呈現出一定的對稱性。當甲板上方以正向來流為主時（即較小的風向角），直升機的縱向周期變距、俯仰角隨風向、風速的變化不大，因為此時艦船高層建筑的阻擋減弱了抵達艦面懸停處的正向來流速度，故直升機的縱向平衡特性變化較小。大側風下的俯仰角變化是由該風況下機身復雜的氣動特性引起的。算例直升機機身在大側風下提供了較大的前向力，故需重力向后傾斜產生反向分量以維持縱向力平衡，因而直升機產生了較大的抬頭姿態。

圖14 縱向周期變距隨風向和風速變化曲線Fig.14 Variation curves of longitudinal cyclic pitch angle with wind direction and wind speed

圖15 俯仰角隨風向和風速變化曲線Fig.15 Variation curves of pitch angle with wind direction and wind speed

圖16展示了算例直升機尾槳總距在不同風向、風速下的平衡計算結果。由圖16可知，結果曲線也呈現出反對稱的特點。右側風時，尾槳總距隨著風速的增大逐漸增大；而左側風時，則相反。上述現象的出現與側風吹來時尾槳拉力的變化密不可分。當甲板上方以右側風為主時，尾槳槳葉剖面迎角的減小致使其拉力不足以維持航向力矩平衡，故需增大尾槳總距；再者，尾槳拉力還需平衡機身、垂尾產生的較大偏航力矩（向右），進而需進一步地增大尾槳總距。

圖16 尾槳總距隨風向和風速變化曲線Fig.16 Variation curves of tail rotor collective pitch with wind direction and wind speed

4 結 論

本文利用Simulink平臺建立了耦合艦面流場的直升機艦面起降飛行動力學仿真模型，并使用該平臺的線性分析工具箱對算例無人直升機進行了配平仿真，得到了該機近艦面懸停時的平衡特性：

（1）較小的來流風向、風速時，直升機氣動面受正向來流影響不大，致使其操縱量與姿態角隨風向、風速的變化較小，但旋翼總距受0°風向角時甲板下沖氣流的影響較大，導致其隨風速的增大基本不變。

（2）隨著風向、風速的增大，尾流側向速度分量對旋翼、機身和垂尾的氣動力、力矩產生了較大影響，直升機操縱量、姿態角隨之產生了較大的變化，特別是尾槳操縱與滾轉角。

（3）直升機橫向周期變距、尾槳總距及滾轉角隨著來流風向、風速的變化趨勢，呈現出反對稱的特點，而其余操縱及姿態的變化趨勢則表現出一定的對稱性。

由此可見，本文建立的艦面起降飛行動力學Simulink模型，能夠較好地捕捉到算例無人直升機近艦面懸停時艦面流場對其平衡特性的影響，可以用于開展艦面起降仿真及風限圖計算相關的研究。