吳宏宇， 王春潔,2， 丁建中， 王巨濤， 滿劍鋒， 羅 敏(. 北京航空航天大學(xué) 機(jī)械工程及自動(dòng)化學(xué)院, 北京 009； 2. 虛擬現(xiàn)實(shí)技術(shù)與系統(tǒng)國家重點(diǎn)實(shí)驗(yàn)室, 北京 009； . 北京空間飛行器總體設(shè)計(jì)部, 北京 00094)

在嫦娥三號任務(wù)中，我國利用腿式著陸器將探測機(jī)器人送達(dá)月球表面，順利完成月球探測任務(wù)。隨著空間探測技術(shù)的進(jìn)步，我國將開展其它星球的探測計(jì)劃，而腿式著陸器作為一種可靠性較高的著陸緩沖裝置，有必要持續(xù)對其進(jìn)行研究[1]。出于減少物理樣機(jī)試驗(yàn)成本的目的，對著陸器的著陸沖擊過程進(jìn)行仿真分析顯得非常必要[2-3]。為了初步檢驗(yàn)新型著陸緩沖機(jī)構(gòu)的可靠性，并降低仿真建模與計(jì)算的復(fù)雜度，本文建立著陸器的單腿動(dòng)力學(xué)仿真模型以模擬著陸沖擊過程，著陸器與著陸面的碰撞利用非線性阻尼彈簧模型模擬[4-5]。針對著陸器的一種特定構(gòu)型，為了提高仿真模型計(jì)算結(jié)果的準(zhǔn)確性，需要根據(jù)有限次的物理樣機(jī)試驗(yàn)結(jié)果，對碰撞模型的參數(shù)進(jìn)行修正。

目前，模型修正技術(shù)廣泛應(yīng)用于機(jī)械、建筑等領(lǐng)域。文獻(xiàn)[5]應(yīng)用重復(fù)仿真的方法，實(shí)現(xiàn)月球著陸器動(dòng)力學(xué)模型的修正，修正后模型仿真得到的著陸器整機(jī)質(zhì)心距翻倒墻最小距離與試驗(yàn)值更加貼近。文獻(xiàn)[6-10]通過構(gòu)造不同類型的響應(yīng)面，實(shí)現(xiàn)汽車制動(dòng)盤和橋梁等結(jié)構(gòu)有限元模型的修正，修正后的模型預(yù)測精度顯著提高。文獻(xiàn)[11]運(yùn)用響應(yīng)面法，實(shí)現(xiàn)數(shù)控機(jī)床主軸系統(tǒng)熱模型的多目標(biāo)修正。文獻(xiàn)[12-13]應(yīng)用Kriging代理模型，實(shí)現(xiàn)拱橋結(jié)構(gòu)和五層框架建筑結(jié)構(gòu)有限元模型的修正，驗(yàn)證了該方法的合理性。

有鑒于此，為了提升計(jì)算效率，同時(shí)保證模型修正的準(zhǔn)確性，本文應(yīng)用估計(jì)方差最小的無偏估計(jì)模型Kriging代理模型，針對著陸器停穩(wěn)后主體下降距離與投放高度之差、輔助支柱最大緩沖行程兩項(xiàng)著陸性能指標(biāo)值，實(shí)現(xiàn)著陸器單腿動(dòng)力學(xué)仿真模型參數(shù)修正。將修正后的參數(shù)代入仿真模型驗(yàn)證，兩項(xiàng)著陸性能指標(biāo)的仿真結(jié)果與物理樣機(jī)試驗(yàn)結(jié)果更加吻合。修正后的參數(shù)也可為著陸器整機(jī)動(dòng)力學(xué)仿真模型的碰撞參數(shù)提供參考。

1 新型著陸器動(dòng)力學(xué)建模

