基于SST-DDES方法的孤立气膜孔流动研究

王 鹏

(中国航发沈阳发动机研究所，沈阳110015)

1 引言

航空发动机气冷涡轮中的冷气/主流掺混过程直接影响涡轮的流场结构、流动损失和冷却效果。其中的复杂机理及规律已成为先进气冷涡轮设计中的一个重要问题和技术瓶颈，同时也是进一步深入挖掘气冷涡轮效率潜能的重要研究方向。

目前，研究手段和方法是制约冷气掺混研究的关键因素。对于真实复杂气动热力环境下的气膜冷却研究和设计优化，参数繁多且相互制约，相应流动机理的揭示还需要有详细可靠的流场数据支撑。用于设计优化结果验证的实验手段常常受到诸多限制，而数值手段的适用性则主要取决于其精度，目前制约冷气掺混数值计算精度的主要因素之一是湍流的预测。气膜冷却的各向异性湍流流场中存在着复杂的时空多尺度涡系结构，数值模拟必须考虑到所有相干结构的综合作用才能更精确地预测掺混流动，而这恰恰是RANS(雷诺平均)方法的力不能及之处。LES(大涡模拟)虽然可以较好地解决这一问题，但当雷诺数较大或几何较复杂时，对计算资源要求很高，且对壁面附近流动的模拟存在固有困难，要在气冷涡轮的系统性研究中应用存在很大困难。

在这种情况下，DES[1](分离涡模拟)类RANS/LES混合计算方法受到了广泛重视。其通过对流场特征的判断，在以耗散为主的区域采用RANS，以大涡输运为主的区域采用LES，并实现两者的自动过渡和切换，兼顾了计算的精度和效率[2]。然而DES类方法自建立以来主要用于外流中的大尺度分离及噪声等问题[3-6]，对内流中小尺度分离尤其是涡轮叶片气膜孔附近复杂流动的DES类研究相对较少[7-11]，且研究中采用的大都为DES97等原始的DES类方法。这些方法易产生雷诺应力模化不足即MSD(模型雷诺应力损耗)问题，进而带来网格诱导分离等非物理问题，而近年来改进的DES类方法如延迟的DES方法(DDES)等还少有应用。为此，本文采用基于SST湍流模型的DDES方法，从涡轮叶片气膜冷却孔附近流动的本质即横流射流问题入手，开展射流与主流掺混过程中的流动机理和损失机理研究。

2 数值方法

湍流模拟中采用基于SST湍流模型的DDES方法[12-13]。该DDES方法的构造基础是Strelets[3]将DES的基本思想与SST湍流模型相结合所得到的SST-DES方法，即将SST模型的湍动能输运方程耗散项改写为与网格尺度Δ和湍流长度尺度Lt相关的DES形式：

该SST-DES方法与早期的SA-DES等方法一样，存在着MSD问题和相应的网格诱导分离问题。为减小此类风险，Menter等基于SST模型的第一混合函数F1和第二混合函数F2对其进行改进：

式中：系数CDES取0.61。本文取FSST=F2，以最为严格地避免发生网格诱导分离。图1为改进前后FDES函数中第一项随混合函数变化特性对比。可见，改进后原SST-DES方法具有了延迟特征，成为一种DDES方法。

在基于有限体积法的数值计算中，空间离散格式的精度是难点之一。对于DES类方法，由于在不同流场区域同时存在着RANS和LES模式，若全场采用中心型格式则人工粘性系数将引入较多不确定性，若全场采用迎风型格式则耗散较大进而抑制小尺度结构的生成和发展。因此，引入了具有加权混合思想的混合型空间离散方法，在不同流动区域采用不同的空间离散格式[14-15]。对于编号为ip的空间某点处的通量Φip可表示为：

式中：σ为与流场有关的混合函数，Φip,U为中心型格式计算的通量，Φip,C为迎风型格式计算的通量。

