[摘要]采用带预处理情形下的大涡模拟方法，结合高阶WENO格式，针对椭圆斜管湍流横向射流进行了数值模拟。计算结果显示在选取截面上，三个方向上平均速度分布与试验值较为接近。对计算结果进行分析，描述流场中发卡涡拟序结构及其运动演化过程，阐述发卡涡对流向速度分布的影响，讨论流场中雷诺应力与发卡涡运动演化之间的关系。
[关键词]预处理 斜管湍流横向射流 大涡模拟 发卡涡拟序结构
中图分类号：TB1文献标识码：A文章编号：1671－7597（2009）0420111－02
一、引言
射流在工程中应用非常广泛，而大多数射流都是湍流的。了解湍流射流中大尺度拟序结构的产生、发展和破碎等动力学演化过程，对有效地控制湍流射流将起到很重要的作用。大涡模拟方法（LES）是对流场中大尺度运动结构进行直接模拟，而对小尺度运动结构采用模型近似。与常用的雷诺平均计算方法（RANS）相比，大涡模拟方法能够更好地捕捉流场中的主要流动结构。通常在低马赫数流动并且伴有传热情况下，流场中密度是变化的，因而，可压缩情形下的控制方程是不适用的[2]。同时大多数的可压缩流动计算格式在解决低速流动问题时，其稳定性和收敛性将变得恶化。采用预处理方法能够比较好地解决这一问题；此外，在运用预处理技术后，能够对流场中的不同流动区域（可压缩流域和不可压缩流域）采用同一个程序代码进行计算，在提高格式的计算精度、稳定性和收敛性的同时，避免了针对不同的流动状况采用不同计算方法的复杂性，使得问题大大简化。
二、控制方程
程序采用预处理情形下[1,2]的大涡模拟计算方法，首先对Navier-Stokes方程进行形式的过滤，得到曲线坐标系下，无量纲通量形式的控制方程为：
三、预处理技术
为了改善低马赫数流动情况下计算格式的稳定性和计算结果的收敛性，本文采用了文献[1,2]中的预处理技术。预处理情形下控制方程的通量形式改写为：
方程进一步可以改写为：（3）
四、亚格子模型
计算中采用Smagorinsky模型作为亚格子应力计算模型。具体如下：
算网格尺度，本文中取C为常数0.01。在近壁区域内，Smagorinsky模型不具备正确的渐进行为，粘性系数不趋近于零而为一有限值，因此在亚格子模型中引入近壁面修正，修正函数形式为：
（5）
其中 为亚格子模型中修正后的常数。
五、计算物理模型
本文选取的算例为椭圆斜管湍流横向射流，射流管与水平方向的夹角为35度[4]，计算域网格如图1所示：
主流区域流向长度为9.2D，法向高度约为1.6D，展向宽度为3.2D，射流管高度为0.8D，其中D为射流管直径，原点取为射流管出口与主流椭圆交界面的中心，总网格量870000。流动特征马赫数为 ，吹气比约为1.0，基于主流进口速度和射流管直径的雷诺数 。空间上采用高阶WENO格式进行插值得到网格界面上的守恒变量，然后运用Roe格式求得界面上的无粘数值通量；时间上采用三步Runge-Kutta法伪时间步显示推进[5]，精度为三阶。
六、计算结果演示及分析
图2为处数值模拟计算结果与试验值[4]对比图。从图2上可以看出，数值模拟计算结果得到的平均速度分布曲线与试验结果比较接近，其差别可能是由于计算网格、进口真实湍流边界条件等方面的因素引起的，但是数值计算结果基本上能够反映出给定截面三个方向上平均速度的大致变化趋势。
（a）流向；（b）法向；（c）展向
图3显示的是主流流场中不同时刻的等涡面图，从图3（a）～（j）不同时间序列的流场等涡面图可以看出发卡涡的生成、发展以及破碎过程。如图3（a），流场中靠近的位置上出现了发卡涡A，此时，A正处于生成阶段，发卡涡的结构尚不清晰；图3（b）中发卡涡A已经移动到
的位置上，此时可以比较清晰地看出发卡涡A的涡头和涡腿；图3（c），A移动至处，可以看出上游有一个新的发卡涡生成；图3（d），A处于的位置，并且逐渐从上游新生成的发卡涡中脱离出来，此时的发卡涡已经渐渐发展得比较完整，可以明显地分辨出典型的发卡涡结构特征；而在图3（f）时刻，A呈现出破碎的趋势，此时A逐渐步入其发展的末期；图3（g），A的位置为 ，此时A接近完全破碎的状态。至此，发卡涡A经历了从生成，发展，直至破碎的三个阶段。可以看出，在此刻A的下游前一时刻的发卡涡已经发生了破碎，流场呈现出完全湍流的流动状态；而在A的上游不断伴随有新的发卡涡生成，流场中形成的发卡涡都在重复着上述的三个过程。同时，注意到在流场中有部分的发卡涡在流动过程中较早地发生了破碎现象。
图3（a）～（j）主流流场不同时刻等涡面图
图4（a）显示的为y=0.315D截面上流向速度等值线图，可以看出，在靠近中间区域内流向速度值相对较小，图示最小值为0.48Ue，产生这种低速条带的原因是流体在绕着发卡涡涡腿旋转的同时，会对近壁处的主流低速流体产生卷吸作用，使其离开壁面，这就使得垂直涡腿之间的流体流向速度相对较低，如图4（a）和图5（b）。而在远下游的位置上，发卡涡扭曲并相继发生破碎，主流流体与射流管流体相互掺混作用程度较强，总体来讲流向速度将会有所提高，但仍然会形成如图4（c）所示的低速区域。以上原因也有助于解释图2（a）中流向速度分布曲线图上存在的速度亏损现象。
图4（a）～（c） 不同位置上流向速度等值线图：
（a）y=0.315D；（b）x=1.522D；（c）x=2.534D。
（c）x=2.939D（d）x=3.950D
（e）x=4.961D（f） x=5.972D
图5 （a）～（f）主流流场流向不同位置处雷诺应力等值线图
图5（a）～（f）显示的为主流流场流向不同截面处雷诺应力等值线图。从图上可以看出靠近射流管出口附近应力值较大，对应此处射流流体与主流流体相互作用程度较为剧烈，如图5（a）。而往下游雷诺应力值变得相对较小，并且在粘性扩散的作用下，应力影响区域有向周围两边扩展的趋势，其影响范围相对较大，如图5（e）和图5（f）。结合图3（a）～（j）等涡量图分析，在靠近射流管出口附近，新生成的发卡涡的强度较大；而在远下游，之前形成的发卡涡结构发生破碎进而形成完全湍流的状态，涡的相对强度变弱，而破碎后其影响区域变大，这一过程与图5（a）～（f）中各图是相吻合的。
七、结论
1．本文主要工作在于将大涡模拟计算方法与预处理技术和高阶WENO格式相结合，针对椭圆斜管湍流横向射流进行了数值模拟。在给定界面上，计算结果所得到的流向、法向和展向平均速度的分布曲线与试验值较为吻合，并且观察到流场中出现了发卡涡拟序结构；
2．针对发卡涡运动演化过程作了简要地描述和分析。计算结果较清晰地反映出了发卡涡生成、发展并最终破碎的运动演化过程。发卡涡拟序结构的存在影响了主流流向速度分布，在流场中形成了低速条带结构；
3．发卡涡的发展与雷诺应力之间存在着相互关系。在发卡涡新生成阶段雷诺应力较强，且影响区域较小；而在发卡涡逐渐发展直至破碎的过程中，雷诺应力变得相对较弱，但其影响区域却变大。
参考文献：
[1]Xiao Wang. A preconditioning Algorithm for Turbo machinery Viscous Flow Simulation，Doctor Degree Dissertation of Philosophy in computational Engineering in the Engineering Research Center of Mississippi State University，December 2005.
[2]W.R.Briley，L.K.Taylor，D.L.Whitfield. High-resolution viscous flow simulations at arbitrary Mach number，Journal of Computational Physics 184(2003）79-105.
[3]Chi-Wang Shu，Essentially Non-Oscillatory and Weighted Essentially Non-Oscillatory Schemes for Hyperbolic Conservation Laws，NASA/CR-97-206253, ICASE Report No.97-65.
[4]Mayank Tyagi，Sumanta Acharya. Large Eddy Simulation of Film Cooling Flow From an Inclined Cylindrical Jet，ASME Vol.125，October 2003，Page734-742.
[5]N.Alkishriwi,M.Meinke，W.Schroder，A large-eddy simulation method for low Mach number flows using preconditioning and multigrid.Computers&Fluids 35(2006）Page 1126-1136.
[6]Chunhua Sheng and Xiao Wang.Characteristic Variables Boundary Conditions for arbitrary Mach number algorithm frame，AIAA 2003.
推荐访问： 湍流 射流 预处理 横向 数值
