安朝，霍貴璽，孟楊，謝長川，楊超

北京航空航天大學(xué) 航空科學(xué)與工程學(xué)院，北京 100191

無人機(jī)（Unmanned Aerial Vehicle，UAV）起源于工業(yè)時(shí)代，隨著航空需求的多樣化發(fā)展，無人機(jī)技術(shù)不斷進(jìn)步，逐步成為航空器的重要分支。相比有人機(jī)，無人機(jī)具有成本低、人員零傷亡、持續(xù)作戰(zhàn)能力強(qiáng)等優(yōu)勢，在軍事偵察、快速打擊、物資運(yùn)輸及通訊中繼等方面發(fā)揮重要作用，受到各國人員的廣泛關(guān)注［1-2］。近年來，輕質(zhì)材料結(jié)構(gòu)、計(jì)算機(jī)及控制技術(shù)的飛速發(fā)展為無人機(jī)在民用及軍事中的廣泛應(yīng)用創(chuàng)造了條件。高空長航時(shí)無人機(jī)（High-Altitude Long-Endurance，HALE）、無人機(jī)集群是頗具潛力的主要發(fā)展方向。高空長航時(shí)無人機(jī)以臨近空間太陽能無人機(jī)為代表，飛行高度高、續(xù)航能力強(qiáng)、任務(wù)范圍大、生存能力強(qiáng)，但靈活性較差，搭載的任務(wù)載荷類型單一。特別是高空長航時(shí)無人機(jī)存在幾何非線性氣動(dòng)彈性問題［3-5］，顫振速度低，分析難度高。無人機(jī)集群載荷種類多，智能化程度高，分布式任務(wù)能力強(qiáng)，能夠完成大量單機(jī)系統(tǒng)無法完成或者由群體更有效完成的任務(wù)。近年來無人機(jī)集群技術(shù)不斷進(jìn)步，正在向去中心化、多任務(wù)、低成本方向發(fā)展，但由于體積受限和機(jī)動(dòng)靈活的控制需要，無人機(jī)集群存在航時(shí)短、巡航高度低的不足?？紤]到任務(wù)需求及復(fù)合平臺(tái)化發(fā)展的方向，融合高空長航時(shí)無人機(jī)及無人機(jī)集群的優(yōu)勢是新概念無人機(jī)可預(yù)見的發(fā)展方向之一。

多體組合式無人機(jī)（Multi-body Combined Aircraft，MCA）是一類由多個(gè)單體固定翼無人機(jī)在翼尖或其他位置組合形成的新概念飛行器，組合形式包括但不限于平行翼尖并聯(lián)組合、平面式網(wǎng)狀組合等，連接形式包括但不限于剛性連接、鉸接或柔性連接等。多體組合式無人機(jī)概念中，以翼尖位置處鉸接形成的多體組合式無人機(jī)兼具高空長航時(shí)無人機(jī)和固定翼無人機(jī)集群的優(yōu)勢，具有極大發(fā)展?jié)摿Α?/p>

組合式無人機(jī)概念最早出現(xiàn)于1931 年，美國海軍基于F9C-2“獵鷹”式戰(zhàn)斗偵察機(jī)，利用機(jī)翼掛鉤與飛艇組合成可分離回收組合體。20 世紀(jì)40 年代，美國空軍嘗試依托固定翼飛機(jī)機(jī)翼翼尖附加存放燃油的飛行模塊來增加飛機(jī)航程。隨后直至1956 年，美國空軍多次開展運(yùn)輸機(jī)與小型靶機(jī)、戰(zhàn)斗機(jī)的空中對接、組合飛行試驗(yàn)，對接形式為翼尖固接。試驗(yàn)中發(fā)現(xiàn)組合飛行后飛行器產(chǎn)生復(fù)雜的“拍打運(yùn)動(dòng)”導(dǎo)致難以完成飛行控制。在發(fā)生2 次飛行事故后，研究中止［6-8］。

