动力喷流效应混合方法数值模拟研究

龚志斌, 李杰, 单继祥, 张恒

(1.中国工程物理研究院 总体工程研究所, 四川 绵阳 621900; 2.西北工业大学 航空学院, 陕西 西安 710072)

发动机喷流是动力推进、气动噪声、燃烧以及低可探测技术的关键部分,对其流场结构更好的理解和更为准确的预测手段对于发展高效、静音和环保飞行器具有重要的意义。对短距起降(STOL)飞行器而言,动力喷流和增升装置之间的气动干扰尤为强烈,流动现象更为复杂,噪声问题更加突出,对动力喷流流场进行准确模拟,可为STOL构型全机气动设计和噪声预测奠定良好基础。

当前航空领域对于喷流流场的预测主要采用RANS方法,但经过几十年发展,RANS方法仍然不能够为喷流流场提供可靠的预测,标准湍流模型仍然无法准确预测剪切层增长速率和不受剪切层影响的喷流无黏速势核心长度。高温条件下,喷流流动中湍流Prantl数从近场到远场变化剧烈,喷流附近区域Prantl数低,掺混现象明显,模拟起来更加困难。

国外不少学者采用LES方法对喷流湍流流场进行数值模拟研究。Pope[1]指出LES方法可以很好地模拟自由剪切流动中的动量、热量和质量交换过程。喷流发动机喷口壁面所产生的湍流影响通常被忽略或者假定很小,但Brunet[2]指出它们对于喷流自由剪切层的影响非常显著,当发动机轴向长度减小或者涵道比增大时喷口壁面的湍流影响还会增大。在应用 LES方法求解包含发动机壁面的湍流问题时,为了模拟壁面附近很小尺度的运动,边界层内所需的网格量几乎与DNS方法相当,准确模拟边界层流动所需的计算资源也很庞大。因此,LES方法主要用于简单的喷流流动和低雷诺数情形,还不能应用于工程实际中的复杂发动机喷流流动。此外,LES方法近壁面及可压缩亚格子应力模式的发展还不够完善,对大多数包含壁面的高雷诺数喷流流动模拟能力不足。同时LES方法喷流预测结果的调查统计[3]表明,喷流中心线上的速度和湍流强度等信息的预测非常不理想,其主要原因在于计算网格质量不高且求解方法的格式离散精度不足。

近年来RANS/LES混合方法逐步发展并应用于动力喷流流场的数值模拟,其基本思想是在喷管壁面附近采用高效的RANS方法来对边界层内的小尺度湍流进行模拟,下游喷流中的大尺度湍流结构则采用LES方法进行模拟,基准RANS湍流模式和LES区域亚格子应力模型在很大程度上决定了RANS/LES混合方法特性。结合RANS和LES方法的各自优点发展而来的混合方法是兼顾模拟精度和计算效率的最优选择,它可以认为是LES方法处理高雷诺数喷流流动的可行方法。Shur和Spalart等人[4]针对复杂喷管喷流流动开展了RANS/LES混合方法数值模拟研究,Cai等[5]基于该方法结合高阶格式讨论了网格分布、边界条件处理等对喷流流动数值模拟的影响。Brunet等[6]基于Zonal-DES方法,针对较为实际的机翼/挂架/短舱构型,开展了跨声速情况下的动力喷流数值模拟研究,基于简单的SA、SST湍流模型和复杂的DRSM湍流模型的混合方法都取得了与试验较为吻合的结果。国外还有一些其他研究人员,如Grinstein和Fureby等[7]认为不使用亚格子模式的隐式LES方法(ILES方法,又称NLES方法)能够较好地捕捉复杂喷流湍流结构,Eastwood等[8-9]采用RANS/NLES方法对简单轴对称喷管和发动机复杂喷流进行模拟并得到良好效果。从公开的文献资料来看,国内采用RANS/LES混合方法对动力喷流流场数值模拟研究成果仍相对较少。

本文拟开展动力喷流效应混合方法分析及应用研究。基于自由发展喷管流动和ARN2喷管流动算例,分别讨论网格密度和空间离散格式精度对混合方法数值模拟结果的影响。在此基础上,对较为实际的涡扇发动机喷流流场进行混合方法模拟分析。

