L1點(diǎn)月球天梯繞地運(yùn)行過程動(dòng)力學(xué)

王曉慧，毛李恒

（北京航空航天大學(xué)宇航學(xué)院，北京100191）

1 引言

天梯概念最早由 Artsutanov于1960年提出［1］，然而他的設(shè)想并沒有引起足夠的關(guān)注；1975年，美國人Pearson［2］也獨(dú)立提出了相關(guān)的概念，并將其發(fā)表在了期刊上，人們開始關(guān)注這一新型的運(yùn)輸系統(tǒng)；1979年，Pearson提出了月球天梯的概念，根據(jù)Pearson的描述，月球天梯是穿過地月L1點(diǎn)（或L2點(diǎn)）的一種纜繩結(jié)構(gòu)，一端連著月球表面，另一端連著配重以保持平衡，攀爬器沿著纜繩上下移動(dòng)來運(yùn)送貨物［3］。由于材料的限制，在當(dāng)時(shí)看來建造一座月球天梯是不可能的。1991年，日本科學(xué)家飯島澄男研制出了世界上第一根碳納米管，這種材料的強(qiáng)度是鋼的100倍，密度卻只有鋼的 1／6［4?5］。 碳納米管的出現(xiàn)讓月球天梯的建造變成可能。

目前對于天梯的動(dòng)力學(xué)研究大都針對地球天梯。Noboru Takeichi忽略了纜繩的柔性和彈性，利用連續(xù)體模型對地球天梯纜繩釋放過程的動(dòng)力學(xué)問題開展了研究，并給出了天梯空間站的控制策略［6］；Paul Williams和 Pamela Woo 都對攀爬器的運(yùn)行給地球天梯帶來的影響進(jìn)行了研究，Paul還根據(jù)天梯的動(dòng)力學(xué)響應(yīng)，對攀爬器的運(yùn)行進(jìn)行了優(yōu)化［7?8］；Vladimir 等人利用多點(diǎn)模型，考慮地球引力和大氣阻力，對地球天梯纜繩斷裂后的墜毀過程開展了動(dòng)力學(xué)研究［9］。但目前對于月球天梯的動(dòng)力學(xué)研究還比較少見。

月球天梯動(dòng)力學(xué)問題與繩系衛(wèi)星動(dòng)力學(xué)問題在本質(zhì)上類似，國內(nèi)外學(xué)者對繩系衛(wèi)星的動(dòng)力學(xué)問題開展過研究。Krupa等人利用變分原理建立了繩系衛(wèi)星連續(xù)體模型的偏微分方程，進(jìn)行了動(dòng)力學(xué)和控制方面的研究［10］；于紹華以一組偏微分方程建立的數(shù)學(xué)模型和距離速率控制算法來描述有質(zhì)繩系的受控運(yùn)動(dòng)，并利用遞推算法求得了不可拉長繩系的駐形（即繩系的最終穩(wěn)定形狀）［11］；余本嵩利用一種改進(jìn)的有限差分法對粘彈性有質(zhì)繩系的釋放過程進(jìn)行了數(shù)值模擬［12］。

目前專門針對月球天梯繞地運(yùn)行過程中偏離平衡位置的動(dòng)力學(xué)研究還沒有看到。雖然已有不少針對地球天梯和繩系衛(wèi)星的動(dòng)力學(xué)研究，但這些研究很少同時(shí)考慮纜繩質(zhì)量、柔性和彈性，而且這些研究都只考慮地球引力對繩系的影響，月球天梯則需同時(shí)考慮地球引力和月球引力的影響，這就導(dǎo)致了這些研究方法無法直接用于月球天梯運(yùn)行過程的動(dòng)力學(xué)研究。本文采用離散體模型，同時(shí)考慮纜繩的質(zhì)量、柔性和彈性，建立起月球天梯繞地運(yùn)行過程動(dòng)力學(xué)模型，并利用Newmark法對動(dòng)力學(xué)模型進(jìn)行求解，以得到月球天梯的動(dòng)力學(xué)響應(yīng)。