21 世紀(jì)后，組合飛行概念再次被各國研究人員關(guān)注，其應(yīng)用目標(biāo)由翼尖拖帶擴(kuò)展段轉(zhuǎn)變至同構(gòu)飛行器組合提高航程航時(shí)與多任務(wù)處理能力。2002 年，Magill 和Durham［9］在風(fēng)洞試驗(yàn)中采集數(shù)據(jù)，對翼尖對接飛行、近距離編隊(duì)飛行和攜帶飛行等飛行模式進(jìn)行了研究。2014—2016 年，Montalvo 等［10-12］對翼尖相連及首尾相連的飛行器進(jìn)行分析研究，提出“元飛行器”概念，基于升力線方法仿真發(fā)現(xiàn)飛機(jī)連接時(shí)，飛行動(dòng)力學(xué)特性如短周期模態(tài)、荷蘭滾模態(tài)等隨連接的數(shù)量和連接方式變化而變化。K?the 和Luckner［13］基于Kane 方法及片條理論進(jìn)行了翼尖鉸接組合式無人機(jī)動(dòng)力學(xué)建模與穩(wěn)定性分析，發(fā)現(xiàn)鉸接連接時(shí)飛行系統(tǒng)先天帶有不穩(wěn)定飛行模態(tài)。2017—2020 年，K?the 等［14］和Cracau［15］基于零極點(diǎn)配置法設(shè)計(jì)飛行控制系統(tǒng)并結(jié)合常規(guī)平直翼構(gòu)型無人機(jī)進(jìn)行了飛行試驗(yàn)驗(yàn)證。Cooper和Rothhaar［16］基于翼尖渦模型及通用非線性氣動(dòng)模型給出了多架無人機(jī)空中組合對接過程的動(dòng)力學(xué)模型及控制方案。安朝等［17］基于片條理論及Newton-Euler 方程建立多機(jī)組合的飛行力學(xué)模型，指出了該類組合式飛行器存在相對滾轉(zhuǎn)不穩(wěn)定飛行模態(tài)并設(shè)計(jì)飛行控制律。楊延平等［18］在綜述文章中介紹了集群組合式柔性無人機(jī)概念，從總體、氣動(dòng)、結(jié)構(gòu)等6 個(gè)方面梳理和分析了此類飛行器亟待解決的關(guān)鍵技術(shù)，具體包括合作博弈的組合體無人機(jī)總體設(shè)計(jì)、柔性多機(jī)組合非定常數(shù)值仿真及對接技術(shù)、多智能體的協(xié)同集群運(yùn)動(dòng)控制技術(shù)等。武建明等［19］探究了翼尖連接對于太陽能飛機(jī)降低能耗和增加飛機(jī)續(xù)航能力的影響，并導(dǎo)出了翼尖連接飛機(jī)的Osward 因子、誘導(dǎo)阻力系數(shù)及功率消耗模型。馬培洋等［20］首先基于升力線理論分析了鏈?zhǔn)浇M合式飛行器氣動(dòng)特性，梳理了其在氣動(dòng)性能、續(xù)航能力、巡航高度和巡航速度上的優(yōu)勢，進(jìn)行了初步飛行試驗(yàn)驗(yàn)證。周洲等［21］采用計(jì)算流體力學(xué)方法建立氣動(dòng)數(shù)據(jù)庫，基于擬坐標(biāo)形式的Lagrange 方程建立動(dòng)力學(xué)模型并分析了飛行力學(xué)穩(wěn)定性問題。

多體組合式無人機(jī)在動(dòng)力學(xué)特性分析、飛行控制系統(tǒng)設(shè)計(jì)、結(jié)構(gòu)形式及應(yīng)用方式等與以往無人機(jī)系統(tǒng)設(shè)計(jì)有明顯不同。當(dāng)前多體組合式無人機(jī)設(shè)計(jì)研究側(cè)重于動(dòng)力學(xué)特性建模與初步試驗(yàn)，鮮有氣動(dòng)布局參數(shù)影響研究。此類新概念飛行器在飛行動(dòng)力學(xué)特性特別是穩(wěn)定性方面與傳統(tǒng)飛行器有顯著區(qū)別，氣動(dòng)布局設(shè)計(jì)過程中要考慮其飛行動(dòng)力學(xué)特殊性。本文針對翼尖鉸接多體組合式無人機(jī)構(gòu)型，基于狀態(tài)空間渦格法提出一種考慮多機(jī)氣動(dòng)耦合效應(yīng)的氣動(dòng)導(dǎo)數(shù)計(jì)算方法，結(jié)合飛行動(dòng)力學(xué)模型開展布局參數(shù)影響研究，重點(diǎn)關(guān)注布局參數(shù)對飛行力學(xué)穩(wěn)定性影響，進(jìn)行不同氣動(dòng)布局參數(shù)下的氣動(dòng)特性與飛行力學(xué)穩(wěn)定性對比分析。翼尖鉸接多體組合式無人機(jī)概念如圖1 所示。

圖1 翼尖鉸接多體組合式無人機(jī)概念Fig.1 Concept of MCA with hinge joint wing-tip

1 氣動(dòng)建模方法

1.1 狀態(tài)空間渦格法

基于狀態(tài)空間渦格法［22-24］建立考慮氣動(dòng)耦合效應(yīng)的多體組合式無人機(jī)氣動(dòng)特性分析與氣動(dòng)導(dǎo)數(shù)計(jì)算方法，重點(diǎn)關(guān)注多體組合式無人機(jī)組合狀態(tài)下，某單體無人機(jī)運(yùn)動(dòng)參數(shù)的變化對組合式無人機(jī)整體及其他單體無人機(jī)氣動(dòng)特性的影響?；趧萘骼碚摰姆嵌ǔu格法在飛行器氣動(dòng)特性分析及氣動(dòng)彈性計(jì)算領(lǐng)域得到了廣泛的應(yīng)用。非定常渦格法相較于片條理論、計(jì)算流體力學(xué)方法及氣動(dòng)數(shù)據(jù)庫方法，兼顧計(jì)算精度與求解效率優(yōu)勢。

定義氣動(dòng)求解坐標(biāo)系如下：x軸沿來流方向，y軸水平向右，z軸由右手定則確定。將氣動(dòng)面如機(jī)翼中弧面等沿弦向和展向劃分為若干四邊形氣動(dòng)網(wǎng)格，中弧面形狀可以考慮機(jī)翼彎度及預(yù)扭轉(zhuǎn)等因素影響。在每個(gè)中弧面氣動(dòng)網(wǎng)格的1/4 弦線處布置渦線段，四段等強(qiáng)度渦線首尾相連構(gòu)成一個(gè)渦環(huán)，所有中弧面氣動(dòng)網(wǎng)格上的渦環(huán)構(gòu)成物面渦。在機(jī)翼中弧面后方布置渦線段構(gòu)成尾渦。選取中弧面氣動(dòng)網(wǎng)格3/4 弦線中點(diǎn)為控制點(diǎn)，物面渦及尾渦產(chǎn)生的誘導(dǎo)速度在控制點(diǎn)處須滿足Neumann 邊界條件。物面渦及尾渦分布示意如圖2 所示［24］。

圖2 渦分布示意圖［24］Fig.2 Illustration of vortex distribution［24］

