劉東旭，謝長(zhǎng)川，洪冠新

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

作為一種特定的新概念飛行器系統(tǒng)，翼尖鉸接復(fù)合飛行技術(shù)具體是以單架次分別進(jìn)行起飛爬升和下降著陸，在巡航過(guò)程中通過(guò)特定設(shè)計(jì)的翼尖鉸柱接合機(jī)構(gòu)實(shí)現(xiàn)飛機(jī)復(fù)合飛行，這種技術(shù)能夠有效降低飛機(jī)的油耗率，提升飛機(jī)的航程和航時(shí)等續(xù)航性能指標(biāo)。對(duì)于遠(yuǎn)距離長(zhǎng)航時(shí)飛機(jī)，這種技術(shù)的優(yōu)點(diǎn)具體體現(xiàn)在：可以增加展弦比、翼型彎度，同時(shí)飛機(jī)還會(huì)受到彼此翼尖渦流的上洗抬升作用以增加低速巡航性能；多架飛行器可以在苛刻條件下以單一架次起降，在空中對(duì)接后再執(zhí)行飛行任務(wù)和作戰(zhàn)任務(wù)，這樣就降低了飛機(jī)的起降環(huán)境要求。從經(jīng)濟(jì)性方面考慮，氣動(dòng)性能的提升意味著復(fù)合飛行系統(tǒng)可以節(jié)省燃油的布置空間和攜帶質(zhì)量，在總體設(shè)計(jì)工作流程方面提高飛機(jī)的設(shè)計(jì)靈活性。作為一種提升續(xù)航性能的飛行技術(shù)，翼尖對(duì)接復(fù)合技術(shù)經(jīng)常被用來(lái)和空中加油技術(shù)進(jìn)行比較。翼尖對(duì)接復(fù)合技術(shù)在一定程度上使作戰(zhàn)飛機(jī)可以免于空中加油機(jī)的輔助燃油保障，從而簡(jiǎn)化作戰(zhàn)出動(dòng)飛機(jī)的機(jī)種，改善機(jī)隊(duì)的出動(dòng)任務(wù)規(guī)劃方法和起降保障作業(yè)的流程工作。

翼尖鉸接復(fù)合飛行［1］的概念最初由德國(guó)科學(xué)家在二戰(zhàn)期間提出，此后美國(guó)空軍和業(yè)內(nèi)各研究機(jī)構(gòu)制定了Tom-Tom項(xiàng)目、Ficon項(xiàng)目為主的軍方試驗(yàn)計(jì)劃［2］驗(yàn)證這個(gè)技術(shù)理論，同時(shí)設(shè)計(jì)研發(fā)了多個(gè)翼尖鉸接復(fù)合飛行概念驗(yàn)證機(jī)型號(hào)。曾在美國(guó)康奈爾大學(xué)和德雷伯實(shí)驗(yàn)室工作的Edgar教授通過(guò)使用渦格法、CFD法等數(shù)值仿真方法［3］和風(fēng)洞測(cè)試模擬仿真方法［4］開(kāi)展復(fù)合機(jī)翼的氣動(dòng)特性研究，驗(yàn)證了這類復(fù)合飛行技術(shù)的氣動(dòng)性能改善效果。美國(guó)空軍從事的“Ficon”計(jì)劃設(shè)計(jì)了一個(gè)由戰(zhàn)略轟炸機(jī)組成的進(jìn)攻單元，由1架B-29轟炸機(jī)和2架F-84 D／E，戰(zhàn)斗機(jī)以翼尖鉸接形式復(fù)合飛行；而在風(fēng)洞測(cè)試、地效實(shí)驗(yàn)以及幾次成功的對(duì)接實(shí)驗(yàn)等驗(yàn)證工作之后，“Ficon”計(jì)劃的飛行試驗(yàn)遭遇了嚴(yán)重失敗。由于缺少協(xié)作飛行策略，其中一架飛機(jī)從對(duì)接的飛機(jī)脫離，該項(xiàng)目組在對(duì)接后僅幾秒內(nèi)就墜海失事［2］。此后翼尖鉸接復(fù)合飛行技術(shù)的動(dòng)力學(xué)研究集中在帶轉(zhuǎn)動(dòng)鉸約束條件的多剛體系統(tǒng)動(dòng)力學(xué)問(wèn)題。多剛體系統(tǒng)指的是多個(gè)剛體按照確定的方式相互聯(lián)系所組成的機(jī)械系統(tǒng)。多剛體系統(tǒng)動(dòng)力學(xué)是在經(jīng)典力學(xué)基礎(chǔ)上產(chǎn)生的新學(xué)科分支。處理由多個(gè)剛體組成的系統(tǒng)，原則上可以利用傳統(tǒng)的經(jīng)典力學(xué)方法，即以Newton-Euler方程為代表的矢量力學(xué)方法或以拉格朗日方法為代表的分析力學(xué)方法。同時(shí)以此為基礎(chǔ)，Roberson-Wittenburg方法主要是利用圖論的概念及數(shù)學(xué)工具描述多剛體系統(tǒng)的結(jié)構(gòu)，以鄰接剛體之間的相對(duì)位移作為廣義坐標(biāo)，導(dǎo)出適合于任意多剛體系統(tǒng)的普遍動(dòng)力學(xué)方程，利用費(fèi)舍爾的增廣體概念對(duì)方程的系數(shù)矩陣A作出物理解釋［5］。在此基礎(chǔ)上，佛羅里達(dá)大學(xué)的Chakravarthy針對(duì)變后掠翼飛行器［6］，上海交通大學(xué)的安繼光［7］、北京航空航天大學(xué)的樂(lè)挺和王立新教授［8］等針對(duì)Z形折疊翼面飛行器［7-8］、非對(duì)稱斜翼飛行器［9-10］和變后掠翼飛行器［11-12］開(kāi)展了動(dòng)力學(xué)建模和仿真工作，Thomas［13-14］、Li［15］與Yan［16］等針對(duì)類似形式的多種多剛體系統(tǒng)飛行器也開(kāi)展了建模仿真研究，為多剛體系統(tǒng)飛行器動(dòng)力學(xué)建模和仿真方法提供了參考。對(duì)于這種特定的翼尖鉸接復(fù)合飛行器，Brandon和Montalvo［17］給出了按單體飛機(jī)添加約束的傳統(tǒng)牛頓歐拉動(dòng)力學(xué)建模方法，同時(shí)使用增益調(diào)度方法完成了控制特性研究，又對(duì)對(duì)接飛行軌跡展開(kāi)了探討，但對(duì)于多個(gè)飛行器組成的復(fù)合系統(tǒng)基于單個(gè)飛機(jī)的氣動(dòng)數(shù)據(jù)和配平方案，沒(méi)有建立多機(jī)復(fù)合的方法給出多機(jī)的配平方案。

