修正PANS模型在非定常空化流动数值计算中的应用与评价

2018-11-30胡常莉曹友铨王学德

船舶力学 2018年11期

胡常莉，曹友铨，王学德

（南京理工大学能源与动力工程学院，南京210094）

0 引言

空化是发生在水力机械中一种常见的流动现象，往往会产生强烈的噪声、振动和空蚀[1]。如今，用数值计算的方法来研究空化现象已经比较广泛，尤其是对云状空化发展过程中出现的大空泡团脱落的模拟已经取得了一定进展。Shyy[2]，Coutier-Delgosha[3]等人的研究表明：湍流模型对于空化流动的预测有十分重要的影响。目前广泛采用的是基于雷诺平均N-S方程的涡黏模型，由于发生空化时流动中湍动能产生项和耗散项间的不平衡，采用标准的k-ε模型不能很好地模拟云状空化[4]。近年来发展了许多混合RANS/LES湍流模型，多位学者研究表明这种混合模型可以较好地应用到空化流动的计算中。黄彪等人[5]应用Spalart[6]提出的分离涡模型（DES），求解了非定常空化流动，达到了较好的效果。Wu等人[2]则将滤波器模型（FBM）[7]应用到了Clark-Y水翼的云状空化计算中，结果可以捕捉到非定常流动细节。另外，多位研究者将Girimaji[11]提出的基于标准k-ε模型的局部时均化模型（PANS），应用到不同绕流体的空化流动数值计算中，主要讨论了PANS模型中不同控制参数fk取值对计算结果的影响[8-10]。

由于Girimaji[11]提出的基于标准k-ε模型的PANS模型历史较短且有不足之处，近年来有学者对其进行了各种修正及发展，如Song等人[12]用基于k-w模型的PANS模型计算了圆柱绕流的流场并取得了较好的结果；AlaaElmiligui等人[13]借助于混合模型的思想对PANS模型的控制参数fk值进行了修正，应用PAB3D求解器研究了圆柱绕流问题，Hu等人[14]针对空化流动的特点提出了一种基于密度修正的PANS模型并取得较好的结果，Huang等人[15]基于当地网格尺度和湍流尺度对PANS模型进行了修正，并将其应用到后台阶绕流流动中。

在前人的研究基础之上，本文基于AlaaElmiligui等人[13]的研究，对标准PANS模型的控制参数fk值进行了修正，通过CFX二次开发技术加入到求解器中，模拟了绕Clark-Y型水翼的云状空化流动，基于数值及实验结果，评价并讨论了该模型的特性。

1 计算模型和数值方法

1.1 基本方程

采用均质平衡流模型，则Favre平均的N-S方程为：

式中：下标i和j分别代表坐标方向，ρm，u和p分别为混合介质的密度、速度和压强，μ和 μt分别为混合介质的层流和湍流粘性系数。

1.2 空化模型

空化流动计算中，选用Kubota空化模型[16]，其表达式如下：

式中：RB为简化气泡半径；pv为汽化压强；ρl和ρv分别为液体的密度和蒸汽的密度；αnuc为气核的体积分数；Fe和Fc分别是蒸发和凝结常数项。Kubota空化模型重点考虑了空化初生和发展时空泡体积变化的影响，适于模拟空化的非定常特性。

许多实验表明湍动能对空化产生重要的影响[17]，本文采用文献[17]中提出的方法来计算湍动能k对当地汽化压强的影响：

汽化压强采用下式计算：

式中：psat和k分别表示当地饱和蒸汽压强和流场的当地湍动能。

1.3 PANS模型及其修正

PANS模型的湍动能ku和耗散率εu的输运方程分别为：

湍动粘度：

其它常数取值分别为：

PANS模型的两个控制参数分别定义为[12]：

在高雷诺数的流动中，fε值通常取1，当fk=1时，说明湍流控制方程复原到RANS模型；当fk=0时，表示数值计算过程没有湍流模型的引入，为直接求解的方式。

本研究中，对于fk的取值，考虑到流场中各处的湍流尺度不同，采用（13）式[13]进行计算：

1.4 数值计算方法

1.4.1 计算网格和边界条件

计算采用Clark-y型水翼，弦长C=0.07 m。图2给出了计算区域及其边界条件。由图2所示，采用速度入口，压力出口，流动区域上下边界为自由滑移壁面条件，翼型表面采用绝热、无滑移固壁条件。攻角设定为8°，空化数设定为0.8，流速设定为10 m/s，雷诺数为7×105。图3为翼型周围网格分布图，对其近壁区域进行了网格加密，近壁面y+值为20～80之间，满足壁面函数要求。并且为了提高网格的质量，翼型前端的区域采用C型结构化网格划分，这样可以较好地匹配翼型头部的形状。

1.4.2 无量纲数的定义

计算中的主要无量纲参数为空化数σ、雷诺数Re、升力系数分别定义为：

式中：P∞、U∞和Pv分别为距实验段上游入口210 mm处参考断面上的平均静压强、断面平均速度和汽化压强，Fy是水翼所受到的升力。

2 结果分析及讨论

2.1 基于空化流场，对修正PANS模型的特性分析

