CN110516342A - 一种基于OpenFOAM平台的螺旋桨可压缩空化流动数值预测方法 - Google Patents

Abstract

本发明公开了一种基于OpenFOAM平台的螺旋桨可压缩空化流动数值预测方法，该方法首先建立螺旋桨模型及三维计算流域并划分网格，然后建立由可压缩空化流动控制方程组、热力学模型、可压缩空化模型和可压缩湍流模型组成的可压缩空化流动计算流体动力学模型；其次进行边界条件和计算设置；再进行可压缩空化流场数值计算；最后获得螺旋桨的非定常可压缩空化流场特性，通过求解可压缩空化流动计算流体动力学模型得到螺旋桨流场速度分布、压力分布、相分布和温度分布。本发明能够实现对螺旋桨空化流动诱导压力脉动、压力脉冲、激波等现象的精确模拟，为建立螺旋桨空化破坏预测及控制提供依据。

Description

一种基于OpenFOAM平台的螺旋桨可压缩空化流动数值预测 方法

技术领域

本发明涉及一种基于OpenFOAM平台的螺旋桨可压缩空化流动数值预测方法，属于计算流体动力学、空化与多相流技术领域。

背景技术

空化是高速水动力学的核心关键科学问题，不可避免发生于水力机械、水下兵器等运行工况下，诱导复杂流体动力。空化问题是水力与海洋工程装备最严酷的服役环境，制约着我国水力与海洋工程装备水动力创新发展，深入对空化流动的认识是进一步探究空化不稳定性机理、预测及抑制空化破坏的重要基础。随着船舶向大型化和高速化发展，作为海洋船舶推进系统的重要组成的螺旋桨及舵的空化及空化破坏问题日益严峻，精确地预测螺旋桨空化流场特性是高速水动力学的重点问题。

实验是空化研究的主要方法，然而受可接触性以及测量手段有限的制约，空化流场内部结构精细测量存在困难。对于螺旋桨等大型水力机械，存在着原型与模型之间的尺寸效应、可重复性等问题，对其规律的解释高度依赖于经验，并且实验室工作无法对材料剥蚀现象进行充分长时间测量。因此采用计算流体动力学(CFD)数值计算方法对实验进行补充具有重要意义。目前，主流的空化流动数值计算多基于商业软件，如Ansys-CFX，Ansys-Fluent，等，采用不可压缩数值计算方法，不可压缩空化流动计算方法可以对非定常空化流动现象，如附着型空穴生长，回射流产生及发展，空泡断裂脱落，以及流体动力等进行预测。然而，进一步精确预测螺旋桨空化流体动力特性，建立空化破坏准则，亟需对空化诱导压力脉动、压力脉冲、激波等现象的精确模拟，这些现象传统的不可压缩计算方法无法进行预测捕捉。理论研究表明，空化流场由于含气率的存在，空化介质中局部声速急剧下降，可降至3-5m/s，空化流动具有高度的压缩性。为深入空化流动机理，精确预测空穴溃灭机制以及其诱导激波压力脉冲，必须综合考虑空化汽液介质可压缩性效应，采用可压缩数值计算方法。由于目前商业软件多基于不可压缩框架，在可压缩空化计算方法上存在不足，OpenFOAM作为大型开源计算流体力学类库，具有良好的可拓展性，便于科研人员进行开发，可针对所解决的问题定制自己的求解器，并且已逐渐应用于空化流动中，因此，基于开源软件OpenFOAM进行可压缩空化流动求解器的开发具有必要性和可行性。

发明内容

有鉴于此，本发明提供了一种基于OpenFOAM平台的螺旋桨可压缩空化流动数值预测方法，能够实现对螺旋桨空化流动诱导压力脉动、压力脉冲、激波等现象的精确模拟，为建立螺旋桨空化破坏预测及控制提供依据。

一种基于OpenFOAM平台的螺旋桨可压缩空化流动数值预测方法，该方法的实现包括以下步骤：

步骤一：建立螺旋桨模型及三维计算流域并划分网格；

步骤二：建立由可压缩空化流动控制方程组、热力学模型、可压缩空化模型和可压缩湍流模型组成的可压缩空化流动计算流体动力学模型；

步骤三：进行边界条件和计算设置；

步骤四：进行可压缩空化流场数值计算；

步骤五：获得螺旋桨的非定常可压缩空化流场特性，通过求解可压缩空化流动计算流体动力学模型得到螺旋桨流场速度分布、压力分布、相分布和温度分布。

