面向稀薄流非线性本构预测的机器学习方法

1、面向稀薄流非线性本构预测的机器学习方法摘要：稀薄非平衡流域内连续介质假设已经失效，主要围绕Boltzmann方程及模型方程对稀薄非平衡流开展理论与 计算研究，统一气体动理论格式(UGKS)是其中一种代表性方法在稀薄非平衡流数值模拟中，Navier-Stokes (N-S)方 程连续介质假设已经失效，不能有效描述流场非平衡特征 UGKS方法虽然计算精度高，但速度空间离散导致计算效率 低下，多维高速条件下数值计算难以开展基于数据驱动的思想，在N-S方程与UGKS方法的研究基础上发展出了一 种稀薄非平衡流非线性本构关系求解方法(DNCR)%该方法以N-S与UGKS求解器获得的流场数值模拟计算结果作为

2、 训练数据集，基于流场特征参数采用极端随机树算法生成机器学习模型，对预测流场中线性黏性应力项与热流项进行非 线性修正，并耦合非线性本构关系求解宏观守恒方程得到目标状态稀薄非平衡流动数值解针对DNCR方法中所采用 的机器学习方法-极端随机树模型，通过二维顶盖驱动方腔流算例对高维非线性建模涉及的特征参数选取、参数调优开 展了相关验证工作，选取若干典型状态对极端随机树模型的泛化性能开展研究，并评估了相关模型与方法的计算精度与 计算效率%关键词：稀薄非平衡流；本构关系；数据驱动；极端随机树；模型参数钱学森根据努森数Kn定义流体的稀薄程 度，将流动区域划分为连续流域(Kn 0.01 )、 滑移流域(0.

3、 01 Kn 0. 1 )、过渡流域(0.1 Kn 10 )，稀薄非 平衡流域包含除连续流域以外的3大流域(1) %连续流域符合流体力学连续介质假设条件， 连续介质流体力学是其理论基础，其中最具有代 表性的控制方程为Navier-Stokes (N-S)方程，形 成了以广义牛顿定律和傅里叶热传导定律为基础 的基于 NSF (Navier-Stokes-Fourier)关系的连续 流数值模拟计算方法，此方法的发展为其他各种 方法的建立提供了宝贵经验稀薄非平衡流域中 气体稀薄属性逐渐凸显，例如物面出现了十分显 著的速度滑移与温度跳跃现象此时在物面附近 的努森层内，连续介质假设已经失效，准确地说是

4、基于3大守恒方程(质量守恒、动量守恒和能量守 恒)推导出来的黏性应力项与热流项已经不能再 简单使用低阶宏观物理量(速度、温度)的梯度线 性表征，即线性本构关系已经不再适用于精准描 述稀薄非平衡流动问题对稀薄非平衡流问题的研究主要围绕稀薄气 体流动控制方程Boltzmann方程开展，它是 分子气体动力学的基本方程，可以不受努森数的 限制对连续流到自由分子流整个流域进行统一描述。Boltzmann方程是一个复杂的七维积分/微 分方程，大部分研究均是对其直接或者间接求解, 发展形成了多种数值计算方法与理论，统一气体 动理论格式(Unified Gas-Kinetic Scheme, UGKS)是其中

5、一种代表性方法。UGKS方法 使用 Boltzman 方程的 Bhatnagar!Gross-Krook (BGK)&-类模型方程作为控制方程，用有限个离 散速度替代整个速度空间，在一定条件下均能给 出较为准确的流动特征，具备求解各流域多尺度 问题的能力30 UGKS方法采用BGK类模型方 程基于当地积分解对通量进行构造，将分子运动 与碰撞过程自然地耦合在一起，物理意义更加明 确，突破了 DSMC方法由解耦处理所带来的计算 时间步长小于分子的平均碰撞时间、计算单元网 格尺度小于分子平均自由程的限制。但由于 UGKS方法在求解过程中依赖速度分布函数对 速度空间进行离散，这对计算与存储的要求非常

