可压缩槽道流动的大涡模拟研究

王锁柱，杨天鹏，杨依峰，秦绪国，苏 伟

可压缩槽道流动的大涡模拟研究

王锁柱，杨天鹏，杨依峰，秦绪国，苏 伟

（北京航天长征飞行器研究所，北京，100076）

可压缩湍流边界层是高速飞行器研究面临的基础科学问题，对其流动机理的认识对于高速飞行器气动力热设计具有重要意义。采用大涡模拟方法对可压缩槽道流动进行了数值模拟研究，并基于Favré过滤建立了适用于可压缩流动的大涡模拟方法，将混合亚格子模型延拓到了动态形式，形成了适用于可压缩流动的动态混合模型。通过分析亚格子模型、网格疏密和计算格式对数值模拟结果的影响，结果表明所建立的大涡模拟方法在较少网格数下能够获得可靠的计算结果。数值模拟还获得了可压缩槽道流动时间转捩的过程，捕捉到了大尺度旋涡结构和近壁速度条带结构的演化过程。

可压缩湍流；槽道流动；大涡模拟

0 引 言

可压缩湍流现象广泛存在于高速飞行器的边界层流动中，具有多尺度、高频脉动、强非线性等特点，对其流动机理的研究对于高速飞行器设计具有重要意义。随着计算流体力学（Computational Fluid Dynamics, CFD）和计算机技术的迅速发展，数值模拟方法已经成为可压缩湍流问题的重要研究手段。

大涡模拟（Large Eddy Simulation，LES）是综合考虑了工程应用的要求以及计算能力的局限而发展起来的一种介于雷诺平均方法（Reynolds Averaged Navier-Stokes，RANS）和直接数值模拟（Direct Numerical Simulation，DNS）之间的数值方法。一方面，LES比RANS易于建立更为普适通用的湍流模型，能够获得真实的瞬态流场信息，可用于预测分离、转捩等复杂流动。另一方面，有效的亚格子模型可使LES所需网格总数比DNS少一到两个数量级，从而其计算量要明显少于DNS。亚格子模型是大涡模拟的核心问题之一，对于大涡模拟的准确性和计算效率起关键作用。Smagorinsky模型[1]是最早的亚格子模型，其简单易行，计算稳定性和鲁棒性也较好，但包含非普适的经验常数，同时耗散性过大。随后发展的尺度相似模型[2]能够较准确地表达大尺度和小尺度间的动量输运关系，其缺陷是耗散严重不足，数值计算存在稳定性问题。综合Smagorinsky模型和尺度相似模型的优点可形成混合模型[3]，既有正确的亚格子动量输运，也有足够的亚格子耗散。Germano等[4]提出了著名的动态模型，该模型本身不提出新的模型，而是建立在给定的基准模型上，通过动态方法确定模型中的系数，如动态Smagorinsky模型、动态混合模型[5]等。

可压缩槽道流动作为一种可压缩壁湍流，不仅包含了壁面效应，还可体现可压缩效应的影响，是一种典型的可压缩湍流流动。Gamet等[6]对马赫数为0.2的槽道流动进行了DNS研究，由于马赫数较低，流场特征与不可压湍流非常接近。Lenormand等[7]采用LES对马赫数为0.5和1.5的槽道湍流进行了模拟，研究了不同亚格子模型对计算结果的影响。Mossi和Sagaut[8]也对马赫数为0.5和1.5的槽道湍流进行了模拟，分析了不同激波捕获格式对LES结果的影响。李新亮等[9]采用DNS研究了马赫数为0.8的可压缩槽道湍流场，分析了压缩性效应对近壁相干结构的影响。

本文采用大涡模拟方法对马赫数为0.5的可压缩槽道流动进行了数值模拟研究，分析了亚格子模型、网格疏密和计算格式对大涡模拟结果的影响，并给出了可压缩槽道流动时间转捩过程中大尺度流向涡和近壁速度条带结构的演化过程。

1 物理模型和数值方法

1.1 物理模型

槽道流动示意如图1所示，两个平行平板间充满粘性可压缩流体沿流向运动。设笛卡尔坐标系的坐标轴、和方向分别平行于流向、法向和展向，为槽道半高。本文取槽道的半高=1，=1对应槽道的中心线。

图1 槽道流动示意

在不可压缩槽道流动的数值模拟中，为了保证流向的周期性，通常假设流动是由均匀压力梯度驱动的，即以流向压差平衡上下壁面的摩擦。但在可压缩槽道流动的模拟中，若仍假设由平均压力梯度驱动，则流场中必然诱导出密度和温度梯度，进而影响整个流场参变量的分布，导致流向的周期条件无法维持。本文采用Coleman等[10]提出的可压缩槽道流动的驱动模型，认为流体是在均匀体积力的作用下流动，即通过在流向动量方程中添加体积力项1对壁面摩擦力进行平衡，以保证流场变量满足流向的周期性条件，同时在能量方程也需要引入一项11。为避免计算体积力时出现数值不稳定，本文采用Lenormand等[11]发展的可压缩流动修正算法来更新每一时间步的体积力。

