高空气球外形计算及通用算法设计

周江华，金伟城,2，李智斌

(1. 中国科学院空天信息创新研究院， 北京 100094； 2. 中国科学院大学， 北京 100049；3. 山东科技大学 电气与自动化工程学院， 山东 青岛 266590)

高空气球利用充入密度小于空气的气体所产生的浮力飞行，是目前可在临近空间工作的几种飞行器之一，具有价格低、工作周期短、吊舱易回收、飞行时间可以根据需要灵活控制等优点，为在高层大气中进行研究提供了可靠的低成本平台。就其结构形式来说，高空气球又分为底部排气管与大气连通的零压式气球以及全封闭的超压式气球[1]。

高空气球的外形设计以及升空过程中外形的变化是浮空器设计以及浮空器飞行动力学研究的基础。气球外形的第一个理论性成果是Upson假设仅存在纵向张力的情况下提出的“自然形”气球[2]。由于其以Σ作为球形控制因子，因此也称为Σ-shape。现今的球形设计大多直接采用“自然形”或以“自然形”为基础：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-Ⅱ工具，对转化得到的最优控制问题进行求解。计算结果与文献数据对比表明了所提方法的有效性。随后分析球形数据，设计了通用算法。该算法虽需通过迭代计算球形，增加了计算时间，但取消了调参过程，且相较于多重打靶法，计算代价已大大减小，便于工程应用。