返回艙著水沖擊過載近似模型建模與分析

房紅軍 楊雷

(1 航天東方紅衛(wèi)星有限公司，北京 100094) (2 中國空間技術(shù)研究院載人航天總體部，北京 100094)

載人飛船返回艙可以采取海上回收方式回收，返回艙著水時，艙體及乘員將承受沖擊過載[1]。實踐和分析表明，返回艙著水過程出現(xiàn)的沖擊過載可能對乘員造成傷害[2]。阿波羅12號在與水面碰撞時姿態(tài)偏離5.5°～7.5°，導(dǎo)致15倍重力加速度(15gn)的著水沖擊過載(“阿波羅”所有飛行任務(wù)中最嚴(yán)重的)，艙壁上的一臺攝像機(jī)脫落，將一名航天員頭部砸傷[3-4]。高出允許范圍的沖擊載荷會威脅到航天員的安全，關(guān)系到任務(wù)的成敗，因此在飛船設(shè)計階段需要對該過程的特性有足夠的了解。

國內(nèi)外對返回艙著水沖擊問題進(jìn)行了大量理論和試驗研究。在分析返回艙著水沖擊過載的統(tǒng)計特性、惡劣著水工況及設(shè)計安全著水初始條件的過程中，往往需要獲取大量不同工況條件下的著水沖擊過載。其主要方法有理論公式、物理試驗和數(shù)值計算[5]。目前，解析形式的理論公式難以求解返回艙這種以復(fù)雜外形和姿態(tài)入水的問題；而面對數(shù)千次以上的工況量，也難以通過物理試驗完成。較為常用的手段是采用有限元等數(shù)值計算方法，基于通用瞬態(tài)動力學(xué)分析軟件，求解返回艙著水沖擊這一固、液、氣耦合問題[5]。例如，NASA的星座計劃中，使用LS-DYNA瞬態(tài)動力學(xué)分析程序計算獵戶座(Orion)飛船乘員艙的著水沖擊過載，進(jìn)而結(jié)合蒙特卡羅仿真分析得到“統(tǒng)計上可行”的著水條件，用于飛行任務(wù)設(shè)計[6]。然而，基于有限元等數(shù)值計算方法求解返回艙著水問題，對計算資源的要求非常高。此外，隨著艙體模型復(fù)雜程度的增加，所需計算時間會成倍增長，這極大限制了可分析的工況量。

本文采用克里金(Kriging)方法建立返回艙著水沖擊過載的近似模型，建模過程所需的樣本工況通過拉丁超立方抽樣(LHS)方法獲取，樣本工況下的沖擊過載通過有限元模型計算。該近似模型能夠快速、準(zhǔn)確地計算返回艙著水沖擊過載，可有效提高返回艙著水條件設(shè)計及著水安全性分析的效率。

1 近似模型建模

近似模型(Metamodel)是根據(jù)原問題若干采樣點擬合出來的，是反映原問題的輸入與輸出之間關(guān)系的近似數(shù)學(xué)模型，在某些領(lǐng)域也稱為代理模型(Surrogate Model)或響應(yīng)曲面模型，是一種數(shù)據(jù)擬合(插值或逼近)[7-9]。常用的近似模型有響應(yīng)曲面模型、徑向基函數(shù)模型和Kriging模型。

返回艙著水沖擊過載問題如圖1所示，返回艙以某一速度、姿態(tài)著陸在某一狀態(tài)的海浪上，產(chǎn)生一定的沖擊過載?？紤]到返回艙及海浪真實運動過程的復(fù)雜性，以及著水時間短暫，將問題進(jìn)行一定的簡化：僅考慮返回艙的平動速度，忽略轉(zhuǎn)動速度；海浪采用簡化模型描述。

返回艙著水沖擊過載近似模型建模的任務(wù)是建立返回艙著水狀態(tài)與著水沖擊過載之間的近似關(guān)系，建模流程如圖2所示。近似模型建模的依據(jù)是若干樣本工況以及樣本工況下返回艙的著水沖擊過載，其中，樣本工況是經(jīng)過挑選的若干返回艙著水狀態(tài)(即著水工況)，為提高建模效率，同時保證準(zhǔn)確度，本文通過試驗設(shè)計方法中的LHS方法選擇樣本工況。綜合考慮實現(xiàn)難度、結(jié)果準(zhǔn)確度等因素，采用有限元模型計算出樣本工況下返回艙的著水沖擊過載。根據(jù)樣本工況以及對應(yīng)的沖擊過載，采用Kriging方法即可建立著水沖擊過載的近似模型。

圖1 返回艙著水問題示意Fig.1 Capsule water landing scenario

圖2 返回艙著水沖擊過載近似模型建模過程

1.1 近似模型輸入輸出分析

