基于改进型多项式混沌展开的固体火箭发动机药柱低温点火不确定性量化分析

李阳天，李海滨，韦广梅，翁洁鑫

(1.内蒙古财经大学 统计与数学学院，内蒙古 呼和浩特 010070；2.内蒙古工业大学 理学院，内蒙古 呼和浩特 010051；3.内蒙古动力机械研究所，内蒙古 呼和浩特 010010)

0 引言

不确定性量化是指根据系统输入、外部环境和系统本身的不确定性，在系统内部传播机制下，对系统输出的不确定性进行量化的过程，它能够合理地考虑不确定性的影响[1]。其中，多项式混沌展开(PCE)方法应用广泛。何琨等[2]基于稀疏多项式混沌展开(SPCE)方法对微电网概率潮流进行了不确定性分析；Ni等[3]基于PCE方法对桥梁结构的位移响应进行了不确定性分析；Peng等[4]基于PCE方法对输入数据不足情况下的复合层合板进行了不确定性分析；Nguyen等[5]基于PCE方法对可变层化等离子体内的电磁传播特性进行了不确定性分析。然而，映射法求解PCE系数的过程中，不恰当的PCE阶数与高斯节点数易引起过拟合与欠拟合问题。对此，本文提出了改进型PCE方法，以确定自适应PCE阶数与自适应非均匀网格求解PCE系数。

固体火箭发动机在药柱低温状态下点火，点火后的发动机受力情况恶劣[6-7]：一方面，受到一定预应变的影响；另一方面，在点火瞬间需承受燃气压力的作用。温度载荷与内压载荷的双重作用下，材料参数的随机性所引起的应变响应不确定性更容易破坏药柱的结构完整性。刘中兵等[8]研究了推进剂弹性模量、泊松比和药柱m数(药柱外径与内径之比)等参数对药柱在低温点火条件下结构完整性的影响；邓康清等[9]分析了自由填装式固体火箭发动机药柱在低温点火条件下的结构完整性。但是，以上研究仅关注到推进剂稳定工作所产生的内压载荷与外部环境温度所引起的预应变作用，忽略了点火药被点燃瞬间所形成的冲击载荷影响。尤其对于小型固体火箭发动机，由于其初始自由容积较小，点火药与推进剂燃烧所引起的真实压力为推进剂稳定工作状态下工作压力的数倍[10]。同时，在固体火箭发动机生产、加工、贮存与工作过程中，输入参数的不确定性对药柱的结构响应影响巨大[11]。工程实践中，药柱结构响应的不确定性研究也一直受到广泛关注：Yaman等[12]基于实验研究了影响推进剂药柱燃烧速率的重要因素；Griego等[13]以铝颗粒的燃烧速率为输出响应，研究了输入参数对输出响应的灵敏度，并进行了不确定性量化分析。

因此，基于所提改进型PCE方法分析材料参数随机性对药柱输出响应的影响，获得输出响应的概率分布，对固体火箭发动机不确定性分析计算精度的提高及其在工程中的应用，具有积极的参考价值。本文主要介绍了PCE方法和PCE后处理的基本原理；基于自适应PCE阶数、自适应均匀网格和自适应非均匀网格概念，给出了改进型PCE方法的基本步骤；以固体火箭发动机低温点火过程中的不确定性分析为例，基于所提方法分析材料参数随机性对药柱输出响应的影响，获得输出响应的概率分布，验证了所提方法相对传统方法的优越性。

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]，然而分子需要运用蒙特卡洛法或高斯积分获得。高斯积分所用高斯节点值与积分权值通过文献[18]可查，6组典型高斯厄米特积分节点的节点值和积分权值如表1所示。其中，N为高斯节点数，ξq和wq(q=1,2,…，N)分别为一维高斯厄米特积分节点的节点值与积分权值。(8)式可以表示为

表1 高斯- 厄米特积分典型权值和积分点

(9)

式中：nm表示第m维变量上的第n个高斯节点；ξnm和wnm分别表示第m维变量上第n个高斯节点(即第nm个高斯节点)的节点值与积分权值；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*为自适应高斯节点数，R为1个常数，其初始值R0=0，ΔR为R的增量。

1.3.1 均匀网格法

均匀网格法基于全因子数值积分技术(FFNI)[23]，其具体步骤为：

1)建立均匀网格，即d维高斯节点的确定。将所得1维节点与权值进行直接张量积操作，得到d维节点与权值，其中，各维节点数相等N1=N2=…=Nd=K.

