吳宏宇，王春潔，丁宗茂，丁建中，董 洋，滿劍鋒

(1. 北京航空航天大學(xué)機(jī)械工程及自動(dòng)化學(xué)院，北京 100083； 2. 北京航空航天大學(xué)虛擬現(xiàn)實(shí)技術(shù)與系統(tǒng)國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 100083;3. 清華大學(xué)機(jī)械工程系, 北京 100084)

0 引 言

隨著科學(xué)技術(shù)的進(jìn)步，人類的研究空間已不再局限于地球，探索月球、火星、乃至距地球更遠(yuǎn)星球上的奧秘已成為當(dāng)今科學(xué)界的研究熱點(diǎn)之一。而探測(cè)任務(wù)中，將機(jī)器人或宇航員安全送往外星球表面是探測(cè)任務(wù)的關(guān)鍵環(huán)節(jié)[1]。縱觀世界各國(guó)的深空探測(cè)任務(wù)，常用腿式著陸器搭載探測(cè)設(shè)備軟著陸外星球表面，如美國(guó)的阿波羅號(hào)載人登月艙、海盜2號(hào)火星探測(cè)器、鳳凰號(hào)火星探測(cè)器及中國(guó)的嫦娥三號(hào)月球探測(cè)器等。腿式著陸器作為一種簡(jiǎn)單、可靠的著陸緩沖裝置，需要持續(xù)開(kāi)展其研究工作，而研究工作的重點(diǎn)內(nèi)容包括著陸器軟著陸過(guò)程的動(dòng)力學(xué)仿真與著陸器的性能優(yōu)化。

關(guān)于腿式著陸器軟著陸過(guò)程仿真計(jì)算與性能優(yōu)化的研究，已有較多學(xué)者開(kāi)展了相關(guān)的工作。針對(duì)軟著陸過(guò)程仿真分析的研究，文獻(xiàn)[2-6]利用多體動(dòng)力學(xué)分析軟件ADAMS建立了腿式著陸器軟著陸過(guò)程的動(dòng)力學(xué)仿真模型，分析了著陸器在典型工況下的軟著陸性能，其中文獻(xiàn)[5-6]通過(guò)對(duì)比仿真結(jié)果與試驗(yàn)數(shù)據(jù)證明了建模方法的合理性；文獻(xiàn)[7]基于基本力學(xué)原理，研究了著陸器沖擊動(dòng)力學(xué)相似性技術(shù)，提出了等效低重力工況的新方法，并利用物理樣機(jī)試驗(yàn)驗(yàn)證了方法的有效性；文獻(xiàn)[8-10]基于Monte Carlo法研究了不確定著陸工況下的月球和火星著陸器的軟著陸性能；文獻(xiàn)[11-12]則考慮多種工況因子，重點(diǎn)針對(duì)腿式著陸器的抗翻倒能力與載荷緩沖性能繪制了穩(wěn)定性邊界；文獻(xiàn)[13-14]建立了腿式著陸器軟著陸過(guò)程的動(dòng)力學(xué)方程，并對(duì)方程進(jìn)行求解，研究了影響著陸器軟著陸性能的關(guān)鍵因素。針對(duì)腿式著陸器優(yōu)化設(shè)計(jì)的研究，文獻(xiàn)[15]基于響應(yīng)面法與單一惡劣著陸工況，優(yōu)化了月球著陸器的緩沖機(jī)構(gòu)；文獻(xiàn)[16]同時(shí)考慮多組惡劣著陸工況，直接利用動(dòng)力學(xué)仿真模型參與優(yōu)化迭代計(jì)算，優(yōu)化了某型著陸器的緩沖機(jī)構(gòu)；文獻(xiàn)[17]除考慮多組惡劣著陸工況外，同時(shí)考慮某型著陸器的兩種著陸模式，結(jié)合響應(yīng)面法對(duì)某型著陸器的緩沖機(jī)構(gòu)進(jìn)行了優(yōu)化設(shè)計(jì)。

上述研究的集中點(diǎn)在于腿式著陸器軟著陸過(guò)程的仿真方法與特殊工況下的著陸器性能優(yōu)化方法，但鮮有文章報(bào)道著陸姿態(tài)不確定條件下的著陸器優(yōu)化方法。由于控制延遲、誤差等原因，著陸器與星球表面接觸時(shí)姿態(tài)具有隨機(jī)性，故而將著陸姿態(tài)的不確定性引入著陸器的優(yōu)化設(shè)計(jì)過(guò)程中非常必要。本文結(jié)合響應(yīng)面法與分級(jí)優(yōu)化設(shè)計(jì)方法，以某型著陸器為例，重點(diǎn)研究了不確定著陸姿態(tài)下的著陸器緩沖機(jī)構(gòu)優(yōu)化方法。

