稀疏发现湍流边界层壁面函数公式

张 晖,叶 涛,王新光,万 钊

(1. 西南科技大学计算机科学与技术学院,四川绵阳621000;2.中国空气动力研究与发展中心计算空气动力研究所,四川绵阳621000)

1 引言

计算流体力学技术已经在飞行器设计、流体机械设计等领域广泛应用,能够实现复杂外形高超声速流动的数值模拟,成为风洞实验的重要补充。数值模拟技术若要准确地获得壁面处的物理信息分布,需要在壁面边界层中布置非常细密的网格,这样会极大地增加迭代收敛的计算步数,使迭代计算的稳定性下降,同时还会带来严重的数值刚性问题。工程上一般采用壁面函数的方法来获得壁面热流和速度等物理信息,壁面函数的使用在很大程度上放宽了近壁面网格密度,降低了对网格的依赖性,加速了收敛过程,在处理边界层计算问题时,壁面函数发挥着重要作用[1]。

边界层沿壁面垂直方向可以分成三个子层:粘性底层,过渡层,对数律层,不同层具有不同的流动状态。由于过渡层的复杂性,标准壁面函数通常设计为两层模型,分别给出了粘性底层和对数律层的壁面规律,并未对过渡层做出任何处理,常常会产生速度和温度的不连续。而统一壁面函数与实验值相差较大,无法准确地得到壁面附近的物理量变化情况,仅适用于理论研究中观察其趋势,并不适用于工程应用。

机器学习能够发现数据中存在的规律,近年来,数据驱动发现控制方程的机器学习方法引起了广泛的关注。由于计算量较小,在多种机器学习算法中,稀疏回归方法在发现各种问题的控制方程方面表现出巨大的潜力。使用稀疏回归可以从预先定义的大候选库中识别构成控制方程的少量项,得到一个简约模型。Brunton[2]提出了一个非线性动力学稀疏辨识框架来发现动态系统的控制方程。Rudy[3]在Brunton的基础上提出了顺序阈值岭回归算法,产生了一个非线性动力学的偏微分方程函数辨识框架。在这之后,大量的文献开始研究稀疏回归在控制方程数据驱动发现中的应用[4-13]。本工作之前的成果表明稀疏回归能够用于仿真数据的控制方程获取并从中得到壁面函数[14]。除了稀疏回归方法,其它技术也被用于执行数据驱动的控制方程的发现,如高斯过程和神经网络[15-17]。Maziar[18]等提出了一个利用高斯过程来发现控制方程的框架。Long[19]等在其之前的工作基础[20]上提出了一种数字和符号混合的深度神经网络,通过学习复杂系统的动力学来揭示隐藏的偏微分方程模型。这些方法方法不会对偏微分方程中的无关项进行惩罚,因此不能产生简约的模型。

本文基于弹性网络法提出了一种新的稀疏回归算法,称为顺序阈值弹性网络(Sequentially Threshold Elastic Net, STEN),并首次将稀疏回归方法用于壁面实验数据识别并得到了适用于工程的壁面函数。STEN结合了LASSO和岭回归的优势,在保证稳定性的同时也具有可快速计算的优点。将提出的STEN算法和传统的稀疏回归LASSO算法分别用于壁面实验数据拟合。计算结果表明,LASSO和STEN两种算法在壁面数据拟合中均表现出色,且STEN算法比LASSO算法相比结果更准确,效果更好。

2 壁面函数

湍流壁面函数是一种使用预定函数作为边界条件来解决近壁面复杂现象的方法。壁面函数允许在壁面附近使用相对粗糙的网格单元,从而节省计算时间。应用简化的雷诺平均Navier-Stokes方程,可以将不同变量(速度、温度等)表示为关于壁面距离的函数,并将边界层中的流动信息作为第一个节点的边界条件来进行数值模拟,这些函数称为壁面函数。

边界层不同层的分子粘性和湍流粘性作用不同,其中过渡层分子粘性和湍流粘性的作用相当,流动状态比较复杂,很难用一个公式或定律描述,且其厚度较小,工程上一般忽略过渡层的影响。所以,标准壁面函数常设计为两层模型。采用两层模型,易造成物理信息分布的不连续性,因此有学者提出了统一壁面函数。Crocco-Busemann提出了在对数层和粘性底层内统一有效的可压缩流动边界层壁面规律,其中温度统一壁面函数如式(1)所示[21-23]

T/Tw=1+βu+-Γ(u+)2

(1)

其中

u+为无量纲速度,μw为壁面粘性系数,qw为壁面热流,ρw为壁面密度,kw为壁面热传导系数,uτ为摩擦速度,Cp为压力系数,r为壁面律常数,Tw为壁面温度,T为温度。

但计算流体学界对统一壁面函数这一经验公式的正确性一直存在争议,本文从实验数据中获取方程以验证该方程的准确性。除了温度壁面函数,工程上还会用到速度壁面函数,本文主要研究温度壁面函数,有关速度壁面函数以及标准壁面函数等其它内容可以参考文献21到23。

