江 昊 王伯福 卢志明

(上海大学,力学与工程科学学院,上海市应用数学和力学研究所,上海 200072)

引言

非线性动力系统广泛存在于自然界及工业界中,对非线性动力系统的建模是亟待解决的问题.非线性动力系统(如电力网络[1]、金融学[2]、神经科学[3]等)的控制方程大多是非线性偏微分方程,从传统的守恒定律和物理原理等第一性原理出发去推导一些复杂系统的控制方程是困难的.现代流体力学虽然已经具备了较完备的理论方程和模型,但仍然有很多复杂问题的理论研究尚不完备,如带化学反应流、多相流、非牛顿流、稀薄流等,这些系统的动力学方程或模型的理论推导难以实现[4].近些年来,随着计算机硬件以及计算科学的发展,基于数据驱动方法对非线性复杂系统建模得到了研究者的青睐[5-13].其中稀疏识别方法[14]可以简化复杂的数学模型,得到项数较少的简单模型来表征系统的动力学特性,并可以综合考虑过拟合和模型精确性之间的矛盾.稀疏识别方法已经被广泛用于复杂非线性系统控制方程的识别[14-21].1996 年,Tibshirani[14]首次提出最小绝对收缩和选择算子(least absolute shrinkage and selection operator,LASSO)方法,并且把LASSO 方法应用在Cox模型的稀疏识别中[15].2016 年,Brunton 等[18]通过稀疏识别在含有噪声的非线性动力系统的数据中识别出了常微分方程.2017 年,Rudy 等[19]提出了偏微分方程函数识别(partial differential equations functional identification of nonlinear dynamics,PDE-FIND)方法可以通过测量或者数值模拟的时空数据来稀疏识别偏微分方程,正确地识别出了圆柱绕流流场的涡量输运方程.胡军等[20]将贝叶斯稀疏识别方法识别结果与PDE-FIND 方法识别结果相比较,表明了贝叶斯稀疏识别方法对偏微分方程具有非常强的稀疏恢复能力,但相比PDE-FIND 方法对噪声更加敏感.

近年来,数据驱动方法在流场中的应用颇受研究者[4,22-27]的重视.2019 年,Chang 和Zhang[23]使用LASSO 方法从数据出发,识别出了单相地下水流动方程和污染物输运方程.2020 年,Raissi 等[24]使用深度学习,对流场可视化的图片进行学习,并反演了流场的压力和速度分布场,可以帮助得到实验难以测量的数据.Zhang 和Ma[25]采用PDE-FIND 方法,对基于直接模拟蒙特卡罗(direct simulation Monte Carlo,DSMC) 方法模拟产生的数据进行稀疏识别,反演了流体的对流、扩散以及涡量输运方程,从数据驱动角度清晰地展示了玻尔兹曼方程和Navier-Stokes 方程之间的相关性.Schmelzer 等[26]采用SINDy 框架在数值模拟的数据中是识别出了代数雷诺应力模型.张亦知等[27]构建了基于数据驱动的湍流模型修正框架并在槽道流中进行了验证.

本文采用基于数据驱动的稀疏识别方法(PDEFIND 方法和LASSO 方法)对圆柱绕流、顶盖驱动方腔流、RB 对流和三维槽道湍流的控制方程进行稀疏识别,研究了不同流场,不同控制方程以及不同区域数据对于流场控制方程稀疏识别的影响.

1 流动控制方程的识别方法

非线性系统的控制方程大多是非线性偏微分方程(组),例如Burgers 方程(1),科尔特弗-德弗里斯(korteweg-de Vries,kdV)方程(2)

这些非线性偏微分方程可以写成通用形式

其中,u表示时空变量,ut表示变量u的时间导数,下标x的个数表示变量u的空间n阶导数,μ 表示系统的参数,N(·)表示时空变量u(x,t)及其时空各阶导数的非线性函数.在非线性偏微分方程的稀疏识别中,将包含变量u及其导数相关的项组成矩阵U,将其他可能出现的项组成矩阵Q.将U,Q合并成矩阵Θ

对于离散的数据点,U和Q分别为(N×Mu) 和的二维矩阵,其中N=n×m为矩阵中每项对应的时空点数(n为空间网格总数,m为时间步数),Mu和Mq分别为U和Q的项数.将非线性偏微分方程表示为以下线性方程

