跨域無人平臺水面垂直起飛動態特性數值模擬

2022-03-19 08:40:04邱磊鄭巢生

中國艦船研究 2022年1期

邱磊，鄭巢生

1 海裝重大專項裝備項目管理中心，北京 100071

2 中國船舶科學研究中心船舶振動噪聲重點實驗室，江蘇無錫 214082

0 引言

水下?空中跨域兩棲無人平臺結合了飛機、水面船和水下航行器的優點，靈活、機動且能自主運行，既可在水下潛航、協同作業并在必要時提供中繼服務，也可貼近海面飛行，執行快速搜索、運輸及中轉等任務[1]。

近年來，國內外越來越多的學者對跨域無人平臺進行了相關研究，同時，隨著計算機應用技術的發展，數值模擬作為一種有效的研究手段，也越來越多地被應用到跨域無人平臺研究領域。吝科等[2]采用數值模擬手段，針對升力型航行器開展了飛行時的空氣動力學和潛航時的水動力學特性研究。齊贊強[3]針對設計的一種新構型傾轉四旋翼無人機，開展了傾轉過渡狀態下螺旋槳對機身、機翼和尾翼的干擾數值模擬研究。吳江等[4]采用CFD軟件，對某型傾轉旋翼機的螺旋槳拉力及功率等性能進行了計算。鄧見等[5]采用數值模擬方法研究了仿飛魚跨介質無人平臺在水下游動和水面滑行階段的水動力性能，解釋了飛魚滑跑階段的合理性。廖保全等[6]采用Fluent軟件，針對一種可改變外形的水?空跨介質航行器進行了空中和水中2種外形的氣動力/水動力特性數值模擬，結果表明，通過改變外形航行器能夠同時滿足水下航行和空中飛行的要求。魏洪亮等[7]提出了水下發射航行體跨介質動態載荷預報方法，通過模擬出水過程中航行體逐步由水到氣的跨介質卸載作用，獲得了航行體軸向速度、俯仰角度、俯仰角速度和彎矩等關鍵參數隨時間的變化規律，模擬結果顯示與試驗數據吻合較好。譚駿怡等[8]針對一種可變體的水下構型，采用Fluent軟件模擬了該構型在典型工況下的出水過程，結果表明，航行器傾斜跨越水?空界面時兩側的流場和載荷會出現不對稱的劇烈變化，俯仰角越大，在水下的部分其流場受動范圍越小，跨越出水部分受影響越大；與零攻角出水相比，航行器帶攻角傾斜出水會導致表面所受流體作用力出現高頻率、大幅度的反向震蕩，影響出水的穩定性。杜特專等[9]針對跨介質航行器三角型截面舵翼，采用Fluent軟件對其運動過程中的空化流動與結構振動耦合進行模擬，分析了不同攻角下的空泡形態以及舵翼與空化流動的相互作用，結果表明，當來流攻角為2°～6°時，航行器及舵翼幾乎被包在超空泡內部；當來流攻角為8°時，舵翼的自由端會穿透空泡界面，使其所受水動力比小攻角條件下的大一個量級，振動特征也更加復雜。賈力平等[10]采用基于自由液面捕捉法的FINE/Marine軟件，模擬分析了DLR-F4飛機模型在水上迫降過程中的流場和水動力特性，并與已有實驗結果進行了對比，結果顯示該模擬方法可以很好地預測跨介質航行器的運動特性和水動力特性。譚駿怡等[11]針對一種特殊的雙層半環形閉合翼構型，采用Fluent軟件對該構型跨越水?空介質后變體過程中的不同狀態氣動力進行了數值模擬，結果表明，機翼回收過程中的氣動參數要大于定常狀態，展開過程中的氣動參數則小于定常狀態，且變體速率越大，非定常效應越明顯，而產生這種差異的原因來自于流場遲滯的影響。