3 稀疏辨识算法

本文采用稀疏辨识算法来处理回归问题。根据收集到的数据设计回归模型的函数候选库,并使用稀疏辨识方法从数据中发现方程。稀疏辨识很好地避免了过拟合,同时还能通过压缩系数来降低函数的复杂度。

假设自变量和因变量间存在如下规律

y=f(x)

(2)

其中f代表一个非线性函数,根据收集到的数据可以得到以下矩阵关系

(3)

构造一个n列的候选非线性函数的函数库Μ(x)如式(4)所示。Μ(x)可以由多种不同的函数项构成,候选库的设计根据领域知识由领域专家确定。

(4)

Μ(x)每一列代表式(2)右边的一个候选函数项,在许多物理系统中,这些非线性函数项只有少数在f中是活跃的,因此可以通过稀疏回归来确定系数η=[η1,η2,η3,…,ηm]的稀疏矢量

y=Μ(x)η

(5)

η中的每一列代表一个稀疏向量。其中η=[η1,η2,η3,…,ηm]是大小为m的系数向量,m是库Μ(x)中的候选项个数。将根据自变量数据x设计的候选函数库Μ(x)和收集到的y带入到(5)中,通过求解系数矩阵η来获得函数形式

(6)

若η中某一特定项的系数不为零,则其相对应的函数项在本模型中是活跃的。函数f通常只包含几个重要的项,使得它在可能的函数空间中表示出强稀疏性。

针对稀疏求解获得η,可采用最小二乘法,岭回归、LASSO算法和顺序阈值最小二乘法。最小二乘法会辨识所有函数项,使得所有系数均为非零值,可能会导致过拟合;LASSO算法压缩系数能力较强,但稳定性有所减弱;岭回归稳定性高但压缩系数的能力较低。而弹性网络算法结合了LASSO和岭回归两种算法的优点,既能够压缩系数也能够保证稳定性。本文基于弹性网络算法提出了一种新的稀疏求解算法,称为顺序阈值弹性网络法(Sequentially Threshold Elastic Net, STEN)来求解稀疏向量。该方法是以缩小变量集为思想的压缩估计方法,它通过设计一个惩罚函数来压缩变量的系数并使部分函数项回归系数为零,以此来进行变量选择。STEN的代价函数为

(7)

其中λ1和λ2分别对L1和L2惩罚函数的权重进行正则化。方程(7)中的λ1和λ2通过对惩罚项施加更多的权重来控制稀疏性的数量,从而导致系数的收缩和避免过拟合。可以使用一个正则化权重参数α和混合参数λ来修正方程(7)

(8)

利用λ参数控制L1和L2的凸组合,当λ=1时为LASSO算法,而λ=0为岭回归算法。正则化的权重α与惩罚项相乘,可以决定模型所需的稀疏强度,α=0代表一般的最小二乘法,α越大则稀疏性越强。为保证结果的稀疏性,对方程(8)结果添加一个阈值,对小于阈值的系数设为0,再将剩余系数进行递归计算直到非零系数收敛。该算法计算效率高,能够通过少量迭代快速收敛得到稀疏解。

4 壁面数据挖掘

4.1 实验准备

本文所用数据来自于壁面网格实验,整个网格实验在Ubuntu16.04系统下结合OpenFOAM5.0版本软件进行,对5种不同马赫数的平板边界层进行了数值模拟,详细的流动情况见表1。本工作针对高马赫数下壁面函数的挖掘,因此选取的均为马赫数大于5的算例,且表中算例是平板边界层中的经典算例[24]。马赫数的增加会增强湍流结构的变化,不同的马赫数算例可以得到多样的数据,这使得挖掘出的壁面函数更具有普适性。

表1 算例信息

Crocco-Busemann温度统一壁面函数与实验数据存在较大差异,其结果如图1所示。由图1可以看出,Crocco-Busemann函数并不适合用于工程应用获得温度信息。为了获得适用于工程的温度壁面函数,对温度实验数据进行数据挖掘。

图1 Crocco-Busemann函数与实验值对比图

4.2 实验数据挖掘

Crocco-Busemann温度壁面函数计算结果虽然与实验数据存在较大差异,但其整体趋势是符合实验结果的,同时相关理论研究也肯定了其价值。本文参考Crocco-Busemann函数形式设计了新的函数形式,如式(9)所示。在无量纲速度为零的时候,流体还未进行流动,此时壁面边界层的温度就是壁面本身的温度,壁面边界层不具有热量梯度不会产生热传递。因此在无量纲速度为零时,边界层的温度与壁面温度一致,基于此理由,将常数项设置为1,与理论公式一致。

T/Tw=1+B*βu++A*Γ(u+)2

(9)

结合式(9)的函数形式对壁面实验数据进行挖掘得到系数A和B,进而得到适用于工程的统一壁面函数。

