馬小兵，任宏道，蔡義坤

(1.北京航空航天大學(xué) 可靠性與系統(tǒng)工程學(xué)院，北京100191;2.北京航空航天大學(xué) 可靠性與環(huán)境工程技術(shù)國(guó)防科技重點(diǎn)實(shí)驗(yàn)室，北京100191)

高超聲速飛行器在大氣環(huán)境中的高速飛行具有嚴(yán)酷的氣動(dòng)加熱問(wèn)題［1］，導(dǎo)致飛行器舵面前緣、進(jìn)氣道等部位產(chǎn)生局部高溫，威脅飛行器的結(jié)構(gòu)安全.因此進(jìn)行高超聲速飛行器高溫結(jié)構(gòu)可靠性分析具有重要意義.目前，高溫結(jié)構(gòu)可靠性研究主要分為對(duì)結(jié)構(gòu)熱響應(yīng)的確定性分析和結(jié)構(gòu)可靠性分析兩部分.其中確定性分析的難點(diǎn)在于氣動(dòng)熱和結(jié)構(gòu)熱響應(yīng)的計(jì)算.目前主要的研究方式是通過(guò)數(shù)值模擬來(lái)進(jìn)行氣動(dòng)加熱相關(guān)的熱-溫度-結(jié)構(gòu)耦合分析［2-3］.在結(jié)構(gòu)可靠性分析方面，傳統(tǒng)的研究思路一般是先建立出結(jié)構(gòu)的極限狀態(tài)方程，然后利用解析法或數(shù)值計(jì)算方法求解結(jié)構(gòu)的可靠度等可靠性指標(biāo).由于復(fù)雜結(jié)構(gòu)的極限狀態(tài)方程一般沒(méi)有解析表達(dá)式或難以得到顯式表達(dá)式，一些研究工作將代理模型技術(shù)應(yīng)用在結(jié)構(gòu)可靠性分析中.代理模型技術(shù)是以擬合精度為約束，利用近似技術(shù)對(duì)離散數(shù)據(jù)進(jìn)行回歸或插值的數(shù)學(xué)模型.目前常用的代理模型有響應(yīng)面、人工神經(jīng)網(wǎng)絡(luò)、徑向基函數(shù)、Kriging插值和支持向量機(jī)等［4-8］.但是這些方法并不能很好地反映結(jié)構(gòu)可靠度的時(shí)變特征［9］，用這些方法計(jì)算得到的可靠度實(shí)際上是結(jié)構(gòu)在特定時(shí)刻的可靠度.當(dāng)結(jié)構(gòu)可靠度隨時(shí)間變化范圍很小時(shí)，采用這些方法具有較好的實(shí)用性.然而對(duì)于高超聲速飛行器的高溫結(jié)構(gòu)問(wèn)題，由于存在快速的升溫過(guò)程，其承載條件、幾何尺寸、材料參數(shù)等都會(huì)隨時(shí)間快速、大幅度地變化，這必將導(dǎo)致結(jié)構(gòu)可靠度在工作過(guò)程中也有較大變化.傳統(tǒng)的結(jié)構(gòu)可靠性模型在解決這類(lèi)時(shí)變問(wèn)題時(shí)需要針對(duì)結(jié)構(gòu)典型時(shí)刻的狀態(tài)分別建模并計(jì)算可靠度，這樣處理只能得到結(jié)構(gòu)在離散時(shí)間點(diǎn)上的可靠度且計(jì)算效率低下［10-11］.

本文通過(guò)對(duì)響應(yīng)面法進(jìn)行擴(kuò)展，提出了一種時(shí)變響應(yīng)面法.在綜合考慮結(jié)構(gòu)熱響應(yīng)及響應(yīng)量閾值時(shí)變效應(yīng)的基礎(chǔ)上建立了結(jié)構(gòu)時(shí)變極限狀態(tài)函數(shù)，給出了結(jié)構(gòu)時(shí)變可靠度計(jì)算方法，并通過(guò)算例驗(yàn)證了該方法的有效性.