為研究翼尖鉸接復(fù)合飛行器的動(dòng)力學(xué)特性，本文以2架中空長(zhǎng)航時(shí)無(wú)人機(jī)通過(guò)翼尖鉸接機(jī)構(gòu)組成的復(fù)合飛行器系統(tǒng)作為研究對(duì)象，考慮復(fù)合飛行過(guò)程中系統(tǒng)整體和內(nèi)部的運(yùn)動(dòng)特點(diǎn)，利用 Newton-Euler方法、Lagrange方法和Roberson-Wittenburg方法開(kāi)展系統(tǒng)動(dòng)力學(xué)建模工作，使用CFD方法建立面向復(fù)合運(yùn)動(dòng)特征參數(shù)和復(fù)合飛行舵面部件的雙機(jī)復(fù)合飛行器系統(tǒng)的氣動(dòng)力數(shù)據(jù)庫(kù)。在此基礎(chǔ)上研究復(fù)合飛行器的本征動(dòng)力學(xué)特性和動(dòng)力學(xué)響應(yīng)，研究結(jié)果可以為復(fù)合飛行器的動(dòng)力學(xué)特性和配平策略的設(shè)計(jì)提供參考依據(jù)。

1 復(fù)合動(dòng)力學(xué)建模

1.1 復(fù)合飛行器系統(tǒng)

算例采用的復(fù)合飛行器系統(tǒng)是由2架常規(guī)尾翼布局中空長(zhǎng)航時(shí)無(wú)人機(jī)通過(guò)翼尖鉸接機(jī)構(gòu)實(shí)現(xiàn)復(fù)合飛行，三維實(shí)體設(shè)計(jì)如圖1所示。

翼尖鉸接機(jī)構(gòu)主要由鉸柱約束機(jī)構(gòu)組成，如圖2所示。翼尖鉸接約束機(jī)構(gòu)可以限制2架飛機(jī)相互之間的3個(gè)相對(duì)平動(dòng)自由度和2個(gè)相對(duì)轉(zhuǎn)動(dòng)自由度，具體是前后、上下和左右方向上的相對(duì)平動(dòng)運(yùn)動(dòng)以及俯仰相對(duì)轉(zhuǎn)動(dòng)、偏航相對(duì)轉(zhuǎn)動(dòng)，同時(shí)允許2架飛機(jī)之間存在相對(duì)滾轉(zhuǎn)運(yùn)動(dòng)。

圖2 翼尖鉸接機(jī)構(gòu)Fig.2 Wingtip joint mechanism

1.2 復(fù)合系統(tǒng)動(dòng)力學(xué)和運(yùn)動(dòng)學(xué)方程

翼尖鉸接復(fù)合飛行器的研究屬于帶轉(zhuǎn)動(dòng)鉸約束條件的無(wú)根樹系統(tǒng)動(dòng)力學(xué)問(wèn)題，可以根據(jù)自由度和約束形式確定列寫方程的數(shù)量。

在歐美制單剛體機(jī)體坐標(biāo)軸系下，分析飛機(jī)的運(yùn)動(dòng)學(xué)與動(dòng)力學(xué)問(wèn)題。2架飛機(jī)分別各自擁有6個(gè)自由度，具體為x、y、z軸的3個(gè)平動(dòng)自由度和3個(gè)轉(zhuǎn)動(dòng)自由度。如圖3所示，在翼尖一側(cè)添加鉸鏈約束后，約束了2架飛機(jī)的x、y、z軸3個(gè)軸向自由度和偏航、俯仰2個(gè)轉(zhuǎn)動(dòng)自由度。因此，2架飛機(jī)共有6＋6－3－2＝7種自由度，具體為：整體x、y、z軸方向的3個(gè)軸向自由度、1個(gè)俯仰轉(zhuǎn)動(dòng)自由度、1個(gè)偏航轉(zhuǎn)動(dòng)自由度和單獨(dú)的2個(gè)滾轉(zhuǎn)自由度。根據(jù)多剛體系統(tǒng)一般動(dòng)力學(xué)原理，建立7自由度力學(xué)模型，解決7自由度的動(dòng)力學(xué)問(wèn)題，需要按照2架飛機(jī)的7個(gè)自由度建立7個(gè)動(dòng)力學(xué)方程才可求解，而每個(gè)方向的動(dòng)力學(xué)方程都會(huì)需要具體分析才能完整建立。

圖3 復(fù)合飛行器系統(tǒng)運(yùn)動(dòng)分析Fig.3 Composite aircraft system’s motion analysis

