馮維明，李 源，苗 楠

(山東大學(xué) 工程力學(xué)系，濟(jì)南 250062)

基于傅立葉平均法下的連續(xù)小推力動(dòng)力學(xué)分析①

馮維明，李 源，苗 楠

(山東大學(xué) 工程力學(xué)系，濟(jì)南 250062)

通過(guò)將小推力展開(kāi)為偏近點(diǎn)角的傅立葉級(jí)數(shù)，并對(duì)高斯攝動(dòng)方程在一個(gè)軌道周期上的平均，將原方程的推力轉(zhuǎn)化為僅由14個(gè)傅立葉系數(shù)表示的控制變量。仿真計(jì)算表明，平均化后的高斯方程使計(jì)算量與牛頓積分相比顯著減少，且對(duì)小推力而言有足夠的精度。對(duì)利用平均化后的高斯方程計(jì)算軌道根數(shù)時(shí)產(chǎn)生誤差的原因進(jìn)行了研究，并進(jìn)一步分析小推力的范圍和小推力近似表達(dá)式對(duì)上述誤差的影響，為今后小推力下非開(kāi)普勒軌道動(dòng)力學(xué)分析提供了理論依據(jù)和參數(shù)。

傅立葉級(jí)數(shù);連續(xù)小推力;平均法;高斯方程;動(dòng)力學(xué)分析;非開(kāi)普勒軌道

0 引言

小推力推進(jìn)系統(tǒng)為許多星際航行和地球軌道任務(wù)提供了一種高效的新選擇，但是這種系統(tǒng)卻對(duì)最優(yōu)化控制提出了新的挑戰(zhàn)［1］對(duì)于某些特殊的小推力作用的情況下，軌道轉(zhuǎn)移的最優(yōu)控制問(wèn)題已經(jīng)出現(xiàn)了解析法和數(shù)值近似解的方法［2-3］，但是對(duì)于一般的小推力問(wèn)題，需要對(duì)每一個(gè)初始條件和推力值進(jìn)行完整的數(shù)值積分，軌道參數(shù)的確定往往對(duì)這些變量很敏感。因此，普通的最優(yōu)控制法則對(duì)于確定螺旋軌道的數(shù)十甚至上百圈軌道是十分困難的［4］。

在一些特殊情形的軌道轉(zhuǎn)移中，數(shù)值解法已經(jīng)得到了挖掘。目前較多的解法應(yīng)用了變分法或直接優(yōu)化法［5-7］，以確定特定約束條件下的小推力的最優(yōu)控制律。此外，文獻(xiàn)［8］運(yùn)用李亞普諾夫反饋控制的方法，用以解決開(kāi)環(huán)軌道最優(yōu)時(shí)間控制和最優(yōu)軌道問(wèn)題。然而實(shí)踐證明，把平均法和其他方法結(jié)合的思路在解決對(duì)初始軌道和推力變量預(yù)測(cè)敏感的問(wèn)題上是很有效的［9-10］，但所有這些方法都局限于特定的推力和軌道參數(shù)范圍。

本文討論了一種有效解決航天器在小推力下軌道動(dòng)力學(xué)問(wèn)題的新方法。推力加速度的各分量以偏近點(diǎn)角展開(kāi)為傅立葉級(jí)數(shù)，而高斯變分方程則在一個(gè)軌道周期上進(jìn)行平均并且通過(guò)正交條件進(jìn)行簡(jiǎn)化，從而定義出了一組久期方程。這組久期方程是一個(gè)含有14個(gè)推力傅立葉系數(shù)的函數(shù)(無(wú)論初始傅立葉級(jí)數(shù)的階數(shù)是多少)。因此，該方法將普遍形式的推力分布情況簡(jiǎn)化為僅有14個(gè)參數(shù)的形式，這使計(jì)算量與牛頓積分相比大大減少。并以此分析了影響傅立葉系數(shù)表示的平均久期方程求解精度的主要因素。

1 數(shù)學(xué)建模

考慮一個(gè)在軌航天器，其質(zhì)量對(duì)于其環(huán)繞的中心體可忽略不計(jì)。將航天器認(rèn)為是一個(gè)質(zhì)點(diǎn)。它受一個(gè)大小和方向可能隨時(shí)間變化的連續(xù)推力加速度的作用。其軌道可用牛頓運(yùn)動(dòng)方程來(lái)表述:

力加速度矢量F可以沿著徑向FR、法向FW和周向FS分解，即

式中r和w分別為徑向和法向的單位矢量。

牛頓方程可分解為拉格朗日行星方程，高斯形式的拉格朗日行星方程可表示為

式中 a為半長(zhǎng)軸;e為偏心率;i為軌道傾角;Ω為升交點(diǎn)赤經(jīng);ω為近地點(diǎn)幅角;f為真近點(diǎn)角;E為偏近點(diǎn)角，而ε1+∫ndt=l為平黃經(jīng)。

平近點(diǎn)角是平黃經(jīng)和近地點(diǎn)幅角的差值:

根據(jù)傅立葉理論，對(duì)于任何一個(gè)在(0，L)上只存在有限個(gè)跳躍間斷點(diǎn)的分段光滑的函數(shù)f(θ)，都可以表示為一系列周期性延伸的傅立葉級(jí)數(shù)，這個(gè)級(jí)數(shù)無(wú)限逼近于函數(shù)本身，當(dāng)跳躍間斷點(diǎn)存在時(shí)，傅立葉級(jí)數(shù)將逼近于左右極限的平均值。因此，這種表示方法可以應(yīng)用于幾乎任何一種普通的小推力航天器的控制律中。對(duì)于給定的任意加速度矢量F，其各個(gè)分量都可以展開(kāi)為任意時(shí)間間隔上的傅立葉級(jí)數(shù)。實(shí)際上，傅立葉級(jí)數(shù)可按時(shí)間展開(kāi)，也可按隨時(shí)間變化的軌道參數(shù)展開(kāi)，如真近點(diǎn)角、偏近點(diǎn)角或平近點(diǎn)角?？紤]一個(gè)軌道周期L=2π，以θ來(lái)表示這一任意的軌道參數(shù):

現(xiàn)在開(kāi)始一階平均化分析，先假設(shè)一個(gè)加速度矢量，它可在一個(gè)軌道周期(L=2π)中表示，這個(gè)加速度矢量的數(shù)值足夠小，以使軌道的大小和形狀在一圈中的改變不是很明顯。因此，可將高斯方程相對(duì)于平近點(diǎn)角在一個(gè)軌道周期上平均來(lái)得到平均軌道根數(shù)。

式中 γ代表任意的軌道根數(shù)。

在對(duì)推力加速度矢量的各分量進(jìn)行傅立葉展開(kāi)時(shí)，所使用的軌道參數(shù)的選擇是很重要的，有的參數(shù)會(huì)使由此得出的久期方程將變得冗長(zhǎng)而復(fù)雜(如按真近點(diǎn)角展開(kāi)成傅立葉級(jí)數(shù)的話(huà))，而有的參數(shù)會(huì)使久期方程會(huì)被簡(jiǎn)化。經(jīng)比較，加速度矢量的每一個(gè)分量都是按偏近點(diǎn)角進(jìn)行傅立葉級(jí)數(shù)展開(kāi)是最佳選擇，并且各軌道參數(shù)的平均化也是按偏近點(diǎn)角進(jìn)行，由dM=(1-ecosE)dE，式(13)可寫(xiě)為

將式(11)代入式(4)～式(9)中，并由式(14)可得到關(guān)于偏近點(diǎn)角的平均化高斯方程。由傅立葉級(jí)數(shù)的正交性，可消掉推力加速度傅立葉展開(kāi)級(jí)數(shù)中第二階以上及個(gè)別第二階傅立葉系數(shù)。因此，無(wú)論初始推力加速度傅立葉級(jí)數(shù)的階數(shù)如何，軌道根數(shù)a、e、i、Ω、ω和ε1的平均變化率只與這14個(gè)傅立葉系數(shù)有關(guān)，即(k=0，1，2)(k=0，1，2)(k=0，1，2)，(k=1，2)(1，2)則高斯平均方程為

2 計(jì)算分析

為在計(jì)算時(shí)數(shù)據(jù)的簡(jiǎn)便起見(jiàn)，將方程做歸一化處理。在此，將重力參數(shù)μ歸一化常數(shù)，μ=3.986×105km3/s2=3.986×1014m3/s2，因此相當(dāng)于各變量的長(zhǎng)度(以m為單位)除以。初始的軌道根數(shù)由表1給出。

表1 初始軌道根數(shù)Table 1 initial orbit elements

