基于大规模并行计算的结冰翼型失速流场特性数值模拟研究

李立, 武君胜, 梁益华, 田增冬

(1.西北工业大学 计算机学院, 陕西 西安 710072;2.中国航空工业集团公司西安航空计算技术研究所 预研部, 陕西 西安 710068)

飞机结冰安全性是运输类飞机设计评估中必须关注的重要课题。美国联邦航空管理局(FAA)和中国民用航空总局发布的适航条款中均有结冰安全性评估的相关条款和内容[1-2]。

在过去数年,随着计算机科学、计算数学等相关学科发展,计算流体力学(CFD)技术得到长足发展,采用CFD方法开展飞机结冰安全性的预评估已越来越普遍。CFD已逐步成为评估飞机飞行性能[3-7],及深入理解结冰流动机理的强有力工具[8-10]。目前,业界形成的普遍共识是,尽管在实际工作过程中,精细CFD仿真存在非常耗时的问题,但与试验相比,投资回报巨大。而且,一般来说,开展全尺寸(或接近全尺寸)结冰试验的费用非常昂贵,并且在飞机设计完成之前几乎不可能开展[11]。

目前,可用于飞行器气动性能评估的CFD方法主要包括直接数值模拟(DNS)、大涡模拟(LES),以及雷诺平均Naiver-Stokes方程(RANS)方法等。其中,RANS方法由于其实现了在计算资源及预测准确性上的有效平衡,在过去数十年间,得到广泛青睐和应用,对工程设计产生了巨大影响[12]。然而,众所周知,以现代工程湍流模型为基石的RANS方法往往无法准确预测具有大范围分离的复杂湍流流场,存在模型部分失效的问题[13]。DNS和LES方法在理论上可应用于上述情况,但由于需要巨大的计算成本,迄今还难以直接应用到复杂几何和全尺寸工程问题的计算。

为此,近年来人们致力于发展RANS-LES混合方法,以弥补RANS与LES(甚至DNS)在预测大范围分离复杂流场时的能力差距。RANS-LES混合方法是包括混合了不同RANS和LES模型的一大类方法,主要针对复杂湍流问题[4-5,8-10]。RANS-LES混合方法的独到之处是,将计算代价降低到非常接近时域非定常RANS(即URANS)计算的程度,但仍可以提供比RANS更准确的方法来解析关键流场区域的高度分离湍流流场结构。目前,RANS-LES混合方法正成为在RANS模型失效情况下、替代用于复杂工程设计的一种有前景方法。

综合当前技术发展来看,构造RANS-LES混合方法的方式主要包括2种:①从RANS出发、自下而上的方法,属于早期普遍采用的方法。典型代表是基于Spalart-Allmaras(SA)方程湍流模型发展的脱体涡模拟(DES)方法,自提出以来即得到广泛跟踪和研究[14-15];②从LES出发、自上而下的方法,典型代表是壁面模化大涡模拟(WMLES)方法,这是近年逐步引起重视的一大类方法。在美国NASA和波音公司发布的CFD 2030发展愿景中,WMLES方法被认为是“面向复杂整机应用的下一代LES方法”[13]。

然而,RANS-LES混合方法作为一种精细CFD方法,与常规RANS计算相比,不单必须采用长时非定常计算(以建立稳定的瞬态变化流场),而且,由于需要解析空间流场的复杂漩涡结构(特别是三维横向涡结构),所需的计算网格数量往往需要增加数十倍甚至上百倍。这意味着计算过程中,引入大规模并行计算以加速计算过程非常关键。

由于结冰翼型失速流场预测的复杂性,目前针对该问题,采用RANS-LES混合方法进行精细分析的工作并不多见。近年来比较具有代表性的工作包括ALAM等[8]提出的一种变种的DDES方法(即DHRL)、XIAO等[9-10]基于自主程序发展的以SST为基底模式的IDDES方法等。这些方法均是采用自底向上的传统方法进行RANS-LES方法构造。本文针对该问题,提出发展结合大规模并行计算和壁面模化大涡模拟(WMLES)的有效计算策略和方法。通过针对双角冰结冰翼型GL305/944的后失速流场特性数值模拟研究,并与以两方程SST模型为基底模式的RANS计算及IDDES计算开展细致的结果对比,证实了本文发展方法的有效性,可为基于数值计算的飞机结冰安全性评估提供有效技术手段。

