高空氣球外形計算及通用算法設計

周江華，金偉城,2，李智斌

(1. 中國科學院空天信息創新研究院， 北京 100094； 2. 中國科學院大學， 北京 100049；3. 山東科技大學 電氣與自動化工程學院， 山東 青島 266590)

高空氣球利用充入密度小于空氣的氣體所產生的浮力飛行，是目前可在臨近空間工作的幾種飛行器之一，具有價格低、工作周期短、吊艙易回收、飛行時間可以根據需要靈活控制等優點，為在高層大氣中進行研究提供了可靠的低成本平臺。就其結構形式來說，高空氣球又分為底部排氣管與大氣連通的零壓式氣球以及全封閉的超壓式氣球[1]。

高空氣球的外形設計以及升空過程中外形的變化是浮空器設計以及浮空器飛行動力學研究的基礎。氣球外形的第一個理論性成果是Upson假設僅存在縱向張力的情況下提出的“自然形”氣球[2]。由于其以Σ作為球形控制因子，因此也稱為Σ-shape?，F今的球形設計大多直接采用“自然形”或以“自然形”為基礎：Yajima在文獻[3]中指出“自然形”氣球的極限形式為歐拉體，而歐拉體作為一種無加強筋的超壓球是各種超壓球的設計基礎[4-5]。

雖然“自然形”氣球外形方程具有簡單的形式，但無法獲得解析解。現一般采用打靶法進行數值求解。對完全膨脹的設計球形求解：Smalley利用打靶法以Σ表的形式給出了Σ在0～1之間的部分球形數據[6]；姜魯華采用Gill法并以計算機程序的方式給出0～1之間任意Σ值的氣球外形[7]。對于Σ≥1時的氣球外形，并未見到有相關文獻給出參考數據，且2.2節結果顯示，Σ從1.5變化到2.0，其初始角度只改變了1°左右，可見此時初值的敏感度較高，打靶法求解困難。

上升段球形分析對升空過程中的動力學分析以及工程應用有重要的參考價值。Farley認為升空段氦氣泡為水滴形，頂部比較平坦，將其用圓球阻力系數近似并不合適，因此Farley用圓錐體近似升空段球形并給出參考阻力系數(Cd≈0.8)[8]。若能通過求解升空球形，獲得不同高度下氣泡最大橫截面積更準確的估計，便可更為精準地估算升空時的氣動阻力。上升段球形分析另一個重要的應用是估計升空段球頂壓差。高空氣球利用球頂排氣閥來調節升空段的浮力，而浮力排出率與球頂壓差相關，零壓氣球通常不安裝壓差傳感器，需要依賴事前估計。通過求解升空球形得到的球頂壓差比采用圓球近似更為準確[9]。

然而高空氣球為由極薄聚乙烯(PolyEthylene, PE)膜材料構成的軟式充氣結構，隨著升空高度的變化，氣球體積不斷膨脹，導致其外形輪廓的變化非常之大。相較于平飛段(滿充)外形，升空段球體大部分情況下底部呈聚束狀，初始角度幾乎為零且需一同猜測初始角度及零壓面高度兩個參數。Smalley[10]將底部聚束視為負載，用打靶法求解外形，但計算不同高度的球形都需調整模型。Baginski等[11]采取解方程的思路：基于MATLAB的fsolve函數，利用多重打靶法求得底部張角及零壓面高度，再借助龍格庫塔法獲得氣球外形，但氣球高度越低，分段數越多，計算量越大。楊燕初等[12]基于MATLAB的fmincon函數，結合多重打靶法求解外形，計算時長依賴于初值的選取。

試驗飛行需要充足的高空氣球升空過程中的形體數據，原有的計算方法需花費大量的時間和精力，因而需要一種簡單快速的計算方法來求解“自然形”氣球外形方程。

高斯偽譜法因其具有收斂性好以及初值敏感度低的優點，已經成為一種較為常用的解決最優控制問題的算法，而基于高斯偽譜法的高斯偽譜法最優化軟件包(Gauss Pseudospectral OPtimization Software,GPOPS-Ⅱ)也成為解決最優控制問題的有力工具之一。在此基礎之上，本文提出將“自然形”外形方程轉化為最優控制描述形式；然后用GPOPS-Ⅱ 工具求解的思路，以期精確、簡捷得到氣球完全膨脹以及升空狀態下的外形數據；同時設計通用算法，以實現無須調參便可求解任意Σ、任意高度下的“自然形”球形的目標。

1 “自然形”氣球外形方程

“自然形”氣球也稱為Σ外形氣球?？捎扇缦铝鶄€常微分方程表示[6,13]：

(1)

(2)

圖1 Σ外形氣球示意圖Fig.1 Schematic diagram of Σ-shape balloon

對于Σ氣球，設計球形及升空段球形的計算，本質上是求解一個兩點邊值問題：

2 數值求解

GPOPS-Ⅱ是用于解決非線性最優控制問題的通用MATLAB軟件包。該軟件包利用可變階的高斯偽譜法，將最優控制問題離散并轉換為非線性規劃問題，通過非線性規劃問題求解器對非線性規劃問題進行求解[14-15]。

由于最優控制問題本身就是一個兩點邊值問題，因此，GPOPS-Ⅱ可用于氣球外形的計算。但需對氣球模型方程組做適當的處理，將其轉換為最優控制問題。

2.1 GPOPS-Ⅱ求解工具介紹

GPOPS-Ⅱ基于高斯偽譜法，其求解最優控制問題的思路是將連續的最優控制問題離散化為非線性規劃問題，然后通過Karush-Kuhn-Tucker條件求解[16-17]。Benson證明了高斯偽譜法滿足協調映射定理，即利用高斯偽譜法離散最優控制問題后得到的非線性規劃問題的Karush-Kuhn-Tucker條件與最優一階必要條件的離散形式完全等價[18]，這使得高斯偽譜法同時具備了直接法的快速性與間接法的精確性。

在GPOPS-Ⅱ中，需要設定變量的初值、終值以及變量范圍三組參數。得益于高斯偽譜法對初值、終值敏感度低的優點，初值、終值只需猜測一次便可，在求解不同球形時無須調整，因而對于不同的球形，只需猜測變量范圍即可。相較于多重打靶法每次計算需猜測數十個參數[11-12]，此方法可將其減少至五個左右；同時一組參數可以計算一段Σ(或高度)區間的外形，與打靶法相比具有通用性。

2.2 設計球形的求解

2.2.1 微分方程組的處理

(3)

(4)

2.2.2 設計球形數值解

文獻[6]用打靶法給出了部分參考結果，以其計算結果作為初值并用四階龍格庫塔法驗算，發現終值符合要求。表1給出了本文與文獻[6]計算結果比較，可以看出，兩者的計算結果基本一致，說明了本文計算方法的有效性。同時表2補充了Σ≥1時的部分參考數據，從表中的數據也可發現：Σ從1.5變化到2.0，其初始角度只改變

表1 高斯偽譜法計算結果與Smalley計算結果對比

表2 部分Σ≥1的計算結果

圖2 大Σ外形與歐拉體外形的對比圖Fig.2 Comparison of the shape of large sigma and Euler′s elastica

2.3 升空狀態球形求解

2.3.1 升空微分方程組處理

(5)

(6)

2.3.2 升空狀態球形數值解

由于當τb≥10之后，θ0基本等于零，此時無法用四階龍格庫塔法驗證，文獻[10-11]分別采用打靶法、多重打靶法，得到的計算結果基本一致，因此從側面證實了計算結果的可靠性。表3給出了本文與Baginski計算結果[11]的比較，可以看出，兩者的計算結果基本一致，說明了本文計算方法的有效性。表4給出了Σ=2.0時，氣球升空過程參考計算結果。圖3給出了Σ=2時的“自然形”氣球從地面(海拔為零)到升限高度(海拔為20 km)的歸一化升空外形曲線。從表4以及圖3可以得到球形隨高度變化規律：隨著氣球升高，零壓面高度減小，氣球高度減小，最大半徑增加，在16～20 km處球形變化劇烈。

表3 Σ=0.4的計算結果

表4 Σ=2.0時，氣球升空過程參考計算結果

圖3 Σ=2.0時，部分高度球形曲線Fig.3 Balloon profile for partial hight values while Σ=2.0

3 通用算法設計

表5 平飛狀態參數取值

算法1 球形計算通用算法

其中，算法1中21組升空狀態計算參數見表6,m0=500。

表6 升空狀態計算參數取值

根據測試，算法1對常用Σ以及任意整數τb均能完成計算，且平均用時20 s左右，而用文獻[12]的方法至少需4 min(相同計算機)，且計算時長對初值的選取更敏感。本算法在計算代價上已大大減小，且無須猜測參數，便于快速獲得大量所需球形數據。

4 結論

本文對“自然形”氣球外形方程求解進行了系統性分析，將原問題進行無量綱化處理后轉換為最優控制形式，借助GPOPS-Ⅱ工具，對轉化得到的最優控制問題進行求解。計算結果與文獻數據對比表明了所提方法的有效性。隨后分析球形數據，設計了通用算法。該算法雖需通過迭代計算球形，增加了計算時間，但取消了調參過程，且相較于多重打靶法，計算代價已大大減小，便于工程應用。