一种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}(\gamma -\sin \gamma \cos \gamma )+\frac{\pi }{3}{(1-\cos \gamma )}^2(2+\cos \gamma )=f·(\pi \frac{H}{R}+\pi \frac{4}{3})---\left(3\right)$其中，H为贮箱高度；R为贮箱半径；γ为推进剂液面半张角；f为充液比，即当前推进剂的体积与贮箱容积之比；$Y_{\mathrm{mc}}=\frac{\frac{H}{R}(\gamma -\sin \gamma \cos \gamma )·\frac{4}{3}\frac{{\sin }^3\gamma }{2\gamma -\sin 2\gamma }}{\frac{H}{R}(\gamma -\sin \gamma \cos \gamma )+\frac{\pi }{3}{(1-\cos \gamma )}^2(2+\cos \gamma )}+\frac{\frac{\pi }{3}{(1-\cos \gamma )}^2(2+\cos \gamma )\frac{3}{4}\frac{{(1+\cos \gamma )}^2}{2+\cos \gamma }}{\frac{H}{R}(\gamma -\sin \gamma \cos \gamma )+\frac{\pi }{3}{(1-\cos \gamma )}^2(2+\cos \gamma )}---\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 \beta +\tilde{y}\cos \mu ---\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 \alpha -\tilde{x}\cos \mu ---\left(11\right)$其中，$\tilde{x}=(x-x_c)/D;$通过公式(10)可得到推力器产生的Z轴控制力矩M_z；$\frac{M_z}{F·D}=-\tilde{y}\cos \alpha +\tilde{x}\cos \beta ---\left(10\right)$通过公式(11)可得到位保推力矢量$\overrightarrow{F}=[F\cos \alpha ,F\cos \beta ,F\cos \mu ]---\left(11\right)$步骤5：在步骤2操作的基础上，通过公式(12)计算东西位保所需速度增量ΔV_EW；$\left\{\begin{array}{c}\frac{\sin (0.5{\omega }_s{T}_{\mathrm{EW}}+\delta )}{\sin \left(0.5{\omega }_s{T}_{\mathrm{EW}}\right)}=\frac{e_r}{e_k}\\ \Delta V_{\mathrm{EW}}=V_s{e}_k\sin \delta \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{\Delta m}=m_0[1-e^{{-\mathrm{\Delta V}}_{\mathrm{EW}}/\left(I_{\mathrm{sp}}\mathrm{\eta 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{\lambda M}}_{\mathrm{Tx}}{t}_x\\ M_{y}t={\mathrm{\lambda M}}_{\mathrm{Ty}}{t}_y\\ M_{z}t={\mathrm{\lambda 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{\Delta V}}_{\mathrm{EW}}/\left(I_{\mathrm{sp}}\mathrm{\eta 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{\Delta m}}^{\prime \prime }=m_0[1-e^{-\mathrm{\Delta V}/\left(I_{\mathrm{sp}}\mathrm{\eta 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·\rho }=\pi {\mathrm{HR}}^{\prime 2}-\frac{2}{3}\pi R^{\prime 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^{\prime }-R}{R}\left|\right)\le ϵ---\left(22\right)$其中，ε为容差。