月球著陸器著陸腿非線性有限元建模與仿真

梁東平 柴洪友 曾福明

(北京空間飛行器總體設(shè)計(jì)部，北京100094)

月球著陸器著陸緩沖分系統(tǒng)的主要功能是緩沖著陸器月面軟著陸時(shí)的沖擊載荷，保證結(jié)構(gòu)和有效載荷不因著陸沖擊而破壞，并吸收著陸過(guò)程中的動(dòng)能，使著陸器安全著陸到月面上.由于著陸沖擊試驗(yàn)成本高、周期長(zhǎng)、難度大，只能進(jìn)行有限次數(shù)的特定工況地面模擬驗(yàn)證試驗(yàn)，因此必須開展大量著陸沖擊動(dòng)力學(xué)仿真分析，驗(yàn)證著陸緩沖分系統(tǒng)的性能，獲得相關(guān)性能參數(shù)，從而為著陸緩沖分系統(tǒng)的優(yōu)化設(shè)計(jì)和試驗(yàn)驗(yàn)證提供指導(dǎo)，并為著陸器結(jié)構(gòu)和儀器設(shè)備的設(shè)計(jì)提供輸入條件［1-4］.

著陸緩沖分系統(tǒng)著陸沖擊仿真大多數(shù)是基于多剛體動(dòng)力學(xué)模型，多剛體模型比較簡(jiǎn)單，計(jì)算效率高，但沒(méi)有考慮結(jié)構(gòu)柔性［5-7］.為研究結(jié)構(gòu)柔性對(duì)著陸腿性能的影響，通常采用部件模態(tài)綜合的方法.通過(guò)有限元分析軟件生成著陸腿的模態(tài)中性文件，并導(dǎo)入多體動(dòng)力學(xué)軟件中，建立柔性體模型［8］.由于著陸沖擊過(guò)程中涉及材料、幾何和接觸等非線性問(wèn)題，采用多體動(dòng)力學(xué)建模方法很難將這些非線性因素考慮在內(nèi).例如考慮結(jié)構(gòu)柔性時(shí)內(nèi)外套筒接觸面法向力很復(fù)雜，而多體動(dòng)力學(xué)模型內(nèi)外套筒之間往往采用滑動(dòng)副的約束形式，因此很難準(zhǔn)確地計(jì)算著陸沖擊過(guò)程中內(nèi)外套筒之間的摩擦力.著陸腿多柔性體動(dòng)力學(xué)模型的柔性信息通常是以沖擊前著陸腿的構(gòu)型和邊界條件給出的，在著陸沖擊過(guò)程中，著陸腿支柱內(nèi)外套筒會(huì)產(chǎn)生壓縮行程，著陸腿構(gòu)型和邊界條件是不斷變化的，用著陸沖擊前的構(gòu)型和邊界條件生成的模態(tài)中性文件來(lái)分析著陸沖擊過(guò)程也會(huì)產(chǎn)生一定的模型誤差.此外，多體動(dòng)力學(xué)模型中著陸腿足墊與月壤的接觸力只能根據(jù)具體工況采用經(jīng)驗(yàn)公式計(jì)算，通用性較差.

本文以非線性有限元理論為基礎(chǔ)，充分考慮著陸器著陸腿在著陸沖擊過(guò)程中的幾何、材料和接觸等非線性因素，建立了詳細(xì)的著陸腿有限元模型，并采用ABAQUS/Explicit作為求解器進(jìn)行著陸沖擊非線性動(dòng)力學(xué)分析，并與試驗(yàn)數(shù)據(jù)進(jìn)行了對(duì)比.仿真結(jié)果表明，非線性有限元?jiǎng)恿W(xué)模型能較準(zhǔn)確地預(yù)測(cè)著陸腿的著陸緩沖性能.

1 著陸腿的構(gòu)型及緩沖原理

