基于人工神经网络的非结构网格尺度控制方法1)

王年华 , 鲁 鹏 * 常兴华 张来平 邓小刚

* (中国空气动力研究与发展中心空气动力学国家重点实验室,四川绵阳 621000)

† (西南科技大学信息工程学院,四川绵阳 621010)

** (重庆文理学院智能制造工程学院,重庆 402160)

†† (国防科技创新研究院无人系统技术研究中心,北京 100071)

*** (军事科学院,北京 100091)

引言

网格生成是计算流体力学(computational fluid dynamics,CFD)数值计算的第一步,也是未来CFD六大重要研究领域之一[1-2].在现代CFD 应用过程中,自动生成复杂构型的高质量网格(包括网格自适应)依然是一个重大挑战性问题.自动化程度和网格质量是网格生成过程中最重要的两个问题[3-4].据统计,网格生成通常占据整个计算周期大约60%的人力时间,高度自动化的网格生成方法无疑可以很大程度节约CFD 计算周期内的人工成本.

网格空间尺度分布控制在网格生成中至关重要,对于网格质量和流动求解精度影响较大.通常希望在几何曲率大、流动梯度大等重点区域加密网格,而在非重点区域则希望网格尽可能稀疏并均匀过渡.传统的非结构网格尺度控制方法主要有:函数指定法、插值类方法、背景网格法和根据流场特征进行自适应等.

函数指定法适用于简单问题的全场网格尺度控制,比如可以指定线性函数控制翼型的网格尺度分布;或用于局部网格尺度控制,比如可以指定几何级数或指数函数对局部网格进行加密[5].

采用背景网格法[6-9]进行网格尺度控制,需要先用规则的矩形结构网格或采用较稀疏的非结构网格覆盖全计算域,在计算域内分布一定的“点源”“线源”或者“面源”等局部网格分布控制参数,将这些源视作离散的“热源”,求解热传导方程(泊松方程),得到的稳态解即为全计算域的尺度控制参数分布[6].或者根据局部几何特征(曲率、狭缝、窄边等信息)确定局部网格尺度,再在背景网格上求解网格尺度满足的梯度限制方程,将局部网格尺度光滑到全场,得到背景网格上的网格尺度分布[7-9].在生成网格时,根据控制空间中某点在矩形背景网格中所处的单元,通过背景网格的尺度分布插值得到该点处的网格尺度.

除此之外,还有一些学者提出一些其他形式的背景网格法.如Deister 等[10]提出由最大最小尺度及最大曲率角对几何进行栅格化,计算得到局部网格尺度,并存储在自适应的背景笛卡尔网格上.Quadros等[11-12]提出采用几何体离散骨架的几何临近信息、特征尺寸、边界曲率来测量几何复杂度,并根据几何复杂度生成点源,根据点源确定网格尺度分布,最终将网格尺度存储在叉树结构的笛卡尔背景网格上.Ruiz-Girones 等[13]提出通过在背景网格上求解一种新的非线性方程来控制四边形网格尺度,等等.

径向基函数(radial basis function,RBF)[14-15]可用于数据插值,在已知边界上的网格尺度后,可用RBF 方法将边界网格尺度插值到内场,采用贪婪算法还可以一定程度提高RBF 插值的效率,是一种简单高效的插值类尺度控制方法.

根据流动特征物理量的梯度量等判据进行网格自适应[16-17]也是控制网格尺度分布的一种有效方法,能够根据流场变化控制网格疏密,在梯度大的区域生成更密的网格,能够更精细地捕捉流动特征,具有更高的计算精度和适应性,是来重要发展方向之一[18].

近年来,基于人工神经网络的深度学习方法在工业社会甚至流体力学领域得到广泛研究和应用[19-21].大数据驱动的人工智能方法成为除理论分析、数值计算和实验技术以外一种新的研究范式,为各个领域带来了新的研究思路和方法.

在网格生成领域,经过多年的工程实践,已经积累了大量各种类型的网格数据,这些数据包含了网格生成规则及技术人员在网格生成方面的知识和经验,是天然的机器学习训练样本.通过机器学习对网格生成规则进行学习,可以简化传统算法,提高网格生成效率[22].而网格尺度分布的控制也需要技术人员根据对流动问题的分析和经验合理确定,采用机器学习方法对网格尺度分布进行控制,有望减少人工工作量和对人工经验的依赖.