推導(dǎo)過(guò)程分為2部分：①采用整體法推導(dǎo)復(fù)合飛行器系統(tǒng)整體的六自由度動(dòng)力學(xué)方程和運(yùn)動(dòng)學(xué)方程；②采用分離法推導(dǎo)復(fù)合飛行器系統(tǒng)內(nèi)部的相對(duì)轉(zhuǎn)動(dòng)動(dòng)力學(xué)方程。

1.2.1 系統(tǒng)整體動(dòng)力學(xué)方程和運(yùn)動(dòng)學(xué)方程

根據(jù)復(fù)合飛行的運(yùn)動(dòng)受力狀態(tài)分析，重新定義一種適用于復(fù)合飛行運(yùn)動(dòng)的機(jī)體坐標(biāo)系——復(fù)合機(jī)體坐標(biāo)系Oxfyfzf，如圖4所示。

圖4 翼尖鉸接復(fù)合飛行器機(jī)體坐標(biāo)系Fig.4 Airframe coordinate system of wingtip-jointed composite aircraft

原點(diǎn)位于兩復(fù)合飛機(jī)機(jī)體的鉸接點(diǎn)O。Oxf軸在復(fù)合飛機(jī)對(duì)稱平面內(nèi)，平行于各單體飛機(jī)機(jī)身軸線，指向前；Ozf軸亦在對(duì)稱平面內(nèi)，垂直于Oxf軸，指向下；Oyf軸垂直于對(duì)稱平面，指向右。這樣選取對(duì)應(yīng)的坐標(biāo)系可以直接表示復(fù)合飛機(jī)整體的運(yùn)動(dòng)和姿態(tài)，便于飛行器駕駛員或控制系統(tǒng)操縱監(jiān)視過(guò)程；單一架次飛機(jī)的氣動(dòng)力和氣動(dòng)力矩也可直接通過(guò)矩陣變換以整體形式表示，重力和發(fā)動(dòng)機(jī)推力的表示方法也不繁瑣。同時(shí)，復(fù)合飛機(jī)整體質(zhì)心C位于2架飛機(jī)質(zhì)心C1和C2連線的中點(diǎn)，OC可由2飛機(jī)鉸接角度表出。

參照以上復(fù)合機(jī)體坐標(biāo)系對(duì)應(yīng)定義機(jī)體復(fù)合機(jī)體氣流坐標(biāo)系Oxfayfazfa、復(fù)合機(jī)體航跡坐標(biāo)系Oxfkyfkzfk和常規(guī)大地坐標(biāo)系Oxeyeze等。

根據(jù)一般剛體動(dòng)力學(xué)方法，2架飛機(jī)各自在復(fù)合機(jī)體坐標(biāo)系Oxfyfzf下的速度可以表示為

vi＝vo＋˙rOCi＋ωi×rOCii＝1，2（1）

式中：vo為O點(diǎn)的速度；rOCi為第i個(gè)剛體質(zhì)心C點(diǎn)相對(duì)O點(diǎn)的矢經(jīng)。

式中：mi為各飛機(jī)的質(zhì)量；vi為各飛機(jī)的飛行速度；ω為復(fù)合機(jī)體坐標(biāo)系相對(duì)于慣性坐標(biāo)系的轉(zhuǎn)動(dòng)角速度；ICi為各飛機(jī)相對(duì)整體質(zhì)心的轉(zhuǎn)動(dòng)慣量張量；ωir為各飛機(jī)相對(duì)復(fù)合機(jī)體坐標(biāo)系的轉(zhuǎn)動(dòng)角速度；S為飛行器系統(tǒng)整體相對(duì)質(zhì)心的靜矩。

根據(jù)已有數(shù)學(xué)關(guān)系，得到復(fù)合飛行器多剛體系統(tǒng)整體動(dòng)力學(xué)方程的標(biāo)量形式

式中：Lef為復(fù)合機(jī)體坐標(biāo)系向大地坐標(biāo)系轉(zhuǎn)換的轉(zhuǎn)換矩陣；φ、θ、ψ為大地坐標(biāo)系和復(fù)合機(jī)體坐標(biāo)系之間的3個(gè)歐拉角；u、v、w為復(fù)合機(jī)體坐標(biāo)系上x、y、z三個(gè)坐標(biāo)軸方向上的速度；p、q、r在飛行動(dòng)力學(xué)學(xué)科中默認(rèn)為滾轉(zhuǎn)角速度、俯仰角速度、偏航角速度。

1.2.2 系統(tǒng)內(nèi)部動(dòng)力學(xué)方程和運(yùn)動(dòng)學(xué)方程

在飛機(jī)復(fù)合飛行的過(guò)程中，存在多剛體相對(duì)轉(zhuǎn)動(dòng)，如圖5所示，φ12為兩復(fù)合飛機(jī)的相對(duì)滾轉(zhuǎn)角（右側(cè)飛機(jī)滾轉(zhuǎn)角為φ1，左側(cè)飛機(jī)滾轉(zhuǎn)角為φ2）?？傻孟到y(tǒng)內(nèi)部的相對(duì)轉(zhuǎn)動(dòng)動(dòng)力學(xué)方程為

圖5 翼尖鉸接復(fù)合飛行器內(nèi)部運(yùn)動(dòng)分析Fig.5 Internal motion analysis of wingtip-jointed composite aircraft

式中：n為復(fù)合的飛機(jī)架次；g為重力加速度；mk、Vkx、φk、zk分別為第k架飛機(jī)的質(zhì)量、速度、滾轉(zhuǎn)角、z向坐標(biāo)。

將動(dòng)勢(shì)代入再考慮復(fù)合鉸接角度和飛行器姿態(tài)的變換關(guān)系φ1＝φ－φ12后，將式（11）展開(kāi)得

