基于自适应和投影Wiener混沌的圆筒实验不确定度量化*
2019-06-05王瑞利
梁 霄，王瑞利
（1.山东科技大学数学学院，山东 青岛 266590；2.北京应用物理与计算数学研究所，北京 100089）
爆轰是极其复杂的物理化学过程，发生在极短的时间尺度(10～20 μs)和极小的空间尺度((1～9)×10-4m)，是爆炸力学的重点和难点。理论、实验、数值模拟是研究爆轰的三种途径，且这三类方法各有千秋，又相辅相成。实验能直接反应真实物理现象，探索爆轰的发生机理。但受环境和仪器的限制，往往只能表征初始状态和最终结果.特别是部分炸药由于感度过高，又导致实验人员处于危险的操作环境。数值模拟可以全时空的重复模拟爆轰的发生过程，是研究爆轰的一条安全经济方法。同时数值模拟可以与理论、实验相互印证，重视实验结果，并给出了没有实验结果的预测[1-6]。然而由于爆轰的极端复杂性以及认知的缺陷，导致多物理爆轰数学模型特征鲜明。首先，模型结构及其运算格式高度复杂，尤其高精度数值模拟的研发和运算成本高昂。其次，目前状态方程（EOS）和反应率方程均为唯象建模，含有大量经验参数，没有物理意义，根据经验将其限定在某个范围。最后，测量技术受外界因素的干扰，以及物质本身固有的不确定性也会导致待测物理量的随机波动。以上原因导致真实的爆轰模型是一个众多不确定因素耦合在一起的非线性耦合Euler方程。要发展高可信度的爆轰过程数值模拟程序，必须对炸药起爆、描述反应区的反应率函数、状态方程中的众多不确定性参数进行量化。
目前，受制于爆轰模型隐藏的大量不确定度和受限于爆轰多物理过程的计算成本，多随机因素扰动的爆轰建模与模拟的不确定度量化（uncertainty quantification，UQ）在国际上仍然是一个极具挑战性的课题，且实质性的结果不多。然而，量化和评估不确定因素对系统输出结果的影响，直接关系到模型的稳健性和数值模拟的可靠性。
研究爆轰模型的不确定度量化，普遍的做法是使用Monte Carlo方法(MC)。MC不需要任何假设，操作简单，不依赖于模型的几何结构，可以认为是系统在不确定意义下的精确解。然而，MC较慢的收敛速度(O(N-1/2))，导致采样点数量介于103～106之间的MC方法很难应用于需要较高运行机时的高精度爆轰格式，这是MC方法在处理高维不确定爆轰模型的局限性之一。Wiener提出的多项式混沌（polynomial chaos，PC）方法在满足特定条件下可以替代MC方法，Ghanem又将PC和有限元结合，从此PC在工程领域应用广泛[7-9]。按照求解系数方式的不同，PC又分为嵌入多项式混沌(IPC)和非嵌入多项式混沌(NIPC)。IPC需要不断地更改程序，而爆轰软件的开发往往是一个团队多年工作的结晶，因此此方法求解爆轰模型并不现实。复杂工程问题中常用的高维不确定度NIPC方法有求积法、回归方法、采样法等。但是鉴于高维不确定系统截断长度和维数灾难，NIPC很难直接应用于爆轰模型。
而从截断长度和取点规模来看，上述方法会遇到如下问题。求积法遇到的一个问题是维数灾难问题。简而言之，对于一个含有20个不确定性参数的爆轰系统，每个参数选择5个求积点，则一次积分需要运行系统520≈1014次，对于4阶展开需要运行系统约1020次。对于高精度复杂多物理爆轰过程这是不现实的。至于采样法，包含传统的MC采样和替代性的拉丁超立方体采样等，目前采样规模仍很庞大，不能直接适用于复杂高精度爆轰系统。对于回归方法，由于爆轰计算成本较高，采样点的个数小于切断长度，回归方法使用采样将其转化为一个欠定线性方程组，然后使用压缩感知方法求解相应的欠定线性方程组。此方法大幅度减少了采样的规模，Huan等[10-11]已经将其应用计算高度复杂的未来高超音速燃烧发动机数值模拟。但是这种方法至今没有经过严格数学推导的误差收敛速度和数值精度，理论上不具有完备性，使用存在未知的风险。
本文中，用自适应基函数和投影的多项式混沌方法处理含有高维不确定度的爆轰问题，并以圆筒模型为例，给出UQ结果。基本思想是：将高维随机变量做旋转变换，得到一个新的随机向量组，在此随机向量组的基础上，重构Hermite多项式，进而利用Cameron-Martin公式重新展开，同时选取部分基函数生成相应的线性子空间。最后，将展开式在子空间上投影。 一个明显的优势是可以减少截断长度，例如对于20个不确定性参数的爆轰系统，标准的4次展开PC截断长度为(20+4)!/(20!4!)-1=10 624，通过自适应基变换和投影变换，截断长度为(1+4)!/(1!4!)-1=4，从而在一定程度上提高了计算速度，缓解了‘维数灾难’，并且此方法的误差有精确的误差公式。
然而，高维不确定度爆轰模型自身独特的特征导致其还必须处理如下问题。首先，JWL状态方程中的参数是相关的，不具备独立同分布(independent identical distribution, IID)。其次，自适应基函数在Wiener-Hilbert空间是完备的，而基于经验的参数未必服从正态分布。再次，以往的假设参数服从均匀分布，但是均匀分布的强间断性，令其转化为正态分布有一点难度。
圆筒实验是研究爆轰波传播驱动过程的基准实验。在实验中，研究人员将药柱置于金属圆筒内，通过多普勒类光学速度计、狭缝扫描相机获取金属圆筒外界面位移曲线和壳体膨胀速度，完整检验爆轰模型的合理性。将本文所述方法应用于圆筒实验为基准实验的UQ研究。也可以为后续众多复杂物理过程的数值模拟提供条件，并能推广到其他武器物理过程的UQ。
1 爆轰流体力学模型及不确定度
1.1 爆轰流体力学方程组
爆轰流体模型为守恒原理的Euler方程耦合化学反应的常微分方程以及复杂非线性状态方程，在拉格朗日坐标系下的质量守恒、动量守恒、能量守恒、状态方程和反应率方程分别为：
式中：ρ、u、E、e、p分别表示密度、速度、总能、内能与压力，uu是指并矢张量，u·u是指内积，E=e+1/2u·u，0≤F≤1为爆轰产物的燃烧函数。
1.2 反应率方程
采用Wilkins反应率模型研究爆轰, 其形式为：
式中：F1为CJ比容燃烧函数，F2为时间燃烧函数,nb可调参数。有：
式中：v=1/ρ表示比容，v0是初始比容，vJ=γv0/(γ+1)是 CJ比容，γ是多方指数（理想气体常数），tb是起爆时间，指爆轰波到达计算网格的时间，通过惠更斯原理计算[12]。ΔL=rbΔR/DJ，ΔR表示网格宽度，DJ是爆速，rb可调。
1.3 物态方程
爆炸产物常采用JWL状态方程，其形式为：
式中：相对体积V=v/v0, 唯象参数A、B、R1、R2、ω可调且相关, 利用爆轰波阵面上的守恒关系以及CJ爆轰条件，得：
式中：pJ、VJ是炸药CJ状态下的爆压、比容，Q是爆热。
实际应用中，用燃烧函数F将炸药和爆轰产物的状态方程连接起来p=FpEOS计算反应区。
1.4 计算方法
本文的模型计算采用自主研发、具有完全知识产权的爆轰流体力学软件LAD2D。软件的基本算法结构如下：代码的空间离散化基于任意多边形非结构网格，将拉氏自相容有限体积方法及Landshoff一次黏性aL、冲击波黏性、von Neumann-Richtmyer二次黏性aNR、子网格压力、人工热流等抑制非物理振荡的方法，以及多介质滑移计算推广到任意多边形非结构网格，构造了一套自适应AMR任意多边形非结构网格、流体大变形强适应的高精度有限体积格式，求解全耦合质量、动量、能量守恒方程和化学反应率方程。并采用了网格大变形处理的邻域可变技术。时间离散化基于预设条件双时间步。利用全显式多步格式构造系统的隐式解，使其具有4阶时间精度，并与预设条件耦合。
1.5 爆轰模型中的不确定度
关于不确定度的描述，目前并没有统一的说法。这里不确定度有两种类型：唯象不确定性参数和不确定性物理量。不确定度量化就是考虑两种不确定度对系统响应量的影响。表1给出了爆轰流体力学模型的输入不确定度。
表1 爆轰流体力学中的不确定度来源Table 1 Sources of uncertainty in detonation hydrodynamics
表中，B[ας，βς，aς，bς]代表参数ς服从具有双参数ας和βς、取值范围在[aς，bς]的Beta分布。做线性变换Z=(ς-aς)/(bς-aς)，则Z满足标准Beta分布B[ας，βς]。由于测量的不精确性以及认知的局限性，爆轰流体力学中的唯象参数，由于没有物理意义，因此无法通过实验标定。这类参数的取值范围取决于工程设计的接受程度。特别指出，以往的研究假设这类唯象参数满足简单易行的均匀分布，而本文中假设参数满足Beta分布，原因在于均匀分布的密度函数的强间断性导致其不易转化为正态分布。此外，N(μρ,σρ2)表示初始密度ρ服从期望μρ、方差σρ2的正态分布。ρ与其余唯象参数不同，它具有明确的物理意义，产生于晶体的错位以及炸药凝结过程中颗粒的不均匀性。由于样本容量大，μρ和σρ2的选取来源于统计观测数据。炸药选取JOB-9003, 根据实验统计结果μρ=1.845 g/cm3,σρ=0.005 g/cm3。
图1 不确定度的概率密度函数Fig.1 Probability density function of uncertainty
Beta分布函数的参数的确定取决于专家建议和工程经验, 本文中参数具体取值如下：
本文中不确定度的概率密度函数如图1所示。
2 基于自适应基函数的模型简化
2.1 Wiener-Hilbert多项式混沌
令代表完备概率空间, 其中Ω为样本空间，是Ω上的σ代数，P是定义在上的概率测度，本文使用 Gauss测度。θ∈Ω为基本事件。代表Ω上的分量服从标准正态分布的随机向量，且分量独立同分布。本文中选取d=11。令H代表Gauss空间, F(H)代表H上的σ代数。H：n：代表n次齐次∫Wiener噪声集合。同时令L2(Ω)=L2(Ω, F (H),P)代表Ω上的∫平方可积空间，赋予内积
上的随机函数，且满足式(1)～(12)，代表密度、压力、速度分量、总能、管壁位置等响应量。均方可积函数U(t,x,ξ)∈L2([0,T]× O ×Ω,([0,T]× O ×Ω), dt×μ×P)，其中 F ([0,T]× O ×Ω)代表[0,T]× O ×Ω 上的 σ 代数，μ是Rk上的Lebesgue测度，dt×μ×P代表([0,T]× O ×Ω)上的乘积测度，且满足：
则根据Cameron-Martin定理[13-15],U(t,x,ξ)具有如下形式：
其中α=(α1,α2, …,αd)表示指标集, Ip={α|α1+α2+···+αd≤p}，且：
式中：Hα(ξ)代表d维 Hermite多项式。
在实际应用中，式(13)需要截断有限次：
且根据Cameron-Martin定理，几乎处处，。
2.2 基变换
A表示Rd上的正交变换，满足AAT=I，I表示Rd上的单位矩阵，定义：
因此η也是H的一组基，且H：n：也可以由η生成，且令：
因此：
式中
2.3 基于投影法的模型简化
令则U(A)(t,x,η)在VI上 的U(A,I)(t,x,η)投影定义如下：
另外，U(A)(t,x,η)在VI上的投影，又表示为：
从而：
将U限制在UI={Uγ,γ ∈ I}。
2.4 非嵌入方法求解高斯自适应基函数的系数
通过第2.3节的叙述，可以看到，模型简化的关键在于A和 I 的选取。本文中使用高斯自适应方法选取A和 I。若U在高斯空间H的投影已知，则A构造如下。
首先令
式中：ei=(0, ···, 1, ···, 0)，第i个位置为 1，其余位置全为 0。
不难看出η1包含U的全部高斯统计量。η的其余分量通过Gram-Schmidt方法求解, 从而解出A。对于 I 的选取，令其包含η1的小于等于P的指标, 即 I=ℰ1， ℰ1是 Ip的子集, 第i个位置不为0，其余位置全为0，因此ei是 ℰ1的单位向量。此时
且：
误差为：
误差量的期望为0，且与高斯量正交。
式(22)的系数通过非嵌入方法求解，方法如下：
式中：η(r)=(η1(r), 0, ···, 0)，η(r)和wr分别为求积点和权重，s为求积点的个数。且U(t,x,ξ)满足式 (1)～(12),并满足关系式 ξ=A-1η=ATη。
2.5 逆累积分布函数和Rosenblatt变换
由第2.1节知，自适应基函数理论成立的前提是{ξ1,ξ2, ···,ξd}为一列独立同分布的标准正态随机变量组。然而，这并不适用于本文的研究对象。由随机变量组（见表1），{ξ1,ξ2, ···,ξn}不仅含有非正态分布随机变量，而且即使服从正态分布也不是标准正态分布。由第1.3节式(10)～(12)知，即给定R1、R2、ω可计算对应的A、B，即A、B可由R1、R2、ω给出。因此A、B、R1、R2、ω并非相互独立。因此，必须对随机变量组做适当改进，方可使用第2.1～2.4节方法。
使用逆累积分布函数变换将一般Beta分布化成标准正态分布。具体步骤如下：
首先，设X、Y分别满足累积分布函数FX(x)、FY(y)，利用概率等交换原则，令FX(x)=FY(y)，则
假设X～ N (0, 1)，则：
若Y～ N (μ,σ2), 则X=(Y-μ)/σ服从标准正态分布。
其次，使用Rosenblatt变换[16]，将一列相关随机变量化为服从正态分布的独立随机变量组。利用概率等交换原则：
其中：
表示条件概率。从而：
则{Y1,Y2, ···,Yn}服从标准正态分布独立同分布。
进而，利用Rosenblatt变换，可将不确定参数化为一列服从标准正态分布的独立同分布随机变量。
3 圆筒实验的不确定度量化
图2 圆筒实验装置示意图Fig.2 Schematic diagram of cylinder test
圆筒实验是确定炸药爆轰产物JWL状态方程和评估炸药做功能力的标准化实验，应用广泛。实验原理是将炸药放入等壁厚的铜质圆筒中，从圆筒的一端将其引爆，利用高速转镜式扫描相机记录筒壁在爆轰产物驱动下的膨胀过程。圆筒实验布局见图2，炸药尺寸为 ∅ 25.0 mm×305 mm，炸药为JOB-9003炸药，实验测得爆速为8 712 m/s；圆筒内径25.0 mm，壁厚2.5 mm，材料为紫铜。狭缝扫描采用平行光后照明技术，狭缝位置距离起爆端200 mm。计算格式采用方法见第1.4节。图3是针对JOB-9003炸药25.0 mm圆筒实验测试结果，给出了距起爆端200 mm狭缝处，筒壁在爆轰产物驱动下的膨胀过程，径向壁位置与径向壁速度随时间的演化过程。
图3 JOB-9003圆筒实验结果Fig.3 Experimental results of JOB-9003 in cylinder test
利用本文所述的方法，给出了圆筒实验的不确定度量化结果，可以从多角度观测输入不确定度对输出结果特别是速度和管壁位置的影响。具体内容如图4～8所示。图4～5给出了管壁位置和速度随时间变化的期望值及标准差。可以看出，爆轰波在25 μs到达管壁，并继续向前传播，并在30 μs附近速度到达峰值，管壁在惯性作用下继续向前运动，同时速度峰值最大的点也是速度标准差最大的点。从图5可以看出，本文的计算格式在强间断处带有明显的波后震荡特征，这可能也是激波位置标准差最大的部分原因。而管壁位置的期望和标准差在激波达到后呈现一定的类线性关系。
图6给出了圆筒管壁位置和速度随时间变化的置信区间的计算结果，区间宽度满足3σ法则，并将置信区间涂黄。计算结果看出，爆轰波过后，置信区间宽度变大。波前标准差为0，因此区间宽度也是0。在图7中计算结果和实验结果对比发现，实验样本均落在置信区间内，符合预期。为观测方便，图8为图7的局部放大，进一步验证了观测结论。同时，也确认了本文方法的有效性。也说明，本文中参数在特定的有限时间内可以在一定范围内变动，均可符合精度要求。同时，进一步观测试验数据，发现爆轰波到达后，速度的实验数据也有明显的震荡，这与标准差的震荡吻合。说明我们的算法是可信的。
图4 管壁位置的期望与标准差Fig.4 Expectation and standard deviation of cylindrical wall position
由于本文为含有11个不确定性参数的爆轰系统，标准的4次展开PC截断长度为(11+4)!/(11!4!)-1=1 365，通过自适应基变换和投影变换，截断长度为(1+4)!/(1!4!)-1=4，从而在一定程度上提高了计算速度，缓解了维数灾难。
图5 管壁速度的期望与标准差Fig.5 Expectation and standard deviation of cylindrical wall velocity
图6 管壁位置和速度的置信区间Fig.6 Confidence intervals of position of cylindrical wall position and velocity
图7 实验数据和计算结果的置信区间Fig.7 Confidence intervals of experimental and simulation results
4 结论与展望
通过Wiener-Hilbert空间中基函数的旋转变换和投影变换，明显减少了截断长度，从而减轻了计算量。并将此方法应用于爆轰圆筒实验，给出了不确定度量化和传播结果。结果以期望、标准差、置信区间表示，并与实验数据做比较，发现实验数据均落在置信区间内。同时，观测到波后速度震荡剧烈与计算结果中强间断出数值模拟标准差达到峰值吻合。从而确认了模型的有效性。
下一步，将本文应用到其他模型并与现有的结果对比, 同时开展如下工作。
（1）本文中研究局限在Wiener-Hilbert空间，导致随机变量必须服从正态分布。对于非正态分布变量通过等概率原则和Rosenblatt变换将其转化成正态分布变量，这增加了计算量。而众所周知的是，在PC理论中均匀分布对应勒让德多项式，Beta分布对应Jacob多项式等。能否在其他完备的Banach空间中研究此自适应投影问题，从而不需要转化，减少计算量？若把唯象参数当做认知不确定度，如何处理爆轰系统的不确定度量化？这仍是一个具有挑战性的课题。
（2）波后震荡与实验数据吻合，没必要消除，但是能否使用高精度格式改进模型算法，减少波后震荡?需进一步提高模型的精确性。
(3)本文中主要研究参数和物理量不确定的爆轰系统的不确定度量化，然而爆轰系统的复杂性导致状态方程和反应率方程的选择也是不确定的，如何研究不同形式状态方程和反应率方程对系统的影响，即模型形式不确定度的量化？这也是一个具有挑战性的课题。
(4)仅研究炸药状态方程和反应率方程中的不确定参数是不够的。在爆轰实验中，无论是光学扫描方法还是激光干涉仪，仍然有很多不确定因素值得研究，如测量不确定度等。需要进一步深入分析实验中的不确定度来源，量化不确定度的传播，提高数值模拟的可信度和模型的预测能力。因此，实验不确定度是下一步研究的课题。