綜上所述，目前跨域無人平臺的數值模擬研究大多聚焦于水下水動力或空中氣動力性能，較少考慮平臺自身在流體作用下的自由運動，且由于跨域無人平臺作業任務多樣，導致其運動模態也復雜多變，其中，作為多模態運動的重要環節——無人平臺從水面到空中的跨越過程，對于無人平臺的跨域功能性和安全性尤為重要，相關的跨域過程動態特性數值模擬研究相對較少。為此，本文擬采用黏流CFD方法結合重疊網格技術和多自由度DFBI運動模型，針對跨域無人平臺從水面垂直起飛至空中跨域過程的動態特性開展數值模擬研究，分析其運動及動力學特性，為后續跨域無人平臺優化設計及控制提供有力的評估手段。

1 研究對象及數值模擬方法

1.1 研究對象

本文選取一種跨域無人平臺作為研究對象，如圖1所示。該平臺總長L=3.0 m，翼展長B=2.8 m。無人平臺在空中飛行時自身重力通過兩側機翼產生的升力克服，飛行姿態通過垂直尾翼與水平尾翼的配合來控制，當無人平臺在水面起飛或降落時，推進器豎直向上產生升力，以克服無人平臺起飛或降落時的重力。本文主要開展無人平臺從水面垂直起飛至空中這一跨域過程的動態特性數值模擬研究。

圖1 一種跨域無人平臺Fig.1 A trans-media unmanned vehicle

1.2 數值方法

1.2.1 控制方程

本文中無人平臺在初始時刻平穩地浮于水面，此時需考慮自由液面的影響，可選用流體體積法（VOF）模擬自由面的水、氣兩相流，即流動中包含水、空氣2種流體，其中水的密度ρwater=998.2 kg/m3，空氣的密度ρair=1.225 kg/m3。定義空氣在計算網格上的空間占比為α，當α=1時網格內全部為空氣，α=0時網格內全部為水，當0<α<1時網格內既包含氣體，同時也包含水。

水、氣兩相流的密度 ρ和黏性 μ為空氣和水的線性組合：

為了簡化問題，鑒于流速比較低，假定空氣和水都為不可壓縮流體。流動的控制方程如下所示。

不可壓縮方程：

式中：V為速度矢量；p為流體壓力；g為重力加速度。

湍流模型選取SSTk-ω模型[12]，運動模型采用六自由度DFBI模型[13]。

1.2.2 數值求解方法

流動的數值求解使用有限體積法。連續性方程和動量守恒方程中的對流項離散使用二階迎風格式，擴散項使用二階中心差分格式，流場中的物理量梯度計算使用基于單元的Green-Gauss方法。離散方程求解利用SIMPLE方法和Gauss-Seidel迭代，同時，求解過程中使用多重網格技術加速迭代的收斂。壓力和速度計算松馳因子分別取0.2和0.5。計算采用隱式非定常方法，計算時間步長選取螺旋槳旋轉3°所用的時間。

1.2.3 初始邊界條件

無人平臺計算域如圖2所示。由于無人平臺為左右對稱布局，在不考慮左、右側向運動及橫滾的情況下，為了減少計算量，采用對稱邊界取一半進行計算。為研究無人平臺在水面起飛時的運動及動力特性，采用重疊網格技術結合多自由度DFBI模型方法，計算中不考慮風、浪、流對無人平臺起飛特性的影響，入口速度定義為0 m/s。自由面初始位置為排水體積為60 kg時的吃水深度。

圖2 數值模擬計算域Fig.2 Computational domain for numerical simulation

1.2.4 計算網格

采用切割體單元網格對計算域進行網格劃分，網格總數為432萬，其中背景網格160萬，重疊區域網格104萬，空氣螺旋槳旋轉區域網格168萬；對自由面、平臺表面、螺旋槳表面及附近，以及背景網格中原理樣機的運動路徑進行網格加密，如圖3和圖4所示。

圖3 計算域總體網格Fig.3 Total grids of computational domain

圖4 重疊區域網格Fig.4 Local grids of overset zone

2 數值模擬結果及分析

2.1 數值方法驗證

在開始跨域無人平臺水面垂直起飛動態特性模擬之前，需要開展數值方法驗證，特別是空氣螺旋槳氣動力數值計算方法的驗證。