1 数值方法

1.1 流场控制方程

在三维笛卡尔坐标系下,Navier-Stokes方程写成守恒形式如下

(1)

Q=(ρ,ρu,ρv,ρw,ρE)T

ρUiw+pδi3,(ρE+p)Ui)T

式中:ρ,U=(u,v,w)T,E,p,τij,qi分别表示密度、速度、总内能、压强、剪切应力张量和热流通量。剪切应力张量按照Boussinesq假设,可表示为

(2)

压强可根据状态方程计算,对理想气体有

(3)

式中:γ为比热容比,对空气,γ=1.4。

对方程(1),本文采用非结构混合网格策略和有限体积法进行数值求解。目前,利用非结构网格进行有限体积离散主要有2种方式[16]:①直接采用网格单元作为控制体,称为格心(cell-centered)格式;②以网格单元顶点为中心,采用对偶网格方法建立控制体,称为格点(cell-vertex)格式。本文选用格点格式。

值得指出,对RANS-LES混合方法来说,由于涉及流场中小涡结构的捕捉,数值格式的耗散控制是必须解决的一个重要问题。为此,在实际数值求解中,对空间离散,本文对流项引入一种全速域自适应低耗散AUSM+格式进行计算[17],具有对边界层及全流场计算自动保持低耗散及高性能的特点;黏性项则采用标准中心差分格式进行计算。对时间迭代,为了获得瞬态流场解,采用改进的2阶精度后向差分(backward differentiation formula,BDF)格式及双时间步方法进行非定常计算。

1.2 湍流模型

为了实现湍流模拟,除了WMLES,同时采用以两方程SST模型为基底模式的RANS和IDDES方法进行对比计算[18-19]。本文WMLES模型的具体定义公式为

(4)

式中:K=0.4为卡门常数;CSMAG=0.18为LES亚格子模型常数,而

fD=1-exp(-y+/A+)A+=26

(5)

为普朗特混合长度模型系数,Δ为计算网格的长度尺度,简单定义为

Δ=max(Δx,Δy,Δz)

(6)

依据模型(4)式,当Ky≪CSMAGΔ,即在在黏性固壁附近,会使用与代数RANS模型相当的普朗特混合长度假设;而当Ky>CSMAGΔ时,表示远离黏性固壁,则采用经典Smagorinsky亚格子LES模型。由此不难分析,上述WMLES模型本质上可以看做是经典Smagorinsky亚格子LES模型和普朗特混合长度理论的一种混合模型[19]。采用上述方法,事实上建立了一种避免在黏性子层物理解析壁湍流的方法,从而建立了一类典型的、自顶向下的RANS-LES混合方法。

以两方程SST模型为基底模式的IDDES模型完整计算公式可在文献[18]中找到。作为一种典型的脱体涡(DES)模型,其主要设计思想是将SST模型中湍动能方程的RANS长度尺度替换为IDDES混合长度尺度,从而使该模型表现为RANS和LES方法混合的一种行为。

修正后的湍动能方程为

(7)

式中,IDDES混合长度尺度,lIDDES,定义为

(8)

fd=1-tanh(8rd)3是一个范围从0.0到1.0的混合函数,用来实现模型在RANS(1.0)和LES(0.0)之间的自动切换。

1.3 基于大规模并行的加速计算

需要指出,无论是从RANS出发的DES类方法,还是从LES出发的WMLES方法,本质上都是一种基于时间准确性(瞬态)计算的方法。因而,当使用WMLES或DES方法进行CFD仿真时,不仅必须采用长时非定常计算,以建立稳定的瞬态变化流场,并且,由于模拟过程中,空间流场复杂漩涡结构的准确解析对该类型模拟非常关键,导致即便对二元机翼计算,所需的计算网格数量均相较RANS计算有数十倍甚至上百倍增加。也就是说,仅从计算网格规模考虑,如按二维RANS普遍采用数十万网格规模的网格开展计算来估算,采用DES或WMLES方法,计算网格普遍在数百万至数千万网格规模以上。

