基于径向基函数和Delaunay图映射的高效高鲁棒性的非结构网格变形方法

2024-01-08王昊达崔晓春

气体物理 2023年6期

王昊达, 刘南, 张颖, 崔晓春

(中国航空工业空气动力研究院, 辽宁沈阳 110034)

引言

在使用数值模拟方法求解复杂流场时, 须对控制方程进行空间离散和时间离散, 然后求解离散后的线性方程组得到计算网格上的流场变量。在对空间进行离散时往往须将计算域切分为若干计算网格。当计算域的边界运动或者变形时, 须利用网格变形方法生成新的网格。网格变形效率和鲁棒性是限制其应用的关键因素, 尤其是在气动外形优化、气动弹性计算等涉及大规模网格变形的问题中[1-2]。

根据采用的计算模型的不同, 将网格变形方法可大致分为两类: 物理模型法、代数插值法[1-2]。物理模型法将网格中的点、线和单元视为可变形的实体, 考虑计算网格节点之间的连接关系, 可以有效地避免网格单元交叉。物理模型法主要包含: 弹簧类方法和弹性体方法等。弹簧类方法的应用比较广泛, 首先由Batina[3]提出, 其基本思想是将计算域内两节点间的连接线看作一根拉压弹簧, 将整个网格看作一个弹簧系统, 通过求解该弹簧系统的静平衡方程求解节点位移。Farhat等[4]随后提出了一种扭转弹簧法, 用以防止网格单元的交叉。弹性体法最早是由Tezduyar[5]提出的, 作为弹簧类法的延伸, 其核心是将整个计算域内的网格看作一整个弹性体, 通过对线性弹性方程组的求解来获得各网格节点的位移。由于弹性体法要耗费大量时间迭代求解大型的控制方程组, 所以其变形效率很低, 因此工程上不适用于网格数量很多的复杂外形以及三维问题。代数插值法是指用某种数学插值方法将物面边界的变化插值到内部节点上, 其优势在于不需要计算域内部网格节点间的连接信息, 只需要网格节点坐标。主要包含: 超限插值法(transfinite interpolation, TFI)、反距离加权函数插值法(inverse distance weighted, IDW)、径向基函数插值法以及基于Delaunay图映射方法。TFI方法最早由Gaitonde等[6]提出, 虽然其变形效率高、网格质量好但仅能应用于结构网格。Liu等[7]提出的DGM方法, 通过构建背景网格与计算网格的映射关系求解变形后的网格节点坐标。该方法效率较高但不适用于大变形。与之相反, Rendall等[8-9]提出的RBF方法以及Witteveen等[10]提出的IDW方法在大变形时能够保证较高的网格质量但计算效率很低[11]。

目前的各类网格变形方法很难同时满足网格变形效率和网格变形质量的需求[2], 因此本文通过网格聚合算法自动建立粗化网格作为背景网格; 通过RBF方法建立背景网格和物面计算网格之间的代数关系; 并通过DGM方法建立背景网格和所有计算网格之间的映射关系, 结合RBF方法变形能力强和DGM方法变形效率高的优点, 建立一种高效、高鲁棒性的网格变形方法。最后通过不同算例验证本方法的网格变形能力和效率, 以及不同网格粗化参数对网格变形能力的影响。

1 方法描述

1.1 RBF&DGM网格变形方法

利用网格聚合方法结合RBF和DGM方法的优点, 建立RBF&DGM网格变形方法。方法的主要思路是: 通过网格聚合算法对计算网格进行粗化并自动生成一套背景网格, 然后结合物面计算网格点、计算域边界点以及4个极大点构建Delaunay控制体, 最后建立Delaunay控制体与计算网格点之间的几何映射关系。当物面边界发生变化时, 通过RBF插值方法更新背景网格, 物面计算网格通过对应的变化规律进行更新, 空间计算网格结合更新后的背景网格点坐标通过建立的几何映射关系进行更新。由于网格粗化是从物面边界出发一层层推进的, 因此可以更好地保证物面附近黏性网格的正交性, 并且除计算边界网格点、物面边界网格点外, 背景网格中还存在空间计算网格点, 因此通过背景网格建立的几何映射关系更为精确, 能更好地传递物面边界的变化规律。为方便程序使用, 本文针对二维网格以及三维网格均直接构造三维Delaunay控制体进行几何映射关系的建立。

