网格生成

网格生成（通用8篇）

网格生成 篇1

摘要：目前计算机工作站的低成本和可用性使其成为高性能计算的有吸引力的解决方案。网络技术的显著进步使得能够来实现高性能的全球计算。据此, 介绍了在动画生成中使用网格计算组成的系统中的性能结果。利用多节点的计算机结构及使用Globus网格管理软件对网格计算的性能和损耗进行定量评估和分析。实验结果表明, 多节点的网格计算系统使得生成的时间显著减少, 同时在动画生成过程中的CPU负载大大减少。

关键词：网格计算,动画生成,监控

1 介绍

随着信息技术的发展, 大量的数据需要被处理。单台计算机不再能满足处理这么大量数据处理的需要。因此, 许多解决方案被提出, 其中一个就是网格并行计算。网格计算是一种管理计算机集合以实现共同任务的技术。与超级计算机相比, 网格技术是松散耦合的, 并且使用网格管理工具包, 很容易构建具有强计算能力的可扩展和安全的网格。网格技术的另一个优点是经济实惠。例如, 公司可以在自己的机器中实现其网格基础架构, 并在机器空闲时向其分配任务。

本文设计了一个用于动画生成的网格架构。第2节将介绍设计目标, 第3节将讨论网格计算结构设计细节, 然后第4节将描述硬件配置及讨论实验结果, 最后在第5节进行总结。

2 设计目标

基于具有一定数量的计算机节点的系统设计, 把pov格式文件存储在系统中的每台计算机里。每台计算机能够从pov文件生成高清晰度图像, 而生成的图像将对网格中的所有计算机可见。在生成图像之后, 任何计算机应该能够将图像逐帧转换为高清晰度动画。但在这种情况下, 在图像或者动画生成期间, CPU和内存资源将会不足。为了解决这个问题, 可通过在该系统中使用网格计算, 并且可以通过添加多个服务器节点以减少CPU和存储器的占用率, 以实现加速处理数据的目标。其中, 网格上的每个节点都应该安装Globus Toolkit工具包。

3 网格计算结构设计

3.1 可扩展性

Globus工具包为网格管理提供了一些组件。Globus资源分配和管理GRAM是使得用户能够定位、提交、监视和取消远程作业的组件之一。此外, 监控和发现服务MDS可提供有关节点的状态和可用性的信息, Grid管理器可以使用这些信息来选择使用资源。

当前网格资源管理有其弱点。公共采用的架构是2层层次资源管理。然而, 在大规模网格网络中使用该架构, 一些资源将变得不可访问。此外, 使用2层架构使得网络QoS较差。

在新的框架中, 我们提供分层资源管理器。资源管理器有三种类型: (1) 管理一个特定资源的个人资源管理器 (IRM) ; (2) 管理集群中的资源的集群资源管理器 (CRM) ; (3) 网格资源管理器 (GRM) , 用于管理整个网格网络中的资源。在4层网格网络中, IRM用在最低层级1中, 中间层采用CRM, 在顶层框架使用GRM进行资源管理。

3.2 安全性

网格的安全性通过以下手段维护: (1) 防火墙:部署防火墙以保护网格免受恶意攻击。防火墙通过检查数据包的源IP地址和目的IP地址来管理网络。恶意和可疑数据包将被防火墙从网关外部过滤。防火墙可以限制从Internet到网格的访问, 它还可以限制从网格计算系统到外部Internet的访问尝试。 (2) Globus工具包:Globus工具包适用于我们的企业动画生成网格系统的设计。它提供了一个称为“Globus安全基础设施” (GSI) 的安全标准和模块。Globus安全基础设施在网格计算环境中的计算机之间提供了秘密的、防篡改的、可委托的通信支持。我们使用非对称加密RSA用于GSI以实现安全和可认证的通信。 (3) 证书:网格上的每个用户和服务都有一个已识别的证书。证书包含主题名称、公钥、证明公钥属于主题的证书颁发机构 (CA) 和CA的数字签名的信息。 (4) 相互认证:GSI使用安全套接字层 (SSL) 进行相互认证。SSL使用1024和2048位密钥长度的RSA算法。当两个单元的网格相互通信时, 他们将首先验证对方的第三方CA。在双重验证成功之后, 然后建立连接。 (5) 密码通信:缺省情况下, GSI不保证双方之间的加密通信。如果请求机密通信, GSI可以提供用于加密和解密的共享密钥。对于我们网格的设计, GSI的这两个特性都用于保护节点和服务器之间的通信。

3.3 外部集成

在我们的实验网格系统中, 需要允许集成外部资源。为了解决这个问题, 我们部署了开放网格服务架构OGSA。OGSA通过分布式异构和动态网格环境提供服务和资源的集成, 无论是在外部资源共享或服务提供方面。我们的设计有一些要求: (1) 全局名称空间:为了容易地访问外部数据和资源, 网格系统应当能够在不考虑位置或复制的安全约束下透明地与其他节点交互。 (2) 元数据服务:我们必须认识到调用和跟踪外部资源是很重要的。我们需要访问和管理跨管理域的实体元数据的权限。 (3) 场地自治性:获取资源的机制需要符合地方控制和政策。 (4) 资源使用数据:这是在网格网络上集成和交换外部资源使用数据的机制和模式。

3.4 监控

数据收集和分布机制的尺度是非常重要的。一个监测机制应建立以监测网格系统的当前性能。当pov文件的大小, 或者动画的质量增加时, 网格系统监控机制应该能够检测相应每个节点的性能。因此, 服务器将能够注意到潜在的资源缺乏并确定节点的数量。通常, 以下特性性能监视信息是系统或程序产生的数据的最重要的部分。

4 硬件配置及实验结果

在我们的实验网格环境中, 我们部署五台计算机, 包括一个服务器和四个客户端节点组成网格网络。在服务器端, 我们部署了英特尔至强服务器。这样的处理器可以满足动画和POV射线软件实现的计算。4GB DDR3内存RAM可为网格计算设备提供高吞吐量。西部数据2TB硬盘为30分钟的动画存储提供足够的空间。AMD HD5630还提供足够的图处理。对于客户端, 主要工作是网格计算, 因此我们上面部署的设备满足要求。

通过对动画的生成演算得到的实验结果, 我们发现4节点下比2节点或1点所消耗的计算时间有所减少, 同时CPU占用率也有所减少, 可见多节点对动画演算是有作用的。

5 结论

并行计算正在成为当代计算机科学中流行和重要的概念。在网格计算的环境中, 如何有效地监视节点信息和应用性能是网格计算中的一个难点问题。本文介绍了一种基于Globus toolkit的网格计算系统组成及共同的动画演算, 并介绍了网格技术研究和设计。设计了基于多节点和多级分布式结构的运算模型。

未来我们的想法是实现自动计算节点选择算法, 并将其集成到我们的集群和网格计算平台中。每个站点中的计算节点被选择为实时提供的运算信息, 软件包还可以针对计算节点的不同速度处理器, 意味着哪个所选计算节点应当执行什么计算在异构集群和网格协定环境中。

参考文献

[1]I.A.Klimonov, V.D.Korneev, V.M.Sveshnikov.Parallelization technologies for solving three-dimensional boundary value problems on quasi-structured grids using the cpu+gpu hybrid computing environment[J].Numerical Methods and Programming, Advanced Computing, 2016, (17) :65-71.

[2]S.Iturriaga, S.Nesmachnow, F.Luna, E.Alba.A parallel local search in cpu/gpu for scheduling independent tasks on large heterogeneous computing systems[J].Journal of Supercomputing, 2015, 71 (2) :648-672.

[3]Z.Junbo, W.Jian-Syuan, P.Yi, L.Tianrui.A parallel matrixbased method for computing approximations in incomplete information systems[J].IEEE Transactions on Knowledge and Data Engineering, 2015, 27 (2) :326-339.

网格生成 篇2

复杂外形三维非结构粘性直角网格的生成

