CN110610042A - 基于阵面推进的非结构网格结构化处理方法 - Google Patents
基于阵面推进的非结构网格结构化处理方法 Download PDFInfo
- Publication number
- CN110610042A CN110610042A CN201910853374.0A CN201910853374A CN110610042A CN 110610042 A CN110610042 A CN 110610042A CN 201910853374 A CN201910853374 A CN 201910853374A CN 110610042 A CN110610042 A CN 110610042A
- Authority
- CN
- China
- Prior art keywords
- grid
- array
- nextfronts
- edge
- unstructured
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
Landscapes
- Image Generation (AREA)
Abstract
本发明公开了一种基于阵面推进的非结构网格结构化处理方法，目的是解决非结构网格结构化处理中存在的局部方向与应用问题无关、三角形非结构网格正交性差、计算量较大的问题。技术方案是生成二维非结构网格；根据非结构网格相关数据结构定义并初始化非结构网格单元的关键元素数组和局部方向数组；采用阵面推进确定非结构网格的关键元素和局部方向，输出关键元素数组和局部方向数组，实现非结构网格结构化处理。采用本发明可高效生成二维三角形、四边形或三角形‑四边形混合非结构网格的局部方向，提高局部方向正交性能，且其中一个局部方向与壁面法向方向接近，有效提高火箭壁面温度梯度的精度，且避免了穷举法的冗余计算。
Description
技术领域
本发明涉及非结构网格结构化处理方法，尤指采用阵面推进进行非结构网格的局部方向生成，使非结构网格具备规则性和方向性的非结构网格结构化处理方法。
背景技术
计算机分析和处理工程应用的过程中常需使用网格对连续物理空间进行离散。随着工程应用研究问题的复杂度增加，空间离散所需的网格量剧增，提高工程应用的计算效率和计算精度成为国内外研究的热点问题。网格处理是解决这些问题的重要方法之一，针对不同网格类型有不同的网格处理方法。
目前网格主要分为结构网格和非结构网格两大类。结构网格形状规则，网格单元沿着笛卡尔坐标系的坐标轴方向规则排列。采用结构网格容易获得较高的计算精度和计算效率，但结构网格的网格生成、网格变形、网格加密等都非常困难且耗时，需要大量人工干预，大大增加了工程应用的时间和人力开销。非结构网格没有排列和形状的限制，采用非结构网格进行网格生成、网格变形、网格加密等基本不需要人工参与，可直接使用商业软件自动实现，更适用于复杂的工程应用。但非结构网格的计算精度和计算效率较低，一方面，非结构网格依靠网格生成软件自动生成，网格质量得不到保证，可能影响计算精度；另一方面，非结构网格不具备规则性和方向性，无法选择对称的模板单元进行截断误差抵消，影响了计算精度和计算效率。
非结构网格结构化处理是一种有效提高非结构网格计算精度和计算效率的方法，其主要内容是通过为非结构网格生成局部方向，使其具有如结构网格一样的规则性和方向性。结构化处理后的非结构网格可提高工程应用的性能，例如将结构化处理后的非结构网格用于火箭壁面温度的梯度重构，梯度重构过程可沿局部方向进行，可减少温度梯度求解的计算量并提高火箭温度梯度的精度，使得火箭温度测点布置更加合理。
现有研究中关于非结构网格结构化处理的研究较少。专利号为201510577643.7的中国专利“结构网格动网格技术用于非结构网格流场求解器的方法”，涉及到了通过生成初始外形结构网格，并建立非结构网格与结构网格映射实现非结构网格结构化处理。但该方法仍需要生成结构网格，适用范围较窄，无法针对已生成好的非结构网格进行结构化处理，无法避免结构网格生成困难的问题。Sozer在文献[Emre S等，Gradient calculationmethods on arbitrary polyhedral unstructured meshes for Cell-Centered CFDsolvers，52th AIAA Aerospace Sciences Meeting，2014。Emre S等，单元中心型CFD求解器中任意多边形非结构网格的梯度计算方法，52届美国航空航天协会会议，2014]中提出了曲线梯度重构法，实现了基于几何外形的二维非结构网格结构化处理。该方法选取雅克比矩阵行列式值最大的两组邻接网格单元构成非结构网格单元的局部方向，实现了非结构网格的结构化。曲线梯度重构法中的非结构网格结构化处理方法(简称背景技术1)总体流程如图1所示，主要包括四个步骤，其中第三步为采用最大化雅克比矩阵行列式值确定非结构网格的关键元素和局部方向，即进行非结构网格结构化处理，第三步的流程图如图2所示。背景技术1的具体步骤如下：
第一步：生成二维非结构网格。
1.1使用专业软件将工程应用需要计算的二维空间进行网格化处理，输出二维非结构网格相关文件。常用的专业软件包括ANSYS公司的Fluent软件[http://www.ansys.com，版本18.0以上]等商业软件、美国国家航空航天局的FUN3D软件[Biedron R等，FUN3D Manual:13.2，2017，版本13.0以上]等In-house软件、以及OpenFOAM[https://www.openfoam.com.版本3.0以上]、SU2[https://su2code.github.io，版本5.0以上]等开源软件。工程应用需要计算的空间在这里指二维计算空间，如火箭截面轮廓的外部空间。二维非结构网格相关文件包括网格单元文件、网格点文件、网格单元的邻接单元文件、网格点的邻接单元文件等。
1.2根据非结构网格相关文件初始化相关数据结构。
1.2.1令网格单元文件为网格单元数组cs，其中Nc为网格单元总数，Ckc为第kc个网格单元。网格单元采用网格点编号数组表示，Ckc＝{kp1,...,kpw,...,kpW}表示第kc个网格单元由第kp1,...,kpw,...,kpW个网格点构成，是一个W边形网格单元，W为正整数。如C1＝{1,4,6}表示第1个网格单元由第1、4、6个网格点构成，是一个三角形网格单元，其它任意形状的二维非结构网格也可类似表示；
1.2.2令网格点文件为网格点数组pts，其中Npt为网格点总数，Ptkp为第kp个网格点。网格点采用笛卡尔坐标数组表示，Ptkp＝(x,y)表示第kp个网格点的坐标为(x,y)，x,y为实数。如Pt1＝{1.0,2.0}表示第1个网格点的坐标为(1.0,2.0)；
1.2.3令网格单元的邻接单元文件为数组cCells，其中CCkc表示网格单元Ckc的邻接网格单元编号组成的数组。如CC1＝{2,3,5,7,11}表示与网格单元C1相邻接的单元为C2,C3,C5,C7,C11；
1.2.4令网格点的邻接单元文件为数组pCells，其中PCkp表示网格点Ptkp的邻接网格单元编号组成的数组。如PC1＝{1,2,4,5,7,8}表示与网格点Pt1相邻接的网格单元为C1,C2,C4,C5,C7,C8；
第二步：根据非结构网格相关数据结构，定义并初始化非结构网格单元的关键元素数组和局部方向数组。
2.1定义关键元素数组其中Elkc为网格单元Ckc的关键元素编号组成的数组。关键元素定义为确定局部方向的网格单元，每个二维非结构网格单元都有4个关键元素，如El1＝{2,4,7,9}表示C1的4个关键元素分别为C2、C4、C7、C9这四个网格单元。初始化均为{0,0,0,0}，从而与有意义的单元编号区别以避免歧义；
2.2根据keyEl定义局部方向数组locDir，其中LDkc表示网格单元Ckc的局部方向向量组成的数组。LDkc包含4个数据项，分别表示Ckc的4个局部方向向量，由Elkc确定，数组A的第ii个数据项用A[ii]表示，ii为正整数：
·LDkc[1]为网格单元Elkc[3]中心指向网格单元Elkc[1]中心构成的向量；
·LDkc[2]为网格单元Elkc[4]中心指向网格单元Elkc[2]中心构成的向量；
·LDkc[3]为网格单元Elkc[1]中心指向网格单元Elkc[3]中心构成的向量；
·LDkc[4]为网格单元Elkc[2]中心指向网格单元Elkc[4]中心构成的向量；
由上述定义易知LDkc[1]＝-LDkc[3],LDkc[2]＝-LDkc[4],1≤kc≤Nc。由于局部方向需要在关键元素确定后才能确定，故初始化时的每一项都初始化为零向量；
第三步：采用最大化雅克比矩阵行列式值确定非结构网格的关键元素和局部方向，即进行非结构网格结构化处理，如图2所示，方法是：
3.1根据网格单元数组cs和网格点数组pts计算cs中网格单元中心坐标，得到网格单元中心坐标数组center，其中Ctkc＝{xkc,ykc}表示网格单元Ckc中心的横坐标xkc和纵坐标ykc，中心坐标根据cs和pts和获取网格单元的顶点坐标，并通过多边形中心坐标公式直接计算得到。以三角形网格单元为例，假设Ckc(1≤kc≤Nc)为三角形网格单元，Ckc的中心坐标Ctkc＝{xkc,ykc}计算公式为：
其中Ckc[1]、Ckc[2]、Ckc[3]为Ckc的三个顶点的编号， 为Ckc[1]、Ckc[2]、Ckc[3]这三个顶点横坐标， 为Ckc[1]、Ckc[2]、Ckc[3]这三个顶点纵坐标；
3.2初始化网格单元索引i＝1；
3.3初始化最大雅克比的值MaxJacobi＝0，初始化邻接单元个数CCnum＝0；
3.4以i为网格单元的邻接单元文件cCells的索引，获取网格单元Ci的邻接单元编号，构成Ci的邻接网格单元编号数组CCi，将CCnum更新为CCi的数据项个数；
3.5初始化索引a＝1；
3.6初始化索引b＝a+1；
3.7初始化索引c＝a；
3.8初始化索引d＝c+1；
3.9以a,b,c,d为CCi索引获取Ci的四个邻接单元的编号，分别为CCi[a]、CCi[b]、CCi[c]、CCi[d]；
3.10以CCi[a]、CCi[b]、CCi[c]、CCi[d]为网格单元中心坐标数组center的索引获取这四个邻接单元的单元中心坐标，为方便表示，将这四个邻接单元的单元中心坐标分别记为(ax,ay)、(bx,by)、(cx,cy)、(dx,dy)；
3.11根据雅克比矩阵行列式计算公式，计算这四个邻接单元的雅克比矩阵行列式的值Jacobi：
Jacobi＝(ax-bx)*(cy-dy)-(ay-by)*(cx-dx)；
3.12若Jacobi>MaxJacobi，转3.13；否则转3.14；
3.13更新MaxJacobi＝Jacobi，更新Eli＝{CCi[a],CCi[c],CCi[b],CCi[d]}；
3.14更新d＝d+1，若d<CCnum+1，转3.9；否则转3.15；
3.15更新c＝c+1，若c<CCnum+1，转3.8；否则转3.16；
3.16更新b＝b+1，若b<CCnum+1，转3.7；否则转3.17；
3.17更新a＝a+1，若a<CCnum+1，转3.6；否则转3.18；
3.18根据Eli，按照第2.2步局部方向的定义，计算LDi；
3.19更新i＝i+1，若i<Nc+1，转3.3；否则说明得到了局部方向数组locDir，转3.20；
3.20将关键元素数组keyEl作为关键元素文件输出，将局部方向数组locDir作为局部方向文件输出，实现非结构网格结构化处理。
关键元素文件及局部方向文件即为结构化的非结构网格，将其用于火箭表面温度梯度计算等工程应用中时，可有效提高工程应用的计算效率等方面性能。
背景技术1可有效减少计算量，提高应用的计算效率，但其存在三个方面问题：
1.仅依靠最大化雅克比选择局部方向，局部方向只与网格单元的几何外形有关，而与应用问题无关，对提高火箭表面温度梯度精度的效果有限。
2.对于常用二维三角形非结构网格，局部方向正交性能较差。如对于网格质量较好的正三角形网格，背景技术1生成的局部方向夹角为60°，正交性能差。
3.局部方向采用穷举法获得，计算量大，效率不高。
因此，常用的非结构网格结构化处理方法(背景技术1)仍存在若干问题，目前还没有较好的解决这些问题的方法。
发明内容
本发明要解决的技术问题是：针对非结构网格结构化处理中存在的局部方向与应用问题无关、三角形非结构网格正交性差、计算量较大的问题，提出一种采用阵面推进进行非结构网格结构化处理的方法，可高效生成二维三角形、四边形或三角形-四边形混合非结构网格的局部方向，提高局部方向正交性能，且其中一个局部方向与壁面法向方向接近，有效提高火箭壁面温度梯度的精度。
本发明具体技术方案为：
第一步：生成二维非结构网格。
采用与背景技术1中第一步相同的方法输出非结构网格相关文件，以及相关数据结构初始化。
第二步：根据非结构网格相关数据结构，定义并初始化非结构网格单元的关键元素数组和局部方向数组。
2.1定义关键元素数组其中Elkc为Ckc的关键元素编号组成的数组。本发明中Ckc的关键元素是用于确定局部方向的Ckc的网格边或网格点。具体而言，三角形网格单元的关键元素为3条网格边和1个网格点，如当C1为三角形网格且El1＝{2,4,7,9}时，表示C1的4个关键元素分别为C1的网格边E2、E4、E7以及C1的网格点Pt9；四边形网格单元的关键元素为4条网格边，如当C1为四边形网格且El1＝{2,4,7,9}时，表示C1的4个关键元素分别为构成C1的网格边E2、E4、E7、E9。初始化均为{0,0,0,0}，从而与有意义的单元编号区别以避免歧义；
2.2根据定义局部方向数组locDir，其中LDkc表示网格单元Ckc的局部方向向量组成的数组。LDkc包含4个数据项，分别表示Ckc的4个局部方向向量，由Elkc确定：
·LDkc[1]为Elkc[3]边中点指向Elkc[1]边中点构成的向量；
·若Ckc为四边形，LDkc[2]为Elkc[4]边中点指向Elkc[2]边中点构成的向量；若Ckc为三角形，LDkc[2]为Elkc[4]网格点指向Elkc[2]边中点构成的向量；
·LDkc[3]为Elkc[1]边中点指向Elkc[3]边中点构成的向量；
·若Ckc为四边形，LDkc[4]为Elkc[2]边中点指向Elkc[4]边中点构成的向量；若Ckc为三角形，LDkc[4]为Elkc[2]边中点指向Elkc[4]网格点构成的向量。
由上述定义易知LDkc[1]＝-LDkc[3],LDkc[2]＝-LDkc[4],1≤kc≤Nc。由于局部方向需要在关键元素确定后才能确定，故的每一项都初始化为零向量；
第三步：采用阵面推进确定非结构网格的关键元素和局部方向，即进行非结构网格结构化处理，方法是：
3.1定义并初始化阵面推进相关数据结构，方法是：
3.1.1定义并初始化阵面数组fronts。fronts为动态数组，数据项个数记为Nf，每一项表示已进行访问的网格点的编号。fronts初始化为位于壁面边界(壁面边界指物体表面，如火箭表面等)的网格点的编号数组，Nf初始化为位于壁面边界的网格点的总数；
3.1.2定义并初始化下一阵面数组nextFronts。nextFronts为动态数组，表示需要加入fronts的网格点的编号，nextFronts的数据项个数记为Nnf。nextFronts初始化为空集，Nnf初始化为0；
3.1.3初始化阵面层数组frontsLevel，其中Flkp表示网格点Ptkp的阵面层值。若kp在数组fronts中，初始化Flkp＝0，否则令Flkp＝-1；
3.1.4定义阵面层变量n，表示frontsLevel中所有数据项的最大值，n初始化为0；
3.1.5初始化fronts的三个索引：令头索引head＝0、令尾索引tail＝Nf、令最终尾索引finalTail＝Nf；
3.2以head为fronts索引，获取网格点的编号front[head]，为方便表述，记front[head]＝j,1≤j≤Npt；
3.3以j为网格点的邻接单元数组pCells的索引，获取所有与网格点Ptj相邻的网格单元编号组成的数组PCj，记PCj的数据项总数为PCnum；
3.4初始化网格单元索引i＝1；
3.5以i为PCj的索引，获取第i个与Ptj相邻的网格单元编号PCj[i]，为方便表述，令o＝PCj[i],1≤o≤Nc。若Elo[1]＝0，则Co的关键元素和局部方向还没有确定，转3.6；若Elo[1]≠0，转3.9；
3.6根据网格单元数组cs判断Co的形状，若Co包含三个数据项，说明Co为三角形，转3.7；否则Co包含四个数据项，转3.8；
3.7确定三角形网格单元Co的关键元素，得到Elo，根据Elo计算得到局部向量LDo，并更新nextFronts，具体流程如下：
3.7.1根据网格单元数组cs获得Co三个顶点的编号Co[1]、Co[2]、Co[3]；
3.7.2以Co[1]、Co[2]、Co[3]为阵面层数组frontsLevel索引，获得对应阵面层值
3.7.3对比按照阵面层值从小到大的顺序，将Co[1]、Co[2]、Co[3]重排序，重排序后的Co[1]、Co[2]、Co[3]记为kp1、kp2、kp3，即
3.7.4若转3.7.5；否则转3.7.6；
3.7.5更新Elo[4]＝kp3，Elo[2]为kp1、kp2组成的网格边的编号，Elo[1]为kp1、kp3组成的网格边的编号，Elo[3]为kp2、kp3组成的网格边的编号，转3.7.7；
3.7.6更新Elo[4]＝kp1，Elo[2]为kp2、kp3组成的网格边的编号，Elo[1]为kp1、kp2组成的网格边的编号，Elo[3]为kp1、kp3组成的网格边的编号，转3.7.7；
3.7.7根据Elo，按照第2.2步中局部方向的定义，计算得到LDo，即：
·LDo[1]为Elo[3]边中点指向Elo[1]边中点构成的向量；
·LDo[2]为Elo[4]网格点指向Elo[2]边中点构成的向量；
·LDo[3]为Elo[1]边中点指向Elo[3]边中点构成的向量；
·LDo[4]为Elo[2]边中点指向Elo[4]网格点构成的向量。
3.7.8更新nextFronts，具体步骤为：
3.7.8.1若且kp1不在nextFronts中，将kp1加入nextFronts，Nnf＝Nnf+1；
3.7.8.2若且kp2不在nextFronts中，将kp2加入nextFronts，Nnf＝Nnf+1；
3.7.8.3若且kp3不在nextFronts中，将kp3加入nextFronts，Nnf＝Nnf+1，转3.9；
3.8确定四边形网格单元Co的关键元素，得到Elo，根据Elo计算得到局部向量LDo，并更新nextFronts，具体流程如下：
3.8.1根据网格单元数组cs获得Co四个顶点的编号Co[1]、Co[2]、Co[3]、Co[4]；
3.8.2以Co[1]、Co[2]、Co[3]、Co[4]为frontsLevel索引，获得对应阵面层值
3.8.3对比按照阵面层值从小到大的顺序，将Co[1]、Co[2]、Co[3]、Co[4]重排序，将重排序后的Co[1]、Co[2]、Co[3]、Co[4]记为kp1、kp2、kp3、kp4，即
3.8.4更新Elo[1]为kp1、kp2组成的网格边的编号，Elo[3]为kp3、kp4组成的网格边的编号。若kp1、kp3构成Co的网格边，Elo[2]更新为kp1、kp3组成的网格边的编号，Elo[4]更新为kp2、kp4组成的网格边的编号；否则若kp1、kp4构成Co的网格边，Elo[2]更新为kp1、kp4组成的网格边的编号，Elo[4]更新为kp2、kp3组成的网格边的编号；
3.8.5根据Elo，按照第2.2步中局部方向的定义，计算得到LDo，即
·LDo[1]为Elo[3]边中点指向Elo[1]边中点构成的向量；
·LDo[2]为Elo[4]边中点指向Elo[2]边中点构成的向量；
·LDo[3]为Elo[1]边中点指向Elo[3]边中点构成的向量；
·LDo[4]为Elo[2]边中点指向Elo[4]边中点构成的向量。
3.8.6更新nextFronts，具体步骤为：
3.8.6.1若且kp1不在nextFronts中，将kp1加入nextFronts，Nnf＝Nnf+1；
3.8.6.2若且kp2不在nextFronts中，将kp2加入nextFronts，Nnf＝Nnf+1，转3.9；
3.9令i＝i+1，若i<PCnum+1，转3.5；否则转3.10。
3.10令head＝head+1，若head<tail+1，转3.2；否则转3.11；
3.11根据nextFronts更新fronts、frontsLevel、n、tail及finalTail，具体步骤为：
3.11.1将nextFronts的数据项加入到fronts中，更新Nf＝Nf+Nnf；
3.11.2更新frontsLevel中数据项Flkp＝n+1，其中kp为nextFronts的数据项；
3.11.3更新n＝n+1；
3.11.4更新nextFronts为空集，令Nnf＝0；
3.11.5更新tail＝finalTail＝Nf；
3.12若head<finalTail+1，转3.2；否则所有网格单元的关键元素和局部方向都已确定，得到关键元素数组keyEl和局部方向数组locDir，转3.13；
3.13将关键元素数组keyEl作为关键元素文件输出，将局部方向数组locDir作为局部方向文件输出，实现非结构网格结构化处理。
与现有技术相比，本发明采用阵面推进的方式，从壁面边界向内部区域按序进行局部方向生成，其中一个局部方向与壁面法向方向接近，可有效提高壁面附近温度梯度的精度；本发明对于三角形网格单元也可获得正交性能好的局部方向；本发明的阵面推进过程高效，避免了穷举法的冗余计算。具体原因为：
1)第3.1.1步中，以壁面网格点初始化fronts，第3步通过head、tail和finalTail索引实现了从壁面边界到内部区域的推进过程，frontsLevel记录了推进的顺序，基于frontsLevel确定的局部方向中存在一个方向与壁面法向方向接近。壁面法向方向温度变化快，梯度大，沿着该方向进行温度梯度重构可有效提高温度梯度的精度，使得火箭温度测点的布置更加合理。
2)第2.1和2.2步中对关键元素和局部方向进行了重定义，使得局部方向正交性能更好，对于二维非结构应用中使用广泛的三角形网格，本发明的局部方向也能获得较好的正交性能。
3)第3步通过head、tail和finalTail索引实现的阵面推进过程，没有冗余计算，效率较高。
综上所述，本发明提出了一种采用阵面推进进行非结构网格结构化处理的方法，将该方法处理后的结构化非结构网格用于火箭表面温度梯度计算应用，可在提高计算效率的基础上，进一步提高温度梯度的精度，使得火箭温度测点的布置更加合理。
附图说明
图1是背景技术1的总流程图。
图2是背景技术1第三步的流程图。
图3是本发明的总流程图。
图4是本发明第三步的流程图。
具体实施方式
图3是本发明的总流程图。如图3所示，本发明包括以下步骤：
第一步：生成二维非结构网格。
采用与背景技术1中第一步相同的方法输出非结构网格相关文件，以及相关数据结构初始化。
第二步：根据非结构网格相关数据结构，定义并初始化非结构网格单元的关键元素数组和局部方向数组。
2.1定义关键元素数组其中Elkc为Ckc的关键元素编号组成的数组。本发明中Ckc的关键元素是用于确定局部方向的Ckc的网格边或网格点。具体而言，三角形网格单元的关键元素为3条网格边和1个网格点，如当C1为三角形网格且El1＝{2,4,7,9}时，表示C1的4个关键元素分别为C1的网格边E2、E4、E7以及C1的网格点Pt9；四边形网格单元的关键元素为4条网格边，如当C1为四边形网格且El1＝{2,4,7,9}时，表示C1的4个关键元素分别为构成C1的网格边E2、E4、E7、E9。初始化均为{0,0,0,0}，从而与有意义的单元编号区别以避免歧义；
2.2根据定义局部方向数组其中LDkc表示网格单元Ckc的局部方向向量组成的数组。LDkc包含4个数据项，分别表示Ckc的4个局部方向向量，由Elkc确定：
·LDkc[1]为Elkc[3]边中点指向Elkc[1]边中点构成的向量；
·若Ckc为四边形，LDkc[2]为Elkc[4]边中点指向Elkc[2]边中点构成的向量；若Ckc为三角形，LDkc[2]为Elkc[4]网格点指向Elkc[2]边中点构成的向量；
·LDkc[3]为Elkc[1]边中点指向Elkc[3]边中点构成的向量；
·若Ckc为四边形，LDkc[4]为Elkc[2]边中点指向Elkc[4]边中点构成的向量；若Ckc为三角形，LDkc[4]为Elkc[2]边中点指向Elkc[4]网格点构成的向量。
由上述定义易知LDkc[1]＝-LDkc[3],LDkc[2]＝-LDkc[4],1≤kc≤Nc。由于局部方向需要在关键元素确定后才能确定，故的每一项都初始化为零向量；
第三步：采用阵面推进确定非结构网格的关键元素和局部方向，即进行非结构网格结构化处理，方法如图4所示：
3.1定义并初始化阵面推进相关数据结构，方法是：
3.1.1定义并初始化阵面数组fronts。fronts为动态数组，数据项个数记为Nf，每一项表示已进行访问的网格点的编号。fronts初始化为位于壁面边界的网格点的编号数组，Nf初始化为位于壁面边界的网格点的总数；
3.1.2定义并初始化下一阵面数组nextFronts。nextFronts为动态数组，表示需要加入fronts的网格点的编号，nextFronts的数据项个数记为Nnf。nextFronts初始化为空集，Nnf初始化为0；
3.1.3初始化阵面层数组frontsLevel，其中Flkp表示网格点Ptkp的阵面层值。若kp在数组fronts中，初始化Flkp＝0，否则令Flkp＝-1；
3.1.4定义阵面层变量n，表示frontsLevel中所有数据项的最大值，n初始化为0；
3.1.5初始化fronts的三个索引：令头索引head＝0、令尾索引tail＝Nf、令最终尾索引finalTail＝Nf；
3.2以head为fronts索引，获取网格点的编号front[head]，为方便表述，记front[head]＝j,1≤j≤Npt；
3.3以j为网格点的邻接单元数组pCells的索引，获取所有与网格点Ptj相邻的网格单元编号组成的数组PCj，记PCj的数据项总数为PCnum；
3.4初始化网格单元索引i＝1；
3.5以i为PCj的索引，获取第i个与Ptj相邻的网格单元编号PCj[i]，为方便表述，令o＝PCj[i],1≤o≤Nc。若Elo[1]＝0，则Co的关键元素和局部方向还没有确定，转3.6；若Elo[1]≠0，转3.9；
3.6根据网格单元数组cs判断Co的形状，若Co包含三个数据项，说明Co为三角形，转3.7；否则Co包含四个数据项，转3.8；
3.7确定三角形网格单元Co的关键元素，得到Elo，根据Elo计算得到局部向量LDo，并更新nextFronts，具体流程如下：
3.7.1根据网格单元数组cs获得Co三个顶点的编号Co[1]、Co[2]、Co[3]；
3.7.2以Co[1]、Co[2]、Co[3]为阵面层数组frontsLevel索引，获得对应阵面层值
3.7.3对比按照阵面层值从小到大的顺序，将Co[1]、Co[2]、Co[3]重排序，重排序后的Co[1]、Co[2]、Co[3]记为kp1、kp2、kp3，即
3.7.4若转3.7.5；否则转3.7.6；
3.7.5更新Elo[4]＝kp3，Elo[2]为kp1、kp2组成的网格边的编号，Elo[1]为kp1、kp3组成的网格边的编号，Elo[3]为kp2、kp3组成的网格边的编号，转3.7.7；
3.7.6更新Elo[4]＝kp1，Elo[2]为kp2、kp3组成的网格边的编号，Elo[1]为kp1、kp2组成的网格边的编号，Elo[3]为kp1、kp3组成的网格边的编号，转3.7.7；
3.7.7根据Elo，按照第2.2步中局部方向的定义，计算得到LDo，即：
·LDo[1]为Elo[3]边中点指向Elo[1]边中点构成的向量；
·LDo[2]为Elo[4]网格点指向Elo[2]边中点构成的向量；
·LDo[3]为Elo[1]边中点指向Elo[3]边中点构成的向量；
·LDo[4]为Elo[2]边中点指向Elo[4]网格点构成的向量。
3.7.8更新nextFronts，具体步骤为：
3.7.8.1若且kp1不在nextFronts中，将kp1加入nextFronts，Nnf＝Nnf+1；
3.7.8.2若且kp2不在nextFronts中，将kp2加入nextFronts，Nnf＝Nnf+1；
3.7.8.3若且kp3不在nextFronts中，将kp3加入nextFronts，Nnf＝Nnf+1，转3.9；
3.8确定四边形网格单元Co的关键元素，得到Elo，根据Elo计算得到局部向量LDo，并更新nextFronts，具体流程如下：
3.8.1根据网格单元数组cs获得Co四个顶点的编号Co[1]、Co[2]、Co[3]、Co[4]；
3.8.2以Co[1]、Co[2]、Co[3]、Co[4]为frontsLevel索引，获得对应阵面层值
3.8.3对比按照阵面层值从小到大的顺序，将Co[1]、Co[2]、Co[3]、Co[4]重排序，将重排序后的Co[1]、Co[2]、Co[3]、Co[4]记为kp1、kp2、kp3、kp4，即
3.8.4更新Elo[1]为kp1、kp2组成的网格边的编号，Elo[3]为kp3、kp4组成的网格边的编号。若kp1、kp3构成Co的网格边，Elo[2]更新为kp1、kp3组成的网格边的编号，Elo[4]更新为kp2、kp4组成的网格边的编号；否则若kp1、kp4构成Co的网格边，Elo[2]更新为kp1、kp4组成的网格边的编号，Elo[4]更新为kp2、kp3组成的网格边的编号；
3.8.5根据Elo，按照第2.2步中局部方向的定义，计算得到LDo，即
·LDo[1]为Elo[3]边中点指向Elo[1]边中点构成的向量；
·LDo[2]为Elo[4]边中点指向Elo[2]边中点构成的向量；
·LDo[3]为Elo[1]边中点指向Elo[3]边中点构成的向量；
·LDo[4]为Elo[2]边中点指向Elo[4]边中点构成的向量。
3.8.6更新nextFronts，具体步骤为：
3.8.6.1若且kp1不在nextFronts中，将kp1加入nextFronts，Nnf＝Nnf+1；
3.8.6.2若且kp2不在nextFronts中，将kp2加入nextFronts，Nnf＝Nnf+1，转3.9；
3.9令i＝i+1，若i<PCnum+1，转3.5；否则转3.10。
3.10令head＝head+1，若head<tail+1，转3.2；否则转3.11；
3.11根据nextFronts更新fronts、frontsLevel、n、tail及finalTail，具体步骤为：
3.11.1将nextFronts的数据项加入到fronts中，更新Nf＝Nf+Nnf；
3.11.2更新frontsLevel中数据项Flkp＝n+1，其中kp为nextFronts的数据项；
3.11.3更新n＝n+1；
3.11.4更新nextFronts为空集，令Nnf＝0；
3.11.5更新tail＝finalTail＝Nf；
3.12若head<finalTail+1，转3.2；否则所有网格单元的关键元素和局部方向都已确定，得到关键元素数组keyEl和局部方向数组locDir，转3.13；
3.13将关键元素数组keyEl作为关键元素文件输出，将局部方向数组locDir作为局部方向文件输出，实现非结构网格结构化处理。
Claims (6)
1.一种基于阵面推进的非结构网格结构化处理方法，其特征在于包括以下步骤：
第一步：生成二维非结构网格，方法是：
1.1 将工程应用需要计算的二维空间进行网格化处理，输出二维非结构网格相关文件，工程应用需要计算的空间在指二维计算空间，二维非结构网格相关文件包括网格单元文件、网格点文件、网格单元的邻接单元文件、网格点的邻接单元文件；
1.2 根据非结构网格相关文件初始化相关数据结构，方法是：
1.2.1 令网格单元文件为网格单元数组cs，其中Nc为网格单元总数，Ckc为第kc个网格单元；网格单元采用网格点编号数组表示，Ckc＝{kp1,...,kpw,...,kpW}表示第kc个网格单元由第kp1,...,kpw,...,kpW个网格点构成，是一个W边形网格单元，W为正整数；
1.2.2 令网格点文件为网格点数组pts，其中Npt为网格点总数，Ptkp为第kp个网格点；网格点采用笛卡尔坐标数组表示，Ptkp＝(x,y)表示第kp个网格点的坐标为(x,y)，x,y为实数；
1.2.3 令网格单元的邻接单元文件为数组cCells，其中CCkc表示网格单元Ckc的邻接网格单元编号组成的数组；
1.2.4 令网格点的邻接单元文件为数组pCells，其中PCkp表示网格点Ptkp的邻接网格单元编号组成的数组；
第二步：根据非结构网格相关数据结构，定义并初始化非结构网格单元的关键元素数组和局部方向数组，方法是：
2.1 定义关键元素数组其中Elkc为Ckc的关键元素编号组成的数组，Ckc的关键元素是用于确定局部方向的Ckc的网格边或网格点；三角形网格单元的关键元素为3条网格边和1个网格点，四边形网格单元的关键元素为4条网格边，初始化均为{0,0,0,0}；
2.2 根据keyEl定义局部方向数组locDir，其中LDkc表示网格单元Ckc的局部方向向量组成的数组，LDkc包含4个数据项，分别表示Ckc的4个局部方向向量，由Elkc确定：
LDkc[1]为Elkc[3]边中点指向Elkc[1]边中点构成的向量；
若Ckc为四边形，LDkc[2]为Elkc[4]边中点指向Elkc[2]边中点构成的向量；若Ckc为三角形，LDkc[2]为Elkc[4]网格点指向Elkc[2]边中点构成的向量；
LDkc[3]为Elkc[1]边中点指向Elkc[3]边中点构成的向量；
若Ckc为四边形，LDkc[4]为Elkc[2]边中点指向Elkc[4]边中点构成的向量；若Ckc为三角形，LDkc[4]为Elkc[2]边中点指向Elkc[4]网格点构成的向量；
将的每一项都初始化为零向量；
第三步：采用阵面推进确定非结构网格的关键元素和局部方向，即进行非结构网格结构化处理，方法是：
3.1 定义并初始化阵面推进相关数据结构，方法是：
3.1.1 定义并初始化阵面数组fronts，fronts为动态数组，数据项个数记为Nf，每一项表示已进行访问的网格点的编号；fronts初始化为位于壁面边界即物体表面的网格点的编号数组，Nf初始化为位于壁面边界的网格点的总数；
3.1.2 定义并初始化下一阵面数组nextFronts，nextFronts为动态数组，表示需要加入fronts的网格点的编号，nextFronts的数据项个数记为Nnf。nextFronts初始化为空集，Nnf初始化为0；
3.1.3 初始化阵面层数组frontsLevel，其中Flkp表示网格点Ptkp的阵面层值，若kp在数组fronts中，初始化Flkp＝0，否则令Flkp＝-1；
3.1.4 定义阵面层变量n，表示frontsLevel中所有数据项的最大值，n初始化为0；
3.1.5 初始化fronts的三个索引：令头索引head＝0、令尾索引tail＝Nf、令最终尾索引finalTail＝Nf；
3.2 以head为fronts索引，获取网格点的编号front[head]，记front[head]＝j,1≤j≤Npt；
3.3 以j为网格点的邻接单元数组pCells的索引，获取所有与网格点Ptj相邻的网格单元编号组成的数组PCj，记PCj的数据项总数为PCnum；
3.4 初始化网格单元索引i＝1；
3.5 以i为PCj的索引，获取第i个与Ptj相邻的网格单元编号PCj[i]，令o＝PCj[i],1≤o≤Nc，若Elo[1]＝0，则Co的关键元素和局部方向还没有确定，转3.6；若Elo[1]≠0，转3.9；
3.6 根据网格单元数组cs判断Co的形状，若Co包含三个数据项，说明Co为三角形，转3.7；否则Co包含四个数据项，转3.8；
3.7 确定三角形网格单元Co的关键元素，得到Elo，根据Elo计算得到局部向量LDo，并更新nextFronts，转3.9；
3.8 确定四边形网格单元Co的关键元素，得到Elo，根据Elo计算得到局部向量LDo，并更新nextFronts，转3.9；
3.9 令i＝i+1，若i<PCnum+1，转3.5；否则转3.10。
3.10 令head＝head+1，若head<tail+1，转3.2；否则转3.11；
3.11 根据nextFronts更新fronts、frontsLevel、n、tail及finalTail，具体步骤为：
3.11.1 将nextFronts的数据项加入到fronts中，更新Nf＝Nf+Nnf；
3.11.2 更新frontsLevel中数据项Flkp＝n+1，其中kp为nextFronts的数据项；
3.11.3 更新n＝n+1；
3.11.4 更新nextFronts为空集，令Nnf＝0；
3.11.5 更新tail＝finalTail＝Nf；
3.12 若head<finalTail+1，转3.2；否则所有网格单元的关键元素和局部方向都已确定，即得到关键元素数组keyEl和局部方向数组locDir，转3.13；
3.13 将关键元素数组keyEl作为关键元素文件输出，将局部方向数组locDir作为局部方向文件输出，实现非结构网格结构化处理。
2.如权利要求1所述的基于阵面推进的非结构网格结构化处理方法，其特征在于1.1步所述将工程应用需要计算的二维空间进行网格化处理的方法是使用专业软件，专业软件指版本18.0以上的Fluent软件、版本13.0以上的FUN3D软件、版本3.0以上的OpenFOAM、版本5.0以上的SU2。
3.如权利要求1所述的基于阵面推进的非结构网格结构化处理方法，其特征在于3.7步所述确定三角形网格单元Co的关键元素Elo和局部方向LDo，并更新nextFronts的流程为：
3.7.1 根据网格单元数组cs获得Co三个顶点的编号Co[1]、Co[2]、Co[3]；
3.7.2 以Co[1]、Co[2]、Co[3]为阵面层数组frontsLevel索引，获得对应阵面层值
3.7.3 对比按照阵面层值从小到大的顺序，将Co[1]、Co[2]、Co[3]重排序，重排序后的Co[1]、Co[2]、Co[3]记为kp1、kp2、kp3，即
3.7.4 若转3.7.5；否则转3.7.6；
3.7.5 更新Elo[4]＝kp3，Elo[2]为kp1、kp2组成的网格边的编号，Elo[1]为kp1、kp3组成的网格边的编号，Elo[3]为kp2、kp3组成的网格边的编号，转3.7.7；
3.7.6 更新Elo[4]＝kp1，Elo[2]为kp2、kp3组成的网格边的编号，Elo[1]为kp1、kp2组成的网格边的编号，Elo[3]为kp1、kp3组成的网格边的编号，转3.7.7；
3.7.7 根据Elo，按照第2.2步中局部方向的定义，计算得到LDo；
3.7.8 更新nextFronts，具体步骤为：
3.7.8.1 若且kp1不在nextFronts中，将kp1加入nextFronts，Nnf＝Nnf+1；
3.7.8.2 若且kp2不在nextFronts中，将kp2加入nextFronts，Nnf＝Nnf+1；
3.7.8.3 若且kp3不在nextFronts中，将kp3加入nextFronts，Nnf＝Nnf+1。
4.如权利要求3所述的基于阵面推进的非结构网格结构化处理方法，其特征在于3.7.7步所述根据Elo计算得到LDo的方法是：
LDo[1]为Elo[3]边中点指向Elo[1]边中点构成的向量；
LDo[2]为Elo[4]网格点指向Elo[2]边中点构成的向量；
LDo[3]为Elo[1]边中点指向Elo[3]边中点构成的向量；
LDo[4]为Elo[2]边中点指向Elo[4]网格点构成的向量。
5.如权利要求1所述的基于阵面推进的非结构网格结构化处理方法，其特征在于3.8步所述确定四边形网格单元Co的关键元素Elo和局部方向LDo，并更新nextFronts的流程如下：
3.8.1 根据网格单元数组cs获得Co四个顶点的编号Co[1]、Co[2]、Co[3]、Co[4]；
3.8.2 以Co[1]、Co[2]、Co[3]、Co[4]为frontsLevel索引，获得对应阵面层值
3.8.3 对比按照阵面层值从小到大的顺序，将Co[1]、Co[2]、Co[3]、Co[4]重排序，将重排序后的Co[1]、Co[2]、Co[3]、Co[4]记为kp1、kp2、kp3、kp4，即
3.8.4 更新Elo[1]为kp1、kp2组成的网格边的编号，Elo[3]为kp3、kp4组成的网格边的编号；若kp1、kp3构成Co的网格边，Elo[2]更新为kp1、kp3组成的网格边的编号，Elo[4]更新为kp2、kp4组成的网格边的编号；否则若kp1、kp4构成Co的网格边，Elo[2]更新为kp1、kp4组成的网格边的编号，Elo[4]更新为kp2、kp3组成的网格边的编号；
3.8.5 根据Elo，按照第2.2步中局部方向的定义，计算得到LDo；
3.8.6 更新nextFronts，具体步骤为：
3.8.6.1 若且kp1不在nextFronts中，将kp1加入nextFronts，Nnf＝Nnf+1；
3.8.6.2 若且kp2不在nextFronts中，将kp2加入nextFronts，Nnf＝Nnf+1。
6.如权利要求5所述的基于阵面推进的非结构网格结构化处理方法，其特征在于3.8.5步所述根据Elo计算得到LDo的方法是：
LDo[1]为Elo[3]边中点指向Elo[1]边中点构成的向量；
LDo[2]为Elo[4]边中点指向Elo[2]边中点构成的向量；
LDo[3]为Elo[1]边中点指向Elo[3]边中点构成的向量；
LDo[4]为Elo[2]边中点指向Elo[4]边中点构成的向量。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910853374.0A CN110610042A (zh) | 2019-09-10 | 2019-09-10 | 基于阵面推进的非结构网格结构化处理方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910853374.0A CN110610042A (zh) | 2019-09-10 | 2019-09-10 | 基于阵面推进的非结构网格结构化处理方法 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN110610042A true CN110610042A (zh) | 2019-12-24 |
Family
ID=68891131
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201910853374.0A Pending CN110610042A (zh) | 2019-09-10 | 2019-09-10 | 基于阵面推进的非结构网格结构化处理方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN110610042A (zh) |
Cited By (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111223185A (zh) * | 2019-12-30 | 2020-06-02 | 杭州数孪科技有限公司 | 三维模型的网格单元信息展示方法，其装置及电子设备 |
| CN112989673A (zh) * | 2021-04-16 | 2021-06-18 | 中国空气动力研究与发展中心计算空气动力研究所 | 一种基于二叉树模型的离散网格点快速重规整化方法 |
| CN114444274A (zh) * | 2022-01-05 | 2022-05-06 | 中国空气动力研究与发展中心计算空气动力研究所 | 一种从非结构网格中重构原结构网格的方法、介质及装置 |
| CN114549793A (zh) * | 2022-04-28 | 2022-05-27 | 中国空气动力研究与发展中心计算空气动力研究所 | 一种从二维非结构网格重构结构网格的方法、介质及装置 |
| CN116227043A (zh) * | 2023-05-10 | 2023-06-06 | 中国空气动力研究与发展中心计算空气动力研究所 | 一种飞行器数值模拟方法、系统、设备及计算机存储介质 |
-
2019
- 2019-09-10 CN CN201910853374.0A patent/CN110610042A/zh active Pending
Cited By (9)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111223185A (zh) * | 2019-12-30 | 2020-06-02 | 杭州数孪科技有限公司 | 三维模型的网格单元信息展示方法，其装置及电子设备 |
| CN111223185B (zh) * | 2019-12-30 | 2023-10-27 | 苏州数设科技有限公司 | 三维模型的网格单元信息展示方法，其装置及电子设备 |
| CN112989673A (zh) * | 2021-04-16 | 2021-06-18 | 中国空气动力研究与发展中心计算空气动力研究所 | 一种基于二叉树模型的离散网格点快速重规整化方法 |
| CN112989673B (zh) * | 2021-04-16 | 2022-03-15 | 中国空气动力研究与发展中心计算空气动力研究所 | 一种基于二叉树模型的离散网格点快速重规整化方法 |
| CN114444274A (zh) * | 2022-01-05 | 2022-05-06 | 中国空气动力研究与发展中心计算空气动力研究所 | 一种从非结构网格中重构原结构网格的方法、介质及装置 |
| CN114444274B (zh) * | 2022-01-05 | 2023-05-02 | 中国空气动力研究与发展中心计算空气动力研究所 | 一种从非结构网格中重构原结构网格的方法、介质及装置 |
| CN114549793A (zh) * | 2022-04-28 | 2022-05-27 | 中国空气动力研究与发展中心计算空气动力研究所 | 一种从二维非结构网格重构结构网格的方法、介质及装置 |
| CN116227043A (zh) * | 2023-05-10 | 2023-06-06 | 中国空气动力研究与发展中心计算空气动力研究所 | 一种飞行器数值模拟方法、系统、设备及计算机存储介质 |
| CN116227043B (zh) * | 2023-05-10 | 2024-03-12 | 中国空气动力研究与发展中心计算空气动力研究所 | 一种飞行器数值模拟方法、系统、设备及计算机存储介质 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN110610042A (zh) | 基于阵面推进的非结构网格结构化处理方法 | |
| Kedward et al. | Efficient and exact mesh deformation using multiscale RBF interpolation | |
| Du et al. | Meshfree, probabilistic determination of point sets and support regions for meshless computing | |
| CN112862972B (zh) | 一种表面结构网格生成方法 | |
| Lehto et al. | A radial basis function (RBF) compact finite difference (FD) scheme for reaction-diffusion equations on surfaces | |
| Ronchi et al. | The “cubed sphere”: A new method for the solution of partial differential equations in spherical geometry | |
| Lin et al. | Constructing iterative non-uniform B-spline curve and surface to fit data points | |
| Zhang et al. | Dynamic load balancing based on constrained kd tree decomposition for parallel particle tracing | |
| CN110633149B (zh) | 均衡非结构网格单元计算量的并行负载均衡方法 | |
| CN109726441B (zh) | 体和面混合gpu并行的计算电磁学dgtd方法 | |
| CN110516316B (zh) | 一种间断伽辽金法求解欧拉方程的gpu加速方法 | |
| CN107403466A (zh) | 基于全局加密的超大规模非结构网格生成方法 | |
| CN110663064A (zh) | 用于矢量图形和图像处理的并行化流水线 | |
| Lewis et al. | Scattered data interpolation and approximation for computer graphics | |
| Li et al. | PSNet: Fast data structuring for hierarchical deep learning on point cloud | |
| CN111489447A (zh) | 一种适用于格子Boltzmann方法的直角网格自适应建模方法 | |
| Mueller‐Roemer et al. | Ternary sparse matrix representation for volumetric mesh subdivision and processing on GPUs | |
| Ansari et al. | Mesh partitioning and efficient equation solving techniques by distributed finite element methods: A survey | |
| Rahn et al. | Scalable distributed-memory external sorting | |
| CN106484532B (zh) | 面向sph流体模拟的gpgpu并行计算方法 | |
| Feng et al. | Cellular topology optimization on differentiable Voronoi diagrams | |
| Ferretti | The cell method: A purely algebraic computational method in physics and engineering | |
| Ling et al. | Vina-FPGA: A hardware-accelerated molecular docking tool with fixed-point quantization and low-level parallelism | |
| Sharma et al. | Multiple K means++ clustering of satellite image using Hadoop MapReduce and Spark | |
| WO2023216915A1 (zh) | 一种基于图形处理器的直升机流场数值模拟系统及方法 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| WD01 | Invention patent application deemed withdrawn after publication | ||
| WD01 | Invention patent application deemed withdrawn after publication |
Application publication date: 20191224 |