設(shè)物面渦個(gè)數(shù)為M，強(qiáng)度列向量Γb維數(shù)為M×1，氣動(dòng)面展向網(wǎng)格數(shù)為Ns，弦向網(wǎng)格數(shù)為Nc，最靠近氣動(dòng)面后緣的一排尾渦強(qiáng)度列向量Γw0維數(shù)為Ns×1，其他尾渦網(wǎng)格數(shù)為Nwl，其他尾渦強(qiáng)度列向量Γwl維數(shù)為Nwl×1。設(shè)氣動(dòng)面控制點(diǎn)處法向下洗速度為M×1 維列向量w，結(jié)合Biot-Savart 定律，控制點(diǎn)處Neumann 邊界條件可寫為

式中：w中各分量wi=為來流速度列向量，n0i為第i個(gè)氣動(dòng)網(wǎng)格中控制點(diǎn)處法向量；Kb、Kw0、Kwl分別為附著渦、第1 排尾渦及其他尾渦的誘導(dǎo)系數(shù)矩陣。

根據(jù)后緣Kutta 條件，后緣點(diǎn)環(huán)量為0，脫尾渦的過程保持渦強(qiáng)守恒，最后一排物面渦渦強(qiáng)與第1 排尾渦渦強(qiáng)相等［25］：

式中：C1為含0 和1 的常數(shù)提取矩陣；Δt為間隔時(shí)間。由式（2）可得：

根據(jù)Kelvin 定理進(jìn)一步規(guī)定尾渦脫出后保持渦強(qiáng)不變：

式中：C2、C3同樣為表征渦對應(yīng)關(guān)系的常數(shù)提取矩陣，只包含0 和1。

綜合方程（1）、（3）、（4）可得渦格法狀態(tài)空間方程形式為

式中：尾渦強(qiáng)度列向量Γw=[Γw0Γwl]T；Aa、Ba為狀態(tài)空間系數(shù)矩陣，矩陣內(nèi)元素值僅與氣動(dòng)面幾何形狀及渦的劃分有關(guān)：

根據(jù)Kutta-Joukowski 定理可推導(dǎo)得到氣動(dòng)力輸出方程，將氣動(dòng)力分為定常部分Fst及非定常部分Funst。氣動(dòng)力定常部分Fst作用于渦環(huán)中點(diǎn)，可表示為

式中：eΓ為前緣切向渦方向矢量；Γ為渦量數(shù)值；ρ為來流空氣密度。氣動(dòng)力非定常部分Funst可表示為

式中：fa為氣動(dòng)力列向量；Ca、Da為輸出方程系數(shù)矩陣。

綜上，方程（5）、（10）構(gòu)成狀態(tài)空間渦格法狀態(tài)空間模型，狀態(tài)變量為尾渦強(qiáng)度列向量Γw，輸入變量為氣動(dòng)網(wǎng)格控制點(diǎn)處法洗速度w及其導(dǎo)數(shù)列向量，輸出變量為fa。

1.2 氣動(dòng)導(dǎo)數(shù)計(jì)算方法

導(dǎo)數(shù)氣動(dòng)力模型在飛行動(dòng)力學(xué)中應(yīng)用廣泛，其表達(dá)簡潔便于動(dòng)力學(xué)仿真及控制系統(tǒng)設(shè)計(jì)。傳統(tǒng)的飛行動(dòng)力學(xué)導(dǎo)數(shù)氣動(dòng)力模型往往基于試驗(yàn)數(shù)據(jù)獲得，成本較高，渦格法計(jì)算精度滿足要求同時(shí)計(jì)算效率高，在獲得渦格法狀態(tài)空間形式后，氣動(dòng)力具有解析表達(dá)式，直接推導(dǎo)氣動(dòng)導(dǎo)數(shù)表達(dá)的解析形式成為求解氣動(dòng)導(dǎo)數(shù)的新思路。同時(shí)，多體組合式無人機(jī)存在多機(jī)之間的氣動(dòng)耦合，組合后任一單體無人機(jī)氣動(dòng)導(dǎo)數(shù)受到其他組合無人機(jī)狀態(tài)變量的影響，解析形式的氣動(dòng)導(dǎo)數(shù)求解方法在考慮多機(jī)氣動(dòng)耦合特性方面更加便利。

考慮渦格法狀態(tài)空間形式中的定常部分，控制方程可表達(dá)為

對于剛體飛機(jī)，氣動(dòng)網(wǎng)格控制點(diǎn)處法向量n0i在機(jī)體坐標(biāo)系中的投影保持不變。飛行器某點(diǎn)運(yùn)動(dòng)速度v可表示為［22］

式中：V0為飛行速度矢量；ω為飛行角速度矢量；r0為控制點(diǎn)至質(zhì)心矢徑。對于第i個(gè)氣動(dòng)網(wǎng)格控制點(diǎn)，法洗速度為

投影到體軸系上的剛體氣動(dòng)力列向量fAr為

式中：G為投影矩陣。全機(jī)氣動(dòng)載荷的合外力列向量FAr與合外力矩列向量MAr可通過結(jié)點(diǎn)氣動(dòng)載荷求和得到：

式中：求和矩陣Φt、Φr具有如下形式：

