基于数据驱动转捩模型的翼型动态失速气动力计算

李金瑛, 戴玉婷,2, 杨 超

(1. 北京航空航天大学航空科学与工程学院, 北京 100083; 2. 天目山实验室, 浙江杭州 310023)

引 言

动态失速是机翼做大幅度俯仰运动时产生的非线性、 非定常气动现象。低Reynolds数下, 飞行器大幅俯仰运动引起分离, 诱导转捩, 并可能进一步诱发失速颤振等气动弹性失稳现象[1-2]。国内外对动态失速开展了大量实验与仿真研究[3-5]。Spentzos等[6]采用非定常Reynolds平均N-S(unsteady Reynolds-averaged Navier-Stokes, URANS)方程及k-ωSST二方程湍流模型对三维机翼进行了数值仿真计算。Wang等[7]尝试了将URANS与分离涡(detached eddy simulation, DES)方法结合用于动态失速仿真。Kim等[8]采用了大涡模拟(large-eddy simulation, LES)方法对NACA0012翼型进行了动态失速仿真。

传统湍流模型包括Spalart-Allmaras一方程湍流模型和k-ε、k-ω、k-ωSST二方程湍流模型等。其中, Menter[9]提出的k-ωSST模型具有良好的综合性能, 但无法对转捩过程进行准确计算。现有转捩过程计算方法主要有3种: 基于线性稳定性理论的eN方法[10-12], 采用低Reynolds数湍流模型的仿真方法[13], 及采用引入转捩间歇因子(γ)的湍流模型。其中, eN方法仅适用于自然转捩的计算; 低Reynolds数湍流模型忽略了转捩区物理特性, 易造成误差[14]。采用转捩间歇因子控制转捩发展是目前主流的计算方式。国内外发展了多种γ的计算方法及湍流模型[15-18]: Dhawan等[19]提出了γ的代数计算方法。Suzen等[20]结合前人工作[21-22]提出了耦合RANS的γ输运方程。Menter等[14]提出了γ-Reθ四方程湍流模型, 并将其简化为SST-γ三方程湍流模型[23], 该模型在RANS方程中对转捩过程模拟精度较高, 但增加了一个输运微分方程, 计算效率低于k-ωSST二方程模型。目前, 须发展同时具备高精度与高效率的转捩湍流模型。

利用数据驱动建模可以提高湍流模拟的效率。当前数据驱动建模方法包括直接辨识气动力系数的黑箱模型、 耦合求解湍流方程的灰箱模型及辨识优化湍流方程参数的白箱模型。在白箱模型方面, Singh等[24-25]使用机器学习修正了Spalart-Allmaras湍流模型参数; Yang等[26]采用随机森林与人工神经网络对k-ω-γ-Ar四方程湍流模型参数进行了修正; Zafar等[27]基于递推神经网络与卷积神经网络建立了具有高泛化能力的转捩湍流模型。此外, 国内外学者通过将神经网络与CFD湍流方程耦合求解[28-30]建立了多种灰箱模型。

本文基于SST-γ三方程湍流模型的流场计算结果, 训练深度神经网络模型预测转捩间歇因子, 建立流场参数与γ的关系, 据此修正k-ωSST二方程湍流模型, 建立将深度神经网络嵌入CFD工具的SST-γ-machine learning(下文简称为SST-γ-ML)耦合二方程湍流模型。利用该耦合二方程湍流模型对NACA0012翼型进行动态失速流场仿真。

1 耦合二方程湍流模型

1.1 耦合湍流模型建模

本节基于k-ωSST与数据驱动的转捩间歇因子模型建立SST-γ-ML耦合二方程湍流模型, 在不增加输运微分方程个数的前提下, 使耦合二方程湍流模型具备模拟转捩过程的能力。

SST-γ-ML耦合二方程湍流模型的控制方程为

其中,S为Reynolds应力,Ω为涡量幅值,μ为动力黏度,μt为湍流黏度。其他常数与函数为

其中,ReT为湍流Reynolds数

(1)

