陳 偉，馮高鵬，拜云山

(中國(guó)工程物理研究院 總體工程研究所，四川 綿陽(yáng) 621900)

小型變體無人機(jī)動(dòng)力學(xué)建模及配平分析

陳 偉，馮高鵬，拜云山

(中國(guó)工程物理研究院 總體工程研究所，四川 綿陽(yáng) 621900)

為了更好地分析變體無人機(jī)(Unmanned Aerial Vehicle，UAV)在變體過程中的動(dòng)態(tài)響應(yīng)，針對(duì)一種小型變后掠翼無人機(jī)，通過靜矩來描述變體過程中重心位置的變化，推導(dǎo)了其六自由度非線性運(yùn)動(dòng)方程。變體飛機(jī)相對(duì)于固定翼飛機(jī)的運(yùn)動(dòng)學(xué)方程多出變體過程中引起的慣性力和慣性力矩項(xiàng)以及轉(zhuǎn)動(dòng)慣量和慣性積的微分項(xiàng)。利用準(zhǔn)定常假設(shè)來簡(jiǎn)化計(jì)算變體過程中的氣動(dòng)力，認(rèn)為在變體的任一瞬時(shí)，全機(jī)的氣動(dòng)力取決于該瞬時(shí)飛機(jī)的靜態(tài)構(gòu)型和飛行狀態(tài)。從計(jì)算結(jié)果可知，在一定飛行速度和飛行高度范圍內(nèi)，影響氣動(dòng)參數(shù)的主要變量為機(jī)翼后掠角和迎角。為了確保變體飛機(jī)在變體過程中具有較理想的靜穩(wěn)定度，通過配平給出了不同構(gòu)型時(shí)飛機(jī)重心在機(jī)體縱軸上的合理位置，有助于變體飛機(jī)控制系統(tǒng)設(shè)計(jì)。

變體無人機(jī)；運(yùn)動(dòng)方程；動(dòng)態(tài)響應(yīng)；配平

0 引 言

變體飛機(jī)在變體過程中，氣動(dòng)外形的變化必然引起飛機(jī)的轉(zhuǎn)動(dòng)慣量、質(zhì)心位置、氣動(dòng)焦點(diǎn)位置等一系列機(jī)體特性的改變。近年來，一些學(xué)者針對(duì)特定的變體飛機(jī)進(jìn)行了動(dòng)力學(xué)建模與分析[1-3]。陳偉等對(duì)變后掠翼飛機(jī)進(jìn)行了建模與分析[4]。樂挺等對(duì)Z型翼變體飛機(jī)進(jìn)行了多體動(dòng)力學(xué)建模，分析了質(zhì)心位置移動(dòng)和氣動(dòng)外形變化對(duì)變體飛機(jī)動(dòng)態(tài)特性的影響，并研究了非定常效應(yīng)對(duì)機(jī)翼折疊過程中氣動(dòng)特性的影響，認(rèn)為當(dāng)機(jī)翼折疊角速度較小時(shí)，可以忽略非定常氣動(dòng)力效

應(yīng)[5-6]。Obradovic B等采用擴(kuò)展剛體動(dòng)力學(xué)方法得到了一種機(jī)翼可上下擺動(dòng)的變體飛機(jī)的動(dòng)力學(xué)方程[7-8]。Seigler T M等建立了固定于機(jī)身的機(jī)體坐標(biāo)軸系，采用笛卡兒矢量描述變體飛機(jī)上各點(diǎn)相對(duì)于地面坐標(biāo)軸系的坐標(biāo)，利用達(dá)朗貝爾原理推導(dǎo)了變體飛機(jī)的一般動(dòng)態(tài)特性[9-11]。

文中在Seigler T M等研究的基礎(chǔ)上，推導(dǎo)了一種小型變體無人機(jī)的非線性運(yùn)動(dòng)方程。利用準(zhǔn)定常假設(shè)來簡(jiǎn)化計(jì)算變體過程中的氣動(dòng)力，為了確保變體飛機(jī)在變體過程中具有較理想的靜穩(wěn)定度，通過配平給出了不同構(gòu)型時(shí)飛機(jī)重心在機(jī)體縱軸上的合理位置，為變體飛機(jī)控制系統(tǒng)設(shè)計(jì)提供了基礎(chǔ)。

1 變體飛機(jī)變體過程中的動(dòng)力學(xué)建模

變體飛機(jī)在變體過程中，質(zhì)心位置會(huì)發(fā)生改變，而質(zhì)心速度不變。研究飛機(jī)上某固定點(diǎn)的速度要比質(zhì)心的速度更有意義[10]。建立與機(jī)身固定的機(jī)體坐標(biāo)軸系Sb-Oxyz，如圖1所示。