首先，需要確定近似模型的輸入?yún)?shù)及這些參數(shù)的范圍。對于圖1所示的返回艙著水問題，可通過返回艙的6個絕對運動參數(shù)(3個平動速度和3個姿態(tài)角)和海浪的6個運動參數(shù)描述，即工況條件。但是，設(shè)計變量越多，近似模型的維度就越高，需要的樣本工況數(shù)量會急劇增加。進(jìn)行進(jìn)一步的簡化[10]，將海浪等效為接觸點處的無限大平面，并忽略海浪的水平速度，選擇返回艙相對海浪的5個運動參數(shù)作為近似模型的輸入，即法向速度、切向速度、俯仰角、偏航角和滾動角。5個相對運動參數(shù)的具體值可以根據(jù)運動學(xué)關(guān)系由艙體和海浪的絕對運動參數(shù)轉(zhuǎn)換得到[10]，分析工程中返回艙實際的著水需求和特性，設(shè)計各參數(shù)的范圍(即近似模型的適用范圍)如表1所示。

表1 近似模型輸入?yún)?shù)

工程中關(guān)注的返回艙著水沖擊過載，包括乘員、結(jié)構(gòu)、儀器設(shè)備等位置的數(shù)據(jù)，考慮到本文采用剛性體有限元模型計算樣本工況下返回艙的著水沖擊過載，因此只選取質(zhì)心處加速度的峰值作為近似模型的輸出，其他位置的過載數(shù)據(jù)可根據(jù)運動學(xué)關(guān)系換算得到。沖擊過載的另一個重要特性是增長率，本文也選取質(zhì)心處過載到達(dá)峰值的時間作為近似模型輸出，以反映過載到達(dá)峰值的速度。

1.2 樣本工況抽樣

在確定了近似模型的輸入?yún)?shù)及范圍后，需要通過抽樣方法得到具體的樣本工況，進(jìn)而通過有限元方法計算響應(yīng)輸出。近似模型是根據(jù)樣本點數(shù)據(jù)對原問題其他位置的數(shù)據(jù)進(jìn)行擬合得到的，因此，樣本工況的選擇對近似模型的準(zhǔn)確程度至關(guān)重要，本文采用LHS方法。LHS方法屬于充滿空間的設(shè)計方法，設(shè)計點具有隨機(jī)性，適用于基于計算機(jī)仿真試驗的近似模型構(gòu)造[11]。為獲取更加優(yōu)良的樣本工況，本文采用最大最小距離準(zhǔn)則[12]下的對稱LHS設(shè)計方法，在表1的范圍內(nèi)，抽取300個樣本工況。

1.3 近似模型

F(x)=M(x)+Z(x)=∑Pφ(x)+Z(x)

(1)

式中：M(x)為趨勢項；Z(x)為均值為零的平穩(wěn)隨機(jī)函數(shù)，且該隨機(jī)函數(shù)的協(xié)方差E [Z(w)Z(x)]滿足式(2),只與點的距離有關(guān)，各向同性；P為隨機(jī)變量；φ(x)為一普通基底函數(shù)。

E [Z(w)Z(x)]=σ2R(θ,‖w-x‖)

(2)

式中：w為變量；σ2為F(x)的方差；R(θ,‖w-x‖)為待定的基底函數(shù)，其中，θ為該基底函數(shù)的待定參數(shù)。

Kriging方法就是在隨機(jī)函數(shù)F(x)的線性模型中求解其最小方差線性無偏估計。本文選擇M(x)為隨機(jī)多項式，根據(jù)最小方差要求和無偏要求，采用Lagrange乘子法，可得

f*(x)=m(x)+rTR-1(Y-m)

(3)

式中：未知點x與樣本工況點{si}之間的相關(guān)列陣r=(R(θ,x,si))N×1，N為樣本工況點總數(shù)，i為樣本工況點序號；樣本工況點{si}的協(xié)方差矩陣R=(R(θ,sj,sk))N×N，j，k為樣本工況點序號；樣本工況點{si}的響應(yīng)值Y=(f(si))N×1；m(x)為隨機(jī)過程F(x)的期望函數(shù)，一般未知，見式(4)；樣本工況點{si}上m(x)的函數(shù)值列陣m=(m(si))N×1。

m(x)=E (F(x))=E (∑Pφ(x))=φTβ

(4)

式中：本文選擇的基底函數(shù)列陣φ=(φ(x))B×1=(xξ)B×1，ξ為多項式階數(shù)，B為基底函數(shù)列陣維數(shù)；β為隨機(jī)向量(P)B×1的期望。

同樣，可根據(jù)最小方差及無偏的要求，估計為

β*=(ΦTR-1Φ)-1ΦTR-1Y

(5)

式中：樣本工況點的基底函數(shù)值矩陣Φ=(φ(si))N×B。

