CN110852020B - 一种基于时间尺度的旋转湍流涡粘阻尼的计算方法 - Google Patents

Abstract

本发明涉及工程计算水力学，提供了一种基于时间尺度的旋转湍流涡粘阻尼的计算方法，包括：根据旋转湍流的特征频率获取解析时间尺度、模化时间尺度和相对螺旋度；根据解析时间尺度、模化时间尺度和相对螺旋度，获取旋转湍流的涡粘阻尼。本发明提供的基于时间尺度的旋转湍流涡粘阻尼的计算方法，涡粘阻尼的计算无需联系网格空间尺度，有利于保证数值计算的高效性和鲁棒性；新的时间尺度经验表达式能直接适应强瞬态、高脉动、大分离的流动工况；明确考虑了水力机械内显著存在的能量反级串输运现象。该计算方法在应用中能在适应流动特点的同时更好地兼顾计算精度与计算成本，可为水力机械内旋转湍流的高效解析提供新思路。

Description

一种基于时间尺度的旋转湍流涡粘阻尼的计算方法

技术领域

本发明涉及工程计算水力学，特别涉及一种基于时间尺度的旋转湍流涡粘阻尼的计算方法。

背景技术

近年来，在经典涡粘模型的基础上开发可兼顾计算精度与计算效率的湍流模型一直是工程计算水力学领域的热点工作。为保证工程流动的高效解析，MIT的Emilio教授利用时间尺度计算涡粘阻尼来调节湍动粘度，进而拓展湍流模型的解析能力。该方法在激活阻尼函数时无需联系网格的空间尺度，可降低数值计算对网格数量、质量及区域配置方案的要求，进而提高数值计算的效率和鲁棒性。

然而，在水轮机和水泵等诸多水力机械中，因受强旋转、大曲率、多壁面等因素的影响，内部旋转湍流常表现出强瞬态、高脉动、大分离等特点。在此种湍流的数值计算中，Emilio原始方法无法直接应用，存在以下显著的问题：(1)原方法中基于平面湍流特征而获得的平均时间尺度的经验表达式不合适旋转湍流的解析，且均化处理本身即无法适应强瞬态、大分离的流动工况；(2)原方法未明确反映旋转湍流中显著存在的能量反级串输运现象。以上不足均可导致模型无法准确预测旋转湍流的发展过程。

针对上述问题，本专利在继承原方法的高效性和稳健性的基础上，首次构造了旋转湍流特征时间尺度的经验表达式，并明确考虑了能量反级串输运现象，进而形成了一种新的旋转湍流涡粘阻尼计算方法，该方法能直接适应水力机械内旋转湍流的特点，可为水力机械中旋转湍流的高效解析提供新思路。

发明内容

(一)要解决的技术问题

鉴于上述技术缺陷和应用需求，本申请提出一种基于时间尺度的旋转湍流涡粘阻尼的计算方法，旨在解决现有方法无法直接适应水力机械内旋转湍流特点的问题，为水力机械程中旋转湍流的高效解析奠定科学基础。

(二)技术方案

为解决上述问题，本发明提供一种基于时间尺度的旋转湍流涡粘阻尼的计算方法，包括：

步骤一：引入涡母矩阵[S_C]以捕捉旋转湍流中局部水体的特征频率；

步骤二：根据涡母矩阵确定旋转湍流的解析时间尺度T_r；

步骤三：根据比耗散率确定旋转湍流的模化时间尺度T_m；

步骤四：引入相对螺旋度H_n以实现能量反级串输运现象的自适应追踪；

步骤五：根据解析时间尺度T_r、模化时间尺度T_m和相对螺旋度H_n确定旋转湍流涡粘阻尼T_D以保证混合模式激活的光滑性和有效性：

步骤六：在应用时，需以旋转湍流数值计算中经典的Baseline k-ω模型为基底，根据涡粘阻尼T_D确定旋转湍流高效解析所需的湍动粘度M_T，并对Baseline k-ω模型中的湍动能k与比耗散率ω的输运方程进行修正，将修正后的方程与雷诺方程结合形成控制方程组，利用有限体积法对控制方程组进行离散后即可直接用于水力机械内旋转湍流流场的数值计算；

