賈晨輝 駱無意 柳嘉潤 呂新廣 李新明 鞏慶海 張 雋
北京航天自動控制研究所,宇航智能控制技術國家級重點實驗室,北京 100854
可重復使用運載火箭是近年來新興的一項航天飛行控制技術。美國當地時間2015年12月21日晚,特斯拉公司創始人埃隆·馬斯克創辦的太空探索技術公司(以下簡稱SpaceX)發射的“獵鷹9”火箭在佛羅里達州卡納維拉爾角成功實現第一節火箭軟著陸,從而開創了火箭從太空直接垂直回收的歷史[1]。近年來,火箭的垂直回收技術受到了國內外多個高校、研究機構的廣泛關注,該技術不僅極大降低了火箭的發射成本,使火箭具備了可重復使用的能力,而且可以使火箭帶回更多的飛行數據供研究人員參考、分析,為提升后續火箭的各項性能提供數據基礎。
“孔雀”系列飛行器為一款可重復使用的垂直起降技術驗證平臺,可用于驗證運載火箭垂直回收、非程序制導及智能飛行控制等各項關鍵技術[2]。。該飛行器目前已服役3年,參加了多次飛行驗證試驗,實現了多項關鍵技術的突破。
為支持航空航天飛行器的飛行控制研究,各研究機構建立了多種飛行器的數學模型供航空航天研究人員或愛好者進行控制仿真研究。如NASA公布了F-18 HARV戰斗機數學模型[3],以及一些學者公布了一系列超音速飛行器的數學模型[4-5]等。但由于運載火箭垂直回收技術研究剛剛起步,目前鮮有對外發布的可重復使用運載火箭數學模型供航天工作者進行運載火箭垂直回收等技術的控制仿真分析與研究工作。
本文建立孔雀飛行器的數學模型用于為垂直起降飛行器控制技術研究人員提供一套可進行控制仿真研究的被控對象模型,該數學模型根據孔雀飛行器的設計過程及飛行試驗過程的參數建立,作為孔雀飛行試驗前的數學仿真和半實物仿真依據,實現孔雀飛行器各項控制算法的驗證,亦可作為可重復使用運載火箭飛行控制仿真研究的被控對象模型使用。
孔雀飛行器模型的常用符號如表1所示。
1)飛行器為剛體,忽略彈性振動;
2)平坦靜止大地,忽略地球曲率和地球轉動的影響;
3)忽略燃油消耗導致的轉動慣量、慣性積、質心位置變化;
4)忽略發動機轉子轉動導致的陀螺效應;
5)忽略發動機擺動角加速度導致的附加力矩。
孔雀飛行器外形及主要艙段如圖1所示。本文中所介紹的模型為不包含涵道風扇與柵格舵(僅計算二者重量)的飛行器模型。

表1 孔雀飛行器模型的常用符號

圖1 孔雀飛行器外形圖
采用俄羅斯體制的“北-天-東”坐標系,定義符合右手定則的空間直角坐標系,其定義見參考文獻[6]。
地面系與箭體系之間的關系可以用俯仰角φ、偏航角ψ和滾轉角γ這3個姿態角來確定。姿態角按照從地面系到箭體系為3-2-1的轉序定義。
2.3.1 風模型


(1)
考慮飛行器處于水平定常風場中,即僅考慮水平方向的風。風向角θw定義為:按右手定則,從地面系OdXd軸開始,繞OdYd軸轉動到風速矢量的角度。地面系下的風速按式(2)計算:
(2)
在一次仿真中風向角θw為常數。
風速大小vw設置為隨高度變化的仿真量,其表達式為:
vw=kwindyvw1≤vw≤vw2
(3)
其中,kwind為風速相對于高度變化的比例系數;vw1,vw2為仿真中風速的上下限幅值。

(4)
2.3.2 飛行器運動方程
飛行器運動方程的狀態為:地面系位置[x,y,z]T;地面系速度[vx,vy,vz]T;姿態角[φ,ψ,γ]T;箭體軸轉動角速度[ωx,ωy,ωz]T;此外,發動機耗油導致質量變化。微分方程如式(4)所示,其中:

(5)
式中,Tt→d為箭體系到地面系的轉換矩陣;[FA,FN,FZ]T為氣動力在箭體系下的分量表示;[Pxt,1,Pyt,1,Pzt,1]T;[Pxt,2,Pyt,2,Pzt,2]T分別為1號、2號發動機推力矢量在箭體系的分量表示;mg為重力;Qeng,1,Qeng,2分別為1號、2號發動機的耗油率,單位為kg/s。
氣動力矩、推力力矩等的合外力矩在箭體系的分量[Mxt,Myt,Mzt]T可表示為:
(6)
式中,[L,M,N]T為箭體系下的氣動力矩;[lx,1,ly,1,lz,1]T,[lx,2,ly,2,lz,2]T分別為1號、2號發動機轉軸在箭體系相對于飛行器質心的坐標;?表示矢量的叉乘,其定義為:
(7)
飛行器地速大小v按下式計算:
(8)
地面系下的空速按下式計算:
(9)
箭體系下的空速按下式計算:
(10)
空速大小按下列2式計算均可:
(11)
(12)
飛行器攻角α和側滑角β根據下2式計算:

(13)
(14)
其中,α∈[-π,π]。
2.3.3 氣動特性
氣動力和力矩需要計算相應的氣動系數。記:
CA,CN,CZ,Cmx,Cmy和Cmz分別為軸向力系數、法向力系數、側向力系數、滾轉力矩系數、偏航力矩系數和俯仰力矩系數,正方向為其所在(或所繞)箭體系坐標軸的正方向。
氣動系數是空速、攻角和側滑角非線性函數,其標稱值以多維數表的形式給出:
Ck=Ck(u,α,-β) (k=A,N,Z,mx,my,mz)
箭體系下氣動力為:
Fk=qSrefCk(k=A,N,Z)
(15)
氣動力矩系數的參考點為飛行器頂點。因此,箭體系下繞質心的氣動力矩為:

(16)
式中,Sref為氣動參考面積;Lref為氣動參考長度;Xcg為質心位置,從飛行器頂點向下為正;q為動壓:
(17)
其中,ρ為大氣密度;u為空速大小。計算時,氣動力系數、氣動力矩系數、質心位置和大氣密度均應考慮偏差。
2.3.4 發動機推力矢量
發動機的擺動由伺服電機驅動實現。每臺發動機由4個電機驅動(可等效為2個舵機),每臺發動機的2個舵機互不耦合,分別在正交的2個方向產生擺角,控制發動機向任意方向偏轉。2臺發動機呈“一”字形安裝。每臺發動機有俯仰、偏航2個轉軸。其后視圖如圖2所示。

圖2 孔雀飛行器發動機安裝與擺角定義示意圖
結合發動機擺角的定義,2臺發動機的推力矢量分別為:
(18)
其中,δins為安裝角,即擺角為零時發動機推力線與箭體軸線夾角,以使得推力線由平行于箭體X軸偏向質心方向。
2.4.1 發動機推力變化的動態特性
“孔雀”飛行器發動機采用JETCAT P550 Pro,發動機推力變化特性可表示為如下形式:
PC=f(PWMeng)
(19)
其中推力指令是發動機控制PWM信號PWMeng的函數,由數表插值得到,PWMeng是一個正整數。若某一自變量的值超出插值表,則以數表邊界值代替。
根據地面測試數據,將發動機動態特性簡化為一個二階線性環節,用于計算實際推力。
(20)
式中,PC為推力指令;P為發動機實際推力大小;Keng為發動機穩態增益;ωeng為發動機自然頻率;ξeng為發動機阻尼比。
2.4.2 驅動發動機擺動的伺服電機動態特性
將電動舵機的動態特性簡化為二階傳遞函數、間隙特性與零位的串聯。
二階傳遞函數特性如下:
(21)
式中,δ*C(s)為擺角指令;δ*TF(s)為經過傳遞函數后的舵機擺角。下標*表示1,2,3,4。Kδ為伺服穩態增益,ωδ為伺服自然頻率,ξδ為伺服阻尼比。
本文所列各飛行器相關參數來源于設計過程中的計算數值或根據飛行試驗數據所測得的值。
孔雀飛行器模型參數見表2:

表2 孔雀飛行器總體參數標稱值
分別在起落架收起上升段、無柵格舵起落架收起下降段2種狀態下,根據試驗數據進行擬合。飛行器在上升段與下降段,在起落架收起與起落架放下的狀態下,其氣動特性均有不同。受篇幅限制,本文未列出模型仿真過程中所采用的氣動參數,對此部分感興趣的讀者可通過作者郵箱進行索取。
發動機參數的標稱值見表3。

表3 發動機參數標稱值
伺服電機參數的標稱值見表4。其中伺服間隙半寬度大小為通過試驗實測。
本仿真實例采用一套制導控制律對模型進行控制,目標為將“孔雀”飛行器由起點(0m,0m,0m)起飛,在終點(40m,0m,0m)處著陸,在飛行過程中以上升速度(Y向速度)10m/s經過交班點(40m,320m,0m),在經過交班點之后開始減速上升,至最高點后進行加速下降,當速度達到-7m/s時進行減速下降直至著陸。受篇幅限制,本文未列出模型仿真過程中所采用的制導與控制算法,對此部分感興趣的讀者可通過作者郵箱進行索取。

表4 伺服電機參數標稱值
飛行器運動方程與制導控制律均通過MFC環境進行編程,對于模型中的微分方程,采用固定步長的四階龍格-庫塔法進行數值積分。仿真開始時間為ts=0.0s,積分步長為Δt=0.001s,控制周期ΔtC=0.01s。圖3~4分別為仿真過程中飛行器的位置、速度曲線。可見按照所設定的制導控制律,可完成飛行器按照預定軌跡在47.32s左右飛至交班點,并垂直降落至指定位置。圖5~6為真實飛行試驗中采用相同的制導控制律得出的真實飛行試驗曲線。從曲線的趨勢可以看出,該模型能夠較為相似地反映實際飛行控制狀態。

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

圖4 飛行控制仿真姿態曲線

圖5 實際飛行試驗中飛行器位置速度曲線

圖6 實際飛行試驗中飛行器姿態曲線
本文介紹了“孔雀”飛行器垂直起降飛行驗證平臺的數學模型,以微分方程的形式描述了該模型的運動學方程、氣動特性和發動機及伺服機構等執行機構特性,并基于本模型采用一套制導控制律進行了飛行控制半實物仿真。該數學模型可作為被控對象用于運載火箭垂直回收技術等飛行控制技術的仿真研究。