进一步地，所述步骤一中的三维计算流域采用上游入口边界距螺旋桨盘面为0.6-0.9倍螺旋桨直径，下游出口边界距螺旋桨盘面为4-6倍螺旋桨直径，外部圆柱壁面边界距螺旋桨桨轴为3倍螺旋桨直径。网格划分分为内域和外域，外域采用高质量的结构化网格划分，考虑到螺旋桨叶面倾斜大，几何形状复杂，内域采用混合网格。

进一步地，所述步骤二中可压缩空化流动控制方程组包括可压缩连续性方程、动量方程、能量方程以及相体积分数输运方程，如下所列：

其中，ρ、U、p、e和K分别为密度、速度、压力、内能和动能，μ为介质动力粘度，I为单位张量，α为相体积分数，α_eff＝α₁κ₁/c_v,1+α_vκ_v/c_v,v是有效热扩散系数，κ₁,κ_v,c_v,l和c_v,v分别是水相和蒸汽相的热导率和比热，m⁺和m^-为空化源项，反映汽液相间质量传输过程，下标m，l和v代表水汽混相介质、液相和蒸汽相；表面张力σ＝0.0728；方程(1)-(4)采用水相和蒸汽相状态方程进行热力学闭合。

进一步地，所述步骤二中的热力学模型由水相和蒸汽相状态方程组成，水相采用Tait状态方程：

其中，p_l,sat(T)和ρ_l,sat(T)为温度决定的水介质饱和蒸汽压力和饱和密度，B＝3.06×10⁸Pa，N＝7.1；

蒸汽相采用理想气体状态方程：

p_v＝ρ_vR_vT_v (6)

其中，R_v＝461.6J/(kg.K)；

水介质饱和蒸汽压力采用温度修正公式：

其中，θ＝1-T/T_c，下标“c”表示水介质临界状态参数，p_c＝22.064MPa，ρ_c＝322kg/m3，T_c＝647.14K；

水介质饱和蒸汽密度采用温度修正公式：

ρ_l(T)＝ρ_c(1+1.99206θ^1/3+1.10123θ^2/3-0.512506θ^5/3-1.75263θ^16/3-45.4485θ^43/3-6.75615×10⁵θ^110/3) (8)

饱和水蒸汽密度采用温度修正的公式：

ρ_v,sat(T)＝ρ_cexp(-2.02957θ^2/6-2.68781θ^4/6-5.38107θ^8/6-17.3151θ^18/6-44.6384θ^37/6-64.3486θ^71/6) (9)。

进一步地，所述步骤二中的可压缩空化模型为：

ρ＝αρ₁+(1-α)ρ₂ (10)

由Wallis声速公式得混相介质声速为，

因此，该可压缩空化模型的蒸发速率和凝结速率分别为

当p<p_v时，

当p>p_v时，

其中，K_l＝472.26J/kgK，T₀＝3836.94K，ρ_c＝1944.61MPa，C_c和C_e分别为空化模型凝结系数和蒸发系数，选择C_c＝C_e＝0.1。

进一步地，所述步骤二中的可压缩湍流模型采用WALE大涡模拟数值计算方法，在可压缩空化湍流中，该湍流模型亚格子应力为：

其中S_ij为可解尺度的应变率张量，μ_t为亚格子涡粘性系数；

其中L_s为亚格子混合长度，k为冯卡门常数，d为距螺旋桨桨叶壁面最近距离，Δx,Δy,Δz分别为计算单元尺寸，C_s＝0.5为WALE常数。

进一步地，所述步骤三中，设置三维流域边界条件，进口处采用速度边界条件，出口采用压力边界条件，壁面采用无滑移边界条件，螺旋桨运动采用OpenFOAM中刚体运动类型进行控制，运动方式为使其按照设定的转动轴与角度速度运动，内部动网格与外部静止网格之间的滑移面采用OpenFOAM提供的周期性任意网格界面方法；采用PISO算法进行速度-压力耦合求解，时间项采用二阶精度的隐式算法，时间步选择保证流场中最大库郎数小于0.4，空间项采用二阶精度算法，梯度项采用对计算单元进行非正交修正；在相体积分数输运方程求解中采用基于二阶差分并且有界的算法。

有益效果：

