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

年 月第 卷第 期西 北 工 业 大 学 学 报 ：收稿日期：基金项目：陕西省重点研发计划（）资助作者简介：李立（），西北工业大学博士研究生，主要从事计算流体力学仿真及并行计算技术研究。通信作者：武君胜（），西北工业大学教授，主要从事应用软件、计算方法、科学计算可视化及复杂领域软件工程研究。：基于大规模并行计算的结冰翼型失速流场特性数值模拟研究李立，武君胜，梁益华，田增冬西北工业大学 计算机学院，陕西 西安；中国航空工业集团公司西安航空计算技术研究所 预研部，陕西 西安 摘 要：结冰安全性评估是民用飞机适航的重要工作内容。翼面结冰将引起机翼前缘外形及边界层状态变化，并诱导大范围分离，进而导致飞行器升力面性能急剧降低，甚至带来严重的飞行安全问题。针对常规方法难以有效准确预测结冰翼型后失速流场空气动力学特性的问题，发展了一种结合大规模并行计算和壁面模化大涡模拟（）的有效数值计算方法，成功用于双角冰结冰翼型的后失速流场特性的数值模拟研究，取得了满意效果。数值模拟研究中，计算状态选取马赫数，雷诺数，攻角，对应了该翼型在风洞试验中后失速附近的流动状态。作为对比，同时给出雷诺平均（）方法及改进的时间延迟脱体涡模拟（）方法的计算结果，并与试验结果进行了综合比较。结果表明，是一种适于计算大范围分离流动的有效方法，针对结冰翼型后失速流场的数值预测，可大幅提高预测精度；针对文中的 结冰翼型，能相对准确地预测总体气动力、压力平顶长度和压力恢复，以及角状冰引起的剪切层失稳，且预测的升力系数相对误差仅为，远小于 方法的。关 键 词：结冰翼型；失速；壁面模化大涡模拟；并行计算；飞行安全中图分类号：文献标志码：文章编号：（）飞机结冰安全性是运输类飞机设计评估中必须关注的重要课题。美国联邦航空管理局（）和中国民用航空总局发布的适航条款中均有结冰安全性评估的相关条款和内容。在过去数年，随着计算机科学、计算数学等相关学科发展，计算流体力学（）技术得到长足发展，采用 方法开展飞机结冰安全性的预评估已越来越 普 遍。已 逐 步 成 为 评 估 飞 机 飞 行 性能，及 深 入 理 解 结 冰 流 动 机 理 的 强 有 力 工具。目前，业界形成的普遍共识是，尽管在实际工作过程中，精细 仿真存在非常耗时的问题，但与试验相比，投资回报巨大。而且，一般来说，开展全尺寸（或接近全尺寸）结冰试验的费用非常昂贵，并且在飞机设计完成之前几乎不可能开展。目前，可用于飞行器气动性能评估的 方法主要包括直接数值模拟（）、大涡模拟（），以及雷诺平均 方程（）方法等。其中，方法由于其实现了在计算资源及预测准确性上的有效平衡，在过去数十年间，得到广泛青睐和应用，对工程设计产生了巨大影响。然而，众所周知，以现代工程湍流模型为基石的 方法往往无法准确预测具有大范围分离的复杂湍流流场，存在模型部分失效的问题。和 方法在理论上可应用于上述情况，但由于需要巨大的计算成本，迄今还难以直接应用到复杂几何和全尺寸工程问题的计算。西 北 工 业 大 学 学 报第 卷为此，近年来人们致力于发展 混合方法，以弥补 与（甚至）在预测大范围分离复杂流场时的能力差距。混合方法是包括混合了不同 和 模型的一大类方法，主要针对复杂湍流问题，。混合方法的独到之处是，将计算代价降低到非常接近时域非定常（即）计算的程度，但仍可以提供比 更准确的方法来解析关键流场区域的高度分离湍流流场结构。目前，混合方法正成为在 模型失效情况下、替代用于复杂工程设计的一种有前景方法。综合当前技术发展来看，构造 混合方法的方式主要包括 种：从 出发、自下而上的方法，属于早期普遍采用的方法。典型代表是基于（）方程湍流模型发展的脱体涡模拟（）方法，自提出以来即得到广泛跟踪和研究；从 出发、自上而下的方法，典型代表是壁面模化大涡模拟（）方法，这是近年逐步引起重视的一大类方法。在美国 和波音公司发布的 发展愿景中，方法被认为是“面 向 复 杂 整 机 应 用 的 下 一 代 方法”。然而，混合方法作为一种精细 方法，与常规 计算相比，不单必须采用长时非定常计算（以建立稳定的瞬态变化流场），而且，由于需要解析空间流场的复杂漩涡结构（特别是三维横向涡结构），所需的计算网格数量往往需要增加数十倍甚至上百倍。这意味着计算过程中，引入大规模并行计算以加速计算过程非常关键。由于结冰翼型失速流场预测的复杂性，目前针对该问题，采用 混合方法进行精细分析的工作并不多见。近年来比较具有代表性的工作包括 等提出的一种变种的 方法（即）、等基于自主程序发展的以 为基底模式的 方法等。这些方法均是采用自底向上的传统方法进行 方法构造。本文针对该问题，提出发展结合大规模并行计算和壁面模化大涡模拟（）的有效计算策略和方法。通过针对双角冰结冰翼型 的后失速流场特性数值模拟研究，并与以两方程 模型为基底模式的 计算及 计算开展细致的结果对比，证实了本文发展方法的有效性，可为基于数值计算的飞机结冰安全性评估提供有效技术手段。数值方法 流场控制方程在三维笛卡尔坐标系下，方程写成守恒形式如下 （）（，）（，（）（，）式中：，（，），分别表示密度、速度、总内能、压强、剪切应力张量和热流通量。剪切应力张量按照 假设，可表示为 （）（）式中：为黏性系数；，为湍动能。压强可根据状态方程计算，对理想气体有 （）（）（）式中：为比热容比，对空气，。对方程（），本文采用非结构混合网格策略和有限体积法进行数值求解。目前，利用非结构网格进行有限体积离散主要有 种方式：直接采用网格单元作为控制体，称为格心（）格式；以网格单元顶点为中心，采用对偶网格方法建立控制体，称为格点（）格式。本文选用格点格式。值得指出，对 混合方法来说，由于涉及流场中小涡结构的捕捉，数值格式的耗散控制是必须解决的一个重要问题。为此，在实际数值求解中，对空间离散，本文对流项引入一种全速域自适应低耗散 格式进行计算，具有对边界层及全流场计算自动保持低耗散及高性能的特点；黏性项则采用标准中心差分格式进行计算。对时间迭代，为了获得瞬态流场解，采用改进的 阶精度后向差分（，）格式及双时间步方法进行非定常计算。第 期李立，等：基于大规模并行计算的结冰翼型失速流场特性数值模拟研究 湍流模型为了实现湍流模拟，除了，同时采用以两方程 模型为基底模式的 和 方法进行对比计算。本文 模型的具体定义公式为（），（）（）式中：为卡门常数；为 亚格子模型常数，而 （）（）为普朗特混合长度模型系数，为计算网格的长度尺度，简单定义为 （，）（）依据模型（）式，当 ，即在在黏性固壁附近，会使用与代数 模型相当的普朗特混合长度假设；而当 时，表示远离黏性固壁，则采用经典 亚格子 模型。由此不难分析，上述 模型本质上可以看做是经典 亚格子 模型和普朗特混合长度理论的一种混合模型。采用上述方法，事实上建立了一种避免在黏性子层物理解析壁湍流的方法，从而建立了一类典型的、自顶向下的 混合方法。以两方程 模型为基底模式的 模型完整计算公式可在文献中找到。作为一种典型的脱体涡（）模型，其主要设计思想是将 模型中湍动能方程的 长度尺度替换为 混合长度尺度，从而使该模型表现为 和 方法混合的一种行为。修正后的湍动能方程为（）（）（）（）式中，混合长度尺度，定义为（）（），（）（）是一个范围从 到 的混合函数，用来实现模型在（）和（）之间的自动切换。基于大规模并行的加速计算需要指出，无论是从 出发的 类方法，还是从 出发的 方法，本质上都是一种基于时间准确性（瞬态）计算的方法。因而，当使用 或 方法进行 仿真时，不仅必须采用长时非定常计算，以建立稳定的瞬态变化流场，并且，由于模拟过程中，空间流场复杂漩涡结构的准确解析对该类型模拟非常关键，导致即便对二元机翼计算，所需的计算网格数量均相较 计算有数十倍甚至上百倍增加。也就是说，仅从计算网格规模考虑，如按二维 普遍采用数十万网格规模的网格开展计算来估算，采用 或 方法，计算网格普遍在数百万至数千万网格规模以上。因此，为了克服上述情况带来的单机难以开展数百万至数千万网格规模非定常计算的困难，在计算过程中，引入基于大规模并行计算的方法非常关键。一方面，既解决无法计算的问题，另一方面又能充分发挥超级计算机的效能，提高全过程计算的效率。具体思路是，针对本文采用的非结构网格有限体积法，引入基于 的图分割算法实现计算域的自动划分，进而可非常方便地采用传统基于 的区域分解算法进行流场控制方程的并行求解。是分区软件 的并行实现版本，采用基于图的多级分区算法自动、高效实现网格的自动剖分。基于图的多级分区算法一般分 步完成（见图）：图的粗化，通过网格单元聚团，自动划分形成大粒度区域；初始化分区，对聚团后的大粒度区域进行初始分区；多级优化分区，对分区结果进行逐层优化，形成最终的网格分区。图 网格分区：多级图划分算法采用上述多级图划分算法，不仅可以很好地保持子区域的连通性，而且由于优化阶段可以采用较好的算法，不仅可实现任意分区，并较好地做到各计算域的负载平衡，而且能够尽可能地有效减少子区域之间的边界单元数目，大幅度降低边界通信量，提西 北 工 业 大 学 学 报第 卷高计算效率。问题定义及计算网格生成 问题定义选取的模型问题是，带双角冰冰型的 结冰翼型（即）的后失速流场模拟。翼型是典型的商务公务机翼型，在世界范围内得到了广泛的研究。风洞试验研究表明，相对其他冰形，结冰时间在 的双角冰（编号为：）对翼型气动性能的影响最为明显。选取的计算状态为：，（）对应了该结冰翼型在物理风洞中后失速附近的一个状态。由于涉及前缘冰形诱导引起的大范围分离流动，采用常规方法进行数值预测非常困难。图 给出了 结冰翼型的几何构型。其中，在上表面标记了 个典型位置以方便通过监测速度剖面数据监测空间流场。标记的位置分别位于：，。其中，是翼型的弦长。由于 和 模拟必须采用三维计算，在实际计算中，均根据二元风洞试验，假设三维沿展向方向形状保持不变。图 结冰翼型几何模型示意图 计算网格生成网格生成对 混合方法的应用非常关键。针对本文的应用，为了构造一个特别适合的计算网格，首先遵循 的 网格生成基本原则，按照对流场的认知，对空间流场区域提前进行了规划，以在不同区域设计获得更合理的网格分布及间距。如图 所示，规划的主要空间区域包括 区（）、区（）、焦点区（），分离区（）。在上述空间网格区域划分的规划中，区的特点是，尽管占据了大部分空间区域，但由于在该区域通常主要表现为无黏流特征，不涉及湍流或涡流计算，因而在该区域采用一般的网格分辨率，就能满足网格设计要求；区主要是边界层附近相应区域，在这一区域，计算网格需要满足 计算对边界层网格的一般要求，尤其在壁面法向上需保证黏性子层到对数律层的网格分辨率；区、区均属于 区，在这一区域需要设计满足 计算需求的网格，充分保证网格分辨率，以更精细地捕捉空间涡流场结构。图 结冰翼型流场区域的划分按照上述规划原则，图 给出了实际采用计算网格在 平面的示意图。在具体的计算网格生成中，采用了纯六面体网格生成策略，以保持其计算效率和计算准确性高的优势。同时，本文采用纯六面体网格的另一个考虑是，相对四面体网格，空间网格分布相对来说更加容易控制。具体计算网格生成中，计算网格的远场边界位于离机翼大约 倍弦长的位置，以充分规避远场的影响。核心区第一层网格到壁面的法向距离约为，并且为保证核心区内具有较充分网格密度，网格增长率取为。由于采用三维非定常计算，计算网格展向宽度取为，对应的展向网格分布数为。最终生成的整个计算域，总的计算网格单元数约为 万。这与文献采用 方法计算时给出的 万网格规模的密网格规模基本相当，也与文献采用 方法计算时给出的 万网格规模的基准网格规模基本密度基本相当。本文采用 万网格规模的主要考虑有：充分保证网格生成中 区和 区的网格密度；基于文献网格规模经验，避免开展基于更大网格规模的网格无关性验证。实际上，文献的研究表明，采用 方法，基于 万网格规模给出的计算结果，尽管在流场结构捕捉精细程度上有细微差异，但总体与 万网格规模的计算结果基本相当。这为本文计算中网格规模的确定提供了参考依据。计算中，边界条件为：远场采用 无反射第 期李立，等：基于大规模并行计算的结冰翼型失速流场特性数值模拟研究远场边界条件；物面采用绝热无滑移物面边界条件；对称面采用周期性边界条件。图 计算网格示意图 计算结果及分析正如前文多次提及的，和 模拟必须采用三维非定常计算，以获得充分发展的三维湍流结构。为了加快计算，本文确定的基本计算测策略为：首先利用三维稳态 计算获得一个合适的初始场，然后利用 混合方法（本文为 和）计算大约 个周期的流场特征时间以建立较为充分发展的湍流非定常流场；之后，继续计算 个流场周期的特征时间，以便对流场数据进行统计获得平均场。流场特征时间 定义为翼型弦长 与自由来流速度 的比值，即 。在利用 混合方法进行非定常计算时，如不做特殊说明，物理时间步均取为（对应无量纲物理时间步长 约为），足够小的时间步长充分保证了计算域中 数小于 的收敛性条件成立。值得指出，通过多少个流场特征时间来建立充分发展流场，并采用多少个流场特征时间来计算平均场，是 混合方法这类非定常计算方法在实际应用中常涉及的一个基本问题。目前，鲜有文献针对这一问题展开讨论，也没有形成公认的普适性确定原则。理论上，建立流场的计算时间越长，流场发展越充分；计算平均场的计算时间越长，平均场相对误差越小。但相应地，计算代价越大。本文计算中，采用 个流场特征时间来建立流场，并采用 个流场特征时间来计算平均场，主要是从计算的经济性及对该问题本身频谱特性的分析来综合确定的。表 给出利用本文采用 种不同方法得到的升力系数结果与风洞试验结果及文献计算结果的对比。其中，文献 公布的 方法是 方法的一种变种，优于传统 方法。表 总体气动特性比较方法升力系数相对误差 风洞试验 由表中的数据对比可以看出，、和 计算得到的升力系数均小于试验结果，而 稍大于试验结果。相比之下，方法最接近试验值。从该对比计算结果可以看到，对于本算例，计算无法准确地预测升力系数，在所有结果中，相对误差最大，达到。同时，还可以发现，通过本文 和 方法计算获得的升力系数非常接近，且均优于文献方法提供的计算结果。图 给出本文 计算给出的翼型上方空间流场典型涡结构示意图，其中，涡结构采用速度梯度张量的二阶不变量表示，并用流向速度进行染色得到。从图中可以看出，本文计算得到了结冰翼型后失速流场的充分发展湍流结构，由前缘双角冰诱导产生的三维涡街结构清晰可见，并且，沿展向的三维涡旋结构的涡脱落情况也清晰地得到解析。图 典型瞬态流场的空间涡结构示意图图 为不同方法得到的壁面平均压力系数对比。从结果来看，几乎不能预测上壁面平均压力系数，其中 和 都能给出与试验数据非常吻合的结果。相比之下，和 均能较准确地预测上表面平顶区域在 西 北 工 业 大 学 学 报第 卷之间。文献中的 可以给出比 更好的结果，但在这方面不如 和。图 不同方法得到的物面压力系数分布对比由图中的计算结果可知，在上表面流场区域，流场的加速（靠近角冰外缘）导致计算得到的吸力峰值远大于试验峰值，分区 方法计算得到的上表面压力平顶与试验值基本重合，较好地描述了压力恢复过程，与试验值吻合最好；与分区 方法计算结果总体一致，但压力恢复过程略为陡峭；计算得到的上表面压力平顶明显低于试验值。在下表面弱分离流动区域，、分区 和 结果趋势基本一致。这一结果进一步证实，方法在该类型大范围分离流动模拟中一定程度存在模型能力不足的问题。图 给出不同方法得到的平均流向速度场的对比。所有的模拟都给出了一个分离泡主导的空间流场。但是，以分离再附位置的预测作为判断依据，与试验相比，过度预测了分离泡，而 和 则没有。本文中，和 预测的分离线再附位置分别位于 和 处，与试验给出的 非常接近。图 不同方法得到的流向平均速度对比 图 给出不同方法的 个监测点的平均流向速度剖面数据与试验结果对比。从图中结果可以看出，在 位置，所有的计算都能观察到较强的回流，并且 和 都明显低估了速度峰值。除了前 个位置（，），计算均无法得到与其他位置一致的速度分布，而 和 在各个站位都可以给出与试验数据较为一致的速度分布计算结果。相比之下，在几乎所有的站位均能给出与试验结果趋势一致、符合程度最优的速度分布预测结果。图 进一步给出了、与通过 试验获得的脉动速度均方根（）对比。从图中可以看出，和 都能准确预测流场核心区的速度脉动。在 时，通过 预测了脉动速度的最大均方根为 ，而在 时，通过 预测了脉动速度的最大均方根为 ，与试验中 处脉动速度的最大均方根 非常接近。相比之下，本文 优于同一套计算网格 给出的计算结果。第 期李立，等：基于大规模并行计算的结冰翼型失速流场特性数值模拟研究图 不同方法得到的各站位平均流向速度剖面数据对比图 不同方法得到的脉动速度均方根（）云图对比图 给出不同方法给出的所有 个监测点位置脉动速度均方根结果与试验数据的对比。从图中结果可以看出，在双角冰附近，即 处，所有计算结果均显示出强烈的湍流脉动。与试验数据相比，所有结果都高估了脉动速度的均方根峰值。相比之下，如果观察所有 个监控站位，可以给出更好的脉动速度均方根预测结果。在 的前 个位置，数值计算结果与试验结果符合良好。对于 和 的剩余位置，尽管对脉动速度预测不足，但计算结果整体趋势与试验一致。西 北 工 业 大 学 学 报第 卷图 不同方法得到的各站位脉动速度均方根剖面数据对比 结 论本文针对结冰翼型失速特性计算，提出发展了一种结合壁面模化大涡模拟（）和大规模并行计算的有效数值方法，并将其成功应用于 结冰翼型后失速流场计算，获得了满意的结果。给出了包括总体气动力、速度曲线、平均速度场，脉动速度均方根等在内的较为详细的计算结果，并与试验数据和文献公布的 方法计算结果（一种 变种方法）进行了综合比较。作为计算的对比验证，同时给出了 和 计算的对比结果。研究得到以下结论：）由于必须采用三维非定常计算，或 计算中，引入大规模并行计算非常必要。）针对双角冰结冰翼型后失速流场特性计算，与其他方法相比，方法非常有效。针对 结冰翼型，本文 方法能够较为准确地预测总体气动力、压力平顶长度和压力恢复，以及角状冰引起的剪切层失稳，并能较好地预测分离和再附位置、速度脉动等关键参数；计算给出的升力系数相对误差仅为 ，远小于 和 的和。相比之下，方法对该类型涉及大范围分离的问题并不总是有效。）与 和 方法相比，内在综合了 和 的能力，针对翼型结冰诱导的大范围分离流动，也不失为一种有效的方法。第 期李立，等：基于大规模并行计算的结冰翼型失速流场特性数值模拟研究参考文献：，：，中国民用航空局 中国民用航空规章第 部：运输类飞机适航标准 ，：，魏扬，徐浩军，薛源，等 机翼前缘结冰对大飞机操稳特性的影响 北京航空航天大学学报，（）：，（）：（），：，：，：，：，：陈迎春，张美红，张淼，等 大型客机气动设计综述 航空学报，（）：，（）：（），：，：阎超，于剑，徐晶磊，等 模拟方法的发展成就与展望 力学进展，（）：，（）：（）李立，麻蓉，梁益华 一类全速域低耗散 格式的构造及性能分析第 届全国计算流体力学会议，厦门，（），：，赵钟，张来平，何磊，等 适用于任意网格的大规模并行 计算框架 计算机学报，（）：西 北 工 业 大 学 学 报第 卷，：，（）：（），（）：，：，；，：，（），（），；，：；（）；引用格式：李立，武君胜，梁益华，等 基于大规模并行计算的结冰翼型失速流场特性数值模拟研究 西北工业大学学报，（）：，（）：（）（：），