新型著陸器由主體結(jié)構(gòu)和四套對稱分布在主體周圍的緩沖機(jī)構(gòu)組成。主體結(jié)構(gòu)作為動(dòng)力學(xué)分析中的負(fù)載。四套緩沖機(jī)構(gòu)構(gòu)型與尺寸完全相同，均由主支柱、輔助支柱、足墊和緩沖桿組成，用以吸收著陸器與著陸面碰撞時(shí)產(chǎn)生的沖擊能量，其具體構(gòu)型如圖1所示。緩沖機(jī)構(gòu)中，主支柱起支撐作用；輔助支柱利用內(nèi)外筒之間的緩沖元件吸收沖擊能量；足墊用以增大著陸器與著陸面的接觸面積，防止著陸時(shí)緩沖機(jī)構(gòu)過度下陷；緩沖桿與主體固連，當(dāng)其受到?jīng)_擊時(shí)會發(fā)生彎曲變形，吸收沖擊能量，防止主體受力過大。

(a)側(cè)視圖(b)軸測圖

圖1 著陸緩沖機(jī)構(gòu)

Fig.1 Landing gear

為了驗(yàn)證緩沖機(jī)構(gòu)的著陸緩沖性能，有必要進(jìn)行著陸沖擊過程的仿真分析，考慮四套緩沖機(jī)構(gòu)完全相同，為了降低仿真分析的復(fù)雜度以提升分析效率，本文建立具有單套緩沖機(jī)構(gòu)的著陸器動(dòng)力學(xué)仿真模型進(jìn)行著陸沖擊的模擬[2]。仿真模型中主體上僅安裝一套緩沖機(jī)構(gòu)，著陸器整體只能沿重力方向運(yùn)動(dòng)，將模型中的足墊距著陸面的豎直高度(h)進(jìn)行參數(shù)化，可模擬多組著陸沖擊過程。模型示意圖，如圖2所示。

圖2 仿真模型示意圖

本文應(yīng)用多體動(dòng)力學(xué)分析軟件ADAMS建立單腿動(dòng)力學(xué)仿真剛體模型，并實(shí)現(xiàn)沖擊過程的仿真計(jì)算。其中，輔助支柱內(nèi)部的緩沖元件通過在輔助支柱內(nèi)外筒之間施加與內(nèi)、外筒相對位移(D)有關(guān)的力(FD)加以模擬，D與FD關(guān)系曲線如圖3所示。

圖3 緩沖元件的力學(xué)特性

考慮單腿動(dòng)力學(xué)仿真模型中緩沖桿發(fā)生側(cè)向彎曲的可能性很小，本文采用如圖4所示的方法對緩沖桿進(jìn)行剛體等效。

圖4 緩沖桿等效方法

首先，參照緩沖桿的實(shí)際質(zhì)量特性將其分為兩段剛性桿，桿1與主體之間添加轉(zhuǎn)動(dòng)副，即桿1僅可繞轉(zhuǎn)軸o轉(zhuǎn)動(dòng)，桿1和桿2之間添加移動(dòng)副，即桿2只能沿桿1的軸線做直線運(yùn)動(dòng)；之后，建立測量函數(shù)以實(shí)時(shí)測量桿1相對于轉(zhuǎn)軸o的轉(zhuǎn)角(θ)，并在桿1的轉(zhuǎn)軸位置施加與轉(zhuǎn)角(θ)相關(guān)的轉(zhuǎn)矩(M)，以模擬緩沖桿受彎時(shí)產(chǎn)生的緩沖彎矩，在桿2上施加與θ相關(guān)、沿桿1軸線方向的運(yùn)動(dòng)(S)，以保證等效模型受力時(shí)自由端位置變化與真實(shí)情況相符；之后，在有限元軟件Dytran中建立緩沖桿模型，施加轉(zhuǎn)矩，測得θ與M，θ與S關(guān)系曲線如圖5、圖6所示，導(dǎo)出數(shù)據(jù)文件，以Spline函數(shù)的形式賦給ADAMS模型中的M與S，最終實(shí)現(xiàn)緩沖桿剛體等效模型的建立。

圖5 θ-M關(guān)系曲線

圖6 θ-S關(guān)系曲線

著陸器與著陸面實(shí)際的碰撞過程十分復(fù)雜，為了便于進(jìn)行仿真計(jì)算，將足墊與著陸面的相互作用力分解為法向碰撞力(Fn)和切向摩擦力(Ff)，利用非線性阻尼彈簧碰撞模型描述Fn，表達(dá)式如下

(1)

