楊建忠，盧 勇，楊士斌

（中國民航大學(xué)適航學(xué)院，天津 300300）

飛機(jī)舒適性、機(jī)動性等性能的提高及生產(chǎn)和使用成本的降低使得飛機(jī)結(jié)構(gòu)的相對重量更輕，表現(xiàn)出越來越高的柔性，從而降低了彈性模態(tài)的頻率[1-3]。 因此，飛機(jī)剛體運動與氣動彈性響應(yīng)之間會出現(xiàn)耦合，導(dǎo)致柔性飛機(jī)受氣動力的影響與剛體飛機(jī)不同，如果按照傳統(tǒng)剛體假設(shè)去設(shè)計飛行控制，可能會與真實飛行有較大差異，甚至產(chǎn)生不利影響危及飛機(jī)的正常飛行，因此必須考慮剛體和彈性體耦合的飛行動力學(xué)模型。

在建立飛機(jī)剛彈耦合模型時，大多基于機(jī)翼進(jìn)行小變形假設(shè)：Waszak 等[4]采用平均軸系法將剛體飛行動力學(xué)與線性結(jié)構(gòu)動力學(xué)耦合，將彈性和剛性自由度解耦，使兩者之間的耦合僅通過廣義力產(chǎn)生；Meirovitch 等[5-6]采用模態(tài)分支法將飛機(jī)視為由多個彈性體構(gòu)成，剛體運動由飛機(jī)重心所在體軸系描述，彈性運動則由各彈性體建立與其固連的各參考系之間的運動關(guān)系表示。平均軸系法與模態(tài)分支法的結(jié)構(gòu)動力學(xué)都是線性的，因此，當(dāng)機(jī)翼等彈性結(jié)構(gòu)存在較為明顯的變形時，用這些方法分析飛機(jī)的運動是不適用的。

當(dāng)機(jī)翼作大變形運動時，可以將機(jī)翼視為非線性梁，用非線性結(jié)構(gòu)動力學(xué)來描述運動。絕大多數(shù)非線性梁模型分為3 種：本征梁模型[7-9]、基于位移的非線性梁模型[10]、基于應(yīng)變的非線性梁模型[11-13]。 由于這3種模型選用的自變量不同導(dǎo)致各有不同的優(yōu)缺點。 本征梁模型自變量個數(shù)過多，基于位移的非線性梁模型建模過程較為繁瑣且收斂較慢，而基于應(yīng)變的非線性梁模型雖然未知數(shù)求解量較多，但具有較低的階數(shù)。

目前，非線性梁模型主要運用于高空長航時無人機(jī)建模，而用該模型建立民用運輸類飛機(jī)剛彈耦合模型的研究較少。本研究采用基于應(yīng)變的非線性梁模型建立民用運輸類飛機(jī)的剛彈耦合模型；通過研究從剛體到剛彈耦合模型的轉(zhuǎn)換獲得與飛行動力學(xué)有關(guān)的等效參數(shù)；基于大變形的剛彈耦合模型與剛體模型響應(yīng)的差異采用非線性動態(tài)過程的時域仿真進(jìn)行分析。該方法可應(yīng)用于大柔性飛機(jī)的飛行控制設(shè)計。

1 大變形機(jī)翼飛機(jī)運動方程

1.1 基于應(yīng)變的非線性梁模型

基于應(yīng)變的非線性梁模型考慮了慣性坐標(biāo)系（G）、機(jī)體坐標(biāo)系（B）、局部坐標(biāo)系（w）和輔助坐標(biāo)系（b），如圖1 所示。 慣性坐標(biāo)系用于定義飛機(jī)相對于地球的位置；機(jī)體坐標(biāo)系用于定義飛機(jī)剛體的自由度；局部坐標(biāo)系用于定義梁每個節(jié)點的位置和方向；輔助坐標(biāo)系用于定義連接到梁上的剛性單元的位置和方向，如連接到機(jī)翼上的發(fā)動機(jī)。 需要注意各坐標(biāo)系軸的定義與研究傳統(tǒng)剛體飛機(jī)的不同，如機(jī)體坐標(biāo)系中Bx垂直于機(jī)身對稱軸線指向右，By沿機(jī)身對稱軸線指向前，Bz根據(jù)右手螺旋定則定義。

圖1 坐標(biāo)系統(tǒng)Fig.1 Coordinate system

根據(jù)該坐標(biāo)系統(tǒng)，飛機(jī)的運動可以表示為

