CN112182489B - 一种基于偏微分方程求解的二维高阶网格生成方法 - Google Patents

Abstract

本发明公开了一种基于偏微分方程求解的二维高阶网格生成方法。当前，生成高质量、高精度的高阶网格仍是一个开放性难题。本发明首先生成覆盖飞机翼形二维几何模型的线性网格，建立线性网格边界和飞机翼形二维几何模型边界的映射关系；然后，将线性网格进行升阶，以升阶线性网格边界偏移矢量为约束，建立控制内部网格点偏移矢量场的偏微分方程，并采用高阶边界元数值方法求解该偏微分方程；最后，将升阶线性网格的所有网格点沿各自对应的偏移矢量进行偏移，偏移后的每个线单元用高阶几何形式表示，得到高阶网格。本发明有效提高了获得偏移矢量场的精度和速度，实现了基于飞机翼形模型快速生成网格数目尽量少的高质量的高阶网格单元。

Description

一种基于偏微分方程求解的二维高阶网格生成方法

技术领域

本发明涉及数值模拟领域前处理中的高精度网格生成过程，具体涉及一种基于PDE求解的二维高阶网格生成方法。

背景技术

网格生成是有限元法、有限体积法和有限差分法等数值模拟技术中的前处理过程，该过程将连续的几何区域剖分成有限个基本几何形体的组合，这些基本几何形体被称为网格单元。常用的单元类型有三角形单元、四边形单元、四面体单元和六面体单元等。网格单元的数目和质量对数值计算的精度和效率有直接影响。好的网格需要用尽量小的网格规模和自由度获得尽量高的数值计算精度。为了满足高精度复杂模型仿真需求，工业上往往采用大规模网格生成，导致时间消耗和计算资源开销非常巨大，几何精度上仍有部分损失。对这类问题，比较好的选择是采用以高阶网格作为输入的高阶数值方法进行数值模拟。

相对于低阶数值方法，高阶数值方法能有效减小计算误差、提升计算精度。另外，在相同计算精度情况下，其对网格规模的要求低于低阶数值方法，从而提升计算效率。为充分发挥高阶数值方法的精度优势，需要将低阶数值方法中应用的直边网格替换为曲边网格，并使曲边网格在几何边界上逼近几何，对应高阶数值方法，这类网格称为高阶网格，网格单元称为高阶单元。

然而，高阶网格的生成却并非易事。在高阶网格中，高阶单元在几何上需由二次、三次甚至更高次曲线、曲面来表示，这无疑增加了其生成难度。此外，相比于线性网格，生成的高阶网格中很容易出现扭曲单元和相交单元等低质量网格单元，进而影响后期数值分析。当前，生成高质量、高精度的高阶网格仍是一个开放性难题。现有的大多高阶网格生成方法依赖于全局迭代求解，需要消耗较多的时间资源和计算资源，难以适应于工业中的大规模高阶网格生成需求。

发明内容

针对二维几何模型的高阶网格生成需求，本发明提供一种基于偏微分方程求解的二维高阶网格生成方法，根据飞机翼形二维几何模型快速生成网格数目尽量少的高质量的高阶网格单元，减少时间和计算资源。

本发明采用的技术方案是：输入飞机翼形二维几何模型，生成覆盖该几何模型的线性网格，并建立线性网格边界和飞机翼形二维几何模型边界的一一映射关系；然后，将线性网格进行升阶，以升阶线性网格边界偏移矢量为约束，建立控制内部网格点偏移矢量场的偏微分方程，并采用高阶边界元数值方法求解该偏微分方程；最后，将升阶线性网格的所有网格点沿各自对应的偏移矢量进行偏移，偏移后的每个线单元用高阶几何形式表示，得到高阶网格。具体包括以下步骤：

步骤1、根据飞机翼形二维几何模型生成覆盖飞机翼形二维几何模型的线性网格，建立线性网格边界到飞机翼形二维几何模型边界的映射关系，具体如下：