6、高，多维与高速条件下计算效率较低！随着计算机理论与技术的迅猛发展，大数据 时代启发了人类思考问题新的思维，基于数据驱 动和机器学习的研究方式也应运而生！研究人员 开始通过机器学习的方法解决流体力学领域中难 以解决的问题，目前在气动优化设计&沁、湍流问 题项、非定常气动力建模16-20等领域上已经开 展了很多卓有成效的研究工作！以N-S方程和UGKS方法为理论基础，提出 了一种适用于稀薄非平衡流数值模拟的基于数据 驱动非线性本构关系(Data-driven Nonlinear Constitutive Relations, DNCR)的数值计算方 法。基于流场特征参数采用机器学习方法对N-S 方

7、程线性的黏性应力项与热流项进行非线性离散 重构，理论适用范围可以覆盖较大来流努森数条 件，通过耦合非线性本构关系求解宏观守恒方程 得到待预测状态稀薄非平衡流动数值解！与传统 机器学习建模方法摒弃基本物理规律直接对数据 进行训练与预测的思想相比，DNCR方法未抛弃 物理规律而是对流体本构关系进行修正后耦合守 恒方程迭代求解，同时保留了数据建模与物理建 模的优点，方法物理意义更加清晰明确。在DNCR方法中，采用的具体机器学习算法 为极端随机树zl-zz ( Extremely Randomized Trees)算法。在机器学习算法中，针对所研究的 具体问题选取低冗余、高代表性的特征参数对于 机器学

8、习算法的泛化性能和运行效率具有重要影 响。本文拟通过二维顶盖驱动方腔流算例对高维 非线性建模涉及的特征参数选取、参数调优开展 相关验证与研究工作，选取若干典型状态对极端 随机树模型的泛化性能开展研究，并评估相关模 型与方法的计算精度与计算效率！1数据驱动非线性本构关系计算方法DNCR基于N-S方程和UGKS方法的数值 模拟计算结果作为流场样本训练数据，采用机器 学习方法构建热流/应力项与流场特征参数的高 维复杂非线性回归关系模型，对N-S方程本构关 系进行非线性修正。通过耦合数据驱动的非线性 本构关系求解宏观量守恒方程得到稀薄非平衡流 数值解！DNCR方法计算流程如图1所示，欧代表在 N-S流

9、场数值模拟结果中提取的特征参数，作为 机器学习模型的输入量，为标记参数，作为机 器学习模型的输出量，即机器学习模型建立了 欧 与之间的复杂高维回归关系，与非线性本构 物理建模函数的区别在于其本构关系没有具体数 学表达式，而是全流域存在离散当地映射=? ,所形成的离散本构映射关系在全流场计算 域适用！而模型泛化性能取决于训练集特征代表 性与特征参数选取,理论上若特征参数选择未涉 及到任何与空间网格与当地坐标值直接相关的参 数且训练集包含相似的非平衡流动特征，机器学 习模型就具备一定的迁移预测能力！训练输入:特征参数预测输入:喘参数耦合非线性本构关系求解N-S方程得训练输入:特征参数预测输入:喘参