其中,Θ是人为给定的一个过完备的候选项的库,即Θ库中包含所有可能属于该非线性偏微分方程的项.ξ是稀疏向量,表示的是非线性偏微分方程系数矩阵,该稀疏向量ξ的非稀疏项等于该非线性偏微分方程的系数,则偏微分方程可以写成线性加权和的形式

流体力学中非线性偏微分方程的稀疏识别就是通过流场的部分空间位置的时间序列数据,来找到稀疏系数矩阵ξ使得方程(6)满足该流场动力学行为.时间序列数据可以通过实验或者数值模拟来获得,变量的时空导数通过多项式插值法或有限差分法获得.对于不含有噪声的数据通常采用有限差分法计算获得,但对于有噪声的数据,多项式插值法计算的导数比有限差分法计算的导数更精确[19].

如果Θ内所有的项都属于该流场系统的控制方程,此时可以使用最小二乘法确定系数矩阵

然而,在过完备的库Θ中通过最小二乘法求解系数矩阵ξ,会引起过拟合问题,解出的系数矩阵ξ中所有元素都是非零,这样所识别出的非线性偏微分方程很复杂,比流场精确控制方程多出很多项,不符合流体动力学规律.

对于偏微分方程的稀疏识别,一种比较直观的方法就是将所有项可能的组合都进行计算,比较所有项组合的误差,得到最优的项组合,数学表示为加入L0正则化最小二乘法

但是求解该问题会遇到NP-hard 问题,即把所有的组合都进行计算和比较,这对于计算资源的消耗是巨大的.Tibshirani[14]提出的LASSO 算法是引入L1正则化的最小二乘法

该方法是对于L0正则化最小二乘法的一个凸优化问题,避免了L0正则化的NP-hard 问题,图1 是LASSO方法的流程图.但LASSO 方法对列相关性很高的数据矩阵的稀疏识别比较困难,例如对于KdV 方程的识别[19].

图1 LASSO 方法流程图Fig.1 Flow chart of LASSO scheme

Rudy 等[19]提出了一个新的稀疏识别方法:PDEFIND 方法.该算法是在加入L2正则化的最小二乘法(1) 中加入阈值筛选,使用Ridge 算法[28]求解L2正则化的最小二乘法,将小于阈值的系数所对应的项筛去,对剩余的项重复递归Ridge 算法求解和阈值筛选,直到所有项的系数都高于阈值,此时返回稀疏系数矩阵ξ,该过程被称为STRidge 算法

显然,阈值的选取对于STRidge 算法很重要,不同的阈值会导致识别的偏微分方程的稀疏程度不同,因此为了找到最佳的阈值,Rudy 等[19]提出了一个训练阈值的算法(TrainSTRidge 算法).TrainSTRidge 算法和STRidge 算法相互耦合一起形成了PDE-FIND方法.图2 是PDE-FIND 方法的流程图.

图2 PDE-FIND 方法流程图Fig.2 Flow chart of PDE-FIND scheme

2 识别结果

2.1 圆柱绕流涡量输运方程的识别

Rudy 等[19]使用PDE-FIND 方法对圆柱绕流的涡量输运方程进行了稀疏识别,在未添加噪声的流场数据中识别出的涡量输运方程的平均系数误差为1%±0.2%.本文同样也对圆柱绕流问题的涡量输运方程进行了识别.

采用LBM 方法对雷诺数为300 的圆柱绕流进行了直接数值模拟,计算网格为1000×250,圆柱圆心位置为(201,133),半径为26,来流最大速度为u=0.1.图3 给出了t=104时流场的涡量云图.如图3 所示将圆柱绕流后面的区域分为x∈ (250,500),x∈ (500,750),每个区域随机选择8000 个空间点并追踪记录60 个时间步长的速度数据和压力数据,这样每个区域就共获得了48 万个数据样本,这对于稀疏识别所需的样本是充足的,且采样的方法不会影响识别的结果[19].采用5 次切比雪夫多项式插值[29]计算导数,导数阶数最高保留到二阶,变量次数最高保留到二阶,非线性项最高保留到四阶,将这些数据用来构建稀疏识别所需Θ库(4).表1 给出了PDE-FIND 方法和LASSO方法对于x∈(250,500) 圆柱绕流数据的稀疏识别结果.

图3 t=104 时圆柱绕流的涡量云图Fig.3 Vorticity contour of flow past a cylinder at t=104

表1 Re=300,x ∈(250,500)圆柱绕流控制方程的识别Table 1 Identification of vorticity transport equation for flow past a cylinder with Re=300 in x ∈(250,500)

