CN111079353A

CN111079353A - 应用于复杂流体分析的快速均匀网格划分的方法及装置

Info

Publication number: CN111079353A
Application number: CN201911304744.1A
Authority: CN
Inventors: 甘华权; 李海艳; 骆继明; 黄运保
Original assignee: Guangdong University of Technology
Current assignee: Guangdong University of Technology
Priority date: 2019-12-17
Filing date: 2019-12-17
Publication date: 2020-04-28
Anticipated expiration: 2039-12-17
Also published as: CN111079353B

Abstract

应用于复杂流体分析的快速均匀网格划分的方法，包括如下步骤：步骤一：定义区域轮廓，实现初始多区域边界的建模；步骤二：输入所述定义区域轮廓的轮廓信息，判断其封闭轮廓环是否由两条或以上的边构成，若是，则进入步骤三；若否，则先将所述封闭轮廓环分割成四段，再进入步骤三；步骤三：确定所述封闭轮廓环的四条边界曲线，再实现U向边界和V向边界的定义。应用于复杂流体分析的快速均匀网格划分的装置，包括第一定义模块、输入判断模块、分割模块、第二定义模块、划分模块、拟合模块，生成模块。本发明能够对复杂轮廓区域进行有效网格划分的方法，以保证划分网格的质量进而保证计算流体力学中的数值计算的精度要求。

Description

应用于复杂流体分析的快速均匀网格划分的方法及装置

技术领域

本发明涉及网格划分技术领域，尤其涉及应用于复杂流体分析的快速均匀网格划分的方法及装置。

背景技术

计算流体力学是通过计算机数值计算和图像显示，对包含有流体流动和热传导等相关物理现象的系统所做的分析。其基本思想是把时间域及空间域上连续的物理量的场，如速度场和压力场，用一系列有限个离散点上的变量值的集合来替代，通过一定的原则和方式建立起关于这些离散点上场变量之间关系的代数方程组，然后求解方程得到场变量的近似值。

计算流体力学的任务是流体力学的数值模拟，建立数学模型后寻求高效率，高准确度的计算方法。计算方法不仅包括数学方程的离散化及求解方法，还包括计算网格的建立，边界条件的处理，其中网格的划分属于计算流体力学计算中的预处理部分。网格的生成是将计算域划分为较小的，不重叠的子域或单元(网格)，单元的数目决定了计算流体力学的计算精度，而网格的细密程度决定精度等，并且网格划分是计算流体力学中的关键技术之一，网格划分实际就是将结构进行离散化，从而得到有限单元及节点，进而通过求解用近似解代替实际解。

网格剖分主要分为映射网格划分和自由网格划分两类。其中自由网格划分又分为节点连接法、拓扑分解法、几何分解法、栅格法等。

映射网格划分是对规整模型的一种规整网格划分方法，映射网格要求面或体的形状是规则的，也就是说它们必须遵循一定的规则。其基本思路为在简单区域内采用某种映射函数构造简单区域的边界点和内点，并按某种规则连接结点构成网格单元。也就是根据形体边界的参数方程，利用映射函数，把参数空间内单元正方形或单元三角形(对于三维问题是单元立方体或单元四面体)的网格映射到欧氏空间，从而生成实际的网格。

节点连接法根据给定的精度对曲面的边界及参数域边界进行离散,生成一组边界节点,再根据给定初始密度对参数域内部布置一定数量节点。然后根据波前法或Delaunay三角剖分得到网格，将参数域所得的三角形网格一一映射到空间曲面,形成参数域与空间域一一对应的拓扑连接关系.对空间曲面网格的每个三角形进行校验,如果该三角形对曲面的逼近误差超过所需精度,则对该三角形进一步网格加密。

拓扑分解法假设最后的网格顶点全部由目标边界顶点组成，那么可以用一种三角化算法将目标用尽量少德尔三角形完全分割覆盖。这些三角形主要由目标拓扑结构决定，这样目标的复杂拓扑结构被分解成简单的三角形拓扑结构。

几何分解法凡在产生结点的同时也确定结点间连接关系的方法均称为几何分解法，其节点与单元同步生成，与拓扑分解法完全不同的是在实体分解过程中，考虑了所生成的单元形状及大小，确保生成的单元质量尽量的好。

