蘭宇馨,岳曉奎

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

?

利用時域配點法的衛(wèi)星相對運動周期軌道求解研究

蘭宇馨,岳曉奎

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

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

時域配點法； 衛(wèi)星相對運動； 非線性代數(shù)方程組； 牛頓法； 周期解; 數(shù)值方法；J2攝動； C-W方程

0 引言

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

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

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

f2ncos(nωft)).

(1)

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

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

f2ncos(nωftj)).

(2)

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

(3)

為便于表達，定義矩陣

(4)

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

[f0f1… f2N]T=

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

由式(1)，f(t)對時間的一階導數(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)

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

當存在J2項時，衛(wèi)星Sj的相對運動在LVLH系中可表示為

(11)

圖1 ECI與LVLH坐標系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為隨相對位置xr，yr，zr變化的函數(shù)，且

(18)

此處：

(19)

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

考慮衛(wèi)星LVLH系中的相對運動由三個方向的分量組成，且頻率各異，假設

對模型進行模擬計算，可得在設計極限荷載和錨桿預拉力作用下基礎各部位受力及位移變形情況，錨桿等效應力及豎向位移圖、高強灌漿層第一主應力及第三主應力云圖、基礎混凝土第一主應力及第三主應力云圖如圖5—圖10所示。根據(jù)應力分布云圖，可以得到高強灌漿層下部附近區(qū)域、下錨板上、下部附近區(qū)域等應力，相應區(qū)域混凝土第一主應力變化曲線如圖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)相同，但其各自相對x、y、z軸有不同的相對運動頻率ωxr，ωyr，ωzr。為確定式(20)～(25)中未知的系數(shù)和頻率，須從式(12)～(17)進行推導。

(26)

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

(27)

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

(28)

式中：

Q=[XrYrZrVxVyVz]T；

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

只要迭代向量Q確定，就能用式(26)給出的關系計算相對運動模型近似的周期解X(t)。

2 TDC方法初值選取

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

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

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

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

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

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

(29)

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

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

3 計算結果與分析

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

在不同偏心率和軌道傾角條件下，1個周期內采用和未采用TDC修正的相對運動軌跡仿真結果分別如圖2～5所示。

圖2 e=0.02，i=π/3時采用和未采用 TDC修正的相對運動軌跡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時采用和未采用 TDC修正的相對運動軌跡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時采用和未采用 TDC修正的相對運動軌跡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時采用和未采用 TDC修正的相對運動軌跡Fig.5 Relative orbit with and without correction of TDC searching scheme under e=0.1 and i=π/3

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

4 結束語

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

[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—)，男，碩士生，主要研究方向為飛行器設計、非合作目標捕獲和軌道力學。

V412.41

A

10.19328/j.cnki.1006-1630.2016.05.015