陳昭岳， 劉莉， 陳樹(shù)霖， 崔穎

(北京理工大學(xué) 宇航學(xué)院 飛行器動(dòng)力學(xué)與控制教育部重點(diǎn)實(shí)驗(yàn)室， 北京 100081)

0 引言

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

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

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

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

1 基于Chebyshev多項(xiàng)式展開(kāi)的探測(cè)器著陸動(dòng)響應(yīng)區(qū)間分析方法

1.1 Chebyshev多項(xiàng)式方法

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

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

原函數(shù)f(x)用n維k階截?cái)郈hebyshev多項(xiàng)式近似表示為

(1)

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

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

(2)

(3)

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

(4)

式中：θjk表示插值點(diǎn)，

(5)

m由多項(xiàng)式截?cái)嚯A數(shù)k確定，一般取m=k+1.

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

(6)

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

(7)

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

(8)

(7)式可進(jìn)一步變換為

(9)

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

1.2 月球探測(cè)器仿真分析模型

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

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

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

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

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

1.3 探測(cè)器著陸動(dòng)響應(yīng)Chebyshev多項(xiàng)式區(qū)間分析計(jì)算流程

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

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

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

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

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

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

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

2.1 探測(cè)器著陸分析工況確定

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

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

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

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

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

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

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

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

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

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

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

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

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

2.3 Chebyshev模型截?cái)嚯A數(shù)影響分析

圖5 測(cè)試點(diǎn)1不確定響應(yīng)結(jié)果對(duì)比Fig.5 Comparison of interval bounds of dynamic response of Point 1

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

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

Chebyshev多項(xiàng)式的近似誤差可以表示為

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

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

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

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

3 結(jié)論

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

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