2 月球天梯繞地運(yùn)行過程動(dòng)力學(xué)模型

月球天梯主要由5部分組成，分別是纜繩、配重、攀爬器、空間站和月面平臺(tái)［13］，其基本結(jié)構(gòu)如圖1所示。

2.1 定義坐標(biāo)系

圖1 月球天梯基本結(jié)構(gòu)Fig．1 Basic architecture of lunar space elevator

假設(shè)月球公轉(zhuǎn)軌道為圓軌道。為了方便動(dòng)力學(xué)方程的描述和結(jié)果的展示，本文自行定義了兩個(gè)坐標(biāo)系：地心?白道慣性坐標(biāo)系OE?xyz和月表旋轉(zhuǎn)坐標(biāo)系 OL?xyz（圖 2）。 地心?白道慣性坐標(biāo)系OE?xyz的原點(diǎn)位于地心，x?OE?y 平面與白道面重合，x軸指向初始時(shí)刻的月心位置，y軸指向該時(shí)刻的月球公轉(zhuǎn)速度方向，z軸由右手定則確定。動(dòng)力學(xué)方程中描述的矢量均為該慣性坐標(biāo)系下的矢量。 月表旋轉(zhuǎn)坐標(biāo)系 OL?xyz的 x?OL?y 平面同樣與白道面重合，原點(diǎn)位于地月連線與月表的交點(diǎn)處，x軸位于地心和月心的連線上且始終指向月心，y軸垂直于x軸指向月球公轉(zhuǎn)速度方向，z軸由右手定則確定。月表旋轉(zhuǎn)坐標(biāo)系OL?xyz主要用于觀察月球天梯配重的振蕩軌跡。則位置矢量r在兩個(gè)坐標(biāo)系下的轉(zhuǎn)換關(guān)系如式（1）所示。

圖2 坐標(biāo)系定義Fig．2 Coordinate system definition

式中，rL表示t時(shí)刻物體在OL?xyz坐標(biāo)系下的位置矢量，rE表示同一時(shí)刻該物體在OE?xyz坐標(biāo)系的位置矢量，rOL表示該時(shí)刻OL點(diǎn)在OE?xyz坐標(biāo)系的位置矢量，ω表示月球公轉(zhuǎn)角速度。

2.2 纜繩離散化

本文采用離散模型來模擬月球天梯纜繩。離散方法是將纜繩看成一系列由無質(zhì)量彈簧連接的質(zhì)點(diǎn)，彈簧可伸長但不可彎曲［14］。纜繩的一端連著月表，另一端連著配重（圖3）。該離散模型既能表現(xiàn)纜繩的柔性，也能體現(xiàn)彈性，同時(shí)也將纜繩的質(zhì)量考慮其中。離散的質(zhì)點(diǎn)越多，則離散纜繩模型越接近真實(shí)模型。纜繩振蕩過程中，離散點(diǎn)個(gè)數(shù)保持不變。

圖3 月球天梯纜繩離散模型Fig．3 Discrete model of lunar space elevator

2.3 動(dòng)力學(xué)方程

假設(shè)配重的質(zhì)量為mc，纜繩的截面積為A，彈性模量為E，密度為ρ，每個(gè)離散單元的長度為l，質(zhì)量為 m0，易知有式（2）：

根據(jù)牛頓第二定律，質(zhì)點(diǎn) i（i＝2，3，…，n）的動(dòng)力學(xué)方程為式（3）：

式中，gi表示質(zhì)點(diǎn)i所處位置的引力加速度，包括地球引力和月球引力；Ti表示質(zhì)點(diǎn)i?1和質(zhì)點(diǎn)i＋1對質(zhì)點(diǎn)i的合拉力，可分別表示為式（4）、（5）：

