Uncertainty quantification of cylindrical test through Wiener chaos with basis adaptation and projection
-
摘要:
由于炸药爆轰现象的复杂性和人们对它的认知缺陷,其表征爆轰流体力学过程的物理数学模型具有较强的不确定性,要降低基于爆轰建模与模拟的数值结果做出决策的风险,量化和评估不确定输入对爆轰系统输出结果的影响尤为重要。本文中针对具有高维随机变量的爆轰问题的不确定度量化,使用自适应基函数的Wiener混沌方法、耦合旋转变换和投影方法,减少截断空间的长度。针对输入变量相关性,使用Rosenblatt变换使其相互独立。针对不符合标准正态分布的变量使用等概率原则,将它化为标准正态分布。最后,使用自主研发的具有完全知识产权的爆轰数值模拟软件LAD2D, 研究了具有高维不确定参数的圆筒实验的不确定度量化,给出期望、标准差、置信区间等统计信息,所得问题与实验数据比对,从而确认了模型的有效性。
-
关键词:
- 圆筒模型 /
- 不确定度量化 /
- 自适应基函数 /
- JWL状态方程 /
- Rosenblatt变换
Abstract:The mathematical-physical model used to describe the detonation dynamics has many uncertain factors due to the complexity and lack of knowledge for detonation phenomenon. Quantifying and assessing the impact of input uncertainties on output of detonation systems has a direct influence on reducing the risk based on the numerical model and simulation results for detonation. The Wiener chaos based on adapted basis is used to deal with the uncertainty quantification of high-dimensional random variables for detonation simulation. The rotation transformation and projection method is used to reduce the length of truncation number. Rosenblatt transformation is used to transform the set of dependent random variables into independent random variables. The equality of probability principle is used to change the non-Gaussian random variables into standard random variables. Uncertainty quantifications of the cylinder test with high dimensional input uncertainties are studied. The statistical informations such as mean, standard deviations, and confidence intervals are presented. The simulation results coincide with the experimental data, and the accuracy of the model is validated.
-
Key words:
- cylinder test /
- uncertainty quantification /
- adapted basis /
- JWL EOS /
- Rosenblatt transform
-
爆轰是极其复杂的物理化学过程,发生在极短的时间尺度(10~20 μs)和极小的空间尺度((1~9)×10−4 m),是爆炸力学的重点和难点。理论、实验、数值模拟是研究爆轰的三种途径,且这三类方法各有千秋,又相辅相成。实验能直接反应真实物理现象,探索爆轰的发生机理。但受环境和仪器的限制,往往只能表征初始状态和最终结果. 特别是部分炸药由于感度过高,又导致实验人员处于危险的操作环境。数值模拟可以全时空的重复模拟爆轰的发生过程,是研究爆轰的一条安全经济方法。同时数值模拟可以与理论、实验相互印证,重视实验结果,并给出了没有实验结果的预测[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方程耦合化学反应的常微分方程以及复杂非线性状态方程,在拉格朗日坐标系下的质量守恒、动量守恒、能量守恒、状态方程和反应率方程分别为:
∂ρ∂t+∇⋅(ρu)=0 (1) ∂(ρu)∂t+∇⋅(ρuu)+∇p=0 (2) ∂(ρE)∂t+∇⋅(ρEu)+∇⋅(ρu)=0 (3) p={p(ρ,e)非炸药p(ρ,e,F)炸药 (4) dFdt=R(p,e,F) (5) 式中:ρ、u、E、e、p分别表示密度、速度、总能、内能与压力,uu是指并矢张量,u·u是指内积,E=e+1/2u·u,0≤F≤1为爆轰产物的燃烧函数。
1.2 反应率方程
采用Wilkins反应率模型研究爆轰, 其形式为:
F=[max(F1,F2)]nb (6) 式中:F1为CJ比容燃烧函数,F2为时间燃烧函数, nb可调参数。有:
F1={0v≥v0(v0−v)/(v0−vJ)vJ<v<v01v≤vJ (7) F2={0t≤tb(t−tb)/ΔLtb<t<tb+ΔL1t≥tb+ΔL (8) 式中:v=1/ρ表示比容,v0是初始比容,vJ=γv0/(γ+1)是CJ比容,γ是多方指数(理想气体常数),tb是起爆时间,指爆轰波到达计算网格的时间,通过惠更斯原理计算[12]。ΔL=rbΔR/DJ,ΔR表示网格宽度,DJ是爆速,rb可调。
1.3 物态方程
爆炸产物常采用JWL状态方程,其形式为:
p=A(1−ωR1V)e−R1V+A(1−ωR2V)e−R2V+ωeV (9) 式中:相对体积V=v/v0, 唯象参数A、B、R1、R2、ω可调且相关, 利用爆轰波阵面上的守恒关系以及CJ爆轰条件,得:
Aexp(−R1VJ)+Bexp(−R2VJ)+CVJ−(ω+1)=pJ (10) AR1exp(−R1VJ)+BR2exp(−R2VJ)+C(1+ω)VJ−(ω+2)=ρ0D2J (11) AR1exp(−R1VJ)+BR2exp(−R2VJ)+CωVJ−ω=ρ0Q+12pJ(1−VJ) (12) 式中: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符号 不确定度 描述 概率分布 ξ1 nb 可调参数 nb\simfont~B[αnb,βnb,anb,bnb] ξ2 rb 可调参数 rb\simfont~B[αrb,βrb,arb,brb] ξ3 R1 JWL EOS系数 R1\simfont~B[αR1,βR1,aR1,bR1] ξ4 R2 JWL EOS系数 R2\simfont~B[αR2,βR2,aR2,bR2] ξ5 ω JWL EOS系数 ω\simfont~B[αω,βω,aω,bω] ξ6 aNR N-R人为黏性系数 aNR\simfont~B[αaNR,βaNR,aaNR,baNR] ξ7 aL Landshoff人为黏性系数 aL\simfont~B[αaL,βaL,aaL,baL] ξ8 γ Gruneissen系数 γ\simfont~B[αγ,βγ,aγ,bγ] ξ9 σ 体积起爆阈值 σ\simfont~B[ασ,βσ,aσ,bσ] ξ10 ρ TNT初始密度 N(μρ, σ2ρ) ξ11 tb 起爆时间 tb\simfont~B[αtb,βtb,atb,btb] 表中,
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。Beta分布函数的参数的确定取决于专家建议和工程经验, 本文中参数具体取值如下:
αnb=2,βnb=2,anb=1.6,bnb=3.0;αrb=2,βrb=2,arb=0.9,brb=1.9;αrb=2,βrb=5,aR1=4.6,bR1= 5.2;αR2=2,βR2=5,aR2=1.2,bR2=1.6;αω=2,βω=5,aω=0.2,bω=0.5;αaNR=2,βaNR=2,aaNR=0.1,baNR= 0.3;αaL=2,βaL=2,aaL=0.05,baL=0.07;αγ=2,βγ=2,aγ=3.0,bγ=3.2;ασ=2,βσ=3,aσ=1.02,bσ= 1.04;αtb=5,βtb=2,atb=2.5,btb=23.2。 本文中不确定度的概率密度函数如图1所示。
2. 基于自适应基函数的模型简化
2.1 Wiener-Hilbert多项式混沌
令
(Ω,F,P) 代表完备概率空间, 其中Ω为样本空间,F 是Ω上的σ代数,P是定义在F 上的概率测度,本文使用Gauss测度。θ∈Ω为基本事件。ξ={ξ1, ξ2, ···, ξd}代表Ω上的分量服从标准正态分布的随机向量,且分量独立同分布。本文中选取d=11。令H代表Gauss空间,F (H)代表H上的σ代数。H:n:代表n次齐次Wiener噪声集合。同时令L2(Ω)= L2(Ω,F (H), P)代表Ω上的平方可积空间,赋予内积<fg>=∫Ωf(θ)g(θ)dP(θ) ,则构成一个Hilbert空间。由基础概率论,知<f>=∫Ωf(θ)dP(θ)=Ef ,E 表示期望算子, 且满足L2(Ω)=⊕nH:n:。U(t, x, ξ)为(Ω,
F ,P )上的随机函数,且满足式(1)~(12),代表密度、压力、速度分量、总能、管壁位置等响应量。均方可积函数U(t, x, ξ)∈L2([0, T]×O ×Ω,F ([0, T]×O ×Ω), dt×μ×P ),其中F ([0, T]×O ×Ω)代表[0, T]×O ×Ω上的σ代数,μ是Rk 上的Lebesgue测度,dt×μ×P 代表F ([0, T]×O ×Ω)上的乘积测度,且满足:∫T0E(‖ 则根据Cameron-Martin定理[13-15], U(t, x, ξ)具有如下形式:
U(t,{{x}},{{\xi }}) = \sum\limits_{{{\alpha }} \in {\mathcal{I}_{\rm{p}}}} {{U_{\rm{\alpha }}}(t,{{x}}){\varPsi _{\rm{\alpha }}}({{\xi }})} (13) 其中α=(α1, α2, …, αd)表示指标集,
\mathcal{I}_{\rm{p}} ={α|α1+α2+···+αd≤p},且:{\varPsi _{\rm{\alpha }}} = \frac{{{H_{\rm{\alpha }}}({{\xi }})}}{{\sqrt {{{\alpha }}!} }},\;\quad{{\alpha }}! = \prod\limits_{i = 1}^d {{\alpha _i}!} ,\;\quad{\left\| {{H_{\rm{\alpha }}}({{\xi }})} \right\|_{{{{L}}^2}({\rm{\Omega}})}} = {{\alpha }}! 式中:Hα(ξ)代表d维Hermite多项式。
在实际应用中,式(13)需要截断有限次:
{U^P}(t,{{x}},{{\xi }}) = \sum\limits_{|{{\alpha }}| \leqslant P} {{U_{\rm{\alpha }}}(t,{{x}}){\varPsi _{\rm{\alpha }}}} \quad\quad|{{\alpha }}| = {\alpha _1} + {\alpha _2} + \cdots + {\alpha _{{d}}} (14) 且根据Cameron-Martin定理,几乎处处,
\mathop {\lim }\limits_{P \to \infty } {U^P}(t,{{x}},{{\xi }}) = U(t,{{x}},{{\xi }}) 。2.2 基变换
A表示
{R}^{\rm d} 上的正交变换,满足AAT=I,I表示{{{R}}^{\rm{d}}} 上的单位矩阵,定义:{{\eta }} = {{A\xi }} (15) 因此η也是H的一组基,且H:n:也可以由η生成,且令:
{\varPsi _{\rm{\alpha }}}^{({{A}})}({{\xi }}) = {\varPsi _{\rm{\alpha }}}({{\eta }}),\quad{U^{({{A}})}}(t,{{x}},{{\eta }}) = U(t,{{x}},{{\xi}}) (16) 因此:
{U^{({{A}})}}(t,{{x}},{{\eta }}) = \sum\limits_{{{\beta }} \in {\mathcal{I}_{\rm{P}}}} {{U_{\rm{\beta }}}^{({{A}})}(t,{{x}}){\varPsi _{\rm{\beta }}}({{\eta }})} (17) 式中:
{U_{\rm{\beta }}}^{({{A}})}(t,{{x}}) = \displaystyle\sum\limits_{{{\alpha }} \in {\mathcal{I}_{{P}}}} {{U_{\rm{\alpha }}}(t,{{x}}) < {\varPsi _{\rm{\alpha }}}{\varPsi _{\rm{\beta }}}^{({{A}})} > } 。2.3 基于投影法的模型简化
令
{{V}_{\mathcal{I}}}=\mathcal{L}({{\varPsi }_{\rm{\beta }}}^{\left({A} \right)}),{\beta }\in \mathcal{I}, \mathcal{I}\in {{\mathcal{I}}_{\rm{P}}} ,则U(A)(t, x, η)在{V}_{\mathcal{I}} 上的{{U}^{\left({A},\mathcal{I} \right)}}\left(t,{x},{\eta } \right) 投影定义如下:{U^{({{A}},\mathcal{I})}}(t,{{x}},{{\eta }}) = \sum\limits_{{{\beta }} \in \mathcal{I}} {\sum\limits_{{{\alpha }} \in {\mathcal{I}_{{P}}}} {{U_{\rm{\alpha }}}(t,{{x}}){\varPsi _{\rm{\beta }}}({{\eta }}) < {\varPsi _{\rm{\alpha }}}{\varPsi _{\rm{\beta }}}^{({{A}})} > } } (18) 另外,U(A)(t, x, η)在
{V}_{\mathcal{I}} 上的投影,又表示为:{U^{({{A}},\mathcal{I})}}(t,{{x}},{{\eta }}) = \sum\limits_{{{\gamma }} \in \mathcal{I}} {U_{\bf{\gamma }}^\mathcal{I}(t,{{x}}){\varPsi _{\bf{\gamma }}}({{\xi }})} 从而:
U_{\bf{\gamma }}^\mathcal{I}(t,{{x}}) = \sum\limits_{{{\gamma }} \in \mathcal{I}} {\sum\limits_{{{\beta }} \in \mathcal{I}} {\sum\limits_{{{\alpha }} \in {\mathcal{I}_{{P}}}} {{U_{\bf{\alpha }}}(t,{{x}}) < {\varPsi _{\bf{\alpha }}}{\varPsi _{\bf{\beta }}}^{({{A}})} > < {\varPsi _{\bf{\beta }}}{\varPsi _{\bf{\gamma }}} > } } } (19) 将U限制在
{U^{\mathcal{I}}}=\left\{ {{U}_{\rm{\gamma }}},\gamma \in \mathcal{I} \right\} 。2.4 非嵌入方法求解高斯自适应基函数的系数
通过第2.3节的叙述,可以看到,模型简化的关键在于A和
\mathcal{I} 的选取。本文中使用高斯自适应方法选取A和\mathcal{I} 。若U在高斯空间H的投影已知,则A构造如下。首先令
{\eta _1} = \sum\limits_{i = 1}^d {{U_{{{\bf{e}}_{\rm{i}}}}}(t,{{x}})} {\xi _{\rm{i}}} (20) 式中:ei=(0, ···, 1, ···, 0),第i个位置为1,其余位置全为0。
不难看出η1包含U的全部高斯统计量。η的其余分量通过Gram-Schmidt方法求解, 从而解出A。对于
\mathcal{I} 的选取,令其包含η1的小于等于P的指标, 即\mathcal{I}=\mathcal{E}_1 ,\mathcal{E}_1 是\mathcal{I}_{\rm{p}} 的子集, 第i个位置不为0,其余位置全为0,因此ei是\mathcal{E}_1 的单位向量。此时\begin{split} {U^{({{A}})}}(t,{{x}},{{\eta }}) =& U_0^{({{A}})} + {U_{{{{e}}_1}}}(t,{{x}}){\eta _1} + \sum\limits_{1 < |{{\beta }}| \leqslant P} {{U_{\bf{\beta }}}^{({{A}})}(t,{{x}}){\varPsi _{\bf{\beta }}}({{\eta }})} = \\ & U_0^{({{A}})} + {U_{{{{e}}_1}}}(t,{{x}}){\eta _1} + \sum\limits_{1 < |{{\beta }}| \leqslant P,{{\beta }} \in {\varepsilon _1}} {{U_{\bf{\beta }}}^{({{A}})}(t,{{x}}){\varPsi _{\bf{\beta }}}({{\eta }})} + \sum\limits_{1 < |{{\beta }}| \leqslant P,{{\beta }} \notin {\varepsilon _1}} {{U_{\bf{\beta }}}^{({{A}})}(t,{{x}}){\varPsi _{\bf{\beta }}}({{\eta }})} \end{split} (21) 且:
{U^{({{A}},\mathcal{I})}}(t,{{x}},{{\eta }}) = U_0^{({{A}})} + \sum\limits_{{{\beta }} \in {\varepsilon _1}} {{U_{\bf{\beta }}}^{({{A}})}(t,{{x}}){\varPsi _{\bf{\beta }}}({{\eta }})} (22) 误差为:
{U^{({{A}})}}(t,{{x}},{{\eta }}) - {U^{({{A}},\mathcal{I})}}(t,{{x}},{{\eta }}) = \sum\limits_{1 < |{{\beta }}| \leqslant P,{{\beta }} \notin {\varepsilon _1}} {{U_{\bf{\beta }}}^{({{A}})}(t,{{x}}){\varPsi _{\bf{\beta }}}({{\eta }})} (23) 误差量的期望为0,且与高斯量正交。
式(22)的系数通过非嵌入方法求解,方法如下:
{U_{\rm{\beta }}}^{({{A}})} \approx \sum\limits_{r = 1}^s {{U^{({{A}},\mathcal{I})}}(t,{{x}},{{{\eta }}^{(r)}})} {\varPsi _{\rm{\beta }}}({{{\eta }}^{(r)}}){w_{\rm{r}}}\quad{{\beta }} \in {\varepsilon _1} (24) 式中:η(r)=(η1(r), 0, ···, 0),η(r)和wr分别为求积点和权重,s为求积点的个数。且U(t, x, ξ)满足式(1)~(12), 并满足关系式
{{\xi }} = {{{A}}^{ - 1}}{{\eta }} = {{{A}}^T}{{\eta }} 。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=FX−1FY(y)。
假设X~
\mathcal{N} (0, 1),则:x=\varPhi _{{X}}^{-1}{{F}_{{Y}}}\left( y \right) 若Y~
\mathcal{N} (μ,σ2), 则X=(Y-μ)/σ服从标准正态分布。其次,使用Rosenblatt变换[16],将一列相关随机变量化为服从正态分布的独立随机变量组。利用概率等交换原则:
\begin{aligned} {{F}_{{{Y}_{1}}}}\left({{y}_{1}} \right)&={{F}_{{{X}_{1}}}}\left({{x}_{1}} \right)\\ {{F}_{{{Y}_{2}}}}\left({{y}_{2}}|{{y}_{1}} \right)&={{F}_{{{X}_{2}}}}\left({{x}_{2}}|{{x}_{1}} \right)\\ &\;\;\vdots \\ {{F}_{{{Y}_{{n}}}}}\left({{y}_{{n}}}|{{y}_{1}},{{y}_{2}},\cdots ,{{y}_{{n-1}}} \right)&={{F}_{{{X}_{{n}}}}}\left({{x}_{{n}}}|{{x}_{1}},{{x}_{2}},\cdots ,{{x}_{{n-1}}} \right)\end{aligned} 其中:
{{F}_{{{Y}_{\rm{n}}}}}\left({{y}_{{n}}}|{{y}_{1}},{{y}_{2}},\cdots ,{{y}_{{n-1}}} \right)=P\left\{ {{Y}_{{n}}}\text{≤} {{y}_{{n}}}|{{Y}_{1}}={{y}_{1}},{{Y}_{2}}={{y}_{2}},\cdots ,{{Y}_{{n-1}}}={{y}_{{n-1}}} \right\} 表示条件概率。从而:
\begin{aligned} & {{y}_{1}}=F_{{{Y}_{1}}}^{-1}\left( {{F}_{{{X}_{1}}}}\left( {{x}_{1}} \right) \right) \\ & {{y}_{2}}={{F}_{{{Y}_{2}}}}\left( {{F}_{{{X}_{2}}}}\left( {{x}_{2}}|{{x}_{1}} \right)|{{y}_{1}} \right) \end{aligned} (25) {{y}_{{n}}}=F_{\rm{Y_n}}^{-1}\left({{F}_{{{X}_{{n}}}}}\left({{x}_{{n}}}|{{x}_{1}},{{x}_{2}},\cdots ,{{x}_{{n-1}}} \right)|{{y}_{1}},{{y}_{2}},\cdots ,{{y}_{{n-1}}} \right) 则{Y1, Y2, ···, Yn}服从标准正态分布独立同分布。
进而,利用Rosenblatt变换,可将不确定参数化为一列服从标准正态分布的独立同分布随机变量。
3. 圆筒实验的不确定度量化
圆筒实验是确定炸药爆轰产物JWL状态方程和评估炸药做功能力的标准化实验,应用广泛。实验原理是将炸药放入等壁厚的铜质圆筒中,从圆筒的一端将其引爆,利用高速转镜式扫描相机记录筒壁在爆轰产物驱动下的膨胀过程。圆筒实验布局见图2,炸药尺寸为
\varnothing ∅ 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狭缝处,筒壁在爆轰产物驱动下的膨胀过程,径向壁位置与径向壁速度随时间的演化过程。利用本文所述的方法,给出了圆筒实验的不确定度量化结果,可以从多角度观测输入不确定度对输出结果特别是速度和管壁位置的影响。具体内容如图4~8所示。图4~5给出了管壁位置和速度随时间变化的期望值及标准差。可以看出,爆轰波在25 μs到达管壁,并继续向前传播,并在30 μs附近速度到达峰值,管壁在惯性作用下继续向前运动,同时速度峰值最大的点也是速度标准差最大的点。从图5可以看出,本文的计算格式在强间断处带有明显的波后震荡特征,这可能也是激波位置标准差最大的部分原因。而管壁位置的期望和标准差在激波达到后呈现一定的类线性关系。
图6给出了圆筒管壁位置和速度随时间变化的置信区间的计算结果,区间宽度满足3σ法则,并将置信区间涂黄。计算结果看出,爆轰波过后,置信区间宽度变大。波前标准差为0,因此区间宽度也是0。在图7中计算结果和实验结果对比发现,实验样本均落在置信区间内,符合预期。为观测方便,图8为图7的局部放大,进一步验证了观测结论。同时,也确认了本文方法的有效性。也说明,本文中参数在特定的有限时间内可以在一定范围内变动,均可符合精度要求。同时,进一步观测试验数据,发现爆轰波到达后,速度的实验数据也有明显的震荡,这与标准差的震荡吻合。说明我们的算法是可信的。
由于本文为含有11个不确定性参数的爆轰系统,标准的4次展开PC截断长度为(11+4)!/(11!4!)−1=1 365,通过自适应基变换和投影变换,截断长度为(1+4)!/(1!4!)−1=4,从而在一定程度上提高了计算速度,缓解了维数灾难。
4. 结论与展望
通过Wiener-Hilbert空间中基函数的旋转变换和投影变换,明显减少了截断长度,从而减轻了计算量。并将此方法应用于爆轰圆筒实验,给出了不确定度量化和传播结果。结果以期望、标准差、置信区间表示,并与实验数据做比较,发现实验数据均落在置信区间内。同时,观测到波后速度震荡剧烈与计算结果中强间断出数值模拟标准差达到峰值吻合。从而确认了模型的有效性。
下一步,将本文应用到其他模型并与现有的结果对比, 同时开展如下工作。
(1)本文中研究局限在Wiener-Hilbert空间,导致随机变量必须服从正态分布。对于非正态分布变量通过等概率原则和Rosenblatt变换将其转化成正态分布变量,这增加了计算量。而众所周知的是,在PC理论中均匀分布对应勒让德多项式,Beta分布对应Jacob多项式等。能否在其他完备的Banach空间中研究此自适应投影问题,从而不需要转化,减少计算量?若把唯象参数当做认知不确定度,如何处理爆轰系统的不确定度量化?这仍是一个具有挑战性的课题。
(2)波后震荡与实验数据吻合,没必要消除,但是能否使用高精度格式改进模型算法,减少波后震荡?需进一步提高模型的精确性。
(3)本文中主要研究参数和物理量不确定的爆轰系统的不确定度量化,然而爆轰系统的复杂性导致状态方程和反应率方程的选择也是不确定的,如何研究不同形式状态方程和反应率方程对系统的影响,即模型形式不确定度的量化?这也是一个具有挑战性的课题。
(4)仅研究炸药状态方程和反应率方程中的不确定参数是不够的。在爆轰实验中,无论是光学扫描方法还是激光干涉仪,仍然有很多不确定因素值得研究,如测量不确定度等。需要进一步深入分析实验中的不确定度来源,量化不确定度的传播,提高数值模拟的可信度和模型的预测能力。因此,实验不确定度是下一步研究的课题。
-
表 1 爆轰流体力学中的不确定度来源
Table 1. Sources of uncertainty in detonation hydrodynamics
符号 不确定度 描述 概率分布 {\xi _1} {n_{\rm{b}}} 可调参数 {n_{\rm{b}}}\simfont\text{~} {B}[{\alpha _{{{n}_{\rm{b}}}}},{\beta _{{{n}_{\rm{b}}}}},{a_{{{n}_{\rm{b}}}}},{b_{{{n}_{\rm{b}}}}}] {\xi _2} {r_{\rm{b}}} 可调参数 {r_{\rm{b}}}\simfont\text{~} {B}[{\alpha _{{{r}_{\rm{b}}}}},{\beta _{{{r}_{\rm{b}}}}},{a_{{{r}_{\rm{b}}}}},{b_{{{r}_{\rm{b}}}}}] {\xi _3} {R_1} JWL EOS系数 {R_1}\simfont\text{~} {B}[{\alpha _{{\rm{R}_1}}},{\beta _{{\rm{R}_1}}},{a_{{\rm{R}_1}}},{b_{{\rm{R}_1}}}] {\xi _4} {R_2} JWL EOS系数 {R_2}\simfont\text{~} {B}[{\alpha _{{\rm{R}_2}}},{\beta _{{\rm{R}_2}}},{a_{{\rm{R}_2}}},{b_{{\rm{R}_2}}}] {\xi _5} \omega JWL EOS系数 \omega \simfont\text{~} {B}[{\alpha _{{\omega}} },{\beta _{{\omega}} },{a_{{\omega}} },{b_{{\omega}} }] {\xi _6} {a_{\rm{NR}}} N-R人为黏性系数 {a_{\rm{NR}}}\simfont\text{~} {B}[{\alpha _{{{a}_{\rm{NR}}}}},{\beta _{{{a}_{\rm{NR}}}}},{a_{{{a}_{\rm{NR}}}}},{b_{{{a}_{\rm{NR}}}}}] {\xi _7} {a_{\rm{L}}} Landshoff人为黏性系数 {a_{\rm{L}}}\simfont\text{~} {B}[{\alpha _{{{a}_{\rm{L}}}}},{\beta _{{{a}_{\rm{L}}}}},{a_{{{a}_{\rm{L}}}}},{b_{{{a}_{\rm{L}}}}}] {\xi _8} γ Gruneissen系数 \gamma \simfont\text{~} {B}[{\alpha _{{\gamma}} },{\beta _{{\gamma}} },{a_{{\gamma}} },{b_{{\gamma}} }] {\xi _9} σ 体积起爆阈值 \sigma \simfont\text{~}{B}[{\alpha _\sigma },{\beta _\sigma },{a_\sigma },{b_\sigma }] {\xi _{10}} \rho TNT初始密度 {N}\left({{\mu }_{{\rho}}},\ \sigma _{{\rho}} ^{2} \right) {\xi _{11}} {t_{\rm{b}}} 起爆时间 {{t}_{\rm{b}}}\simfont\text{~}{B}[{{\alpha }_{{{{t}}_{\rm{b}}}}},{{\beta }_{{{{t}}_{\rm{b}}}}},{{a}_{{{{t}}_{\rm{b}}}}},{{b}_{{{{t}}_{\rm{b}}}}}] -
[1] 章冠人, 陈大年. 凝聚炸药起爆动力学[M]. 北京: 国防工业出版社, 1991. [2] 孙锦山. 凝聚炸药非理想爆轰的数值模拟 [J]. 力学进展, 1995, 25(1): 127–133. DOI: 10.6052/1000-0992-1995-1-J1995-043SUN Jinshan. Numerical modeling of non-ideal detonation in condensed explosives [J]. Advances in Mechanics, 1995, 25(1): 127–133. DOI: 10.6052/1000-0992-1995-1-J1995-043 [3] 王瑞利, 江松. 多物理耦合非线性偏微分方程与数值解不 [J]. 中国科学: 数学, 2015, 45(6): 723–738. DOI: 10.1360/N012014-00115WANG Ruili, JIANG Song. Mathematical methods for uncertainty quantification in nonlinear multi-physics systems and their numerical simulations [J]. Scientia Sinica: Mathematica, 2015, 45(6): 723–738. DOI: 10.1360/N012014-00115 [4] 王成, SHU Chiwang. 爆炸力学高精度数值模拟研究进展 [J]. 科学通报, 2015, 60(10): 882–898. DOI: 10.1360/N972014-00936WANG Cheng, SHU Chiwang. Progresses in high-resolution numerical simulation of explosion mechanics [J]. Chinese Science Bulletin, 2015, 60(10): 882–898. DOI: 10.1360/N972014-00936 [5] 梁霄, 王瑞利. 爆轰流体力学模型敏感度分析与模型确认 [J]. 物理学报, 2017, 66(11): 116401. DOI: 10.7498/aps.66.116401LIANG Xiao, WANG Ruili. Sensitivity analysis and validation of detonation computational fluid dynamics model [J]. Acta Physica Sinica, 2017, 66(11): 116401. DOI: 10.7498/aps.66.116401 [6] 梁霄, 王瑞利. 混合不确定度量化方法及其在计算流体动力学迎风格式中的应用 [J]. 爆炸与冲击, 2016, 36(4): 505–519. DOI: 10.11883/1001-1455(2016)04-0509-07LIANG Xiao, WANG Ruili. Mixed uncertainty quantification and its application in upwind scheme for computational fluid mechanics [J]. Explosion and Shock Waves, 2016, 36(4): 505–519. DOI: 10.11883/1001-1455(2016)04-0509-07 [7] 汤涛, 周涛. 不确定性量化的高精度数值方法和理论 [J]. 中国科学: 数学, 2015, 45(7): 891–928. DOI: 10.1360/N012014-00218TANG Tao, ZHOU Tao. Recent developments in high order numerical methods for uncertainty quantification [J]. Scientia Sinica Mathematica, 2015, 45(7): 891–928. DOI: 10.1360/N012014-00218 [8] GHANEM R, SPANOS R. Stochastic finite elements: A spectral approach[M]. New York: Springer, 1991. [9] 梁霄, 王瑞利. 爆炸波中的混合不确定度量化方法 [J]. 计算物理, 2017, 34(5): 574–582. DOI: 1001-246X(2017)05-0574-09LIANG Xiao, WANG Ruili. Mixed uncertainty quantification of blast wave problem [J]. Chinese Journal of Computational Physics, 2017, 34(5): 574–582. DOI: 1001-246X(2017)05-0574-09 [10] HUAN X, GERACI G, SAFTA C, et al. Multifidelity statistical analysis of large eddy simulation in Scramjet computations[C]//AIAA non-deterministic Approaches Conference. 2018: 1−11. DOI: 10.2514/6.2017-1089. [11] HUAN X, SAFTA C, SARGSYAN K, et al. Compressive sensing with cross-validation and stop-sampling for sparse polynomial chaos expansions [J]. SIAM/ASA Journal on Uncertainty Quantification, 2018, 6(2): 907–936. DOI: 10.1137/17M1141096. [12] 王瑞利, 林忠, 魏兰, 等. 利用前沿推进法计算复杂炸药结构的起爆时间 [J]. 高压物理学报, 2015, 29(4): 286–292. DOI: 10.11858/gywlxb.2015.04.008WANG Ruili, LIN Zhong, WEI Lan, et al. Calculating the initiation time of the explosive with complex structure using the advancing front technique [J]. Chinese Journal of High Pressure Physics, 2015, 29(4): 286–292. DOI: 10.11858/gywlxb.2015.04.008 [13] TIPIREDDY R, GHANEM R. Basis adaptation in homogeneous chaos spaces [J]. Journal of Computational Physics, 2014, 259(3): 304–317. DOI: 10.1016/j.jcp.2013.12.009. [14] TSILIFIS P, GHANEM R. Reduced Wiener chaos representation of random fields via basis adaptation and projection [J]. Journal of Computational Physics, 2017, 341: 102–120. DOI: 10.1016/j.jcp.2017.04.009. [15] CAMERON R, MARTIN W. The orthogonal development of nonlinear functional in series of Fourier-Hermite functional [J]. Annals Mathematics, 1947, 48(2): 385–392. DOI: 10.2307/1969178. [16] ROSENBLATT M. Remarks on a multivariate transformation [J]. Annals of Mathematical Statistics, 1952, 23: 470–472. DOI: 10.1007/978-l-4419-8339-8_8. 期刊类型引用(4)
1. 梁霄,王朔,王瑞利,陈江涛. 基于时间演化响应面模型的JWL状态方程参数标定方法. 山东大学学报(理学版). 2024(04): 117-126 . 百度学术
2. 王瑞利,胡星志,王朔. 建模与模拟中验证和确认及不确定度量化关键方法. 数学建模及其应用. 2021(02): 1-16+107 . 百度学术
3. 梁霄,陈江涛,王瑞利. 高维参数不确定爆轰的不确定度量化. 兵工学报. 2020(04): 692-701 . 百度学术
4. 梁霄,陈江涛,王瑞利,胡星志. 非接触水下爆炸下舰船冲击环境的不确定度量化. 中国舰船研究. 2020(06): 128-136 . 百度学术
其他类型引用(0)
-