1.2 大涡模拟方法

引入Favré过滤后，除密度和压强仍采用直接空间过滤表示外，其余流场变量，如速度场、温度等均取Favré过滤量，则可导出基于Favré过滤后的N-S方程组为

其中，模型系数、I可基于瞬态流场动态计算得到。

采用有限差分法对大涡模拟控制方程进行离散求解。对流项经通量分裂后采用五阶迎风型紧致格式[12]，粘性项则采用六阶中心型紧致格式[13]，时间离散采用满足TVD特性的三阶Runge-Kutta方法。

1.3 初始条件和边界条件

对于充分发展的槽道湍流，流向和展向是均匀同性的。因此可以沿流向和展向采用周期边界条件，法向上下壁面采用等温无滑移条件，即=w和w=0。

槽道初始的速度场、温度场按二维层流Poiseuille流动的解析解给出，密度场则设为均匀分布。为了诱导初始层流流动顺利转捩为充分发展的湍流流动，在流向速度中需要加入一定幅值的随机扰动。具体初始条件如下：

2 结果分析

本文对马赫数为0.5的亚声速槽道流动的时间转捩过程进行了大涡模拟研究。计算域沿流向、法向和展向的长度分别为2π、2和4π/3，基于体积平均密度、体积速度、槽道半高和壁面粘性系数的雷诺数为3000。计算网格在流向和展向均匀分布，在法向采用双曲拉伸对壁面网格加密。槽道流动在流向随机扰动的作用下发生了明显的时间转捩，从层流转变为完全湍流。本文对处于完全湍流状态下的流场进行统计平均，由于流向和展向的均匀性，统计除了按时间平均，也沿流向和展向进行空间平均。

2.1 亚格子模型影响

为了考察不同亚格子模型对计算结果的影响，本文采用4种不同亚格子模型分别对可压缩槽道流动进行了数值模拟。表1给出了不同亚格子模型计算得到的典型平均量的对比。表中f为壁面摩擦系数、τ为摩擦雷诺数、u为壁面摩擦速度、c为槽道中心处速度、b为体积速度。DMM代表本文所发展的动态混合模型，SM代表Smagorinsky模型，DSM代表动态Smagorinsky模型，HSSM代表混合模型。表中还给出了几组对比数据，NAH代表Niederschulte等[14]的实验测量结果，KMM代表Kim等[15]的不可压缩DNS数据，Dean为经验公式给出的结果[16]。可以发现，SM和DSM计算得到的壁面摩擦系数和壁面摩擦速度相对于DMM误差稍大，而HSSM的结果与KMM很接近。HSSM计算得到的槽道中心速度较为准确，其他3种模型的结果都稍大。

图2给出了采用不同亚格子模型计算得到的流向平均速度分布。在全局坐标轴下，4种模型得到的速度型相差不大，SM和DSM的速度型在靠近槽道中心时偏高，而HSSM的结果则偏低。在壁面坐标轴+下，4种模型在粘性底层得到的结果都较为一致，而在对数区SM和DSM的速度分布出现了上抬，因为该两种模型所计算得到的壁面摩擦速度稍小。

表1 不同亚格子模型计算的平均量对比

Tab.1 Mean Flow Variables for Different SGS Models

算例CfReτuτucuc/uτub/uτ NAH8.29×10-31790.06441.1517.8615.54 KMM8.18×10-31800.06371.1618.2015.63 Dean8.29×10-3——1.16—— DMM8.32×10-31960.06371.1818.5215.92 SM7.96×10-31910.06241.1818.9116.28 DSM7.76×10-31880.06181.1819.0916.53 HSSM8.19×10-31950.06281.1618.4715.94

图2 不同亚格子模型计算的流向平均速度

图3给出了采用不同亚格子模型计算得到的速度脉动均方根分布。由图3可以看出，SM和DSM得到的流向速度脉动均方根的峰值要比DMM和KMM的结果大，峰值位置也有所偏离；而HSSM的峰值偏小。对于法向速度脉动均方根，相对于DMM的结果，SM、DSM和HSSM得到的峰值依次减小，峰值位置也向槽道中心偏离。展向速度脉动均方根的分布与法向类似。可以看出，本文所建立的动态混合模型在总体上要优于其他3种模型，得到的数值模拟结果也更为准确。

图3 不同亚格子模型计算的速度脉动均方根

2.2 网格疏密影响

