CN112446126B - 超空泡航行体尾拍运动状态的模拟方法 - Google Patents

Abstract

本发明公开了一种超空泡航行体尾拍运动状态的模拟方法，为了给出有普遍意义的超空泡航行尾拍运动状态的模拟方法，基于超空泡航行体的纵平面运动方程和正态分布下的3σ原则，在航行体尾部空泡截面上进行建模，为研究航行体尾拍运动提供更为普遍的方法，旨在通过理论计算的手段预报超空泡航行体的尾拍状态，为分析航行体尾拍运动提供依据。

Description

超空泡航行体尾拍运动状态的模拟方法

技术领域

本发明涉及一种超空泡航行体尾拍运动仿真的计算方法，特别是一种超空泡航行体尾部空泡截面上的尾拍状态的模拟方法。

背景技术

超空泡航行体在水下运动时，由空气和水蒸汽所组成的空泡会从空化器延伸至整个航行体表面，阻隔航行体与水的接触，摩擦阻力大大减小。水下超空泡航行体在一定距离内可实现鱼雷侵彻拦截、毁伤舰船、打击蛙人等，具有十分重要的应用。

当超空泡航行体在水下运动时，头部的空化器会诱导产生流体空化，所产生的空化泡，将超空泡航行体完全包裹起来，此时仅航行体头部受力，航行体受力不均衡，航行体尾部会连续撞击空泡壁面，由此来达到自身平衡的运动称为尾拍运动。尾拍运动是航行体在水下运动时其在水下运动时的主要方式。因其运动方式复杂，且伴随有强烈的流固耦合效应，使超空泡航行体水下运动的研究简化为二维纵平面上的运动。导致模型在很大程度上的不精确，研究不具备普遍意义，目前尚没有公开的三维情况下研究超空泡航行体尾拍运动的例子。

发明内容

本发明的目的是提供一种超空泡航行体尾拍运动的模拟方法，旨在通过理论计算的手段预报超空泡航行体的尾拍状态，为分析航行体尾拍运动提供依据。

实现本发明的技术解决方案为：超空泡航行体尾拍运动的模拟方法，包括以下步骤：

步骤1、建立3σ正态分布模型：

其中f(x)为正态分布的概率密度函数，σ为方差，μ为前一次的超空泡射弹尾拍角度，k为前一次的超空泡射弹尾拍角度；

步骤3、建立σ-κ线性半经验估计方程：基于σ近似解，给出线性半经验的σ估计方程，定义无量纲数κ：

得到σ-κ线性半经验估计方程：

其中κ为每一个次超空泡航行体尾拍位置处的更新值，κ₀为计算初值；

步骤4、确定初始尾拍位置：初始位置的选择通过生成在[0，360]上的随机数确定，服从均匀分布；

步骤5：确定初始尾拍角速度：超空泡航行体的初始尾拍角速度，受不确定性的扰动而产生，因此超空泡航行体初始尾拍角速度的给定是任意的；

步骤6、确定尾部空泡截面：超空泡航行体尾部空泡面积的计算公式

其中p₀为大气压强，ρ为水的密度，g重力加速度，h航行体航行深度，S₀为空化器截面积，t为航行时间，τ为航行体航行至在计算空泡截面处的时刻，

为空泡的初始扩张速度，v₀为航行体初速度，C_x0＝0.82为空化器阻力系数，m为航行体质量，L为航行体长度，k’为系数微弱依赖于空化数σ＝Δp/(0.5ρv²)，取k’＝π。

步骤7：确定超空泡航行体尾拍位置：由初始位置的均匀分布所确定的初始超空泡航行体的尾拍位置，来计算下一次尾拍的正态分布的期望值μ_i＝k_i+180(i＝1,2,3…n)，其中k_i相对于下一次超空泡航行体的尾拍位置计算为确定值，而μ_i为对应的超空泡航行体的尾拍期望。

对每一次的μ_i加以限定，限定后记作μ_i*，限定计算方法为：

其中n为左取整函数，表达式如下：

建立极坐标系，则每一次超空泡航行体的尾拍位置由μ_i*来计算：

当在作第二时间步计算时，由于期望为μ_i的正态分布结果，存在负值的可能性，通过超空泡航行体尾部空泡截面的几何分析结果，将负值角度通过正定限定在[0,360]之间，使计算具备角度的物理意义。具体实现方法：

1:当|μ_i+1|>360计算2,否则计算7；

2:n＝[|μ_i+1|]；

3:ξ_i＝|μ_i+1|-360n；

4:μ_i+1＝ξ_i sign(μ_i+1)；

5:如果sign(μ_i+1)＝-1则μ_i+1＝360-ξ_i,否则计算6；