進(jìn)一步整理可得采用分離法得到的多剛體內(nèi)部相對(duì)轉(zhuǎn)動(dòng)參數(shù)和復(fù)合整體參數(shù)之間存在的運(yùn)動(dòng)關(guān)系表達(dá)式：

建立復(fù)合飛行多體動(dòng)力學(xué)模型得到最終表達(dá)形式，7個(gè)自由度坐標(biāo)方向下的零階量分別為（xf、yf、zf、φf(shuō)、θf(wàn)、ψf、φ12）。根據(jù)飛行器動(dòng)力學(xué)基本原理，在復(fù)合機(jī)體坐標(biāo)系Oxfyfzf下，得到復(fù)合系統(tǒng)整體和內(nèi)部動(dòng)力學(xué)方程（見(jiàn)式（14）），T為發(fā)動(dòng)機(jī)推力。

對(duì)應(yīng)有運(yùn)動(dòng)學(xué)方程和轉(zhuǎn)動(dòng)動(dòng)力學(xué)方程以及系統(tǒng)內(nèi)部運(yùn)動(dòng)學(xué)方程為

2 CFD方法建立氣動(dòng)力數(shù)據(jù)庫(kù)

由于翼尖鉸接復(fù)合飛行器系統(tǒng)飛行運(yùn)動(dòng)區(qū)別于單機(jī)飛行的特殊運(yùn)動(dòng)形式，相關(guān)的CFD數(shù)值仿真和風(fēng)洞模擬仿真試驗(yàn)［5］證明復(fù)合飛行器系統(tǒng)采用2架飛機(jī)以翼尖鉸接的形式實(shí)現(xiàn)復(fù)合飛行時(shí)，2架飛機(jī)的流場(chǎng)會(huì)產(chǎn)生氣動(dòng)耦合效應(yīng)，從而氣動(dòng)特性相比單機(jī)飛行時(shí)存在顯著變化，同時(shí)2架飛機(jī)之間的相對(duì)滾轉(zhuǎn)造成了飛行器系統(tǒng)的氣動(dòng)參數(shù)和氣動(dòng)舵面控制之間高度耦合無(wú)法解耦，因此需要使用CFD方法建立復(fù)合飛行氣動(dòng)力數(shù)據(jù)庫(kù)。

建立氣動(dòng)力數(shù)據(jù)庫(kù)時(shí)要考慮復(fù)合運(yùn)動(dòng)特征參數(shù)（主要是復(fù)合系統(tǒng)內(nèi)部的相對(duì)滾轉(zhuǎn)角φ12）和2架飛機(jī)各自的舵面部件偏轉(zhuǎn)角度（δa1，δa2）、（δe1，δe2）和（δr1，δr2）對(duì)氣動(dòng)特性的影響。考慮復(fù)合運(yùn)動(dòng)特征參數(shù)的原因是：①不同大小的相對(duì)滾轉(zhuǎn)角會(huì)對(duì)復(fù)合飛行器氣動(dòng)參數(shù)產(chǎn)生影響；②便于研究關(guān)聯(lián)飛行性能指標(biāo)的復(fù)合運(yùn)動(dòng)特征參數(shù)最優(yōu)值。考慮復(fù)合飛行時(shí)舵面部件偏轉(zhuǎn)角度的原因是：①相比單機(jī)常規(guī)飛行時(shí)，副翼、升降舵和方向舵舵面偏轉(zhuǎn)角度對(duì)復(fù)合飛行器的影響機(jī)理更加復(fù)雜；②復(fù)合機(jī)體配平偏轉(zhuǎn)影響。風(fēng)洞模擬試驗(yàn)［5］證明復(fù)合飛行時(shí)機(jī)翼存在的內(nèi)外側(cè)翼面壓差使內(nèi)外側(cè)機(jī)翼承受不相等的氣動(dòng)升力，可能要采用不同于單機(jī)飛行時(shí)的雙機(jī)舵面差動(dòng)配平方案，需要研究這種配平舵偏帶來(lái)的復(fù)合飛行氣動(dòng)特性變化；③便于精確研究復(fù)合運(yùn)動(dòng)特征參數(shù)的飛行性能指標(biāo)最優(yōu)值。

在建立面向復(fù)合運(yùn)動(dòng)特征參數(shù)和復(fù)合飛行舵面部件的氣動(dòng)數(shù)據(jù)庫(kù)之前，由多剛體飛行器系統(tǒng)的相關(guān)氣動(dòng)耦合效應(yīng)研究［10-11］仿真結(jié)果中可知，復(fù)合飛行過(guò)程里非定常效應(yīng)對(duì)復(fù)合飛行動(dòng)態(tài)特性的影響較小，同時(shí)在復(fù)合飛機(jī)相對(duì)滾轉(zhuǎn)角角速度較小時(shí)，能夠忽略非定常氣動(dòng)力效應(yīng)，即忽略所有氣動(dòng)力與相對(duì)滾轉(zhuǎn)角速度˙φ12的氣動(dòng)導(dǎo)數(shù)影響，因此可以利用以上氣動(dòng)力準(zhǔn)定常假設(shè)來(lái)簡(jiǎn)化計(jì)算機(jī)翼相對(duì)滾轉(zhuǎn)過(guò)程中復(fù)合飛行器系統(tǒng)的氣動(dòng)力。假設(shè)工況動(dòng)壓為k，翼面積為S，氣動(dòng)幾何平均弦長(zhǎng)為ˉc，復(fù)合機(jī)翼展長(zhǎng)為b，則復(fù)合飛行過(guò)程中氣動(dòng)力和氣動(dòng)力矩表達(dá)式為