1 時(shí)變極限狀態(tài)函數(shù)的建立方法

1.1 時(shí)變響應(yīng)面法

響應(yīng)面法是一種統(tǒng)計(jì)學(xué)綜合試驗(yàn)技術(shù)，用于求解復(fù)雜系統(tǒng)輸出與輸入之間關(guān)系的近似表達(dá)式.傳統(tǒng)的響應(yīng)面法有線(xiàn)性、二階和高階等形式.線(xiàn)性形式的響應(yīng)面法近似能力較差，而高階形式不僅計(jì)算量大，還可能導(dǎo)致在樣本點(diǎn)區(qū)域的外部出現(xiàn)不穩(wěn)定等問(wèn)題［12］，因此一般常用的是帶交叉項(xiàng)的二次響應(yīng)面法，其形式如下:

式中，y表示結(jié)構(gòu)響應(yīng)量;xi(i=1，2，…，n)表示基本輸入變量;αi(i=0，1，…，n)和 αij(i，j=1，2，…，n)均為模型待定系數(shù).

在高速飛行過(guò)程中，高超聲速飛行器舵面前緣、發(fā)動(dòng)機(jī)進(jìn)氣道等局部結(jié)構(gòu)的溫度迅速升高，結(jié)構(gòu)的溫度場(chǎng)、熱應(yīng)力場(chǎng)等響應(yīng)量都隨時(shí)間變化.式(1)只能針對(duì)結(jié)構(gòu)特定時(shí)刻的響應(yīng)量建模，不能反映響應(yīng)量的時(shí)變特征.為了針對(duì)結(jié)構(gòu)在工作過(guò)程中的響應(yīng)量建模，考慮在傳統(tǒng)的響應(yīng)面函數(shù)后面增加與時(shí)間有關(guān)的函數(shù)f(t)，即建立形式如下的響應(yīng)函數(shù):

只要結(jié)構(gòu)在工作過(guò)程中響應(yīng)量隨時(shí)間的變化曲線(xiàn)是近似光滑可導(dǎo)的，就能找到一個(gè)合適的函數(shù)f(t)來(lái)擬合響應(yīng)量隨時(shí)間的變化趨勢(shì).對(duì)于高溫結(jié)構(gòu)來(lái)說(shuō)，不僅其熱響應(yīng)是一個(gè)動(dòng)態(tài)過(guò)程，其熱響應(yīng)對(duì)于輸入?yún)?shù)的靈敏度也是隨時(shí)間變化的，對(duì)式(2)所示的響應(yīng)函數(shù)來(lái)說(shuō):

其中αij=αji.由上式可見(jiàn)響應(yīng)函數(shù)對(duì)基本變量的導(dǎo)數(shù)與時(shí)間無(wú)關(guān)，不能反映結(jié)構(gòu)熱響應(yīng)對(duì)參數(shù)靈敏度隨時(shí)間的動(dòng)態(tài)變化.為了考慮時(shí)間與結(jié)構(gòu)輸入?yún)?shù)的耦合效應(yīng)同時(shí)又控制響應(yīng)面模型的復(fù)雜程度，可建立如下形式的響應(yīng)函數(shù):

據(jù)此:

上式反映了結(jié)構(gòu)響應(yīng)量對(duì)參數(shù)靈敏度的時(shí)變特性.其中f(xi，t)的具體形式需要結(jié)合結(jié)構(gòu)熱物理方程與樣本數(shù)據(jù)規(guī)律來(lái)綜合確定.

1.2 結(jié)構(gòu)響應(yīng)量閾值的時(shí)變特征