6:μ_i+1＝ξ_i；

7:如果sign(μ_i+1)＝-1则μ_i+1＝360+ξ_i,否则计算8；

8:μ_i+1＝μ_i+1

步骤8：确定超空泡航行体两次尾拍间隔：尾拍前后角速度的变化关系式为：

Q⁺＝Q^-[1-G(1+e)]

其中：

Q^-为超空泡航行体尾拍前的超空泡航行体绕空化器的旋转角速度，Q⁺为超空泡航行体尾拍后的超空泡航行体旋转绕空化器的旋转角速度，x_cm为超空泡航行体质心到尾部的长度。

e为补偿系数，其表达式为：

其中U^-为超空泡航行体的速度，α为超空泡航行体尾拍角，即航行体尾拍时轴线与空泡壁面的夹角:

设前一次尾拍位置坐标为(x₁,y₁)，后一次尾拍位置为(x₂,y₂)，则两次尾拍的平面距离为：

超空泡航行体尾部的线速度为：

v_L＝|Q⁺|L

由此得两次尾拍的时间间隔为；

本发明与现有技术相比，其显著优点为：1)本发明将超空泡航行体尾拍运动作为一种随机运动，仿真模拟结果更具有普遍性2)本发明得到的航行体尾拍结果超出了试验可观测范围具有弥补试验水平不足的优点，减少了试验研究的成本3)本发明计算所使用的计算方法充分考虑了影响超空泡航行体尾拍运动的因素，不受条件限制。

附图说明

图1为本发明方法流程图。

图2为超空泡弹丸尾部空泡截面图。

图3为正态分布模型。

图4为改进3σ-正态分布模型。

图5为超空泡弹丸尾部空泡刻度等分图。

图6(a)为300m/s初速度下的尾拍状态1

图6(b)为300m/s初速度下的尾拍状态2

图6(c)为300m/s初速度下的尾拍状态3

具体实施方式

下面结合附图和具体实施例对本发明方案作进一步说明。

为达到本发明目的，现给出如下假设：

(1)初始扰动对超空泡航行体尾拍位置的影响在航行体尾部空泡截面内服从均匀分布；

(2)忽略航行体重力的影响，两次尾拍之间航行体尾部运动轨迹为直线；

(3)超空泡航行体的尾拍位置k是服从参数为μ，σ²的正态分布，记作k～N(μ_i,σ²)，连续尾拍过程下μ_i+1＝k+180,σ＝f(ρ,R)，μ_i为前一次尾拍位置的角度值，ρ为空泡的曲率半径，R为航行体尾部截面的曲率半径，其中标准差σ与空泡的曲率半径成正比，与航行体尾部截面的曲率半径成反比。

(4)超空泡航行体尾拍导致的空泡变形被忽略，航行体尾部空泡截面上的航行体截面积作为点处理。

步骤1、建立3σ正态分布模型：基于正态分布自身概率分布两侧递减特征，结合超空泡航行体尾拍运动方向依照在空泡内部角度分布概率递减的规律，这里假设超空泡航行体的尾拍运动撞击方向随机变量θ服从正态分布，即θ～N(μ,σ²)，如图1，其中μ为均值，σ²为方差。由于正态分布概率密度函数自变量x的取值范围为实数域R，而在超空泡航行体尾拍运动的问题中，尾拍方向随机变量θ的选择范围为[0°,360°]，如附图5所示对应的正态分布概率密度函数自变量x取值范围为[0，360]。因此传统的正态分布模型(图3)在这里是不适用的，需采用与实际问题相匹配的3σ-正态分布模型(图4)。

根据正态分布的3σ原则，若θ～N(μ,σ²)，则θ落在(μ-3σ，μ+3σ)之间的概率为0.9973，而落在该区间以外的概率仅有0.0027，即θ几乎不可能落在(μ-3σ，μ+3σ)之外。如图2，将超空泡航行体尾部空泡截面用角度从[0°,360°]进行等分，即给定了3σ区间，不考虑该区间以外的情形，将正态分布概率密度函数自变量x应用在[0，360]上。

在区间[k,k+360]上满足正态分布概率密度函数积分值大于0.9973。其数学表达式为：

步骤2、化简3σ正态分布模型：

对3σ正态分布模型不等式进行换元，令t＝(x-μ)/σ得：

令m＝2^0.5t/2得:

由于正态分布均值μ位于区间上[k,k+360]，因此有不等式k<μ<k+360成立，且满足：

将μ的表达式带入上式并作如下处理：

化简可得标准差σ在数学意义上的取值范围为：

0＜σ≤59.84

因此得到3σ-正态分布模型为：

步骤3、建立σ-κ线性半经验估计方程：基于σ近似解，给出线性半经验的σ估计方程，根据假设(3)定义无量纲数κ：

