0 引言

月球著陸器研發(fā)過程中需要準(zhǔn)確地預(yù)測(cè)月面軟著陸沖擊載荷，從而為結(jié)構(gòu)設(shè)計(jì)提供載荷輸入條件。著陸器著陸沖擊載荷條件的確定有兩種方法。一種方法是將所有可能的輸入?yún)?shù)的離散值進(jìn)行排列組合，求出所有組合工況條件下結(jié)構(gòu)所受到的最大沖擊載荷作為設(shè)計(jì)載荷。采用這種方法的缺點(diǎn)是求出的設(shè)計(jì)載荷可能過于保守。對(duì)于著陸沖擊載荷分析這樣的多工況問題，這種方法存在一個(gè)難以克服的困難，即所需的仿真量過大。如果以 3 個(gè)著陸線速度、3 個(gè)著陸角速度、3 個(gè)著陸姿態(tài)角和著陸面坡度 10 個(gè)變量作為仿真輸入?yún)?shù)，每個(gè)變量取 3 個(gè)值，可能的組合工況數(shù)高達(dá) 310個(gè)，進(jìn)行這樣大規(guī)模工況的仿真是不現(xiàn)實(shí)的。另一方法叫蒙特卡洛方法，即根據(jù)各參數(shù)的概率密度函數(shù)隨機(jī)選擇輸入?yún)?shù)，根據(jù)選擇的輸入?yún)?shù)樣本計(jì)算著陸沖擊載荷，然后采用統(tǒng)計(jì)方法確定著陸沖擊載荷可能的最大值，并給出相應(yīng)的置信水平，這種方法可以大大降低所需的仿真工況數(shù)。蒙特卡洛法在 Surveyor，Apollo，Viking等行星探測(cè)器的設(shè)計(jì)中已經(jīng)得到應(yīng)用[1-7]，但由于受到當(dāng)時(shí)計(jì)算速度的限制，多采用解析的剛體模型計(jì)算著陸沖擊載荷。隨著計(jì)算機(jī)計(jì)算速度的不斷提高和有限元方法的發(fā)展成熟，使得采用有限元柔性體模型進(jìn)行著陸沖擊載荷的分析計(jì)算成為了可能。

 

本文通過參數(shù)化建模系統(tǒng)建立著陸器有限元模型，并采用具有強(qiáng)大非線性瞬態(tài)計(jì)算功能的 Abaqus/Explicit 作為求解器。通過與試驗(yàn)數(shù)據(jù)的對(duì)比驗(yàn)證和修正非線性有限元模型。根據(jù)

輸入?yún)?shù)的概率密度函數(shù)隨機(jī)生成 100 個(gè)著陸工況，然后利用經(jīng)試驗(yàn)驗(yàn)證的有限元模型對(duì)這100 個(gè)著陸工況開展著陸沖擊仿真分析，并通過 3 法則確定著陸器著陸沖擊的設(shè)計(jì)載荷。

 

1 顯式非線性有限元?jiǎng)恿W(xué)分析理論

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

 

 

 

其中 為加速度列陣， 為質(zhì)量矩陣， 是單元應(yīng)力引起的內(nèi)部節(jié)點(diǎn)力列陣， 為外載荷列陣。

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

 

 

 

其中 為速度列陣； 為位移列陣；上標(biāo)(i)表示增量步的次數(shù)，(i-1/2)和(i+1/2)是指增量步的中間值。這樣，在增量步開始時(shí)提供了滿足動(dòng)力學(xué)平衡條件的加速度。得到了加速度后，在時(shí)間上“顯式地”前推速度和位移。顯式積分算法中，不需要形成整體的剪切剛度矩陣。由于是顯式的前推模型的狀態(tài)，所以不需要迭代和收斂準(zhǔn)則。顯式方法中可以很容易地模擬接觸條件和其他一些極度不連續(xù)的情況，通過調(diào)整節(jié)點(diǎn)加速度來平衡在接觸時(shí)的外力和內(nèi)力。這些特點(diǎn)使得顯式方法特別適合于求解著陸沖擊瞬態(tài)動(dòng)力學(xué)問題。雖然顯式動(dòng)力學(xué)求解方法的效率很高，但它卻是條件穩(wěn)定的，最大穩(wěn)定時(shí)間積分步長(zhǎng) 與單元的特征尺寸和波速 相關(guān)：

 

 

 

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

 

 

 

