基于高精度格式的翼型失速特性脱体涡数值模拟

2023-08-28钟伯文

科技创新与应用 2023年24期

王瑞，钟伯文*

（1.南昌航空大学飞行器工程学院，南昌 330063；2.江西省飞行器设计与气动仿真重点实验室，南昌 330063）

机翼在大攻角状态下的气动特性因湍流分离流动的存在呈现高度的非定常特征。当攻角超过某临界值时，机翼的升力会急剧减小，阻力会大幅度增大，导致飞机处于失速状态而引发飞行事故。因此，精确评估飞机的失速特性对飞机的设计和飞机的性能优化具有重要意义。

通常，机翼翼型在小攻角状态的流动接近定常状态，可以使用传统的雷诺平均方程（Reynolds-averaged Navier Stokes Equation，RANS）进行较为精确的数值模拟。然而，当攻角一旦增大至接近甚至超过某临界值，翼型外流场会出现大幅度的分离，传统的RANS 方法无法精确模拟这类湍流的大分离非定常流动。

除了RANS 方法，还有直接数值模拟（Direct Numerical Simualtion，DNS）和大涡模拟（Large Eddy Simulation，LES）[1]等湍流模拟方法。从理论上讲，DNS方法能够准确模拟这类分离流动，但由于其对网格量的要求很高，计算耗费巨大，无法满足飞机工程型号的应用需求。LES 方法采用空间过滤的方式过滤湍流的小尺度脉动，直接模拟大尺度涡结构，而小尺度湍流脉动的影响通过亚格子模型来模拟，LES 方法相对于DNS可以大幅度降低计算量，但在近壁面仍然对网格密度有很高的要求[2]，使得LES 难于获得工程型号的实际应用。波音公司的Spalart 等[3]根据LES 和RANS 两种方法的特点和优点，提出了RANS/LES 混合方法，即脱体涡模拟（Detached Eddy Simulation ，DES），在近壁区采用RANS 方法，在远离壁面区域采用LES 方法，使得其方法在飞机工程型号的应用成为可能。其后，还形成了不同的改进方法，如延迟脱体涡模拟方法[4]（Delayed-DES，DDES）、DDES 增强版本[5]（Improved-DDES，IDDES）及分区域脱体涡模拟方法[6]（Zonal-DES，ZDES）等。

目前DES 类方法已经对翼型大分离流动问题进行了广泛研究和应用[7-8]。Shur 等[9]将DES 方法应用于NACA0012 翼型低速状态的失速特性预测，并且与URANS 进行了对比，得到了不错的模拟效果。本文作者Zhong 等[10]将DES 方法与URANS、隐式大涡模拟（Implicit LES，ILES）进行了对比，体现了DES 方法对预测NACA0012 翼型失速特性的较强能力。Li[11]将BL 模型（Baldwin -Lomax）与SA 模型（Spalart-Allmaras）结合构造了基于SA+BL 模型的DES 方法，并对NACA 系列翼型进行了失速特性的数值模拟。

分离流动的精确模拟不仅受到湍流模拟方法的影响，同时也与数值格式密切相关。Ghosal[12]和Chow 等[13]评估了数值误差对LES 方法的影响，结果都发现数值格式对LES 计算结果的影响占比甚至超过了亚格子模型。Kawai 等[14]将基于B-L 模型（Baldwin-Lomax）的RANS/LES 混合方法与高阶紧致差分格式（Compact Differencing Scheme，CMPT）结合模拟了NACA64A-006 薄翼型失速特性，模拟结果表明了相同网格密度下高阶CMPT与传统迎风格式相比可以得到更可靠的计算结果。

综上所述，翼型的失速现象所表现的分离流动十分复杂，精确高效的模拟对湍流模拟方法和数值格式具有很高的要求。尽管DES 方法已经在翼型失速特性的预测中得到了应用，但是在数值格式、方法的适用性及模拟精度等方面，仍然存在进一步的提升空间，需要深入地探索研究。

本文采用脱体涡模拟（DES）方法，并将HLLC 低耗散近似黎曼求解器与五阶WENO 格式结合对NACA0015 翼型的失速特性和相应的流场进行了数值模拟研究，并与二阶Jameson 中心差分格式（Centered Differencing Scheme，CDS）[15]的数值模拟结果分别从失速特性、压力系数分布及非定常湍流结构3 个方面进行对比分析。

1 数值模拟方法

1.1 HLLC 格式

Toro 等[16]为解决HLL 格式在间断面会产生较大的数值耗散现象提出了HLLC 格式。HLLC 格式本身是一阶精度并且构造简单，具有低耗散、高分辨率的特性，与Roe 格式相比不需要进行熵的修正来保持物理数值解，适合跨音速和低速计算。

HLLC 通量和相关的黎曼求解器都包含了黎曼问题的所有波，不需要将方程线性化，而且不通过固定熵值的方法来解决低密度问题。图1 给出了黎曼问题示意图，Toro 等[16]改进的两态HLLC 格式可以准确地解决孤立激波、间断波和剪切波，并且在黎曼扇形内部结构中保留了初始的正密度和压力。

