PREDICTION OF REYNOLDS STRESS ANISOTROPIC TENSOR BY NEURAL NETWORK WITHIN WIDE SPEED RANGE
-
摘要: 基于Pope修正的有效黏度假设, 张量基神经网络(tensor based neural network, TBNN)构建了从雷诺平均方程湍流模型(RANS)的平均应变率张量和平均旋转率张量到高精度数值解的雷诺应力各向异性张量的映射. 将高精度数值解用于TBNN的训练, 从而使TBNN根据RANS求解的湍动能、湍流耗散率和速度梯度预测其雷诺应力各向异性张量, 并与对应的高精度数值模拟结果以及风洞实验结果对比以评估TBNN的预测能力. 本工作将TBNN的预测能力从低速域拓展至高超声速工况, 分别对低速槽道流、低速 NACA0012 翼型以及高超声速平板边界层3种工况进行了小样本的训练并成功预测, 并以槽道流训练的TBNN较好地预测了低速平板边界层, 验证了模型的泛化能力. 对于外推的低速槽道流算例, TBNN预测的结果在y + >5的区域与直接数值模拟(DNS)以及实验的误差均在10%以内, 预测结果揭示了TBNN对雷诺应力各向异性张量的良好预测能力; 对于翼型的预测效果尽管相较于槽道流略有下降, 但近壁关键区域较RANS结果仍有显著提升; 对于高超声速平板, TBNN在边界层内展现出了良好的预测能力, 在y + >5的区域与DNS的误差同样在10%以内. 基于Pope本构关系的TBNN方法在平板的高超声速工况下仍能较准确预测边界层内的雷诺应力各向异性张量, 方法在宽速域下的预测能力具有较好的表现, 且模型泛化能力亦得到了验证.
-
关键词:
- 张量基神经网络 /
- 雷诺应力各向异性张量 /
- 湍流建模 /
- 宽速域 /
- 小样本
Abstract: Tensor based neural network (TBNN) is constructed based on Pope’s effective viscosity hypothesis, and it’s used to produces a mapping from the mean strain rate tensor, mean rotation rate tensor calculated by Reynolds averaged Navier-Stokes (RANS) to the high resolution Reynolds stress anisotropy tensor. The high resolution data is used to train TBNN, then TBNN will give prediction results of Reynolds stress anisotropic tensor from the RANS result. The prediction of TBNN will be compared with high resolution numerical simulation and wind tunnel results to evaluate the prediction ability of TBNN. This work expands the predictive ability of TBNN from the low speed domain to hypersonic conditions. Small sample training is performed on low speed channel flow, NACA0012 and hypersonic boundary layer and the prediction accuracy is satisfactory. In addition, the TBNN trained with channel flow accurately predicts the boundary layer of the low-speed flat plate, which verifies the generalization ability of the model. For the extrapolation channel flow at low-speed, TBNN can predict the Reynolds stress anisotropy tensor well in the range of y + > 5, the error between direct numerical simulation (DNS), experiment and TBNN is inside 10%. Although the prediction accuracy of the low-speed airfoil is slightly lower than that of the channel flow, the cloud images predicted in the key area have significant improvement compared with RANS. For the hypersonic boundary layer, TBNN shows good predictive ability in the boundary layer, and the error between TBNN and DNS is also within 10% in the range of y + > 5. Although Pope’s constitutive law is proposed for most incompressible flows, TBNN can still predict the Reynolds stress anisotropy tensor under hypersonic conditions. The predictive ability of this method in a wide speed range is confirmed and the generalization ability of the model has also been verified. -
引 言
湍流模型是为了封闭Navier-Stokes (N-S)方程中的雷诺应力项而额外补充的方程, 主要目的是构建时均流动、空间位置与雷诺应力张量或湍流涡粘之间的数学关系式[1]. 雷诺平均(Reynolds averaged Navier-Stokes, RANS)方程是当前湍流研究中使用最广泛的模型之一. RANS框架的核心是确定流场的平均量, 并针对脉动量对于平均场本身的影响进行建模; 而脉动场则是通过雷诺应力张量的散度嵌入到RANS方程中, RANS模型正是需要对此张量进行建模[2]. 最早的有效黏度假设由Boussinesq提出, 该假设也是使用最为广泛的有效黏度假设, 通常使用的两方程RANS模型(例如k−ε, k−ω模型)通过Boussinesq涡黏假设完成雷诺应力的封闭. 但Tracey等[3]指出湍流涡黏模型的主要误差来源是无法解释雷诺应力的各向异性, 即Boussinesq涡黏假设认为湍流黏性系数是各向同性的标量. 此外, 有诸多学者根据推导出的雷诺应力输运方程进行求解, 例如雷诺应力方程模型(RSM), 闫超等[4]指出RSM存在雷诺应力方程的建模困难、数值刚性问题较严重、计算量较大的问题. 不同于RSM的思路, Pope[5]针对各向异性提出了新的有效黏度假设. 该假设的有效性已经得到了充分的验证[6], 但由于其在N-S求解器中较差的鲁棒性, 并未得到广泛应用. 这也意味着研究者往往需要从DNS, LES等高分辨率数据中获得更加准确的雷诺应力各向异性张量.
随着机器学习方法的不断发展, 在雷诺应力的封闭中结合机器学习的研究方法日益增加. 主要的研究方向分为两种[7]: 采用机器学习方法缩小RANS解与高精度数据或实验值的偏差; 或是基于高精度数据或实验值直接构建某些湍流变量的代理模型. 缩小偏差的方法既可以是改变模型的控制方程形式, 也可以构建针对RANS模型偏差函数进行叠加修正. 代表性工作包括Tracey等[8]构建了替代SA模 型控制方程中源项的神经网络模型; Wu等[9]针对 RANS 模型计算结果和高分辨率DNS数据之间的雷诺应力偏差进行建模, 提高了原有模型的准确性. 相比于缩小偏差的方法, 直接构建代理模型的方法同样颇具亮点. Zhu等[10-11]通过径向基函数神经网络以及深度神经网络构建了涡黏系数νt与平均流动变量之间的映射关系, 并将得到的映射关系与N-S求解器耦合用于封闭湍流模型. 而Ling等[12]针对张量不变性分别在随机森林和神经网络中进行了测试. 随机森林作为决策树的集合, 是将多个决策树的预测组合成一个模型, 其结构并不像神经网络一样容易改变, 难以将Pope的本构关系嵌入到随机森林中. 而神经网络由于体系结构的灵活性在张量不变性的处理上别具优势. 为了在RANS模型的基础上, 利用机器学习方法提高雷诺应力各向异性的计算精度, Ling等[7]采用具有伽利略不变性的输入特征并提出张量基神经网络(TBNN)架构. 诸多成果表明机器学习方法在湍流模型化工作中具有良好的前景[13].
对于TBNN, Ling等[7]使用TBNN预测了低速的管流、周期山流动. Fang等[2]针对槽道流提出了新的神经网络模型, 使其在低速槽道流的预测效果优于TBNN. 而该研究也指出, 其新的神经网络未能嵌入旋转不变性, 且仅适用于特定流动. 张珍等[14]也将TBNN与一个预测涡黏系数的神经网络嵌入到了RANS求解器中, 并对低速周期山流动进行了预测, 提升了RANS求解的精度. 由于TBNN基于Pope的有效黏度假设构建, 而该假设是Pope针对绝大部分不可压缩流提出的, 因此TBNN的研究工作往往集中在低速的领域.
本文在蒙特利尔大学和谷歌开发的Theano深度学习框架(https://pypi.org/project/Theano/)上完成了神经网络的编译. 以Pope[5]修正的有效黏度模型作为理论基础, 并基于Ling等[7]搭建的TBNN内核构建了神经网络模型. 本文通过TBNN构建了从雷诺平均方程(RANS)湍流模型的湍动能、湍流耗散率和速度梯度到高精度数值解的雷诺应力各向异性张量的映射. 除了对低速槽道流、低速NACA0012翼型进行了预测以外, 还将TBNN的应用范围从低速拓展到了高超声速, 较精确地预测了马赫6的平板边界层的雷诺应力各向异性张量, 并通过一组低速平板进行了模型泛化能力的进一步验证. 验证了TBNN在宽速域流动下对雷诺应力各向异性的预测能力. 对于本文涉及的3个TBNN模型, 经由槽道流训练得的模型标记为TBNN-C, 翼型对应模型为TBNN-N, 高超声速平板对应模型为TBNN-H.
1. 计算方法
1.1 有效黏度假设
时均形式的N-S方程为
$$ \frac{{\partial \left( {\rho \overline {{u_i}} } \right)}}{{\partial t}} + \frac{{\partial \left( {\rho \overline {{u_i}{u_j}} } \right)}}{{\partial {x_j}}} = - \frac{{\partial \bar p}}{{\partial {x_i}}} + \frac{\partial }{{\partial {x_j}}}\left( {\mu \frac{{\partial \overline {{u_i}} }}{{\partial {x_j}}} - \rho \overline {u_i^\prime u_j^\prime } } \right) $$ (1) 式中ρ为密度,
$ \overline {{u_i}} $ 为时均形式下的各方向速度,$ \bar p $ 为时均的压力,$ \rho \overline {u_i^\prime u_j^\prime } $ 被称作雷诺应力. 雷诺应力作为对称的二阶张量, 其对角线上的分量$ \overline {u_1^\prime u_1^\prime } $ ,$ \overline {u_2^\prime u_2^\prime } $ 和$ \overline {u_3^\prime u_3^\prime } $ 称为正应力, 非对角线分量为切应力. 而湍动能k(x,t)则被定义为雷诺应力张量迹线的一半, 即$$ k \equiv \frac{1}{2}\overline {u_i^\prime u_i^\prime } = \frac{1}{2}\left( {\overline {u_1^\prime u_1^\prime } + \overline {u_2^\prime u_2^\prime } + \overline {u_3^\prime u_3^\prime } } \right) $$ (2) 主应力和切应力的不同取决于坐标系的选择, 而根据雷诺应力内在的差别可以将其区分为各向同性和各向异性
$$ \overline {u_i^\prime u_j^\prime } = \frac{2}{3}k{\delta _{ij}} + {a_{ij}} $$ (3) 式中
$ k{\delta _{ij}} $ 为各向同性张量,$ {a_{ij}} $ 为各向异性张量.为封闭RANS方程, 需要给出雷诺应力与平均流场间的关系, 这一封闭的过程可基于有效黏度假设完成. 最早的有效黏度假设由Boussinesq提出, 该假设也是使用最为广泛的有效黏度假设. 基于该假设, 雷诺应力可被定义为
$$ \overline {u_i^\prime u_j^\prime } = \frac{2}{3}k{\delta _{ij}} - {\mu _{{\text{eff}}}}\left( {{U_{i,j}} + {U_{j,i}}} \right) $$ (4) 式中k为湍动能,
$ {\delta _{ij}} $ 为克罗内克符号,$ {\mu _{{\text{eff}}}} $ 为有效黏度, Ui,j和Uj,i为速度梯度. 发现Boussinesq假设是基于各向同性假设的, 这也使得该假设无法准确地捕捉各向异性.为了更好的捕获雷诺应力各向异性张量, Pope[5]进一步修正了有效黏度模型. 将各向异性张量通过湍动能进行归一化
$$ {b_{ij}} = {a_{ij}}/(2k) $$ (5) Pope[5] 根据Caley–Hamilton理论推导了归一化的各向异性张量b与基张量之间的本构关系的一般形式
$$ {\boldsymbol{b}} = \sum\limits_{n = 1}^{10} {{g^{(n)}}\left( {{\lambda _1},{\lambda _2} \ldots ,{\lambda _5}} \right)} {{\boldsymbol{T}}^{(n)}} $$ (6) 式中
${{\boldsymbol{T}}^{(1)}}, {{\boldsymbol{T}}^{(2)}},\cdots ,{{\boldsymbol{T}}^{(10)}}$ 为基张量,${\lambda _1}, {\lambda _2}\cdots ,{\lambda _5}$ 为张量不变量, 二者都是无量纲化后的张量S和R相关的函数. S为平均应变率张量, R为平均旋转率张量. 不同于三维流场中的10个基张量, 该本构关系在二维流场中仅需要4个基张量, 具体表达为$$ \left. \begin{array}{l} {\boldsymbol{b}}={\displaystyle \sum _{n=1}^{4}{g}^{(n)}\left({\lambda }_{1},{\lambda }_{2}\cdots ,{\lambda }_{5}\right)}{{\boldsymbol{T}}}^{(n)}\\ {{\boldsymbol{T}}}^{(1)}={\boldsymbol{S}}\\ {{\boldsymbol{T}}}^{(2)}={\boldsymbol{SR}}-{\boldsymbol{RS}}\\ {{\boldsymbol{T}}}^{(3)}={{\boldsymbol{S}}}^{2}-\dfrac{1}{3}{\boldsymbol{I}}\cdot{{{\rm{Tr}}}}\left({{\boldsymbol{S}}}^{2}\right)\\ {{\boldsymbol{T}}}^{(4)}={{\boldsymbol{R}}}^{2}-\dfrac{1}{3}{\boldsymbol{I}} \cdot {{{\rm{Tr}}}}\left({{\boldsymbol{R}}}^{2}\right) \end{array} \right\} $$ (7) 对于张量不变量的构造, Johnson[15]枚举了对称和反对称张量的7个相关不变量: Tr(S), Tr(S2), Tr(S3), Tr(R2), Tr(R2S), Tr(R2S2)和Tr(R2SRS2). 而Pope将本构关系中的张量不变量确定为
$$\left. \begin{array}{l} {\lambda _1} = {{{\rm{Tr}}}}\left( {{{\boldsymbol{S}}^2}} \right) \hfill \\ {\lambda _2} = {{{\rm{Tr}}}}\left( {{{\boldsymbol{R}}^2}} \right) \hfill \\ {\lambda _3} = {{{\rm{Tr}}}}\left( {{{\boldsymbol{S}}^3}} \right) \hfill \\ {\lambda _4} = {{{\rm{Tr}}}}\left( {{{\boldsymbol{R}}^2}{\boldsymbol{S}}} \right) \hfill \\ {\lambda _5} = {{{\rm{Tr}}}}\left( {{{\boldsymbol{R}}^2}{{\boldsymbol{S}}^2}} \right) \hfill \end{array} \right\} $$ (8) 无量纲化后的张量S和R表达式为
$$ {S_{ij}} = \frac{k}{{2\varepsilon }}\left( {\frac{{\partial {{\bar u}_i}}}{{\partial {x_j}}} + \frac{{\partial {{\bar u}_j}}}{{\partial {x_i}}}} \right) $$ (9) $$ {R_{ij}} = \frac{k}{{2\varepsilon }}\left( {\frac{{\partial {{\bar u}_i}}}{{\partial {x_j}}} - \frac{{\partial {{\bar u}_j}}}{{\partial {x_i}}}} \right) $$ (10) 由于张量不变量和基张量的存在, 任何满足该有效黏度模型的本构关系的张量b都会自动满足伽利略不变性.
1.2 深度神经网络
全连接神经网络的大致结构如图1所示. 每一个圆圈代表神经元, 其中包含两步计算: 第一步为输入向量X和权值向量W及偏置b的线性运算a = WX + b, 第二步为激活函数y=h(a)的非线性运算. 常用的激活函数h(x)包括Tanh, Sigmoid, ReLU和阶跃函数等. 本文沿用了TBNN中采用了leaky ReLU激活函数[6].
通常采用带有梯度下降的反向传播方法来训练神经网络, 本文选择了Adam方法. 该方法通过迭代调整神经网络中各节点的权重, 从而在各向异性张量的预测值和真实值之间提供最低的均方误差, 以使模型更适合训练数据. 此外隐藏层的层数、每层的宽度也会显著影响神经网络的性能. 神经网络的层数和节点数较多时能容纳更复杂的数据, 但也容易造成过拟合.
1.3 张量基神经网络
TBNN目标是确定式(6)中的
$ {g^{(n)}}\left( {{\lambda _1},{\lambda _2}, \cdots ,{\lambda _5}} \right) $ 函数. 一旦确定了函数, 则可以使用公式(6)求解张量b. 对于机器学习而言, 将相同的流动在不同坐标系下的流场用于机器学习时, 机器学习模型可能会产生不同的预测, 因此在使用TBNN时保证神经网络的伽利略不变性是必要的.基于Pope的有效黏度模型, 将5个张量不变量
$ {\lambda _1},{\lambda _2}, \cdots ,{\lambda _5} $ 作为TBNN的输入变量, 标量函数$ {g^{(n)}} $ 作为输出变量, 从而避免在神经网络中引入张量运算, 保证了其伽利略不变性. 在完成了TBNN模型的构建后, 结合基张量进行式(7)的运算, 即可得到归一化的雷诺应力各向异性张量b. TBNN的结构图如图2所示.此外为了衡量训练效果, 引入均方根误差作为损失函数
$$ \begin{split} &{RMSE}=\\& \sqrt{\dfrac{1}{4{N}_{\text{data}}}\displaystyle\sum _{{m}=1}^{{N}_{\text{data}}}\left[\displaystyle\sum _{i=1}^{3}\displaystyle\sum _{j=1}^{i}{\left({b}_{ij,m,\text{predicted}}-{b}_{ij,m,{\rm{DNS}}}\right)}^{2}\right]} \end{split} $$ (11) 学习率采用学习率衰减(learning rate decay), 初始学习率为0.01, 最小学习率为1 × 10−6.
并根据Banerjee等[16]提出的湍流不变量与特征值方法的各向异性特性的表述, 对b和其特征值
$ {\xi _1} \geqslant {\xi _2} \geqslant {\xi _3} $ 添加约束, 从而对预测好的各向异性进行后处理$$\left. \begin{split} & - \frac{1}{3} \leqslant {b_{ii}} \leqslant \frac{2}{3},\quad \quad \quad \quad {\xi _1} \geqslant (3\left| {{\xi _2}} \right| - {\xi _2})/2 \hfill \\ & - \frac{1}{2} \leqslant {b_{ij}} \leqslant \frac{1}{2}\quad {\text{for}}\;j \ne i,\;\;\quad \quad {\xi _1} \leqslant \frac{1}{3} - {\xi _2} \hfill \end{split} \right\}$$ (12) 2. 数据集
Ling等[7]的工作对多种基础流动进行了研究, 诸如管流(Reb=3500)、周期山流动等. 张珍等[14]使用EVNN和TBNN的组合神经网络进一步研究了周期山流动. 为了验证基于Pope有效黏度模型的TBNN在宽速域和小样本的前提下的适用范围, 本文选取了低速槽道流、低速翼型以及高超声速平板三类算例展开研究. 由于获得的高精度数据均对应二维, 本文采用翼型、槽道流以及高超平板的二维算例, 以精确匹配高精度DNS和LES数据的二维工况.
针对低速工况, 选择NACA0012翼型以及槽道流作为研究算例. 槽道流的DNS数据来自Moser等[17-19], 用于TBNN训练的为Reτ = 395工况下湍流充分发展区域的一条截线上136个点的流场数据, Reτ = 590工况对应的136个点将作为外推的测试集验证TBNN预测能力; 针对Reτ = 590, Krogstad等[20]在风洞实验中通过热线风速仪测量了该工况的雷诺应力各向异性张量, 相关实验中的18个流场点可作为测试集的验证, 与TBNN预测结果和DNS结果形成充分的对比, 进一步衡量TBNN预测能力. 翼型来自Vinuesa等[21-22]通过LES求解的NACA0012翼型, 工况为Rec = 4×105, AoA = 0º, 不可压缩. 数据包含了流场上翼面的1800个数据点, 其中80%的点作为训练数据, 余下20%的点将通过散点图和云图的形式对TBNN的预测效果进行验证, 对比测试中尚不包含实验结果.
高超声速流动涉及了强激波、强逆压梯度、强压缩效应等诸多因素[23], 因此Pope[5]针对大部分不可压缩流动提出的本构关系在高超是否适用需要重新判断. 本文选择一组Ma = 6的平板作为研究算例. 该算例来自Sun等[24]最新计算的高超声速平板, 以边界层内的流场为研究对象, 并于该平板上选取了8条截线共计1540个流场点进行研究, 其中的6条截线用于训练, 相关DNS计算方法[25]以及平板计算结果[24]见以下文章. 而相较于低速槽道流的实验, 高超声速工况下的雷诺应力难以使用热线风速仪进行测量. 其它测量手段如纳米示踪平面激光散射技术(NPLS)[26]可以获得超声速、高超声速条件下的雷诺应力分布, 但由于本文需要进行边界层内不同站位法向方向雷诺应力各向异性分量的定量对比, 故仅采用DNS计算结果与本文计算结果进行对照.
本文中用于TBNN的训练集包含两部分. 第一部分为RANS求解所得结果, 其中湍流模型采用k−ε, SST模型进行求解. 根据流场的湍动能、湍流耗散率和速度梯度构造无量纲化后的张量S和R, 进而计算出5个张量不变量作为神经网络输入. 第二部分为高精度的雷诺应力数据, 经求解得雷诺应力各向异性张量, 作为TBNN的输出. 在经过训练集完成训练后, 便可运用TBNN模型对筛选出的测试集进行预测.
在完成RANS计算后, 可基于RANS求解结果并根据公式(12)对雷诺应力各向异性张量进行初步的反推[7]
$$ {a_{ij}} = - 2{\mu _t}\frac{{{S_{ij}}}}{{2k}} = \left\{ \begin{gathered} - {C_\mu }\frac{k}{\varepsilon }{S_{ij}},\;\;\;\; k - \varepsilon \;{\text{model}} \hfill \\ - \frac{{{S_{ij}}}}{\omega },\;\;\;\;k - \omega \;{\text{model}} \hfill \\ \end{gathered} \right. $$ (13) 通过此式亦可对Boussinesq涡黏假设在各向异性的预测能力进行初步的判断.
3. 数值模拟结果
本文主要研究雷诺应力各向异性张量的预测, 需要完成RANS求解、训练集构建以及TBNN的训练和预测.
3.1 低速槽道流
低速槽道流的RANS计算由k−ε模型完成. 选择Reτ = 395的工况作为训练集, 并将训练好的TBNN模型命名为TBNN-C. 将TBNN-C用于Reτ = 590的工况进行外推算例的预测. 开源的DNS数据中[17-19]给出了完全发展后的一个截线, 共有136个流场点; 两种工况在RANS计算的网格细节如表1所示.
表 1 槽道流RANS计算的网格信息Table 1. Mesh information for channel flow in RANSRe Mesh size Lx Ly Δx + Δyc + Reτ = 395 1024 × 257 8πδ 2δ 10.0 6.5 Reτ = 590 1024 × 257 8πδ 2δ 9.7 7.2 表1中Lx为计算域在x方向的长度, δ为槽道的半高, Ly = 2δ为y方向的长度. Δx + 为x方向两个网格点的距离差, Δyc + 为y方向中心位置两个网格点的距离差. 此处针对RANS所用网格的进行了网格无关性验证. 在确保了网格不会对求解精度造成影响后, 选取其中湍流充分发展的位置作截线, 与DNS数据中的流场点进行匹配. 由于槽道流的对称性, 只取下半部分作为研究对象, 网格尺寸为1024 × 257. 在保证RANS求解的速度型与DNS相近之后, 将RANS求得的速度梯度、湍动能和湍流耗散率以及DNS中的雷诺应力张量输入到TBNN-C中. 隐藏层的层数分别设置为2, 4和6层并对训练结果的均方根误差进行对比, 最终选择6层, 每层20个节点.
该工况的实验由Krogstad等[20]在回流式风洞中完成, 通过两个平行的光滑平板形成了槽道, 试验段长5 m, 入口面积1.35 m × 0.10 m, 通过热线风速仪测量速度脉动, 从而获得雷诺应力. 图3为Reτ = 590的工况下, 雷诺主应力的各分量随着y + 的变化. y + 的计算公式如式(14)所示, 其中
$ {u_\tau } $ 为摩擦速度(即湍流壁面附近的黏性速度量级),$ v $ 为运动黏度![]()
图 3 Reτ = 590雷诺应力各向异性分量随y + 值的分布Figure 3. Reynolds normal stress anisotropy components vs y + at Reτ = 590$$ {y^ + } = \frac{{{u_\tau }y}}{v} $$ (14) 区别于雷诺应力张量, 从图3结果我们可以发现, TBNN-C的预测结果在y + 较小的区域(即y + <5的黏性子层区域)有较大偏离, 但在过渡子层以及完全湍流区域TBNN-C的预测效果与DNS以及风洞实验的结果相近, 预测效果良好, 能够揭示雷诺正应力的各向异性部分在槽道流中的发展规律. 在y + >5的范围内, TBNN-C预测的bii与DNS结果的整体误差仅有10%左右, 与实验相同点位上TBNN-C的预测误差均不超过10%; 而与TBNN相比, RANS的预测结果存在量级上的差异, 也说明了基于Boussinesq涡黏假设的湍流k−ε湍流模型无法准确捕捉雷诺应力各向异性张量.
此外, 针对槽道流的TBNN-C模型仅根据Reτ = 395的136个流场点进行训练, 在样本极小的情况下, 训练完成的TBNN-C模型针对Reτ = 590的不同流场外推预测依旧展现出极佳的泛化能力. 这也说明对于槽道流这种外形简单的研究对象, Pope[5]提出的本构关系能够充分诠释其流动机理, 从而极大的提升了TBNN模型的泛化能力.
3.2 低速NACA0012翼型
本文所研究的NACA翼型为NACA0012. 由于NACA0012在零度攻角下的对称性, 仅研究其上翼面. 大涡模拟方法(LES)作为研究复杂湍流问题的重要工具, 其求解结果拥有较高的精度[27-28], 可作为神经网络的训练标签. 在开源的LES数据中, 共有1800个流场点[21-22], 皆位于机翼上翼面近壁区. 翼型RANS计算选择的湍流模型为SST模型. 在完成RANS计算后, 将其中的流场点提取并与LES流场点一一匹配. 为验证RANS计算精度, 图4中对比了RANS解得的速度场云图以及LES的1800个数据点插值所得速度场云图, RANS的计算结果能够初步的反映LES结果中的速度变化趋势.
![]()
图 4 RANS与LES的速度场对比(续)Figure 4. Comparison of velocity field for RANS and LES (continued)在LES数据中选择1440个流场点作为训练集对TBNN进行训练, 训练完成后的TBNN模型命名为TBNN-N. 将额外的360个流场点作为测试集输入到TBNN-N模型中, 检验神经网络对于翼型的预测效果. 测试集测点分布如图5所示.
在对隐藏层的层数进行了测试后, 选择三层隐藏层, 每层20个节点进行训练. 训练步数为8000步, 最终训练完成时训练集的均方根误差为0.052, 训练的样本数为1440, 样本数量较小. 图6和图7给出了测试集的360个流场点的各向异性分量b11, b22的分布, 并选取关键区域的流场点插值绘制云图. 插值方法为Matlab中的V4双调和样条插值(biharmonic spline interpolation).
![]()
图 6 雷诺正应力各向异性分量b11结果Figure 6. The result of Reynolds normal stress anisotropy component b11对于NACA0012翼型, 基于RANS求解的平均速度场几乎完全无法捕捉雷诺应力各向异性分量, 而基于RANS结果预测的TBNN-N各向异性分布在个别流场点的预测较差, 但在重点关注区域TBNN-N能够较好地预测雷诺应力各向异性. 同样在低速工况下, TBNN-N针对翼型的预测相比于槽道流精度有所下降, 但在小样本的前提下TBNN-N取得的预测效果尚可接受, 较RANS结果显著提升. 要对翼型进行精准预测, 可考虑进一步扩大样本范围, 以提升TBNN的预测精度.
![]()
图 7 归一化雷诺正应力各向异性分量b22云图分布Figure 7. The result of normalized Reynolds normal stress anisotropy component b223.3 高超声速平板
针对低速工况的模拟方法面对高超声速往往会产生不适用性[29-30], 故通过高超声速平板来衡量TBNN在高超声速下预测的精确度. 当前所拥有的平板DNS数据为二维算例[24], 来流参数中马赫数为6, 单位雷诺数为12000, 在湍流区域共有8条截线. 选取截线位置边界层内的数据用于本算例, 其中B11距平板前缘点51.68%处, 各截线位置以及边界层厚度如图8所示.
![]()
图 8 平板DNS数据各截线位置以及边界层厚度Figure 8. Position of transversals and boundary layer thickness of DNS data在8条截线中, 选择B13和B16作为测试集, 其余6组截线数据作为训练集对TBNN进行训练, 训练完成后的模型命名为TBNN-H. 训练集中共有1155个流场点, 用于预测的B13和B16截线分别有187和199个流场点.
高超声速平板的RANS计算由Fluent完成, 使用可压缩k−ε模型作为湍流模型. 根据边界层厚度在RANS结果中完成B11到B18的各截线位置的匹配. RANS计算的网格尺寸为3869 × 320, 在完成网格无关性验证后, 保证相同站位RANS与DNS 边界层厚度一致, 对比结果如图9所示. 此外, 如图9还给出了边界层内近壁区域的网格加密.
根据DNS计算结果, 截线B13的边界层厚度δ = 6.77 mm, B16的边界层厚度δ = 7.91 mm. 在边界层厚度相同的情况下, 相比于DNS计算结果, RANS的速度发展较慢. 以此RANS结果为基础的TBNN-H依旧取得了极佳的预测效果, 如图10和图11所示. (由于流场点过多, 进行数据点显示控制以美观曲线, 每两个点显示一个点)
![]()
图 10 截线B13的雷诺主应力各向异性分量的分布Figure 10. Reynolds normal stress anisotropy components on transversals B13针对高超平板的TBNN-H使用了4个隐藏层, 每层18个神经元. 在训练完成后, 训练集的均方根误差(root mean square error, RMSE)维持在0.01左右, 而训练好的TBNN-H对于测试集的预测表现良好, RMSE约为0.015左右. 从结果中可以发现, 在y + 较小的区域(即y + <5的黏性子层区域)TBNN-H的预测结果有一定偏离; 与槽道流相同, TBNN-H在过渡子层以及完全湍流区域TBNN-H的预测效果与DNS的结果相近, 预测效果良好. 这也说明了TBNN的预测在一定程度上可能会受到k−ε模型的固有限制, 因此在湍流边界层的黏性子层的预测效果较差.
![]()
图 11 截线B16的雷诺主应力各向异性分量的分布Figure 11. Reynolds normal stress anisotropy components on transversals B16此外, 从图10和图11的结果中我们可以发现TBNN-H对于高超声速平板边界层依旧具有良好的预测能力, 尽管Pope的本构关系针对不可压缩流动提出, 但基于该本构关系构造的TBNN-H依旧能够对强压缩性的高超声速平板流动雷诺应力各向异性张量进行预测. 而这也意味着通过高超声速样本构建的TBNN模型可以用来提高可压缩模型的针对性, 为高超声速湍流模型的定制化提供方法基础[1].
3.4 模型泛化能力验证
神经网络模型的泛化能力至关重要, 本节将对训练好的TBNN模型能否应用于与训练算例不同的算例展开讨论. 在3.1节中, TBNN-C模型由Reτ = 395的低速槽道流完成训练, 并对Reτ = 590的槽道流工况进行了较好的预测, 初步验证了TBNN的外推能力. 为了进一步验证其泛化能力, 本节选取了一组Reθ = 1100的低速平板作为验证算例, 该平板的DNS结果由Jimenez等[31]计算. 将Reτ = 395的低速槽道流作为训练数据, 通过训练好的TBNN-C模型对该低速平板进行预测.
验证算例的开源DNS数据[31]给出了边界层厚度为δ99 = 2.7568的对应截线上的数据点, 其中边界层内的流场点共有130个. 该算例的RANS计算由standard k−ε模型完成, 网格参数见表2.
表 2 验证算例RANS计算的网格信息Table 2. Mesh information for validation case in RANSMomentum thickness Reθ Mesh size Lx Ly △x + Reθ = 1100 1226 × 359 224.58 12.17 0.183 表中Lx为计算域在x方向的长度, Ly为y方向的长度, Δx + 为x方向两个网格点的距离差. RANS计算的y方向的网格分布与DNS计算保持一致. 此处针对RANS所用网格的进行了网格无关性验证. 在确保了网格不会对求解精度造成影响后, 选取其中湍流充分发展的位置作截线, 与DNS数据中的流场点进行匹配以便于最终进行对比.
将RANS算得的平板边界层的湍动能、湍流耗散率以及速度梯度输入到TBNN-C模型中, 得到预测的结果. 同时将预测结果与DNS结果以及通过公式(13)解得的RANS结果进行对比, 如图12所示.
![]()
图 12 Reθ = 1100平板边界层的雷诺主应力各向异性分量的分布Figure 12. Reynolds normal stress anisotropy components on Reθ = 1100 boundary layer相比于3.1节的结果, 基于槽道流的136个流场点训练的TBNN-C模型对于不同工况的平板边界层的预测精度尽管稍有下降, 但与DNS结果相比误差仍不超过20%, 依旧较为准确的预测了雷诺应力各向异性张量, 验证了TBNN模型的泛化能力. 而使用上述模型预测翼型以及高超声速平板边界层, 结果显示模型的预测的精度大幅下滑, 在一定程度上说明了模型性能对于训练数据的依赖. 尽管TBNN是基于Pope的本构关系构建, 但依旧会受到神经网络特性的固有限制.
4. 结 论
本文宽速域下TBNN对雷诺应力各向异性张量的预测结论如下.
(1) 无论是低速还是高超声速, 基于Boussinesq有效黏度假设的RANS模型均难以准确地捕捉雷诺应力各向异性张量.
(2) 基于Pope[5]提出的有效黏度假设构造出的TBNN-C对于槽道流的预测结果较好, 在训练样本仅有136个流场点的极小样本的情况下, 依旧可以准确地预测不同雷诺数下槽道流的雷诺应力各向异性张量. TBNN预测结果与DNS误差在10%左右, 部分点位与风洞实验结果误差亦在10%以内, 泛化能力较好.
(3) TBNN-N对于NACA0012翼型的预测在个别流场点存在偏差. 在选取了流场中1440个点作为训练集之后, TBNN对于同流场额外的360个点的预测结果尚可. 在小样本的训练前提下, TBNN-N能够对NACA0012的关键区域进行较为准确的预测. 与槽道流相比, TBNN对于复杂流场需要增大样本量以提升准确性.
(4) 对于高超声速平板, 在以湍流区域部分位置进行训练后, TBNN-H能够较好地给出流场中其他位置的雷诺应力各向异性张量分布. 尽管Pope[5]有效黏度假设是针对不可压缩流提出的, TBNN-H依旧可以在小样本的前提下对高超声速平板边界层内的流场进行精准的预测.
(5) 在低速槽道流和低/高超声速平板中, TBNN在黏性子层(即y + <5区域)表现较差, 猜测原因可能是训练集的RANS部分来自于k−ε模型求解, TBNN的预测能力会受到RANS结果的固有限制.
(6) 低速槽道流训练的TBNN-C模型能够较为精确的预测低速平板算例, 模型的泛化能力得到了验证. 而相同模型对于翼型和高超声速平板的预测出现的预测精度下滑则在一定程度上说明了模型性能对于训练数据的依赖.
本文在宽速域下, 对TBNN预测能力进行了充分的验证. 相比于Pope针对不可压缩工况提出的本构关系, 神经网络凭借出色的从数据中提取信息的能力在映射关系中学习到了可压缩相关的信息, 使得基于Pope本构关系的TBNN能够对高超声速工况进行较好的预测. 尽管神经网络的“黑匣子”性质使得TBNN难以像本构关系一样具有可解释性, 但本工作对于高超声速的雷诺应力各向异性的求解仍具有参考价值.
在下一步的工作中, 将针对更多的高超声速复杂工况展开研究, 同时增加三维的相关预测. 在进一步完成TBNN在复杂高超声速工况的适用性验证后, 将通过单向耦合的方式将训练好的TBNN模型与可压缩湍流模型以及N-S求解器耦合, 尝试改善对于流动的预测精度. 而耦合对应的稳定性与收敛性的问题将是研究的重点与难点. 此外, 还将通过在神经网络中添加物理约束, 构建物理驱动的张量基神经网络(physics informed tensor based neural network), 进一步增强TBNN在小样本的泛化能力.
致 谢
感谢Julia Ling于https://github.com/tbnn/tbnn开源的TBNN内核代码; 感谢Ricardo Vinuesa提供的Naca0012翼型LES数据和帮助; 感谢Robert D. Moser提供的槽道流DNS数据以及Jimenez提供的低速平板DNS数据.
-
![]()
图 3 Reτ = 590雷诺应力各向异性分量随y + 值的分布
Figure 3. Reynolds normal stress anisotropy components vs y + at Reτ = 590
![]()
图 4 RANS与LES的速度场对比(续)
Figure 4. Comparison of velocity field for RANS and LES (continued)
![]()
图 6 雷诺正应力各向异性分量b11结果
Figure 6. The result of Reynolds normal stress anisotropy component b11
![]()
图 7 归一化雷诺正应力各向异性分量b22云图分布
Figure 7. The result of normalized Reynolds normal stress anisotropy component b22
![]()
图 8 平板DNS数据各截线位置以及边界层厚度
Figure 8. Position of transversals and boundary layer thickness of DNS data
![]()
图 10 截线B13的雷诺主应力各向异性分量的分布
Figure 10. Reynolds normal stress anisotropy components on transversals B13
![]()
图 11 截线B16的雷诺主应力各向异性分量的分布
Figure 11. Reynolds normal stress anisotropy components on transversals B16
![]()
图 12 Reθ = 1100平板边界层的雷诺主应力各向异性分量的分布
Figure 12. Reynolds normal stress anisotropy components on Reθ = 1100 boundary layer
表 1 槽道流RANS计算的网格信息
Table 1 Mesh information for channel flow in RANS
Re Mesh size Lx Ly Δx + Δyc + Reτ = 395 1024 × 257 8πδ 2δ 10.0 6.5 Reτ = 590 1024 × 257 8πδ 2δ 9.7 7.2 
下载: 导出CSV
表 2 验证算例RANS计算的网格信息
Table 2 Mesh information for validation case in RANS
Momentum thickness Reθ Mesh size Lx Ly △x + Reθ = 1100 1226 × 359 224.58 12.17 0.183
下载: 导出CSV
-
[1] 张伟伟, 寇家庆, 刘溢浪. 智能赋能流体力学展望. 航空学报, 2021, 42(4): 524689 (Zhang Weiwei, Kou Jiaqing, Liu Yilang. Prospect of artificial intelligence empowering fluid mechanics. Acta Aeronautica et Astronautica Sinica, 2021, 42(4): 524689 (in Chinese) [2] Fang R, Sondak D, Protopapas P, et al. Neural network models for the anisotropic Reynolds stress tensor in turbulent channel flow. Journal of Turbulence, 2019, 21(9-10): 525-543
[3] Tracey B, Duraisamy K, Alonso J. Application of supervised learning to quantify uncertainties in turbulence and combustion modeling//AIAA Aerospace Sciences Meeting, 2013-0259
[4] 阎超, 屈峰, 赵雅甜等. 航空航天CFD物理模型和计算方法的述评与挑战. 空气动力学学报, 2020, 38(5): 829-857 (Yan Chao, Qu Feng, Zhao Yatian, et al. Review of development and challenges for physical modeling and numerical scheme of CFD in aeronautics and astronautics. Acta Aerodynamica Sinica, 2020, 38(5): 829-857 (in Chinese) [5] Pope SB. A more general effective-viscosity hypothesis. Journal of Fluid Mechanics, 1975, 72: 331-340
[6] Gatski T, Speziable C. On explicit algebraic stress models for complex turbulent flows. Journal of Fluid Mechanics, 1993, 254: 59-78
[7] Ling J, Kurzawski A, Templeton J. Reynolds averaged turbulence modelling using deep neural networks with embedded invariance. Journal of Fluid Mechanics, 2016, 807: 155-166
[8] Tracey BD, Duraisamy K, Alonso JJ. A machine learning strategy to assist turbulence model development//Proceedings of the 53rd AIAA Aerospace Sciences Meeting, Kissimmee, Florida, 2015
[9] Wu J, Xiao H, Paterson E. Physics-informed machine learning approach for augmenting turbulence models: A comprehensive framework. Physical Review Fluids, 2018, 3(7): 074602
[10] Zhu LY, Zhang WW, Kou J. Machine learning methods for turbulence modeling in subsonic flows around airfoils. Physics of Fluids, 2019, 31(1): 015105
[11] Zhu LY, Zhang WW, Sun X, et al. Turbulence closure for high Reynolds number airfoil flows by deep neural networks. Aerospace Science and Technology, 2021, 110: 106452
[12] Ling J, Jones R, Templeton J. Machine learning strategies for systems with invariance properties. Journal of Computational Physics, 2016, 318: 22-35
[13] 袁先旭, 陈坚强, 杜雁霞等. 国家数值风洞工程中的CFD基础科学问题研究进展. 航空学报, 2021, 42(9): 625733 (Yuan Xianxu, Chen Jianqiang, Du Yanxia, et al. Research progress on fundamental CFD issues in national numerical windtunnel project. Acta Aeronautica et Astronautica Sinica, 2021, 42(9): 625733 (in Chinese) [14] 张珍, 叶舒然, 岳杰顺等. 基于组合神经网络的雷诺平均湍流模型多次修正方法. 力学学报, 2021, 53(6): 1532-1542 (Zhang Zhen, Ye Shuran, Yue Jieshun, et al. A combined neural network and multiple modification strategy for Reynolds-Averaged Navier-Stokes turbulence modeling. Chinese Journal of Theoretical and Applied Mechanics, 2021, 53(6): 1532-1542 (in Chinese) doi: 10.6052/0459-1879-21-073 [15] Johnson R. The Handbook of Fluid Dynamics, CRC Press, 1998
[16] Banerjee S, Kraal R, Durst F, et al. Presentation of anisotropy properties of turbulence invariants versus eigenvalue approaches. Journal of Turbulence, 2007, 8: 1-27
[17] Moser RD, Kim J, Mansour NN. Direct numerical simulation of turbulent channel flow up to Reτ=590. Physics of Fluids, 1999, 11: 943-945
[18] Kim J, Moin P, Moser R. Turbulence statistics in fully developed channel flow at low Reynolds number. Journal of Fluid Mechanics, 1987, 177: 133-166
[19] Lee M, Moser R. Direct numerical simulation of turbulent channel flow up to. Journal of Fluid Mechanics, 2015, 774: 395-415
[20] Krogstad P, Andersson H, Bakken O, et al. An experimental and numerical study of channel flow with rough walls. J. Fluid Mech., 2005, 530: 327-352
[21] Tanarro Á, Vinuesa R, Schlatter P. Effect of adverse pressure gradients on turbulent wing boundary layers. Journal of Fluid Mechanics, 2019, 883: A8
[22] Vinuesa R, Negi P, Atzori M, et al. Turbulent boundary layers around wing sections up to Rec=1000,000. International Journal of Heat and Fluid Flow, 2018, 72: 86-99
[23] 陈坚强, 袁先旭, 涂国华等. 高超声速边界层转捩的几点认识. 中国科学:物理学 力学 天文学, 2019, 49: 114701 (Chen Jianqiang, Yuan Xianxu, Tu Guohua, et al. Recent progresses on hypersonic boundary-layer transition. Sci. Sin.-Phys. Mech. Astron., 2019, 49: 114701 (in Chinese) doi: 10.1360/SSPMA-2019-0071 [24] Sun D, Guo QL, Yuan XX, et al. Decomposition formula for the wall heat flux of a compressible boundary layer. arXiv preprint, 2106.10968
[25] Sun D, Chen JQ, Yuan XX, et al. On the wake structure of a micro-ramp vortex generator in hypersonic flow. Physics of Fluids, 2020, 32: 126111
[26] 易仕和, 陈植, 朱杨柱等. (高)超声速流动试验技术及研究进展. 航空学报, 2015, 36(1): 98-119 Yi Shihe, Chen Zhi, Zhu Yangzhu, et al. Progress on experimental techniqgues and studies of hypersonic/supersonic flows. Acta Aeronautica et Astronautica Sinica, 2015, 36(1): 98-119 (in Chinese)
[27] 谢晨月, 袁泽龙, 王建春等. 基于人工神经网络的湍流大涡模拟方法. 力学学报, 2021, 53(1): 1-16 (Xie Chenyu, Yuan Zelong, Wang Jianchun, et al. Artificial neural network-based subgrid-scale models for large-eddy simulation of turbulence. Chinese Journal of Theoretical and Applied Mechanics, 2021, 53(1): 1-16 (in Chinese) [28] 时北极, 何国威, 王士召. 基于滑移速度壁模型的复杂边界湍流大涡模拟. 力学学报, 2019, 51(3): 754-766 (Shi Beiji, He Guowei, Wang Shizhao. Large-eddy simulation of flows with complex geometries by using the slip-wall model. Chinese Journal of Theoretical and Applied Mechanics, 2019, 51(3): 754-766 (in Chinese) doi: 10.6052/0459-1879-19-033 [29] 向星皓, 陈坚强, 袁先旭等. C-γ-Reθ高超声速三维边界层转捩预测模型研究. 航空学报, 2021, 42(9): 625711 Xiang Xinghao, Chen Jianqiang, Yuan Xianxu, et al. Study on C-γ-Reθ model for hypersonic three-dimensional boundary layer transition prediction. Acta Aeronautica et Astronautica Sinica, 2021, 42(9): 625711 (in Chinese)
[30] 童福林, 李新亮, 唐志共. 激波与转捩边界层干扰非定常特性数值分析. 力学学报, 2017, 49(1): 93-104 (Tong Fulin, Li Xinliangy, Tang Zhigong. Numerical analysis of unsteady motion in shock wave/transitional boundary layer interaction. Chinese Journal of Theoretical and Applied Mechanics, 2017, 49(1): 93-104 (in Chinese) doi: 10.6052/0459-1879-16-224 [31] Jimenez J, Hoyas S, Simens MP, et al. Turbulent boundary layers and channels at moderate Reynolds numbers. Journal of Fluid Mechanics, 2021, 657: 335-360
-
期刊类型引用(1)
1. 刘春，乐天呈，施斌，朱遥. 颗粒离散元法工程应用的三大问题探讨. 岩石力学与工程学报. 2020(06): 1142-1152 . 
百度学术
其他类型引用(0)
计量
- 文章访问数: 1471
- HTML全文浏览量: 479
- PDF下载量: 161
- 被引次数: 1
