标准k-ε 模型[1]作为各向同性假设的经典湍流涡黏模型,因其强大的通用性与可靠性,在工程界的流体流动数值模拟中得到广泛的应用[2-4]。而无论是现在常用于科学研究的直接数值模拟(direct numerical simulation,DNS),还是介于两者之间的大涡模拟(large eddy simulation,LES),由于其所消耗的计算资源较多,目前还难以进行广泛的工程实践应用。

但由于标准k-ε 模型在对雷诺应力各个分量的计算中,所假定的湍流黏度为各向同性的标量,使得其计算模型仅适用于雷诺数较高的情况。而无法对流体分子黏性力起较大影响的壁面黏性底层及其过渡区域进行准确地计算。对此,其处理方法是,不对黏性影响比较明显的区域(黏性底层和过渡区)进行求解,而是用一组半经验的公式(壁面函数)将壁面上的物理量与湍流核心区内的相应物理量联系起来,这就是壁面函数法。结合所谓的标准壁面函数(standard wall function,SWF)[1,5],不但避免了壁面附近流体分子黏性的影响,而且在划分网格的时候,也不需要在壁面区加密,只需要把第一个节点布置在边界层的对数律成立的区域内(y  30),即配置在边界层的湍流充分发展区域。而对较密的边界网格,则采用增强型壁面函数(Enhanced wall treatment,EWT)(双层模型)[5]。以上处理方法在工程界对壁面边界层较薄的高雷诺数流动模拟带来了便利,但对雷诺数较低的湍流流动,特别是模拟壁面附近的流动,会存在较大的误差。

为了能够让基于k-ε 模型的数值计算能从较高Re数的湍流主流区域一直延伸到固体壁面上,对壁面区域进行准确求解,需要对标准k-ε 模型进行修正,使得修正后的模型适应不同的雷诺数区域。所得到的模型即为所谓的低雷诺数k-ε 湍流模型,见式(1)—(5)。

附壁剪切流动是建筑物通风及工业干燥与冷却等过程中的常见流动现象。本文采用Fluent14.0中增强型壁面函数(EWT)的标准k-ε 模型[5](SKE)及6种低雷诺数模型：AB模型[6],LB模型[7],LS模型[8],YS模型[9],AKN模型[10],CHC模型[11],对三维附壁剪切流进行了数值模拟。根据湍流的壁面限定条件对这6种低雷诺数模型进行了理论分析;并对模拟结果与实验数据进行对比验证;其后,对它们所得的模拟结果是否符合湍流壁面流动特性的情况进行了分析;从这3方面对相关模型进行了分析与验证。以期为工业数值模拟中低雷诺数模型的选取提供较为可靠的参考与依据。

近半个多世纪以来,诸多研究者提出了超过20种不同的低雷诺数k-ε 湍流模型[9-13],前期很多模型中的fμ 等阻尼函数中存在y,使得回流区脱离与附着点处存在奇点的问题。其后,很多模型采用Kolmogorov速度uε 替代y 表达式中的壁面摩擦速度uτ,克服了模型这方面的问题。Fluent中选用的6种低雷诺数模型都不存在这方面的问题。Fluent是目前应用最广泛的通用CFD软件,但迄今为止,还没有文献对其所包含的这6种模型就三维附壁剪切流流场进行统一的比较与分析。

连续性方程：

\(\frac{\partial {{U}_{i}}}{\partial {{x}_{i}}}=0\) (1)

动量方程：

\(\frac{\text{D}{{U}_{i}}}{\text{D}t}=-\frac{1}{\rho }\frac{\partial P}{\partial {{x}_{i}}}+\frac{\partial }{\partial {{x}_{i}}}(\nu \frac{\partial {{U}_{i}}}{\partial {{x}_{j}}}-\overline{{{u}_{i}}{{u}_{j}}})\) (2)

式中线性的涡黏性模型采用应力和应变关系如下：

