论文部分内容阅读
摘要： 针对六面体网格自动划分的难度远高于四面体网格的问题，用投影法对简单形状的初始网格进行投影变换得到每个块体的实际网格，用节点合并算法和再分割技术实现不同块体之间的不同密度网格的过渡，从而生成复杂三维几何形体的全六面体网格.算例表明：该方法生成的网格质量很好，易于实施，适用性广.
关键词： 全六面体网格； 投影法； 超限插值； 再分割技术； 网格过渡； 爆炸模拟
中图分类号： O241.82 文献标志码： A
Full hexahedral meshing method for 3D model based on
projection method
REN Huilong， GONG Shenglong， CAI Yongchang
（Department of Geotechnical Engineering， Tongji University， Shanghai 200092，China）
Abstract： As to the issue that the automatic hexahedral meshing is far more difficult than the automatic tetrahedral meshing， the projection method is used to transform the initial mesh of simple shapes into the actual mesh of each block， the mesh transition among blocks with different mesh density is implemented by node merging algorithm and resegmentation technique. The numerical examples show that the mesh quality generated by the method is good and the method is easy to bring into effect with wide applicability.
Key words： full hexahedral mesh； projection method； transfinite interpolation； resegmentation technique； mesh transition； explosion simulation
收稿日期： 2013-03-29 修回日期： 2013-05-24
基金项目： 西部交通建设科技项目（2011ZB04）
作者简介： 任辉龙（1989—），男，江西宜春人，硕士研究生，研究方向为岩土工程稳定性分析，（E-mail）hl.ren@hotmail.com
0 引 言
目前，在三维有限元模型网格生成中普遍采用四面体单元或六面体单元.与四面体单元相比，六面体单元具有更好的力学性能，达到同样精度所需的六面体单元和节点数远少于四面体单元，故六面体单元在有限元分析领域有较高的价值，是一种十分重要的单元.[1]
很多专家认为，在诸如锻压成形和爆炸冲击等三维有限元仿真领域中，六面体应该是首选单元.在某些情况下，如当用有限体积法和边界适应坐标系统求解复杂形体的控制守恒方程时，只能采用六面体单元进行有限元分析.[2-3]从技术角度讲，六面体网格自动划分的难度远远高于四面体网格.目前，有代表性的全六面体网格生成方法有映射法、扫描法、基于栅格法、扩展的AFT方法、多子区域法和投影法等.
映射法的优点是速度快、网格质量高、密度可控制，但同时也有2个严重的问题：复杂三维域自动分解问题和网格密度过渡问题.近年来，一些研究者采用中面法[4]进行子域自动分解，但对于复杂的几何形体难以实现全自动分区.扫描法是通过对离散化的基本单元形体进行旋转、扫描和拉伸等操作获得高维网格的一种方法.该方法难度较低，容易实现，但只适合于形状简单的三维物体.单元转换法将其他单元转化为六面体单元，可以实现全自动化[5-7]，但这种方法的边界拟合能力弱，生成的网格质量较差.左旭等[8] 改进单元转换法，将十节点曲边四面体转换为六面体，并采用非线性约束优化算法大幅提高六面体网格的单元质量.基于栅格法能实现网格生成的自动化，生成速度也非常快，但其最大弱点是边界单元的质量较差，网格密度很难控制.AFT方法的扩展主要有编须算法[9]和粘贴算法[10].编须算法基于空间缠绕连续集[11]概念，将AFT方法扩展到三维空间，生成的六面体网格质量是所有算法中最好的，但程序实现非常困难；粘贴算法可能会留下一些孔洞，这些未被划分的区域只能用四面体填充.
LI等[12]提出一种带梯度的全六面体网格划分方法，利用单元分割模板得到不同块体之间节点的约束关系，建立整数规划数学模型，解得不同块体之间的单元划分数量关系.该方法可以实现不同网格密度的合理过渡，但求取约束关系过于复杂.SU等[13]利用单元分割模板提出一种基于复合投影法建立多块体的全六面体网格划分方法，但该方法涉及到三棱柱单元的再划分.
本文基于投影法，利用单元分割技术实现不同块体的单元节点协调过渡，建立疏密有致的全六面体网格模型.
1 投影法
投影法的理论基础是射影几何，通过投影法将点投影到曲面上，将普通形状变换成具有相同拓扑性质的曲面形状.[14-15]复杂的三维模型都能分解为简单块体，而简单块体又可以通过简单形状经过布尔运算和投影得到.投影法生成六面体网格正是基于这样一种思想，将复杂几何模型分解为简单块体，提取所有简单块体的边界几何信息；生成具有初始六面体网格的简单实体，删除某些区域，使其与简单块体具有相同的拓扑结构，将具有初始网格的简单实体表面节点投影到简单块体的边界表面，内部网格节点通过超限插值分布在简单块体的内部，形成具有六面体网格的简单块体，最后将简单块体组装起来形成全六面体的复杂几何模型. 基于投影法划分网格是一种自顶向下的建模方法，与普通商业软件先建立实体模型再进行网格划分的流程相反，而正是这种自顶向下的建模方法可以建立任意复杂几何形状的全六面体网格模型.
基于投影法的六面体网格生成算法可以分解成以下步骤：（1）形成具有初始网格的初始结构；（2）将初始网格的自由节点投影到模型表面；（3）用最小变形角法匹配网格边界；（4）用最小扭转角法匹配模型顶点；（5）用超限插值技术获得网格内部节点坐标；（6）用再分割技术和节点合并算法形成整体网格模型.
1.1 超限插值法
超限插值法是最重要的映射方法，它将一个单位立方体映射到任意曲棱线的六面体区域，见图1.
图 1 基于曲棱线的六面体超限插值
Fig.1 Hexahedral transfinite interpolation based on ridge
由文献[16]有
P（ξ，η，γ）=（1-η）（1-γ）f1（ξ）+（1-η）γf2（ξ）+ηγf3（ξ）+η（1-γ）f4（ξ）+
（1-ξ）（1-γ）f5（η）+（1-ξ）γf6（η）+ξγf7（η）+ξ（1-γ）f8（η）+
（1-ξ）（1-η）f9（γ）+（1-ξ）ηf10（γ）+ξηf11（γ）+
ξ（1-η）f12（γ）+C（ξ，η，γ） （1）
C（ξ，η，γ）=-2（（1-ξ）（1-η）（1-γ）v（0，0，0）+（1-ξ）（1-η）γv（0，0，1）+
（1-ξ）η（1-γ）v（0，1，0）+（1-ξ）ηγv（0，1，1）+ξ（1-η）（1-γ）v（1，0，0）+
ξ（1-η）γv（1，0，1）+ξη（1-γ）v（1，1，0）+ξηγv（1，1，1）） （2）
式中：P（ξ，η，γ）={x（ξ，η，γ），y（ξ，η，γ），z（ξ，η，γ）}为六面体区域中任意一插值点；fi（s）（i=1，2，…，12，s=ξ，η，γ）为边界曲线函数； v（i，j，k）={x（i，j，k），y（i，j，k），z（i，j，k）}，（i，j，k=0，1）为8个顶点；ξ∈[0，1]；η∈[0，1]；γ∈[0，1].
六面体区域内插值点的坐标（x，y，z）与12条边界曲线、8个顶点及在立方体中相应坐标（ξ，η，γ）有关，而在生成六面体网格时只需要区域内有限的插值点，因此，对于边界曲线只需要有限节点的坐标值.设线1，5和9上节点分别为m+1，n+1和p+1，取ξ=（i-1）/m， i=1，2，…，m+1
η=（j-1）/n， j=1，2，…，n+1
γ=（k-1）/p， k=1，2，…，p+1（3）则可求得曲棱线六面体内任意节点的插值坐标.
1.2 形成整体协调模型
在简单块体的组装过程中，会面临不同块体对应面上节点的协调性问题.当单元节点不协调时，可以通过添加节点耦合方程获得整体网格模型.本文通过共用节点法形成全六面体网格模型.
在组装简单块体时，建立交界面的局部坐标，见图2，ζ和η方向的单元个数分别为n和m.需要合并节点的2个交界面分别称为主面和从面，不妨称网格数目较多的交界面为主面，单元数目为nm×mm；网格数目较少的交界面为从面，单元数目为ns×ms.如果块体交界面上节点满足一定的关系，就可通过节点合并算法和再分割模板形成整体模型.
图 2 交界面上节点局部坐标
Fig.2 Local coordinates of interface node
1.2.1 节点合并算法
节点合并算法指若主面和从面在ζ和η方向的单元数目对应相等，则可以建立主面和从面上节点的一一对应，实现对应节点的合并.
1.2.2 再分割技术
再分割技术 [14]指通过按一定的规则再分割已有的六面体单元，在适当位置添加相应节点，再经过适当组合，形成多个六面体单元的方法.通过再分割技术，可以使主面和从面上ζ和η方向单元节点数目相等，再利用节点合并算法就可实现简单块体的组装，形成协调的全六面体网格有限元模型.本文采用几种较简单、方便实施的再分割模板，单元再分割模板见图3.
图 3 单元再分割模板
Fig.3 Element resegmentation templates
1.2.3 不同网格密度过渡
交界面上ζ和η方向的节点数目相同是实现网格协调的充要条件，然而当交界面上主面与从面之间的网格数量满足一定的比例条件时，利用如图3所示的再分割模板，使得交界面上主、从面的节点数量相等.这些比例条件有：（1）当nm∶ns=1∶1，mm∶ms=1∶1时，通过节点合并技术就可以实现单元协调过渡.（2）当nm∶ns=4∶2，mm∶ms=1∶1或nm∶ns=1∶1，mm∶ms=4∶2时，用节点合并和再分割法过渡，见图4（a）.（3）当nm∶ns=4∶2，mm∶ms=4∶2时，用图4（b）方法过渡.（4）当nm∶ns=3∶1，mm∶ms=1∶1或nm∶ns=1∶1，mm∶ms=3∶1时，用图4（c）方法过渡.（5）当nm∶ns=3∶1，mm∶ms=3∶1时，用图4（d）方法过渡.（6）当nm∶ns=3∶1，mm∶ms=4∶2时，用图4（e）方法过渡.
因此，对于复杂三维模型，首先将三维模型分解为简单块体，简单块体的实质就是拓扑形状较简单，简单块体通过之间的共用面形成整体三维模型，利用共用面就可以快速匹配交界面.确定初始网格的数目，使得交界面上的网格数目满足比例条件（1）～（6）中的一条，就可以实现块体之间不同网格密度的过渡，满足节点协调性要求.已有的有限元软件主要通过逐步扩大过渡实现节点协调性，但这只是比例条件（1）的特殊情况.当简单块体之间有2个或多个相邻的交界面，或从面有其他单元时，则采用图4中5种过渡方法之一，就会在过渡面块体之外产生节点不协调的单元，见图5.（a）过渡方法1 （b）过渡方法2
（c）过渡方法3
（d）过渡方法4
（e）过渡方法5
图 4 交界面上不同密度网格过渡
Fig.4 Transition of mesh with different density at interface
图 5 单元节点不协调网格
Fig.5 Unconformable mesh of element nodes
不协调单元产生的原因是从面内的单元进行再分割，而从面周围的单元没有进行分割，产生一层节点不协调的单元.分析图4中5种过渡方法，从面周围若有其他单元，再分割的单元与周围的单元是不协调的.因此，对于节点不协调的单元需进行再分割.对于图5中节点不协调的单元，可采用图3（f）和3（g）模板进行再分割，实现单元节点的协调.
1.2.4 六面体模型修复
假设简单块体已通过投影法和过渡方式1～5组合生成八节点六面体复杂模型，模型由n个节点和m个单元组成.节点为（Ni，xi，yi，zi），i=1，2，…n.单元为（Eid，Aj，Bj，Cj，Dj，Ej，Fj，Gj，Hj），j=1，2，…，m；Aj，Bj，…，Hj为节点编号，顺序见图6.
图 6 八节点实体单元
Fig.6 8-node solid elements
单元的8个节点组成6个面，依次为（A，D，C，B），（A，E，H，D），（A，B，F，E），（B，C，G，F），（C，D，H，G）和（E，F，G，H），右手手指沿面上节点顺序，大拇指方向为面的法向方向，如图6中面EFGH的法向方向为n，n=（EF×EH）/|EF×EH|（4）式中：EF=（xF-xE，yF-yE，zF-zE）；EH=（xH-xE，yH-yE，zH-zE）.
提取所有六面体单元的6个面，总共有6m个面，Si=（Ni1，Ni2，Ni3，Ni4），i=1，2，…，6m.单元之间的交界面上节点不协调，则该界面不协调，见图7，设交界面S1=（1，4，3，2），S2=（1，2，6，5）：（1）具有公共节点（1，2）；（2）面法向量n1与n2反向；（3）面积|S1|<|S2|；（4）面积较小的S1面重心落在面积较大的S2面内.若4个条件同时满足，则S1和S2是节点不协调的交界面，由S2定位需要重新分割的单元.通过图7所示的交界面的判定方法，可以找出4种节点不协调单元的交界面，见图8.
图 7 不协调交界面
Fig.7 Unconformable interface
图 8 节点不协调单元交界面
Fig.8 Unconformable element interface of nodes
图8中，情况1和2是以2∶4为过渡方式的节点不协调单元的交界面，情况3和4是以1∶3为过渡方式的节点不协调单元的交界面.基于对不协调交界面的分析，提出消除上述节点不协调的算法.
步骤1 读入全六面体网格模型的节点数据和单元数据.每个节点包括有节点号和x，y，z方向的坐标；每个单元由8个节点、6个面组成，其中每个面由4个节点构成，形成每个节点的面数据.
步骤2 找出需要分割的六面体单元：比较所有单元的面，找出属于图8（a）～8（d）的不协调面，找出需要分割的六面体单元.
步骤3 对步骤2中找出的六面体单元，按模板6或7的方式进行再分割，形成新的六面体单元，并删除原单元.
步骤4 通过节点合并算法，合并图8（a）～8（d）中交界面上的节点，更新节点数据和单元数据.
步骤5 按一定格式输出全六面体网格模型.
按该步骤编制MATLAB程序，对存在部分节点不协调的六面体网格模型进行修复，形成协调的全六面体网格复杂模型.
2 全六面体网格建模实例
若几何模型较规则，按照该方法能确保单元质量；若几何模型过于复杂，按照该方法就不可避免地会产生少量畸变单元.虽然这些不良单元数量不多，占总单元数比例也不是很大，但会影响计算精度，甚至无法进行计算.因此，检测网格单元是否已经畸变必不可少，可依次检查最大和最小顶角以及纵横比确定.
（1）最大顶角，即单元表面相邻两边夹角的最大值，其警戒值为155°.
（2）最小顶角，即单元表面相邻两边夹角的最小值，其警戒值为18°.
（3）纵横比.过单元表面投影图形的四边中点做一个矩形，其长边与短边的比值就是纵横比.该比值的警戒值为5.
实例1 具有初始网格的简单实体向边界曲面投影，得到简单块体的六面体网格，利用单元的再分割和简单块体的组装形成整体六面体网格模型，实例1的六面体网格见图9.通过单元质量检查，最大顶角为151°，最小顶角为26°，纵横比为1～2.3，均小于警戒值.
图 9 实例1的六面体网格
Fig.9 Hexahedral mesh of example 1
实例2 浅埋连拱隧道爆破数值模型的建模.爆破模拟可以利用高性能炸药单元或等效炮轰载荷的方法.第一种方法的理论更可靠，该方法能很好地模拟爆源近区的冲击响应问题，也可以研究中远区的爆破震动传播问题；但第一种方法存在数值模型建立困难的问题.这是因为该方法要求炸药材料单元和炮孔周围网格划分得十分细小，而地下工程爆破震动数值模拟的计算模型又很大，如果没有合理的网格过渡，就会使模型单元数量太多，超过普通计算机的处理能力.
合理的数值建模是爆破模拟中的一个难题.由于精细模型建立的复杂性，数值模型只能降低单元质量，即采用过于刚硬的四面体单元或采用扫掠拉伸方法建立的六面体单元，但采用扫掠拉伸方法建立的部分六面体单元过于细长，纵横比超过20；或者只能缩小模型规模，简化为平面应变模型，无法真正建立疏密有致、又很好控制单元质量和单元总数的模型. 利用本文提出的方法，选取计算区域：140 m（x方向）×80 m（z方向）×72 m（y方向），隧道掌子面爆破施工，前、后洞掌子面部分区域采用边长为20 cm的六面体网格，经过多面过渡，依次过渡边长为60， 100， 200和400 cm的单元.该隧道爆破数值模型见图10，所有六面体单元的长宽比均不超过4，严格控制单元质量，模型单元总数193 256个，节点总数202 817个.通过单元质量检查，模型的最小顶角为16°，低于18°的单元少于总数0.5%，最大顶角为154°，纵横比为1.0～6.3，超过警戒值5的单元总数少于0.3%，单元质量总体合格.如果按20 cm边长的六面体网格进行精细划分，网格数量会超过1 000万个，而通过本文的方法将单元数量控制在20万个以内，普通的计算机也能计算.在文献[17]中，利用该模型模拟分析浅埋连拱隧道爆破过程，见图11～12，监测结果与模拟结果一致，振速合矢量波形很接近，特征节点最大振速为15.951 cm/s， 实测为15.726 cm/s，说明六面体网格模型有效，网格质量很好.
（a）开挖区网格划分
（b）隧道整体网格模型
图 10 浅埋连拱隧道爆破数值模型
Fig.10 Numerical explosive model of shallow buried multi-arch tunnel
图 11 测点处振速合矢量
Fig.11 Resultant vibration velocity at measurement point
图 12 特征节点合速度
Fig.12 Resultant velocity of feature node
3 结束语
将复杂三维几何模型分为简单块体，基于投影法，对简单形状初始网格进行投影变换，得到每个块体的实际网格，利用节点合并算法和再分割技术，实现不同块体之间的不同网格密度过渡，生成整体全六面体单元模型.算例表明，本文所提出的方法简单易行，生成的网格质量很好，可以有效解决爆炸与冲击数值模拟中模型建立问题，适合于各种要求六面体单元数值计算.在以后的工作中，将考虑如何使用更多形式分割模板及其组合形式实现不同密度网格过渡，以及寻找三维几何模型简单块体的自动划分方法.
参考文献：
[1] 金晶， 吴新跃. 有限元网格划分相关问题分析研究[J]. 计算机辅助工程， 2005， 14（2）： 75-78.
JIN Jing， WU Xinyue. Analysis of finite element mesh problems[J]. Comput Aided Eng， 2005， 14（2）： 75-78.
[2] 关振群， 宋超， 顾元宪， 等. 有限元网格生成方法研究的新进展[J]. 计算机辅助设计与图形学学报， 2003， 15（1）： 1-14.
GUAN Zhenqun， SONG Chao， GU Yuanxian， et al. Recent advances of research on finite element mesh generation methods[J]. J Computer-aided Des & Compu Graphics， 2003， 15（1）： 1-14.
[3] LANDER S F， STEFFAN H. Method to generate complex computational meshes efficiently[J]. Int J Commun Numer Methods Eng， 1994， 10（5）： 373-384.
[4] PRICE M， ARMSTRONG C. Hexahedral mesh generation by medial surface subdivision： Part II： solids with flat and concave edges[J]. Int J Numer Methods Eng， 1997， 40（1）： 111-136.
[5] 关振群， 单菊林， 顾元宪. 基于转换模板的三维实体全六面体网格生成方法[J]. 计算力学学报， 2005， 22（1）： 32-37.
GUAN Zhenqun， SHAN Julin， GU Yuanxian. All hexahedral mesh generation method for 3D solid model based on extended transform templates[J]. Chin J Comput Mech， 2005， 22（1）： 32-37.
[6] 陈沸镔， 谢步瀛. 建筑结构模型的四边形网格生成算法[J]. 计算机辅助工程， 2010， 19（1）： 17-21.
CHEN Feibin， XIE Buying. Quadrilateral mesh generation algorithm of building structure model[J]. Comput Aided Eng， 2010， 19（1）： 17-21.
[7] 徐敏艳， 陈华书， 范秦寅， 等. 多块网格划分技术在CFD仿真中的应用[J]. 计算机辅助工程， 2012， 21（4）： 65-68.
XU Minyan， CHEN Huashu， FAN Qinyin， et al. Application of multi-block meshing technology in CFD simulation[J]. Comput Aided Eng， 2012， 21（4）： 65-68. [8] 左旭， 卫原平， 陈军， 等. 三维六面体有限元网格自动划分中的一种单元转换优化算法[J]. 计算力学学报， 1999， 16（3）： 343-348.
ZUO Xu， WEI Yuanping， CHEN Jun， et al. An element transformation and optimization method in 3D automatic hexahedral meshes generation[J]. Chin J Comput Mech， 1999， 16（3）： 343-348.
[9] TAUTGES T J， BLACKER T， MITCHELL S. The whisker-weaving algorithm： a connectivity based method for constructing all hexahedral finite element meshes[J]. Int J Numer Methods Eng， 1996， 39（19）： 3327-3349.
[10] BLACKER T D， MEYERS R J. Seams and wedges in plastering： a 3D hexahedral mesh generation algorithm[J]. Eng Computers， 1993， 9（2）： 83-93.
[11] MURDOCH P， BENZLEY S E. The spatial twist continuum[C]//Proc 4th Int Meshing Roundtable， Albuquerque， 1995： 243-251.
[12] LI H， CHENG G. New method for graded mesh generation of all hexahedral finite elements[J]. Computers & Structures， 2000， 76（6）： 729-740.
[13] SU Y， LEE H K， KUMAR A S. Automatic hexahedral mesh generation for multi-domain composite models using a hybrid projective mesh-based method[J]. Computer-Aided Des， 2004， 36（3）： 203-215.
[14] GORDON W J， THIEL L C. Transfinite mappings and their application to mesh generation[J]. Appl Math & Comput， 1982（10）： 171-233.
[15] MIDDLECOFF J F， THOMAS P D. Direct control of the mesh point distribution in meshes generated by elliptic equations[J]. AIAA J， 1980， 18（6）： 652-656.
[16] COOK W A. Body oriented （natural） co-ordinates for generating three-dimensional meshes[J]. Int J Numer Methods Eng， 1974， 8（1）： 27-43.
[17] 任辉龙， 段群苗， 蔡永昌. 浅埋连拱隧道爆破的数值模拟[J]. 爆破， 2012， 29（4）： 70-75.
REN Huilong，DUAN Qunmiao，CAI Yongchang. Numerical simulation on explosion of shallow buried multi-arch tunnel[J]. Blasting， 2012， 29（4）： 70-75.
（编辑 于杰）
关键词： 全六面体网格； 投影法； 超限插值； 再分割技术； 网格过渡； 爆炸模拟
中图分类号： O241.82 文献标志码： A
Full hexahedral meshing method for 3D model based on
projection method
REN Huilong， GONG Shenglong， CAI Yongchang
（Department of Geotechnical Engineering， Tongji University， Shanghai 200092，China）
Abstract： As to the issue that the automatic hexahedral meshing is far more difficult than the automatic tetrahedral meshing， the projection method is used to transform the initial mesh of simple shapes into the actual mesh of each block， the mesh transition among blocks with different mesh density is implemented by node merging algorithm and resegmentation technique. The numerical examples show that the mesh quality generated by the method is good and the method is easy to bring into effect with wide applicability.
Key words： full hexahedral mesh； projection method； transfinite interpolation； resegmentation technique； mesh transition； explosion simulation
收稿日期： 2013-03-29 修回日期： 2013-05-24
基金项目： 西部交通建设科技项目（2011ZB04）
作者简介： 任辉龙（1989—），男，江西宜春人，硕士研究生，研究方向为岩土工程稳定性分析，（E-mail）hl.ren@hotmail.com
0 引 言
目前，在三维有限元模型网格生成中普遍采用四面体单元或六面体单元.与四面体单元相比，六面体单元具有更好的力学性能，达到同样精度所需的六面体单元和节点数远少于四面体单元，故六面体单元在有限元分析领域有较高的价值，是一种十分重要的单元.[1]
很多专家认为，在诸如锻压成形和爆炸冲击等三维有限元仿真领域中，六面体应该是首选单元.在某些情况下，如当用有限体积法和边界适应坐标系统求解复杂形体的控制守恒方程时，只能采用六面体单元进行有限元分析.[2-3]从技术角度讲，六面体网格自动划分的难度远远高于四面体网格.目前，有代表性的全六面体网格生成方法有映射法、扫描法、基于栅格法、扩展的AFT方法、多子区域法和投影法等.
映射法的优点是速度快、网格质量高、密度可控制，但同时也有2个严重的问题：复杂三维域自动分解问题和网格密度过渡问题.近年来，一些研究者采用中面法[4]进行子域自动分解，但对于复杂的几何形体难以实现全自动分区.扫描法是通过对离散化的基本单元形体进行旋转、扫描和拉伸等操作获得高维网格的一种方法.该方法难度较低，容易实现，但只适合于形状简单的三维物体.单元转换法将其他单元转化为六面体单元，可以实现全自动化[5-7]，但这种方法的边界拟合能力弱，生成的网格质量较差.左旭等[8] 改进单元转换法，将十节点曲边四面体转换为六面体，并采用非线性约束优化算法大幅提高六面体网格的单元质量.基于栅格法能实现网格生成的自动化，生成速度也非常快，但其最大弱点是边界单元的质量较差，网格密度很难控制.AFT方法的扩展主要有编须算法[9]和粘贴算法[10].编须算法基于空间缠绕连续集[11]概念，将AFT方法扩展到三维空间，生成的六面体网格质量是所有算法中最好的，但程序实现非常困难；粘贴算法可能会留下一些孔洞，这些未被划分的区域只能用四面体填充.
LI等[12]提出一种带梯度的全六面体网格划分方法，利用单元分割模板得到不同块体之间节点的约束关系，建立整数规划数学模型，解得不同块体之间的单元划分数量关系.该方法可以实现不同网格密度的合理过渡，但求取约束关系过于复杂.SU等[13]利用单元分割模板提出一种基于复合投影法建立多块体的全六面体网格划分方法，但该方法涉及到三棱柱单元的再划分.
本文基于投影法，利用单元分割技术实现不同块体的单元节点协调过渡，建立疏密有致的全六面体网格模型.
1 投影法
投影法的理论基础是射影几何，通过投影法将点投影到曲面上，将普通形状变换成具有相同拓扑性质的曲面形状.[14-15]复杂的三维模型都能分解为简单块体，而简单块体又可以通过简单形状经过布尔运算和投影得到.投影法生成六面体网格正是基于这样一种思想，将复杂几何模型分解为简单块体，提取所有简单块体的边界几何信息；生成具有初始六面体网格的简单实体，删除某些区域，使其与简单块体具有相同的拓扑结构，将具有初始网格的简单实体表面节点投影到简单块体的边界表面，内部网格节点通过超限插值分布在简单块体的内部，形成具有六面体网格的简单块体，最后将简单块体组装起来形成全六面体的复杂几何模型. 基于投影法划分网格是一种自顶向下的建模方法，与普通商业软件先建立实体模型再进行网格划分的流程相反，而正是这种自顶向下的建模方法可以建立任意复杂几何形状的全六面体网格模型.
基于投影法的六面体网格生成算法可以分解成以下步骤：（1）形成具有初始网格的初始结构；（2）将初始网格的自由节点投影到模型表面；（3）用最小变形角法匹配网格边界；（4）用最小扭转角法匹配模型顶点；（5）用超限插值技术获得网格内部节点坐标；（6）用再分割技术和节点合并算法形成整体网格模型.
1.1 超限插值法
超限插值法是最重要的映射方法，它将一个单位立方体映射到任意曲棱线的六面体区域，见图1.
图 1 基于曲棱线的六面体超限插值
Fig.1 Hexahedral transfinite interpolation based on ridge
由文献[16]有
P（ξ，η，γ）=（1-η）（1-γ）f1（ξ）+（1-η）γf2（ξ）+ηγf3（ξ）+η（1-γ）f4（ξ）+
（1-ξ）（1-γ）f5（η）+（1-ξ）γf6（η）+ξγf7（η）+ξ（1-γ）f8（η）+
（1-ξ）（1-η）f9（γ）+（1-ξ）ηf10（γ）+ξηf11（γ）+
ξ（1-η）f12（γ）+C（ξ，η，γ） （1）
C（ξ，η，γ）=-2（（1-ξ）（1-η）（1-γ）v（0，0，0）+（1-ξ）（1-η）γv（0，0，1）+
（1-ξ）η（1-γ）v（0，1，0）+（1-ξ）ηγv（0，1，1）+ξ（1-η）（1-γ）v（1，0，0）+
ξ（1-η）γv（1，0，1）+ξη（1-γ）v（1，1，0）+ξηγv（1，1，1）） （2）
式中：P（ξ，η，γ）={x（ξ，η，γ），y（ξ，η，γ），z（ξ，η，γ）}为六面体区域中任意一插值点；fi（s）（i=1，2，…，12，s=ξ，η，γ）为边界曲线函数； v（i，j，k）={x（i，j，k），y（i，j，k），z（i，j，k）}，（i，j，k=0，1）为8个顶点；ξ∈[0，1]；η∈[0，1]；γ∈[0，1].
六面体区域内插值点的坐标（x，y，z）与12条边界曲线、8个顶点及在立方体中相应坐标（ξ，η，γ）有关，而在生成六面体网格时只需要区域内有限的插值点，因此，对于边界曲线只需要有限节点的坐标值.设线1，5和9上节点分别为m+1，n+1和p+1，取ξ=（i-1）/m， i=1，2，…，m+1
η=（j-1）/n， j=1，2，…，n+1
γ=（k-1）/p， k=1，2，…，p+1（3）则可求得曲棱线六面体内任意节点的插值坐标.
1.2 形成整体协调模型
在简单块体的组装过程中，会面临不同块体对应面上节点的协调性问题.当单元节点不协调时，可以通过添加节点耦合方程获得整体网格模型.本文通过共用节点法形成全六面体网格模型.
在组装简单块体时，建立交界面的局部坐标，见图2，ζ和η方向的单元个数分别为n和m.需要合并节点的2个交界面分别称为主面和从面，不妨称网格数目较多的交界面为主面，单元数目为nm×mm；网格数目较少的交界面为从面，单元数目为ns×ms.如果块体交界面上节点满足一定的关系，就可通过节点合并算法和再分割模板形成整体模型.
图 2 交界面上节点局部坐标
Fig.2 Local coordinates of interface node
1.2.1 节点合并算法
节点合并算法指若主面和从面在ζ和η方向的单元数目对应相等，则可以建立主面和从面上节点的一一对应，实现对应节点的合并.
1.2.2 再分割技术
再分割技术 [14]指通过按一定的规则再分割已有的六面体单元，在适当位置添加相应节点，再经过适当组合，形成多个六面体单元的方法.通过再分割技术，可以使主面和从面上ζ和η方向单元节点数目相等，再利用节点合并算法就可实现简单块体的组装，形成协调的全六面体网格有限元模型.本文采用几种较简单、方便实施的再分割模板，单元再分割模板见图3.
图 3 单元再分割模板
Fig.3 Element resegmentation templates
1.2.3 不同网格密度过渡
交界面上ζ和η方向的节点数目相同是实现网格协调的充要条件，然而当交界面上主面与从面之间的网格数量满足一定的比例条件时，利用如图3所示的再分割模板，使得交界面上主、从面的节点数量相等.这些比例条件有：（1）当nm∶ns=1∶1，mm∶ms=1∶1时，通过节点合并技术就可以实现单元协调过渡.（2）当nm∶ns=4∶2，mm∶ms=1∶1或nm∶ns=1∶1，mm∶ms=4∶2时，用节点合并和再分割法过渡，见图4（a）.（3）当nm∶ns=4∶2，mm∶ms=4∶2时，用图4（b）方法过渡.（4）当nm∶ns=3∶1，mm∶ms=1∶1或nm∶ns=1∶1，mm∶ms=3∶1时，用图4（c）方法过渡.（5）当nm∶ns=3∶1，mm∶ms=3∶1时，用图4（d）方法过渡.（6）当nm∶ns=3∶1，mm∶ms=4∶2时，用图4（e）方法过渡.
因此，对于复杂三维模型，首先将三维模型分解为简单块体，简单块体的实质就是拓扑形状较简单，简单块体通过之间的共用面形成整体三维模型，利用共用面就可以快速匹配交界面.确定初始网格的数目，使得交界面上的网格数目满足比例条件（1）～（6）中的一条，就可以实现块体之间不同网格密度的过渡，满足节点协调性要求.已有的有限元软件主要通过逐步扩大过渡实现节点协调性，但这只是比例条件（1）的特殊情况.当简单块体之间有2个或多个相邻的交界面，或从面有其他单元时，则采用图4中5种过渡方法之一，就会在过渡面块体之外产生节点不协调的单元，见图5.（a）过渡方法1 （b）过渡方法2
（c）过渡方法3
（d）过渡方法4
（e）过渡方法5
图 4 交界面上不同密度网格过渡
Fig.4 Transition of mesh with different density at interface
图 5 单元节点不协调网格
Fig.5 Unconformable mesh of element nodes
不协调单元产生的原因是从面内的单元进行再分割，而从面周围的单元没有进行分割，产生一层节点不协调的单元.分析图4中5种过渡方法，从面周围若有其他单元，再分割的单元与周围的单元是不协调的.因此，对于节点不协调的单元需进行再分割.对于图5中节点不协调的单元，可采用图3（f）和3（g）模板进行再分割，实现单元节点的协调.
1.2.4 六面体模型修复
假设简单块体已通过投影法和过渡方式1～5组合生成八节点六面体复杂模型，模型由n个节点和m个单元组成.节点为（Ni，xi，yi，zi），i=1，2，…n.单元为（Eid，Aj，Bj，Cj，Dj，Ej，Fj，Gj，Hj），j=1，2，…，m；Aj，Bj，…，Hj为节点编号，顺序见图6.
图 6 八节点实体单元
Fig.6 8-node solid elements
单元的8个节点组成6个面，依次为（A，D，C，B），（A，E，H，D），（A，B，F，E），（B，C，G，F），（C，D，H，G）和（E，F，G，H），右手手指沿面上节点顺序，大拇指方向为面的法向方向，如图6中面EFGH的法向方向为n，n=（EF×EH）/|EF×EH|（4）式中：EF=（xF-xE，yF-yE，zF-zE）；EH=（xH-xE，yH-yE，zH-zE）.
提取所有六面体单元的6个面，总共有6m个面，Si=（Ni1，Ni2，Ni3，Ni4），i=1，2，…，6m.单元之间的交界面上节点不协调，则该界面不协调，见图7，设交界面S1=（1，4，3，2），S2=（1，2，6，5）：（1）具有公共节点（1，2）；（2）面法向量n1与n2反向；（3）面积|S1|<|S2|；（4）面积较小的S1面重心落在面积较大的S2面内.若4个条件同时满足，则S1和S2是节点不协调的交界面，由S2定位需要重新分割的单元.通过图7所示的交界面的判定方法，可以找出4种节点不协调单元的交界面，见图8.
图 7 不协调交界面
Fig.7 Unconformable interface
图 8 节点不协调单元交界面
Fig.8 Unconformable element interface of nodes
图8中，情况1和2是以2∶4为过渡方式的节点不协调单元的交界面，情况3和4是以1∶3为过渡方式的节点不协调单元的交界面.基于对不协调交界面的分析，提出消除上述节点不协调的算法.
步骤1 读入全六面体网格模型的节点数据和单元数据.每个节点包括有节点号和x，y，z方向的坐标；每个单元由8个节点、6个面组成，其中每个面由4个节点构成，形成每个节点的面数据.
步骤2 找出需要分割的六面体单元：比较所有单元的面，找出属于图8（a）～8（d）的不协调面，找出需要分割的六面体单元.
步骤3 对步骤2中找出的六面体单元，按模板6或7的方式进行再分割，形成新的六面体单元，并删除原单元.
步骤4 通过节点合并算法，合并图8（a）～8（d）中交界面上的节点，更新节点数据和单元数据.
步骤5 按一定格式输出全六面体网格模型.
按该步骤编制MATLAB程序，对存在部分节点不协调的六面体网格模型进行修复，形成协调的全六面体网格复杂模型.
2 全六面体网格建模实例
若几何模型较规则，按照该方法能确保单元质量；若几何模型过于复杂，按照该方法就不可避免地会产生少量畸变单元.虽然这些不良单元数量不多，占总单元数比例也不是很大，但会影响计算精度，甚至无法进行计算.因此，检测网格单元是否已经畸变必不可少，可依次检查最大和最小顶角以及纵横比确定.
（1）最大顶角，即单元表面相邻两边夹角的最大值，其警戒值为155°.
（2）最小顶角，即单元表面相邻两边夹角的最小值，其警戒值为18°.
（3）纵横比.过单元表面投影图形的四边中点做一个矩形，其长边与短边的比值就是纵横比.该比值的警戒值为5.
实例1 具有初始网格的简单实体向边界曲面投影，得到简单块体的六面体网格，利用单元的再分割和简单块体的组装形成整体六面体网格模型，实例1的六面体网格见图9.通过单元质量检查，最大顶角为151°，最小顶角为26°，纵横比为1～2.3，均小于警戒值.
图 9 实例1的六面体网格
Fig.9 Hexahedral mesh of example 1
实例2 浅埋连拱隧道爆破数值模型的建模.爆破模拟可以利用高性能炸药单元或等效炮轰载荷的方法.第一种方法的理论更可靠，该方法能很好地模拟爆源近区的冲击响应问题，也可以研究中远区的爆破震动传播问题；但第一种方法存在数值模型建立困难的问题.这是因为该方法要求炸药材料单元和炮孔周围网格划分得十分细小，而地下工程爆破震动数值模拟的计算模型又很大，如果没有合理的网格过渡，就会使模型单元数量太多，超过普通计算机的处理能力.
合理的数值建模是爆破模拟中的一个难题.由于精细模型建立的复杂性，数值模型只能降低单元质量，即采用过于刚硬的四面体单元或采用扫掠拉伸方法建立的六面体单元，但采用扫掠拉伸方法建立的部分六面体单元过于细长，纵横比超过20；或者只能缩小模型规模，简化为平面应变模型，无法真正建立疏密有致、又很好控制单元质量和单元总数的模型. 利用本文提出的方法，选取计算区域：140 m（x方向）×80 m（z方向）×72 m（y方向），隧道掌子面爆破施工，前、后洞掌子面部分区域采用边长为20 cm的六面体网格，经过多面过渡，依次过渡边长为60， 100， 200和400 cm的单元.该隧道爆破数值模型见图10，所有六面体单元的长宽比均不超过4，严格控制单元质量，模型单元总数193 256个，节点总数202 817个.通过单元质量检查，模型的最小顶角为16°，低于18°的单元少于总数0.5%，最大顶角为154°，纵横比为1.0～6.3，超过警戒值5的单元总数少于0.3%，单元质量总体合格.如果按20 cm边长的六面体网格进行精细划分，网格数量会超过1 000万个，而通过本文的方法将单元数量控制在20万个以内，普通的计算机也能计算.在文献[17]中，利用该模型模拟分析浅埋连拱隧道爆破过程，见图11～12，监测结果与模拟结果一致，振速合矢量波形很接近，特征节点最大振速为15.951 cm/s， 实测为15.726 cm/s，说明六面体网格模型有效，网格质量很好.
（a）开挖区网格划分
（b）隧道整体网格模型
图 10 浅埋连拱隧道爆破数值模型
Fig.10 Numerical explosive model of shallow buried multi-arch tunnel
图 11 测点处振速合矢量
Fig.11 Resultant vibration velocity at measurement point
图 12 特征节点合速度
Fig.12 Resultant velocity of feature node
3 结束语
将复杂三维几何模型分为简单块体，基于投影法，对简单形状初始网格进行投影变换，得到每个块体的实际网格，利用节点合并算法和再分割技术，实现不同块体之间的不同网格密度过渡，生成整体全六面体单元模型.算例表明，本文所提出的方法简单易行，生成的网格质量很好，可以有效解决爆炸与冲击数值模拟中模型建立问题，适合于各种要求六面体单元数值计算.在以后的工作中，将考虑如何使用更多形式分割模板及其组合形式实现不同密度网格过渡，以及寻找三维几何模型简单块体的自动划分方法.
参考文献：
[1] 金晶， 吴新跃. 有限元网格划分相关问题分析研究[J]. 计算机辅助工程， 2005， 14（2）： 75-78.
JIN Jing， WU Xinyue. Analysis of finite element mesh problems[J]. Comput Aided Eng， 2005， 14（2）： 75-78.
[2] 关振群， 宋超， 顾元宪， 等. 有限元网格生成方法研究的新进展[J]. 计算机辅助设计与图形学学报， 2003， 15（1）： 1-14.
GUAN Zhenqun， SONG Chao， GU Yuanxian， et al. Recent advances of research on finite element mesh generation methods[J]. J Computer-aided Des & Compu Graphics， 2003， 15（1）： 1-14.
[3] LANDER S F， STEFFAN H. Method to generate complex computational meshes efficiently[J]. Int J Commun Numer Methods Eng， 1994， 10（5）： 373-384.
[4] PRICE M， ARMSTRONG C. Hexahedral mesh generation by medial surface subdivision： Part II： solids with flat and concave edges[J]. Int J Numer Methods Eng， 1997， 40（1）： 111-136.
[5] 关振群， 单菊林， 顾元宪. 基于转换模板的三维实体全六面体网格生成方法[J]. 计算力学学报， 2005， 22（1）： 32-37.
GUAN Zhenqun， SHAN Julin， GU Yuanxian. All hexahedral mesh generation method for 3D solid model based on extended transform templates[J]. Chin J Comput Mech， 2005， 22（1）： 32-37.
[6] 陈沸镔， 谢步瀛. 建筑结构模型的四边形网格生成算法[J]. 计算机辅助工程， 2010， 19（1）： 17-21.
CHEN Feibin， XIE Buying. Quadrilateral mesh generation algorithm of building structure model[J]. Comput Aided Eng， 2010， 19（1）： 17-21.
[7] 徐敏艳， 陈华书， 范秦寅， 等. 多块网格划分技术在CFD仿真中的应用[J]. 计算机辅助工程， 2012， 21（4）： 65-68.
XU Minyan， CHEN Huashu， FAN Qinyin， et al. Application of multi-block meshing technology in CFD simulation[J]. Comput Aided Eng， 2012， 21（4）： 65-68. [8] 左旭， 卫原平， 陈军， 等. 三维六面体有限元网格自动划分中的一种单元转换优化算法[J]. 计算力学学报， 1999， 16（3）： 343-348.
ZUO Xu， WEI Yuanping， CHEN Jun， et al. An element transformation and optimization method in 3D automatic hexahedral meshes generation[J]. Chin J Comput Mech， 1999， 16（3）： 343-348.
[9] TAUTGES T J， BLACKER T， MITCHELL S. The whisker-weaving algorithm： a connectivity based method for constructing all hexahedral finite element meshes[J]. Int J Numer Methods Eng， 1996， 39（19）： 3327-3349.
[10] BLACKER T D， MEYERS R J. Seams and wedges in plastering： a 3D hexahedral mesh generation algorithm[J]. Eng Computers， 1993， 9（2）： 83-93.
[11] MURDOCH P， BENZLEY S E. The spatial twist continuum[C]//Proc 4th Int Meshing Roundtable， Albuquerque， 1995： 243-251.
[12] LI H， CHENG G. New method for graded mesh generation of all hexahedral finite elements[J]. Computers & Structures， 2000， 76（6）： 729-740.
[13] SU Y， LEE H K， KUMAR A S. Automatic hexahedral mesh generation for multi-domain composite models using a hybrid projective mesh-based method[J]. Computer-Aided Des， 2004， 36（3）： 203-215.
[14] GORDON W J， THIEL L C. Transfinite mappings and their application to mesh generation[J]. Appl Math & Comput， 1982（10）： 171-233.
[15] MIDDLECOFF J F， THOMAS P D. Direct control of the mesh point distribution in meshes generated by elliptic equations[J]. AIAA J， 1980， 18（6）： 652-656.
[16] COOK W A. Body oriented （natural） co-ordinates for generating three-dimensional meshes[J]. Int J Numer Methods Eng， 1974， 8（1）： 27-43.
[17] 任辉龙， 段群苗， 蔡永昌. 浅埋连拱隧道爆破的数值模拟[J]. 爆破， 2012， 29（4）： 70-75.
REN Huilong，DUAN Qunmiao，CAI Yongchang. Numerical simulation on explosion of shallow buried multi-arch tunnel[J]. Blasting， 2012， 29（4）： 70-75.
（编辑 于杰）