下载:
-
太阳系边际探测旨在了解外太空物质分布、星系化学演化、日光层结构、附近星系介质等，探测太阳系边际的任务需要较长的转移时间和先进的推进系统。目前，只有两个航天器“旅行者1号”（Voyager 1）和“旅行者2号”（Voyager 2）[1-3]已到达太阳系边际。尽管“旅行者1号”和“旅行者2号”已取得了重大发现，但是探测器上携带的是30多年的科学仪器载荷，无法完成众多重要的观测任务，需考虑新的太阳系边际探测任务。采用化学推进或核能推进技术需要消耗足够多的燃料才有可能到达太阳系边际。因此，太阳系边际探测的一个重要的要求是使用性能更优的推进系统来执行相应的任务，以使燃料消耗尽可能少、任务周期尽可能短。
连续小推力航天器的发展，为深空远距离探测任务提出了一种更为优越的解决方案。在各类连续小推力航天器中，电动帆是一种新的推进概念[4-5]，利用太阳风的动压产生连续的推力，并不需要消耗燃料。电动帆可以在整个任务周期内提供连续推力，非常适合长距离的深空探测。在电动帆结构中，航天器绕着一个对称轴旋转，利用离心力展开许多细长的导电金属链。金属链被保持在高正电位，电子束大致沿着自旋轴射出，产生的静电场扰动了入射太阳风质子的轨迹，从而产生太阳风等离子体流到金属链的动量转移。对于电动帆这种连续小推力推进系统，在最优转移轨迹设计问题中，设计变量的急剧增加使得小推力飞行任务的设计在数值计算上具有极大的挑战性。
本文以电动帆作为主要推进系统，实现对逃逸太阳系转移轨迹的快速生成。由于电动帆的推力特性，其在飞行过程中产生连续的小推力，其任务设计需要一种近似飞行轨迹和任务成本的方法。在转移轨迹中，航天器需要满足运动学方程和边界条件，对于这个问题，间接或直接优化方法都需要一个合理的初始近似值，并且通常需要大量计算，导致在任务的初始设计阶段非常困难。因此，快速的初始轨迹设计对任务分析非常重要。
近年来，学者提出了使用形状法实现连续推力航天器转移轨迹快速生成。Petropoulos等[6]提出了形状法，使用指数正弦函数将航天器的飞行轨迹描述出来，然后通过优化相关系数满足相应的约束条件。之后，又有很多作者提出了多种形状函数[7-9]。Wall和Conway[8]等使用一个逆多项式来求解极坐标系中平面轨迹的径向距离，Xie等[9]提出一种考虑两个共面椭圆轨道之间轨迹约束的低推力轨迹近似，Peloni等[10]提出通过组合指数项和正弦项以近似航天器交会任务中的转移轨迹，Taheri和Abdelkhalik[11-13]提出使用有限快速傅里叶级数（Fast Flourier Series，FFS）形状法获得了低推力航天器的二维和三维转移轨迹；Huo等利用贝塞尔曲线[14]设计了电动帆探测小行星的转移轨迹，对轨迹终点位置和速度做了具体的约束，控制设计难度相对较小。
本文利用贝塞尔曲线快速设计电动帆逃逸太阳系的三维运动轨迹，只对轨迹终点的能量进行约束，通过将贝塞尔形函数方法获得的结果作为高斯伪谱方法（Gauss Pseudospectral Method，GPM）的初始解进行计算，并与几种方法做了对比，最后对贝塞尔形函数方法获得的结果对直接解算器初始解的适用性进行了评估。
-
电动帆在飞行过程中无燃料消耗，在该研究中以电动帆的转移时间作为性能指标，求解三维连续轨迹优化问题。
电动帆的推力模型和轨道动力学通常在两个笛卡尔坐标系下描述，即轨道参考系
${O_o}{x_o}{y_o}{z_o}$ 和日心黄道惯性系${O_i}{x_i}{y_i}{z_i}$ 。如图1所示，轨道坐标系的原点${O_o}$ 定义在电动帆的质心处，${z_o}$ 轴是沿太阳与电动帆连线的方向，${y_o}$ 轴垂直于${z_o}$ 轴和黄道面法线，形成右手系。日心黄道惯性系的定义：原点${O_i}$ 在太阳质心，${x_i}$ 轴定义在日分点方向，${z_i}$ 轴沿黄道面法线和黄道北极方向，${y_i}$ 轴形成右手系。最后，在贝塞尔方法中，圆柱坐标$(\rho ,\theta ,z)$ 也用来描述电动帆的位置。
图 1 参考坐标系
Figure 1. Reference coordinate system
在初步任务分析中，电动帆的推力矢量常用俯仰角
${\alpha _n}$ 和时钟角$\sigma $ 来描述。如图2所示，俯仰角${\alpha _n} \in [0, \text{π} /2]$ 是电动帆标称平面相对于${z_o}$ 轴的倾角，时钟角$\sigma \in [ - \text{π} , \text{π} ]$ 在${O_{{o}}}{x_{{o}}}{y_{{o}}}$ 平面上，是${x_o}$ 轴与电动帆自旋轴分量之间的夹角。$\alpha $ 是推力锥角，定义为推进加速度的方向向量与太阳航天器方向（${z_{{o}}}$ 轴）之间的角。推力矢量位于由${x_{{o}}}$ 轴和自旋轴矢量构成的平面内，因此${x_{{o}}}$ 轴与${O_{{o}}}{x_{{o}}}{y_{{o}}}$ 平面上推力矢量分量之间的夹角等于时钟角$\sigma $ 。
图 2 参考坐标系和推进加速度的特征角
Figure 2. Reference coordinate system and characteristic angle of propulsion acceleration
在圆柱坐标系中，
$\rho $ 为电动帆径向距离，$\theta $ 为方位角，$z$ 为电动帆的高度。在圆柱坐标系中描述电动帆的运动学方程为$$ \left\{ {\begin{array}{*{20}{l}} {\ddot \rho - \rho {{\dot \theta }^{\text{2}}} + \mu \rho /{s^3} = {a_\rho }} \\ {\rho \ddot \theta + 2\dot \rho \dot \theta = {a_\theta }} \\ {\ddot z + \mu z/{s^3} = {a_z}} \end{array}} \right. $$ (1) 其中：
${{s}} = \sqrt {{\rho ^2} + {z^2}}$ 为电动帆中心与太阳之间的距离；$\mu $ 为太阳的引力参数；${a_\rho }$ 、${a_\theta }$ 和${a_z}$ 为电动帆推力加速度的分量。在轨道参考系
${O_o}{x_o}{y_o}{z_o}$ 中描述的解析精确推力模型为$$ {\left[ {\begin{array}{*{20}{c}} {{a_{ox}}} \\ {{a_{oy}}} \\ {{a_{oz}}} \end{array}} \right]_{\rm refined}} = \frac{{{a_c}{r_ \oplus }\kappa }}{{2s}}\left[ {\begin{array}{*{20}{c}} {\cos {\alpha _n}\sin {\alpha _n}\cos \sigma } \\ {\cos {\alpha _n}\sin {\alpha _n}\sin \sigma } \\ {{{\cos }^2}{\alpha _n} + 1} \end{array}} \right] $$ (2) 其中：
${r_ \oplus }$ =1 AU是一个参考距离；${a_c}$ 是电动帆的特征加速度；$\kappa \in [0,1]$ 是推力系数，可以通过安装在电动帆上的电子枪来调节。根据轨道参考系与圆柱坐标系的几何关系，可得到在圆柱坐标系中推进加速度的3个分量为
$$ \begin{aligned}[b] {\left[ {\begin{array}{*{20}{l}} {{a_\rho }} \\ {{a_\theta }} \\ {{a_z}} \end{array}} \right]_{\rm refined}} = & \left[ {\begin{array}{*{20}{l}} {\cos \varphi }&0&{\sin \varphi } \\ 0&1&0 \\ { - \sin \varphi }&0&{\cos \varphi } \end{array}} \right]\left[ {\begin{array}{*{20}{l}} {{a_{ox}}} \\ {{a_{oy}}} \\ {{a_{oz}}} \end{array}} \right] \hfill =\\ & \left[ {\begin{array}{*{20}{l}} {\cos \varphi {a_{ox}} + \sin \varphi {a_{oz}}} \\ {{a_{oy}}} \\ { - \sin \varphi {a_{ox}} + \cos \varphi {a_{oz}}} \end{array}} \right] \hfill \end{aligned} $$ (3) 其中，
$$ \sin \varphi = {\rho \mathord{\left/ {\vphantom {\rho {\sqrt {{\rho ^2} + {z^2}} }}} \right. } {\sqrt {{\rho ^2} + {z^2}} }}\text{，}\cos \varphi = {z \mathord{\left/ {\vphantom {z {\sqrt {{\rho ^2} + {z^2}} }}} \right. } {\sqrt {{\rho ^2} + {z^2}} }}。 $$ 根据式（2）中在轨道参考系
${O_o}{x_o}{y_o}{z_o}$ 中描述的解析精确推力模型，获得关于${a_{ox}}$ 和${a_{oy}}$ 的方程为$$ a_{ox}^2 + a_{oy}^2 = {\left( {\frac{{{a_c}{r_ \oplus }\kappa }}{{2s}}} \right)^2}{\cos ^2}{\alpha _n}{\sin ^2}{\alpha _n} $$ (4) 同时，由式（2）可得：
${\cos ^2}{\alpha _n} = {{2s{a_{oz}}} \mathord{\left/ {\vphantom {{2s{a_{oz}}} {({a_c}{r_ \oplus }\kappa )}}} \right. } {({a_c}{r_ \oplus }\kappa )}} - 1$ ，将其代入式（4）中，可以得到关于推力系数的二次方程为$$ 2{\left( {\frac{{{a_c}{r_ \oplus }}}{{2s}}} \right)^2}{\kappa ^2} - 3\left( {\frac{{{a_c}{r_ \oplus }{a_{oz}}}}{{2s}}} \right)\kappa + (a_{ox}^2 + a_{oy}^2{\text{ + }}a_{oz}^2) = 0 $$ (5) 由式（4）可知，当
${\alpha _n} = {0^ \circ }$ 时，$\kappa = {a_{oz}}s/{a_c}{r_ \oplus }$ ，于是可得$$ \kappa = \frac{{s(3{a_{oz}} - \sqrt {a_{oz}^2 - 8a_{ox}^2 - 8a_{oy}^2} }\;)}{{2{a_c}{r_ \oplus }}} $$ (6) 在一般的三维轨迹设计问题中，电动帆在出发点和到达点的位置和速度需要满足以下12个边界条件为
$$ \left\{ \begin{gathered} \rho (\tau = 0) = {\rho _i},\rho (\tau = 1) = {\rho _f} \hfill \\ {\rho'}(\tau = 0) = T{{\dot \rho }_i},{\rho'}(\tau = 1) = T\dot \rho \hfill \\ \theta (\tau = 0) = {\theta _i},\theta (\tau = 1) = {\theta _f} \hfill \\ {\theta'}(\tau = 0) = T{{\dot \theta }_i},{\theta'}(\tau = 1) = T{{\dot \theta }_f} \hfill \\ {\text{z}}(\tau = 0) = {z_i}{\text{,z}}(\tau = 1) = {z_f} \hfill \\ {{\text{z}}'}(\tau = 0) = T{{\dot z}_i},{{\text{z}}'}(\tau = 1) = T{{\dot z}_f} \hfill \end{gathered} \right. $$ (7) 其中：T是电动帆的总飞行时间，下标“i”和“f”分别表示电动帆转移轨迹的初始条件和最终条件；0≤
$\tau {\text{ = }}t/T$ ≤1是无量纲化时间，t是电动帆飞行的瞬时时间。对于逃离太阳系问题，从数学的观点来看，这个问题相当于规定了电动帆在逃离时刻的机械能[15]为
$$ {\varepsilon _f} \triangleq \frac{{v_f^2}}{2} - \frac{\mu }{{{s_f}}} = \frac{{V_\infty ^2}}{2} $$ (8) 其中：
${v_f}$ 为逃离时电动帆的飞行速度；${V_\infty }$ 为电动帆的双曲线过渡速度。当${V_\infty }{\text{ = 0}}$ 时，可认为电动帆逃离太阳系。 -
在贝塞尔方法中[14]，
$\rho $ 、$\theta $ 、z在时间域中展开。为了避免重复，这里只介绍坐标$\rho $ 的处理过程（$\theta $ 、$z$ 可以近似处理）$$\left\{ \begin{gathered} \rho (\tau ) = \sum\limits_{j = 0}^{{n_\rho }} {{B_{\rho ,j}}(\tau )} {P_{\rho ,j}} \hfill \\ \rho '(\tau ) = \sum\limits_{j = 0}^{{n_\rho }} {B_{\rho ,j}'(\tau )} {P_{\rho ,j}} \hfill \\ \rho ''(\tau ) = \sum\limits_{j = 0}^{{n_\rho }} {B_{\rho ,j}{''}(\tau )} {P_{\rho ,j}} \hfill \\ \end{gathered} \right.$$ (9) 其中：
${n_\rho }$ 是贝塞尔曲线函数的阶数；${P_{\rho ,j}}$ 是未知贝塞尔系数，${B_{\rho ,j}}(\tau )$ 是状态变量的贝塞尔基函数，$\rho '$ 和$\rho ''$ 分别为$\rho $ 关于$\tau $ 的一阶和二阶导数。$$ {B_{\rho ,j}}(\tau ) = \frac{{{n_\rho }!}}{{j!({n_\rho } - j)!}}{\tau ^j}{(1 - \tau )^{{n_\rho } - j}} ,\cdots \quad ,j \in [0,{n_\rho }] $$ (10) $B_{\rho ,j}'$ 和$B_{\rho ,j}{''}$ 分别为$B_{\rho ,j}^{}$ 关于$\tau $ 的一阶和二阶导数，其具体形式可以查阅参考文献[14]。当$\tau = 0$ 和$\tau = 1$ 时，由式（10）可以得到贝塞尔基函数在边界处的特征值为$$ \left\{ \begin{gathered} {B_{\rho ,j}}(\tau = 0) = \left\{ \begin{gathered} 1\;\;\;\;, \quad j = 0 \hfill \\ 0\;\;\;\;, \quad j \in [1,\,{n_\rho }] \hfill \\ \end{gathered} \right. \hfill \\ {B_{\rho ,j}}(\tau = 1) = \left\{ \begin{gathered} 0\;\;\;\;, \quad j \in [0,\,{n_\rho } - 1] \hfill \\ 1\;\;\;\;, \quad j = {n_\rho } \hfill \\ \end{gathered} \right. \hfill \\ B_{\rho ,j}'(\tau = 0) = \left\{ \begin{gathered} - {n_\rho }\;\;\;\;, \quad j = 0 \hfill \\ {n_\rho }\quad , \;\;j = 1 \hfill \\ 0\;\;\;\;, \quad j \in [2,\,{n_\rho }] \hfill \\ \end{gathered} \right. \hfill \\ B_{\rho ,j}'(\tau = 1) = \left\{ \begin{gathered} 0\;\;\;\;, \;\;\;\; j \in [0,\,{n_\rho } - 2] \hfill \\ - {n_\rho }\;\;, \quad j = {n_\rho } - 1 \hfill \\ {n_\rho }\;\;\;\;,\quad \;\;j \in {n_\rho } \hfill \\ \end{gathered} \right. \hfill \\ \end{gathered}\right. $$ (11) 由式（7）、（9）、（10）和（11）可以得到方程
$$\left\{ \begin{array}{cc}\begin{array}{l}{\rho }_{i}=\rho (\tau =0)={P}_{\rho ,0}\\ {\rho }_{f}=\rho (\tau =1)={P}_{\rho ,{n}_{\rho }}\\ T{\dot{\rho }}_{i}=\rho \text{'}(\tau =0)={n}_{\rho }({P}_{\rho ,1}-{P}_{\rho ,0})\\ T{\dot{\rho }}_{f}=\rho \text{'}(\tau =1)={n}_{\rho }({P}_{\rho ,{n}_{\rho }}-{P}_{\rho ,{n}_{\rho }-1})\end{array}& \end{array} \right.$$ (12) 于是，通过方程（12）可以获得4个贝塞尔系数
${P_{\rho ,0}}$ 、${P_{\rho ,1}}$ 、${P_{\rho ,{n_\rho } - 1}}$ 和${P_{\rho ,{n_\rho }}}$ 为$$ \left\{ \begin{gathered} {P_{\rho ,0}} = {\rho _i} \hfill \\ {P_{\rho ,1}} = {\rho _i} + T{{\dot \rho }_i}/{n_\rho } \hfill \\ {P_{\rho ,{n_\rho } - 1}} = {\rho _f} - T{{\dot \rho }_f}/n \hfill \\ {P_{\rho ,{n_\rho }}} = {\rho _f} \hfill \end{gathered} \right. $$ (13) 其中：
${\rho _i}$ 和${\rho _f}$ 分别为电动帆出发位置和到达位置的$\rho $ 值。为了沿飞行轨迹提供有效的约束，需要在一些离散点对飞行轨迹进行评估。飞行轨迹的离散点采用勒让德–高斯分布，即离散点分布选取为
$m$ 次勒让德多项式的根为$$ {\tau _1}{\text{ = }}0 < {\tau _2} < \cdots < {\tau _{m - 1}} < {\tau _m} = 1 $$ (14) 因此，3个坐标(
$\rho ,\theta ,z$ )可以表示为矩阵形式。以坐标$\rho $ 为例，$\rho $ 及其相关导数可以表示为$$ \left\{\begin{aligned} &{[\rho ]_{m \times 1}} = {[{B_\rho }]_{m\, \times \,({n_\rho } + 1)}}{[{P_\rho }]_{({n_\rho } + 1){\kern 1pt} \times {\kern 1pt} 1}} \hfill \\ & {[\rho ']_{m \times 1}} = {[B_\rho ']_{m\, \times \,({n_\rho } + 1)}}{[{P_\rho }]_{({n_\rho } + 1){\kern 1pt} \times {\kern 1pt} 1}} \hfill \\ & {[\rho '']_{m \times 1}} = {[B_\rho {''}]_{m\, \times \,({n_\rho } + 1)}}{[{P_\rho }]_{({n_\rho } + 1){\kern 1pt} \times {\kern 1pt} 1}} \hfill \end{aligned}\right. $$ (15) 其中：
$$ \begin{aligned}[b] {[{P_\rho }]_{({n_\rho } + 1){\kern 1pt} \times {\kern 1pt} 1}} = & [{P_{\rho ,0}}\quad {P_{\rho ,1}}\quad [{X_\rho }]_{({n_\rho } - 3){\kern 1pt} \times {\kern 1pt} 1}^{\rm{T}} \hfill \\ & {P_{\rho ,{n_\rho } - 1}}\quad {P_{\rho ,{n_\rho }}}{]^{\rm{T}}} \hfill \end{aligned} $$ (16) $$ [{X_\rho }]_{({n_\rho } - 3){\kern 1pt} \times {\kern 1pt} 1}^{} = {[{P_{\rho ,2}}\;\;\; \cdots \;\;\;{P_{\rho ,{n_\rho } - 2}}]^{\rm{T}}} $$ (17) 其中：
$ {P_{\rho ,0}},\,{P_{\rho ,1}},\,{P_{\rho ,{n_\rho } - 1}},\,{P_{\rho ,{n_\rho }}} $ 均可以由式（13）获得，而$[{X_\rho }]_{({n_\rho } - 3){\kern 1pt} \times {\kern 1pt} 1}^{\rm{T}}$ 是需要优化的未知贝塞尔系数。由这些坐标的矩阵形式可以将推力加速度的各个分量的方程写为
$$ \left\{ \begin{gathered} {[{a_\rho }]_{m \times 1}} = {a_\rho }\left( {{{\left[ \rho \right]}_{m \times 1}},{{\left[ z \right]}_{m \times 1}},{{\left[ {\theta '} \right]}_{m \times 1}},{{\left[ {\rho ''} \right]}_{m \times 1}}} \right) \hfill \\ {\left[ {{a_\theta }} \right]_{m \times 1}} = {a_\theta }\left( {{{\left[ \rho \right]}_{m \times 1}},{{\left[ {\rho '} \right]}_{m \times 1}},{{\left[ {\theta '} \right]}_{m \times 1}},{{\left[ {\theta ''} \right]}_{m \times 1}}} \right) \hfill \\ {\left[ {{a_z}} \right]_{m \times 1}} = {a_z}\left( {{{\left[ \rho \right]}_{m \times 1}},{{\left[ z \right]}_{m \times 1}},{{\left[ {z''} \right]}_{m \times 1}}} \right) \hfill \end{gathered} \right. $$ (18) 因此，推力系数也可以使用矩阵形式为
$$ {[\kappa ]}_{m\times 1}=\frac{[s](3(\mathrm{sin}[\phi ]·[{a}_{\rho }]+\mathrm{cos}[\phi ]·[{a}_{z}])-\sqrt{[D]})}{2{a}_{c}{r}_{\oplus }} $$ (19) 其中，
$$ \begin{aligned}[b] {[D]_{m \times 1}} &= {\left( {\sin [\varphi ] · \left[ {{a_\rho }} \right] + \cos [\varphi ] · \left[ {{a_z}} \right]} \right)^2}-\\ & 8{\left[ {{a_\theta }} \right]^2} - 8{\left( {\cos [\varphi ] · \left[ {{a_\rho }} \right] - \sin [\varphi ] · \left[ {{a_z}} \right]} \right)^2} \geqslant 0 \end{aligned}$$ (20) 于是，可以将电动帆转移过程中的所需求解的连续轨迹优化问题转化为非线性规划问题，即
$$ \begin{array}{*{20}{c}} {\mathop {\min }\limits_{\left[ {{X_\rho }} \right],\;\left[ {{X_\theta }} \right],\;\left[ {{X_z}} \right],\;T} \;\;\;T} \\ {{\rm s.t.}\;\;\;{{[\kappa ]}_{m \times 1}} \geqslant 0,\;{{[\kappa ]}_{m \times 1}} \leqslant 1,\;{{[D]}_{m \times 1}} \geqslant 0,\;{{[{a_{o{\kern 1pt} z}}]}_{m \times 1}} \geqslant 0} \end{array} $$ (21) 其中：
$[{X_\rho }]、\;[{X_\theta }]、\;[{X_z}]$ 分别是3个坐标的未知贝塞尔系数。推力系数
$\kappa \in [0,1]$ ，因此可以由式（6）、（19）和（20）获得式（21）的4个约束条件。在研究的问题中，飞行时间在迭代计算过程中是不断进行优化的，电动帆的飞行时间T既是优化变量，又是性能指标。通过对式（9）~（13）的分析可以看到，当
${n_\rho } = {n_\theta } = {n_z} = 3$ 时，贝塞尔曲线函数的所有变量均可以由边界条件（7）计算获得。可使用三阶贝塞尔曲线函数（${n_\rho } = {n_\theta } = {n_z} = 3$ ）对三个坐标($\rho ,\;\theta ,\;z$ )进行初始化。仍以坐标$\rho $ 为例，其初始化结果为$$\begin{aligned}[b] \rho (\tau ) =& {(1 - \tau )^3}{P_{\rho ,0}} + 3\tau {(1 - \tau )^2}{P_{\rho ,1}}+\\ & 3{\tau ^2}(1 - \tau ){P_{\rho ,2}} + {\tau ^3}{P_{\rho ,3}} \end{aligned} $$ (22) 其中：4个贝塞尔函数系数
${P_{\rho ,0}}$ 、${P_{\rho ,1}}$ 、${P_{\rho ,2}}$ 和${P_{\rho ,3}}$ 均可以由式（7）和（13）来确定。在三维轨迹设计中，贝塞尔曲线函数的阶数分别为
$ {n}_{\rho }=\text{8}，{n}_{\theta }\text{=8}，{n}_{z}=8 $ ，离散点数m = 120。${a_c}$ = 2×10−3 m/s2，发射日期为2 459 022.5。FFS方法的参数设置与贝塞尔方法一致。GPM的轨迹优化中，勒让德高斯点的数目为120。采用贝塞尔方法对初始轨迹进行设计后，通过GPM得到进一步优化的轨迹。图3是使用贝塞尔方法、FFS方法和GPM方法获得电动帆逃离太阳系优化转移轨迹的结果。图4是使用贝塞尔方法、FFS方法和GPM方法获得电动帆飞行位置随时间的变化结果，图4（a）为电动帆在x方向的位置，图4（b）为电动帆在y方向的位置，图4（c）为电动帆在z方向的位置。图5是使用贝塞尔方法、FFS方法和GPM方法获得电动帆飞行速度随时间的变化结果，图5（a）为电动帆在x方向的飞行速度，图5（b）为电动帆在y方向的飞行速度，图5（c）为电动帆在z方向的飞行速度。图6是在轨道坐标系中使用贝塞尔方法获得的总的推进加速度及其分量随时间的变化曲线。图7是使用贝塞尔方法、FFS方法和GPM方法获得的电动帆俯仰角
${\alpha _n}$ 随时间变化的结果。
图 3 使用贝塞尔方法、FFS方法和GPM方法获得的逃离太阳系优化转移轨迹结果
Figure 3. The results of the trajectories of escape from the solar system obtained by the Bezier method，the FFS method and the GPM method
通过图3可以看到，3种方法获得的转移轨迹相近，电动帆在相同的空间位置达到可以逃离太阳系的飞行速度。由图4可以看到，3种方法获得的电动帆到达终点的坐标位置基本相同。由图7可以看到，3种方法获得的电动帆俯仰角
${\alpha _n}$ 均在[$0, \text{π} /2$ ]范围内，满足约束要求。
图 4 使用贝塞尔方法、FFS方法和GPM方法获得的电动帆飞行位置随时间的变化结果
Figure 4. The position of the electric sail obtained by the Bezier method，the FFS method and the GPM method
图 5 使用贝塞尔方法、FFS方法和GPM方法获得的电动帆飞行速度随时间的变化结果
Figure 5. The velocity of the electric sail obtained by the Bezier method，the FFS method and the GPM method
图 6 在轨道坐标系中，用贝塞尔方法获得的推进加速度随时间的变化
Figure 6. In orbit coordinate system，the propulsion acceleration obtained by the Bezier method
图 7 使用贝塞尔方法、FFS方法和GPM方法获得的电动帆俯仰角
${\alpha _n}$ 结果Figure 7. The results of the pitch angle
${\alpha _n}$ of the electric sail obtained by the Bezier method，the FFS method and the GPM method -
利用贝塞尔和FFS方法设计了电动帆逃离太阳系的转移轨迹，并将贝塞尔方法获得的结果作为初值代入GPM得到进一步优化轨迹。3种方法得到的飞行时间和计算时间如表1所示。在相同条件下，贝塞尔方法的计算时间10.642 8 s，电动帆逃离太阳系所需的飞行时间0.741 4 a；FFS方法使用18.199 3 s的计算时间获得转移轨迹，转移过程电动帆需要飞行0.760 2 a；GPM获得的飞行时间0.7091年，计算时间694.506 7 s。
表 1 贝塞尔方法、FFS方法和GPM的结果比较
Table 1. Comparison of the Bezier method，the FFS method and the GPM
方法 飞行时间/a 计算时间/s 贝塞尔 0.7414 10.6428 FFS 0.7602 18.1993 GPM 0.7091 694.5067 仿真结果表明，和FFS方法相比，贝塞尔方法可以使用更短的计算时间获得性能指标更好的转移轨迹。通过将贝塞尔方法获得的结果作为GPM的初始解，贝塞尔方法和GPM计算得到的飞行时间仅相差约4.56%，而贝塞尔方法的计算时间仅约为GPM计算时间的1.53%。因此贝塞尔方法在进行电动帆逃离太阳系转移轨迹的快速计算上具有良好的性能，非常适合在初始任务设计阶段，求解大量仿真算例的情形。
-
本文提出了一种利用贝塞尔方法设计电动帆逃离太阳系的三维转移轨迹的方法，且对轨迹终点位置和速度不做具体的约束，该方法假定电动帆的飞行轨迹具有贝塞尔曲线函数的形式。为了验证贝塞尔方法的先进性，将其获得的结果作为GPM的初始值进行后续计算，并将结果进行比较。数值结果表明，贝塞尔方法可以在较短的时间内设计出一个合理的三维初始轨迹。这对于在初步任务设计阶段快速评估大量电动帆飞行场景具有重要意义。
Initial Trajectory Design of the Electric Sail Escaping from Solar System Based on the Bezier Curve Method
-
摘要: 探测太阳系边际需要寻找高效的推进系统和设计合理逃离太阳系的轨迹。提出了利用贝塞尔形函数方法快速生成电动帆逃离太阳系的三维轨迹，为更精确的轨迹优化算法提供合适的初始解。通过将贝塞尔形函数方法获得的结果作为高斯伪谱法（Gauss Pseudospectral method，GPM）的初始解进行计算，并对贝塞尔形函数方法获得的结果对直接解算器初始解的适用性进行了评估。仿真结果表明，贝塞尔形函数方法可以在短时间内为直接优化求解器设计出合理的电动帆逃离太阳系的三维初始轨迹。对于在太阳系边际探测初步任务设计阶段快速评估大量电动帆飞行场景具有参考意义。Abstract: The detection of the solar system boundary has gradually attracted people's attention. To explore the solar system boundary, we need to find an efficient propulsion system and design a reasonable trajectory to escape from the solar system. Therefore, this paper proposes to use the Bezier shape-based method to quickly generate the three-dimensional trajectory of the electric sail escaping from the solar system, which provides an appropriate initial solution for a more accurate trajectory optimization algorithm. By taking the results obtained by the Bezier shape-based method as the initial solution of Gauss pseudospectral method（GPM）, the applicability of the results obtained by the Bezier shape-based method to the initial solution of the direct solver is evaluated. The simulation results show that the Bezier shape-based method can design a reasonable three-dimensional initial trajectory of the electric sail escaping from the solar system for the direct optimization solver in a short time. This is of great significance for the rapid evaluation of a large number of electric sail flight scenarios in the preliminary mission design stage.Highlights
● A fast algorithm is proposed to generate the trajectories of escaping solar system by using the shape-based method. ● The transfer trajectory of the electric sail escaping from the solar system is studied. ● The applicability of the results obtained by the Bezier shape-based method to the initial solution of the direct solver is evaluated. -
图 2 参考坐标系和推进加速度的特征角
Fig. 2 Reference coordinate system and characteristic angle of propulsion acceleration
图 3 使用贝塞尔方法、FFS方法和GPM方法获得的逃离太阳系优化转移轨迹结果
Fig. 3 The results of the trajectories of escape from the solar system obtained by the Bezier method，the FFS method and the GPM method
图 4 使用贝塞尔方法、FFS方法和GPM方法获得的电动帆飞行位置随时间的变化结果
Fig. 4 The position of the electric sail obtained by the Bezier method，the FFS method and the GPM method
图 5 使用贝塞尔方法、FFS方法和GPM方法获得的电动帆飞行速度随时间的变化结果
Fig. 5 The velocity of the electric sail obtained by the Bezier method，the FFS method and the GPM method
图 6 在轨道坐标系中，用贝塞尔方法获得的推进加速度随时间的变化
Fig. 6 In orbit coordinate system，the propulsion acceleration obtained by the Bezier method
图 7 使用贝塞尔方法、FFS方法和GPM方法获得的电动帆俯仰角
${\alpha _n}$ 结果Fig. 7 The results of the pitch angle
${\alpha _n}$ of the electric sail obtained by the Bezier method，the FFS method and the GPM method表 1 贝塞尔方法、FFS方法和GPM的结果比较
Table 1 Comparison of the Bezier method，the FFS method and the GPM
方法 飞行时间/a 计算时间/s 贝塞尔 0.7414 10.6428 FFS 0.7602 18.1993 GPM 0.7091 694.5067 
下载: 导出CSV
-
[1] BURLAGA L F,BURLAGA L F,NESS N F,et al. Intermittency and q-Gaussian distributions in the magnetic field of the very local interstellar medium(VLISM)Observed by Voyager 1 and Voyager 2[J]. The Astrophysical Journal Letters,2020,901(1):1-6. doi: 10.3847/1538-4357/abaa48 [2] GURNETT D A,KURTH W S. Plasma densities near and beyond the heliopause from the Voyager 1 and 2 plasma wave instruments[J]. Nature Astronomy,2019,3(11):1024-1028. doi: 10.1038/s41550-019-0918-5 [3] GURNETT D A,KURTH W S,STONE E C,et al. A foreshock model for interstellar shocks of Solar origin:Voyager 1 and 2 observations[J]. The Astronomical Journal,2021,161(1):1-11. [4] 张恒,文浩,罗操群. 旋转电动帆推进性能指标分析[J]. 动力学与控制学报,2019,17(3):281-287. doi: 10.6052/1672-6553-2019-018ZHANG H,WEN H,LUO C Q. Analysis of rotating electric solar sail for spacecraft propulsion[J]. Journal of Dynamics and Control,2019,17(3):281-287. doi: 10.6052/1672-6553-2019-018 [5] 霍明英,齐乃明,刘宇飞,等. 电动帆改进推力模型及深空探测性能分析[J]. 动力学与控制学报,2018,16(2):151-156. doi: 10.6052/1672-6553-2018-003HUO M Y,QI N M,LIU Y F,et al. Refined thrust model of electric sail and analysis of deep space detection performance[J]. Journal of Dynamics and Control,2018,16(2):151-156. doi: 10.6052/1672-6553-2018-003 [6] PETROPOULOS A,LONGUSKI J. Shape-based algorithm for the automated design of low-thrust gravity assist trajectories[J]. Journal of Spacecraft and Rockets,2004,41(5):787-796. doi: 10.2514/1.13095 [7] PASCALE P,VASILE M. Preliminary design of low-thrust multiple gravity-assist trajectories[J]. Journal of Spacecraft and Rockets,2006,43(5):1069-1076. [8] WALL B,CONWAY B. Shape-based approach to low-thrust rendezvous trajectory design[J]. Journal of Guidance,Control,and Dynamics,2009,32(1):95-101. [9] XIE C,ZHANG G,ZHANG Y. Simple shaping approximation for low-thrust trajectories between coplanar elliptical orbits[J]. Journal of Guidance,Control,and Dynamics,2015,38(12):2448-2455. [10] PELONI A,DACHWALD B,CERIOTTI M. Multiple near-Earth asteroid rendezvous mission:solar-sailing options[J]. Advances in Space Research,2018,62(8):2084-2098. doi: 10.1016/j.asr.2017.10.017 [11] 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 [12] TAHERI E,ABDELKHALIK O. Fast initial trajectory design for low thrust restricted-three-body problems[J]. Journal of Guidance,Control,and Dynamics,2015,38(11):2146-2160. [13] TAHERI E,ABDELKHALIK O. Initial three-dimensional low-thrust trajectory design[J]. Advances in Space Research,2016(57):889-903. [14] HUO M Y,MENGALI G,QUARTA A,et al. Electric sail trajectory design with bezier curve-based shaping approach[J]. Aerospace Science and Technology,2019(88):126-135. [15] QUARTA A A,MENGALI G. Electric sail mission analysis for outer solar system exploration[J]. Journal of Guidance,Control,and Dynamics,2010,33(3):740-755. -
点击查看大图
计量
- 文章访问数: 575
- HTML全文浏览量: 229
- PDF下载量: 64
- 被引次数: 0