式中：b 表示慣性坐標(biāo)系下的廣義位移；PB、θB表示機(jī)體坐標(biāo)系相對于慣性坐標(biāo)系的位置和方向；表示b對時間求導(dǎo)，可得機(jī)體坐標(biāo)系的速度β，再對時間求導(dǎo)，可得加速度；VB=[U V W]T表示Bx、By、Bz軸方向上的速度；ωB= [P Q R]T表示繞Bx、By、Bz軸轉(zhuǎn)動的角速度。

基于應(yīng)變的非線性梁理論中結(jié)構(gòu)變形用應(yīng)變表示為

式中：εx是梁單元中沿著梁展向的變形；κx、κy、κz分別表示梁單元相對于局部坐標(biāo)系中wx、wy、wz軸的扭轉(zhuǎn)。

在圖1 坐標(biāo)系統(tǒng)中Pw表示局部坐標(biāo)系相對于機(jī)體坐標(biāo)系的位置，局部坐標(biāo)系相對于慣性坐標(biāo)系的位移Pa表示為

由式（5）可以得到梁單元中s 節(jié)點t 時刻位置和的方向矢量表示

1.2 剛彈耦合模型的運動方程

根據(jù)虛功原理，飛機(jī)內(nèi)力產(chǎn)生的虛功和外力產(chǎn)生的虛功之和為0。 通過計算內(nèi)、外虛功，可以推導(dǎo)出剛彈耦合模型的運動方程表示為

式中：MFF、MFB、MBF、MBB表示廣義質(zhì)量矩陣的子矩陣；CFF、CFB、CBF、CBB表示廣義阻尼矩陣的子矩陣；KFF表示結(jié)構(gòu)剛度矩陣；RF和RB表示廣義力，通過計算施加到每個結(jié)構(gòu)節(jié)點的氣動力，可以獲得重力和發(fā)動機(jī)推力。由于廣義質(zhì)量矩陣不是對角陣，可以看到飛機(jī)剛體運動β 與結(jié)構(gòu)彈性運動ε 存在慣性耦合。

方程（9）中等號左邊的子矩陣為

式中：M 表示彈性結(jié)構(gòu)質(zhì)量矩陣，其僅取決于彈性結(jié)構(gòu)的慣性和質(zhì)量而不取決于應(yīng)變；MRB表示廣義剛體機(jī)身質(zhì)量矩陣，其取決于剛體機(jī)身慣性和質(zhì)量；C 表示結(jié)構(gòu)阻尼矩陣，此處，C=ζKFF，ζ 為阻尼比；CRB是廣義剛體機(jī)身阻尼矩陣，取決于剛體機(jī)身質(zhì)量和角速度；Hhb是以ωB的叉乘矩陣為子矩陣的對角矩陣與Jhb相乘的矩陣。

方程（9）中的廣義力表示為

式中：Fpt、Mpt表示點力和點力矩；BF、BM、N 表示影響系數(shù)，由數(shù)值積分得到；Fdist、Mdist表示分布力和分布力矩；G 表示重力；表示廣義外力， 由作用于剛體機(jī)身的外力和外力矩組成。

描述柔性飛機(jī)的運動學(xué)方程為

式中：θ、φ、ψ 分別表示飛機(jī)的俯仰角、滾轉(zhuǎn)角、偏航角；x、y、H 表示慣性坐標(biāo)系中飛機(jī)3 個方向的位移。

2 飛機(jī)模型參數(shù)

2.1 飛機(jī)幾何參數(shù)

原始飛機(jī)幾何數(shù)據(jù)來源于Silva[14]的研究。Silva 提出的柔性飛行器模型有3 種，飛機(jī)幾何數(shù)據(jù)選擇的是第1 種模型Conf1。由于在剛彈耦合模型和剛體模型之間的比較側(cè)重于飛機(jī)的剛體運動，Conf1 模型可以方便地對具有適度結(jié)構(gòu)柔性的飛機(jī)進(jìn)行建模和分析。如果使用具有更高結(jié)構(gòu)柔性的另外兩種飛機(jī)模型，其柔性可能會顯著影響飛機(jī)的操縱性和穩(wěn)定性，導(dǎo)致飛機(jī)剛體運動受到嚴(yán)重影響。Conf1 模型中飛機(jī)長33 m，機(jī)身截面半徑1.5 m，彈性部件幾何參數(shù)如表1 所示。飛機(jī)視圖如圖2 所示。

表1 彈性部件幾何參數(shù)Tab.1 Geometric parameters of elastic components

圖2 飛機(jī)視圖Fig.2 View of aircraft

2.2 飛機(jī)氣動參數(shù)分布