其中：xi、yi、zi為第i個(gè)氣動(dòng)網(wǎng)格控制點(diǎn)在體軸系中的坐標(biāo)。在渦格法剛體氣動(dòng)力表達(dá)式（14）中，系數(shù)矩陣Aa、Ba、Ca、Da僅與氣動(dòng)面幾何形狀及渦的劃分有關(guān)，而與飛機(jī)運(yùn)動(dòng)速度相關(guān)的變量及配平變量等都包含在法洗速度表達(dá)式（13）中。將剛體運(yùn)動(dòng)角速度列向量ω=[p q r]T及速度列向量V0=[-V0-V0β-V0α]T代入式（14）中，可求得剛體氣動(dòng)力FAr和氣動(dòng)力矩MAr對配平變量的各項(xiàng)導(dǎo)數(shù)如下：

其中：α、β、p、q、r分別為飛行器攻角、側(cè)滑角、滾轉(zhuǎn)角速度、俯仰角速度、偏航角速度；cp、cq、cr為列向量，各分量形式為

當(dāng)n個(gè)單體無人機(jī)進(jìn)行組合后，令整體氣動(dòng)面區(qū)域?yàn)棣福髯訜o人機(jī)對應(yīng)氣動(dòng)面區(qū)域分別為Ωk，k=1，2，…，n，設(shè)求和矩陣Φtk、Φrk為

式中：I3×3及所在位置對應(yīng)Ωi在整體氣動(dòng)面模型Ω中的位置。無人機(jī)組合狀態(tài)下第k個(gè)無人機(jī)氣動(dòng)力FArk和氣動(dòng)力矩MArk對第j個(gè)無人機(jī)配平變量的各項(xiàng)導(dǎo)數(shù)如下：

式中：k，j=1，2，…，n。

2 飛行動(dòng)力學(xué)建模

多體組合式無人機(jī)飛行力學(xué)特性與傳統(tǒng)無人機(jī)飛行器存在顯著差異，其氣動(dòng)布局設(shè)計(jì)應(yīng)考慮飛行力學(xué)特性影響。本文基于Newton-Euler方程建立動(dòng)力學(xué)模型［26］，組合式無人機(jī)中第k個(gè)無人機(jī)相對慣性坐標(biāo)系的位置由坐標(biāo)rk=(xk，yk，zk)及歐拉角φk=(φk，θk，ψk)表示，其剛體運(yùn)動(dòng)學(xué)方程為

式中：上標(biāo)(k)表示機(jī)體坐標(biāo)系下變量；mk及為單體無人機(jī)質(zhì)量及慣量矩陣；為單體無人機(jī)所受外力及外力矩；=(pk，qk，rk)為無人機(jī)角速度向量；LEB為體軸系至慣性系轉(zhuǎn)換矩陣；S(·)為外積算子。結(jié)合飛行力學(xué)理論，單體無人機(jī)動(dòng)力學(xué)模型為

考慮n個(gè)單體無人機(jī)進(jìn)行翼尖鉸接組合，僅允許相鄰單體無人機(jī)沿機(jī)翼翼尖弦向方向自由轉(zhuǎn)動(dòng)，即允許單體無人機(jī)間的相對滾轉(zhuǎn)，約束關(guān)系不考慮摩擦、阻尼及剛度。約束關(guān)系可表達(dá)為

式中：u=。引入Lagrange 乘子λ表達(dá)約束關(guān)系，多體組合式無人機(jī)飛行動(dòng)力學(xué)方程為

式中：含下標(biāo)矩陣表示對應(yīng)的Jacobi 矩陣；t為時(shí)間；Ξ=diag(Ξ1，Ξ2，…，Ξn)；β=diag(β1，β2，…，βn)。氣動(dòng)力模型選擇氣動(dòng)導(dǎo)數(shù)模型，在規(guī)定好未知配平變量及配平自由度后，直接求解該非線性方程組即可給出配平狀態(tài)解。

將無人機(jī)組合系統(tǒng)動(dòng)力學(xué)方程（37）記做：

引入變量z=[uT，σT]T，=λ，對方程（39）進(jìn)行小擾動(dòng)線化處理，線化小擾動(dòng)方程可表示為

式中：、、分別為方 程ξ關(guān)于引入變量z的Jacobi 矩陣。小擾動(dòng)方程（40）的狀態(tài)空間形式可記做：

系統(tǒng)矩陣體現(xiàn)了多體組合式無人機(jī)系統(tǒng)的穩(wěn)定性特征，系統(tǒng)矩陣特征值在復(fù)平面左側(cè)則意味著飛行動(dòng)力學(xué)系統(tǒng)保持穩(wěn)定，反之則不能保持穩(wěn)定，通過求解系統(tǒng)矩陣特征值可完成對無人機(jī)系統(tǒng)的飛行力學(xué)穩(wěn)定性分析。

3 多體組合式無人機(jī)分析

3.1 單體無人機(jī)模型

多體組合式無人機(jī)算例由多個(gè)相同的單體常規(guī)布局固定翼無人機(jī)鉸接組成，單體無人機(jī)模型設(shè)計(jì)圖如圖3 所示，設(shè)計(jì)參數(shù)如表1 所示。其機(jī)身設(shè)計(jì)參考美國密歇根大學(xué)的X-HALE 無人機(jī)［27］。為了體現(xiàn)大展弦比機(jī)翼設(shè)計(jì)的同時(shí)盡量減小翼尖變形效應(yīng)的影響，該無人機(jī)模型采用雙機(jī)身結(jié)構(gòu)，全機(jī)設(shè)計(jì)質(zhì)量8 kg，展弦比約為11.1，共分為3 段，其中左機(jī)翼、右機(jī)翼及中部機(jī)翼展長相等均為1.0 m，機(jī)翼至尾翼距離（前緣）為1.02 m。由于設(shè)計(jì)飛行速度較低，無人機(jī)采用矩形平直翼構(gòu)型，其中左、右兩側(cè)機(jī)翼后緣布置副翼，左副翼及右副翼分別貫穿左機(jī)翼及右機(jī)翼，控制面弦長占機(jī)翼弦長的1/3。機(jī)身與尾翼之間通過由機(jī)身延伸的尾撐桿連接。

