引用本文
孙逸凡, 朱炜, 吴玉新, 祁海鹰. Gao-Yong湍流模型对边界层转捩的适用性研究[J]. 清华大学学报(自然科学版), 2023, 63(4): 642-648.  
SUN Yifan, ZHU Wei, WU Yuxin, QI Haiying. The applicability study of Gao-Yong turbulence model to boundary layer transitions[J]. Journal of Tsinghua University (Science and Technology), 2023, 63(4): 642-648.  
Gao-Yong湍流模型对边界层转捩的适用性研究
孙逸凡, 朱炜, 吴玉新, 祁海鹰    
清华大学 能源与动力工程系，北京 100084
收稿日期：2022-11-04
基金项目：国家科技重大专项(Y2019-I-0022-0021)
作者简介：孙逸凡(1999—)，男，博士研究生
通讯作者：祁海鹰，教授，E-mail: hyqi@tsinghua.edu.cn
摘要：准确预测边界层转捩对深入理解湍流边界层发展和飞行器翼型优化工业应用等具有重要意义。为提升Gao-Yong(G-Y)湍流模型对平板边界层发展及转捩特性预测的准确性与适用性，该文提出了采用G-Y漂移位矢Reynolds数(ReT)表征当地湍流脉动强度的方法，建立实度系数(Cs)与ReT间的函数关系，解决G-Y湍流模型中Cs取值间断、分布离散的问题。基于开源OpenFOAM软件开发了改进G-Y湍流模型的求解器，对经典平板边界层和转捩实验开展计算流体力学(computational fluid dynamics，CFD)数值模拟并验证。研究结果表明，相比k-ω等传统雷诺平均(Reynolds averaged Naviers-Stokes，RANS)模型，改进后的G-Y湍流模型不仅能准确模拟边界层流动特征，还能合理预测边界层发展过程及转捩位置，准确性较高。结果表明，经过改进的G-Y湍流模型可进一步用于对边界层转捩特性的研究。
关键词：数值模拟    G-Y湍流模型    平板边界层    转捩    
The applicability study of Gao-Yong turbulence model to boundary layer transitions
SUN Yifan, ZHU Wei, WU Yuxin, QI Haiying    
Department of Energy and Power Engineering, Tsinghua University, Beijing 100084, China
Abstract: Objective Transition is one of the most significant progressions in fluid mechanics. Accurate boundary layer transition prediction is also essential for a complete understanding of the production of turbulent boundary layers and the optimization of airfoil shapes for industrial applications. Turning turbulence models can accurately forecast transition by including an empirical turning criterion, a low Reynolds number viscous flow field characteristic, or an intermittency component for flow in the transition zone. As for the conventional turbulence models using Reynolds averaging method, such as the k-ε and the k-ω models, accurately forecasting this phenomenon is tough since the transition process is poorly characterized or described. Methods A functional correction relation for the G-Y model is suggested and put to the test in the calculations to increase its precision and usefulness for forecasting boundary layer development and transitional features across flat plates. The coefficient of substance Cs had the problem of discontinuous distribution and intermittent values in the original model. This work provided a correlation formula between Cs and the drift vector Reynolds number ReT to solve this problem. The drift vector Reynolds number ReT is built using the G-Y model's drift velocity and vector to describe the strength of local turbulent pulsations. With the aid of physical understanding, the Cs-ReT relationship is produced, resulting in a continuous distribution of the real degree coefficient with ReT variation, close to 0 at the location of mild turbulent pulsation and close to 2/3 at the location of strong turbulent pulsation. The turbulent viscosity is bound by the size of the real degree coefficient at the point of weak turbulent pulsation by solving the discrete distribution's initial problem in this manner. Results A solver for the modified G-Y turbulence model was made using the open-source program OpenFOAM. The CFD numerical simulations and validation were then performed for the tests involving the T3A and T3B transitions as well as the conventional flat plate boundary layer. The G-Y model's computational findings showed that: (1) The formula of Cs and ReT considerably increased the G-Y model's accuracy for boundary layer instances and gave the model the ability to forecast transitions. The G-Y model accurately predicted the transitional positions of the two boundary layer experiments, T3A and T3B. (2) Results of the modified G-Y model were in good agreement with the experimental and theoretical values for the distribution of surface friction coefficients before and after the turn. The relative error of the friction coefficients in the segment just entering turbulence was only 3%. (3) The G-Y model accurately replicated the transition from laminar to turbulent flow, which caused the velocity profile of the boundary layer to change from fullness to loss. Conclusions The findings show that the improved G-Y model has clear advantages over the established k-ω model in terms of its ability to accurately simulate the boundary layer transition process of a flat plate. This model can be used to investigate the properties of boundary layer transitions in greater depth.
Key words: numerical simulation    G-Y turbulence model    flat plate boundary layer    turbulent transition    