栅格法又分为栅格叠合法和四(八)叉树法，栅格叠合法是将栅格置放于物体上，既可在栅格的规则点处布置节点，也可在栅格单元中随机布置结点，容易得到均匀网格，然后对栅格与物体求交，将在物体外的栅格去掉，物体内的栅格保留作为网格。在边界处的栅格需要变形，也就是要调整节点的位置以满足边界条件的要求，最终产生的网格内部单元完全一样，只是边界单元不同。栅格越密网格质量将越好。

叉树法网格生成的过程是完全面向几何特征的，它通过一系列几何操作，同时从几何模型数据库中取得分析物体的几何信息，以及网格控制参数信息等，来完成网格划分。

然而，映射法它不是完全面向几何特征的，很难完成自动化，尤其对3D区域。它是通过低维点生成高维单元，使得对单元的定位是十分困难的，使得对局部的控制能力下降。且各映射块之间的网格密度相互影响程度很大。

节点连接法主要依赖于Delaunay三角剖分法，其容易出现所谓的退化现象或产生所谓的Sliver单元，对三维的情形，基于最大二维最小角判断的对角线交换规则不在成立，而基于外接圆判据的Delaunay三角化一般也不再能保证生成的网格质量。且Delaunay三角化网格自动生成方法的事前(a-priori)网格控制能力不强。

拓扑分解法生成的网格一般相当粗糙，必须与其他方法相结合，通过网格加密等过程，才能生成合适的网格，且单一的拓扑分解法因只依赖于几何体的拓扑结构使网格剖分不理想，有时甚至很差。

几何分解法其自动化程度低，不利于复杂网格的生成，而栅格法所生成网格与所选择的初始栅格及其取向有关,网格边界单元质量差。

传统的网格划分方法，对形状是规则的轮廓区域具有很好的适用性，对复杂的轮廓区域划分的网格很差且还需要用户进行繁琐的调整，而对于计算流体力学中由多种物质组成的复杂流体的区域来说这些方法均不能保证有效的划分网格，也无法保证网格的质量。

发明内容

本发明的目的在于提出应用于复杂流体分析的快速均匀网格划分的方法及装置，能够对复杂轮廓区域进行有效网格划分的方法，以保证划分网格的质量进而保证计算流体力学中的数值计算的精度要求。

为达此目的，本发明采用以下技术方案：

应用于复杂流体分析的快速均匀网格划分的方法，包括如下步骤：

步骤一：定义区域轮廓，实现初始多区域边界的建模；

步骤二：输入所述定义区域轮廓的轮廓信息，判断其封闭轮廓环是否由两条或以上的边构成，若是，则进入步骤三；若否，则先将所述封闭轮廓环分割成四段，再进入步骤三；

步骤三：确定所述封闭轮廓环的四条边界曲线，再实现U向边界和V向边界的定义；

步骤四：对所述U向边界和V向边界进行划分,并确定划分点的数目及划分比例；

步骤五：用B样条曲线对各轮廓边界曲线进行拟合，得到轮廓曲线，再对拟合后的轮廓曲线生成B样条曲面；

步骤六：通过所述U向边界和V向边界划分的划分点的参数，利用等参线自动生成比例均匀网格。

进一步，还包括：步骤七：调整所述网格，确保边的端点是网格顶点。

进一步，步骤五具体的过程为：用B样条曲线对各轮廓边界曲线进行拟合，得到轮廓曲线，用拟合后的各条轮廓曲线的U向节点矢量和V向节点矢量作为B样条曲面的节点矢量，通过插值于轮廓曲线和跨界导矢得到相对应的B样条曲面，再通过逼近内部数据点得到对应的准均匀B样条曲面。

进一步，在步骤五中，采用能量优化法对所述数据点进行B条样曲线拟合。

进一步，在步骤五中，采用带光顺项的最小二乘拟合方法对所述数据点进行B条样曲线拟合。

应用于复杂流体分析的快速均匀网格划分的装置，包括：

第一定义模块，用于定义区域轮廓，实现初始多区域边界的建模；

输入判断模块，用于输入所述定义区域轮廓的轮廓信息，判断其封闭轮廓环是否由两条或以上的边构成；

分割模块，用于将所述封闭轮廓环分割成四段；

第二定义模块，用于确定所述封闭轮廓环的四条边界曲线，再实现U向边界和V向边界的定义；

划分模块，用于对所述U向边界和V向边界进行划分,并确定划分点的数目及划分比例；

拟合模块，用于用B样条曲线对各轮廓边界曲线进行拟合，得到轮廓曲线，再对拟合后的轮廓曲线生成B样条曲面；

生成模块，用于通过所述U向边界和V向边界划分的划分点的参数，利用等参线自动生成比例均匀网格。

进一步，还包括：

调整模块，用于调整所述网格，确保边的端点是网格顶点。

进一步，所述拟合模块用于用B样条曲线对各轮廓边界曲线进行拟合，得到轮廓曲线，再对拟合后的轮廓曲线生成B样条曲面的具体的过程为：用B样条曲线对各轮廓边界曲线进行拟合，得到轮廓曲线，用拟合后的各条轮廓曲线的U向节点矢量和V向节点矢量作为B样条曲面的节点矢量，通过插值于轮廓曲线和跨界导矢得到相对应的B样条曲面，再通过逼近内部数据点得到对应的准均匀B样条曲面。

进一步，所述拟合模块采用能量优化法对所述数据点进行B条样曲线拟合。

进一步，所述拟合模块采用带光顺项的最小二乘拟合方法对所述数据点进行B条样曲线拟合。

本发明根据上述内容提出应用于复杂流体分析的快速均匀网格划分的方法及装置。

本发明采用一种基于B样条曲线/曲面的网格划分的方法，从而使得对任意复杂曲面的进行网格划分，都可以得到质量很好的网格。

本发明通过保证组成边的端点必须是网格顶点，从而使得边界网格的质量得到了很好的保证。

本发明通过拟合生成均匀B样条曲面后再采用等参线从而自动生成比例均匀的网格，不管多复杂封闭轮廓都可以快速生成比例均匀的网格，且时间与复杂程度是线性正相关的关系，并不会因为图形复杂而导致运算时间指数级增加。

附图说明

图1是本发明其中一个实施例的方法流程图；

图2是本发明其中一个实施例的网格划分过程图；

图3是本发明其中一个实施例的网格划分过程图；

图4是本发明其中一个实施例的网格划分过程图。

具体实施方式

下面结合附图并通过具体实施方式来进一步说明本发明的技术方案。

如图1-4所示，应用于复杂流体分析的快速均匀网格划分的方法，包括如下步骤：

步骤一：定义区域轮廓，实现初始多区域边界的建模；本实施例是利用自主的三维CAD系统InteSolid中的二维草图定义模块，快速定义区域轮廓，实现初始多区域边界的建模。

步骤二：输入在步骤一定义的所述区域轮廓的轮廓信息，其边集会自动按照顺或逆时针首尾相接，构成封闭轮廓环，判断所述封闭轮廓环是否由两条或以上的边构成，若是，则进入步骤三；若为一条边构成，如圆，则先将所述封闭轮廓环分割成四段，再进入步骤三进行操作。

步骤三：确定所述封闭轮廓环的四条边界曲线，再实现U向边界和V向边界的定义(可退为空，即为某一边由边退化为点，例如一个四边形的某一边退化为一个点了，此时四边形就变成三角形了)；如图2所示,选取左边的边零0和边一1为第一边界,上边的边二2和边三3为第二边界,右边的边四4为第三边界,下边的边五5为第四边界。然后定义第一边界的起点为U方向的起始点，第二边界的起点为V方向的起始点，从而实现U向边界及V向边界的定义。

步骤四：对所述U向边界和V向边界进行划分,可按均匀划分或一定比例划分，开始点位于环向的起点，并可定义区域内材料的密度，材料ID等信息。即首先定义划分比例，然后再定义划分点的个数，最后自由定义区域内材料的密度，材料ID等信息。如图1所示,第一边界与第三边界逻辑上互为对边，第二边界与第四边界逻辑上互为对边，故假设设定第一边界的划分点的个数为20个,划分比例为1,即U向边界划分点的个数为20个且U向边界起始点是划分起点，第二边界的划分点的个数为10,划分比例也为1,即V向边界划分点的个数为10个且V向起始点是划分起点，则可由求和公式

(其中:q为比例,l为其对应边长,a0为第一段的长度,n为点分割后曲线的段数)确定点划分后每段曲线的长度或每个点在对应边上的位置。