在保证非结构网格逻辑关系不变的条件下,采用光顺、投影、分层加密的方法,解决直角网格物面不贴体的`难点,建立三维非结构贴体直角网格生成系统,以适应粘性流场和湍流的计算,同时对自适应网格的流场计算进行研究.对航空航天领域几类复杂外形飞行器进行了三维非结构贴体直角网格的生成.

作 者：李盾 Li Dun 作者单位：航天科技集团公司第十一研究院,北京,100074刊 名：计算机与数字工程 ISTIC英文刊名：COMPUTER AND DIGITAL ENGINEERING年，卷(期)：35(3)分类号：V2 TP3关键词：直角网格 贴体 自适应 复杂外形

网格生成 篇3

非结构化网格对复杂计算域有较强适应性,被广泛地应用在有限元分析和计算流体力学中。在多种网格生成算法中,Delaunay三角化法具有“最小内角最大”和“空外接圆”性质, 能够进行任意多连通域网格的自动生成,因此是目前主流的非结构网格生成方法。然而单纯地对某一模型进行网格剖分并不能满足工程技术需求。因为完整的数值模拟是由几何造型、网格剖分、数值计算和数据后处理四部分组成,网格剖分是几何造型与数值计算的桥梁。首先它要继承几何层上的拓扑信息保证块与块之间衔接部分的信息融合。并且网格单元自身能够搜索其周围单元的拓扑信息以方便数值计算过程。因此需要建立一套完整的网格生成体系来满足实际工程应用。而基于对象的编程技术正是为大规模软件开发应运而生的。本文采用面向对象的设计思想,结合课题组已有的几何造型体系,开发了一套三维非结构化网格生成软件,为实际工程应用提供技术支持。

曲面网格是三维非结构网格生成的基础,本文采用映射法实现曲面Delaunay三角剖分。为控制曲面网格品质,在网格剖分过程中引入黎曼度量加以控制[7],并给出了黎曼空间下的外接圆准则判定标准。另外三维Delaunay三角化过程中需要保证模型边界的完整性。本文采用文献[4,5]中的方法,首先通过一致边界恢复进行辅助插点,之后移动或分解边界上的辅助点,实验表明该方法对复杂模型有很强的适应力。另外,因为浮点运算和边界恢复带来的影响,网格插点过程中并不能完全保证形成空腔的正确性。针对此问题引入一个空腔检查机制控制空腔的形成,保证程序的健壮性。

1 网格生成算法

a) 根据实现过程,把生成Delaunay网格的各种算法分为三类:分治算法,三角网生长法,逐点插入法。在此采用Bowyer[1]-Watson[2]的逐点插入法,该算法可以描述为:1)首先确定一个插入点Q,并在原始三角面(四面体)集合T中搜索外接圆(外接球)包含插入点的三角形集合S。2)确定区域S的外边界集合E,删除集合S中的三角面(四面体)并形成一个空腔。3)将插入点Q与边界集合E连接,形成新的Delaunay网格。4)重复步骤1)到3),直到没有插入点为止。

b) 以上是算法的逐点插入过程,完整的Delaunay三角化过程如下:

1) 输入模型边界点集合,二维情况下确定一个四边形能够完全覆盖边界点区域,并将该四边形剖分为两个初始三角形;三维情况下确定一个六面体,并将该六面体剖分为五个初始四面体。

2) 利用Bowyer-Watson算法将边界点逐一插入。

3) 保证边界的完整性,删除计算域外的网格单元。

4) 利用Weatherill[3]的重心插点法进行计算域自动插点。

5) 分别用拓扑法和Laplace法对网格进行优化、光顺。

经过前人的研究和完善,目前已经能够对二维平面生成完美的Delaunay三角网格。但曲面网格的生成和三维网格边界恢复问题则仍是网格生成的研究重点,下面将针对这两方面进行介绍。

1.1 曲面网格生成

利用Delaunay三角化法进行曲面网格生成时需要通过映射的方式,即在二维参数域中进行Delaunay三角剖分,之后将剖分结果变换到三维物理域。然而映射法在处理简单平滑曲面时可以保证映射后的网格品质,但对于曲率变化剧烈的复杂曲面则会使映射后的网格产生很大的扭曲,网格品质下降。本文将在平面参数域进行三角剖分时引入一个黎曼度量矩阵加以控制,并给出了黎曼空间下的外接圆准则判定标准,以此保证最终曲面网格的品质。

有空间曲面∑,其参数表达式为

${\rm P}(u,v)=\left(\begin{array}{l}x(u,v)\\ y(u,v)\\ z(u,v)\end{array}\right),(u,v)\in \Omega (1)$

式中,Ω是∑的参变量平面。设在任意位置存在一阶导数:$\stackrel{\rightharpoonup }{{\rm P}}{}_u=\frac{\partial \stackrel{\rightharpoonup }{{\rm P}}}{\partial u}‚\stackrel{\rightharpoonup }{{\rm P}}{}_v=\frac{\partial \stackrel{\rightharpoonup }{{\rm P}}}{\partial v}$,两个矢量可以定义出P点的切平面。设Ω上有一段直线AB,则对应于∑上的一条曲线记为ab,AB用参数表示为A+tAB,t∈[0,1]。则ab长度为

$l(a,b)={\displaystyle {\int }_0^1\sqrt{\stackrel{\rightharpoonup }{{\displaystyle AB}}{}^{{\rm T}}{\rm M}(A+tAB)\stackrel{\rightharpoonup }{{\displaystyle AB}}\text{d}t}}(2)$

式中,

$\stackrel{\rightharpoonup }{{\displaystyle AB}}=\left(\begin{array}{l}\Delta u\\ \Delta v\end{array}\right)$

;

${\rm M}(A+tAB)=\left(\begin{array}{l}a(t)b(t)\\ b(t)c(t)\end{array}\right)$

是点A+tAB处的黎曼度量矩阵,且是正定矩阵。其中,$a(t)=\stackrel{\rightharpoonup }{{\rm P}}{}_u(t)\cdot \stackrel{\rightharpoonup }{{\rm P}}{}_v(t)$;$b(t)=\stackrel{\rightharpoonup }{{\rm P}}{}_u(t)\cdot \stackrel{\rightharpoonup }{{\rm P}}{}_v(t)$;$c(t)=\stackrel{\rightharpoonup }{{\rm P}}{}_v(t)\cdot \stackrel{\rightharpoonup }{{\rm P}}{}_v(t)$。式(2)可以转换为

$l(a,b)={\displaystyle {\int }_0^1\sqrt{a(t)\cdot \text{Δ}u{}^2+2b(t)\cdot \text{Δ}u\cdot \text{Δ}v+c(t)\cdot \text{Δ}v{}^2\text{d}t}}$。 (3)

式(3)可以用梯形公式近似计算,并将其定义为Ω上两点的距离。这个度量反映了通常曲面∑上的曲线段长度,用该度量对在Ω进行三角剖分, 并转换到∑上, 则可以得到曲面上满意的非结构网格。

进行Delaunay三角化过程中外接圆准则是算法的关键之一,需要准确计算三角面的外接圆圆心和半径。通常二维平面下能够准确地计算出三角面的圆心,但在三维曲面下每一点的度量矩阵并不相同,很难准确计算三角面的圆心。因此本文采用一种近似的求三角面圆心的方法。设Z是参变量平面Ω上的一点,则在黎曼空间域(Ω,M(Z))内两点的距离记为lz。对于三角形P1P2P3采用以下公式计算圆心O:

并通过公式(5)判断插入点P是否满足外接圆准则。

lP(O,P)<lP(O,P1) (5)

1.2 边界恢复

因为通过Delaunay三角剖分生成初始网格后会出现边界丢失的现象,所以要进行边界恢复过程。二维情况下的边界恢复问题可以通过边交换的方式方便地解决,但在三维区域下,情况则要复杂的多。通常解决的方法有两种:一致边界恢复和约束边界恢复。一致边界恢复需要在边界上添加辅助点来完成边界恢复。但添加辅助点破坏了原有单元尺寸场的定义,且可能在其附近形成品质较差、甚至退化的网格单元。约束边界恢复是通过边、面交换操作和在一些局部区域添加辅助点来实现的。但在文献[8]中给出一些无法解决的例子,这些例子在边界上常有一些狭长的三角单元,或存在局部区域单元尺寸剧烈过度的现象。本文采用先在边界上添加辅助点完成一致边界恢复,再移动或分解边界上的辅助点,以达到剔除它们的目的。一致边界恢复方法在[3]中有详细的说明,这里将主要介绍如何将辅助点分解的过程。

图1给出了丢失边上辅助点分解的过程。图中边AC是丢失边,E是添加的辅助点,三角面ABC和ACD共享边AC,G是BC边上的辅助点,F是三角面ABC的辅助点。则分解E点的过程为:

1) 确定三角面ABC和ACD的法向量$\stackrel{\rightharpoonup }{n}{}_1‚\stackrel{\rightharpoonup }{n}{}_2$,令$\stackrel{\rightharpoonup }{n}=(\stackrel{\rightharpoonup }{n}{}_1+\stackrel{\rightharpoonup }{n}{}_2)/2$。沿$\stackrel{\rightharpoonup }{n}$扰动E,扰动距离为ε1,则$E{}_1=E+\stackrel{\rightharpoonup }{n}\times \epsilon {}_1$。并保证E1对E点所连接的四面体组成的球体可视。

2) 沿$\stackrel{\rightharpoonup }{n}$相反方向扰动E,扰动距离为ε2,则$E{}_2=E-\stackrel{\rightharpoonup }{n}\times \epsilon {}_2$。同样保证E2对E点所连接的四面体组成的球体可视。

3) 注意到三维多边形AFGCD可以是由多边形ACD、AFGC组成,分别三角化两组多边形则可以得到丢失边AC。利用分解后的三角形分别与E1,E2两点生成四面体。

由于AC上只有一个辅助点E,E被分解后,丢失边便直接恢复了。如果AC边上还有其他辅助点,则继续按上面的方式进行分解,最终恢复丢失边。

丢失面上的辅助点分解与丢失边的方法类似,仅是在求解扰动方向$\stackrel{\rightharpoonup }{n}$时,可以直接通过辅助点所在丢失面的法矢量确定。

1.3 空腔判断准则

由于浮点计算误差和边界恢复过程带来的影响,插点过程中形成的空腔并不能完全保证其正确性。因此必须保证插入点对空腔边界单元的可视,下面给出了保证空腔正确性的方法:

1) 设空腔集合为T,找出空腔的外边界三角面集合F。通过Fi和插入点P计算出混合积Si,再由Fi和其所在空腔的四面体的对顶点计算出混合积Ki。

2) 逐一计算Si×Ki,

如果Si×Ki>0,则继续比较;

如果Si×Ki<0,从集合T中移除四面体Ti;

如果Si×Ki=0,则将三角面Fi连接的另外一个四面体加入到集合T。

3) 如果集合T经过步骤2)没有改变,则结束空腔正确性判断;否则转到1)重新进行判断。

2 基于对象的网格生成软件体系

2.1 网格生成软件体系

目前网格生成技术多种多样,算法也趋于成熟,但在解决工程实际问题中仍面临一些问题。首先,网格生成的前提是几何模型,而如何从模型中获得合适的几何信息则是首要的问题。大多CAD软件的几何、拓扑信息由于其自身的特点,不能直接运用到网格生成中去。需要对于几何信息进行特殊的清理与修补[6],这是一件繁琐的工作。另外,CAD软件很难通过参数化构建几何模型进行实体外形设计。因此在课题组已有的几何造型体系下,搭建一个网格生成系统将极大地方便几何信息与网格信息之间的融合;并且能够通过参数化设计快速生成几何模型以及计算网格,方便下游的数值计算模拟,缩短设计周期。

网格生成体系的完整结构如图2所示,包括几何造型和网格生成两部分。几何造型体系具有生成标准图形及NURBS曲线、曲面的功能。几何实体由几何点、几何线、几何曲面组成,通过它们之间建立的拓扑关系,能够互相进行上下层几何信息搜索。网格生成部分可以进行二维平面、三维曲面和三维实体的非结构化网格生成。针对几何层上的点、线、面建立对应的网格信息存储模块,即网格点(Grid_Vertex)、线网格(Grid_Line)、面网格(Grid_Face)、体网格(Grid_Solid) 四部分。其中线网格中包含曲线上的网格点集合、网格边集合信息;面网格中包含曲面上的网格点集合、网格边集合以及三角面集合信息;体网格中记录网格点集合、四面体集合信息。通过网格信息存储模块与几何层上的点、曲线、曲面建立的关系,便实现了网格信息部分与几何信息部分的有机融合,继承了几何层上的拓扑关系,为网格生成体系搭建了几何拓扑基础。

2.2 基于对象的数据存储结构

通过面向对象技术将几何层中的点、线、面和网格信息存储模块中的网格点、线网格、面网格、体网格建立相应的类,并且在各自类中记录上下层对象的指针,以此将几何信息与网格信息联系起来。需要说明的是体网格的上层没有几何体的信息,因此它通过记录组成体的几何面指针数组来实现拓扑信息的融合。

网格生成过程中的基本网格单元包括网格点(Grid_Vertex)、网格边(Grid_Edge)、三角面(Grid_Triangle)和四面体(Grid_Tetrahedron),同样建立相应的类来表示。因为考虑到Delaunay三角化过程中需要进行频繁的单元拓扑信息搜索和网格的删除、生成工作,例如在确定插入点形成的空腔时就需要通过一个基四面体进行相邻单元的树型搜索来完成。因此建立完整的网格单元拓扑结构可以方便地进行单元之间的拓扑搜索,提高Delaunay三角化效率。网格单元之间的拓扑结构如图3所示。其中,网格点类通过指针链表记录点所连接的边;网格边类同时记录组成点的指针数组以及边所连接三角面的指针链表;三角面类记录组成的点、边的指针数组以及两侧的四面体单元指针;四面体类则同时记录组成的点、边、面数组指针。如此每种类型的单元均能方便地对上下层单元进行搜索。

3 算例

定义表征网格单元品质的扭曲率为

$SR=\max (\frac{\alpha {}_{\max }}{120°}-0.5,1.0-\frac{\alpha {}_{\min }}{60°})$(6)

二维情况下,α为网格单元中任意相邻网格边的夹角,因此当SR→0.0,网格品质越好。三维情况下α为任意相邻三角面的的夹角,正四面体的二面角约为70.5°,因此当SR→0.083时,网格品质越好。

例1:参数曲面表达式为z=sin(πx)sin(πy),其中-1≤x,y≤1。图4显示了参数域下的三角剖分结果,可以看出三角网格呈现各向异性的特征。图5是由参数域映射到曲面空间域的三角网格,网格为各向同性。平均网格品质为0.2438,网格品质在0.1以下的单元占90.12%,三角形夹角大多在50°～70°之间,网格品质较好。

例2:带狭缝的长方六面体非结构网格。表面三角网格总数为7318,丢失的三角网格数为24,占总网格的0.33%。三维四面体网格数为43530,平均品质为0.32,网格品质在0.2以下的占80.62%。网格恢复过程如图6所示。

例3:几何模型为上表面是Volcano参数曲面,下半部分为一立方体。表面三角网格总数为5234,初始网格生成后丢失面个数为48,占总网格的0.92%。三维四面体网格数为42365,平均网格品质为0.2921,网格品质在0.2以下的占78.23%。网格恢复过程如图7所示。

4 总结

通过以上算例表明该网格生成软件很好地解决了曲面网格和边界恢复等网格生成中的关键问题,对复杂三维实体有较强的适应性和通用性,网格品质较高。采用面向对象的程序设计思想成功地将课题组已有的几何造型体系和网格生成体系进行了有机的结合,为工程应用提供了方便。

参考文献

[1] Bowyer A. Computing Dirichlet tessellations[J]. The Computer Journal, 1981(24):162-166.

[2] Watson D F. Computing the n-dimensional Delaunay tessellation with application to Voronoi polytopes[J]. The Computer Journal, 1981(24):167-172.

[3] Weatherill N P, Hassan O. Efficient three-dimensional Delaunay triangulation with automatic point creation and imposed boundary constraints[J]. International Journal For Numerical Methods in Engineering, 1994(37):2005-2039.

[4] George P L, Borouchaki H, Saltel E. ‘Ultimate’ robustness in meshing an arbitrary polyhedron[J]. International Journal for Numerical Methods in Engineering, 2003(58):1061-1089.

[5]Qiang Du,Desheng Wang.Constrained boundary recovery for three dimensional Delaunay triangulations[J].International Jour-nal for Numerical Methods in Engineering,2004(61):1471-1500.

[6]Gopalsamy S,TzuYi Y.A Geometry engine for CAD/Grid inte-gration[J].41st Aerospace Sciences Meeting and Exhibit,2003.

[7]Borouchaki H,George P L,Hecht F,et al.Delaunay mesh gen-eration governed by metric specifications[J].Part I.Algorithms.Finite Element Analysis and Design,1997,25(1-2):61-83.

网格生成 篇4

一种改进的四边界插值网格生成方法及其应用

通过改进一种四边界插值方法,并引入提出的网格正交控制及加权平均光顺措施,给出了一种非常有效地代数网格生成方法.应用该方法生成的网格贴体性、疏密性及均匀性都较好,且网格生成非常迅速.应用该方法对各种复杂气动外形进行了较高质量网格生成,表明了此方法的`有效性、实用性.

