楊建忠，盧 勇，楊士斌

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

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

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

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

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

1 大變形機(jī)翼飛機(jī)運(yùn)動(dòng)方程

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)系用于定義梁每個(gè)節(jié)點(diǎn)的位置和方向；輔助坐標(biāo)系用于定義連接到梁上的剛性單元的位置和方向，如連接到機(jī)翼上的發(fā)動(dòng)機(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ī)的運(yùn)動(dòng)可以表示為

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

基于應(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é)點(diǎn)t 時(shí)刻位置和的方向矢量表示

1.2 剛彈耦合模型的運(yùn)動(dòng)方程

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

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

方程（9）中等號(hào)左邊的子矩陣為

式中：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表示點(diǎn)力和點(diǎn)力矩；BF、BM、N 表示影響系數(shù)，由數(shù)值積分得到；Fdist、Mdist表示分布力和分布力矩；G 表示重力；表示廣義外力， 由作用于剛體機(jī)身的外力和外力矩組成。

描述柔性飛機(jī)的運(yùn)動(dòng)學(xué)方程為

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

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

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

原始飛機(jī)幾何數(shù)據(jù)來源于Silva[14]的研究。Silva 提出的柔性飛行器模型有3 種，飛機(jī)幾何數(shù)據(jù)選擇的是第1 種模型Conf1。由于在剛彈耦合模型和剛體模型之間的比較側(cè)重于飛機(jī)的剛體運(yùn)動(dòng)，Conf1 模型可以方便地對具有適度結(jié)構(gòu)柔性的飛機(jī)進(jìn)行建模和分析。如果使用具有更高結(jié)構(gòu)柔性的另外兩種飛機(jī)模型，其柔性可能會(huì)顯著影響飛機(jī)的操縱性和穩(wěn)定性，導(dǎo)致飛機(jī)剛體運(yùn)動(dòng)受到嚴(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ī)氣動(dòng)參數(shù)分布

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

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

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

式中：M、Mx分別表示1/4 弦長處計(jì)算的力矩、1/2 弦長處計(jì)算的力矩，由于機(jī)翼彈性軸位于機(jī)翼1/2 弦長處， 而氣動(dòng)力和力矩是在1/4 弦長處的氣動(dòng)中心計(jì)算的，所以1/4 弦長處的力和力矩必須移到1/2 弦長處；ρ 是大氣密度；α 是局部迎角；t 時(shí)刻的局部迎角α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ù)表示的含義與飛行動(dòng)力學(xué)中表示的含義不同。 例如，在飛行動(dòng)力學(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ī)翼單元俯仰力矩的影響。

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

式中：CLαw、CLαh分別表示機(jī)翼升力線斜率和平尾升力線斜率，Sw、Sh分別表示機(jī)翼面積和平尾面積；CLαi表示機(jī)翼第i 個(gè)單元的升力線斜率；Ai表示第i 個(gè)單元的面積；?表示下洗。表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 彈性部件其他氣動(dòng)參數(shù)Tab.5 Aerodynamic parameters of other elastic components

3 仿真結(jié)果分析

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

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

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

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

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

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

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

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

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

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

4 結(jié)語

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