發射載具的重力轉彎與姿態控制

國 立 交 通 大 學

電控工程研究所

碩 士 論 文

發射載具的重力轉彎與姿態控制

Gravity Turn and Attitude Control of Launch Vehicle

研 究 生：李政興

指導教授：林清安 教授

中華民國九十八年九月

發射載具的重力轉彎與姿態控制

Gravity Turn and Attitude Control of Launch Vehicle

研 究 生：李政興 Student：Jeng-Shing Li 指導教授：林清安 教授 Advisor：Prof. Ching-An Lin

國 立 交 通 大 學 電控工程研究所

碩 士 論 文

A Thesis

Submitted to Institute of Electrical and Control Engineering College of Electrical and Computer Engineering

National Chiao Tung University in partial Fulfillment of the Requirements

for the Degree of Master in

Electrical and Control Engineering

September 2009

Hsinchu, Taiwan, Republic of China

中華民國九十八年九月

發射載具的重力轉彎與姿態控制

學生：李政興 指導教授：林清安 教授

國立交通大學 電控工程研究所

中文摘要

發射載具通常是以垂直或接近垂直的角度發射升空，在酬載入軌前除 了高度及速度的需求外還要作適當的轉彎將飛行路徑角控制到接近零 度。載具的轉彎通常利用推力向量控制系統來達成。在大氣層的飛行由於 動態壓力相當大，控制的目標之一是防止轉彎造成的側向力矩太大導致載 具結構損壞。因此載具在大氣層的初期飛行常採所謂重力轉彎的方式，即 保持載具的速度方向與機身一致。本論文探討重力轉彎控制系統的設計。 載具氣動力模式的時變及非線性的特性，使其不方便直接用來作控制系統 設計。本論文提一個建立近似線性模式的方式，並根據線性模式設計以推 力向量為致動的控制系統。垂直發射的載具在進行重力轉彎前必須透過姿 態控制改變俯仰角。本文以幾個數值例子，探討姿態控制及重力轉彎控制 系統的設計。模擬結果顯示，所提出的設計能維持載具的穩定飛行並達到 姿態及重力轉彎的設計目標。 i

Gravity Turn and Attitude Control of Launch Vehicle

Student：Jeng-Shing Li Advisor：Prof. Ching-An Lin

Institute of Electrical and Control Engineering

National Chiao Tung University

Abstract

To put a payload into orbit, a launch vehicle must meet altitude, speed and flight path angle specifications. Usually the flight path angle is near zero, which means the vehicle must turn from its vertical or near vertical launch attitude. The turning of vehicle is usually achieved by a thrust vector control system. Due to the large dynamic pressure during atmospheric flight, one of the control objective during this phase of flight is to prevent the induced bending moment from damaging the vehicle structure. Hence most vehicles adopt the so called gravity turn during the early atmospheric part of flight, that is, to maintain zero angle o attack and zero sideslip angle. This thesis investigates gravity turn control system design. Nonlinear and time-varying nature of aerodynamic model makes it inconvenient for control system design. This thesis proposes a method that construct approximate linear model and use it for control system design. Vertically launched vehicles need attitude control to achieve a certain pitch angle before gravity turn. There numerical examples are used to demonstrate the design attitude control system and gravity turn control system. Simulation results show that the designed control systems yield stable flight and achieve attitude and gravity turn design objectives.

誌 謝

首先我要感謝我的指導教授 林清安博士兩年來細心的

教導以及指引研究的方向，並在學習過程中給於協助以及鼓

勵，除了讓我對於專業領域有更深一層的了解外，也學習到

許多做研究的態度與方法。

另外要感謝實驗室裡的學長與學弟，不僅在研究方面提

供寶貴的經驗以及建議，在生活方面也彼此照顧，他們的鼓

勵給了我很大的助力。

最後要感謝我的家人以及朋友，求學的過程中他們不斷

地支持以及鼓勵，是我精神上最大的支柱。

iii

目錄

中文摘要……….i 英文摘要………ii 誌謝………...iii 目錄………...iv 表目錄………v 圖目錄………...vi 第一章、緒論……….1 第二章、火箭模型………2 2.1 座標系統定義………..2

2.2 火箭運動方程式………..4 2.2.1 座標轉換矩陣………...4 2.2.2 尤拉方程式………...5 2.2.3 四元素方程式………...5 2.2.4 火箭六自由度運動方程式………...6 2.3 氣動力模式………..8 2.4 火箭推力向量控制系統……….11 2.4.1 推力方程式……….12 2.4.2 推力力矩方程式……….12 第三章、斜角發射重力轉彎………14 3.1 重力轉彎方程式………15 3.2 重力轉彎軌跡………...15 3.3 火箭系統轉移函數………...17 3.4 火箭重力轉彎控制器………...20 3.5 火箭的物理性質………...26 3.5.1 火箭介紹………26 3.5.2 火箭的自由飛行軌跡………28 3.6 火箭重力轉彎的軌跡………30 第四章、垂直發射控制………...33 4.1 垂直發射………...33 4.2 火箭姿態控制器………...34 4.3 火箭垂直發射控制的軌跡………...36 第五章、模擬結果與比較………..44 第六章、結論………..48 參考資料………..49 iv

表目錄

表 2.1 座標係定義………2 表 2.2 氣動力係數表………..10 表 3.1 火箭的質量中心位置與轉動慣量………..27 表 5.1 三種不同情況的控制模擬………..44 v

圖目錄

圖 2.1 地心慣性座標………3 圖 2.2 發射座標………3 圖 2.3 地面座標………3 圖 2.4 體座標………3 圖 2.5 火箭的速度分量以及速度之間夾角的定義圖………....8 圖 2.6 火箭推力向量系統………...11 圖 3.1 火箭的自由體圖………..14 圖 3.2 俯仰角的角速度………..16 圖 3.3 俯仰角………..16 圖 3.4 火箭系統方塊圖………..17 圖 3.5 火箭的近似二階轉移函數………..20 圖 3.6 系統的控制方塊圖………..20 圖 3.7 近似的二階函數………..21 圖 3.8a 俯仰角θ………...22 圖 3.8b 角速度q………..22 圖 3.9 時變量………..23 圖 3.10a 俯仰角θ………24 圖 3.10b 角速度q………24 圖 3.11 近似後的控制器方塊圖……….25 圖 3.12 增益值k₂………25 圖 3.13 火箭的外型………26 圖 3.14 火箭的轉動慣量，質心位置，推力分布，質量………27 圖 3.15 馬赫數，推力分布，動態壓力，火箭的加速度………29 圖 3.16 火箭的尤拉角，攻角以及側滑角，角速度，高度對距離……….29 圖 3.17 馬赫數，推力分布，動態壓力，質量……….30 圖 3.18 尤拉角，攻角以及側滑角，角速度，以及高度對距離………..…..31 圖 3.19 推力向量控制系統的角度………32 圖 4.1 火箭的推力分布………..34 圖 4.2 姿態控制器的方塊圖………..35 圖 4.3 增益值k₄……….……….35 圖 4.4 設計的俯仰角與實際模擬的俯仰角………..36 圖 4.5 馬赫數，推力分布，動態壓力，質量……..………37 圖 4.6 尤拉角，攻角以及側滑角，角速度，以及高度對距離………..38 圖 4.7 推力向量控制系統的角度………..39 圖 4.8 設計的俯仰角與實際模擬的俯仰角………..40 圖 4.9 馬赫數，推力分布，動態壓力，質量………..41 vi

圖 4.10 尤拉角，攻角以及側滑角，角速度，高度對距離...42 圖 4.11 推力向量控制系統的角度………..43 圖 5.1 俯仰角，馬赫數，高度對距離，推力分布………..44 圖 5.2 角速度 q，攻角，側滑角，動態壓力………45 圖 5.3 推力向量控制系統角度，質量，加速度………..46 vii

第一章 緒論

載具火箭從發射到太空中，大致上會分成四個階段。一開始幾乎於垂直地面 發射，離開超越大氣層之後，火箭會開始作比較大的轉向，最後會飛行到外太空， 然後繞地球運轉。而本文主要是針對載具火箭從地表發射，到第一節火箭脫離的 這段時間內，共為 63 秒，來進行火箭的控制與軌跡的探討，目的是讓在載具火 箭能夠穩定的升空，並且能讓載具火箭能夠做些許的轉向，以減少載具火箭離開 大氣層後所要轉向的角度。 載具火箭從地面上發射最後進入軌道繞地球運行，其中飛行的過程最常採用 的方法就是重力轉彎；重力轉彎的飛行方法會讓火箭的機身幾乎與速度方向重 合，也是就攻角等於零度，如此一來，在氣動力影響相當大的氣層下飛行，會明 顯降低氣動力對火箭的扭力，增加飛行時的穩定度。然而飛行最終的目的是進入 繞地球運轉的軌道，也就是俯仰角要等於九十度，但並不是所有的重力轉彎的軌 跡都可以達到這個目的，必須根據火箭本身的推力與重量比、燃料的燃燒時間等 計算出火箭所該具備的速度以及初始俯仰角，才能使火箭用重力轉彎的方法進入 軌道中[8]。然而計算出來後的初始俯仰角並不能從現實情況中實現，也就是發 射角度過大；另一方面重力轉彎的軌跡並不是最佳化的軌跡[6]，假使重力轉彎 並不能達到所要求的軌跡，應該如何修正？甚至在氣動力相當大的氣層中，火箭 是否真的能順利的做重力轉彎？本文在第二章中介紹火箭的模型，火箭的運動方 程式，氣動力模型，火箭推力向量控制系統，第三章中介紹斜角發射重力轉彎的 方法與重力轉彎控制器的設計，第四章介紹火箭垂直發射與姿態控制器的設計， 第五章為模擬結果以及討論，第六章為結論。 1