為了便于計算，選擇準(zhǔn)定常氣動力模型計算飛機(jī)受到的氣動力載荷。 機(jī)翼剖面如圖3 所示，圖中b 為機(jī)翼半弦長，d 為機(jī)翼的氣動軸系統(tǒng)到翼弦中點的距離。為速度V 的分量，L 表示升力，D 表示阻力，M表示力矩，aOz表示零升力氣動坐標(biāo)系的z 軸，aIz表示氣流坐標(biāo)系的z 軸。

圖3 機(jī)翼剖面圖Fig.3 Profile of airfoil

根據(jù)亞音速不可壓二維流的準(zhǔn)定常空氣動力學(xué)理論[15]，可得到機(jī)翼單元上氣動系數(shù)的氣動力和力矩，表示如下

式中：M、Mx分別表示1/4 弦長處計算的力矩、1/2 弦長處計算的力矩，由于機(jī)翼彈性軸位于機(jī)翼1/2 弦長處， 而氣動力和力矩是在1/4 弦長處的氣動中心計算的，所以1/4 弦長處的力和力矩必須移到1/2 弦長處；ρ 是大氣密度；α 是局部迎角；t 時刻的局部迎角αt=arctan是零升力迎角；CLα、CD0、Cm0是穩(wěn)定性導(dǎo)數(shù)，分別表示升力線斜率、零升力矩系數(shù)和阻力系數(shù)；δi是飛機(jī)上舵面的偏轉(zhuǎn)量（i=a、e、r 分別表示副翼、升降舵、方向舵）；CLδi、Cmδi是控制導(dǎo)數(shù)。 需要注意的是，穩(wěn)定性導(dǎo)數(shù)和控制導(dǎo)數(shù)表示的含義與飛行動力學(xué)中表示的含義不同。 例如，在飛行動力學(xué)中，穩(wěn)定性導(dǎo)數(shù)CLα表示迎角變化對飛機(jī)升力的影響，控制導(dǎo)數(shù)Cmδe表示升降舵偏轉(zhuǎn)對飛機(jī)重心位置附近俯仰力矩的影響；在此處，穩(wěn)定性導(dǎo)數(shù)CLα表示迎角變化對翼型截面所在當(dāng)?shù)貦C(jī)翼單元升力的影響，控制導(dǎo)數(shù)Cmδe表示升降舵偏轉(zhuǎn)對翼型截面1/4 弦長處翼型截面當(dāng)?shù)厮跈C(jī)翼單元俯仰力矩的影響。

氣動力系數(shù)通過Tornado 代碼計算，該代碼是Tomas Melin 采用渦格法開發(fā)的亞聲速機(jī)翼氣動特性分析程序。在Silva 的模型中，氣動力系數(shù)是飛機(jī)上產(chǎn)生的力和力矩的總系數(shù)。在采用基于應(yīng)變的非線性梁模型建立的剛彈耦合模型中，氣動力系數(shù)是沿機(jī)翼展長分布的。 要調(diào)整沿機(jī)翼展長方向上氣動力系數(shù)的分布，氣動力系數(shù)的分布結(jié)果必須與用剛體方法建模產(chǎn)生的力和力矩效果相同。以飛機(jī)升力系數(shù)CLα為例，計算公式為

式中：CLαw、CLαh分別表示機(jī)翼升力線斜率和平尾升力線斜率，Sw、Sh分別表示機(jī)翼面積和平尾面積；CLαi表示機(jī)翼第i 個單元的升力線斜率；Ai表示第i 個單元的面積；?表示下洗。表2～表5 給出了飛機(jī)在速度為224.6 m/s、高度為10 000 m 平飛條件下等效的穩(wěn)定性導(dǎo)數(shù)和控制導(dǎo)數(shù)分布值（其中CLαv表示垂尾升力線斜率）。

表2 機(jī)翼等效導(dǎo)數(shù)分布Tab.2 Equivalent derivative distribution of wing

表3 平尾等效導(dǎo)數(shù)分布Tab.3 Equivalent derivatives distribution of horizontal tail

表4 垂尾等效導(dǎo)數(shù)分布Tab.4 Equivalent derivative distribution of vertical tail

表5 彈性部件其他氣動參數(shù)Tab.5 Aerodynamic parameters of other elastic components

3 仿真結(jié)果分析

3.1 滾轉(zhuǎn)機(jī)動

在1 s 時給定-5°～5°周期為2 s 的副翼偏轉(zhuǎn)輸入，如圖4（a）所示，滾轉(zhuǎn)角速度隨時間變化的曲線如圖4（b）所示。 由圖4 可看出，副翼偏轉(zhuǎn)為正時，滾轉(zhuǎn)角速度為負(fù)值。這是因為此時右側(cè)機(jī)翼產(chǎn)生升力增加，左側(cè)機(jī)翼產(chǎn)生升力減小，飛機(jī)向左側(cè)滾轉(zhuǎn)。通過剛彈耦合和剛體兩種模型下滾轉(zhuǎn)角速度的峰值差，計算出副翼偏轉(zhuǎn)時剛彈耦合模型副翼效率比剛體模型低18.77%。