式中： 為彈性模量； 為泊松比； 為密度。

 

2 著陸器有限元模型及參數(shù)化建模系統(tǒng)

2.1 著陸器有限元模型

采用 ABAQUS/CAE 建立的著陸器著陸沖擊有限元模型及坐標(biāo)系如圖 1 所示，主要由著陸器模擬結(jié)構(gòu)、四套著陸腿（包括鋁蜂窩緩沖元件、拉桿緩沖元件、足墊等）和著陸面組成。著陸腿主支柱外筒及輔助支柱外筒頂端通過萬向節(jié)(UJOIN)連接單元連接；輔助支柱內(nèi)筒底

端與主支柱外筒之間、主支柱內(nèi)筒與足墊之間通過球鉸(JOIN+ROTATION)連接單元連接。著陸器足墊采用離散剛體建模。模擬結(jié)構(gòu)主要由鋼板和槽鋼組成，鋼板部分網(wǎng)格較粗，故采用全積分殼單元建模，而槽鋼結(jié)構(gòu)采用梁?jiǎn)卧?。著陸腿?nèi)外套筒采用殼單元建模。鋁蜂窩緩沖元件采用實(shí)體單元建模，采用 Crushable Foam 體積硬化模型作為鋁蜂窩緩沖器的材料模型，具體參數(shù)見文獻(xiàn)[8]。著陸面模型由一個(gè) 7m×7m×1m 實(shí)體單元集合組成，采用Drucke-Prager/Cap 模型作為模擬月壤的本構(gòu)模型，模型參數(shù)見文獻(xiàn)。

 

 

圖 1 著陸器有限元模型

Fig.1 Finite element model of lunar lander

 

2.2 著陸沖擊仿真分析及試驗(yàn)驗(yàn)證

著陸仿真分析的工況見表 1。著陸器的偏航角定義為繞 X 軸轉(zhuǎn)動(dòng)的角度，滾轉(zhuǎn)角為繞 Z軸轉(zhuǎn)動(dòng)的角度，俯仰角為繞 Y 軸轉(zhuǎn)動(dòng)的角度，三者均按右手定則方向?yàn)檎?。試?yàn)工況著陸器的初始角速度、俯仰角、滾轉(zhuǎn)角、偏航角都近似為 0。為了與試驗(yàn)結(jié)果進(jìn)行對(duì)比，模型中重力場(chǎng)取 9.8 m/s2。試驗(yàn)時(shí)模擬著陸器垂直方向的速度是在重力加速度作用下作自由落體運(yùn)動(dòng)獲得的，為了節(jié)省計(jì)算時(shí)間，將模擬著陸器置于一個(gè)十分接近著陸面的位置，然后賦予初始速度來等效其在自由落體下的運(yùn)動(dòng)效果。著陸面模型底部和四周固支。

 

表 1 著陸工況

Tab.1 Landing scenarios

 

 

著陸緩沖裝置主支柱接頭載荷分析值與試驗(yàn)結(jié)果對(duì)比見表 2,分析值與試驗(yàn)實(shí)測(cè)值的最大相對(duì)誤差約為 8.97%，平均相對(duì)誤差約為 5.18%，分析精度滿足工程要求。仿真結(jié)果與試驗(yàn)測(cè)量值的誤差來源主要有以下幾個(gè)方面：1）真實(shí)緩沖器的壓縮載荷與有限元模型中設(shè)定的設(shè)計(jì)值有偏差；2）著陸腿足墊與模擬月壤之間的摩擦系數(shù)與模型中設(shè)定的固定摩擦系數(shù)存在偏差；3）力傳感器的測(cè)量值本身也存在一定的測(cè)量誤差；4）有限元模型做了一些理想化的假設(shè)和簡(jiǎn)化，與真實(shí)試驗(yàn)樣機(jī)之間也存在一定的模型誤差。

 

表 2 分析與試驗(yàn)結(jié)果對(duì)比

