利用時(shí)域配點(diǎn)法的衛(wèi)星相對(duì)運(yùn)動(dòng)周期軌道求解研究

蘭宇馨,岳曉奎

(西北工業(yè)大學(xué) 航天學(xué)院,陜西 西安 710072)

?

利用時(shí)域配點(diǎn)法的衛(wèi)星相對(duì)運(yùn)動(dòng)周期軌道求解研究

蘭宇馨,岳曉奎

(西北工業(yè)大學(xué) 航天學(xué)院,陜西 西安 710072)

為更精確而有效地求解衛(wèi)星相對(duì)運(yùn)動(dòng)周期性軌道，提出了一種新的數(shù)值求解方法，該法是基于時(shí)域配點(diǎn)(TDC)方法研究存在J2項(xiàng)的非線性相對(duì)運(yùn)動(dòng)模型的周期解。在TDC方法中，預(yù)先將期望得到的周期解假定為傅里葉級(jí)數(shù)展開，再將傅里葉級(jí)數(shù)展開形成的近似解代入原非線性方程中，得到剩余誤差函數(shù)并令其等于0，其本質(zhì)是一種平衡剩余誤差的方法?？赏ㄟ^此方法用C-W方程或T-H方程的相對(duì)運(yùn)動(dòng)軌道的初始條件生成邊界軌道。TDC方法較一般的數(shù)值方法更精確。數(shù)值仿真表明：應(yīng)用TDC方法后，能獲得閉合的相對(duì)運(yùn)動(dòng)軌跡。研究對(duì)解決相對(duì)運(yùn)動(dòng)軌道維持中節(jié)省燃料和編隊(duì)飛行動(dòng)力學(xué)有重要的意義。

時(shí)域配點(diǎn)法； 衛(wèi)星相對(duì)運(yùn)動(dòng)； 非線性代數(shù)方程組； 牛頓法； 周期解; 數(shù)值方法；J2攝動(dòng)； C-W方程

0 引言

近年來，衛(wèi)星在地球軌道上的相對(duì)運(yùn)動(dòng)研究引起了越來越多的關(guān)注。對(duì)兩個(gè)在軌衛(wèi)星的相對(duì)運(yùn)動(dòng)問題，已建立了多種模型，并將其用于不同空間任務(wù)。首先得到發(fā)展、最著名的相對(duì)運(yùn)動(dòng)數(shù)學(xué)模型是Clohessy-Wiltshire(C-W)方程，在解決交會(huì)問題時(shí)有很高的實(shí)用價(jià)值[1]。但C-W方程使用了一些假設(shè)條件，引起模型存在誤差，這些誤差是不能被忽略的：C-W方程提供的初值條件，僅在圓形參考軌道、地球?yàn)榍蛐?、線性重力加速度同時(shí)滿足時(shí)才有效。由于周期相對(duì)運(yùn)動(dòng)軌道對(duì)相對(duì)運(yùn)動(dòng)的維持和航天器編隊(duì)飛行來說十分重要，有必要建立更精確的相對(duì)運(yùn)動(dòng)模型以尋找周期性軌道。大量研究致力于不斷擴(kuò)展C-W方程[2-3]。但當(dāng)考慮非線性微分重力、J2攝動(dòng)和大偏心率的影響時(shí)，這些研究所用的初始條件就無法使用。XU，WANG提出了一種衍生的衛(wèi)星相對(duì)運(yùn)動(dòng)的動(dòng)力學(xué)模型，該模型考慮了J2攝動(dòng)、重力勢(shì)的非線性項(xiàng)，以及偏心率的影響[2]。因動(dòng)力學(xué)方程中不存在任何近似，故這是一種精確的帶J2項(xiàng)的非線性相對(duì)運(yùn)動(dòng)模型。本文采用了上述動(dòng)力學(xué)模型，并結(jié)合TDC方法研究周期性相對(duì)運(yùn)動(dòng)軌道。TDC方法由代紅華等提出，其本質(zhì)是一種平衡剩余誤差的方法，該算法已成功用于解決多種非線性動(dòng)力學(xué)問題[3-6]。在TDC方法中，期望得到的周期解被預(yù)先假定為傅里葉級(jí)數(shù)展開，并將傅里葉級(jí)數(shù)展開形成的近似解代入原非線性方程中，得到剩余誤差函數(shù)[7-8]。再令剩余誤差函數(shù)在選定的配點(diǎn)范圍內(nèi)等于0，就獲得了一組非線性代數(shù)方程，這些方程中都含有傅里葉系數(shù)作為變量，并均可用成熟的方法求解，如經(jīng)典的牛頓迭代(Newton-Raphson)方法。在以前的研究中，TDC方法主要用于求解范德波振蕩方程(Van der Pol oscillator)和達(dá)芬方程(Duffing oscillator)的周期解，本文對(duì)TDC方法進(jìn)行拓展，以近似求解周期性相對(duì)運(yùn)動(dòng)軌道的初值。在軌道近似為圓形且無J2攝動(dòng)時(shí)，這些初始條件能由C-W方程直接獲得，但在攝動(dòng)的影響下，就必須對(duì)C-W方程給出的初值進(jìn)行修正，可通過用TDC方法對(duì)初值進(jìn)行調(diào)整使其配適攝動(dòng)的影響。為此，本文對(duì)用TDC方法的衛(wèi)星相對(duì)運(yùn)動(dòng)周期軌道求解進(jìn)行了研究。