本文方法主要分为前处理和网格变形两步, 前处理的作用是自动生成背景网格并建立Delaunay控制体以及Delaunay控制体与计算网格之间的几何映射关系, 网格变形的作用是利用上述几何映射关系和物面计算网格点的变化规律更新计算网格, 流程分别如图1、图2所示, 图中下标f表示背景网络中的物面网格节点, 下标s表示背景网格中的空间网格节点。

图1 前处理步骤Fig. 1 Pre-processing steps

图2 网格变形步骤Fig. 2 Grid deformation steps

前处理步骤如下:

1) 通过网格聚合方法建立粗网格作为背景网格;

2) 利用RBF方法建立背景网格中物面和空间网格点变形之间的关系;

3) 利用物面网格点、计算域边界点、 RBF控制点、 4个极大点建立Delaunay控制体;

4) 构建Delaunay控制体与所有计算网格点之间的映射关系, 即得到每个空间网格点对应的4个背景网格点, 以及相应的权重系数。

网格变形步骤如下:

1) 根据输入获取物面变形量;

2) 利用RBF方法, 根据物面网格点变形得到背景网格中空间网格点的变形, 更新背景网格;

3) 利用前处理步骤中建立的几何映射关系, 根据物面网格点变形和背景网格点变形计算得到所有空间网格点变形, 更新计算网格。

1.2 网格聚合方法

混合网格的聚合方法较结构网格复杂很多, 本文参考Berglind[12]提出的方法。在该算法中, 强制将前沿阵面从物面边界处逐层向外扩展, 通过设置每个网格节点的优先级来完成(该优先级表征每个节点与壁面的最小连接数), 以确保每一个附面层网格单元每一层顶点都具有相同的优先级。

如图3所示, 所有节点的优先级由以下算法来确定:

图3 网格节点优先级Fig. 3 Priority numbering of all vertices

1) 设置所有网格节点的优先级为p=-1;

2) 设置所有物面边界网格节点的优先级p=0;

3) 创建与所有优先级为p的网格节点相邻的边构成的前沿阵面;

4) 将所有与非负优先级节点相邻节点的优先级加1,p=p+1;

5) 返回第3)步, 直到所有网格节点的优先级均非负。

优先级用于种子点选择的标准。以下选择标准按给定顺序, 从前沿阵面的候选节点中选择种子点:

1) 具有最低优先级数p的网格节点。

2) 与已聚合的控制体间连接数最大的网格节点。

为了保证聚合后的粗网格物面附近的正交性, 本算法中规定了不同区域的不同定向聚合策略: 在附面层网格区域中(二维为四边形网格单元、三维为棱柱形网格单元), 网格聚合的方向只在法向方向而切向不聚合。在附面层网格以外的区域, 长宽比将作为定向聚合的判据。长宽比被定义为三维V/S1.5, 二维V/L2。其中V是控制体的体积,S是表面积,L是控制体积的周向长度。

粗化比R用于定义与粗网格控制体聚合的细网格控制体的数量。Rp、Rt分别表示具有棱柱状元素和四面体的区域的粗化比, 二者在大多数情况下会有所不同。因此本文通过全局平均粗化比来控制网格粗化程度。

在以下的描述中, “点”特指对偶网格的控制单元。具体网格聚合算法可概述如下:

1) 在初始前沿阵面中根据标准选择一个“种子点”(初始聚合网格)。

2) 检查种子点是否位于边界层内。

3) 若种子点位于边界层内, 则只与优先级比先前聚集的节点高一个单位的相邻细网格点沿法向聚合。若种子点位于边界层之外, 则从其相邻的控制体中聚合使其长宽比最大化的点, 生成一个粗网格控制单元。

4) 更新该粗网格控制单元的相邻细网格点列表, 重复步骤3)直至达到所需的网格粗化比。

5) 更新前沿阵面, 将与最后被聚合的细网格点相邻的未被聚合的点加入到前沿阵面。

6) 若仍有细网格点未完全聚合, 则重复步骤1), 直至所有细网格点均被聚合。

7) 将所有的孤立点与和它相邻的一个粗网格控制体进行融合(清除孤立点)。

对于同一初始网格, 本文方法可以进行多次网格聚合算法对初始网格进行粗化, 该过程是通过给定网格粗化层数l确定的。

