¹浙江理工大学，计算机科学与技术学院(人工智能学院)，浙江 杭州

²中国人民解放军92341部队，河南 洛阳

收稿日期：2023年4月24日；录用日期：2023年7月4日；发布日期：2023年7月12日

摘要

边界层网格的质量是影响计算流体力学中粘性计算精度的关键因素，如何在给定任意复杂模型下生成高质量的边界层网格仍是计算流体力学中的瓶颈。针对边界层网格的质量问题，本文提出一种基于局部特征尺寸场的边界层网格生成算法。该方法将表面网格顶点的最大推进高度作为局部特征尺寸场，并基于局部特征尺寸场构造顶点密度场，在密度场引导下进行保质量与保特征的表面三角网格重网格化，最终在新的表面三角网格上生成边界层网格。实验表明本文方法能在给定任意复杂模型下鲁棒地生成高质量的边界层网格。

关键词

局部特征尺寸场，重网格化，边界层网格

Boundary Layer Mesh Optimization Algorithm Based on Local Feature Size

Linqi Chen¹, Ying Zhou^2*

¹School of Computer Science and Technology (School of Artificial Intelligence), Zhejiang Sci-Tech University, Hangzhou Zhejiang

²92341 Troops of the Chinese People’s Liberation Army, Luoyang Henan

Received: Apr. 24^th, 2023; accepted: Jul. 4^th, 2023; published: Jul. 12^th, 2023

ABSTRACT

The quality of the boundary layer mesh is a important factor affecting the accuracy of viscous calculations in computational fluid dynamics, and how to generate a high quality boundary layer mesh for a given arbitrary complex model is still a bottleneck in computational fluid dynamics. To address the boundary layer mesh quality optimization problem, this paper proposes a boundary layer mesh generation algorithm based on the local feature size. The method takes the maximum advance height of the surface mesh vertices as the local feature size field, and constructs the vertex density field based on the local feature size, then remesh the surface triangular mesh under the guidance of the density field to preserve the quality and features, and finally generates the boundary layer mesh on the new surface triangular mesh. Experiments show that the method in this paper can robustly generate high quality boundary layer meshes for a given arbitrary complex model.

Keywords:Local Feature Size, Remesh, Boundary Layer Mesh

Copyright © 2023 by author(s) and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY 4.0).

http://creativecommons.org/licenses/by/4.0/

1. 引言

近年来，随着算力的发展，计算流体力学在航天航空，汽车等各行各业中扮演着重要的角色。计算模拟仿真的核心是用计算机来求解各种物理微分，积分方程，在具有复杂拓扑和几何的曲面或者实体上求解偏微分方程，通常采用有限元法，有限体积法或者等几何分析方法，而这些方法都强烈依赖于曲面实体的网格生成。边界层网格是一种混合网格，因其能更好地兼顾计算精度和网格易用性 [1] ，所以在计算流体力学中边界层网格被广泛采用。

前沿推进法是应用最为广泛的边界层网格生成方法之一，Lu等人 [2] 提出一种结构多块边界层网格生成方法，该方法首先提取表面网格的几何特征，基于几何特征构造边界层网格的框架，然后在框架内通过TFI技术生成边界层网格。曹杰等人 [3] 提出一种基于离散中轴面的前沿分块层进算法，首先基于四面体单元的外心构建三角网格内部的离散中轴面，再根据表面与中轴面的距离将边界面分为两类分别生成边界层网格。王硕等人 [4] 提出了基于径向基函数的边界层网格生成算法，通过计算表面分区好的网格点的法向构造支撑线，利用径向基函数插值构造网格的轮廓线并生成边界层网格。为解决鲁棒性与网格质量问题，甘洋科等人 [5] 通过收缩推进距离的策略来防止局部自交，并且通过建立空间八叉树来提前检测全局相交问题，再对收缩后的推进方向采用拉普拉斯平滑来改善局部网格单元的质量。刘田田等人 [6] 提出一种基于双前沿推进的边界层网格生成算法，并结合多种网格技术处理局部几何特征以保证边界层网格的质量。Ye等人 [7] 提出基于三种策略的边界层网格生成算法，该算法进一步提高了边界层网格的质量，并且能够鲁棒地生成拓扑合法的边界层网格。