高速飛行器在飛行過(guò)程中材料強(qiáng)度、結(jié)構(gòu)屈曲系數(shù)等承載指標(biāo)即響應(yīng)量閾值也是隨時(shí)間變化的.為了得到結(jié)構(gòu)的時(shí)變極限狀態(tài)函數(shù)，還需要建立結(jié)構(gòu)響應(yīng)量閾值的時(shí)變模型.

通過(guò)試驗(yàn)，測(cè)試出結(jié)構(gòu)在不同溫度下的強(qiáng)度等響應(yīng)量閾值參數(shù).利用統(tǒng)計(jì)方法對(duì)樣本數(shù)據(jù)進(jìn)行擬合，得到結(jié)構(gòu)響應(yīng)量閾值與溫度的近似函數(shù):y′=f(T) (6)式中，y′為響應(yīng)量閾值;T為溫度.

由仿真計(jì)算結(jié)果或試驗(yàn)測(cè)試數(shù)據(jù)可得到結(jié)構(gòu)在飛行過(guò)程中的溫度變化數(shù)據(jù)，對(duì)樣本數(shù)據(jù)進(jìn)行擬合，得到結(jié)構(gòu)溫度與時(shí)間的函數(shù):

由式(6)、式(7)可以得到結(jié)構(gòu)響應(yīng)量閾值與時(shí)間的近似函數(shù)，即

1.3 結(jié)構(gòu)時(shí)變可靠度的計(jì)算

結(jié)構(gòu)的時(shí)變極限狀態(tài)函數(shù)等于結(jié)構(gòu)時(shí)變響應(yīng)量閾值函數(shù)與時(shí)變響應(yīng)量函數(shù)之差，由式(4)和式(8)得到:

這里，極限狀態(tài)函數(shù)g(X，t)是一個(gè)隨機(jī)過(guò)程.為了得到結(jié)構(gòu)時(shí)變可靠度等可靠性指標(biāo)，首先計(jì)算其均值μg和標(biāo)準(zhǔn)差σg，其中:

由于結(jié)構(gòu)響應(yīng)量y與響應(yīng)量閾值y′相互獨(dú)立，則

材料的響應(yīng)量閾值(如拉伸強(qiáng)度等)的變異系數(shù)用cv表示，則

當(dāng)f(xi，t)不是線(xiàn)性函數(shù)時(shí)，σy的求解比較困難，可以采用Monte Carlo抽樣的方法求得近似解.一般情況下令f(xi，t)=α′ixit即可以得到較好的擬合精度.時(shí)變響應(yīng)量y的方差可表示為

結(jié)構(gòu)可靠度計(jì)算公式為

由于均值和標(biāo)準(zhǔn)差均是時(shí)間的函數(shù)，因此求解得到的可靠度也是時(shí)間的連續(xù)函數(shù).利用時(shí)變響應(yīng)面法可以計(jì)算得到飛行過(guò)程中熱結(jié)構(gòu)在任意時(shí)刻的可靠度.

2 案例分析

飛行器迎流結(jié)構(gòu)通常是飛行過(guò)程中氣動(dòng)加熱最嚴(yán)酷的部位.下面以某高超聲速飛行器的迎流結(jié)構(gòu)為例，對(duì)上述高溫結(jié)構(gòu)可靠性分析的時(shí)變響應(yīng)面法進(jìn)行說(shuō)明.案例分析只考慮熱載荷，針對(duì)熱強(qiáng)度破壞問(wèn)題進(jìn)行結(jié)構(gòu)時(shí)變可靠性分析，此時(shí)結(jié)構(gòu)響應(yīng)量為熱應(yīng)力，響應(yīng)量閾值為材料的熱強(qiáng)度.

2.1 有限元建模及仿真計(jì)算

根據(jù)結(jié)構(gòu)設(shè)計(jì)參數(shù)及載荷條件建立有限元模型，采用完全耦合算法仿真計(jì)算得到結(jié)構(gòu)的熱溫度場(chǎng)和熱應(yīng)力場(chǎng).結(jié)果表明前緣部位溫度最高、熱應(yīng)力最大，圖1為結(jié)構(gòu)前緣點(diǎn)的應(yīng)力時(shí)間歷程.

