γ－Reθ转捩模型的标定与应用

牟 斌，江 雄，肖中云，陈作斌

（中国空气动力研究与发展中心，四川 绵阳 621000）

0 引 言

转捩边界层流动研究在CFD工程应用的许多方面非常重要，如层流翼型、风轮机、轮船船体、涡轮机叶栅等流动特性研究方面。这些流动的雷诺数范围确定了其层流区域和湍流区域往往具有相同量级，全湍流计算在阻力等方面带来的偏差不可忽视，必须考虑转捩的影响。

在过去的几十年中，全湍流的计算取得了辉煌的成就，并发展起以并行计算、非结构网格为代表的现代CFD技术；而转捩的数值模拟则严重滞后。Menter［1］将其归结为两方面的原因：一、转捩机制复杂，包含了自然转捩、bypass转捩、分离诱导转捩，同时湍流边界层在顺压梯度下可能再层流化，导致很难用一个简单的数学模型描述如此多的效应。二、转捩流动中线性和非线性效应相关，RANS平均消除了线性扰动的增长，因此传统的雷诺平均方程不能描述转捩流动。

在工程应用中，考虑转捩效应的经典方法是经验关联，就是把自由流的湍流度和当地的压力梯度等通过实验数据关联到转捩动量厚度雷诺数，如Abu－Ghannam和Shaw的关联等。遗憾的是以前的关联往往是非局部变量的函数，需要作诸如边界层厚度判定，驻点、分离点查找，沿流线积分等，这些算法受到外形、网格拓扑结构的制约，应用的范围很窄。

为了与现代 CFD 兼容，Menter与 Langtry［1－5］设计了γ－Reθ转捩模型。该模型采用经验关联的方法，并不模拟转捩的物理机理。其成功之处在于用应变率雷诺数触发转捩，并通过求解额外的两个输运方程避免了非局部变量的采用，从而形成与现代CFD相容的框架。在 Menter［2－5］早期的文章中，由于商业保密的原因，两个重要的经验关联函数没有给出，国外许多学者为了应用，依据平板实验结果提出了各种in－house公式［6－8］，并在各种 CFD 工程应用中获得成功；国内空气动力中心张玉伦［13－14］，西工大张晓东［15］等也提出了各自标定的经验关系式。虽然2009年Menter及Langtry［1］给出了一组关联函数，然而，转捩模型对不同的代码及标定过程中的微小差异非常敏感，因此这一组关联函数在CFX中能够得到很好的结果，在另一个解算器中的表现可能就会差一些［9］。为了较好模拟转捩效应，有必要在代码中实现转捩模型并标定关联函数。

本文首先在内部软件mbns3d上实现了转捩模型，依据零压力梯度平板转捩实验数据确定了一组关联函数，然后针对A翼型做了验证计算。随后数值模拟了NLR7301两段翼型，研究了转捩、远场边界条件及远场范围对精确计算阻力的影响，获得了有益结论。

1 γ－Reθ 转捩模型简介

（1）控制方程

无量纲的γ－Reθ输运方程的守恒形式为：

式中各项的详细解释可见文献［5］。

（2）关联函数

模型中存在三个关联函数：

对于关联函数Reθt，不同文献中的形式略有不同，本文采用Langtry［5］建议的形式。下面的工作是确定其余两个关联函数。

2 γ－Reθ 关联函数的确定

标定工作在内部软件mbns3d软件上进行，mbns3d采用多块结构网格（含拼接、对接、重叠），有限体积求解，具备多重网格、预处理、并行计算功能。

标定计算中，对流项离散选用三阶精度的Roe's FDS格式，粘性项二阶中心格式离散，AFADI方法求解离散方程。为了适应低速情况和加速收敛，采用了预处理技术和多重网格技术。

标定依据文献［5］提供的 T3A，T3B，T3A－和Schubauer＆Klebanof等平板试验数据，试验条件见表1。

表1 平板转捩试验的进口条件Table 1 Inlet condition for the flat plate test cases

平板长度选为2．0，来流远场边界距平板前缘0．22，法向远场边界距壁面0．22，网格为H型，壁面第一层网格法向距离1．0×10－5，平板前缘第一网格的x向距离1．5×10－3。网格规模为273×97，其中壁面上有241网格点，满足四重多重网格计算要求。