输入飞机翼形二维几何模型，(采用现有低阶网格生成方法)生成覆盖飞机翼形二维几何模型的较为稀疏的线性三角形网格，记线性三角形网格边界处的线单元为ΔL_j,j∈[1,N],N为线单元数；记飞机翼形二维几何模型边界为将线性三角形网格边界处每个线单元上的网格点投影到最靠近该线单元的飞机翼形二维几何模型边界上，从而将飞机翼形二维几何模型边界划分成N段，记为几何边界段/>j∈[1,N]；建立线单元ΔL_j到几何边界段/>的映射/>

步骤2、对覆盖飞机翼形二维几何模型的线性网格中每个网格单元进行升阶，得到升阶三角形网格，具体如下：

在每个网格单元的线单元上以及内部均添加额外的网格点，并通过插值变换将每个线单元用高阶几何形式表示，得到升阶三角形网格；插值变换如下：设以m次函数N_α(ε)作为插值基函数，且每个线单元均有m+1个网格点，则线单元的高阶几何形式为：

式中，为第α个网格点在i方向上的值；p为线单元上网格点通式表达，i表示二维上的x或y分量，每个线单元有m+1个网格点，N_α(ε)表示第α个网格点对应的插值基函数，ε为插值基函数的参数，取值范围是[-1,1]。

步骤3、建立偏移矢量场最优化模型和约束条件，具体如下：

以升阶三角形网格边界的线单元偏移矢量向升阶三角形网格内部推进时满足最小变化为目标，建立偏移矢量场最优化模型：

约束条件为：

其中，Ω为升阶三角形网格的内部区域，u_i(p)为位于升阶三角形网格上的网格点p(p为网格点表达参量)处的偏移矢量在i方向的分量，为u_i(p)梯度的2范数，和/>分别为由升阶三角形网格边界得到的狄利克雷边界和诺依曼边界，n为边界上网格点p处的外法向量，u_i,0(p)为升阶三角形网格边界的网格点p到飞机翼形二维几何模型边界的偏移矢量沿i方向的分量，/>为u_i,0(p)在n方向上的偏导。

偏移矢量场最优化模型的最优解等价为如下拉普拉斯边值问题的解：

式(4)中，Δu_i(p)为u_i(p)的拉普拉斯算子。

通过线单元ΔL_j到几何边界段的映射关系，并采用步骤二中的插值变换将升阶三角形网格边界的每个线单元偏移矢量的约束条件用高阶几何形式表示，得到：

p^α表示第α个网格点；

步骤4、求解得到升阶三角形网格的偏移矢量场，将升阶三角形网格的所有网格点沿各自对应的偏移矢量进行偏移，偏移后的每个线单元用高阶几何形式表示，得到高阶网格，具体过程如下：

根据公式(4)，升阶三角形网格的任意一个内部点p₁处的偏移矢量沿i方向的分量通过如下边界积分公式求解：

其中，ΔL为线性三角形网格边界处线单元的集合，点p₂为位于升阶三角形网格边界线单元上的点的一种参量表达,表示点p₂邻域的微分，G(p₂,p₁)为Δu_i(p)＝0方程的基本解，即/>n2为点p₂处的外法向量。

同时，由公式(4)，当升阶三角形网格边界线单元上每个点的u_i(p₂)和中有一个值已知时(u_i(p₂)可由公式(5)求得)，另一个值通过如下积分公式得到：

其中，点p₃为升阶三角形网格边界线单元上的点的另一种参量表达，假设升阶三角形网格边界的各个线单元之间光滑连接，则常量 表示点p₃邻域的微分，/>n3为点p₃处的外法向量。公式(7)的求解采用高阶边界元数值方法，结合公式(5)，得到如下离散形式：

