李陽天，李海濱，韋廣梅，翁潔鑫

(1.內(nèi)蒙古財經(jīng)大學(xué) 統(tǒng)計與數(shù)學(xué)學(xué)院，內(nèi)蒙古 呼和浩特 010070；2.內(nèi)蒙古工業(yè)大學(xué) 理學(xué)院，內(nèi)蒙古 呼和浩特 010051；3.內(nèi)蒙古動力機械研究所，內(nèi)蒙古 呼和浩特 010010)

0 引言

不確定性量化是指根據(jù)系統(tǒng)輸入、外部環(huán)境和系統(tǒng)本身的不確定性，在系統(tǒng)內(nèi)部傳播機制下，對系統(tǒng)輸出的不確定性進行量化的過程，它能夠合理地考慮不確定性的影響[1]。其中，多項式混沌展開(PCE)方法應(yīng)用廣泛。何琨等[2]基于稀疏多項式混沌展開(SPCE)方法對微電網(wǎng)概率潮流進行了不確定性分析；Ni等[3]基于PCE方法對橋梁結(jié)構(gòu)的位移響應(yīng)進行了不確定性分析；Peng等[4]基于PCE方法對輸入數(shù)據(jù)不足情況下的復(fù)合層合板進行了不確定性分析；Nguyen等[5]基于PCE方法對可變層化等離子體內(nèi)的電磁傳播特性進行了不確定性分析。然而，映射法求解PCE系數(shù)的過程中，不恰當(dāng)?shù)腜CE階數(shù)與高斯節(jié)點數(shù)易引起過擬合與欠擬合問題。對此，本文提出了改進型PCE方法，以確定自適應(yīng)PCE階數(shù)與自適應(yīng)非均勻網(wǎng)格求解PCE系數(shù)。

固體火箭發(fā)動機在藥柱低溫狀態(tài)下點火，點火后的發(fā)動機受力情況惡劣[6-7]：一方面，受到一定預(yù)應(yīng)變的影響；另一方面，在點火瞬間需承受燃?xì)鈮毫Φ淖饔?。溫度載荷與內(nèi)壓載荷的雙重作用下，材料參數(shù)的隨機性所引起的應(yīng)變響應(yīng)不確定性更容易破壞藥柱的結(jié)構(gòu)完整性。劉中兵等[8]研究了推進劑彈性模量、泊松比和藥柱m數(shù)(藥柱外徑與內(nèi)徑之比)等參數(shù)對藥柱在低溫點火條件下結(jié)構(gòu)完整性的影響；鄧康清等[9]分析了自由填裝式固體火箭發(fā)動機藥柱在低溫點火條件下的結(jié)構(gòu)完整性。但是，以上研究僅關(guān)注到推進劑穩(wěn)定工作所產(chǎn)生的內(nèi)壓載荷與外部環(huán)境溫度所引起的預(yù)應(yīng)變作用，忽略了點火藥被點燃瞬間所形成的沖擊載荷影響。尤其對于小型固體火箭發(fā)動機，由于其初始自由容積較小，點火藥與推進劑燃燒所引起的真實壓力為推進劑穩(wěn)定工作狀態(tài)下工作壓力的數(shù)倍[10]。同時，在固體火箭發(fā)動機生產(chǎn)、加工、貯存與工作過程中，輸入?yún)?shù)的不確定性對藥柱的結(jié)構(gòu)響應(yīng)影響巨大[11]。工程實踐中，藥柱結(jié)構(gòu)響應(yīng)的不確定性研究也一直受到廣泛關(guān)注：Yaman等[12]基于實驗研究了影響推進劑藥柱燃燒速率的重要因素；Griego等[13]以鋁顆粒的燃燒速率為輸出響應(yīng)，研究了輸入?yún)?shù)對輸出響應(yīng)的靈敏度，并進行了不確定性量化分析。