图1 黎曼问题示意图

HLLC 格式求解的格式如下

1.2 WENO 格式

WENO 格式（Weighted Essentially Non-oscillatory）是Liu 等[17]在ENO 算法（Essentially Non-oscillatory）的基础上发展的。近年来WENO 算法在复杂流动的计算得到了广泛应用[18]。Jiang 等[19]建议采用IS 光滑因子度量函数的光滑性，构成了WENO-JS 格式。

五阶WENO-JS 格式的三阶数值通量在3 个模板中可以表示为

光滑因子表示为

非线性权重因子定义如下

式中：k 的取值范围为[0，1，2]；ε 为防止分母为0 的充分小正数，本文取10-6；P 为用于控制ENO 的行为，本文取4。

最终的重构函数形式如下

1.3 Jameson 中心差分格式

Jameson 中心差分格式用于空间离散，具体形式如下

该格式允许杂波振荡，为了解决这个问题，使用耗散通量模拟了人工耗散，得到的最终Jameson 中心差分格式如下

1.4 时间推进格式

以偏微分方程y′=f（x，y），y（x，y）=y0为例，在点x=x0+h 处，y 和h 为矢量。五阶Runge-Kutta 格式的具体形式如下

式中：ν 为Runge-Kutta 阶数，本文取5。

2 数值模拟分析

本文采用NACA0015 翼型作为研究算例，展长设置为0.2c，c 为翼型弦长。计算网格采用多块结构化网格，网格块划分如图2 所示，流场计算域为15c，网格点数分布为：流向点数647 个、法向点数115 个、展向点数81 个，网格单元总数为6×106。实验数据来源LM Glasfiber 风洞[20]，计算的模拟流动状态为：自由来流马赫数为0.078，基于弦长的雷诺数为1.6×106。计算程序中，为了加速收敛和减少计算耗费，非定常时间推进采用双时间步隐式迭代法。时间步长采用无量纲形式，时间步长Δt=0.01，对应的物理时间步长为3.5×10-5s。

图2 多块网格划分

2.1 失速特性预测分析

本文对比了HLLC+五阶WENO 格式与CDS 的计算结果。图3 和图4 分别给出了不同数值格式得到的升力系数曲线（CL）和阻力系数（CD）曲线。图3 结果表明，2 种数值格式都能准确预测失速的临界攻角及小攻角的升力系数。但是，在近失速段（攻角在8.13～13°范围内），CDS 预测得到的近失速升力系数与实验值相比偏小，13°攻角的升力系数预测偏差最大，而HLLC+五阶WENO 格式得到的升力系数误差较小，增长趋势和实验值相同。在过失速段（攻角在14.3～19°范围内），HLLC+五阶WENO 格式能够准确预测升力系数，而CDS 的预测的整体数值偏小。

图3 不同数值格式的升力系数曲线

图4 不同数值格式的阻力系数曲线

根据图4 阻力系数预测结果，当攻角小于17°，2种格式得到的数值结果基本一致，在临界攻角附近与实验值相比偏大。当攻角为17°时，CDS 预测得到的阻力系数与实验相比偏大，HLLC+五阶WENO 格式的结果更接近实验值，而在19°攻角状态下，2 种差分格式的预测结果都偏大。

2.2 压力系数分布预测

图5 给出了在15.66°攻角下时间平均后的翼型表面压力系数（CP）分布模拟结果。结果表明，在翼型上表面，两种格式得到的上表面压力分布系数转折点基本一致，并且与实验值相同，因此2 种格式都准确地预测了流动的分离点位置。CDS 预测翼型表面压力系数分布的整体数值偏小，这导致了该攻角下升力系数的数值结果小于实验值。HLLC+五阶WENO 格式预测得到的压力峰值和分离区的压力系数偏大，而其他位置的数值与实验值比较接近，因此在该攻角下预测的升力系数比实验值略大。

图5 AOA=15.66°时间平均的翼型表面压力系数分布

表1 给出了15.66°攻角的平均升力系数对比情况。根据升力系数相对误差的对比结果，HLLC+五阶WENO 格式得到的升力系数预测的相对误差较小，因此该格式得到的大攻角压力系数分布的数值模拟结果更接近实验结果。

表1 AOA=15.66°平均升力系数对比

2.3 大攻角湍流结构分析

图6 给出了在19°攻角状态下的瞬时湍流结构等值面图（Q 准则）。结果显示，2 种格式都能捕捉到大尺度湍流涡结构。HLLC+五阶阶WENO 格式的耗散性相对较低、模拟精度较高，能够捕捉到小尺度湍流涡结构，能够更为准确地模拟大分离湍流运动，从而更准确地预测翼型的失速特性。而CDS 的模拟精度低、耗散性大，无法捕捉较小尺度的湍流涡结构。