結(jié)合式(4)和式(5)，可得到m(x)的估計m*(x)，以及m*=(m*(si))N×1。

根據(jù)式(4)和式(5)，估計隨機(jī)過程的方差為

(6)

本文取隨機(jī)過程F(x)的自相關(guān)函數(shù)為Gauss型，即

(7)

式中：q為變量元素序號。

對于式(7)中的未知參數(shù)θ，采用極大似然估計的方法，引入似然函數(shù)，結(jié)合最優(yōu)化方法求解。本文結(jié)合Kriging建模工具箱DACE[15]進(jìn)行建模，該工具箱使用模式搜索法確定最優(yōu)θ，則最終給出原響應(yīng)函數(shù)的Kriging近似模型為

(8)

2 近似模型誤差與效率分析

作為算例，本文采用Matlab工具箱Model-Based Calibration進(jìn)行LHS，抽取返回艙著水的300個樣本工況和20個檢驗工況；采用LS-DYNA軟件建立返回艙著水有限元模型，計算300個樣本工況下返回艙的著水沖擊過載，對計算結(jié)果進(jìn)行截斷及濾波處理后，采用Matlab工具箱DACE(Kriging近似模型建模工具箱)建立著水沖擊過載峰值及峰值時間的Kriging近似模型。

2.1 近似模型誤差分析

本文采用獨立的檢驗點檢驗近似模型的誤差，以返回艙質(zhì)心處縱向(X軸)過載為例，圖3和圖4分別給出了20個檢驗點處近似模型的預(yù)測結(jié)果與有限元模型計算結(jié)果的對比?？梢姡诟鞣N工況下，近似模型對艙體過載峰值以及到達(dá)過載峰值的時間的擬合程度均較高。

(9)

(10)

(11)

式中：Nc為檢驗點的數(shù)量；t=1，2，…，Nc。

圖3 返回艙質(zhì)心X軸向過載峰值Fig.3 Centroid peak loads along X axis of capsule

圖4 返回艙質(zhì)心X軸過載到達(dá)峰值時間Fig.4 Time of centroid peak loads along X axis of capsule

在20個檢驗工況上，近似模型對過載峰值預(yù)測的偏差統(tǒng)計如表2所示?？梢姡祷嘏撡|(zhì)心過載響應(yīng)近似模型的最大偏差在15 m/s2以下，均方根偏差與平均偏差普遍在7 m/s2以下，而對返回艙安全構(gòu)成威脅的沖擊過載峰值一般在100 m/s2以上[16]，這樣的偏差是在工程可接受的范圍內(nèi)。結(jié)合圖3可判斷，本文建立的返回艙質(zhì)心處著水沖擊過載的近似模型，在全局和局部的準(zhǔn)確程度均比較高。

表2 近似模型的偏差

2.2 近似模型計算效率分析

完成近似建模后，可對沖擊過載進(jìn)行快速預(yù)測。以返回艙著水沖擊問題的蒙特卡羅仿真分析為例，隨機(jī)抽樣5000個工況，分別使用有限元模型和近似模型進(jìn)行計算。表3給出了2種模型計算耗時的對比，其中近似模型的計算耗時考慮了樣本工況和檢驗工況的計算過程耗時。本文采用高性能計算系統(tǒng)進(jìn)行有限元模型計算，單次計算耗時30 min；采用普通計算機(jī)(2.93 GHz處理器×2)進(jìn)行近似模型的計算。

由表3可見，在抽樣量為5000時，結(jié)合少量的有限元計算結(jié)果建立近似模型后進(jìn)行蒙特卡羅仿真，所需計算時間僅是完全采用有限元模型的6.4%，效率提高約16倍；并且，當(dāng)抽樣量增加時，有限元模型的蒙特卡羅仿真計算耗時的增長是線性的，而近似模型的仿真耗時基本不會增長。

3 結(jié)論

本文利用LHS方法抽取了返回艙著水樣本工況，采用有限元模型計算了樣本工況下著水沖擊過載結(jié)果，以樣本工況和工況結(jié)果為輸入，選用Kriging方法建立返回艙著水沖擊過載的近似模型。經(jīng)過與有限元模型進(jìn)行準(zhǔn)確度和計算效率對比，表明這種近似模型可以較為高效、準(zhǔn)確地完成返回艙著水沖擊過載計算。

(1)通過合理地選擇樣本工況，利用少量的有限元計算，采用Kriging方法，能夠建立返回艙著水沖擊過載(峰值與峰值時間)的近似模型，并較為準(zhǔn)確地對返回艙著水沖擊過載峰值進(jìn)行預(yù)測。在本文抽取的20個檢驗工況上，峰值過載的最大偏差不超過15 m/s2，平均偏差不超過7 m/s2。

(2)近似模型能夠?qū)Ψ祷嘏摰闹疀_擊過載進(jìn)行高效計算，對于抽樣量為5000的蒙特卡羅仿真，采用近似模型與有限元模型相比，效率提高約16倍。