1、本发明的可压缩空化预测方法，通过引入空化汽液介质状态方程以及空化介质能量方程，综合考虑了空化密度脉动、温度脉动及其耦合作用对空化流动特性及动力的影响，该可压缩空化流动模型更加全面地反映了空化流动物理本质，可提高空化预测精度。

2、本发明的可压缩空化数值预测方法通过考虑空化汽、液及其混相介质压缩性效应对相间质量传输过程的影响，对空化模型进行了压缩性效应修正，该可压缩空化模型可进一步提高空化模拟精度，可获得准确的螺旋桨空化流场特性。

3、本发明的可压缩空化预测方法基于OpenFOAM开源软件平台开发的可压缩空化流动求解器，综合考虑空化气液介质压缩性效应，以及空化介质物质属性受温度变化影响，以及高效率动网格技术，能够提高对螺旋桨空化流场及动力特性预测的精度。

附图说明

图1是本发明一种基于OpenFOAM平台的螺旋桨可压缩空化流动数值预测方法的流程图；

图2是本发明螺旋桨可压缩空化流动数值计算三维流域示意图。

具体实施方式

下面结合附图并举实施例，对本发明进行详细描述。

如附图1所示，本发明提供了一种基于OpenFOAM平台的螺旋桨可压缩空化流动数值预测方法，该方法的实现步骤如下：

步骤一：建立螺旋桨模型及三维计算流域并划分网格

采用建模软件进行螺旋桨建模，建立三维螺旋桨计算流域，并划分网格。如附图2所示，三维计算流域采用上游入口边界距螺旋桨盘面为0.7倍螺旋桨直径，下游出口边界距螺旋桨盘面为5倍螺旋桨直径，外部圆柱壁面边界距螺旋桨桨轴为3倍螺旋桨直径。网格划分分为内域和外域，外域采用高质量的结构化网格划分，考虑到螺旋桨叶面倾斜大，几何形状复杂，内域采用混合网格。

步骤二：建立可压缩空化流动计算流体动力学模型

为精确预测螺旋桨空化流动诱导流体动力特性，基于OpenFOAM平台开发的可压缩空化流动求解器采用可压缩空化流动数值计算方法，可压缩空化流动计算流体动力学模型由可压缩空化流动控制方程组、热力学模型、可压缩空化模型和可压缩湍流模型组成，可压缩空化流动控制方程组包括可压缩连续性方程、动量方程、能量方程以及相体积分数输运方程，如下所列：

其中，ρ、U、p、e和K分别为密度、速度、压力、内能和动能，μ为介质动力粘度，I为单位张量，α为相体积分数，α_eff＝α₁κ₁/c_v,1+α_vκ_v/c_v,v是有效热扩散系数，κ₁,κ_v,c_v,l和c_v,v分别是水相和蒸汽相的热导率和比热，m⁺和m^-为空化源项，反映汽液相间质量传输过程，下标m，l和v代表水汽混相介质、液相和蒸汽相。考虑到大尺度空穴溃灭过程末期空泡界面急剧减小，表面张力影响增大，在数值计算方法中考虑表面张力，σ＝0.0728。方程(1)-(4)采用水相和蒸汽相状态方程进行热力学闭合。

热力学模型由水相和蒸汽相状态方程组成，水相采用Tait状态方程：

其中，p_l,sat(T)和ρ_l,sat(T)为温度决定的水介质饱和蒸汽压力和饱和密度，B＝3.06×10⁸Pa，N＝7.1。

蒸汽相采用理想气体状态方程：

p_v＝ρ_vR_vT_v (6)

其中，R_v＝461.6J/(kg.K)。

水介质饱和蒸汽压力采用温度修正的公式：

其中，θ＝1-T/T_c，下标“c”表示水介质临界状态参数，p_c＝22.064MPa，ρ_c＝322kg/m3，T_c＝647.14K。

水介质饱和蒸汽密度采用温度修正的公式：

ρ_l(T)＝ρ_c(1+1.99206θ^1/3+1.10123θ^2/3-0.512506θ^5/3-1.75263θ^16/3-45.4485θ^43/3-6.75615×10⁵θ^110/3) (8)

饱和水蒸汽密度采用温度修正的公式：

ρ_v,sat(T)＝ρ_cexp(-2.02957θ^2/6-2.68781θ^4/6-5.38107θ^8/6-17.3151θ^18/6-44.6384θ^37/6-64.3486θ^71/6) (9)