式中：K，e，C為分別表示著陸面的剛度系數(shù)、非線性指數(shù)和阻尼系數(shù)；δ表示著陸面變形量。

Ff利用庫倫摩擦方程描述，表達(dá)式如下

Ff=μFn

(2)

式中：μ為足墊與著陸面間的摩擦因數(shù)。

根據(jù)參考文獻(xiàn)[5,14]并結(jié)合工程經(jīng)驗(yàn)設(shè)定K，e，C的初值，并確定鋁合金質(zhì)試驗(yàn)足墊與著陸面的摩擦因數(shù)μ，各參數(shù)取值如表1所示。

表1 參數(shù)初值

2 沖擊試驗(yàn)與仿真計(jì)算

2.1 沖擊試驗(yàn)

為了檢驗(yàn)仿真模型的準(zhǔn)確性，需進(jìn)行著陸器單腿樣機(jī)的沖擊試驗(yàn)。試驗(yàn)裝置參照仿真模型建立，主要由一套緩沖機(jī)構(gòu)、吊籃及配重、導(dǎo)軌、起吊裝置、釋放機(jī)構(gòu)和著陸面組成，如圖7所示。緩沖機(jī)構(gòu)與裝有配重的吊籃相連接，吊籃相當(dāng)于著陸器主體通過釋放機(jī)構(gòu)與起吊裝置相連，可沿導(dǎo)軌豎直滑動(dòng)。在地球重力試驗(yàn)環(huán)境下，通過改變吊籃的投放高度，可實(shí)現(xiàn)多次著陸沖擊試驗(yàn)。試驗(yàn)裝置裝有傳感器，以記錄緩沖機(jī)構(gòu)的各項(xiàng)性能指標(biāo)。

圖7 著陸沖擊試驗(yàn)裝置正視圖

試驗(yàn)序號投放高度/mHexpi/mDexpi/m10.30.0850.02020.50.1230.03230.70.1590.04440.90.1930.055510.2110.06161.10.2250.06671.20.2390.072

2.2 仿真計(jì)算及誤差分析

表3 仿真結(jié)果

為了定量地檢驗(yàn)仿真結(jié)果相較于試驗(yàn)結(jié)果的準(zhǔn)確程度，選取第1、3、4、5、7組試驗(yàn)結(jié)果與仿真結(jié)果，相應(yīng)計(jì)算出兩項(xiàng)性能指標(biāo)的相對誤差之和，如式(3)、式(4)所示。另外兩組試驗(yàn)結(jié)果用來驗(yàn)證后續(xù)模型修正結(jié)果的正確性。

(3)

(4)

模型修正前，f1，f2計(jì)算結(jié)果如表4所示。

表4 相對誤差之和

3 碰撞模型參數(shù)修正

3.1 Kriging代理模型的建立

在著陸器構(gòu)型參數(shù)確認(rèn)無誤的情況下，碰撞模型中剛度系數(shù)(K)，阻尼系數(shù)(C)及非線性指數(shù)(e)的取值將極大地影響仿真結(jié)果。上述參數(shù)的初值通常根據(jù)工程經(jīng)驗(yàn)給出，與真實(shí)的著陸面參數(shù)存在一定偏差，因此本文選取上述參數(shù)作為模型修正參數(shù)。根據(jù)工程經(jīng)驗(yàn)并結(jié)合參考文獻(xiàn)[5,14]，確定各參數(shù)的范圍如表5所示。

表5 模型修正參數(shù)取值范圍

為了提高計(jì)算效率，本文基于Kriging代理模型實(shí)現(xiàn)仿真模型的修正。采用優(yōu)化拉丁超立方實(shí)驗(yàn)設(shè)計(jì)在表4所示各參數(shù)的取值范圍內(nèi)抽取40個(gè)樣本點(diǎn)，仿真計(jì)算出各樣本點(diǎn)對應(yīng)的著陸器性能指標(biāo)值，再由式(3)、式(4)計(jì)算出f1與f2。將修正參數(shù)值作為輸入，f1與f2作為響應(yīng)，構(gòu)造Kriging代理模型。

本文構(gòu)造的Kriging代理模型由多項(xiàng)式回歸函數(shù)和隨機(jī)過程函數(shù)組成，其表示如下。