式（4）中，μE為地球引力常數(shù)，μM為月球引力常數(shù)；rM為月心位置向量。 公式（5）中，li－1，i為質(zhì)點(diǎn) i－1 到質(zhì)點(diǎn) i的實(shí)際距離，li，i＋1為質(zhì)點(diǎn) i到質(zhì)點(diǎn)i＋1的實(shí)際距離，表達(dá)為式（6）：

需要注意的是，公式（3）只描述了質(zhì)點(diǎn)2到質(zhì)點(diǎn)n的動(dòng)力學(xué)方程，不包含質(zhì)點(diǎn)1和纜繩末端重物的動(dòng)力學(xué)方程。質(zhì)點(diǎn)1（配重）的動(dòng)力學(xué)方程與公式（3）相比只有拉力T1發(fā)生如式（7）所示變化：

質(zhì)點(diǎn)n＋1固連于月球表面，繞地球做勻速圓周運(yùn)動(dòng)，故質(zhì)點(diǎn)n＋1的動(dòng)力學(xué)方程為式（8）：

這樣，我們就得到了所有n＋1個(gè)質(zhì)點(diǎn)的動(dòng)力學(xué)方程?，F(xiàn)將公式（3）和公式（8）的等式右端統(tǒng)一用fi（i＝ 1，2，…，n，n ＋1） 表示，則可得到整個(gè)離散纜繩系統(tǒng)的動(dòng)力學(xué)方程為式（9）：

其中各項(xiàng)滿足式（10）：

式中，m2＝m3＝… ＝mn＋1＝m0。

3 動(dòng)力學(xué)方程求解方法

得到動(dòng)力學(xué)方程（10）后，需選擇合適的方法對其進(jìn)行求解。月球天梯纜繩的總長通常超過104km，而且其振蕩周期很可能達(dá)到十幾天甚至幾十天，這就意味著該動(dòng)力學(xué)問題是一個(gè)大自由度、長歷程的問題。因此，求解算法和積分步長的選擇對于該問題的求解非常關(guān)鍵。如果積分步長太短，求解時(shí)間會(huì)變得很長，每調(diào)整一次參數(shù)都需付出很大的時(shí)間成本；如果積分步長太長，則會(huì)帶來較大的誤差，甚至導(dǎo)致算法不收斂。

經(jīng)過比較，本研究選用Newmark法對方程進(jìn)行求解。Newmark法是線性加速法的推廣，它所采用的假定如式（11） ～（12）［15］：

公式（11）和公式（12）中，δ和α的關(guān)系為如式（13）：

δ的幾何意義如圖4所示。

Newmark法的優(yōu)點(diǎn)在于其穩(wěn)定性與系統(tǒng)的故有頻率無關(guān)，當(dāng)其中的兩個(gè)參數(shù)δ＝0.5、α＝0.25時(shí)，算法是無條件穩(wěn)定的［15］，這將節(jié)約大量調(diào)整參數(shù)的時(shí)間；另外，Newmark法至少具有2階精度，而該動(dòng)力學(xué)問題的響應(yīng)頻率又比較小，因此取較大的積分步長也能得到較高的計(jì)算精度［16］，容易在計(jì)算精度和求解時(shí)長之間找到平衡。

動(dòng)力學(xué)方程（10）中只含質(zhì)量矩陣M，不含剛度矩陣K和阻尼矩陣C，則Newmark法的求解步驟如下：

1）形成質(zhì)量矩陣 M0， 計(jì)算初始值 R0、R˙0、R¨0；

2）選取 δ、α、時(shí)間步長 Δt，計(jì)算積分常數(shù)：

4）計(jì)算t＋Δt時(shí)刻的有效載荷：

5）計(jì)算t＋Δt時(shí)刻的位移：

6）計(jì)算t＋Δt時(shí)刻的位移、速度和加速度：

動(dòng)力學(xué)方程求解流程如圖5所示。

4 動(dòng)力學(xué)求解及結(jié)果分析

圖5 動(dòng)力學(xué)方程求解流程Fig．5 The process of dynamic equation solving