2)考虑非正态变量的相关性，根据Nataf变换的逆变换方法[24]，将每个高斯节点从标准随机空间(ξ空间)变换到原随机空间(x空间)：x1m=F1(ξ1m),…,xnm=Fn(ξnm),…,xkm=FK(ξKm)，其中，Fn(ξnm)表示将第m维变量中的第n个高斯节点ξnm，从标准随机空间变换到原随机空间。

3)将样本点x=[x1,…,xm,…,xd]代入实际响应模型，得到各高斯节点处的真实响应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，各维随机变量高斯节点数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，各维随机变量初始高斯节点数N=3，利用均匀网格法求解PCE系数。

②基于(10)式计算已得PCE模型均值μ与标准差σ，将结果代入(16)式[25]，得E1、E2.如果E1、E2满足(16)式，则自适应高斯节点数N*=N；否则，设N=N+1，返回步骤①，开始新运算，直至(16)式得到满足。

(16)

式中：根据文献[25]，e1、e2为区间(0，0.05]的常数。

3)自适应非均匀网格的确定。

①根据所得自适应PCE阶数P*和自适应均匀网格(N*)d，利用均匀网格法求解PCE系数。基于(13)式，得到各变量的总灵敏度指数ST1,…,STi,…,STd.

②基于(17)式，计算变量ξ1,…,ξi,…,ξd上的高斯节点数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为区间(0，0.03]的常数。

⑤基于自适应PCE阶数与自适应非均匀网格，求解PCE系数和输出响应统计矩。

2 固体火箭发动机低温点火不确定性量化

2.1 问题描述

某固体火箭发动机由壳体、绝缘层和六角星形药柱构成，且药柱材料为异佛尔酮二异氰酸酯/端羟基聚丁二烯(IPDI-HTPB)丁羟基推进剂。考虑到发动机几何构型和载荷的对称性，建立药柱结构的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 药柱横截面

以该固体发动机低温点火为例：点火前，药柱经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阶数的确定

设均匀网格方法中各变量的高斯节点数为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，基于均匀网格方法计算各变量高斯节点数分别为N=3,4,…时的PCE模型，求解其模型期望与标准差，计算结果如图5 所示。N=4时，E1=0.002 4<0.05，E2=0.033 7<0.05；由(16)式或图5可判断，N*=4为均匀网格方法中自适应高斯节点数，即4×4网格为所求自适应均匀网格。由于3×3均匀网格的重复使用，该过程的实际计算成本为4×4=16次真实响应模型运算。图5中PCE(P=3)表示3阶PCE代理模型方法。

图5 均匀网格自适应高斯节点数N*的确定

2.3.3 全局灵敏度分析

基于自适应PCE阶数P*=3和自适应均匀网格(N*)2=4×4的PCE模型对药柱的线膨胀系数与泊松比进行全局灵敏度分析，计算结果如图6 所示。

图6 输入变量灵敏度指数

2.3.4 自适应非均匀网格确定

基于1.3.2节中的步骤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个高斯节点，输入变量PRXY上取5个高斯节点，构成“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方法相近且远小于其他方法，表明药柱内表面最大等效应变响应的尖锐程度误差减小。综上所述，本文所提方法的不确定性分析结果优于其他方法。

图8 不同方法所得药柱内表面应变响应的PDF

图9 不同方法所得药柱内表面应变响应的CDF

表4 药柱内表面的最大等效应变的前4阶统计矩

3 结论

传统PCE求解方法，包括回归法和映射法(如FFNI方法)在求解PCE系数的过程中，对各输入变量等量采样，并未考虑各输入变量对输出响应结果影响程度的不同。不恰当的PCE阶数与采样点数会在PCE模型求解的过程中，引起过拟合与欠拟合问题。为提高PCE模型求解精度，本文提出了改进型PCE方法。该方法改变了传统映射法与回归法将输入变量同等对待的习惯，求解基于各输入变量全局灵敏度指数的自适应非均匀网格。进而基于此非均匀网格，采用映射法求解PCE系数。

本文基于固体火箭发动机低温点火过程中的不确定性材料参数(药柱的线膨胀系数与泊松比)对于药柱内表面最大应变响应灵敏度结果的不同，建立了4×5的自适应非均匀网格，以此求解PCE系数。结果表明，本文所提方法较传统方法，能够在服从效率要求的条件下，明显提高了计算精度，证明其具有良好的工程应用价值。但本文所提方法还只是在输入变量较少的低维问题中证明其有效，对于输入变量较多的高维问题，还有待进一步深入研究。