CN104050338A - 一种geo卫星小推力器推进剂消耗耦合分析方法 - Google Patents

Abstract

本发明涉及一种GEO卫星小推力器推进剂消耗耦合分析方法，属于卫星平台总体设计技术领域。本发明方法考虑了不同阶段推进剂消耗量与贮箱半径之间的耦合关系，推进剂的质心位置根据当前卫星所在阶段进行更新，明显改进了小推力器所提供的南北位保和东西位保推力大小和控制力矩的准确性。在推进剂消耗计算过程中，估算了采用小推力器卸载干扰力矩部分所需推进剂，使得推进剂消耗估算模型与实际工程任务相符，提高了模型准确性与实际工程能力。同时针对静止轨道阶段位保过程和转移轨道阶段推进剂在贮箱内位置不同，提出了不同质心计算方法，提高了位保推力和控制力矩计算模型精度，有益于提高计算卫星全寿命周期内推进剂消耗量计算精度。

Description

一种GEO卫星小推力器推进剂消耗耦合分析方法

技术领域

本发明涉及一种GEO卫星小推力器推进剂消耗耦合分析方法，属于卫星平台总体设计技术领域。 

背景技术

由于太阳和月球引力作用会造成地球同步轨道(Geostationary Orbit，GEO)卫星的南北倾角摄动，而地球形状非标准球形会造成GEO卫星的东西摄动。因此，在卫星寿命期内，GEO卫星通常要通过南北轨道保持(North-South Station Keeping，NSSK)和东西轨道保持(East-West Station Keeping，EWSK)来抵消轨道摄动。国内外众多专家学者针对GEO卫星的NSSK和EWSK组合策略提出了不同的策略。现在航天工业越来越多地选用电/离子推力器进行轨道位置保持，主要是由于它高效的推进剂效率和高比冲，并且可提供适当的低加速度。但是，考虑到电推力器中氙气和氦气的成本以及卫星自身的鲁棒性和适应能力，现有的主流卫星平台依然会选用化学小推力器承担部分卫星姿态控制和轨道保持的任务，或者作为电推力器的备份，如：SS/L-1300平台，LS3000平台，以及我国的DFH-4平台。 

通常，在卫星发射阶段，化学推进剂大约占卫星总重的50～60％，而到达同步轨道后，用于轨道保持和姿态控制的推进剂重量大约占卫星总重的30％左右。因此，在保证完成航天任务的前提下，减少推进剂携带量可以使得卫星携带更多的有效载荷，从而提高卫星性能。不同于电推力器可以调整推力方向，化学小推力器的构型与安装角度不可调整，并决定了推力器的控制效率。 

由于化学推力器所能提供的推力与控制力矩不仅仅与推力器的组合方式、安装位置和安装角度有关，还与当前卫星的质心位置有关。卫星的质心位置主要由三部分组成：有效载荷、平台和推进剂。其中，有效载荷和平台在构型一定的情况下，其质量和质心位置固定；而推进剂质量和质心位置在整个寿命期内都是不断变化。推进剂主要分为两部分，一部分用于轨道转移阶段主发动机提供转移轨道所需要的推力和小推力器姿态调整，另一部分用于静止轨道阶段小推力器进行位置保持和姿态调整。卫星变轨推进剂消耗量不仅与速度增量有 关，同时与卫星质量有关。另一方面，推进剂当前质心位置不仅与质量有关，同时还与贮箱尺寸相关。然而，小推力器在进行位置保持和姿态调整所需推进剂受当前卫星质心位置影响。因此，静止轨道阶段推进剂消耗量、转移轨道阶段推进剂消耗量和贮箱尺寸三者之间是相互影响相互耦合的过程。 

发明内容

本发明的目的是为了解决GEO卫星小推力器在转移轨道阶段和静止轨道阶段推进剂消耗量与贮箱尺寸三者之间的耦合关系问题，提出一种GEO卫星小推力器推进剂消耗耦合分析方法。 

本发明的目的是通过下述技术方案实现的。 

一种GEO卫星小推力器推进剂消耗耦合分析方法，具体实现步骤如下： 

步骤1：定义卫星机械坐标系和卫星本体坐标系，给出卫星的基本设计参数，并给定静止轨道阶段推进剂消耗量预估值(用符号M_p表示)和贮箱半径预估值(用符号R表示)；所述卫星的基本设计参数包括：卫星寿命、平台尺寸、南北位保周期、东西位保周期、主发动机推力、主发动机比冲、小推力器推力、小推力器比冲、推力器工作效率、贮箱与平台间隙、有效载荷质量、有效载荷质心Z方向高度、平台结构质量、推进剂密度。 

所述卫星机械坐标系的原点位于卫星与运载火箭机械分离面内，原点固联在卫星上，并与星箭对接环圆心重合，用符号X、Y和Z分别表示卫星机械坐标系的三个轴；其中，Z轴垂直于卫星与运载火箭的对接环面，其正方向从原点指向对地板，X轴的正方向从原点指向卫星东板，Y轴正方向从原点指向卫星南板。 

所述卫星本体坐标系的原点在卫星质心，用符号x_b、y_b和z_b分别表示卫星本体坐标系的三个轴；x_b、y_b和z_b分别平行于卫星机械坐标系X轴、Y轴和Z轴。 