設(shè)控制率為周向加速度FS=1×10-8×(1-sinE)、徑向加速度FR=1×10-8×sinE、法向加速度FW=1×10-8×(1+cosE)，上述推力的加速度分量為歸一化后的值，如航天器為1 000 kg，推力的最大值幅值為2.208 N，可算作小推力范圍。用龍格-庫(kù)塔法進(jìn)行積分來(lái)估算軌道根數(shù)的變化情況。該情況下由牛頓方程和平均變分方程求解出的軌道根數(shù)分別用實(shí)線(xiàn)和虛線(xiàn)在圖1中給出。

圖1 航天器受到3個(gè)方向推力加速度作用時(shí)軌道根數(shù)吻合情況Fig.1 Osculating orbital elements of spacecraft subject to acceleration in three directions

由圖1可看到，所有的軌道根數(shù)均給出了非常吻合的結(jié)果。即用傅立葉系數(shù)表示的平均久期方程的求解結(jié)果和直接對(duì)拉格朗日方程進(jìn)行積分得到的結(jié)果具有緊密的一致性。

再考慮一個(gè)簡(jiǎn)單的周向分步加速度的問(wèn)題，F(xiàn)R=FW=0，F(xiàn)S按如圖2中所示變化。每個(gè)周期中包含2次點(diǎn)火和2個(gè)滑行弧，初始軌道同表1所示。切向分步推力加速度作用時(shí)軌道根數(shù)吻合情況見(jiàn)圖3。

圖2 周向分步加速度Fig.2 Step tangential acceleration

圖3 切向分步推力加速度作用時(shí)軌道根數(shù)吻合情況Fig.3 Osculating orbital elements of spacecraft subject to step tangential acceleration

由圖3知，平均久期方程和牛頓方程給出的結(jié)果具有完全相同的趨勢(shì)，數(shù)值上也相差不大。

另外，基于傅立葉平均化后的高斯方程最顯著優(yōu)勢(shì)之一就是計(jì)算量大規(guī)模減少，尤其在進(jìn)行非開(kāi)普勒軌道優(yōu)化計(jì)算時(shí)，計(jì)算精度和計(jì)算耗時(shí)是困擾著眾多學(xué)者的問(wèn)題，考慮在表2所列初始條件下，歸一化后的周向推力為FS=5×10-9(實(shí)際推力加速度為7.36×10-5m/s2)，分別用牛頓方程和平均變分方程求解出的軌道根數(shù)及誤差列入表2。

推力表達(dá)式為如下形式:

表2 連續(xù)小推力作用下牛頓方程和平均方程10圈后的計(jì)算結(jié)果比較Table 2 Comparison of the calculated results determined by the two methods on low-thrust continuous controls after running 10 laps

航天器運(yùn)行時(shí)間為65 973.4 s，約運(yùn)行10圈。由表2計(jì)算結(jié)果看出，由平均化后的高斯方程計(jì)算得到的10圈后的軌道的6個(gè)根與牛頓方程精確積分所得到的結(jié)果相比，除近地點(diǎn)幅角ω相對(duì)誤差略大些，其他相對(duì)誤差非常之小。利用牛頓法進(jìn)行精確積分計(jì)算耗時(shí)(CPU耗時(shí))為13.556 407 s，而利用基于傅立葉平均化后的高斯方程計(jì)算耗時(shí)為0.450 808 s，后者比前者快30倍。

3 影響計(jì)算精度的主要因素

平均久期方程能夠精確有效解決小推力螺旋軌道問(wèn)題，且與通常的牛頓問(wèn)題的積分方法相比能大大減少計(jì)算量。但每一種近似計(jì)算方法都有其局限性，對(duì)于平均久期方程也是如此，這也是很多文獻(xiàn)所未曾涉及到的問(wèn)題，因此很有必要對(duì)利用該方法進(jìn)行軌道計(jì)算可能導(dǎo)致的精度進(jìn)行分析。

3.1 推力幅值大小的影響

由該方法的限制條件可知，推力幅值過(guò)大，會(huì)導(dǎo)致周期內(nèi)軌道變形太大，從而使平均法的結(jié)果產(chǎn)生累計(jì)誤差，于是，關(guān)心在保證精度條件下，推力的極限值是多少。為討論推力增大對(duì)計(jì)算結(jié)果帶來(lái)的影響，研究推力加速度大小在切向分步加速度形式(如圖2示)的控制率中對(duì)計(jì)算結(jié)果的影響。令周向加速度FS的幅值分別增大為 1 ×10-7(歸一化后)的 1.0、1.5、2.0和2.5倍，限于篇幅僅討論對(duì)軌道形狀影響最大的根數(shù)解半長(zhǎng)軸a和偏心率e，圖4給出其吻合情況(幅值為1×10-7見(jiàn)圖3)。

