在直升机、固定翼等飞行器设计领域，大分离流动和旋涡主导流动是经常会遇到的流动现象，如飞机大迎角飞行时机翼后部形成分离涡、直升机旋翼脱出的尾迹涡等，这些脱体涡结构在形成、发展、破碎等过程中，存在着复杂的流动机理，在整个流场中起着主导作用。因而，开展湍流大分离流动问题研究，具有实际的工程应用价值，而在数值计算技术方面，建立高效、高精度的数值模拟技术，形成能够有效模拟高雷诺数大分离流动问题的方法，将是一项十分有意义的工作。

雷诺平均Navier-Stokes (Reynolds average Navier-Stokes, RANS) 方程耦合湍流模型 (如SA一方程模型、SST k-ω两方程湍流模型)，在处理高雷诺数湍流问题中得到了广泛且成功的应用，但是对于湍流大分离流动问题，该方法面临不小的技术障碍。目前，大涡模拟技术^[1](Large eddy simulation，LES) 被认为是处理存在涡结构流动问题的有效方法，但是由于其网格需求量大、计算耗时，在工程应用方面依然未能得到广泛的推广应用。为解决工程实际的需要，针对这一情况，Spalart^[2]提出了一种脱体涡模拟 (Detach ed eddy simulation, DES) 技术，并将其应用于SA湍流模型，其核心思想是在近壁面附着涡主导区域仍采用RANS的湍流模型，在离开物面附近的脱体涡主导区域采用大涡模拟 (Large eddy simulation, LES)。之后, Strelets^[3]结合两方程SST湍流模型发展了SST-DES模式，围绕典型大分离流动基础算例，开展了SST-DES、SA-DES以及URANS结果的详细比对。Ment er等人^[4]改进了Strelets发展的SST-DES模式，给出了一种SST-DDES实现模式。Sp al art^[5]对DES技术优缺点与发展方向做了综述性报告。国内陈龙^[6]，刘学强^[7]等也对DES技术展开了应用研究。

网格自适应技术^[8](Adaptive mesh refinement，AMR) 通过设置自适应判据，能够动态地跟踪、捕捉、保持流场特征。因此，AMR结合DES技术，基于流场的涡特征，能够在旋涡集中区域布置合适的网格尺度，从而成为处理大分离流动问题的有效方法。文献[9]中介绍了一种自适应混合笛卡尔网格 (Adaptive hybrid Cartesian grid，AHCG) 方法，具有网格生成迅速、自动化程度高等优点。同时，文献[9]介绍了基于流场特征 (如旋涡、激波等) 的网格自适应技术，以提高流场模拟精度。本文在此基础上，提出一种改进的基于流场涡特征的网格自适应技术，针对高雷诺数湍流大分离流动问题，利用格心型有限体积方法，采用SST k-ω湍流模型的DDES方法模拟湍流，并引入隐式双时间步的低速预处理技术，数值求解了Navier-Stokes方程，并以此为技术手段，模拟了低速大迎角三角翼绕流及大迎角椭球绕流问题，计算结果展示了由于分离引起的旋涡生成、发展及破碎过程，所得的定量计算数据与实验结果相吻合，表明SST-DDES模型耦合低速预处理技术可以较好地模拟高雷诺数大分离流动问题。

1 数学模型与数值离散 1.1 流动控制方程与数值解法

采用Navier-Stokes方程作为黏性大分离流动的控制方程，其积分守恒型形式如下

$ \frac{\partial }{{\partial t}}\int_\mathit{\Omega } {\boldsymbol{W}d\mathit{\Omega }} + \oint\limits_{\partial \mathit{\Omega }} {({\boldsymbol{F}_c}{\boldsymbol{F}_v})} {\rm{d}}S = 0 $

(1)

式中：W为守恒变量; F_c和F_v分别为对流通量和黏性通量，具体形式见文献[10]。

采用格心型有限体积法求解方程式 (1)，求解过程如下，具体细节见文献[9]：

(1) 空间项离散：对流项使用HLLC^[11]格式求解，并通过解的线性重构^[10]获得二阶精度，黏性项采用二阶中心型格式离散。

(2) 时间项离散：采用双时间步推进方法^[12]，在每一个物理时间步长内，采用LU-SG S隐式迭代方法^[13]，并利用当地时间步长加速收敛。

(3) 湍流模型：采用SST k-ω湍流模型^[14]。