步骤2：在步骤1操作的基础上，通过公式(1)和公式(2)分别计算卫星整个寿命期内南北位保次数与两次南北位保间隔期内东西位保次数。 

$n_{\mathrm{SN}}=\frac{T_{\mathrm{life}}}{T_{\mathrm{SN}}}---\left(1\right)$

其中，n_SN为卫星整个寿命期内南北位保次数；T_life为卫星寿命；T_SN为南北位保周期。 

$n_{\mathrm{EW}}=\frac{T_{\mathrm{SN}}}{T_{\mathrm{EW}}}---\left(2\right)$

其中，n_EW为两次南北位保间隔期内东西位保次数；T_EW为东西位保周期。 

步骤3：在步骤2操作的基础上，通过公式(3)得到液面半张角(用符号γ表示)，用公式(4)计算静止轨道期间第k次南北位置保持时贮箱内推进剂质心位置(用符号Y_mc表示)。其中，k为正整数，1≤k≤n_SN。k的初始值设置为1。 

$\frac{H}{R}(\mathrm{\γ}-\sin \mathrm{\γ}\cos \mathrm{\γ})+\frac{\mathrm{\π}}{3}{(1-\cos \mathrm{\γ})}^2(2+\cos \mathrm{\γ})=f\·(\mathrm{\π}\frac{H}{R}+\mathrm{\π}\frac{4}{3})---\left(3\right)$

其中，H为贮箱高度；R为贮箱半径；γ为推进剂液面半张角；f为充液比，即当前推进剂的体积与贮箱容积之比。 

$Y_{\mathrm{mc}}=\frac{\frac{H}{R}(\mathrm{\γ}-\sin \mathrm{\γ}\cos \mathrm{\γ})\·\frac{4}{3}\frac{{\sin }^3\mathrm{\γ}}{2\mathrm{\γ}-\sin 2\mathrm{\γ}}}{\frac{H}{R}(\mathrm{\γ}-\sin \mathrm{\γ}\cos \mathrm{\γ})+\frac{\mathrm{\π}}{3}{(1-\cos \mathrm{\γ})}^2(2+\cos \mathrm{\γ})}+\frac{\frac{\mathrm{\π}}{3}{(1-\cos \mathrm{\γ})}^2(2+\cos \mathrm{\γ})\frac{3}{4}\frac{{(1+\cos \mathrm{\γ})}^2}{2+\cos \mathrm{\γ}}}{\frac{H}{R}(\mathrm{\γ}-\sin \mathrm{\γ}\cos \mathrm{\γ})+\frac{\mathrm{\π}}{3}{(1-\cos \mathrm{\γ})}^2(2+\cos \mathrm{\γ})}---\left(4\right)$

其中，Y_mc为贮箱内推进剂质心位置。 

步骤4：在步骤3操作的基础上，通过公式(5)至公式(11)计算静止轨道期间第k次位置保持时位保推力矢量和X轴、Y轴和Z轴控制力矩。 

用符号表示单个推力器的位保推力矢量，推力矢量在卫星本体坐标系中方位角和仰角分别用符号和θ表示；用符号表示卫星质心在机械坐标系中的位置矢量，用符号表示推力器在卫星机械坐标系中的位置矢量， 用符号表示推力器在卫星本体坐标系中的位置矢量，

用符号M_x、M_y和M_z分别表示推力器产生的卫星机械坐标系X轴、Y轴和Z轴控制力矩；用符号α、β和μ分别表示推力矢量与卫星本体坐标系的x_b、y_b和z_b三个轴的夹角，则有如下关系： 

其中，和θ分别为卫星本体坐标系中方位角和仰角。 

μ＝arccos(sinθ)    (7) 

通过公式(8)可得到推力器产生的X轴控制力矩M_x。 

$\frac{M_x}{F\·D}=-\tilde{z}\cos \mathrm{\β}+\tilde{y}\cos \mathrm{\μ}---\left(8\right)$

其中，F和D分别为推力矢量和位置矢量的模，F和D为人为设定值；  $\tilde{z}=(z-z_c)/D;\tilde{y}=(y-y_c)/D.$

通过公式(9)可得到推力器产生的Y轴控制力矩M_y。 

$\frac{M_y}{F\·D}=\tilde{z}\cos \mathrm{\α}-\tilde{x}\cos \mathrm{\μ}---\left(9\right)$

其中， $\tilde{x}=(x-x_c)/D.$

通过公式(10)可得到推力器产生的Z轴控制力矩M_z。 

$\frac{M_z}{F\·D}=-\tilde{y}\cos \mathrm{\α}+\tilde{x}\cos \mathrm{\β}---\left(10\right)$

通过公式(11)可得到位保推力矢量

$\stackrel{\→}{F}=[F\cos \mathrm{\α},F\cos \mathrm{\β},F\cos \mathrm{\μ}]---\left(11\right)$

步骤5：在步骤2操作的基础上，通过公式(12)计算东西位保所需速度增量(用符号ΔV_EW表示)。 