空氣螺旋槳在市場上有多種產品可供選型，并且都有廠商提供的拉力試驗測試數據。本文首先選用了某型二葉商用空氣螺旋槳，如圖5所示，該螺旋槳槳葉直徑0.660 4 m。由于沒有該槳的詳細幾何型值，先對所選購的空氣螺旋槳進行了幾何逆向掃描，獲得了精確的槳葉幾何模型，然后，針對該幾何模型進行了數值建模和氣動力特性數值計算分析，其中空氣螺旋槳表面網格如圖6所示。計算中，采用的是SSTk-ω湍流模型結合非定?；凭W格方法。表1所示為廠商提供的二葉商用空氣螺旋槳拉力試驗數據和氣動力數值計算結果比較。

圖5 二葉商用空氣螺旋槳Fig.5 The two-bladed commercial air propeller

圖6 空氣螺旋槳氣動力計算網格Fig.6 The surface grids of air propeller for aerodynamic performance calculation

表1 某型二葉商用空氣螺旋槳性能試驗值與數值計算值Table 1 The experimental and numerical values of aerodynamic performance of a two-bladed commercial air propeller

從表1中可以看出，計算值與廠商數據吻合較好，說明目前的數值預報方法可以較好地模擬空氣螺旋槳的氣動力特性，通過數值計算，能夠作為空氣螺旋槳氣動力性能評估的依據。

2.2 計算工況

在開始水面起飛跨域過程的動態數值模擬之前，針對設計的三葉空氣螺旋槳進行了氣動力特性計算研究，其中槳葉直徑D=1.2 m。為了提供單槳294 N的拉力，使無人平臺以1 m/s的速度上升，計算得到空氣螺旋槳不同轉速下的拉力與扭矩如表2所示。

表2 不同轉速下單個三葉空氣螺旋槳的拉力、扭矩與功率Table 2 The pull force, torque and power of a three-bladed air propeller at different rotation speeds

為滿足294 N的拉力，通過插值計算，初步判斷三葉空氣螺旋槳槳葉直徑D=1.2 m時轉速需達到2 352 r/min。但這僅是針對單個空氣螺旋槳在敞開環境中的評估結果，沒有考慮空氣螺旋槳與無人平臺之間的相互干擾。因此，本文后續計算中，將考慮空氣螺旋槳與無人平臺之間的相互干擾，并放開無人平臺的運動自由度，模擬無人平臺從水面靜止狀態到起飛的動態過程。

2.3 計算結果與分析

計算時，首先選取空氣螺旋槳的轉速為2 352 r/min，針對無人平臺起飛過程中的運動特性與動力學特性進行計算，其自由度示意圖如圖7(a)所示，即放開上下、前后、俯仰3個自由度。計算坐標原點位于無人平臺首部，位移及速度計算參考點位于無人平臺重心處，如圖7(b)所示。

圖7 自由度方向及參考點示意圖Fig.7 The schematic diagram of direction of freedom degree and reference points

分析無人平臺開始從水面起飛的0.9 s內上升位移隨時間的變化情況，如圖8所示。從中可以看到，在螺旋槳轉速N=2 352 r/min時無人平臺開始是向上移動的，但約0.5 s后開始下降，說明在2 352 r/min轉速下空氣螺旋槳無法提供順利拉起無人平臺（一半為30 kg）所需的拉力。

圖8 N=2 352 r/min時的垂向位移分量Fig.8 The vertical displacement components at N=2 352 r/min

考察無人平臺開始從水面起飛的0.9 s內單槳拉力，如圖9所示。從中可以看到，在螺旋槳轉速N=2 352 r/min下空氣螺旋槳的拉力T基本穩定在219 N，確實無法順利拉起重量為30 kg的無人平臺（此處取一半計算）。

圖9 N=2 352 r/min時的單個螺旋槳拉力Fig.9 The single air propeller pull force at N=2 352 r/min

