月球探測器著陸動響應區間不確定性分析

陳昭岳， 劉莉， 陳樹霖， 崔穎

(北京理工大學 宇航學院 飛行器動力學與控制教育部重點實驗室， 北京 100081)

0 引言

月球探測器成功軟著陸是完成月球探測任務的重要條件。當探測器降落在月球表面時，對月球表面的撞擊會形成較大的沖擊，探測器攜帶的電子儀器設備能否經受得住這一沖擊是軟著陸能否成功的關鍵。探測器軟著陸動力學分析對探測器結構和緩沖系統的設計具有重要指導意義。

只考慮確定參數的著陸過程動力學分析，國內外學者已進行了大量的工作[1-6]。探測器在實際著陸過程中的著陸速度、著陸姿態等都會存在不確定性，相對確定值存在區間波動，現有的文獻對含不確定性因素探測器著陸過程的動力學響應分析較少，并且主要使用蒙特卡洛分析方法。Merchant等[7]將蒙特卡洛方法引入著陸動力學分析中，對“阿波羅”飛船登月艙的沖擊載荷進行了不確定性分析。陳金寶等[8]以4腿“懸臂梁式”月球探測器為研究對象，對132種著陸工況進行仿真分析，研究了著陸姿態和著陸速度對探測器著陸性能的影響。

實際分析中，探測器著陸狀態參數只能獲得其不確定性區間邊界信息，而區間方法適用于這類只含有邊界信息的動力學問題。探測器模型龐大，結構復雜，具有強非線性特點，而目前常用的動力學區間分析方法如Taylor級數法[9-11]、Taylor模型法[12-15]計算過程復雜，且對于強非線性問題無法滿足精度要求。基于Chebyshev多項式的區間分析方法可以有效地處理這類非線性問題，相對于傳統Taylor方法，Chebyshev多項式方法在處理非線性問題時，能更有效地壓縮區間算法的包裹效應，具有更高的求解精度和計算效率。Wu等[16-19]提出了利用Chebyshev多項式展開構造替代模型求解區間動力學問題，并將區間Chebyshev多項式分析方法應用于不確定條件下汽車懸架動態響應分析。Xia等[20]將Chebyshev多項式方法應用于不確定時域系統動態響應分析中。Li等[21]利用稀疏回歸和Chebyshev多項式，提出了高效的非線性動力學系統區間響應分析方法。劉堅等[22]針對多層穿孔板超材料聲學透射率分析，提出了一種區間Chebyshev展開——蒙特卡洛模擬分析方法，該方法在區間變量條件下，高效地分析多層穿孔板超材料聲學透射率的傳輸不確定特性。陳劍等[23]應用Chebyshev區間方法求解汽車動力總成懸置系統的固有頻率和解耦率隨元件剛度值波動的范圍。

本文在Abaqus計算機輔助工程(CAE)軟件平臺上建立了一個全柔性探測器軟著陸非線性有限元模型，并計算探測器傾斜著陸工況下關鍵點的動力學響應。之后基于Chebyshev多項式方法，提出了不確定條件下探測器著陸動響應區間分析流程，對比分析了Chebyshev多項式方法與蒙特卡洛仿真分析方法分析結果。

1 基于Chebyshev多項式展開的探測器著陸動響應區間分析方法

1.1 Chebyshev多項式方法

根據Weierstrass逼近理論，任何在閉區間內連續的函數都可均勻地近似為多項式，并且達到所期望精度要求。

對于變量x的函數f(x)，其輸入參數含有n個不確定性區間變量，且不確定區間范圍分別為xi∈[ai,bi]，i=1，2，…，n，利用k階Chebyshev多項式用于近似原始函數f(x)，以提高其近似精度。

原函數f(x)用n維k階截斷Chebyshev多項式近似表示為

(1)

式中：l為n維Chebyshev級數Ci1,i2,…,in(x)的下標i1,i2,…,in中包含0的個數；Ci1,i2,…,in(x1,x2,…,xn)表示k階Chebyshev級數，

Ci1,i2,…,in(x1,x2,…,xn)=cosi1θ1cosi2θ2…cosinθn，

(2)

(3)

對(3)式采用Gauss-Chebyshev積分，可得

(4)

式中：θjk表示插值點，

(5)

m由多項式截斷階數k確定，一般取m=k+1.

如果函數f(x)具有k+1階連續導數，則Chebyshev多項式與原始函數之間的近似誤差為

(6)

將(1)式中的變量x替換為區間變量[x]，就得到了f(x)的Chebyshev多項式區間表示：

(7)

Ci1,i2,…,in([x])=
cosi1[θ]cosi2[θ]…cosin[θ]=[-1,1].

(8)

(7)式可進一步變換為

(9)

(9)式即為最終得到的區間上下界。

1.2 月球探測器仿真分析模型

本文在Abaqus/CAE平臺上建立用于分析計算的全柔性探測器非線性分析模型。探測器模型主要由著陸器中心體、月球車、緩沖系統組成[24]。