为精确预测空化流动压缩性效应对空化流场结构的影响，需建立可压缩空化模型。在可压缩空化流动中，流动-热力学特性耦合效应显著，为精确预测空化流体动力学特性，需对相间质量交换过程进行压缩性效应和热力学参数修正。在现有Saito空化模型中，该模型基于平板蒸发冷凝理论建立，考虑了温度变化对空化汽液相间质量交换的影响。为考虑空化压缩性效应对汽液相间质量交换速率的作用，本发明以Saito空化模型为基础，同时考虑空化汽液及其混相介质密度变化对相间质量交换过程影响，建立了可压缩空化模型。

ρ＝αρ₁+(1-α)ρ₂ (10)

由Wallis声速公式得混相介质声速为，

因此，该可压缩空化模型的蒸发速率和凝结速率分别为

当p<p_v时，

当p>p_v时，

其中，K_l＝472.26J/kgK，T₀＝3836.94K，ρ_c＝1944.61MPa，C_c和C_e分别为空化模型凝结系数和蒸发系数，计算中，选择C_c＝C_e＝0.1。

为精确模拟空化流动可压缩湍流效应，需采用高精度可压缩湍流模型。为捕捉螺旋桨空化流动细节，采用WALE大涡模拟数值计算方法。

亚格子应力为：

其中S_ij为可解尺度的应变率张量，μ_t为亚格子涡粘性系数。

其中L_s为亚格子混合长度，k为冯卡门常数，d为距螺旋桨桨叶壁面最近距离，Δx,Δy,Δz分别为计算单元尺寸，C_s＝0.5为WALE常数。

步骤三：边界条件及计算设置

在OpenFOAM软件中，设置三维流域边界条件，进口处采用速度边界条件，出口采用压力边界条件，壁面采用无滑移边界条件，螺旋桨运动采用OpenFOAM中刚体运动类型(solidBodyMotionFvMesh)进行控制，运动方式为rotatitingMotion使其按照设定的转动轴与角度速度运动，内部动网格与外部静止网格之间的滑移面采用OpenFOAM提供的周期性任意网格界面(cyclicAMI，cyclic Arbitrary Mesh Interface)方法。采用PISO算法进行速度-压力耦合求解，时间项采用二阶精度的隐式算法，时间步选择保证流场中最大库郎数小于0.4，空间项采用二阶精度算法，梯度项采用对计算单元进行非正交修正。在相体积分数输运方程求解中，采用基于二阶差分并且严格保证有界的算法。

步骤四：进行可压缩空化流场数值计算

基于OpenFOAM平台开发的可压缩空化流动求解器，以步骤五中所设置的边界条件和计算方法，进行非定常可压缩空化流动数值计算，得到螺旋桨空化流动三维可压缩数值计算结果。

步骤五：获得螺旋桨的非定常可压缩空化流场特性，通过求解可压缩空化流动计算流体动力学模型得到螺旋桨流场速度分布、压力分布、相分布和温度分布。即实现对螺旋桨可压缩空化流动现象及其流体动力特性的高精度数值预测。对螺旋桨空化流场进行后处理，获得非定常空穴形态分布及其演化、螺旋桨叶片不同测点压力脉动曲线等螺旋桨动力特性、非定常空穴结构演化与螺旋桨运动关系以及螺旋桨推力等水力性能。

综上所述，以上仅为本发明的较佳实施例而已，并非用于限定本发明的保护范围。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (7)

1.一种基于OpenFOAM平台的螺旋桨可压缩空化流动数值预测方法，其特征在于，该方法的实现包括以下步骤：

步骤一：建立螺旋桨模型及三维计算流域并划分网格；

步骤二：建立由可压缩空化流动控制方程组、热力学模型、可压缩空化模型和可压缩湍流模型组成的可压缩空化流动计算流体动力学模型；

步骤三：进行边界条件和计算设置；

步骤四：进行可压缩空化流场数值计算；

步骤五：获得螺旋桨的非定常可压缩空化流场特性，通过求解可压缩空化流动计算流体动力学模型得到螺旋桨流场速度分布、压力分布、相分布和温度分布。

2.如权利要求1所述的基于OpenFOAM平台的螺旋桨可压缩空化流动数值预测方法，其特征在于，所述步骤一中的三维计算流域采用上游入口边界距螺旋桨盘面为0.6-0.9倍螺旋桨直径，下游出口边界距螺旋桨盘面为4-6倍螺旋桨直径，外部圆柱壁面边界距螺旋桨桨轴为3倍螺旋桨直径。网格划分分为内域和外域，外域采用高质量的结构化网格划分，考虑到螺旋桨叶面倾斜大，几何形状复杂，内域采用混合网格。