y(x)=fT(x)β+z(x)

(5)

式中:x表示修正參數(shù)向量；y(x)表示Kriging代理模型的響應(yīng)函數(shù)，f(x)表示多項(xiàng)式回歸函數(shù)向量；β表示回歸系數(shù)向量；z(x)表示均值為0、方差為σ2的隨機(jī)過程函數(shù)，代理模型詳細(xì)構(gòu)建過程詳見文獻(xiàn)[15]，在此不做贅述。

利用式(6)所示的決定系數(shù)R2和式(7)所示的均方根誤差相對值RMSE檢驗(yàn)Kriging代理模型對仿真結(jié)果的擬合精度。

(6)

(7)

表6 Kriging代理模型誤差分析

3.2 參數(shù)修正

xL

(8)

式中：x=(K，C，e)T，代表設(shè)計(jì)變量；x(L)與x(U)分別表示設(shè)計(jì)變量的下界與上界。

本文利用第二代非劣排序遺傳算法(NSGA-Ⅱ)實(shí)現(xiàn)優(yōu)化計(jì)算，算法參數(shù)設(shè)置如表7所示。

表7 優(yōu)化參數(shù)設(shè)置

綜上，單腿動(dòng)力學(xué)仿真模型修正流程如圖8所示。

經(jīng)過優(yōu)化計(jì)算，得到帕累托最優(yōu)解集如表8所示。以第1組解為例進(jìn)行分析，將該組解代入動(dòng)力學(xué)仿真模型，計(jì)算得到f1與f2的值，并與修正前結(jié)果對比，如表9所示。

表8 帕累托最優(yōu)解集

修正后模型仿真計(jì)算出著陸器的兩項(xiàng)性能指標(biāo)值如表10所示。

表9 修正前后參數(shù)對比

表10 修正后模型仿真結(jié)果

針對用于修正模型的5組試驗(yàn)高度，將模型修正前后的仿真結(jié)果與試驗(yàn)結(jié)果做柱狀圖計(jì)較如圖9、圖10所示。

圖9 主體下降距離與投放高度之差結(jié)果對比

圖10 輔助支柱最大緩沖行程結(jié)果對比

針對模型修正過程中未用到的兩組試驗(yàn)高度，將模型修正前后的仿真結(jié)果與試驗(yàn)結(jié)果做柱狀圖計(jì)較如圖11、圖12所示。

圖11 主體下降距離與投放高度之差結(jié)果對比

圖12 輔助支柱最大緩沖行程結(jié)果對比

針對圖11、圖12所示的兩組試驗(yàn)高度，計(jì)算得到著陸器主體下降距離與投放高度之差平均仿真誤差由修正前的29.4%降低到5.9%，輔助支柱最大緩沖行程平均仿真誤差由修正前的24.7%降低到8.4%。

4 結(jié) 論

(1) 本文建立了新型著陸器的單腿動(dòng)力學(xué)仿真模型，著陸器與著陸面的碰撞采用非線性阻尼彈簧模型模擬。

(2) 為了高效地評估仿真模型的精度，針對模型中的碰撞參數(shù)，結(jié)合5組物理樣機(jī)的試驗(yàn)結(jié)果，構(gòu)造出著陸器兩項(xiàng)性能指標(biāo)相對誤差之和的Kriging代理模型。

(3) 以最小化代理模型計(jì)算出的性能指標(biāo)相對誤差之和作為目標(biāo)，應(yīng)用NSGA-Ⅱ算法尋優(yōu)，實(shí)現(xiàn)碰撞參數(shù)修正。將修正后的參數(shù)代入動(dòng)力學(xué)仿真模型，另取2組試驗(yàn)結(jié)果驗(yàn)證，模型計(jì)算出的著陸器主體下降距離與投放高度之差誤差由修正前的29.4%降低到5.9%，輔助支柱最大緩沖行程誤差由修正前的24.7%降低到8.4%。

[1] 楊建中,滿劍鋒,曾福明,等.“嫦娥三號”著陸緩沖機(jī)構(gòu)的研究成果及其應(yīng)用[J].航天返回與遙感，2014,35(6):20-27.