\(\overline{{{u}_{i}}{{u}_{j}}}=\frac{2}{3}k{{\delta }_{ij}}-{{\nu }_{t}}(\frac{\partial {{U}_{i}}}{\partial {{x}_{j}}}+\frac{\partial {{U}_{j}}}{\partial {{x}_{i}}})\) (3)

\(\frac{\text{D}k}{\text{D}t}=\frac{\partial }{\partial {{x}_{j}}}[(\nu +\frac{{{\nu }_{t}}}{{{\sigma }_{k}}})\frac{\partial k}{\partial {{x}_{j}}}]+{{P}_{k}}-\tilde{\varepsilon }-D\) (4)

\(\begin{align} \frac{\text{D}\tilde{\varepsilon }}{\text{D}t}=\frac{\partial }{\partial {{x}_{j}}}[(\nu +\frac{{{\nu }_{t}}}{{{\sigma }_{\varepsilon }}})\frac{\partial \tilde{\varepsilon }}{\partial {{x}_{j}}}]+{{C}_{\varepsilon 1}}{{f}_{1}}\frac{1}{{{T}_{t}}}{{P}_{k}}- \\ \quad \quad {{C}_{\varepsilon 2}}{{f}_{2}}\frac{{\tilde{\varepsilon }}}{{{T}_{t}}}+E \\\end{align}\) (5)

式中：fμ、f1、f2、\(\tilde{\varepsilon }\)、D和E随给出的模型而异,其中fμ、f1和f2是当地湍流雷诺数的阻尼函数;\(\tilde{\varepsilon }\)为修正后的各向同性耗散率,与 ε 有关;D、E是附加的源项;Pk为湍流动能产生项;Tt为湍流时间尺度[14]。湍流涡黏系数在低雷诺数k-ε 模型中定义为ν t = Cμ fμ k2/ε,Cμ 是常数。

低雷诺数区域的特点是分子黏性对湍流动量的直接影响,湍流雷诺数Ret(= k2/(νε ))小于150时,分子黏性影响较为显著,通常称为低雷诺数流动[15]。低雷诺数k-ε 湍流模型使用阻尼函数( fμ、f1、f2)以及附加项(D、E),使得控制方程在靠近墙壁的低雷诺数的区域变得合法,从而可直接对黏性底层进行求解,而不需要象标准k-ε 模型那样采用壁面函数去“桥接”层流底层。其中fμ 降低了靠近壁面的湍流黏度,是各种低雷诺数k-ε 模型最重要的阻尼函数;f1与湍流动能产生项Pk相乘,增强了壁面附近的湍流耗散率,所有模型中,仅YS模型使用了阻尼函数f1,其余模型的f1皆可看成常数1;f2与耗散方程的消失项相乘,将各向同性湍流假设在低雷诺数区域的偏差考虑了进去[16]。各模型的阻尼函数表达式见表1。

附加项D出现在湍流动能方程中,是由于LS模型采用了0neumann边界条件。Jones和Launder[17]曾论述了求解耗散方程时采用零边界条件的优势。因此,附加项D就等于壁面处湍流耗散率的大小,但该方法也并非被所有低雷诺数模型采用,见表2。E项加入耗散方程中,是部分低雷诺数模型[8-9]为得到近壁面较合理的湍流动能k的分布。部分模型(AB、AKN)对其部分模型常数Cε 1、Cε 2、Cμ、σ k、

表1 各湍流模型中阻尼系数的表达式 Tab. 1 The damping functions used in the Low Re. turbulence model

表2 模型常数及随着边界条件的改变各个模型源项(D、E)的表达式 Tab. 2 The model constants, and additional terms (D and E) in different wall boundary conditions