控制方程组由k与ω的输运微分方程及γ的神经网络代数方程组成。此方程组中, 前两个输运方程即构成由γ修正的k-ωSST二方程湍流模型。γ在SST-γ三方程模型中由输运微分方程求出, 说明其与流场中U、 压力p、k、ω等流场参数具有复杂的非线性关系, 因此, 采用神经网络代替输运微分方程预测γ, 可提高计算效率, 该神经网络具体结构在1.2节详细介绍。

SST-γ-ML耦合二方程湍流模型的训练与嵌入策略如图1所示。红色虚线框内表示了训练过程, 选取对转捩过程呈现相对准确的SST-γ三方程模型计算的流场参数和转捩间歇因子作为神经网络模型的输入输出数据集。

图1 模型的训练与使用Fig. 1 Training and use of model

神经网络与CFD平台的耦合方式有松耦合与紧耦合。松耦合在平台间传递流场数据, 紧耦合将神经网络模型直接嵌入CFD平台。对非定常流动, 如动态失速过程, 每个时间步间流场变化大, 须在每个时间步对CFD仿真计算与神经网络模型进行交互, 松耦合并不适用。本文采用紧耦合方式建立数据驱动的转捩模型, 如图1绿色线框所示。

基于SST-γ-ML耦合二方程湍流模型建模主要分5个步骤, 包括训练神经网络、 冻结导出、 嵌入耦合模型、 模型编译与CFD仿真计算。为将训练完成的神经网络完整嵌入OpenFOAM开源CFD软件中, 使用protocol buffer(pb)格式文件导出神经网络模型。

在CFD计算中对每个流体网格单元均独立调用神经网络模型, 预测每个流体网格单元的转捩间歇因子。虽然该建模方式计算效率低于直接预测全流场转捩间歇因子, 但大幅简化了神经网络结构, 训练集数据量降低2～3个数量级。此外, 它对流场不同网格单元具有泛化性, 能够直观反映出流场参数与不同位置γ之间的关系。

1.2 神经网络模型和输入参数设计

本文基于反向传播的全连接结构人工神经网络建立流场参数与转捩间歇因子之间的预测模型。反向传播(back-propagation, BP)神经网络的参数及样本数据预处理方法具体如下:

采用Adam算法作为误差反向传播的优化算法, 均方误差(mean squared error, MSE)作为误差计算算法; 神经网络隐藏层层数为4层, 每层分别含有[40, 20, 10, 10]个神经元, 批处理数为256, 学习率为0.001, 衰减率为0.001; 在输入层到第1层隐藏层间加入批归一化层, 从第4层隐藏层到输出层选择Sigmoid作为激活函数, 以控制转捩间歇因子输出范围为0～1。隐藏层间激活函数选取Tanh以扩大近壁面数据梯度; 训练前对样本数据进行打乱预处理, 并采用多个时间步流场参数作为训练集数据。

通过对流场参数与转捩间歇因子关系的物理分析及优化设计, 选取神经网络模型的流场输入参数组合为

(2)

其中,Ux为平行来流方向速度,Uy为垂直来流方向速度,dω为网格单元与壁面的距离。

6个输入参数中,p,μt/ρ,Ux,Uy直接从流场中提取。k/ω与湍流Reynolds数ReT呈正相关(式(1)),k表征湍流运动分量具有的动能

流动从层流转捩到湍流的过程中, 湍流动能逐渐增大, 与壁面间的剪切应力逐渐减小, 对应k增大,ω减小, 即输入参数k/ω增大。

其中, dUy/dy为垂直来流方向速度导数。比例系数λθL表征流场压力梯度, 壁面距离平方项与比例系数呈线性关系, 说明转捩程度对壁面距离非常敏感。在接近壁面的位置, 极小距离下会发生显著的转捩发展, 故采用平方项输入参数。

2 数据驱动的转捩模型验证

以标准T3A平板为研究对象, 在Re=6.12×105下进行稳态转捩仿真计算, 验证数据驱动的耦合二方程湍流模型的有效性。

T3A平板长度为1.7 m, 流场高度1 m, 平板前缘半径为0.75 mm, 流场主要区域为平板前区域(0～0.04 m)、 前缘区域与平板区域(0.04～1.70 m), 平板流场网格如图2所示。在平板前缘附近进行网格加密, 网格量为26 820, 网格平均Y+为0.53。平板表面(wall)采用无滑移边界, 流场顶部及平板前区域的流场底部(above)采用滑移边界, 远场设置自由来流速度U∞、 湍流度Tu∞与黏度比Rμ=μt/μ, 出口压力为0。