当航行体尾部空泡截面积与航行体尾部截面相等时κ＝1，此时超空泡航行体尾拍不会发生，等价于标准差为零的正态分布，则此时标准差与无量纲数κ为一一对应的关系，初始超空泡航行体尾拍动能最大对应的航行体尾部空泡截面积最大，应对应标准差的最大值59.84，此时对应的无量纲数为κ₀，由标准差σ和无量纲数κ组成两组对应关系(0，1)，(59.84，κ₀)，得到σ-κ线性半经验估计方程：

其中κ为每一个次超空泡航行体尾拍位置处的更新值；

步骤4、确定初始尾拍位置：因尾拍初始位置受多因素影响，不确定性往往较强，初始尾拍时航行体尾部空泡由云状空化和层状空化转变为超空化过程中，航行体尾部空化面会经历一个由不完整到完整的过程。尾拍起始位置出现在尾部空泡截面各个方向的概率是一致的。因此初始位置的选择通过生成在[0，360]上的随机数确定，服从均匀分布；

步骤5：确定初始尾拍角速度：超空泡航行体的初始尾拍角速度，受不确定性的扰动而产生，因此超空泡航行体初始尾拍角速度的给定是任意的；

步骤6、确定尾部空泡截面

鉴于Logvinovich空泡扩张方程中压强所用为无穷远处压强，并非空泡附近压强，则有：p＝p₀+ρgh；

设τ为截面x的形成时刻，对Logvinovich空泡扩张方程积分可以得到：

设t_m＝t，n≤i≤m，对(29)式使用复化Simpson公式进行离散，令f＝(t-u)(p-p_v)则：

将式(28)带入式(30)可得超空泡航行体尾部空泡面积的计算公式

其中p₀为大气压强，ρ为水的密度，g重力加速度，h航行体航行深度，S₀为空化器截面积，t为航行时间，τ为航行体航行至在计算空泡截面处的时刻，

为空泡的初始扩张速度，v₀为航行体初速度，C_x0＝0.82为空化器阻力系数，m为航行体质量，L为航行体长度，k为系数微弱依赖于空化数σ＝Δp/(0.5ρv²)，取k＝π。

步骤7：确定超空泡航行体尾拍位置：由初始位置的均匀分布所确定的初始超空泡航行体的尾拍位置，来计算下一次尾拍的正态分布的期望值μ_i＝k_i+180(i＝1,2,3…n)，其中k_i相对于下一次超空泡航行体的尾拍位置计算为确定值，而μ_i为对应的超空泡航行体的尾拍期望。

由于k_i和μ_i的取值范围均为[0,360]，因此当在μ_i-1上加180以后，取值会超出该范围，会失去实际的超空泡航行体的尾拍位置的物理意义，因此需要对每一次的μ_i加以限定，限定后记作μ_i*，限定计算方法为：

其中n为左取整函数，表达式如下：

以航行体尾部空泡截面中心为极点O空泡半径为极轴r建立极坐标系，则每一次超空泡航行体的尾拍位置由μ_i*来计算：

当在作第二时间步计算时，由于期望为μ_i的正态分布结果，存在负值的可能性，通过超空泡航行体尾部空泡截面的几何分析结果，将负值角度通过正定限定在[0,360]之间，使计算具备角度的物理意义。具体实现方法：

1:当|μ_i+1|>360计算2,否则计算7；

2:n＝[|μ_i+1|]；

3:ξ_i＝|μ_i+1|-360n；

4:μ_i+1＝ξ_i sign(μ_i+1)；

5:如果sign(μ_i+1)＝-1则μ_i+1＝360-ξ_i,否则计算6；

6:μ_i+1＝ξ_i；

7:如果sign(μ_i+1)＝-1则μ_i+1＝360+ξ_i,否则计算8；

8:μ_i+1＝μ_i+1

步骤8：确定超空泡航行体两次尾拍间隔：尾拍前后角速度的变化关系式为：

Q⁺＝Q^-[1-G(1+e)]

其中：

Q^-为超空泡航行体尾拍前的超空泡航行体绕空化器的旋转角速度，Q⁺为超空泡航行体尾拍后的超空泡航行体旋转绕空化器的旋转角速度，x_cm为超空泡航行体质心到尾部的长度。

e为补偿系数，其表达式为：

其中U^-为超空泡航行体的速度，α为超空泡航行体尾拍角，即航行体尾拍时轴线与空泡壁面的夹角:

设前一次尾拍位置坐标为(x₁,y₁)，后一次尾拍位置为(x₂,y₂)，则两次尾拍的平面距离为：

超空泡航行体尾部的线速度为：

v_L＝|Q⁺|L