$\left\{\begin{array}{c}\frac{\sin (0.5{\mathrm{\ω}}_s{T}_{\mathrm{EW}}+\mathrm{\δ})}{\sin \left(0.5{\mathrm{\ω}}_s{T}_{\mathrm{EW}}\right)}=\frac{e_r}{e_k}\\ \mathrm{\Δ}{V}_{\mathrm{EW}}=V_s{e}_k\sin \mathrm{\δ}\end{array}\right.---\left(12\right)$

其中，ω_s为地球公转角速度；T_EW为东西位保周期；δ为轨道偏移角度；e_r为偏心率摄动圆半径，e_r＝1.5a_s/(ω_sω_xa)；ω_x为地球自转角速度；a_s为太阳辐射压加速度；a为静止轨道半长轴；e_k为偏心率容许值；V_s为静止轨道速度。 

步骤6：在步骤5操作的基础上，静止轨道阶段，小推力器提供推力，通过公式(13)计算东西位保所需推进剂质量(用符号Δm表示)。 

$\mathrm{\Δm}=m_0[1-e^{-\mathrm{\Δ}{V}_{\mathrm{EW}}/\left(I_{\mathrm{sp}}\mathrm{\ηg}\right)}]---\left(13\right)$

其中，m₀为推力器开始工作时卫星的初始质量；I_sp为推力器比冲，η为推力器工作效率，I_sp和η为人为设定值；g为重力加速度，为常量。 

步骤7：在步骤5操作的基础上，通过公式(14)至公式(16)计算卸载东西位保干扰力矩所需的推进剂质量(用符号Δm_c表示)。 

$\left\{\begin{array}{c}{M}_{x}t=\mathrm{\λ}{M}_{\mathrm{Tx}}{t}_x\\ M_{y}t=\mathrm{\λ}{M}_{\mathrm{Ty}}{t}_y\\ M_{z}t=\mathrm{\λ}{M}_{\mathrm{Tz}}{t}_z\end{array}\right.---\left(14\right)$

其中，t为小推力器工作时间，可通过公式(15)得到；λ为卸载干扰力矩的效率，为人为设定值；M_Tx、M_Ty和M_Tz分别表示X轴、Y轴和Z轴干扰力矩的大小；t_x表示卸载X轴干扰力矩的工作时间；t_y表示卸载Y轴干扰力矩的工作时间；t_z表示卸载Z轴干扰力矩的工作时间。 

$F\·t/\left(I_{\mathrm{sp}}g\right)=m_0[1-e^{-{\mathrm{\ΔV}}_{\mathrm{EW}}/\left(I_{\mathrm{sp}}\mathrm{\ηg}\right)}]---\left(15\right)$

Δm_c＝FΔt/(I_spg)                     (16) 

其中，Δt＝t_x+t_y+t_z。 

步骤8：在步骤5操作的基础上，通过公式(17)至公式(18)计算静止轨道期间第k次南北位置保持时所需的速度增量(用符号ΔV_SN表示)。 

其中，Δi为轨道倾角变化量；Ω_mT为终了月球升交点黄经，Ω_m0为初始月球升交点黄经；Ω_m为月球升交点黄经，Ω_m＝259.183°-0.05295°·t′，t′为相对于1900年1月1日12时的儒略日。 

ΔV_SN＝ΔiπV_s/180                         (18) 

步骤9：在步骤8操作的基础上，计算静止轨道期间第k次南北位置保持时所需推进剂质量(用符号Δm′表示)和卸载干扰力矩所需推进剂质量(用符号Δm′_c表示)，此过程与步骤6计算东西位保时所需推进剂质量和步骤7计算卸载东西位保干扰力矩所需推进剂质量的方法相同。 

步骤10：在步骤9操作的基础上，计算当前贮箱内推进剂质量(用符号M′_p表示)，M′_p＝M_p-Δm-Δm_c-Δm′-Δm′_c，并使用M′_p更新贮箱内推进剂质量M_p，M_p＝M′_p。判断k＜n_SN是否成立，如果成立，则使k值增1，并重复执行步骤3至步骤9的操作；否则，计算静止轨道阶段推进剂消耗总量(用符号M_GEO表示)，并执行步骤11的操作； 

所述推进剂消耗总量M_GEO为静止轨道阶段第1次至第n_SN次位保时Δm、Δm_c、Δm′、Δm′_c四个量之和。 

步骤11：在步骤10操作的基础上，通过公式(19)计算转移轨道阶段第l次轨道转移时贮箱内推进剂质心高度(用符号Z_mc表示)。其中，l为正整数，1≤l≤l_trans。l的初始值设置为转移轨道次数l_trans，l_trans为人为设定值，3≤l_trans≤8。 

$Z_{\mathrm{mc}}=\frac{{8\mathrm{RH}}_0-{3H}_0^2}{4(3R-H_0)}---\left(19\right)$

其中，R为贮箱半径；H₀为推进剂液面高度。 

步骤12：在步骤11操作的基础上，轨道转移阶段，主发动机提供推力进行 轨道转移，在第l次轨道转移时所需的速度增量为已知量的情况下，通过公式(20)计算第l次轨道转移时所需推进剂质量(用符号Δm″表示)。 

${\mathrm{\Δm}}^{\′\′}=m_0[1-e^{-\mathrm{\ΔV}/\left(I_{\mathrm{sp}}\mathrm{\ηg}\right)}]---\left(20\right)$

其中，ΔV为在第l次轨道转移时所需的速度增量。 

步骤13：在步骤11操作的基础上，计算第l次轨道转移时X轴、Y轴和Z轴控制力矩，此过程与步骤4中计算静止轨道期间第k次位置保持时X轴、Y轴和Z轴控制力矩的方法相同。 

步骤14：在步骤13操作的基础上，计算第l次轨道转移时卸载干扰力矩所需推进剂质量(用符号Δm″_c表示)，此过程与步骤7计算卸载东西位保干扰力矩所需推进剂质量的方法相同。 

步骤15：在步骤14操作的基础上，判断l的值是否不为1，如果不为1，则使l值减1，并重复执行步骤11至步骤14的操作；否则，计算转移轨道阶段推进剂消耗总量(用符号M_trans表示)，并执行步骤16的操作； 

所述转移轨道阶段推进剂消耗总量M_trans为转移轨道阶段第1次至第l_trans次轨道转移时Δm″和Δm″_c之和。 

步骤16：在步骤15操作的基础上，根据静止轨道阶段推进剂总消耗量M_GEO和转移轨道阶段推进剂总消耗量M_trans，采用公式(21)计算所需贮箱半径(用符号R′表示)。 

$\frac{M_{\mathrm{GEO}}+M_{\mathrm{trans}}}{N\·\mathrm{\ρ}}=\mathrm{\π}{\mathrm{HR}}^{\′2}-\frac{2}{3}{\mathrm{\πR}}^{\′3}---\left(21\right)$

其中，N为贮箱个数；ρ为推进剂密度。 

步骤17：在步骤16操作的基础上，判断静止轨道阶段推进剂总消耗量M_GEO与所需贮箱半径R′是否满足收敛条件式(22)，如果满足，则流程终止，得到卫星寿命期内消耗推进剂的总质量为静止轨道阶段推进剂总消耗量M_GEO和转移轨道阶段推进剂总消耗量M_trans的总和；如果不满足，则令M_p＝M_GEO和R＝R′，重复执行步骤3至步骤17的操作。 

$\max \left(\right|\frac{M_{\mathrm{GEO}}-M_P}{M_P}|,|\frac{R^{\′}-R}{R}\left|\right)\≤\mathrm{\ϵ}---\left(22\right)$

其中，ε为容差。 

有益效果 

本发明与现有技术相比的优点在于： 

(1)本发明方法充分考虑了不同阶段推进剂消耗量与贮箱半径之间的耦合关系，推进剂的质心位置不再局限于采用寿命中期质心位置估算，而是时时根据当前卫星所在阶段进行更新，明显改进了小推力器所提供的南北位保和东西位保推力大小和控制力矩的准确性。 

(2)本发明在推进剂消耗计算过程中，估算了采用小推力器卸载干扰力矩部分所需推进剂，使得推进剂消耗估算模型与实际工程任务相符，提高了模型准确性与实际工程能力。 

(3)本发明方法针对静止轨道阶段位保过程和转移轨道阶段推进剂在贮箱内位置不同，提出了不同质心计算方法，提高了位保推力和控制力矩计算模型精度，有益于提高计算卫星全寿命周期内推进剂消耗量计算精度。 

附图说明

图1为本发明具体实施方式中8个小推力器布局构型示意图； 

图2为本发明具体实施方式中卫星寿命期内南北位保所需速度增量示意图。 

图3为本发明具体实施方式中寿命期内东西位保每个周期内所需推进剂消耗量示意图； 

图4为本发明具体实施方式中寿命期内南北位保每个周期内所需推进剂消耗量示意图。 

具体实施方式

下面结合附图和具体实施例对本发明做详细说明。 

本实施例中，卫星平台共有8个小推力器，通过各自组合实现4个方位的定点位置保持和三轴姿态控制，8个小推力器(1A、2A、3A、4A、1B、2B、3B、4B)安装于长方体形状的卫星平台的8个顶角，安装位置如图1所示。8个小推力器安装角度如表1所示，小推力器的组合方式如表2所示。 

本实施例中卫星的变轨次数为5次，每次所需的速度增量分别为[461.2，404，584，170，102](m/s)。变轨期间，由于推力偏心产生的X，Y，Z三轴干扰力矩大小为[6，6，0.02](N·m)。 

表1 小推力器安装角度 

表2 小推力器组合方式 

采用本发明计算推进剂消耗具体步骤如下： 

步骤1：定义卫星机械坐标系和卫星本体坐标系，给出卫星的基本设计参数，如表3所示，并给定静止轨道阶段推进剂消耗量预估值M_p＝2000kg和贮箱半径预估值R＝450mm。 

表3 卫星基本设计参数 

卫星机械坐标系的原点位于卫星与运载火箭机械分离面内，原点固联在卫星上，并与星箭对接环圆心重合，用符号X、Y和Z分别表示卫星机械坐标系的三个轴；其中，Z轴垂直于卫星与运载火箭的对接环面，其正方向从原点指向对地板，X轴的正方向从原点指向卫星东板，Y轴正方向从原点指向卫星南板。 

卫星本体坐标系的原点在卫星质心，用符号x_b、y_b和z_b分别表示卫星本体坐标系的三个轴；x_b、y_b和z_b分别平行于卫星机械坐标系X轴、Y轴和Z轴。 

步骤2：在步骤1操作的基础上，计算得到卫星整个寿命期内南北位保次数n_SN为391次，两次南北位保间隔期内东西位保次数n_EW为28次。 

步骤3：在步骤2操作的基础上，通过公式(3)得到液面半张角γ，用公式(4)计算静止轨道期间第k次南北位置保持时贮箱内推进剂质心位置Y_mc。当k＝1时，贮箱内推进剂质心位置Y_mc为[0，0，1755](mm)。 

步骤4：在步骤3操作的基础上，通过公式(5)至公式(11)计算静止轨道期间第k次位置保持时位保推力矢量和X轴、Y轴和Z轴控制力矩。 

此时能提供的向南位保推力为17.4N，控制力矩为：+Mx＝26Nm，-Mx＝-23Nm，+My＝23Nm，-My＝-23Nm，+Mz＝5.6Nm，-Mz＝-6.5Nm。 

步骤5：在步骤2操作的基础上，通过公式(12)计算东西位保所需速度增量ΔV_EW为0.0045m/s。 

步骤6：在步骤5操作的基础上，静止轨道阶段，小推力器提供推力，通过公式(13)计算东西位保所需推进剂质量Δm。 

步骤7：在步骤5操作的基础上，通过公式(14)至公式(16)计算卸载东西位保干扰力矩所需的推进剂质量Δm_c。 

此时得到在南北位保间隔期内由于东西位保和卸载干扰力矩所消耗的推进剂为0.79kg。 

步骤8：在步骤5操作的基础上，通过公式(17)至公式(20)计算静止轨道期间第k次南北位置保持时所需的速度增量ΔV_SN。当k＝1时，ΔV_SN＝1.7m/s. 

步骤9：在步骤8操作的基础上，计算静止轨道期间第k次南北位置保持时所需推进剂质量Δm′和卸载干扰力矩所需推进剂质量Δm′_c。当k＝1时，Δm′＝2.57kg，Δm′_c＝1.24kg。 

步骤10：在步骤9操作的基础上，计算当前贮箱内推进剂质量M′_p，M′_p＝M_p-Δm-Δm_c-Δm′-Δm′_c＝1995.4kg，并使用M′_p更新贮箱内推进剂质量M_p，M_p＝M′_p。判断k＜n_SN是否成立，如果成立，则使k值增1，并重复执行步骤3至步骤9的操作；否则，计算静止轨道阶段推进剂消耗总量M_GEO，并执行步骤11的操作； 

经过391次位保，此时卫星寿命结束，得到静止轨道阶段推进剂消耗总量M_GEO为1277kg。 

步骤11：在步骤10操作的基础上，通过公式(19)计算转移轨道阶段第l次轨道转移时贮箱内推进剂质心高度Z_mc。l的初始值设置为5，此时质心位置为[0，0，1770](mm)。 

步骤12：在步骤11操作的基础上，轨道转移阶段，主发动机提供推力进行轨道转移，在第l次轨道转移时所需的速度增量为已知量的情况下，通过公式(20)计算第l次轨道转移时所需推进剂质量Δm″。当前整星质量为3277kg。根据设计需求，第5次轨道转移所需速度增量为102m/s，轨道转移所需时间为1204s，经过计算得到第5次轨道转移所需推进剂质量为116.4kg。 

步骤13：在步骤11操作的基础上，计算第5次轨道转移时X轴、Y轴和Z轴控制力矩为：+Mx＝26Nm，-Mx＝-23Nm，+My＝23Nm，-My＝23Nm，+Mz＝5.6Nm，-Mz＝-6.5Nm； 

步骤14：在步骤13操作的基础上，计算第5次轨道转移时卸载干扰力矩所需推进剂质量Δm″_c＝2.6kg。 

步骤15：在步骤14操作的基础上，判断l的值是否不为1，如果不为1，则使l值减1，并重复执行步骤11至步骤14的操作；否则，计算转移轨道阶段推进剂消耗总量M_trans，并执行步骤16的操作； 

l＝1时，得到卫星在轨道转移阶段消耗推进剂总量M_trans为2347kg。 

步骤16：在步骤15操作的基础上，根据静止轨道阶段推进剂总消耗量M_GEO和转移轨道阶段推进剂总消耗量M_trans，采用公式(21)计算所需贮箱半径R′。本实例卫星平台采用4个贮箱并排的方式，单个贮箱上下两端分别是半球形，中间段为圆柱体，4个贮箱并排放在对接环上。在静止轨道阶段和转移轨道阶段总消耗推进剂为3624kg，计算所需贮箱半径为436mm。 

步骤7：设定收敛容差ε＝0.01，根据收敛条件公式计算相对误差为0.57，不满足收敛条件；则根据更新后的静止轨道阶段推进剂消耗量和贮箱半径，返回步骤2，重新计算，直至满足收敛条件。最后，计算出来静止轨道阶段推进剂消耗量为924kg，转移轨道阶段推进剂消耗量为2094kg，贮箱半径为393mm。 

步骤17：在步骤16操作的基础上，判断静止轨道阶段推进剂总消耗量M_GEO与所需贮箱半径R′是否满足收敛条件式(22)，如果满足，则流程终止，得到卫星寿命期内消耗推进剂的总质量为静止轨道阶段推进剂总消耗量M_GEO和转移轨道阶段推进剂总消耗量M_trans的总和；如果不满足，则令M_p＝M_GEO和R＝R′，重复执行步骤3至步骤17的操作。 

设定收敛容差ε＝0.01，根据收敛条件公式计算相对误差为0.57，不满足收敛条件；则根据更新后的静止轨道阶段推进剂消耗量和贮箱半径，返回步骤3，重新计算，直至满足收敛条件。最后，计算出来静止轨道阶段推进剂消耗量为924kg，转移轨道阶段推进剂消耗量为2094kg，贮箱半径为393mm。 

图2给出了卫星在南北位保所需的速度增量，主发动机在5次轨道转移消耗推进剂分别为[656，510，601，177，104](kg)，小推力器在轨道转移阶段提供姿态控制力矩消耗推进剂依次为[14.8，11.5，13.6，4，2.4](kg)。在静止轨道阶段，推进剂消耗量为924kg，其中，东西位保与姿态控制共消耗推进剂153kg，南北位保与姿态控制共消耗推进剂771kg。在寿命期内，卫星东西位保每个周期内所需推进剂消耗量如图3所示，南北位保每个周期内所需推进剂消耗量如图4所示。 

根据具体的实施案例，本发明方法基本实现了预期的发明目的，能够有效的计算GEO卫星推进剂消耗量，本发明有助于提高卫星小推力器布局设计质量。 

Claims (1)

1.一种GEO卫星小推力器推进剂消耗耦合分析方法，其特征在于：其具体实现步骤如下：

步骤1：定义卫星机械坐标系和卫星本体坐标系，给出卫星的基本设计参数，并给定静止轨道阶段推进剂消耗量预估值M_p和贮箱半径预估值R；所述卫星的基本设计参数包括：卫星寿命、平台尺寸、南北位保周期、东西位保周期、主发动机推力、主发动机比冲、小推力器推力、小推力器比冲、推力器工作效率、贮箱与平台间隙、有效载荷质量、有效载荷质心Z方向高度、平台结构质量、推进剂密度；

所述卫星机械坐标系的原点位于卫星与运载火箭机械分离面内，原点固联在卫星上，并与星箭对接环圆心重合，用符号X、Y和Z分别表示卫星机械坐标系的三个轴；其中，Z轴垂直于卫星与运载火箭的对接环面，其正方向从原点指向对地板，X轴的正方向从原点指向卫星东板，Y轴正方向从原点指向卫星南板；

所述卫星本体坐标系的原点在卫星质心，用符号x_b、y_b和z_b分别表示卫星本体坐标系的三个轴；x_b、y_b和z_b分别平行于卫星机械坐标系X轴、Y轴和Z轴；

步骤2：在步骤1操作的基础上，通过公式(1)和公式(2)分别计算卫星整个寿命期内南北位保次数与两次南北位保间隔期内东西位保次数；

$n_{\mathrm{SN}}=\frac{T_{\mathrm{life}}}{T_{\mathrm{SN}}}---\left(1\right)$

其中，n_SN为卫星整个寿命期内南北位保次数；T_life为卫星寿命；T_SN为南北位保周期；

$n_{\mathrm{EW}}=\frac{T_{\mathrm{SN}}}{T_{\mathrm{EW}}}---\left(2\right)$

其中，n_EW为两次南北位保间隔期内东西位保次数；T_EW为东西位保周期；

步骤3：在步骤2操作的基础上，通过公式(3)得到液面半张角γ，用公式(4)计算静止轨道期间第k次南北位置保持时贮箱内推进剂质心位置Y_mc；其中，k为正整数，1≤k≤n_SN；k的初始值设置为1；

$\frac{H}{R}(\mathrm{\γ}-\sin \mathrm{\γ}\cos \mathrm{\γ})+\frac{\mathrm{\π}}{3}{(1-\cos \mathrm{\γ})}^2(2+\cos \mathrm{\γ})=f\·(\mathrm{\π}\frac{H}{R}+\mathrm{\π}\frac{4}{3})---\left(3\right)$

其中，H为贮箱高度；R为贮箱半径；γ为推进剂液面半张角；f为充液比，即当前推进剂的体积与贮箱容积之比；

$Y_{\mathrm{mc}}=\frac{\frac{H}{R}(\mathrm{\γ}-\sin \mathrm{\γ}\cos \mathrm{\γ})\·\frac{4}{3}\frac{{\sin }^3\mathrm{\γ}}{2\mathrm{\γ}-\sin 2\mathrm{\γ}}}{\frac{H}{R}(\mathrm{\γ}-\sin \mathrm{\γ}\cos \mathrm{\γ})+\frac{\mathrm{\π}}{3}{(1-\cos \mathrm{\γ})}^2(2+\cos \mathrm{\γ})}+\frac{\frac{\mathrm{\π}}{3}{(1-\cos \mathrm{\γ})}^2(2+\cos \mathrm{\γ})\frac{3}{4}\frac{{(1+\cos \mathrm{\γ})}^2}{2+\cos \mathrm{\γ}}}{\frac{H}{R}(\mathrm{\γ}-\sin \mathrm{\γ}\cos \mathrm{\γ})+\frac{\mathrm{\π}}{3}{(1-\cos \mathrm{\γ})}^2(2+\cos \mathrm{\γ})}---\left(4\right)$

其中，Y_mc为贮箱内推进剂质心位置；

步骤4：在步骤3操作的基础上，通过公式(5)至公式(11)计算静止轨道期间第k次位置保持时位保推力矢量和X轴、Y轴和Z轴控制力矩；

用符号表示单个推力器的位保推力矢量，推力矢量在卫星本体坐标系中方位角和仰角分别用符号和θ表示；用符号表示卫星质心在机械坐标系中的位置矢量，用符号表示推力器在卫星机械坐标系中的位置矢量，用符号表示推力器在卫星本体坐标系中的位置矢量，

用符号M_x、M_y和M_z分别表示推力器产生的卫星机械坐标系X轴、Y轴和Z轴控制力矩；用符号α、β和μ分别表示推力矢量与卫星本体坐标系的x_b、y_b和z_b三个轴的夹角，则有公式(5)至公式(7)表示的关系；

其中，和θ分别为卫星本体坐标系中方位角和仰角；

μ＝arccos(sinθ)                               (7)

通过公式(8)可得到推力器产生的X轴控制力矩M_x；

$\frac{M_x}{F\·D}=-\tilde{z}\cos \mathrm{\β}+\tilde{y}\cos \mathrm{\μ}---\left(8\right)$

其中，F和D分别为推力矢量和位置矢量的模，F和D为人为设定值； $\tilde{z}=(z-z_c)/D;\tilde{y}=(y-y_c)/D;$

通过公式(9)可得到推力器产生的Y轴控制力矩M_y；

$\frac{M_y}{F\·D}=\tilde{z}\cos \mathrm{\α}-\tilde{x}\cos \mathrm{\μ}---\left(11\right)$

其中， $\tilde{x}=(x-x_c)/D;$

通过公式(10)可得到推力器产生的Z轴控制力矩M_z；

$\frac{M_z}{F\·D}=-\tilde{y}\cos \mathrm{\α}+\tilde{x}\cos \mathrm{\β}---\left(10\right)$

通过公式(11)可得到位保推力矢量

$\stackrel{\→}{F}=[F\cos \mathrm{\α},F\cos \mathrm{\β},F\cos \mathrm{\μ}]---\left(11\right)$

步骤5：在步骤2操作的基础上，通过公式(12)计算东西位保所需速度增量ΔV_EW；

$\left\{\begin{array}{c}\frac{\sin (0.5{\mathrm{\ω}}_s{T}_{\mathrm{EW}}+\mathrm{\δ})}{\sin \left(0.5{\mathrm{\ω}}_s{T}_{\mathrm{EW}}\right)}=\frac{e_r}{e_k}\\ \mathrm{\Δ}{V}_{\mathrm{EW}}=V_s{e}_k\sin \mathrm{\δ}\end{array}\right.---\left(12\right)$

其中，ω_s为地球公转角速度；T_EW为东西位保周期；δ为轨道偏移角度；e_r为偏心率摄动圆半径，e_r＝1.5a_s/(ω_sω_xa)；ω_x为地球自转角速度；a_s为太阳辐射压加速度；a为静止轨道半长轴；e_k为偏心率容许值；V_s为静止轨道速度；

步骤6：在步骤5操作的基础上，静止轨道阶段，小推力器提供推力，通过公式(13)计算东西位保所需推进剂质量Δm；

$\mathrm{\Δm}=m_0[1-e^{{-\mathrm{\ΔV}}_{\mathrm{EW}}/\left(I_{\mathrm{sp}}\mathrm{\ηg}\right)}]---\left(13\right)$

其中，m₀为推力器开始工作时卫星的初始质量；I_sp为推力器比冲，η为推力器工作效率，I_sp和η为人为设定值；g为重力加速度，为常量；

步骤7：在步骤5操作的基础上，通过公式(14)至公式(16)计算卸载东西位保干扰力矩所需的推进剂质量Δm_c；

$\left\{\begin{array}{c}{M}_{x}t={\mathrm{\λM}}_{\mathrm{Tx}}{t}_x\\ M_{y}t={\mathrm{\λM}}_{\mathrm{Ty}}{t}_y\\ M_{z}t={\mathrm{\λM}}_{\mathrm{Tz}}{t}_z\end{array}\right.---\left(14\right)$

其中，t为小推力器工作时间，可通过公式(15)得到；λ为卸载干扰力矩的效率，为人为设定值；M_Tx、M_Ty和M_Tz分别表示X轴、Y轴和Z轴干扰力矩的大小；t_x表示卸载X轴干扰力矩的工作时间；t_y表示卸载Y轴干扰力矩的工作时间；t_z表示卸载Z轴干扰力矩的工作时间；

$F\·t/\left(I_{\mathrm{sp}}g\right)=m_0[1-e^{{-\mathrm{\ΔV}}_{\mathrm{EW}}/\left(I_{\mathrm{sp}}\mathrm{\ηg}\right)}]---\left(15\right)$

Δm_c＝FΔt/(I_spg)    (16)

其中，Δt＝t_x+t_y+t_z；

步骤8：在步骤5操作的基础上，通过公式(17)至公式(18)计算静止轨道期间第k次南北位置保持时所需的速度增量ΔV_SN；

其中，Δi为轨道倾角变化量；Ω_mT为终了月球升交点黄经，Ω_m0为初始月球升交点黄经；Ω_m为月球升交点黄经，Ω_m＝259.183°-0.05295°·t′，t′为相对于1900年1月1日12时的儒略日；

ΔV_SN＝ΔiπV_s/180    (18)

步骤9：在步骤8操作的基础上，计算静止轨道期间第k次南北位置保持时所需推进剂质量Δm′和卸载干扰力矩所需推进剂质量Δm′_c，此过程与步骤6计算东西位保时所需推进剂质量和步骤7计算卸载东西位保干扰力矩所需推进剂质量的方法相同；

步骤10：在步骤9操作的基础上，计算当前贮箱内推进剂质量M′_p，M′_p＝M_p-Δm-Δm_c-Δm′-Δm′_c，并使用M′_p更新贮箱内推进剂质量M_p，M_p＝M′_p；判断k＜n_SN是否成立，如果成立，则使k值增1，并重复执行步骤3至步骤9的操作；否则，计算静止轨道阶段推进剂消耗总量M_GEO，并执行步骤11的操作；

所述推进剂消耗总量M_GEC为静止轨道阶段第1次至第n_SN次位保时Δm、Δm_c、Δm′、Δm′_c四个量之和；

步骤11：在步骤10操作的基础上，通过公式(19)计算转移轨道阶段第l次轨道转移时贮箱内推进剂质心高度Z_mc；其中，l为正整数，1≤l≤l_trans；l的初始值设置为转移轨道次数l_trans，l_trans为人为设定值，3≤l_trans≤8；

$Z_{\mathrm{mc}}=\frac{{8\mathrm{RH}}_0-{3H}_0^2}{4(3R-H_0)}---\left(19\right)$

其中，R为贮箱半径；H₀为推进剂液面高度；

步骤12：在步骤11操作的基础上，轨道转移阶段，主发动机提供推力进行轨道转移，在第l次轨道转移时所需的速度增量为已知量的情况下，通过公式(20)计算第l次轨道转移时所需推进剂质量Δm″；

${\mathrm{\Δm}}^{\′\′}=m_0[1-e^{-\mathrm{\ΔV}/\left(I_{\mathrm{sp}}\mathrm{\ηg}\right)}]---\left(20\right)$

其中，ΔV为在第l次轨道转移时所需的速度增量；

步骤13：在步骤11操作的基础上，计算第l次轨道转移时X轴、Y轴和Z轴控制力矩，此过程与步骤4中计算静止轨道期间第k次位置保持时X轴、Y轴和Z轴控制力矩的方法相同；

步骤14：在步骤13操作的基础上，计算第l次轨道转移时卸载干扰力矩所需推进剂质量Δm″_c，此过程与步骤7计算卸载东西位保干扰力矩所需推进剂质量的方法相同；

步骤15：在步骤14操作的基础上，判断l的值是否不为1，如果不为1，则使l值减1，并重复执行步骤11至步骤14的操作；否则，计算转移轨道阶段推进剂消耗总量M_trans，并执行步骤16的操作；

所述转移轨道阶段推进剂消耗总量M_trans为转移轨道阶段第1次至第l_trans次轨道转移时Δm″和Δm″_c之和；

步骤16：在步骤15操作的基础上，根据静止轨道阶段推进剂总消耗量M_GEO和转移轨道阶段推进剂总消耗量M_trans，采用公式(21)计算所需贮箱半径R′；

$\frac{M_{\mathrm{GEO}}+M_{\mathrm{trans}}}{N\·\mathrm{\ρ}}=\mathrm{\π}{\mathrm{HR}}^{\′2}-\frac{2}{3}\mathrm{\π}{R}^{\′3}---\left(21\right)$

其中，N为贮箱个数；ρ为推进剂密度；

步骤17：在步骤16操作的基础上，判断静止轨道阶段推进剂总消耗量M_GEO与所需贮箱半径R′是否满足收敛条件式(22)，如果满足，则流程终止，得到卫星寿命期内消耗推进剂的总质量为静止轨道阶段推进剂总消耗量M_GEO和转移轨道阶段推进剂总消耗量M_trans的总和；如果不满足，则令M_p＝M_GEO和R＝R′，重复执行步骤3至步骤17的操作；

$\max \left(\right|\frac{M_{\mathrm{GEO}}-M_P}{M_P}|,|\frac{R^{\′}-R}{R}\left|\right)\≤\mathrm{\ϵ}---\left(22\right)$

其中，ε为容差。