圖1 坐標(biāo)系

其中，Sg-Ogxgygzg為地面坐標(biāo)軸系。變體飛機(jī)由μm個(gè)相互運(yùn)動(dòng)的部分組成，各個(gè)運(yùn)動(dòng)部分質(zhì)量為mi，質(zhì)心為ci。ri和Ri分別為ci在機(jī)體坐標(biāo)軸系和地面坐標(biāo)系中的坐標(biāo)；R0為機(jī)體坐標(biāo)軸系原點(diǎn)在地面坐標(biāo)系中的坐標(biāo)?；谝陨纤⒌臋C(jī)體坐標(biāo)軸系，SeiglerTM推導(dǎo)了變體飛機(jī)的力方程和力矩方程[10]。

(1)

(2)

(3)

⌒為運(yùn)算符，將矢量轉(zhuǎn)換成反對(duì)稱矩陣，有：

(4)

(5)

其中，Ix、Iy和Iz為繞機(jī)體坐標(biāo)軸系的轉(zhuǎn)動(dòng)慣量；Ixy,Iyz和Izx為慣性積。

Fb為飛機(jī)受到的合外力在機(jī)體坐標(biāo)軸系中的矢量。

(6)

其中，Tbg為地面坐標(biāo)軸系到機(jī)體坐標(biāo)軸系的轉(zhuǎn)換矩陣，Tba為氣流坐標(biāo)軸系到機(jī)體坐標(biāo)軸系的轉(zhuǎn)換矩陣，具體坐標(biāo)系定義及轉(zhuǎn)換矩陣表達(dá)式見文獻(xiàn)[12]；g為重力加速度；T為推力，與機(jī)體x軸方向平行；L為升力，Y為側(cè)力，D為阻力；Mb為Fb產(chǎn)生的力矩在機(jī)體坐標(biāo)軸系中的矢量：

(7)

2 六自由度非線性運(yùn)動(dòng)方程

基于圖2所示的小型變體無人機(jī)進(jìn)行研究。圖2中給出了機(jī)翼后掠角為0°、30°和45°時(shí)的構(gòu)型。機(jī)翼后掠角可以在0°～45°之間任意變化，與MFX-2變體飛機(jī)的機(jī)翼變形機(jī)構(gòu)類似[13]。該變體飛機(jī)在機(jī)翼后掠角變化過程中，機(jī)翼寬度會(huì)按照一定的比例隨著機(jī)翼后掠角一起變化，即機(jī)翼后掠角確定，機(jī)翼氣動(dòng)弦長(zhǎng)和機(jī)翼面積確定，可認(rèn)為變體自由度為1。在機(jī)翼后掠角最大時(shí)，機(jī)翼寬度最小。為了方便敘述，文中用機(jī)翼后掠角表示飛機(jī)構(gòu)型，如機(jī)翼后掠角為0°時(shí)，簡(jiǎn)稱0°構(gòu)型。表1給出了變體飛機(jī)在0°構(gòu)型、30°構(gòu)型和45°構(gòu)型時(shí)的外形參數(shù)。飛機(jī)質(zhì)量m=25kg。

圖2 變體飛機(jī)

表1中，Sw為機(jī)翼參考面積；cA為機(jī)翼的平均幾何弦長(zhǎng)；b為機(jī)翼展長(zhǎng)。變體飛機(jī)從0°構(gòu)型向45°構(gòu)型變化過程中，飛機(jī)的機(jī)翼參考面積和展長(zhǎng)逐漸減小，機(jī)翼的參考面積最大可改變36%，機(jī)翼展長(zhǎng)最大可改變25.7%，機(jī)翼平均幾何弦長(zhǎng)最大可改變15%。

表1 變體飛機(jī)的不同構(gòu)型參數(shù)

在機(jī)翼后掠角變化過程中，飛機(jī)質(zhì)心和氣動(dòng)焦點(diǎn)沿機(jī)體x軸移動(dòng)。設(shè)左機(jī)翼質(zhì)心在機(jī)體坐標(biāo)軸系中的坐標(biāo)為r1=r1xi+r1yj+0k，質(zhì)量為m1；右機(jī)翼質(zhì)心坐標(biāo)為r2=r2xi+r2yj+0k，質(zhì)量為m2；機(jī)身加機(jī)尾的質(zhì)心坐標(biāo)為r3=r3xi+0j+0k，質(zhì)量為m3。因?yàn)轱w機(jī)左右對(duì)稱，則有m1=m2,r1x=r2x,r1y=-r2y,Ixy=Iyz=0。飛機(jī)的轉(zhuǎn)動(dòng)慣量張量I為:

(8)

飛機(jī)的靜矩S為:

(9)

其中，Sx≈2m1r1x+m3r3x,Sy=Sz=0，有:

(10)

將式(8)和式(10)代入式(1)和式(2)得到該小型變體飛機(jī)的力和力矩方程，按照文獻(xiàn)[12]中推導(dǎo)固定翼飛機(jī)非線性運(yùn)動(dòng)方程的步驟，推導(dǎo)出變體飛機(jī)六自由運(yùn)動(dòng)方程:

xg=Vcosγcosχ

(11)

yg=Vcosγsinχ

(12)

(13)

(14)

(15)

(16)

(17)

(18)

(19)

(20)

(21)

(22)

其中，xg、yg、zg為飛機(jī)在地面坐標(biāo)軸系的位置分量；V為飛行速度；h為飛行高度；χ為航跡方位角；γ為航跡傾斜角；μ為航跡滾轉(zhuǎn)角；α為迎角；β為側(cè)滑角；θ為俯仰角；φ為滾轉(zhuǎn)角；ψ為偏航角；FIx、FIy、FIz為變體過程引起的慣性力在氣流坐標(biāo)軸系中的分量，有:

(23)

MIx、MIy、MIz為變體過程引起的慣性力矩在機(jī)體坐標(biāo)軸系中的分量，有:

(24)

從變體飛機(jī)運(yùn)動(dòng)方程來看，變體飛機(jī)相對(duì)于固定翼飛機(jī)的運(yùn)動(dòng)學(xué)方程多出了變體過程中引起的慣性力和慣性力矩項(xiàng)以及轉(zhuǎn)動(dòng)慣量和慣性積的微分項(xiàng)。在機(jī)翼后掠角變化過程中，飛機(jī)重心沿機(jī)體x軸移動(dòng)，靜矩Sx會(huì)隨之變化。當(dāng)靜矩Sx不等于零時(shí)，重力相對(duì)于機(jī)體坐標(biāo)軸系原點(diǎn)會(huì)產(chǎn)生力矩，在力矩方程組中分別體現(xiàn)在Sxgcosθsinφ、-Sxgcosθcosφ、Sxgcosθsinφ三項(xiàng)中。另外，在機(jī)翼后掠角變化過程中，飛機(jī)氣動(dòng)外形和氣動(dòng)焦點(diǎn)會(huì)發(fā)生改變，飛機(jī)的氣動(dòng)參數(shù)會(huì)發(fā)生較大變化。

3 不同構(gòu)型的配平狀態(tài)

利用準(zhǔn)定常假設(shè)來簡(jiǎn)化計(jì)算變體過程中的氣動(dòng)力，認(rèn)為在變體的任一瞬時(shí)，全機(jī)的氣動(dòng)力取決于該瞬時(shí)飛機(jī)的靜態(tài)構(gòu)型和飛行狀態(tài)。

表2給出了在飛行速度V=25m/s、飛行高度h=1 000m、迎角α=0°時(shí)變體飛機(jī)不同構(gòu)型的部分氣動(dòng)參數(shù)。

表2 變體飛機(jī)不同構(gòu)型氣動(dòng)參數(shù)

其中，CL0為零迎角升力系數(shù)；CDα為阻力系數(shù)相對(duì)迎角的導(dǎo)數(shù)；Cm0為零迎角俯仰力矩系數(shù)；CYβ為側(cè)力系數(shù)對(duì)側(cè)滑角的導(dǎo)數(shù)；Clβ為滾轉(zhuǎn)力矩系數(shù)對(duì)側(cè)滑角的導(dǎo)數(shù)；Cnβ為偏航力矩系數(shù)對(duì)側(cè)滑角的導(dǎo)數(shù)。