因此,为了克服上述情况带来的单机难以开展数百万至数千万网格规模非定常计算的困难,在计算过程中,引入基于大规模并行计算的方法非常关键。一方面,既解决无法计算的问题,另一方面又能充分发挥超级计算机的效能,提高全过程计算的效率。具体思路是,针对本文采用的非结构网格有限体积法,引入基于ParMetis的图分割算法实现计算域的自动划分,进而可非常方便地采用传统基于MPI的区域分解算法进行流场控制方程的并行求解[20]。

ParMetis是分区软件Metis的并行实现版本,采用基于图的多级分区算法自动、高效实现网格的自动剖分[21]。基于图的多级分区算法一般分3步完成(见图1):①图的粗化,通过网格单元聚团,自动划分形成大粒度区域;②初始化分区,对聚团后的大粒度区域进行初始分区;③多级优化分区,对分区结果进行逐层优化,形成最终的网格分区。

图1 网格分区:多级图划分算法[21]

采用上述多级图划分算法,不仅可以很好地保持子区域的连通性,而且由于优化阶段可以采用较好的算法,不仅可实现任意分区,并较好地做到各计算域的负载平衡,而且能够尽可能地有效减少子区域之间的边界单元数目,大幅度降低边界通信量,提高计算效率。

2 问题定义及计算网格生成

2.1 问题定义

选取的模型问题是,带双角冰冰型的GL305结冰翼型(即GL305/944)的后失速流场模拟。GL305翼型是典型的商务公务机翼型,在世界范围内得到了广泛的研究[10]。风洞试验研究表明,相对其他冰形,结冰时间在22.5 min的双角冰(编号为:944)对翼型气动性能的影响最为明显[22]。

选取的计算状态为:

M∞=0.12,Re=3.5×106,α=6°

(9)

对应了该结冰翼型在物理风洞中后失速附近的一个状态。由于涉及前缘冰形诱导引起的大范围分离流动,采用常规方法进行数值预测非常困难[10-11]。

图2给出了GL305/944结冰翼型的几何构型。其中,在上表面标记了9个典型位置以方便通过监测速度剖面数据监测空间流场。标记的位置分别位于:x/c=0.12,0.15,0.20,0.40,0.45,0.50,0.55,0.60,0.75。其中,c是翼型的弦长。

图2 GLC-305/944结冰翼型几何模型示意图

由于IDDES和WMLES模拟必须采用三维计算,在实际计算中,均根据二元风洞试验,假设三维沿展向方向形状保持不变。

2.2 计算网格生成

网格生成对RANS-LES混合方法的应用非常关键。针对本文的应用,为了构造一个特别适合的计算网格,首先遵循Spalart的DES网格生成基本原则[22],按照对流场的认知,对空间流场区域提前进行了规划,以在不同区域设计获得更合理的网格分布及间距。如图3所示,规划的主要空间区域包括Euler区(ER)、RANS区(RR)、焦点区(FR),分离区(DR)。在上述空间网格区域划分的规划中,ER区的特点是,尽管占据了大部分空间区域,但由于在该区域通常主要表现为无黏流特征,不涉及湍流或涡流计算,因而在该区域采用一般的网格分辨率,就能满足网格设计要求;RANS区主要是边界层附近相应区域,在这一区域,计算网格需要满足RANS计算对边界层网格的一般要求,尤其在壁面法向上需保证黏性子层到对数律层的网格分辨率;FR区、DR区均属于LES区,在这一区域需要设计满足LES计算需求的网格,充分保证网格分辨率,以更精细地捕捉空间涡流场结构。

图3 GLC-305/944结冰翼型流场区域的划分

按照上述规划原则,图4给出了实际采用计算网格在X-Y平面的示意图。在具体的计算网格生成中,采用了纯六面体网格生成策略,以保持其计算效率和计算准确性高的优势。同时,本文采用纯六面体网格的另一个考虑是,相对四面体网格,空间网格分布相对来说更加容易控制。

图4 计算网格示意图