圖1 結(jié)構(gòu)前緣點(diǎn)應(yīng)力時(shí)間歷程Fig.1 Stress-time relation of the leading edge structure

2.2 結(jié)構(gòu)可靠性分析

2.2.1 試驗(yàn)設(shè)計(jì)

結(jié)構(gòu)的基本輸入變量包括材料參數(shù)、結(jié)構(gòu)幾何形狀、載荷等.這些變量都存在一定程度的分散性.為了簡(jiǎn)化分析過(guò)程，減少計(jì)算量，通過(guò)多次仿真計(jì)算，篩選出對(duì)該結(jié)構(gòu)熱應(yīng)力影響最大的3個(gè)基本變量:熱流密度系數(shù)k、材料拉伸模量E以及線(xiàn)膨脹系數(shù)α.假設(shè)這3個(gè)基本變量均服從正態(tài)分布，且變異系數(shù)分別為 0.1，0.1 和0.01.

常用的響應(yīng)面試驗(yàn)設(shè)計(jì)方法有中心復(fù)合設(shè)計(jì)(CCDS)和 Box-Behnken 設(shè)計(jì)(BBD)［13-14］.本案例屬于三水平設(shè)計(jì)問(wèn)題，采用BBD設(shè)計(jì)方法.試驗(yàn)設(shè)計(jì)得到13組樣本輸入數(shù)據(jù)，分別進(jìn)行有限元仿真計(jì)算.由于結(jié)構(gòu)前緣中心點(diǎn)處溫度最高、熱應(yīng)力最大，因此選取該單元為研究對(duì)象.圖2是0～60 s內(nèi)各組仿真計(jì)算所得到的熱應(yīng)力時(shí)間歷程.

圖2 各組樣本仿真計(jì)算的應(yīng)力時(shí)間歷程Fig.2 Stress-time relation of each simulation test

2.2.2 建立熱應(yīng)力的時(shí)變響應(yīng)模型

利用13組仿真計(jì)算所得到的樣本數(shù)據(jù)，可以建立時(shí)變響應(yīng)面模型.由于熱應(yīng)力數(shù)據(jù)在0～20 s和20～60 s這兩段時(shí)間內(nèi)的變化趨勢(shì)差異較大，為了簡(jiǎn)化響應(yīng)面形式并提高計(jì)算精度，本文將熱應(yīng)力時(shí)間歷程分為0～20 s和20～60 s兩個(gè)階段，分別建立熱應(yīng)力的時(shí)變響應(yīng)面.通過(guò)分析仿真計(jì)算所得的數(shù)據(jù)特征可知該問(wèn)題中結(jié)構(gòu)響應(yīng)量(熱應(yīng)力)對(duì)輸入?yún)?shù)的靈敏度的變化趨勢(shì)較為簡(jiǎn)單，因此可令式(4)中的 f(xi，t)= αixit，f(t)=b0t+b1t2，從而確定如下形式的響應(yīng)方程:

采用逐步回歸法對(duì)有限元仿真計(jì)算得到的數(shù)據(jù)進(jìn)行回歸擬合，得到結(jié)構(gòu)在0～20 s和20～60 s的熱應(yīng)力函數(shù)如式(17)和式(18)所示:

2.2.3 材料強(qiáng)度的時(shí)變特征分析

采用有限元仿真計(jì)算得到迎流結(jié)構(gòu)前緣點(diǎn)的溫度時(shí)間歷程如圖3所示.

對(duì)超高溫陶瓷材料的強(qiáng)度-溫度數(shù)據(jù)［15］進(jìn)行多項(xiàng)式回歸可得到強(qiáng)度與溫度的近似函數(shù)關(guān)系為