湍流是流体复杂的流动形式，是经典物理学的难题，对传热、摩擦、边界层发展等问题有重要影响。20世纪以来，随着计算机技术及计算流体动力学(computational fluid dynamics，CFD)的快速发展，人们有了探究湍流这一物理问题的机会。实际上，传统湍流模式的理论并未很好地解决这一问题，很多流动问题也难以得到令人满意的答案，对湍流的研究仍然具有很高的理论意义和应用价值。因此本文开展了对平板流动边界层转捩与发展预测问题的研究。

转捩问题由于其强烈的非定常、非线性特性，一直是流体力学领域的核心问题之一。传统的雷诺平均(Reynolds averaged Naviers-Stokes，RANS)方法缺少对转捩过程的认识和描述，经典湍流模型如一方程模型、二方程模型和Reynolds应力模型等，皆难以对这一现象进行准确预测，为了在模拟中准确反映转捩区流动特征，学者们建立了相应的转捩湍流模型，类似经典湍流模型的分类方法，可主要分为零方程、单方程和多方程3类模型。零方程转捩模型从大量实验数据中总结出各种情况下判断转捩的经验关系，并直接应用或融入现存代数模型后，间接应用于流场分析中。其中，常被使用的转捩经验准则包括Abu-ghannam和Shaw[1]的经验准则和Walker和Gostelow[2]的最小长度理论等。直接使用单、多方程经典湍流模型预测转捩的方法仍有待改进[3-4]，所以人们将低Reynolds数黏性流场特征或转捩区流动的间歇因子引入湍流模型，以反映转捩边界层非平衡黏性流动现象；目前涉及的模型有k方程、k-ε方程、RST(Reynolds stress transport)模型和k-ε-γ方法[5]等多种形式。

高歌和熊焰[6]建立的Gao-Yong(G-Y)湍流模型不同于传统RANS方法。该模型通过侧偏平均处理将湍流脉动分解为平均流和漂移流；在保留了湍流脉动的一阶统计信息的基础上，通过严谨的推导获得了对称的漂移流连续方程和动量方程。该模型在充分考虑湍流实际物理过程的基础上，通过湍流的多尺度现象引入级数形式的机械能方程[7]。这样漂移流的连续性方程、动量方程和机械能方程3个方程共同使G-Y湍流模型封闭，且不含任何的经验系数和壁面函数。该模型自1999年提出，经过了诸多算例的检验，在可压、不可压、二维流动和三维流动中均展现出较高的精度，具有较强的适应性与生命力，但对边界层流动现象的预测能力仍有待探究。

为此，本文对G-Y在平板边界层转捩与发展中的模拟精度进行研究。首先对模型中的实度系数提出了修正关系式，通过与当地湍流强度建立关系，将原本离散的实度系数值变为连续分布，使模型具备预测转捩的能力，模拟精度有所提高。在此基础上，本文通过标准平板边界层流动算例和T3A、T3B平板边界层转捩实验工况对改进的G-Y湍流模型开展验证，同时比较了k-ω湍流模型和k-ω SSTLM(k-ω shear stress transport Langtry-Menter)低Re湍流模型对此类问题的预测性能，并对改进G-Y湍流模型的转捩预测精度进行分析。

1 数值模型与方法 1.1 控制方程

G-Y湍流模型的湍流参量为漂移速度(${\mathit{\boldsymbol{\tilde u}}_i}$)、漂移压力($\tilde p$)和漂移位矢(lj)。其中，${\mathit{\boldsymbol{\tilde u}}_i}$和$\tilde p$为通过侧偏平均得到的速度、压力脉动量。N-S(Navier-Stokes)方程的非线性使其在高Re时会有无穷多解；与传统RANS方法的系综平均不同，侧偏平均方法将所有的脉动值按照一定的划分标准分成2组，第一组有NⅠ个解，第二组有NⅡ个解，2组脉动速度的系综平均表示如下：