式中：CD0為復(fù)合機(jī)體零升阻力系數(shù)；CDα和CDV分別為阻力系數(shù)對(duì)α和V的無(wú)因次導(dǎo)數(shù)；CL0為迎角為0時(shí)的升力系數(shù)；CLV、CLq、CLα分別為升力系數(shù)對(duì)V、q、α的無(wú)因次導(dǎo)數(shù)；Cm0為零升俯仰力矩系數(shù)；Cmα、CmV分別為俯仰力矩系數(shù)對(duì)α、V的無(wú)因次導(dǎo)數(shù)；Cmq和Cm˙α分別為俯仰阻尼導(dǎo)數(shù)和洗流時(shí)差導(dǎo)數(shù)；CDδa1δa2＝fDδa（δa1，δa2）、CDδe1δe2＝fDδe（δe1，δe2）、CDδr1δr2＝fDδr（δr1，δr2）等為復(fù)合機(jī)體舵面操縱導(dǎo)數(shù)，其中（δa1，δa2）、（δe1，δe2）和（δr1，δr2）為復(fù)合飛行器系統(tǒng)的舵面偏角二維數(shù)組，由雙機(jī)復(fù)合飛行時(shí)的舵面配平方案決定；y＝fDδa（x1，x2），y＝fDδe（x1，x2）和y＝fDδr（x1，x2）等函數(shù)關(guān)系由CFD仿真得到的離散數(shù)據(jù)點(diǎn)采用二元非線性回歸插值獲得。

具體地，建立氣動(dòng)數(shù)據(jù)庫(kù)的CFD仿真工作內(nèi)容、氣動(dòng)參數(shù)和建模方法具體如表1所示。

表1 氣動(dòng)力數(shù)據(jù)庫(kù)獲取方法Table 1 Aerodynamic database acquisition method

建立氣動(dòng)力數(shù)據(jù)庫(kù)時(shí)，復(fù)合飛行器系統(tǒng)的靜導(dǎo)數(shù)和操縱導(dǎo)數(shù)計(jì)算可以認(rèn)為是定常流動(dòng)，采用穩(wěn)態(tài)流動(dòng)的CFD數(shù)值模擬方法，計(jì)算過(guò)程具體如下：通過(guò)復(fù)合飛行器系統(tǒng)三維實(shí)體建模得到定常流動(dòng)計(jì)算所需的網(wǎng)格文件和邊界條件，采用ANSYS R15.0 Fluent軟件進(jìn)行數(shù)值模擬迭代求解；復(fù)合飛行器巡航速度為聲速，因此采用壓力基求解器進(jìn)行穩(wěn)態(tài)流動(dòng)迭代求解。湍流模型采用航空外流場(chǎng)計(jì)算常用的對(duì)受到逆壓梯度作用下的邊界層模擬效果較好的Spalart-Allmaras模型，求解算法使用Simple算法；計(jì)算入口邊界條件采用壓力遠(yuǎn)場(chǎng)，采用二階迎風(fēng)格式對(duì)控制方程進(jìn)行離散推進(jìn)，并最終通過(guò)監(jiān)測(cè)復(fù)合機(jī)翼上下表面壓力平均值的穩(wěn)定性作為計(jì)算結(jié)果的收斂判定，復(fù)合飛行器系統(tǒng)非結(jié)構(gòu)網(wǎng)格具體如圖6所示。建立氣動(dòng)力數(shù)據(jù)庫(kù)時(shí)，復(fù)合飛行器系統(tǒng)的靜導(dǎo)數(shù)和操縱導(dǎo)數(shù)計(jì)算可以認(rèn)為是定常流動(dòng)，采用穩(wěn)態(tài)流動(dòng)的CFD數(shù)值模擬方法，計(jì)算入口邊界條件采用壓力遠(yuǎn)場(chǎng)，采用二階迎風(fēng)格式對(duì)控制方程進(jìn)行離散推進(jìn)，并最終通過(guò)監(jiān)測(cè)復(fù)合機(jī)翼上下表面壓力平均值的穩(wěn)定性作為計(jì)算結(jié)果的收斂判定。動(dòng)導(dǎo)數(shù)CLq、CL˙α、Cmq、Cm˙α等主要由Zaero軟件的非定常氣動(dòng)力計(jì)算方法獲得。

圖6 復(fù)合飛行器系統(tǒng)ANSYS非結(jié)構(gòu)網(wǎng)格劃分Fig.6 Unstructured grid generation of composite aircraft system by ANSYS

經(jīng)過(guò)基本氣動(dòng)參數(shù)、靜導(dǎo)數(shù)、動(dòng)導(dǎo)數(shù)、各舵面偏轉(zhuǎn)角度與復(fù)合運(yùn)動(dòng)特征參數(shù)即相對(duì)滾轉(zhuǎn)姿態(tài)的CFD仿真計(jì)算后，就可以建立得到考慮氣動(dòng)舵面偏轉(zhuǎn)和復(fù)合運(yùn)動(dòng)特征參數(shù)的復(fù)合飛行氣動(dòng)數(shù)據(jù)庫(kù)。

3 動(dòng)力學(xué)特性

3.1 仿真方法

