[發明專利]一種GEO衛星小推力器推進劑消耗耦合分析方法

在審

• 申請號: 201410320109.3
• 申請日: 2014-07-07
• 公開（公告）號: CN104050338A
• 公開（公告）日: 2014-09-17
• 發明（設計）人: 劉莉;彭磊;龍騰;郭曉松;史人赫
• 申請（專利權）人: 北京理工大學
• 主分類號: G06F17/50
• 分類號: G06F17/50
• 代理公司: 暫無信息
• 代理人: 暫無信息

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)分別計算衛星整個壽命期內南北位保次數與兩次南北位保間隔期內東西位保次數；

nSN=TlifeTSN---(1)]]>

其中，n_SN為衛星整個壽命期內南北位保次數；T_life為衛星壽命；T_SN為南北位保周期；

nEW=TSNTEW---(2)]]>

其中，n_EW為兩次南北位保間隔期內東西位保次數；T_EW為東西位保周期；

步驟3：在步驟2操作的基礎上，通過公式(3)得到液面半張角γ，用公式(4)計算靜止軌道期間第k次南北位置保持時貯箱內推進劑質心位置Y_mc；其中，k為正整數，1≤k≤n_SN；k的初始值設置為1；

HR(γ-sinγcosγ)+π3(1-cosγ)2(2+cosγ)=f·(πHR+π43)---(3)]]>

其中，H為貯箱高度；R為貯箱半徑；γ為推進劑液面半張角；f為充液比，即當前推進劑的體積與貯箱容積之比；

Ymc=HR(γ-sinγcosγ)·43sin3γ2γ-sin2γHR(γ-sinγcosγ)+π3(1-cosγ)2(2+cosγ)+π3(1-cosγ)2(2+cosγ)34(1+cosγ)22+cosγHR(γ-sinγcosγ)+π3(1-cosγ)2(2+cosγ)---(4)]]>

其中，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；

MxF·D=-z~cosβ+y~cosμ---(8)]]>

其中，F和D分別為推力矢量和位置矢量的模，F和D為人為設定值；z~=(z-zc)/D;y~=(y-yc)/D;]]>

通過公式(9)可得到推力器產生的Y軸控制力矩M_y；

MyF·D=z~cosα-x~cosμ---(11)]]>

其中，x~=(x-xc)/D;]]>

通過公式(10)可得到推力器產生的Z軸控制力矩M_z；

MzF·D=-y~cosα+x~cosβ---(10)]]>

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

F→=[Fcosα,Fcosβ,Fcosμ]---(11)]]>

步驟5：在步驟2操作的基礎上，通過公式(12)計算東西位保所需速度增量ΔV_EW；

sin(0.5ωsTEW+δ)sin(0.5ωsTEW)=erekΔVEW=Vseksinδ---(12)]]>

其中，ω_s為地球公轉角速度；T_EW為東西位保周期；δ為軌道偏移角度；e_r為偏心率攝動圓半徑，e_r＝1.5a_s/(ω_sω_xa)；ω_x為地球自轉角速度；a_s為太陽輻射壓加速度；a為靜止軌道半長軸；e_k為偏心率容許值；V_s為靜止軌道速度；

步驟6：在步驟5操作的基礎上，靜止軌道階段，小推力器提供推力，通過公式(13)計算東西位保所需推進劑質量Δm；

Δm=m0[1-e-ΔVEW/(Ispηg)]---(13)]]>

其中，m₀為推力器開始工作時衛星的初始質量；I_sp為推力器比沖，η為推力器工作效率，I_sp和η為人為設定值；g為重力加速度，為常量；

步驟7：在步驟5操作的基礎上，通過公式(14)至公式(16)計算卸載東西位保干擾力矩所需的推進劑質量Δm_c；

Mxt=λMTxtxMyt=λMTytyMzt=λMTztz---(14)]]>

其中，t為小推力器工作時間，可通過公式(15)得到；λ為卸載干擾力矩的效率，為人為設定值；M_Tx、M_Ty和M_Tz分別表示X軸、Y軸和Z軸干擾力矩的大小；t_x表示卸載X軸干擾力矩的工作時間；t_y表示卸載Y軸干擾力矩的工作時間；t_z表示卸載Z軸干擾力矩的工作時間；

F·t/(Ispg)=m0[1-e-ΔVEW/(Ispηg)]---(15)]]>

Δ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；

Zmc=8RH0-3H024(3R-H0)---(19)]]>

其中，R為貯箱半徑；H₀為推進劑液面高度；

步驟12：在步驟11操作的基礎上，軌道轉移階段，主發動機提供推力進行軌道轉移，在第l次軌道轉移時所需的速度增量為已知量的情況下，通過公式(20)計算第l次軌道轉移時所需推進劑質量Δm″；

Δm′′=m0[1-e-ΔV/(Ispηg)]---(20)]]>

其中，Δ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′；

MGEO+MtransN·ρ=πHR′2-23πR′3---(21)]]>

其中，N為貯箱個數；ρ為推進劑密度；

步驟17：在步驟16操作的基礎上，判斷靜止軌道階段推進劑總消耗量M_GEO與所需貯箱半徑R′是否滿足收斂條件式(22)，如果滿足，則流程終止，得到衛星壽命期內消耗推進劑的總質量為靜止軌道階段推進劑總消耗量M_GEO和轉移軌道階段推進劑總消耗量M_trans的總和；如果不滿足，則令M_p＝M_GEO和R＝R′，重復執行步驟3至步驟17的操作；

max(|MGEO-MPMP|,|R′-RR|)≤ϵ---(22)]]>

其中，ε為容差。