本文通过网格聚合算法得到粗化后的背景网格点如图4所示。

(a) Original grid

1.3 RBF网格变形方法

RBF插值函数形式如下[13]

式中,s(x)是x处径向基函数值,φ(‖·‖)是基函数, 一般选择Gauss或Wendland′s C2函数,i代表RBF中心,xi是各中心点坐标。对于网格变形问题,p(x)=0。令函数值等于物面计算网格点位移, 并写成矩阵和向量形式

Δxb=Cbbax
Δyb=Cbbay
Δzb=Cbbaz

(1)

其中x方向(y和z方向与之类似)如下

因此,Cbb矩阵如下

其中径向基函数为

φbibj=φ(‖xbi-xbj‖)

式中, 下标b代表流场物面计算网格点,Nb为网格点数, 令函数值等于空间点位移则有

Δxv=Avbax
Δyv=Avbay
Δzv=Avbaz

(2)

其中Avb矩阵如下

式中, 下标v代表网格聚合后的背景网格点,Nv为网格点数。利用式(1)和(2)得到物面计算点和背景网格点位移之间的关系, 如下

(3)

1.4 DGM网格变形方法

首先在平面或空间上给定的包括流场边界点和网格聚合后的点在内的一组点的基础上, 对整个计算域进行唯一的Delaunay图三角化, 具体算法可参考Leatham[14], 从而完成计算域Delaunay控制体的构建。随后建立计算网格节点与Delaunay控制体节点之间的几何映射关系。几何映射关系通过相对面积(二维)和相对体积(三维)坐标建立, 如图5所示。结合背景网格, 通过DGM方法建立的Delauny控制体如图6所示。

(a) Two-dimensional case

(a) Triangle

对于二维情况如图5(a)所示, 任意网格节点O位于三角形单元△MNQ中, 由节点O和三角形组成的3个三角形△MON, △NOQ, △QOM的面积分别为S1,S2,S3, 令

对于三维情况如图5(b)所示, 任意网格节点O位于四面体单元MNPQ中, 4个四面体OMNQ,OMQP,ONPQ,OMNP的体积分别为V1,V2,V3,V4, 令

式中,Si表示各顶点对应三角形的面积,Vi表示各顶点对应的四面体体积,S表示三角形单元面积,V表示四面体单元体积,λ表示表征几何映射关系的定位参数。当背景网格更新后, 根据定位参数可实现计算网格节点坐标的更新。

对于二维网格有

对于三维网格有

(4)

式中,XO为任意计算网格节点坐标,X′为更新后的背景网格点坐标。

通过式(4)可知, 当完成背景网格的更新后, 通过λ即可计算出任意计算网格节点坐标, 完成计算网格的更新。

2 方法测试

通过算例对本文所建立的基于RBF&DGM的网格变形方法进行测试, 并与RBF、 DGM方法进行对比, 验证本文方法的网格变形能力、网格变形效率。同时改变不同网格粗化参数, 研究不同网格粗化参数下生成的背景网格对网格变形质量的影响。

2.1 网格质量验证标准

本文采用尺寸斜交度量函数fsize-skew定量衡量网格质量[15], 计算公式如下

fsize=min(τ, 1/τ),τ=α/w

对于三角形单元

对于四面体单元

式中,fsize为相对尺寸度量函数, 用于检验目标网格单元相比于参考网格单元的尺寸是否合适, 选取变形前网格单元作为参考单元, 则w为参考网格单元面积的2倍,α为目标单元面积的2倍。当fsize=1时, 表示目标单元和参考单元面积相等, 当fsize=0时, 表示目标单元已经退化。

fskew为形状斜交度量函数, 用于检测网格单元形状的扭曲程度, 对于三角形单元选取正三角形作为参考网格单元, 对于四面体单元选取正四面体作为参考网格单元,λ通过构造张量度量矩阵获得。当fskew=1时, 表示目标单元为正三角形或正四面体,fskew=0时, 表示目标单元已退化。

fsize-skew为尺寸斜交度量函数, 其取值范围在0～1之间。当单元退化或具有负面积/负体积时,fsize-skew值为0, 当fsize-skew为1时, 说明此时变形后的网格为理想的正三角形或正多面体。

2.2 二维矩形算例