$ {\mathit{\boldsymbol{\bar u}}_{'{\rm{Ⅰ}}}} = \mathop {\lim }\limits_{{N_{\rm{Ⅰ}}} \to \infty } \frac{1}{{{N_{\rm{Ⅰ}}}}}\sum\limits_{k = 1}^{{N_{\rm{Ⅰ}}}} {{\mathit{\boldsymbol{u}}_{{\rm{'Ⅰ}}k}}} , $ (1)
$ {\mathit{\boldsymbol{\bar u}}_{'{\rm{Ⅱ}}}} = \mathop {\lim }\limits_{{N_{\rm{Ⅱ}}} \to \infty } \frac{1}{{{N_{\rm{Ⅱ}}}}}\sum\limits_{k = 1}^{{N_{\rm{Ⅱ}}}} {{\mathit{\boldsymbol{u}}_{{\rm{'Ⅱ}}k}}} . $ (2)

其中：${\mathit{\boldsymbol{\bar u}}_{'{\rm{Ⅰ}}}}$和${\mathit{\boldsymbol{\bar u}}_{'{\rm{Ⅱ}}}}$为2组脉动速度的系综平均，“－”为系综平均；“′”为脉动量。

由于所有脉动速度的总系综平均为零，故有

$ \begin{array}{c} {{\mathit{\boldsymbol{\bar u}}}_{'}} = \mathop {\lim }\limits_{N \to \infty } \frac{1}{N}\sum\limits_{k = 1}^N {{\mathit{\boldsymbol{u}}_{'k}}} = \\ \mathop {\lim }\limits_{{N_{\rm{Ⅰ}}} \to \infty } \frac{{{N_{\rm{Ⅰ}}}}}{N}\frac{1}{{{N_{\rm{Ⅰ}}}}}\sum\limits_{k = 1}^{{N_{\rm{Ⅰ}}}} {{\mathit{\boldsymbol{u}}_{'Ⅰk}}} + \mathop {\lim }\limits_{{N_{{\rm{Ⅱ}}}} \to \infty } \frac{{{N_{{\rm{Ⅱ}}}}}}{N}\frac{1}{{{N_{{\rm{Ⅱ}}}}}}\sum\limits_{k = 1}^{{N_{{\rm{Ⅱ}}}}} {{\mathit{\boldsymbol{u}}_{'Ⅱk}}} = 0. \end{array} $ (3)

其中，$\mathop {\lim }\limits_{{N_{\rm{Ⅰ}}} \to \infty } \frac{{{N_Ⅰ}}}{N}{\mathit{\boldsymbol{u}}_{'{\rm{Ⅰ}}}}$与$\mathop {\lim }\limits_{{N_{\rm{Ⅱ}}} \to \infty } \frac{{{N_Ⅱ}}}{N}{\mathit{\boldsymbol{u}}_{'{\rm{Ⅱ}}}}$分别为2组的加权脉动速度，也即漂移速度。漂移压力定义与速度相同。

由式(3)可知，漂移量是对称的，任取其中一组漂移速度、压力进行推导，即可得到G-Y湍流模型的漂移速度连续性方程与动量方程。这2个方程与求解漂移位矢lj的机械能方程一同使得模型封闭，详细推导见文[8]。此处直接给出稳态不可压条件下G-Y湍流模型的平均流连续方程、平均流动量方程、漂移流连续方程、漂移流动量方程、漂移流机械能方程，依次表示如下：

