计算机模拟升降法试验随机序列产生与统计检验

董海平; 蔡瑞娇; 严楠; 张天飞

doi:10.11883/1001-1455(2004)01-0049-5

计算机模拟升降法试验随机序列产生与统计检验

doi: 10.11883/1001-1455(2004)01-0049-5

计量
- 文章访问数: 2197
- HTML全文浏览量: 89
- PDF下载量: 62
- 被引次数: 0
出版历程
- 刊出日期: 2004-01-15

摘要

摘要: 阐述了计算机模拟升降法试验中随机序列产生的基本原理。从参数分布、均匀性和不相关性三个方面对该随机序列进行了统计检验 ,并以该随机序列为基础 ,通过计算机模拟升降法试验 ,求得了感度分布标准差估计值的纠偏因子表 ,实例证明该纠偏因子表是有效的。

HTML全文

随着数值计算方法的进步和计算能力的提升，数值仿真技术为解决水下爆炸机理分析、结构响应机制规律揭示等科学技术难题提供了有力的手段。在水下爆炸数值仿真研究中，研究人员已对水介质状态方程、人工黏性系数、网格密度等影响计算精度的因素开展了较为细致的探讨^[1-3]。

在进行水下爆炸问题计算时，需设定水介质的初始参数以设置流场初始压力。已有研究中，大多采用以下两种方式：一是认为水介质密度变化很小，可忽略不计，仅通过改变内能设置流场初始压力^[4-6]；二是认为流体静压力的增大导致水的密度增大，仅通过调整比容来设置初始压力^[7-8]。当采用前一种方式时，冲击波超压峰值随水深的增大略有减小^[6]，后一种方式下的变化趋势则相反^[7]。可以看出，初始参数的设置方式会带来水下爆炸载荷计算结果的差异，而对此的细致分析尚未见报道。

除直接设置水介质初始参数外，也可通过施加重力、采用程序提供的特定功能等方式实现流场压力初始化^[9-10]。这些方式简化了操作，但程序仍需提供水介质的初始参数使其具有一定的压力。以LS-DYNA提供的关键字INITIAL_EOS_ALE为例，其说明文档指出，程序将采用一种迭代方法根据用户提供的压力计算介质的初始内能和比容^[11]。简而言之，该类方法其本质上也是对水介质的初始参数进行设置。

水介质初始参数设置本质上是对水介质参数随水深变化的规律的反映。本文中将首先从常用的状态方程中选定符合参考状态时参数的水介质状态方程；然后从热力学角度分析现有两种设置方式对应的热力学过程，并给出第3种参数设置方式，同时以关键字INITIAL_EOS_ALE为例，对程序给出的初始化结果进行分析；随后采用LS-DYNA程序进行一维球形装药水下爆炸仿真研究，细致分析前三种初始参数设置方式对水下爆炸载荷特性的影响，并与已有研究成果进行对比；最后确定适当的初始参数设置方式。

1. 水介质状态方程

数值仿真软件中常见的水介质状态方程为Mie-Grüneisen状态方程和Polynomial状态方程。Mie-Grüneisen状态方程以冲击绝热线作为参考线，最终形式为：