中心體的主承力全部為復合材料夾層制成的板殼結構和起支撐作用的梁結構，所以建模為一個純彈性的殼體和梁結構[25]；緩沖系統包括主著陸腿、輔著陸腿和足墊[26-27]，著陸腿外筒與中心體之間、輔著陸腿內筒與足墊之間的連接均設置為球鉸。有效載荷(雷達、相機、桅桿和慣性組件等)的質量較小，將有效載荷的質量均勻分布到鄰近的板殼結構上。

著陸腿由主緩沖腿和輔助緩沖腿組成，并采用筒式結構，其中鋁蜂窩緩沖器嵌入緩沖腿中。根據參考文獻[24]，結合實際情況，分別使用梁單元與殼單元對緩沖腿內外筒進行建模；著陸腿緩沖過程包含較強非線性，計算較為復雜，本文未考慮著陸腿緩沖載荷行程的參數不確定性，其緩沖載荷行程采用Abaqus連接器模型進行簡化模擬，即主輔緩沖腿的載荷行程曲線簡化成由彈性段和塑性段構成的理想化臺階模式。

月壤[28]建模為立方體實體模型，材料本構模型選擇Cap Drucker-Prager模型用于分析。著陸沖擊仿真有限元模型和響應點如圖1所示。

圖1 探測器有限元分析模型Fig.1 Analysis model of lander

1.3 探測器著陸動響應Chebyshev多項式區間分析計算流程

探測器著陸過程的非線性反映在材料非線性和邊界條件的非線性中，材料非線性指的是月球土壤材料以及緩沖支柱滑動副的彈塑性本構，非線性邊界條件指的是足墊和月壤間的接觸。

由于探測器的非線性特性，在軟著陸動力學過程中每個時間增量步內都要重新計算整個系統的切線剛度矩陣，計算時間長且可重復性差。傳統的區間分析方法，如Taylor級數法和Taylor模型法均為插入式方法，每一時刻需要單獨調用求解器，并重新計算剛度、質量矩陣的偏導數，計算效率低。

Chebyshev區間分析是一種非插入式方法，無需重復調用求解器[16]，因此探測器著陸分析中利用Chebyshev方法可以利用有限元分析求解全過程著陸動響應，而不必每一時刻調用ODE求解器求解當前時刻響應，從而提高編程和計算效率。

本文基于探測器非線性著陸動力學分析模型和Chebyshev區間分析方法提出了探測器著陸動響應Chebyshev區間分析計算流程。采用Chebyshev區間多項式求解帶有區間參數的探測器著陸動力學分析主要包括3個步驟：1)建立探測器非線性著陸動力學分析模型；2)輸入不確定性參數和相應區間，進行參數插值，設定總仿真分析時間ttot，計算每個時刻當區間參數取值為插值點時，關鍵點動力學響應；3)計算每個時刻對應的Chebyshev系數，并構造Chebyshev區間多項式，利用區間算法計算解區間。

Chebyshev區間多項式求解帶有區間參數的探測器著陸動力學分析的實現過程如圖2所示。

圖2 Chebyshev模型法探測器動響應分析流程Fig.2 Flow chart of Chebyshev method for dynamic analysis of lander

2 月球探測器著陸動響應區間分析

2.1 探測器著陸分析工況確定

影響探測器著陸動響應的主要因素為探測器的著陸速度和著陸姿態傾角[3]。探測器著陸速度直接決定探測器的總動能，進而影響探測器動力學響應；著陸過程中著陸傾角有助于保持穩定性，不確定條件下的著陸傾角，對著陸時探測器受到不確定的沖擊作用造成影響。

由于目前我國檢評規范中的評價指標均采用百分制形式，為了與其它指標統一值域，故對橫向裂縫狀況指數TCCI進行百分化處理，以利于結合其他性能評價指標對路面性能進行綜合評價。在此，采用“專家評分技術”，結合TCCI、TCS和TCL、TWR的發展規律分析結果，主客觀相結合，對橫向裂縫狀況指數進行百分制處理。

為分析上述參數對著陸動力學響應的影響，本文考慮著陸姿態角和著陸速度的不確定性，參數區間中值分別設定為4 m/s(等價于地面試驗中從0.87 m高處自由落體)和15°傾斜角，不確定性參數區間半徑分別取5%和15%兩種工況，對應小區間參數范圍和較大區間參數范圍。目前基于傳統Taylor模型法的區間不確定性分析方法主要適用于區間半徑在小不確定性區間的情況，而當不確定性區間較大時，傳統方法需要將不確定性區間分割為數個小區間，分別求解響應結果區間，并進行區間集合，計算效率會顯著降低。2種工況中的不確定性變量變化范圍如表1所示。

表1 不確定性變量區間范圍

著陸動響應有限元分析步長選取0.5 ms，分析時間選取0.15 s，區間分析中使用4 階Chebyshev級數。關鍵點選取探測器中心體的頂板中心點和太陽翼角點。分別利用Chebyshev多項式方法和蒙特卡洛仿真方法計算得到響應點的加速度響應區間，以蒙特卡洛方法得到的數據結果作為準確值。