该混合型空间离散方法具有以下特点：在LES模式流动区域趋近于0，表现为低耗散的中心型格式，利于解析小尺度流动结构；在壁面及远场无旋区域趋近于1，表现为耗散较大的迎风型格式，利于抑制数值振荡[15]。为进一步避免DES计算中发生网格诱导分离风险，对σ作如下附加限制：

式中：BF1和BF2为混合因子，本文选择第二混合函数，即BF1=0，BF2=1。为避免在中高CFL数情况下中心型格式带来数值振荡，对σ又施加了基于CFL数的限制器：

式中：CFLmax、CFLEXP 分别取5.0和1.0。

3 物理模型和计算设置

Ajersch等[16]对射流雷诺数Rejet为4 700的方孔横流射流开展了大量实验研究，获取了0.5、1.0和1.5三种吹风比[17]下的平均流参数和湍流脉动参数在流场不同位置处的空间分布，可用于数值模拟方法在此类流动中适用性的校验。本文以其物理模型为计算模型，采用SST-DDES方法研究了真实跨声气冷涡轮气膜冷却中较典型的0.5吹风比状态，同时采用相同的物理模型开展了基于SST湍流模型的RANS和URANS计算，以进行对比分析。

实验中有6个均布方孔，其排列方向和射流方向均与横流流动方向垂直，可实现多孔流动的周期性。数值模拟中选取一个周期作为计算域，图2为具体的计算几何模型。定义流向为x方向，展向为y方向，法向为z方向，射流孔中心为坐标原点。令射流孔边长D为特征长度，与实验模型一致取为12.7 mm，计算域流向长度共为51.0D，法向为25.0D，展向为3.0D。射流孔位置为 x=(-0.5～0.5)D，y=(-0.5～0.5)D，射流孔入口位置为 z=-5.0D，主流入口位置为x=-10.5D，出口位置为x=40.5D。

主流进口边界条件为给定总温、速度分布和湍流度分布，具体分布形式按照实验中的实测结果给出[18]，见图3；出口边界给定大气压作为背压。主流端壁为无滑移壁面，两侧面为对称边界，上边界为滑移壁面。射流进口给定总温、速度和中等湍流度，侧壁为无滑移壁面。主流和射流进口总温均为293.15 K。

计算域网格总数为300万，分27块，全场均为结构化网格。图4为计算域网格分布示意图，各处壁面的 y+均保证在1以下。DDES计算和URANS计算中的物理时间步长取0.000 15 s，虚拟时间步数取10步。非定常时间离散采用二阶欧拉后差，空间离散以准二阶迎风格式为基础，能量方程求解中采用热能方程形式。

4 结果分析

4.1 时均流场分析

图5为 y/D=0截面上x/D=0、1.0、8.0位置处流向速度U沿法向的分布，其中DDES和URANS的结果为时间平均值。速度采用射流出口速度Vjet进行无量纲处理，法向距离采用射流孔边长进行无量纲处理。由于复杂流动集中在壁面附近，因此主要针对法向约4.0D范围内的结果进行分析。

在x/D=0即射流孔中心位置，不同层次的数值模拟方法所预测的壁面附近流向速度分布几乎无差异，且都与实验结果吻合较好。这主要是由于该位置处于射流与主流相互作用的起始阶段，平均流的压力梯度和速度梯度占主导，流动的非定常性和掺混程度都较弱，RANS就可以实现较好的预测。从分布趋势看，射流对主流的影响主要集中在法向2.0D以下。值得注意的是，在法向2.0D以上区域，计算值大于实验值。这主要是由于考虑到计算资源的限制，对非研究重点但实际中存在的主流上壁面边界层不予模拟，进而使得主流的流向速度在该位置偏高，但这并不影响对下壁面附近流动特征的研究。在法向(0.2～0.4)D范围内，流向速度沿法向的梯度突然减小，实验由于测量点有限未能捕捉到该现象，该现象对应主流与射流之间的强剪切过程。