YANG Jianzhong, MAN Jianfeng, ZENG Fuming, et al. The landing process of Mars phoenix lander entry descent, and landing[J].Spacecraft Recovery & Remote Sensing，2014,35(6):20-27.

[2] 曾福明,楊建中,滿劍峰，等.月球著陸器著陸緩沖機(jī)構(gòu)設(shè)計(jì)方法研究[J].航天器工程,2011,20(2):46-51.

ZENG Fuming, YANG Jianzhong, MAN Jianfeng, et al.Study of design method of landing gear for lunar lander[J].Spacecraft Engineering，2011,20(2):46-51.

[3] 陳金寶,聶宏,萬峻麟.深空探測著陸器數(shù)字化設(shè)計(jì)及著陸性能影響因素[J].航空學(xué)報(bào),2014,35(2):541-554.

CHEN Jinbao, NIE Hong, WAN Junlin.Digital design and landing performance influence factors of deep space lander[J].Acta Aeronautica et Astronautica Sinica，2014,35(2):541-554.

[4] NOHMIN M, MIYAHARA A. Modeling for lunar lander by mechanical dynamics software[R]. AIAA-2005-6416, 2005.

[5] 吳建云,王春潔,汪瀚.月球著陸器模型參數(shù)修正[J].北京航空航天大學(xué)學(xué)報(bào)，2013,39(10):1366-1375.

WU Jianyun, WANG Chunjie, WANG Han. Model parameters updating of Lunar lander[J].Journal of Beijing University of Aeronautics and Astronautics，2013,39(10):1366-1375.

[6] GUO B Y, ZHANG W J, LU Q H,et al.A successive selection method for finite element model updating[J].Mechanical Systems and Signal Processing，2016,70/71:320-333.

[7] NIU J, ZONG G H, CHU F P,et al. Damage identification method of girder bridges based on finite element model updating and modal strain energy[J]. Science China Technological Sciences，2013,28(3):210-226.

[8] ZHOU L R，YAN G R，OU J P.Response surface method based on radial basis functions for modeling large-scale structures in model updating[J]. Computer-Aided Civil and Infrastructure Engineering，2015,58(4):701-711.

[9] REN W X,FANG S E,DENG M Y. Finite element model updating in structural dynamics by using the response surface method[J].Journal of Mechanical Engineering，2011,137(4):248-257.

[10] CHAKRABORTY S,SEN A. Adaptive response surface based efficient finite element model updating[J]. Finite Elements in Analysis & Design，2014，80(3):33-405.

[11] 鄧小雷,傅建中,夏晨輝,等.數(shù)控機(jī)床主軸系統(tǒng)熱模型參數(shù)多目標(biāo)修正方法[J].機(jī)械工程學(xué)報(bào),2014,50(15):121-126.

DENG Xiaolei, FU Jianzhong, XIA Chenhui, et al. Multi-objective correction method for thermal model parameters of CNC machine tool spindle system[J].Journal of Mechanical Engineering，2014,50(15):121-126.

[12] 胡俊亮,顏全勝,鄭恒斌,等.基于Kriging模型的鋼管混凝土連續(xù)梁拱橋有限元模型修正[J].振動(dòng)與沖擊,2014,33(14):33-39.

HU Junliang, YAN Quansheng, ZHENG Hengbin,et al. CFST arch/continuous beam bridge FEM model updating based on Kriging model[J].Journal of Vibration and Shock，2014,33(14):33-39.

[13] JIN S S, JUNG H J. Sequential surrogate modeling for efficient finite element model updating[J].Computers and Structures，2016,168: 30-45.

[14] 吳曉君,鐘世英,凌道盛,等.數(shù)著陸器足墊垂直沖擊模型試驗(yàn)研究[J].巖土力學(xué),2012,33(4):1045-1050.

WU Xiaojun, ZHONG Shiying, LING Daosheng, et al. Model test study of vertical impact of space lander footpad[J].Rock and Soil Mechanics，2012,33(4):1045-1050.

[15] WANG J T, WANG C J, ZHAO J P. Structural dynamic model updating based on Kriging model using frequency response date[J]. Journal of Vibroengineering，2016,18(6):3484-3498.