作 者：杨旭东 乔志德 朱兵 作者单位：西北工业大学,翼型研究中心,陕西,西安,710072刊 名：西北工业大学学报 ISTIC EI PKU英文刊名：JOURNAL OF NORTHWESTERN POLYTECHNICAL UNIVERSITY年，卷(期)：20(1)分类号：V211.3关键词：网格生成 四边界插值 正交控制 加权平均

网格生成 篇5

关键词：纹理特征,网格细分,二叉树,三维地形

(一) 引言

地形三维可视化计算近年来一直是相关领域的热点研究问题, 在三维仿真和虚拟地形环境中占有十分重要的地位。关于虚拟环境中三维地形的生成大多数是采用简化算法, 在保留地形细部特征的同时来提高其显示效率, 在简化过程中, 需要对消除裂缝、T型连接和简化模型误差度量等进行计算。本文给出了一种无需考虑消除裂缝和连接计算的三维地形生成方法, 依据地形纹理特征曲线从低精度网格地形向高精度网格细化, 使得平坦的地形区域用少量的三角形网格绘制, 起伏较大的地形区域用较多的三角形网格显示, 也形成一种多分辨地形模型结构, 以提高显示质量和效率。

(二) 方法依据和思路

在三维地形场景中, 虽然地形表现为三维模型, 但表示地形的DEM数据可以等价于由高度场生成的二维灰度纹理图像, 亮度变化小的地方, 表明此区域较为平坦, 亮度变化大的地方, 此地形可能山脉、谷等区域, 而“灰度图像边缘”正可以用来描述纹理图像局部区域亮度变化显著的部分。若能在三维空间中细化“边缘特征曲线”周围的网格正是本方法的关键思路所在, 使得边缘特征曲线周围网格密集, 增强其几何特征。