根據(jù)第2節(jié)中推導(dǎo)得到的7自由度復(fù)合機(jī)體動(dòng)力學(xué)方程全量模型，在MATLAB／Simulink平臺(tái)上搭建復(fù)合飛行器系統(tǒng)動(dòng)力學(xué)仿真模塊，如圖7所示。其中飛機(jī)多體系統(tǒng)狀態(tài)模塊主要是計(jì)算復(fù)合飛行器系統(tǒng)的增廣體靜矩、增廣體張量和增廣體質(zhì)量；飛行器動(dòng)力學(xué)模型計(jì)算模塊主要負(fù)責(zé)根據(jù)飛機(jī)多體狀態(tài)參數(shù)、運(yùn)動(dòng)姿態(tài)參數(shù)和舵面輸入來(lái)確定動(dòng)力學(xué)方程中力和力矩的計(jì)算；多剛體系統(tǒng)非線性方程組求解器主要負(fù)責(zé)計(jì)算下一步計(jì)算過(guò)程中的運(yùn)動(dòng)參數(shù)和姿態(tài)參數(shù)；控制系統(tǒng)反饋模塊主要負(fù)責(zé)根據(jù)選取好的復(fù)合飛行器配平方案，給定對(duì)應(yīng)舵面偏轉(zhuǎn)參數(shù)對(duì)應(yīng)的數(shù)值。迭代求解器采用ODE45格式。

圖7 翼尖鉸接復(fù)合飛行器系統(tǒng)動(dòng)力學(xué)仿真模塊Fig.7 Dynamic simulation module of wingtip-jointed composite aircraft system

使用全量方程搭建復(fù)合飛行器系統(tǒng)的動(dòng)力學(xué)仿真平臺(tái)，可以開(kāi)展未配平、準(zhǔn)配平和全配平方案下的動(dòng)力學(xué)響應(yīng)研究。未配平方案指復(fù)合飛行器各控制舵面（主要有油門舵偏角、升降舵偏角、方向舵偏角和副翼偏角等）輸入量為任意值，不需要使復(fù)合飛行器整體運(yùn)動(dòng)和內(nèi)部相對(duì)運(yùn)動(dòng)動(dòng)力學(xué)方程平衡的舵面控制輸入方案；準(zhǔn)配平方案指復(fù)合飛行器各控制舵面需要滿足復(fù)合飛行器系統(tǒng)整體運(yùn)動(dòng)動(dòng)力學(xué)方程平衡，不需要滿足內(nèi)部相對(duì)運(yùn)動(dòng)動(dòng)力學(xué)方程平衡的舵面控制輸入方案；全配平方案指復(fù)合飛行器系統(tǒng)各舵面控制輸入量需要同時(shí)滿足復(fù)合飛行器系統(tǒng)整體運(yùn)動(dòng)和內(nèi)部相對(duì)運(yùn)動(dòng)動(dòng)力學(xué)方程的平衡。

3.2 準(zhǔn)配平方案和全配平方案的動(dòng)力學(xué)響應(yīng)

在算例工況高度5000 m，馬赫數(shù)0.282（速度95.91 m／s），相對(duì)滾轉(zhuǎn)角0°和迎角2.27°條件下，研究準(zhǔn)配平方案和全配平方案下的雙機(jī)復(fù)合飛行器系統(tǒng)算例的動(dòng)力學(xué)響應(yīng)。

準(zhǔn)配平方案具體為：為保持復(fù)合飛行器整體處于定直平飛狀態(tài)，使用升降舵偏角、油門桿位移來(lái)平衡復(fù)合飛行器系統(tǒng)任務(wù)飛行高度、速度，同時(shí)副翼輸入來(lái)平衡復(fù)合機(jī)體系統(tǒng)的滾轉(zhuǎn)力矩，方向舵來(lái)配平復(fù)合飛行器系統(tǒng)的偏航力矩。按以上策略在算例工況下實(shí)現(xiàn)準(zhǔn)配平方案如表2所示，保證復(fù)合飛行器整體動(dòng)力學(xué)方程中所有初始狀態(tài)導(dǎo)數(shù)（˙xi除外）小于10－8，則檢查的所有配置都滿足準(zhǔn)配平要求；仿真開(kāi)始后飛機(jī)控制舵面在整體運(yùn)動(dòng)方程中始終保持系統(tǒng)初始穩(wěn)定狀態(tài)的輸入。

全配平配平方案具體為：為保持復(fù)合飛行器整體處于定直平飛狀態(tài)，使用升降舵偏角、油門桿位移來(lái)平衡復(fù)合飛行器系統(tǒng)任務(wù)飛行高度、速度。使用副翼輸入來(lái)平衡復(fù)合機(jī)體系統(tǒng)的滾轉(zhuǎn)力矩同時(shí)平衡系統(tǒng)內(nèi)部的相對(duì)滾轉(zhuǎn)力矩，方向舵來(lái)配平復(fù)合飛行器系統(tǒng)的偏航力矩。按以上策略在算例工況下實(shí)現(xiàn)全配平方案如表2所示，保證復(fù)合飛行器整體動(dòng)力學(xué)方程與內(nèi)部相對(duì)轉(zhuǎn)動(dòng)動(dòng)力學(xué)方程中所有初始狀態(tài)導(dǎo)數(shù)（˙xi除外）小于10－8，則檢查的所有配置都滿足全配平要求。仿真開(kāi)始后飛機(jī)控制輸入在整體運(yùn)動(dòng)方程暫不調(diào)整，始終保持系統(tǒng)初始穩(wěn)定狀態(tài)的輸入，仿真時(shí)間10 s，仿真步長(zhǎng)間隔0.01 s，仿真結(jié)果如圖8所示。

表2 復(fù)合飛行配平方案Table 2 Composite flight trimming strategy