1 TDC方法在周期性軌道問題中實(shí)現(xiàn)原理

為實(shí)現(xiàn)TDC方法，考慮一含部分展開項(xiàng)的傅里葉級(jí)數(shù)構(gòu)成的函數(shù)，用于近似方程的周期解，該函數(shù)可表示為

f2ncos(nωft)).

(1)

式中：N為近似用到的諧波數(shù)；ωf為周期運(yùn)動(dòng)假定的頻率；fi(i=1，2，…，2N)為諧波系數(shù)；f0為f(t)的初值[9]。顯然，TDC方法求得的結(jié)果本身就是周期性變化的。

在f(t)的1個(gè)周期T內(nèi)，采集K個(gè)等分區(qū)間點(diǎn)處的函數(shù)值f(tj)(j=1，2，…，K)，則由式(1)可得

f2ncos(nωftj)).

(2)

可發(fā)現(xiàn)，2N+1個(gè)系數(shù)fi配置點(diǎn)處的f(tj)存在轉(zhuǎn)換關(guān)系

(3)

為便于表達(dá)，定義矩陣

(4)

則，如能得到K個(gè)等分區(qū)間點(diǎn)處的函數(shù)值f(tj)，就能確定諧波系數(shù)fi，即

[f0f1… f2N]T=

E-1[f(t1) f(t2) … f(tK)]T.

由式(1)，f(t)對(duì)時(shí)間的一階導(dǎo)數(shù)可寫為

cos(nωft)-f2nsin(nωft)).

(5)

f2nsin(nωftj)).

(6)

[f1f2… f2N]T.

(7)

式(7)可更簡潔地表示為

(8)

式中：

(9)

由式(3)可知：fj能通過f(tj)表示。將式(3)代入式(8)，可得

(10)

用上述轉(zhuǎn)換方程推導(dǎo)非線性相對(duì)運(yùn)動(dòng)模型的TDC代數(shù)方程。圖1中包含了2個(gè)笛卡爾坐標(biāo)系。其中：地心慣性坐標(biāo)系(ECI)由單位向量[xyz]表示；衛(wèi)星軌道坐標(biāo)系(LVLH)依附于參考衛(wèi)星S0。

當(dāng)存在J2項(xiàng)時(shí)，衛(wèi)星Sj的相對(duì)運(yùn)動(dòng)在LVLH系中可表示為

(11)

圖1 ECI與LVLH坐標(biāo)系Fig.1 ECI and LVLH coordinate frames

(12)

(13)

(14)

zrωzωx-(ζr-ζ)sinisinθ-

r((ηr)2-(η)2)+Fx；

(15)

(ωz)2-(ωx)2)+zrαx-

(ζr-ζ)sinisinθ+Fy；

(16)

zr((ηr)2-(η)2)-(ζr-ζ)cosi+Fz.

(17)

式中：Fx，F(xiàn)y，F(xiàn)z為衛(wèi)星Sj在衛(wèi)星LVLH系中的控制力；ζr，ηr為隨相對(duì)位置xr，yr，zr變化的函數(shù)，且

(18)

此處：

(19)

式(11)與式(12)～(17)是等價(jià)的。通過增加方程數(shù)，避免了原方程中出現(xiàn)的二階導(dǎo)數(shù)項(xiàng)。

