电动帆改进推力模型及深空探测性能分析*

霍明英 齐乃明 刘宇飞 曹世磊 叶炎茂

(1.哈尔滨工业大学 航天工程系,哈尔滨 150001) (2.中国航天科技集团公司钱学森空间技术实验室，北京 100094)

引言

电动帆是由芬兰学者Janhunen于2004年提出的一种无质损的推进方式[1].与太阳帆不同的是,电动帆的动力来源不是太阳光压,而是太阳风带电粒子的动能冲力,其飞行原理示意图如图1所示.电动帆由很多根长而细的柔性金属链所组成,这些金属链通过飞行器自旋展开.空间飞行器上的太阳能电子枪向外喷射电子,使金属链始终保持在高度的正电位,这些带电的金属链会排斥太阳风正电粒子,从而利用太阳风的动能冲力推动空间飞行器驶向目标方向.由于电动帆能够利用太阳风的动能冲力飞行而不需要消耗推进剂,因此电动帆非常适用于长期的空间飞行任务[2,3].

欧盟第七框架协议资助电动帆项目组数百万欧元开展理论研究及原型样机研制.美国航空航天局(NASA)马歇尔航天中心也正在实施HERTS(太阳风静电高速运输系统)计划,计划将电动帆应用于太阳风层顶的探测.电动帆的飞行原理已经经过空间飞行测试(ESTCube-1卫星),其原型样机的测试将在未来的几年内完成(Aalto-1卫星[4,5],如图2所示).目前,国内已有航天科技集团五院502所[6,7]、北京理工大学[8]和哈尔滨工业大学[9,11]等正在开展相关的研究工作.

近年来,电动帆推力矢量数学模型方面的研究得到了相关研究团体的关注.在电动帆的传统推力模型中,忽略了电动帆姿态对推力标量的影响,并假设电动帆的推进角为入射角的一半[12].针对此问题,本文推导得出一种解析形式的改进推力模型,并与最新的多项式拟合改进推力模型[13]进行了对比.另外,为了重新评估电动帆在深空探测中的性能,以地球至火星飞行任务为算例,分别采用传统推力模型和解析改进推力模型进行了电动帆轨迹优化仿真,并以不同的特征加速度情况进行对比分析.

图1 电动帆原理示意图Fig.1 Conceptual sketch of an electric sail

图2 电动帆原型样机Aalto-1试验卫星Fig.2 Aalto-1 prototype test satellite of electric sail

1 电动帆推力模型及对比

1.1 坐标系及角度定义

本文所涉及到的参考系共三个,分别是日心黄道参考系、轨道参考系和体参考系,如图3所示.日心黄道参考系os-xsyszs的原点为太阳中心,正xs轴指向历元J2000.0时刻平春分点方向,正zs轴垂直于J2000.0时刻黄道面并指向黄道北极方向,ys轴与xs轴和zs构成右手系.轨道参考系oo-xoyozo的原点位于电动帆航天器质心,正zo轴为太阳-电动帆航天器的矢量方向,yo轴与zo轴和日心黄道参考系中的zs轴垂直,方向指向飞行运动方向,xo轴与yo轴和zo构成右手系.电动帆体参考系ob-xbybzb的原点位于电动帆航天器质心,正xb轴沿着1号金属链方向,正zb轴垂直于电动帆工作面,并指向自旋角速度矢量方向,yb与xb和zb构成右手系.

如图3所示,太阳风入射角αn为太阳风入射方向zo与电动帆回转体轴zb之间的夹角,推进锥角α为推进加速度矢量a与电动帆回转体轴zb之间的夹角,推进钟角δ为推进加速度矢量在ooxoyo平面内分量与xo轴之间的夹角.电动帆航天器在轨道参考系下的姿态可以通过三个角进行描述,即φ,θ和ψ三个姿态角,由轨道参考系转换至体参考系的姿态转换顺序为x(φ)-y(θ)-z(ψ),姿态旋转矩阵为Abo(φ,θ,ψ).

图3 参考系及角度示意图Fig.3 Reference frame and angles

1.2 传统推力模型

在意大利比萨大学Mengali教授提出的电动帆传统推力模型[12]中,忽略了电动帆姿态对推力标量的影响,即假设在电动帆工作面与太阳风入射方向不垂直时,推力幅值不变.实际上电动帆与太阳帆一样,不只推力的方向由帆体姿态所决定,推力大小也一定程度上取决于帆体相对太阳风速度方向的姿态,原因是当电动帆帆体平面相对太阳风粒子运动方向产生角度变化时,太阳风粒子与带电金属量的动量交互效率将发生改变,从而最终影响推力的大小.另外,传统推力模型中还假设推进锥角α近似等于太阳风入射角αn的一半,即α≃αn/2.传统推力模型中电动帆推力矢量a在轨道参考系oo-xoyozo下可写为:

(1)

其中a⊕为电动帆的特征加速度,即电动帆距离太阳r⊕=1AU处所能产生的最大加速度值；κ∈[0,1]为电动帆推力开关系数,可以通过电子枪调整金属链电压来调整电动帆整体的推力.