3.如权利要求1或2所述的基于OpenFOAM平台的螺旋桨可压缩空化流动数值预测方法，其特征在于，所述步骤二中可压缩空化流动控制方程组包括可压缩连续性方程、动量方程、能量方程以及相体积分数输运方程，如下所列：

其中，ρ、U、p、e和K分别为密度、速度、压力、内能和动能，μ为介质动力粘度，I为单位张量，α为相体积分数，α_eff＝α₁κ₁/c_v,1+α_vκ_v/c_v,v是有效热扩散系数，κ₁,κ_v,c_v,l和c_v,v分别是水相和蒸汽相的热导率和比热，m⁺和m^-为空化源项，反映汽液相间质量传输过程，下标m，l和v代表水汽混相介质、液相和蒸汽相；表面张力σ＝0.0728；方程(1)-(4)采用水相和蒸汽相状态方程进行热力学闭合。

4.如权利要求3所述的基于OpenFOAM平台的螺旋桨可压缩空化流动数值预测方法，其特征在于，所述步骤二中的热力学模型由水相和蒸汽相状态方程组成，水相采用Tait状态方程：

其中，p_l,sat(T)和ρ_l,sat(T)为温度决定的水介质饱和蒸汽压力和饱和密度，B＝3.06×10⁸Pa，N＝7.1；

蒸汽相采用理想气体状态方程：

p_v＝ρ_vR_vT_v (6)

其中，R_v＝461.6J/(kg.K)；

水介质饱和蒸汽压力采用温度修正公式：

其中，θ＝1-T/T_c，下标“c”表示水介质临界状态参数，p_c＝22.064MPa，ρ_c＝322kg/m3，T_c＝647.14K；

水介质饱和蒸汽密度采用温度修正公式：

ρ_l(T)＝ρ_c(1+1.99206θ^1/3+1.10123θ^2/3-0.512506θ^5/3-1.75263θ^16/3-45.4485θ^43/3-6.75615×10⁵θ^110/3) (8)

饱和水蒸汽密度采用温度修正的公式：

ρ_v,sat(T)＝ρ_cexp(-2.02957θ^2/6-2.68781θ^4/6-5.38107θ^8/6-17.3151θ^18/6-44.6384θ^37/6-64.3486θ^71/6) (9)。

5.如权利要求4所述的基于OpenFOAM平台的螺旋桨可压缩空化流动数值预测方法，其特征在于，所述步骤二中的可压缩空化模型为：

ρ＝αρ₁+(1-α)ρ₂ (10)

由Wallis声速公式得混相介质声速为，

因此，该可压缩空化模型的蒸发速率和凝结速率分别为

当p<p_v时，

当p>p_v时，

其中，K_l＝472.26J/kgK，T₀＝3836.94K，ρ_c＝1944.61MPa，C_c和C_e分别为空化模型凝结系数和蒸发系数，选择C_c＝C_e＝0.1。

6.如权利要求5所述的基于OpenFOAM平台的螺旋桨可压缩空化流动数值预测方法，其特征在于，所述步骤二中的可压缩湍流模型采用WALE大涡模拟数值计算方法，在可压缩空化湍流中，该湍流模型亚格子应力为：

其中S_ij为可解尺度的应变率张量，μ_t为亚格子涡粘性系数；

其中L_s为亚格子混合长度，k为冯卡门常数，d为距螺旋桨桨叶壁面最近距离，Δx,Δy,Δz分别为计算单元尺寸，C_s＝0.5为WALE常数。

7.如权利要求6所述的基于OpenFOAM平台的螺旋桨可压缩空化流动数值预测方法，其特征在于，所述步骤三中，设置三维流域边界条件，进口处采用速度边界条件，出口采用压力边界条件，壁面采用无滑移边界条件，螺旋桨运动采用OpenFOAM中刚体运动类型进行控制，运动方式为使其按照设定的转动轴与角度速度运动，内部动网格与外部静止网格之间的滑移面采用OpenFOAM提供的周期性任意网格界面方法；采用PISO算法进行速度-压力耦合求解，时间项采用二阶精度的隐式算法，时间步选择保证流场中最大库郎数小于0.4，空间项采用二阶精度算法，梯度项采用对计算单元进行非正交修正；在相体积分数输运方程求解中采用基于二阶差分并且有界的算法。