蒙特卡洛分析中子樣本數最少應在數百，本文中，蒙特卡洛分析中各變量的分布為均勻分布，子樣數為400，可以認為子樣數已足夠。

2.2 探測器著陸動力學響應區間分析

在工況1 和工況2兩種參數不確定區間條件下，考慮著陸姿態角和豎直著陸速度不確定性的探測器著陸動響應區間分析與蒙特卡洛分析結果分別如圖3和圖4所示。

圖3 測試點1不確定響應邊界Fig.3 Interval bound of dynamic response of Point 1

圖4 測試點2不確定響應邊界Fig.4 Interval bound of dynamic response of Point 2

月球探測器著陸動響應不確定性分析中，最關注的是探測器加速度峰值的不確定分析精度。在加速度峰值附近，工況1中Chebyshev區間動響應分析結果與蒙特卡洛分析結果最大相對誤差為12%，工況2中最大相對誤差為16.8%. 從中可以得出，采用Chebyshev模型分析得到的探測器著陸動力學響應區間上界、下界與蒙特卡洛分析得到的區間上界、下界基本吻合，通過Chebyshev 區間多項式方法獲得的分析結果可以完全包含系統的真實解，且結果區間沒有出現被過度放大的包裹效應，波動現象不強。

工況2中不確定性參數的區間半徑超過15%，屬于大不確定性區間問題，但采用Chebyshev模型法獲得的區間結果仍能夠完全包含系統的真實解，結果區間的波動現象并不劇烈，結果區間被放大的較少。

現有算法針對大不確定性問題，均會出現一定程度的計算誤差增大，屬于各區間算法固有的缺點，本算法當區間半徑超過15%時，最大相對誤差為16.8%，仍保留較高的計算精確度。

2.3 Chebyshev模型截斷階數影響分析

圖5 測試點1不確定響應結果對比Fig.5 Comparison of interval bounds of dynamic response of Point 1

圖6 測試點2不確定響應結果對比Fig.6 Comparison of interval bounds of dynamic response of Point 2

由圖5分析對比結果：當k=4時，Chebyshev方法得到的區間響應上界和下界可以包絡蒙特卡洛方法求得響應邊界，且結果的最大相對誤差為15%；k=4和k=9時求得的區間響應上界和下界非常接近，最大相對誤差小于5%，證明本文選取的截斷階數為4，滿足計算精度的要求。

Chebyshev多項式的近似誤差可以表示為

隨著多項式階數k的增大，Chebyshev方法的誤差項ek(f)將迅速減小，因此當選取足夠大階數k時，本方法的分析誤差可忽略。本節分別選取4階和9階Chebyshev進行分析，由誤差表達式可以看出，k取4階和9階時，近似誤差均已足夠小，因此本算例中階數對分析誤差影響可以忽略。

探測器單次著陸仿真所需時間約40 min，構建Chebyshev模型所用時間相對于仿真時間可以忽略不計。

根據算例分析結果，對于低階問題，采用構造較低階數的Chebyshev多項式的方法，如本文構造4階Chebyshev多項式，已滿足精度要求。本文算例構造的Chebyshev多項式僅需要采樣調用有限元模型進行25次仿真計算，計算成本很小。與之對比，蒙特卡洛方法通常需要至少進行數百次的分析，在本算例中，蒙特卡洛方法的子樣本數為400；兩種分析方法結果誤差不超過15%，而Chebyshev方法的采樣數比蒙特卡洛方法小一個數量級，因此有更高的計算效率。

值得注意的是，隨著時間步的推進，Chebyshev區間多項式方法的區間參數預測結果誤差有放大現象，并且多維區間參數下計算次數仍然呈指數型增加，因此不建議在長時間、多維區間參數動力學仿真分析中使用。

3 結論

1)本文構建月球探測器著陸動力學模型，并基于Chebyshev多項式理論，提出了探測器著陸動響應Chebyshev區間分析計算流程，可用于不確定條件下的探測器著陸動響應分析。首先構建月球探測器著陸動力學模型，然后根據輸入的不確定性參數和相應區間進行參數插值，計算每個時刻當區間參數取值為插值點時，關心點的動力學學響應，利用截斷Chebyshev多項式方法計算探測器動力學響應區間結果。

2)采用Chebyshev區間計算方法與蒙特卡洛仿真分析方法對月球探測器著陸動力學響應進行對比分析，結果表明，與蒙特卡洛方法相比，Chebyshev區間方法可以準確、高效地得到探測器著陸動響應上界和下界，區間結果能夠完全包含系統的真實解，結果區間的波動現象并不劇烈，結果區間被放大的較少。對Chebyshev區間方法中多項式級數的選擇進行分析，結果表明恰當的多項式階數選取不僅使分析得到的區間結果滿足精度要求，且結果區間放大較小。