本文从网格质量、效率、灵活性和自动化程度4 个方面综合分析各类网格尺度控制方法的优缺点.为了克服传统背景网格插值法效率低、自动化程度不高等方面的不足,本文从效率和自动化程度角度提出两种网格尺度控制方法,首先将RBF 插值方法应用于网格尺度控制,采用贪婪算法对RBF 插值参考点序列进行精简,实现高效的RBF 网格尺度分布控制方法.进一步将提出一种采用人工神经网络进行非结构网格尺度控制的方法,通过引入相对壁面距离和相对网格尺度,初步确定合理的神经网络输入输出参数,建立人工神经网络训练模型,采用商业软件生成二维圆柱和二维翼型非结构三角形网格作为训练样本,通过训练和学习建立起相对壁面距离和相对网格尺度之间的映射关系,进而实现不同密度的二维圆柱、不同二维翼型在不同远场大小情况下的网格尺度分布控制.

1 基于背景网格法的尺度控制方法回顾

背景网格法可采用规则的笛卡尔直角结构网格、非规则的非结构网格或者自适应笛卡尔网格[9],各种类型的背景网格各有优缺点,由于规则结构网格求解和插值效率较高而得到广泛使用.

在背景网格上布置“点源”、“线源”或者“面源”等局部网格分布控制参数,通过求解热传导方程,得到稳态解即为全流场的尺度控制参数分布[5].该方法能够生成分布均匀的网格,且能更灵活地考虑流动局部特征对网格分布的影响.

1.1 背景网格法

以二维问题为例,网格尺度S满足如式(1)所示的稳态热传导方程(泊松方程),以确保尺度源项的作用在物理空间内光滑分布.

式中G为源项,定义如式(2)所示

式中下标“i”,“j”代表背景网格节点,N为热源总数,Ψn为第n个源的强度因子,函数In及Jn分别定义如下

式中S和f为在各源处的尺度,r为背景网格节点到源的距离,|ln|为线源长度.

文献[5]还提出了考虑方向性的网格密度控制方法,可以通过改变源的强度函数来实现,具体可参考文献,本文不再赘述.

因此,采用背景网格法进行网格尺度控制的步骤为:

(1)在矩形背景网格上采用中心差分离散泊松方程;

(2)采用Gauss-Seidel 超松弛迭代求解离散的泊松方程,得到背景网格上的网格尺度分布;

(3)在网格生成过程中,根据当地位置在背景网格中插值得到当地网格尺度,用于控制网格生成过程.

1.2 背景网格法网格生成实例

本节采用二维圆柱、NACA0012 翼型和30P30N多段翼型作为考核算例,对背景网格法进行实例测试.本文算例中采用的三角形网格生成算法为作者发展的基于ANN (artificial neural network)的阵面推进法[22],该方法在传统阵面推进法的基础上,通过引入ANN 进行生成模式判断和新点预测,减少了相交性判断次数,网格生成效率提高30%,其网格生成步骤可简要概括为:

(1)将几何边界离散成初始阵元;

(2)从最小阵面出发,自动选择网格模板点,人工神经网络根据网格模板判断生成模式并预测新点坐标;

(3)根据生成模式、新点坐标和局部网格尺度生成新网格单元;

(4)判断新单元是否合适,合适则更新数据结构;

(5)回到步骤2,直至所有面变成非活跃面,整个计算域被网格填满.

图1～图3 给出了3 个算例点源设置示意图和生成的网格,图中蓝色圆点即为点源所在位置.3 个算例人工设置的点源数量分别为12,24 和44.图中背景网格仅作为示意,实际背景网格节点数量需要根据边界网格尺度进行调整,以确保对网格尺度分布场的分辨率和插值精度.比如对于在前后缘网格尺度较小的翼型,背景网格必须足够密才可以有效反映出空间网格尺度的变化.

图1 圆柱算例点源设置及网格生成情况Fig.1 Nodal source settings and corresponding triangular mesh over a 2D cylinder

图2 NACA0012 算例点源设置及网格生成情况Fig.2 Nodal source settings and corresponding triangular mesh over NACA0012 airfoil

图3 30P30N 算例点源设置及网格生成情况Fig.3 Nodal source settings and corresponding triangular mesh over 30P30N airfoil

在本文三个算例中,背景网格规模分别为51 × 51,301 × 301 和401 × 401.由于是在背景网格上迭代求解泊松方程,因此背景网格的网格数量直接决定了迭代求解的效率.

图中结果显示,在人工设置合适的点源参数后,可以在背景网格上得到恰当的网格尺度分布.

要说明的是本文暂未采用线源,实际上线源可以看作按线段排列的具有一定强度分布的点源集合,因此在处理本文的简单问题时,只采用点源进行网格尺度的控制,用以说明背景网格法的优缺点.

2 基于RBF 方法的尺度控制方法