表1中的FSTI为是平板前缘的湍流强度，由于网格外边界距平板前缘0．22m，因此外边界的湍流强度实际上是依据下式估计［5］：

其中β＝0．09，β＊＝0．0828，x是下游距进口的流向距离，V是平均速度。

关联函数中Reθc和Flength是耦合性很强的函数，一个确定，另一个也随之确定。为了简便，这里先选定Reθc函数。其选取有一定的原则，保证转捩能够发生；在低雷诺数转捩区尽量大，在高雷诺数转捩区适当小。这里采用以下形式：

下面的主要任务是确定Flength函数，先选取一个平板算例，按以下几个步骤进行迭代：

（1）以常数Flength进行计算，并与试验结果比较，记下与试验吻合较好的Flength值；

（2）Flength值不变，代入常数Reθc进行计算，当获得与（1）相同结果时记录此时Reθc值；

（3）对比Reθc的表达式获得此时的eθt值；

最后以拟合公式对四个算例进行验证计算，如有必要，可对Reθc表达式进行适当调整。平板算例中，Flength函数只有在当地值下才对转捩进行作用，因此，拟合公式几乎不需要调整。

表2 平板算例关联函数值Table 2 Correlations for Reθcand Flengthas function ofeθt

表2 平板算例关联函数值Table 2 Correlations for Reθcand Flengthas function ofeθt

Case ~Reθt ~Reθc Flength T3B 107．136 102．2 37．01 T3A 205．816 178 31．74 T3A － 714．516 508．6 0．46 S＆K 927．066 685．1 0．40

图1 壁面摩擦阻力系数与试验的比较Fig．1 Comparisons of skin friction coefficient distribution between CFD calculations and wind tunnel tests

3 验证算例与讨论

3．1 A－aerofoil

A－aerofoil是验证转捩模型的标准算例之一。目前存在两套试验数据，ONERA F1风洞迎角13．1°和F2风洞迎角13．3°的试验数据，这里采取了与文献［5］相同的比较方式。计算参数：M＝0．15，Re＝2．07×106，α＝13．1°，13．3°。湍流度 FSTI＝0．2%，粘性比RT＝10。计算网格（图4）481×97（流向×法向），远场延伸100倍弦长，网格第一层1．0×10－5C，满足y＋≤1的要求。

数值模拟显示，翼型上表面在0．12C处出现层流分离，且导致分离诱导转捩，向下游发展为湍流，并在翼型后缘出现湍流分离。转捩模型计算得到升力为1．5741，阻力 0．0180，试验值分别 为 1．562、0．0208。全湍流1．4470、0．0282，升力偏小，阻力偏大。与文献［5］中全湍流阻力差别较大。全湍流结果在后缘与试验结果吻合更好些。

湍流项和转捩方程离散均采用二阶精度。图5比较了压力分布，湍流计算压力系数比试验明显偏低。图6中比较了湍流及转捩变量用一阶精度离散的摩阻分布，采用转捩模型大大提高了摩阻计算精度，同时摩阻分布显示一阶精度离散导致转捩位置略微提前。

图5 壁面压力系数曲线比较Fig．5 Comparison of wall pressure coefficient curve with expriment

图6 摩阻系数曲线比较Fig．6 Comparison of wall friction coefficient curve with experiment

3．2 NLR7301两段翼型

NLR7301两段翼型的风洞试验是在70年代末期在NLR3m×2m低速风洞中完成的，试验结果包括了总体气动特性、压力分布、典型站位的速度型分布等多种数据［9］。该翼型的襟翼偏角为20°、来流马赫数0．185，试验雷诺数2．51×106、主翼／襟翼重叠区均为5．3%c；试验构型包含了两种不同的缝隙宽度，一种为2．6%c，另一种为1．3%c。本文计算了主翼／襟翼的缝隙宽度2．6%c的构型，其中c为弦长。