实验过程中结合序列最小优化思想(Sequential minimal optimization, SMO)对STEN和LASSO算法进行优化,将初始问题不断地分解为子问题,再对子问题进行求解分析,最终使得全部变量都满足条件。根据SMO算法思想对方程(9)进行控制变量处理,控制A和B其中一个值为定值,对另外一个值进行回归计算。

首先控制修正系数A为定值,定值的选择上延续原函数的取值,原函数中A的值为-1,方程表示如下

T/Tw=1+B*βu+-Γ(u+)2

(10)

使用平方误差最小法计算得到各算例的B值如表2所示。

表2 控制修正系数A后各算例的B值

在流体力学中,各个物理量之间都存在可以用函数来表示的关系,而各个物理量都可以通过基本物理量来计算得到。因此选取壁面函数研究中的基本物理量作为自变量,分别是:马赫数Ma,壁面剪切应力τw,壁面密度ρw和壁面温度Tw。根据自变量设计了88项函数项来组成候选函数库,库中包含常数项、一次项、二次项、三次项、平方根项,三角函数项等等。函数侯选库设计完成后,使用STEN和LASSO算法进行回归分析,STEN算法在正则化权重参数α=0.09和混合参数λ=0.8时可以得到统一结果

(11)

把回归得到的B的统一形式带入到方程(10)中进行计算,并将计算结果与实验结果和理论函数进行比较,部分算例的对比图如图2所示,图中只展示了Case3算例的效果图,其余算例效果与此相似。由图2可知,两种稀疏算法在控制A情况下所得的统一函数效果均不好。其误差虽然较原理论函数有所降低,但与实验结果之间的误差仍然较大,并不能很好地匹配实验值。

图2 控制A得到的统一函数效果对比

控制变量B取原函数值,而对修正系数A进行计算来获取统一壁面函数,所得到地函数形式为

T/Tw=1+βu++A*Γ(u+)2

(12)

原函数中B值为常数1,选取与原函数一致的B值,能够得到如式(12)表示的函数表达式。使用平方误差最小法计算得到各算例的A值如表3所示

表3 控制修正系数B后各算例的A值

使用STEN和LASSO算法对表3的结果进行回归分析,得到修正系数A的统一形式:

(13)

将所得A值带入到(12)中进行计算并与实验数据进行对比。由图3可知,两种算法的结果均与实验数据匹配度较高,相对于原温度壁面函数效果均有所提升。图3仅展示了Case3算例的效果图,其余算例效果与此相似。

图3 控制B得到的统一函数效果对比

从所得结果可知,B取理论公式值,计算控制变量A所得结果对原壁面函数改进最为明显,与实验数据的匹配度较高。进一步对控制变量A所得结果进行误差分析,选取误差最小的作为最终回归结果并进行验证。

4.3 误差分析及结果验证

本文的工作是从壁面数据中抽取出一个统一的壁面函数,更关注壁面函数的计算结果与真实值数据之间的整体匹配度。因此对壁面实验数据使用不同稀疏算法得到的统一壁面函数进行误差分析,使用均方误差来衡量误差大小。误差计算结果如表4所示,虽然STEN和LASSO算法所得结果与实验值误差均较小,但相比之下,STEN算法误差更小,效果较LASSO算法更好。

表4 均方误差结果表

根据误差分析结果选取误差更小的STEN算法所得结果为最终结果。虽然式(13)的A表达式有表达出一定的物理意义,若0.0013这样的小数能用物理闭合系数来表示则能够更好的体现出其与物理量的相关性。其中常数-1.0021进行四舍五入为-1,对0.0013系数使用壁面常见闭合系数进行表示可以得到,常见的物理闭合系数有κ,B,Cp,r,Prl,Prt等。

(14)

将A的最终表达式(14)带入到(12)中得到最终确定的统一壁面函数形式为式(15)。

(15)

将该结果在各算例上进行验证,对各个算例的效果进行展示,如图4所示。

图4 所得函数与实验数据对比效果图

式(15)为本文根据原壁面函数形式结合SMO算法优化后STEN回归算法对实验数据,进行回归分析得到的最终函数形式。所得结果不仅能够较好地匹配实验数据,使得后期工程应用中关于壁面的处理成本降低,同时也能提高数值模拟的效率。

5 总结

本文提出了一种基于弹性网络的稀疏回归算法STEN,并将其对壁面实验数据进行数据挖掘得到控制方程。STEN与传统的LASSO稀疏回归算法相比,二者在壁面实验数据的挖掘中效果均较好。误差分析结果表明,STEN算法较LASSO算法误差更小,准确度更高,STEN算法对于壁面数据更具适用性。通过对实验数据进行拟合得到了统一壁面函数,其与实验数据匹配度高,对比理论壁面函数效果有较大提升。将本文研究所得的统一壁面函数应用于工程应用中可以大大缩短数值模拟的时间,提高工程效率。

在接下来的工作中将进一步分析挖掘出的统一壁面函数背后的物理机制,为该公式的工程应用提供更多的理论支撑。