因此，基于所提改進型PCE方法分析材料參數(shù)隨機性對藥柱輸出響應(yīng)的影響，獲得輸出響應(yīng)的概率分布，對固體火箭發(fā)動機不確定性分析計算精度的提高及其在工程中的應(yīng)用，具有積極的參考價值。本文主要介紹了PCE方法和PCE后處理的基本原理；基于自適應(yīng)PCE階數(shù)、自適應(yīng)均勻網(wǎng)格和自適應(yīng)非均勻網(wǎng)格概念，給出了改進型PCE方法的基本步驟；以固體火箭發(fā)動機低溫點火過程中的不確定性分析為例，基于所提方法分析材料參數(shù)隨機性對藥柱輸出響應(yīng)的影響，獲得輸出響應(yīng)的概率分布，驗證了所提方法相對傳統(tǒng)方法的優(yōu)越性。

1 不確定性量化方法

1.1 PCE方法

PCE方法是一種將模型輸出映射在隨機輸入的正交隨機多項式基函數(shù)之上的概率統(tǒng)計方法。隨機變量u(θ)的PCE可以表示為

(1)

(2)

式中：ξ為由P個高斯隨機變量(ξi1,…,ξiP)組成的隨機向量。(1)式也可表示為

(3)

式中：bi和Ψi[ξ(θ)]分別與ai1,…,aiP和ΓP(ξi1,…,ξiP)一一對應(yīng)?；煦缍囗検降恼恍钥梢员磉_為

(4)

式中：δij表示克羅內(nèi)克符號函數(shù)。根據(jù)計算需要，PCE要求截斷為

(5)

式中：NP為截斷展開的項數(shù)，其值取決于隨機變量維數(shù)d和多項式展開階數(shù)P，NP可以表示為

(6)

映射法求解PCE系數(shù)：基于混沌多項式的正交性，映射法可以作為一種非侵入式方法被用于計算響應(yīng)的展開系數(shù)[14-15]。(3)式兩邊同時乘以Ψj(ξ)，并取期望值，得

(7)

由于多項式基函數(shù)的正交性，(7)式可以表示為

(8)

式中：分母可以以一種解析解的形式查表得到[16-17]，然而分子需要運用蒙特卡洛法或高斯積分獲得。高斯積分所用高斯節(jié)點值與積分權(quán)值通過文獻[18]可查，6組典型高斯厄米特積分節(jié)點的節(jié)點值和積分權(quán)值如表1所示。其中，N為高斯節(jié)點數(shù)，ξq和wq(q=1,2,…，N)分別為一維高斯厄米特積分節(jié)點的節(jié)點值與積分權(quán)值。(8)式可以表示為

表1 高斯- 厄米特積分典型權(quán)值和積分點

(9)

式中：nm表示第m維變量上的第n個高斯節(jié)點；ξnm和wnm分別表示第m維變量上第n個高斯節(jié)點(即第nm個高斯節(jié)點)的節(jié)點值與積分權(quán)值；xnm為ξnm在原隨機空間(x空間)中的變換值。

1.2 PCE后處理

1.2.1 統(tǒng)計特性計算

根據(jù)厄米特多項式的展開特性，可以基于多項式展開系數(shù)直接計算輸出響應(yīng)量前4階統(tǒng)計矩的解析解。輸出響應(yīng)量u的均值μu、標(biāo)準(zhǔn)差σu、偏度δu、和峰度κu的解析表達式[19-20]分別為

(10)

1.2.2 全局靈敏度計算

根據(jù)厄米特多項式的展開特性，可以基于多項式展開系數(shù)直接計算各輸入變量的全局靈敏度指數(shù)。

輸入變量ξi的1階靈敏度指數(shù)[21]為

(11)

式中：分子為(1)式中所有只涉及單變量ξi項的方差之和；分母為輸出變量uPC(ξ)的總方差；多變量指數(shù)集li可表示為

li1,i2,…,is=

(12)

式中：li1,i2,…,is為從PCE表達式中選擇只含有s個變量ξi1，…，ξis的展開項所得的指數(shù)集；s為指數(shù)集所選變量個數(shù)。

輸入變量ξi的總靈敏度指數(shù)可表示為

(13)

式中：指數(shù)集Ii通過選取PCE中所有出現(xiàn)變量ξi的項所構(gòu)成[16]，

Ii1,i2,…,is={α:αk>0 ?k=1,…,d，

k∈(i1,i2,…,is)}.

(14)

1.3 改進型PCE方法