σε 也进行了修正,其值相对于标准k-ε 模型有所改变,具体见表2。模型常数的意义及影响在文献[16]进行了描述。对此,不同研究者具有不同的看法,多数研究者认为,标准k-ε 湍流模型对高雷诺数湍流的模拟经过了近一个世纪的实验与数值模拟等实践工作的反复检验,具有较好的通用性,在雷诺数相对较高的主流区,往往能得到较好的计算结果。故而,在远离壁面的高雷诺数区域,低雷诺数k-ε 湍流模型的方程形式应该还原成标准k-ε 湍流模型[9]。也有少数研究者认为,对于不同形态的流动,模型常数具有其特殊性,需根据不同流动形态进行修改[10]。

不同计算模型所得的结果,往往是这些阻尼函数、源项以及常数的综合作用。

本文应用标准k-ε 湍流模型及其6种低雷诺数模型,共7种模型,对三维壁面剪切流动进行了数值模拟。采用FLUENT14.0作为解算器。验证数据采用文献[18]试验数据。

流动域由矩形槽截面和光滑壁面组成,进口处被一平板分隔成两个进口(进口1：u = 2m/s,进口2：umax = 8.4m/s),平板距壁面b = 5mm,工作介质为空气(ρ = 1.225kg/m3,μ = 1.7894 × 10-5kg/(m•s))。流道长410mm,宽150mm,高100mm,取样线l1-8 = 30mm位于对称面上,如图1所示,单位为mm。流道网格绘制如图2所示。高速贴壁入口的雷诺数Reslot = ubρ /μ = 2875。

图1 流道几何尺寸简图 Fig. 1 Diagram flow geometry size

图2 流道网格尺寸图(网格数361×301×147) Fig. 2 Diagram flow grid size (grid number 361×301×147)

对模型的第一层网格的y 值而言,不同文献有着不同建议值[19-20],但本文根据从低雷诺数模型的物理原则出发,在黏性层流壁面的黏性底层(y ; 5)之内,都设置3—5层网格,根据各模型的最终收敛结果,壁面处第一层网格的y 在0.3~0.8之间。其中心截面的网格形态见图2,总网格数量为 1597万。计算过程中,对于压力与速度的耦合采用simple算法[21]。对方程离散采用二阶迎风差分格式,所有计算参数的收敛残值小于10-5。但为了保证收敛,达到收敛精度后每计算2000步,对取样线处前后计算结果的速度场进行比对,保证其波动变化值小于1/‰,方认为其收敛。每种模型算例的计算在曙光TC5000并行计算集群上的单个节点进行,每个节点的配置为双CPU(八核2.6G),计算时间约为4天(96h)。

图3为6种低雷诺数k-ε 湍流模型和标准k-ε 湍流模型的沿流向的速度计算值与实验数据之间的对比。其横坐标进行无量纲化,为纵坐标Y与高速射流进口2的间距之比值,坐标原点位于壁下面上。总体来看,对于标准k-ε 模型来说,由于其各向同性的假定和壁面函数的应用,对射流初期的湍流动能k与湍流耗散率 ε 有着较高预测,见图4、5,射流衰减明显高于其他模型。从速度场的计算准确度上来看,几乎所有低雷诺数湍流模型,相比标准k-ε 模型都有所改进。

在6种低雷诺数模型中,LS模型射流衰减最

图3 速度场与试验的对比 Fig. 3 Comparison of velocity field between the simulation results and the experimental result

图4 湍流动能 Fig. 4 Turbulent kinetic energy

快,在速度梯度较大的近壁面处,湍流动能k值与湍流耗散率 ε 值与其他模型相比明显过高,见图4与5。这与文献[14,22]曾分别采用Fluent12与Fluent13对该模型进行数值试验所得结论相似。

Patel根据量级分析发现,近壁面湍流的限制条

件拥有以下关系[16,23]：当y → 0,\(-\overline{{u}'{v}'}\propto {{y}^{3}}\),k  y2,

ν t  y3,ε → ε w,fμ → y-1。