方法实现基本步骤及思路:

(1) 采用图像处理技术, 对地形高程灰度图像进行边缘检测, 从中提取图像的边缘曲线, 作为下一步网格细分的依据。

(2) 将含有地形特征曲线的图像作为纹理, 映射到一个低精度平面网格上。

(3) 根据三角边与它的纹理特征曲线相交情况来细分网格, 细分后的地形网格存储在二叉树数据层次结构中。

(4) 根据高程图像灰度值为地形网格的每个顶点赋予高度值, 最后绘制二叉树结构中所有叶子数据就是地形最终显示结果。

(三) 方法实现过程

1. 网格模型的定义

三角网格除了是由三维空间中的三角片通过边和顶点连接而成的分片线性曲面外, 我们将纹理图像信息作为曲面描述数据定义在网格模型中, 便于细分网格时作为参考对象, 以及增强显示效果, 所以, 这里网格M可以定义为一个三元组, M= (K, V, U) 。

其中, U={u1, …, um}, uiЄR2表示像素点在二维纹理图像中的坐标。

V={v1, …, vm}, viЄR3表示网格M的顶点在三维空间中的位置。

K为表征网格拓扑结构的单纯复形, 一个单纯复形包含一组单形, 其中{i, j}称为边, 记为L;{i, j, k}称为三角形面, 记为T, 三角形每条边最多只能包含在两个三角片中, 并且一个顶点和一个纹理坐标相对应。

2. 网格三角面拆分过程

网格细分的主要目的是为了增加网格细节层次;在地形网格模型中, 这种细节层次主要表现在地形的起伏特征, 根据地形起伏特征自动调整三角网格绘制数量, 在起伏大的区域, 三角网格数量多, 平坦区域, 三角网格数量少;首先预生成一个规则初始地形网格M0, 将含有地形纹理特征的图像映射到网格M0中, 经过网格n次细分, 生成M=M0→M1→M2…→Mn序列, 设T是M0经第i次细分的网格模型Mi中一个三角形, 对T细分过程如图1所示。

因为每个三角形网格顶点都有一个纹理坐标 (ui, uj) 相对应, 首先沿着三角形的绘制方向判断边和特征曲线是否相交, 那么新的交点与该边对应的顶点连线, 将三角形拆分成两个;在图1中的虚线表示地形纹理图像特征曲线, 图1 (a) 中的v1、v2和v3是三角形与特征曲线交点, 首先边{1, k}将三角形T拆分成T1、T2两个三角形, 然后边{2, 1}将T2拆分成T3、T4, 边{3, 1}将T1拆分成T5、T6。将新的顶点v1、v2和v3映射回三维空间, 依据地形高度场赋予新的高度值。一个三角形可能被拆分成4个、3个、2个三角形或不被拆分4种情况, 这样完成一个三角形的细分, 其它三角形网格采用同样的方法, 最终完成整个地形网格模型的细分, 即Mi→Mi+1。

由上述三角面拆分过程, 可以得到图2所示的二叉树细分层次结构。图中的某一层Mx的叶子或没有子树的根是该层网格模型最终绘制结果, 层次越深, 细化程度越高, 树型分支和叶子数量增多, 导致绘制的三角形网格数量的增多。

网格地形的简化过程, 是细化过程的逆过程, 而且无需判断三角形边的拆分情况, 只需将某一层Mx的三角形网格的子树和叶子删除, 就可以达到简化的目的。

3. 细分规则

由于存在提取的边缘线可能与三角形边有多个交点且细分后的三角形法线方向确定等问题, 约束了细分网格效果, 所以, 为了得到较为理想的生成结果, 需要对网格按照一定的规则细分和绘制。

(1) 为了使细分前和后三角面法线方向一致, 避免绘制三角面空洞现象, 将细分后的三角面绘制方向与上一级的三角面绘制方向应该保持一致, 如:三角面T的绘制方向为vivjvk, 那么拆分后的三角面T1、T2的一种绘制方向分别为viv1vk和v1vjvk。

(2) 在三角面细分过程中, 当新增加的拆分顶点v′与该边的端点很接近时, 造成细分后的某一三角形与细分前的三角形很接近, 分裂前后几乎没有发生变化, 即图3中的T2≌T, 所以应给出ωЄ (0, 0.5) , 当|v'-vi|<ω|vj-vi|或|v'-vj|<ω|vj-vi|时, 则新增加的顶点v′为无效点, 此时三角面T的边{i, j}不被拆分。

(3) 三角面的边可能与边缘线存在多个相交点, 为了保证三角面尽可能地均匀地被拆分, 取三角形边中点最近的交点作为该边的拆分点, 即:拆分点P要满足up={|up-0.5× (ui+uj) |}min。

(4) 按照三角形绘制方向来拆分三角形, 如图1中三角形T的绘制方向为vivjvk, 那它的拆分顺序是边vivj, vjvk, vkvi。

(5) 拆分顶点应与最近被拆分的三角形顶点连线来拆分该三角形, 如图1中 (b) 的拆分顶点v2与最近被拆分的三角形T2顶点v1连线来拆分T2, 而不是与三角形T的顶点vi的连线。

4. 数据结构描述和实现过程

根据上述细分过程, 地形网格细分数据可以采用结构规范、搜索较快的二叉树存储, 将每个网格三角形拆分为两个三角形, 拆分前的三角形作为二叉树的根节点, 拆分后的两个三角形作为根节点的子节点, 若子节点不能够再拆分, 则这两个三角形称之为叶子;这些三角形逻辑上就像一组相连的邻居 (左右邻居) , 在细分完地形后, 就建立了三角网格二叉树型结构, 树的叶节点保存了图形绘制缓冲区中的三角形, 针对二叉树描述地形网格特点, 本文如下描述网格二叉树结构:

随着细分层次的增加, 三角形顶点的数量也不断增加且不固定, 因此顶点可以采用链表结构进行存储, 而变量m_vertexIndices就存储了链表队列中的位置;levelNUM代表了三角形拆分过程中所在的层次, 相当于图2中的M层次位置, 便于细化和简化的判定;canDivide是为了提高细分速度, 用来存储该三角形是否能再次细分, 若不能再细分, 就不用再判断二维纹理中三角形各边与特征曲线相交状态。

首先提取由地形的高度场得到灰度图的特征曲线图形, 将基本网格坐标映射到特征曲线的二维平面中, 那么在二维平面中也构成了网格, 判断网格与特征曲线相交情况, 产生新的顶点和三角形网格, 特征曲线密集的地方, 网格细分程度越高, 从而得到的三维地形中, 在地形平坦的地方以最少的存储空间, 起伏变化大的地方用较多的存储空间;如图3所示。

(四) 方法测试

测试目的:采用基于图像纹理特征网格细分的三维地形生成方法对显示效率、质量的影响和关系。

测试手段:

使用C++语言, 以OpenGL作为底层的图形接口, 在C++Builder6.0编译环境中进行测试。建立了CLand (地形生成类) 、CPicture (图像处理类) 、BST (二叉树类) 、CMath (数学计算类) 、CView (窗口显示类) 等主要类。在PⅣ1.8GHz, 256MB内存, Geforce3 64MB显卡机器上, 在800×600窗口, 水平视角450, 误差限2个像素情况下, 定义初始的规则网格数据M0为30×30×2=1800个三角片;原始网格定义为511×511×2=522242个三角片 (数目是根据高程灰度纹理分辨率定义) , 并赋予高度值;细分程度定义为细分后的三角形数目与原始网格三角形数目的百分比;误差定义为原始网格和细分网格的双边Hausdorff距离与原始网格的包围盒对角线长度的百分比, 其测试数据如表1所示, 绘制效果图4所示。