10、数预测输出:标记参数机器学习模型图1 预测输出:标记参数机器学习模型Fig. 1 Schematic of DNCRDNCR方法的一个重要特点是训练与预测 过程相对独立，图1中红线与绿线分别表示了训 练与预示流程。模型训练过程通过对包含不同典 型流动特征的基础流场数据集开展训练，获得欧 与(&之间的复杂回归关系。而预测过程则首先 采用N-S求解器对待预测状态开展初步计算，获 得待预测当地特征值=，然后采用已训练生成的 回归关系=? (&得到待预测流场当地标记值 ,最终通过时间推进方式求解质量、动量与能 量守恒方程，耦合非线性本构关系完成计算，守恒 方程为子气体Ar作为模拟气体，其动力黏性系数(

11、采用逆幂律公式)24*计算：(3)单原子气体Ar物性参数取= 2.272 X10%Pa( s,O0 = 300 K,+ = 1. 667,Fr = 0. 667,E=208. 16 J/(kg . K),其中(o 为 Ar 在温度 O# 下 的动力黏性系数，O0为Sutherland温度常量，+ 为比热比，Fr为普朗特数，$为气体常数，逆幂律 幂次取e = 0. 72。计算状态如表1所示。对于一维激波结构，如前所述将UGKS数值)Q +)E,)F,)El ( )3, _ B)8)9)8)9式中：；为守恒量E、F为无黏通量E, #3,为黏性 通量。控制方程中黏性应力张量与热流项表达式为(1)模拟

12、结果中的=8、，8、磬、萼数据提取出来作为/ 4)u I,8，(61 ( ( 8 4(61(9 4(61)+/ 4)u)9DNCR计算程序的读入数据。计算收敛后对比 DNCR、N-S、UGKS密度速度U、压强F、温度 O等基本物理量，如图2图5所示。可以看出 DNCR的计算结果与N-S相比更加接近UGKS 的计算结果，表明了本文所提出的本构关系非线 性修正方案的正确性与可行性。,89)u)8参数,89)u)8参数数值来流马赫数Ma8. 0激波上游静压/Pa73 477激波上游静温/K288. 15表1 一维Ar激波结构计算状态Table 1 Calculation state of one-d

13、imensionalArgon shock structure图2激波结构密度分布Fig. 2 Density distributions in shock wave,98 = ,89,(61 ( (9 4(61 + )8 4 (61 ) +/)u . )k()9 )8_ /4 k |，99 (61 . 9 4(61 TOC o 1-5 h z / 4)v)0)0 =8 = 8(61 -()8 4(61 )+()8)()0)0=9 = =9心-()9 l(6l)+()9)(2)式中：，为黏性应力项(为黏性系数；O为温度； U、K为速度； 为热传导系数=为热流项；下标 (ag表示待预测流场当地标

14、记值。式(1)与式(2) 随着时间推进，相关热流、应力张量与流场梯度量 向标记值趋近，最终完成计算收敛，获得待预测流 场定常解。当物体或扰动源的速度大于扰动信息的传播 时，扰动无法向上游传播所形成的强压缩间断称 为激波。一维激波结构是最简单、最基本的非平 衡流动现象之一，是研究本构关系与非平衡流动 的典型算例，可以用来对模型进行验证230N-S、UGKS、DNCR相同网格下，采用单原7SEM图3激波结构速度分布7SEMFig. 3 Velocity distributions in shock wave图4激波结构压强分布Fig. 4 Pressure distributions in sho

15、ck wave图5激波结构温度分布Fig. 5 Temperature distributions in shock wave2训练数据集生成本文选取图6所示顶盖驱动方腔流动作为测 试算例，计算网格如图7所示，计算网格为61X 61均匀网格，速度离散为28X28O选取5组稀薄 非平衡流动状态生成训练数据集，计算状态如表 2所示。图6顶盖驱动方腔流示意图Fig. 6 Schematic of lid-driven cavity flow1.00.60.40.200.20.40.60.81.0 x图7二维顶盖驱动方腔流计算网格1.00.60.40.200.20.40.60.81.0 xFig. 7

16、 Grid of two!d imensiona l lid-driven cavity flow表2二维方腔计算状态参数数值方腔尺寸/m参数数值方腔尺寸/m1上平板移动速度/(m s-1)50努森数Kn0. 5/0. 7/1. 0/1. 3/1/5来流静温/K273壁面温度/K2733机器学习模型在二维顶盖驱动方腔流问题中，标记参数(&共包括: TOC o 1-5 h z .)u2)T . )u2)(&共包括:(=! (,3 ! ! ! ! !)8)9)9其中23为张量指示符#本文回归问题建模使用 并行集成学习算法的典型代表$极端随机树! 该算法的基学习器采用分类与回归树Z5 (Classi

