孔雀垂直起降飛行器數(shù)學(xué)仿真模型*

賈晨輝 駱無(wú)意 柳嘉潤(rùn) 呂新廣 李新明 鞏慶海 張 雋

北京航天自動(dòng)控制研究所，宇航智能控制技術(shù)國(guó)家級(jí)重點(diǎn)實(shí)驗(yàn)室，北京 100854

0 引言

可重復(fù)使用運(yùn)載火箭是近年來(lái)新興的一項(xiàng)航天飛行控制技術(shù)。美國(guó)當(dāng)?shù)貢r(shí)間2015年12月21日晚，特斯拉公司創(chuàng)始人埃隆·馬斯克創(chuàng)辦的太空探索技術(shù)公司(以下簡(jiǎn)稱(chēng)SpaceX)發(fā)射的“獵鷹9”火箭在佛羅里達(dá)州卡納維拉爾角成功實(shí)現(xiàn)第一節(jié)火箭軟著陸，從而開(kāi)創(chuàng)了火箭從太空直接垂直回收的歷史[1]。近年來(lái)，火箭的垂直回收技術(shù)受到了國(guó)內(nèi)外多個(gè)高校、研究機(jī)構(gòu)的廣泛關(guān)注，該技術(shù)不僅極大降低了火箭的發(fā)射成本，使火箭具備了可重復(fù)使用的能力，而且可以使火箭帶回更多的飛行數(shù)據(jù)供研究人員參考、分析，為提升后續(xù)火箭的各項(xiàng)性能提供數(shù)據(jù)基礎(chǔ)。

“孔雀”系列飛行器為一款可重復(fù)使用的垂直起降技術(shù)驗(yàn)證平臺(tái)，可用于驗(yàn)證運(yùn)載火箭垂直回收、非程序制導(dǎo)及智能飛行控制等各項(xiàng)關(guān)鍵技術(shù)[2]。。該飛行器目前已服役3年，參加了多次飛行驗(yàn)證試驗(yàn)，實(shí)現(xiàn)了多項(xiàng)關(guān)鍵技術(shù)的突破。

為支持航空航天飛行器的飛行控制研究，各研究機(jī)構(gòu)建立了多種飛行器的數(shù)學(xué)模型供航空航天研究人員或愛(ài)好者進(jìn)行控制仿真研究。如NASA公布了F-18 HARV戰(zhàn)斗機(jī)數(shù)學(xué)模型[3]，以及一些學(xué)者公布了一系列超音速飛行器的數(shù)學(xué)模型[4-5]等。但由于運(yùn)載火箭垂直回收技術(shù)研究剛剛起步，目前鮮有對(duì)外發(fā)布的可重復(fù)使用運(yùn)載火箭數(shù)學(xué)模型供航天工作者進(jìn)行運(yùn)載火箭垂直回收等技術(shù)的控制仿真分析與研究工作。

本文建立孔雀飛行器的數(shù)學(xué)模型用于為垂直起降飛行器控制技術(shù)研究人員提供一套可進(jìn)行控制仿真研究的被控對(duì)象模型，該數(shù)學(xué)模型根據(jù)孔雀飛行器的設(shè)計(jì)過(guò)程及飛行試驗(yàn)過(guò)程的參數(shù)建立，作為孔雀飛行試驗(yàn)前的數(shù)學(xué)仿真和半實(shí)物仿真依據(jù)，實(shí)現(xiàn)孔雀飛行器各項(xiàng)控制算法的驗(yàn)證，亦可作為可重復(fù)使用運(yùn)載火箭飛行控制仿真研究的被控對(duì)象模型使用。

1 符號(hào)說(shuō)明與基本假設(shè)

1.1 符號(hào)說(shuō)明

孔雀飛行器模型的常用符號(hào)如表1所示。

1.2 基本假設(shè)

1)飛行器為剛體，忽略彈性振動(dòng)；

2)平坦靜止大地，忽略地球曲率和地球轉(zhuǎn)動(dòng)的影響；

3)忽略燃油消耗導(dǎo)致的轉(zhuǎn)動(dòng)慣量、慣性積、質(zhì)心位置變化；

4)忽略發(fā)動(dòng)機(jī)轉(zhuǎn)子轉(zhuǎn)動(dòng)導(dǎo)致的陀螺效應(yīng)；

5)忽略發(fā)動(dòng)機(jī)擺動(dòng)角加速度導(dǎo)致的附加力矩。