其中，涡母矩阵[S_C]由下式得出：

[S_C]＝[R]²+[D]²

其中，[R]为旋转率张量，[D]为应变率张量；

解析时间尺度T_r由下式得出：

其中，λ_i(S_C)为涡母矩阵[S_C]的实特征值，f₁为解析时间尺度的经验常数；

模化时间尺度T_m由下式得出：

其中，ω为旋转湍流的比耗散率，f₂为模化时间尺度的经验常数；

相对螺旋度H_n由下式得出：

其中，V为速度矢量，

为速度旋度；

旋转湍流涡粘阻尼T_D由下式得出：

其中，d₁第一旋转湍流特征常数，d₂第二旋转湍流特征常数，d₃第三旋转湍流特征常数，d₄第四旋转湍流特征常数；

在应用时，需以旋转湍流数值计算中经典的Baseline k-ω模型为基底，根据下式：

其中，ρ为水流密度，k为旋转湍流湍动能，ω为旋转湍流比耗散率。

进一步地，第一旋转湍流特征常数d₁＝1.5708，第二旋转湍流特征常数d₂＝1.7183，第三旋转湍流特征常数d₃＝0.7274，第四旋转湍流特征常数d₄＝0.8653。

进一步地，解析时间尺度的经验常数f₁＝0.7071。

进一步地，模化时间尺度的经验常数f₂＝0.0845。

(三)有益效果

本发明提供的基于时间尺度的旋转湍流涡粘阻尼的计算方法，根据旋转湍流的解析时间尺度T_r、模化时间尺度T_m和相对螺旋度H_n，获取旋转湍流的涡粘阻尼T_D。涡粘阻尼的计算无需联系网格空间尺度，有利于保证数值计算的高效性和鲁棒性；新的时间尺度经验表达式能直接适应强瞬态、高脉动、大分离的流动工况；明确考虑了水力机械内显著存在的能量反级串输运现象。该计算方法在应用中能在适应流动特点的同时更好地兼顾计算精度与计算成本，可为水力机械内旋转湍流的高效解析提供新思路。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作一简单地介绍，显而易见地，下面描述中的附图是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为水翼绕流算例计算域示意图。

图2为泰勒-库埃特流动算例计算域示意图。

图3为水翼绕流算例监测位置处法向速度脉动值对比图。

图4为泰勒-库埃特流动算例监测位置处周向速度脉动值对比图。

图5为所述基于时间尺度的旋转湍流涡粘阻尼的计算方法及其应用的实现流程图。

具体实施方式

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

如图1和如图2所示分别为NACA0009水翼绕流算例和泰勒-库埃特流动算例的计算域示意图，此二算例均是水力机械中旋转湍流模型开发时的经典算例。在水翼绕流算例中，水翼弦长C＝100mm，尾缘厚度h＝3.22mm，展向宽度B＝150mm，上游来流流速V₀＝20m/s；在泰勒-库埃特流动算例中，圆筒内径R₁＝100mm，圆筒外径R₂＝200mm，内壁旋转速度U₀＝0.08m/s。对于此二算例，所述的基于时间尺度的旋转湍流涡粘阻尼的计算方法均可按如下过程实施。

本发明实施例提供一种基于时间尺度的旋转湍流涡粘阻尼的计算方法，包括：

步骤一：引入涡母矩阵[S_C]以捕捉旋转湍流中局部水体的特征频率；

步骤二：根据涡母矩阵确定旋转湍流的解析时间尺度T_r；

步骤三：根据比耗散率确定旋转湍流的模化时间尺度T_m；

步骤四：引入相对螺旋度H_n以实现能量反级串输运现象的自适应追踪；

步骤五：根据确定解析时间尺度T_r、模化时间尺度T_m和相对螺旋度H_n确定旋转湍流涡粘阻尼T_D以保证混合模式激活的光滑性和有效性：

步骤六：在应用时，需以旋转湍流数值计算中经典的Baseline k-ω模型为基底，根据涡粘阻尼T_D确定旋转湍流高效解析所需的湍动粘度M_T并对Baseline k-ω模型中的湍动能k与比耗散率ω的输运方程进行修正，将修正后的方程与雷诺方程结合形成控制方程组，利用有限体积法对控制方程组进行离散后即可直接用于水力机械内旋转湍流流场的数值计算，以获得水力机械内水流的速度和压力等流动参数。