考慮衛(wèi)星LVLH系中的相對(duì)運(yùn)動(dòng)由三個(gè)方向的分量組成，且頻率各異，假設(shè)

對(duì)模型進(jìn)行模擬計(jì)算，可得在設(shè)計(jì)極限荷載和錨桿預(yù)拉力作用下基礎(chǔ)各部位受力及位移變形情況，錨桿等效應(yīng)力及豎向位移圖、高強(qiáng)灌漿層第一主應(yīng)力及第三主應(yīng)力云圖、基礎(chǔ)混凝土第一主應(yīng)力及第三主應(yīng)力云圖如圖5—圖10所示。根據(jù)應(yīng)力分布云圖，可以得到高強(qiáng)灌漿層下部附近區(qū)域、下錨板上、下部附近區(qū)域等應(yīng)力，相應(yīng)區(qū)域混凝土第一主應(yīng)力變化曲線如圖11—圖13所示。

xr2ncos(nωxrt)；

(20)

yr2ncos(nωyrt)；

(21)

zr2ncos(nωzrt)；

(22)

vx2ncos(nωxrt)；

(23)

vy2ncos(nωyrt)；

(24)

vz2ncos(nωzrt).

(25)

其中xr，yr，zr，vx，vy，vz的形式與f(t)相同，但其各自相對(duì)x、y、z軸有不同的相對(duì)運(yùn)動(dòng)頻率ωxr，ωyr，ωzr。為確定式(20)～(25)中未知的系數(shù)和頻率，須從式(12)～(17)進(jìn)行推導(dǎo)。

(26)

根據(jù)式(10)，式(16)的左邊可轉(zhuǎn)換為

(27)

將前文的方程代入式(26)，得式(12)～(17)的時(shí)域配點(diǎn)法代數(shù)方程組

(28)

式中：

Q=[XrYrZrVxVyVz]T；

Ji(Q)=(R(Q′i)-R(Q))/δqi.

只要迭代向量Q確定，就能用式(26)給出的關(guān)系計(jì)算相對(duì)運(yùn)動(dòng)模型近似的周期解X(t)。

2 TDC方法初值選取

大量研究表明：主星軌道的偏心率會(huì)對(duì)相對(duì)運(yùn)動(dòng)軌道的形狀產(chǎn)生顯著影響，不僅體現(xiàn)在面內(nèi)運(yùn)動(dòng)，而且包括異面運(yùn)動(dòng)。此外，主星軌道的高度也會(huì)影響相對(duì)運(yùn)動(dòng)軌跡，主要取決于空氣粘滯效應(yīng)的差異。本文主要關(guān)注主星軌道的偏心率和軌道傾角對(duì)相對(duì)運(yùn)動(dòng)的影響。

為能使牛頓迭代順利進(jìn)行，需確定迭代的初值。牛頓迭代法存在的固有局限性——迭代的初值須保證能足夠接近精確解，因此本文采用成熟的相對(duì)運(yùn)動(dòng)模型，為TDC方法提供合適的初始迭代向量。

TDC方法的特點(diǎn)是能直接估計(jì)非線性系統(tǒng)的周期解，這是該法的優(yōu)勢(shì)，但同時(shí)也會(huì)帶來周期運(yùn)動(dòng)頻率ω和周期運(yùn)動(dòng)幅值A(chǔ)的初值的選取問題(此處：A即為諧波系數(shù))。如ω，A初值估計(jì)與真實(shí)解相差太大，牛頓迭代法很有可能不會(huì)向真實(shí)解逼近，即配點(diǎn)應(yīng)在一個(gè)接近真實(shí)周期解的函數(shù)上進(jìn)行，并盡量在1個(gè)周期內(nèi)配點(diǎn)。

在二體相對(duì)運(yùn)動(dòng)問題中，C-W，T-H方程能提供無攝動(dòng)引力場(chǎng)中的相對(duì)運(yùn)動(dòng)周期解的初始化條件S0=[r(t0)v(t0)]，當(dāng)有攝動(dòng)存在時(shí)，該初始條件能提供一條近似周期的軌跡(用RK4方法就能獲得)。然后用這條軌跡估計(jì)出ω，并在該軌跡上配點(diǎn)。