對于PCE代理模型映射求解過程中的過擬合與欠擬合問題，Suryawanshi等[22]提出基于參數(shù)的相對重要性建立非均勻網(wǎng)格，利用映射法求解PCE系數(shù)可以減少PCE系數(shù)求解的計算量和提高計算精度，但并未說明自適應(yīng)PCE階數(shù)和自適應(yīng)非均勻網(wǎng)格的確定方法。本文提出了改進型PCE方法，即基于全局靈敏度指數(shù)的非均勻網(wǎng)格方法，映射求解PCE系數(shù)，所提非均勻網(wǎng)格方法以PCE系數(shù)映射求解中的均勻網(wǎng)格法為基礎(chǔ)，流程圖如圖1所示。圖1中P*為自適應(yīng)PCE階數(shù)，N*為自適應(yīng)高斯節(jié)點數(shù)，R為1個常數(shù)，其初始值R0=0，ΔR為R的增量。

1.3.1 均勻網(wǎng)格法

均勻網(wǎng)格法基于全因子數(shù)值積分技術(shù)(FFNI)[23]，其具體步驟為：

1)建立均勻網(wǎng)格，即d維高斯節(jié)點的確定。將所得1維節(jié)點與權(quán)值進行直接張量積操作，得到d維節(jié)點與權(quán)值，其中，各維節(jié)點數(shù)相等N1=N2=…=Nd=K.

2)考慮非正態(tài)變量的相關(guān)性，根據(jù)Nataf變換的逆變換方法[24]，將每個高斯節(jié)點從標(biāo)準(zhǔn)隨機空間(ξ空間)變換到原隨機空間(x空間)：x1m=F1(ξ1m),…,xnm=Fn(ξnm),…,xkm=FK(ξKm)，其中，F(xiàn)n(ξnm)表示將第m維變量中的第n個高斯節(jié)點ξnm，從標(biāo)準(zhǔn)隨機空間變換到原隨機空間。

3)將樣本點x=[x1,…,xm,…,xd]代入實際響應(yīng)模型，得到各高斯節(jié)點處的真實響應(yīng)u，其中，xm=[x1m,…,xnm,…,xKm]T,u=[u1,…,un,…,uK]T。

4)將ξ=[ξ1,…,ξn,…,ξK]T、w=[w1,…,wn,…,wK]T和u=[u1,…,un,…,uK]T代入(9)式，得到PCE系數(shù)bj(j=1,2,…,NP)，其中，ξn=[ξn1,…,ξnm,…,ξnd]，wn=[wn1,…,wnm,…,wnd]。

1.3.2 改進型PCE方法

如圖1所示，改進型PCE方法計算過程為：

1)自適應(yīng) PCE階數(shù)P*的確定。

① 設(shè)初始PCE階數(shù)P=2，各維隨機變量高斯節(jié)點數(shù)K=3，利用均勻網(wǎng)格法求解PCE系數(shù)。

② 基于(10)式計算已得PCE模型方差σ2，將結(jié)果代入(15)式[25]：

(15)

式中：根據(jù)文獻[25]，設(shè)常數(shù)e0=5.

如果E0滿足(15)式，則自適應(yīng)PCE階數(shù)P*=P；否則，設(shè)P=P+1，返回步驟①，開始新運算，直至(15)式得到滿足。

2)自適應(yīng)均勻網(wǎng)格的確定。

①設(shè)自適應(yīng)PCE階數(shù)P*=P，各維隨機變量初始高斯節(jié)點數(shù)N=3，利用均勻網(wǎng)格法求解PCE系數(shù)。

②基于(10)式計算已得PCE模型均值μ與標(biāo)準(zhǔn)差σ，將結(jié)果代入(16)式[25]，得E1、E2.如果E1、E2滿足(16)式，則自適應(yīng)高斯節(jié)點數(shù)N*=N；否則，設(shè)N=N+1，返回步驟①，開始新運算，直至(16)式得到滿足。

(16)

式中：根據(jù)文獻[25]，e1、e2為區(qū)間(0，0.05]的常數(shù)。

3)自適應(yīng)非均勻網(wǎng)格的確定。

①根據(jù)所得自適應(yīng)PCE階數(shù)P*和自適應(yīng)均勻網(wǎng)格(N*)d，利用均勻網(wǎng)格法求解PCE系數(shù)?；?13)式，得到各變量的總靈敏度指數(shù)ST1,…,STi,…,STd.