由準(zhǔn)配平方案下復(fù)合飛行動(dòng)力學(xué)響應(yīng)仿真結(jié)果說(shuō)明，復(fù)合飛行器系統(tǒng)動(dòng)力學(xué)仿真在4.58 s終止，原因是相對(duì)滾轉(zhuǎn)角達(dá)到180°，即2架復(fù)合飛機(jī)通過(guò)翼尖鉸接機(jī)構(gòu)最終在4.58 s以相對(duì)滾轉(zhuǎn)運(yùn)動(dòng)形式相撞。不同于單機(jī)飛行動(dòng)力學(xué)響應(yīng)，雙機(jī)復(fù)合飛行會(huì)存在相對(duì)滾轉(zhuǎn)力矩造成內(nèi)部相對(duì)滾轉(zhuǎn)運(yùn)動(dòng)，需要通過(guò)偏轉(zhuǎn)副翼、升降舵和方向舵實(shí)現(xiàn)系統(tǒng)內(nèi)部動(dòng)力學(xué)方程中的相對(duì)滾轉(zhuǎn)力矩。差動(dòng)雙機(jī)升降舵會(huì)造成多余的俯仰力矩，使得復(fù)合飛行器配平問(wèn)題更加復(fù)雜，差動(dòng)方向舵偏角會(huì)造成內(nèi)部相對(duì)偏航力矩，從而對(duì)飛行器的翼尖鉸接結(jié)構(gòu)強(qiáng)度和機(jī)構(gòu)對(duì)接精度產(chǎn)生更高要求。因此，全配平方案主要使用差動(dòng)副翼舵面偏角來(lái)配平內(nèi)部相對(duì)滾轉(zhuǎn)力矩，其中造成運(yùn)動(dòng)發(fā)散的系統(tǒng)內(nèi)部相對(duì)滾轉(zhuǎn)力矩（l2－l1）可以由面向復(fù)合運(yùn)動(dòng)特征參數(shù)和舵面部件控制參數(shù)的氣動(dòng)力數(shù)據(jù)庫(kù)得到。仿真結(jié)果如圖8所示，馬赫數(shù)、迎角、俯仰角和相對(duì)滾轉(zhuǎn)角的發(fā)散速度均小于準(zhǔn)配平方案仿真結(jié)果，同時(shí)包括馬赫數(shù)、迎角、俯仰角和相對(duì)滾轉(zhuǎn)角在內(nèi)的整體和系統(tǒng)內(nèi)部運(yùn)動(dòng)零階參數(shù)在仿真時(shí)間10 s內(nèi)的變化幅值均小于0.3%。

圖8 翼尖鉸接復(fù)合飛行器動(dòng)力學(xué)響應(yīng)仿真結(jié)果Fig.8 Dynamic response simulation results of wingtip-jointed composite aircraft

3.3 本征動(dòng)力學(xué)特性分析

在全配平方案狀態(tài)下，馬赫數(shù)、迎角、俯仰角和相對(duì)滾轉(zhuǎn)角在內(nèi)的整體和系統(tǒng)內(nèi)部運(yùn)動(dòng)零階參數(shù)在仿真時(shí)間10 s內(nèi)的數(shù)值趨于穩(wěn)定，以此為基準(zhǔn)運(yùn)動(dòng)狀態(tài)，使用小擾動(dòng)假設(shè)［18］下開(kāi)展全量非解耦的線性化工作，具體包括飛行器系統(tǒng)力和力矩的線化和飛行力學(xué)方程組的線化，再將線性化結(jié)果整理成狀態(tài)空間方程有

求解飛行器系統(tǒng)運(yùn)動(dòng)狀態(tài)方程的系統(tǒng)矩陣A的特征值和特征向量可得雙機(jī)復(fù)合運(yùn)動(dòng)模態(tài)，通過(guò)特征值和特征向量的分析，并與單機(jī)運(yùn)動(dòng)模態(tài)對(duì)照可得如表3所示的雙機(jī)復(fù)合運(yùn)動(dòng)模態(tài)分布規(guī)律。相比單機(jī)運(yùn)動(dòng)模態(tài)，雙機(jī)復(fù)合運(yùn)動(dòng)模態(tài)多出了2個(gè)發(fā)散的相對(duì)滾轉(zhuǎn)模態(tài)，如表3中命名的“復(fù)合運(yùn)動(dòng)模態(tài)1”和“復(fù)合運(yùn)動(dòng)模態(tài)2”。雙機(jī)飛行時(shí)常規(guī)運(yùn)動(dòng)模態(tài)和單機(jī)飛行時(shí)的模態(tài)特性變化規(guī)律可以根據(jù)各模態(tài)特征值幅值和相位值分析得出為比較雙機(jī)復(fù)合飛行整體系統(tǒng)模態(tài)和單機(jī)飛行時(shí)的變化規(guī)律，給出了各運(yùn)動(dòng)模態(tài)參數(shù)指標(biāo)如表3所示，具體為ωn、ε、ω、T、t1／2、N1／2、t2等特性指標(biāo)的數(shù)值變化；為便于分析2個(gè)發(fā)散模態(tài)的運(yùn)動(dòng)特點(diǎn)，在這里分別給出5個(gè)模態(tài)的特征向量：復(fù)合運(yùn)動(dòng)模態(tài)1的特征矢量：［1.67×10－80.00252.19×10－80.0001－7.88×10－11－1.03×10－10－2.55×10－91.16×10－70.99990.0019］T；復(fù)合運(yùn)動(dòng)模態(tài)2的特征矢量：［0.0024－0.00090.00010.00022.02×10－91.19×10－10－5.6×10－58.42×10－50.15530.9878］T，其中復(fù)合運(yùn)動(dòng)模態(tài)1的特征矢量：忽略數(shù)量在10－4以下的參數(shù)數(shù)值，則相關(guān)參數(shù)之比α∶θ∶φ12∶˙φ12＝0.0025∶0.0001∶1∶0.0019。復(fù)合運(yùn)動(dòng)模態(tài)2的特征矢量：忽略數(shù)量在10－4以下的參數(shù)數(shù)值，則相關(guān)參數(shù)之比為V∶α∶q∶θ∶φ12∶˙φ12＝0.0024∶－0.0009∶0.0001∶0.0002∶0.1553∶1。