各參數(shù)的取值為：地月距離R＝384400 km，月球半徑 RM＝1738 km，地球引力常數(shù) μE＝3.986×1014m3／s2，月球引力常數(shù)μM＝4.902 77 ×1012m3／s2，月球公轉(zhuǎn)角速度 ω ＝2.6617×10－6rad／s，纜繩截面積A＝20 mm2，（碳納米管材質(zhì)）纜繩彈性模量E ＝1000 GPa［5］，纜繩初始應(yīng)力 σ ＝ 40 GPa，纜繩密度 ρ＝1300 kg／m3，纜繩總長度 L＝200 000 km，纜繩單個(gè)離散單元長度l＝100 km，纜繩初始時(shí)刻偏離地月連線的角度θ分別取3°和5°，積分步長Δt＝0.7 s，該積分步長是綜合考慮計(jì)算精度和計(jì)算時(shí)間得到的，一方面，通過初步計(jì)算發(fā)現(xiàn)，Δt取0.1 s和0.7 s結(jié)果沒有明顯差別；另一方面，當(dāng)Δt取0.7 s時(shí)，仿真10 d所需的計(jì)算時(shí)間約為1 h（CPU： Intel Core i7?4790； 軟件： Matlab 2016a），該計(jì)算效率可接受。

配重質(zhì)量可根據(jù)以上數(shù)據(jù)計(jì)算得到。配重受力如圖6所示。

圖6 配重受力分析Fig．6 The force condition of counterweight

圖中，F(xiàn)E為地球引力，F(xiàn)L為月球引力，T為纜繩拉力，F(xiàn)為離心力。受力平衡方程為式（14）：

將各力的具體表達(dá)式及數(shù)據(jù)代入公式（14）求得配重質(zhì)量mc＝8.12×104t。

為了進(jìn)行對比分析，同時(shí)對不考慮纜繩彈性時(shí)的動(dòng)力學(xué)問題進(jìn)行了研究。本文的處理方法為：將纜繩的彈性模量設(shè)置為碳納米管的1000倍（即1000 TPa），此時(shí)，若纜繩應(yīng)力為40 GPa，則纜繩的變形不會(huì)超過0.004%，可近似認(rèn)為沒有彈性。

因?yàn)槔|繩所受的力均在白道平面內(nèi)，故只考慮纜繩在白道平面內(nèi)的運(yùn)動(dòng)，如圖7。

圖7 θ角定義Fig．7 Definition of θ

當(dāng)纜繩初始偏離角分別為3°和5°時(shí)，配重的振幅如圖8所示。從圖中可以看出，配重偏離平衡位置后將做周期性振蕩，月球天梯并不會(huì)墜毀。產(chǎn)生周期性振蕩的原因是當(dāng)配重偏離地月連線時(shí)，地球引力和月球引力會(huì)產(chǎn)生一個(gè)垂直于地月連線的分力，配重在該分力的作用下持續(xù)振蕩。配重的振幅基本保持不變，若考慮纜繩彈性，則偏離3°時(shí)的振幅約為1.125×104km，偏離5°時(shí)的振幅約為1.864×104km，二者的振蕩周期基本一致，約為6.7 d。若忽略纜繩彈性，則偏離3°時(shí)的振幅約為1.047×104km，偏離5°時(shí)的振幅約為1.743×104km，振蕩周期約為7.0 d。忽略纜繩彈性時(shí)振幅稍小，原因是此時(shí)纜繩沒有彈性變形。

圖8 纜繩振幅Fig．8 The amplitude of counterweight

配重的振蕩速度如圖9所示，若考慮纜繩彈性變形，配重偏離3°時(shí)的最大振蕩速度約為135 m／s，偏離5°時(shí)的最大振蕩速度約為208 m／s。若不考慮纜繩彈性變形，配重偏離3°時(shí)的最大振蕩速度約為110 m／s，偏離5°時(shí)的最大振蕩速度約為182 m／s。不考慮纜繩彈性變形時(shí)最大振蕩速度稍小，原因是此時(shí)纜繩的振幅小，回復(fù)力做功少。