1 著陸器仿真模型建立

本文以圖1所示的某型著陸器為研究對(duì)象，著陸器由主體和四套完全相同的緩沖機(jī)構(gòu)組成。緩沖機(jī)構(gòu)均勻分布在主體四周，由主支柱、輔助支柱、足墊和緩沖桿組成，輔助支柱內(nèi)外筒間裝有吸收沖擊能量的緩沖元件[16-17]。

采用ADAMS軟件建立非0高度關(guān)機(jī)[17]軟著陸模式下的著陸器軟著陸過(guò)程仿真模型，模型從著陸器與星球表面接觸瞬間開(kāi)始進(jìn)行仿真分析，且各構(gòu)件均簡(jiǎn)化為剛體或利用剛體等效建模。

考慮大變形構(gòu)件緩沖桿主要通過(guò)彎曲變形的方式緩沖載荷，其等效建模方式如圖2所示。首先，建立兩段軸線共線的剛體，其質(zhì)量特性參考真實(shí)緩沖桿賦予。剛體1利用萬(wàn)向節(jié)副與主體在O點(diǎn)聯(lián)接，并限制剛體2只能沿軸線移動(dòng)。緩沖桿最初狀態(tài)時(shí)，軸線重合于X軸，當(dāng)其末端受力導(dǎo)致其彎曲時(shí)，則等效于O點(diǎn)產(chǎn)生反作用彎矩M，將M投影在Y與Z軸上，得到兩個(gè)彎矩分量MY,MZ。反之，在O點(diǎn)處對(duì)剛體1施加兩個(gè)彎矩分量MY,MZ可以合成彎矩M。同時(shí)，為了模擬末端位置的變化，在剛體2上施加運(yùn)動(dòng)S。結(jié)合圖2可知，M,S均為O點(diǎn)與末端連線同X軸夾角α的函數(shù)。利用測(cè)量函數(shù)，求得緩沖桿末端相對(duì)于O點(diǎn)的坐標(biāo)(x,y,z)，進(jìn)而可得α的表達(dá)式如式(1)所示，并可得到M,S的表達(dá)式如式(2)和式(3)所示。

(1)

M=f(α)

(2)

S=g(α)

(3)

建立緩沖桿的有限元模型，末端施加彎矩，可仿真計(jì)算得到α-M，α-S函數(shù)關(guān)系f，g如圖3、圖4所示。

僅考慮受彎，緩沖桿的俯視圖為一條線段，如圖5所示。結(jié)合圖5，根據(jù)幾何關(guān)系，可計(jì)算出M在Y軸、Z軸的分量MY，MZ的表達(dá)式如式(4)、式(5)所示。

(4)

(5)

進(jìn)而，利用式(3)～(5)結(jié)合圖2所示的等效模型可模擬緩沖桿的受力彎曲情況。

輔助支柱內(nèi)緩沖元件的吸能過(guò)程通過(guò)在內(nèi)外筒之間施加與壓縮行程D相關(guān)的緩沖力FD來(lái)模擬，F(xiàn)D與D的函數(shù)關(guān)系如圖6所示。

足墊和星球表面之間的垂向力和法向力分別采用非線性阻尼彈簧模型和庫(kù)倫摩擦模型來(lái)模擬[4,16]。

2 著陸工況參數(shù)及軟著陸性能判據(jù)

2.1 著陸工況參數(shù)

設(shè)著陸器與星球表面接觸瞬間的整機(jī)質(zhì)量為800 kg，豎直速度vx為3.7 m/s、水平速度與三軸角速度為0，星球表面的坡度為10°、動(dòng)摩擦系數(shù)為0.4[9,17]，不考慮著陸器足墊撞擊巖石或落入凹坑，此時(shí)著陸工況參數(shù)主要為著陸器的姿態(tài)角，而通過(guò)調(diào)整俯仰角θf(wàn)與偏航角θp即可獲得評(píng)判第2.1節(jié)提出的軟著陸性能所需的著陸工況?？紤]著陸器整機(jī)的對(duì)稱性，設(shè)定θf(wàn)與θp的取值范圍如表1所示。

表1 著陸姿態(tài)角取值范圍Table 1 Range of parameters of landing attitude angle

θf(wàn)與θp的幾何含義如圖7所示，圖中xM，yM，zM表示著陸器的機(jī)械坐標(biāo)系。