第二章 火箭模型

在本章中，我們將介紹火箭發射所需要用到的座標系以及方程式。2.1 節介 紹座標系的定義，2.2 節介紹火箭運動方程式，2.3 節介紹氣動力模式，2.4 節介 紹火箭推力向量控制系統。 2.1 座標係統定義[1] 為了方便描述火箭在飛行中的動態以及姿態，定義三個座標系統如表 2.1： 座標名稱 原點 座標定義 地心慣性座標 （S_I） 地心位置 X_I :指向格林威治子午線以及赤道面 : I Z 指向北極 : I Y 符合右手座標系統 發射座標 （S_L） 發射點的位 置 L X ：發射點所在的地表平面垂直向上 L Z ：發射點所在的地表平面指向要發射的方向 L Y ：符合右手座標系統 地面座標 （S_G） 火箭所在的 地平面位置 G X ：火箭所在的地表平面指向北方 G Y ：火箭所在的地表平面指向東方 G Z ：火箭所在的地表平面指向地下 體座標 （S_B） 火箭的質心 位置 B X ：為火箭的長軸方向，指向前方 B Y ：垂直XB在火箭的左右平面上，指向右方 B Z ：符合右手座標系統 表 2.1 座標係定義 地心慣性座標（S_I ） 地心慣性座標如圖 2.1 所示，是為了方便描述火箭在空間中的飛行軌跡，像 是高度，方位，飛行距離等。 發射座標（S_L） 發射座標作為慣性座標如圖 2.2 所示，是為了方便描述火箭與發射點位置之 間的關係並藉此定義出尤拉角。 2

地面座標（S_G） 又叫做 NED 座標係，如圖 2.3 所示，其座標係的Z 為指向地下，而不是指_G 向地心是因為所使用的地球模型並不一定為球體。 體座標（S_B） 體座標，如圖 2.4 所示是為了方邊描述火箭的的姿態，例如尤拉角以及攻角 及侧滑角。 圖 2.1 地心慣性座標 圖 2.1 地心慣性座標 圖 2.2 發射座標 圖 2.4 體座標 圖 2.3 地面座標 3

2.2 火箭運動方程式 2.2.1 座標轉換矩陣 尤拉角(Euler angle)大多用來描述火箭的姿態，為體座標對於發射座標的相 對角度，以( , , )φ ψ θ 表示。依序旋轉發射座標的三個座標軸，使得發射座標平行 於體座標，而為了方便推導旋轉過程，我們假設兩個座標系的原點是一致的；旋 轉的過程如下： 1 繞X 軸旋轉_L φ 角，使得旋轉後的Y ′ 投影在_L Y_B−Z_B的平面上為 ，此時經一次 旋轉後的座標軸為 B Y L L L X′ −Y′− ′ 。 Z 2 繞Z′ 軸旋轉_L ψ 角，使得旋轉後的Y ′′ 與_L 重合，此時經過兩次旋轉後的座標軸 為 B Y L L L X′′−Y′′−Z′′。 3 繞Y ′′_L 軸旋轉θ角，使得旋轉後的X ′′′ 與_L X 重合，_B Z ′′′ 與_L Z 重合，則此時經過三_B 次旋轉後的發射座標軸X_L′′′− −Y_L′′′ Z_L′′′將平行於體座標軸。 從上述的旋轉關係中，可以依序推得每個旋轉中發射座標跟體座標的轉換關 係，如(2.1)-(2.3)式所示： ₁ 1 0 0 0 cos sin 0 sin cos L L L L L L L L L X X X Y Y T Y Z Z Z φ φ φ φ ′ ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ ⎢ _{′ =}⎥ ⎢ ⎥ ⎢ ⎥₌ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ′ ⎢ ⎥ ⎢ − ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ ⎣ ⎦ ⎣ ⎦ ⎣ ⎦ (2.1) ₂ cos sin 0 sin cos 0 0 0 1 L L L L L L L L L X X X Y Y T Y Z Z Z ψ ψ ψ ψ ′′ ′ ′ ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ ⎢ _′′⎥ ⎢_{= −} ⎥ ⎢ _′⎥₌ ⎢ _′⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ′′ ′ ′ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ ⎣ ⎦ ⎣ ⎦ ⎣ ⎦ (2.2) ₃ cos 0 sin 0 1 0 sin 0 cos L B L L B L L B L L L L X X X Y Y Y T Y X Z Z Z θ θ θ θ ′′′ − ′′ ′′ ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ ⎢ _′′′⎥ ⎢₌ ⎥ ⎢₌ ⎥ ⎢ _′′⎥₌ ⎢ _′′ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ′′′ ′′ ′′ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ ⎣ ⎦ ⎣ ⎦ ⎣ ⎦ ⎣Z ⎦ ⎥ ⎥ L L L (2.3) 因此我們可以得到發射座標跟體座標之間的轉換矩陣為 3 2 1 B L B L BL B L X X X Y T T T Y T Y Z Z Z ⎡ ⎤ ⎡ ⎤ ⎡ ⎢ ⎥₌ ⎢ ⎥₌ ⎢ ⎢ ⎥ ⎢ ⎥ ⎢ ⎢ ⎥ ⎢ ⎥ ⎢ ⎣ ⎦ ⎣ ⎦ ⎣ ⎤ ⎥ ⎥ ⎥⎦ (2.4) 4

2.2.2 尤拉方程式 假設角速度為 5 B ω= pX_B+qY_B+rZ (2.5) 根據 2.2.1 所推導出發射座標跟體座標之間的轉換關係，將角速度改寫成 L B L X Y Z ω φ= ′ +θ +ψ ′′ (2.6) 其中 sin cos L B B Z′′ = − θX + θZ (2.7)

cos sin cos cos cos sin sin

L L L B B B X′ = ψX′′− ψY′′= ψ θX − ψ θZ − ψY B B (2.8) 代入整理後可得

(cos cos cos sin sin ) ( sin cos ) ( cos cos sin ) ( sin ) ( cos sin cos )

L B L B B B B B B B X Y Z X Z Y Y X Z X Y Z ω φ θ ψ φ ψ θ ψ θ ψ θ ψ θ θ φ ψ θ ψ θ φ ψ θ φ ψ θ ψ θ ′ ′′ = + + = − − + + − + = − + − + + − + 與(2.5)式比較係數後，可以得到

cos cos sin _{cos cos 0 sin}

sin sin 1 0

cos sin 0 cos cos sin cos

p _p q q r r φ ψ θ ψ θ _{ψ θ θ φ} φ ψ θ ψ θ ψ θ θ φ ψ θ ψ θ ⎧ = − ⎡ ⎤ ⎡ − ⎤⎡ ⎤ ⎪ _{⎢ ⎥ ⎢ ⎥⎢ ⎥} = − + ⇒ = − ⎨ _{⎢ ⎥ ⎢ ⎥⎢ ⎥} ⎪ _{= − + ⎢ ⎥ ⎢ ⎥⎢ ⎥} ⎣ ⎦ ⎣ _{⎦ ⎣ ⎦} ⎩ ψ ⎡ ⎤ ⎥ ⎢ ⎥ ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ (2.9) 經過整理，可以得到

sec cos 0 sec sin

tan cos 1 tan sin

sin 0 cos p q r φ ψ θ ψ θ θ ψ θ ψ θ ψ θ θ ⎡ ⎤ ⎡ ⎤ ⎢ ⎥ ⎢₌ ⎢ ⎥ ⎢ ⎢ ⎥ ⎢_⎣ − ⎥_⎦ ⎣ ⎦ (2.10) 上式即為尤拉方程式。 2.2.3 四元素方程式[10] 我們可以用四個係數 、 、 、 來定義尤拉角e₀ e₁ e₂ e₃ φ、ψ 、φ，如(2.11)至(2.14) 式所示：

0 cos cos cos sin sin sin 2 2 2 2 2 e 2 ψ θ φ ψ θ = + φ (2.11)

1 cos cos sin sin sin cos

2 2 2 2 2

e

2

ψ θ φ ψ θ

= − φ (2.12)

2 cos sin cos sin cos sin

2 2 2 2 2

e

2

ψ θ φ ψ θ

= − φ (2.13)

3 sin cos cos cos sin sin

