一種適用于彈性飛機(jī)飛行仿真的補(bǔ)丁方法

王培涵，吳志剛，，楊超，孫曉旭

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

2．海鷹航空通用裝備有限責(zé)任公司，北京 100070

飛行仿真技術(shù)是以飛行器的運(yùn)動(dòng)為研究對(duì)象，面向復(fù)雜系統(tǒng)的仿真分析方法[1-2］，是計(jì)算飛行力學(xué)研究的重要領(lǐng)域。傳統(tǒng)飛行力學(xué)主要研究低速飛行器的運(yùn)動(dòng)特性和安全性，并逐步衍生出有控飛行力學(xué)、計(jì)算飛行力學(xué)等諸多領(lǐng)域[3-4］。在此基礎(chǔ)上發(fā)展出的飛行仿真技術(shù)，大多基于剛性假設(shè)，在運(yùn)動(dòng)響應(yīng)分析、飛行控制系統(tǒng)設(shè)計(jì)等方面展開(kāi)研究。

隨著航空技術(shù)的飛速發(fā)展，現(xiàn)代飛機(jī)呈現(xiàn)出輕結(jié)構(gòu)、大柔性和低阻尼等特點(diǎn)[5］，氣動(dòng)彈性問(wèn)題成為不可忽視的重要因素，使得飛行仿真中的剛性假設(shè)在特定情況下不再適用。彈性飛機(jī)飛行仿真的關(guān)鍵，在于剛體運(yùn)動(dòng)與彈性變形之間的耦合關(guān)系，需考慮結(jié)構(gòu)的彈性效應(yīng)、不同模態(tài)間的氣動(dòng)力耦合等諸多因素，情況較復(fù)雜。國(guó)內(nèi)外學(xué)者大多將剛體自由度和彈性自由度統(tǒng)一考慮，選取適當(dāng)?shù)淖鴺?biāo)系構(gòu)建彈性飛機(jī)運(yùn)動(dòng)方程，包括平均體軸系、準(zhǔn)坐標(biāo)系及瞬態(tài)坐標(biāo)系等[6］。

目前，平均體軸系在彈性飛機(jī)飛行仿真研究中應(yīng)用較為廣泛。Waszak 等[7-8］通過(guò)拉格朗日方程和片條理論構(gòu)建彈性飛機(jī)的運(yùn)動(dòng)方程，實(shí)現(xiàn)了基于剛彈耦合的彈性飛機(jī)飛行仿真，為后續(xù)研究工作提供了理論基礎(chǔ)。Schmidt 等[9］將彈性自由度添加到剛性飛機(jī)動(dòng)力學(xué)模型中，統(tǒng)一考慮剛體和彈性模態(tài)廣義坐標(biāo)，并以B-1 和高速商用運(yùn)輸機(jī)為例，深入研究了基于剛彈耦合的彈性飛機(jī)飛行仿真。Saltari 等[10］推導(dǎo)了全耦合方程，構(gòu)建了配平狀態(tài)下小擾動(dòng)線(xiàn)化氣動(dòng)彈性的狀態(tài)空間方程，系統(tǒng)地討論了慣性和氣動(dòng)耦合效應(yīng)對(duì)穩(wěn)定性、動(dòng)態(tài)響應(yīng)的影響。Liu 等[11］基于非線(xiàn)性耦合動(dòng)力學(xué)方程對(duì)彈性飛機(jī)進(jìn)行了仿真，研究了幾何非線(xiàn)性對(duì)彈性飛機(jī)運(yùn)動(dòng)方程的影響。Zú?iga等[12］針對(duì)彈性飛行器進(jìn)行了基于剛彈耦合的飛行仿真，深入探討了彈性效應(yīng)對(duì)短周期和長(zhǎng)周期模態(tài)的影響。孫若琳[13］基于氣動(dòng)力降階模型開(kāi)展了彈性飛機(jī)飛行仿真研究，師妍等[14-15］進(jìn)一步研究適用于飛行仿真的非線(xiàn)性、非定常氣動(dòng)力降階建模，實(shí)現(xiàn)了低速飛機(jī)在強(qiáng)陣風(fēng)作用下的飛行仿真分析。

由于剛彈耦合動(dòng)力學(xué)方程需統(tǒng)一考慮剛體和彈性的自由度，有學(xué)者提出了將剛體和彈性模態(tài)廣義坐標(biāo)分開(kāi)考慮，以簡(jiǎn)化彈性飛機(jī)飛行仿真的建模過(guò)程。Looye[16］結(jié)合氣動(dòng)模型，提出了僅考慮動(dòng)力學(xué)增量對(duì)剛性飛機(jī)運(yùn)動(dòng)方程影響的研究思路。Kier 與Looye[17］合作，在小擾動(dòng)線(xiàn)化模型中，利用該方法改進(jìn)非線(xiàn)性運(yùn)動(dòng)方程的動(dòng)態(tài)機(jī)動(dòng)載荷分析模型，形成了統(tǒng)一的機(jī)動(dòng)和陣風(fēng)載荷建模分析方法。Castrichini 等[18］應(yīng)用該思想，研究半氣動(dòng)彈性鉸鏈裝置對(duì)氣動(dòng)彈性、飛行動(dòng)力學(xué)耦合效應(yīng)的影響。Chen 等[19］以此思路為基礎(chǔ)，設(shè)計(jì)了非線(xiàn)性動(dòng)態(tài)飛行仿真（Dynamic Flight Simulation，DFS）工具，結(jié)合飛行動(dòng)力學(xué)和氣動(dòng)彈性力學(xué)，實(shí)現(xiàn)了基于小擾動(dòng)線(xiàn)化方程的彈性飛機(jī)飛行仿真。而后，美國(guó)空軍研究實(shí)驗(yàn)室在StS-RtS 仿真程序中應(yīng)用該飛行仿真工具，提高了載荷計(jì)算精度。Tuegel 根據(jù)DFS 便于集成的特點(diǎn)，在計(jì)算飛機(jī)壽命周期的數(shù)字孿生技術(shù)中，應(yīng)用該方法考慮氣動(dòng)彈性對(duì)飛機(jī)疲勞載荷的影響，提高了壽命周期的計(jì)算精度[20］。

與剛性飛機(jī)的飛行仿真相比，基于剛彈耦合的彈性飛機(jī)建模方法相對(duì)復(fù)雜，使得動(dòng)力學(xué)特性、控制系統(tǒng)設(shè)計(jì)等的研究和驗(yàn)證工作需要重新開(kāi)發(fā)，不利于其在實(shí)際工程中的應(yīng)用。分開(kāi)考慮剛體和彈性自由度，通過(guò)引入氣動(dòng)彈性效應(yīng)將原有剛性飛機(jī)飛行仿真推廣到彈性飛機(jī)，具有便于集成的特點(diǎn)。但該方法基于小擾動(dòng)線(xiàn)化模型，未推廣至非線(xiàn)性6 自由度全量方程。而實(shí)際工程應(yīng)用中，全量方程的剛性飛機(jī)飛行仿真占比較大，該方法的應(yīng)用受到了限制。另外，這種分開(kāi)考慮的思想在國(guó)內(nèi)外研究中相對(duì)較少，具有深入分析的價(jià)值。

針對(duì)以上這些問(wèn)題，本文提出了一種適用于彈性飛機(jī)飛行仿真的補(bǔ)丁方法。結(jié)合剛彈耦合動(dòng)力學(xué)方程的特點(diǎn)，拆分廣義氣動(dòng)力、提取出與彈性效應(yīng)相關(guān)的部分疊加到剛性飛機(jī)運(yùn)動(dòng)方程中，以實(shí)現(xiàn)彈性飛機(jī)運(yùn)動(dòng)方程的求解；同時(shí)將彈性振動(dòng)導(dǎo)致的測(cè)量信號(hào)增量疊加到剛性飛機(jī)運(yùn)動(dòng)響應(yīng)中，實(shí)現(xiàn)了彈性飛機(jī)測(cè)量信號(hào)的模擬。這兩部分共同構(gòu)造了反映氣動(dòng)彈性效應(yīng)的“補(bǔ)丁模塊”。本文的補(bǔ)丁方法通過(guò)置入該模塊，將6 自由度全量方程剛性飛機(jī)的飛行仿真推廣到彈性飛機(jī)的飛行仿真。該方法能充分利用原有剛性飛機(jī)仿真模型，簡(jiǎn)化了彈性飛機(jī)的建模和驗(yàn)證過(guò)程，便于穩(wěn)定性分析和控制系統(tǒng)設(shè)計(jì)等工作開(kāi)展。

1 飛行仿真補(bǔ)丁方法總體思路

剛性飛機(jī)[2，21］和基于剛彈耦合的彈性飛機(jī)飛行仿真[7-8］，分別是采用6 自由度剛性飛機(jī)運(yùn)動(dòng)方程和剛彈耦合動(dòng)力學(xué)方程來(lái)實(shí)現(xiàn)的，如圖1 和圖2所示。

圖2 基于剛彈耦合的彈性飛機(jī)飛行仿真Fig.2 Flight simulation of flexible aircraft based on rigid-elastic coupling equations

相比于6 自由度剛性飛機(jī)運(yùn)動(dòng)方程，圖2 中的剛彈耦合動(dòng)力學(xué)方程的特點(diǎn)在于引入了n個(gè)彈性自由度，方程總個(gè)數(shù)變?yōu)?+n。該變化使得飛行仿真從剛性飛機(jī)推廣到彈性飛機(jī)時(shí)，飛行動(dòng)力學(xué)模塊的仿真程序需要做較大幅度改動(dòng)。

通過(guò)分析剛彈耦合方程，考慮將廣義氣動(dòng)力拆分，把飛機(jī)的彈性效應(yīng)整理為氣動(dòng)彈性狀態(tài)空間方程、剛體氣動(dòng)力的增量等幾部分，并將它們組合成“補(bǔ)丁模塊”。將補(bǔ)丁模塊與已有的剛性飛機(jī)飛行仿真程序相結(jié)合，實(shí)現(xiàn)了彈性飛機(jī)的飛行仿真，飛行仿真方法的總體思路如圖3 所示。圖3 中的“剛性飛機(jī)運(yùn)動(dòng)方程”模塊[22-23］已是成熟技術(shù)，這里不作贅述，本文重點(diǎn)討論“補(bǔ)丁模塊”的構(gòu)建。

圖3 基于“補(bǔ)丁模塊”的彈性飛機(jī)飛行仿真示意圖Fig.3 Schematic diagram of flight simulation of flexible aircraft based on patch module method

2 彈性飛機(jī)動(dòng)力學(xué)模型

在彈性飛機(jī)的飛行仿真中，選擇平均體軸系，可減弱或消除剛體和彈性自由度的慣性耦合，飛機(jī)的剛彈耦合方程可拆分成剛體運(yùn)動(dòng)和彈性振動(dòng)方程，其中的氣動(dòng)彈性效應(yīng)可通過(guò)彈性振動(dòng)方程和廣義氣動(dòng)力表征。

在平均體軸系下，應(yīng)用拉格朗日方程推導(dǎo)彈性飛機(jī)的剛彈耦合動(dòng)力學(xué)方程[7-9］，表達(dá)式為

式（1）為剛體動(dòng)力學(xué)方程：M為飛機(jī)的質(zhì)量；Ixx、Iyy、Izz為 飛 機(jī) 繞 體 軸 系3 個(gè) 軸 轉(zhuǎn) 動(dòng) 的 轉(zhuǎn) 動(dòng) 慣量；Ixy、Ixz、Iyz為交叉慣性積；U、V、W為飛行速度在體軸系3 個(gè)方向的投影；p、q、r為角速度在體軸系 上 的3 個(gè) 分 量；φ、θ、ψ為 歐 拉 角；Fx、Fy、Fz和Mx、My、Mz為除重力外，飛機(jī)所受的外力和外力矩在體軸系上的分量。式（2）為彈性振動(dòng)方程：Mi、ξi、μi、ωi和Qξi分別為第i階彈性模態(tài)的廣義質(zhì)量、廣義坐標(biāo)、阻尼比、無(wú)阻尼自然頻率和廣義力；n為彈性模態(tài)的階數(shù)。整理式（1）和式（2）的右端項(xiàng)，得到廣義力的表達(dá)式為

飛機(jī)所受的外力，包括重力、發(fā)動(dòng)機(jī)推力、舵面偏轉(zhuǎn)引起的慣性力，以及空氣動(dòng)力。重力已在式（1）左端給出，廣義力的表達(dá)式為

式 中：T為 發(fā) 動(dòng) 機(jī) 推 力 在 廣 義 力 中 的 體 現(xiàn)；Mqδ為舵面耦合慣性質(zhì)量為舵面偏轉(zhuǎn)角的二階導(dǎo)數(shù)；fA為廣義氣動(dòng)力，由剛體運(yùn)動(dòng)和彈性振動(dòng)耦合而成。

從式（1）和式（2）可以看出，在平均體軸系下，剛體與彈性自由度在慣性和剛度上不存在耦合，它們的耦合主要體現(xiàn)在廣義氣動(dòng)力上，即為氣動(dòng)耦合[10］。

3 廣義氣動(dòng)力拆分

根據(jù)剛彈耦合動(dòng)力學(xué)方程，廣義氣動(dòng)力由剛體運(yùn)動(dòng)和彈性振動(dòng)耦合產(chǎn)生。因此，需定義表征剛體運(yùn)動(dòng)和彈性振動(dòng)的變量。其中，剛體運(yùn)動(dòng)由狀態(tài)變量η進(jìn)行表征，表達(dá)式為η=

式中：X、Y為體軸系沿x、y軸方向的位移；H為高度的變化；α、β分別為迎角和側(cè)滑角。

彈性振動(dòng)由彈性模態(tài)廣義坐標(biāo)ξ定義，表達(dá)式為

不失一般性，式（4）中的廣義氣動(dòng)力可表示為非線(xiàn)性形式，其表達(dá)式為

式中：δ為舵面偏轉(zhuǎn)角；fη和fξ為剛體氣動(dòng)力、彈性模態(tài)廣義氣動(dòng)力，其與η、ξ、δ存在非線(xiàn)性的函數(shù)關(guān)系。

在6 自由度全量方程的剛性飛機(jī)飛行仿真中，氣動(dòng)力一般采用非線(xiàn)性形式。而對(duì)于一般的彈性飛機(jī)，在線(xiàn)彈性小變形的假設(shè)條件下，彈性模態(tài)產(chǎn)生的廣義氣動(dòng)力可以表示成線(xiàn)性形式。故將廣義氣動(dòng)力拆分成非線(xiàn)性部分和線(xiàn)性增量的形式，其表達(dá)式為

式 中：h1(η，δ) 為 剛 性 飛 機(jī) 非 線(xiàn) 性 氣 動(dòng) 力；h2(η，δ)為剛體運(yùn)動(dòng)對(duì)彈性模態(tài)廣義氣動(dòng)力的影響，同樣為非線(xiàn)性形式；Δfη(ξ)和Δfξ(ξ)分別表示剛體氣動(dòng)力增量和彈性模態(tài)廣義氣動(dòng)力增量，它們可由偶極子格網(wǎng)法[24］計(jì)算、并通過(guò)有理函數(shù)擬合技術(shù)[25-26］轉(zhuǎn)化為時(shí)域形式，得到的表達(dá)式為

式中：R為對(duì)角陣，其對(duì)角線(xiàn)元素為氣動(dòng)力滯后根；Eξ為彈性模態(tài)廣義坐標(biāo)導(dǎo)數(shù)的系數(shù)矩陣。

拆分后的廣義氣動(dòng)力中，h1(η，δ)和h2(η，δ)為非線(xiàn)性部分，Δfη和Δfξ為線(xiàn)性部分。

3.1 剛性飛機(jī)非線(xiàn)性氣動(dòng)力h1

在原有6 自由度全量方程的剛性飛機(jī)飛行仿真中，剛性飛機(jī)的氣動(dòng)力與狀態(tài)變量η的非線(xiàn)性關(guān)系，一般不采用表達(dá)式的形式，而是根據(jù)氣動(dòng)力插值表等離散形式的數(shù)據(jù)庫(kù)[2，21］，通過(guò)線(xiàn)性插值得到。這些數(shù)據(jù)庫(kù)通常由CFD 計(jì)算、風(fēng)洞試驗(yàn)或飛行試驗(yàn)等方式獲得。以法向力系數(shù)Cn為例，給出氣動(dòng)力插值表示例如表1 所示。

表1 氣動(dòng)力Cn插值表示例Table 1 Examples of aerodynamic Cninterpolation

在飛行仿真過(guò)程中，根據(jù)每一時(shí)間步的飛行狀態(tài)，例如迎角α和側(cè)滑角β，可從氣動(dòng)力插值表中得到當(dāng)前的氣動(dòng)力數(shù)據(jù)，用于剛性飛機(jī)運(yùn)動(dòng)方程求解。

3.2 彈性模態(tài)廣義力h2 計(jì)算與修正

與h1有所不同，彈性模態(tài)廣義氣動(dòng)力的非線(xiàn)性部分h2與彈性模態(tài)振型和氣動(dòng)力分布有關(guān)，一般無(wú)法通過(guò)試驗(yàn)獲得，采用數(shù)值計(jì)算得到。若采用CFD 方法計(jì)算，時(shí)間較長(zhǎng)，數(shù)據(jù)量也非常龐大，難以用于飛行仿真，因此本文采用面元法[27］計(jì)算。

在氣動(dòng)力計(jì)算時(shí)，需要給定飛機(jī)6 個(gè)剛體運(yùn)動(dòng)模態(tài)，包括3 個(gè)剛體平動(dòng)模態(tài)Tx、Ty、Tz和3 個(gè)繞飛機(jī)質(zhì)心的轉(zhuǎn)動(dòng)模態(tài)Rx、Ry、Rz，得到關(guān)于剛體模態(tài)的氣動(dòng)力影響系數(shù)矩陣。剛體模態(tài)與狀態(tài)變量η之間的關(guān)系[28］為

需要說(shuō)明的是：式（12）、式（13）實(shí)際上假設(shè)了φ、θ等為小量，不適合大轉(zhuǎn)動(dòng)等情況。最終可以得到由η和δ產(chǎn)生的剛體氣動(dòng)力和彈性模態(tài)廣義氣動(dòng)力，這里記為和，它們的線(xiàn)性形式為

式中：為動(dòng)壓；Aηη、Aηξ、Aδη、Aδξ為 對(duì)應(yīng)于η和δ的氣動(dòng)力影響系數(shù)矩陣，上標(biāo)η和δ為矩陣對(duì)應(yīng)的變量，下標(biāo)η或ξ為對(duì)應(yīng)剛體氣動(dòng)力或彈性模態(tài)廣義氣動(dòng)力；hη0和hξ0分別為初始下洗角產(chǎn)生的氣動(dòng)力分量。

值得注意的是：在實(shí)際飛行仿真中使用的剛性飛機(jī)氣動(dòng)力h1(η，δ)具有非線(xiàn)性特性，與h1在數(shù)值上不匹配，需要對(duì)和修正。

假定由面元法計(jì)算的氣動(dòng)力分布在比例上是準(zhǔn)確的，則修正思路為：將與h1相比較求得修正因子，再將修正因子應(yīng)用于，從而得到飛行仿真中需要的h2(η，δ)。

以迎角α產(chǎn)生的氣動(dòng)力Cn為例，的氣動(dòng)力修正如圖4 所示。由面元法計(jì)算得到的零舵偏時(shí)剛體氣動(dòng)力在體軸系的z向分量為初始下洗角產(chǎn)生的貢獻(xiàn)為hz0。

首先，在剛性飛機(jī)非線(xiàn)性氣動(dòng)力h1(η，δ)中，獲得零迎角、零舵偏時(shí)的法向力系數(shù)Cn0，根據(jù)Cn0與hz0的關(guān)系，計(jì)算截距的修正系數(shù)μ，即

式中：S為參考面積。

系數(shù)μ用于修正初始下洗角產(chǎn)生的氣動(dòng)力hη0和hξ0，然后對(duì)氣動(dòng)力影響系數(shù)進(jìn)行修正。在飛行仿真的t時(shí)刻，從h1(η，δ)得到零舵偏時(shí)的非線(xiàn)性法向力系數(shù)則修正系數(shù)λ1(t)為

同理，由面元法計(jì)算得到的零迎角、舵偏δ產(chǎn)生的剛體法向氣動(dòng)力為從h1(η，δ)得到零迎角、舵偏δ對(duì)應(yīng)的法向力系數(shù)增量ΔCn，則修正系數(shù)λ2(t)為

將 得 到 的 修 正 系 數(shù)λ1(t)、λ2(t)和μ應(yīng) 用 于，得到飛行仿真中的h2(η，δ)，其表達(dá)式為

4 補(bǔ)丁模塊構(gòu)建

對(duì)比剛性飛機(jī)的氣動(dòng)力，式（8）增加了表征彈性效應(yīng)的分量h2(η，δ)、Δfη和Δfξ。這些分量分別構(gòu)成氣動(dòng)彈性狀態(tài)空間方程、彈性模態(tài)廣義氣動(dòng)力和剛體氣動(dòng)力的增量模塊，用于實(shí)現(xiàn)彈性飛機(jī)的氣動(dòng)耦合，再與測(cè)量信號(hào)的彈性增量模塊結(jié)合，共同構(gòu)成補(bǔ)丁模塊。

4.1 氣動(dòng)彈性狀態(tài)空間方程

根據(jù)式（2）～式（8），忽略發(fā)動(dòng)機(jī)推力對(duì)彈性振動(dòng)的影響，整理彈性振動(dòng)方程，表達(dá)式為

式中：M、C、K分別為廣義質(zhì)量、廣義阻尼、廣義剛度的矩陣形式；Mξδ為彈性模態(tài)與舵面偏轉(zhuǎn)模態(tài)的耦合慣性質(zhì)量。方程（19）左端各系數(shù)矩陣的表達(dá)式為

舵面耦合慣性質(zhì)量Mqδ為6+n行的矩陣，各行分別對(duì)應(yīng)式中的各廣義力。其中的后n行對(duì)應(yīng)矩陣Mξδ，兩者之間的關(guān)系為

式中：Mηδ為剛體模態(tài)與舵面偏轉(zhuǎn)模態(tài)的耦合慣性質(zhì)量。

結(jié)合式（10）、式（11）和式（19），構(gòu)建氣動(dòng)彈性狀態(tài)空間方程，表達(dá)式為

各系數(shù)矩陣的表達(dá)式為

根據(jù)式（22）構(gòu)建氣動(dòng)彈性的狀態(tài)空間方程模塊，輸出彈性模態(tài)廣義坐標(biāo)及其導(dǎo)數(shù)ξ,以及氣動(dòng)力滯后根的狀態(tài)變量xa用于后續(xù)計(jì)算，如圖5所示。

圖5 氣動(dòng)彈性狀態(tài)空間方程模塊Fig.5 Module of aeroelastic state-space equation

4.2 彈性模態(tài)廣義氣動(dòng)力

氣動(dòng)彈性狀態(tài)空間方程需要輸入彈性模態(tài)廣義氣動(dòng)力的非線(xiàn)性部分h2(η，δ)。3.2 節(jié)給出了該部分的計(jì)算和修正方法，表達(dá)式如式（18）所示。因此，輸入狀態(tài)變量η和舵面偏轉(zhuǎn)角δ，即可構(gòu)建的彈性模態(tài)廣義氣動(dòng)力模塊如圖6 所示。

圖6 彈性模態(tài)廣義氣動(dòng)力模塊Fig.6 Module of generalized aerodynamic forces based on elastic modes

4.3 剛體氣動(dòng)力增量

利用氣動(dòng)彈性狀態(tài)空間方程的輸出，根據(jù)式（9），計(jì)算剛體氣動(dòng)力的彈性效應(yīng)增量Δfη，如圖7 所示。根據(jù)式（8），將該增量作用于原有剛性飛機(jī)運(yùn)動(dòng)方程的非線(xiàn)性氣動(dòng)力h1(η，δ)中，實(shí)現(xiàn)剛體氣動(dòng)力的彈性效應(yīng)修正。上述3 部分構(gòu)成廣義氣動(dòng)力的彈性效應(yīng)模塊，可將剛性飛機(jī)的運(yùn)動(dòng)方程拓展到彈性飛機(jī)的剛彈耦合動(dòng)力學(xué)方程。

圖7 剛體氣動(dòng)力的增量模塊Fig.7 Module of increments of rigid aerodynamic forces

4.4 測(cè)量信號(hào)彈性增量

彈性飛機(jī)上的傳感器可同時(shí)感受剛體運(yùn)動(dòng)和彈性振動(dòng)信息，需要在原有剛性飛機(jī)運(yùn)動(dòng)響應(yīng)的基礎(chǔ)上，疊加彈性振動(dòng)引起的測(cè)量信號(hào)增量。

在剛性飛機(jī)飛行仿真中，傳感器測(cè)量飛機(jī)飛行的加速度、角速度和歐拉角。加速度的表達(dá)式為[29］

式中：x、y、z表示測(cè)量點(diǎn)在結(jié)構(gòu)坐標(biāo)系中的坐標(biāo)，一般選取慣導(dǎo)和翼尖等需要關(guān)注的位置。

在此基礎(chǔ)上，疊加彈性振動(dòng)引起的測(cè)量信號(hào)增量，包括加速度、角速度和歐拉角的彈性效應(yīng)增量，其表達(dá)式[10］為

式中：Δax、Δay、Δaz為加速度增量在體軸系3 個(gè)方向的投影；Δp、Δq、Δr為角速度增量在體軸系3個(gè) 方 向 的 分 量；Δφ、Δθ、Δψ為 歐 拉 角 的 增 量；分別為第i階彈性模態(tài)振型的平動(dòng)和轉(zhuǎn)動(dòng)分量。

引 入 彈 性 模 態(tài) 廣 義 坐 標(biāo) 及 其 導(dǎo) 數(shù)ξ、ξ?、ξ?，即可計(jì)算測(cè)量信號(hào)的彈性增量如圖8 所示。將該部分與廣義氣動(dòng)力的彈性效應(yīng)模塊結(jié)合，構(gòu)成表征氣動(dòng)彈性效應(yīng)的“補(bǔ)丁模塊”。將該模塊置入剛性飛機(jī)飛行仿真中，即可實(shí)現(xiàn)基于“補(bǔ)丁模塊”的彈性飛機(jī)飛行仿真，如圖9 所示。

圖8 測(cè)量信號(hào)的彈性增量模塊Fig.8 Module of incremental measurement signals due to elastic vibrations

圖9 基于“補(bǔ)丁模塊”的彈性飛機(jī)飛行仿真Fig.9 Flight simulation of flexible aircraft based on patch module method

5 算例驗(yàn)證

以某大展弦比無(wú)人機(jī)[15］為研究對(duì)象，如圖10所示，驗(yàn)證基于“補(bǔ)丁模塊”的彈性飛機(jī)飛行仿真分析方法的可行性和準(zhǔn)確性。根據(jù)飛機(jī)有限元模型，應(yīng)用模態(tài)截?cái)喾╗8］，選取前20 階彈性模態(tài)進(jìn)行驗(yàn)證。主要彈性模態(tài)頻率如表2 所示，彈性模態(tài)振型如圖11 和圖12 所示。

圖10 某大展弦比無(wú)人機(jī)示意圖Fig.10 Unmanned aerial vehicle with large aspect ratio

圖12 機(jī)翼反對(duì)稱(chēng)一彎模態(tài)示意圖Fig.12 Mode of the first antisymmetric bend

表2 主要彈性模態(tài)頻率Table 2 Frequencies of main elastic modes

圖11 機(jī)翼對(duì)稱(chēng)一彎模態(tài)示意圖Fig.11 Mode of the first symmetric bend

針對(duì)該算例，研究基于“補(bǔ)丁模塊”的彈性飛機(jī)飛行仿真分析，給出對(duì)應(yīng)的配平方法，并對(duì)飛機(jī)的機(jī)動(dòng)響應(yīng)、閉環(huán)反饋、氣動(dòng)力非線(xiàn)性和氣動(dòng)伺服彈性穩(wěn)定性，從不同方面對(duì)該飛行仿真方法驗(yàn)證。

5.1 配平方法

在飛行仿真工作開(kāi)始之前，首先需要確定初始條件，即針對(duì)基于“補(bǔ)丁模塊”的彈性飛機(jī)飛行仿真分析方法，研究飛機(jī)的配平問(wèn)題。由于“補(bǔ)丁模塊”在剛性飛機(jī)飛行仿真的基礎(chǔ)上添加，該補(bǔ)丁方法保留了原有剛性飛機(jī)的配平結(jié)果。但該結(jié)果在彈性飛機(jī)在飛行仿真初期不能保持穩(wěn)定，需重新給出針對(duì)彈性飛機(jī)的配平方法。

根據(jù)飛行仿真特點(diǎn)，給出2 種可行的配平方法，分別是剛彈耦合和閉環(huán)反饋輔助配平。剛彈耦合配平是直接應(yīng)用剛彈耦合動(dòng)力學(xué)方程，同時(shí)考慮剛體和彈性自由度進(jìn)行配平。閉環(huán)反饋輔助配平是構(gòu)建如圖13 所示的迎角反饋回路輔助配平，a=0.036 為舵機(jī)參數(shù)；δ為升降舵偏轉(zhuǎn)角；δc為舵機(jī)偏轉(zhuǎn)指令信號(hào)；K為反饋系統(tǒng)的控制律。

圖13 迎角閉環(huán)反饋系統(tǒng)Fig.13 Close systems of angle of attack feedback

在高度6 km、飛行速度32.8 m/s 的飛行狀態(tài)下，進(jìn)行彈性飛機(jī)的飛行仿真，得到2 種配平方法的響應(yīng)曲線(xiàn)，并與原有剛性飛機(jī)配平結(jié)果在彈性飛機(jī)中的響應(yīng)曲線(xiàn)進(jìn)行對(duì)比如圖14 所示。

從圖14 中可以看出，由于彈性自由度的存在，以原有剛性飛機(jī)配平結(jié)果為初始條件時(shí)，飛機(jī)響應(yīng)存在明顯的振蕩收斂現(xiàn)象。但z方向（豎直方向）加速度和俯仰角速度無(wú)法收斂到0，導(dǎo)致俯仰角發(fā)散，彈性飛機(jī)無(wú)法保持穩(wěn)定，故引入針對(duì)彈性飛機(jī)的配平方法是必要的。

圖14 不同配平方法在彈性飛機(jī)中的響應(yīng)曲線(xiàn)Fig.14 Response of flexible aircraft with different trim methods

在2 種配平方法中，剛彈耦合配平能夠使得飛機(jī)處于穩(wěn)定狀態(tài)，不需要經(jīng)過(guò)收斂過(guò)程。而閉環(huán)反饋輔助配平方法，通過(guò)迎角反饋使得加速度和俯仰角速度在1 s 內(nèi)振蕩收斂至0，使得俯仰角趨于穩(wěn)定，滿(mǎn)足配平要求。故這2 種配平方法都是行之有效的，可根據(jù)實(shí)際情況選擇。

在工程實(shí)際中，可選擇閉環(huán)反饋輔助配平的方法，減少“補(bǔ)丁模塊”對(duì)原有剛性飛機(jī)飛行仿真的修改；而對(duì)于簡(jiǎn)單模型，選取剛彈耦合配平方法相對(duì)簡(jiǎn)單，易于實(shí)現(xiàn)。出于簡(jiǎn)化考慮，本文選取剛彈耦合配平方法進(jìn)行后續(xù)研究工作。

5.2 俯仰機(jī)動(dòng)響應(yīng)

在剛彈耦合配平的基礎(chǔ)上，研究基于“補(bǔ)丁模塊”彈性飛機(jī)飛行仿真對(duì)機(jī)動(dòng)響應(yīng)的影響，分為俯仰機(jī)動(dòng)響應(yīng)、偏航和滾轉(zhuǎn)機(jī)動(dòng)響應(yīng)這2 部分。在俯仰機(jī)動(dòng)響應(yīng)中，以如圖15 所示的舵面偏轉(zhuǎn)信號(hào)作為升降舵激勵(lì)，對(duì)比剛性飛機(jī)飛行仿真、基于剛彈耦合的彈性飛機(jī)飛行仿真和基于“補(bǔ)丁模塊”的彈性飛機(jī)飛行仿真，驗(yàn)證該補(bǔ)丁方法的可行性和準(zhǔn)確性。俯仰機(jī)動(dòng)響應(yīng)曲線(xiàn)如圖16 所示。

圖15 舵面偏轉(zhuǎn)信號(hào)示意圖Fig.15 Deflection of control surfaces

由圖16 的響應(yīng)曲線(xiàn)可知，在縱向飛行仿真中，氣動(dòng)彈性效應(yīng)使得飛機(jī)機(jī)動(dòng)響應(yīng)的超調(diào)量和峰值時(shí)間發(fā)生變化，導(dǎo)致彈性飛機(jī)的俯仰機(jī)動(dòng)響應(yīng)明顯區(qū)別于剛性飛機(jī)。

圖16 俯仰機(jī)動(dòng)響應(yīng)曲線(xiàn)Fig.16 Response of pitch maneuver

對(duì)比2 種彈性飛機(jī)飛行仿真方法的結(jié)果，發(fā)現(xiàn)基于“補(bǔ)丁模塊”彈性飛機(jī)飛行仿真與基于剛彈耦合飛行仿真具有一致的仿真結(jié)果。該結(jié)果表明“補(bǔ)丁模塊”能夠有效地引入氣動(dòng)彈性效應(yīng)，將飛行仿真由剛性飛機(jī)推廣到彈性飛機(jī)，驗(yàn)證了該方法在6 自由度全量方程剛性飛機(jī)飛行仿真俯仰通道上應(yīng)用的可行性和準(zhǔn)確性。

5.3 偏航和滾轉(zhuǎn)機(jī)動(dòng)響應(yīng)

針對(duì)彈性飛機(jī)在6 自由度全量方程中的偏航和滾轉(zhuǎn)通道，研究“補(bǔ)丁模塊”引入氣動(dòng)彈性效應(yīng)對(duì)偏航和滾轉(zhuǎn)機(jī)動(dòng)響應(yīng)的影響效果。分別在方向舵和副翼上，作用如圖15 所示的舵偏激勵(lì)，研究表征橫航向動(dòng)力學(xué)特性的機(jī)動(dòng)響應(yīng)。其中，偏航機(jī)動(dòng)響應(yīng)曲線(xiàn)如圖17 所示；滾轉(zhuǎn)機(jī)動(dòng)響應(yīng)曲線(xiàn)如圖18 所示。比較偏航和滾轉(zhuǎn)的機(jī)動(dòng)響應(yīng)曲線(xiàn)可知，氣動(dòng)彈性效應(yīng)對(duì)飛機(jī)橫航向機(jī)動(dòng)響應(yīng)具有顯著影響，在滾轉(zhuǎn)機(jī)動(dòng)的翼尖加速度響應(yīng)中該現(xiàn)象尤為明顯。

圖17 偏航機(jī)動(dòng)響應(yīng)曲線(xiàn)Fig.17 Response of yaw maneuver

圖18 滾轉(zhuǎn)機(jī)動(dòng)響應(yīng)曲線(xiàn)Fig.18 Response of roll maneuver

基于“補(bǔ)丁模塊”彈性飛機(jī)的飛行仿真與基于剛彈耦合的飛行仿真基本重合，能夠很好地反映飛機(jī)的氣動(dòng)彈性效應(yīng)，驗(yàn)證了“補(bǔ)丁模塊”引入彈性效應(yīng)在偏航和滾轉(zhuǎn)通道上應(yīng)用的可行性和準(zhǔn)確性。

5.4 閉環(huán)反饋

與基于剛彈耦合彈性飛機(jī)的飛行仿真相比，“補(bǔ)丁模塊”的優(yōu)勢(shì)在于，對(duì)彈性飛機(jī)進(jìn)行控制系統(tǒng)設(shè)計(jì)、穩(wěn)定性分析和驗(yàn)證工作時(shí)，僅需在原有剛性飛機(jī)飛行仿真模型中進(jìn)行處理，減小了彈性自由度對(duì)其他設(shè)計(jì)和驗(yàn)證程序模塊的影響。

因此，以具有簡(jiǎn)單控制系統(tǒng)的閉環(huán)反饋為例，驗(yàn)證補(bǔ)丁模塊對(duì)閉環(huán)系統(tǒng)響應(yīng)的影響。選取同樣的無(wú)人機(jī)模型作為研究對(duì)象，設(shè)計(jì)簡(jiǎn)單的俯仰角速度反饋如圖19 所示。

圖19 俯仰角速度閉環(huán)反饋系統(tǒng)Fig.19 Close systems of pitch rate feedback

選取如圖20 所示的階躍信號(hào)作為升降舵機(jī)偏轉(zhuǎn)指令，研究閉環(huán)系統(tǒng)在舵偏擾動(dòng)下的俯仰機(jī)動(dòng)特性如圖21 所示。比較不同仿真方法的響應(yīng)曲線(xiàn)可知，基于“補(bǔ)丁模塊”的彈性飛機(jī)飛行仿真與基于剛彈耦合的彈性飛機(jī)飛行仿真具有相同結(jié)果，能夠精確地反映閉環(huán)系統(tǒng)動(dòng)響應(yīng)特性，特別是翼尖加速度的振蕩現(xiàn)象，驗(yàn)證了該方法在閉環(huán)系統(tǒng)響應(yīng)仿真中應(yīng)用的可行性和準(zhǔn)確性，為彈性飛機(jī)控制系統(tǒng)設(shè)計(jì)和驗(yàn)證工作提供了新的思路。

圖20 舵機(jī)偏轉(zhuǎn)指令信號(hào)Fig.20 Command signal of control surface deflections

圖21 閉環(huán)系統(tǒng)俯仰機(jī)動(dòng)響應(yīng)曲線(xiàn)Fig.21 Response of pitch maneuver in close system

5.5 非線(xiàn)性氣動(dòng)力驗(yàn)證

在6 自由度全量方程的剛性飛機(jī)飛行仿真中，氣動(dòng)力一般采用非線(xiàn)性形式。因此，需要驗(yàn)證基于“補(bǔ)丁模塊”彈性飛機(jī)飛行仿真方法對(duì)氣動(dòng)力非線(xiàn)性情況的適用性。針對(duì)該大展弦比無(wú)人機(jī)算例，由CFD 技術(shù)分析得到法向力系數(shù)與迎角的關(guān)系如圖22 所示。

圖22 法向力系數(shù)隨迎角的變化曲線(xiàn)Fig.22 Curve of aerodynamic coefficient of normal force with angle of attack

根據(jù)該非線(xiàn)性法向力系數(shù)曲線(xiàn)，應(yīng)用3.2 節(jié)對(duì)彈性模態(tài)廣義氣動(dòng)力的非線(xiàn)性部分進(jìn)行計(jì)算與修正，使得彈性模態(tài)廣義氣動(dòng)力中與剛體運(yùn)動(dòng)相關(guān)部分，呈現(xiàn)出同樣的非線(xiàn)性變化趨勢(shì)。

以5.2 節(jié)中的俯仰機(jī)動(dòng)響應(yīng)為基礎(chǔ)，將振幅為7°、圖15 所示的激勵(lì)信號(hào)作用到升降舵上，使得響應(yīng)中的法向力系數(shù)達(dá)到非線(xiàn)性段。以第一階彈性模態(tài)廣義氣動(dòng)力系數(shù)為例，俯仰機(jī)動(dòng)響應(yīng)過(guò)程中，彈性模態(tài)廣義氣動(dòng)力系數(shù)隨迎角的非線(xiàn)性變化情況如圖23 所示。

圖23 第一階彈性模態(tài)廣義氣動(dòng)力系數(shù)隨迎角變化曲線(xiàn)Fig.23 Curve of aerodynamic coefficient of the first elastic mode with angle of attack

應(yīng)用線(xiàn)性和非線(xiàn)性氣動(dòng)力的剛性飛機(jī)飛行仿真、線(xiàn)性氣動(dòng)力的彈性飛機(jī)飛行仿真，以及非線(xiàn)性氣動(dòng)力的基于剛彈耦合和“補(bǔ)丁模塊”的彈性飛機(jī)飛行仿真，進(jìn)行升降舵偏轉(zhuǎn)激勵(lì)下的俯仰機(jī)動(dòng)響應(yīng)如圖24 所示。根據(jù)仿真結(jié)果可知，非線(xiàn)性氣動(dòng)力改變了飛機(jī)動(dòng)響應(yīng)的變化趨勢(shì)，慣導(dǎo)和翼尖處加速度、俯仰角速度響應(yīng)的超調(diào)量明顯降低，翼尖加速度存在明顯的振蕩現(xiàn)象；俯仰角的響應(yīng)由于非線(xiàn)性氣動(dòng)力出現(xiàn)了不同的變化趨勢(shì)。這些現(xiàn)象表明非線(xiàn)性氣動(dòng)力對(duì)剛性和彈性飛機(jī)的機(jī)動(dòng)響應(yīng)具有顯著影響。

圖24 氣動(dòng)力非線(xiàn)性時(shí)的俯仰機(jī)動(dòng)響應(yīng)曲線(xiàn)Fig.24 Response of pitch maneuver with nonlinear aerodynamics

對(duì)比基于剛彈耦合和“補(bǔ)丁模塊”的彈性飛機(jī)飛行仿真響應(yīng)曲線(xiàn)可知，在非線(xiàn)性氣動(dòng)力的條件下，“補(bǔ)丁模塊”能夠有效地引入氣動(dòng)彈性效應(yīng)和非線(xiàn)性效果，與求解剛彈耦合方程具有相同的效果。因此，基于“補(bǔ)丁模塊”的彈性飛機(jī)飛行仿真能夠通過(guò)非線(xiàn)性氣動(dòng)力修正，有效地模擬彈性飛機(jī)的動(dòng)響應(yīng)，簡(jiǎn)化了彈性飛機(jī)飛行仿真的建模和分析流程。

5.6 氣動(dòng)伺服彈性穩(wěn)定性

在彈性飛機(jī)閉環(huán)系統(tǒng)中，氣動(dòng)彈性和飛行控制系統(tǒng)存在一定的耦合現(xiàn)象，形成了氣動(dòng)伺服彈性力學(xué)[30-33］。氣動(dòng)伺服彈性穩(wěn)定性便是其中的一個(gè)重要領(lǐng)域。

針對(duì)相同的算例模型，構(gòu)建滾轉(zhuǎn)角速度閉環(huán)反饋系統(tǒng)如圖25 所示，δ為副翼偏轉(zhuǎn)角。調(diào)節(jié)控制系統(tǒng)的增益系數(shù)至K=-3，使得彈性飛機(jī)的閉環(huán)系統(tǒng)處于臨界穩(wěn)定狀態(tài)。此時(shí)，以如圖20 所示的階躍信號(hào)作為舵機(jī)輸入，采用不同的方法進(jìn)行飛行仿真分析，得到副翼偏轉(zhuǎn)擾動(dòng)下的動(dòng)響應(yīng)曲線(xiàn)如圖26 所示。

圖25 滾轉(zhuǎn)角速度閉環(huán)反饋系統(tǒng)Fig.25 Close systems of rolling rate feedback

圖26 氣動(dòng)伺服彈性臨界穩(wěn)定條件下的動(dòng)響應(yīng)曲線(xiàn)Fig.26 Dynamic response under condition of aeroservoelasticity critical stability

由響應(yīng)曲線(xiàn)結(jié)果可知，剛性飛機(jī)在副翼偏轉(zhuǎn)擾動(dòng)下處于穩(wěn)定狀態(tài)，而彈性飛機(jī)由于引入氣動(dòng)彈性效應(yīng)，氣動(dòng)彈性和控制系統(tǒng)存在一定的耦合，使得飛機(jī)處于臨界穩(wěn)定狀態(tài)。當(dāng)飛機(jī)受到副翼偏轉(zhuǎn)擾動(dòng)時(shí)，z方向（豎直方向）加速度和滾轉(zhuǎn)角速度存在明顯的振蕩現(xiàn)象。該振蕩現(xiàn)象的頻率6.40 Hz，接近表2 中的機(jī)翼反對(duì)稱(chēng)一彎頻率7.06 Hz。表明該臨界穩(wěn)定現(xiàn)象由機(jī)翼反對(duì)稱(chēng)一彎及其他模態(tài)耦合產(chǎn)生。

對(duì)比不同方法的仿真結(jié)果可知，基于“補(bǔ)丁模塊”的彈性飛機(jī)飛行仿真能很好地反映氣動(dòng)彈性效應(yīng)導(dǎo)致的臨界穩(wěn)定現(xiàn)象，與基于剛彈耦合的彈性飛機(jī)飛行仿真具有一致的響應(yīng)結(jié)果。

因此，“補(bǔ)丁模塊”能良好地反映彈性飛機(jī)的氣動(dòng)彈性效應(yīng)，驗(yàn)證了該飛行仿真方法在氣動(dòng)伺服彈性穩(wěn)定性分析中應(yīng)用的可行性和準(zhǔn)確性。同時(shí)，由于該方法能充分利用原有剛性飛機(jī)的飛行仿真模型，以及在此基礎(chǔ)上設(shè)計(jì)的控制系統(tǒng)，豐富了氣動(dòng)伺服彈性穩(wěn)定性的分析和驗(yàn)證方法。

6 結(jié) 論

1）建立了一種適用于彈性飛機(jī)飛行仿真的補(bǔ)丁方法。該方法通過(guò)廣義氣動(dòng)力的拆分、測(cè)量信號(hào)彈性增量的疊加構(gòu)建“補(bǔ)丁模塊”，將該模塊置入6 自由度全量方程的剛性飛機(jī)飛行仿真中，實(shí)現(xiàn)基于“補(bǔ)丁模塊”的彈性飛機(jī)飛行仿真。

2）以某大展弦比無(wú)人機(jī)為例，給出了基于“補(bǔ)丁模塊”彈性飛機(jī)飛行仿真方法的配平方法，驗(yàn)證了該飛行仿真方法針對(duì)飛機(jī)機(jī)動(dòng)響應(yīng)的可行性和準(zhǔn)確性，顯示了該方法能夠考慮氣動(dòng)力的非線(xiàn)性，可應(yīng)用于氣動(dòng)伺服彈性穩(wěn)定性分析。

3）基于“補(bǔ)丁模塊”彈性飛機(jī)飛行仿真能夠在引入氣動(dòng)彈性效應(yīng)的同時(shí)，充分利用原有剛性飛機(jī)飛行仿真的整體框架，便于將動(dòng)力學(xué)特性研究、控制律設(shè)計(jì)與驗(yàn)證等工作直接推廣到彈性飛機(jī)上，簡(jiǎn)化了彈性飛機(jī)的建模過(guò)程，為彈性飛機(jī)的飛行仿真提供了新思路。