涡母矩阵[S_C]由下式得出：

[S_C]＝[R]²+[D]²

其中，[R]为旋转率张量，[D]为应变率张量；

解析时间尺度T_r由下式得出：

其中，λ_i(S_C)为涡母矩阵[S_C]的实特征值，f₁为解析时间尺度的经验常数；

模化时间尺度T_m由下式得出：

其中，ω为旋转湍流的比耗散率，f₂为模化时间尺度的经验常数；

相对螺旋度H_n由下式得出：

其中，V为速度矢量，

为速度旋度；

旋转湍流涡粘阻尼T_D由下式得出：

其中，d₁第一旋转湍流特征常数，d₂第二旋转湍流特征常数，d₃第三旋转湍流特征常数，d₄第四旋转湍流特征常数。

以旋转湍流数值计算中经典的Baseline k-ω模型为基底，根据下式：

其中，ρ为水流密度，k为旋转湍流湍动能，ω为旋转湍流比耗散率。

如图5所示，该基于时间尺度的旋转湍流涡粘阻尼的计算方法包括如下步骤：

旋转湍流特征参数计算：按下式引入涡母矩阵[S_C]以捕捉旋转湍流中局部水体的特征频率，按下式引入相对螺旋度H_n以实现能量反级串输运现象的自适应追踪。

[S_C]＝[R]²+[D]²

式中，[R]为旋转率张量，[D]为应变率张量，V为速度矢量，

为速度旋度，<>为矢量夹角算符。

旋转湍流时间尺度计算：利用涡母矩阵和湍流比耗散率确定转湍流的解析时间尺度T_r和模化时间尺度T_m。

式中，λ_i(S_C)为涡母矩阵[S_C]的实特征值，f₁为解析时间尺度的经验常数；ω为旋转湍流的比耗散率，f₂为模化时间尺度的经验常数。优选的，解析时间尺度的经验常数的取值宜为0.7071。ω为旋转湍流的比耗散率，f₂为模化时间尺度的经验常数。优选的，模化时间尺度的经验常数f₂的取值宜为0.0845。

旋转湍流涡粘阻尼计算：利用旋转湍流时间尺度和特征参量确定高效解析所需的转湍流涡粘阻尼T_D。以保证混合模式激活的其光滑性和有效性。

式中，ln[]为自然对数函数，ch()为双曲余弦函数。优选的，旋转湍流特征常数d₁～d₄的取值依次为1.5708、1.7183、0.7274和0.8653。

输运方程修正：以旋转湍流数值计算中经典的Baseline k-ω模型为基底，根据下式

确定旋转湍流高效解析所需的湍动粘度。式中，ω为旋转湍流的比耗散率，ρ为旋转湍流的水流密度，k为旋转湍流的湍动能。利用所述湍动粘度M_T对Baseline k-ω模型中的湍动能k与比耗散率ω的输运方程进行修正，以获得高效解析所需的湍动能k与比耗散率ω的输运方程。

新方法应用：将修正后的方程与雷诺方程结合形成控制方程组，利用有限体积法对控制方程组进行离散后即可直接用于水力机械内旋转湍流流场的数值计算以获得水力机械内水流的速度和压力等流动参数。

需要说明的是，在实际计算过程中，可根据实际需要调整各步骤的顺序，在此不做限定。

本实施例中，采用高质量六面体网格对计算域进行空间离散，水翼绕流算例中总网格数约为397万，泰勒-库埃特流动算例中总网格数约为26万。将本专利提出的新方法引入后，采用全隐式耦合求解技术进行瞬态计算，时间步长均可满足工程计算流体力学中的CFL条件。另外，为更好地显示新方法的优点，同时取目前水力机械旋转湍流计算中常用的SST k-ω方法和背景技术中所述的Emilio原始法进行对比。