根据表1中公式,进行同样的量级分析可知,对于低雷诺数k-ε 模型中最关键的阻尼函数fμ,在壁面附近当y → 0时,不同的低雷诺数模型fμ 具有不同的特性,即对于AB、AKN、CHC、YS模型,fμ → y-1;对于LB、LS,fμ → 常数。而在所有模型中,当远离壁面时,fμ → 1.0。

根据以上分析,LB、LS模型都不符合壁面限定条件。对此,众多研究者曾有过争议,一些研究者认为[16],在近壁面处由于湍流黏性对流动所起的作用较小,fμ 的量级大小对流动计算所起的作用很小,故而对流场的计算准确性影响不大。但近几十年以来,特别是对黏性层流底层外边界处湍流涡团猝发现象的认识,人们发现壁面限定条件对壁面传热、传质的计算准确性影响较大[11,24]。

而模拟所得的速度场的结果往往是阻尼函数、源项、常数综合作用的结果,仅仅通过算例中速度场模拟的准确性去判断模型的优劣是较为片面的,需结合湍流的物理本质及其流动现象,特别是湍流壁面的流动特征,对其他湍流参数的计算结果进行综合分析。

尽管AB与CHC模型符合fμ 的壁面限定条件,但在壁面速度梯度较大的射流初期(x/b = 1~20),在

图5 湍流耗散率 Fig. 5 Turbulent dissipation rate

黏性底层k值明显偏小,甚至趋向于零,见图4(每图的右上角为其近壁面处的放大图)。而根据现有的大量研究表明,贴壁湍流的湍流涡始发于边界层中层流底层的外边界处[25],这会使得此处存在一个湍流动能的峰值。AB与CHC模型在壁面附近的k值计算结果与实际流动物理特性并不相符。其余模型的计算结果均符合这一湍流壁面流动的物理特征,在壁面附近出现明显的小峰,见图4。而另一大峰为不同流速流体相遇后产生的强剪切作用所引起的。

同时,由于近壁面湍流局部平衡条件的限定,AB与CHC模型近壁处的湍流耗散率 ε 也偏低,见图5。对贴壁处湍流动能k以及相应的湍流耗散率 ε 的低估,也使得贴壁射流在后期衰减与实验值相比较慢,见图3。

由图3可知,LB、YS、AKN模型对速度场的模拟结果相对较准确,其中LB模型的准确性总体来说相对较高。但由于LB模型不符合壁面限定行为条件,因此无法预测高普朗特数和施密特数流动的传热问题[11]。Nagano和Tagawa曾在1990年提出了NT模型,使得其模型符合壁面限定条件,得到了较好的计算结果[26-27],但由于NT模型方程中包含了摩擦速度uτ,从而使得在uτ = 0处,雷诺应力消失,出现奇点。对此,K. Abe、T. Kondoh与Y. Nagano等人对NT模型提出了修改,采用了Kolmogorov速度替代摩擦速度,来计算近壁面区域以及低雷诺数效应,提出了AKN模型[10,28]。同时采用二维模型模拟结果结合实验数据,对常数Cε 1、Cε 2、Cμ、σ k、σε 中Cε 1、Cε 2常数以及湍流普朗特数 σ k、σε 进行了修改,使得其二维计算结果与试验数据更相符。

而本文在对三维附壁剪切流所进行模拟的基础上,将标准k-ε 模型的常数(见表2)用于AKN模型,进行了重新计算,计算结果见图6—8中的

图6 AKN1、AKN、YS与LB模型的速度场 Fig. 6 Velocity field of AKN1, AKN, YS and LB model

图7 AKN1、AKN、YS与LB模型的湍流动能 Fig. 7 k of AKN1, AKN, YS and LB model