步骤五：用B样条曲线对各轮廓边界曲线进行拟合，得到轮廓曲线，再对拟合后的轮廓曲线生成B样条曲面。

采用B样条曲线对各条轮廓边界曲线进行拟合,得到轮廓曲线，用拟合后的各条轮廓曲线的U向节点矢量和V向节点矢量作为B样条曲面的节点矢量，通过插值于其曲线和跨界导矢得到相对应的B样条曲面，再通过逼近内部数据点得到对应的准均匀B样条曲面。如图2，以第一边界为例，以U为节点矢量，采用能量优化法(或带光顺项的最小二乘拟合方法)对数据点进行B样条曲线拟合，拟合时要求插值于两端点，并逼近内部数据点，由此可得到以U为节点矢量的边界轮廓的拟合3次准均匀B样条曲线，然后检查拟合误差，若存在不满足精度要求的点，则以U为节点矢量，重新对数据点进行拟合，直到满足精度要求为止。其余边界也类似的得到拟合的3次准均匀B样条曲线。最后用拟合的3次均匀B样条曲线的U向节点矢量和V向节点矢量作为B样条曲面的节点矢量，通过插值于其曲线和跨界导矢得到相对应的B样条曲面，再通过逼近内部数据点得到更好准均匀B样条曲面。

定义了U向边界和V向边界后，而U向节点矢量＝{u₀，u₁，u₂，u₃，...，u_n+p+1}，相当于在定义了U向边界后，U向节点矢量就是在U边界内取值，即点集，可以理解为一条直线参数化后，在其内取n+p+1个点组成U向节点矢量，V向节点矢量也同理。

具体是，由步骤四得到U向边界及V向边界划分的划分点，然后根据划分点的U、V的值(通过步骤三定义了U向边界和V向边界，把曲线参数化，步骤四的U向边界得分划分点在U上有相对应的U、V值，同理V边界划分点也如此，其中u∈[0,1]，v∈[0,1])，利用四边域B样条曲面插值方法，使样条曲面上等参数线与生成U向网格或V向网格线近似重合，进而可利用等参线自动生成比例均匀网格，以满足二维辐射流体力学模拟需要。如图2，U方向划分点为20个，V方向划分点为10个，则等参线自动生成比例均匀网格数目为20*10＝200个。由于划分点是通过U向边界及V向边界均匀划分或一定比例划分得到的并且采用的是均匀B样条曲面，从而保证了利用等参线自动生成的网格成比例均匀。

步骤七：调整所述网格，对于边界处确保边的端点必须是网格顶点，先按以上步骤一至步骤六生成网格，同时考虑到区域边界上一阶不连续点的约束，可自动将临近网格均匀微调，使其通过这些不连续点。即先以组成的边计算每条边的离散网格数，再以各组成边分别进行比例划分，确保组成边端点必须是网格顶点。对比如图3和图4，可知原来在Ⅰ、Ⅱ处由步骤四中得到的划分点并没有与边的端点重合，使得其边零0、边一1的交点处的网格划分质量不够好，故通过自动调整划分点使得分别与网格顶点跟边零0、边一1处的交点重合、边二2、边三3处的交点重合，从而保证了生成网格边界单元的质量。

本方法通过边集自动成环及对边界轮廓曲线的边界选定，从而实现U向边界及V向边界的定义。

本方法通过用自定义划分点的个数和划分比例来对U向边界和V向边界进行均匀划分或一定比例划分来实现对网格数量的控制。

本方法通过划分点的参数利用插值方法，使样条曲面上等参数线与生成网格的U向网格线或V向网格线近似重合的方法，进而做到了利用等参线自动生成比例均匀的网格，从而实现了高质量的网格划分。

本方法可应用于任意复杂的区域轮廓，并且可以很好的保证网格质量，并通过略微调整即可保证网格边界单元的质量，很好的保证了计算流体力学中的数值计算的精度要求。并且本方法拟合生成均匀B样条曲面后，采用等参线一次自动生成的网格成比例均匀，时间效率高。即使对于很复杂的区域轮廓来说，仍时间效率高。

在步骤五中，采用能量优化法或带光顺项的最小二乘拟合方法对所述数据点进行B条样曲线拟合。其中能量法属于整体优化方法，即在给定的容差条件下，使曲线的应变能函数E达到最小，通过优化问题进行求解从而实现对曲线的光顺，这两种方法生成的曲线具有较好的数值稳定性、具有较好的光顺性，