如图3所示为水翼绕流算例中X＝C+h位置处监测线上法向速度脉动值分布，其中包含各方法计算值与实验值的对比；如图4所示为泰勒-库埃特流动算例中半径方向上周向脉动速度空间平均值分布，其中包含各方法计算值与直接数值模拟值(DNS)的对比。分析可知，相比于现有的SST k-ω方法和Emilio原方法，本专利提出的新方法的数值模拟结果与实验值或DNS值更加接近，即能更加精确地获得旋转湍流中的流动信息，更适用于水力机械内的流动计算。

综上所述，本专利提出的一种基于时间尺度的旋转湍流涡粘阻尼的计算方法的优点在于：(1)涡粘阻尼的计算无需联系网格空间尺度，有利于保证数值计算的高效性和鲁棒性；(2)新的时间尺度的经验表达式能直接适应强瞬态、高脉动、大分离的流动工况；(3)计算方法明确考虑了水力机械内显著存在的能量反级串输运现象。因此，新方法在应用中能在适应流动特点的同时更好地兼顾计算精度与计算成本，可为水力机械内旋转湍流的高效解析提供新思路。

以上所描述的装置实施例仅仅是示意性的，其中所述作为分离部件说明的单元可以是或者也可以不是物理上分开的，作为单元显示的部件可以是或者也可以不是物理单元，即可以位于一个地方，或者也可以分布到多个网络单元上。可以根据实际的需要选择其中的部分或者全部模块来实现本实施例方案的目的。本领域普通技术人员在不付出创造性的劳动的情况下，即可以理解并实施。

最后应说明的是：以上实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的精神和范围。

Claims (4)

1.一种基于时间尺度的旋转湍流涡粘阻尼的计算方法，其特征在于，包括：

步骤一：引入涡母矩阵[S_C]以捕捉旋转湍流中局部水体的特征频率；

步骤二：根据涡母矩阵确定旋转湍流的解析时间尺度T_r；

步骤三：根据比耗散率确定旋转湍流的模化时间尺度T_m；

步骤四：引入相对螺旋度H_n以实现能量反级串输运现象的自适应追踪；

步骤五：根据确定解析时间尺度T_r、模化时间尺度T_m和相对螺旋度H_n确定旋转湍流涡粘阻尼T_D以保证混合模式激活的光滑性和有效性：

步骤六：在应用时，需以旋转湍流数值计算中经典的Baseline k-ω模型为基底，根据涡粘阻尼T_D确定旋转湍流高效解析所需的湍动粘度M_T，并对Baseline k-ω模型中的湍动能k与比耗散率ω的输运方程进行修正，将修正后的方程与雷诺方程结合形成控制方程组，利用有限体积法对控制方程组进行离散后即可直接用于水力机械内旋转湍流流场的数值计算；

涡母矩阵[S_C]由下式得出：

[S_C]＝[R]²+[D]²

其中，[R]为旋转率张量，[D]为应变率张量；

解析时间尺度T_r由下式得出：

其中，λ_i(S_C)为涡母矩阵[S_C]的实特征值，f₁为解析时间尺度的经验常数；

模化时间尺度T_m由下式得出：

其中，ω为旋转湍流的比耗散率，f₂为模化时间尺度的经验常数；

相对螺旋度H_n由下式得出：

其中，V为速度矢量，

为速度旋度；

旋转湍流涡粘阻尼T_D由下式得出：

其中，d₁第一旋转湍流特征常数，d₂第二旋转湍流特征常数，d₃第三旋转湍流特征常数，d₄第四旋转湍流特征常数；

以旋转湍流数值计算中经典的Baseline k-ω模型为基底，根据下式：

其中，ρ为水流密度，k为旋转湍流湍动能，ω为旋转湍流比耗散率。

2.根据权利要求1所述的基于时间尺度的旋转湍流涡粘阻尼的计算方法，其特征在于，第一旋转湍流特征常数d₁＝1.5708，第二旋转湍流特征常数d₂＝1.7183，第三旋转湍流特征常数d₃＝0.7274，第四旋转湍流特征常数d₄＝0.8653。

3.根据权利要求1所述的基于时间尺度的旋转湍流涡粘阻尼的计算方法，其特征在于，解析时间尺度的经验常数f₁＝0.7071。

4.根据权利要求1所述的基于时间尺度的旋转湍流涡粘阻尼的计算方法，其特征在于，模化时间尺度的经验常数f₂＝0.0845。

Images