$ \frac{{\partial {{\mathit{\boldsymbol{\bar u}}}_i}}}{{\partial {x_i}}} = 0, $ (4)
$ \frac{{\partial {{\mathit{\boldsymbol{\bar u}}}_i}}}{{\partial t}} + {\mathit{\boldsymbol{\bar u}}_j}\frac{{\partial {{\mathit{\boldsymbol{\bar u}}}_i}}}{{\partial {x_j}}} = - \frac{1}{\rho }\frac{{\partial \bar p}}{{\partial {x_i}}} + v\frac{{{\partial ^2}{{\mathit{\boldsymbol{\bar u}}}_i}}}{{\partial {x_j}\partial {x_j}}} - \frac{{\partial {{\mathit{\boldsymbol{\bar \sigma }}}_{ij, {\rm{T}}}}}}{{\partial {x_j}}}, $ (5)
$ \frac{{\partial {{\mathit{\boldsymbol{\tilde u}}}_i}}}{{\partial {x_i}}} = 0, $ (6)
$ \begin{array}{c} \frac{{\partial {{\mathit{\boldsymbol{\tilde u}}}_i}}}{{\partial t}} + {{\mathit{\boldsymbol{\tilde u}}}_j}\frac{{\partial {{\mathit{\boldsymbol{\tilde u}}}_i}}}{{\partial {x_j}}} + {{\mathit{\boldsymbol{\tilde u}}}_j}\frac{{\partial {{\mathit{\boldsymbol{\bar u}}}_i}}}{{\partial {x_j}}} = \\ - \frac{1}{\rho }\frac{{\partial \tilde p}}{{\partial {x_i}}} + v\frac{{{\partial ^2}{{\mathit{\boldsymbol{\tilde u}}}_i}}}{{\partial {x_j}\partial {x_j}}} + \frac{{\partial {{\mathit{\boldsymbol{\tilde \sigma }}}_{ij, {\rm{T}}}}}}{{\partial {x_j}}} - v\frac{{{\partial ^2}{{\mathit{\boldsymbol{\bar u}}}_i}}}{{\partial {x_j}\partial {x_j}}} - \frac{{\partial {{\mathit{\boldsymbol{\bar \sigma }}}_{ij, T}}}}{{\partial {x_j}}}, \end{array} $ (7)
$ \begin{array}{c} \mathop \sum\limits_{n = 1}^\infty \frac{{\left( {{{\left( {{l_j}\frac{\partial }{{\partial {x_j}}}} \right)}^n}{{\mathit{\boldsymbol{\bar u}}}_i}} \right) \cdot {{\mathit{\boldsymbol{\tilde u}}}_i}}}{{n!}} = \\ {\mathit{\boldsymbol{l}}_i} \cdot \left( { - \frac{1}{\rho }\frac{\partial }{{\partial {x_j}}}\left( {{{\mathit{\boldsymbol{\bar u}}}_i}{{\mathit{\boldsymbol{\tilde u}}}_j} + \tilde p{\mathit{\boldsymbol{\delta }}_{ij}}} \right) + v\frac{{{\partial ^2}{{\mathit{\boldsymbol{\tilde u}}}_i}}}{{\partial {x_j}\partial {x_j}}} + } \right.\\ \left. {\frac{{\partial {{\mathit{\boldsymbol{\tilde \sigma }}}_{ij, {\rm{T}}}}}}{{\partial {x_j}}} - v\frac{{{\partial ^2}{{\mathit{\boldsymbol{\bar u}}}_i}}}{{\partial {x_j}\partial {x_j}}} - \frac{{\partial {{\mathit{\boldsymbol{\bar \sigma }}}_{ij, {\rm{T}}}}}}{{\partial {x_j}}}} \right). \end{array} $ (8)

其中：${\mathit{\boldsymbol{\bar \sigma }}_{ij, {\rm{T}}}}$与${\mathit{\boldsymbol{\tilde \sigma }}_{ij, {\rm{T}}}}$分别为平均流湍流应力和漂移流湍流应力，i、j为Einstein张量表示法，T为湍流，“~”为G-Y模型的侧偏平均；δij为Kronecker符号；t为时间；ρ为密度；υ为层流黏性；n为阶数；x为坐标方向。

G-Y湍流模型在计算湍流应力时，首先经式(9)计算各向异性涡黏系数张量ξij。ξij在平均流方向上建立的正交各向异性坐标系[7]，将其经过坐标变换，可得到(x, y, z)坐标系下的湍流黏性系数。

$ {\mathit{\boldsymbol{\xi }}_{ij}} = {\rm{Cs}} \cdot \rho {\mathit{\boldsymbol{\tilde u}}_i}{\mathit{\boldsymbol{l}}_j}. $ (9)

其中：ξij为涡黏系数张量；Cs为G-Y湍流模型提出的实度系数，表征湍流脉动受限制的大小。