②基于(17)式，計算變量ξ1,…,ξi,…,ξd上的高斯節(jié)點數(shù)N1,…,Ni,…,Nd，進而得到一個d維網(wǎng)格。如果所得網(wǎng)格滿足(18)式，則保留此網(wǎng)格；否則，設(shè)R=R+ΔR，重新基于(17)式，求解新網(wǎng)格，直至(18)式得到滿足。

Ni=3+[STi×R]，

(17)

N1(Rj)×…×Nd(Rj)>N1(Rj-1)×…×Nd(Rj-1)，

(18)

(19)

式中：[·]表示取整運算；min{STi}表示集合{ST1,…,STi,…,STd}中的最小值。

③基于所得網(wǎng)格，求解PCE系數(shù)。

④基于(10)式計算已得PCE模型均值μ與標(biāo)準(zhǔn)差σ，將結(jié)果代入(20)式[25]，得E3、E4.如果E3、E4滿足(20)式，則此網(wǎng)格為自適應(yīng)非均勻網(wǎng)格；否則，設(shè)R=R+ΔR，返回步驟②，開始新運算，直至(20)式得到滿足。

(20)

式中：根據(jù)文獻[25]，e3、e4為區(qū)間(0，0.03]的常數(shù)。

⑤基于自適應(yīng)PCE階數(shù)與自適應(yīng)非均勻網(wǎng)格，求解PCE系數(shù)和輸出響應(yīng)統(tǒng)計矩。

2 固體火箭發(fā)動機低溫點火不確定性量化

2.1 問題描述

某固體火箭發(fā)動機由殼體、絕緣層和六角星形藥柱構(gòu)成，且藥柱材料為異佛爾酮二異氰酸酯/端羥基聚丁二烯(IPDI-HTPB)丁羥基推進劑。考慮到發(fā)動機幾何構(gòu)型和載荷的對稱性，建立藥柱結(jié)構(gòu)的1/12周期對稱性模型，如圖2所示。幾何參數(shù)：藥柱長度L=1 600×10-3m、殼體厚度Rc=4×10-3m、絕緣層厚度Ri=2.7×10-3m、藥柱外半徑r=250×10-3m、藥柱內(nèi)半徑ri=32.5×10-3m、肉厚D=130×10-3m、星角數(shù)S=6、星槽圓弧半徑Pg=5×10-3m、星邊夾角B=47°、星角系數(shù)ε=0.8、星跟圓弧半徑Pb=11×10-3m.

圖2 藥柱橫截面

以該固體發(fā)動機低溫點火為例：點火前，藥柱經(jīng)4 d時間由零應(yīng)力溫度(58 ℃)降至-40 ℃，藥柱承受-40 ℃溫度載荷單獨作用；藥柱在-40 ℃自身溫度條件下點火后，藥柱承受溫度載荷預(yù)應(yīng)變與內(nèi)壓載荷的聯(lián)合疊加作用。設(shè)點火后的0.15 s時刻，藥柱燃燒室內(nèi)壓(點火藥與推進劑燃燒所引起)達到初始壓力峰值且壓力峰值為44 MPa(假設(shè)點火藥與推進劑燃燒充分)[10]；點火后藥柱內(nèi)表面溫度為2 500 ℃.

工程實踐表明，藥柱低溫點火過程中，藥柱內(nèi)孔表面是最危險部位[26-27]，此時以藥柱內(nèi)表面伸長率為失效判據(jù)[28-29]。結(jié)合von Mises等效應(yīng)變理論，取0.15 s時刻藥柱內(nèi)表面最大等效應(yīng)變值作為研究對象，研究藥柱材料參數(shù)的隨機性對輸出響應(yīng)的影響。建立三維黏彈性有限元分析模型如圖3所示。

圖3 藥柱三維有限元模型

2.2 有限元建模

藥柱低溫點火過程中，受溫度載荷與內(nèi)壓載荷雙重疊加作用。文獻[30]指出：在溫度載荷作用下，主要影響結(jié)構(gòu)完整性的是藥柱的泊松比PRXY與線膨脹系數(shù)ALPX；在內(nèi)壓載荷作用下，主要影響結(jié)構(gòu)完整性的是藥柱泊松比及推進劑初始模量，但只有較小的推進劑初始模量對應(yīng)變響應(yīng)作用明顯。藥柱低溫點火過程中的初始模量不屬于較小初始模量范疇，故本文選擇推進劑線膨脹系數(shù)與泊松比作為研究的不確定性因素，其分布信息如表2所示。