$p = \left\{ \begin{array}{*{20}{l}} \dfrac{{{\rho _0}c_0^2\mu \left[ {1 + \left( {1 - \displaystyle\frac{{{\gamma _0}}}{2}} \right)\mu - \dfrac{a}{2}{\mu ^2}} \right]}}{{{{\left[ {1 + \mu -\displaystyle \sum\limits_{i = 1}^3 {{S_i}\dfrac{{{\mu ^i}}}{{{{\left( {1 + \mu } \right)}^{i - 1}}}}} } \right]}^2}}} + \left( {{\gamma _0} + a\mu } \right){e_V}&\mu \text{≥} {\text{0}} \\ {\rho _0}c_0^2\mu + \left( {{\gamma _0} + a\mu } \right){e_V}&\mu {{\text{＜} 0}} \end{array} \right.$

(1)

式中：c₀为参考状态下的声速；ρ₀为参考状态下的密度；μ为压缩比，μ = ρ/ρ₀−1；γ₀为Mie-Grüneisen系数；a为体积修正系数；S₁、S₂和S₃为试验拟合系数；e_V为体积内能增量。

从Mie-Grüneisen状态方程的推导过程可以看出，该方程忽略了参考压力p₀项，内能项e_V实际上是相对于参考状态下的内能增量^[12]。在参考状态下，μ = 0，e_V = 0，进而求得p = 0，这与真实参考压力p₀ = 101325 Pa不符。现有仿真计算中，大多将e_V设置为一个非零值，以使p = p₀^[8]。

Polynomial状态方程一般是对Mie-Grüneisen状态方程或者试验数据拟合得到^{[4, 8]}，常见形式为：

$p = \left\{ \begin{array}{*{20}{l}} {A_1}\mu + {A_2}{\mu ^2} + {A_3}{\mu ^3} + \left( {{B_0} + {B_1}\mu } \right){\rho _0}{e_M}&\mu \text{≥} {\text{0}} \\ {T_1}\mu + {T_2}{\mu ^2} + {B_0}{\rho _0}{e_M}&\mu \text{＜} 0 \end{array} \right.$

(2)

$p = \left\{ \begin{array}{*{20}{l}} {a_1}\mu + {a_2}{\mu ^2} + {a_3}{\mu ^3} + \left( {\displaystyle\sum\limits_{i = 0}^3 {{b_i}{u^i}} } \right){\rho _0}{e_M}&\mu \text{≥} 0 \\ {a_1}\mu + \left( {{b_0} + {b_1}\mu } \right){\rho _0}{e_M}&\mu \text{＜} 0 \end{array} \right.$

(3)

$p = \left\{ \begin{array}{*{20}{l}} {C_0} + {C_1}\mu + {C_2}{\mu ^2} + {C_3}{\mu ^3} + \left( {\displaystyle\sum\limits_{i = 4}^6 {{C_i}{\mu ^{i - 4}}} } \right){e_V}&\mu \text{≥} 0 \\ {C_0} + {C_1}\mu + {C_3}{\mu ^3} + \left( {{C_4} + {C_5}\mu } \right){e_V}&\mu \text{＜} 0 \hfill \\ \end{array} \right.$

(4)

式中：A₁、A₂、A₃、B₀、B₁、T₁和T₂均为状态方程系数，e_M为质量内能增量，a₁、a₂、a₃、b₀、b₁、b₂和b₃均为状态方程系数，C₀、C₁、C₂、C₃、C₄、C₅和C₆均为状态方程系数。

与Mie-Grüneisen状态方程类似，式（2）～（3）也存在与实际压力不符的问题。式（4）在C₀ = p₀时，可满足μ = 0，e_V = 0以及p = p₀。因此，宜选择Polynomial状态方程式（4）作为水介质状态方程，具体参数见表1^[13]。

表 1 水介质状态方程参数

Table 1. EOS parameters of water

C₀/Pa	C₁/GPa	C₂/GPa	C₃/GPa	C₄	C₅	C₆
101325	2.2	9.54	14.57	0.28	0.28	0

下载: 导出CSV

| 显示表格

2. 初始参数设置

2.1 设置方式及状态参数变化

根据水介质状态方程，设置初始压力时可改变参量μ和e_V。现今大多数仿真研究是通过改变e_V（方式Ⅰ）来进行流场初始压力设置^[4-6]。从热学角度看，方式Ⅰ中不同水深下的水介质参数按等容过程（dV=0）变化，此时外界做功为零，流场的内能增量全部由外界传导的热量提供（dE_i = dQ）。这表明采用方式Ⅰ时，流体静压力的增大单纯源于外界传热，与重力等外力做功等因素无关，这与实际深水环境以及加压模拟深水环境严重不符。

少数研究人员通过调整μ（方式Ⅱ）来进行初始压力设置^[7-8]。此时，内能变化为零（dE_i = 0），即等内能形式，外界对流体做功与外界向流体传导的热量之和为零（dE_i = dQ − pdV = 0），水介质被压缩后向外放出热量。为分析温度变化趋势，引入含温度T的内能表达式^[14]：

$\dfrac{{{e_V}}}{{{p_0}}} = 10\;160\dfrac{{\Delta T}}{{{T_0}}} + 0.63 \times {10^5} \times \dfrac{{\left( {1 - {{{\rho _0}}}/{\rho }} \right)\left( {0.71 - {{{\rho _0}}}/{\rho }} \right)}}{{{{\left( {{{{\rho _0}}}/{\rho }} \right)}^{{4 /3}}}}}\left( {1 - 2\dfrac{\rho }{{{\rho _0}}}{{\rm{e}}^{ - {{\left( {{\rho }/{{{\rho _0}}}} \right)}^2}}}} \right)$

(5)

式中：T₀ = 273 K。当0 ≤ μ ≤ 0.025时，T随μ的增大而逐渐增大（见图1）。在实际深水环境中，随着水深的增大，水温逐渐减小^[15]；在通过加压构造的模拟深水环境中，水温与环境温度基本相同，即温度可认为保持不变。因此，从温度变化趋势的角度来看，方式Ⅱ与实际情况有一定的差距。

图 1 不同设置方式下的μ、e_V和ΔT值

Figure 1. Values of μ, e_V and ΔT in different setting modes

下载: 全尺寸图片幻灯片

热力学研究结果表明，温度到处相同是重力场的热平衡条件^[16]。据此，本文中提出第3种初始参数设置方式，即假设水介质参数随水深按照等温过程（dT = 0，方式Ⅲ）变化。对于该方式，ΔT = 0，将式（5）代入式（4），可得等温过程对应的状态方程。显然，对于实际深水环境，方式Ⅲ只是一种粗糙的近似，不过理论上优于方式Ⅰ和方式Ⅱ；对于模拟深水爆炸试验，近似程度较高。

除上述设置方式外，数值仿真软件往往也具备根据重力和水深信息或者设定的流体静压力计算水介质初始参数的功能，只是用户往往无法直接获知其遵循的规律。现以关键字INITIAL_EOS_ALE为例，通过读取不同水深下的流场压力初始化结果来初步了解LS-DYNA程序中μ和e_V随水深的变化规律（方式Ⅳ），结果如图1所示，可以看出，方式Ⅳ对应的热力学过程接近方式Ⅱ。

图1中给出了4种设置方式在不同水深条件时的压缩比μ、内能增量e_V和温度变化ΔT的计算结果。方式Ⅱ和方式Ⅳ的计算结果几乎重合，二者均与方式Ⅲ较接近；方式Ⅰ的计算结果与其他3种相去甚远。采用方式Ⅱ、Ⅲ、Ⅳ时，μ随水深增加而增大；采用方式Ⅰ时μ则始终为零；在同一水深处，μ的值从小到大依次为：方式Ⅰ、方式Ⅳ、方式Ⅱ、方式Ⅲ。采用方式Ⅰ和方式Ⅳ时，e_V随水深增加而增大，采用方式Ⅱ时e_V始终为零，采用方式Ⅲ时e_V随水深的增加而减小；在同一水深处，e_V的值从小到大依次为：方式Ⅲ、方式Ⅱ、方式Ⅳ、方式Ⅰ。采用方式Ⅰ、方式Ⅱ和方式Ⅳ时，ΔT随水深增加而增大；采用方式Ⅲ时ΔT始终为零；在同一水深处，ΔT的值从小到大依次为：方式Ⅲ、方式Ⅱ、方式Ⅳ、方式Ⅰ。

2.2 流体可压缩性

流体的可压缩性决定流体内微弱扰动波的传播速度，即流体内的声音传播速度，因此常用声速c来表示流体的可压缩性。声速c的表达式为：

${c^2} = {\left( {\frac{{{\text{d}}p}}{{{\text{d}}\rho }}} \right)_S}$

(6)

式中：下标S表示等熵过程。

当初始状态为( μ_∞ , e_V∞ )时，等熵压缩状态下的压力可表示为：

$p = {C_0} + {C_1}\mu + {C_2}{\mu ^2} + {C_3}{\mu ^3} + \left( {{C_4} + {C_5}\mu + {C_6}{\mu ^2}} \right)\left( {{e_{V\infty }} + \int\nolimits_{{\mu _\infty }}^\mu {\frac{p}{{{{\left( {1 + \mu } \right)}^2}}}{\text{d}}\mu } } \right)$

(7)

于是，声速可表示为：

${c^2} = \frac{{{\text{d}}p}}{{{\text{d}}\rho }} = \frac{{{\text{d}}p}}{{{\text{d}}\mu }}\frac{{{\text{d}}\mu }}{{{\text{d}}\rho }} = \frac{{{C_1} + 2{C_2}\mu + 3{C_3}{\mu ^2}}}{{{\rho _0}}} + \frac{{{C_5} + 2{C_6}\mu }}{{{\rho _0}}}\left( {{e_{V\infty }} + \int\nolimits_{{\mu _\infty }}^\mu {\frac{p}{{{{\left( {1 + \mu } \right)}^2}}}{\text{d}}\mu } } \right) + \frac{{{C_4} + {C_5}\mu + {C_6}{\mu ^2}}}{{{\rho _{0}}}}\frac{p}{{{{\left( {1 + \mu } \right)}^2}}}$

(8)

当μ = μ_∞时，声速为：

${c^2} = \frac{{{C_1} + 2{C_2}{\mu _\infty } + 3{C_3}\mu _\infty ^2}}{{{\rho _0}}} + \frac{{{C_5} + 2{C_6}{\mu _\infty }}}{{{\rho _0}}}{e_{V\infty }} + \frac{{{C_4} + {C_5}{\mu _\infty } + {C_6}\mu _\infty ^2}}{{{\rho _0}}}\frac{{{p_\infty }}}{{{{\left( {1 + {\mu _\infty }} \right)}^2}}}$

(9)

图2中给出了4种设置方式在不同水深条件时的声速。从图2中可以看出，随着水深的增大，4种方式下的水介质声速逐渐增大，即可压缩性均出现减弱，但减弱幅度不同，由小到大的排序为：方式Ⅰ、方式Ⅳ、方式Ⅱ、方式Ⅲ。可见，声速的变化与压缩比的变化情形一致。

图 2 不同设置方式下的声速变化

Figure 2. Acoustic velocity in different modes

下载: 全尺寸图片幻灯片

3. 初始参数设置对水下爆炸载荷的影响

3.1 数值仿真模型

建立了一维球对称水下爆炸数值仿真模型，并应用LS-DYNA数值仿真软件进行计算。其中，炸药采用MAT_HIGH_EXPLOSIVE_BURN材料模型以及JWL状态方程，参数取自文献[17]。水介质的状态方程采用式（4），参数在表1中列出。

炸药选用1 kg TNT，水深范围0 ～ 5 km，由于方式Ⅱ和方式Ⅳ对应的初始参数极为接近，因此只对比前3种方式，各自对应的μ和e_V值取自图1。为保证计算精度，网格尺寸为装药半径的1/100；为减小边界反射冲击波对气泡脉动的影响，水域半径取237.2 m，大于水深为0 m时声波在一个脉动周期内传播距离的1/2。

3.2 冲击波载荷特性

选取水深为5 km，爆距R = 55R₀（R₀为装药半径），计算得到了3种方式下的冲击波压力-时间曲线，如图3所示，对应的冲击波到达时间t_a、超压峰值ΔP_m和正压持续时间t_c列于表2中。可以看出，3条曲线并不完全重合，除方式Ⅰ外，其他2种方式的仿真结果相差较小；方式Ⅲ、Ⅱ、Ⅰ的t_a及t_c值依次减小，而ΔP_m依次增大。该规律与流体可压缩性结果基本一致，即流体可压缩性越弱，冲击波到达时间越短，冲击波超压峰值越大，正压持续时间越短。

图 3 H=5 km及R=55R₀时3种方式下的冲击波压力-时间曲线

Figure 3. Shock wave pressure-time curves in three modes when H=5 km and R=55R₀

下载: 全尺寸图片幻灯片

表 2

$H=5\;{\rm {km}}$ 及

$R=55R_0$ 时3种方式下的

$t_{\rm a}$ 、

$\Delta P_{\rm m}$ 及

$t_{\rm c}$

Table 2. Values of

$t_{\rm a}$ ,

$\Delta P_{\rm m}$ and

$t_{\rm c}$ in three modes when

$H=5\;{\rm {km}}$ and

$R=55R_0$

方式	t_a/ms	ΔP_m/MPa	t_c/ms
Ⅰ	1.7745	12.080	0.2676
Ⅱ	1.6578	13.603	0.2587
Ⅲ	1.6518	13.692	0.2581

下载: 导出CSV

| 显示表格

图4中给出了R = 55R₀时3种设置方式在不同水深下获得的冲击波超压峰值、冲量I和能流密度e_s，其中，I和e_s分别由下式给出；

图 4 R=55R₀时不同水深下3种方式所对应的冲击波载荷参数

Figure 4. Values of shock wave load in three modes when R=55R₀

下载: 全尺寸图片幻灯片

$I = \int\nolimits_0^{{t_{\text{c}}}} {\left( {p - {p_\infty }} \right){\text{d}}t}$

(10)

${e_{\text{s}}} = \frac{1}{{{\rho _\infty }{c_\infty }}}\int\nolimits_0^{{t_{\text{c}}}} {{{\left( {p - {p_\infty }} \right)}^2}{\text{d}}t}$

(11)

式中：e_s为冲击波能流密度，ρ_∞、c_∞分别为未扰动流体的密度和声速，为便于比较，统一取ρ_∞ = ρ₀=1000 kg/m³，c_∞ =c₀ = 1483.2 m/s。

表3中给出了冲击波载荷参数在5 km水深处相对于0 m的变化幅度。可以看出，随着水深的增大，采用方式Ⅰ时ΔP_m逐渐减小，其他2种方式则逐渐增大，且变化幅度略大于方式Ⅰ。对于I和e_s，3种方式给出的仿真结果均随水深的增大而减小，方式Ⅱ和Ⅲ的变化幅度略小于方式Ⅰ。

表 3 冲击波载荷在5 km水深处相对于0 m的变化幅度

Table 3. Changing amplitudes of shock wave load when depth changes from 0 m to 5 km

方式	变化幅度/%
方式	ΔP_m	I	e_s
Ⅰ	−1.202	−74.304	−32.896
Ⅱ	11.255	−73.227	−21.877
Ⅲ	11.981	−73.171	−21.247

下载: 导出CSV

| 显示表格

3.3 气泡脉动特性

分析H=5 km时3种设置方式下气泡的脉动曲线，如图5所示，相应的气泡最大半径R_max和第1次脉动周期T列于表4。可以看出，除方式Ⅰ外，其余两条曲线基本重合；方式Ⅲ、方式Ⅱ、方式Ⅰ对应的R_max和T依次增大。显然，该规律与水介质压缩比的变化一致。由于气泡脉动与水介质的惯性密切相关，当体积不变时，水介质密度的增大导致其惯性增大，在密度变化较小时惯性增幅有限，因此，3种方式下气泡脉动曲线的差别并不明显。

图 5 H=5 km时3种方式下气泡的半径-时间曲线

Figure 5. Bubble radius-time curves in three modes when H=5 km

下载: 全尺寸图片幻灯片

表 4

$H=5\;{\rm{km}}$ 时3种方式下对应的

$R_{\max}$ 和

$T$

Table 4. Values of

$R_{\max}$ and

$T$ in three modes when

$H=5\;{\rm{km}}$

方式	R_max/mm	T/ms
Ⅰ	178.722	1.7129
Ⅱ	180.618	1.7445
Ⅲ	180.731	1.7460

下载: 导出CSV

| 显示表格

4. 对比分析

4.1 冲击波载荷对比分析

自20世纪60年代起，研究人员开展了较多的水深对冲击波载荷特性影响的研究。Baum等^[18]在模拟深水爆炸容器中进行了深水爆炸试验，当1 < R/R₀ < 7时，ΔP_m在0 ～ 4 km的模拟水深范围内基本保持不变。Vanzant等^[19]给出了14.75 g Pentolite炸药在更大爆距范围内的试验结果，本文对其试验数据进行了线性拟合，结果表明，ΔP_m随水深的增大而略有增大，且在较大爆距处的变化更明显，如图6所示。钟帅等^[20]开展了0～150 m模拟水深水下爆炸试验，指出ΔP_m也随水深增大而略有增加。Slifko^[21]、Xiao等^[22]也进行了类似研究，不过试验在大洋中进行，爆源和测点并不位于同一深度，而且爆距较大，与本文分析工况不符。

图 6 Vanzant等^[19]的试验结果及线性拟合曲线

Figure 6. Test results from Vanzant et al^[19] and linear fitting results

下载: 全尺寸图片幻灯片

理论上，Baum等^[18]从相似性和量纲理论出发，在浅水经验公式的基础上，给出了考虑流体静压力时ΔP_m的计算式：

$\Delta {P_{\text{m}}} = A{\left( {\frac{{{W^{{1 \mathord{\left/ {\vphantom {1 3}} \right. } 3}}}}}{R}} \right)^\alpha }{\left( {\frac{{{p_\infty }}}{B} + 1} \right)^{{{\left( {N - \alpha } \right)} \mathord{\left/ {\vphantom {{\left({N - \alpha }\right)} N}} \right. } N}}}$

(12)

式中：A和α为常数，当6 ≤ R/R₀ < 12时，A = 44.1 MPa，α = 1.5；当12 ≤ R/R₀ < 240时，A = 52.4 MPa，α = 1.13；B = 304.9 MPa为Tait状态方程中的常数；N为几何指数，对于球面冲击波，N = 3。

式（12）表明，ΔP_m随水深的增大而缓慢增大；同时随着爆距的增大，α逐渐减小，流体静压力对ΔP_m的影响逐渐增大。

从前述研究中可以看到，随着水深的增大，采用方式Ⅱ和方式Ⅲ时ΔP_m会逐渐增大。从图7中可以看到，5 km水深条件下当7R₀ ≤ R≤ 65R₀时，随着爆距的增大，水深对ΔP_m的影响也在逐渐增大。图8为H = 5 km时3种方式下ΔP_m计算结果相对于式（12）的误差。综合来看，采用方式Ⅱ和方式Ⅲ获得的计算结果与试验值及理论值均较为吻合，方式Ⅲ的相对误差更小。

图 7 ΔP_m在5 km水深处相对于在0 m处的变化幅度

Figure 7. Changing amplitudes of ΔP_m when depth changes from 0 m to 5 km

下载: 全尺寸图片幻灯片

图 8 H=5 km时ΔP_m相对于计算式(8)的误差

Figure 8. Relative errors of ΔP_m between simulation results and calculation formulas when H=5 km

下载: 全尺寸图片幻灯片

4.2 气泡脉动参数对比分析

Cole^[23]基于不可压缩理论模型推导了R_max和T与水深和装药量的关系，Swift等^[24]基于试验数据给出了计算气泡脉动参数的经验公式：

${R_{{\text{max}}}} = 3.36\frac{{{W^{{1 \mathord{\left/ {\vphantom {1 3}} \right. } 3}}}}}{{{{\left( {H + 10.3} \right)}^{{1 \mathord{\left/ {\vphantom {1 3}} \right. } 3}}}}}$

(13)

$T = 2.11\frac{{{W^{{1 \mathord{\left/ {\vphantom {1 3}} \right. } 3}}}}}{{{{\left( {H + 10.3} \right)}^{{5 \mathord{\left/ {\vphantom {5 6}} \right. } 6}}}}}$

(14)

Baum等^[18]、Slifko^[21]和梁浩哲等^[25]借助试验进一步研究了气泡脉动参数随水深的变化规律，其中，Baum等^[18]验证了R_max ∝ p_∞^−1/3及T ∝ p_∞^−5/6在大水深（H > 1 km）下的适用性。

图9中给出了3种方式下气泡脉动参数仿真结果相对于经验值的误差。可以看出，在同一水深下，方式Ⅱ、Ⅲ误差较小，与经验值更接近，且方式Ⅲ的相对误差最小。

图 9 气泡脉动参数相对于经验值的误差

Figure 9. Relative errors of R_max and T between simulation and empirical formula

下载: 全尺寸图片幻灯片

4.3 设置方式合理性探讨

从热力学角度看，2.1节中已指出，采用方式Ⅰ时，流体静压力的升高单纯由外界传热引起，这与真实情形不符；采用方式Ⅱ和方式Ⅳ时，温度随水深的增大而略有增大，与实际情况有一定差距。采用方式Ⅲ时，水介质温度不随水深变化，与模拟深水试验条件吻合，与实际深水环境的吻合性较方式Ⅰ、Ⅱ和Ⅳ要好。从与已有成果的对比分析看，方式Ⅱ、Ⅲ的结果与已有成果吻合较好；相较而言，采用方式Ⅲ时载荷值的相对误差最小。综上分析认为，初始参数设置宜选用方式Ⅲ，方式Ⅱ和方式Ⅳ次之，在讨论水深影响时不建议选用方式Ⅰ。