径向基函数是一种常用的插值函数,常用于网格变形的插值[14-15,23-24],如图4 所示即为用NACA-0012 翼型的S 型启动来模拟鱼的游动过程,采用RBF 插值方法将物面的变形插值到空间来生成动网格.RBF 插值方法也可以用于网格尺度的插值控制,在已知边界网格尺度的情况下,求得插值系数矩阵,可将边界尺度插值到整个计算域.

图4 RBF 网格变形方法Fig.4 Mesh deformation controlled by RBF interpolation

2.1 基于贪婪算法的RBF 插值方法

RBF 插值方法是比较成熟的插值方法,对于一个变量场,如网格变形情况下的位移场,或网格尺度控制情况下的网格密度场,可以用RBF 插值公式表示为

其中,N为参考点数目,f(r) 为某待求点的函数值(网格尺度),r为待求点的位置矢量,ri为参考点的位置矢量,φ 为RBF 基函数,‖r-ri‖ 为参考点与待求点之间的欧氏距离.wi为第i个参考点的权重系数.

以网格尺度分布控制为例,权重系数满足

式中Sp为各个参考点的网格尺度,参考点通常为边界点,也即边界网格尺度.基函数的类型包括全域型基函数和紧支型基函数,具体可以参考文献[14-15,23-24],本文不再赘述.

权重系数的求解通常需要求解系数矩阵的逆,在参考点数目较大时,会导致求解效率较低,因此可以引入贪婪算法对参考点进行精简.

基于贪婪算法的RBF 方法的步骤可简要归纳为:

(1)初始参考点集为空集,随机取一边界点作为初始参考点;

(2)根据参考点,采用RBF 插值计算其他边界点处的网格尺度;

(3)由于参考点数目太少,RBF 插值得到的边界网格尺度与给定值存在一定误差,找出误差最大的点;

(4)若误差最大的点不为已有参考点,则将该点做为新参考点加入参考点序列,否则重新任选一点加入参考点序列;

(5)重复步骤(2)～ (4),直到最大误差或者参考点数目满足要求,确定最终参考点序列;

(6)根据最终参考点序列进行RBF 插值,得到空间所有位置的网格尺度分布.

2.2 RBF 方法网格生成实例

本节仍然采用二维圆柱、NACA0012 翼型和30P30N 多段翼型作为考核算例,对RBF 方法进行实例测试.算例中所采用的三角形网格生成算法仍为基于ANN 的阵面推进法,具体可参考文献[22].

RBF 基函数取为Wendland’s C0,紧支半径取为计算域范围大小的1/4.网格尺度控制对尺度插值误差的要求不高,本文以尺度的最大相对误差不超过5%为标准,进行参考点序列精简,表1 给出了采用贪婪算法对参考点进行精简的结果,结果显示对于翼型算例,精简序列后参考点数量至少减少一半,RBF 插值的耗时也减少了近一半,在保证了插值效果的同时提高了插值效率.

表1 RBF 方法参考点的数目及插值耗时Table 1 Number of reference nodes and time consumption on interpolation

图5～ 图7 给出了精简后的参考点的位置及根据精简参考点进行尺度控制生成的非结构网格.结果显示RBF 方法能够根据边界网格尺度插值得到空间的网格尺度,插值得到的网格尺度分布过渡均匀.由于初始参考点的选择是随机的,参考点的分布具有一定的随机性,但是均能够保证插值误差满足要求.

图5 圆柱算例精简后的参考点及网格生成情况Fig.5 Reference nodes corresponding triangular mesh over a 2D cylinder

图6 NACA0012 算例精简后的参考点及网格生成情况Fig.6 Reference nodes and corresponding triangular mesh over NACA0012 airfoil

图7 30P30N 算例精简后的参考点及网格生成情况Fig.7 Reference nodes and corresponding triangular mesh over 30P30N airfoil

3 基于人工神经网络的尺度控制方法

基于人工神经网络的深度学习具有较强的非线性拟合能力,能够通过现有样本数据的训练识别出数据中隐含的非线性映射关系.网格尺度分布实际上是一个关于几何特征、流场特征的非线性映射,如图8 所示,几何特征(物面曲率、狭缝、细小结构等)可以直接影响网格尺度分布,几何特征也可以决定流动特征(梯度量等),从而间接决定网格分布.不同几何外形的网格尺度分布不同,同一几何外形在不同来流条件(流场特征)情况下,网格尺度分布也不同.

图8 网格分布控制与几何特征和流场特征的关系Fig.8 Relationship between mesh size control,geometry,and flow features

流场特征与几何外形、来流条件和边界条件相关,目前采用机器学习方法进行流场预测是一个热点研究问题[25-26].而根据流场特征确定网格分布则可以根据网格自适应的相关准则进行,也可以采用机器学习的方法进行[27].本文初步考虑几何外形对网格分布的影响,采用ANN 建立几何特征与网格尺度分布之间的关系.

3.1 ANN 输入输出参数

Chedid 和Najiar 等[28]曾尝试用ANN 建立起网格密度与几何特征之间的关系,如图9 所示,神经网络输入考虑了空间点gi到边界点的最小距离 δ1,次小距离 δ2及其对应的夹角a*和a**及到夹角对应两边的投影距离j1a,j1b,j2a和j2b共8 个输入参数,输出该空间点处的网格密度,定义为一定大小区域内节点的数量.该方法能够较为全面地反映空间点与边界之间的几何关系(距离及投影距离),以及空间点对应的边界处的几何特征(夹角、曲率),但计算量较大不利于提高尺度控制效率,同时由于无法网格局部加密控制,失去了对网格尺度控制的灵活性,因此优势并未得到发挥.

图9 文献[28]中神经网络的输入参数Fig.9 Input parameters for the artificial neural network in Ref.[28]

为提高效率,本文初步选择最小壁面距离wdist作为输入参数,网格尺度Sp作为输出参数.为提高ANN 的泛化性,其输入输出参数通常需要进行归一化操作,本文采用计算域的大小Lr_d,远场边界网格尺度Lr_f和物面边界网格尺度Lr_w对输入输出进行分别归一化.

同时,由于从物面到远场,网格尺度变化较大,网格尺度相对值可能在103量级以上,为尽量缩小输入输出的值域范围,提高ANN 训练效果,本文对输入输出参数进行开根号.此外,在物面附近网格尺度变化快,需要同时考虑物面参考值和远场参考值进行归一化,具体输入输出参数形式如表2 所示.

表2 ANN 输入输出模型Table 2 Parameter model for artificial neural network

3.2 训练方法及参数

本文基于Matlab 神经网络训练工具设计全连接的人工神经网络,网络含有1 个输入层,1 个隐藏层,1 个输出层.输入层含有1 个神经元,输出层含有1 个神经元,隐藏层的神经元数量为10 个,激活函数采用Sigmoid 函数,损失函数为均方误差函数,训练方法采用Levenberg-Marquardt 反向传播方法,该方法比常规的梯度下降反向传播算法训练效率更高,具体可以参考文献[29].神经网络的结构如图10所示.

图10 基于Matlab 的人工神经网络训练工具Fig.10 Artificial neural network training tool based on Matlab

选择二维圆柱网格和NACA0012 翼型网格作为训练样本,如图11 所示,网格面个数分别为4295,3950,每个网格面对应一组样本数据点.网格训练样本按70%,15%和15%的比例随机划分为训练集、测试集和验证集,图12 给出了在翼型三角形网格训练集上的Loss 值及在验证集上的预测精度收敛历程.结果显示:经过120 次迭代后,Loss 值及预测误差均下降到了0.001 58 左右.

图11 网格分布训练样本网格Fig.11 Sample grids for ANN training

图12 训练Loss 值和精度收敛历程Fig.12 Convergence of loss and accuracy on sample grids

3.3 ANN 预测结果

基于前述的ANN 训练结果,分别预测生成了不同密度的二维圆柱网格、NACA0012 密网格、RAE2822 翼型以及三段翼型,在物面附近和远场均取得了较好的效果,同时还能适应不同远场大小的情况,如图13 所示.算例中所采用的三角形网格生成算法仍为基于ANN 的阵面推进法,具体可参考文献[22].

图13 ANN 模型各向同性网格预测结果Fig.13 Mesh size controlled by ANN model for isotropic triangular grids

图13 ANN 模型各向同性网格预测结果(续)Fig.13 Mesh size controlled by ANN model for isotropic triangular grids (continued)

同时,本文还将ANN 应用于各向异性混合网格尺度控制,各向异性四边形采用层推进逐层推进生成[30-32],而各向同性三角形采用基于ANN 的阵面推进生成[22].层推进的推进方向、多方向推进数量均采用ANN 预测,同时在凹角处考虑局部推进步长,避免网格相交,具体方法细节可以参考文献[30].为与各向同性网格尺度控制相一致,本文将层推进生成的最后一层网格作为虚拟物面,用于计算最小壁面距离,作为ANN 控制网格尺度的输入参数.另外,层推进的网格尺度控制仍然采用指定物面第一层网格高度和增长率的方式给定.图14 给出了NA0012翼型和30P30N 三段翼型的生成结果.由图中结果可见,网格分布均匀合理,网格质量满足要求,说明本文发展的ANN 方法可应用于各向异性混合网格的尺度控制.