在采用了精確的J2相對(duì)運(yùn)動(dòng)模型后，再利用上述初始條件，就可得逐漸發(fā)散的無邊界相對(duì)運(yùn)動(dòng)軌道。假設(shè)軌道的一個(gè)周期為T，通過采集此周期內(nèi)無邊界軌道上平均分布的點(diǎn)，就能得到TDC方法代數(shù)方程組的初始值Q0。

通常，時(shí)域配點(diǎn)的個(gè)數(shù)應(yīng)與近似函數(shù)中的諧波系數(shù)的個(gè)數(shù)相同。由于ωxr，ωyr，ωzr未知，為求解式(28)，需要3個(gè)額外的約束。本文在起始點(diǎn)引入約束狀態(tài)：xr(t0)=p1，yr(t0)=p2，zr(t0)=p3(p1，p2，p3是預(yù)先定義的)，則三個(gè)額外的約束狀態(tài)可表示為

(29)

式(28)含配置點(diǎn)的位置和速度，式(29)可展開為

只要確定了配置點(diǎn)并假定頻率，就能通過式(1)、(3)中的轉(zhuǎn)換關(guān)系得到一條閉合的軌道。顯然，該軌道開始并不能滿足式(28)，但通過牛頓法或其它迭代方法，配點(diǎn)處的位置、速度和頻率就能得以修正，最終使剩余誤差達(dá)到最小。當(dāng)由配點(diǎn)選擇造成的殘余誤差足夠小時(shí)，相應(yīng)的軌道即可用于近似相對(duì)運(yùn)動(dòng)的周期解。

3 計(jì)算結(jié)果與分析

首先用普通的數(shù)值積分方法(RK4)求解C-W方程，得到一個(gè)發(fā)散的軌道(無非常明顯的周期解)，該軌道在配點(diǎn)后可用TDC進(jìn)行修正，其過程就是牛頓法的迭代過程。

在不同偏心率和軌道傾角條件下，1個(gè)周期內(nèi)采用和未采用TDC修正的相對(duì)運(yùn)動(dòng)軌跡仿真結(jié)果分別如圖2～5所示。

圖2 e=0.02，i=π/3時(shí)采用和未采用 TDC修正的相對(duì)運(yùn)動(dòng)軌跡Fig.2 Relative orbit with and without correction of TDC searching scheme under e=0.02 and i=π/3

圖3 e=0.05，i=π/3時(shí)采用和未采用 TDC修正的相對(duì)運(yùn)動(dòng)軌跡Fig.3 Relative orbit with and without correction of TDC searching scheme under e=0.05 and i=π/3

圖4 e=0.05，i=0時(shí)采用和未采用 TDC修正的相對(duì)運(yùn)動(dòng)軌跡Fig.4 Relative orbit with and without correction of TDC searching scheme under e=0.05 and i=0

圖5 e=0.1，i=π/3時(shí)采用和未采用 TDC修正的相對(duì)運(yùn)動(dòng)軌跡Fig.5 Relative orbit with and without correction of TDC searching scheme under e=0.1 and i=π/3

由圖2～4可知：使用TDS后，相對(duì)運(yùn)動(dòng)軌跡就會(huì)由原來發(fā)散的狀態(tài)變?yōu)榉忾]的曲線，當(dāng)e=0.02，i=π/3時(shí)，T=7 123 s；當(dāng)e=0.05，i=π/3時(shí)，T=7 130 s；當(dāng)e=0.05，i=π/3時(shí)，T=7 130 s。由圖5可知：當(dāng)e=0.1，i=π/3時(shí)，此時(shí)軌道偏心率的值已非常大，但用TDC修正后仍獲得了較理想的周期性相對(duì)運(yùn)動(dòng)軌道，此時(shí)T=7 157 s。

4 結(jié)束語

TDC是一種求解非線性動(dòng)力學(xué)方程周期解的新方法，已成功用于求解平面二自由度翼型的振顫問題[6]。本文討論了在J2項(xiàng)攝動(dòng)存在時(shí)，用TDC求解在軌衛(wèi)星的周期性相對(duì)運(yùn)動(dòng)。通過將這些周期解假設(shè)成為傅里葉級(jí)數(shù)的形式，推導(dǎo)出TDC的動(dòng)力學(xué)方程組，再用牛頓-拉夫森迭代方法對(duì)方程進(jìn)行求解。TDC的優(yōu)點(diǎn)是能處理許多復(fù)雜方程組中的非線性項(xiàng)，使用較簡便，結(jié)果重復(fù)性較好，發(fā)散速度很慢。理論上TDC方法還能用于解決其他的非線性問題，但由于該方法本身的特性(即在開始就假設(shè)方程組的解均為周期變化的解)，目前還局限于求解已知或有可能為周期解的方程組。若方程組本身并不存在周期解，該方法就無法適用。后續(xù)可對(duì)擴(kuò)展TDC方法的適用范圍進(jìn)行研究。隨著許多問題的動(dòng)力學(xué)模型越來越精確，TDC求得周期解的精度也會(huì)越來越高，其在相對(duì)運(yùn)動(dòng)軌道維持中的燃料節(jié)省和保存，以及編隊(duì)飛行動(dòng)力學(xué)的研究中有廣闊的應(yīng)用前景。