图2 T3A平板计算网格及边界条件Fig. 2 Mesh and boundary conditions of T3A flat plate

采用如表1所示入口边界条件, 分别采用SST-γ-ML耦合二方程湍流模型,k-ωSST二方程湍流模型与SST-γ三方程湍流模型计算T3A标准平板表面湍流度与摩阻系数, 结果如图3所示。

表1 T3A平板入口条件

(a) Turbulence intensity

结果表明, 耦合二方程湍流模型对T3A平板定常流动条件下的转捩位置及对应摩阻系数预测准确, 与实验值[31]及三方程湍流模型结果吻合较好。耦合二方程湍流模型摩阻系数预测结果相对三方程湍流模型结果误差为1%, 验证了数据驱动的SST-γ-ML耦合二方程湍流模型在转捩预测中的有效性。

3 数值算例与结果讨论

3.1 数值模型

以NACA0012二维翼型为研究对象, 应用第2节验证的SST-γ-ML耦合二方程湍流模型在Re=1.35×105下分别进行静态小迎角及动态失速过程的气动力计算。

NACA0012二维翼型弦长c=0.15 m, 翼型流场网格如图4所示。在翼型动态失速气动力计算中, 旋转域与固定域采用滑移网格以耦合不连续边界, 适用于旋转几何体数值仿真。将全流场划分为内部旋转区域(红圈内)与外部固定区域(红圈外)。内部区域为596×100的O形结构网格, 最大Y+为0.93; 外流场网格量为9 779的混合网格。为确保尾迹区流动的充分发展, 压力出口边界设置在后缘下游20倍弦长位置。

图4 NACA0012计算网格及边界条件Fig. 4 Mesh and boundary conditions of NACA0012

远场设置U∞=14 m/s,Tu∞=0.08%, 采用PIMPILE瞬态求解器进行压力速度解耦计算, 计算时间步长固定为1×10-5s。对SST-γ三方程湍流模型, 控制方程中对流项U、k与ω的散度计算均采用线性迎风离散格式,γ的散度计算采用1阶迎风离散格式。对SST-γ-ML耦合二方程湍流模型,γ无须进行控制方程计算, 其他参数离散格式与三方程湍流模型相同。

3.2 稳态小迎角

本节应用SST-γ-ML耦合二方程湍流模型对NACA0012二维翼型进行小迎角稳态气动力计算。

应用SST-γ三方程湍流模型与SST-γ-ML耦合二方程模型对2°～8°迎角流场进行计算, 耦合二方程模型中神经网络训练集迎角及数值仿真计算测试集迎角如表2所示。

表2 神经网络训练集与测试集迎角设置

耦合二方程湍流模型与三方程湍流模型在各迎角下升力系数计算结果如图5所示。

图5 不同迎角下升力计算结果Fig. 5 Lift coefficient at different angles of attack

耦合二方程湍流模型对未经训练的α=5°, 6°, 7°内插迎角升力系数计算准确: 与三方程模型计算结果对比, 平均相对误差为2.60%; 与实验结果对比[32], 平均相对误差为8.81%。未经训练的α=2°外插迎角结果, 对比实验结果相对误差为19.63%, 对比三方程结果相对误差为23.81%。

仅采用4°, 8°迎角数据训练的耦合二方程湍流模型在6°迎角下计算翼型表面压力系数曲线如图6所示。对未经训练的6°迎角情况, 耦合二方程湍流模型预测表面压力系数结果与三方程模型准确对应。对4°与8°迎角算例, 耦合二方程湍流模型同样可以准确预测压力系数。表明该数据驱动的转捩模型具有较好的迎角泛化能力。

图6 6°迎角下压力系数预测效果Fig. 6 Prediction of pressure coefficient at α=6°

3.3 动态失速

本节应用SST-γ-ML耦合二方程湍流模型对NACA0012翼型进行动态失速过程计算。设定翼型做给定正弦规律的俯仰运动

α(t)=10+15sin(18.67t)