AKN1图标所示。其结果表明,在采用标准k-ε 模型的模型常数的计算结果AKN1中,其速度场的准确度明显优于采用原常数的AKN模型与YS模型,也略优于LB模型,见图6。事实证明,文献[10,28]中采用二维计算得到的结果来修正模型常数会带来一定的误差。从另一方面也证明,经过了长期检验的标准k-ε 湍流模型的常数具有较好的通用性和客观性。由图7、8可知,AKN1的湍流动能k与湍流耗散率 ε 的结果与LB相似。在射流初期(X/b = 1,X/b = 5)LB模型的湍流动能k相对较高,这使得此处的射流衰减略微快于AKN1。而YS模型由于受到阻尼函数f1的影响,使得其近壁面处的湍流动能k与湍流耗散率 ε 较大。但在剪切区域,结果显示采用YS模型使得总体的湍流动能偏小,见图7,使得其射流衰减较慢。

1）大部分低雷诺数k-ε 模型相对于标准k-ε 模型来说,在对三维壁面剪切流的模拟中,其计算准确性都有所改进;但其中,LB、LS模型不符合壁面限定条件,而AB与CHC模型的k值计算结果与实际湍流壁面流动的流动特性不符。

2）在对速度场的模拟中,LB、YS、AKN模型准确度较高。其中YS与AKN模型符合壁面限定条件,且计算结果基本符合实际湍流壁面流动的流动特性。

3）在对三维附壁剪切流的模拟过程中,采用标准k-ε 模型所用常数的AKN模型的计算结果最为准确,且其模型也相对较为合理。

图8 AKN1、AKN、YS与LB模型的湍流耗散率 Fig. 8 ε of AKN1, AKN, YS and LB model

[1] Launder B E,Spalding D B.The numerical computation of turbulent flows[J].Computer Methods in Applied Mechanics and Engineering,1974,3(2):269-289.

[2] 许昌,韩星星,王欣,等.基于改进致动盘和拓展k-ε 湍流模型的风力机尾流数值研究[J].中国电机工程学报,2015,35(8):1954-1961. Xu Chang,Han Xingxing,Wang Xin,et al.Study of wind turbine wake modeling based on a modified actuator disk model and extended k-ε turbulence model[J].Proceedings of the CSEE,2015,35(8):1954-1961(in Chinese).

[3] 梁健,舒立春,胡琴,等.风力机叶片雨淞覆冰的三维数值模拟及试验研究[J].中国电机工程学报,2017,37(15):4430-4436. Liang Jian,Shu Lichun,Hu Qin,et al.3-D numerical simulation and experiment on glaze ice accretion of wind turbine blade[J].Proceedings of the CSEE,2017,37(15):4430-4436(in Chinese).

[4] 韩志敏,刘坐东,王景涛,等.结构对圆形楞涡流发生器CaSO4污垢特性和综合换热特性的影响[J].中国电机工程学报,2017,37(10):2920-2926. Han Zhimin,Liu Zuodong,Wang Jingtao,et al.Effect of structure on the circular vortex generators in CaSO4 fouling characteristics and comprehensive heat transfer characteristics[J].Proceedings of the CSEE,2017,37(10):2920-2926(in Chinese).

[5] Chmielewski M,Gieras M.Three-zonal wall function for k-ε turbulence models[J].Computational Methods in Science and Technology,2013,19(2):107-114.

[6] Abid R.Evaluation of two-equation turbulence models for predicting transitional flows[J].International Journal of Engineering Science,1993,31(6):831-840.

[7] Lam C K G,Bremhorst K.A modified form of the k-ε model for predicting wall turbulence[J].Journal of Fluids Engineering,1981,103(3):456-460.

[8] Launder B E,Sharma B I.Application of the energy-dissipation model of turbulence to the calculation of flow near a spinning disc[J].Letters in Heat and Mass Transfer,1974,1(2):131-137.

[9] Yang Z,Shih T H.New time scale based k-epsilon model for near-wall turbulence[J].AIAA Journal,1993,31(7):1191-1198.

[10] Abe K,Kondoh T,Nagano Y.A new turbulence model for predicting fluid flow and heat transfer in separating and reattaching flows-I.Flow field calculations[J].International Journal of Heat and Mass Transfer,1994,37(1):139-151.