[1] CLOHESSY W H. Terminal guidance system for satellite rendezvous[J]. Aerospace Sci, 1960, 27(9): 653-658.

[2] CARTER T E, BRIENT J. Optimal bounded-thrust space trajectories based on linear equations[J]. Journal of Optimization Theory & Applications, 1991, 70(2): 299-317.

[3] MELTON R G. Time-explicit representation of relative motion between elliptical orbits[J]. Journal of Guidance Control & Dynamics, 2000, 23(4): 604-610.

[4] XU G, WANG D. Nonlinear dynamic equations of satellite relative motion around an oblate earth[J]. Journal of Guidance Control and Dynamics, 2008, 31(5): 1521-1524.

[5] DAI H H, SCHNOOR M, ATLURI S N. A simple collocation scheme for obtaining the periodic solutions of the Duffing equation, and its equivalence to the high dimensional harmonic balance method: subharmonic oscillations[J]. Computer Modeling in Engineering & Sciences, 2012, 84(5): 459-497.

[6] DAI H, YUE X, YUAN J, et al. A time domain collocation method for studying the aeroelasticity of a two dimensional airfoil with a structural nonlinearity[J]. Journal of Computational Physics, 2014, 270(3): 214-237.

[7] YUE X K, DAI H H, LIU C S. Optimal scale polynomial interpolation technique for obtaining periodic solutions to the Duffing oscillator[J]. Nonlinear Dynamics, 2014, 77(4): 1455-1468.

[8] DAI H H, YUE X K, YUAN J P. A time domain collocation method for obtaining the third superharmonic solutions to the Duffing oscillator[J]. Nonlinear Dynamics, 2013, 73(73): 593-609.

[9] WANG W, YUAN J, ZHAO Y, et al. Fourier series approximations toJ2-bounded equatorial orbits[J]. Mathematical Problems in Engineering, 2014: 1-8.

[10] KASDIN N J, GURFIL P. Hamiltonian modelling of relative motion[J]. Annals of the New York Academy of Sciences, 2004, 1017(1): 138-157.

Study on a Time Domain Collocation Method to Search for Periodic Orbits of Satellite Relative Motion

LAN Yu-xin, YUE Xiao-kui

(School of Astronautics, Northwestern Polytechnical University, Xi’an 710072, Shaanxi, China)

To study for the periodic orbits of satellite relative motion, a new numerical approach was proposed, based on using the time domain collocation (TDC) method to search for the periodic solutions of an exactJ2nonlinear relative model. In TDC method, the desired periodic solution was preassumed as a truncated Fourier series first, then the approximate solution of the Fourier series was substituted into the equations of nonlinear system, resulting in a residual error function, which was essentially one of the weighted residual methods. The initial conditions for periodic relative orbits of the Clohessy-Wiltshire (C-W) equations or Tschauner-Hempel (T-H) equations could be refined with this approach to generate nearly bounded orbits. Numerical simulations show that the presented TDC searching scheme is more precise. It is obvious that the TDC searching scheme gives the closed orbit of relative motion. Consequently, it is believed that this approach has potential applications on fuel saving relative orbit maintenance and spacecraft formation flying.

Time domain collocation method; Satellite relative motion; Nonlinear algebraic equations; Newton method; Periodic solutions; Numerical simulationslJ2perturbation; C-W equations

1006-1630(2016)05-0095-06

2016-07-20；

2016-09-12

國家自然基金資助(11502203)

蘭宇馨(1991—)，男，碩士生，主要研究方向?yàn)轱w行器設(shè)計(jì)、非合作目標(biāo)捕獲和軌道力學(xué)。

V412.41

A

10.19328/j.cnki.1006-1630.2016.05.015