表1 單體無人機(jī)設(shè)計(jì)參數(shù)Table 1 Design parameters of single UAV

圖3 單體無人機(jī)模型Fig.3 Model of single UAV

3.2 氣動(dòng)導(dǎo)數(shù)計(jì)算及驗(yàn)證

通過對單體無人機(jī)氣動(dòng)導(dǎo)數(shù)的計(jì)算驗(yàn)證狀態(tài)空間渦格法氣動(dòng)力模型的正確性。氣動(dòng)網(wǎng)格沿機(jī)翼氣動(dòng)面展向劃分90 個(gè)網(wǎng)格，弦向劃分8 個(gè)網(wǎng)格，共720 個(gè)氣動(dòng)網(wǎng)格單元；沿尾翼氣動(dòng)面展向劃分46 個(gè)網(wǎng)格，弦向劃分8 個(gè)網(wǎng)格，共368 個(gè)氣動(dòng)網(wǎng)格單元，全機(jī)共計(jì)1 088 個(gè)物面渦氣動(dòng)單元，劃分如圖4 所示。本文在計(jì)算氣動(dòng)導(dǎo)數(shù)時(shí)不考慮機(jī)身影響。尾渦弦向劃分網(wǎng)格數(shù)量給定為20 倍物面渦弦向網(wǎng)格數(shù)量，即沿來流方向劃分160 個(gè)網(wǎng)格，共21 760 個(gè)氣動(dòng)網(wǎng)格單元。計(jì)算工況設(shè)定飛行速度20 m/s，迎角1°，大氣密度1.225 kg/m3。取矩中心指定為距離機(jī)翼前緣0.067 5 m（機(jī)翼弦長的1/4 處）?；跔顟B(tài)空間渦格法計(jì)算得到的氣動(dòng)壓強(qiáng)分布結(jié)果如圖5 所示，圖中ΔCp表示氣動(dòng)壓強(qiáng)數(shù)值。圖6 給出了機(jī)翼展向升力分布情況，圖7 給出了尾翼展向升力分布情況?？梢姍C(jī)翼及尾翼升力分布大致滿足橢圓分布，符合大展弦比機(jī)翼的升力分布規(guī)律。由于存在尾翼物面渦及尾渦的誘導(dǎo)，機(jī)翼中對應(yīng)尾翼位置的升力存在一定程度的下降。

圖4 單體無人機(jī)氣動(dòng)網(wǎng)格劃分Fig.4 Aerodynamics grid of single UAV

圖5 單體無人機(jī)氣動(dòng)壓強(qiáng)分布Fig.5 Aerodynamic pressure distribution of single UAV

圖6 機(jī)翼展向升力分布Fig.6 Lift distribution of wing in spanwise

圖7 尾翼展向升力分布Fig.7 Lift distribution of tail in spanwise

采用本文方法計(jì)算單體無人機(jī)模型氣動(dòng)導(dǎo)數(shù)，在氣動(dòng)力分析軟件ZONAIR［28］中建立相同的氣動(dòng)模型，氣動(dòng)網(wǎng)格劃分與渦格法狀態(tài)空間形式保持一致，給定相同的計(jì)算工況條件，2 種方法的計(jì)算結(jié)果如表2 所示。可以看出本文建立的氣動(dòng)導(dǎo)數(shù)計(jì)算方法與ZONAIR 軟件的氣動(dòng)導(dǎo)數(shù)計(jì)算結(jié)果一致性良好，重要?dú)鈩?dòng)導(dǎo)數(shù)如CLα、CMyα及CLq等最大偏差在3%左右，與升力相關(guān)的氣動(dòng)導(dǎo)數(shù)偏差不超過1%。另外，由于算例模型中構(gòu)型對稱，橫縱向氣動(dòng)導(dǎo)數(shù)解耦，諸如CMxα、CLβ及CLp等橫縱向耦合氣動(dòng)導(dǎo)數(shù)在2 種方法計(jì)算中結(jié)果均為0。計(jì)算結(jié)果驗(yàn)證了本文氣動(dòng)力建模方法及氣動(dòng)導(dǎo)數(shù)計(jì)算方法的有效性和正確性。

表2 氣動(dòng)導(dǎo)數(shù)計(jì)算結(jié)果Table 2 Results of aerodynamic derivatives

3.3 組合式無人機(jī)飛行力學(xué)分析

以2 架單體無人機(jī)構(gòu)建組合式無人機(jī)為例，進(jìn)行氣動(dòng)特性分析與氣動(dòng)布局參數(shù)影響分析。組合式無人機(jī)中將2 個(gè)單體固定翼無人機(jī)翼尖沿展向鉸鏈連接，只允許各無人機(jī)間出現(xiàn)沿機(jī)身x軸線的相對滾轉(zhuǎn)運(yùn)動(dòng)，其他相對自由度被鉸鏈約束。雙機(jī)組合示意圖如圖8 所示。

圖8 雙機(jī)組合示意圖Fig.8 Illustration of combination of two UAVs

基于前述飛行力學(xué)方程可求解按照指定配平狀態(tài)求解配平變量［26］。一般而言，飛行器實(shí)際飛行滿足對稱構(gòu)型飛行姿態(tài)。圖9 及圖10 分別給出了1 號機(jī)具有10°滾轉(zhuǎn)角及2 號機(jī)具有-10°滾轉(zhuǎn)角情況下的氣動(dòng)網(wǎng)格劃分情況及壓強(qiáng)分布情況。計(jì)算工況仍取飛行速度20 m/s，迎角1°，大氣密度1.225 kg/m3。