圖9 纜繩振蕩速度Fig．9 The hunting speed of counterweight

配重在三個(gè)振蕩周期里的運(yùn)動(dòng)軌跡（月表旋轉(zhuǎn)坐標(biāo)系OL?xyz下）如圖10所示。我們發(fā)現(xiàn)，若考慮纜繩彈性（如圖10（a）和（b）），配重并不是做簡單的單擺運(yùn)動(dòng)，而是在做y方向振蕩的同時(shí)，還伴隨著x方向的振蕩。初始偏離角越大，x方向的振蕩范圍也越大，初始偏離角為3°時(shí)x方向的最大振蕩范圍約為4381 km，初始偏離角為5°時(shí)約為5275 km。如果忽略彈性（如圖10（c）），則配重的振蕩軌跡接近一段圓弧，與預(yù)期相符。

配重在振蕩過程中纜繩應(yīng)力的變化情況如圖11所示?？紤]彈性時(shí)，纜繩的最大應(yīng)力約為60 GPa，相對于40 GPa的初始應(yīng)力來講，增大了約50%。這一結(jié)果提醒我們，月球天梯纜繩的工作應(yīng)力不應(yīng)該過大，否則在振蕩過程中很可能導(dǎo)致纜繩斷裂。不考慮彈性時(shí)，最大應(yīng)力均沒有超過45 GPa，主要原因是此時(shí)沒有x軸方向的振蕩。

圖10 配重運(yùn)動(dòng)軌跡Fig．10 Path of counterweight

圖11 纜繩應(yīng)力Fig．11 Tether Stress

以上結(jié)果都是不加任何控制的自由振蕩結(jié)果，我們同時(shí)對施加一定反向推力的振蕩過程開展了研究?，F(xiàn)假設(shè)配重的初始偏離角為5°，考慮纜繩彈性，在配重上施加一個(gè)10 000 N的持續(xù)推力抑制振蕩，該推力由配重上的反推發(fā)動(dòng)機(jī)提供，方向始終垂直于地月連線。最終得到配重的振幅變化過程如圖12所示。從圖中可以看出，大約經(jīng)過3天配重就可回到平衡點(diǎn)附近，由于推力的存在，配重并沒有完全回到平衡點(diǎn)。振蕩速度的變化如圖13所示，振蕩過程中的最大速度約為102 m／s。

圖12 纜繩振幅（受控條件下）Fig．12 The amplitude of counterweight（under con?trol）

纜繩的應(yīng)力變化過程以及配重的運(yùn)動(dòng)軌跡分別如圖14和圖15所示。纜繩的最大應(yīng)力約為46 GPa。

圖14 纜繩應(yīng)力（受控條件下）Fig．14 Tether Stress （under control）

圖15 配重運(yùn)動(dòng)軌跡（受控條件下）Fig．15 Path of counterweight （under control）

5 結(jié)論

1）月球天梯偏離平衡位置后不會(huì)墜毀，而是會(huì)在一定范圍內(nèi)振蕩，如果不加控制，振幅基本保持不變，且振蕩周期與初始偏離角基本無關(guān)。

2）月球天梯在振蕩過程中，纜繩應(yīng)力最多會(huì)增大50%左右，因此，纜繩的工作應(yīng)力不應(yīng)太大，否則可能導(dǎo)致纜繩斷裂。

3）通過對比發(fā)現(xiàn)，若不考慮纜繩彈性，則月球天梯的振幅較小，振蕩形式較簡單，振蕩過程中的最大應(yīng)力較小，很可能會(huì)低估振蕩的危險(xiǎn)性，因此，開展天梯動(dòng)力學(xué)仿真時(shí)，考慮纜繩彈性是必要的。

4）如果在配重上施加一定的控制力，月球天梯能夠在較短時(shí)間內(nèi)回到平衡位置附近。受控條件下纜繩應(yīng)力及振蕩速度都較小。

本文求得的動(dòng)力學(xué)響應(yīng)能為未來月球天梯的運(yùn)行和控制提供一定的參考。