為了能夠順利從水面起飛，根據空氣螺旋槳提供的實際升力并結合表2中數據進行換算，預估空氣螺旋槳轉速至少達2 840 r/min時才能拉起無人平臺。因此，選取2 840 r/min作為空氣螺旋槳轉速，重新針對無人平臺起飛過程中的運動特性與動態特性進行計算，其中自由度示意圖仍如圖6所示，即放開上下、前后、俯仰3個自由度。

2.3.1 起飛特性

選擇無人平臺開始從水面起飛1.6 s的數據對起飛特性進行比較分析。分析起飛運動特性時，主要針對無人平臺上升的位移隨時間的變化，以及向上的速度分量隨時間的變化進行比較。

圖10所示為無人平臺開始從水面起飛運動時向上的位移分量。從中可以看到，當螺旋槳轉速N=2 840 r/min時無人平臺向上移動的速度更快，且被順利拉起升離水面。

圖10 N=2 840 r/min時的垂向位移分量Fig.10 The vertical displacement components at N=2 840 r/min

圖11所示為無人平臺開始從水面起飛運動時向上的位移速度。從中可以看到，當螺旋槳轉速N=2 840 r/min時無人平臺主要存在2個較明顯的加速區和1個減速區，其中0～0.2 s為第1個加速區，0.2～0.72 s為減速區，0.72 s之后為第2個加速區。

圖11 N=2 840 r/min時的垂向速度分量Fig.11 The vertical velocity components at N=2 840 r/min

分析起飛過程中空氣螺旋槳拉力與無人平臺在垂直方向的整體受力情況分別如圖12和圖13所示。其中，無人平臺的整體受力包含螺旋槳拉力、無人平臺重力和起飛過程氣動力。

圖12 N=2 840 r/min時的空氣螺旋槳拉力Fig.12 The air propeller pull force at N=2 840 r/min

圖13 N=2 840 r/min時的無人平臺垂向受力Fig.13 The vertical components of force on unmanned vehicle at N=2 840 r/min

由圖12和圖13可以看到，在第1個上升加速區，隨著無人平臺逐漸上升，浮力減小，整體垂向受力逐漸減小至0。在上升減速區，空氣螺旋槳拉力保持穩定，無人平臺整體垂向受力為負值，究其原因，是因為隨著上升速度的增加，無人平臺的上升阻力（氣動力）也隨之增加。在第2個上升加速區，受無人平臺與空氣螺旋槳氣動力耦合作用的影響，螺旋槳拉力開始波動，對比上升減速區瞬間（t=0.4 s）和第2個上升加速區瞬間（t=0.85 s），空氣螺旋槳的表面壓力系數分布如圖14和圖15所示。其中，壓力系數CP的定義為：

圖14 t=0.4 s時螺旋槳表面壓力系數分布Fig.14 The pressure coefficient distribution on air propeller surface at t=0.4 s

圖15 t=0.85 s時螺旋槳表面壓力系數分布Fig.15 The pressure coefficient distribution on air propeller surface at t=0.85 s

式中：P為當地壓力；P∞為遠前方壓力。

對比圖14和圖15可以看到，與上升減速區相比，在第2個上升加速區空氣螺旋槳的表面壓力存在明顯的非穩定低壓與高壓集中區，這也是空氣螺旋槳拉力劇烈波動的原因。

此外，對比上升減速區瞬間（t=0.4 s）和第2個上升加速區瞬間（t=0.85 s），得到無人平臺主體上表面的壓力系數分布分別如圖16和圖17所示。同樣，從圖中可以看到，與上升減速區相比，在第2個上升加速區無人平臺主體的上表面壓力明顯增加且更加不均勻，這也將導致無人平臺整體垂向受力波動明顯。

圖16 t=0.4 s時無人平臺上表面壓力系數分布Fig.16 The pressure coefficient distribution on the upper surface of unmanned vehicle at t=0.4 s

圖17 t=0.85 s時無人平臺上表面壓力系數分布Fig.17 The pressure coefficient distribution on the upper surface of unmanned vehicle at t=0.85 s