从表1 可得出,PDE-FIND 方法和LASSO 方法都正确地识别出了涡量输运方程,并且系数相同.两种算法对于涡量输运方程的稀疏识别的能力一致,且平均系数误差为3.21%.涡量输运方程是一个线性偏微分方程,PDE-FIND 方法和LASSO 方法对线性偏微分方程的稀疏识别能力一致.对于x∈(500,750) 区域的数据,两种算法也都识别出了正确的控制方程,但所得平均系数误差有所增加,为8.75%.这是因为越靠近圆柱区域的流场信息越丰富,因而识别出的控制方程系数越精确.这与文献[20]的结论是吻合的,即选取具有丰富动力学行为的时空演化数据可以提高对动力学控制方程的识别的准确性.

2.2 顶盖驱动方腔流控制方程的识别

本节对二维顶盖驱动方腔流的NS 方程组进行稀疏识别.相比于上一节涡量输运方程的识别,原始变量形式的NS 方程组是一组偏微分方程组

其中,方程(12)是连续性方程,方程(13)是动量方程.二维顶盖驱动方腔流场数据使用基于有限差分的直接数值模拟得到.方腔宽高比为2,方腔顶部的壁面以大小为1 的速度匀速向左运动,流动雷诺数为100.模拟计算网格数为200×100,图4 给出了t=40 时流场的压力云图和流线图.

图4 t=40 时顶盖驱动方腔流的压力云图和速度矢量图Fig.4 Pressure contour and velocity vector of lid-driven cavity flow at t=40

注意到连续性方程(12),不含有时间偏导项,将连续性方程可改写为

将x看成时间t,就可以采用稀疏识别方法来识别连续性方程了.

在整个方腔区域选取数据,采样方法、导数计算方法以及Θ库的构建与2.1 节相同.表2 给出了PDE-FIND 方法和LASSO 方法对于Re=100 顶盖驱动方腔流控制方程的稀疏识别结果.从表2 可得出,对于顶盖驱动方腔流的数据,PDE-FIND 方法正确地识别出了连续性方程和NS 方程的各项,且平均系数误差为2.77%,而LASSO 方法只能正确地识别出连续性方程,没有识别出动量方程,例如在u方程中未选择出uux,可能是因为LASSO 算法对含强非线性项的方程的识别会失效,具体的原因将在下一节进行讨论.

表2 Re=100,顶盖驱动方腔流控制方程的识别Table 2 Identification of vorticity transport equation for lid-driven cavity flow with Re=100

2.3 Rayleigh-B´enard 对流控制方程的识别

下面考虑识别Rayleigh-B´enard(RB)对流系统中的控制方程,其物理模型是一个下底板加热,上底板冷却的对流系统,系统中充满流体介质,由于上下底板温度不一样,导致不同高度位置的流体之间存在密度差,在浮力的驱动下系统中流体产生运动.在OberbeckBoussinesq (OB) 近似下的无量纲控制方程组如下

其中,方程(15)是连续性方程,方程(16)是含有热浮力项的NS 方程,方程(17)是热输运方程,Rayleigh 数(Ra)和Prandtl 数(Pr)是RB 系统的重要控制参数.

RB 对流流场通过Zhang[30]所发展的基于有限差分的直接数值模拟得到.不失一般性,模拟方形区域的RB 对流,Pr数固定为1,计算了4 个不同的Ra数的工况.Ra数分别取106,107,108,109,对应的计算网格分别是128×128,512×512,720×720,1024×1024.当Ra数从106增大到109时,RB 对流系统会从层流状态转变为强非线性的湍流流动状态,整个系统的速度及温度脉动剧烈,会出现局部梯度较大的数据,对于系统控制方程的识别可能造成一定的困难.图5 是t=80 时,不同Ra数的RB 对流系统的温度云图和速度矢量图.

图5 Pr=1,不同Ra 数的RB 对流系统的温度云图和速度矢量图Fig.5 Temperature contour and velocity vector of RB convection for different Rayleigh number with Pr=1

在整个计算区域选取数据,采样方法、导数计算方法以及Θ库的构建与2.1 节相同.