其中，q_k1为位于线单元ΔL_k1上的点，k1∈[1,N]，q_j为位于ΔL_j上的点，nq_j为q_j处的外法向量，表示点q_j邻域的微分，/>为第j个边界线单元上的第α个网格点。当q_k1和q_j重合或处于同一个线单元时，采用高斯积分求解。公式(8)表示一个由N*m个未知数和N*m个方程组成的方程组。通过求解该方程组可得到位于升阶三角形网格边界线单元上的边值u_i(p)或/>随后，通过公式(6)求解位于升阶三角形网格内部的网格点偏移矢量，从而获得整个升阶三角形网格的偏移矢量场。

将升阶三角形网格的所有网格点沿各自对应的偏移矢量进行偏移，并采用步骤二中的插值变换将升阶三角形网格的所有网格点偏移后的每个线单元用高阶几何形式表示，得到高阶网格。

本发明的有益效果是：

本发明通过建立覆盖飞机翼形二维几何模型区域的高阶偏移矢量场物理模型，并采用高阶边界元法进行数值求解，有效提高获得偏移矢量场的精度和速度，再根据偏移矢量场更新所有网格点位置，插值拟合生成高阶网格，实现了基于飞机翼形模型快速生成网格数目尽量少的高质量的高阶网格单元。本发明面向高精度数值模拟需求，可以分析二维翼形几何结构，快速计算生成对应的高阶网格，也适用其他复杂二维几何模型。

附图说明

图1为本发明中飞机翼形二维几何模型的边界示意图；

图2为本发明中覆盖飞机翼形二维几何模型的线性三角形网格示意图；

图3为图2中线性三角形网格升阶得到的升阶三角形网格示意图；

图4为本发明针对图3中边界线单元的偏移矢量场示意图；

图5为本发明针对图3中升阶三角形网格的偏移矢量场图；

图6(a)为图5中A部分的局部放大图；

图6(b)为图5中B部分的局部放大图；

图7为本发明针对图1中飞机翼形二维几何模型生成的高阶网格示意图；

图8(a)为图7中C部分的局部放大图；

图8(b)为图7中D部分的局部放大图。

具体实施方式

下面结合附图对本发明作进一步说明。

一种基于偏微分方程求解的二维高阶网格生成方法，具体步骤如下：

步骤1、根据飞机翼形二维几何模型生成覆盖飞机翼形二维几何模型的线性网格，建立线性网格边界到飞机翼形二维几何模型边界的映射关系，具体如下：

如图1所示，输入飞机翼形二维几何模型，(采用现有低阶网格生成方法)生成覆盖飞机翼形二维几何模型的较为稀疏的线性三角形网格，如图2所示，记线性三角形网格边界处的线单元为ΔL_j,j∈[1,N],N为线单元数；记飞机翼形二维几何模型边界为将线性三角形网格边界处每个线单元上的网格点投影到最靠近该线单元的飞机翼形二维几何模型边界上，从而将飞机翼形二维几何模型边界划分成N段，记为几何边界段/>j∈[1,N]；建立线单元ΔL_j到几何边界段/>的映射/>

步骤2、对覆盖飞机翼形二维几何模型的线性网格中每个网格单元进行升阶，得到升阶三角形网格，具体如下：

在每个网格单元的线单元上以及内部均添加额外的网格点，并通过插值变换将每个线单元用高阶几何形式表示，得到升阶三角形网格，如图3所示；插值变换如下：设以m次函数N_α(ε)作为插值基函数，且每个线单元均有m+1个网格点，则线单元的高阶几何形式为：

式中，为第α个网格点在i方向上的值；p为线单元上网格点通式表达，i表示二维上的x或y分量，每个线单元有m+1个网格点，N_α(ε)表示第α个网格点对应的插值基函数，ε为插值基函数的参数，取值范围是[-1,1]。

步骤3、建立偏移矢量场最优化模型和约束条件，具体如下：

以升阶三角形网格边界的线单元偏移矢量向升阶三角形网格内部推进时满足最小变化为目标，建立偏移矢量场最优化模型：

约束条件为：