2.2 軟著陸性能判據(jù)

本文主要針對(duì)著陸器的軟著陸穩(wěn)定性，提出以下三點(diǎn)軟著陸性能判據(jù)：

1)著陸器抗翻倒能力判據(jù)。定義包含著陸器任意相鄰兩足墊中心點(diǎn)的豎直平面為翻倒平面，著陸器的抗翻倒能力與著陸器在著陸過(guò)程質(zhì)心距翻倒平面的最小距離LD成正相關(guān)，且LD應(yīng)始終大于0，否則認(rèn)為著陸器翻倒[18]。

2)輔助支柱緩沖行程約束判據(jù)?？紤]輔助支柱內(nèi)部緩沖元件的承載能力有限，需要保證輔助支柱在正常工作狀態(tài)下的最大緩沖行程DM不大于80 mm。

3)著陸器底面抗損壞能力判據(jù)。著陸器的底面抗損壞能力與著陸器尾噴管底面中心點(diǎn)在著陸過(guò)程中距星球表面的最小距離HM成正相關(guān)，且HM過(guò)小不利于探測(cè)車的釋放，在此約束HM不小于350 mm。

3 著陸器軟著陸性能分析與緩沖機(jī)構(gòu)優(yōu)化

3.1 不確定著陸姿態(tài)下的軟著陸性能分析

本文基于ADAMS模型，采用Monte Carlo法分析著陸器在不確定著陸姿態(tài)下的軟著陸性能，且認(rèn)為各種著陸姿態(tài)對(duì)軟著陸性能的影響同等重要。首先利用描述性采樣方法在表1所示的姿態(tài)角取值范圍內(nèi)等概率抽取樣本點(diǎn)，之后將樣本點(diǎn)代入仿真模型進(jìn)行仿真計(jì)算，最后對(duì)仿真模型輸出結(jié)果進(jìn)行統(tǒng)計(jì)，以定量分析著陸器的軟著陸性能。本文抽取2000個(gè)樣本點(diǎn)進(jìn)行仿真分析，根據(jù)2000組仿真結(jié)果繪制性能指標(biāo)值的散點(diǎn)圖如圖8所示，并得到表2所示的軟著陸性能分析結(jié)果。

圖8中，方形散點(diǎn)表示HM小于350 mm的樣本點(diǎn)。

表2 Monte Carlo模擬分析結(jié)果Table 2 Analysis results of Monte Carlo simulation

分析結(jié)果顯示，著陸器底面抗損壞能力不理想，有必要對(duì)著陸器進(jìn)行優(yōu)化設(shè)計(jì)以提升其軟著陸性能。

3.2 緩沖機(jī)構(gòu)優(yōu)化設(shè)計(jì)

3.2.1優(yōu)化設(shè)計(jì)變量

為保證著陸器在運(yùn)載火箭中的安裝要求，輔助支柱與主體的連接點(diǎn)位置、主支柱與緩沖桿的連接點(diǎn)位置、輔助支柱與主支柱的連接點(diǎn)相對(duì)于足墊的位置均難以調(diào)整[17]。由緩沖機(jī)構(gòu)構(gòu)型可知，著陸器足墊的中心點(diǎn)至主體中軸線的水平距離dz、主支柱長(zhǎng)度LM與輔助支柱緩沖元件的緩沖力大小可作為優(yōu)化設(shè)計(jì)變量。dz，LM的幾何表示如圖9所示，緩沖元件緩沖力的大小利用其增減倍數(shù)kz來(lái)表達(dá)。

結(jié)合圖9，可直接分析出dz，LM和kz對(duì)著陸器軟著陸性能有以下三方面的影響：

1)在其他兩個(gè)參數(shù)不變的前提下，dz的增加會(huì)直接加大質(zhì)心距翻倒平面的初值，減小主體的高度，故dz與著陸器的抗翻倒能力成正相關(guān)、與底面抗損壞能力成負(fù)相關(guān)。

2)在其他兩個(gè)參數(shù)不變的前提下，LM的增加會(huì)直接提升著陸器主體的高度，故dz與著陸器底面抗損壞能力成正相關(guān)，而主體提高的同時(shí)著陸器質(zhì)心勢(shì)必增高，故LM的增加會(huì)間接導(dǎo)致著陸器抗翻倒能力下降。

在其他兩個(gè)參數(shù)不變的前提下，kz的增加會(huì)直接導(dǎo)致輔助支柱緩沖行程與主體下降距離減小，不難得知kz與著陸器的抗翻倒能力和輔助支柱最大緩沖行程成負(fù)相關(guān)、與底面抗損壞能力成正相關(guān)。