图14 ANN 模型各向异性混合网格预测结果Fig.14 Mesh size controlled by ANN model for anisotropic hybrid grids

4 背景网格法、RBF 方法与ANN 方法的比较

本文分别采用3 种方法生成了几个典型几何外形的非结构网格,其网格生成质量、生成效率、尺度控制灵活性、自动化程度等情况存在一定差别,本节对3 种方法进行对比分析.

从网格质量的角度,在网格尺度较小的部位,如翼型前后缘、狭缝等,背景网格法要求加密背景网格,以分辨最小网格尺度,而RBF 方法和ANN 方法不依赖于背景网格,因此在尺度较小的部位得到的尺度分布控制效果比采用矩形背景网格要好.

从自动化程度和灵活性角度,采用在背景网格上设置点源控制网格分布,对分布的控制最为灵活,能够根据几何特征和流场特征预先设置点源来改变网格分布.而背景网格法需要人工设置点源,本文算例中,在远场和翼型附近人工设置数十个点源,每个点源人为给定强度、位置等参数,需要一定的经验,自动化程度较低.而RBF 方法和ANN 方法只需要给定离散的边界节点,就可以得到空间网格尺度分布,自动化程度相对较高,但是灵活性较低.

从网格生成效率的角度,背景网格法需要在笛卡尔网格上迭代求解泊松方程,在背景网格很密时,求解效率较低.而RBF 方法虽然需要求解矩阵的逆,但在引入贪婪算法之后,参考点数量减少,RBF 方法的效率得到提高.ANN 方法虽然需要求解最小壁面距离,但壁面距离的值并不需要十分精确,因此可以采用一些近似求法,也可以达到更高的效率.

表3 给出了3 种方法在生成网格过程中尺度控制所耗费的时间,由于每种方法生成的网格单元数存在一定差异,因此表中也给出了网格单元数量.

表3 3 种方法控制网格尺度耗时对比Table 3 Efficiency comparison of the three methods

背景网格法耗时基本由背景网格的规模及点源的数量决定.在本文三个算例中,背景网格规模为51 × 51,301 × 301 和401 × 401,点源的数量大致为20 个,因此在采用背景网格法时,翼型算例耗时明显增加.

由表3 中数据可见,在背景网格数量较少时,背景网格法的效率较高,但是随着背景网格数量和点源数量的增加,背景网格法效率明显下降.相较于传统背景网格法,RBF 方法和ANN 方法的耗时明显减少,耗时仅为背景网格法耗时的1/10～ 1/5,网格尺度控制效率相应提高了5～ 10 倍.而且,在外形相对复杂的情况下,ANN 方法展现了更好的控制效率.可以预见,在三维复杂外形情况下,ANN 的控制效果和效率会更高.

5 总结与展望

本文提出了一种采用ANN 进行网格尺度分布控制的方法,初步确定了合理的神经网络输入输出参数,基于Matlab 建立人工神经网络模型,采用商业软件生成二维圆柱和NACA0012 翼型非结构三角形网格作为网格密度训练样本,通过训练和学习建立起壁面距离和网格尺度的映射关系,实现了不同密度的二维圆柱和不同二维翼型的网格尺度预测.同时,发展了基于RBF 方法的网格尺度控制方法,采用贪婪算法对插值参考点序列进行精简,将RBF 插值效率提高了一倍.

与传统背景网格法相比,基于ANN 方法和RBF方法的网格尺度预测效率提高5～10 倍,有助于进一步提高网格生成效率.最后将基于ANN 的网格尺度控制方法拓展应用于各向异性混合网格的尺度控制,得到的网格质量满足要求,证明了方法的实用性.

展望ANN 方法在尺度分布控制领域的应用前景,在以往的CFD 研究和工程实践中,已经积累了大量各种类型的网格,其中包含了网格尺度分布信息.对已有的网格数据按照几何外形进行分类,训练得到可以处理不同类别几何外形的神经网络,在需要对新的几何外形进行网格尺度控制时,只需选择已经训练好的该类别的神经网络即可高效快速完成尺度控制.进一步还可以考虑流场特征(来流条件),对神经网络的适用范围进行细分,使得其能够同时考虑几何外形和流场特征进行尺度控制.

可以进一步考虑物面几何曲率、局部特征尺寸(窄边、狭缝)等参数,提高神经网络的泛化性,同时研究通过神经网络进行流场预测,并根据流场特征进行网格自适应.