為分析無人平臺起飛過程中自由液面的形態以及周圍流場的變化，選取上升減速區瞬間（t=0.4 s）和第2個上升加速區瞬間（t=0.85 s）的自由液面波高分布及無人平臺下表面壓力系數分布情況，分別如圖18～圖21所示。從中可以看到，無人平臺從水面起飛出水過程中，自由液面的波高Hw分布不均勻性較明顯，但當無人平臺完全升離水面后，自由液面的波高Hw分布不均勻性明顯降低。分析圖19和圖21可知，自由液面變化對無人平臺下表面的壓力系數分布影響較小，這主要是因為無人平臺采用的是小水線面船體結構，下表面與水的接觸面較小。

圖18 t=0.4 s時自由液面波高分布Fig.18 The wave height distribution on free surface at t =0.4 s

圖19 t=0.4 s時無人平臺下表面壓力系數分布Fig.19 The pressure coefficient distribution on the lower surface of unmanned vehicle at t =0.4 s

圖20 t=0.85 s時自由液面波高分布Fig.20 The wave height distribution on free surface at t =0.85 s

圖21 t=0.85 s時無人平臺下表面壓力系數分布Fig.21 The pressure coefficient distribution on the lower surface of unmanned vehicle at t =0.85 s

2.3.2 平移特性

無人平臺從水面起飛時的水平位移分量如圖22所示。從圖中可以看到，當螺旋槳轉速N=2 840 r/min時，無人平臺一直保持向前運動。

圖22 N=2 840 r/min時的水平位移分量Fig.22 The horizontal displacement components at N=2 840 r/min

無人平臺從水面起飛時的水平速度分量如圖23所示。從圖中可以看到，當螺旋槳轉速N=2 840 r/min時，無人平臺主要存在2個較明顯的加速區，其中0.2～0.72 s為第1個加速區，0.72 s之后為第2個加速區，且第2個加速區的加速度更大。但與圖11所示的垂向速度相比，無人平臺的水平速度值偏小，達到一個數量級，即無人平臺的主運動為垂向上升運動。

圖23 N=2 840 r/min時的水平速度分量Fig.23 The horizontal velocity components at N=2 840 r/min

分析無人平臺在水平方向的受力如圖24所示。從中可以看到，對應第2個平移加速區的無人平臺水平受力均值大于第1個平移加速區。

圖24 N=2 840 r/min時平臺的水平受力Fig.24 The horizontal components of force on unmanned vehicle at N=2 840 r/min

2.3.3 俯仰及載荷特性

圖25所示為跨域無人平臺在起飛過程中俯仰角的變化曲線，該角度定義為原理樣機繞橫軸的旋轉角度，機頭抬高、機尾下降為正，反之為負。由圖中可以看到，在螺旋槳轉速N=2 840 r/min下，無人平臺在0.72 s之前俯仰角變化較小，在0.72 s之后俯仰角變化較大，存在快速低機頭現象。分析圖26所示的無人平臺俯仰力矩可知，在0.72 s之前俯仰力矩變化較小，在0.72 s之后俯仰力矩變化較大，且波動明顯。這表明原理樣機極有可能失穩，因為此時未加入控制程序，所以在該狀態下無人平臺最終會失穩，從而無法順利垂直起飛。

圖25 N=2 840 r/min時無人平臺的俯仰角度Fig.25 The pitch angle of unmanned vehicle at N=2 840 r/min

圖26 N=2 840 r/min時無人平臺的俯仰力矩Fig.26 The pitch moment of unmanned vehicle at N=2 840 r/min

圖27所示為螺旋槳轉速N=2 840 r/min時無人平臺快速低頭瞬間（t=1.6 s）的原理樣機表面的壓力分布。從中可以看到，空氣螺旋槳高速的下洗氣流沖擊在無人平臺上表面，使得以橫軸為界，至機頭的總體壓力大于至機尾的總體壓力，出現了將機頭向下壓的狀態。隨著無人平臺“低頭”，在總體下壓載荷作用下無人平臺的俯仰角變化呈指數形式上升，最終導致無人平臺失穩及墜落。因此，為保證無人平臺順利升空，在水面至空中的垂直起飛階段必須加入手動或自動的控制程序，以實時調整推進器傾轉角度。