圖4 不同切向分步推力加速度作用時(shí)軌道根數(shù)吻合情況Fig.4 Osculating orbital elements of spacecraft subject to different step tangential acceleration

為更清晰地表現(xiàn)牛頓方程與平均方程的計(jì)算結(jié)果誤差，現(xiàn)將不同幅值的推力作用下航天器運(yùn)行20圈2種方法計(jì)算得到的半長(zhǎng)軸a、偏心率e及相對(duì)誤差列入表3。

由圖4和表3可看到，隨著加速度幅值的增大，a和e的誤差都在增大，a的誤差增幅較小，而e的誤差顯著增大。在FS達(dá)到2.5×10-7時(shí)，偏心率的2種求解結(jié)果已有顯著偏離，此時(shí)對(duì)應(yīng)質(zhì)量為1 000 kg的航天器，推力大小約為18.4 N，約20圈后，偏心率的結(jié)果已經(jīng)完全不可靠了。另外，由圖4也可看到，即使推力較大時(shí)，在前12圈之內(nèi)，平均久期方程和牛頓方程給出的結(jié)果仍然保持著良好的一致性。由表3可見(jiàn)，半長(zhǎng)軸變化很大，因此，若在一定圈數(shù)內(nèi)完成變軌任務(wù)，平均方程仍有足夠的可靠度，即小推力的上限應(yīng)根據(jù)變軌任務(wù)確定。

表3 推力加速度大小對(duì)a、e計(jì)算誤差的影響Table 3 Influence of acceleration on the accuracy of a and e

3.2 近似推力表達(dá)式的影響

由式(15)～式(20)可知，平均化后的高斯方程僅與14個(gè)傅立葉系數(shù)有關(guān)，推力的傅立葉級(jí)數(shù)高階項(xiàng)會(huì)因平均過(guò)程而被消去，平均化后非零的傅立葉系數(shù)是否能精確的描述原推力形式將成為關(guān)鍵。在本文第1個(gè)算例中，由傅立葉系數(shù)表達(dá)式(12)不難得到各推力加速度矢量分量的傅立葉表達(dá)式中的系數(shù)為=1 ×10-8=1 ×10-8=-1 ×10-8=1 ×10-8和=1×10-8，而其他系數(shù)均為零。將上述非零系數(shù)帶入式(11)可得

高斯方程經(jīng)平均化后得到的推力加速度分量傅立葉級(jí)數(shù)近似表達(dá)式與原初始推力表達(dá)式完全一致，因此傅立葉級(jí)數(shù)近似表達(dá)式對(duì)計(jì)算精度沒(méi)有影響。下面再來(lái)考慮3個(gè)推力分量皆為分步加速度的形式，即

在圖5中可表示加速度分量近似表達(dá)式與原表達(dá)式的差異。圖5中用點(diǎn)表示的線(xiàn)是高斯方程平均后3個(gè)加速度分量變化曲線(xiàn)，注意到此時(shí)徑向推力FR已變成常數(shù)，為分段推力的平均值;而周向推力FS和法向推力FW近似為正弦函數(shù)形式，從圖形上看比較接近原推力形式。這也就是在前例中(僅有周向推力FS)計(jì)算精度較高的原因之一。平均方程與解析解得到的計(jì)算結(jié)果如圖6所示。

圖5 分段切向加速度的傅立葉級(jí)數(shù)表示Fig.5 Fourier series for step circumferential acceleration

圖6 三向分步推力加速度作用時(shí)軌道根數(shù)吻合情況Fig.6 Osculating orbital elements of spacecraft subject to step circumferential acceleration in three directions

圖6表示的是航天器運(yùn)行了10圈的根數(shù)時(shí)程圖。長(zhǎng)半軸a和偏心率e平均方程的解與牛頓積分法的解吻合的較好，軌道傾角i，升交點(diǎn)赤經(jīng)Ω，近地點(diǎn)幅角ω在10圈后誤差開(kāi)始增大，較明顯的為Ω的值，平均結(jié)果開(kāi)始偏離高斯方程積分結(jié)果。

4 結(jié)論