(4) 边界与初始条件：壁面处采用无滑移边界条件，远场边界处采用考虑低速预处理计算的远场特征边界条件^[10]。初始流场为均匀流场，SST k-ω湍流模型来流值的设置方法见文献[10]。

(5) 网格生成：采用文献[9]描述的AHCG方法生成计算网格，并根据流场特征进行网格自适应。

1.2 双时间步LU-SGS低速预处理技术

当采用可压缩Navier-Stokes方程求解低速流动问题时，由于方程刚性的影响，会导致计算收敛速度慢且难以获得准确的流场结果，因此，需在双时间步推进方法^[12]的基础上，引入预处理矩阵Γ

$ \mathit{\boldsymbol{ \boldsymbol{\varGamma} }}\frac{{\partial ({\mathit{\Omega }_I}{\boldsymbol{Q}_I})}}{{\partial \tau }} + \frac{{\partial ({\mathit{\Omega }_I}{\boldsymbol{W}_I})}}{{\partial t}} =-{\boldsymbol{R}_I} $

(2)

式中：R，Ω表示残差项和控制体体积，下标I表示控制体序号。将虚拟时间步中的原始变量Q改成守恒变量W后，式 (2) 变为

$ \mathit{\boldsymbol{ \boldsymbol{\varGamma} M}}{^{-1}}\frac{{\partial ({\mathit{\Omega }_I}{\mathit{\boldsymbol{W}}_I})}}{{\partial \tau }} + \frac{{\partial ({\Omega _I}{\mathit{\boldsymbol{W}}_I})}}{{\partial t}} =-{\mathit{\boldsymbol{R}}_I} $

(3)

式中：雅克比矩阵M=əQ/əW。双时间步方法中物理时间采用二阶向后欧拉格式离散，伪时间步采用一阶向后欧拉格式离散，离散后的方程式 (3) 变为

$ \left[{(\frac{\mathit{\boldsymbol{ \boldsymbol{\varGamma} }}}{{\Delta {\tau _I}}} + \frac{{3\mathit{\boldsymbol{M}}}}{{2\Delta t}}{\rm{ }}){\rm{ }}{\mathit{\boldsymbol{M}}^{-1}}\mathit{\Omega } + \frac{{\partial {\mathit{\boldsymbol{R}}_I}}}{{\partial {\mathit{\boldsymbol{W}}_I}}}} \right]\Delta W_I^* = -{(\mathit{\boldsymbol{R}}_I^*)^l} $

(4)

令

$ {\mathit{\boldsymbol{ \boldsymbol{\varGamma} }}_u} = \mathit{\boldsymbol{ \boldsymbol{\varGamma} }} + \frac{{3\Delta {\tau _I}}}{{2\Delta {t_I}}}\mathit{\boldsymbol{M}}, \frac{1}{b} = 1 + \frac{{3\Delta {\tau _I}}}{{2\Delta {t_I}}} $

(5)

式中：Δτ和Δt分别表示虚拟时间步长和真实时间步长。因此，非定常预处理矩阵Γ_u可写成