2 被控對(duì)象模型

2.1 孔雀飛行器外形

孔雀飛行器外形及主要艙段如圖1所示。本文中所介紹的模型為不包含涵道風(fēng)扇與柵格舵(僅計(jì)算二者重量)的飛行器模型。

表1 孔雀飛行器模型的常用符號(hào)

圖1 孔雀飛行器外形圖

2.2 坐標(biāo)系

采用俄羅斯體制的“北-天-東”坐標(biāo)系，定義符合右手定則的空間直角坐標(biāo)系，其定義見(jiàn)參考文獻(xiàn)[6]。

地面系與箭體系之間的關(guān)系可以用俯仰角φ、偏航角ψ和滾轉(zhuǎn)角γ這3個(gè)姿態(tài)角來(lái)確定。姿態(tài)角按照從地面系到箭體系為3-2-1的轉(zhuǎn)序定義。

2.3 被控對(duì)象模型

2.3.1 風(fēng)模型

(1)

考慮飛行器處于水平定常風(fēng)場(chǎng)中，即僅考慮水平方向的風(fēng)。風(fēng)向角θw定義為：按右手定則，從地面系OdXd軸開(kāi)始，繞OdYd軸轉(zhuǎn)動(dòng)到風(fēng)速矢量的角度。地面系下的風(fēng)速按式(2)計(jì)算：

(2)

在一次仿真中風(fēng)向角θw為常數(shù)。

風(fēng)速大小vw設(shè)置為隨高度變化的仿真量，其表達(dá)式為：

vw=kwindyvw1≤vw≤vw2

(3)

其中，kwind為風(fēng)速相對(duì)于高度變化的比例系數(shù)；vw1，vw2為仿真中風(fēng)速的上下限幅值。

(4)

2.3.2 飛行器運(yùn)動(dòng)方程

飛行器運(yùn)動(dòng)方程的狀態(tài)為：地面系位置[x,y,z]T；地面系速度[vx,vy,vz]T；姿態(tài)角[φ,ψ,γ]T；箭體軸轉(zhuǎn)動(dòng)角速度[ωx,ωy,ωz]T；此外，發(fā)動(dòng)機(jī)耗油導(dǎo)致質(zhì)量變化。微分方程如式(4)所示，其中：

(5)

式中，Tt→d為箭體系到地面系的轉(zhuǎn)換矩陣；[FA,FN,FZ]T為氣動(dòng)力在箭體系下的分量表示；[Pxt,1,Pyt,1,Pzt,1]T；[Pxt,2,Pyt,2,Pzt,2]T分別為1號(hào)、2號(hào)發(fā)動(dòng)機(jī)推力矢量在箭體系的分量表示；mg為重力；Qeng,1，Qeng,2分別為1號(hào)、2號(hào)發(fā)動(dòng)機(jī)的耗油率，單位為kg/s。

氣動(dòng)力矩、推力力矩等的合外力矩在箭體系的分量[Mxt,Myt,Mzt]T可表示為：

(6)

式中，[L,M,N]T為箭體系下的氣動(dòng)力矩；[lx,1,ly,1,lz,1]T，[lx,2,ly,2,lz,2]T分別為1號(hào)、2號(hào)發(fā)動(dòng)機(jī)轉(zhuǎn)軸在箭體系相對(duì)于飛行器質(zhì)心的坐標(biāo)；?表示矢量的叉乘，其定義為：

(7)

飛行器地速大小v按下式計(jì)算：

(8)

地面系下的空速按下式計(jì)算：

(9)

箭體系下的空速按下式計(jì)算：

(10)

空速大小按下列2式計(jì)算均可：

(11)

(12)

飛行器攻角α和側(cè)滑角β根據(jù)下2式計(jì)算：

(13)

(14)

其中，α∈[-π，π]。

2.3.3 氣動(dòng)特性

氣動(dòng)力和力矩需要計(jì)算相應(yīng)的氣動(dòng)系數(shù)。記：

CA,CN,CZ,Cmx,Cmy和Cmz分別為軸向力系數(shù)、法向力系數(shù)、側(cè)向力系數(shù)、滾轉(zhuǎn)力矩系數(shù)、偏航力矩系數(shù)和俯仰力矩系數(shù)，正方向?yàn)槠渌?或所繞)箭體系坐標(biāo)軸的正方向。

氣動(dòng)系數(shù)是空速、攻角和側(cè)滑角非線性函數(shù)，其標(biāo)稱(chēng)值以多維數(shù)表的形式給出：