圖9 多體組合式無人機(jī)氣動(dòng)網(wǎng)格劃分Fig.9 Aerodynamics grid of MCA

圖10 多體組合式無人機(jī)氣動(dòng)壓強(qiáng)分布Fig.10 Aerodynamic pressure distribution of MCA

圖11 及圖12 分別給出了多體組合式無人機(jī)機(jī)翼及尾翼展向升力分布計(jì)算結(jié)果，可以看出機(jī)翼及尾翼升力分布同樣大致滿足橢圓分布，符合大展弦比機(jī)翼的氣動(dòng)特性規(guī)律。由于存在雙尾翼物面渦及尾渦的誘導(dǎo)，機(jī)翼中對應(yīng)尾翼位置的升力存在一定程度的下降。

圖11 多體組合式無人機(jī)機(jī)翼展向升力分布Fig.11 Lift distribution of wing of MCA in spanwise

圖12 多體組合式無人機(jī)尾翼展向升力分布Fig.12 Lift distribution of tail of MCA in spanwise

基于該狀態(tài)進(jìn)行飛行力學(xué)分析，求解配平變量，計(jì)算飛行力學(xué)模態(tài)。模態(tài)求解結(jié)果如表3 所示。整機(jī)飛行力學(xué)模態(tài)中除長周期模態(tài)外，縱向穩(wěn)定性模態(tài)均收斂，橫向穩(wěn)定性模態(tài)中的滾轉(zhuǎn)模態(tài)及荷蘭滾模態(tài)亦收斂。除了傳統(tǒng)飛行力學(xué)模態(tài)外，組合構(gòu)型因自由度增加在特征值分析中多出現(xiàn)一個(gè)發(fā)散的飛行力學(xué)復(fù)合運(yùn)動(dòng)模態(tài)1 及一個(gè)收斂的飛行力學(xué)復(fù)合運(yùn)動(dòng)模態(tài)2。基于復(fù)合模態(tài)1 中的1 號機(jī)滾轉(zhuǎn)角對動(dòng)力學(xué)方程特征模態(tài)做歸一化，即令Δφ1=1，兩階復(fù)合運(yùn)動(dòng)模態(tài)特征向量中角度量分析結(jié)果如表4 所示，角速度量分析結(jié)果如表5 所示。發(fā)散的復(fù)合運(yùn)動(dòng)模態(tài)1 呈現(xiàn)相對向上的“V 字”形態(tài)，其特征時(shí)間為0.42 s，對飛行力學(xué)穩(wěn)定性有顯著影響，飛行器先天不穩(wěn)定，這是多體組合式無人機(jī)在飛行力學(xué)特性方面與常規(guī)無人機(jī)飛行器的顯著差異。復(fù)合運(yùn)動(dòng)模態(tài)1 示意圖如圖13 所示。

表3 飛行力學(xué)模態(tài)特征根求解結(jié)果Table 3 Results of flight dynamics modes

表4 飛行力學(xué)模態(tài)特征向量角度量Table 4 Angular components in flight dynamics modes

表5 飛行力學(xué)模態(tài)特征向量角速度量Table 5 Angular velocity components in flight dynamics modes

圖13 發(fā)散的飛行力學(xué)模態(tài)示意圖Fig.13 Illustration of divergent flight dynamics mode

4 氣動(dòng)布局參數(shù)影響

與傳統(tǒng)飛行器不同，組合式無人機(jī)中單體無人機(jī)滾轉(zhuǎn)角變化會(huì)帶來整機(jī)氣動(dòng)布局的差異，需要針對單體無人機(jī)滾轉(zhuǎn)角變化開展參數(shù)影響研究。同時(shí)，后掠角、機(jī)翼尾翼距離及上反角是典型的飛行器氣動(dòng)布局設(shè)計(jì)參數(shù)，單體無人機(jī)改變上述氣動(dòng)布局參數(shù)同樣會(huì)帶來組合式無人機(jī)整機(jī)的氣動(dòng)布局變化。本文針對相關(guān)參數(shù)變化對飛行力學(xué)特性的影響，特別是穩(wěn)定性的影響作簡要討論。

4.1 單機(jī)滾轉(zhuǎn)角

本文研究中假設(shè)雙機(jī)飛行中保持“上折”姿態(tài)，即1 號機(jī)滾轉(zhuǎn)角為正，2 號機(jī)滾轉(zhuǎn)角為負(fù)。為方便討論，后文中2 號機(jī)滾轉(zhuǎn)角值取絕對值。改變單體無人機(jī)滾轉(zhuǎn)角，逐步從1°增加至15°，圖14～圖18給出了與升力及俯仰力矩相關(guān)的重要?dú)鈩?dòng)導(dǎo)數(shù)計(jì)算結(jié)果。從左至右分別為組合式無人機(jī)整機(jī)氣動(dòng)導(dǎo)數(shù)計(jì)算結(jié)果，組合式無人機(jī)1 號機(jī)氣動(dòng)導(dǎo)數(shù)計(jì)算結(jié)果及組合式無人機(jī)2 號機(jī)氣動(dòng)導(dǎo)數(shù)計(jì)算結(jié)果。橫坐標(biāo)為1 號機(jī)滾轉(zhuǎn)角，縱坐標(biāo)為2 號機(jī)滾轉(zhuǎn)角。定義氣動(dòng)導(dǎo)數(shù)展示結(jié)果中的45°對角線為兩機(jī)滾轉(zhuǎn)角一致線，在此線上，1 號機(jī)與2 號機(jī)滾轉(zhuǎn)角絕對值相同。整機(jī)氣動(dòng)導(dǎo)數(shù)CLα隨兩機(jī)相對滾轉(zhuǎn)角增大而減小且沿兩機(jī)滾轉(zhuǎn)角一致線對稱分布；1 號機(jī)與2 號機(jī)的CLα沿一致線非對稱分布且無人機(jī)自身滾轉(zhuǎn)角對自身氣動(dòng)導(dǎo)數(shù)影響最大，這與展弦比逐漸減小時(shí)的氣動(dòng)力特性規(guī)律相似。類似規(guī)律也出現(xiàn)在CLp、CLq、CMyα、CMyq計(jì)算中。

