传统的雷诺平均模拟(RANS)方法已经在工程领域得到广泛应用，然而大量的研究表明RANS方法对大分离流动的模拟能力明显不足。受有限的计算资源所影响，直接数值模拟(DNS)和大涡模拟(LES)目前还无法在工程领域广泛应用，因此计算效率相对较高的RANS/LES混合方法引起了研究者们的广泛关注。

在众多RANS/LES混合方法中，DES方法以其简单的构造和较广的适用范围被众多研究者们所关注。Spalart等^[1]在1997年提出了基于一方程模型的脱体涡模拟(DES)方法，这一方法也是后来众多DES类改进方法^[2-3]的雏形。理想状态下RANS/LES方法对附着流动采用RANS模式，对分离流动则采用LES模式。然而实际应用中RANS和LES区域之间的转换并不是立即完成的，存在一个既不是纯粹RANS也不是纯粹LES的过渡区域，这就是所谓的“灰区”现象。灰区问题会导致一些非物理的流动现象，比如模型应力耗散(MSD)现象^[4]，它会引起网格诱导分离(GIS)^[5]和对数层不匹配(LLM)^[6]等问题。

为了解决MSD引起的GIS问题，Spalart提出了改进版本的延迟DES方法(DDES)^[2]。Shur等^[3]在此基础上又提出了改进的DDES方法(IDDES)，其目的是进一步解决模型应力耗散引起的对数层不匹配的问题。另一种思路是分区域的DES方法^[7]，这种方法的思想是分别指定流场中的RANS和DES计算区域，通过预先设定交界面来避免MSD和GIS问题。

作为经典的大分离流动算例，对后台阶流动的数值模拟能够很好地评估各种湍流模拟方法对流动分离的预测能力。原始的DES方法是基于一方程SA模型所构造的，本文采用基于两方程SST湍流模型的多种DES类混合方法对后台阶分离流动进行了数值计算研究，通过对计算结果的分析对比，验证各方法的数值模拟能力，从而为利用DES类方法预测工程实际中的复杂流动分离问题提供参考。

1 数值方法及湍流模型

本文工作是基于中心有限体积方法，时间离散采用二阶全隐式LU-SGS-τTS算法，空间离散采用基于五阶WENO插值的Roe格式。采用两方程k-ω SST模型作为RANS/LES混合方法的基础湍流模型。

1.1 基于k-ω SST模型的DES类混合方法

原始的SST模型^[8]对近壁区流动采用Wilcox提出的k-ω模型，而在边界层边缘及自由剪切层则采用k-ε模型，两者通过混合函数进行过渡。为了构造DES类混合方法，这里需引入混合长度尺度l, 则湍动能方程可以写成如下形式

(1)

式中，ρ为密度, k为湍动能, σ_k=0.5, μ为分子运动黏度, 湍流涡黏性, τ_ij是雷诺应力张量。F₂为模型中的一个混合函数

(2)

混合长度尺度l根据所需的DES方法可以取不同的值, 对于原始的SST模型。

Menter等^[5]在2003年提出了一种基于SST模型中混合函数的DDES方法(以下简称DDES-2003), 混合长度尺度的构造如下

(3)

式中，C_DES=0.78·F₁+0.61·(1-F₁), 亚格子尺度Δ_max=max(Δx, Δy, Δz), Δx, Δy, Δz分别为x、y、z 3个方向上的网格步长。混合函数F₁为

(4)

这里交叉扩散项, d为到当地最近壁面的距离。混合函数F_SST可以取不同的值:当F_SST=0时(3)式即为SST模型; 当取混合函数F₁或F₂时(3)式为DDES模型, 本文中F_SST=F₂。

Spalart等^[2]在2006年提出了基于一方程模型的DDES方法(以下简称DDES-2006), 为将该方法也推广到两方程SST模型中, Gritskevich等^[9]对模型中参数进行了修正, 混合长度尺度的构造如下:

(5)

f_d为延迟函数

(6)

(7)

式中，v_t为湍流黏性, κ为Karman常数取值为0.41, S为应变变化率张量, Ω为旋度张量。

当采用IDDES方法^{[2, 9]}时, 混合长度l取值为

(8)

式中亚格子长度Δ=min[C_wmax(d, Δ_max), Δ_max], d是网格单元到壁面的距离, C_w=0.15。l_SST、C_DES和Δ_max的定义与(3)式中一致。经验混合函数的定义如下:

(9)

当f_e=0时, IDDES以DDES模式运行; 当f_e> 0且时, IDDES方法在近壁区以WMLES模式运行。函数f_e的具体构成如下:

(10)

(9) 式和(10)式中的各模型常数取值如下:

1.2 基于五阶WENO插值的Roe格式

对无黏通量项离散采用基于迎风型通量差分裂的Roe格式。在三维贴体正交坐标系下穿过网格单元界面的无黏通量可表示为:

(11)

采用Roe格式对N-S方程空间项进行离散求解时, 首先求得网格单元体交界面两侧变量的值q_L和q_R。对网格交界面左右两侧的变量通过插值得到。对交界面处变量插值的精度就决定了整个空间离散格式的精度。本文采用基于改进的WENO方法对交界面两边状态参量进行重构的WENO-Roe五阶格式^[10-11]。这种空间离散格式的精度已经在之前对超声速底部流动^[12]的研究中得到验证。

以一维标量模型方程为例, 五阶精度的边界外推值形式为:

(12)

式中，通过第k(k=0, 1, 2)个模板S_k=(i-k, i-k+1, i-k+2)中的单元平均值外推得到。即:

(13)