Ck=Ck(u,α,-β) (k=A,N,Z,mx,my,mz)

箭體系下氣動(dòng)力為：

Fk=qSrefCk(k=A,N,Z)

(15)

氣動(dòng)力矩系數(shù)的參考點(diǎn)為飛行器頂點(diǎn)。因此，箭體系下繞質(zhì)心的氣動(dòng)力矩為：

(16)

式中，Sref為氣動(dòng)參考面積;Lref為氣動(dòng)參考長(zhǎng)度;Xcg為質(zhì)心位置，從飛行器頂點(diǎn)向下為正；q為動(dòng)壓：

(17)

其中，ρ為大氣密度;u為空速大小。計(jì)算時(shí)，氣動(dòng)力系數(shù)、氣動(dòng)力矩系數(shù)、質(zhì)心位置和大氣密度均應(yīng)考慮偏差。

2.3.4 發(fā)動(dòng)機(jī)推力矢量

發(fā)動(dòng)機(jī)的擺動(dòng)由伺服電機(jī)驅(qū)動(dòng)實(shí)現(xiàn)。每臺(tái)發(fā)動(dòng)機(jī)由4個(gè)電機(jī)驅(qū)動(dòng)(可等效為2個(gè)舵機(jī))，每臺(tái)發(fā)動(dòng)機(jī)的2個(gè)舵機(jī)互不耦合，分別在正交的2個(gè)方向產(chǎn)生擺角，控制發(fā)動(dòng)機(jī)向任意方向偏轉(zhuǎn)。2臺(tái)發(fā)動(dòng)機(jī)呈“一”字形安裝。每臺(tái)發(fā)動(dòng)機(jī)有俯仰、偏航2個(gè)轉(zhuǎn)軸。其后視圖如圖2所示。

圖2 孔雀飛行器發(fā)動(dòng)機(jī)安裝與擺角定義示意圖

結(jié)合發(fā)動(dòng)機(jī)擺角的定義，2臺(tái)發(fā)動(dòng)機(jī)的推力矢量分別為：

(18)

其中，δins為安裝角，即擺角為零時(shí)發(fā)動(dòng)機(jī)推力線與箭體軸線夾角，以使得推力線由平行于箭體X軸偏向質(zhì)心方向。

2.4 執(zhí)行機(jī)構(gòu)模型

2.4.1 發(fā)動(dòng)機(jī)推力變化的動(dòng)態(tài)特性

“孔雀”飛行器發(fā)動(dòng)機(jī)采用JETCAT P550 Pro，發(fā)動(dòng)機(jī)推力變化特性可表示為如下形式：

PC=f(PWMeng)

(19)

其中推力指令是發(fā)動(dòng)機(jī)控制PWM信號(hào)PWMeng的函數(shù)，由數(shù)表插值得到，PWMeng是一個(gè)正整數(shù)。若某一自變量的值超出插值表，則以數(shù)表邊界值代替。

根據(jù)地面測(cè)試數(shù)據(jù)，將發(fā)動(dòng)機(jī)動(dòng)態(tài)特性簡(jiǎn)化為一個(gè)二階線性環(huán)節(jié)，用于計(jì)算實(shí)際推力。

(20)

式中，PC為推力指令;P為發(fā)動(dòng)機(jī)實(shí)際推力大小;Keng為發(fā)動(dòng)機(jī)穩(wěn)態(tài)增益;ωeng為發(fā)動(dòng)機(jī)自然頻率；ξeng為發(fā)動(dòng)機(jī)阻尼比。

2.4.2 驅(qū)動(dòng)發(fā)動(dòng)機(jī)擺動(dòng)的伺服電機(jī)動(dòng)態(tài)特性

將電動(dòng)舵機(jī)的動(dòng)態(tài)特性簡(jiǎn)化為二階傳遞函數(shù)、間隙特性與零位的串聯(lián)。

二階傳遞函數(shù)特性如下：

(21)

式中，δ*C(s)為擺角指令;δ*TF(s)為經(jīng)過(guò)傳遞函數(shù)后的舵機(jī)擺角。下標(biāo)*表示1,2,3,4。Kδ為伺服穩(wěn)態(tài)增益，ωδ為伺服自然頻率，ξδ為伺服阻尼比。

3 模型參數(shù)