2 2 2 2 2 e 2 ψ θ φ ψ θ = + φ (2.14) 透過四元素與三個尤拉角之間的關係，我們可以用一個四元素方程式來表示 四元素與角速度之間的關係，如(2.15)所示 0 1 2 3 1 0 3 2 2 3 0 1 3 2 1 0 1 2 e e e e p e e e e q e e e e r e e e e − − − ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ ⎢ ⎥ ⎢ ₋ ⎥ ⎢ ⎥ ⎢ ⎥₌ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ − ⎥ ⎢ ⎥ ⎢ ⎥ ⎢₋ ⎥ ⎣ ⎦ ⎣ ⎦ ⎣ ⎦ (2.15) (2.15)式中四元素的初始值可以從(2.11)式至(2.14)式中計算尤拉角的初始值 可得，而尤拉角可以透過四元素計算出來如(2.16)式所示： 1 2 3 0 1 2 2 2 2 0 1 2 3 1 0 3 1 2 1 1 3 0 2 2 2 2 2 0 1 2 3 2( ) tan sin ( ) 2( ) tan e e e e e e e e e e e e e e e e e e e e φ ψ θ − − − ⎛ + ⎞ = _{⎜ ⎟} − − + ⎝ ⎠ = − ⎛ + ⎞ = _{⎜ ⎟} + − − ⎝ ⎠ (2.16) 我們可以從(2.10)式中看到角速度與尤拉角的一次微分中之間的關係存在著 secψ 函數。在本文中在火箭的偏航角方面我們希望控制在零，若在模擬中存在 著secψ 函數，可能會出現錯誤的結果。因此模擬過程中，我們將用四元素方程 式來計算火箭角速度與四元素之間的關係，再利用(2.16)式計算出火箭的尤拉角。 2.2.4 火箭六自由度運動方程式 火箭的六自由度運動方程式包含了平移運動方程式(translational equation)以 及旋轉運動方程式(rotational equation)。 6

平移運動方程式[1] 描述物體質心運動的位置、速度以及加速度的改變；由於 座標為慣性座 標軸，因此在 座標三個座標軸上各別的加速度 ，速度 ，以及位移量 I S I S a_I v_I x 之_I 間的關係就如同兩個積分器： xI x yI yI v a v a v a = = = I zI zI I I I zI (2.17a) 以及 I x I y x v y v z v = = = (2.17b) 如果可以得知火箭在 座標下所受的總力 ，包括火箭的推力，氣動力以 及重力等，根據牛頓第二運動定律 I S F_I I F =ma_I，再透過上述的平移運動方程式，可 以計算出火箭在飛行中的位置，速度以及加速度。 旋轉運動方程式[1] 假設火箭的質心在幾何中心，且火箭的外型以及質量對稱於S_B座標系的xy 平面以及 xz 平面，因此Ixy =Iyz =Ixz =0，則火箭的力矩跟角速度的關係為: ( ) ( ) ( ) y z x x x y z x y y x y z z z qr I I M p I I M rp I I q I I pq I I M r I I − = + − = + − = + (2.18) 其中M 、_x M_y、M 為在z SB座標下的力矩，p、 、q r為在SB座標下的角速度， x I 、I_y、I 為火箭的在z SB座標下的轉動慣量。 若能知道推力所產生的力矩，再加上氣動力所產生的力矩，根據旋轉運動方 程式，可以以計算出在 座標下的角速度，再經由尤拉方程式便可以算出火箭 在飛行的姿態，也就是三個尤拉角。 B S 7

2.3 氣動力模型[1] 在火箭飛行的過程中，為了計算氣動力，我們得先知道馬赫數、攻角以及風 向角，而這些角的定義如圖 2.5 所示。 圖 2.5 火箭的速度分量以及速度之間夾角的定義圖 V 為火箭的速度方向，u 為火箭在 x 方向的速度分量， v 為在y方向的速度 分量， w 為在 方向的速度分量。 z 攻角定義為 1 tan w u α ₌ − ⎛ ⎞ ⎜ ⎟ ⎝ ⎠ (2.19) 側滑角定義為 tan 1 v u β ₌ − ⎛ ⎞ ⎜ ⎟ ⎝ ⎠ (2.20) 全攻角定義為 2 2 1 tan c w v u α ₌ − ⎛ + ⎞ ⎜⎜ ⎝ ⎠⎟⎟ (2.21) 8

氣動力所產生的力跟力矩方程式如(2.22)-(2.23)所示：

2 2

1 1 1

2 2 2

AxB fx AyB fy AzB fz

2 F = ρV C S F = ρV C S F = ρV C S (2.22) 2 2 1 1 1 2 2 2

AxB mx AyB my AzB

2 mz

M = ρV C SL M = ρV C SL M = ρV C SL (2.23)

其中L為參考長度(reference length)， 為參考面積(reference area)，S ρ 為大氣密 度，V 為火箭速度，而方程式中的各個係數C_fx、C_fy、C_fz、 、 、 定 義如(2.24)以及(2.25)所示 mx C C_my Cmz ( , ) ( , ) ( , ) fx A fy Y fz N C C M C C M C C M c α β α = − = − = − (2.24)

[

[

]

0 ( , ) ( , ) ( , ) ( , ) mx my m CG N mz n CG Y C C C M X C M C C M X C M

]

α α β β = = + = + (2.25) (2.24)為力的部分，(2.25)為力矩的部分，其中M 為馬赫數，X_CG為火箭的質 心位置到火箭頂端的距離；C_A、C_Y、C_N、C_m、C_n則由表 2.2 作線性計算可得。 所以要計算氣動力，必須要知道攻角、側滑角以及馬赫數，經由表 2.2 作線 性計算求出 、 、 、 、 等係數，再代入(2.24)及(2.25)式中計算氣動 力的力跟力矩的六個係數，最後再由(2.22)及(2.23)式計算出氣動力所產生的力跟 力矩。 A C C_Y C_N C_m C_n 9

C

N M

α

0.5 0.8 1.2 1.5 2.0 3.0 4.0 0 0 0 0 0 0 0 0 2 0.736 0.791 0.811 0.774 0.726 0.607 0.548 5 1.906 2.141 2.327 2.223 2.089 1.710 1.494

Cm

M

α

0.5 0.8 1.2 1.5 2.0 3.0 4.0 0 0 0 0 0 0 0 0 2 -8.203 -8.962 -9.149 -8.626 -8.090 -6.539 -5.809 5 -21.403 -24.420 -26.787 -25.280 -23.682 -18.742 -15.978

C

A M

α

c 0.5 0.8 1.2 1.5 2.0 3.0 4.0 0 0.769 0.835 2.063 1.994 1.710 1.233 1.010 2 0.769 0.835 2.063 1.994 1.710 1.233 1.010 5 0.896 0.877 2.198 2.097 1.730 1.241 1.010

C

Y M

β

0.5 0.8 1.2 1.5 2.0 3.0 4.0 0 0 0 0 0 0 0 0 2 -0.559 -0.611 -0.713 -0.558 -0.414 -0.390 -0.305 5 -1.394 -1.539 -1.838 -1.420 -1.055 -0.992 -0.855

Cn

M

β

0.5 0.8 1.2 1.5 2.0 3.0 4.0 0 0 0 0 0 0 0 0 2 6.885 7.647 9.013 6.772 4.646 4.223 2.996 5 17.167 19.242 23.246 17.168 11.755 10.580 8.508 表 2.2 氣動力係數表 10

2.4 火箭推力向量控制系統 火箭推力控制系統所採用的座標係為S_B座標係，其定義如圖 2.6 所示： F₁ F1y F_1x F_1z δ_y1 δ_z1 F₂ F_2x δ_y2 F₃ F3y F_3x F_3z δ_y3 W W L F_2z

X

Y X₁ X₂ X₃ Z₁ Z₃ CG F3y δ_z2 δ_z3

L

Δ

Z₂ Z F₁ F1y F_1x F_1z δ_y1 δ_z1 F₂ F_2x δ_y2 F₃ F3y F_3x F_3z δ_y3 W W L F_2z

X

Y X₁ X₂ X₃ Z₁ Z₃ CG F3y δ_z2 δ_z3

L

Δ

Z₂ Z 圖 2.6 火箭推力向量控制系統 : CG 為火箭發射前以及發射後在 X 軸上的質心位置。 : L Δ 火箭發射後的質心位置與火箭發射前的質心位置在 X 軸方向上的距離。 : L 推力器的頂部與質心在 X 軸方向上的距離。 : W 兩翼的推力器與中間推力器的距離。 : i F 推力器的推力向量，可以被分解成三個方向的向量F_ix、F_iy、F_iz，i= ∼1 3。 : yi δ 為F_i在X − 平面上的投影跟Y X 之間的夾角，_i i= ∼1 3。 : zi δ 為F_i在X − 平面上的投影跟Z X 之間的夾角，_i i= ∼1 3。 11

2.4.1 推力方程式

根據上述的圖示與定義，推力器的推力向量( )F₁ 與其分量之間的關係為

1 1 2 2 2

1 1 1

1 1

tan y (2.26a) tan z (2.26b)

y y x x F F F F F F F F δ = δ = = 2 1x+ 1y+ 1z (2.26c)