Tab.2 Comparison of simulation and test results

 

 

2.3 參數(shù)化建模系統(tǒng)

月球著陸器有限元建模涉及幾何、材料和接觸等非線性因素，模型比較復(fù)雜，同時(shí)由于著陸的工況很多，建模的工作量很大。為進(jìn)行蒙特卡洛仿真分析，開發(fā)了月球著陸器著陸沖擊參數(shù)化非線性有限元建模系統(tǒng)。以 Abaqus/CAE 作為建模平臺(tái)，通過 Abaqus 腳本接口進(jìn)行二次開發(fā)，實(shí)現(xiàn)探測(cè)器參數(shù)化有限元建模。

 

參數(shù)化建模系統(tǒng)的框架見圖 2。參數(shù)化建模系統(tǒng)采用 Excel 2007 作為參數(shù)輸入的圖形用戶界面（GUI），通過 VBA 宏(Macro)生成著陸器著陸工況、著陸面條件及著陸器模型的參數(shù)文件。Python 建模腳本導(dǎo)入建模參數(shù)文件后，通過 Python 解釋器編譯后進(jìn)入到Abaqus/CAE 內(nèi)核中執(zhí)行，進(jìn)入到內(nèi)核中的命令將轉(zhuǎn)換為 INP 文件，再經(jīng)過 Abaqus/Explicit顯式求解器進(jìn)行分析，最后得到輸出數(shù)據(jù)庫(kù)（ODB）文件。Python 建模腳本是參數(shù)化建模系統(tǒng)的核心，它包括有數(shù)據(jù)文件、材料、零件、裝配體、相互作用、連接單元、初始條件、接觸、分析步、邊界條件、載荷、分析作業(yè) 12 個(gè)模塊。

 

 

圖 2 參數(shù)化建模系統(tǒng)框架

Fig.2 Framework of parameterized modeling system

 

3 基于蒙特卡洛法的著陸沖擊載荷分析

本節(jié)將采用蒙特卡洛方法來確定著陸器著陸沖擊設(shè)計(jì)載荷。選用三個(gè)線速度(v x, v y, v z)、三個(gè)角初速度(vYaw, vPitch, vRoll)、三個(gè)姿態(tài)角(yaw, pitch, roll)和月面坡度(slope)共 10個(gè)變量作為分析的輸入?yún)?shù)。坐標(biāo)系定義見圖 1。這些參數(shù)的分布概率密度函數(shù)如表 3 所示。除了偏航角和月面坡度外，其余參數(shù)均服從正態(tài)分布。偏航角(yaw)在 0~45o之間服從均勻分布，0o對(duì)應(yīng) 1-2-1 著陸模式，即一個(gè)著陸腿先觸地，然后 Y 軸上的兩條著陸腿同時(shí)觸地，+Z軸上的著陸腿最后觸地；45o對(duì)應(yīng) 2-2 著陸模式，即兩條腿同時(shí)先觸地，另兩條腿最后同時(shí)著陸。由于采用模擬著陸器模型進(jìn)行仿真分析，著陸器結(jié)構(gòu)是關(guān)于中心軸對(duì)稱的，因此偏航角在 0~45o均勻分布與在 0-360o間均勻分布是等效的。假設(shè)月面坡度服從系數(shù)為 2.113 是 3自由度的卡方(chi-square)分布，即

 

 

 

x 的概率密度函數(shù)為：

 

 

 

式中自由度 n 取 3。

 

表 3 輸入?yún)?shù)概率密度函數(shù)

Tab.3 Probability density function of input parameters

 

 

 

根據(jù)輸入?yún)?shù)的概率密度函數(shù)，隨機(jī)生成 100 個(gè)著陸工況，提交高性能計(jì)算機(jī)進(jìn)行計(jì)算。輸出變量主要包括：著陸腿主支柱與著陸器結(jié)構(gòu)接頭處三個(gè)方向上的分力的峰值載荷 Fx，F(xiàn)y，F(xiàn)z；著陸腿輔助支柱與著陸器結(jié)構(gòu)接頭處三個(gè)分方向的分力的峰值載荷 fx，fy，fz。假設(shè)接頭處三個(gè)方向的載荷峰值服從正態(tài)分布，則根據(jù) 3法則確定著陸器接頭處的設(shè)計(jì)載荷。根據(jù) 3法則，設(shè)計(jì)載荷 F 可通過下式確定：

 

 