Cs并非经验系数，而是由lj的活动域受限比例、${\mathit{\boldsymbol{\tilde u}}_i}$的活动域受限比例和微元体承受湍流应力的面积比这3个物理因素决定。在不同情况下，Cs的取值为：三维流动外流，1；二维流动外流，2/3；三维边界，3/16；二维边界，1/8[9]。此外，文[10-11]也采用通过速度梯度构造判别式等方法，区分Cs不同值的作用域。Cs采用不同取值的原因是考虑各种物理因素对湍流黏性的影响，区别不同区域中湍流应力受影响的程度。但这种取值方式也导致了Cs的分布存在间断的情况。例如，在靠近壁面的位置上，Cs大小为1/8，但在距离壁面一定距离后，Cs的值会一跃变为2/3，如图 1所示。红线为由于Cs只有2种值而带来的间断。

图 1 Cs的间断分布示意图
图选项

为解决Cs取值不连续的问题，本文提出了G-Y湍流模型的ReT，并建立Cs与ReT的修正关系式，表示如下：

$ R{e_{\rm{T}}} = \frac{{\left| {{{\mathit{\boldsymbol{\tilde u}}}_i}{\mathit{\boldsymbol{l}}_i}} \right|}}{v}, $ (10)
$ {\rm{Cs}} = \frac{2}{3}\left( {1 - {{\rm{e}}^{ - R{e_{\rm{T}}}/390}}} \right){\rm{.}} $ (11)

其中ReT为湍流模型的漂移位矢Re。无量纲量ReT用以表征当地湍流脉动的强弱。

式(11)为Cs与ReT修正关系式，函数曲线如图 2所示，Cs与ReT均为无量纲。从图中可以看出，该函数使Cs变为随ReT变化的连续分布，在湍流脉动较弱的地方接近0，在湍流脉动较强的地方则接近红色虚线的值0.667(约为2/3)，而2/3也是原G-Y湍流模型中Cs的最大值。Cs与ReT的修正关系式借助物理认识建立，在湍流脉动较弱的地方通过Cs对湍流脉动进行限制，解决了原来分布离散的问题，使Cs的值在流场中连续分布。

图 2 Cs-ReT函数关系示意图
图选项

1.2 数值方法

传统湍流模型在计算时直接求解湍流参量的微分方程，而改进G-Y湍流模型则使用SIMPLE算法求解湍流场中耦合的漂移速度与漂移压力。本文充分利用OpenFOAM软件开源的特性，基于软件原有的不可压simpleFoam求解器，开发了一套改进G-Y湍流模型求解器。图 3为求解器计算流程图。

图 3 求解器计算流程
图选项

为分析改进后G-Y湍流模型对平板边界层转捩预测特性，本文分别采用Re=6×106的标准平板边界层流动算例和T3A、T3B平板边界层转捩实验算例进行模型验证，并比较了k-ω模型和k-ω SSTLM低Re数模型对T3A、T3B实验的预测性能。

在上述的模拟中均采用二维矩形网格，其几何结构、边界分布和网格加密方式均类似，如图 4所示。Rumsey和Gatski[12]在2003年通过网格比较研究，得到一种适用于在Re=6×106的标准平板边界层流动中进行湍流模型评估的默认网格，流向方向网格均布，宽度为0.028 033 m；在法向方向网格向壁面加密，相邻网格拉伸比为1.18，第一层网格高度为1×10-6 m，无量纲高度Y+约为0.2；当网格高度达到0.020 833 m后，网格在法向等间距分布。3种不同湍流模型均在网格下开展数值模拟。

图 4 平板边界层网格及计算域示意图
图选项

2 结果与分析 2.1 Re=6×106的标准平板边界层流动算例

为检验改进后G-Y湍流模型计算准确性，本文首先对平板边界层的经典算例进行了计算。算例各项设置与文[12]相同：平板长L=1.000 m，计算域高H=1.000 m，Re=6×106，网格数为65×97。计算区域进口取在平板前缘1/3 m处。右、上边界均为出口，平板为无滑移壁面，平板前部过渡段为对称边界。

图 5为本文分别使用改进前、后的G-Y湍流求解器得到的表面摩擦系数Cf分布。由图 5可知，原模型直接应用于平板边界层时，计算得到的Cf偏大，在入口处Cf出现了比较大的值。而在模型中加入了Cs-ReT修正关系后，G-Y湍流模型的计算结果显著改善。改进后的模型预测了来流无湍流下的自然转捩过程：从入口流过一段较短距离后，边界层由层流转变为湍流流动，Cf逐渐与湍流经验公式保持一致，误差较小。