在x/D=1.0位置，流动的非定常性和掺混程度加强，不同数值模拟方法预测结果间的差异开始显现。相比而言，DDES方法预测结果与实验结果吻合得更好。主要差异集中在法向0.6D以下范围内，URANS和RANS方法均过高地预测了射流尾迹中回流流动的程度，特别是RANS方法。在法向约0.8D位置，也存在一个流向速度法向梯度突变的现象，这与射流内侧边缘的剪切流动相对应。

在x/D=8.0位置，经过一定距离的流动掺混后，DDES的预测精度明显更高，主要体现在法向1.2D以下射流主体区域，URANS和RANS预测的流向速度均明显偏小。

图 6为 y/D=-0.5截面上 x/D=0、3.0、5.0位置处展向速度W的分布。在x/D=0位置，各方法的预测结果仍基本无差异，且均与实验结果吻合较好。该位置处于射流孔侧面边缘，射流从孔内喷出后，孔内的边界层突然变为自由剪切层，在K-H不稳定性的作用下，剪切层失稳进而卷起一对反向旋转的旋涡。图中展向速度在近壁面附近的速度峰即对应该旋涡的起始发展过程。

在x/D=3.0和x/D=5.0位置，反旋涡对进一步发展，强度减弱且相互融合，DDES预测结果与实验结果吻合得更好，URANS和RANS均过高预测了反旋涡对上半区的强度。在反旋涡对下半区，各方法均预测出了近壁面附近的速度峰，表征了反旋涡对的强烈下洗过程，但实验结果受测量手段影响未能捕捉到该速度峰。值得注意的是，在x/D=3.0位置，实验结果在法向0.8D附近存在两个较弱的速度峰，URANS和RANS均未捕捉到该特征，而DDES实现了对该特征的捕捉。整体而言，DDES在展向速度方面的预测精度相对较高，但在局部也与实验结果存在一定差异。这除了表明DDES方法自身还需进一步改进外，也与计算域受计算资源限制仅截取了单个流动周期有关。

图7为y/D=-1.0截面x/D=0位置处法向速度V的分布。在射流影响下，展向1.0D位置处的主流也具有显著的法向速度，在法向1.0D附近达到峰值；DDES的预测结果同样与实验结果吻合得较好，其他两种方法的预测结果均偏小。

图8为y/D=0截面x/D=8.0位置处无量纲湍动能的分布。由前文可知，在法向1.2D附近区域，∂U/∂z较大，URANS和RANS高估了该区域的湍动能峰值，DDES则较好地预测了湍动能的峰值水平。但应注意，实验测量结果中在壁面附近还存在另一个局部湍动能峰值，而各方法计算结果均未捕捉到该峰值，这和计算预测的速度场等与实验结果存在一定差异有关。

图9为y/D=0截面x/D=3.0、8.0位置处雷诺切应力的分布。在法向1.0D附近即射流迎风面与主流发生强剪切的区域，URANS和RANS捕捉的局部应力峰值显著过高，这与图8中计算涡粘系数的湍动能局部峰值预测过高相对应，而DDES则实现了相对更为合理的捕捉。对于该现象，Andreopoulos等[19]的研究认为，RANS方法的涡粘模型在与∂U/∂z正负相反的流动区域理应具有一定的预测能力，但当存在流线汇聚或发散以及∂W/∂y的影响较显著等情况时其预测能力较弱，该区域附近的流动恰好由于主流绕射流的流动而存在这些情况。

图10进一步给出了DDES方法预测的射流孔出口法向速度(以Vjet进行无量纲)等值线与实验结果对比。在所研究的0.5吹风比状态下，射流孔出口流动受主流强烈影响，从射流迎风侧到背风侧法向速度逐渐增加，最高值约1.7Vjet，上游半区的梯度相对较大。图11为DDES方法预测的射流孔出口平面流线分布与实验结果中速度矢量分布的对比，绝对速度沿流向逐渐增加反映了射流在流向的偏转程度，同时流动也向两侧偏转，且其偏转程度随距射流孔中心距离的增加而增加。从图10和图11的对比结果看，DDES方法预测的主要特征均与实验结果吻合较好。