17、fication and Regression Tree, CART) #3.1 CART与极端随机树算法x与y分别为特征输入和标记输出变量，并 且丫为连续变量，给定训练数据集：d，)(8j, 9) !(8- , 9-) !！(85 , 95 ) + #回归树是用于研究回归问题的决策树模型， 决策树是一种树形结构的基本分类与回归算法# 在机器学习中决策树是一个预测模型!代表的是 对象属性(即特征参数)与对象值(即标记值)之间 的一种映射关系# 一个回归树对应着输入空间 (即特征空间)的一个划分以及在划分的单元上的 输出值。假设将训练数据集所在的输入空间划分 为M个单元:R,R-,*,Rm，这里

18、的单元可以理 解为从最顶端往下依次划分构建的M个分支,每 个分支只对应一个输出值#递归地将每个区域划 分为两个子区域并决定每个子区域的输出值c , 构建二叉决策树#1)选择最优切分变量3与切分点s :特征空 间包含多个特征参数即特征变量，决策树的构建 就是每层进行二叉树的划分，划分时采用某特征 变量能够使误差最小，此变量即为最优切分变量, 切分点指的是切分变量进行左右划分的值，小于 此值划分为左子树，大于此值划分为右子树，切分 时某切分点能够使误差最小即为最优切分点#通 过求解：min min , (92 C )- (C1 8+ $13,$)min , (92 c-)-(4)c-8 + $2

19、(3 ,s)遍历变量3 !对固定切分变量3扫描切分点 $ !选择使式(4)最小的3 ,s) #-)用选定的(j ,$)划分区域并决定相应的输 出值:$1 (j $) =)8 | 8() $)$- (j ,$) =)8 | 833 * $)C+ = 1, 92,8 + $+ , + = 1,2(5)5+ 8+Rm(j ,$)继续对两个子区域调用步骤1)、2) !直至 满足条件#将输入空间划分成M个区域R1R-,Rm !生成决策树：Mf(8)=U , 8 + Rm(6),1极端随机树在生成CART过程中选择划分 属性不再是选择最优属性而是完全随机选择#采 用训练数据集进行并行训练得到 #选取低冗余

20、、高代表性的特征参数对机 器学习算法的泛化性能和运行效率具有重要影 响#首先从物理机理上考虑,选择表征流场稀薄 非平衡特征的参数,例如基于当地流场梯度的局 部努森数k)gll5,P,T)作为主要的流场特征参 数!然后尝试加入其他参数并最终选取适当参数 集作为流场特征参数建立特征空间#本文采用方差过滤准则进行特征选择!方差 越小!表示该特征的数据差异越小，可以认为这个 特征对区分样本贡献不大!导致预测的效果越差! 因此可以剔除此特征或降低特征权重#另外考虑 到真实物理问题中不同物理量之间数据的量级可 能存在较大差异!量级较大导致方差较大!为了避 免这一不利因素!使用标准差与极差的商值消除 量级影

21、响进行特征选择#max(X) max(X) min(X)式中：X , )81,8-,，8)代表某一特征参数；) 代表样本点的数量；(为X , )81,8-,，8)的 算数平均值#标准差与极差做商的目的在于消除数据量纲的影响！采用二维顶盖驱动方腔流的数值模拟结果进行特征选择工作！流场中初始待选特征参数包括p、U、K、P、T、)p )P )P )8 p、U、K、P、T、)p )P )P )8 )x )y )9、K)5 Kn*、K)O方向上各特征参数的标准差情况如图8所示，特首先从稀薄非平衡特征上考虑，确定了基于当地流场梯度的局部努森数KnGd(p，P，T)作为 主要流场特征参数。由图8可知F2、F

