网格生长动力学模型的建立和求解及其混合黏性网格生成中的应用1)

任一鹏 刘 枫 林 其 胡华伟

∗(北京宇航系统工程研究所，北京100076)

†(中国空气动力研究与发展中心，四川绵阳621000)

∗∗(成都铭峰新源科技有限公司，成都610042)

真实飞行器外形复杂，所涉及的流动结构和现象对数值模拟中所需的计算网格提出了较高的要求。一方面，飞行器本身构型多样，各种部件无论从尺度还是位置都存在着多种变化和组合，将传统的参数化的网格生成方法应用其中存在重大障碍；另一方面，真实飞行条件多变，所产生的流动结构性质各异，对计算网格的要求也不尽相同。在流动的数值模拟中，一个熟知的原则就是网格分布应该尽可能符合当地的流动状态。因此，建立能够适应复杂外形的网格生成方法是必要的。

国内外一些学者在复杂外形的网格生成方法[1-2]上进行了诸多探索。在过去的二十年中，许多重要的网格生成技术得到建立和发展，比如多块对接/拼接结构网格、重叠/嵌套网格、非结构(四面体/金字塔)网格等。近年来混合网格(包含有多种形态的网格)越来越受到关注，得到了较快的发展。Zhang等[3-4]建立了基于各向异性四面体凝聚的混合网格生成方法，并成功运用于大雷诺数的复杂飞行器数值模拟中；Esquieu[5]和Kallinderis等[6]建立了三棱柱/四面体混合网格生成算法，并运用于喷气式客机绕流流场的计算；Timothy[7]建立了三棱柱/金字塔/四面体/六面体混合网格生成方法；Kannan等[8]和Wang[9]建立了自适应笛卡尔网格生成算法。众所周知，边界层网格的质量对计算结果有重大影响。三棱柱网格作为一种各向异性网格可以任意沿其法向拉伸或者压缩而不影响其形状，能够满足边界层计算的需要。然而，传统的三棱柱网格生成算法(层推进算法及其改进[10-11])在凹凸相邻或相交区域都遇到了较大的困难，因此难以运用于复杂三维外形的网格生成中。本文试图建立网格生长动力学模型，通过迭代求解动态优化调整网格点生长法向和长度，形成边界层网格生成算法，然后结合各向同性四面体网格的阵面推进Delaunay网格细化算法，建立一种新型的混合黏性网格生成方法。

1 网格的类型和性质

根据网格的几何特征，可以将网格类型分为各向同性网格和各向异性网格，如图1所示。众所周知，不同流动结构对网格尺度的要求是不同的，在流动变化剧烈的方向需要加密网格。因此，各向同性网格适合布置在流动梯度较小、物理量分布较光滑的区域；各向异性网格由于网格几何特征决定了网格方向性的差异，不同几何方向可以存在较大的尺寸比，可以用于流动梯度较大、物理量分布存在间断的区域。

图1 网格类型图

2 网格生长动力学模型

2.1 网格生长受力分析

如图2所示，根据生长边界的几何特性的不同，一般可以分为凸壳生长、非凸壳生长以及二者的组合。这三种生长类型都试图将网格的法向“光顺化”，由法向变化剧烈的最底层网格逐渐过渡到法向光滑的顶层网格。这种“光顺化”在二维表现为圆弧，在三维情况表现为球面。

图2 网格生长示意图

假设网格点之间存在力的相互作用，并进一步假设作用力满足胡克定律，即

式中，k为弹性系数，为网格点之间的矢径。

不失一般性，分别考虑凸壳和非凸壳中法向变化最剧烈的角点受力情况，如图3和图4所示。显然，在最底层网格中，除角点外的其他网格点所受合力为0。凸壳中角点受力沿着网格生长的方向，合力使得网格生长加速；非凸壳中角点受力沿着网格生长相反的方向，合力会抑制网格生长。这就使得凸壳中的角点在光顺化条件下，生长总距离较其他网格点长，网格生长速度较快；在非凸壳中反之。一般来说，网格生成成败的关键在于“角点”的处理，网格交叉、网格扭曲最严重的区域也在于此。

图3 凸壳角点受力分析

图4 非凸壳角点受力分析

2.2 网格生长动力学模型

如图5所示，建立O−xyz笛卡尔网格坐标系。由上节分析可知，对于无论凸壳或是非凸壳情况，网格点上的合力是网格生长的动力，这种合力作用随着网格生长而改变。由牛顿第二定律，即有

图5 网格坐标系

不妨假设在生长过程中网格点作用合力不变，对式(2)积分，显然有

其中，c0，c1为常数。随着时间推进，|r|呈抛物线快速增长或者减小。这种生长模式简单直接，然而由于网格尺度增加或缩减速度过快，对于复杂外形极易出现网格交叉和负体积，导致网格生成失败。根据对图1和图2的分析，如果引入某种阻尼，对网格生长速度进行动态抑制，则有

式中，Kv为阻尼函数。初始条件为

式中，r0为初始位置向量，n0为初始法向量。

对整个网格生长阵面而言，则有

式(6)是一个N维二阶常微分非线性系统。

2.3 模型求解

动力学系统公式(6)存在强烈的非线性且相互耦合，难以直接求得解析解。然而，我们知道在每一层网格阵面进行推进，特别是生成边界层网格时，由于边界层尺度η远小于飞行器的特征尺度L，不妨假设在每一层网格推进过程中，各网格点生长过程相互独立，则有

式(7)是一个二阶非齐次常微分方程。不失一般性，去掉方程中的下标i，式(7)对应的齐次方程为

式(8)对应的特征方程是

方程(9)有两个不相等的实数解M1,2=0,−Kv，则齐次方程(8)的通解形式为

不难找出

是方程(7)的一个特解。根据常微分方程解析理论，非齐次方程的通解为对应齐次方程的通解与非齐次方程的一个特解之和，则方程(7)的通解为

代入初始条件(5)，则有

将方程(13)代入方程(12)中，则方程(7)的解为

将初始条件式(5)代入式(3)中，则有

对比式(14)和式(15)，不难看出两种不同的生长模式对网格生长高度和方向的影响。下面分别考虑凸壳和非凸壳中网格生长随时间的变化模式。在凸壳(如图3)中，令则图6和图7分别给出了凸壳中模式一(阻尼模式)和模式二(非阻尼模式)中网格生长高度随时间的变化。在非凸壳(如图4)中，令则图8和图9分别给出了非凸壳中模式一(阻尼模式)和模式二(非阻尼模式)中网格生长高度随时间变化。可以看出，阻尼模式会抑制网格生长高度的过快增长，从而增强网格生成的鲁棒性。

图6 凸壳角点网格生长高度r x随时间变化

图7 凸壳角点网格生长高度r y随时间变化

图8 非凸壳角点网格生长高度r x随时间变化

图9 非凸壳角点网格生长高度r y随时间变化

2.4 法向量求解与阻尼函数选取

法向量存在于网格推进阵面的面心处，而网格格点的法向量有

S i为格点p相邻网格阵面的面积，n i F为相邻网格阵面的法向量。然而，网格格点法向量与周围格点法向量密切相关且相互作用，本文取周围格点法向量的平均值，则有

式(16)等价于线性系统

其中

线性系统(17)对应的放大因子显然小于1，因此该线性系统可以通过Jacobi迭代或者Gauss–Seidel迭代求解。迭代求解初值按照式(16)求得。

根据式(4)的定义，阻尼函数的选取应该与网格生长距离l的模成反比，即

3 各向异性网格生成算法

上节建立和求解了网格生长动力学模型，运用该模型可以得到各向异性网格生成算法，归结起来可以分为以下几步：

步骤1：计算每层网格的预计生长高度。每层网格预计生长高度按指数增加计算，即

式中，α为层增长率，β为预设网格层数，h0为第一层网格高度。

步骤2：计算网格推进阵面上各个网格点的法向矢量。根据式(16)计算法向矢量初值，求解式(17)得到各个网格点的法向矢量。

步骤3：计算当前时间步中网格推进阵面上各个网格点的阻尼函数Kv和所受合力F。分别根据式(1)和式(19)计算合力F和阻尼函数Kv。

步骤4：根据式(14)计算每个时间步中生长矢量r。

步骤5：判断当前推进阵面中任意网格点的生长高度是否达到该网格点的预计生长高度。如果是，则转入下一层网格生长；如果否，则返回步骤2继续时间推进，直到当前推进阵面中任意网格点的生长高度达到该网格点的预计生长高度。

4 各向同性四面体生成算法

各向同性四面体生成算法，国内外已经有众多学者进行了广泛而深入的研究，已经相对成熟。各类算法的主要区别在于网格生成质量和效率的差异。阵面推进法的特点是搜索量大，每次搜索只能完成一个网格的生成，网格生成效率较低，生成网格不一定能满足外接球准则即Delaunay准则；而Delaunay方法因其生成网格的效率较高，而且网格质量高，天然满足Delaunay准则，然而需要高质量的边界恢复算法才能保形。本文基于计算几何的基本原理，建立了基于可视面的快速四面体初始化算法，然后建立了阵面推进Delaunay网格细化算法，并结合网格优化算法，发展了一套任意多面体的快速非结构网格生成算法。该算法一方面克服了阵面推进法网格生成效率低的缺点，另一方面克服了Delaunay方法在边界恢复中不保形的缺陷。详细算法描述及定义参考文献[12-13]。

其主要步骤包括：

步骤1：通过建立的基于可视面的初始化四面体算法，将网格剖分区域进行初始化，形成初始四面体阵列；

步骤2：然后通过运用改进的最差面细化算法将网格细化，其间运用Delaunay面交换算法优化网格质量；

步骤3：求解基于顶点弹簧的线性系统进一步提高网格质量。

5 算例验证

5.1 双椭球算例

双椭球算例是高超声速风洞试验的标准模型，外形由两个不同形状的椭球叠加而成。图10给出了双椭球外形的表面网格，共由14102个三角形构成。三棱柱网格共10层，网格量为141020。图11～图13分别给出了头部、拐角和底部转角处的局部网格，可以看出网格无交叉且光顺性较好。

图10 表面网格

图11 头部网格

图13 底部转角网格

5.2 空天飞机算例

空天飞机算例是翼身融合体模型，外形中有数个凹凸部件相连的区域，对网格生成算法的鲁棒性有很高的要求。图14给出了空天飞机外形的表面网格，共由80606个三角形构成。三棱柱网格共2层，网格量为1612120。图15～图17分别给出了不同截面的局部网格，可以看出网格光顺性较好。

图12 拐角网格

图14 表面网格

图15 截面x/L=0.25网格

图17 截面y/L=0.2网格

图16 截面x/L=0.9网格

6 结论

(1)本文提出并求解的网格生长动力学模型，可以用于黏性网格的边界层网格生成中，具有较高的可靠性；

(2)基于各向异性网格生成算法和各向同性四面体网格生成算法的网格生成策略，可以有效地运用于飞行器的网格生成中。