式中： 為變量的平均值； 為變量的標(biāo)準(zhǔn)差。

 

采用 3法則的前提是研究的變量服從正態(tài)分布。為了驗(yàn)證 100 個(gè)接頭峰值載荷樣本服從正態(tài)分布，通過 Kolmogorov-Smirnov 檢驗(yàn)來驗(yàn)證輸出的接頭載荷峰值樣本服從正態(tài)分布。Kolmogorov-Smirnov 檢驗(yàn)是檢驗(yàn)單一樣本是否來自某一特定分布的方法，它是以樣本數(shù)據(jù)的累積頻數(shù)分布 F(x)與特定理論分布 F 0(x)比較，若兩者間的差距很小，則推論該樣本取自某特定分布族。即對(duì)于假設(shè)檢驗(yàn)問題：

 

 

 

結(jié)論：當(dāng)實(shí)際觀測(cè) D>D時(shí)，則接受 H1，反之則接受 H0假設(shè)。對(duì)主支柱和輔助支柱三個(gè)方向上的峰值載荷樣本進(jìn)行置信水平為 0.05 的 Kolmogorov- Smirnov 檢驗(yàn)表明，6組樣本均服從正態(tài)分布。主支柱接頭處 X 向峰值載荷載的概率密度函數(shù)和概率分布函數(shù)的經(jīng)驗(yàn)值和理論值見圖 3 和圖 4，峰值載荷樣本的經(jīng)驗(yàn)分布函數(shù)與理論正態(tài)分布函數(shù)的匹配性很好。

 

 

圖 3 峰值載荷頻率直方圖和理論正態(tài)概率密度函數(shù)圖

 

圖 4 峰值載荷經(jīng)驗(yàn)分布函數(shù)和理論正態(tài)分布函數(shù)圖

 

根據(jù) 3法則得到主支柱與輔助支柱接頭處峰值載荷樣本均值、標(biāo)準(zhǔn)差和 3值如表 26所示, 峰值載荷超出 3值的概率僅為 0.13%。

 

表 4 接頭峰值載荷均值、標(biāo)準(zhǔn)差及 3值

 

 

為確定蒙特卡洛法的有效性，對(duì)表 5 中 12 個(gè)典型著陸工況進(jìn)行仿真分析。分析工況的垂直速度和水平速度均為著陸器著陸初始速度的上限。分析得到的主支柱接頭 x 向峰值載荷見表 6，12 個(gè)工況中各條腿的 x 向峰值載荷最大值均小于蒙特卡洛法給出的 3值，證明了蒙特卡洛方法得到 3值是合理的。

 

表 5 典型著陸工況

Tab.5 Typical landing scenarios

 

 

表 6 著陸緩沖裝置主支柱接頭處 x 向峰值載荷

Tab.6 X-direction Maximum Loads at the join of main strut of langding gear

 

 

4 結(jié)束語

在通過試驗(yàn)數(shù)據(jù)驗(yàn)證分析模型正確性的基礎(chǔ)上，利用Abaqus/CAE提供的基于Python腳本的二次開發(fā)接口，設(shè)計(jì)開發(fā)了月球著陸器著陸參數(shù)化建模系統(tǒng)，并利用該系數(shù)對(duì)著陸器著陸過程中的沖擊載荷進(jìn)行蒙特卡洛分析，通過Kolmogorov-Smirnov檢驗(yàn)證明著陸緩沖裝置接頭處的峰值載荷服從正態(tài)分布，并通過3法則確定了最大著陸沖擊載荷。通過對(duì)典型極限著陸工況的仿真分析證明了通過蒙特卡洛方法的有效性。

 

本文提出的基于蒙特卡洛法的著陸器著陸沖擊載荷分析方法可用于確定著陸器著陸過程中使用載荷，從而為著陸器主結(jié)構(gòu)的設(shè)計(jì)提供輸入條件。