图 5 改进前、后摩擦系数结果对比
图选项

Cs在Y=0.001 m和X=0.800 m的分布如图 6所示。u0/U∞表示当地速度与来流最大速度的比值。由图 6a可知，在转捩前，Cs一直保持较小的值，模型计算得到的湍流黏性因此较低，流动状态接近层流。在当地Rex约为2×105时，Cs逐渐增大，并在Rex为1×106时，逐渐接近红色虚线的值0.667(约为2/3)，标志着转捩的完成，边界层流动开始进入湍流阶段。由图 6b可知，Cs在壁面处为0，随着法向高度的增加，Cs迅速增大到最大值0.667(约为2/3)，随后在进入远场后逐渐降至0。通过与边界层内速度分布的对比发现，Cs较大且对湍流模型计算起显著作用的区域对应速度边界层的内部；而在湍流作用很小的位势流区域，Cs的值很小，湍流黏性较低，符合实际边界层的物理规律。

图 6 实度系数分布
图选项

2.2 T3A、T3B平板边界层转捩实验

T3A、T3B边界层转捩实验来自ERCOFTAC的T3系列平板边界层实验[13]，被广泛应用于转捩建模以及评价模型的转捩模拟精度。本文选取这2个零压梯度平板边界层转捩实验进行模拟，结合k-ω模型和k-ω SSTLM低Re模型，对比考察改进后的G-Y湍流模型对转捩过程的模拟精度。

在T3A实验中，进口来流速度为5.4 m/s，湍流强度为3%；T3B进口来流速度为9.4 m/s，湍流强度为6%。根据实验装置的实际大小，该算例的计算域大小设置为1.900 m×0.220 m，其中平板长L =1.700 m，计算域高H =0.200 m，网格总数为92×64。除上边界为滑移壁面外，其余边界条件设置均与上一算例相同。

各模型的Cf值计算结果如图 7所示。从图中可以看出，改进后的G-Y湍流模型展现出较好的转捩模拟能力。对2个算例，模型都较准确地预测出了转捩的发生位置，T3A实验转捩位置与实验测量值(黑点)几乎重合。与此同时，G-Y湍流模型模拟的层流和湍流段Cf值误差都很小，与实验、理论值(黑线)都十分吻合。值得注意的是，G-Y湍流模型对转捩过程的模拟与实际存在些许误差，具体表现为模拟转捩的发生较实际稍快，并且Cf在刚进入湍流段时略偏高(T3A相对误差为3.30%，T3B相对误差为3.25%)，但此偏差很快消除。这有可能是Cs- ReT修正关系式的函数形式不够准确，模型对当地湍流强度ReT较敏感导致的。因此，如何选择合适的函数也是后续工作中有待完善的地方。

图 7 改进G-Y湍流模型的Cf计算结果
图选项

k-ω模型[14]作为一个典型的二方程RANS模型，尽管计算得到的Cf误差较小，但无法模拟出转捩过程。这主要是因为传统RANS方法是针对充分发展的高Re湍流模型建立的，而转捩过程包含从层流到湍流变化范围很大的Re，因此k-ω这类传统RANS模型并不具备模拟低Re下转捩的能力。

作为k-ω SST模型在低Re下的修正，k-ω SSTLM模型[15]通过引入间歇因子和动量Re，具备了准确预测边界层转捩发生与发展的能力。在T3A实验的模拟中，k-ω SSTLM模型的Cf计算值在湍流段与实验值几乎一致，在转捩段仅与实验值相差约2×10-4。

图 8展示了经改进的G-Y湍流模型计算结果中，不同流向位置上的边界层速度变化规律。其中，Rex=0.162× 105位于层流区域，Rex=1.782×105位于转捩段，Rex=5.382×105位于湍流区域，对应的流向距离分别为0.045，0.495和1.495 m。从图中可以看出，G-Y湍流模型的模拟结果较好地复现了边界层由层流向湍流转捩过程中，边界层内的速度型曲线由饱满向亏损的变化过程，符合实际的物理规律。

图 8 改进G-Y湍流模型计算的边界层内速度分布
图选项

3 结论

本文采用ReT表征当地湍流脉动强弱的方法，建立了Cs与ReT的函数关系，进而改进了G-Y湍流模型，并通过开源软件OpenFOAM实现了数值求解。通过对平板边界层转捩的经典算例进行验证并分析，主要得到以下结论：