应用于复杂流体分析的快速均匀网格划分的装置，包括：

分割模块，用于将所述封闭轮廓环分割成四段；

进一步，还包括：调整模块，用于调整所述网格，确保边的端点是网格顶点。

进一步，所述拟合模块采用能量优化法或对带光顺项的最小二乘拟合方法所述数据点进行B条样曲线拟合。

本装置通过边集自动成环及对边界轮廓曲线的边界选定，从而实现U向边界及V向边界的定义。

本装置通过用自定义划分点的个数和划分比例来对U向边界和V向边界进行均匀划分或一定比例划分来实现对网格数量的控制。

本装置通过划分点的参数利用插值装置，使样条曲面上等参数线与生成网格的U向网格线或V向网格线近似重合的装置，进而做到了利用等参线自动生成比例均匀的网格，从而实现了高质量的网格划分。

本装置可应用于任意复杂的区域轮廓，并且可以很好的保证网格质量，并通过略微调整即可保证网格边界单元的质量，很好的保证了计算流体力学中的数值计算的精度要求。并且本装置拟合生成均匀B样条曲面后，采用等参线一次自动生成的网格成比例均匀，时间效率高。即使对于很复杂的区域轮廓来说，仍时间效率高。

本发明的所述能量优化法，即要求在满足插值条件的前提下，使得曲线或曲面的能量达到最小，而使生成的曲线具有较好的光顺性。

在U向节点矢量＝{u₀，u₁，u₂，u₃，...，u_n+p+1}上的p次B样条曲线为

其中，V_j(j＝0,1,...,n)是控制顶点，N_j,p(t)(j＝0,1,...,n)是定义在U上的p次B样条基函数，由下式确定：

规定0/0＝0

由于实际应用中使用最多的是三次B样条曲线，故p取3，此外，在U向节点矢量中，通常取u₀＝u₁＝....＝u_p＝0，u_n+1＝u_n+2＝......＝u_n+p+1＝1,则在B样条曲线的断点处有：P(0)＝V₀,P(1)＝V₁,根据各条轮廓边界曲线得到m+1个型值点P_j(j＝0,1,...,n)，并设P_j对应的参数是

且

则插值于P_j(j＝0,1,...,n)且以U为节点矢量的3次B样条曲线P(t)为：

根据插值条件，可得m+1个方程：

令

V＝[V₀,V₁,V₂,...,V_n]^T

P＝[P₀,P₁,P₂,...,P_n]^T

则

可以写成矩阵的形式：

AV＝P

其中

为(m+1)行(m+1)列的矩阵。

一般将B样条曲线P(t)看成一弹性样条，其能量函数为：E＝∫k²ds由于其涉及到曲率和弧长，一般计算中，通常简化为

将

代入其中可得：

上式是V＝[V₀,V₁,V₂,...,V_n]^T的二次式，可用矩阵表示为：

E＝V^TMV

式中

是一个(n+1)阶的方阵。

采用能量法插值就是求曲线P(t)插值于P_j(j＝0,1,...,n)并且使得E＝V^TMV最小，即在线性约束条件AV＝P下，使E＝V^TMV达到最小。故引入C＝(c₀,c₁,...,c_m)^T,将约束问题转化无约束问题：

分别对V和C求导并令其等于0，可得

MV+A^TC＝0

AV-P＝0

写成分块矩阵的形式

求解线性方程组即可得到控制顶点向量V＝[V₀,V₁,V₂,...,V_n]^T。当M和AM^-1A^T均可逆，方程组有唯一解。

且

代入AV-P＝0中得

将C＝-(AM^-1A^T)^-1P代入V＝-M^-1A^TC,得

V＝M^-1A^T(AM^-1A^T)^-1P

根据上式，即可得到插值与型值点P_j(j＝0,1,...,n)并且使能量函数最小的B样条曲线P(t)的控制顶点V＝[V₀,V₁,V₂,...,V_n]^T。

本发明的所述带光顺项的最小二乘拟合

给定U向节点矢量＝{u₀＝u₁＝u₂＝u_3＝0，...,u_i,...，u_n+1＝u_n+2＝u_n+3＝u_n+4＝1}和m+1个型值点P_j(j＝0,1,...,n)，并设P_j对应的参数是

