2. 航天飞行动力学技术国家重点实验室, 陕西 西安 710072
20世纪60年代开始，国内外学者对太阳帆航天器进行了深入的研究和分析，除了实现无推进剂消耗的轨道转移之外，太阳帆航天器还能实现一些特殊的非开普勒轨道任务，如悬浮轨道任务、人工拉格朗日点任务等[1]。
McInnes[2]首次阐述了利用太阳帆构造广义二体问题的思想，在该文献中，作者使用正、反2种方法研究了航天器在广义引力场中的运动特性，通过调节太阳帆产生的推力，实现共面圆轨道之间的轨道转移和中心引力场逃逸。然而，如果单纯采用太阳帆推进模式，由于其产生的推力很小，轨道转移常常需要很长时间。针对这个问题，Mengali[3]基于McInnes的广义二体问题模型，提出了太阳帆和电推进结合的混合推进方法，定量分析了太阳帆推力对电推进燃料消耗的影响。
近年来，针对航天器连续推力机动轨道设计问题，为了提高计算效率，快速地给轨道优化算法提供可行有效的优化初值，国内外学者提出了基于形状的轨道设计方法(简称形状方法)，如正弦指数法[4-5]、逆多项式法[6-7]、对数螺旋线法[8-9]、速度曲线法[10]、傅里叶级数法[11-14]等。形状方法是轨道设计的反方法，是先通过一系列参数表示的轨道特性曲线来近似推力作用下的轨道，然后反推航天器沿这一形状运动所需的推力[9]。但是目前为止，未见到将基于形状的轨道设计方法用于太阳帆和电推进结合的混合推进模式的航天器轨道设计中。
本文依据现有研究成果，引入太阳帆/电推进结合的混合推力模式，融合形状方法的基本思想，研究了混合推力模式下的基于傅里叶级数展开的航天器机动轨道设计方法，并分析了混合推力对傅里叶级数形状轨道和推力产生的影响。首先，根据太阳帆/电推进混合推力模型，建立航天器在极坐标系下的运动模型；然后，依据傅里叶级数展开式，将轨道设计问题转化为关于傅里叶系数的非线性规划问题，并应用Matlab中的fmincon函数求解该非线性规划问题。另外，本文改进原有的形状轨道设计方法，研究了航天器轨道转移的时间最优问题。在考虑航天器任务约束的情况下，基于上述方法设计了混合推力作用下的“地球-火星”转移轨道。仿真分析了太阳光压对形状轨道设计方法的影响，验证了应用混合推力设计形状轨道的可行性，并且在一定范围内可以得到燃料消耗更优的航天器连续推力转移轨道。
1 太阳帆推力模型假设太阳帆的表面为平面，且太阳帆的法向与太阳光方向一致，则作用在太阳帆上的太阳光压加速度可以表示为[1]：
|
(1) |
式中，τ为太阳帆的反射系数; As为太阳帆的面积; S0为一个天文单位处的太阳能通量; AU为一个天文单位长度, 即地球到太阳的平均距离; c为光速; m为太阳帆航天器的质量; rs为太阳帆航天器到太阳的距离。
由于太阳光压力和太阳引力的大小都与太阳帆到太阳的距离的平方成反比, 对于给定的太阳帆, 太阳光压力与太阳引力的比值为常数。定义光压因子η为作用在太阳帆上的光压力和引力的比值, 则η与太阳帆的面质比有关, 比值越大, η越大。太阳帆受到的光压力除了与帆的面质比有关外, 还与太阳帆的姿态有关, 则太阳帆受到的光压加速度可以表示为:
|
(2) |
式中，μs为太阳的引力常数, β为太阳帆的法线方向与太阳光方向的夹角。
2 航天器轨道机动模型假设不考虑中心天体自转及所有摄动力的影响, 在极坐标系下, 考虑太阳光压力, 航天器绕中心天体运动的单位质量运动方程可以描述为
|
(3) |
式中，r为航天器到中心天体中心的距离; θ为极角; u为航天器径向速率; v为航天器周向速率; ae为电推进器产生的推力加速率; α为加速度转向角。
由于本文需要考虑太阳光压力, 航天器绕太阳运动, 如图 1所示, 则(3)式中的μ特指太阳引力常数, r特指航天器到太阳的距离。假设太阳帆的法向与太阳光方向一致, 即β=0, 则航天器动力学模型可表示为:
|
| 图 1 太阳帆/电推进混合推力下轨道机动示意图 |
|
(4) |
研究表明, 实际的最优推力方向是随时间变化的, 沿切向施加推力非常接近于最优推力方向。因此, 本文假设电推进推力加速度ae的方向与速度V的方向平行, 即ae转向角与飞行航迹角ξ满足α=ξ+nπ, 其中n=0、1, 如图 2所示。
|
| 图 2 加速度方向假设 |
由(4)式的第三式和第四式, 可得:
|
(5) |
由切向推力假设, 可得推力加速度转向角α与飞行航迹角ξ有如下关系:
|
(6) |
将(4)式的第一、二式和(6)式代入(5)式, 则(5)式可以改写为:
|
(7) |
根据傅里叶级数展开式, 距离r和角度θ可以近似为:
|
(8) |
式中，x代表(r, θ); nx为(r, θ)傅里叶级数展开式的项数; (ax0, axi, bxi)为傅里叶系数; T为航天器总飞行时间。
将(8)式代入运动方程(7)式, 以此将运动方程式(7)转化为只包含傅里叶系数和时间t的非线性代数方程:
|
(9) |
由此可见, 混合推力下的航天器转移轨道设计问题转化为了关于傅里叶系数的非线性规划问题。
在混合推力的作用下, 航天器的转移轨道满足(4)式, 由(4)式可得航天器电推进推力加速度ae大小的表达式为:
|
(10) |
式中, cosα=±cosξ, 即
|
(11) |
针对时间最优问题, 本节引入航天器推力广义加速率, 定义航天器绕中心天体运行1 rad由电推进产生的速率增量, 即
|
(12) |
广义加速度
利用(4)式中的第二式, 则航天器动力学模型可以转换为
|
(13) |
由(13)式的第二式和第三式, 可得:
|
(14) |
由切向推力假设, 可得
|
(15) |
将(13)式的第一式和(15)式代入(14)式, 可得:
|
(16) |
(16)式各项同除以v2, 可得:
|
(17) |
考虑切向加速度假设, 则(17)式可以改写为:
|
(18) |
化简可得:
|
(19) |
根据傅里叶级数展开式, 距离r和周向速率v可以近似为:
|
(20) |
式中，x代表(r, v); Θ=θf-θi+2Nrevπ, (θi, θf)分别为航天器始末位置的极角, Nrev为航天器绕中心天体运动的总圈数。
将(20)式代入运动方程式(19), 以此将运动方程式(19)转化为只包含傅里叶系数和角度θ的非线性代数方程。
|
(21) |
由此可见, 混合推力下的航天器转移轨道设计问题转化为了关于傅里叶系数的非线性规划问题。
航天器的机动轨道需满足(13)式, 由(13)式可得广义加速度
|
(22) |
式中，cosα=±cosξ, 即
|
(23) |
一般轨道设计的性能指标是时间最省或者燃料消耗最少, 即:
|
(24) |
在实际航天任务中, 发动机产生的推力大小是有限的。为此, 本文考虑了电推力器的推力约束:
|
(25) |
式中, amax为给定的最大推力。
4 仿真分析 4.1 燃料最优问题考虑航天器从地球转移到火星进行火星探测。本算例取距离单位为1 DU=1 AU, 时间单位为2π TU=1 year。该算例轨道转移的边界条件和输入参数如表 1所示, 表中η∈[0, 0.05]代表光压因子的取值范围。
| 参数 | 数值 |
| ri | 1 DU |
| θi | 0 rad |
| rf | 1.523 4 DU |
| θf | 3.548 rad |
 | 0 DU/TU |
 | 1 rad/TU |
 | 0 DU/TU |
 | 0.531 8 rad/TU |