將(2.20a)改寫成F_1y =F_1xtanδ_y1 (2.20b)改寫成F_1z =F_1xtanδ_z1 代入(2.20c)可得 2 2 2 2 2 2 1 1x 1xtan y1+ tan1x 1 F =F +F δ F δ_z (2.27) 經整理可以得到 1 1 2 2 1 1 1 1 tan tan x y z F F δ δ = + + (2.28) 將(2.22)式代回(2.20a)跟(2.20b)就可以得到 1 1 2 2 1 1 1 tan 1 tan tan z z y z F F δ₁ δ δ = + + (2.29) 1 1 2 2 1 1 1 tan 1 tan tan y y y z F F δ ₁ δ δ = + + (2.30) 接著令 2 2 1 for 1 ~ 3 1 tan tan i yi zi i δ δ Δ = = + + (2.31) 因此推力器的推力方程式為 1 1 2 2 3 3 1 1 1 2 2 2 3 3 3 1 1 1 2 2 2 3 3 3

tan tan tan tan tan tan

TxB TyB y y y TzB z z z F F F F F F F F F F F F δ δ δ δ δ δ = Δ + Δ + Δ = − Δ − Δ − Δ = − Δ − Δ − Δ (2.32) 12

2.4.2 推力力矩方程式 根據圖 2.6 所示與定義，我們可以列出推力器產生的力矩關係式(2.33)-(2.35)： M_TxB =F W_1z −F W_3z (2.33) 1 ( ) 2 ( ) 3 ( TyB z z z ) M = −F L+ Δ −L F L+ Δ −L F L+ Δ (2.34) L 1 3 1 ( ) 2 ( ) 3 ( TzB x x y y y ) M =F W−F W +F L+ Δ +L F L+ Δ +L F L+ Δ (2.35) L 將之前所推導出的(2.28)-(2.30)式代入(2.33)-(2.35)中，並用(2.31)式的定義，經整 理後，可以將(2.33)-(2.35)改寫成： M_TxB = ΔF_{1 1}tanδ_z1W− ΔF_{3 3}tanδ_z3W (2.36)

1 1tan 1( ) 2 2tan 2( ) 3 3tan 3( )

TyB z z z x

M = − ΔF δ L+ Δ − ΔL F δ L+ Δ − ΔL F δ L+ ΔL (2.37)

1 1 3 3 1 1tan 1( ) 2 2tan 2( ) 3 3tan 3( )

TzB y y y

M = ΔF W− ΔF W + ΔF δ L+ Δ + ΔL F δ L+ Δ + ΔL F δ L+ ΔL (2.38)

第三章 斜角發射重力轉彎

火箭從地面升空之後，在平流層(0~11km)以及對流層(11~50km)之中飛行， 由於在平流層裡的空氣密度很大，加上火箭本身的速度，會使得動態壓力 (dynamic pressure)非常的大;而在對流層中飛行，雖然空氣密度沒有平流層來的 大，但火箭仍然持續的在加速中，動態壓力依舊相當的大，不過已經有隨著空氣 密度迅速的減少而使得動態壓力逐漸下降的趨勢。在動態壓力這麼大的環境下飛 行，只要產生些許的攻角，就會在體座標的 方向產生相當大的力矩，造成火 箭控制上的困難。因此在空氣密度大的氣層中，要讓火箭能夠穩定的升空，就必 須控制攻角維持在幾乎為零的角度。然而就實際層面來看，在火箭飛行的過程 中，所能量測得到的資訊有限，而攻角就是其中一項相當難以量測的數據；因此 要控制攻角為零，必須藉由控制其他的參數來達成，而在這裡，我們採用了重力 轉彎。 B Y 重力轉彎的方法會讓速度方向盡量與火箭的體座標的X 軸本身重合，也就 是會使得攻角趨近於零，一旦攻角趨近於零，氣動力的影響就會降低很多；而重 力的影響在垂直部分，會使得火箭作緩慢的轉向。 mg , T V θ cos mg θ sin mg θ D Local horizon To center of Earth mg , T V θ cos mg θ sin mg θ D Local horizon To center of Earth 圖 3.1 火箭的自由體圖 14