图12为y/D=0截面射流孔附近DDES方法预测的时均流场流线分布。在主流作用下，射流还未到射流孔出口便开始发生弯曲和流管收缩，流速相应增加，射流主体逐渐偏转，最终与主流方向趋于平行，在射流迎风侧和背风侧均存在回流区。由于吹风比较小，随着向下游的发展，射流速度的法向分量迅速减小，即射流的穿透能力较弱。

图13为射流孔背风侧边缘即x/D=0.5截面DDES方法预测的时均流线分布。在y/D=±0.5的射流孔两侧区域卷起了两个旋向相反且具有良好对称性的旋涡，其核心区域的流速高于主流流速，且该反旋涡对在该位置上游便已形成并逐渐变大。图14为x/D=2.0截面时均流场流线分布。反旋涡对离开射流孔区域后向下游继续发展，尺寸不断增加，相互之间的影响逐渐加强，各自的强度则逐渐减弱，旋涡核心区的流速已显著低于主流区。

图15为射流孔附近主流壁面的极限流线分布。射流的存在对主流流动带来阻碍，主流绕射流主体流动，可近似看作方柱扰流流动，在前缘存在马蹄涡系，在扰流后方存在回流区。但射流不同于固体，自身可以变形且带来自由剪切层，使得该区域涡系结构更加复杂。

4.2 瞬态流场分析

针对某时刻瞬态流场中射流孔附近的流场结构进行了基于Q法则的旋涡识别，Q值取320 000，同时采用流向涡量对涡系结构进行渲染，见图16。由图可知，瞬态流场中的大涡拟序结构，主要包括射流前缘稳定的马蹄涡、射流侧壁产生的反旋涡对以及射流后逐渐形成并发展的发卡涡结构等。在本文较低吹风比状态下未见高吹风比状态下存在的从射流前缘开始随射流向下游移动的反向涡头[20]。从大涡结构上的渲染状态，可清楚看到反旋涡对的反向旋转特征及前缘马蹄涡两个分支的反向旋转特征，下游的发卡涡无论是平行涡腿还是垂直涡腿也都具有两侧分支旋转方向相反的特征，同一侧的马蹄涡分支和反旋涡对分支的旋转方向也彼此相反。流向涡较强的区域主要集中在反旋涡对和发卡涡的平行涡腿等位置。

图17进一步给出了采用展向涡量和法向涡量渲染的涡系结构。马蹄涡和发卡涡的涡头等位置具有相对较高的展向涡量，且方向相同。法向涡量较强的区域则主要集中在发卡涡的垂直涡腿等位置，且也具有两侧方向相反的特征。图18为某瞬时x/D=1.0截面的熵分布，可看到流动损失主要来源于包含反旋涡对在内的射流背风侧下游区域。

5 结论

采用基于SST湍流模型的DDES法，对典型吹风比下的孤立气膜孔流动进行了数值模拟，并与实验测量结果及RANS/URANS结果进行了详细对比。分析了包括各速度分量、湍动能、雷诺应力和大涡拟序结构等在内的流场信息。研究结果表明：

(1)在吹风比0.5的条件下，横流射流的入射能力较弱，进入主流后迅速偏转，对下游壁面形成覆盖，其影响区域主要集中在法向2.0D范围内；同时，射流还存在一定的展向扩展，但在方孔边长为相邻孔距离1/3的条件下，相邻射流孔流动不存在明显干涉。

(2)掺混流场中存在着射流前缘马蹄涡、产生于侧壁且具有强烈下洗特征的反旋涡对以及包络其上的发卡涡等大涡拟序结构，且这些大涡结构直接影响平均流特征和气膜覆盖效果，同时也是掺混损失的主要来源。