22、3、F5即u、 k、O标准差较大，第1组特征参数选取u、k、首先从稀薄非平衡特征上考虑，确定了基于Kn5、Kn*、Kno，为了对比不同参数选取对预测 结果的影响，选取 5、u、k、*、O、|5、|*、|5、|*、 )8 )8 )y )y所选特征参数对模型泛化性能的影响，数值越小Kn5、Kn*、Kno作为第-组特征参数。本文采用 余弦相似度-6的另一形式为评价指标用于评估 所选特征参数对模型泛化性能的影响，数值越小Cos_similarity = 1 一手些严(9)式中：Ypredct代表预测的标记值，Yactual代表真实的 标记值，均可看作形如6 = *yi，y-，9n,的向 量，y,代表某

23、网格点标记值，n代表计算网格点 的数量。采用Kn = 0.7与Kn = 1. 3两组状态对机器 学习模型进行训练，对Kn = 1.。状态进行预测并 与Kn = 1.0的UGKS数据进行对比，计算各预 测量余弦相似度，第1组特征值(6特征值)与第2 组特征值(12特征值)结果对比如图9所示。由 图9可以看出，虽然训练集中特征参数U、K、O、图8特征参数标准差对比Fig. 8 Comparison of standard deviation ofcharacteristic parameters图9余弦相似度对比Fg.9 Comparson of cos ne s mlartyKnp、Kn*、Kn

24、o的数据差异更大，但极端随机树方 法在DNCR中采用12特征值时泛化性能更优， 在本文后续工作中特征参数将选用第图9余弦相似度对比Fg.9 Comparson of cos ne s mlarty数:P、U、K、P、T、5-、)8)8Kn*、Knr !)5 )5 )y在研究具体问题时，不同机器学习算法中都 有许多的参数需要人为设定，即使是同一种算法， 面向不同回归问题的参数配置也有所区别，最终 得到的模型性能表现往往也存在十分明显的 差别!对于极端随机树算法来说，主要参数包括框 架参数：基学习器的数量n_estimators、度量分裂 的标准criterion、是否采用袋外样本来评估模型 的好

25、坏oob_score；决策树参数：决策树的每个节 点随机选择的最大特征数max_features、叶子结 点上应有的最少样例数min_sampleseaf、决策 树最大深度max _depth等参数0本文重点对 n_est imators与 max .features两个参数进行研 究，其余参数均采取默认设置，研究方法为单一变 量法!本文采用可决系数&27评价参数调节后的模 型优劣情况!可决系数表示一个随机变量与多个 随机变量关系的数字特征，即可决系数值越大越 接近于1时说明特征参数对标记参数的解释程度 越高!res$ -tot(10)式中：表示可决系数; SSre?与SStot分别表示 残差

26、平方和与总平方和；9 + Y，D为Y均值；9 (10)采用Kn，0.7与Kn，1.3两组状态对机器 学习模型进行训练，对Kn = 1.0状态进行预测。 选取 n _estimators 在1,1 000(范围内,R- 随 n_estimators变化趋势如图10所示&可以看出 R-最终趋近于一个比较稳定的值，拟合较好且随 着基学习器的数量增加未出现严重过拟合现象& 如图11所示，模型的训练时间与基学习器数目基 本满足线性关系，考虑到模型训练的时间成本，本 文中选取 n_estimators = 300&图10R-随基学习器数目变化曲线Fig. 10 R- variation curve wit

27、h n_estimators图11训练时间随基学习器数目变化曲线Fig. 11 Trani time variation curve with n_estimators根据前文中确定使用的12个特征值,max_ features在1,12(范围中取值&为避免偶然误 差，训练1 000次后取平均值，得到R2值随maxfeatures 变化趋势如图12所示&可以看出在 DNCR方法中极端随机树的预测准确性随着 max_features变化大致呈现递增的趋势,max features的大小对模型训练预测时间的影响可以 忽略不计，在本文中选择max_features = 12&其 余模型参数及其他待预

28、测物理量的模型参数调优 过程与之类似&图12 R2随最大特征数变化曲线Fig. 12 R2 var i at i on curve wi th max_features.DNCR计算精度与极端随机树泛化性 能分析采用Kn = 0. 7与Kn = 1. 3两组状态对机器 学习模型进行训练，选取一组中间状态Kn = 1.0A )T .A )T 上一 进仃预测，(=2、) 、 一、 标)8)8记值预测结果如图13所示&为更好地反映预测值与真实值的对比情况, 采用均方误差(Mean Square Error,MSE)和均方 根误差(Root Mean Square Error, RMSE)可以衡 量预