为了考察网格疏密对计算结果的影响，对网格加密后进行数值模拟，亚格子模型均采用动态混合模型。表2给出了网格特征参数对比，n、n、n分别表示流向、法向和展向的网格点数，Δx和Δz是壁面坐标系下的流向、展向网格间距，min(Δy)是壁面坐标系下法向离壁面第1层的网格间距，max(Δy)则是槽道中心的网格间距。=0.5为原始网格，=0.5为对法向网格进行加密，0.5为对展向网格进行加密。

表3给出了网格加密后计算得到的典型平均量对比。由表3可以看出，法向网格加密后得到的壁面摩擦系数和壁面摩擦速度略小于加密前的量，而展向网格加密后的结果则变大。法向或展向加密计算得到的槽道中心速度均为1.15，略小于加密前的结果。

表2 网格参数对比

Tab.2 Comparison of Grid Characteristics

算例nxnynzΔx+min(Δy+)max(Δy+)Δz+ KMM192160129120.054.47 Ma=0.541656530.751.0512.312.81 Ma=0.5-y411296530.770.516.212.82 Ma=0.5-z416512931.521.0812.66.57

表3 不同网格计算的平均量对比

Tab.3 Mean Flow Variables for Different Grids

算例CfReτuτucuc/uτub/uτ NAH8.29×10-31790.06441.1517.8615.54 KMM8.18×10-31800.06371.1618.2015.63 Dean8.29×10-3——1.16—— Ma=0.58.32×10-31960.06371.1818.5215.92 Ma=0.5-y8.18×10-31960.06271.1518.3415.91 Ma=0.5-z8.59×10-32000.06421.1517.9115.53

图4给出了网格加密后计算得到的流向平均速度分布。在全局坐标系下，网格加密后的速度型相差不大，只是在靠近槽道中心处偏低，与KMM的误差也很小。在壁面坐标系下，法向加密后的分布与加密前的结果较为一致，而展向加密后的速度型在对数区略有降低。

图4 不同网格计算的流向平均速度

图5给出了网格加密后计算得到的速度脉动均方根分布。由图5可以看出，法向加密后对流向速度脉动均方根的分布影响不大，而展向加密后分布的峰值有所减小，但峰值位置无变化。对于法向速度脉动均方根，法向加密后得到的分布基本不变，展向加密后在近壁面处与KMM符合得更好，而在靠近槽道中心时误差稍大。网格加密后展向速度脉动均方根的分布与法向相似。根据以上分析，综合考虑数值模拟精度和计算量的影响，本文采用的网格参数已可满足槽道流动大涡模拟的需求。

图5 不同网格计算的速度脉动均方根

2.3 计算格式影响

为了考察普通迎风型格式是否可用于进行大涡模拟计算，本文尝试采用目前较流行的几种迎风型格式开展了数值模拟研究，包括Harten-Yee TVD格式、AUSMDV格式和Roe格式。但在同样的网格参数下，即使将流向速度扰动幅值增大到0.5，流动也未能发生转捩，始终保持层流状态。显然，以上几种迎风型格式的数值耗散过大，不适用于湍流机理的大涡模拟研究。

2.4 流场结构分析

大涡模拟的一大优点在于能够获得大量瞬态流场，可以从中提取大尺度相干结构随时间或空间的演化过程，有利于对流动现象的深化认识。速度梯度张量的第二不变量是流向涡和展向涡的综合体现，本文采用法则来识别旋涡结构，即主要观察速度梯度张量第二不变量的等值面。

图6给出了计算过程中3个典型时刻=0.85的等值面分布。由图6可以看出，在计算的初始阶段（=5时）流场中几乎不存在大尺度相干结构；=50时，流场结构非常紊乱，存在大量不规则小涡结构；随着时间的进一步推移，开始出现大尺度旋涡结构，到=100时流场中已存在大量的准流向涡，这些准流向涡与壁面和流向均呈一定角度，并在向下游运动的过程中不断演化发展。

图6 不同时刻大尺度相干结构

图7给出了对应时刻瞬态流场=0.04法向截面的流向速度云图。由图7可知=5时流向速度较小，分布也较均匀；当=50时，流向速度出现大幅值的扰动，但未出现明显的展向速度条带结构；=100时，沿展向开始出现明显的高速区和低速区，并且呈条带状沿展向交错分布。这样的近壁速度条带结构是壁湍流重要的拟序结构。

图7 不同时刻法向截面流向速度云图

3 结束语

本文采用大涡模拟方法对可压缩槽道流动进行了数值模拟研究，分析了亚格子模型、网格疏密以及计算格式对数值模拟结果的影响。结果表明，本文采用的大涡模拟方法和发展的动态混合模型能够适用于可压缩湍流的数值模拟，采用较少的网格可以获得准确的结果。数值模拟还获得了可压缩槽道流动时间转捩的过程，捕捉到了大尺度旋涡结构和近壁速度条带结构的演化过程。

[1] Smagorinsky J. General circulation experiment with the primitive equations[J]. Monthly Weather Review, 1963, 91(3): 99-165.

