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

(1.內蒙古財經大學 統計與數學學院，內蒙古 呼和浩特 010070；2.內蒙古工業(yè)大學 理學院，內蒙古 呼和浩特 010051；3.內蒙古動力機械研究所，內蒙古 呼和浩特 010010)

0 引言

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

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

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

1 不確定性量化方法

1.1 PCE方法

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

(1)

(2)

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

(3)

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

(4)

式中：δij表示克羅內克符號函數。根據計算需要，PCE要求截斷為

(5)

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

(6)

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

(7)

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

(8)

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

表1 高斯- 厄米特積分典型權值和積分點

(9)

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

1.2 PCE后處理

1.2.1 統計特性計算

根據厄米特多項式的展開特性，可以基于多項式展開系數直接計算輸出響應量前4階統計矩的解析解。輸出響應量u的均值μu、標準差σu、偏度δu、和峰度κu的解析表達式[19-20]分別為

(10)

1.2.2 全局靈敏度計算

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

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

(11)

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

li1,i2,…,is=

(12)

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

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

(13)

式中：指數集Ii通過選取PCE中所有出現變量ξi的項所構成[16]，

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

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

(14)

1.3 改進型PCE方法

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

1.3.1 均勻網格法

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

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

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

3)將樣本點x=[x1,…,xm,…,xd]代入實際響應模型，得到各高斯節(jié)點處的真實響應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系數bj(j=1,2,…,NP)，其中，ξn=[ξn1,…,ξnm,…,ξnd]，wn=[wn1,…,wnm,…,wnd]。

1.3.2 改進型PCE方法

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

1)自適應 PCE階數P*的確定。

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

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

(15)

式中：根據文獻[25]，設常數e0=5.

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

2)自適應均勻網格的確定。

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

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

(16)

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

3)自適應非均勻網格的確定。

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

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

Ni=3+[STi×R]，

(17)

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

(18)

(19)

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

③基于所得網格，求解PCE系數。

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

(20)

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

⑤基于自適應PCE階數與自適應非均勻網格，求解PCE系數和輸出響應統計矩。

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

2.1 問題描述

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

圖2 藥柱橫截面

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

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

圖3 藥柱三維有限元模型

2.2 有限元建模

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

表2 不確定性參數概率分布信息

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

(21)

式中：T為溫度。

經Prony級數擬合，藥柱的松弛模量可以表示為

(22)

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

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

表3 IPDI-HTPB推進劑的Prony級數

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

2.3.1 自適應PCE階數的確定

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

圖4 自適應PCE階數P* 的確定

2.3.2 自適應均勻網格的確定

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

圖5 均勻網格自適應高斯節(jié)點數N*的確定

2.3.3 全局靈敏度分析

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

圖6 輸入變量靈敏度指數

2.3.4 自適應非均勻網格確定

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

圖7 自適應非均勻網格確定

2.4 結果分析

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

所得結果的概率密度函數(PDF)與累積分布函數(CDF)如圖8、圖9所示。在表4中，由于FFNI、SRSM和所提方法計算成本相近，在精度方面對所得結果進行對比。結果顯示：本文所提方法的均值誤差與3階SRSM方法相近且遠小于其他方法，表明藥柱內表面最大等效應變響應的平均值誤差減??；本文所提方法的標準差誤差遠小于其他方法，表明藥柱內表面最大等效應變響應的離散性誤差減小；本文所提方法的偏度誤差與FFNI方法相近且遠小于其他方法，表明藥柱內表面最大等效應變響應的非對稱程度誤差減??；本文所提方法的峰度誤差與FFNI方法相近且遠小于其他方法，表明藥柱內表面最大等效應變響應的尖銳程度誤差減小。綜上所述，本文所提方法的不確定性分析結果優(yōu)于其他方法。

圖8 不同方法所得藥柱內表面應變響應的PDF

圖9 不同方法所得藥柱內表面應變響應的CDF

表4 藥柱內表面的最大等效應變的前4階統計矩

3 結論

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

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