本算例测试网格变形方法在大变形下的适应能力。采用长a=1 m, 宽b=0.5 m的矩形作为物面边界, 正方形作为网格计算域, 网格数量9 100, 网格种类为非结构网格且不生成黏性网格, 原始网格如图7所示。

(a) Global view

绕矩形中心做旋转变形, 旋转角度θ以Δθ=10°, 从10°增加至90°。对3种方法变形后的网格进行质量评估, 分别对比最小质量系数和平均质量系数, 结果如图8所示。

(a) Minimum quality coefficient

由图8可知: 1) 由于矩形网格没有生成黏性网格, 以正三角形作为斜交度量函数的参考网格单元, 3种方法下变形得到的网格的最小质量系数和平均质量系数均较高; 2) 随着旋转角度的增加, 最小网格质量系数曲线快速下降, 平均网格质量曲线变化较缓, 仅有小部分物面附近的网格因旋转变形而受到挤压, 大部分计算网格并未受到过多影响; 3) 本文方法与RBF方法的最小、平均质量系数曲线几乎重合, 说明本方法具有与RBF方法相当的网格变形能力。

图9～图11分别给出了物面旋转角度θ=90°大变形下, 3种方法的网格变形情况。从图中可以看出, 本文方法和RBF方法可以很好地将物面的变形传递至空间网格中, 而DGM方法在大变形时对物面变形的传递能力较差, 因此只有少部分计算网格受到物面变形的影响, 其平均网格质量并未与其他两种方法相差过大。

(a) Global view

2.3 网格粗化参数影响

背景网格的生成是本文方法的基础, 不同的网格粗化层数、网格粗化比会影响到背景网格的密度, 进而影响几何映射关系建立的精准度。因此本节将研究网格粗化层数l、网格粗化比R对网格变形质量的影响。

图12给出了矩形旋转θ=90°后网格质量随网格层数l的变化曲线, 网格变形质量以RBF方法作为参考, 径向基函数半径r取7.5。

图12 不同l下的网格变形质量Fig. 12 Grid deformation quality with different l

网格聚合算法共设置网格粗化层数5层, 随着网格层数的增加, 背景网格密度逐渐变小, 分别取每一层粗化后的网格作为背景网格并建立Delaunay控制体及几何映射关系。由图12可以看出: 随着l的不断增大, 网格平均质量几乎与RBF方法一致, 网格最小质量呈现先上升后下降的趋势, 当选取第5层粗化网格时, 最小网格质量已经接近于0, 说明此时已经出现网格单元交叉情况。

图12的结果表明, 背景网格的层数不宜过大也不宜过小, 应选取较为适中的背景网格密度所对应的网格层数l。当l为1时, 此时背景网格密度较高, 背景网格单元尺寸较小, 建立的几何映射关系过于精细, 导致过度传递物面的旋转变形, 使得物面附近的计算网格受挤压更严重。因此, 最小网格质量小于RBF方法。当l为5时, 此时背景网格密度较低, 背景网格单元尺寸较大, 建立的几何映射关系较为粗糙, 导致无法传递物面的旋转变形, 使得出现网格单元交叉的情况。因此在选取网格粗化层数时, 应当选择背景网格密度适中的粗化网格层数建立背景网格。

但单一地通过网格层数选取合适的背景网格密度有时是很难实现的, 因此还可以通过调节网格粗化比R来进一步选取合适的背景网格密度。选定第3层粗化网格研究R对网格质量的影响, 见图13。

图13 不同R下的网格变形质量Fig. 13 Grid deformation quality with different R

粗化比R表示与粗网格控制单元聚合的细网格控制单元的数量,R越大说明网格粗化能力越强、背景网格的密度越小, 图14中展示了不同R下的背景网格。由图13可以看出, 随着R的增加平均网格质量先增大后减小, 最小网格质量快速减小, 且当R大于2.5时最小网格质量约等于0, 出现网格单元交叉的情况。

(a) R=2

图13的结果表明, 在固定网格层数的情况下, 当R<2.75时, 随着R增大, 生成的背景网格密度逐渐变小, 导致对物面变形的传递能力逐渐减小, 最终出现网格单元交叉情况, 但空间网格受物面变形影响较小, 因此平均网格质量逐渐增加。但当R>2.75时, 此时背景网格密度过大, 导致物面附近第1层网格高度变大, 受物面变形影响导致更多的空间网格出现网格单元交叉情况, 因此平均网格质量开始逐渐减小。