1 数值模拟方法

在有限体积法基础上,对三维可压缩非定常N-S方程进行求解。对无黏项采用三阶MUSCL格式,黏性项采用二阶中心差分格式进行空间离散;时间推进采用隐式LU-SGS格式。为了提高空间离散格式精度,采用Jiang等[10]提出的有限体积五阶WENO-JS格式进行左右状态变量重构来实现。

采用Shur等[11]的方法,在k-ωSST湍流模型的基础上,结合DDES方法和LES壁面模型的特点,重新定义亚格子尺度,解决边界层附近对数边界层不匹配问题,加快RANS方法到LES方法的转换,实现了SST-IDDES方法。其中湍动能输运方程中引入长度尺度lIDDES,耗散项为

Dk=ρk1.5/lIDDES

Δ=min[max(CwΔmax,Cwd,Δmin),Δmax)]

Cw=0.15,Δmax和Δmin分别为网格单元三向最大和最小尺度,d为网格单元与壁面距离。

发动机进排气边界条件有不同的构造方法。其中一种将进气口作为流场出口,给定质量流量,或者是静压和总温,排气口作为流场入口,给定质量流量和总温条件。但是这种方法在整个网格边界面上需要指定均匀的质量流量,会影响计算过程收敛性和计算结果精度。本文根据自由来流条件、进气口质量流量和喷口总温总压确定发动机工作状态,以总压比(Po/P∞)和总温比(To/T∞)统一定义进排气动力边界条件,其中Po和To为当地总压和总温,P∞和T∞为自由来流静压和静温。在风扇入口指定质量流量比,在已知自由来流条件下,由能量守恒关系式和等熵关系式,可迭代求出风扇入口马赫数,进而得到风扇入口总压比和总温比,内外涵道尾喷出口的总压比和总温比在已知涵道比的情况下可以类似得到,详细推导参见文献[12-13]。

2 自由发展喷管流动计算验证分析

Panda等[14]在NASA的Glenn研究中心对自由喷管流动进行了试验测量,文献[15]给出了相应的喷管直径、试验环境和喷管出口流动条件。计算过程中除喷管出口给定排气边界外,喷管其他物面边界给定滑移壁面边界。试验条件下出口马赫数为1.4,对应喷管出口边界总压比Poj/P∞=3.192、总温比Toj/T∞=1.012 4,其中Poj和Toj为喷管出口当地总压和总温。自由来流马赫数为0,为了使得计算稳定,给定来流马赫数为0.01。计算过程中给定雷诺数1×106。

为准确计算喷口下游的动力喷流,对喷口下游15D(D为喷管出口直径)范围内的流场网格进行了适当加密,生成了粗、细2套计算网格,喷流区域内粗网格3个方向上最大尺寸为0.08D,网格规模约为700万;细网格流向最大尺寸为0.04D,网格规模约为2 100万,网格空间截面如图1所示。

在非定常RANS方法计算结果的基础上采用结合了三阶MUSCL插值格式的SST-IDDES方法进行续算,无量纲时间步长Δt取为0.000 5。

图1 喷管网格拓扑结构

图2给出了喷管下游中心线上的时均速度与试验结果的对比。X/D(喷管出口下游位置X与喷管直径D之比)小于8的范围内,时均速度接近于波动范围内的平均值。X/D大于8时,时均速度结果与试验存在偏差,网格密度增加有利于结果改善。从图3不同位置速度分布结果来看,粗网格情况喷口下游X/D=2和X/D=4处速度峰值与试验接近,但是计算得到的喷流直径比试验情况小,当地剪切层速度梯度很大。网格加密后,X/D=2处速度分布没有变化,X/D=4和X/D=6位置速度分布结果得到较大改善。

图2 时均流场喷管下游中心线上速度计算结果

图4给出了喷流流向截面内的瞬态马赫数云图。从粗网格结果来看,在喷口下游保持了较长一段相对稳定的流动状态,然后高速喷流与周围静止气流在黏性的作用下,产生了剪切层失稳。当网格加密后,剪切层失稳更早出现。结合图4结果来看,喷口下游出现的较长流动稳定状态与物理情况不符,网格加密在一定程度上能够更加逼近实际情况。图6给出了以马赫数渲染的流场瞬态Q等值面图,网格密度的影响体现得更为明显,在细网格的基础上可获得喷流流场更多、更丰富的三维小尺度湍流结构。