[11] Chang K C,Hsieh W D,Chen C S.A modified low-Reynolds-number turbulence model applicable to recirculating flow in pipe expansion[J].Journal of Fluids Engineering,1995,117(3):417-423.

[12] 陶文铨. 数值传热学[M].2版.西安:西安交通大学出版社,2001:362-370. Tao Wenquan.Numerical heat transfer[M].2nd ed.Xi’an:Xi’an Jiaotong University of Press,2001:362-370(in Chinese).

[13] 张兆顺,崔桂香,许春晓.湍流理论与模拟[M].北京:清华大学出版社,2005:219-222. Zhang Zhaoshun,Cui Guixiang,Xu Chunxiao.Theory and modeling of turbulence[M].Beijing:Tsinghua University of Press,2005:219-222(in Chinese).

[14] Gorji S,Seddighi M,Ariyaratne C,et al.A comparative study of turbulence models in a transient channel flow[J].Computers & Fluids,2014,89:111-123.

[15] 范维澄,万跃鹏.流动及燃烧的计算模型[M].合肥:中国科学技术大学出版社,1992:124-125. Fan Weicheng,Wan Yuepeng.Calculation models of flow and combustion[M].Hefei:University of Science & Technology China Press,1992:124-125(in Chinese).

[16] Patel V C,Rodi W,Scheuerer G.Turbulence models for near-wall and low Reynolds number flows-a review[J].AIAA Journal,1985,23(9):1308-1319.

[17] Jones W P,Launder B E.The prediction of laminarization with a two-equation model of turbulence[J].International Journal of Heat and Mass Transfer,1972,15(2):301-314.

[18] Sato Y,Hishida K,Maeda M,et al.Effect of dispersed phase on modification of turbulent flow in a wall jet[J].Journal of Fluids Engineering,1996,118(2):307-315.

[19] Jagadeesh P,Murali K.Application of low-Re turbulence models for flow simulations past underwater vehicle hull forms[J].Journal of Naval Architecture and Marine Engineering,2005,2(1):41-54.

[20] Rathore S K,Das M K.Comparison of two low-Reynolds number turbulence models for fluid flow study of wall bounded jets[J].International Journal of Heat and Mass Transfer,2013,61:365-380.

[21] Patankar S V,Heat N.Numerical heat transfer and fluid flow[M].New York:Hemisphere,1980:146-153.

[22] Mathur A,He S.Performance and implementation of the Launder-Sharma low-Reynolds number turbulence model[J].Computers & Fluids,2013,79:134-139.

[23] Hsieh W D,Chang K C.Calculation of wall heat transfer in pipe expansion turbulent flows[J].International Journal of Heat and Mass Transfer,1996,39(18):3813-3822.

[24] Nešić S,Postlethwaite J,Bergstrom D J.Calculation of wall-mass transfer rates in separated aqueous flow using a low Reynolds number k-ε model[J].International Journal of Heat and Mass Transfer,1992,35(8):1977-1985.

[25] Versteeg H K,Malalasekera W.An introduction to computational fluid dynamics-the finite volume method[M].Harlow,Essex:Longman Scientific & Technical,1995:400-401.

[26] Nagano Y,Tagawa M,Niimi M,et al.An Improvement of the k-ε turbulence model:the limiting behavior of wall and free turbulence,and the effect of adverse pressure gradient[J].Transactions of the Japan Society of Mechanical Engineers Series C,1989,55(512):1008-1015.

[27] Nagano Y,Tagawa M.An improved k-epsilon model for boundary layer flows[J].Journal of Fluids Engineering,1990,112(1):33-39.

[28] Abe K,Kondoh T,Nagano Y.A new turbulence model for predicting fluid flow and heat transfer in separating and reattaching flows-II.Thermal field calculations[J].International Journal of Heat and Mass Transfer,1995,38(8):1467-1481.