| Nrev | 1 |
| nr | 3 |
| nθ | 4 |
| T | 13.433 TU |
| η | [0, 0.05] |
| amax | 0.05 DU/TU2 |
由于对于不同的光压因子, 航天器转移轨道相差不大, 所以仅以η=0.03为例, 图 3为光压因子η取0.03时的航天器转移轨迹; 图 5整合上述不同光压因子η取值下的电推进推力曲线; 图 7为航天器完成“地球-火星”轨道转移的特征速度随光压因子变化的曲线。
|
| 图 3 η=0.03时的形状轨道 |
|
| 图 4 η取值不同时相应的推力曲线 |
|
| 图 5 特征速度随η变化的曲线 |
|
| 图 7 η取值不同时相应的推力曲线 |
由图 3可见, 在太阳帆/电推进混合推力作用下, 应用傅里叶级数法设计航天器转移轨道是可行的。对比文献[14]仿真结果可知, 傅里叶级数的项数取值对推力曲线的影响较大。
在η=0的情况下, 由图 5通过对比有无推力大小约束2种情况的推力曲线可知:在有推力约束的情况下, 傅里叶级数法得到的推力曲线满足给定的推力约束。由图 5可知, 当光压因子η增大时, 在航天器飞行的前期, 电推进施加的最大推力加速度相应变小, 但是在后期, 其推力加速度相应增大。
由图 7可知, 在太阳帆/电推进混合推力作用下, 应用傅里叶级数法的航天器转移轨道的特征速度并不会随着光压因子η的增大而一直减小, 会存在一个极点, 决定该极点的因素除了始末边界条件, 还包括傅里叶展开式的项数、总飞行时间等。在本算例的初始条件下, 当η=0.03时, 特征速度最小。
4.2 时间最优问题考虑航天器从地球轨道转移到火星轨道进行火星探测。该算例轨道转移的边界条件如表 1所示, 输入参数如下:Nrev=1;nr=3;nv=3;η∈[0, 0.03];amax=0.03 DU/TU2。
然而在飞行时间自由情况下, 表 1给出的边界条件不能直接应用, 需要对其重新计算得到(r, v)及其关于极角θ的一阶导数, 如(26)式所示。
|
(26) |
图 6为不同光压因子取值下的航天器轨迹; 图 7为不同光压因子取值下的电推进推力曲线; 图 8为航天器完成轨道转移的飞行时间随光压因子变化的曲线。
|
| 图 6 η取值不同时相应的形状轨道 |
|
| 图 8 飞行时间随η变化的曲线 |
由图 6和图 7可知，针对带推力约束的航天器飞行时间最优问题，傅里叶级数法得到的推力曲线满足给定的推力约束。由此说明该方法可为飞行时间最优问题的轨迹优化提供初值。
由图 8可知，不同于燃料最优问题，航天器飞行时间会随着光压因子的增大而呈线性增加。由图 6和图 7可知，光压因子增大，电推进在飞行前期推力增大，导致航天器所经历的轨迹增长，从而飞行时间增加。
5 结论本文借鉴了当前太阳帆在轨道设计应用中的部分研究成果，引入太阳帆/电推进混合推力模型，结合傅里叶级数法的思路，研究了混合推力作用下的傅里叶级数法，并基于该方法设计了混合推力作用下的“地球-火星”转移轨道。另外，本文研究了基于傅里叶级数法的航天器轨道转移的时间最优问题。仿真算例验证了在具有推力约束的情况下应用傅里叶级数法设计混合推力转移轨道的可行性，得到了一些有意义的结论。
| [1] |
龚胜平. 太阳帆航天器动力学与控制研究[D]. 北京: 清华大学, 2009. Gong Shengping. Study on Dynamics and Control of Sail-Craft[D]. Beijing, Tsinghua University, 2009(in Chinese) http://cdmd.cnki.com.cn/Article/CDMD-10003-2009083598.htm |
| [2] | Mclnnes C R. Orbits in a Generalized Two-Body Problem[J]. Journal of Guidance Control and Dynamics, 2003, 26(5): 743-749. DOI:10.2514/2.5129 |
| [3] | Mengali G, Quarta A A. Trajectory Design with Hybrid Low-Thrust Propulsion System[J]. Journal of Guidance Control and Dynamics, 2007, 30(2): 419-426. DOI:10.2514/1.22433 |
| [4] | Petropoulos A E, Longuski J M. Shape-Based Algorithm for Automated Design of Low-Thrust, Gravity-Assist Trajectories[J]. Journal of Spacecraft and Rockets, 2004, 41(5): 787-796. DOI:10.2514/1.13095 |
| [5] | Izzo D. Lambert's Problem for Exponential Sinusoids[J]. Journal of Guidance Control and Dynamics, 2006, 29(5): 1242-1245. DOI:10.2514/1.21796 |
| [6] | Wall B J, Conway B A. Shape-Based Approach to Low-Thrust Rendezvous Trajectory Design[J]. Journal of Guidance Control and Dynamics, 2009, 32(1): 95-101. DOI:10.2514/1.36848 |
| [7] | Wall B J. Shape-Based Approximation Method for Low-Thrust Trajectory Optimization[C]//AIAA/AAS Astrodynamics Specialist Conference and Exhibit, Honolulu, USA, 2008 |
| [8] |
郑莉莉, 袁建平, 朱战霞. 基于对数螺旋线的非开普勒轨道设计[J]. 宇航学报, 2010, 31(9): 2075-2081.
Zheng Lili, Yuan Jianping, Zhu Zhanxia. Logarithmic Spiral-Based Non-Keplerian Orbit Design[J]. Journal of Astronautics, 2010, 31(9): 2075-2081. (in Chinese) |
| [9] |
郑莉莉. 连续推力作用下的机动轨道设计与优化[D]. 西安: 西北工业大学, 2011. Zheng Lili. Orbit Design and Optimization for Spacecraft Maneuver under Continuous Thrust[D]. Xi'an, Northwestern Polytechnical University, 2011(in Chinese) |
| [10] | Gondelach D J, Noomen R. Hodographic-Shaping Method for Low-Thrust Interplanetary Trajectory Design[J]. Journal of Spacecraft and Rockets, 2015, 52(3): 728-738. DOI:10.2514/1.A32991 |
| [11] | Taheri E, Abdelkhalik O. Shaped-Based Approximation of Constrained Low-Thrust Space Trajectory Using Fourier Series[J]. Journal of Spacecraft and Rockets, 2012, 49(3): 535-545. |
| [12] | Abdelkhalik O, Taheri E. Approximate on-off Low-Thrust Space Trajectories Using Fourier Series[J]. Journal of Spacecraft and Rockets, 2012, 49(5): 962-965. DOI:10.2514/1.A32307 |
| [13] |
王雪峰, 方群, 孙冲. 一种航天器空间机动轨道的改进形状设计方法[J]. 宇航学报, 2015, 36(11): 1242-1247.
Wang Xuefeng, Fang Qun, Sun Chong. An Improved Shape-Based Method for Spacecraft 3-D Maneuver Orbit Design[J]. Journal of Astronautics, 2015, 36(11): 1242-1247. DOI:10.3873/j.issn.1000-1328.2015.11.004 (in Chinese) |
| [14] | Wang X, Fang Q, Sun C, Yuan J. An Improved Shape-Based Method for Interplanetary Low-Thrust Orbit Design[C]//6th European Conference for Aeronautics and Space Sciences, Krakow, Poland, 2015 |
2. National Key Laboratory of Aerospace Flight Dynamics, Xi'an 710072, China
