基于Laplacian能量的CAD模型四六面体网格拟合算法研究

计算设计与工程杂志，卷。号12（2014）88~95www.jcde.org基于Laplacian能量Yusuke Imai*、Hiroyuki Hiraoka和Hiroshi Kawaharada112-8551东京都文京区春日1-13-27中央大学精密机械系（2013年9月16日接收;2013年10月24日修订;2013年11月1日接受摘要由于使用实际产品进行性能测试的成本很高，因此制造商使用成本较低的计算机辅助设计仿真来实现此功能。在本文中，我们建议使用六面体网格，这是更准确的比四面体网格，有限元分析。我们提出了自动六面体网格生成尖锐的功能，以精确地表示目标形状的相应功能我们的六面体网格生成使用基于体素的算法。在我们以前的工作中，我们使用拉普拉斯能量最小化将体素的表面拟合到目标表面我们在拟合中使用法向量来保留尖锐的特征。然而，这种方法不能精确地表示凹尖特征在这个建议中，我们改进了我们以前的拉普拉斯能量最小化，增加了一项，依赖于多个正常的向量，而不是使用正常的向量。此外，我们强调一个凸/凹表面子集，以表示凹尖锐的功能。关键词：CAD模型;六面体网格;尖锐特征;拟合算法;多法向量1. 介绍在制造业中，计算机模拟的成本低于测试实际原型。因此，大多数制造商运行仿真，这需要体积网格。直到十年前，模拟过程还从使用计算机辅助设计（CAD）软件制作的然而，实际产品和CAD模型之间的形状存在差异，这是由于生产压力机中的回弹等制造因素造成的。即使制造商使用与CAD模型中形状相同的金属模具，通过冲压加工等工艺获得的零件也不具有与CAD模型相同的因此，实际产品与CAD模型不同。对于真实的模拟，CAD模型必须与产品完全相同。如今，仿真过程从实际产品扫描的点云开始这个过程被称为逆向工程。从这样的点云生成的四面体/六面体网格六面体网格很重要，因为它们在精确分析方面优于四面体网格（见图1（a））。因此，在本文中，我们使用六-*通讯作者。联系电话：+81-3-3817-1821，传真：+81-3 -3817-1820电子邮件地址：imai@lcps.mech.chuo-u.ac.jp© 2014 CAD/CAM工程师协会Techno-Press doi：10.7315/JCDE. 2014. 009六面体体积网格（见图1（b）），其单元仅为六面体单元（称为全六面体网格）并考虑其表面。在两个物体之间的结构分析中，峰值应力发生在接触区域附近或周围这些区域通常是尖锐的特征。尖锐特征通常是对象的尖点部分（诸如边缘或点）。因此，为了获得精确的模拟结果，表面网格必须表示尖锐的特征。在本文中，我们考虑一个基于体素的六面体网格生成算法[1，2]。此外，有限元的体积网格雅可比矩阵是一个三重标量积a（ b × c），其中a、b和c是与单元的角顶点相邻的向量（边）。我们的最终目标是自动六面体网格生成没有负雅可比矩阵。在本文中，我们讨论了六面体网格的四边形表面（我们在以前的研究中研究的算法是基于体素的[1，2]）。因此，我们的输入是目标表面网格和作为体素网格表面的四边形网格一般来说，我们可以将边界拟合六面体网格的方法分类为基于体素的[3，4]，推进前沿[5-7]，晶须波动[8]，循环消除[9]，中间基于轴的[10]和扫描/映射方法[11，12]，尽管存在其他类型。六面体网格生成算法可以是全自动或半自动的，但没有一种方案能保证所有的雅可比矩阵都是正的。另Y. Imai等人/Journal of Computational Design and Engineering，Vol.号12（2014）88~9589（一）（b）第（1）款图1.四面体/六面体网格：（a）四面体网格，（b）六面体网格。一方面，有一个计划，保证正雅可比- ans确实存在四面体网格。在考虑如何实现正雅可比矩阵之前，我们生成了代表目标曲面的全六面体网格的四边形曲面本文在文献[1，2，13]的基础上，提出了一种具有尖锐特征的自动拟合算法2. 以前的六面体网格划分首先，我们总结了以前的六面体网格生成算法[1，2，13]。基本算法可以分解如下。1. 输入目标表面网格（参见图2（a）[14]）。2. 使用Polymender [15]生成体素以包裹目标表面（参见图2（b））。3. 提取体素的边界面。4. 拟合体素的边界表面（参见图2（c）[1，2]和图2（d）[13]中的拟合表面）。5. 确定内部顶点的位置。6. 应用后处理。7. 输出六面体网格。（一）（b）第（1）款（c）第（1）款90Y. Imai等人/Journal of Computational Design and Engineering，Vol.号12（2014）88~95（d）其他事项图2.齿轮：（a）目标表面（三角形网格），（b）齿轮体素，（c）无尖锐特征的四边形网格，（d）具有尖锐特征。Y. Imai等人/Journal of Computational Design and Engineering，Vol.号12（2014）88~9591M我我我我min{c}（vk+1-1（k+1）k+11i我我|(i,j)|J（i，j） ∈E- （vk-12∑（i，j） ∈Evk）<$我|（i，j）|J+c vk+1- v k2 + c vk+1- wk+ 122i i3i i- （k+1）（ihk+1+c4{∑（h，t）（h）我J（一）图3.带反向部分的螺栓。这里，c = 1，c =1，c =2，c=0.01。12234在上面的步骤4中，我们将体素的边界表面拟合到目标表面网格。我们称之为“装配部分”。我们提出了一种拟合算法，它不代表尖锐的特征（见图2（c））[1，2]。因此，我们提出了一种改进的拟合算法[13]。修改后的算法表示凸尖特征（见图2（d））。但表面拉普拉斯平滑用于防止顶点在目标表面的尖锐特征附近变得拥挤然而，普通的拉普拉斯平滑消除了尖锐的特征。因此，我们使用特征保持拉普拉斯平滑，而不是普通的拉普拉斯平滑。上述最小化中的第一项是连续拉普拉斯算子之间第二项是vk+1和vk之间的差。第三个任期是不能代表凹的尖锐特征。因此，我们i ivk+1和wk+1的差。 第四个任期是现在提出了一种新的拟合算法，法向量和凸/凹加重。多法向量k，h = 1，2，3 <$;nwk+1，t = 1，2，3 <$。 我们描述伊赫它3. 拟合算法多法向量和不匹配率，第4款.在这一项中，多法线向量k，h =在本节中，我们将解释新的拟合算法。让v0，i = 1，2，v是体素边界的位置1，2，3，n是vk+1的一组线性公式它(not单位向量我第0步的顶点（初始位置）。新的拟合算法是迭代的（k= 0 ~T）。算法一：1. 确定vk+1，i = 1，2，k，通过最小化以下公式确定：tors）。在以前的研究中，我们使用了类似的最小化（方程。(2))[13]第10段。1imin{c（vk+1-（k+1）汀江 二次 能源 功能 [16]得双曲余切值.v k表示在步骤k处边界顶点v i的位置，并且wk+1表示输入表面上与v k对应的点。k+11i我我|(i,j)|J（i，j）∈E22. 将表面拉普拉斯平滑应用于第（k+1）个- （vk-1∑（i，j）∈Evk）<$面如果k G T，则将k增加1并返回步骤1。设E为边界曲面的边界边我|（i，j）|J+C vk+1- v k2 +c vk+1- wk+1体素， （i，j）在...之间占优势vi 而vj|是连接到v i的边的数量，|be the number of edgesconnecting to v i, and nv i2i i3i i2和nwk+1是第k个上vk+c4 {k+1-（nk+1 n wk+1）2}}（2）ik+1i我我我在目标表面上的边界表面和Wi，- 是的 表示向量的范数。 一个巨大的-注意向量a和b之间的内积。 T是用户定义的最大迭代次数（本文中T = 30）。二次能量函数相应的第四项Eq. (2)是类似地表示法向量nwk+1和nwk+1之间的不匹配率的二次公式（非线性）。怎么--问题如下：ii曾经，Eq.的第四项的值（2）不改变用-k+1代替k+1因此，让一个正常我我vv92Y. Imai等人/Journal of Computational Design and Engineering，Vol.号12（2014）88~95我我我我我我我我我我ih我k+1 kih我我我我我我ih×1我我我我向量n是方程的相应第四项的最小值（二）、那么，-n也是一个极小值，因为相应的第四项是二次的。图3显示了结果网格与反向部分。因此，我们将c4设置为非常小。在这种情况下，所得网格不代表凸/凹尖锐特征（参见图5（c））。图5（a）和5（b）示出了目标表面和螺栓的体素。因此，我们定义了一个新的第四项方程。（1）作为线性公式。在一种方法[16]中，wk+1是上最接近vk的点。连接到vk的三角形（面）可能不平行于连接到wk+1的每个三角形（面）的法向量。全六面体网格的边界表面的面是四边形的。我们考虑一个有四个三角形和四个法向量的四边形面因此，我们将vk的多重法向量设置为此类三角形对vk的法向量。我们将这个多法线向量的几个并行向量组合成一个向量，并重复应用此操作。因此，该多法线向量没有子平行法向量的集合例如，如果所有三角形我我目标曲面网格。然而，在其他文献[1，2]中，w k+1是目标曲面网格上最近的顶点。 这两个定义的wk +1是明显不同的。我们通过上采样增加了目标表面网格上的顶点数量[1，2]. 我们选择细分三角形进行上采样。 与循环细分方法一样，我们将目标曲面上的每个三角形细分为四个三角形。 我们将这种细分应用了五次。设Sj，j = 0，1，2，…相邻于vk的元素在一个平面上，则vk的多法向量只有一个元素是该平面的法向量。我们设k，h = 1，2，3，k为这个收缩的多重法向量. 我们对Sj应用类似的过程。类似地，我们将多重法向量ns_j，t = 1，2，3，…，设为包含s_j的三角形的法向量。 设A k，B j是在上述条件下，目标曲面上的顶点（包括上采样）让是目标表面上Sj处的单位法向量。本文中的wk+1的定义是：α∑（h，t）（knsjt）w arg min.{v-sj（1 -）}过程4.1 凸/凹加重在本小节中，我们在第二部分（考虑法向量）中为第3ijiKi（三）以表示凹的尖锐特征。具体来说，在算法1的步骤1中确定wk+1之前，我们对边界的子集添加凸/凹加重此定义基于数据库之间的不匹配率和NSJ.此时，α = 0.8。wk+1的旧定义是[1，2]：w k +1最小值为10000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000全六面体网格的表面为了表示尖锐特征，连接对应于尖锐特征的尖锐边界边缘的vk+1必须在对应的尖锐边缘上。因此，重要的是，相应的wk+1在尖锐的边缘上因此，我们认为，iji我K+1在最小化之前（决定wi），我们强调W k+1的另一个旧定义 是[13]：wk+1是一个简单的函数。{vk-sj（1 -αnvknsj）}（5）J候选顶点sj的法向量nsj更接近于wk，并且顶点sj倾向于被选择为wk+1。全六面体网格边界曲面的凸/凹子集图4示出了凸/凹加速度。黑线表示目标表面。红线表示全六面体网格的边界曲面在这种情况下，红色顶点不在黑线的凸角和凹角上，因此边界曲面ii不代表尖锐的特征。新的拟合算法由两部分组成。在第一部分中，我们应用拟合算法与方程。（一）. 我们为了表示尖锐的特征，我们首先计算Laplaci-边界表面的vk让CCk成为向量-设c4 = 0.0，T = 2。在第二部分中，我们应用了ikik算法与Eq. (1)用c= 0.01，将K重置为0，平行于Δih的平均值的tor4T = 30（ K = 0 ~ 30）。其范数等于多法向量的平均值与拉普拉斯向量之间的内积。如果内积为负，CCk的方向相反，4. 表示凹尖特征的拟合算法我站点的平均值为0.01k，h = 1，2，3，0.01k，Ak。设cl为在本节中，我们将解释新拟合算法中的多法向量和凸/凹强调初始体素的边缘的长度，并且D是用户-定义的系数（在本文中，我们将其设置为0.5）。我们的凸/凹加重是vk += D × cl ×首先，我们解释了不匹配率和相应的点w k+1 关于V K 基于多法线向量。 中CCkK/3（）下一页2我，其中i = 1，2，3，如果K是3的倍数。我我在以前的研究[13]中，我们定义了平均值（aver-k）和平均值（aver-k +1）。如果K不是3的倍数，我们就不应用重音。我我步骤K时三角形法向量的年龄 和步骤K+1，分别。在这里，我们意识到，即使阿斯克在不依赖于法向量的第一部分v T- w T 几乎为0，因为第三个我我我与nwk+1相同的方向，每个在等式中， (1)几乎变成了0。我我我NY. Imai等人/Journal of Computational Design and Engineering，Vol.号12（2014）88~9593我我我我ihihihihih我我我伊赫伊iji图4.凸/凹加重。表1.法向量匹配表。克I1克I2克I3…克尔克伊什基MaxNSJ10.100.96-0.020.380.96NSJ20.75-0.100.11-0.290.82……nsjBj0.050.800.360.550.80Max0.750.960.99…0.759.64接下来，我们解释Eq的第四项（一）. wk+1是{sj}j=1，2，…的元素.因此，我们有一个多法向量nwk+1和针对wk+1的法向量匹配表。让它我Ck+1 是nwk+1的法向量的数目。有一个我我在第二部分中，没有凸/凹重音，多重法向量nwk+1，t = 1，2，3，n，Ck+1.它我对于所有K，wk+1不变，因为内部在等式的右侧最小化的项（3）总是为了将第四项定义为vk，我们定义∑（h，t）（nk+1 <$nwk+1）为我爱它在相同的sj下几乎为0。因此，我们使用了nwk +1（线性公式）和nwk +1之间的凸/凹伊赫它强调将vk与wk+1分开。在这个重音之后-在法向量匹配表中的彩色单元格处，我我因此，wk +1在拐角（尖锐特征）上的概率增加。4.2 多法线间的匹配率wk+1。（在蓝色单元格中，我们将内积加到和上两次）。设k+1是对应于以下的两倍三角形面积：k，Mk（= Ak + Ck+1）是行数，在本节中，我们解释了多个之间的匹配率伊伊伊K+1t-法向量k，h = 1，2，3，k，Ak和Ak，列的法线向量匹配表中。NSJTih，t = 1，2，3，λ，Bj.ii在最小化过程中，k+1是常数。在第二部分之后我们进行后处理，以提高质量的为了做到这一点，我们计算所有的内积是-在多重法向量之间，h = 1，2，3，k，Ak和得到全六面体网格。伊赫伊ns jt，t = 1，2，3，n，B j（见表1）。内积的值域是[-1，1]，因为nik和ns jt是单位向量。首先，我们得到每一行和每一列的最大值。 我们将∑（h，t）（kns jt）定义为最大值之和（表1右下角的单元格）。 在表1中，我们显示了匹配。红色/绿色单元格是行/列的最大单元格。蓝色单元格是行和列的最大单元格右下角的单元格显示行和列最大值的总和。我们使用该值作为在k，h = 1，2，3，k，Ak和nsjt，t = 1，2，3，k，Bj之间的匹配率。设 Nk（=Ak+Bj）为表1中的行数和列数缩放的匹配率∑（h，t）（k_k_ns_jt）/N_k为[-1，1]。4.3 基于多法向量的拟合算法在本节中，我们解释了Eq。（1）Eq.（3）详细。由方程式（3），我们使用多法向向量nk和ns之间的缩放匹配率。本文中wk+1的定义反映了目标表面的细节5. 结果在本节中，我们展示了我们方法的实验结果。图5（d）显示了通过我们的方法计算的全六面体网格（螺栓模型）的边界四边形表面网格图6（a）、6（b）、6（c）和6（d）分别是目标表面、使用先前方法[13]的轴承模型的边界四边形表面网格、体素和通过我们的方法计算的轴承模型的全六面体网格螺栓和轴承模型的凹尖特征被成功地表示。类似地，我们的方法可以表示目标表面网格的凸尖锐特征图7中的平滑特征模型有一个逐渐减小的凸尖特征图7（a）、7（b）和7（c）分别是目标表面、体素和我们的方法针对平滑特征模型的结果我们实验的计算时间只有几秒钟。Letn^b是体素的边界的矢量的数量，并且N^b是目标表面的矢量（包括我们解决94Y. Imai等人/Journal of Computational Design and Engineering，Vol.号12（2014）88~95（一）（b）第（1）款（c）第（1）款（d）其他事项图5.螺栓：（a）目标表面，（b）螺栓体素，（c）使用先前方法[13]的全六面体网格的边界四边形表面，（d）使用我们的方法的全六面体网格。Y. Imai等人/Journal of Computational Design and Engineering，Vol.号12（2014）88~9595（一）（b）第（1）款（c）第（1）款（d）其他事项图6.轴承：（a）目标表面，（b）边界四边形表面没有，尖锐的特点，（c）体素，（d）边界四边形表面使用我们的方法96Y. Imai等人/Journal of Computational Design and Engineering，Vol.号12（2014）88~95（一）（b）第（1）款（c）第（1）款图7.平滑特征：（a）目标表面，（b）体素，（c）我们方法的结果。Y. Imai等人/Journal of Computational Design and Engineering，Vol.号12（2014）88~9597在拟合部分的线性系统，计算时间为O（n^2）。我们如何解决这个问题是因为这个线性系统的矩阵是稀疏的。所有算法的计算时间都是O（n^N^）。因此，使用octree，我们将cree视为实际计算时间。6. 结论提出了一种新的基于多法向量和凸/凹强调的拟合算法，以表示目标表面上的凸/凹尖锐特征将该算法用于六面体网格的自动生成，得到了具有凸/凹尖锐特征的全六面体网格边界曲面在未来的工作中，我们打算确定我们的算法的最佳系数除了表面网格，我们想扩大这一研究的内部顶点。致谢这项工作得到了JSPS KAKENHI资助号24760124的支持。引用[1] Kawaharada H，Hiraoka H.数值模拟中体积细分的边界表示In：Proceedings of the ACDDE; 2011 Aug 27-29;Shanghai，China; p.三比八[2] 作者：Kawaharada H，Sugihara K.基于细分的六面体网格生成。计算工程2011; 16（2）：12-15.[3] Schneiders R，Schindler R，Weiler F.基于八叉树的六面体单元网格生成。在：第五届国际网格- ING圆桌会议; 1996年10月10日至11日;匹兹堡，宾夕法尼亚州;p.205-215[4] Loic M.一种基于八叉树的六面体网格划分新方法在：第10届国际网格圆桌会议; 2001年10月10日至11日;纽波特海滩，加利福尼亚州; p.209-221[5] Saten ML，Owen SJ，Blacker TD.无约束铺层与剖分：一种全六面体网格生成的新思想在：第14届国际网格圆桌会议; 2005年9月11日至14日;圣地亚哥，加利福尼亚州; p.399-416[6] 黄文忠，陈晓梅，陈文忠.第二部分介绍了智能局部法网格自动生成方法在六面体网格自动生成中的应用。日本机械工程师学会会刊系列2001; 67（655）：496-502.[7] 布莱克TD，迈耶斯RJ。三维六面体网格生成算法。工程与计算机1993; 2（9）：83-93.[8] Tautges TJ，Blacker T，Mitchell SA.晶须编织算法：一种基于连通性的全六面体有限元网格构建方法。国际工程数值方法1996; 39：3327-3349.[9] 米勒-汉内曼六面体网格剖分的逐次对偶循环消元法。In ： The 7th International Mesh- ing Roundtable; 1998Oct 26-28; Dearborn，MI; p.365-378.98Y. Imai等人/Journal of Computational Design and Engineering，Vol.号12（2014）88~95[10] 谢佛A，艾齐安M，拉波波特A，柏科维尔M.基于嵌入式Voronoi图的六面体网格生成。在：第七届国际啮合圆桌会议; 1998年10月26日至28日;亲爱的出生，MI; p.347-364[11] Jason S，Mitchell SA，Knupp P，White D.多扫描自动化方法。In：第九届国际啮合圆桌会议; 2000年10月2日至5日;新奥尔良，CA; p.七十七比八十七[12] Xevi R，Sarrate J.扫描网格的自动和通用最小二乘投影程序In ： The 15th Interna- tional Meshing Roundtable; 2006Sep 17-20; Birmingham，AB; p.487-506[13] 作者：Moriya S，Hiraoka H，Kawaharada H.使用改进的拉普拉斯能量最小化的尖锐特征六面体网格的复制。第三届亚洲设计与数字工程会议论文集; 2012年12月6日至8日;日本北海道;论文号100062。[14] [互联网]第二届中国国际汽车工业展览会达姆施塔特：德国; [引自2012年7月11日]。可从以下网址获得：http://shapes.aimatshape.net/[15]Polymender[Internet]. 西雅图;[引用2012年7月11日]。可从以下网址获得：http://www1.cse.wustl.edu/~taoju/code/polymender.htm[16] MehraR，ZhouQ，LongJ，ShefferA，Gooch A，Mitra NJ.抽象的人造形状。 ACMTransactionsonGraphics。2009;28（5）：1-10.