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

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

0 引言

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

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

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

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

1 基于Chebyshev多項式展開的探測器著陸動響應(yīng)區(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ù)導(dǎo)數(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]。

中心體的主承力全部為復(fù)合材料夾層制成的板殼結(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)外筒進行建模；著陸腿緩沖過程包含較強非線性，計算較為復(fù)雜，本文未考慮著陸腿緩沖載荷行程的參數(shù)不確定性，其緩沖載荷行程采用Abaqus連接器模型進行簡化模擬，即主輔緩沖腿的載荷行程曲線簡化成由彈性段和塑性段構(gòu)成的理想化臺階模式。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2.2 探測器著陸動力學(xué)響應(yīng)區(qū)間分析

在工況1 和工況2兩種參數(shù)不確定區(qū)間條件下，考慮著陸姿態(tài)角和豎直著陸速度不確定性的探測器著陸動響應(yīng)區(qū)間分析與蒙特卡洛分析結(jié)果分別如圖3和圖4所示。

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

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

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

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

由圖5分析對比結(jié)果：當k=4時，Chebyshev方法得到的區(qū)間響應(yīng)上界和下界可以包絡(luò)蒙特卡洛方法求得響應(yīng)邊界，且結(jié)果的最大相對誤差為15%；k=4和k=9時求得的區(qū)間響應(yīng)上界和下界非常接近，最大相對誤差小于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ù)預(yù)測結(jié)果誤差有放大現(xiàn)象，并且多維區(qū)間參數(shù)下計算次數(shù)仍然呈指數(shù)型增加，因此不建議在長時間、多維區(qū)間參數(shù)動力學(xué)仿真分析中使用。

3 結(jié)論

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

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