转捩模型的应用就是为了提高阻力计算精度，文献［10］指出，高升力翼型的远场边界距离对气动特性、尤其是阻力特性影响显著，这要求对远场边界进行修正，或者扩大远场边界。然而对于本文使用的代码，远场边界取多大合适，这里首先采用环量修正远场边界条件在不同大小远场的网格下进行验证计算。

环量修正远场边界条件其实质为：由于翼型的存在，远场的自由来流受到了扰动，不再是均匀流动。远场的扰动可以近似用位于翼型压心点涡代替，由于阻力一般比升力小2个量级，因此阻力代表的点源不予考虑。点涡诱导速度为：

ξ、η为气流坐标系下坐标。从上式可以看到，诱导速度与升力成正比，与网格远场范围成反比，因此在升力较大（如高升力翼型），而远场不够远的范围下，诱导速度必须考虑。ξ、η表达式为：

代码中实际应用的远场条件通过变化得到：

压力、密度可以由等熵条件求得。

计算基准网格如图7所示，为2004年中国空气动力研究与发展中心与中国航空计算技术研究所联合举办的“CFD统一算例研讨”活动发布的标准网格。采用多块对接结构网格，网格分7块，网格单元数为15万，主翼布置457个点，襟翼布置249个点，尾迹区89个，壁面第一层网格距离1．0×10－5c，远场边界延伸12c。为研究网格外边界的影响，这里在基准网格基础上又做了7套网格进行对比计算。其他网格均是在基准网格基础上应用双曲型生成方法推进得到，网格径向间距放大因子取为1．1，网格生成过程始终保证了大网格内部网格与小网格完全一致，计算迎角为6°。

计算中采用Roe格式，应用4重多重网格及预处理，隐式离散用AFADI求解。翼型湍流度取 ，粘性比RT＝5。在网格外边界较大时，湍流强度可能衰减厉害，导致计算中转捩位置偏后，甚至可能出现气动力振荡。因此本算例应用了Spalart［11］建议的环境源项法来控制湍流度衰减，本文具体做法是使用物面距离来判断，即y＞10，控制源项作用，反之则使其自由衰减。也可以使用x坐标来控制。计算中应用了环量修正［12］来作对比计算。

图7 NLR7301计算基准网格及拓展网格Fig．7 Base grid and big grid of NLR7301

计算结果见表3、表4，对于无环量修正情况，阻力随着远场扩大单调减少，升力单调增加，达到6000C后收敛。从120C的结果来看，阻力距离最后收敛值还差了27count。应用标准网格计算的结果与最后收敛值相比大了接近100%。加入远场环量修正后，网格外边界达到120C后气动力就不再变化，且最终收敛值与无环量修正时完全一致。表中CDp为压力贡献，CDv为摩擦贡献。可以看到摩阻贡献在12C网格下就可以满足需要，阻力的主要差别在于压力的阻力分量。

因此，在高升力翼型设计中，计算网格远场边界至少应达到100C以上，尽量应用环量修正远场边界。

表3 无环量修正时计算的升力、阻力Table 3 Lift and drag without vortex correction

表4 有环量修正时计算的升力、阻力Table 4 Lift and drag with vortex correction

下面以300C的网格计算，加环量修正，研究转捩对阻力的影响，计算方法及条件同上，分别计算了6°、10．1°、13．1°状态。计算结果见表5。应用转捩模型后，阻力最大误差4%，全湍流计算最大误差达34．3%。本算例中，转捩计算提高阻力计算精度的同时，也高估了升力，这是以后应当考虑的问题。

表5 计算升力、阻力与试验比较Table 5 Comparison of lift and drag with experiment

表6给出了NLR－7301机翼上翼面、下翼面和襟翼上表面的转捩位置与计算的比较。本文计算中转捩位置的判定是以摩阻开始增加为依据。从表中可以看到，计算主翼上表面和襟翼上表面转捩位置与试验吻合较好，主翼下表面计算转捩位置靠前。图8给出了NLR7301不同迎角下的摩阻分布，从图上可以看到，计算摩阻分布转捩位置6°和13．1°的结果与试验转捩位置很接近，但是在计算中用摩阻系数第一次开始增加来判断不一定合理。以6°为例，摩阻系数从0．625C位置开始增加，但是增长极其缓慢，在0．68C处存在拐点，从这点开始才急剧增长。同时，图8结果表明，应用经验关联转捩模型最重要的是关联函数的标定，由于标定过程均在平板零压力梯度流动中进行，因此在实际外形流动的表现还有待改进。