dz，LM和kz的初值與范圍如表3所示。

表3 設(shè)計(jì)變量范圍及初值Table 3 Initial values and range of design variable

3.2.2響應(yīng)面代理模型

考慮結(jié)合Monte Carlo法，使用ADAMS模型參與緩沖機(jī)構(gòu)的優(yōu)化迭代計(jì)算效率極低。對(duì)此，本文基于響應(yīng)面方法建立LD，DM，HM與dz，LM，kz，θf(wàn)，θp的映射關(guān)系，響應(yīng)面模型的引入將大幅度提高優(yōu)化計(jì)算效率。本文采用二階多項(xiàng)式建立響應(yīng)面模型，包含dz，LM，kz，θf(wàn)，θp五個(gè)參數(shù)的二階多項(xiàng)式表達(dá)式：

φ=β0+β1dz+β2LM+β3kz+β4θf(wàn)+β5θp+

β12dzkz+β13dzθf(wàn)+β14dzθp+β15LMkz+β16LMθf(wàn)+

β17LMθp+β18kzθf(wàn)+β19kzθp+β20θf(wàn)θp

(6)

式中：φ表示LD，DM，HM中任意一項(xiàng)的近似函數(shù)、βi表示多項(xiàng)式的各項(xiàng)系數(shù)。根據(jù)文獻(xiàn)[19]，利用最小二乘法確定βi的取值，進(jìn)而可得響應(yīng)面模型的具體表達(dá)式。

本文根據(jù)均方根相對(duì)誤差值RMSE和決定系數(shù)R2來(lái)判斷響應(yīng)面模型的擬合精度，對(duì)應(yīng)表達(dá)式如式(7)、式(8)所示[19]。

(7)

(8)

利用優(yōu)化拉丁超立方實(shí)驗(yàn)設(shè)計(jì)在表1和表3中的著陸姿態(tài)與緩沖機(jī)構(gòu)構(gòu)型參數(shù)范圍內(nèi)抽取30個(gè)樣本點(diǎn)，利用ADAMS模型計(jì)算出30組LD,DM,HM值如表4所示，進(jìn)而利用以上30組數(shù)據(jù)擬合響應(yīng)面模型。

為避免多項(xiàng)式項(xiàng)數(shù)過(guò)多引起的多項(xiàng)式擺動(dòng)等誤差源，本文引入了響應(yīng)面模型關(guān)鍵項(xiàng)篩選策略，即對(duì)響應(yīng)面模型的多項(xiàng)式進(jìn)行逐項(xiàng)篩選，保留10項(xiàng)來(lái)擬合樣本點(diǎn)，選取RMSE最小的項(xiàng)組合方式來(lái)確定響應(yīng)面模型的最終表達(dá)式如式(9)～式(11)所示。

2.732θpLM-121.739dzkz+1148.830dzLM

(9)

0.685θpLM-120.259dzkz+205.839dzLM

(10)

HM= -21936.723+23353.525dz+13572.952LM-

0.0345θf(wàn)θp+2.612θf(wàn)LM+402.763dzkz-

183.107kzLM

(11)

利用RMSE和R2分析響應(yīng)面模型對(duì)30組樣本點(diǎn)的擬合精度，分析結(jié)果如表5所示。

Dynamic Analysis of Semi-type Floating Offshore Wind Turbine with Failure Conditions Under Metocean Conditions in South China Sea SHI Wei,ZHENG Kan,REN Nianxin (12)

檢驗(yàn)結(jié)果顯示響應(yīng)面模型可準(zhǔn)確描述LD，DM，HM與dz，LM，kz，θf(wàn)，θp的映射關(guān)系。

3.2.3優(yōu)化數(shù)學(xué)模型的建立及求解

上文的分析結(jié)果顯示，著陸器在初始構(gòu)型參數(shù)下的底面抗損壞能力不理想，故本節(jié)將提出優(yōu)化方法計(jì)算dz，LM和kz的最佳取值，以提升著陸器的軟著陸性能。為使各種著陸姿態(tài)在優(yōu)化過(guò)程中起同等重要的作用以保證優(yōu)化結(jié)果的穩(wěn)健性，本文以最大化不確定著陸姿態(tài)下LD與HM的平均值LD-ave與HM-ave為目標(biāo)，以不確定著陸姿態(tài)下DM最大值DM-max與HM的最小值HM-min不超過(guò)許用值作為約束條件建立優(yōu)化數(shù)學(xué)模型，為了防止優(yōu)化后DM-max與HM-min與許用邊界重合，優(yōu)化過(guò)程中將DM許用值降為75 mm,HM許用值提升為360 mm。綜上，得到優(yōu)化數(shù)學(xué)模型如式(12)所示。