具体计算网格生成中,计算网格的远场边界位于离机翼大约50倍弦长的位置,以充分规避远场的影响。核心区第一层网格到壁面的法向距离约为1.0×10-5m,并且为保证核心区内具有较充分网格密度,网格增长率取为1.1。由于采用三维非定常计算,计算网格展向宽度取为0.4c,对应的展向网格分布数为89。最终生成的整个计算域,总的计算网格单元数约为1 710万。这与文献[8]采用DHRL方法计算时给出的1 500万网格规模的密网格规模基本相当,也与文献[10]采用IDDES方法计算时给出的1 750万网格规模的基准网格规模基本密度基本相当。本文采用1 700万网格规模的主要考虑有:①充分保证网格生成中RANS区和LES区的网格密度;②基于文献网格规模经验,避免开展基于更大网格规模的网格无关性验证。实际上,文献[10]的研究表明,采用IDDES方法,基于1 700万网格规模给出的计算结果,尽管在流场结构捕捉精细程度上有细微差异,但总体与3 000万网格规模的计算结果基本相当。这为本文计算中网格规模的确定提供了参考依据。

计算中,边界条件为:远场采用Rieman无反射远场边界条件;物面采用绝热无滑移物面边界条件;对称面采用周期性边界条件。

3 计算结果及分析

正如前文多次提及的,WMLES和IDDES模拟必须采用三维非定常计算,以获得充分发展的三维湍流结构。为了加快计算,本文确定的基本计算测策略为:首先利用三维稳态RANS计算获得一个合适的初始场,然后利用RANS-LES混合方法(本文为WMLES和IDDES)计算大约5个周期的流场特征时间以建立较为充分发展的湍流非定常流场;之后,继续计算3个流场周期的特征时间,以便对流场数据进行统计获得平均场。流场特征时间T定义为翼型弦长c与自由来流速度U∞的比值,即T=c/U∞。在利用RANS-LES混合方法进行非定常计算时,如不做特殊说明,物理时间步均取为Δt=2×10-5s(对应无量纲物理时间步长ΔtU∞/c约为0.000 8),足够小的时间步长充分保证了计算域中CFL数小于1.0的收敛性条件成立。

值得指出,通过多少个流场特征时间来建立充分发展流场,并采用多少个流场特征时间来计算平均场,是RANS-LES混合方法这类非定常计算方法在实际应用中常涉及的一个基本问题。目前,鲜有文献针对这一问题展开讨论,也没有形成公认的普适性确定原则。理论上,建立流场的计算时间越长,流场发展越充分;计算平均场的计算时间越长,平均场相对误差越小。但相应地,计算代价越大。本文计算中,采用5个流场特征时间来建立流场,并采用3个流场特征时间来计算平均场,主要是从计算的经济性及对该问题本身频谱特性[10]的分析来综合确定的。

表1给出利用本文采用3种不同方法得到的升力系数结果与风洞试验结果[22]及文献[8]计算结果的对比。其中,文献[8]公布的DHRL方法是DDES方法的一种变种,优于传统DDES方法。

表1 总体气动特性比较

由表中的数据对比可以看出,RANS、IDDES和DHRL计算得到的升力系数均小于试验结果,而WMLES稍大于试验结果。相比之下,WMLES方法最接近试验值。

从该对比计算结果可以看到,对于本算例,RANS计算无法准确地预测升力系数,在所有结果中,相对误差最大,达到-26.7%。同时,还可以发现,通过本文IDDES和WMLES方法计算获得的升力系数非常接近,且均优于文献[8]方法提供的计算结果。

图5给出本文WMLES计算给出的翼型上方空间流场典型涡结构示意图,其中,涡结构采用速度梯度张量的二阶不变量表示,并用流向速度进行染色得到。从图中可以看出,本文计算得到了结冰翼型后失速流场的充分发展湍流结构,由前缘双角冰诱导产生的三维涡街结构清晰可见,并且,沿展向的三维涡旋结构的涡脱落情况也清晰地得到解析。

图5 典型瞬态流场的空间涡结构示意图

图6为不同方法得到的壁面平均压力系数对比。从结果来看,RANS几乎不能预测上壁面平均压力系数Cp,其中WMLES和IDDES都能给出与试验数据非常吻合的结果。相比之下,WMLES和IDDES均能较准确地预测上表面平顶区域在0.0～0.25之间。文献[8]中的DHRL可以给出比RANS更好的结果,但在这方面不如WMLES和IDDES。

图6 不同方法得到的物面压力系数分布对比