其中，Ω为升阶三角形网格的内部区域，u_i(p)为位于升阶三角形网格上的网格点p(p为网格点表达参量)处的偏移矢量在i方向的分量，为u_i(p)梯度的2范数，和/>分别为由升阶三角形网格边界得到的狄利克雷边界和诺依曼边界，n为边界上网格点p处的外法向量，u_i,0(p)为升阶三角形网格边界的网格点p到飞机翼形二维几何模型边界的偏移矢量沿i方向的分量，/>为u_i,0(p)在n方向上的偏导。

偏移矢量场最优化模型的最优解等价为如下拉普拉斯边值问题的解：

式(4)中，Δu_i(p)为u_i(p)的拉普拉斯算子。

通过线单元ΔL_j到几何边界段的映射关系，并采用步骤2中的插值变换将升阶三角形网格边界的每个线单元偏移矢量的约束条件用高阶几何形式表示(如图4所示)，得到：

p^α表示第α个网格点；

步骤4、求解得到升阶三角形网格的偏移矢量场，将升阶三角形网格的所有网格点沿各自对应的偏移矢量进行偏移，偏移后的每个线单元用高阶几何形式表示，得到高阶网格，具体过程如下：

根据公式(4)，升阶三角形网格的任意一个内部点p₁处的偏移矢量沿i方向的分量通过如下边界积分公式求解：

其中，ΔL为线性三角形网格边界处线单元的集合，点p₂为位于升阶三角形网格边界线单元上的点的一种参量表达,表示点p₂邻域的微分，G(p₂,p₁)为Δu_i(p)＝0方程的基本解，即/>n2为点p₂处的外法向量。

同时，由公式(4)，当升阶三角形网格边界线单元上每个点的u_i(p₂)和中有一个值已知时(u_i(p₂)可由公式(5)求得)，另一个值通过如下积分公式得到：

其中，点p₃为升阶三角形网格边界线单元上的点的另一种参量表达，假设升阶三角形网格边界的各个线单元之间光滑连接，则常量 表示点p₃邻域的微分，/>n3为点p₃处的外法向量。公式(7)的求解采用高阶边界元数值方法，结合公式(5)，得到如下离散形式：

其中，q_k1为位于线单元ΔL_k1上的点，k1∈[1,N]，q_j为位于ΔL_j上的点，nq_j为q_j处的外法向量，表示点q_j邻域的微分，/>为第j个边界线单元上的第α个网格点。当q_k1和q_j重合或处于同一个线单元时，采用高斯积分求解。公式(8)表示一个由N*m个未知数和N*m个方程组成的方程组。通过求解该方程组可得到位于升阶三角形网格边界线单元上的边值u_i(p)或/>随后，通过公式(6)求解位于升阶三角形网格内部的网格点偏移矢量，从而获得整个升阶三角形网格的偏移矢量场，如图5、6(a)和6(b)所示。

将升阶三角形网格的所有网格点沿各自对应的偏移矢量进行偏移，并采用步骤2中的插值变换将升阶三角形网格的所有网格点偏移后的每个线单元用高阶几何形式表示，得到高阶网格，如图7、8(a)和8(b)所示。

Claims (2)

1.一种基于偏微分方程求解的二维高阶网格生成方法，其特征在于：具体包括以下步骤：

步骤1、根据飞机翼形二维几何模型生成覆盖飞机翼形二维几何模型的线性网格，建立线性网格边界到飞机翼形二维几何模型边界的映射关系，具体如下：

输入飞机翼形二维几何模型，生成覆盖飞机翼形二维几何模型的较为稀疏的线性三角形网格，记线性三角形网格边界处的线单元为ΔL_j，j∈[1，N]，N为线单元数；记飞机翼形二维几何模型边界为将线性三角形网格边界处每个线单元上的网格点投影到最靠近该线单元的飞机翼形二维几何模型边界上，从而将飞机翼形二维几何模型边界划分成N段，记为几何边界段/>j∈[1，N]；建立线单元ΔL_j到几何边界段/>的映射/>