2. 方法

2.1. 问题与方法概述

本文以外流场仿真应用为例，目的是生成面向外流场流体计算的高质量三棱柱边界层网格。令 $S=\langle V_s,F_s\rangle$ 表示一个无自交、可定向的流形三角网格，其中 $V_s$ 与 $F_s$ 分别表示网格的顶点与面片集合。本文在三角网格表面与外流场之间生成指定层数n、初始高度h与增长率r的三棱柱边界层网格单元P，在剩余区域内生成四面体网格单元T。

为鲁棒地生成高质量边界层网格，本文提出基于局部特征尺寸场的边界层网格生成算法，该方法基于前沿推进法生成三棱柱网格单元P，并基于Delaunay方法生成四面体网格单元T。具体而言，为防止产生自交单元，首先计算表面三角网格顶点的最大推进高度，将最大推进高度作为局部特征尺寸场，接着基于局部特征尺寸场与密度函数得到顶点的密度场，对密度场进行拉普拉斯平滑以使密度值分布更为均匀，基于密度场进行表面三角网格的重网格化，最终分别基于前沿推进法与Delaunay方法生成三棱柱网格单元P与四面体网格单元T，得到最终的边界层网格。

2.2. 局部特征尺寸场

为获得局部特征尺寸场，需要计算表面三角网格顶点的推进方向以及最大推进高度。表面三角网格顶点的推进方向可通过平均顶点周围的面法向得到：

$n_p=\frac{1}{N_f}{\displaystyle {\sum }_{i=1}^{N_f}{n}_i}$ (1)

其中 $n_p$ 为推进方向， $N_f$ 为点p相邻的总面数， $n_i$ 表示第i个面的法向。

由指定层数n、初始高度h与增长率r可知顶点的总体推进高度H为：

$H=h*\left(r^0+r^1+\cdots +r^{n-1}\right)$ (2)

由于在三角网格及外流场边界之间的一些狭窄区域，过高的推进高度会造成自交单元产生，因此为了得到合法单元，需要对不合理的推进高度采取缩减策略，如图1所示。图1左侧为自交单元，右侧为经过处理后得到的合法单元。

图1. 缩减过大的推进高度

为检测自交单元，可以对单元构建空间八叉树结构，检测单元之间是否发生碰撞。具体而言，从每个表面三角网格顶点发出一条沿着推进反向的射线，并记录下首次于表面三角网格碰撞时所在的三角面片，取三角面片中距离出发点最近的点并计算两点之间的距离l，并以此调整此顶点的最大推进高度 $H_{\max }$ ：