其中,α(t)为t时刻翼型迎角。运动初始迎角为10°, 俯仰运动迎角范围为α∈[-5°, 25°], 俯仰运动周期T=0.337 s。

计算总时长为10T, 取第6～10周期进行相平均, 得到SST-γ三方程湍流模型与SST-γ-ML耦合二方程湍流模型在一个稳定周期内计算的升力系数。实验[32]、 三方程湍流模型、 LES模型[8]及耦合二方程湍流模型所得升力系数滞回曲线如图7所示。

图7 升力系数随迎角变化曲线Fig. 7 Curve of lift coefficient with angle of attack

升力系数的滞回曲线趋势与实验及三方程结果均符合较好, 与三方程模型计算结果相对误差为11.95%。从曲线上看, 升力系数曲线在翼型上仰过程的线性段计算准确; 在20°～25°迎角失速区间, 耦合二方程模型与三方程模型几乎同时达到第1次升力系数峰值, 耦合二方程模型计算的第2次升力系数峰值出现时间略晚于三方程模型; 在下俯过程中, 耦合二方程模型计算所得升力系数相对三方程结果波动小, 且取值偏低。

结合涡量云图, 对比翼型在上仰过程中耦合二方程湍流模型与三方程湍流模型对动态失速过程中典型流动状态的仿真结果, 如图8所示。

图8 翼型上仰过程典型流动状态涡量场: (a) 三方程湍流模型; (b) 耦合二方程湍流模型Fig. 8 Vortex fields of typical flow conditions via upstroke: (a) 3 equation turbulence model; (b) coupled 2 equation turbulence model

在翼型上仰过程中,α∈[-5°, 14°]迎角下, 耦合二方程模型对附着在翼型上表面的湍流边界层及后缘脱落涡仿真准确(图8(a-1; b-1)); 随迎角增大, 耦合二方程模型与三方程模型几乎同时出现前缘涡并迅速增长(图8(a-2, 3; b-2, 3))。前3个典型流动状态上, 耦合二方程模型数值仿真结果与三方程模型均准确对应, 该现象解释了升力系数曲线(图7)在上仰过程的α∈[-5°, 20°]迎角区间耦合二方程与三方程模型曲线的准确对应, 耦合二方程模型曲线的第1个升力系数峰与三方程曲线同时出现。在α∈[20°, 25°]区间, 前缘涡持续生成并脱落, 耦合模型对此期间前缘涡生成与脱落规律、 流动状态与涡量大小等实现了准确仿真, 但速率略慢于三方程仿真结果(图8(a-4; b-4)), 该速率差解释了升力系数曲线(图7)的α∈[20°, 25°]上仰区间里, 第2个升力系数峰的出现时间略晚于三方程曲线的误差原因。总体而言, 耦合二方程模型对翼型上仰过程仿真较为准确。

图9对比了翼型下俯过程中两个湍流模型对涡量场的仿真情况。

图9 翼型下俯过程典型流动状态涡量场: (a) 三方程湍流模型; (b) 耦合二方程湍流模型Fig. 9 Vortex fields of typical flow conditions via downstroke: (a) 3 equation turbulence model; (b) coupled 2 equation turbulence model

在翼型下俯过程中, 交替发生前缘涡与后缘涡的脱落(图9(a-1, 2, 3; b-1, 2, 3)); 在小迎角下, 上翼面尾缘流动再次附着(图9(a-4; b-4)); 此后, 边界层均附着在翼型表面(图9(a-5; b-5))。耦合二方程模型在下俯过程中对典型流动状态的仿真结果与三方程均能准确对应。

4 结论

本文建立了由湍动能k、 单位湍动能耗散率ω的输运方程及预测转捩间歇因子γ的神经网络构成的耦合二方程湍流模型(SST-γ-ML模型), 用于翼型动态失速过程的高精度方程。

1) 利用基于数据驱动的转捩间歇因子预测模型对T3A平板进行了流场仿真, 耦合二方程湍流模型对包含转捩过程的平板表面预测摩阻系数相对于SST-γ三方程模型结果平均误差为1%;

2) 对NACA0012翼型进行了低Reynolds数动态失速过程的流场仿真, 耦合二方程湍流模型预测的非定常气动升力与SST-γ三方程模型结果相比, 误差小于12%。