修正PANS模型可以实现对流场的不同区域采用不同的fk值进行数值计算，即按照不同的湍流尺度，fk值的变化范围为0到1之间。为了研究fk的取值与空化流动结构之间的关系，图4给出了基于特征网格尺度Δ=2L的修正PANS模型得到的空穴形态随时间的演变过程及对应时刻的翼型周围湍动能和fk值分布云图。观察分析各时刻的空穴形态图及fk值、湍动能分布云图可知，在翼型吸力面的空化区域存在较大尺度的涡结构，且云状空化在发展的过程中伴随着大尺度空泡团的脱落现象，因此相比于远流场，翼型周围空化区域的湍动能较大，尤其是翼型尾部的旋涡区域。PANS模型的特点是控制参数fk值越小，流场中可解的湍流尺度越多，因此对于高湍动能区域，需要采用较小的fk值才能模拟出大尺度的湍流结构。由fk值云图可以看出，修正PANS模型恰好满足这个特点，即在翼型周围高湍动能区域对应的fk值较小，而远流场的fk值较大。图4中还给出了由标准k-ε模型得到的空穴形态图，通过对比发现，基于网格尺度Δ=2L的修正PANS模型可以较好地模拟出云状空化阶段翼型尾部的大空泡团脱落现象，这在一定程度上弥补了标准k-ε模型的不足。

为了进一步研究修正PANS模型中控制参数fk的取值与空化流场中湍流尺度之间的关系，图5分别给出了由基于4种网格尺度的修正PANS模型得到的时均fk值及时均湍流尺度分布云图。显然，4种情况下，时均湍流尺度较大的区域对应的fk值均较小，这满足对流场中大尺度湍流结构的求解要求。由时均湍流尺度分布云图可知，在翼型的尾部始终存在较大的湍流尺度结构，而在包裹翼型整个空化区域的剪切层内，由于该区域的耗散较大，导致其湍流尺度较小，从而该区域对应于较大的fk值。通过对比4种不同网格尺度的结果可知，网格尺度越小，对整个流场的求解中fk的取值越小。

2.2 修正PANS模型中特征网格尺度Δ的取值对空化流动计算的影响

由上一节的分析可知，修正PANS模型的控制参数fk的取值不仅与流场当地的湍流尺度有关而且受特征网格尺度的影响。因此，本小节主要讨论特征网格尺度Δ的取值对云状空化流动计算的影响。

图6给出了4种不同网格尺度下的时均蒸汽含量及湍流粘性分布云图。通过对比可知，时均蒸汽含量随着所选取网格尺度的减小而增大，相应地，湍流粘性系数则随着网格尺度的减小而减小。流场中当湍流尺度一定时，网格尺度越小，无量纲化后的湍流尺度λ值越大，结合图1可知，此时fk值较小，流场中会释放出较多的湍流尺度，减小了对湍流粘性的预测，从而减小了对空穴的抑制作用。

云状空化的发展过程具有明显的非定常特性，图7列出了4种不同的网格尺度时，修正PANS模型得到的空穴形态随时间的变化过程并与实验结果进行了对比。由实验结果分析可知，云状空化在一个发展周期内主要经历的形态特征为：翼型前缘出现附着空穴薄层并伴随着尾部大空泡团的脱落—前缘的空穴逐渐发展长大，尾部空泡团脱落溃灭消失—反向射流的形成发展，继而又造成大空泡的脱落。对比数值及实验结果可知，从空穴的整个发展过程来看，4种网格尺度时，修正PANS模型得到的空穴发展过程均与实验较一致，而主要的差异在于空穴发展的最大尺度及反向射流的发展过程。具体地，当t0+38%T时，Δ=0.1L的空穴形态与实验结果吻合较好，结合图8的空穴体积随时间的变化情况可知，随着网格尺度的减小，空穴体积随时间的波动幅度增大，表现为空穴形态的最大尺度较大。另外，当t0+62%T时，与实验结果对比可知，当Δ=0.1L时修正的PANS模型更好地模拟出了反向射流向翼型前缘推进的过程。

为了研究网格尺度Δ值对空化流场湍流流动特征的影响，图9给出了4种不同网格尺度下，图7中t0时刻对应的瞬时速度等值线图。由图7可知，该时刻为云状空化流动的翼型尾部大尺度空泡团的脱落。从图9中可以看出，当Δ值较小时，修正PANS模型计算得到的空泡团内部的湍流速度尺度分布较广且湍流速度脉动明显，随着Δ值的增大，流动速度值的分布范围逐渐减小，湍流速度脉动不明显。

图10分别给出了4种网格尺度时，修正PANS模型得到的翼型升力系数随时间的波动曲线，并将它们一一与实验结果进行了对比，发现与实验相一致的是，4种网格尺度时的升力系数曲线随时间的波动均具有明显的周期性。对比各曲线的波动细节可知，当网格尺度较大时，修正PANS模型不能释放出较多的湍流尺度，不足以模拟流动细节，因而升力系数曲线的波动细节不明显，而随着网格尺度的减小，升力系数曲线的小波动细节逐渐增多，说明此时修正PANS模型求解了较丰富的湍流尺度，与实验结果也较接近。