测试结果:在显示效率不低于30帧/秒时, 其中M3能清楚地描述地形几何特征, 误差较小, 显示效率可以达到32帧/秒, 且地形特征变化明显 (图4的 (a) , (b) 比较) , 从M0~M3误差变化较大, 表明能够迅速收敛于地形细节特征;当显示效率低于30帧/秒时, 地形特征变化并不明显 (图5的 (b) , (c) 比较) , M4~M6误差变化小, 而且显示效率并不理想。所以, 在方法应用中, 通常需要给出x={fFS (Mx) ≥FS}max测试程序, 即在满足显示效率的同时, 最大限度地提高其显示质量, 目的是兼顾3D模型显示效率和显示质量, 其中fFS () 用来测试网格Mx显示效率, FS是预先给定的显示效率常数, 如本文以30帧/秒显示效率作为分界线, x是多层次模型最大细分网格深度, 如本文x=3。在实际程序运行过程中, x会根据给定的FS值, 通过网格细化和简化自动调整其网格绘制深度。

(五) 结论

本文给出了一种基于纹理图像特征网格细分的三维地形生成方法, 方法思路简单且容易实现, 通过提取地形高程特征值作为网格细分依据, 形成多分辨率地形模型结构, 使得显示效率和保留地形细节特征尽可能兼顾。从第3节可知, 由于采用三角边拆分方式, 共享该边的一个三角面被拆分, 另外三角面一定被拆分, 因此, 不存在消除裂缝和连接计算, 在一定程度上减少程序设计的复杂过程。

另外, 该方法还可以结合“视点相关”三维地形显示技术, 在显示效率方面还会进一步提高, 这是今后的主要工作。

参考文献

[1]李偈, 张松海, 刘强.一种基于四叉树结构的动态多分辨率地形模型[J].计算机工程与应用, 2006, 42 (7) :202-204.

[2]杜剑侠, 李凤霞, 战守义.LOD算法研究及其在地形实时显示中的应用.计算机工程与应用, 2005, 41 (13) :211-231.

[3]韩敏, 汤松涛, 李洋.大规模地形实时可视化算法[J].计算机工程, 2008, 34 (13) :270-272.

[4]张立强, 杨崇俊.多进制小波和二叉树实现大规模地形的实时漫游[J].计算机辅助设计与图形学学报, 2005, 17 (1) :467-472.

网格生成 篇6

混合网格由黏性区三棱柱网格及非黏性区四面体网格组成, 很多研究者对此类混合网格进行了相关研究。Pirzadeh[8]起初把层推进法应用于二维情况, 然后把该方法推广到三维情况, 并透过复杂几何体展示了算法的性能。L9hner等[9]采用层推进法对复杂几何体生成边界层网格, 并求解了N-S方程。Sang[10]等采用层推进法在壁面附近生成结构网格, 并同时应用于二维与三维情况, 并通过对机翼等几何模型进行数值计算验证了算法的性能。

在采用层推进法过程中, 计算节点推进矢量是非常关键的一个步骤, 在上述方法中, 节点推进矢量的确定是通过计算与该节点相连的所有面元法向矢量的平均值来确定的, 该方法对于相对平滑的流形能够产生合理的表面法向。然而对于离散复杂拓扑结构这可能违背流形的可视化条件即推进后新节点位置从所有流行面上看是可视的, 比如在离心泵的蜗壳隔舌处以及叶轮进出口处, 采用该方法计算生成的推进矢量可能会与当前的流行面相交, 从而不满足可视化条件。为此, 本文在层推进法的基础上, 采用一种新的计算层推进矢量方法生成三棱柱网格。同时对影响生成网格质量的因素如棱柱的层数和相邻棱柱层之间的比例系数进行相关研究, 在对离心泵中的叶轮及蜗壳划分的网格结果中显示, 棱柱层数越少, 网格质量越高, 且本文算法生成的混合网格比采用已有算法生成混合网格的质量要高。

1 混合网格生成

1.1 算法流程

混合网格的生成包括黏性区三棱柱网格的生成及非黏性区四面体网格的生成。三棱柱网格层推进的方法生成, 具体步骤如下:

(1) 表面网格的生成, 它以非结构化三角形作为三棱柱网格的初始面;

(2) 对物面的每个点计算推进矢量γ, 具体见1.2节;

(3) 对推进矢量进行光顺;

(4) 对物面上的点沿着推进矢量的方向进行推进;

推进步长由下列公式进行控制:

式中:δl为第l层黏性网格的空间分布长度;δ1为黏性网格的第一层的空间分布长度, 根据计算的需要由经验确定, 它必须保证与壁面相邻的第一个节点布置在湍流区, 同时为了较好的捕捉湍流边界层内的变化, 与壁面相邻的第一个节点还应同时位于湍流边界层内, 且湍流边界层内应包含若干节点, 这也要求了网格层数至少要有3层, r代表预先给定的网格压缩因子, 即相邻棱柱层之间的比例系数, 由于三棱柱网格层数及比例系数影响着网格质量及网格划分效率, 因此有必要对其进行测定, 具体见1.3节。当物面上所有的点都向前推进一步后就得到了推进面上的点;

(5) 把这些点按物面上的拓扑关系连接起来就得到一层网格, 并把当前推进面上的点当成物面上的点;

(6) 重复上述过程, 从而得到我们需要的三棱柱网格。

对于非黏性区的四面体网格按文献[11]中的方法生成, 从而生成混合网格。

1.2 节点推进矢量的确定

通常采用计算节点推进矢量的方法是通过计算与该节点相连的所有面元法向矢量的平均值,

式中:γP即为节点P处的推进矢量;ni是与P点相连接的第i个三角形面元的单位法外向矢量;N为与P点相连的三角形面元的个数。

由于离心泵结构复杂, 特别在蜗壳隔舌处以及叶轮进出口处, 采用该方法计算生成的推进矢量很容易发生与当前流行面相交的情况, 即不满足可视化条件。因此, 本文采用一种新的测定节点推进矢量的方法, 其涉及的参量如图1, 计算推进矢量的步骤如下:

(a) 通过寻找面元N1, N2, 确定流形中相交的构成最严重楔形的两个面。

式中:i, j为流行中任意的两个面。

(b) 创建面N1, N2的相交矢量I1。

(c) 构建矢量B位于面N1和N2的平分面上, 且其垂直于I1。

(d) 构建与矢量I1平行的线I, I满足下列条件:

(e) 寻找所有面元与I1相交的点P′t±, t±代表点位于P′的两边;

(f) 计算节点P推进矢量γP。

所有节点的γP值可以在开始进行推进之前求出来, 当推进了一层网格后, 所有节点的γP值要重新计算。

2 比例系数的确定

本文对影响网格质量的棱柱层数及相邻棱柱层之间的比例系数进行研究。图2 (a) 和图2 (b) 分别为叶轮和蜗壳中划分的网格质量随三棱柱层数及相邻棱柱层之间的比例系数的变化情况。从中可以看出, 图2 (a) 和图2 (b) 的网格质量的变化趋势有着很大的相似性, 棱柱层数越少, 网格质量越高, 在棱柱层数较少时, 网格质量随比例系数的增大而减小, 在棱柱层数较多时, 网格质量随着比例系数的增大, 先变小, 后变大。此外, 当比例系数增大到一定程度时, 网格质量在网格层数较多时出现保持不变的情况, 这是由于生成的三棱柱之间发生了重叠, 继续增加棱柱层数对生成的网格已经没有意义。

图3 (a) 、3 (b) 分别为叶轮和蜗壳中网格单元数随棱柱层数及比例系数的变化情况, 从图中可以看出, 两幅图变化趋势相近, 其中在三棱柱层数相同时, 网格单元数基本上是随着比例系数的增加而减少, 而在比例系数相同且较小时, 网格单元数随着棱柱层数的增加而增加。综上所述, 三棱柱层数越少, 网格质量越高, 网格单元数越少, 因此在数值算例中应选取较少的棱柱层数;而在棱柱层数较少时, 随着比例系数减小, 网格质量提高, 但此时网格单元数变多。因此, 基于计算精度和计算效率考虑, 推荐层数为3, 比例系数为1。

3 算例验证

为了验证本文算法生成网格的可靠性, 采用本文算法对叶轮和蜗壳分别划分网格, 通过上述分析, 划分的三棱柱网格层数选取为3, 相邻三棱柱之间的比例系数为1, 并与采用已有计算层推进矢量方法生成网格的质量统计结果进行对比。

图4和图5为采用本文算法与已有算法分别对叶轮及蜗壳划分的网格剖视图及其局部放大图。从局部放大图中可以看出, 采用本文算法时, 除了在生成三棱柱网格质量较差的第一层外, 层推进矢量几乎都垂直于层推进面, 正交性效果较好, 而已有算法的层推进矢量基本上都不垂直于推进面, 这很可能导致网格质量的下降。