表6 NLR－7301两段翼型转捩位置（GA P2．6%）Table 6 Transition location of NLR－7301（GA P2．6%）

图9的升力、阻力收敛历程比较可以看到，转捩与全湍流计算在多重网格迭代2000步后气动力就不再变化，转捩对收敛速度影响很小，适合工程应用。图10给出NLR7301迎角6°湍流强度等值线，图中白线部分为多块网格中的间隙。由于转捩模型采用了输运方程，因此湍流强度在网格间连续、光滑变化，相应的间歇函数在网格块间也光滑变化，这在一定程度上增强了代码的健壮性。

图10 NLR7301迎角6°湍流强度等值线Fig．10 NLR7301contours of turbulence intensity

4 小 结

（1）应用本文方法确定了一组关联函数，在平板测试中摩阻系数与试验吻合较好，说明了本文标定方法的正确性；

（2）A翼型的计算中，转捩结果的壁面压力系数及摩阻系数与实验的吻合程度远远好于全湍流计算，说明本文确定的关联函数在一定程度上能模拟转捩效应；

（3）NLR7301计算研究表明，要提高高升力翼型阻力计算精度，远场边界必须足够大（100C以上），尽量应用环量修正远场边界条件，同时应用转捩模型。

（4）转捩模型还很不完善，标定点太少，需要在今后应用中不断摸索。

致谢：感谢空气动力学国家重点实验室张玉伦、王光学、孟德虹同志的热情帮助。

［1］ 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．

［2］ MENTER F R，LANGTRY R B，LIKKI S R，et al．A correlation based transition model using local variables：part I－model formulation ASME－GT2004－53452［A］．Proceedings of ASME Turbo Expo 2004［C］，Vienna，Austria，2004．

［3］ LANGTRY R B，MENTER F R，LIKKI S R，et al．A correlation based transition model using local variables：part II－model formulation ASME－GT2004－53454［A］．Proceedings of ASME Turbo Expo 2004［C］，Vienna，Austria，2004．

［4］ LANGTRY R B，MENTER F R．Transition modeling for general CFD Applications in aeronautics［R］．AIAA Paper 2005－0522．

［5］ LANGTRY R B．A correlation－based transition model using local variables for unstructured parallelized CFD codes［D］．［Ph．D．thesis］，Stuttgart，2006．

［6］ MISAKA T，OBAYASHI S．A correlation－based transition models to flows around wings［R］．AIAA Paper 2006－918．

［7］ KARL PETTERSSON，SIMONE CRIPPA．Implementation and verification of a correlation based transition prediction method［R］．AIAA 2008－4401．

［8］ PAUL MALAN，KEERATI SULUKSNA．Calibrating the gamma－Re＿theta transition model for commercial CFD［R］．AIAA 2009－1142．

［9］ VAN DEN BERG B．Boundary layer measurements on a two－dimensional wing with flap［R］．NLR TR 79009U，1979．

［10］CHRISTOPHER L RUMSEY，SUSAN X YING．Prediction of high lift：review of present CFD capability［J］．Progress in Aerospace Sciences，2002，38：145－180．

［11］PHILIPPE R SPALART，CHRISTOPHER L RUMSEY，Effective in flow conditions for turbulence models in aerodynamic calculations［J］．AIAA Journal，2007，45（10）：2544－2553．

［12］THOMAS J L，SALAS M D．Far field boundary conditions for transonic lifting solutions to the Euler equations［R］．AIAA－85－0020．

［13］张玉伦，王光学．γ－Reθ转捩模型的标定研究［A］．工程湍流转捩学术研讨会［C］．海口，2010．

［14］孟德虹，张玉伦．γ－Reθ转捩模型在二维低速问题中的应用研究［A］．工程湍流转捩学术研讨会［C］．海口，2010．

［15］张晓东，高正红．关于补充Langtry的转捩模型经验修正式的数值探讨［J］．应用数学和力学，2010，31（5）：544－552．