由图中的计算结果可知,在上表面流场区域,流场的加速(靠近角冰外缘)导致计算得到的吸力峰值远大于试验峰值,分区WMLES方法计算得到的上表面压力平顶与试验值基本重合,较好地描述了压力恢复过程,与试验值吻合最好;IDDES与分区WMLES方法计算结果总体一致,但压力恢复过程略为陡峭;RANS计算得到的上表面压力平顶明显低于试验值。在下表面弱分离流动区域,RANS、分区WMLES和IDDES结果趋势基本一致。这一结果进一步证实,RANS方法在该类型大范围分离流动模拟中一定程度存在模型能力不足的问题。

图7给出不同方法得到的平均流向速度场的对比。所有的模拟都给出了一个分离泡主导的空间流场。但是,以分离再附位置的预测作为判断依据,与试验相比,RANS过度预测了分离泡,而IDDES和WMLES则没有。本文中,IDDES和WMLES预测的分离线再附位置分别位于x/c=0.417 3和x/c=0.415 3处,与试验给出的x/c=0.53非常接近。

图7 不同方法得到的流向平均速度对比

图8给出不同方法的9个监测点的平均流向速度剖面数据与试验结果对比。从图中结果可以看出,在x/c=0.12位置,所有的计算都能观察到较强的回流,并且IDDES和WMLES都明显低估了速度峰值。除了前3个位置(x/c=0.12,0.15,0.20),RANS计算均无法得到与其他位置一致的速度分布,而IDDES和WMLES在各个站位都可以给出与试验数据较为一致的速度分布计算结果。相比之下,WMLES在几乎所有的站位均能给出与试验结果趋势一致、符合程度最优的速度分布预测结果。

图8 不同方法得到的各站位平均流向速度剖面数据对比

图9进一步给出了IDDES、WMLES与通过PIV试验获得的脉动速度均方根(RMS)对比。从图中可以看出,IDDES和WMLES都能准确预测流场核心区的速度脉动。在x/c=0.23时,通过IDDES预测了脉动速度的最大均方根为0.424 5,而在x/c=0.29时,通过WMLES预测了脉动速度的最大均方根为0.409 7,与试验中x/c=0.30处脉动速度的最大均方根0.34非常接近。相比之下,本文WMLES优于同一套计算网格IDDES给出的计算结果。

图9 不同方法得到的脉动速度均方根(RMS)云图对比

图10给出不同方法给出的所有9个监测点位置脉动速度均方根结果与试验数据的对比。从图中结果可以看出,在双角冰附近,即x/c=0.12处,所有计算结果均显示出强烈的湍流脉动。与试验数据相比,所有结果都高估了脉动速度的均方根峰值。相比之下,如果观察所有9个监控站位,WMLES可以给出更好的脉动速度均方根预测结果。在x/c=0.12～0.55的前7个位置,数值计算结果与试验结果符合良好。对于x/c=0.60和x/c=0.75的剩余位置,尽管对脉动速度预测不足,但计算结果整体趋势与试验一致。

图10 不同方法得到的各站位脉动速度均方根剖面数据对比

4 结 论

本文针对结冰翼型失速特性计算,提出发展了一种结合壁面模化大涡模拟(WMLES)和大规模并行计算的有效数值方法,并将其成功应用于GL305/944结冰翼型后失速流场计算,获得了满意的结果。给出了包括总体气动力、速度曲线、平均速度场,脉动速度均方根等在内的较为详细的计算结果,并与试验数据和文献[8]公布的DHRL方法计算结果(一种DDES变种方法)进行了综合比较。作为计算的对比验证,同时给出了RANS和IDDES计算的对比结果。研究得到以下结论:

1) 由于必须采用三维非定常计算,WMLES或DES计算中,引入大规模并行计算非常必要。

2) 针对双角冰结冰翼型后失速流场特性计算,与其他方法相比,WMLES方法非常有效。针对GL305/944结冰翼型,本文WMLES方法能够较为准确地预测总体气动力、压力平顶长度和压力恢复,以及角状冰引起的剪切层失稳,并能较好地预测分离和再附位置、速度脉动等关键参数;计算给出的升力系数相对误差仅为0.47%,远小于RANS和DHRL的-26.7%和-3.03%。相比之下,RANS方法对该类型涉及大范围分离的问题并不总是有效。

3) 与RANS和DHRL方法相比,IDDES内在综合了DDES和WMLES的能力,针对翼型结冰诱导的大范围分离流动,也不失为一种有效的方法。