圖3 結(jié)構(gòu)前緣點(diǎn)溫度時(shí)間歷程Fig.3 Temperature-time relation of the leading edge structure

將圖3中的溫度-時(shí)間數(shù)據(jù)代入式(19)中得到材料的強(qiáng)度-時(shí)間數(shù)據(jù)近似函數(shù)如圖4所示.

圖4 材料強(qiáng)度時(shí)變數(shù)據(jù)Fig.4 Time-varying strength data of the material

2.2.4 可靠度計(jì)算結(jié)果

由式(9)、式(17)、式(18)和圖4中的強(qiáng)度-時(shí)間近似函數(shù)可得結(jié)構(gòu)極限狀態(tài)函數(shù)為

利用式(10)～式(14)計(jì)算得到極限狀態(tài)函數(shù)的均值μg和標(biāo)準(zhǔn)差σg，代入式(15)得到結(jié)構(gòu)時(shí)變可靠度的表達(dá)式:熱結(jié)構(gòu)可靠度R隨時(shí)間的變化規(guī)律見(jiàn)圖5.可見(jiàn)，結(jié)構(gòu)在初始階段的可靠度很高.隨著氣動(dòng)加熱的作用，結(jié)構(gòu)溫度上升，材料強(qiáng)度降低，同時(shí)結(jié)構(gòu)的溫度梯度增大導(dǎo)致熱應(yīng)力增大，因此結(jié)構(gòu)可靠度開(kāi)始顯著下降.在大約55 s后，由于傳熱作用，結(jié)構(gòu)內(nèi)的溫度分布趨于均勻，溫度場(chǎng)梯度開(kāi)始減小，熱應(yīng)力也隨之減小，因此結(jié)構(gòu)可靠度有小幅增大.

圖5 熱結(jié)構(gòu)可靠度R與時(shí)間的關(guān)系Fig.5 Time-varying structural reliability of the high-temperature structure

2.3 結(jié)果對(duì)比

為驗(yàn)證該方法的準(zhǔn)確性，將本文計(jì)算結(jié)果與傳統(tǒng)可靠性分析方法的計(jì)算結(jié)果進(jìn)行對(duì)比.由于傳統(tǒng)的結(jié)構(gòu)可靠性分析模型只能計(jì)算結(jié)構(gòu)在特定時(shí)刻的可靠度，因此本文選取 5，10，15，20，25，30，35，40，45，50，55，60 s 這 12 個(gè)時(shí)刻分別采用二次響應(yīng)面法計(jì)算出可靠度，計(jì)算結(jié)果對(duì)比如圖5散點(diǎn)所示.二者的平均相對(duì)誤差為0.05%.

3 結(jié)論

1)對(duì)結(jié)構(gòu)的熱響應(yīng)建立時(shí)變響應(yīng)面模型并通過(guò)逐步回歸確定得到待定參數(shù)的方法具有較高的精度.本文案例中，采用時(shí)變響應(yīng)面模型對(duì)熱應(yīng)力進(jìn)行擬合的相對(duì)誤差小于0.1%.

2)考慮材料強(qiáng)度退化效應(yīng)的高溫結(jié)構(gòu)時(shí)變響應(yīng)面法可高效地計(jì)算出高溫結(jié)構(gòu)可靠度與時(shí)間的關(guān)系，其在特定任務(wù)時(shí)間下的計(jì)算結(jié)果與傳統(tǒng)計(jì)算方法所得結(jié)果一致.由于本文方法不需要針對(duì)特定時(shí)刻重復(fù)建立響應(yīng)面模型，因此在保證計(jì)算精度的同時(shí)可大幅提高計(jì)算效率.

References)

［1］ Moses P L，Rausch V L，Nguyen L T，et al.NASA hypersonic flight demonstrators-overview，status，and future plans［J］.Acta Astronautica，2004，55(3-9):619-630.

［2］ Joshi O，Pénélope L.Stability analysis of a partitioned fluidstructure thermal coupling algorithm［J］.Journal of Thermophysics and Heat Transfer，2014，28(1):59-67.