综上, 在构建背景网格时, 应当将网格粗化层数与网格粗化比结合考虑, 不同的网格层数对应着不同的最佳网格粗化比。本文建议: 为方便调节参数, 在进行网格粗化时可选取第3层粗化网格作为基础, 若背景网格密度过大则适当减小网格粗化比R, 若背景网格密度过小, 则适当增大网格粗化比R, 进而选取更为合适的背景网格密度。

2.4 二维翼型算例

本算例测试本文方法对黏性网格物面正交性的保持能力。采用NACA64A010翼型, 网格数量为56 000, 计算域为r=30l0(l0为翼型弦长)的圆形域, 网格种类为非结构网格, 在物面处生成25层黏性网格, 原始网格如图15所示。

(a) Global view

绕翼型1/4弦线处作旋转变形, 旋转角度θ以Δθ=10°, 从0°顺时针旋转至90°, 并对3种方法变形后的网格进行质量评估, 分别对最小质量系数和平均质量系数进行对比。翼型旋转变形通过式(5)定义

(5)

由图16可知: 1) 本文所用网格质量系数fsize-skew中的斜交度量函数是以正三角形作为参考网格单元, 黏性网格的存在, 使得3种变形方法变形后的最小网格质量系数和平均网格质量系数均较小。2) 随着旋转角度θ的增加, DGM方法的网格质量曲线快速下降, 但RBF和本文方法的网格质量曲线变化较缓, 表明DGM方法物面附近黏性网格无法随物面运动, 黏性网格单元受挤压影响较大, 而另外两种方法能够使黏性网格跟随物面一起运动。3) 本文方法与RBF方法的网格质量曲线基本相同, 说明本文方法对于黏性网格同样具有较好的网格变形能力。

(a) Minimum quality coefficient

图17～图19分别给出了翼型顺时针旋转至45°时3种方法下的网格变形情况。可以看出: 1) 在物面附近, 本文方法能够与RBF方法保持一致, 较好地保证黏性网格的正交性; 2) 本文方法能够较好地将物面变形传递至空间网格, 因此远离物面的空间计算网格的网格质量也要好于DGM方法并与RBF方法保持一致, 进一步保证了网格质量。

(a) Global view

2.5 三维HIRENASD算例

以HIRENASD翼身组合体的机翼上弯和扭转变形为例验证本方法在三维复杂外形问题中的变形效率和能力。由于DGM方法的网格变形能力较差, 对于HIRENASD翼身组合体在弯曲小变形情况下网格就已出现负体积, 因此本算例仅与RBF方法进行对比。机翼弯曲和扭转测试变形分别如下

式中, Δz(y)为机翼沿z向的弯曲变形位移沿展向分布,θ(y)为机翼绕前缘的扭转角度沿展向分布,y为机翼展向,b为机翼展长,γ为幅度系数。

该构型的计算域共包含网格节点个数为246 043。图20为原始网格拓扑和翼梢处空间网格细节。

(a) Global view

(a) Minimum quality coefficient

设置弯曲变形系数γw=1和扭转变形系数γt=0.5, 弯曲和扭转变形后的网格分别如图21、图22所示, 弯曲变形会导致附面层网格正交性下降, 图23、图24也表明着弯曲变形对网格正交性的影响明显大于扭转变形。由上述几幅图可知, 本文方法和RBF方法的最小质量系数曲线和平均质量系数曲线基本吻合, 表明本文方法在保持网格质量方面与RBF方法基本一致。

表1为本文方法与RBF方法对HIRENASD翼身组合体的网格变形计算时间对比。为避免误差, 分别采用RBF方法以及本文方法执行50次网格变形程序并对计算时间取平均值, 结果如表1所示, 整体网格计算时间由网格前处理时间、网格变形时间两部分构成, 前处理的时间花费是一次性的, 通过前处理得到的Delaunay和RBF映射关系会存储到文件或内存中, 在后续计算中只需反复调用即可, 而网格变形时间的花费是需要反复叠加的, 因此本方法所花费的网格变形时间要远低于RBF方法。在单次计算中本方法的计算效率提升可能并不显著, 但对于大多数需要多次迭代更新网格的计算而言, 本方法的计算效率要远高于RBF方法。且随着网格数量的不断增加, 计算效率提升的效果会更加明显。