$H_{\max }=\{\begin{array}{l}H,l>3H\\ H/3,l\le 3H\end{array}$ (3)

为防止三棱柱单元局部的自交问题，需要计算三棱柱单元的雅克比行列式，若为负值，则发生了局部自交，同样需要缩减推进高度以得到合法的三棱柱单元。本文对发生局部自交的三棱柱单元的推进高度以20%的缩减率进行缩减，每缩减一次则检测是否得到合法单元，若依旧不合法，则继续缩减直到单元合法。

由此得到表面三角网格顶点的最大推进高度场，此最大推进高度场实际上就是局部特征尺寸场，接下来需要基于局部特征尺寸场构建顶点密度场。

2.3. 顶点密度场

局部特征尺寸场反映了三棱柱单元的尺寸厚度，若尺寸厚度过低，而底面三角形的尺寸过大，则会导致三棱柱单元纵横比过小。纵横比是一种衡量网格质量的评价指标，对于三棱柱单元，其计算方式为三棱柱平均厚度与底面三角形平均边长之比，其值越接近1则三棱柱质量越好。因此为了提升棱柱单元的质量，需要通过构建密度场来增加底面三角形的密度，即对表面三角网格进行重网格化。因此，本文计算每个顶点的密度值的方式为：

$\text{Density}\left(v_i\right)=\{\begin{array}{l}10/L,0<H_{\max }\left(v_i\right)<L/10\\ 1/H_{\max }\left(v_i\right),L/10\le H_{\max }\left(v_i\right)\le L\\ 1/L,H_{\max }\left(v_i\right)>L\end{array}$ (4)

其中， $\text{Density}\left(v_i\right)$ 代表顶点 $v_i$ 的密度值， $H_{\max }\left(v_i\right)$ 代表顶点 $v_i$ 的最大推进高度，即局部特征尺寸。L代表输入的表面三角网格的目标平均边长，即重网格化后得到的表面三角网格的平均边长。由平均边长可以推导出重网格化后的目标点数，重网格化后的三角网格应是各向同性的，由各向同性的三角网格可知每个三角形为正三角形，那么每个三角形的面积a为：

$a=\frac{\sqrt{3}}{4}{L}^2$ (5)

由模型的表面积A与每个三角形的面积a可得知总三角形个数 $n_{tri}$ 为：

$n_{tri}=\frac{A}{a}$ (6)

由欧拉公式可知三角形个数为顶点个数的两倍，则 $v_{tri}$ 为：

$v_{tri}=\frac{n_{tri}}{2}$ (7)

若某一点与其邻接点的密度值差距过大，则会造成密度值突变，那么可能导致重网格化后的三角形数量过渡不平滑，会对后续边界层网格生成造成影响，所以需要对密度场进行拉普拉斯平滑。那么对于每个顶点 $v_i$ 的密度值，其新的密度值计算方式为：

$\text{Density}{\left(v_i\right)}^{\prime }=\text{Density}\left(v_i\right)-\frac{1}{n_{v_i}}{\displaystyle {\sum }_j^{n_{v_i}}\text{Density}\left(v_j\right)}$ (8)

其中， $\text{Density}{\left(v_i\right)}^{\prime }$ 代表每个顶点 $v_i$ 经过平滑后的密度值， $\text{Density}\left(v_i\right)$ 代表初始密度值，n代表顶点 $v_i$ 的1-邻域所有邻接点的数量。由此密度场构建完毕，接下来需要基于密度场对三角网格进行重网格化。

2.4. 重网格化

重网格化的目标是对原始三角网格进行重采样以符合密度场要求。本文基于Shen等人 [8] 提出的重网格化算法进行重网格化，该方法需要输入重网格化后的目标点数，密度场以及三角网格的特征线，由3.1节中可知本文已计算得到目标点数和密度场，接下来需要提取三角网格的特征线。本文通过检测相邻三角面片的二面角值来进行特征线的提取，即某条边邻接的两个三角面片的法线夹角大于指定阈值θ，那么将其视为特征边：

$n_1\cdot n_2\le \Vert n_1\Vert \cdot \Vert n_2\Vert \cdot \cos \theta$ (9)

其中， $n_1$ 和 $n_2$ 分别为两邻接三角面片的单位法向，如图2所示。图2左侧为输入三角网格，右侧为提取出来的特征线。

图2. 提取特征线

由此，基于目标点数，密度场以及三角网格的特征线对三角网格进行重网格化后，可得到符合密度场要求的三角网格，接下来需要生成边界层网格。

2.5. 生成边界层网格

本文基于前沿推进法生成边界层网格的三棱柱单元，即从表面三角网格顶点出发，沿着顶点的推进方向推进得到新的顶点，并根据表面三角网格的拓扑关系连接新的顶点得到三棱柱单元。由层数，初始高度以及增长率可得到多层三棱柱单元，如图3所示。

图3. 前沿推进法

本文为得到合法无自交的三棱柱单元，基于局部特征尺寸场进行三棱柱网格的生成，即在保证单元无自交的情况下进行推进高度的调整，并基于局部特征尺寸场得到密度场，在密度场引导下进行表面三角网格的重网格化，最终基于顶点最大推进高度沿着推进方向得到三棱柱单元。

生成得到三棱柱单元后，需要在外流场区域生成四面体网格单元。生成四面体网格的方法有很多 [9] [10] ，本文借助TetGen [10] 开源库生成高效鲁棒的四面体网格单元。其生成算法基于Delaunay算法，从理论上保证了四面体网格的质量。

3. 实验结果与分析

为测试以及评估本文方法，本文通过实验分析得到的结果，实验所用操作系统为windows11，处理器为2.9 GHz IntelCorei5-10400 CPU，内存为16 GB。本文实现该算法为一个单线程的C++程序，其涉及到矩阵运算库，由Eigen [11] 库实现。本文在若干个网格上实现了该算法，并通过一系列质量指标评估结果。

3.1. 边界层网格生成结果

图4(a)，图4(b)分别为airplane模型，scaffold模型与joint模型生成得到的边界层网格，从图中可以看出，各模型生成的边界层网格单元的过渡性较好，具有良好的正交性，三棱柱单元较为规整，四面体单元大部分接近正四面体。表1列出了各个模型的边界层网格的质量评价指标，其中，纵横比指标越接近1越好，扭曲度越小越好，雅克比值越接近1越好。由表1可见，本文生成得到的边界层网格的质量较好。

图4. 边界层网格

表1. 边界层网格质量评价指标

3.2. 重网格化结果分析

本文的重网格化旨在调整表面三角网格的密度以符合密度场要求，目的是为了改善复杂模型狭窄区域的三棱柱网格单元的质量。为了展示重网格化带来的提升，本文对测试模型构造狭窄区域，位于模型与外流场之间。

图5~图7分别展示了assy模型、aries模型与bendix模型重网格化前后得到的边界层网格的切割视

图5. assy模型重网格化前后生成的边界层网格切割视图

图6. aries模型重网格化前后生成的边界层网格切割视图

图7. bendix模型重网格化前后生成的边界层网格切割视图

表2. 边界层网格质量评价指标

图。从图中可以看出，重网格化前模型的狭窄区域的三棱柱单元非常扁平，其纵横比过小，经过重网格化后，狭窄区域的表面三角形的密度增加，因而提升了三棱柱单元的纵横比，从而提高了边界层网格的质量。表2列出了各个模型重网格化前后的质量指标比较。由表2可见，重网格化后能有效提升边界层网格的质量。

4. 总结

本文主要针对计算流体力学中的边界层网格生成问题提出了一种基于局部特征尺寸场的边界层网格生成算法。该算法首先计算三角网格顶点的最大推进高度作为局部特征尺寸场，并基于局部特征尺寸场构造表面三角网格顶点的密度场，为使得密度值过渡平滑，采用了拉普拉斯平滑方法平滑密度场，接着基于密度场对表面三角网格进行重网格化，由于重网格化算法需要提取三角网格的特征线，因此基于二面角的方法对三角网格进行了特征线的提取，最后基于前沿推进法得到三棱柱网格单元，并在外流场区域基于Delaunay法生成四面体网格单元，得到最终的边界层网格。实验结果表明，本文提出的方法能鲁棒地生成边界层网格，并且边界层网格单元拓扑合法，不会有自交单元的产生，此外，结合重网格化方法对三角网格进行密度调整，在新的三角网格上生成边界层网格，能有效提高复杂模型的狭窄区域处的网格单元质量。

文章引用

陈林祺,周 颖. 基于局部特征尺寸场的边界层网格生成算法
Boundary Layer Mesh Optimization Algorithm Based on Local Feature Size[J]. 建模与仿真, 2023, 12(04): 3464-3471. https://doi.org/10.12677/MOS.2023.124317

参考文献