且

带光顺项的最小二乘拟合方法的原理是：除了要求曲线逼近给定的型值点外，还要求曲线的能量函数(可以作为曲线光顺性的度量)尽可能小，即

达到最小，式中，第一项称为误差项；第二项称为光顺项，是曲线光顺性的度量；μ_k(k＝0,1,...,m)和λ均为大于等于0的常数。曲线的能量函数为

且

则

的矩阵形式为：

J＝D+λE＝(AV-P)^TD(AV-P)+λ(V^TMV)

式中，V＝[V₀,V₁,V₂,...,V_n]^T为控制顶点向量；P＝[P₀,P₁,P₂,...,P_n]^T为型值点向量；A为

(m+1)行(m+1)列的矩阵；

M为

(n+1)行(n+1)列的方阵；D＝diag[μ₀,μ₁,...,μ_n]为对角矩阵。

为了使J达到最小，令

得

(A^TDA+λM)V＝A^TDP

可以证明，只要λ＞0，则方程组中的系数矩阵必定为对称正定矩阵，故存在唯一解，且数值稳定性很好。求解方程组就可以得到逼近P_j(j＝0,1,...,m)的B样条曲线P(t)的控制顶点。

以上结合具体实施例描述了本发明的技术原理。这些描述只是为了解释本发明的原理，而不能以任何方式解释为对本发明保护范围的限制。基于此处的解释，本领域的技术人员不需要付出创造性的劳动即可联想到本发明的其它具体实施方式，这些方式都将落入本发明的保护范围之内。

Claims

1.应用于复杂流体分析的快速均匀网格划分的方法，其特征在于：包括如下步骤：

步骤一：定义区域轮廓，实现初始多区域边界的建模；

2.根据权利要求1所述的应用于复杂流体分析的快速均匀网格划分的方法，其特征在于：还包括：

步骤七：调整所述网格，确保边的端点是网格顶点。

3.根据权利要求1所述的应用于复杂流体分析的快速均匀网格划分的方法，其特征在于：

步骤五具体的过程为：用B样条曲线对各轮廓边界曲线进行拟合，得到轮廓曲线，用拟合后的各条轮廓曲线的U向节点矢量和V向节点矢量作为B样条曲面的节点矢量，通过插值于轮廓曲线和跨界导矢得到相对应的B样条曲面，再通过逼近内部数据点得到对应的准均匀B样条曲面。

4.根据权利要求3所述的应用于复杂流体分析的快速均匀网格划分的方法，其特征在于：

在步骤五中，采用能量优化法对所述数据点进行B条样曲线拟合。

5.根据权利要求3所述的应用于复杂流体分析的快速均匀网格划分的方法，其特征在于：

在步骤五中，采用带光顺项的最小二乘拟合方法对所述数据点进行B条样曲线拟合。

6.应用于复杂流体分析的快速均匀网格划分的装置，其特征在于：包括：

分割模块，用于将所述封闭轮廓环分割成四段；

7.根据权利要求6所述的应用于复杂流体分析的快速均匀网格划分的装置，其特征在于：还包括：

调整模块，用于调整所述网格，确保边的端点是网格顶点。

8.根据权利要求6所述的应用于复杂流体分析的快速均匀网格划分的装置，其特征在于：

所述拟合模块用于用B样条曲线对各轮廓边界曲线进行拟合，得到轮廓曲线，再对拟合后的轮廓曲线生成B样条曲面的具体的过程为：用B样条曲线对各轮廓边界曲线进行拟合，得到轮廓曲线，用拟合后的各条轮廓曲线的U向节点矢量和V向节点矢量作为B样条曲面的节点矢量，通过插值于轮廓曲线和跨界导矢得到相对应的B样条曲面，再通过逼近内部数据点得到对应的准均匀B样条曲面。

9.根据权利要求6所述的应用于复杂流体分析的快速均匀网格划分的装置，其特征在于：

所述拟合模块采用能量优化法对所述数据点进行B条样曲线拟合。

10.根据权利要求6所述的应用于复杂流体分析的快速均匀网格划分的装置，其特征在于：

所述拟合模块采用带光顺项的最小二乘拟合方法对所述数据点进行B条样曲线拟合。