表2 不確定性參數(shù)概率分布信息

當(dāng)參考溫度為293.27 K時，表示時間- 溫度移位因子的WLF方程可以表示為

(21)

式中：T為溫度。

經(jīng)Prony級數(shù)擬合，藥柱的松弛模量可以表示為

(22)

式中：E∞為平衡模量；NM為Maxwell模型單元數(shù)；Ei為第i個Maxwell單元的彈性模量；τi為第i個Maxwell單元的松弛時間。

E(t)的各個系數(shù)如表3所示。

表3 IPDI-HTPB推進劑的Prony級數(shù)

2.3 基于改進型PCE方法計算過程

2.3.1 自適應(yīng)PCE階數(shù)的確定

設(shè)均勻網(wǎng)格方法中各變量的高斯節(jié)點數(shù)為K=3，基于均勻網(wǎng)格方法計算PCE階數(shù)P=1,…,5時的PCE模型，求解其方差，計算結(jié)果如圖4所示。當(dāng)P=3時，E0=30.474 5>e0=5；由(15)式或圖4可判斷，P=3為自適應(yīng)PCE階數(shù)，即自適應(yīng)PCE階數(shù)P*=3.該過程的實際計算成本為3×3=9次真實響應(yīng)模型運算。圖4中PCE(3×3)表示基于3×3均勻網(wǎng)格PCE代理模型方法，PCE(4×4)表示基于4×4均勻網(wǎng)格PCE代理模型方法，MCS(1×104)表示基于1×104個樣本的蒙特卡洛仿真(MCS)方法。

圖4 自適應(yīng)PCE階數(shù)P* 的確定

2.3.2 自適應(yīng)均勻網(wǎng)格的確定

設(shè)自適應(yīng)PCE階數(shù)P*=3，基于均勻網(wǎng)格方法計算各變量高斯節(jié)點數(shù)分別為N=3,4,…時的PCE模型，求解其模型期望與標(biāo)準(zhǔn)差，計算結(jié)果如圖5 所示。N=4時，E1=0.002 4<0.05，E2=0.033 7<0.05；由(16)式或圖5可判斷，N*=4為均勻網(wǎng)格方法中自適應(yīng)高斯節(jié)點數(shù)，即4×4網(wǎng)格為所求自適應(yīng)均勻網(wǎng)格。由于3×3均勻網(wǎng)格的重復(fù)使用，該過程的實際計算成本為4×4=16次真實響應(yīng)模型運算。圖5中PCE(P=3)表示3階PCE代理模型方法。

圖5 均勻網(wǎng)格自適應(yīng)高斯節(jié)點數(shù)N*的確定

2.3.3 全局靈敏度分析

基于自適應(yīng)PCE階數(shù)P*=3和自適應(yīng)均勻網(wǎng)格(N*)2=4×4的PCE模型對藥柱的線膨脹系數(shù)與泊松比進行全局靈敏度分析，計算結(jié)果如圖6 所示。

圖6 輸入變量靈敏度指數(shù)

2.3.4 自適應(yīng)非均勻網(wǎng)格確定

基于1.3.2節(jié)中的步驟3，建立非均勻網(wǎng)格。設(shè)自適應(yīng)PCE階數(shù)P*=3，基于所得非均勻網(wǎng)格計算PCE模型，求解其模型期望與標(biāo)準(zhǔn)差，計算結(jié)果如圖7 所示。圖7中：R的取值分別為0、1.041 7、2.083 3、3.125 0、4.166 7；相應(yīng)非均勻網(wǎng)格分別為3×3、4×4、4×5、5×5、6×6.當(dāng)R=2.083 3時，E3=0.001 8<0.03，E4=0.022 4<0.03；由(20)式或圖7可判斷，當(dāng)R=2.083 3時可得自適應(yīng)非均勻網(wǎng)格，即4×5網(wǎng)格為自適應(yīng)非均勻網(wǎng)格。由于3×3和4×4網(wǎng)格的重復(fù)使用，該過程的實際計算成本為4×5=20次真實響應(yīng)模型運算。因此，本文所提方法的總計算成本為9+16+20=45次真實響應(yīng)模型運算。