［3］ 唐健，吳志剛，楊超.考慮結(jié)構(gòu)剛度不確定性的概率顫振分析［J］.北京航空航天大學(xué)學(xué)報(bào)，2014，40(4):569-574.Tang J，Wu Z G，Yang C.Probabilistic flutter analysis with uncertainties in structural stiffness［J］.Journal of Beijing University of Aeronautics and Astronautics，2014，40(4):569-574(in Chinese).

［4］ Bucher C G，Bourgund U.A fast and efficient response surface approach for structural reliability problems［J］.Structural Safety，1990，7(1):57-66.

［5］ 汪保，孫泰.改進(jìn)的Kriging模型的可靠度計(jì)算［J］.計(jì)算機(jī)仿真，2011，28(2):113-116.Wang B，Sun T.Structural reliability computation based on Kriging model［J］.Computer Simulation，2011，28(2):113-116(in Chinese).

［6］ 呂震宙，楊子政，趙潔.基于加權(quán)線(xiàn)性響應(yīng)面法的神經(jīng)網(wǎng)絡(luò)可靠性分析方法［J］.航空學(xué)報(bào)，2006，27(6):1063-1067.Lü Z Z，Yang Z Z，Zhao J.An artificial neural network method for reliability analysis based on weighted linear response surface［J］.Acta Aeronautica et Astronautica Sinica，2006，27(6):1063-1067(in Chinese).

［7］ 趙衛(wèi)，劉濟(jì)科.基于支持向量機(jī)回歸的迭代序列響應(yīng)面可靠度計(jì)算方法［J］.機(jī)械強(qiáng)度，2008，30(6):916-920.Zhao W，Liu J K.Iterative sequence response surface method for reliability computation based on support vector regression［J］.Journal of Mechanical Strength，2008，30(6):916-920(in Chinese).

［8］ 姜潮，黃新萍，韓旭，等.含區(qū)間不確定性的結(jié)構(gòu)時(shí)變可靠度分析方法［J］.機(jī)械工程學(xué)報(bào)，2013，49(10):186-193.Jiang C，Huang X P，Han X，et al.Time-dependent structural reliability analysis method with interval uncertainty［J］.Journal of Mechanical Engineering，2013，49(10):186-193(in Chinese).

［9］ Chakraborty B，Bhar A.Time-varying reliability analysis using a response surface method for a laminated plate under a dynamic load［J］.Journal of Engineering for the Maritime Environment，2014，226(3):260-272.

［10］ Liu J，Li Y.An improved adaptive response surface method for structural reliability analysis［J］.Journal of Central South University，2012，19:1148-1154.

［11］ Rackwitz R.Reliability analysis:a review and some perspectives［J］.Structural Safety，2001，23(4):365-395.

［12］ Gavin H P，Yau S C.High-order limit state functions in the response surface method for structural reliability analysis［J］.Structural Safety，2008，30(2):162-179.

［13］ Ferreira S L C，Bruns R E，F(xiàn)erreira H S，et al.Box-Behnken design:an alternative for the optimization of analytical methods［J］.Analytica Chimica Acta，2007，597(2):179-186.

［14］ 付劍平，陸民燕，阮鐮，等.試驗(yàn)設(shè)計(jì)在軟件可靠性測(cè)試中的應(yīng)用［J］.北京航空航天大學(xué)學(xué)報(bào)，2008，34(12):1379-1383.Fu J P，Lu M Y，Ruan L，et al.Experimental design for software reliability testing［J］.Journal of Beijing University of Aeronautics and Astronautics，2008，34(12):1379-1383(in Chinese).

［15］ Song G M，Zhou Y，Kang S J L.Experimental description of thermomechanical properties of carbon fiber-reinforced TiC matrix composites［J］.Materials ＆ Design，2003，24(8):639-646.