29、测值同真值之间的偏离程度&MSE = 1, (S；.-.)2(11)RMSE = , 9，-s，)2(12)标记值S1S11的MSE、RMSE、R2情况如 表3所示&从表中可以看出MSE与RMSE与标记值本身数据量级有关，并且与其量级相比! MSE与RMSE的值处于较低水平即预测效果较 好，通过观察值也以很好佐证，除S11预测最 差外，其余均能达到0.98以上水平将每个网格点预测得到的标记值读入DN-CR计算程序迭代求解守恒方程，最终得到收敛 待预测K)=1.0时的流场为了直观对比计 算精度，如图14所示，选择丁 =0. 5m截线上 DNCR与NS# UGKS方法的物理量分布进行 对比01 0

30、00 2 000 3 000网格点编号 预测结果-0.04-0.06-0.081000 2 000 3 000 网格点编号寸籍蓝 Z籍京(g) 偌预测结果500.030.020.010-0.01-0.02-0.0301 000 2 000 3 000网格点编号(e)&“.预测结果0.04-0.04-0.0601 000 2 000 3 000网格点编号(c)Aj预测结果1 0009500国 0jg 500-1 00001 000 2 000 3 000网格点编号(f)A罂预测结果网格点编号5080封 -50-150-200网格点编号 (幻针预测结果图13预测情况Fig. 13 Predicti

31、on situation表3标记值误差特性Table 3 Error characteristics of predicted values误差特性S1S2S3S4S5S6S7S8S9S10S11MSE8. 12X10!61.39X10!52. 76X10!81. 76X10!81.06X10 825. 4810.850.9632.9116.74175.54RMSE2.85X10!33.72X10!31.66X10!41.33X10!41.03X10!45.053.290.985.744.0913.40R20.994 30.998 30.992 70.998 30.994 70.989 10.

32、 997 10 998 00 999 80 985 77 0 675 5采用式（13）计算评估精度提升情况:52-1NCR J2 UGKS)-采用式（13）计算评估精度提升情况:52-1NCR J2 UGKS)-(13)52-1:yHgks)-Ts.uos.0图14 Kn - 1.0计算结果对比(y = 0. 5 mFig. 14 Comparison of Kn - 1. 0 calculation results( y - 0. 5 m )式中：若与趋于1,则表明DNCR方法预测值与 UGKS基本一致；若与趋于0,则表示DNCR预测 值与N-S方程计算结果趋于一致&通过计算得到y = 0.

33、 5 m截线上U、P、T的 结果，DNCR相比N-S精度提升分别为92. 23%、 97. 77% 90. 80% &为进一步评估极端随机树的泛化能力，选取 Kn- 0.5与Kn - 1. 5两个超出训练集范围的努 森数条件开展DNCR计算，重点考察DNCR方法 的训练集外延数据泛化性能&其中y=0.5m截线 DNCR与N-SUGKS方法的y方向速度、温度与 压力物理量分布对比分别如图15、图16所示。 DNCR较N-Q方程预测精度提升情况依次为 7.59%、95.71%、71. 20% 和 89. 12%、94.36%、 87.69%,精度计算公式见式（13）,可以看出在 Kn - 0. 5与Kn - 1. 5条件下,DNCR较N-S方程 预测精度提升百分比与Kn - 1. 0状态时相比略有 下降，但仍保持较高水平，即DNCR方法中采用极 端随机树模型具有较好的泛化