圖14 氣動(dòng)導(dǎo)數(shù)CLα 變化Fig.14 Variation of aerodynamics derivatives CLα

圖15 氣動(dòng)導(dǎo)數(shù)CLp 變化Fig.15 Variation of aerodynamics derivatives CLp

圖16 氣動(dòng)導(dǎo)數(shù)CLq 變化Fig.16 Variation of aerodynamics derivatives CLq

圖17 氣動(dòng)導(dǎo)數(shù)CMyα 變化Fig.17 Variation of aerodynamics derivatives CMyα

圖18 氣動(dòng)導(dǎo)數(shù)CMyq 變化Fig.18 Variation of aerodynamics derivatives CMyq

圖19 和圖20 給出了與滾轉(zhuǎn)力矩相關(guān)的重要?dú)鈩?dòng)導(dǎo)數(shù)CMxα及CMxq計(jì)算結(jié)果。CMxα屬于典型的橫縱向耦合氣動(dòng)導(dǎo)數(shù)，對于常規(guī)布局固定翼無人機(jī)而言，CMxα=0。對于具有相對滾轉(zhuǎn)自由度的無人機(jī)組合狀態(tài)而言，當(dāng)兩機(jī)滾轉(zhuǎn)自由度不同時(shí)，整機(jī)構(gòu)型不對稱出現(xiàn)非零CMxα。更重要的，在分析1 號機(jī)及2 號機(jī)各自的氣動(dòng)導(dǎo)數(shù)時(shí)，即使兩機(jī)滾轉(zhuǎn)角相同，由于組合整機(jī)滿足氣動(dòng)力對稱近似的橢圓分布，1 號機(jī)和2 號機(jī)上也會(huì)因?yàn)樽訜o人機(jī)上氣動(dòng)力分布不對稱而具有非零的CMxα，這與傳統(tǒng)構(gòu)型無人機(jī)具有顯著差別，在分析組合飛行狀態(tài)下的飛行動(dòng)力學(xué)時(shí)，必須考慮該類氣動(dòng)耦合效應(yīng)，CMxq計(jì)算展現(xiàn)了相似的規(guī)律。事實(shí)上，從氣動(dòng)導(dǎo)數(shù)計(jì)算結(jié)果看出，幾乎全部氣動(dòng)導(dǎo)數(shù)都存在耦合效應(yīng)，無人機(jī)組合飛行下橫縱向氣動(dòng)力無法解耦。

圖19 氣動(dòng)導(dǎo)數(shù)CMxα 變 化Fig.19 Variation of aerodynamics derivativesCMxα

圖20 氣動(dòng)導(dǎo)數(shù)CMxq 變化Fig.20 Variation of aerodynamics derivativesCMxq

對稱改變各子無人機(jī)滾轉(zhuǎn)角，保證兩單體無人機(jī)滾轉(zhuǎn)角絕對值一致即以對稱“上折”姿態(tài)飛行。圖21 給出了發(fā)散復(fù)合運(yùn)動(dòng)模態(tài)特征值隨滾轉(zhuǎn)角變化的計(jì)算結(jié)果?？梢园l(fā)現(xiàn)，系統(tǒng)穩(wěn)定性隨配平滾轉(zhuǎn)角的增大而逐漸減低，具體表現(xiàn)為特征值的增大及發(fā)散特征時(shí)間的縮短。系統(tǒng)穩(wěn)定性隨配平滾轉(zhuǎn)角的增大而逐漸減低。圖22 給出了俯仰角/偏航角特征模態(tài)分量隨滾轉(zhuǎn)角變化情況，基于復(fù)合模態(tài)1 中的1 號機(jī)滾轉(zhuǎn)角對特征模態(tài)分量進(jìn)行了歸一化處理，1 號機(jī)與2 號機(jī)相同運(yùn)動(dòng)分量數(shù)值幾乎重合?？梢园l(fā)現(xiàn)俯仰角、偏航角分量相對于滾轉(zhuǎn)角而言較小，變化趨勢平緩，發(fā)散的飛行力學(xué)模態(tài)并未由以相對滾轉(zhuǎn)為主導(dǎo)改變至其他模態(tài)運(yùn)動(dòng)形式。

圖21 發(fā)散模態(tài)特征值隨滾轉(zhuǎn)角變化Fig.21 Variation of divergence modal eigenvalue with roll angle

圖22 發(fā)散模態(tài)分量變化Fig.22 Variation of divergence modal eigenvector

4.2 單機(jī)后掠角