为了表征太阳风入射角对推力大小的影响,引入无量纲加速度γ,其表达式为:

γ

(2)

由等式(1)可知,电动帆传统推力模型中无量纲加速度γ1,即假设推力大小不受姿态的影响.

1.3 多项式改进推力模型

为了讨论太阳风入射角αn对推进锥角α和推力大小γ的影响,日本学者Yamaguchi和Yamakawa[13]通过部分实验和数学仿真的方法,以多项式拟合的形式得出了太阳风入射角αn与推进锥角α的关系式(等式3)和太阳风入射角αn与推力大小γ的关系式(等式4).意大利学者Quarta[14]在此模型基础上采用间接优化方法得出了多组优化轨迹:

(3)

(4)

其中bk(k=0,1,…,6)和ck(k=0,1,…,6)是多项式拟合系数,如表1所示.

表1 多项式拟合系数Table 1 Polynomial fitting coefficient

1.4 解析改进推力模型

在本文作者前期的研究[15]中,以电动帆单根带电金属链推力模型为基础,采用有限傅里叶级数加和方法,推导得出了考虑电动帆姿态影响的推力矢量模型,在轨道参考系oo-xoyozo下的表达式为:

(5)

其中φ和θ为体参考系相对轨道系的姿态角(定义见2.1节).由姿态角角度定义可知,太阳风入射角αn的余弦、推进钟角δ的正弦和余弦可以写作如下形式:

=cosφcosθ

(6)

=sinθ/sinαn

(7)

=-cosθsinφ/sinαn

(8)

将等式(6)-等式(8)带入等式(5),可得到以太阳风入射角αn和推进钟角δ描述的推力模型:

(9)

由图3可知,推进锥角α的余弦可写成如下形式:

(10)

化简等式(10)后,可得到推进锥角α与太阳风入射角αn的解析关系式:

(11)

参考γ的定义(等式(2)),可得到无量纲加速度γ与太阳风入射角αn的解析关系式:

(12)

1.5 推力模型对比

推进锥角α与太阳风入射角αn的关系如图4所示,上述的三个推力模型在αn∈[0°,20°]区间内能够很好地重合,说明传统推力模型的推进锥角假设α≃αn/2在小太阳风入射角的情况下是适用的.但随着太阳风入射角的增加,Yamaguchi得到的多项式改进模型和本文得到的解析改进模型都表现出了非常显著的非线性变化行为.本文所得到的推力模型在αn=54.75°时,推进锥角α达到最大值19.47°.

无量纲加速度γ与太阳风入射角αn的关系如图5所示,在两种改进的推力模型中无量纲加速度γ随着太阳风入射角αn的增大而单调减小,并在αn=90°时达到最小γ=0.5.另外,本文还对无量纲加速度径向分量(zo轴分量)和切向分量(ooxoyo平面内分量)进行了对比,如图6和图7所示.两种改进推力模型的无量纲径向加速度γr和无量纲切向加速度γt均小于传统模型的估计值.另外,由图4～图7可知,本文的解析改进推力模型与多项式改进推力模型数值结果很接近,但本文的解析模型形式更简单.

图4 推进锥角α与太阳风入射角αn的关系Fig.4 The relationship between the thrust cone angle α and inclined angle αn

图5 无量纲加速度γ与太阳风入射角αn的关系Fig.5 The relationship between dimensionless acceleration γ and inclined angle αn

图6 径向无量纲加速度γr与太阳风入射角αn的关系Fig.6 The relationship between radial dimensionless acceleration γr and inclined angle αn

2 电动帆地球至火星轨迹优化及对比分析

为了重新评估电动帆在深空探测中的性能,本文以地球至火星的飞行任务为算例,分别采用传统推力模型和解析改进推力模型进行电动帆轨迹优化仿真.在电动帆航天器地球至火星轨迹优化问题中,优化性能指标为飞行时间最优；在边界约束方面,假设电动帆航天器在初始时刻位于地球逃逸抛物线轨迹上,且逃逸剩余能量为零,并考虑了真实星历约束.所采用的方法为一种结合遗传算法和高斯伪谱法的混合优化方法[9],能够在无任何初值猜测的情况下完成电动帆航天器飞行轨迹的优化,避免了单纯高斯伪谱法初值赋值繁琐的问题.

图7 切向无量纲加速度γt与太阳风入射角αn的关系Fig.7 The relationship between tangential dimensionless acceleration γt and inclined angle αn

图8 不同特征加速度下的地球-火星过渡时间Fig.8 Earth-Mars transition time with different characteristic acceleration

以不同特征加速度下电动帆自地球至火星的过渡轨迹进行优化仿真,考虑的特征加速度范围为0.5～1.1mm/s2,每个仿真之间的特征加速度间隔为0.1mm/s2,基于传统推力模型和解析改进推力模型优化得出的飞行时间与特推力模型进行优化仿真,电动帆航天器特征加速度越大,完成过渡所需要的飞行时间均越短.这说明电动帆的推进能力越强,完成任务所需的飞行时间也就越短.