3.1 重力轉彎方程式[7] 15 ) 從圖3.1 中我們可以將作用在火箭上的重力 ( 分解成垂直於火箭體座標 中的 mg B X 軸

(

mgsinθ 與平行於體座標中的

)

X 軸_B

(

mgcosθ ；此外，T 為火箭的推

)

力，D為作用在火箭體座標中X 軸的空氣阻力，_B V 為火箭的速度向量，而θ 為 尤拉角中的俯仰角。假設速度向量和X 軸重合，亦即攻角及側滑角均為零，則_B 我們可以將作用在火箭上的力整理後可以得到下面的式子： cos dV T D mg m dt θ − − = (3.1) sin d mg mV dt θ θ = (3.2) (3.1)式為平行於火箭體座標中X 軸的力方程式；(3.2)式為垂直於火箭體座標中_B B X 軸的力方程式。從(3.2)式中可以看出若火箭一開始的發射角度並不是垂直發 射，即θ ≠0，則在發射之後由於重力的影響，會使火箭產生一個角速度讓火箭 對體座標中的Y_B軸旋轉，就會使得俯仰角 (θ)越來越大。因此我們控制火箭使其 作重力轉彎的方法，便是讓火箭去跟隨重力產生的角速度，將(3.2)改寫成 sin d g dt V θ θ = (3.3) 3.2 重力轉彎軌跡 從(3.3)式中可以看出俯仰角以及速度跟角速度之間的關係。火箭的速度越 大，則重力轉彎的角速度就越小，則火箭所能轉的俯仰角也就越小。雖然我們可 以藉由降低火箭的速度來增加重力轉彎所轉向的角度，但是效果並不明顯。因此 我們在本小節討論火箭一開始的發射角與重力轉彎軌跡之間的關係。我們將計算 出來火箭的速度代入(3.3)式，並將(3.3)中的初始俯仰角以及角速度設為零，根據 (3.3)式的關係，計算火箭一開始在不同的發射角發射後，飛行 32 秒之內的俯仰 角以及角速度，如圖 3.2 以及圖 3.3 所示。其一開始的發射角度分別為 ， ， ， ， 。 2 − −4 6 − −8 −10

0 5 10 15 20 25 30 35 40 -4.5 -4 -3.5 -3 -2.5 -2 -1.5 -1 -0.5 0 time (s) pi tc h an gu lar r a te ( d eg re e/ s ) -2o -4o -6o -8o -10o 圖 3.2 俯仰角的角速度 0 5 10 15 20 25 30 35 40 -35 -30 -25 -20 -15 -10 -5 0 time (s) pi tc h an gl e (deg re e) -2o -4o -6o -8o -10o 圖 3.3 俯仰角 16

圖 3.2 中角速度以及圖 3.3 中的俯仰角在前面的 1.4 秒內皆為零是因為我們 假設火箭在 0 秒之後到 1.4 秒之前是從發射架上發射，直到 1.4 秒只後才離開發 射架，因此在這 1.4 秒期間火箭所受到的重力影響由發射架吸收。從圖 3.2 中可 以看出初始發射角度越大，角速度也越大，這跟(3.3)式的關係一致。而從圖 3.3 中可以看出在越大的初始角度下發射，經由重力轉彎後所產生的俯仰角也越大。 3.3 火箭系統轉移函數 在設計控制器之前，我們必須先了解火箭的系統轉移函數；火箭轉移函數因 為氣動力的部分，使得整個轉移函數變得較為複雜。為了要計算氣動力所產生的 力矩，從(2.25)式中，我們必須根據馬赫數以及攻角來決定係數，再代入(2.23) 式中算出力矩。如圖 3.4 所示，計算出氣動力力矩後，再與推力向量控制系統中 所造成的力矩相加，即為火箭所受到的總力矩，再經由火箭六自由度運動方程式 便可以算出火箭的尤拉角以及角速度。我們曾經嘗試用了一些假設條件來將氣動 力的部分作線性化的處理，來求得火箭的系統轉移函數；然而求得的系統函數為 一個四階的函數，且此四階函數中的係數大多為時變函數，使得控制器在設計上 會有很大的難度，此外火箭在實際模擬的結果也並未完全符合所假設的條件。因 此，我們試著用另一個方法來求得火箭的系統轉移函數。 5 5 s+ eqC

δ

FL eq

δ

TB M Aero dynamics C

α α β

Mach number AB M

φ θ ϕ

6DOF p q r input output1 output2 5 5 s+ eqC

δ

FL eq

δ

TB M Aero dynamics C

α α β

Mach number AB M

φ θ ϕ

6DOF p q r input output1 output2 圖 3.4 火箭系統方塊圖 17

在求系統函數之前，我們先做一些簡單的線性化處理，使得圖 3.4 中的系統 可以依照三個尤拉角而分成三個獨立的系統，以避免有多輸入多輸出的重疊問 題。 假設火箭的尤拉角都趨近於零，那麼我們可以將尤拉方程式(2.10)簡化為: p q φ = θ = ψ = (3.4) r 假設火箭的角速度都趨近於零，則我們可以將力矩方程式(2.18)式簡化為: y x z x y z M M M p q r I I I = = = (3.5) 則原本的六自由度非線性方程式可以簡化成三組互相獨立的雙積分線性方 程組： x x p M p I φ = = (3.6) y y q M q I θ = = (3.7) z z r M r I ψ = = (3.8) 假設δ_yi以及δ_zi很小，則 tanδ_yi ≈ ， tan0 δ_zi ≈ 以及(2.31)式可以近似為0 ，我們可以將推力向量控制系統的力矩方程式(2.36)至(2.38)式改寫成 1 i Δ ≈ 1 3 1 2 3 1 2 3 ( ) ( ( ) TxB z z TyB z z z TzB y y y M FW M FL M FL ) δ δ δ δ δ δ δ δ = − = − − − = + + (3.9) 定義三個方向的等效角 來等效推力向量控制系統中的六 個角度，使得(3.9)式可以改寫成： _ _ _ [δp eq δq eq δr eq]T _ _ _ TxB p eq TyB q eq TzB r eq M FW M FL M FL δ δ δ = = = (3.10) 18

則三個等效角與推力向量控制系統中的六個角度之間的關係為: 1 2 _ 3 _ 1 _ 2 3 y y p eq y q eq z r eq z z D δ δ δ δ δ δ δ δ δ ⎡ ⎤ ⎢ ⎥ ⎢ ⎥ ⎡ ⎤ ⎢ ⎥ ⎢ _{⎥ = ⎢ ⎥} ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ _{⎢ ⎥} ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ (3.11) 其中 0 0 0 1 0 1 0 0 0 1 1 1 1 1 1 0 0 0 D − ⎡ ⎤ ⎢ ⎥ =_⎢ − − − _⎥ ⎢ ⎥ ⎣ ⎦ 雖然兩者之間的轉換矩陣並不是方陣，不過可以觀察出秩等於三，因此兩者 之間的逆轉換矩陣可以寫成： 1 2 _ 3 1 _ 1 _ 2 3 ( ) y y p eq y T T q eq z r eq z z D DD δ δ δ δ δ δ δ δ δ − ⎡ ⎤ ⎢ ⎥ ⎢ ⎥ _{⎡ ⎤} ⎢ ⎥ _{⎢ ⎥} = ⎢ ⎥ _{⎢ ⎥} ⎢ ⎥ _{⎢ ⎥} ⎣ ⎦ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ (3.12) 經由上述的假設以及近似，我們可以將整個系統依照三個尤拉角來分割成三 個獨立的系統，這樣可以將三個尤拉角的控制看作是單輸入單輸出的系統來設 計。 我們假設從推力向量控制系統的角度到火箭的俯仰角之間的系統轉移函數 可以用一個二階系統來近似；然而整個火箭系統為一個時變函數，因此我們希望 用一個標準二階函數來近似系統的暫態響應部分，以及用一個時變函數 來近似 整個系統的時變量。以俯仰角為例，我們嘗試用一個二階函數來近似推力向量控 制系統 Z 軸方向的角度至火箭俯仰角之間的系統轉移函數，也就是近似圖 3.4 中 的 input 至 output1 之間的轉移函數，如圖 3.5 所示。 k 19

input 2 2 2

2

n n n

kw

s

+

ζ

w s

+

w

_ q eqC

δ

output1

θ

input 2 2 2

2

n n n

kw

s

+

ζ

w s

+

w

_ q eqC

δ

output1

θ

圖 3.5 火箭的近似二階轉移函數 因此，我們在輸入的部分輸入適當的推力向量控制系統的角度，然後從整個 系統的模擬中去觀察俯仰角的結果。從俯仰角的暫態響應裡我們可以用最大超越 量以及安定時間來計算近似的二階轉移函數中的ζ 以及 ；然而整個系統為時 變系統，所以模擬出來俯仰角的穩態部分也不會是定值，因此我們在近似的二階 轉移函數中分子部分多乘上了一個時變函數k，用 來近似俯仰角穩態部分的變 化量。因此所近似出來的二階系統函數，除了可以減少控制器的階數，也因為系 統近似到只有一個時變的係數k，也降低了設計控制器的複雜度。 n w k 3.4 火箭重力轉彎控制器 根據 3.3 節的近似，再考慮氣動力的影響，以俯仰角為例，控制方塊圖如圖 3.6 所示[3]。 input 5 5 s+ 1 s 2 k 1 k q q eqC δ _{− θ} q 1 s output FL 1 y I C θ q eq δ ₋ TyB M q q k′ 1 s Aero dynamics α Mach number AyB M input 5 5 s+ 1 s 2 k 1 k q q eqC δ _{− θ} q 1 s output FL 1 y I C θ q eq δ ₋ TyB M q q k′ 1 s Aero dynamics α Mach number AyB M 圖 3.6 系統的控制方塊圖 20

圖 3.6 中由於控制的目的是跟隨重力轉彎的角速度，因此在設計控制器時採 用角速度當作回授，並在外邊的迴路上增加了一個積分器。然而 3.3 節的方法為 近似推力向量系統的角度至俯仰角之間的系統函數，而重力轉彎控制器是控制角 速度，因此我們將近似的二階轉移函數改寫成如圖 3.7 所示： input 2 2 2

2

n n n

kw s

s

+

ζ

w s

+

w

_ q eqC

δ

output1

q

₁

s

θ

( )

G s

input 2 2 2

2

n n n

kw s

s

+

ζ

w s

+

w

_ q eqC

δ

output1

q

₁

s

θ

( )

G s

圖 3.7 近似的二階函數 從圖 3.7 中可以看出，我們將近似的二階系統再加上了個微分器，然後經由 一個積分器使輸出為俯仰角。我們近似的重點跟 3.3 節一樣是放在推力向量系統 的角度與俯仰角之間的關係，只是將近似出來的結果再加上一個微分器來近似成 推力向量控制系統的角度到角速度 q 之間的系統轉移函數。 21 2 2 2 至於輸入部分，我們可以藉由力矩方面來決定推力向量控制系統的角度為 何，才能得到合理的模擬結果來近似出我們所設計的二階系統。根據 2.3 節的氣 動力模式，假設攻角為 0.1 度，馬赫數為 2，以及動態壓力為 ，我們 可以算出氣動力在 方向所造成的力矩為 ；而推力向量系統在 Z 軸方向轉 0.1 度所造成的力矩約為 。因此我們讓推力向量系統的在 Y 軸方向的角度設為 0，Z 軸方向設為 0.1 度，將飛行時的發射角設為 度，來 模擬出整個系統中的 q 跟 4 15 10 N/m× B Y −1.3 10 N/m× 4 4 0.6 10 N/m× 2 − θ，並藉由 q 跟θ的響應圖中的資訊找出安定時間以及 最大超越量，來近似出標準二階系統，其響應圖如圖 3.8a 以及圖 3.8b 所示：

0 5 10 15 20 25 30 -18 -16 -14 -12 -10 -8 -6 -4 -2 pi tc h an gl e (deg re e) time (s) 圖 3.8a 俯仰角θ 0 5 10 15 20 25 30 -0.08 -0.06 -0.04 -0.02 0 0.02 0.04 0.06 0.08 q ( ra d / s ) time (s) 圖 3.8b 角速度 q 22

從圖 3.8a 中可以觀察出θ的穩態不是一個定値，幾乎是一條逐漸遞減的斜直 線，而圖 3.8b 中 q 的穩態幾乎為零。我們根據圖 3.8a 中的安定時間以及最大超 越量來近似圖 3.7 中G s( )裡的二階函數，而G s( )中的 值為近似圖 3.8a 中k θ的穩 態中逐漸遞減的斜直線。 由模擬出來的響應圖中，我們近似系統中安定時間為 15 秒以及最大超越量 為 70%，而根據這兩項數據可以算出： 0.1128 w_n 1.8911 δ = = 因此近似出來的二階系統函數為 2 3.5763 ( ) 0.4267 3.5763 ks G s s s = + + (3.13) 其系統時變量k值就如圖 3.9 所示。 0 5 10 15 20 25 30 35 2.5 3 3.5 4 4.5 5 5.5 time (s) k 圖 3.9 時變量 k 圖 3.10a 以及圖 3.10b 則是近似跟實際模擬的比較圖，紅色的線為實際模擬 的圖，藍色的線為近似的二階系統。 23

0 5 10 15 20 25 30 -18 -16 -14 -12 -10 -8 -6 -4 -2 pi tc h ang le ( de gr ee) time (s) simulation result approximate system 圖 3.10a 俯仰角θ 0 5 10 15 20 25 30 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 q ( rad / s ) time (s) simulation result approximate system 圖 3.10b 角速度 q 24

現在我們可以將整個系統的方塊圖簡化如圖 3.11 所示： 2 k 1 k δz1~ 3z 1 θ s C θ ₁ q s k1 k2

G s

( )

1~ 3 z z δ _{1 θ} s C θ ₁ q s

G s

( )

25 圖 3.11 近似後的控制器方塊圖 圖 3.11 中的G s( )為(3.13)式，則閉迴路轉移函數可以寫成： 1 2 2 2 1 2 (0.4627 3.5763 ) 3.5763 3.5763 c s kk s kk k θ = + + + + (3.14) 從圖 3.11 以及(3.14)式中，我們可以將控制器中 的值設為 ，則利用根 ，來設計 以及 值，就可以設計出控制器中的 兩個增益值 以及 。這樣設計的好處是控制器中的兩個增益值只會有一個是隨 時間變化的參數， 為 ；另一個增益值 就可以用定值來決定。也由於降低了 控制器設計的難度，也就不用再設計圖 中的 3.5763kk k q 2 k n k/ 軌跡法則決定極點與零點的位置 k₁ n 1 k k₂ 即 k₂ k₁ 3.6 k′值。經設計調整後得 ， 為一個隨時間變化的參數，如圖 3.12 1 0.5 k = 所示： 2 k 0 5 10 15 20 25 30 35 0 0.2 0.4 0.6 0.8 1 1.2 1.4 time (s) k2 2 k 圖 3.12 增益值

26 控制器設計完後，輸入方面在俯仰角的部分我們則依照(3.3)式輸入重力轉彎 .5 火箭的物理性質 火箭介紹 的火箭模型分為四節，如圖 3.13 所示。在第一節部份有三個推 火箭第一節有三個推進器，且每個推進器上的燃料為 3648 公斤，因此第一 節的 零，此時 的角速度；至於滾轉角以及偏航角的部分，其控制器也用上述的方法設計；而所 給的輸入方面則都給為零。 3 3.5.1 我們所用到 進器，在第二、三、四節則各為一個推進器。火箭的全長為 15 公尺，總重量包 括燃料質量為 22 公噸，酬載為 50 公斤。 圖 3.13 15 11 6 3 1.7 0.6 0.6 1 1 1 50 (kg) 0.5 2.5 8 22 燃料總重量為 10944 公斤。在飛行的過程中，我們讓燃料的燃燒速率維持為 定值，其火箭的總重量與時間的關係圖如圖 3.14 示。推力分布為從 0 秒的時候 開始燃燒，在 1 秒的時候達到最大值為 5 8.65 10× N，接著持續推力最大值到 31 秒，之後開始下降，在 32 秒時推力降為 燃料燃燒完畢。火箭的推力與 時間關係，以及火箭的重量與時間的關係如圖 3.14 示 全長 (m) 直徑 (m) 重量 (ton) 15 11 6 3 1.7 0.6 0.6 1 1 1 50 (kg) 0.5 2.5 8 22 直徑 (m) 重量 (ton) 全長 (m)

0 20 40 60 80 0 1 2 3x 10 5 time (s) M o m e n t In e rti a (m ) I x I_y I_z 0 20 40 60 80 0 5 10 15 time (s) m a s s c e nt er (m ) CG x CG_y CG_z 0 20 40 60 80 0 5 10x 10 5 time (s) th ru s t ( N ) 0 20 40 60 80 1 1.5 2 2.5x 10 4 time (s) ma s s (m) 圖 3.14 火箭的轉動慣量，質心位置，推力分布，質量 飛行時間t=0s 飛行時間t=32s x CG 10.85m 8.63m y CG 0m 0m z CG 0m 0m x I 2 18156kg−m 6131kg−m2 y I 224390kg−m2 105492kg−m2 z I 2 197690kg−m 100680kg−m2 表 3.1 火箭的質量中心位置與轉動慣量 27

表 3.1 中的CG ，_x CG_y以及CG_z分別表示火箭的質量中心位置；CG 為在火_x 箭的X 軸上質量中心至火箭頂端的距離。_B I 、_x I 、_y I 為火箭的在_z 座標下的轉 動慣量。在飛行時間 到 B S 0 t= s t=32s期間，火箭的質量中心位置以及轉動慣量的 計算，我們將根據表 3.1 的數據做線性內插計算；而在飛行時間 到 期間，由於燃料已燃燒完畢，所以各項數據將維持在飛行時間 的量值。 32 t= s t=63s 32 t= s 火箭的質心位置以及轉動慣量與時間的關係圖如圖 3.14 示。 3.5.2 火箭的自由飛行軌跡 介紹完火箭的物理性質後，在模擬有重力轉彎控制器的控制之下的飛行 前，我們可以先模擬看火箭在沒有控制下的飛行情況。 火箭在第一階段的飛行時間總共為 63 秒，前 32 秒為推進器燃燒燃料使火箭 加速飛行；32 秒至 63 秒推進器已無推力，讓火箭做自由飛行。在本小節中我們 將推力向量系統中的角度皆設為零，來模擬火箭自由飛行的結果，如圖 3.15 及 圖 3.16。 火箭的發射條件為： z 發射角度：−2 z 飛行時間：63秒（前32秒為有推力的飛行，32秒至63秒為自由飛行） z 火箭的推力成梯形分布 28

0 20 40 60 80 0 1 2 3 4 time (s) M ac h num b er 0 20 40 60 80 0 5 10x 10 5 time (s) th rus t (N ) 0 20 40 60 80 0 5 10 15x 10 4 time (s) d y na m ic pr e s s ure (N /m 2 ) 0 20 40 60 80 0 2 4 6 8 time (s) a c c e le ra ti o n (G ) 圖 3.15 馬赫數，推力分布，動態壓力，火箭的加速度 0 20 40 60 80 -6 -4 -2 0 2 time (s) E u ler an gl e ( d eg ree ) roll yaw pitch 0 20 40 60 80 -0.2 0 0.2 0.4 0.6 time (s) (d egree ) alpha beta 0 20 40 60 80 -0.4 -0.2 0 0.2 0.4 time (s) A n g u la r ve lo ci ty ( d e g re e /s) p q r 0 1 2 3 4 0 10 20 30 40 50 range (km) a lti tu d e ( k m ) 圖 3.16 火箭的尤拉角，攻角以及側滑角，角速度，高度對距離 29

從圖 3.16 中可以看到火箭的俯仰角，有些微的震盪現象。而火箭的攻角雖 然大多都維持在 0.1 度以內，但是有著一定頻率的震盪。而角速度方面，也可以 看到震盪的情形。而從圖 2.5 中的高度對距離的關係圖看到火箭在第一階段結束 時，飛行的高度高達 48 公里左右。 3.6 火箭重力轉彎的軌跡 在這節裡我們用3.4 節所設計出來的控制器，在火箭飛行的前 32 秒，控制 火箭作重力轉彎。 火箭的發射條件為： z 發射角度： 2 − z 飛行時間：63秒（前32秒為有推力的飛行，32秒至63秒為自由飛行） z 火箭的推力成梯形分布 0 20 40 60 80 0 1 2 3 4 time (s) M ac h num b er 0 20 40 60 80 0 5 10x 10 5 time (s) th rus t (N ) 0 20 40 60 80 0 5 10 15x 10 4 time (s) d y na m ic pr e s s ure (N /m 2 ) 0 20 40 60 80 0 2 4 6 8 time (s) a c c e le ra ti o n (G ) 圖 3.17 馬赫數，推力分布，動態壓力，質量 30

從圖 3.17 中可看出火箭的馬赫數在前 32 秒不斷的增加，到 32 秒時馬赫數 達到最大值為 3.94，之後便開始下降，至 63 秒時馬赫數為 2.5。火箭的推力成梯 形分布，0 秒時推力為零，在 1 秒時推力達到最大值 8.65╳105 N，然後持續到 31 秒時推力開始下降，至 32 秒時推力降為零。火箭所遭受最大的動態壓力為 1.466╳105 N/m2。火箭的加速度在 1 秒時為 3G，接著不斷地上升至 31 秒時加速 度為 6G。 0 20 40 60 80 -10 -5 0 5 time (s) E u le r a n g le (d e g re e ) 0 20 40 60 80 -0.1 0 0.1 0.2 0.3 time (s) (d e g re e ) 0 20 40 60 80 -0.6 -0.4 -0.2 0 0.2 time (s) A n g u la r v e lo c it y (d e g re e /s ) 0 2 4 6 0 10 20 30 40 50 range (km) al ti tu de ( k m ) roll yaw pitch alpha beta p q r 圖 3.18 尤拉角，攻角以及側滑角，角速度，以及高度對距離 圖 3.18 中可以看出火箭在做重力轉彎的時候，除了前 5 秒火箭還在加速狀 態的時候攻角會顯得有點大，5 秒後火箭的攻角開始變小，之後便慢慢維持在一 個相當小的角度；而側滑角一開始就維持在很小的角度。角速度以及攻角方面可 以很明顯的看出並沒有產生如圖 3.16 中的震盪。因為全程做重力轉彎的關係， 從尤拉角的圖可以看出火箭的俯仰角的改變量並不是很大，最多只轉向到-7.5 度。火箭的高度在 63 秒時可以飛到 47.04 公里，距離則為 4.73 公里。 31

0 20 40 60 80 -5 -4.5 -4 -3.5 -3 -2.5 -2 -1.5 -1 -0.5 0x 10 -3 time (s) th ru s t an gl e s ( y 1 ~y 3 ) (de g re e) y₁ y₂ y 3 0 50 100 -0.005 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 time (s) th ru s t an gl e s z 1 ~z 3 (d e g re e ) z₁ z₂ z 3 圖 3.19 推力向量控制系統的角度 圖 3.19 中所顯示出的角度只有一個，但實際推力向量系統中在兩個軸上的 角度都各有三個，是因為火箭在飛行的過程中，同軸的角度都會維持一樣。因為 火箭做重力轉彎的關係，所以可以看到推力向量系統在 Z 軸的角度很小；而在 Y 軸的角度可以看到比 Z 軸來得更小，也是因為側滑角控制得比攻角來得小。 32

第四章 垂直發射控制

33 在第三章中我們可以知道火箭在越大的初始角度下發射，經由重力轉彎後所 產生的俯仰角也越大。然而若發射的角度太大，則火箭有落地的危險。因此大部 分的火箭的發射角度都採小角度發射，即幾乎是垂直發射。雖然在氣動力影響大 的環境下飛行，重力轉彎是個風險相當低的飛行策略，可惜的是重力轉彎的軌跡 卻不是一般飛行的最佳軌跡。若我們希望能夠在飛行 63 秒以後，火箭的俯仰角 能夠達到 之間，則從圖 3.3 中可以知道火箭的發射角必須在起碼在 左右發射，對於大型火箭這樣的發射角度可能不切實際。因此在本章節中將會介 紹改善的方法，4.1 節介紹垂直發射控制，4.2 節中設計姿態控制器，4.3 節則是 垂直發射控制的模擬。 25 30 − ∼− −8 4.1 垂直發射 我們將火箭採垂直發射，然後再在火箭發射的前 16 秒左右，讓火箭的俯仰 角先偏轉約15 ，之後再讓火箭接著做重力轉彎的飛行；由於俯仰角已經變大， 所以火箭在做重力轉彎的飛行時所偏轉的俯仰角也會變大。在做重力轉彎飛行 時，由於火箭的攻角會非常的小，因此推力向量控制系統所轉的角度並不會太 大；但要讓火箭在前 16 秒先作偏轉俯仰角的動作，則必須完全靠推力向量控制 系統的力量。火箭在做偏轉動作的這段期間裡，攻角會變的很大，使得火箭所受 到的氣動力變大，有可能會讓火箭推力向量控制系統所轉的角度過大。根據(2.22) 以及(2.23)式，氣動力除了跟由攻角決定的係數有關，也跟速度的平方成正比， 因此為了讓火箭能夠比較容易的達到偏轉的動作，我們必須降低火箭在前 16 秒 內的速度。而我們降低火箭速度的做法則是改變火箭燃料的燃燒速度，來改變火 箭的推力分布。

0 10 20 30 40 50 60 70 0 1 2 3 4 5 6 7 8 9x 10 5 T h rus t(N ) time (s) original thrust changed thrust 圖 4.1 火箭的推力分布 如圖 4.1 所示，紅色虛線是原本的推力分布，藍色實線是更改過的推力分布。 我們將火箭推力在前 16 秒的時候降到約為原本的一半，然後 16 秒的時候開始增 大，再 17 秒的時候增加到約 ，維持到 47 秒後開始降低，再 48 秒的 時候推力降至零。雖然推力降低了，但是燃燒時間變長了，不過兩者對時間積分 出來的面積會是一樣的，因為只是改變了燃料燃燒的速率，並沒有改變燃料的總 質量。推力降低的好處是整體的火箭速度降低了，所以在(2.22)以及(2.23)式中的 5 6.32 10 N× 2 1 2ρV 會比原來的小，也就是說所遭受到的氣動力會比較小，比較容易控制。 4.2 火箭姿態控制器 在第三章的重力轉彎控制器的輸入為角速度，因此都以角速度當作回授；而 火箭姿態控制器所要控制的是尤拉角，因此我們用一個雙積分迴路來控制，以角 速度以及尤拉角當作回授，如圖 4.2 所示 34

input k3 k4 G s( ) q q eqC

δ

₋

_θ

θ

1 s output Desired pitch angle C

θ

_q input k3 k4 G s( ) q q eqC

δ

₋

_θ

θ

1 s output Desired pitch angle C

θ

_q 圖 4.2 姿態控制器的方塊圖 圖4.2 中系統函數 的求法跟第三章一樣，然而我們改變了火箭的推力分 布，系統函數會有所不同。經由第三章的方法，發現改變了火箭的推力分布後所 求出來的近似函數在暫態部分改變不多，比較不同的是系統的時變量，即在(3.14) 式中的k值不一樣。則閉迴路轉移函數為 ( ) G s 3 4 2 4 3 3.5763 (0.4267 3.5763 ) 3.5763 3.5763 C kk k s kk s kk k₄ θ θ = + + + + (4.1) 這裡的設計方法跟 3.3 節一樣，經過設計調整後得到 k₃ =0.8， k 為一個隨 時間變化的參數，如圖 4.3 所示。 4 0 2 4 6 8 10 12 14 16 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 time (s) k4 圖 4.3 增益值k₄ 35

從(4.1)式中可以發現轉移函數跟(3.16)一樣，但計算出來的增益並不一樣， 因為將火箭的推力燃料做了調整，讓火箭剛離架後的速度不至於過快，使得系統 函數 裡的 不一樣外；此外轉向控制器的重點在於跟隨角度，所以設計出來 的控制器參數比較著重在上升時間上，希望能盡量跟上所給的輸入。 ( ) G s k 至於輸入部分，在滾轉角以及偏轉角部分輸入皆為零，俯仰角部份則給一個 二次曲線如圖 4.4 所示。 4.3 火箭垂直發射控制的軌跡 模擬結果 I 本節的模擬為先讓火箭藉由4.2 節所設計的姿態控制器讓火箭在前 16 秒的 時候跟隨所設計的俯仰角，如圖 4.4，使其做一個轉向的動作；飛行到 16 秒時則 切換 3.3 節所設計的重力轉彎控制器，讓火箭作重力轉彎的飛行。 z 發射角度： 0 z 飛行時間：63秒（前48秒為有推力的飛行，48秒至63秒為自由飛行） z 火箭的推力為兩段式的階梯狀分布 0 2 4 6 8 10 12 14 16 -16 -14 -12 -10 -8 -6 -4 -2 0 2 pi tc h an gl e ( d egr e e) time (s)

desired pitch angle simulation resullt

圖 4.4 設計的俯仰角與實際模擬的俯仰角 36

圖 4.4 中藍色的線條為我們希望火箭跟隨的俯仰角，紅色的線條為火箭模擬 的俯仰角。我們希望火箭在前 16 秒的時候可以轉向至-15 度，而實際上模擬的結 果也在 16 秒時轉向至-15 度。所設計的俯仰角軌跡與模擬出來的軌跡差異，可能 是推力向量系統中有著約 0.2 秒的延遲時間。 0 20 40 60 80 0 1 2 3 4 time (s) M ac h num b er 0 20 40 60 80 0 2 4 6 8x 10 5 time (s) th rus t (N ) 0 20 40 60 80 0 5 10x 10 4 time (s) d y na m ic pr e s s ure (N /m 2 ) 0 20 40 60 80 0 2 4 6 time (s) a c c e le ra ti o n (G ) 圖 4.5 馬赫數，推力分布，動態壓力，質量 從圖 4.5 中可以看出我們改變了火箭推力的分布，也因此改變了火箭的其他 性質。火箭可以達到的最大馬赫數比改變前來得小，約為 3.5 左右；但所遭受的 最大動態壓力也比改變前來得小，其最大值將近 105 N/m2。加速度方面在第一階 段推力的時候，加速度從 1.6G左右加速到 1.7G；第二階段的推力，使得火箭的 加速度從 2.4G加速到 4.3G。 37

0 20 40 60 80 -40 -20 0 20 time (s) E u le r a n g le (d e g re e ) 0 20 40 60 80 -2 -1 0 1 2 time (s) (d eg re e) 0 20 40 60 80 -2 -1.5 -1 -0.5 0 0.5 time (s) A n g u la r v e lo c it y (d e g re e /s ) 0 5 10 15 0 10 20 30 40 range (km) al ti tu de ( k m ) roll yaw pitch alpha beta p q r 圖 4.6 尤拉角，攻角以及側滑角，角速度，以及高度對距離 從圖 4.6 中可以看出火箭在飛行 63 秒後所轉向的俯仰角可以達到-31.89 度。 攻角部分可以看出火箭在一開始轉向的時候產生的攻角很大，至 16 秒轉向結束 的時候攻角將近 1 度；而在前 16 秒裡，火箭的攻角會變相的原因，是因為跟隨 所設計的俯仰角的關係。一開始火箭為了跟上所設計的俯仰角，因而轉動推力向 量系統的角度；而在接近 16 秒的時候，火箭的推力向量系統的角度必須轉到另 一個方向來使得火箭的俯仰角能夠慢慢的達到我們所設計俯仰角。而之後火箭作 重力轉彎的時候，攻角依舊維持在一個很小的角度。在角速度方面，可以看到在 約 16 秒的時候會有一個明顯的震盪，那是因為火箭在 16 秒的時候開始了第二段 的加速所導致；除此之外，由原本的角度控制器切換成角速度控制器也是造成震 盪的原因。高度方面在 63 秒時最大高度為 30.84 公里，距離為 14.65 公里。 38

0 20 40 60 80 -20 -15 -10 -5 0 5x 10 -3 time (s) thrus t angl es ( y 1 ~y 3 ) (degr ee) y₁ y 2 y₃ 0 20 40 60 80 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 time (s) thr us t angl es z 1 ~z 3 ( degree) z₁ z 2 z₃ 圖 4.7 推力向量控制系統的角度 圖 4.7 可以看出在前 16 秒轉向的時候，因為跟隨角度的關係，導致火箭推 力向量系統在 Z 軸方向，必須轉一個較大的角度來使火箭做一個轉向的動作。 而 16 秒後因為作重力轉彎的關係，所以火箭推力向量系統所轉的角度也就慢慢 的變小。至於在 Y 軸方向，則側滑角依舊維持著相當小的角度，所以在這方向 轉的角度也不大。

模擬結果 II 在本小節的模擬中，我們再次更改火箭的推力分布如圖4.9 所示，然後一樣 在前 16 秒的時候用 4.2 節設計控制器的方法跟隨圖 4.4 中所設計的俯仰角，而 16 秒之後改用第三章的方法讓火箭作重力轉彎。 39

z 發射角度： 0 z 飛行時間：63秒（0秒至63秒皆為有推力的飛行） z 火箭的推力為兩階段的梯形分布 0 2 4 6 8 10 12 14 16 -16 -14 -12 -10 -8 -6 -4 -2 0 pi tc h an gl e ( d egr e e) time (s)

desired pitch angle simulation resullt 圖 4.8 設計的俯仰角與實際模擬的俯仰角 圖 4.8 可以看出在更改火箭的推力分布後，在前 16 秒的時候火箭依舊可以 跟隨著所設計的俯仰角來飛行，只是在一開始俯仰角的下降速度比 4.3.1 節來得 慢一些，而在接近 16 秒的部分，俯仰角也沒有 4.3.1 節來得平緩；而這些差異也 說明了接下來的圖 4.10 中攻角在前 16 秒的分布情形與模擬結果 I 中的圖 4.6 會 有所不同。 40

0 20 40 60 80 0 1 2 3 4 time (s) M ac h num b er 0 20 40 60 80 0 2 4 6x 10 5 time (s) th rus t (N ) 0 20 40 60 80 0 2 4 6 8x 10 4 time (s) d y na m ic pr e s s ure (N /m 2 ) 0 20 40 60 80 0 1 2 3 4 time (s) a c c e le ra ti o n (G ) 圖 4.9 馬赫數，推力分布，動態壓力，質量 圖4.9 中可以看出我們將火箭的推力分布更改成在 63 秒內皆為有推力的梯 形分布，其推力的最大值為 。由於一直都有推力的關係，火箭的馬赫 數會一直上升，在 63 秒的時候達到最大值約 3.24。動態壓力方面則在 50 秒左右 達到最大值為 。在第一階段的加速度部分，火箭的加速度在 1 秒的 時候達到 0.7G 左右，在 16 秒的時候加速到 0.95G；第二階段的部分，在 17 秒 的時候為 1.4G，而在 62 秒時 則加速到 3.1G。 5 4.33 10 N× 4 7.8 10 N/m× 2 41

0 20 40 60 80 -60 -40 -20 0 20 time (s) E u le r a n g le (d e g re e ) 0 20 40 60 80 -1 0 1 2 3 time (s) (d eg re e) 0 20 40 60 80 -3 -2 -1 0 1 time (s) A n g u la r v e lo c it y (d e g re e /s ) 0 5 10 15 20 0 5 10 15 20 range (km) al ti tu de ( k m ) roll yaw pitch alpha beta p q r 圖 4.10 尤拉角，攻角以及側滑角，角速度，高度對距離 圖 4.10 中可以看到在尤拉角部分所轉向的俯仰角比圖 4.6 中所顯示出的結果 還要大，63 秒時的俯仰角角度為-54.54 度。而在攻角方面除了在前 16 秒作轉向 的動作導致會有比較大的攻角產生外，16 秒後便維持在很小的角度。側滑角方 面一樣仍然維持在很小的角度。角速度方面，在 16 秒左右的時候可以看到一個 小震盪，這是因為由姿態控制器切換到重力轉彎控制器之間所造成的。高度對距 離的關係圖可以看到高度最高 15.88 公里，也是因為俯仰角的轉向變大的關係； 63 秒結束所飛行的距離為 17.49 公里。 42

0 20 40 60 80 -0.035 -0.03 -0.025 -0.02 -0.015 -0.01 -0.005 0 0.005 0.01 time (s) th ru s t an gl e s ( y 1 ~y 3 ) (de g re e) y₁ y 2 y 3 0 20 40 60 80 -0.25 -0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 0.2 time (s) th ru s t an gl e s z 1 ~z 3 (d e g re e ) z₁ z 2 z 3 圖 4.11 推力向量控制系統的角度 圖 4.11 可以看出在推力向量系統中的 Z 軸方向在前 16 秒的時候會轉一個比 較大的角度，是為了跟隨所設計的俯仰角；16 秒之後作重力轉彎的關係，角度 便維持著較小的角度。而在 Y 軸的方向，推力向量系統的角度則一樣維持著很 小的角度。而在這個模擬中，可以看到兩個軸在 16 秒至 17 秒時，也就是切換控 制器的時候，都有震盪的產生，尤其是 Z 軸方向，震盪的幅度也比 4.2.2 節的模 擬來得更大，頻率也較高。 43

第五章 模擬結果與比較討論

本章主要是針對第三章以及第四章所做的模擬來進行討論，總共為三個模 擬如下表 5.1 所示，我們將三個模擬結果放在一起，並比較每種情況的優缺點。 發射角度 推力分布 方法 模擬 1 -2 度 0 至 32 秒推力呈梯形分布，32 秒至 63 秒則無推力 0 秒至 32 秒採重力轉彎 控制，32 秒至 63 為自由 飛行 模擬 2 0 度 0 至 16 秒為第一段推力，16 秒 至 48 秒為第二段推力，48 秒至 63 秒則無推力 0 秒至 16 秒作轉向的動 作，16 秒至 48 秒作重力 轉彎，48 秒至 63 秒為自 由飛行 模擬 3 0 度 0 至 16 秒為第一段推力，16 秒 至 63 秒為第二段推力， 0 秒至 16 秒作轉向的動 作，16 秒至 63 秒作重力 轉彎 表 5.1 0 20 40 60 80 -60 -40 -20 0 20 time (s) p it c h ang le ( de gree ) 0 20 40 60 80 0 1 2 3 4 time (s) M a c h nu m b er 0 5 10 15 20 0 10 20 30 40 50 range (km) a lt it ude (k m ) 0 20 40 60 80 0 5 10x 10 5 time (s) th ru s t ( N ) case1 case2 case3 case1 case2 case3 case1 case2 case3 case1 case2 case3 圖 5.1 俯仰角，馬赫數，高度對距離，推力分布 44

圖 5.1 可以看到火箭在模擬 1 中，只做全程的重力轉彎的話，則所轉向的 俯仰角很小，因此火箭飛行的高度會比較高。模擬 2 中我們修改了推力分布，目 的是讓火箭能夠在前 16 秒的時候，能夠比較順利的作一個轉向的動作，而 16 秒後再繼續做重力轉彎。從模擬 2 可以看出先轉向所造成的重力轉彎會比模擬 1 來大得多，因此模擬 2 的飛行高度也從模擬 1 中的 47 公里降至 31 公里。模擬 3 中我們再次修改了推力分布並沿用模擬 2 的方法，可以看出在前 16 秒也有一個 轉向的動作，跟模擬 2 差別不大，但接下來轉向的角度卻又比模擬 2 大得多，這 是因為火箭速度的關係。火箭在 16 秒的時候已經轉了一個約 15 度的俯仰角，也 加快了接下來作重力轉彎的角速度，而從圖 5.1 可以看到模擬 3 的馬赫數又比模 擬 2 來得小，因此速度上的差異便會很明顯的反映在重力轉彎的角速度上，這點 可以從圖 5.2 中看出，也使得模擬 3 所轉的角度會比模擬 2 來大得多。 0 20 40 60 80 -3 -2 -1 0 1 time (s) A ngul ar v el oc it y q ( degr ee/ s ) 0 20 40 60 80 -2 0 2 4 time (s) al pha ( degr ee) 0 20 40 60 80 -0.2 0 0.2 0.4 0.6 time (s) bet a ( degr ee) 0 20 40 60 80 0 5 10 15x 10 4 time (s) dy nam ic pr es s ur e case1 case2 case3 case1 case2 case3 case1 case2 case3 case1 case2 case3 圖 5.2 角速度 q，攻角，側滑角，動態壓力 45