改變單體無人機(jī)機(jī)翼后掠角，逐步從1°增加至15°，多體組合式無人機(jī)飛行動(dòng)力學(xué)特性仍以出現(xiàn)發(fā)散的復(fù)合運(yùn)動(dòng)模態(tài)為主要特征。圖23 給出了單體無人機(jī)后掠布局形式的氣動(dòng)網(wǎng)格劃分情況。圖24 給出了發(fā)散復(fù)合運(yùn)動(dòng)模態(tài)特征值隨單機(jī)后掠角變化的計(jì)算結(jié)果。隨后掠角增加，特征根實(shí)部在保持為正的情況下逐漸減小，特征時(shí)間隨后掠角增大而增大，雙機(jī)對稱滾轉(zhuǎn)的不穩(wěn)定性有所改善。

圖23 單體無人機(jī)后掠布局Fig.23 UAV layout of sweep angle

圖24 發(fā)散模態(tài)特征值隨后掠角變化Fig.24 Variation of divergence modal eigenvalue with sweep angle

4.3 單機(jī)機(jī)翼尾翼間距

改變單體無人機(jī)機(jī)翼至尾翼距離，逐步從0.7 m 增加至1.3 m，圖25 給出了發(fā)散復(fù)合運(yùn)動(dòng)模態(tài)特征值隨單機(jī)機(jī)翼尾翼間距變化的計(jì)算結(jié)果。隨著距離的逐步增大，對應(yīng)特征根保持為正的情況下先減小后增大，特征時(shí)間隨距離增大出現(xiàn)先增大后減小的趨勢，其中距離1.07 m 為極值點(diǎn)。機(jī)翼與尾翼距離可能存在改善飛行穩(wěn)定性的最優(yōu)值。

圖25 發(fā)散模態(tài)特征值隨機(jī)翼尾翼距離變化Fig.25 Variation of divergence modal eigenvalue with distance of wing and tail

4.4 單機(jī)上反角

改變單體無人機(jī)機(jī)翼上反角，逐步從0°增加至15°。圖26 給出了單體無人機(jī)機(jī)翼上反布局形式的氣動(dòng)網(wǎng)格劃分情況，圖27 給出了發(fā)散復(fù)合運(yùn)動(dòng)模態(tài)特征值隨單機(jī)上反角變化的計(jì)算結(jié)果。多體組合式無人機(jī)飛行動(dòng)力學(xué)特性仍以出現(xiàn)發(fā)散的復(fù)合運(yùn)動(dòng)模態(tài)為主要特征。但隨上反角增加，特征實(shí)部在保持為正的情況下變化很小，單體無人機(jī)上反角對于整機(jī)發(fā)散復(fù)合運(yùn)動(dòng)模態(tài)影響很小。

圖26 單體無人機(jī)上反布局Fig.26 UAV layout of dihedral angle

圖27 發(fā)散模態(tài)特征值隨上反角變化Fig.27 Variation of divergence modal eigenvalue with dihedral angle

5 結(jié)論

基于渦格法狀態(tài)空間形式推導(dǎo)了適合于考慮氣動(dòng)耦合效應(yīng)的翼尖鉸接組合式無人機(jī)氣動(dòng)導(dǎo)數(shù)計(jì)算方法，基于氣動(dòng)導(dǎo)數(shù)模型及Newton-Euler 方程建立多體組合式無人機(jī)飛行動(dòng)力學(xué)模型，計(jì)算配平狀態(tài)及飛行力學(xué)穩(wěn)定性。針對翼尖鉸接組合式無人機(jī)特點(diǎn)，進(jìn)行了布局參數(shù)對飛行力學(xué)穩(wěn)定性影響研究，得到以下結(jié)論：

1）狀態(tài)空間渦格法可以給出解析形式的氣動(dòng)導(dǎo)數(shù)計(jì)算方法，同時(shí)針對翼尖鉸接組合式無人機(jī)可以得到考慮多機(jī)間氣動(dòng)耦合效應(yīng)的氣動(dòng)導(dǎo)數(shù)，經(jīng)驗(yàn)證計(jì)算結(jié)果精度好。需要說明的是，機(jī)身對于飛行器橫側(cè)向穩(wěn)定性及氣動(dòng)導(dǎo)數(shù)具有一定影響，本文方法未考慮機(jī)身影響，如何進(jìn)行準(zhǔn)確的機(jī)身影響修正將會(huì)是后續(xù)研究工作的重點(diǎn)。

2）多體組合式無人機(jī)飛行力學(xué)特性與傳統(tǒng)飛行器有較大不同，本文基于Newton-Euler 方程建立多體組合式無人機(jī)飛行動(dòng)力學(xué)模型，考慮鉸接約束對飛行動(dòng)力學(xué)特性的影響。對于多體組合式無人機(jī)，氣動(dòng)導(dǎo)數(shù)及飛行力學(xué)方程橫縱向無法解耦，在配平狀態(tài)下進(jìn)行組合下的動(dòng)力學(xué)系統(tǒng)穩(wěn)定性分析會(huì)出現(xiàn)先天不穩(wěn)定的以相對滾轉(zhuǎn)為主導(dǎo)的復(fù)合運(yùn)動(dòng)模態(tài)，在氣動(dòng)布局設(shè)計(jì)中需要考慮參數(shù)對飛行力學(xué)特性的影響。

3）改變單體無人機(jī)配平狀態(tài)滾轉(zhuǎn)角會(huì)改變整機(jī)氣動(dòng)布局，減小單體無人機(jī)配平滾轉(zhuǎn)角可以改善飛行力學(xué)穩(wěn)定性。對于傳統(tǒng)布局參數(shù)，增加單體無人機(jī)后掠角亦可以改善飛行力學(xué)穩(wěn)定性；上反角影響很??；機(jī)翼與尾翼距離可能存在改善飛行力學(xué)穩(wěn)定的最優(yōu)值。