通过开源的Gmsh网格软件对划分网格的质量进行测量, 其中网格衡量准则是基于体积与所有边长总和的三次方的比值, 表1为划分网格质量的统计结果, 可以看出, 采用本文算法对叶轮及蜗壳划分的网格单元数和网格节点数较已有算法稍少些, 且采用本文算法生成的网格质量要高于已有算法生成的网格质量。

4 结语

(1) 在比例系数相同的情况下, 划分三棱柱网格的层数越少, 网格的整体质量越高, 生成的网格单元数也越少。

(2) 在棱柱层数较少时, 网格质量随着相邻棱柱之间比例系数的增大而减小;在棱柱层数较多时, 网格质量随着比例系数的增大先减小后增大;在棱柱层数相同的情况下, 网格单元数随着比例系数的增大而减小;且比例系数相同且较小时, 网格单元数随着棱柱层数的增加而增加。

(3) 综合上述分析, 选取棱柱层数为3, 比例系数为1。通过与采用已有的计算层推进矢量方法生成网格相比, 得出本文算法生成的网格有着较好的正交性, 且质量较高。

参考文献

[1]Espadafor F J, Villanueva J B, García M T, et al.Experimental and dynamic system simulation and optimization of a centrifugal pump-coupling-engine system.Part 2:System design[J].Engineering Failure Analysis, 2010, 17 (7) :1 551-1 559.

[2]王洋, 王洪玉, 徐小敏, 等.冲压焊接离心泵叶轮有限元计算[J].排灌机械工程学报, 2011, 29 (2) :109-113.

[3]刘厚林, 董亮, 王勇, 等.流体机械CFD中的网格生成方法进展[J].流体机械, 2010, 38 (4) :32-37.

[4]Evazi M, Mahani H.Unstructured-coarse-grid generation using background-grid approach[J].SPE Journal, 2010, 15 (2) :326-340.

[5]ZHANG Laiping, CHANG Xinghua, DUAN Xupeng, et al.Applications of dynamic hybrid grid method for three-dimensional moving/deforming boundary problems[J].Computers&Fluids, 2012, 62:45-63.

[6]Ito Y, Shih A M, Soni B K.Hybrid mesh generation with embedded surfaces using a multiple marching direction approach[J].International Journal for Numerical Methods in Fluids, 2011, 67 (1) :1-7.

[7]Luo H, Spiegel S, Lhner R.Hybrid Grid Generation Method for Complex Geometries[J].AIAA journal, 2010, 48 (11) :2 639-2 647.

[8]Pirzadeh S.Three-dimensional unstructured viscous grids by the advancing-layers method[J].AIAA journal, 1996, 34 (1) :43-49.

[9]Lhner R.Matching semi-structured and unstructured grids for Navier-Stokes calculations[J].AIAA paper, 1993:93-3 348.

[10]SANG Wenmin, LI Fengwei.An Unstructured/structured Multilayer Hybrid Grid Method and Its Application[J].Chinese Journal of Computational Physics, 2004, 21 (4) :345-351.

网格生成 篇7

自适应有限元分h型、p型、h-p型三种。对于h型自适应有限元,最重要的问题有两个:(1)对现有计算结果进行误差估计,以判断是否进行网格调整并给出网格调整的标准;(2)自动生成符合要求的有限元网格。

网格自动生成的研究一直是一个热点,其中关于四面体单元自动生成的研究进步最为明显。目前国际上流行的四面体生成算法为行波法和Delauney法。应用行波法生成的的四面体网格常有一些形状极差的单元,这对有限元求解的精度和效率都有影响,于是通常一个有效的行波法网格自动生成过程应包括两部分,先形成一个拓扑关系正确的初始网格,然后对其进行优化,以得到一个高质量的网格供计算使用[1]。

经过20多年的发展,h型自适应有限元在机械制造、金属成型等工业领域都有成功运用,近年来也有学者尝试在岩土工程等领域引入自适应有限元分析。本文的主要目的就是结合岩土工程等领域的实际情况,讨论利用行波法生成拓扑关系正确的四面体网格的若干关键问题,至于网格优化的问题可参见论文[2]。

行波法首先对三维物体的轮廓面进行三角形离散以形成初始波前,然后生成四面体,将波前逐渐向里推进。当波前上的活动面数为零,四面体充满整个三维域时行波法终止[3]。在这个过程中,有两个问题对行波法的有效性和效率起重要作用,这两个问题分别是不规则曲面的参数表示和几何搜索问题。

物体轮廓面的三角形离散是一个三维问题,通常的办法是先将轮廓面参数化表示为:

然后在参数面(s,t)与三维曲面之间建立双向映射,将问题转化为二维三角形行波法处理[4]。显然不规则曲面的参数化表示是这个过程的关键,本文的第二部分将讨论利用混合样条插值曲面来解决这个问题。

行波法通过几何搜索寻找在一特定的空间区域的几何对象(点、线、面等),这在四面体生成过程中有重要应用。行波法首先将在波前面上确定一三角形作为活动面Δabc。若此时波前上共有n个点,运用几何搜索过程从这n个点中找出位于Δabc某一规定邻域内的点,然后在这些搜索得出的点中确定一点d,形成四面体Tabcd(另一种生成四面体的方法是新生成一点d',形成四面体Tabcd')。由此可知几何搜索是行波法生成每一个单元都必须处理的问题。如果直接对波前上的每一个点依次进行判断,所需的工作量通常较大,行波法的很大一部分计算时间都被几何搜索占用。为了提高行波法的效率,本文第三部分将讨论利用双叉树型结构进行几何搜索的问题。

2不规则曲面的参数化表示

2.1混合样条插值曲面的引入

岩土工程中存在着大量不规则曲面,如地形面、大规模断层面等,如何将这些曲面参数化表示一直是一个困难的问题。计算机辅助几何设计(CAGD)中一些常用的曲面构造方法(如孔斯曲面(Coons Surface))对上述大尺度且几何信息不充分的曲面往往是不适用的。这些方法往往要求许多面片的角点坐标X、切矢Xs、Xt和扭矢Xst等数据,这些要求显然是实际工程所不能满足的。若采用一些近似方法(如Adini处理),又可能导致曲面较严重的变形。在这种情况下,本文将W.J.Gordon[5]所提出的混和样条插值曲面引入岩土工程领域,以期解决上述问题。

2.2混和样条插值曲面的构造

令R:[a,b][cd]为(s—t)平面内一矩形域,定义于该矩形域的连续曲面记为:

将区间I=(a,b)分为π:a=s1≤s2≤…≤sn=b,将区间I'=(c,d)分为π':c=t1≤t2≤…tn'=d,垂直于(s—t)面的系列平面s=si,t=tj与曲面Z的交线记为:

这n+n'个单参数函数给出了曲面Z上的一套曲线网络。

选定两套2p-1次样条基函数{φi(s)},{φj(t)}。

混和样条插值曲面可表示为

式(8)中对曲面式(2)运用的算符L、M为一线性操作,具体意义可见文献[5]。

2.3混和样条插值曲面在岩土工程中的应用

通过以上叙述可知,混和样条插值曲面对原曲面进行的是一种整体拟合,这从根本上有别于常见的曲面片拼接型的曲面构造,适合于对大尺度曲面进行拟合。同时若令p=2,该曲面保有C2的连续性,这是不规则曲面三角离散所必须的。又曲面边界的切矢条件在岩土工程中并无特别意义,在这里将其忽略:

将式(7)带入式(6),得到新的混和样条插值曲面公式:

式(8)中样条基函数可以作为固定模块直接编入程序,根据不同的n或n'直接调用。如果已知曲线网络式(3)上的一些点,利用三次样条插值曲线便可拟合出n+n'条曲线参数方程fj(s)、gi(t),利用公式(8)便可得到曲面式(2)。将曲面参数方程(1)看作定义于同一矩形域的三个独立的曲面方程,依次用上述方法拟合,从而得到不规则曲面的参数表示。

上述简化处理会使曲面临近边界的区域发生一些变形,但与机械、数控操作等领域所研究的精密几何形体的模拟不同,这些误差在岩土工程中是可以接受的。如果选择一个恰当的辅助面,这种影响还会得到更好地控制[6]。

3行波法中的几何搜索问题

3.1双叉树型结构的基本概念

双叉树型结构是一种数据在存储空间非顺序排列的数据结构。它由许多层状结构构成,其中最基本的层状结构为一个节点(母节点)与其下左右两个节点(子节点)相连接,这种结构方式也称为母节点指向两个子节点。在整个树型结构中,除了最上层节点(根节点),每个节点都被且仅被一个节点所指向。为了充分定义某节点,需三个相关数字,首先是存贮于该节点的数据,然后是该节点两个子节点在存贮空间中的位置,后两个数字记为该节点的左连接、右连接。显然左、右连接或者为零,或者为整数。若一个节点左、右连接均为零,即无子节点与其相连,则该节点称为中断点。