（

）

［1］ Artsutanov Yuri.V Kosmos na Electrovoze［N］.Komsomol?skaya Pravda， 1960?07?31.Artsutanov Yuri.Into Space on a Train［N］.Komsomolskaya Pravda， 1960?07?31. （in Russian）

［2］ Pearson J.The orbital tower：a spacecraft launcher using the Earth’s rotational energy［J］.Acta Astronautica， 1975， 2（9?10）： 785?799.

［3］ Pearson J.Anchored lunar satellites for cislunar transportation and communication ［ J］.Journal of the Astronautical Sci?ences， 1979， 27（1）： 39?62.

［4］ Iijima S.Helical microtubules of graphitic carbon［J］.Na?ture， 1991， 354（6348）： 56?58.

［5］ Ruoff R S，Qian D，Liu W K，et al.Mechanical properties of carbon nanotubes：theoretical predictions and experimental measurements［J］.Comptes Rendus Physique， 2003， 4（9）：993?1008.

［6］ Takeichi N.Geostationary station keeping control of a space elevator during initial cable deployment［J］.Acta Astronauti?ca， 2012， 70（1?2）： 85?94.

［7］ Williams P，Ockels W.Climber motion optimization for the tethered space elevator［J］.Acta Astronautica， 2010， 66（9）： 1458?1467.

［8］ Woo P，Misra A K.Dynamics of a partial space elevator with multiple climbers［J］.Acta Astronautica， 2010， 67（7）：753?763.

［9］ Aslanov V S，Ledkov A S，Misra A K， et al.Dynamics of space elevator after tether rupture［J］.Journal of Guidance，Control， and Dynamics， 2013， 36（4）： 986?993.

［10］ Krupa M，Poth W，Schagerl M，et al.Modeling，dynamics and control of tethered satellite systems［J］.Nonlinear Dy?namics， 2006， 43（1?2）： 73?96.

［11］ 于紹華，劉強(qiáng).有分布質(zhì)量繩系的衛(wèi)星系統(tǒng)的動(dòng)力學(xué)［J］. 宇航學(xué)報(bào)，2001， 22（3）：52?61.Yu Shaohua， Liu Qiang.Dynamics of mass?distributed tether system［J］.Journal of Astronautics， 2001， 22（3）： 52?61.（in Chinese）

［12］ 余本嵩，文浩，金棟平.時(shí)變自由度繩系衛(wèi)星系統(tǒng)動(dòng)力學(xué)［J］. 力學(xué)學(xué)報(bào)，2010， 42（5）：926?932.Yu Bensong， Wen Hao， Jin Dongpin.Dynamics of tethered satellite system with a time?varying number of degrees?of?free?dom［J］.Chinese Journal of Theoretical and Applied Mechan?ics， 2010， 42（5）： 926?932.（in Chinese）

［13］ Pearson J， Levin E， Oldson J， et al.The Lunar space eleva?tor［C］ ／／55th International Astronautical Congress， Vancou?ver， Canada， 2004.

［14］ 余本嵩.復(fù)雜太空環(huán)境下柔性繩系衛(wèi)星動(dòng)力學(xué)與控制［D］.南京：南京航空航天大學(xué)，2011.Yu Bensong.Dynamics and Control of Flexible Tethered Sat?ellite in Complex Space Environment［D］.Nanjing： Nanjing University of Aeronautics and Astronautics， 2011.（ in Chi?nese）

［15］ 李東旭.高等結(jié)構(gòu)動(dòng)力學(xué)［M］.長沙：國防科技大學(xué)出版社， 1997： 161?163.Li Dongxu.Advanced Dynamics of Structures［M］.Chang?sha： National University of Defense Technology Press， 1997：161?163.（in Chinese）

［16］ Bathe K J， Wilson E L.Stability and accuracy analysis of di?rect integration methods［J］.Earthquake Engineering and Structural Dynamics， 1997， 1（3）： 283?291.