表4 實(shí)驗(yàn)設(shè)計(jì)結(jié)果Table 4 Result of design of experiment

表5 響應(yīng)面精度檢驗(yàn)結(jié)果Table 5 Accuracy analysis table of response surface models

max{LD-ave,HM-ave} s.t.DM-max≤75HM-min≥360x(L)O≤xO≤x(U)O(12)

本文結(jié)合式(9)與式(11)，利用多元積分法求解LD-ave與HM-ave，鑒于θf(wàn)與θp相互獨(dú)立，多元積分求解如下所示：

(13)

122.948LM-121.739dzkz+1148.830dzLM

(14)

HM-ave=-21936.723+23353.525dz+13572.952LM-

402.763dzkz-183.107kzLM

(15)

本文結(jié)合式(10)與式(11)，在給定dz，LM，kz取值下，利用MIGA算法求解以θf(wàn),θp為設(shè)計(jì)變量的DM-max與HM-min，算法參數(shù)配置如表6所示。

表6 MIGA算法參數(shù)配置Table 6 Parameters configuration of MIGA

對(duì)于整個(gè)多目標(biāo)優(yōu)化數(shù)學(xué)模型，本文利用NSGA-II算法求解，算法參數(shù)配置如表7所示。

表7 NSGA-II算法參數(shù)配置Table 7 Parameters configuration of NSGA-II

綜上，得到不確定著陸姿態(tài)下的著陸器緩沖機(jī)構(gòu)分級(jí)優(yōu)化設(shè)計(jì)流程如圖10所示。

經(jīng)過(guò)優(yōu)化迭代計(jì)算，得到帕累托最優(yōu)解集，并繪制帕累托前沿曲線如圖11所示。

由于著陸器的翻倒更能直接導(dǎo)致探測(cè)任務(wù)的失敗，故在帕累托最優(yōu)解集中選取著陸器抗翻倒能力最強(qiáng)的解作為本文的最優(yōu)解，dz，LM和kz的取值如表8所示。

表8 最優(yōu)解取值Table 8 Value of optimum solution

將最優(yōu)解代入ADAMS模型，再次抽取2000個(gè)樣本點(diǎn)進(jìn)行Monte Carlo模擬，根據(jù)2000組仿真結(jié)果繪制性能指標(biāo)值的散點(diǎn)圖如圖12所示，并得到表9所示的軟著陸性能分析結(jié)果。

將表2、表9的分析結(jié)果做對(duì)比可知：相比于優(yōu)化前，優(yōu)化后著陸器的HM超出許用值的工況數(shù)由73組減為0組，LD與HM的平均值分別提升了3.546%和5.140%，優(yōu)化效果明顯。

4 結(jié) 論

本文在某型著陸器軟著陸動(dòng)力學(xué)仿真模型的基礎(chǔ)上，結(jié)合Monte Carlo法研究了著陸姿態(tài)不確定條件下的著陸器軟著陸性能，并結(jié)合優(yōu)化拉丁超立方實(shí)驗(yàn)設(shè)計(jì)得到樣本點(diǎn)，建立了以著陸器緩沖機(jī)構(gòu)構(gòu)型參數(shù)和著陸姿態(tài)參數(shù)為輸入、軟著陸性能指標(biāo)值為輸出的響應(yīng)面模型；針對(duì)著陸姿態(tài)不確定條件下的著陸器緩沖機(jī)構(gòu)優(yōu)化問(wèn)題，基于響應(yīng)面模型，提出了結(jié)合NSGA-II算法與MIGA算法的分級(jí)優(yōu)化求解方法，并計(jì)算得到了使著陸器軟著陸性能最佳的緩沖機(jī)構(gòu)構(gòu)型參數(shù)；再次結(jié)合Monte Carlo法與動(dòng)力學(xué)仿真模型校驗(yàn)優(yōu)化結(jié)果，優(yōu)化后著陸器的抗翻倒能力和底面抗損壞能力分別提升了3.546%和5.140%，且尾噴管底面中心點(diǎn)距星球表面的最小距離與輔助支柱最大緩沖行程均不超出許用值。

指標(biāo)超出許用值工況數(shù)平均值LD/mm0868.558DM/mm037.605HM/mm0374.169

本文所述的研究方法具有成本低、效率高等特點(diǎn)，也可為其他型號(hào)腿式著陸器的優(yōu)化設(shè)計(jì)提供參考。