對(duì)于該小型變體無人機(jī)來說，由于受到發(fā)動(dòng)機(jī)推力和自身設(shè)計(jì)參數(shù)的限制，飛行速度和飛行高度被限制在一定范圍內(nèi)，在該范圍內(nèi)飛行速度和飛行高度對(duì)氣動(dòng)參數(shù)的影響較小。主要影響氣動(dòng)參數(shù)的變量為機(jī)翼后掠角和迎角。圖3給出了在V=25m/s、h=1 000m、δe=0°、β=0°、q=0 (°)/s時(shí)的升力系數(shù)CL、阻力系數(shù)CD、俯仰力矩系數(shù)Cm隨機(jī)翼后掠角和迎角變化的擬合曲面。在迎角小于8°時(shí)，CL與迎角具有較好的線性關(guān)系，升力線斜率CLα可近似為常值，隨著迎角的增大，CLα逐漸減小最終變?yōu)樨?fù)值，CLα等于零對(duì)應(yīng)的迎角為失速迎角。0°構(gòu)型時(shí)的失速迎角在14°左右，45°構(gòu)型時(shí)的失速迎角要大于16°。CD隨迎角的增大而增大。在機(jī)翼后掠角變化過程中，飛機(jī)重心始終在氣動(dòng)焦點(diǎn)之前，飛機(jī)是縱向靜穩(wěn)定的。

圖3 氣動(dòng)參數(shù)隨機(jī)翼后掠角和迎角變化的擬合曲面

在水平飛行狀態(tài)，對(duì)變體飛機(jī)的不同構(gòu)型進(jìn)行配平。配平代價(jià)函數(shù)為

(25)

在一定的飛行速度V和飛行高度h下，通過計(jì)算配平迎角α0、配平俯仰角θ0、配平推力T0、配平升降舵偏角δe0使得代價(jià)函數(shù)值趨近于0，即飛機(jī)所受到的力和力矩平衡。由于變體飛機(jī)在變體的過程中，氣動(dòng)焦點(diǎn)在機(jī)體x軸上的位置變化較大，為了確保變體飛機(jī)具有較理想的靜穩(wěn)定度，隨著機(jī)翼后掠角的改變應(yīng)當(dāng)通過配重塊將飛機(jī)重心配置在合理位置上，即氣動(dòng)焦點(diǎn)在重心之后的距離與平均氣動(dòng)弦長(zhǎng)的比值在10%附近。這里通過配平給出了不同構(gòu)型時(shí)重心在機(jī)體x軸上相對(duì)建模原點(diǎn)O的位置，如表3所示，同時(shí)表3給出了在飛行高度h=1 000m、飛行速度分別為V=25m/s、V=30m/s、V=40m/s和V=50m/s時(shí)，變體飛機(jī)在0°構(gòu)型、30°構(gòu)型和45°構(gòu)型時(shí)的配平值。

表3 配平狀態(tài)

表3中，由于機(jī)翼面積隨著機(jī)翼后掠角增大而減小，在飛行速度一定時(shí)為了提供相同的平飛升力，配平迎角會(huì)隨著機(jī)翼后掠角增大而增大。由于配平迎角的增大，為了平衡迎角產(chǎn)生的低頭力矩，配平升降舵偏角隨著機(jī)翼后掠角增大而上偏。當(dāng)飛行速度較高時(shí)，由于配平迎角較小，同時(shí)阻力系數(shù)對(duì)迎角的導(dǎo)數(shù)和機(jī)翼面積隨著機(jī)翼后掠角增大而減小，配平推力會(huì)隨著機(jī)翼后掠角增大而減小。當(dāng)V=25m/s時(shí)，45°構(gòu)型的配平迎角達(dá)到了9.204°，在該飛行速度下做機(jī)動(dòng)動(dòng)作時(shí)迎角可能會(huì)達(dá)到失速迎角，表明45°構(gòu)型不適合在低速下飛行。根據(jù)以上分析可知，在低速巡航時(shí)，為了使升阻比最大，增大巡航時(shí)間，可采用0°構(gòu)型或者小后掠角構(gòu)型飛行；在快速進(jìn)入戰(zhàn)區(qū)或逃離戰(zhàn)區(qū)時(shí)，為了減小飛行阻力，可采用45°構(gòu)型飛行。

4 結(jié)束語(yǔ)

文中推導(dǎo)了一種變后掠翼無人機(jī)的非線性運(yùn)動(dòng)方程。利用準(zhǔn)定常假設(shè)來簡(jiǎn)化計(jì)算變體過程中的氣動(dòng)力，認(rèn)為在變體的任一瞬時(shí)，全機(jī)的氣動(dòng)力取決于該瞬時(shí)飛機(jī)的靜態(tài)構(gòu)型和飛行狀態(tài)。為了確保變體飛機(jī)在變體過程中具有較理想的靜穩(wěn)定度，通過配平給出了不同構(gòu)型時(shí)飛機(jī)重心在機(jī)體縱軸上的合理位置。為了確保變體過程中的飛行穩(wěn)定性，所設(shè)計(jì)的飛行控制系統(tǒng)應(yīng)當(dāng)較好地抑制重心位置變化帶來的影響。