而两者之间不同的是,基于电动帆改进推力模型优化得出的飞行时间要大于基于传统模型优化得出的飞行时间.出现这种现象的原因一共有两点:①传统推力模型中忽略了姿态改变对推力加速度大小的影响,假定即使帆体平面与太阳风粒子速度方向不垂直时,动量交换效率依然不会减弱；②传统模型中估计得出的最大推进锥角αmax假定为35°,比改进推力模型计算得出的最大推进锥角αmax=19.47°大,即传统推力模型高估了电动帆的切向推进能力.

3 结论

本文针对电动帆传统推力模型中忽略姿态对推力幅值影响的问题,在原有基础上推导得出了一种解析形式的电动帆改进推力模型,并与最新提出的多项式拟合改进推力模型进行了对比.对比结果显示:①传统推力模型的推进锥角假设在小太阳风入射角的情况下是适用的,但随着太阳风入射角的增大,估计偏差会越来越大；②两种改进推力模型数值结果很接近,但本文的解析改进推力模型形式更简单.另外,为了重新评估电动帆在深空探测中的性能,本文以地球至火星的飞行任务为算例,分别采用传统推力模型和解析改进推力模型进行了电动帆轨迹优化仿真.仿真结果显示,在相同特征加速度情况下采用改进解析推力模型完成任务所需时间大于采用传统推力模型所用时间.分析结果认为,出现上述现象的原因在于传统推力模型中忽略了姿态改变对推力加速度大小的影响,且高估了电动帆所能产生的最大推进角.

1Janhunen P. Electric sail for spacecraft propulsion.JournalofPropulsionandPower, 2004,20(4):763～764

2Janhunen P, Sandroos A. Simulation study of solar wind push on a charged wire: basis of solar wind electric sail propulsion.AnnalesGeophysicae, 2007,25(3):755～767

3Kiprich S, Kurppa R, Janhunen P. Wire-to-wire bonding of lm-diameter aluminum wires for the Electric Solar Wind Sail.MicroelectronicEngineering, 2011,88:3267～3269

4Kestilä A, Tikka T, Peitso P, et al. Aalto-1 nanosatellite-technical description and mission objectives.GeoscientificInstrumentation,MethodsandDataSystems, 2013,2(1):121～130

5Khurshid O, Tikka T, Praks J, et al. Accommodating the plasma brake experiment on-board the Aalto-1 satellite.ProceedingsoftheNationalAcademyofSciencesoftheUnitedStatesofAmerica, 2014,64(2):258～266

6Chen M, Xia G, Wei Y, et al. Charateristics and Stress Analysis of Sheath of Parallel Conducting Tethers for the Electric Sail.ActaPhysicaSinica, 2016,65(209601):1～11

7王昱,魏延明,李永等. 基于粒子群算法的电帆轨迹优化设计. 中国空间科学技术, 2015,3(3):26～34 (Wang Y, Wei Y , Li Y, et al. Bian Bingxiu. Trajectory Optimization of Electric Sail Based on Particle Swarm Algorithm.ChineseSpaceScienceandTechnology, 2015,3(3):26～34 (in Chinese))

8胡权,张景瑞. 光电帆航天器姿态动力学建模. 第十五届全国非线性振动暨第十二届全国非线性动力学和运动稳定性学术会议. 中国, 长沙:2015,05 (Hu Q, Zhang J R. Modeling of attitude dynamics of a photoelectric sail spacecraft. Fifteenth National nonlinear vibration and the Twelfth National Symposium on nonlinear dynamics and motion stability, China, Changsha:2015,05 (in Chinese))

9Huo M, Mengali G, Quarta A A. Optimal planetary rendezvous with an electric sail.AircraftEngineeringandAerospaceTechnology, 2016,88(4):515～522

10 霍明英,彭福军,赵钧等. 电动帆航天器谷神星探测任务轨迹优化.宇航学报, 2015(12):1363～1372 (Huo M Y, Peng F J, Zhao J,et al. Trajectory optimization for ceres exploration with an electric sail.JournalofAstronautics, 2015(12):1363～1372 (in Chinese))

11 齐乃明,霍明英,袁秋帆. 电动帆轨迹优化及其性能分析. 宇航学报, 2013,30(5):634～641 (Qi N M, Huo M Y, Yuan Q F. Trajectory optimization for ceres exploration with an electric sail.JournalofAstronautics, 2013,30(5):634～641 (in Chinese))

12 Mengali G, Quarta A A. Janhunen P. Electric sail performance analysis.JournalofSpacecraftandRocket, 2008,45(1):122～129

13 Yamaguchi K, Yamakawa H. Study on orbital maneuvers for electric sail with on-off thrust control.AerospaceandTechnologyofJapanese, 2013,12:79～88

14 Quarta A A, Mengali G. Minimum-time trajectories of electric sail with advanced thrust model.AerospaceScienceandTechnology, 2016,55:419～430

15 霍明英. 电动帆航天器动力学、控制及轨迹优化研究[博士学位论文]. 哈尔滨:哈尔滨工业大学, 2016 (Huo M Y. Research on Dynamics Control and Trajectory Optimization for Electric Sails[Ph.D Thesis]. Harbin:Harbin Institute of Technology, 2015 (in Chinese))