結(jié)合表3、表4中復(fù)合運(yùn)動(dòng)模態(tài)1和復(fù)合運(yùn)動(dòng)模態(tài)2的特征值和模態(tài)特性指標(biāo)參數(shù)，可以分析得出復(fù)合運(yùn)動(dòng)模態(tài)1是相對(duì)滾轉(zhuǎn)角緩慢發(fā)散運(yùn)動(dòng)主導(dǎo)的同時(shí)伴有微小迎角、俯仰角速度和俯仰角發(fā)散的模態(tài)，倍幅時(shí)間t2為4.813 s；復(fù)合運(yùn)動(dòng)模態(tài)2是相對(duì)滾轉(zhuǎn)角速度快速發(fā)散運(yùn)動(dòng)主導(dǎo)的，同時(shí)伴有緩慢的相對(duì)滾轉(zhuǎn)角與微小的速度、迎角、俯仰角速度和俯仰角發(fā)散運(yùn)動(dòng)的模態(tài)，倍幅時(shí)間t2為0.407 s。復(fù)合飛行時(shí)各模態(tài)指標(biāo)參數(shù)也存在顯著變化：長(zhǎng)周期的無(wú)阻尼自振頻率ωn從0.050 rad／s增加到0.054 rad／s，阻尼比ε從0.154減小到0.080，振蕩頻率從0.0497 rad／s增加到0.0534 rad／s，周期T從126.42 s減小到117.66 s，半幅時(shí)間t1／2從90 s提高到161.16 s，半衰期內(nèi)振蕩次數(shù)N1／2從0.72次增加到1.37次；縱向模態(tài)1的無(wú)阻尼自振頻率ωn從2.3 rad／s減小到2.27 rad／s，半幅時(shí)間t1／2從0.301 s增加到0.305 s；縱向模態(tài)2的無(wú)阻尼自振頻率ωn從0.0552 rad／s減小到0.0486 rad／s，半幅時(shí)間t1／2從12.53 s增加到14.26 s；滾轉(zhuǎn)模態(tài)的無(wú)阻尼自振頻率ωn從4.169 rad／s增加到6.658 rad／s，半幅時(shí)間t1／2從0.1662 s減小到0.1041 s；螺旋模態(tài)的無(wú)阻尼自振頻率ωn從0.022 rad／s減小到0.008 rad／s，半幅時(shí)間 t1／2從31.5 s增加到86.625 s；荷蘭滾模態(tài)的無(wú)阻尼自振頻率ωn從0.0169 rad／s增加到0.163 rad／s，阻尼比ε從0.722減小到0.564，振蕩頻率從0.0118 rad／s增加到0.1346 rad／s，周期T從532.5 s減小到46.7 s，半幅時(shí)間t1／2從56.8 s減小到7.541 s，半衰期內(nèi)振蕩次數(shù)N1／2從0.1067次增加到0.1615次。

表3 運(yùn)動(dòng)模態(tài)分析Table 3 Motion modal analysis

表4 模態(tài)參數(shù)分析Table 4 Modal parameter analysis

由于發(fā)散的復(fù)合合運(yùn)動(dòng)模態(tài)2的倍幅時(shí)間很短，對(duì)飛行穩(wěn)定性影響很大，飛行器在無(wú)控狀態(tài)下無(wú)法穩(wěn)定飛行，這與傳統(tǒng)構(gòu)型飛行器有很大差別，因此為了保障復(fù)合飛行器系統(tǒng)能夠長(zhǎng)期持續(xù)地穩(wěn)定飛行，需要采用合理的控制策略來(lái)保證飛行的穩(wěn)定性。

4 結(jié) 論

1）根據(jù)推導(dǎo)得到的7自由度復(fù)合機(jī)體動(dòng)力學(xué)方程全量模型和建模完成的氣動(dòng)力數(shù)據(jù)庫(kù)的動(dòng)力學(xué)仿真結(jié)果證明類似單機(jī)飛行配平方案的準(zhǔn)配平方案無(wú)法使復(fù)合飛行器系統(tǒng)長(zhǎng)期保持飛行穩(wěn)定，而考慮整體和內(nèi)部動(dòng)力學(xué)方程組和運(yùn)動(dòng)學(xué)方程組平衡的全配平方案下復(fù)合系統(tǒng)的運(yùn)動(dòng)參數(shù)在仿真時(shí)間10 s內(nèi)始終趨于穩(wěn)定，各參數(shù)變化幅度均小于5%。

2）在全配平方案下，使用航空飛行器動(dòng)力學(xué)原理的小擾動(dòng)假設(shè)，對(duì)復(fù)合飛行器動(dòng)力學(xué)方程組和運(yùn)動(dòng)學(xué)方程組開(kāi)展線性化工作，整理得到復(fù)合系統(tǒng)狀態(tài)空間方程和系統(tǒng)本征矩陣，分析得到了系統(tǒng)本征矩陣的特征值和特征向量，得到了雙機(jī)復(fù)合飛行模態(tài)分布規(guī)律，并研究分析了單機(jī)飛行模態(tài)分布的指標(biāo)參數(shù)數(shù)值的變化。

3）研究復(fù)合飛行過(guò)程出現(xiàn)的2個(gè)發(fā)散的復(fù)雜運(yùn)動(dòng)模態(tài)，分別由相對(duì)滾轉(zhuǎn)角和相對(duì)滾轉(zhuǎn)角速度主導(dǎo)，因此需要設(shè)計(jì)合理的控制方法來(lái)保持復(fù)合飛行的穩(wěn)定性。