圖4 副翼偏轉(zhuǎn)響應(yīng)Fig.4 Eeflection response of aileron

圖5 給出了兩側(cè)機(jī)翼的翼根彎矩、 翼根扭矩、機(jī)翼撓度隨時間的變化。 由圖5 可看出，翼根彎矩比翼根扭矩大1 個數(shù)量級，對飛機(jī)滾轉(zhuǎn)操縱影響最大的是機(jī)翼的彎曲。 副翼偏轉(zhuǎn)時，由于剛彈耦合模型中機(jī)翼具有柔性，副翼偏轉(zhuǎn)誘導(dǎo)兩側(cè)機(jī)翼扭轉(zhuǎn)減小了當(dāng)?shù)氐刃в?，另一方面兩?cè)不對稱的彎曲改變了氣動載荷分布，從而造成飛機(jī)的滾轉(zhuǎn)力矩減小，降低了副翼舵效率。

圖5 機(jī)翼彈性響應(yīng)（飛機(jī)滾轉(zhuǎn)）Fig.5 Elastic response of wing（Aircraft roll）

3.2 垂直陣風(fēng)響應(yīng)

根據(jù)CCAR 25.341 條款中規(guī)定的離散突風(fēng)設(shè)計準(zhǔn)則得到

式中：Uw為陣風(fēng)形狀；Uds為設(shè)計陣風(fēng)速度；s 為飛機(jī)進(jìn)入陣風(fēng)區(qū)的距離；Hw為陣風(fēng)梯度，即達(dá)到其峰值速度時與飛機(jī)飛行航跡的平行距離，范圍為9.1～106.7 m；Uref為參考陣風(fēng)速度；Fg為飛行剖面緩和系數(shù)。 由于所選的飛機(jī)是在10 000 m 高度上建模，而10 000 m 高度被認(rèn)為是飛機(jī)的運行高度?；谶@一事實，參數(shù)Fg估計為1.0，計算得到Uds=10.61 m/s 的“1-cos”垂直離散陣風(fēng)。圖6 給出了俯仰角、俯仰角速度隨時間變化的曲線。

圖6 垂直離散陣風(fēng)響應(yīng)Fig.6 Vertical discrete gust response

由圖6 可看出，在受到垂直陣風(fēng)擾動時剛彈耦合模型和剛體模型縱向響應(yīng)基本一致，表明機(jī)翼的剛度對飛機(jī)縱向影響很小。

圖7 給出了右翼根彎矩、右翼根扭矩、左右兩側(cè)機(jī)翼的撓度隨時間的變化。

圖7 機(jī)翼彈性響應(yīng)（垂直陣風(fēng)）Fig.7 Elastic response of wing(Vertical gust)

由圖7 可看出，翼根彎矩比翼根扭矩大1 個數(shù)量級，垂直陣風(fēng)對飛機(jī)影響最大的是機(jī)翼的彎曲。遭遇垂直陣風(fēng)時，剛彈耦合模型預(yù)測的翼根彎矩峰值比剛體模型低8.35%。從能量角度考慮，陣風(fēng)擾動的能量除了轉(zhuǎn)化為剛體俯仰沉浮的動能外，還轉(zhuǎn)化為機(jī)翼和尾翼的彈性能。因而，按照傳統(tǒng)剛體假設(shè)預(yù)測的陣風(fēng)載荷會過于保守，不利于飛機(jī)的減重。

4 結(jié)語

采用基于應(yīng)變的非線性梁模型對柔性飛機(jī)建立了剛彈耦合的飛行動力學(xué)模型，通過時域仿真分析剛彈耦合模型和剛體模型的響應(yīng)，得到了如下結(jié)論：①基于應(yīng)變的剛彈耦合建模方法能夠同時表達(dá)飛機(jī)彈性元件的變形和飛機(jī)的剛體運動，并模擬彈性變形與飛行動力學(xué)參數(shù)的相互影響；②在滾轉(zhuǎn)機(jī)動仿真中，剛彈耦合模型的飛機(jī)能夠模擬副翼舵效率的降低；③在垂直陣風(fēng)響應(yīng)仿真中，剛彈耦合模型與傳統(tǒng)剛體建模相比具有更低的翼根彎矩。 基于應(yīng)變的柔性飛機(jī)剛彈耦合建模方法可用于大柔性飛機(jī)的飛行控制律和陣風(fēng)載荷減緩設(shè)計。