CN112597710B - 一种可压缩空化流动中旋转湍流的数值模拟方法 - Google Patents

Abstract

本发明属于工程计算流体力学领域，尤其是涉及一种可压缩空化流动中旋转湍流的数值模拟方法，其实现过程主要包括：采用可压缩非线性亚格子模型计算涡粘系数和亚格子应力，引入调节因子修正亚格子应力，创建可压缩非线性亚格子模型基类，将可压缩非线性亚格子模型植入可压缩空化流动求解器。本发明考虑了可压缩空化流动中旋转湍流的影响，通过引入调节因子保证了可压缩空化流动计算的稳定性，能够有效降低可压缩空化流动模拟对网格的要求，提高了计算效率。同时，为深入研究可压缩空化流动中可压缩性对旋转湍流的影响提供新思路。

Description

一种可压缩空化流动中旋转湍流的数值模拟方法

技术领域

本发明涉及工程计算流体力学领域，更具体地说，涉及一种可压缩空化流动中旋转湍流的数值模拟方法。

背景技术

空化是高速水动力学中的核心问题。在水中，空化往往会诱发复杂的旋涡结构，使得空化流动被认为是一种旋转湍流，常常表现出非线性、强瞬态、高脉动等特点。同时，空化流动在本质上又是一种存在相变的汽液两相流动，导致空化流动具有强烈的可压缩特性。传统的商用CFD软件多是基于不可压缩框架开发，无法对空化流动的可压缩特性进行准确模拟。开源CFD软件OpenFOAM虽然可以进行可压缩空化流动数值模拟，但现有的湍流模型均为可压缩线性涡粘模型，该类模型虽然结构简单、计算方便，但该类模型在物理上即不成立，也无法充分考虑可压缩空化流动中湍流的强旋转效应。

发明内容

本发明要解决的技术问题在于，提供一种可压缩空化流动中旋转湍流的数值模拟方法，可模拟出空化可压缩性对旋转湍流的影响。

本发明解决其技术问题所采用的技术方案是：构造一种可压缩空化流动中旋转湍流的数值模拟方法，包括以下步骤：

步骤一：利用可压缩非线性亚格子模型计算涡粘系数ν和亚格子应力τ；

步骤二：引入调节因子η修正亚格子应力τ；

步骤三：创建可压缩非线性亚格子模型基类，该基类利用步骤一中的涡粘系数ν和步骤二中的修正亚格子应力τ求解湍流热扩散系数α^th和湍流应力P，并将计算结果传递到可压缩空化流动控制方程组，求解流场的温度T和速度U；

步骤四：将可压缩非线性亚格子模型植入可压缩空化流动求解器。

上述方案中，所述步骤一中的可压缩非线性亚格子模型的滤波方法以滤波对象乘以密度的形式出现，滤波公式如下：

式中，a为滤波对象，ρ为密度，a和ρ上方的短横线“-”表示空间滤波操作，通过空间滤波操作将流场中的物理量分为两部分，一部分称为大尺度可解量，可以直接利用方程进行求解，另一部分称为小尺度不可解量，需要利用物理模型模化。

上述方案中，所述步骤一中的可压缩非线性亚格子模型的涡粘系数ν和亚格子应力τ的计算公式如下：

式中：亚格子应力τ反映了一段时间内流场中小尺度物理量对求解的运动方程的影响；c_s、c₁、c₂均为模型系数，在计算过程中利用最小二乘法动态求解；Δ为滤波尺度用于分离大尺度可解量和小尺度不可解量；下标i,j,k分别代表空间中的x,y,z方向；下标m,n代表自由指标；Π_s为应变率张量S的第二主不变量；应变率张量

和旋转率张量

代表了速度变化和旋转运动导致的流体微团变形作用，二者的定义为：

式中：

代表了i方向的滤波速度在j方向上的变化率。

上述方案中，所述步骤二中的调节因子η和亚格子应力τ修正公式分别定义如下：

η＝min(0.8,max(τ_{ij_mean})) (13)

式中，τ_{ij_mean}为时均亚格子应力，其反映了一段时间内流场中小尺度物理量对求解的运动方程的影响，需要采用物理模型模化；调节因子η设置为一个常数，其数值根据计算流场的特点进行设定。

上述方案中，湍流热扩散系数α^th和湍流应力P求解公式如下：