[1] 陳 偉,盧京潮,王曉光,等.基于backstepping/RHO的變體飛機(jī)控制器設(shè)計(jì)[J].北京航空航天大學(xué)學(xué)報(bào),2014,40(8):1060-1065.

[2] 陳 偉,盧京潮,章衛(wèi)國(guó).基于LQ/RHO多回路技術(shù)的變體飛機(jī)控制器設(shè)計(jì)[J].北京航空航天大學(xué)學(xué)報(bào),2013,39(3):290-294.

[3] 陳 偉,盧京潮,劉志君,等.基于指令濾波器的變體飛機(jī)RHO控制器設(shè)計(jì)[J].飛行力學(xué),2014,32(5):423-426.

[4] 陳 偉,馮高鵬.變體飛機(jī)建模及自適應(yīng)動(dòng)態(tài)面控制[J].測(cè)控技術(shù),2016,35(2):91-95.

[5] 樂 挺,王立新,艾俊強(qiáng).Z型翼變體飛機(jī)的縱向多體動(dòng)力學(xué)特性[J].航空學(xué)報(bào),2010,31(4):679-686.

[6] 樂 挺,王立新,艾俊強(qiáng).變體飛機(jī)設(shè)計(jì)的主要關(guān)鍵技術(shù)[J].飛行力學(xué),2009,27(5):6-10.

[7]ObradovicB,SubbaraoK.Modelingofflightdynamicsofmorphing-wingaircraft[J].JournalofAircraft,2011,48(2):391-402.

[8]ObradovicB,SubbaraoK.Modelingofdynamicloadingofmorphing-wingaircraft[J].JournalofAircraft,2011,48(2):424-435.

[9]SeiglerTM.Dynamicsandcontrolofmorphingaircraft[D].Virginia:theVirginiaPolytechnicInstituteandStateUniversity,2005.

[10]SeiglerTM,NealDA,BaeJS,etal.Modelingandflightcontroloflarge-scalemorphingaircraft[J].JournalofAircraft,2007,44(4):1077-1087.

[11]SeiglerTM,NealDA,InmanDJ.Dynamicmodelingoflarge-scalemorphingaircraft[C]//AIAAstructures,structuraldynamics,andmaterialsconference.Newport:AmericanInstituteofAeronauticsandAstronautics,2006:3668-3678.

[12] 吳森堂.飛行控制系統(tǒng)[M].北京:北京航空航天大學(xué)出版社,2005:20-70.

[13]GandhiN,CooperJ,WardD,etal.Ahardwaredemonstrationofanintegratedadaptivewingshapeandflightcontrollawformorphingaircraft[C]//AIAAguidance,navigation,andcontrolconference.Chicago:AmericanInstituteofAeronauticsandAstronautics,2009.

Dynamics Modeling and Trimming Analysis for Small Morphing UAV

CHEN Wei,FENG Gao-peng,BAI Yun-shan

(Institute of Systems Engineering，China Academy of Engineering Physics，Mianyang 621900,China)

In order to better analyze the dynamic response of morphing Unmanned Aerial Vehicle (UAV),the moving of the center of gravity is described by quiet moment,and the six free nonlinear equations of motion for a small variable sweep UAV is derived.The kinematic equation of morphing aircraft has more items that are inertial force and inertial moment than fixed wing plane.The aerodynamic forces and moments in the process of morphing are calculated with the quasi-steady assumption,the aerodynamic force is depended on the static configuration and flight state of morphing aircraft at the moment of morphing.From the results it can be known that the main variable that affects the aerodynamic parameters ie sweep angle and angle of attack within a certain range of flight speed and altitude.To make sure that the morphing aircraft has an ideal static stability margin,the responsible position of the center of gravity on the longitudinal axis of the body at different configurations is obtained by trimming,which lays a foundation for the control system design for the morphing aircraft.

morphing UAV;equations of motion;dynamic response;trimming

2016-03-26

2016-06-29

時(shí)間：2017-01-04

國(guó)家自然基金聯(lián)合基金(U1430113)；中國(guó)工程物理研究院科學(xué)技術(shù)發(fā)展基金(2014B0101009)

陳 偉(1986-)，男，博士，工程師，研究方向?yàn)橹茖?dǎo)武器控制系統(tǒng)設(shè)計(jì)；馮高鵬，通訊作者，高級(jí)工程師，研究方向?yàn)槲淦飨到y(tǒng)設(shè)計(jì)。

http://www.cnki.net/kcms/detail/61.1450.TP.20170104.1039.068.html

TP273

A

1673-629X(2017)02-0125-05

10.3969/j.issn.1673-629X.2017.02.028