由此得两次尾拍的时间间隔为；

实施例1

下面结合实施例进行更详细的描述。

给定以下初始条件：

a.从水下发射超空泡弹丸，弹丸运动的初速度v₀为200m/s；

b.超空泡弹丸的质量m为0.1kg；

c.超空泡弹丸的长度L为156.81mm；

d.流体介质为水，密度ρ为1000kg/m³；

e.超空泡弹丸最大截面直径D为12.67mm；

f.超空泡弹丸的空化器直径d为4.49mm；

g.重力加速度g为9.8m/s²；

h.大气压强P₀为101325Pa；

i.弹丸飞行深度h为2.3m；

j.对应24℃下饱和蒸气压强P_v为3025Pa；

k.平面上弹丸绕质心的转动惯量为1.46×10^-4N/m；

l.弹丸质心与弹丸尾部距离x_cm为80mm；

m.时间步长Δt为10^-5s；

n.总计算时间t为0.5s。

计算结果如附图6(a)、(b)、(c)。

Claims (3)

1.一种超空泡航行体尾拍运动状态的模拟方法，其特征在于，包括以下步骤：

步骤1、建立3σ正态分布模型：

其中f(x)为正态分布的概率密度函数，σ为方差，μ为前一次的超空泡射弹尾拍角度，k为前一次的超空泡射弹尾拍角度；

步骤2、建立σ-κ线性半经验估计方程：基于σ近似解，给出线性半经验的σ估计方程，定义无量纲数κ：

得到σ-κ线性半经验估计方程：

其中κ为每一个次超空泡航行体尾拍位置处的更新值，κ₀为计算初值；

步骤3、确定超空泡航行体的初始尾拍位置：初始位置的选择通过生成在[0，360]上的随机数确定，服从均匀分布；

步骤4：确定超空泡航行体的初始尾拍角速度：

步骤5、确定尾部空泡截面：超空泡航行体尾部空泡面积的计算公式

其中p₀为大气压强，ρ为水的密度，g重力加速度，h航行体航行深度，S₀为空化器截面积，t为航行时间，τ为航行体航行至在计算空泡截面处的时刻，

为空泡的初始扩张速度，v₀为航行体初速度，C_x0＝0.82为空化器阻力系数，m为航行体质量，L为航行体长度，k’为系数微弱依赖于空化数σ＝Δp/(0.5ρv²)，取k’＝π；

步骤6：确定超空泡航行体尾拍位置：由初始位置的均匀分布所确定的初始超空泡航行体的尾拍位置，来计算下一次尾拍的正态分布的期望值μ_i＝k_i+180i＝1,2,3…n，其中k_i是当前超空泡航行体的尾拍位置角度值，μ_i为对应的超空泡航行体的尾拍期望，当出现负期望值时，需要对期望进行正定；

对每一次的μ_i加以限定，限定后记作μ_i*，限定计算方法为：

其中n为左取整函数，表达式如下：

建立极坐标系，则每一次超空泡航行体的尾拍位置由μ_i*来计算：

步骤7：确定超空泡航行体两次尾拍间隔：尾拍前后角速度的变化关系式为：

Q⁺＝Q^-[1-G(1+e)]

其中：

Q^-为超空泡航行体尾拍前的超空泡航行体绕空化器的旋转角速度，Q⁺为超空泡航行体尾拍后的超空泡航行体旋转绕空化器的旋转角速度，x_cm为超空泡航行体质心到尾部的长度；

e为补偿系数，其表达式为：

其中U^-为超空泡航行体的速度，α为超空泡航行体尾拍角，即航行体尾拍时轴线与空泡壁面的夹角:

设前一次尾拍位置坐标为(x₁,y₁)，后一次尾拍位置为(x₂,y₂)，则两次尾拍的平面距离为：

超空泡航行体尾部的线速度为：

v_L＝|Q⁺|L

由此得两次尾拍的时间间隔为：

。

2.根据权利要求1所述的模拟方法，其特征在于：所述步骤4中超空泡航行体初始尾拍角速度为任意给定。

3.根据权利要求1所述的模拟方法，其特征在于，步骤6中所述超空泡航行体尾拍角度正定的具体实现方法为：

1:当|μ_i+1|>360计算2,否则计算7；

2:n＝[|μ_i+1|]；

3:ξ_i＝|μ_i+1|-360n；

4:μ_i+1＝ξ_isign(μ_i+1)；

5:如果sign(μ_i+1)＝-1则μ_i+1＝360-ξ_i,否则计算6；

6:μ_i+1＝ξ_i；

7:如果sign(μ_i+1)＝-1则μ_i+1＝360+ξ_i,否则计算8；

8:μ_i+1＝μ_i+1。

Images