(3)对于横流射流类复杂掺混流动，与传统的RANS/URANS方法相比，基于RANS/LES混合的SST-DDES方法不仅对时均流场的预测更合理、更准确，还可有效捕捉瞬态大涡拟序结构及其演化过程，更利于流动机理研究和精细化设计。

(4)与实验结果相比，在射流影响区域，SST-DDES方法预测的时均结果仍存在流向速度亏损、展向扩张程度、湍动能和雷诺切应力偏大等不足，SST-DDES方法仍需适当改进以进一步提高对横流射流流动的预测精度。

参考文献：

[1]Spalart P R，Jou W，Strelets M，et al.Comments on the fea⁃sibility of LES for wings and on a hybrid RANS/LES ap⁃proach[J].Advances in DNS/LES，1997，1：4—8.

[2]孙明波，汪洪波，梁寒剑，等.复杂湍流流动的混合RANS/LES 方法研究[J].航空计算技术，2011，41(1)：24—33.

[3]Strelets M.Detached eddy simulation of massively separat⁃ed flows[R].AIAA 2001-0879，2001.

[4]Squires K D，Forsythe J R，Morton S A，et al.Progress on detached-eddy simulation of massively separated flows[R].AIAA 2002-1021，2002.

[5]Forsythe J R，Squires K D，Wurtzler K E，et al.De⁃tached-eddy simulation of fighter aircraft at high alpha[R].AIAA 2002-0591，2002.

[6]Rulik S，Dykas S，Wroblewski W.Modeling of aerodynam⁃ic noise using hybrid SAS and DES methods[R].ASME GT2010-22696，2010.

[7]Roy S，Kapadia S，Heidmann J D.Film cooling analysis us⁃ing DES turbulence model[R].ASME GT2003-38140，2003.

[8]Kim S I，Hassan I G，Zhang X Z.Unsteady simulation of film cooling flow from an inclined cylindrical jet[R].ASME HT2007-32401，2007.

[9]Funazaki K，Kawabata H，Takahashi D，et al.Experi⁃mental and numerical studies on leading edge film cooling performance:effects of hole exit shape and free-stream turbulence[R].ASME GT2012-68217，2012.

[10]梁俊宇，康 顺.基于DES方法的平板孤立冷却孔数值模拟研究[J].工程热物理学报，2011，32(3)：395—398.

[11]翟丽娜.基于DES方法的气膜冷却数值模拟[D].北京：华北电力大学，2012.

[12]Menter F R，Kuntz M.Development and application of a zonal DES turbulence model for CFX-5[R].CFX Valida⁃tion Report CFX-VAL17/0503.

[13]Menter F，Kuntz M.A zonal SST-DES formulation[R].Russia：St.Petersburg Technical University，2003.

[14]Travin A，Shur M，Strelets M，et al.Physical and numeri⁃cal upgrades in the detached-eddy simulation of complex turbulent flows[J].Advances in LES of Complex Flows，2004，65：239—254.

[15]刘 健.采用DES类方法研究基本起落架非定常大范围分离流动[D].北京：清华大学，2011.

[16]Ajersch P，Zhou J M，Ketler S，et al.Multiple-jets in a cross-flow:detailed measurements and numerical simula⁃tions[J].Journal of Turbomachinery，1997，119：330—342.

[17]Holdeman J D，Walker R E.Mixing of a row of jets with a confined cross-flow[J].AIAA Journal，1977，15(2)：243—249.

[18]Zhou J M.A computational and experimental investigation of gas turbine blade film cooling effectiveness[D].Univer⁃sity of British Columbia，1994.

[19]Andreopoulos J，Rodi W.Experimental investigation of jets in a cross-flow[J].Journal of Fluid Mechanics，1984，138：93—127.

[20]刘 斌，任丽芸，叶 建，等.方孔横向射流的大涡模拟研究[J].航空科学技术，2009，(2)：40—43.