1) Cs与ReT的修正关系式使G-Y湍流模型显著提升了对标准平板边界层流动算例的预测精度，并使模型具备了预测转捩的能力。

2) 在转捩前后，改进G-Y湍流模型的表面摩擦系数计算值均与实验值、理论值吻合较好。

3) 与已有的k-ω模型和k-ω SSTLM低Re模型进行比较，改进后的G-Y湍流模型对平板边界层转捩过程具有准确的模拟能力，相比传统二方程湍流模型具有明显的优势。

参考文献
[1]
ABU-GHANNAM B J, SHAW R. Natural transition of boundary layers: The effects of turbulence, pressure gradient, and flow history[J]. Journal of Mechanical Engineering Science, 1980, 22(5): 213-228. DOI:10.1243/JMES_JOUR_1980_022_043_02
[2]
WALKER G J, GOSTELOW J P. Effects of adverse pressure gradients on the nature and length of boundary layer transition[J]. Journal of Turbomachinery, 1990, 112(2): 196-205. DOI:10.1115/1.2927633
[3]
JIANG H D, MOORE J, MOORE J. Low Reynolds number one-equation turbulence modelling for prediction of transitional flows over a flat plate[C]//Proceedings of the 28th Aerospace Sciences Meeting. Reno, USA: AIAA, 1990: 242.
[4]
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. DOI:10.1016/0017-9310(72)90076-2
[5]
VANCOILLIE G, DICK E. A turbulence model for the numerical simulation of the transition zone in a boundary layer[J]. International Journal of Engineering Fluid Mechanics, 1988, 1(1): 28-49.
[6]
高歌, 熊焰. 不可压湍流平均流与拟序流的各向异性模式理论及其数值验证[J]. 航空动力学报, 2000, 15(1): 1-11.
GAO G, XIONG Y. The united modeling theory for both of mean flow and coherent flow of incompressible turbulence[J]. Journal of Aerospace Power, 2000, 15(1): 1-11. (in Chinese)
[7]
GAO G, YONG Y. Partial-average-based equations of incompressible turbulent flow[J]. International Journal of Non-Linear Mechanics, 2004, 39(9): 1407-1419. DOI:10.1016/j.ijnonlinmec.2004.02.002
[8]
GAO G, XU J L. A partial average based study of compressible turbulent flows[J]. International Journal of Mechanic Systems Engineering, 2013, 3(1): 20-35.
[9]
高歌. Gao-Yong理性湍流方程[J]. 推进技术, 2010, 31(6): 666-675.
GAO G. A review on the research development of Gao-Yong equations of turbulent flows[J]. Journal of Propulsion Technology, 2010, 31(6): 666-675. (in Chinese)
[10]
陈广业, 高歌, 尹幸愉. 使用GAO-YONG方程组对不可压转捩/湍流平板边界层的计算[J]. 航空动力学报, 2003, 18(1): 197-202.
CHEN G Y, GAO G, YIN X Y. The calculation of incompressible wall boundary layer using GAO-YONG turbulence equations[J]. Journal of Aerospace Power, 2003, 18(1): 197-202. (in Chinese)
[11]
任鑫, 高歌. 使用GAO-YONG湍流方程组对翼型绕流的计算[J]. 航空学报, 2007, 28(Z1): S28-S34.
REN X, GAO G. The calculation of airfoil flows using GAO-YONG turbulence equations[J]. Acta Aeronautica et Astronautica Sinica, 2007, 28(Z1): S28-S34. DOI:10.3321/j.issn:1000-6893.2007.z1.006 (in Chinese)
[12]
RUMSEY C L, GATSKI T B. Summary of EASM turbulence models in CFL3D with validation test cases[R]. Virginia: Langley Research Center, 2003.
[13]
COUPLAND J. Classic collection database[EB/OL]. (1981-12-19)[2022-11-01]. http://cfd.mace.manchester.ac.uk/ercoftac/.
[14]
WILCOX D C. Formulation of the k-w turbulence model revisited[J]. AIAA Journal, 2008, 46(11): 2823-2838.
[15]
LANGTRY R B, MENTER F R. Correlation-based transition modeling for unstructured parallelized computational fluid dynamics codes[J]. AIAA Journal, 2009, 47(12): 2894-2906.