(1)將推力加速度分量展開(kāi)為偏近點(diǎn)角表示的傅立葉級(jí)數(shù)，而高斯變分方程則在一個(gè)軌道周期上進(jìn)行平均并且通過(guò)正交條件進(jìn)行簡(jiǎn)化。該方法將整個(gè)連續(xù)控制問(wèn)題的參數(shù)減少為14個(gè)(無(wú)論初始傅立葉級(jí)數(shù)的階數(shù)是多少)，這使計(jì)算量與牛頓積分相比大大減少。

(2)通過(guò)平面變軌和空間變軌的計(jì)算結(jié)果分析，用該方法計(jì)算連續(xù)和非連續(xù)小推力下的軌道轉(zhuǎn)移是正確的和準(zhǔn)確的。

(3)對(duì)造成計(jì)算誤差的主要因素進(jìn)行了定性和定量分析，推力超出“小推力”的范圍和傅立葉級(jí)數(shù)下推力的近似表達(dá)式(高斯方程平均后)是影響精度的主要因素，即便如此，在一定范圍內(nèi)平均法計(jì)算結(jié)果仍與牛頓方程的計(jì)算結(jié)果吻合良好。

［1］Gao Y.Advances in low-thrust trajectory optimization and flight mechanics ，dissertation thesis［R］.University of Missouri-Columbia，2003，64-12(B):6178.

［2］Craig A Kluever.Direct approach for computing near-optimal low-thrust earth-orbit transfers［J］.Journal of Spacecraft and Rockets，1998，29(1):45-61.

［3］Akella M R，Broucke R A.Anatomy of the constant radial thrust problem［J］.Journal of Guidance，Control，and Dy namics，2002，25(3):563-570.

［4］李俊峰，龔勝平.非開(kāi)普勒軌道動(dòng)力學(xué)與控制［J］.宇航學(xué)報(bào)，2009，30(1):47-53.

［5］Kluever C A.Optimal low-thrust interplanetary trajectories by direct method techniques［J］.Journal of the Astronautical Sciences，1997，45(3):162-247.

［6］Petropoulos A E.Some analytic integrals of the averaged variational equations for a thrusting spacecraft［R］.Interplanetary Network Progress Rept.42-150，Jet Propulsion Lab.，California Inst.of Technology，Pasadena，CA，Aug.2002:1-29.

［7］John T Betts.Very low-thrust trajectory optimization using a direct SQP method［J］.Journal of Computational and Applied Mathematics，2000，120:27-40.

［8］Gurfill P.Nonlinear feedback control of low-thrust orbital transfer in a central graviational field［J］.Acta Astronautica，2007，60(8-9):631-648.

［9］Jennifer S Hudson，Daniel J Scheeres.Reduction of lowthrust continuous controls for trajectory dynamics［J］.Journal of Guidance，Control，and Dynamics，2009，32(3):780-787.

［10］尚海濱，崔平遠(yuǎn)，欒恩杰.基于平均法的小推力轉(zhuǎn)移軌道優(yōu)化研究［C］//25屆中國(guó)控制論文集.2006.

Dynamic analysis of continuous low-thrust based on fourier average method

FENG Wei-ming，LI Yuan，MIAO Nan
(Department of Engineering Mechanics，Shandong University，Jinan 250062，China)

Each component of the thrust vector was expanded as Fourier series in eccentric anomaly and Gauss variational equations were averaged over one orbit period，then the thrust vector was translated to a variable controlled by fourteen Fourier＇s parameters.Simulation results show that these secular equations are sufficient to accurately determine a low-thrust spiral trajectory with significantly reduced computation as compared with integration of the full Newtonian problem.In addition，error causes of orbit elements witch were calculated by the averaged Gauss equations were studied，and influence of low-thrust range and approximate expressions on the error was further analyzed，providing theoretical basis and parameters for dynamic analysis of the Non-Keplerian orbits of lowthrust.

Fourier series;continuous low-thrust;average method;Gauss equations;dynamic analysis;non-Keplerian orbits

V412

A

1006-2793(2012)03-0285-05

2011-08-18;;

2011-10-10。

國(guó)家863項(xiàng)目。

馮維明(1957—)，男，教授，主要研究方向?yàn)榉蔷€(xiàn)性動(dòng)力學(xué)和軌道動(dòng)力學(xué)。E-mail:fwm@sdu.edu.cn

(編輯:呂耀輝)