月球著陸器多采用腿式著陸緩沖機(jī)構(gòu)，著陸腿圍繞著陸器中心軸均勻分布，如圖1所示.單個(gè)著陸腿緩沖機(jī)構(gòu)由主支柱、輔助支柱和足墊組成.主支柱外筒上端通過(guò)萬(wàn)向節(jié)與著陸器主體連接，主支柱內(nèi)筒下端通過(guò)球鉸與足墊連接.輔助支柱通過(guò)萬(wàn)向節(jié)著陸器主體連接，通過(guò)球鉸與主支柱連接.主、輔支柱內(nèi)都安裝有鋁蜂窩芯子緩沖器，利用蜂窩芯子塑性壓縮變形吸收著陸沖擊過(guò)程中的動(dòng)能，起到緩沖著陸沖擊載荷的作用.主支柱內(nèi)的緩沖器具有單向緩沖功能，輔助支柱緩沖器具有雙向緩沖功能.

圖1 軟著陸機(jī)構(gòu)示意圖

2 顯式非線性有限元分析理論

顯式非線性有限元?jiǎng)恿W(xué)模型的數(shù)值解法，是將結(jié)構(gòu)在空間域進(jìn)行離散，把連續(xù)的微分方程轉(zhuǎn)換成有限階的代數(shù)方程組:

對(duì)求解的時(shí)間歷程在時(shí)間域內(nèi)離散，并采用顯式積分方法進(jìn)行求解.中心差分法是最常用的顯式積分方法，積分步驟如下:

其中，u·為速度列陣;u為位移列陣;上標(biāo)(i)表示增量步的次數(shù)，(i-1/2)和(i+1/2)是指增量步的中間值.這樣，在增量步開始時(shí)提供了滿足動(dòng)力學(xué)平衡條件的加速度.得到了加速度后，在時(shí)間上“顯式地”前推速度和位移.顯式積分算法中，不需要形成整體的剪切剛度矩陣，由于是顯式的前推模型的狀態(tài)，所以不需要迭代和收斂準(zhǔn)則.顯式動(dòng)力學(xué)求解方法是條件穩(wěn)定的，最大穩(wěn)定時(shí)間積分步長(zhǎng)Δts與單元的特征尺寸Le和波速cd相關(guān):

對(duì)于線彈性材料，cd的計(jì)算公式如下:

其中，E為彈性模量;ν為泊松比;ρ為密度.

3 著陸沖擊有限元建模

3.1 著陸腿有限元模型

在ABAQUS/CAE平臺(tái)上建立著陸腿有限元模型，模型主要包括鋁蜂窩緩沖器、主支柱、輔助支柱和足墊，如圖2所示.著陸腿支柱內(nèi)外套筒采用S4R減縮積分殼單元及S3殼單元建模，并對(duì)不同區(qū)域賦予不同的單元截面屬性(section property).足墊采用解析剛體建模.輔助支柱和主支柱之間、主支柱與足墊之間采用球鉸連接單元(connector element)連接.

著陸腿中對(duì)仿真結(jié)果影響不大的細(xì)節(jié)已經(jīng)做了簡(jiǎn)化處理.由于著陸腿一些局部區(qū)域的單元尺寸較小，這些單元控制最大穩(wěn)定時(shí)間增量步長(zhǎng)，對(duì)這些單元較小的區(qū)域采用質(zhì)量縮放技術(shù)，將這些局部區(qū)域單元的質(zhì)量增加5倍，人為地增大材料的密度，由式(4)、式(5)可知，穩(wěn)定時(shí)間步長(zhǎng)可增加約2倍.這些局部區(qū)域質(zhì)量的增大不會(huì)對(duì)著陸腿系統(tǒng)的總體動(dòng)力學(xué)響應(yīng)產(chǎn)生太大影響，但卻可以明顯地提高計(jì)算效率，減少計(jì)算量.

圖2 著陸腿有限元模型

在著陸沖擊過(guò)程中，著陸腿各支柱內(nèi)外套筒接觸面之間以及足墊與著陸面之間的相互作用采用接觸算法來(lái)實(shí)現(xiàn).模型采用ABAQUS/Explicit中的通用接觸算法定義著陸腿構(gòu)件間的接觸，同時(shí)對(duì)不同的接觸區(qū)域定義不同的接觸屬性.采用庫(kù)侖摩擦模型來(lái)描述接觸面間摩擦的相互作用.摩擦系數(shù)定義如下:

其中，μ為摩擦系數(shù);μk為動(dòng)摩擦系數(shù);μs為靜摩擦系數(shù);dc為指數(shù)衰減系數(shù);γ·eq為接觸面間的相對(duì)滑移速度.模型中著陸腿支柱內(nèi)外套筒之間的摩擦參數(shù)取值均由試驗(yàn)標(biāo)定.

3.2 鋁蜂窩緩沖器模型

采用Crushable Foam體積硬化模型作為鋁蜂窩緩沖器的材料模型.根據(jù)鋁蜂窩緩沖器的額定壓縮載荷行程特性，定義材料模型的彈性及硬化參數(shù).初始屈服應(yīng)力由鋁蜂窩的壓縮載荷和鋁蜂窩芯子的橫截面積確定.將弱蜂窩極限壓縮塑性應(yīng)變對(duì)應(yīng)的屈服應(yīng)力設(shè)為一個(gè)大值，當(dāng)弱蜂窩達(dá)到極限壓縮行程后將會(huì)硬化，從而促使強(qiáng)蜂窩開始?jí)嚎s.

為了驗(yàn)證模型的正確性，本文以主支柱二級(jí)鋁蜂窩緩沖器的載荷行程特性進(jìn)行了仿真分析，分析模型如圖3所示.底部剛性板固支，蜂窩緩沖器固定在底板上，頂部剛性板以4 m/s的速度向下運(yùn)動(dòng)，達(dá)到蜂窩的極限壓縮行程時(shí)停止.底部剛性板的支反力(80 Hz濾波)以及鋁蜂窩的壓縮行程曲線如圖4所示.曲線中較平滑的兩段對(duì)應(yīng)的壓縮載荷與設(shè)計(jì)值基本一致.整個(gè)壓縮過(guò)程中鋁蜂窩的總塑性耗散能為7 548 J，與額定值7 520 J相當(dāng).仿真結(jié)果證明，該模型能很好地模擬蜂窩緩沖器的緩沖及吸能特性.

圖3 蜂窩緩沖器有限元模型

圖4 鋁蜂窩緩沖器載荷行程曲線

4 單腿著陸沖擊仿真及試驗(yàn)驗(yàn)證

單腿著陸沖擊試驗(yàn)裝置如圖5所示.吊籃用于模擬著陸腿的真實(shí)安裝情況，吊籃安裝在導(dǎo)軌上，只能沿導(dǎo)軌在豎直方向運(yùn)動(dòng)，吊籃上安裝有配重用于模擬不同的著陸工況.試驗(yàn)時(shí)，由起吊裝置將吊籃及著陸腿吊至指定高度后由釋放裝置釋放，吊籃及著陸腿在重力作用下沖擊到著陸平臺(tái)上.由安裝在吊籃上的加速度儀測(cè)量吊籃中心處的加速度.本文對(duì)表1中的著陸工況進(jìn)行仿真分析，以驗(yàn)證著陸腿模型的正確性.

圖5 單腿著陸沖擊試驗(yàn)示意圖［7］

表1 著陸沖擊試驗(yàn)工況

單腿著陸沖擊有限元仿真模型如圖6所示.模型由吊籃、著陸腿和著陸面組成.由于吊籃為鋼結(jié)構(gòu)，其剛度比著陸腿大，且不是仿真研究所關(guān)注的重點(diǎn)，因此采用離散剛體模型建模，其邊界條件施加在離散剛體的質(zhì)心參考點(diǎn)上，約束參考點(diǎn)上的3個(gè)旋轉(zhuǎn)自由度和y向、z向的平移自由度.著陸腿支柱和吊籃之間采用萬(wàn)向節(jié)連接單元連接.著陸平臺(tái)分為剛性月面和模擬月壤兩種情況.采用Cap Drucke-Prager(CDP)模型作為模擬月壤的材料模型，該模型是在子午線為線性的Drucker-Prager模型上加上一個(gè)帽蓋屈服面形成的，這個(gè)帽蓋屈服面設(shè)置了靜水壓力的上限，因此提供了一個(gè)非彈性硬化機(jī)制用于描述材料的塑性壓縮.模型參數(shù)由JLU2模擬月壤［9］的三軸壓縮試驗(yàn)數(shù)據(jù)及其相關(guān)研究成果標(biāo)定，具體的模型參數(shù)如表2所示.表中，d為材料粘性;β為材料摩擦角;R為帽蓋偏心率;p0為帽蓋初始屈服面位置;r為過(guò)度面半徑參數(shù);K為流動(dòng)應(yīng)力比.模擬月壤的壓縮系數(shù)Cc和膨脹系數(shù)Cs分別取0.02和0.005，月壤的壓縮硬化特性由壓縮系數(shù)、膨脹系數(shù)和初始孔隙率定義.