圖7 自適應(yīng)非均勻網(wǎng)格確定

2.4 結(jié)果分析

基于所得自適應(yīng)PCE階數(shù)P*=3和自適應(yīng)(4×5)非均勻網(wǎng)格(即輸入變量ALPX上取4個高斯節(jié)點，輸入變量PRXY上取5個高斯節(jié)點，構(gòu)成“4×5”的非均勻網(wǎng)格，映射法求解PCE系數(shù))，提出了改進型PCE方法，計算PCE模型，求解其模型期望與標(biāo)準(zhǔn)差。所得結(jié)果與4種傳統(tǒng)方法作比較，其中，包括2階隨機響應(yīng)面(SRSM)方法、3階SRSM方法、采用7×7均勻網(wǎng)格求解3階PCE系數(shù)的FFNI方法和10 000次MCS方法。其中，以MCS方法所得結(jié)果為精確解。

所得結(jié)果的概率密度函數(shù)(PDF)與累積分布函數(shù)(CDF)如圖8、圖9所示。在表4中，由于FFNI、SRSM和所提方法計算成本相近，在精度方面對所得結(jié)果進行對比。結(jié)果顯示：本文所提方法的均值誤差與3階SRSM方法相近且遠(yuǎn)小于其他方法，表明藥柱內(nèi)表面最大等效應(yīng)變響應(yīng)的平均值誤差減?。槐疚乃岱椒ǖ臉?biāo)準(zhǔn)差誤差遠(yuǎn)小于其他方法，表明藥柱內(nèi)表面最大等效應(yīng)變響應(yīng)的離散性誤差減?。槐疚乃岱椒ǖ钠日`差與FFNI方法相近且遠(yuǎn)小于其他方法，表明藥柱內(nèi)表面最大等效應(yīng)變響應(yīng)的非對稱程度誤差減??；本文所提方法的峰度誤差與FFNI方法相近且遠(yuǎn)小于其他方法，表明藥柱內(nèi)表面最大等效應(yīng)變響應(yīng)的尖銳程度誤差減小。綜上所述，本文所提方法的不確定性分析結(jié)果優(yōu)于其他方法。

圖8 不同方法所得藥柱內(nèi)表面應(yīng)變響應(yīng)的PDF

圖9 不同方法所得藥柱內(nèi)表面應(yīng)變響應(yīng)的CDF

表4 藥柱內(nèi)表面的最大等效應(yīng)變的前4階統(tǒng)計矩

3 結(jié)論

傳統(tǒng)PCE求解方法，包括回歸法和映射法(如FFNI方法)在求解PCE系數(shù)的過程中，對各輸入變量等量采樣，并未考慮各輸入變量對輸出響應(yīng)結(jié)果影響程度的不同。不恰當(dāng)?shù)腜CE階數(shù)與采樣點數(shù)會在PCE模型求解的過程中，引起過擬合與欠擬合問題。為提高PCE模型求解精度，本文提出了改進型PCE方法。該方法改變了傳統(tǒng)映射法與回歸法將輸入變量同等對待的習(xí)慣，求解基于各輸入變量全局靈敏度指數(shù)的自適應(yīng)非均勻網(wǎng)格。進而基于此非均勻網(wǎng)格，采用映射法求解PCE系數(shù)。

本文基于固體火箭發(fā)動機低溫點火過程中的不確定性材料參數(shù)(藥柱的線膨脹系數(shù)與泊松比)對于藥柱內(nèi)表面最大應(yīng)變響應(yīng)靈敏度結(jié)果的不同，建立了4×5的自適應(yīng)非均勻網(wǎng)格，以此求解PCE系數(shù)。結(jié)果表明，本文所提方法較傳統(tǒng)方法，能夠在服從效率要求的條件下，明顯提高了計算精度，證明其具有良好的工程應(yīng)用價值。但本文所提方法還只是在輸入變量較少的低維問題中證明其有效，對于輸入變量較多的高維問題，還有待進一步深入研究。