表3 给出了Ra数为109的识别结果,PDE-FIND方法和LASSO 方法都正确识别出了Ra=109的RB对流系统中连续性方程和热输运方程,并且系数一致.而且PDE-FIND 方法仍然正确识别出了RB 对流的动量方程,对应Ra从106到109,识别的控制方程组平均系数误差分别为0.36%,0.06%,0.40%,0.46%.而LASSO 方法对NS 方程仍没有正确识别,在u方程识别中未选择uxx.原因是由变量及变量各阶导数组成的过于完备库中,回归系数趋于相近的项之间会产生分组效应(grouping effect)[31,32].而LASSO 方法对于具有分组效应的候选项的筛选可能会失败,表现为在有分组效应的候选项中只选择一项,而其他的候选项的系数会设为0.在Ra=106,Ra=107,Ra=108的RB 对流控制方程稀疏识别中,也存在同样的现象.

表3 Ra=109,Pr=1RB 对流控制方程识别Table 3 Identification of governing equations for Rayleigh-B´enard convection with Ra=109,Pr=1

2.4 三维槽道湍流控制方程的识别

最后对三维槽道湍流进行了稀疏识别,数据通过离散统一的气动格式(discrete unified gas-kinetic scheme,DUGKS)[33]模拟得到,计算域为4π×2×4π/3,槽道半高为1,计算网格为128×128×128,垂向采用非均匀网格,流向和展向采用均匀网格,壁面摩擦速度定义的雷诺数Reτ=180,此时的槽道流处于充分发展湍流状态.图6 给出了槽道展向中心截面的速度云图,各种尺度的复杂结构清晰可见.在整个槽道中选取数据,采样方法、导数计算方法以及Θ 库的构建与2.1 节相同.

图6 Reτ=180,三维槽道湍流展向中心截面速度云图Fig.6 Velocity contour of 3D turbulent channel flow with Reτ=180 at z=0

表4 给出了三维槽道湍流控制方程的识别结果,PDE-FIND 方法对于正确地识别出了控制方程,并且平均系数误差为0.37%,说明PDE-FIND 可以在在湍流状态的流场数据中识别控制方程.而LASSO方法只识别出了连续性方程,这与前两节的结果一致.说明在过完备候选库中列相关性强弱对于PDEFIND 方法的稀疏识别没有影响,PDE-FIND 方法不存在分组效应.本节还研究了不同y+的数据对识别的影响,取y+∈(0.21,19.3),y+∈(19.3,59.7),y+∈(59.7,182.9),都正确识别出了控制方程,系数平均误差分别为0.42%,0.48%,0.39%.说明在三维槽道湍流控制方程的稀疏识别中,误差与壁面距离几乎无关.

表4 Reτ=180 三维槽道湍流控制方程识别结果Table 4 Identification of governing equations for 3D turbulent channel flow with Reτ=180

3 结论

本文采用PDE-FIND 方法及LASSO 方法对二维圆柱绕流、顶盖驱动方腔流动、RB 对流和三维槽道湍流控制方程进行了稀疏识别,主要得到以下结论:

(1)使用PDE-FIND 方法识别出了流场中的控制方程,得到了正确的方程及精确的系数.

(2) LASSO 方法正确识别出了其中的线性偏微分控制方程如连续性方程,涡量输运和温度输运方程,但对含非线性的偏微分方程识别效果差,这是由于LASSO 方法识别过程中侯选库存在分组效应,降低了识别能力.而在PDE-FIND 方法中是不存在分组效应的,所以PDE-FIND 方法对强非线性方程识别有一定的优势.

(3)丰富的流动结构对于数据驱动的流场控制方程的稀疏识别有重要作用,选择具有丰富流动结构的区域的数据可以提升稀疏识别的有效性和精确性.

从2.3 节和2.4 节结果可知,PDE-FIND 方法对层流和湍流两种流动状态都可以精确的识别出控制方程.这对于复杂流动建模提供了一种新的研究方法,例如将PDE-FIND 方法用于湍流脉动量控制方程和湍流统计量的控制方程的系数识别.未来考虑将稀疏识别算法与k-ε 等湍流模式相结合,提高湍流模式在数值模拟中的精确性.另外,本文的流场数据都是通过精细的直接数值模拟得到的,PDE-FIND 都识别出了正确的方程形式,得到了准确的方程系数,但从实验或者实测数据往往具有噪声,对于含有噪声的数据,PDE-FIND 方法识别能力往往变弱,如何提高PDE-FIND 方法在含有噪声数据的识别能力是值得深入研究的课题.可行的方法之一是结合数据去噪技术和寻找最优的插值点数和多项式插值阶数等,提高高阶导数项的计算精度,从而提高稀疏识别能力.