式中，湍流热扩散系数α^th代表了流体传递热量的能量；Pr为普朗特数，代表了流体的动量扩散能力与其能量扩散能力的对比关系，其值在流场内取1，在壁面上取0.85；湍流应力P代表了流体作湍流运动时所产生的力；ν代表了流体的运动粘度；

代表了i方向的滤波速度在j方向上的变化率。

代表了亚格子应力在j方向上的变化率。

实施本发明的可压缩空化流动中旋转湍流的数值模拟方法，具有以下有益效果：

1、本发明采用一种可压缩非线性亚格子模型模拟可压缩空化流动中的旋转湍流，更加充分的考虑了流场的非定常、脉动和旋转效应；

2、本发明引入调节因子修正亚格子应力，保证了可压缩空化流动计算的稳定性，同时能够有效降低数值模拟对网格的要求，提高了计算效率。

3、本发明开发了一种可压缩非线性亚格子模型基类，该基类包含针对不同空化流动现象的可压缩非线性亚格子模型，可以根据所要计算的空化流动现象选择可压缩非线性亚格子模型。

附图说明

下面将结合附图及实施例对本发明作进一步说明，附图中：

图1为本发明可压缩空化流动中旋转湍流的数值模拟方法的实现流程图；

图2为水翼可压缩空化绕流算例计算示意图；

图3为水翼可压缩空化绕流流场矢量图。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清晰，下面将结合相关附图，对本发明实施例中的技术方案进行完整地描述。特别地，所描述的实施例并非本发明全部的实施例，基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动的前提下所获得的所有其他实施例都属于本发明保护的范围。

在本发明的描述中，所指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所描述的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明的限制。对于本领域的普通技术人员而言，可以根据具体情况理解上述术语在本发明中的具体含义。

如图1所示为NACA66水翼可压缩空化绕流算例的计算域示意图。该算例中，水翼弦长C＝150mm，水翼攻角为6°，空化数σ＝1.25，展向宽度B＝45mm，上游来流流速V₀＝5.33m/s。对于此算例，所述的一种可压缩空化流动中旋转湍流的数值模拟方法可按如下过程实施。

步骤一：计算涡粘系数ν和亚格子应力τ。

采用一种可压缩非线性亚格子模型计算流场的涡粘系数ν和亚格子应力τ，计算公式如下：

式中，亚格子应力τ反映了一段时间内流场中小尺度物理量对求解的运动方程的影响；c_s、c₁、c₂均为模型系数，在计算过程中利用最小二乘法动态求解；Δ为滤波尺度用于分离大瓷都可解量和小尺度不可解量；下标i,j,k分别代表空间中的x,y,z方向；下标m,n代表自由指标；Π_s为应变率张量S的第二主不变量；应变率张量

和旋转率张量

代表了速度变化和旋转运动导致的流体微团变形作用。

步骤二：调节因子η修正亚格子应力τ

在可压缩非线性亚格子模型中引入调节因子η来修正亚格子应力τ。调节因子η和亚格子应力τ的修正公式如下：

η＝min(0.8,max(τ_{ij_mean})) (17)

式中，τ_{ij_mean}为时均亚格子应力，其反映了一段时间内流场中小尺度物理量对求解的运动方程的影响，需要采用物理模型模化；调节因子η设置为一个常数，其数值根据计算流场的特点进行设定。通过C++语言将公式(1)、(2)、(3)和(4)写入OpenFOAM程序，由此建立可压缩非线性亚格子模型。

步骤三：创建可压缩非线性亚格子模型基类

基于OpenFOAM平台创建可压缩非线性亚格子模型基类ComGenNonlinearSGS，该基类可以利用步骤一中的涡粘系数ν和步骤二中的修正亚格子应力τ求解湍流热扩散系数α^th和湍流应力P，并将计算结果传递到可压缩空化流动控制方程组，求解流场的温度T和速度U。湍流热扩散系数α^th和湍流应力P求解公式如下：

式中，湍流热扩散系数α^th代表了流体传递热量的能量；Pr为普朗特数，代表了流体的动量扩散能力与其能量扩散能力的对比关系，其值在流场内取1，在壁面上取0.85；湍流应力P代表了流体作湍流运动时所产生的力；ν代表了流体的运动粘度；