(3)近似模型能夠?qū)Ψ祷嘏摰闹疀_擊過載進(jìn)行準(zhǔn)確高效的預(yù)測，為惡劣工況判斷、安全著水初始條件設(shè)計及著水過程的蒙特卡羅仿真提供了一種有效的解決途徑。

參考文獻(xiàn)(References)

[1] 戚發(fā)軔，朱仁璋，李頤黎.載人航天器技術(shù)[M].2版.北京：國防工業(yè)出版社，2003：17-18

Qi Faren, Zhu Renzhang, Li Yili.Manned spacecraft technology [M]. 2nd ed. Beijing: National Defense Industry Press, 2003: 17-18 (in Chinese)

[2] Lawrence Charles, Littell Justin D, Fasanella Edwin L. Orion crew member injury predictions during land and water landings, NASA/TM-2008-215171 [R]. Washington D.C.： NASA, 2007

[3] Johnston R S, Dietlein L F， Berry C A. Biomedical results of Apollo， NASA-SP-368 [R]. Washington D.C.： NASA, 1975

[4] Woods W D. How Apollo flew to the moon [M]. Berlin： Springer, 2007

[5] 房紅軍,楊雷.載人飛船返回艙著水沖擊問題研究進(jìn)展[J].載人航天,2013,19(6):76-83

Fang Hongjun, Yang Lei. Review on water landing impact problems of manned reentry capsules [J]. Manned Spaceflight, 2013, 19(6): 76-83 (in Chinese)

[6] Bryan K Smith, Richard Bouchard, Chung-Chu Teng, et al. Assessment of ocean wave model used to analyze the Constellation Program (CxP) Orion project crew module water landing conditions， NASA/TM-2009-215752 [R]. Washington D.C.： NASA, 2009

[7] 陳余軍,周志成,曲廣吉.多學(xué)科設(shè)計優(yōu)化技術(shù)在衛(wèi)星設(shè)計中的應(yīng)用[J].航天器工程,2013,22(3):16-24

Chen Yujun, Zhou Zhicheng, Qu Guangji. Application of multidisciplinary design optimization in satellite design [J]. Spacecraft Engineering, 2013, 22(3): 16-24 (in Chinese)

[8] 李恩穎.近似模型優(yōu)化體系關(guān)鍵技術(shù)研究及應(yīng)用[D].長沙：湖南大學(xué),2009

Li Enying. The key of research and application of metamodel-based optimization [D]. Changsha: Hunan University, 2009 (in Chinese)

[9] 穆雪峰,姚衛(wèi)星,余雄慶，等.多學(xué)科設(shè)計優(yōu)化中常用代理模型的研究[J].計算力學(xué)學(xué)報,2005,22(5):608-612

Mu Xuefeng, Yao Weixing, Yu Xiongqing, et al. A survey of surrogate models used in MDO [J]. Chinese Journal of Computational Mechanics, 2005, 22(5)： 608-612 (in Chinese)

[10] Arthur M Whitnah, D B Howes. Statistics concerning the Apollo command module water landing, including the probability of occurrence of various impact conditions, successful impact, and body X-axis loads， NASA-TM-X-2430 [R]. Washington D.C.： NASA, 1971

[11] 游海龍,賈新章,張小波，等.試驗設(shè)計與仿真相結(jié)合構(gòu)造集成電路元模型的方法研究[J].電子學(xué)報,2006,36(6):1159-1162

You Hailong, Jia Xinzhang, Zhang Xiaobo, et al. Study on the method of the design of experiment integrated with simulation for constructing the integrated circuit metamodel [J]. Acta Electronica Sinica, 2006, 36(6): 1159-1162 (in Chinese)

[12] Jared L Deutsch, Clayton V Deutsch. Latin hypercube sampling with multi-dimensional uniformity [J]. Journal of Statistical Planning and Inference, 2012, 142: 763-772

[13] Anoop A Mullur, Achille Messac. Metamodeling using extended radial basis functions: a comparative approach [J]. Engineering with Computers, 2006, 21: 203-217

[14] 吳宗敏.散亂數(shù)據(jù)擬合的模型、方法和理論[M].北京：科學(xué)出版社，2007

Wu Zongmin. Model, method and theory of fitting disheveled data [M]. Beijing: Science Press, 2007 (in Chinese)

[15] S?ren N Lophaven, Hans Bruun Nielsen, Jacob S?n￣der￣gaard. DACE-a Matlab Kriging toolbox， Technical Report IMM-TR-2002-12 [R]. Kongens Lyngby： Technical University of Denmark, 2002

[16] NASA. Human integration design handbook (HIDH)， NASA/SP-2010-3407 [S]. Washington D.C.：NASA, 2010