本文所列各飛行器相關(guān)參數(shù)來(lái)源于設(shè)計(jì)過(guò)程中的計(jì)算數(shù)值或根據(jù)飛行試驗(yàn)數(shù)據(jù)所測(cè)得的值。

3.1 飛行器總體參數(shù)

孔雀飛行器模型參數(shù)見(jiàn)表2：

表2 孔雀飛行器總體參數(shù)標(biāo)稱(chēng)值

3.2 氣動(dòng)系數(shù)

分別在起落架收起上升段、無(wú)柵格舵起落架收起下降段2種狀態(tài)下，根據(jù)試驗(yàn)數(shù)據(jù)進(jìn)行擬合。飛行器在上升段與下降段，在起落架收起與起落架放下的狀態(tài)下，其氣動(dòng)特性均有不同。受篇幅限制，本文未列出模型仿真過(guò)程中所采用的氣動(dòng)參數(shù)，對(duì)此部分感興趣的讀者可通過(guò)作者郵箱進(jìn)行索取。

3.3 發(fā)動(dòng)機(jī)參數(shù)

發(fā)動(dòng)機(jī)參數(shù)的標(biāo)稱(chēng)值見(jiàn)表3。

表3 發(fā)動(dòng)機(jī)參數(shù)標(biāo)稱(chēng)值

3.4 伺服電機(jī)參數(shù)

伺服電機(jī)參數(shù)的標(biāo)稱(chēng)值見(jiàn)表4。其中伺服間隙半寬度大小為通過(guò)試驗(yàn)實(shí)測(cè)。

4 控制仿真實(shí)例

本仿真實(shí)例采用一套制導(dǎo)控制律對(duì)模型進(jìn)行控制，目標(biāo)為將“孔雀”飛行器由起點(diǎn)(0m,0m,0m)起飛，在終點(diǎn)(40m,0m,0m)處著陸，在飛行過(guò)程中以上升速度(Y向速度)10m/s經(jīng)過(guò)交班點(diǎn)(40m,320m,0m)，在經(jīng)過(guò)交班點(diǎn)之后開(kāi)始減速上升，至最高點(diǎn)后進(jìn)行加速下降，當(dāng)速度達(dá)到-7m/s時(shí)進(jìn)行減速下降直至著陸。受篇幅限制，本文未列出模型仿真過(guò)程中所采用的制導(dǎo)與控制算法，對(duì)此部分感興趣的讀者可通過(guò)作者郵箱進(jìn)行索取。

表4 伺服電機(jī)參數(shù)標(biāo)稱(chēng)值

飛行器運(yùn)動(dòng)方程與制導(dǎo)控制律均通過(guò)MFC環(huán)境進(jìn)行編程，對(duì)于模型中的微分方程，采用固定步長(zhǎng)的四階龍格-庫(kù)塔法進(jìn)行數(shù)值積分。仿真開(kāi)始時(shí)間為ts=0.0s，積分步長(zhǎng)為Δt=0.001s,控制周期ΔtC=0.01s。圖3～4分別為仿真過(guò)程中飛行器的位置、速度曲線?？梢?jiàn)按照所設(shè)定的制導(dǎo)控制律，可完成飛行器按照預(yù)定軌跡在47.32s左右飛至交班點(diǎn)，并垂直降落至指定位置。圖5～6為真實(shí)飛行試驗(yàn)中采用相同的制導(dǎo)控制律得出的真實(shí)飛行試驗(yàn)曲線。從曲線的趨勢(shì)可以看出，該模型能夠較為相似地反映實(shí)際飛行控制狀態(tài)。

圖3 飛行控制仿真位置速度曲線

圖4 飛行控制仿真姿態(tài)曲線

圖5 實(shí)際飛行試驗(yàn)中飛行器位置速度曲線

圖6 實(shí)際飛行試驗(yàn)中飛行器姿態(tài)曲線

5 結(jié)論

本文介紹了“孔雀”飛行器垂直起降飛行驗(yàn)證平臺(tái)的數(shù)學(xué)模型，以微分方程的形式描述了該模型的運(yùn)動(dòng)學(xué)方程、氣動(dòng)特性和發(fā)動(dòng)機(jī)及伺服機(jī)構(gòu)等執(zhí)行機(jī)構(gòu)特性，并基于本模型采用一套制導(dǎo)控制律進(jìn)行了飛行控制半實(shí)物仿真。該數(shù)學(xué)模型可作為被控對(duì)象用于運(yùn)載火箭垂直回收技術(shù)等飛行控制技術(shù)的仿真研究。