[2] Bardina J, Ferziger J H, Reynolds W C. Improved subgrid scale models for large-eddy simulation[R]. AIAA Paper 80-1357, 1980.

[3] Meneveau C, Katz J. Scale-invariance and turbulence models for large eddy simulation[J]. Annual Review of Fluid Mechanics, 2000(32): 1-32.

[4] Germano M, et al. A dynamic subgrid-scale eddy viscosity model[J]. Physics of Fluids A, 1991, 3(7): 1760-1765.

[5] Zang Y, Street R L, Koseff J R. A dynamic mixed subgrid-scale model and its application to turbulent recirculating flows[J]. Physics of Fluids A, 1993, 5(12): 3186-3196.

[6] Gamet L, et al. Compact finite difference schemes on non-uniform meshes application to direct numerical simulations of compressible flows[J]. International Journal for Numerical Methods in Fluids, 1999(29): 159-191.

[7] Lenormand E, et al. Subgrid-scale models for large-eddy simulations of compressible wall bounded flows[J]. AIAA Journal, 2000, 38(8): 1340-1350.

[8] Mossi M, Sagaut P. Numerical investigation of fully developed channel flow using shock-capturing schemes[J]. Computers & Fluids, 2003(32): 249-274.

[9] 李新亮, 马延文, 傅德薰. 可压槽道湍流的直接数值模拟及标度律分析[J].中国科学(A辑), 2001, 31(2): 153-164.

Li Xinliang, Ma Yanwen, Fu Dexun. Direct numerical simulation of compressible channel turbulence and scaling law analysis[J]. Science in China (Series A), 2001, 31(2): 153-164.

[10] Coleman G N, Kim J, Moser R D. A numerical study of turbulent supersonic isothermal-wall channel flow[J]. Journal of Fluid Mechanics, 1995(305): 159-183.

[11] Lenormand E, Sagaut P, Ta Phuoc L. Large eddy simulation of subsonic and supersonic channel flow at moderate Reynolds number[J]. International Journal for Numerical Methods in Fluids, 2000(32): 369-406.

[12] 傅德薰, 马延文. 高精度差分格式及多尺度流场特性的数值模拟[J]. 空气动力学学报, 1998, 16(1): 24-35.

Fu Dexun, Ma Yanwen. High order accurate schemes and numerical simulation of multi-scale structures in complex flow fields[J]. Acta Aerodynamica Sinica, 1998, 16(1): 24-35.

[13] Lele S K. Compact finite difference schemes with spectral-like resolution[J]. Journal of Computational Physics, 1992(103): 16-42.

[14] Niederschulte M A, Adrian R J, Hanratty T J. Measurements of turbulent flow in a channel at low Reynolds numbers[J]. Experiments in Fluids, 1990, 9(4): 222-230.

[15] Kim J, Moin P, Moser R. Turbulence statistics in fully developed channel flow at low Reynolds number[J]. Journal of Fluid Mechanics, 1987(177): 133-166.

[16] Dean R B. Reynolds number dependence of skin friction and other bulk flow variables in two-dimensional rectangular duct flow[J]. Journal of Fluids Engineering, 1978, 100(2): 215-223.

Large Eddy Simulation of Compressible Channel Flow

Wang Suo-zhu, Yang Tian-peng, Yang Yi-feng, Qin Xu-guo, Su Wei

(Beijing Institute of Space Long March Vehicle, Beijing, 100076)

The compressible turbulent boundary layer is a basic scientific problem for the research of high-speed aircraft, and the understanding of its flow mechanism is of great significance for the aerodynamic and thermal design of high-speed aircraft. Large eddy simulation (LES) of compressible channel flow is performed. The LES method for compressible flow is established based on Favré filter, and a dynamic mixed subgrid-scale (SGS) model for compressible flow is formed by extending the mixed SGS model to the dynamic form. Through the analysis of SGS model, grid density and numerical scheme, the LES method established in the present work can obtain reliable results with fewer grids. The numerical results also show the transition process of the compressible channel flow, and capture the evolution process of the large-scale vortices and the near-wall velocity streaks.

compressible turbulence; channel flow; large eddy simulation

2097-1974(2023)01-0112-06

10.7654/j.issn.2097-1974.20230122

V211.3

A

2022-11-11；

2022-12-31

王锁柱（1984-），男，博士，高级工程师，主要研究方向为再入飞行器空气动力学设计。

杨天鹏（1993-），男，博士，工程师，主要研究方向为再入飞行器空气动力学设计。

杨依峰（1990-），男，高级工程师，主要研究方向为再入飞行器空气动力学设计。

秦绪国（1982-），男，博士，研究员，主要研究方向为再入飞行器总体设计。

苏 伟（1979-），男，博士，研究员，主要研究方向为再入飞行器空气动力学设计。