$ {\mathit{\boldsymbol{ \boldsymbol{\varGamma} }}_u} = \left[{\begin{array}{*{20}{c}} {\mathit{\Theta '}} & 0 & 0 & 0 & {{\rho _T}}\\ {\mathit{\Theta '}u} & \rho & 0 & 0 & {{\rho _T}u}\\ {\mathit{\Theta '}v} & 0 & \rho & 0 & {{\rho _T}v}\\ {\mathit{\Theta '}w} & 0 & 0 & \rho & {{\rho _T}w}\\ {\mathit{\Theta '}H-1} & {\rho u} & {\rho v} & {\rho w} & {{\rho _T}H + \rho {c_P}} \end{array}} \right] $

(6)

式中：$\mathit{\Theta }\prime = b\mathit{\Theta }-\left( {b-1} \right){\rho _P}, {\rho _P} = \frac{{\partial \rho }}{{\partial P}}{\left| {_T;{\rho _T} = \frac{{\partial \rho }}{{\partial T}}} \right|_P}$; c_P表示定压比热容。

1.3 SST-DDES模型

SST k-ω两方程湍流模型^[14]的积分形式是

$ \frac{\partial }{{\partial t}}\int_\mathit{\Omega } {{\mathit{\boldsymbol{W}}_T}{\rm{d}}\mathit{\Omega }} + \oint\limits_{\partial \mathit{\Omega }} {} ({\mathit{\boldsymbol{F}}_{c, {\rm{ }}T}}-{\mathit{\boldsymbol{F}}_{v, T}}){\rm{d}}S = \int_\mathit{\Omega } {{\mathit{\boldsymbol{Q}}_T}{\rm{d}}\mathit{\Omega }} $

(7)

式中：W_T为守恒变量；F_{c, T}和F_{v, T}分别为对流项和黏性项；Q_T为湍流输运方程的源项，具体形式为

$ \begin{array}{l} {\mathit{\boldsymbol{W}}_T} = \left[\begin{array}{l} \rho K\\ \rho \omega \end{array} \right], {\mathit{\boldsymbol{F}}_{c, {\rm{ }}T}} = \left[\begin{array}{l} \rho KV\\ \rho \omega V \end{array} \right]\\ {\mathit{\boldsymbol{F}}_{v, T}} = \left[\begin{array}{l} {n_x}\tau _{xx}^K + {n_y}\tau _{yy}^K + {n_z}\tau _{zz}^K\\ {n_x}\tau _{xx}^\omega + {n_y}\tau _{yy}^\omega + {n_z}\tau _{zz}^\omega \end{array} \right]\\ {\mathit{\boldsymbol{Q}}_T} = \left[{\begin{array}{*{20}{c}} {P-{\beta ^*}\rho \omega K}\\ {\frac{{{C_\omega }\rho }}{{\mu T}}P-\beta \rho {\omega ^2} + 2(1-{f_1})\frac{{\rho {\sigma _{\omega 2}}}}{\omega }\frac{{\partial K}}{{\partial {x_j}}}\frac{{\partial \omega }}{{\partial {x_j}}}} \end{array}} \right] \end{array} $

以上各式中的具体参数见文献[14]。

针对高雷诺数大分离流动问题，基于SST k-ω两方程湍流模型，定义DES尺度${\tilde l}$

$ \tilde l = \min ({l_{k-\omega }}, {C_{{\rm{DES}}}}\mathit{\Delta }) $

(8)

其中

$ {l_{k-\omega }} = \frac{{{k^{\frac{1}{2}}}}}{{{\beta ^*}\omega }}\;\;\mathit{\Delta = }{\rm{max(}}{\mathit{\Delta }_x}{\rm{, }}{\mathit{\Delta }_y}{\rm{, }}{\mathit{\Delta }_z}{\rm{)}} $

(9)

网格尺度Δ为3个空间方向上的最大网格步长，C_DES为自适应参数

$ {C_{{\rm{DES}}}} = (1-{f_1})C_{{\rm{DES}}}^{k-\varepsilon } + {f_1}C_{{\rm{DES}}}^{k-\omega } $

(10)

式中：f₁由湍流模型给出，C_DES^k-ε=0.61，C_DES^k-ω=0.78。

SST-DDES^[4]模型的实现只需替换k方程中的耗散项，即由D_RNAS^k=β^*ρkω替换为D_DES^k

$ D_{{\rm{DES}}}^k = {\beta ^*}\rho k\omega {F_{DES}} $

(11)

式中：${F_{{\rm{DES}}}} = \max (\frac{{{l_{k-\omega }}}}{{{C_{{\rm{DES}}}}\mathit{\Delta }}}(1-{f_{{\rm{SST}}}}), 1)$，引入了开关函数f_SST=0, f₁, f₂。本文在进行模拟时取f_SST=f₂。

2 算例分析 2.1 大迎角三角翼低速绕流分离问题

大迎角ONERA 70°后掠三角翼低速扰流问题是典型的大分离流动问题，它具有涡的脱落、二次涡的生成、分离再附等复杂流动特征，是湍流大分离流动中很具有挑战性的经典算例，2000年，Mitchell等人^[15]在Onera-F2风洞中对该外形进行了相关实验。计算模型及尺寸如图 1所示，计算条件：迎角α=27°，来流马赫数Ma_∞ =0.069，雷诺数Re=1.56×10⁶。计算网格采用AHCG方法^[9]生成，在定常求解计算过程中，基于流场的涡特征进行4次网格自适应加密，最终形成的计算网格总数目为3 272 356，图 2给出经过自适应后的网格图。无量纲物理时间步长为5×10^-4，每一物理时间步内虚拟迭代步数为1 000步，且每50物理时间步进行一次动态网格自适应。图 3给出了t=2.5时刻，使用SST-DDES模型计算得到三角翼不同剖面的涡量瞬态结果，从该图中可以清晰地看到三角翼背风面上的主涡与二次涡等涡系结构，此时在靠近背风面的后半部分流动会趋于不稳定，涡系处于破碎状态。

依据Morton^[16]提出的涡核破碎点位置的判断方式，图 4展示了本文计算得到涡核破碎点的位置，在涡核破碎点处，流线显示出了逆向运动，约位于x=0.668处。为细致地比较背风面上的流动情况，图 5定量比较了背风面上4个典型站位 (x=500, 600, 700, 8 00 mm) 上的压力系数分布，本文结果在压力峰值上明显优于Morton采用SA-DES的结果^[16]，但是皆与实验值 (图 5中以“EXP”表示) 存在一定的误差，可能与实验中存在洞壁干扰、选定的参考动压不同以及网格敏感性等因素相关。

2.2 大迎角6:1椭球黏性绕流问题

大迎角椭球体黏性绕流问题，虽然几何构型简单，但却包含了十分复杂的流动现象：背风面首先会形成两个对称卷起的主涡结构，并且会有流动再附现象产生，且随着主涡与附面层作用，会伴随着二次涡结构的形成。本文以该算例考核数值方法的有效性，计算条件：迎角α=20°，来流马赫数Ma_∞=0.135，基于长轴直径c的雷诺数Re=4.2×10⁶。

采用自适应混合笛卡尔网格生成计算网格，其中贴体部分的结构网格数目为200× 37×120(流向、法向、周向)，第一层网格的厚度为1.0×10^-5c，外围初始笛卡尔网格数目为316 016。数值求解过程中，经过3次基于旋度的解自适应后，在背风面涡主导区域的笛卡尔网格都得到了加密，经过加密后的笛卡尔网格数目增加到376 594，整体混合网格总数为1 264 594。

图 6通过物面的极流线分布展示了主涡分离线的位置S₁，二次涡分离线的位置S₂与二次涡的再附线的位置R₂，主涡再附线的位置R₁位于z=0对称面附近。表 1将流向x/L=0.77站位处S₁, S₂与R₂位置的本文计算结果与文献值^[17]以及实验值^[18]进行了对比，可以看到S₁与S₂的值与实验吻合很好，R₂的位置与文献值^[17]相近，验证了本节混合笛卡尔网格方法耦合DES技术对于大分离流动问题的模拟能力。

为定量分析背风面流动情况，图 7中给出了x/L=0.77横截面上沿周向的压力系数分布。图 8中给出了该截面上4个不同的典型方位角下3个方向上的无量纲速度随着法向的速度分布。可以看出，在该横截面上，压力分布与实验值^[18]吻合较好，只在二次峰值的大小与位置上存在一点差异。在90°方位角上 (y=0对称面) 的速度分布与实验值吻合良好。u/U_∞ (U_∞为来流速度) 呈现出对数分布趋势，在y_n/L≈3.0×10^-3达到最大值；v/U_∞的值很小，几乎靠近为0，较之于实验值稍稍偏大；w/U_∞的值也在y_n/L≈3.0×10^-3到达最大值，与测量结果完全一致。在120°方位角 (背风面转30°) 上u/U_∞与w/U_∞随法向的变化形状较之于前面一个站位 (90°) 更为陡峭，该站位位于背风面流动分离区域，3个方向的速度型计算结果与实验值吻合很好，证明了本文数值方法能够在流动分离区域正确地捕捉边界层分离。在150°方位角上，该位置处于二次涡结构附近，从图 8(c)中可以看出，u/U_∞的结果与实验值吻合良好，但是v/U_∞与w/U_∞的结果却存在明显的差异。在180°方位角 (z=0对称面) 上，u/U_∞和w/U_∞与实验值吻合良好，v/U_∞较实验值稍稍偏大。

3 结束语

本文为解决高雷诺数湍流大分离低速流动问题，采用了非定常双时间步LU-SGS方法与低速预处理技术，并基于SST k-ω湍流模型研究了延迟脱体涡模拟技术，建立了基于网格自适应技术的数值求解Navier-Stokes方程的方法。文中的两个三维数值算例表明，延迟脱体涡模拟技术耦合低速预处理技术，能够有效模拟低速大分离流动问题，揭示流动机理，并且定量数值模拟结果与实验数据相近，该方法为模拟大分离流动及其他旋涡主导流动问题提供了有益的借鉴和参考。