代表了i方向的滤波速度在j方向上的变化率。

代表了亚格子应力在j方向上的变化率。

步骤四：可压缩非线性亚格子模型植入可压缩空化流动求解器

采用wmake命令编译步骤二中的可压缩非线性亚格子模型和步骤三中的可压缩非线性亚格子模型基类ComGenNonlinearSGS。在主程序的controlDict文件里面采用lib命令将可压缩非线性亚格子模型链接到可压缩空化流动求解器。在turbulenceProperties文件中调用可压缩非线性亚格子模型。

本实施例中采用有限体积法对控制方程进行离散；采用高质量六面体网格对计算域进行空间离散，总网格数约为250万。采用速度进口、压力出口边界条件，计算域壁面设置为无滑移边界条件，计算采用PISO算法进行压力-速度耦合，采用动态时间步长保证最大库朗数小于0.9。如图2所示为某一瞬时NACA66水翼可压缩空化绕流流场矢量图，分析可知，本专利提出的一种可压缩空化流动中旋转湍流的数值模拟方法能够精确捕捉到流场中的湍流涡结构，显示出了丰富的流场细节。另外，该方法中调节因子可以降低可压缩空化流动对网格的要求，提高计算效率。

综合来看，本发明提出的一种可压缩空化流动中旋转湍流的数值模拟方法，采用可压缩非线性亚格子模型，更加充分的考虑了可压缩空化流动中湍流的旋转效应。通过引入调节因子既保证了可压缩空化流动计算的稳定性，又降低了可压缩空化流动对网格的要求，提高了计算效率。该方法为可压缩空化流动中旋转湍流的计算提供了新思路。

上面结合附图对本发明的实施例进行了描述，但是本发明并不局限于上述的具体实施方式，上述的具体实施方式仅仅是示意性的，而不是限制性的，本领域的普通技术人员在本发明的启示下，在不脱离本发明宗旨和权利要求所保护的范围情况下，还可做出很多形式，这些均属于本发明的保护之内。

Claims (1)

1.一种可压缩空化流动中旋转湍流的数值模拟方法，其特征在于，包括以下步骤：

步骤一：利用可压缩非线性亚格子模型计算涡粘系数ν和亚格子应力τ；

步骤二：引入调节因子η修正亚格子应力τ；

步骤三：创建可压缩非线性亚格子模型基类，该基类利用步骤一中的涡粘系数ν和步骤二中的修正亚格子应力τ求解湍流热扩散系数α^th和湍流应力P，并将计算结果传递到可压缩空化流动控制方程组，求解流场的温度T和速度U；

步骤四：将可压缩非线性亚格子模型植入可压缩空化流动求解器；

所述步骤一中的可压缩非线性亚格子模型的滤波方法以滤波对象乘以密度的形式出现，滤波公式如下：

式中，a为滤波对象，ρ为密度；

所述步骤一中的可压缩非线性亚格子模型的涡粘系数ν和亚格子应力τ的计算公式如下：

式中：亚格子应力τ反映了一段时间内流场中小尺度物理量对求解的运动方程的影响；c_s、c₁、c₂均为模型系数，在计算过程中利用最小二乘法动态求解；Δ为滤波尺度用于分离大尺度可解量和小尺度不可解量；下标i,j,k分别代表空间中的x,y,z方向；下标m,n代表自由指标；Π_s为应变率张量S的第二主不变量；应变率张量

和旋转率张量

代表了速度变化和旋转运动导致的流体微团变形作用，二者的定义为：

式中：

代表了i方向的滤波速度在j方向上的变化率；

所述步骤二中的调节因子η和亚格子应力τ修正公式分别定义如下：

η＝min(0.8,max(τ_{ij_mean})) (6)

式中，τ_{ij_mean}为时均亚格子应力，其反映了一段时间内流场中小尺度物理量对求解的运动方程的影响，需要采用物理模型模化；调节因子η设置为一个常数，其数值根据计算流场的特点进行设定；

湍流热扩散系数α^th和湍流应力P求解公式如下：

式中，湍流热扩散系数α^th代表了流体传递热量的能量；Pr为普朗特数；湍流应力P代表了流体作湍流运动时所产生的力；ν代表了流体的运动粘度；

代表了i方向的滤波速度在j方向上的变化率；

代表了亚格子应力在j方向上的变化率。

Images