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

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

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

0 引言

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

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

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

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

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

1.1 Chebyshev多項式方法

根據(jù)Weierstrass逼近理論，任何在閉區(qū)間內(nèi)連續(xù)的函數(shù)都可均勻地近似為多項式，并且達到所期望精度要求。

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

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

(1)

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

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

(2)

(3)

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

(4)

式中：θjk表示插值點，

(5)

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

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

(6)

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

(7)

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

(8)

(7)式可進一步變換為

(9)

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

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

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

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

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

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

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

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

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

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

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

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

Chebyshev區(qū)間多項式求解帶有區(qū)間參數(shù)的探測器著陸動力學分析的實現(xiàn)過程如圖2所示。

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

2 月球探測器著陸動響應區(qū)間分析

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

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

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

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

表1 不確定性變量區(qū)間范圍

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

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

2.2 探測器著陸動力學響應區(qū)間分析

在工況1 和工況2兩種參數(shù)不確定區(qū)間條件下，考慮著陸姿態(tài)角和豎直著陸速度不確定性的探測器著陸動響應區(qū)間分析與蒙特卡洛分析結(jié)果分別如圖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

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

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

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

2.3 Chebyshev模型截斷階數(shù)影響分析

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

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

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

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

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

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

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

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

3 結(jié)論

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

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