3.2双叉树型结构在几何搜索中的应用

设行波法波前上有n个点,将其存贮于树型结构的n个节点上。利用某种规则确定树型结构,使节点与空间某域建立联系,通过对树型结构节点的访问实现几何搜索[7]。

在N维空间RN中定义一超立方体为ai≤xi≤bi(i=1,…,N),并将该立方体上下界记为(a,b)。在树型结构中定义根节点为“0层点”,其对应空间为(a,b)。与树型结构的分叉与层状结构对应,依次用垂直于坐标轴的平面切分空间(a,b)。设节点K位于第L层,其代表的空间子域为(ck,dk),节点K的左连接对应子节点为P,右连接对应子节点为Q,它们分属的子空间为(ckl,dkl)、(ckr,dkr)。

令j=1+mod (l,N) (9)

P对应的子空间为:

c${}_{kl}^i$= c${}_k^i$,d${}_{kl}^i$= d${}_k^i$;i≠j (10)

$c{}_{kl}^i=c{}_k^i,d{}_{kl}^i=\frac{1}{2}(c{}_k^i+d{}_k^i)$; i=j (11)

c${}_{kl}^i$≤x${}_{{\rm P}}^i$≤d${}_{kl}^i$ (12)

Q对应的子空间为:

c${}_{kr}^i$= c${}_k^i$,d${}_{kr}^i$= d${}_k^i$;i≠j (13)

$c{}_{kr}^i=\frac{1}{2}(c{}_k^i+d{}_k^i),d{}_{kr}^i=d{}_k^i$;i=j (14)

c${}_{kr}^i$≤x${}_Q^i$≤d${}_{kr}^i$ (15)

其中i=1,…,n。

由以上定义可知,K节点的所有子节点存贮的数据均位于K所代表的空间(ck,dk)内。若(ck,dk)与搜索空间不相交,则K所含的所有子节点均可排除,不必一一检验。下面具体给出利用双叉树型结构进行几何搜索的算法,该算法略经修改也可实现对树型结构任一节点的访问。

定义一栈数组stack,栈长度m=0,root=根节点存贮位置。

(a) 访问位于root的节点,检查存贮于该节点的数据xk是否在搜索空间(e,f)中,即检查是否满足ei≤x${}_k^i$≤fi(i = 1,N)。

(b) 如果root对应的节点右连接right≠0,且(ckr,dkr)与(e,f)相交,即,令n=n+1,stack(n)=right。

(c) 如果root对应的节点左连接left≠0,root=left go to (a) 。

如果left=0,但n≠0,root=stack(n),n=n-1, go to (a) 。

如果left=0, 且n=0,结束程序。

在四面体生成过程中,行波法的波前不断变化,这就要求从树型结构中不断增加或删除节点。利用上述访问节点的算法可以方便地实现这一点。

最简洁的添加节点的方法是从根节点开始由上至下访问树型结构的节点。如果某节点左(右)连接为0,且需添加点坐标位于该连接所代表的空间区域内,令该连接指向需添加节点即可。

至于删除节点,若被删除点为一个中断点,只需指向该节点的连接为0即可。非中断节点的删除较复杂,需从被删除节点所有的子结构中找出一中断点,令指向被删除节点的连接指向该中断点,指向该中断点的连接为0。

4算例

运用上述方法对一边坡洞室开挖问题进行了有限元网格离散。该边坡有不规则地形面和两条相交断层。为了保证计算精度,根据自适应有限元误差估计理论对网格做了优化。

图2、图3给出的是经过两轮优化后的网格,其网格离散误差为12.3%,有四面体单元13 114个。

摘要：网格自动生成是自适应有限元重要组成部分,其中如何有效地对三维物体轮廓面进行参数化表示一直是一个难点。结合岩土工程的实际情况,利用混合样条插值曲面对地形面、大规模断层面等不规则曲面进行了拟合,成功地解决了这个问题。另外还讨论了利用双叉树型结构进行几何搜索的算法,该算法能提高行波法生成四面体的效率。最后用一边坡洞室开挖的算例验证了编制程序的有效性。

关键词：有限元网格,四面体,行波法

参考文献

[1] Dari E A,Buscaglia G C.Mesh optimization:how to obtain good un-structure 3-D finite element meshes with not-so-good mesh generators.Struct Optim,1994;8,181—188

[2] Zavattieri P D,Dari E A, Buscaglia G C. Optimization strategies in unstructured mesh generation. Int J Numer Meth Eng,1996;39:2055—2071

[3] Peraire J,Peiro J.Adaptive remeshing for three-dimensional compres-sible flow computations,Journal of Computational Physics,1992;103:269—285

[4] Moller P,Hansbo P.On advancing front mesh generation in three di-mensions.Int J Numer Methods Eng,1995;38:3551—3569

[5] Gordon W J. Spline-blended surface interpolation through curve net works. J Math Mech,1969;18,931—952

[6] Xia H X,Chen S H.3-D adaptive FEM in rock slope stability analy-sis.The 10th International Conference on Computer Methods and ad-vances in Geomechanics,Arizona,USA,2001:1635—1639

网格生成 篇8

地下水和溶质在地下的运移过程相当复杂,涉及到气-液-固、热和化学反应、饱和非饱和流动等问题。地下水资源评价和地下水污染防治与管理需要正确了解多相流体中地下水流的运动和溶质的运移,在考虑地下水运动特点的基础上,适当忽略非饱和带水流和气相固相态的运动过程则可大大地将研究问题简化,但很多情况多相流动、饱和非饱和流动不能简化,因此,多相流的数值模拟作为分析水-气-固运动的工具显得尤其重要。

TOUGH2模拟器是多相流模拟软件中最知名、应用最广泛的软件之一,由美国加州大学伯克利分校劳伦斯实验室开发[1,2],它采用积分有限差方法进行空间离散,能进行非等温多组分多相流体数值模拟,可以模拟一维、二维和三维孔隙和裂隙介质中流动,在地下水流、地热、CO2地质封存、天然气水合物等模拟方面应用广泛。其并行版TOUGH2-MP[3]是目前文献中已报道的唯一成功应用于大规模精细模拟(具有几十万到上千万网格的模型)的软件,已成功地应用于多个实际示范和研究工程[4,5]。TOUGH2模拟器的缺点在于它的使用比较困难,需要大量的专业知识和长期的经验才能全面成功应用,本文将针对TOUGH2模拟器应用中最棘手的网格问题提出解决方案,供广大TOUGH2模拟器应用者参考借鉴。

1 已有的网格生成工具评述

TOUGH2从最初20世纪80年代产生到现在,出现了几款网格生成工具,比较代表的程序和软件包括MeshMaker、AMesh、PetraSim、MView和WinGridder等,下面分别进行介绍。

1.1 MeshMaker程序

该网格剖分器最初是嵌在TOUGH2程序中,在后期的TOUGH2版本中,MeshMaker单独成为一个剖分器。它可以将研究区域离散为矩形、同心圆、扇形弧段等网格单元,平面上相邻网格之间的间距可以灵活设置。程序在剖分时还可将不同区域设置为不同的材质代号,方便后期的渗透率等参数的分区赋值,最终程序生成TOUGH2程序运行需要的网格文件,包括ELEME和CONNE数据。正因为MeshMaker程序的这些优点,TOUGH2模拟软件集成了MeshMaker程序,供用户方便进行网格设计。该剖分器对于规则边界和等厚的地质模型具有很好的适用性,但对于不规则边界和复杂的变厚度问题,该剖分器显得无能为力。

1.2 AMesh

它是Haukwa(1999)开发[6]的开放源代码的网格剖分器,能直接生成TOUGH2模拟器运行所需的网格文件。在给定一些点信息或泰森多边形的中心,它可以生成一维、二维和三维网格。缺省的输入文件为IN文件,输出文件包括ELEME、CONNE和SEGMENT文件,前两个输出文件可直接作为TOUGH2模拟器的网格文件,后一个文件可用于后处理的制图。它没有可视化的用户界面,操作起来比较困难。

1.3 PetraSim

PetraSim是用于TOUGH2模拟程序家族的具有图形化界面的商业软件。它以交互式3D环境,包括网格剖分,参数定义和结果展示,使建模者方便应用TOUGH2模拟器的功能。软件集成网格生成、网格和单元格编辑等功能,可自动生成模型的输入文件,操作起来方便,目前版本只能处理规则边界,而且价格昂贵。

1.4 MView

MView软件最初是由Intera公司针对尤卡山(Yucca Mountain)项目开发的TOUGH等系列数值模型前后处理和可视化操作界面,后来扩展到MODFLOW等软件中,集成了钻孔编录、物探、采样分析等野外数据的分析。MView软件目前已成为一款商业化的数值模拟支持系统,包括数据分析和复杂地质数据可视化等功能,包括TOUGH2模拟器的前后处理和可视化工作。