圖6 單腿著陸沖擊仿真模型

表2 模擬月壤材料屬性

著陸沖擊仿真結(jié)果與試驗(yàn)結(jié)果的對(duì)比如表3所示，著陸腿主支柱和輔助支柱的緩沖器的緩沖行程仿真結(jié)果與試驗(yàn)結(jié)果基本一致，兩個(gè)工況的誤差都在工程上可接受的范圍內(nèi).這說(shuō)明非線性有限元模型可以較準(zhǔn)確地預(yù)測(cè)著陸腿的緩沖和吸能特性.工況2的主支柱行程仿真值要比工況1小，這是因?yàn)槟M月壤的塑性變形會(huì)吸收一部分動(dòng)能;而工況2主支柱行程試驗(yàn)值卻比工況1大，原因主要有以下兩方面:蜂窩芯子由于存在制造誤差，工況2所用蜂窩芯子實(shí)際壓縮載荷可能要比工況1所用的蜂窩芯子小;試驗(yàn)中工況1所用鋁蜂窩芯子存在壓偏的現(xiàn)象，這會(huì)增大著陸腿套筒間摩擦力，使得工況1的壓縮行程偏小.

為了研究著陸腿的柔性對(duì)著陸沖擊性能的影響，在工況1條件下［8］，對(duì)著陸腿施加剛性約束進(jìn)行著陸沖擊仿真分析，吊籃中心加速度(80 Hz濾波)時(shí)間歷程如圖7所示.吊籃中心加速度試驗(yàn)峰值約為9.18 g，柔性著陸腿模型和剛性著陸腿模型仿真結(jié)果分別為8.0 g和7.67 g.由于試驗(yàn)過(guò)程中加速度測(cè)量點(diǎn)的測(cè)試結(jié)果受安裝剛度影響較大，而仿真模型的吊籃采用剛體建模，導(dǎo)致加速度試驗(yàn)值和仿真值匹配性不是很好，但總體變化趨勢(shì)基本一致.同時(shí)由圖7可知，柔性腿模型比剛性腿模型的仿真值更接近試驗(yàn)測(cè)試值.

柔性腿和剛性腿主支柱內(nèi)外套筒之間的摩擦力(80 Hz濾波)如圖8所示.柔性腿套筒摩擦力峰值約為5.4 kN，剛性腿峰值約為4.36 kN.剛性腿主支柱套筒間的摩擦力變化比較平緩，柔性腿主支柱摩擦力的波動(dòng)性較大，這主要是主支柱本身，以及輔助支柱的彈性變形引起的.在著陸沖擊前0.15 s過(guò)程中，柔性腿摩擦力曲線大部分都位于剛性腿摩擦力曲線上方，可見(jiàn)著陸腿的柔性變形會(huì)增大套筒間的摩擦力.