ω_k为非线性加权因子, 其定义为

(14)

其中系数, , 为线性因子, ε=10^-6其作用是避免分母为0, 光滑因子β_k的表达式为:

(15)

2 计算设置和网格

数值计算条件是基于Driver和Seegmiller的实验^[13]展开的。自由来流速度为44.2 m/s, 台阶高度为H=12.7 mm, 温度为289 K, 基于台阶高度的雷诺数为Re=36 000。更多实验细节可参见文献[13]。

为了在使台阶上游得到充分发展的湍流边界层, 计算中在选取了一个较长的上游长度(110H)。展向长度取2H, 计算域尺寸如图 1所示。入口为自由来流边界, 出口为插值边界, 展向两侧面为周期性边界条件, 上下表面为绝热无滑移壁面。

计算网格采用结构化多块网格, 台阶上游区域流向取161个网格点, 下游区域流向网格点数为320, 展向网格点数为40, 首层网格尺度为5×10^-6mm, 以保证壁面首层网格y⁺≤1。网格单元总数为2.98×10⁶。网格截面如图 2所示。

非定常计算时物理时间步长为Δt=3×10^-6 s, 对应的无量纲时间步(基于声速和台阶高度H)为0.01, 每个时间步内子迭代步数为5。计算时首先执行定常RANS计算, 将其作为DES计算的初场。为保证平均流场的准确性, 待流场稳定之后, 取最后2×10^-4时间步的计算结果进行统计平均。

3 计算结果分析

在采用DES类方法对后台阶流动进行模拟时，在台阶上游的附体流动区域和通道上壁面采用RANS模式，经过台阶肩部位置后，DDES方法在下游分离区域快速转换成LES模式，而IDDES继承了更多上游的湍流信息，在再附区域的边界层执行WMLES模式。

图 3为3种DES方法计算得到的瞬时湍流结构，图中所展示的是依据Q准则^[14]的等值面图。后台阶分离涡结构的发展是一种典型的Kelvin-Helmholtz不稳定性的表现。剪切层涡结构从台阶肩部脱出形成二维的涡结构，在向下游发展过程中，这些二维的涡逐渐演化成大尺度的三维涡结构。从图中可以看出在再附区域内包含了一些大尺度的发卡涡。由于IDDES方法在下游边界层执行WMLES模式，使其捕捉到了更加丰富的涡结构，尤其是在下游的远场区域，IDDES展示了比其他2种DDES方法更多的小尺度湍流结构。

图 4对比了几种方法计算得到的平均流线图。与2种DDES方法和IDDES方法的计算结果相似，非定常RANS(URANS)计算得到的平均流线图也展示了台阶下游的主分离区域和台阶下方拐角处的二次分离区域，所不同的是URANS所预测的主分离区规模明显要小。

对下游主分离区域再附位置的预测是衡量计算结果好坏的重要指标。图 5比较了下游壁面摩擦力系数分布，摩擦力等于零的位置即为再附位置。在台阶下游分离区内，IDDES和DDES-2006预测的摩擦力系数与实验吻合较好，而URANS结果与实验差距较大。表面摩擦力的预测对网格变化较为敏感，在下游x/H=10以后，由于网格尺度逐渐变大，使得计算结果与实验值有一些差异。计算得到的再附位置见表 1。可以看出，IDDES方法预测的分离再附位置与实验测量结果较好，相比之下2种DDES方法得到的再附位置略靠近上游，而URANS结果与实验值差异较大，这是由于RANS方法不能够很好地捕捉后台阶引起了非定常涡脱落过程导致的，这与平均流线图的比较结果也是一致的。

图 6分别展示了x/H=1、x/H=4、x/H=6、和x/H=10流线位置计算得到的速度分布剖面分布和实验的对比，无量纲时的参考速度为上游x/H=-4位置中线位置处的速度值。可以看出，2种DDES方法和IDDES方法预测的速度分布没有明显的差异，与实验结果吻合较好。由于非定常RANS所预测的再附区域与实验存在差异，使得计算得到的各流向位置速度分布与实验结果差异较大。

高阶量的统计结果也是评价数值模拟精度的重要参考指标，图 7中对4个不同流向位置的流向雷诺应力分布进行了对比。URANS得到的流向雷诺应力与实验值的吻合性较差，而IDDES和DDES计算结果与实验吻合较好，即使在远场区域(x/H=10)依然能够得到较好的预测结果，这也表明DES类混合方法对大分离流动的预测具有一定优势。与速度分布的对比结果相似，IDDES和其他2种DDES方法预测的雷诺应力没有较大差异。

4 结论

本文采用基于两方程k-ω SST湍流模型的IDDES方法和2种不同的DDES方法，结合基于五阶WNEO插值的Roe格式，针对后台阶大分离流动进行了数值模拟，将计算结果与实验值进行对比，结论如下：

1) IDDES和其他2种DDES方法能够很好的模拟后台阶涡脱落的发展过程；由于在下游再附区域附面层采用WMLES模式，IDDES方法在分离区域能够捕捉到更加丰富的涡结构。

2) URANS预测的下游再附区域偏小，再附位置与实验差异较大；2种DDES方法得到的下游再附位置略靠近上游，IDDES方法预测的再附位置与实验吻合较好。

3) DDES和IDDES方法预测的不同流向位置的速度和雷诺应力分布没有明显的差异；URANS结果与实验值吻合较差。

综合上述计算结果，IDDES和DDES方法对后台阶分离流动具有较好的预测效果，结合工程实际问题，将DES类方法应用到更多复杂的分离流动问题中去将是作者下一步研究工作的方向。