房紅軍 楊雷

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

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

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

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

1 近似模型建模

近似模型(Metamodel)是根據(jù)原問題若干采樣點(diǎn)擬合出來的，是反映原問題的輸入與輸出之間關(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)生一定的沖擊過載。考慮到返回艙及海浪真實(shí)運(yùn)動(dòng)過程的復(fù)雜性，以及著水時(shí)間短暫，將問題進(jìn)行一定的簡(jiǎn)化：僅考慮返回艙的平動(dòng)速度，忽略轉(zhuǎn)動(dòng)速度；海浪采用簡(jiǎn)化模型描述。

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

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

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

1.1 近似模型輸入輸出分析

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

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

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

1.2 樣本工況抽樣

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

1.3 近似模型

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

(1)

式中：M(x)為趨勢(shì)項(xiàng)；Z(x)為均值為零的平穩(wěn)隨機(jī)函數(shù)，且該隨機(jī)函數(shù)的協(xié)方差E [Z(w)Z(x)]滿足式(2),只與點(diǎn)的距離有關(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)的線性模型中求解其最小方差線性無偏估計(jì)。本文選擇M(x)為隨機(jī)多項(xiàng)式，根據(jù)最小方差要求和無偏要求，采用Lagrange乘子法，可得

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

(3)

式中：未知點(diǎn)x與樣本工況點(diǎn){si}之間的相關(guān)列陣r=(R(θ,x,si))N×1，N為樣本工況點(diǎn)總數(shù)，i為樣本工況點(diǎn)序號(hào)；樣本工況點(diǎn){si}的協(xié)方差矩陣R=(R(θ,sj,sk))N×N，j，k為樣本工況點(diǎn)序號(hào)；樣本工況點(diǎn){si}的響應(yīng)值Y=(f(si))N×1；m(x)為隨機(jī)過程F(x)的期望函數(shù)，一般未知，見式(4)；樣本工況點(diǎn){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，ξ為多項(xiàng)式階數(shù)，B為基底函數(shù)列陣維數(shù)；β為隨機(jī)向量(P)B×1的期望。

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

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

(5)

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

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

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

(6)

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

(7)

式中：q為變量元素序號(hào)。

對(duì)于式(7)中的未知參數(shù)θ，采用極大似然估計(jì)的方法，引入似然函數(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個(gè)樣本工況和20個(gè)檢驗(yàn)工況；采用LS-DYNA軟件建立返回艙著水有限元模型，計(jì)算300個(gè)樣本工況下返回艙的著水沖擊過載，對(duì)計(jì)算結(jié)果進(jìn)行截?cái)嗉盀V波處理后，采用Matlab工具箱DACE(Kriging近似模型建模工具箱)建立著水沖擊過載峰值及峰值時(shí)間的Kriging近似模型。

2.1 近似模型誤差分析

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

(9)

(10)

(11)

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

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

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

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

表2 近似模型的偏差

2.2 近似模型計(jì)算效率分析

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

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

3 結(jié)論

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

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

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

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

參考文獻(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è)計(jì)優(yōu)化技術(shù)在衛(wèi)星設(shè)計(jì)中的應(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è)計(jì)優(yōu)化中常用代理模型的研究[J].計(jì)算力學(xué)學(xué)報(bào),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] 游海龍,賈新章,張小波，等.試驗(yàn)設(shè)計(jì)與仿真相結(jié)合構(gòu)造集成電路元模型的方法研究[J].電子學(xué)報(bào),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