-
舰船非接触水下爆炸是一系列复杂的非线性多物理过程，迄今仍无法得到全系统、全时空的理解[1-4]。在研究水下爆炸时，主要运用3种方法：试验、理论计算和数值仿真[1, 3, 5-6]。二战后，西方海军强国通过实船水下爆炸试验[7-8]来直观地观测爆炸毁伤效果和评估船船抗冲击性能，以此判断拟定设备装舰的可行性。然而，实船水下爆炸试验也存在固有缺陷，例如成本高昂、过程不可控、对海洋生态环境造成破坏等。数值仿真的优点是安全环保、成本低、过程可控，但建模和模拟（M&S）过程中却包含了大量不确定性因素，使得决策者对M&S方法的预测能力有所顾虑。
不确定度量化（uncertainty quantification，UQ）技术结合了试验与数值这2种方法的优点，运用此技术可提高数值模型的可信度和可靠性。近年来，UQ研究作为一门新兴学科，受到了欧美国家学者的高度关注，被广泛应用于核能[9-11]、安全[12-13]、航空航天[14-18]等重大工程领域。我国在实船水下爆炸试验领域起步较晚，可获得的样本有限，供借鉴的国外公开资料也极度匮乏，使得我国舰船非接触水下爆炸的UQ方法研究具备了广阔的应用前景。然而，有关舰船非接触水下爆炸的UQ研究至今未见相关报道，其中部分原因是水下爆炸的M&S过程复杂而独特，且无法照搬已有的成熟方法。
首先，由于水下爆炸M&S过程中不确定性因素繁多且类型不同，除了有物理量自身固有的波动性和测量技术误差导致的无法消除的不确定度外，还有基于拟合数据的需要，在上述过程中使用了没有物理意义的不确定度且无法通过试验标定的唯象参数； 其次，目前常用的UQ方法成立的前提条件是随机变量服从独立同分布（independent identical distribution, IID），但水下爆炸中的随机变量并不完全服从独立同分布； 再次，部分不确定的物理量要求严格非负，使得概率统计中常见的高斯分布无法直接应用，例如，若假设质量服从正态分布，则理论上样本取值会出现质量为负的非物理情况； 最后，若假设参数服从均匀分布，则容易满足参数的有界性要求，但均匀分布的概率密度函数的强间断性使其很难转化为正态分布。因此，选取合理且符合统计结果的概率分布，对舰船非接触水下爆炸的UQ研究至关重要。
对于UQ方法的选择，很自然地会想到蒙特卡洛（Monte Carlo，MC）方法，但MC方法有着收敛速度慢的缺陷，而能有效代替MC方法的是多项式混沌（polynomial chaos，PC）方法，它也是大规模工程计算中常用的方法[19-23]。然而，在舰船非接触水下爆炸M&S过程中有太多的不确定性因素，使得多元PC方法容易陷入“维数灾难”。简言之，若采用经典的5个求积点方法计算11维随机变量驱动的水下爆炸系统，则需要运行程序
${5^{11}} \approx 4.9 \times {10^7}$ 次，5阶多项式则需要展开（PC截断长度[24]）$( 11 + 5)!/ \left( {11!5!} \right) - 1 = 4\;367$ 次，共需运行程序约$5 \times {10^{11}}$ 次，超出了目前的计算能力。而基于自适应基函数的齐次Wiener混沌方法改进了PC方法，可缓解“维数灾难”问题，其核心思想是通过构造随机基函数的同构酉变换（unitary transformation）得到新随机基函数，待测物理量在新随机基函数展开下的概率集中在低维子空间。鉴于此，本文将重点研究使用基于自适应基函数的齐次Wiener混沌方法处理含高维不确定度的舰船非接触水下爆炸问题。通过设计一个简单的试验装置，给出系统输出量的期望值、标准差、置信区间等统计信息，进而分析、量化和评估不确定性因素对非接触水下爆炸下舰船冲击环境的影响，所得结果可提高数学模型的可靠性、可信度和预测能力，用以为船舶结构设计和设备上舰安装提供依据。
-
因为水下爆炸的压力与时间及位置有关，且不容易通过实测获得，所以需要通过如下经验公式来确定[1-2,25]。
$$p(r,t) = {p_{\rm{m}}}(r)f\left(\frac{t}{{\theta (r)}}\right)$$ (1) 式中：
${p_{\rm{m}}}$ 为峰值压力；t为时间；$r$ 为斜距；$\theta (r)$ 为衰减常数。根据Cole公式[1]，有$${p_{\rm{m}}} = {K_1}{\left(\frac{{\sqrt[3]{W}}}{r}\right)^\alpha }$$ (2) $$\theta (r) = {K_2}{\left(\frac{{\sqrt[3]{W}}}{r}\right)^\beta }\sqrt[3]{W}$$ (3) $$f(x) = a\exp ( - {a_1}x) + b\exp ( - {b_1}x)$$ (4) 式中：
$W$ 为TNT炸药质量； f (x)为经验函数；$a,\;b,\; {a_1},\;{b_1},\; {K_1},\;{K_2},\;\alpha ,\;\beta$ 均为经验参数，取值待定。美国海军舰载设备的抗冲击性能是舰船结构设计及设备装舰的依据。如图1所示，在舰船−水交界面（流固耦合处），入射波分解为2个部分：一是穿透甲板的折射波；二是反射回水中的反射波。反射后的净压力是入射波和反射波的代数运算结果。
Figure 1. Force analysis of experimental setup
令
${p_{\rm{i}}}(t)$ 和${u_{\rm{i}}}(t)$ 分别为入射波导致的压力和位移，${p_{\rm{f}}}(t)$ 和${u_{\rm{f}}}(t)$ 分别为反射波导致的压力和位移，${p_{\rm{t}}}(t)$ 为透射波导致的压力，$x(t)$ 为甲板外法线方向的位移，$m$ 为甲板面密度。由于船体较大，计算面密度时，弹簧系统试验装置的质量忽略不计。根据牛顿运动定律，垂直于甲板的外法线方向有$$m\frac{{{{\rm{d}}^2}x}}{{{\rm{d}}{t^2}}} = {p_{\rm{t}}}(t)$$ (5) $${p_{\rm{t}}}(t) = {p_{\rm{i}}}(t) - {p_{\rm{f}}}(t)$$ (6) $$\frac{{{\rm{d}}x}}{{{\rm{d}}t}} = \frac{{{\rm{d}}{u_{\rm{i}}}}}{{{\rm{d}}t}} - \frac{{{\rm{d}}{u_{\rm{f}}}}}{{{\rm{d}}t}}$$ (7) $${p_{\rm{i}}}(t) = \rho c\frac{{{\rm{d}}{u_{\rm{i}}}}}{{{\rm{d}}t}}$$ (8) $${p_{\rm{f}}}(t) = \rho c\frac{{{\rm{d}}{u_{\rm{f}}}}}{{{\rm{d}}t}}$$ (9) 式中：
$\rho $ 为海水密度，$\rho = 1\;027\;{\rm{kg}}/{{\rm{m}}^3}$ ；c为本地声速，$c = 1\;493$ m/s。其中${p_{\rm{i}}}(t)$ 由式(1)确定。如图1所示，甲板上的弹簧系统试验装置在透射波作用下产生简谐振动。令
$y (t)$ 为试验装置在外法线方向的绝对位移，${\textit{z}}(t) = y(t) - x(t)$ 为试验装置在外法线方向的相对位移。根据牛顿第二定律，则有$$M\frac{{{{\rm{d}}^2}{\textit{z}}}}{{{\rm{d}}{t^2}}} + \mu \frac{{{\rm{d}}{\textit{z}}}}{{{\rm{d}}t}} + k{\textit{z}} = - {p_{\rm{i}}}(t)$$ (10) 式中：
$M$ 为试验装置质量，弹簧质量忽略不计；$\mu $ 为弹簧刚度系数；$k$ 为阻尼系数。根据式(5)～式(10)，可得$$\frac{{{{\rm{d}}^2}{\textit{z}}}}{{{\rm{d}}{t^2}}} + 2\beta \frac{{{\rm{d}}{\textit{z}}}}{{{\rm{d}}t}} + {\omega ^2}{\textit{z}} = - \frac{m}{M}\frac{{{{\rm{d}}^2}x}}{{{\rm{d}}{t^2}}}$$ (11) 式中，
$2\beta = \dfrac{\mu }{M},\;{\omega ^2} = \dfrac{k}{M},$ 其中$\omega $ 为弹簧系统的固有频率。 -
在非接触水下爆炸与舰船相互作用的过程中，存在众多不确定性因素，且可分为2类：一类是不确定的物理量；另一类是不确定的唯象参数（也称“拟合系数”）。表1中，
${\xi _1}\sim {\xi _6}$ 为可通过试验标定的不确定物理量； 表2中，${\xi _7}\sim {\xi _{14}}$ 为无法通过试验标定的不确定唯象参数，。符号 不确定的物理量 对数正态概率分布 ${\xi _1}$ 斜距$r$/m $[2.302,\;0.020]$ ${\xi _2}$ TNT炸药质量$W$/kg $[4.605,\;0.120]$ ${\xi _3}$ 甲板面密度$m$/(kg·m−2) $[4.379,\;0.078\;9]$ ${\xi _4}$ 试验装置质量M/kg $[7.489\;4,\;0.111]$ ${\xi _5}$ 弹簧刚度系数$\mu $/(N·mm) $[16.705,\;0.010]$ ${\xi _6}$ 阻尼系数$k$/(Ns·m−1) $[10.859,\;0.004]$ Table 1. Uncertainty of ship subjected to non-contact underwater explosion（physical quantities）
符号 不确定的唯象参数 Beta概率分布 ${\xi _7}$ $a$ $[2,\;2,\;0.82,\;0.83];0.83]$ ${\xi _8}$ $b$ $[2,\;2,\;0.16,\;0.18]$ ${\xi _9}$ ${a_1}$ $[3,\;4,\;1.2,\;1.4]$ ${\xi _{10}}$ ${a_2}$ $[3,\;4,\;0.175,\;0.185]$ ${\xi _{11}}$ ${K_1}$ $[5,\;4,\;50,\;54]$ ${\xi _{12}}$ ${K_2}$ $[5,\;4,\;0.88,\;0.94]$ ${\xi _{13}}$ $\alpha $ $[3,\;4,\;1.12,\;1.14]$ ${\xi _{14}}$ $\beta $ $[3,\;4,\;- 0.20,\;- 0.16]$ Table 2. Uncertainty of ship subjected to non-contact underwater explosion（empirical parameters）
对于可通过试验标定的
${\xi _1}\sim {\xi _6}$ 这些不确定物理量，假设它们服从对数正态分布。此假设既是基于物理量的统计结果，也是因为对数正态分布可严格保证物理量的非负性，且仅需知道期望值和标准差等低阶统计量即可确定对数正态分布LN$[\tau ,\sigma ]$ 的参数（其中$\tau $ 为对数期望值，$\sigma $ 为对数标准差）。式（12）给出了$\tau ,\sigma $ 计算公式。$$\begin{split} & \tau = \lg \left(\frac{{E{{\left( \xi \right)}^2}}}{{\sqrt {std{{\left( \xi \right)}^2} + E{{\left( \xi \right)}^2}} }}\right)\\&\sigma = \sqrt {\lg \left(\frac{{{E^2}\left( \xi \right) + std{{\left( \xi \right)}^2}}}{{E{{\left( \xi \right)}^2}}}\right)} \end{split} $$ (12) 式中：
$E\left( \xi \right)$ 为随机变量$ \xi $ 的期望值；$std\left( \xi \right)$ 为$ \xi $ 的标准差。以鱼雷等水下武器爆炸为例，由于海浪和暗流作用，此类武器难以停留在固定的位置，且其是在运动中发生爆炸，确定位置的难度更大，这导致表示斜距r的物理量
${\xi _1}$ 具有了不确定性。假设表示斜距r的物理量${\xi _1}$ 满足$E\left( {{\xi _1}} \right) = 10\;{\rm{m}}$ ，$std\left( {{\xi _1}} \right) = 0.2\;{\rm{m}}$ ，由式(12)计算可得${\mu _{{\xi _1}}} = 2.032, {\sigma _{{\xi _1}}} =0.02$ 。本文使用的球形TNT炸药在加工过程中会因存在的孔洞和间隙随机波动，且还会混入杂质颗粒，导致颗粒凝结不均匀，进而使炸药密度分布[9,10,19]呈现随机性。此外，表示TNT炸药质量W的物理量
${\xi _2}$ 满足$E\left( {{\xi _2}} \right) = 100\;{\rm{kg}}$ ，$std\left( {{\xi _2}} \right) = 2\;{\rm{kg}}$ ，由式(12)计算可得${\mu _{{\xi _2}}} = 4.605,{\sigma _{{\xi _2}}} = 0.120$ 。表示甲板面密度m的物理量${\xi _3}$ 和表示试验装置质量M的物理量${\xi _4}$ 因测量因素存在不确定性，导致测得的结果产生波动，数字统计特征满足$E\left( {{\xi _3}} \right) = 80\;{\rm{kg}}$ ，$std\left( {{\xi _3}} \right) = 6\;{\rm{kg}}$ ,$E\left( {{\xi _4}} \right) = 1\;800\;{\rm{kg}}$ ，$std\left( {{\xi _4}} \right) = 200\;{\rm{kg}}$ 。由式（12）计算分别可得${\mu _{{\xi _3}}} = 4.379,\;{\sigma _{{\xi _3}}} = 0.078\;9$ ，${\mu _{{\xi _4}}} = 7.489\;4, \; {\sigma _{{\xi _4}}} = 0.111$ 。表示刚度系数$\mu $ 的物理量${\xi _5}$ 和表示阻尼系数k的物理量${\xi _6}$ 的随机性是材料固有属性所致，数字统计特征分别满足$E\left( {{\xi _5}} \right) = 1.8 \times {10^7}\;{\rm{N/m}}$ ，$std\left( {{\xi _5}} \right) \!= \!1.8 \!\!\times\!\! {10^5}\;{\rm{N/m}}$ ，$E\left( {{\xi _6}} \right) \!\!=\!\! 52\;000\;{\rm{Ns/m}}$ ，$std\left( {{\xi _6}} \right) \!\!=\!\! 200\;{\rm{Ns/m}}$ 。由式(12)计算分别可得${\mu _{{\xi _5}}} = 16.705, \;{\sigma _{{\xi _5}}} = 0.010$ ，${\mu _{{\xi _6}}} = 10.859$ ，${\sigma _{{\xi _6}}} = 0.004$ 。如上所述，唯象参数无法通过试验标定，一般需根据工程经验在一定范围内界定，并假设这些参数服从Beta概率分布。虽然均匀分布也可限定参数的取值范围，但鉴于密度函数的强间断性，使得很难将其转化为正态分布，这也是本文选择Beta概率分布的原因所在。尤其需指出的是，Beta概率分布
$[\chi ,\;\beta ,\;a,\;b]$ 中的参数$\chi,\;\beta$ 决定了概率密度函数(probability density function, PDF)的曲线形状，而$a,\;b$ 则给出了随机变量取值区间的上、下界，且取值时有赖于工程师的经验[1,3,7]。为直观描述不确定度的统计特征，图2给出了随机变量（不确定的物理量）的PDF图像。
Figure 2. Probability density function of random variables
-
二次自适应基函数齐次Wiener混沌方法成立的前提是随机变量必须是满足独立同分布的标准正态随机变量，但由2.1节可知，此条件并未得到满足。本文使用Rosenblatt变换[26]将相关随机变量转化为服从标准正态分布的独立随机变量组。具体步骤为：设
$\{ {X_1},\;{X_2},\;\ldots ,\;{X_n}\} $ 为一列随机变量（其中下标n为随机变量的个数）。令$$ 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) $$ 其中，
$$ \begin{split} & \qquad\qquad F_{Y_{n}}\left(y_{n} | y_{1}, y_{2}, \cdots, y_{n-1}\right)=\\&P\left\{Y_{n} \leqslant y_{n} | Y_{1}=y_{1}, Y_{2}=y_{2}, \cdots, Y_{n-1}=y_{n-1}\right\} \end{split} $$ 表示条件概率，故
$$ y_{n}=F_{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) $$ 则
$\{ {Y_1},{Y_2}, \cdots ,{Y_n}\} $ 服从标准正态分布且相互独立。 -
令
$\varOmega $ 为样本空间，${\cal{F}}$ 为$\varOmega $ 上的$\sigma $ 代数，$P:{\cal{F}} \to [0,1]$ 是定义在${\cal{F}}$ 上的概率测度，Gelfand三元组$(\varOmega ,{\cal{F}},P)$ 为完备概率空间，${L^2}(\varOmega )$ 为$\varOmega $ 上的平方可积空间。 若赋予内积$< fg > = \int_\varOmega {f\left( \theta \right)g\left( \theta \right)} {\rm{d}}P(\theta ),\;{\text{而}}\;f,g \in {L^2}(\varOmega )$ (其中$\theta $ 为Lebesgue意义下的积分元；$f,g \in {L^2}(\varOmega )$ 为广义函数)，则可构成一个Hilbert空间。设
$\left\{ {{\xi _i}} \right\}_{i = 1}^d$ 为表1和表2中的所有随机变量(上标d为维度)，经Rosenblatt变换后，可转化为一列标准的独立高斯随机变量组。弹簧装置的位移满足式(1)～式(11)，其中$\xi = {[{\xi _1},\;{\xi _2},\;\cdots ,\;{\xi _d}]^{\rm{T}}}$ ，而由表1和表2可知$d = 14$ 。根据Cameron-Martin定理[21-23, 27]，有
$${\textit{z}}(t,\xi ) = \sum\limits_{{{\alpha }} \in {\cal{I}}} {{{\textit{z}}_{{\alpha }}}(t){\varPsi _{{\alpha }}}({{\xi }})} $$ (13) 式中：
${{\alpha }} = ({\alpha _1},\;{\alpha _2},\; \cdots ,\;{\alpha _d})$ ，为指标集(Lebesue空间)${\cal{I}} = {{\bf{N}}^d}$ ，其中${\bf{N}}$ 为自然数集合;${{\textit{z}}_{{\alpha }}}(t) \triangleq < {\textit{z}}(t,{{\xi }}){\varPsi _{{\alpha }}}({{\xi }}) >$ （其中${{\textit{z}}_{{\alpha }}}(t)$ 为展开系数），且$$\begin{split} & {\varPsi _{{\alpha }}} = \frac{{{H_{{\alpha }}}({{\xi }})}}{{\sqrt {\alpha !} }} \\& {H_{{\alpha }}}({{\xi }}) = \prod\limits_{i = 1}^d {{h_{{\alpha _i}}}({\xi _i})} \\& {{\alpha }}! = \prod\limits_{i = 1}^d {{\alpha _i}!} \\& {\left\| {{H_{{\alpha }}}({{\xi }})} \right\|_{{L^2}(\varOmega )}} = {{\alpha }}! \end{split} $$ (14) 式中：
${h_{{\alpha _i}}}({\xi _i})$ 为${\alpha _i}$ 次一元Hermite多项式；${H_{{\alpha }}}({{\xi }})$ 为多元Hermite多项式。令
$H = S\!\!pan\left\{ {{\xi _1},\;{\xi _2},\; \cdots ,\;{\xi _d}} \right\}$ ，设${H^{:n:}}$ 为$d$ 维$n$ 次随机Hermite多项式生成的空间。根据Wiener-Ito-Segal同构公式[28]，${L^2}(\varOmega ) = { \oplus _n}{H^{:n:}}$ 。在实际应用中，式(13)需截断有限次，即$${{\textit{z}}^p}(t,{{\xi }}) \triangleq \sum\limits_{{{\alpha }} \in {{\cal{I}}_p}} {{{\textit{z}}_{{\alpha }}}(t){\varPsi _{{\alpha }}}({{\xi }})} $$ (15) 式中：
$p \triangleq {{(n + d)!}}/({{n!d!}} )- 1$ ,为n阶多项式展开次数；${{\cal{I}}_p} = \{ |{{\alpha }}| \leqslant p||{{\alpha }}| = {\alpha _1} + {\alpha _2} + \cdots + {\alpha _d}\} $ ，且在均方意义下$\mathop {\lim }\limits_{p \to \infty } {{\textit{z}}^p} (t,\;{{\xi }}) = {\textit{z}}(t,\;{{\xi }})$ 。 -
${{A}}$ 表示Rd上的酉矩阵（其中R为Euclid空间）, 定义$${{\eta }} = {{A\xi }}$$ (16) 式中，
${{\eta }}$ 为水下爆炸不确定性因素${{\xi }}$ 的线性变换，在数学意义上，其也是$H$ 的一组基, 且${H^{:n:}}$ 也可以由${{\eta }}$ 生成。令弹簧装置的位移由下式表示为$${{\textit{z}}^A}(t,{{\eta }}) \triangleq {\textit{z}}(t,\xi ),\;{\varPsi _{{\alpha }}^A}({{\xi }}) = {\varPsi _{{\alpha }}}({{\eta }})$$ (17) 即
$${\textit{z}}(t,{{\xi }}) = \sum\limits_{{{\alpha }} \in {{\cal{I}}_p}} {{{\textit{z}}_{{\alpha }}}(t){\varPsi _{{\alpha }}}({{\xi }})} = \sum\limits_{{{\beta }} \in {{\cal{I}}_p}} {{{\textit{z}}_{{\beta }}}(t){\varPsi _{{\beta }}}({{\xi }})} $$ (18) 因此，
$$\begin{split} & {{\textit{z}}_{{\beta }}^{{A}}}(t) = \sum\limits_{{{\alpha }} \in {{\cal{I}}_P}} {{{\textit{z}}_{{\alpha }}}(t) < {\varPsi _{{\alpha }}}{\varPsi _{{\beta }}^{{A}} }> } =\\ & \qquad\sum\limits_{{{\alpha }},|{{\alpha }}| = |{{\beta }}|} {{{\textit{z}}_{{\alpha }}}(t) < {\varPsi _{{\alpha }}}{\varPsi _{{\beta }}^{{A}}} > } \end{split} $$ (19) -
令
${V_{\cal{I}}}$ =${\cal{L}}\left( {{\varPsi _\beta }} \right)$ ，${{\beta }} \in {\cal{I}}$ ，${\cal{I}} \subset {{\cal{I}}_p}$ ，则弹簧装置的位移在$V_{\cal I}$ 上的投影定义如下：$${{\textit{z}}^{{{A}},{\cal{I}}}}(t,{ \xi} ) = {{\textit{z}}^{\cal{I}}}(t,{{\eta }}) = \sum\limits_{{{\beta }} \in {\cal{I}}} {\sum\limits_{{{\alpha }} \in {{\cal{I}}_P}} {{{\textit{z}}_{{\alpha }}}(t){\varPsi _{{\beta }}}({{\eta }}) < {\varPsi _{{\alpha }}}{\varPsi _{{\beta }}^{{A}}} > } } $$ (20) 此外，
${{\textit{z}}^A}(t,{{\eta }})$ 在${V_{\cal{I}}}$ 上的投影又可表示为$${{\textit{z}}^{{{A}},{\cal{I}}}}(t,{{\eta }}) = \sum\limits_{{{\gamma }} \in {\cal{I}}} {{\textit{z}}_{{\gamma }}^{\cal{I}}(t){\varPsi _{{\gamma }}}({{\xi }})} $$ (21) 因此
$${\textit{z}}_{{\gamma }}^{\cal{I}}(t) = \sum\limits_{{{\beta }} \in {\cal{I}}} {\sum\limits_{{{\alpha }} \in {{\cal{I}}_P}} {{{\textit{z}}_{{\alpha }}}(t) < {\varPsi _{{\alpha }}}{\varPsi _{{\beta }}^{{A}}} > < {\varPsi _{{\beta }}}{\varPsi _{{\gamma }}} > } } $$ (22) 然后，将
${\textit{z}}$ 限制在${V_{\cal{I}}}$ 上，得到${\textit{z}}_{{\gamma }}^{\cal{I}} = \{ {{\textit{z}}_{{\gamma }}},{{\gamma }} \in {\cal{I}}\} $ 。 -
使用二次自适应方法选取
${{A}}$ 和${\cal{I}}$ , 令$${\textit{z}} = {{\textit{z}}_0} + \sum\limits_{i = 1}^d {{{\hat {\textit{z}}}_i}} {\xi _i} + \sum\limits_{i = 1}^d {\sum\limits_{j = 1}^d {{{\hat {\textit{z}}}_{ij}}} ({\xi _i}{\xi _j} - {\delta _{ij}})} + \sum\limits_{|{{\alpha }}| \geqslant 3} {{{\textit{z}}_{{\alpha }}}(t){\varPsi _{{\alpha }}}({{\xi }})} $$ (23) 式中：z0为位移z的期望值；
${\hat {\textit{z}}_i} = {{\textit{z}}_{{{{e}}_i}}}, \;$ ${\hat {\textit{z}}_{{{e}}_{ij}}} = {{\textit{z}}_{{_{ij}}}}/\sqrt 2$ ，其中ei为${\cal I}_p$ 的子集, 是第i个位置不为1，其余位置全为0的单位向量；${\xi _i}$ ，${\xi _j}$ 分别为表1和表2中的随机变量；${\delta _{ij}}$ 为Kronecker$\delta $ 算子。式（23）的主部可以改写为$${{\textit{z}}_{II}} = {{\textit{z}}_0} - \sum\limits_{i = 1}^d {{{\textit{z}}_{ii}}} + \sum\limits_{i = 1}^d {{{\hat {\textit{z}}}_i}} {\xi _i} + \sum\limits_{i = 1}^d {\sum\limits_{j = 1}^d {{{\hat {\textit{z}}}_{ij}}} {\xi _i}{\xi _j}} $$ (24) 式中，
${\hat {\textit{z}}_{ii}} = {{\textit{z}}_{2{{{e}}_i}}}/\sqrt 2$ 令
$${{S}} = \frac{1}{{\sqrt 2 }}\left( {\begin{array}{*{20}{c}} {{{\textit{z}}_{2{e_1}}}}&{{{\textit{z}}_{{e_1} + {e_2}}}}& \cdots &{{{\textit{z}}_{{e_1} + {e_d}}}} \\ {{{\textit{z}}_{{e_2} + {e_1}}}}&{{{\textit{z}}_{2{e_2}}}}& \cdots &{{{\textit{z}}_{{e_2} + {e_n}}}} \\ \vdots & \vdots & \vdots & \vdots \\ {{{\textit{z}}_{{e_n} + {e_1}}}}&{{{\textit{z}}_{{e_n} + {e_2}}}}& \cdots &{{{\textit{z}}_{2{e_n}}}} \end{array}} \right)$$ (25) 根据式(23), 则有
$${{AS}}{{{A}}^{\rm{T}}} = {{D}}$$ (26) 式中，
${{D}} ={\rm{diag}}\{ {d_1},{d_2}, \cdots ,{d_d}\}$ 。根据线性代数知识，${{A}},{{D}}$ 分别为${{S}}$ 的特征向量及特征值矩阵，其中${{A}}$ 为酉矩阵，则有$${{\textit{z}}^A}({{\eta }}) = {{\textit{z}}_0} + \sum\limits_{i = 1}^d {{b_i}{\eta _i}} + \sum\limits_{i = 1}^d {{d_i}(\eta _i^2 - 1)} + \sum\limits_{|{{\beta }}| \geqslant 3} {{{\textit{z}}_{{\beta }}}(t){\varPsi _{{\beta }}}({{\eta }})} $$ (27) 式中，
${b_i} = \displaystyle\sum\limits_{j = 1}^d {{a_{ij}}{{\textit{z}}_{{{{e}}_j}}}}$ 。令
${\cal I} \buildrel \Delta \over = {\cal E} = \bigcup\limits_{i = 1}^d {{\cal E}_i} $ , 则有$${{\textit{z}}^{A,{\cal{I}}}}({{\eta }}) = {{\textit{z}}_0} + \sum\limits_{i = 1}^d {\sum\limits_{{{\beta }} \in {{{e}}_i}} {{{\textit{z}}}_{{\beta }}^A(t){\varPsi _{{\beta }}}({{\eta }})} } $$ (28) ${{\textit{z}}^A}({{\eta }})$ 与${{\textit{z}}^{A,{\cal{I}}}}({{\eta }})$ 的误差为$${{\textit{z}}^A}({{\eta }}) - {{\textit{z}}^{A,{\cal{I}}}}({{\eta }}) = \sum\limits_{p \geqslant |{{\beta }}| \geqslant 3} {{{\textit{z}}}_{{\beta }}^A(t){\varPsi _{{\beta }}}({{\eta }})} $$ (29) 式(24)的系数通过非嵌入方法求解，即
$${{\textit{z}}_{{\beta }}^{{A}}} \approx \sum\limits_{r = 1}^s {{{\textit{z}}^{{{A}},{\cal{I}}}}(t,{{{\eta }}^{(r)}})} {\varPsi _{{\beta }}}({{{\eta }}^{(r)}}){w_r},\;{{\beta }} \in \varepsilon $$ (30) 式中：
${{{\eta }}^{(r)}}$ 和${w_r}$ 分别为求积点和权重，其中η(r)=(η1(r), η2(r), ···, ηd(r))；s为求积点的个数。${{\textit{z}}^{A,{\cal{I}}}}(t,{{{\eta }}^{(r)}})$ 满足式(1)～式(11), 及${{{\xi }}^{\left( r \right)}} = {{{A}}^{ - 1}}{{{\eta }}^{\left( r \right)}}$ 。 -
运用第2节所述方法，选取
$p = 2$ , 则式(24)与式（31）等价。$$ \begin{split} & {{\textit{z}}^{A,{\cal{I}}}}({{\eta }}) = {{\textit{z}}_0} + \sum\limits_{i = 1}^d {{\textit{z}}_i^A{\eta _i}} + \sum\limits_{i = 1}^d {{\textit{z}}_{ii}^A\frac{{(\eta _i^2 - 1)}}{{\sqrt 2 }}} \\& \qquad\;\; {\textit{z}} _i^A = \sum\limits_{j = 1}^d {{a_{ij}}{{\textit{z}}_{{{{e}}_j}}}} ,\;{\textit{z}}_{ii}^A = \sqrt 2 {d_i} \end{split} $$ (31) 按照第2节的方法，计算得到弹簧系统试验装置在外法线方向的相对位移
${\textit{z}}(t)$ 的期望值和标准差以及置信区间，如图3～图6所示。
Figure 3. Expectation for displacement z in spring system
Figure 4. Standard deviation for displacement z in spring system
Figure 5. Confidence interval for displacement z in spring system
Figure 6. PDF of displacement z in spring system at different times
由图3可见，100 ms时刻后
${\textit{z}}(t)$ 的期望值逐渐趋于0，但仍有轻微振荡，这表明在水下爆炸冲击波作用后，甲板一直处于振动状态。由图4可见，
${\textit{z}}(t)$ 的标准差在40 ms时刻后趋于0，这说明水下爆炸冲击载荷对系统的影响主要集中在初始时刻100 ms内；在100 ms后，爆炸冲击波的影响可忽略不计，这说明受到爆炸冲击波冲击时前10 ms内，舰船振荡最为明显，此时抗爆炸冲击的影响至关重要。$E|{\textit{z}}(t)|$ 在1.87 ms时刻达到最大值12.3 mm，而$std({\textit{z}}(t))$ 在3.86 ms时刻达到最大值2.38 mm。标准差达到峰值的时间比期望值达到峰值的时间延长，这主要是惯性所导致的结果。通过比较图3与图4发现，
${\textit{z}}(t)$ 的标准差振荡比期望值的剧烈，标准差曲线的光滑度也较差。由图5可见，位移
${\textit{z}}(t)$ 的置信区间在10 ms内急剧变宽，在8.65 ms时刻宽度最大，然后逐渐变窄，并趋于0。这说明在初始时刻不确定度较大，不易预测，而这恰恰是舰船防护最重要的时刻。同时，由图5还可见，40 ms时刻后的预测精准，此时爆炸冲击波虽还在波动，但已逐渐消失。为了获得甲板上设备发生振荡的更精准结果，图6选取了6个时刻作为观测点，计算得到了
${\textit{z}}(t)$ 的PDF曲线。由图6可见，PDF曲线大致呈钟形结构，此结果也符合认知。同时，根据上述结果也可得到${\textit{z}}(t)$ 在此时刻的其他数字统计特征，例如峰值和偏度。由图6还可见，随着时间的变化，
${\textit{z}}(t)$ 的变化范围越来越窄，PDF峰值增加，而偏度则交替出现。综上所述，采用本文方法对长时间的动力行为进行预测要比初始阶段容易一些。 -
本文通过设计合适的弹簧系统试验装置，应用概率统计方法研究了中、低强度水下爆炸的不确定性因素对甲板上弹簧系统试验装置的影响。利用基于自适应基函数的齐次Wiener混沌方法，给出了试验装置位移的期望值、标准差、置信区间以及概率密度函数。得到如下主要结论：
1） 甲板受到水下爆炸冲击波的冲击后，一直处于振荡状态。弹簧系统试验装置的期望值和标准差在达到极大值后逐渐趋于0，置信区间也逐渐变窄。标准差的振荡相比期望值大很多，且标准差达到极值的时间落后于期望值。因此，当舰船受到非接触水下武器的攻击时，爆炸冲击初始阶段的毁伤效果最大且不易预测，此时的舰船防护至关重要。
2） 基于自适应基函数的齐次Wiener混沌方法，通过构造随机基函数的同构酉变换，选取合适的投影空间，利用系统响应量在低维空间的结构逼近全系统，可在一定程度上缓解“维数灾难”问题，提高计算效率，节约计算成本，在应用上具有可行性。例如，使用标准多元多项式混沌方法，展开5次多项式，截断长度为
$\left( {{\rm{14 + 5}}} \right)!/14!5! - 1 = 46\;512 - 1 = 46\;511$ ，若使用基于自适应基函数的齐次Wiener混沌方法，截断长度为$\left( {{\rm{1 + 5}}} \right)!/1!5! - 1 = 5$ 次，效率为${\rm{46\;511/5}} \approx {10^4}$ 。本文研究方法可用于指导船舶设计人员预测甲板上物体的振荡范围，判断鱼雷毁伤影响，为舰上人员采取防护措施提供建议，给出舰船加固标准，并判断舰载设备装舰的可行性。基于自适应基函数的齐次Wiener混沌方法还可推广到其他船舶冲击响应研究中。
综上所述，舰船非接触水下爆炸的UQ研究是一个系统工程，需要海洋、工程、数学、物理等各领域专家协同合作，本文仅给出了初步结果。下一步工作拟考虑以下问题：
1） 由于本文尚缺乏实验数据与数值结果的比对，且尚未获取真实的实验数据，所以下一步拟与此领域的专家联合研究水下爆炸试验不确定度的传播和量化，将实验结果与数值结果进行比对，以确认模型的参数。
2） 本文未考虑模型不确定度的影响。事实上，炸药类型不同，峰值压力和衰减常数甚至是拟合函数公式也会不同，即使是同一类型的炸药，也可能会用不同的经验函数表示。因此，研究不同的经验函数对系统输出结果的影响，即模型形式不确定度(model form uncertainty)的量化，将始终是UQ研究的一个重要课题。
-
感谢山东科技大学公派访问学者项目对第一作者在美国南加州大学访学期间的资助。感谢南加州大学Roger Ghanem教授对本文选题提出的建议。
The uncertainty quantification of ship shock environment subjected to non-contact underwater explosion
doi: 10.19693/j.issn.1673-3185.01826
- Received Date: 2019-11-16
- Rev Recd Date: 2020-02-26
- Available Online: 2020-11-10
- Publish Date: 2020-12-30
-
Key words:
- adaptive basis function /
- uncertainty quantification /
- non-contact underwater explosion /
- unitary transformation /
- Rosenblatt transformation /
- homogeneous Wiener chaos
Abstract:
| Citation: | LIANG X, CHEN J T, WANG R L, et al. The uncertainty quantification of ship shock environment subjected to non-contact underwater explosion[J]. Chinese Journal of Ship Research, 2020, 15(6): 128–136 doi: 10.19693/j.issn.1673-3185.01826
|
DownLoad:
