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

第3期
任一鹏等：网格生长动力学模型的建立和求解及其混合黏性网格生成中的应用
407
合当地的流动状态。因此，建立能够适应复杂外形的 网格生成方法是必要的。
国内外一些学者在复杂外形的网格生成方法【1-2] 上进行了诸多探索。在过去的二十年中，许多重要 的网格生成技术得到建立和发展，比如多块对接/拼 接结构网格、重叠/嵌套网格、非结构（四面体/金 字塔）网格等。近年来混合网格（包含有多种形态的 网格）越来越受到关注，得到了较快的发展。Zhang 等GT建立了基于各向异性四面体凝聚的混合网格 生成方法，并成功运用于大雷诺数的复杂飞行器数 值模拟中；Esquieu同和Kallinderis等冏建立了三 棱柱/四面体混合网格生成算法，并运用于喷气式客 机绕流流场的计算；Timothy 1?1建立了三棱柱/金字 塔/四面体/六面体混合网格生成方法；Kannan等冈 和Wangle建立了自适应笛卡尔网格生成算法。众 所周知，边界层网格的质量对计算结果有重大影响。 三棱柱网格作为一种各向异性网格可以任意沿其法 向拉伸或者压缩而不影响其形状，能够满足边界层 计算的需要。然而，传统的三棱柱网格生成算法（层 推进算法及其改进I10-11））在凹凸相邻或相交区域都 遇到了较大的困难，因此难以运用于复杂三维外形 的网格生成中。本文试图建立网格生长动力学模型, 通过迭代求解动态优化调整网格点生长法向和长度, 形成边界层网格生成算法，然后结合各向同性四面 体网格的阵面推进Delaunay网格细化算法，建立一 种新型的混合黏性网格生成方法。
Key words hybrid grid, viscous grid, step self-adaptive, complex geometry, mesh growing dynamics model
真实飞行器外形复杂，所涉及的流动结构和现 象对数值模拟中所需的计算网格提出了较高的要求 。 一方面，飞行器本身构型多样，各种部件无论从尺度 还是位置都存在着多种变化和组合，将传统的参数
第43卷第3期
力学与实践
2021年6月
网格生长动力学模型的建立和求解及其混合黏性 网格生成中的应用"
任一鹏*卫)刘枫十林其t 胡华伟**
*(北京宇航系统工程研究所，北京100076) 卞(中国空气动力研究与发展中心，四川绵阳621000)
**(成都铭峰新源科技有限公司，成都610042)
摘要 真实飞行器外形复杂，涉及的流动结构和现象对数值模拟中所需的网格提出了较高的要求。本文首 先研究了各向异性黏性网格的生成方法，提出了网格生长动力学模型，并建立了松弛算法求解该模型，从而得 到合理的层状网格；提出并建立了各向同性四面体网格的阵面推进Delaunay网格细化算法。最后通过数个算 例对所建立的混合黏性网格生成方法的适用性和可靠性进行考察。结果表明，本文建立的混合黏性网格生成方 法适用于复杂飞行器的网格生成，具有较高的可靠性。
Ren Yipeng, Liu Feng, Lin Qi, et al. Mesh growing dynamics model for hybrid viscous mesh generation method for complex geometry. Mechanics in Engineering, 2021, 43(3): 406-413
iSN(P)
图4非凸壳角点受力分析
2.2网格生长动力学模型 如图5所示，建立O-zyz笛卡尔网格坐标系。
由上节分析可知，对于无论凸壳或是非凸壳情况，网 格点上的合力是网格生长的动力，这种合力作用随 着网格生长而改变。由牛顿第二定律，即有
式中，Kv为阻尼函数。初始条件为
r\t=o = ro
dr
⑸
dt
=n0
1网格的类型和性质
厂各向异性网格 计算网格乂
结构网格 三棱柱网格 笛卡尔网格 六面体网格 四面体网格
各向同性网格
金字塔网格 图1网格类型图
（a）凸壳生长示意图
（b）非凸壳生长示意图 图2网格生长示意图
合。这三种生长类型都试图将网格的法向“光顺化”, 由法向变化剧烈的最底层网格逐渐过渡到法向光滑 的顶层网格。这种“光顺化”在二维表现为圆弧，在 三维情况表现为球面。
假设网格点之间存在力的相互作用，并进一步 假设作用力满足胡克定律，即
根据网格的几何特征，可以将网格类型分为各 向同性网格和各向异性网格，如图1所示。众所周 知，不同流动结构对网格尺度的要求是不同的，在流 动变化剧烈的方向需要加密网格。因此，各向同性网 格适合布置在流动梯度较小、物理量分布较光滑的 区域；各向异性网格由于网格几何特征决定了网格 方向性的差异，不同几何方向可以存在较大的尺寸 比，可以用于流动梯度较大、物理量分布存在间断的 区域。
ቤተ መጻሕፍቲ ባይዱt=o
式中，To为初始位置向量，«0为初始法向量。 对整个网格生长阵面而言，则有
ri
Kvi …• • 0
d2 『2
F2
=
-
dt2
!凤2
_ rn _ _Fn_
0 •…• ■ Kvn, _
rn
式(6)是一个N维二阶常微分非线性系统。
2.3模型求解 动力学系统公式(6)存在强烈的非线性且相互
耦合，难以直接求得解析解。然而，我们知道在每一 层网格阵面进行推进，特别是生成边界层网格时，由 于边界层尺度T］远小于飞行器的特征尺度L,不妨 假设在每一层网格推进过程中，各网格点生长过程 相互独立，则有
化的网格生成方法应用其中存在重大障碍；另一方 面，真实飞行条件多变，所产生的流动结构性质各 异，对计算网格的要求也不尽相同。在流动的数值模 拟中，一个熟知的原则就是网格分布应该尽可能符
本文于2021-01-11收到。 1) 国家自然科学基金资助项目(91016010). 2) 任…鹏，高级工程师，研究方向为飞行器气动与热环境总体设计。E-mail: rtianxiao@ 引用格式：任一鹏，刘枫，林其等.网格生长动力学模型的建立和求解及其混合黏性网格生成中的应用.力学与实践，2021, 43(3): 406-413
关键词 混合网格，黏性网格，步长自适应，复杂外形，网格生长动力学模型
中图分类号：V211.3 文献标识码：A doi: 10.6052/1000-0879-21-017
MESH GROWING DYNAMICS MODEL FOR HYBRID VISCOUS MESH GENERATION METHOD FOR COMPLEX GEOMETRY*1)
|r| = cie_Kvt + C2
(10)
不难找出 F (11)
是方程(7)的一个特解。根据常微分方程解析理论, 非齐次方程的通解为对应齐次方程的通解与非齐次 方程的一个特解之和，则方程(7)的通解为
= cie*t + p + C2 (12)
图6凸壳角点网格生长高度rx随时间变化
代入初始条件(5),则有
式(7)是一个二阶非齐次常微分方程。不失一般性, 去掉方程中的下标式(7)对应的齐次方程为
第3期
任一鹏等：网格生长动力学模型的建立和求解及其混合黏性网格生成中的应用
409
式(8)对应的特征方程是
M2 + KVM = 0
⑼
方程(9)有两个不相等的实数解Mi,2 = 0,-Kv，则 齐次方程(8)的通解形式为
REN Yipeng*，2) LIU Feng+ LIN Qi+ HU Huawei**
* (Beijing Institute of Astronautical System Engineering, Beijing 100076, China) t (China Aerodynamic Research and Development Center, Mianyang 621000, Sichuan, China) ** (Chengdu Mingfeng Newprotech Technology Limited Company, Chengdu 610042, China)
图7凸壳角点网格生长高度Ty随时间变化 图8非凸壳角点网格生长高度rx随时间变化
410
力学与实践
2021年第43卷
其中
"_ ( 0,其他
线性系统(17)对应的放大因子显然小于1,因此 该线性系统可以通过Jacobi迭代或者Gauss-Seidel 迭代求解。迭代求解初值按照式(16)求得。
根据式(4)的定义，阻尼函数的选取应该与网格 生长距离I的模成反比，即
408
力学与实践
2021年第43卷
不妨假设在生长过程中网格点作用合力不变，对 式⑵积分，显然有
卜| = co + + ㊁|戸|产
(3)
其中，co, Cl为常数。随着时间推进，|r|呈抛物线快 速增长或者减小。这种生长模式简单直接，然而由于 网格尺度增加或缩减速度过快，对于复杂外形极易 出现网格交叉和负体积，导致网格生成失败。根据对 图1和图2的分析，如果引入某种阻尼，对网格生 长速度进行动态抑制，则有
2网格生长动力学模型
2.1网格生长受力分析 如图2所示，根据生长边界的几何特性的不同，
一般可以分为凸壳生长、非凸壳生长以及二者的组
耳j = krij
（1）
式中，k为弹性系数，fc= l/||rv||2； rij为网格点之 间的矢径。
不失一般性，分别考虑凸壳和非凸壳中法向变 化最剧烈的角点受力情况，如图3和图4所示。显 然，在最底层网格中，除角点外的其他网格点所受合 力为0。凸壳中角点受力沿着网格生长的方向，合力 使得网格生长加速；非凸壳中角点受力沿着网格生 长相反的方向，合力会抑制网格生长。这就使得凸壳 中的角点在光顺化条件下，生长总距离较其他网格 点长，网格生长速度较快：在非凸壳中反之。一般来 说，网格生成成败的关键在于“角点”的处理，网格 交叉、网格扭曲最严重的区域也在于此。
Abstract A real aircraft has a complex geometry, in its simulation, a high quality mesh is required for its structure and the related flow field. In this paper, first we discuss the grid generation method for the anisotropic viscous mesh and propose a grid growing-iteration model, with a relax algorithm to generate a reasonable layered grid. We also propose an advancing front-Delaunay (AF-D) refine algorithm for the isotropic tetrahedron grid. Finally we study a few cases to see algorithm's feasibility and reliability. The results indicate that the method of the hybrid viscous grid generation for complex geometry is reliable.
将方程(13)代入方程(12)中，则方程⑺的解为
将初始条件式(5)代入式(3)中，则有
r = r0 + 710上+*戸产
(15)
对比式(14)和式(15),不难看出两种不同的生 长模式对网格生长高度和方向的影响。下面分别考 虑凸壳和非凸壳中网格生长随时间的变化模式。在 凸壳(如图 3)中，令 r0 = 0, n0 = (-^/2,^/2/2), 则 F = (-^2/2, v/2/2), Kv = ||r||2o 图 6 和图 7 分别给出了凸壳中模式一(阻尼模式)和模式二(非 阻尼模式)中网格生长高度随时间的变化。在非凸 壳(如图 4)中，令 r0 = 0, n0 = (^/2,^2/2),则 F = (-^/2,-^2/2), Kv = ||r||2» 图 8 和图 9 分别 给出了非凸壳中模式一(阻尼模式)和模式二(非阻 尼模式)中网格生长高度随时间变化。可以看出，阻 尼模式会抑制网格生长高度的过快增长，从而增强 网格生成的鲁棒性。