图4 瞬态速度云图计算结果

图5 流场瞬态Q等值面图

3 ARN2喷管流动计算验证分析

NASA的Glenn研究中心设计了一系列亚声速流动喷管ARN(acoustic reference nozzles),并在此基础上进行了大量的喷流及噪声试验,为CFD数值验证提供了有用的数据,是验证喷流湍流流动的标准算例之一[16]。ARN系列收缩喷管为轴对称外形,其中 ARN2出口直径为5.08 cm。

针对ARN2喷管外形与流动特点,建立与上一节类似的三维轴对称柱形计算域,计算域大小为50D×80D(D为喷管出口直径)。喷管位置与试验保持一致,喷流入口位于X=0处,上游边界距离喷管出口30D,喷管外壁面向上游延伸直到上游边界,下游边界距离喷管出口的距离为50D,周向边界到喷管中心轴线的距离为25D。采用点对接多块结构化网格技术生成流场结构网格,在喷管出口下游20D范围内对网格适当加密,壁面附近网格尺度为1×10-4D,唇口厚度方向布置33个网格点,在出口下游20D范围的关注区域内网格最大尺度不超过0.04D,网格规模为2 200万。

在非定常RANS方法计算结果的基础上采用结合了三阶MUSCL格式和五阶WENO格式的SST-IDDES方法进行续算。在喷管壁面给定无滑移边界条件,延伸的外壁面给定滑移边界条件,在周向和上、下游边界给定黎曼远场边界条件。不考虑温度效应,在喷流入口给定总温比To/T∞=1.0,总压比Po/P∞=1.8,对应喷管出口的马赫数为0.9。给定来流马赫数为0.01,试验雷诺数为2.0×106。

图6 瞬态马赫数云图

图7给出了计算得到的瞬态马赫数云图。喷流入口到喷管出口之间的区域内,流动逐步加速到马赫数0.9,涡量几乎为零,流动是稳定的。采用三阶MUSCL格式计算时,喷管出口下游保持了较长的流动稳定状态,然后才由于黏性和速度差导致剪切层失稳,湍流混合作用增强。采用五阶WENO格式计算时,在唇口下游立即出现较大尺度的流场结构,这些流场结构的大小与唇口厚度处于一个量级,喷流流场中的涡结构更加细致,涡量耗散要小得多。图9给出了瞬态Q等值面结果,相对于三阶MUSCL格式而言,采用五阶WENO格式计算时整个喷流流场小尺度涡结构要更加丰富细致,在唇口附近表现尤为明显。

图7 瞬态Q等值面喷流流场结构图

图8给出了时均流场喷流中心线上的速度分布。在X/D小于5的速势中心区域内,不同离散格式的结果基本相同。X/D大于5时,相对于三阶MUSCL格式而言,采用五阶WENO格式计算得到的速势中心区域要更长,中心线上的速度衰减速率要小得多,量值上与试验结果更为贴近。

图8 喷流中心线上的时均速度与试验结果对比

图9给出了时均流场喷口下游X/D为1,4,10的流向速度沿轴向的分布。速度分布结果与试验较为接近,但不同离散格式的结果还是存在一定差异。在X/D=1位置,流向速度峰值与试验几乎一致,三阶MUSCL格式结果在y/D=±0.5位置处速度存在阶跃,y/D=±0.5时速度很快减小并接近为零,这意味着当地剪切层厚度比实际情况要小得多,五阶WENO格式结果则不存在阶跃,与试验结果较为接近。在X/D=4位置,2种格式均能够获得较好的结果。在X/D=10位置,速度峰值相对于试验结果偏低而速势中心区以外的速度偏高,尤其是三阶MUSCL格式的结果差异增大。

图9 时均流场流向速度分布与试验结果对比

为了定量的表示喷流与周围气流的掺混,定义剪切层厚度参数b为

b=r|u=0.1Uj-r|u=0.9Uj