步骤2、对覆盖飞机翼形二维几何模型的线性网格中每个网格单元进行升阶，得到升阶三角形网格，具体如下：

在每个网格单元的线单元上以及内部均添加额外的网格点，并通过插值变换将每个线单元用高阶几何形式表示，得到升阶三角形网格；插值变换如下：设以m次函数N_α(ε)作为插值基函数，且每个线单元均有m+1个网格点，则线单元的高阶几何形式为：

式中，为第α个网格点在i方向上的值；p为线单元上网格点通式表达，i表示二维上的x或y分量，每个线单元有m+1个网格点，N_α(ε)表示第α个网格点对应的插值基函数，ε为插值基函数的参数，取值范围是[-1，1]；

步骤3、建立偏移矢量场最优化模型和约束条件，具体如下：

以升阶三角形网格边界的线单元偏移矢量向升阶三角形网格内部推进时满足最小变化为目标，建立偏移矢量场最优化模型：

约束条件为：

其中，Ω为升阶三角形网格的内部区域，u_i(p)为位于升阶三角形网格上的网格点p处的偏移矢量在i方向的分量，为u_i(p)梯度的2范数，/>和/>分别为由升阶三角形网格边界得到的狄利克雷边界和诺依曼边界，n为边界/>上网格点p处的外法向量，u_i，0(p)为升阶三角形网格边界的网格点p到飞机翼形二维几何模型边界的偏移矢量沿i方向的分量，/>为u_i，0(p)在n方向上的偏导；

偏移矢量场最优化模型的最优解等价为如下拉普拉斯边值问题的解：

式(4)中，Δu_i(p)为u_i(p)的拉普拉斯算子；

通过线单元ΔL_j到几何边界段的映射关系，并采用步骤2中的插值变换将升阶三角形网格边界的每个线单元偏移矢量的约束条件用高阶几何形式表示，得到：

p^α表示第α个网格点；

步骤4、求解得到升阶三角形网格的偏移矢量场，将升阶三角形网格的所有网格点沿各自对应的偏移矢量进行偏移，偏移后的每个线单元用高阶几何形式表示，得到高阶网格。

2.根据权利要求1所述一种基于偏微分方程求解的二维高阶网格生成方法，其特征在于：步骤4具体过程如下：

根据公式(4)，升阶三角形网格的任意一个内部点p₁处的偏移矢量沿i方向的分量通过如下边界积分公式求解：

其中，ΔL为线性三角形网格边界处线单元的集合，点p₂为位于升阶三角形网格边界线单元上的点的一种参量表达，表示点p₂邻域的微分，G(p₂，p₁)为Δu_i(p)＝0方程的基本解，即/>n2为点p₂处的外法向量；

由公式(5)求得升阶三角形网格边界线单元上每个点的u_i(p₂)，同时，由公式(4)，升阶三角形网格边界线单元上每个点的通过如下积分公式得到：

其中，点p₃为升阶三角形网格边界线单元上的点的另一种参量表达，假设升阶三角形网格边界的各个线单元之间光滑连接，则常量 表示点p₃邻域的微分，n3为点p₃处的外法向量；公式(7)的求解采用高阶边界元数值方法，结合公式(5)，得到如下离散形式：

其中，q_k1为位于线单元ΔL_k1上的点，k1∈[1，N]，q_j为位于ΔL_j上的点，nq_j为q_j处的外法向量，表示点q_j邻域的微分，/>为第j个边界线单元上的第α个网格点；当q_k1和q_j重合或处于同一个线单元时，采用高斯积分求解；通过公式(8)求解得到位于升阶三角形网格边界线单元上的边值u_i(p)或/>随后，通过公式(6)求解位于升阶三角形网格内部的网格点偏移矢量，从而获得整个升阶三角形网格的偏移矢量场；

将升阶三角形网格的所有网格点沿各自对应的偏移矢量进行偏移，并采用步骤2中的插值变换将升阶三角形网格的所有网格点偏移后的每个线单元用高阶几何形式表示，得到高阶网格。