1.5 WinGridder

WinGridder是Pan开发的基于视窗操作系统的网格剖分软件[7],它提供了有效的、友好的界面,集成了网格处理软件,能处理单一连续体和双重连续体。它首先在图形界面下将研究区域划分为二维网格,然后将裂隙和局部区域细化,再拓展为三维网格。输入文件包括域范围、地质层、岩石属性、裂隙和钻孔信息,由基础信息划分网格。软件价格昂贵,而且不容易学。

由上述分析可知,商业软件价格昂贵,而且学习困难。共享软件大多数功能仍不齐全,有必要进行对网格生成进行改进。

2 基于多边形网格的网格生成方法

TOUGH2模拟器的网格可以是任意的,只要按照模型输入文件格式给出网格单元的体积和中心点(图1)。图1(a)和图1(b)显示了典型的TOUGH2模拟器的垂向和平面上网格关系,即必需要知道相邻单元的接触面积和侧面面积等信息。从相邻单元接触面关系来看,很容易想到依据泰森多边形来划分单元,而将研究区域离散为多个三角形网格,然后按照三角形形成泰森多边形,再展开就成为多边形网格,则可满足TOUGH2模拟器网格的特点。TOUGH2模型运行的网格包括ELEME单元和CONNE连结信息,单元可以为5或8字符命名,文件具体格式按固定的输入格式进行如表1所示。

2.1 平面三角网格的生成

对于二维平面三角网格自动生成,从20世纪70年代开始,国内外学者[8,9,10]做了大量的研究,相关技术相当成熟,有Delaunay法、波前推进法、正规法、向量偏移法、有限四(八)叉树法等等。Delaunay法因为方法简便,对于带约束条件的三角剖分更显得灵活自由,同时容易编程实现,适合作为区域三角剖分的方法,然而目前的方法仍不能较好地考虑相邻单元的垂心在单元内,即保证三角形为非钝角三角形。笔者曾对带约束的任意多边形Delaunay三角剖分算法进行了改进,基于结点尺度设计,提出了一种简单的平面区域结点尺度自动确定的方法,在三角单元形成过程中,特别地考虑了三角形状的自动修正和编辑问题[11]。而且大多数有限元软件都有平面三角形网格生成的功能,可完成区域的离散化工作。图2(a)显示了由三角形网格转换为多边形网格的转换。

注:A为字符变量,X为空格,I为整型变量,E为实型数

2.2 多边形网格数据的存储格式

平面上,多边形网格的数据存储信息包括:(1)点信息,为多边形每个顶点的坐标信息,包括x,y,z方向坐标值;(2)每个顶点的侧边结点总数,简记NPSD;(3)每个顶点的侧边结点序号,简记MPSD。而三角形网格的数据存储信息只需要:(1)点信息;(2)组成三角形的各个结点序号及按逆时针排序,即IJK信息。NPSD为一个结点周围的结点总数,MPSD为一个结点按逆时针排列的侧边结点序号,一个多边形由一个NPSD及对应的MPSD组成。NPSD和MPSD的数据信息可由三角网格的IJK结点拓补结构生成。

对于NPSD和MPSD信息的生成,先判断此结点是内结点还是边结点,判断的方法就是按照含此结点的三角形个数与含此三角形的边的个数是否相等,如果两者相等,说明此点为内点;不等则说明此点为边点。如图3(a)所示,对于内结点,先找到一个含此点的边和三角形Ep7,设给定点结点号为O,此边的另一结点号为P1,找到的三角形结点顺序号和OP1序号相应,将边的另一结点号P1作为起始结点序号,对于Ep7单元的结点号P1,设此单元的最后一个结点号为P2,找到与它对边相邻的单元号Ep1,然后将P2作为序列号的第二点;依照这样的方法,找到P3,P4,P5,P6,P7结点号。对于内结点,结点号封闭,此结点的最终Mpsd序号为P1P2P3P4P5P6P7P1。

对于边结点,可采用类似的方法,不同的是先找起始点时,先找到边的一点,如图3(b)所示,有两个边点,P1和P7,对于含O和P7的单元,其单元另一结点号为P6,要求OP7P6逆时针排列,如不是逆时针排列,说明这一点不是要找的边序号起点,从图中可以看出,OP7P6不是逆时针排列,P7不能作为边结点O的起始点;相反,对于P1点,OP1P2逆时针排列,说明P1正好是边结点O的起始点。

平面上生成的多边形网格单元数与结点总数相同,其编号则按结点号编号。垂向上,每层的网格保持相同,按从上层到下层排序,下一层单元编号则在上一层单元编号加上单元总数。而每个单元的平面面积、总体积与相邻单元的接触面积可由图2(b)中的关系计算而成[12],则可生成ELEME和CONNE的单元、面积等信息。

2.3 参数分区及数据生成

模型仍需要介质参数分区文件及辅助的数据生成,可根据介质的水力性质对其进行分区,需要考虑不同岩性的参数分区情况。在介质分区信息生成过程中,设计的分区信息生成的方法是按结点赋值,即在一个封闭的多边形内的结点可统一赋成一个分区号,这些点形成的泰森多边形即为这些点控制的分区区域,离散化的边界只能逼近实际边界,设计的算法是符合TOUGH2模拟器数值计算要求的,如图4所示。

(虚线区域为参数区,灰色区域为格点形成的参数分区)

3 方法的应用

为说明方法的正确性,以北京市平原区地下水流数值模拟为例,将多边形网格生成方法应用于模型中。按照研究区主要水文地质条件以及考虑到实际收集地下水监测信息的代表性[13],将研究区划分为6层水文地质结构,分别为种植土层弱透水层、潜水含水层(砂砾石)、弱透水层(砂粘土、粘土)、第一承压含水层(砂土)、弱透水层(粘土等)、第二承压含水层(砂土),取第四系底为模型底边界。根据已有的钻孔信息,已获知每一模拟层的厚度。

平面上先以北京市平原区为外边界,包括河流、抽水井和观测井点,利用三角网格剖分器或有限元软件FEFLOW[14]将区域离散为多个三角形,然后再通过钝角修正,确保泰森多边形的结点为三角形网络的垂心在三角形之内。平面上离散的泰森多边形为7620个单元,在模型感兴趣重点模拟处进行加密剖分。每个模拟层的网格保持一致,位于地面的多边形网格高程即由已知点的高程信息数据插值而成,插值可借助SURFER软件完成,其它各个模拟层的结点高程也可由上一层结点的高程和该层已知点的厚度信息插值生成,多边形网格实体模型即已完成。

参考北京市多年的地下水研究的水文地质资料可获得每一层水文地质结构的水文地质参数分区信息,以顶层为例,说明参数的生成方法,参数分区如图5所示。结点的数据生成可借助ARCGIS信息平台完成,即对每个分区建立区文件,位于区文件之内的结点赋给对应的岩性分区号;然后对于其它模拟层也采用同样的方法,对于统一的岩性分区给定同样的标识号,在生成ELEME信息时可给出每个结点的岩性分区号,多边形网格和参数分区如图6所示。

(1、2、3、4分别代表种植土、粘泥和黑泥土、细砂土、种植土层和水泥土)

网格及对应参数分区(数字为分区号)

4 结论

TOUGH2模拟器作为强有力的多相流模拟软件,困难的网格生成方法影响了TOUGH2模拟器的方便快捷应用,本研究提出了新的网格生成和参数分区方法,即在三角形网格的基础上,先将三角形转换成多边形网格,然后再利用地理信息系统软件进行参数分区的结点赋值,再利用统计软件SURFER将已知点的高程和厚度信息扩展到模拟层的结点,完成了TOUGH2模拟器运行需要的ELEME和CONNE数据信息。整个过程笔者已在VISUAL C++平台上完成数据生成的集成处理,利用该程序可方便快捷生成网格。本文提出的方法适用于已有三角网格或将研究区域离散为多个三角形的情形,对于四边形或同心圆网格,则需将其转换为三角单元网格方可转换。另外对于三维模型中含孔洞或隧道的剖分问题,需要基于已有的三维网格和剪裁区域进行剪切,本文没有进行深入讨论。

摘要：针对TOUGH2模拟器网格数据较难生成的特点,在总结分析已有的网格生成工具的优缺点基础上,介绍了基于多边形的网格生成方法。该方法先将研究区离散为多个三角形网格,再将其转化为多边形网格。基于水文地质参数分区和地理信息系统软件,可生成各结点的水文地质参数分区信息,再利用统计软件将已知点的高程和厚度信息扩展到模拟层的结点,最终可生成TOUGH2模拟器需要的网格输入文件。以北京市平原区网格生成说明了方法的合理性。这些成果为TOUGH2模拟器的推广应用提供了参考,并将提高TOUGH2的使用效率。