工況2著陸沖擊仿真能量時(shí)間歷程如圖9所示.在著陸沖擊后0.15 s時(shí)刻，系統(tǒng)動(dòng)能變化量、重力做功及人工應(yīng)變能總和為4 135 J，塑性變形和摩擦力耗散能總和約為4 077 J，相差的58 J由數(shù)值求解過(guò)程中引入的粘性耗散掉，由此可見(jiàn)，沖擊過(guò)程中能量的變化符合能量守恒定律.人工應(yīng)變能很小說(shuō)明仿真模型沒(méi)有發(fā)生明顯的沙漏現(xiàn)象，仿真結(jié)果是可靠的.系統(tǒng)著陸時(shí)刻理論初始動(dòng)能為3367 J，仿真分析初始時(shí)刻動(dòng)能為3384 J，這是由于采用了局部質(zhì)量放大技術(shù)引起的，這個(gè)微小的差異對(duì)著陸腿的整體沖擊動(dòng)力學(xué)響應(yīng)的影響不大，但使穩(wěn)定時(shí)間步長(zhǎng)由原來(lái)7.08×10-7s增大到1.46×10-6s，顯著地提高了計(jì)算效率.

表3 仿真結(jié)果與試驗(yàn)結(jié)果對(duì)比

圖7 單腿著陸沖擊仿真模型

圖8 主支柱內(nèi)外套筒間的摩擦力時(shí)間歷程

圖9 工況2能量時(shí)間歷程

5 結(jié)論

1)非線性有限元?jiǎng)恿W(xué)建模方法可以充分考慮著陸腿著陸沖擊過(guò)程中的材料、幾何、接觸等非線性因素，能較準(zhǔn)確地預(yù)測(cè)著陸腿的載荷緩沖和吸能特性，可用于指導(dǎo)著陸腿設(shè)計(jì)及優(yōu)化.

2)著陸腿柔性對(duì)其內(nèi)外套筒之間的摩擦力有較大的影響，分析時(shí)應(yīng)予以考慮;增加著陸腿剛度可以減小摩擦力.

3)采用局部質(zhì)量縮放技術(shù)可以顯著地增大顯式動(dòng)力學(xué)程序的穩(wěn)定時(shí)間增量步長(zhǎng)，提高計(jì)算效率，且不會(huì)對(duì)總體動(dòng)態(tài)響應(yīng)產(chǎn)生太大影響.

References)

［1］ William F R.Apollo experience report-lunar module landing gear subsystem［R］.NASA TN D-6850，1972

［2］ Merchant D H，Sawdy D T.Monte Carlo dynamic analysis for lunar module landing loads［J］.Journal of Spacecraft and Rockets，1970，8(1):48-55

［3］ Doiron H H，George A Z.Apollo lunar module landing dynamics［C］//The 41st Structures，Structural Dynamics，and Materials Conferenceand Exhibit.Washington:AIAA/ASME/ASCE/AHS/ASC，2000

［4］ Andrew L B，James R C，Ralph D L，et al.Planetary landers and entry probes［M］.New York:Cambridge University Press，2007

［5］ Masahiro N，Akira M.Modeling for lunar lander by mechanical dynamics software［C］//Modeling and Simulation Technologies Conference and Exhibit.Washington:AIAA，2005

［6］ Hilchenbach M，Kuchemann O，Rosenbauer H.Impact on a comet:Rosetta lander simulations［J］.Planetary and Space Science，2000，48(5):361-369

［7］朱汪，楊建中.月球探測(cè)器軟著陸機(jī)構(gòu)著陸腿模型與仿真分析［J］.宇航學(xué)報(bào)，2008，29(6):1723-1728 Zhu Wang，Yang Jianzhong.Modeling and simulation of landing leg for lunar landing gear system ［J］.Journal of Astronautics，2008，29(6):1723-1728(in Chinese)

［8］ 曾福明，楊建中，朱汪，等.月球著陸器著陸緩沖性能研究［J］.航天器工程，2010，19(5):43-49 Zeng Fuming，Yang Jianzhong，Zhu Wang，et al.Research on landing impact attenuation performance of lunar lander［J］.Spacecraft Engineering，2010，19(5):43-49(in Chinese)

［9］楊艷靜.模擬月壤和月球車車輪相互作用的數(shù)值仿真和試驗(yàn)研究［D］.北京:中國(guó)空間技術(shù)研究院，2010 Yang Yanjing.Numerical and experimental research of lunar soil stimulant and lunar rover wheel interaction［D］.Beijing:China Academy of Space Technology，2010(in Chinese)