式中，r为与喷流中心线的轴向间距。图10给出了X/D小于8范围内的剪切层厚度分布,速势中心区域基本在此范围。试验结果显示剪切层厚度基本上呈线性增长。三阶MUSCL格式得到的剪切层厚度在唇口附近接近于零,呈类似抛物型增长,在X/D大于4以后呈线性增长且增长速率要明显大于试验结果。五阶WENO格式得到的剪切层厚度比试验结果略高,总体增长趋势基本一致。

图11给出了喷流中心线上的脉动速度(即湍流强度)分布。中心线上的脉动速度在唇口附近最小,往下游逐渐增大,在达到一定峰值后慢慢下降。三阶MUSCL格式和五阶WENO格式计算得到的中心线上脉动速度分布与试验结果趋势上基本一致。三阶MUSCL格式得到的脉动速度峰值出现较早,X/D在5到10的范围内要高于试验结果,当X/D大于10时要明显低于试验结果。五阶WENO格式计算得到的中心线上脉动速度与试验更为接近。

图10 剪切层厚度分布与试验结果对比 图11 喷流中心线上脉动速度与试验结果对比

4 涡扇发动机动力喷流数值模拟分析

下面针对较为实际的涡扇发动机模型,采用SST-IDDES混合方法,对动力喷流效应进行数值模拟分析。计算网格在喷流区域内适当加密,在发动机下游3倍发动机长度范围内网格尺度不超过发动机最大直径的2%,网格总数约为3 500万。取海平面标准大气条件,来流状态为:Ma=0.2,Re=4.66×106,发动机轴线与来流平行。发动机动力边界条件为:质量流量比RMF=1.6,涵道比RB=6.0,风扇入口总压比Po/P∞=1.028 3,总温比To/T∞=1.008;风扇喷流出口总压比Po/P∞=1.8,总温比To/T∞=1.2;涡轮喷流出口总压比Po/P∞=1.6,总温比To/T∞=2.8。在非定常RANS方法计算结果的基础上采用结合了三阶MUSCL格式和五阶WENO格式的SST-IDDES方法进行续算。

涡扇发动机由于存在内外涵道,在出口下游的喷流中会产生内、外剪切层,内侧剪切层源于内涵喷流与外涵喷流之间的黏性作用,外侧剪切层源于外涵喷流与周围环境气流之间的黏性作用,从图12瞬态密度分布和图13瞬态涡量分布结果来看，内侧剪切层相对于外侧剪切层更加稳定,外侧剪切层更早出现失稳现象并诱导内侧剪切层失稳,从而导致内外涵道喷流掺混。图14给出了以马赫数渲染的发动机喷流瞬态Q等值面结果,图中更好地显示出采用五阶WENO格式所解析出的发动机出口下游丰富清晰的三维小尺度湍流结构和喷流的影响范围。

图12 子午面瞬态密度分布云图

图13 子午面瞬态涡量云图

图14 瞬态Q等值面结果

5 结 论

本文从动力喷流效应精细化模拟需求出发,开展了混合方法数值模拟分析与应用研究。在k-ωSST湍流模型基础上建立了IDDES混合方法,在此基础上,针对自由发展喷管流动和ARN2喷管流动开展了模拟验证分析,讨论了网格密度和空间离散格式精度对喷流流场计算结果的影响,并对较为实际的涡扇发动机动力喷流流场进行了数值模拟分析。

自由发展喷管流动混合方法计算结果表明:增加网格密度可有效改善喷流速度分布预测结果,缩短喷口下游非物理流动稳定状态,瞬态流场结构更加清晰合理。ARN2喷管流动混合方法计算结果表明:相比于三阶MUSCL格式,采用五阶WENO格式计算得到的喷流湍流结构更加丰富细腻,喷口下游流场更加合理,中心线上的速度衰减速率要更小且与试验更为吻合,喷流区内的速度分布与试验更为接近,湍流强度尤其是喷口下游附近的湍流强度预测更为准确。

涡扇发动机喷流流场混合方法数值模拟结果表明,外侧剪切层更早出现失稳现象并诱导内侧剪切层失稳,从而导致内、外涵道喷流掺混。提高空间离散格式精度,喷流掺混更早更充分,湍流结构更加细致合理。