六面体网格生成办法
合集下载
相关主题
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
2 α = xη + yη2 ,
β = xξ xη + yξ yη ,
γ = xξ2 + yξ2
和代数法相比，偏微分方程法的计算较为复杂，需要求解偏微分方程，右端强迫函数的选取 也不方便，可以通过强迫函数的选取得到希望的映射网格。
・21・
上海交通大学博士学位论文
3.2.3 超单元映射法
等 参 超 单 元 网 格 生 成 法 是 由 Zienkiewicz 和 Phillips 最早提出的[31,42]，其基本思想是将划分区域分 成更简单的子块（超单元） ，子块被映射为自然坐标系 中的单位直角形（正方形） ，根据每个方向给定的分级 权系数，将自然坐标系中的直角形（正方形）离散并 反向变换回超单元。在三维情况使用超单元生成法时， 要划分的实体用 20 节点六面体单元分成数个子区域 （图3-3） 。在超单元内一点的坐标 x，y 和 z 和自然坐 标之间的关系由下式确定
计算区域 物理区域
ζ η ξ
z y x
图3-1计算区域和物理区域 Fig. 3-1 Computational doma in and physical domain 计算网格 物理网格
x = x(ξ, η, ζ) y = y(ξ, η, ζ) z = z(ξ, η, ζ)
ζ η
z
ξ
y x
图3-2计算区域网格和物理区域网格 Fig. 3-2 Grids in the computational domain and in the physical domain
(3-2)
右端的强迫函数为
N ξ − ξ m −c m ξ −ξ m e − ∑ bi ξ − ξm m=1 i =1 M N η − η m − cm η−ηm Q (ξ ,η ) = − ∑ am e − ∑ bi η − ηm m =1 i =1
P (ξ , η ) = − ∑ a m
M
ξ − ξ m − d i [(ξ −ξ i )2 +(η−ηi )2 ]2 e ξ − ξm 1 η − η m −d i [(ξ −ξi )2 +(η−ηi )2 ]2 e η −ηm
3.2.2 偏微分方程法
类似于代数法，偏微分方程法也是通过映射将曲的边界映射为网格坐标线，只不过这里使用 的映射函数满足某一类微分方程。以二维椭圆偏微分方程为例[41]，满足最大条件的椭圆型微分方 程组可以为
ξ xx + ξ yy = P(ξ , η ) η xx + η yy = Q(ξ , η )
x = ∑ N i xi , y = ∑ N i y i , z = ∑ N i z i
i =1 i=1 i =1 20 20 20
20 5 16 13 12 1
8
19 18
ζ η ξ
7
17
6 15 14 3
4 9
11 10 2
(3-5)
式中，Ni(ξ, η, ζ)是与按曲线坐标系ξ，η和ζ定义的与每 个节点相联系的形函数， ξ， η和ζ的取值范围为−1 到 1。 20 节点六面体单元形函数在各角点定义为
1
(3-3)
在变换空间中，这些函数写成为
αxξξ − 2 βxξη + γxηη = − J 2 (Pxξ + Qxη ) αyξξ − 2 β yξη + γyηη = − J 2 (Pyξ + Qyη )
(3-4)
这里 J 是变换的 Jacobi 矩阵行列式的值，J = xξxη− yξyη，而
(
)
(
)
(3-7)
(
)
如节点的坐标已知，则与特定坐标ξ，η和ζ相对应的点的坐标可由式(3-5)计算出。 沿每个局部自然坐标，给定划分的数目和划分的间距，这些子区域（超单元）被划分为更小 的单元。以 ξ为例，划分时自然坐标按下式增加
ξ i = ξ i −1 + 2 (Wξ )i WξT
(3-8)
Baidu Nhomakorabea
这里WξT = ∑ (Wξ )i ，ξ0 = −1。 超单元映射法计算简单，通过对自然坐标的不同划分，可以得到密度变化的网格，但超单元 描述的边界和划分区域的实际边界存在误差，不能适应复杂的边界。
3.3.1 映射法
当划分区域的几何特性合适时，四边形或六面体映射网格生成法可以得到相当好的网格[56]。 虽然映射法被认为是结构化网格生成法，但许多非结构网格划分软件都提供这种方法。当应用映 射法时，要求划分区域的对边具有相同的子划分。在三维情况下，将要划分区域的拓扑结构的立 方体的每个相对的面必须有相同的表面网格。这种情况在几何形状复杂是不可能达到，或者需要 人工干预将几何区域划分成可映射的子区域并分配适当的间隔。 Yang 等[57,58]使用的方法是将划分区域分成子区域，每个子区域对应于不同的模块，每个模 块中的网格已预先划分，模块中的网格映射为子区域中的实际网格。这里模块中的网格可以是结 构化的，也可以是非结构化的，但要求相邻子区域的交接面上的表面网格具有相同的拓扑结构。 该方法的特点是可以针对子区域的几何特性和物理特性选择相应的映射模块，生成的网格质量 高，缺点是通用性差。提供映射法生成网格的软件有 ANSYS（ANSYS Inc.）、CUBIT （Sandia National Laboratories ）、goemagic Wrap（Raindrop Geomagic Inc.）、Hypermesh （Altair Computing Inc.）等。
・22・
第 3 章 六面体单元生成方法及相关技术
3.3 非结构网格四边形/六面体网格生成法
由于算法的特点，自动非结构网格生成算法主要是三角形或四面体网格划分，有许多文献讨 论了三角形或四面体网格划分方法和相关的技术[43,44,45,46,47,48,49,50,51,52,53,54]， 许多软件也提供了三角 [55] 形或四面体网格划分功能 。非结构网格生成的特点是能适应复杂的边界，需要的人工干预少， 计算费时。 Steven. J. Owen[55]总结了当前学术界和工业界使用的网格划分方法和可以利用的软件。 非结构六面体网格生成方法主要有映射法、转换法、铺层算法和须段编织算法等。不象三角 形网格划分算法向四面体网格划分算法扩展那样直接，从二维四边形网格划分算法向三维六面体 网格算法的扩展并不是直接的。
・20・
第 3 章 六面体单元生成方法及相关技术
x (ξ , η , ζ ) F (ξ ,η , ζ ) = y (ξ , η , ζ ) z (ξ ,η , ζ )
(0 ≤ ξ ≤ 1,
0 ≤ η ≤ 1, 0 ≤ ζ ≤ 1)
(3-1)
多变量的坐标变换式(3-1)使用无限插值(Transfinite Interpolation)来完成。在三维情况下，使 用混合函数和与之相联系的参数（特定点的位置及偏导数）来显式地决定式(3-1)，然后通过对每 个单变量的循环完成无限插值。一般来说，这些混合函数和参数都选取为区域边界处的函数和参 数。文献[40]给出了几种无限插值函数及其参数的选取。 代数法无限插值网格生成法的特点是计算简单，可以采用中间变量的方法方便地控制网格的 密度，对边界简单的区域，可以生成质量较高的网格，但缺点是不适应复杂边界的划分，边界不 规则时生成的网格的质量很差，并可能产生奇异性。可以通过将划分区域分解为子区域的方法， 在子区域上应用，可以在一定程度上克服这些缺点，但不易实现自动划分。
扫描线
目标
源
图3-4 扫描法生成六面体单元 Fig. 3-4 Hex elements generated by sweeping
3.2.1 代数法网格生成
代数法网格生成技术[40]是通过插值函数将理想的直角坐标系表示的计算区域变换为实际物 理区域来实现的（图3-1） 。它的具体做法是将计算区域的直角坐标用均匀的间隔划分为计算区域 的网格，通过变换，计算区域上均匀分布的坐标被映射为物理区域上的坐标（图 3-2） 。网格点数 目的多寡并不影响变换的特性。从计算区域到物理区域的变换是通过一个向量函数完成的
N i (ξ , η , ζ ) =
图3-3 20 节点六面体单元 Fig. 3-3 20-node hexahedral element
1 (1 + ξξi )(1 +ηηi )(1 + ζζi ), i = 1, 2, L, 8 8
(3-6)
在各棱边的中间节点定义为
1 1 − ξ 2 (1 + ηη i )(1 + ζζi ), i = 9, 11, 17, 19 4 1 N i (ξ , η , ζ ) = (1 + ξξi ) 1 − η 2 (1 + ζζi ), i = 10 , 12, 18, 20 4 1 N i (ξ , η , ζ ) = (1 + ξξi )(1 + ηηi ) 1 − ζ 2 , i = 13, 14, 15, 16 4 N i (ξ ,η , ζ ) =
第3章 六面体单元生成方法及相关技术
3.1 引言
有限元法是求解偏微分方程描述的连续体问题的一种近似工程方法。为克服实际连续体问题 难于处理的问题，它将分析区域离散化，将偏微分方程转化为线性方程组，采用数值计算方法求 出连续体问题的近似解。用有限元法进行工程分析的主要过程包括三个阶段： （1）有限元模型的 建立和数据准备； （2）用软件分析计算； （3）分析结果的判断和评定。迅速而合理地划分有限元 网格是完成有限元分析的前提和保证。 有限元网格生成技术发展到现在，已经出现了大量的、不同的实现方法 [22,23,24,25,26,27,28,29,30,31, 32,33,34,35,36,37,38,39] 。网格生成方法一般可以分为两类：结构化网格生成方法和非结构网格生成方法。 所谓结构化网格，严格的讲，是指网格内部节点具有相同的邻接单元的网格，即网格内部节点具 有相同的拓扑结构。结构化网格生成通常使用复合的迭代平滑技术，用边界或物理区域来排列单 元，划分区域的边界不能太复杂，以便能将区域分解为拓扑结构的块。所谓非结构化网格，放松 了对节点的连接要求，允许一个节点属于任意数量的单元。虽然四边形和六面体单元也可以是非 结构网格，但通常非结构网格是指三角形和四面体网格。目前二维区域网格自动划分发展较为成 熟，许多商业软件都提供了成熟的工具。而三维区域的网格自动划分，除个别商业软件提供了四 面体网格的自动划分外，还不成熟，特别是形状较规范的六面体网格的自动生成。 一般来说，由程序自动划分的网格总有一部分的单元的形状不是很好，需要对网格的质量进 行优化。节点序号的的标注方法直接影响刚度矩阵存储容量的大小从而影响计算的效率，因而需 要对网格的节点重新编号，以得到具有较小带宽的网格系统。 本章主要讨论了三维六面体网格生成的方法及其相关技术。
3.2 结构化网格生成
结构化网格生成方法主要分为两类：代数法和 PDE（Partial Differential Equation，偏微分方 程） 方法。 代数法包括超限插值法、 等参映射法和保角映射法； PDE 方法包括椭圆系统方程、 Poisson 系统方程、抛物线系统方程和双曲线系统方程方法。应用最多的是等参映射法。结构化网格的优 点是易于实现，在每个子区域内网格可以得到很好的控制，生成规则的结构化网格，而且能生成 曲面网格。缺点是对不规则的形体有时生成的网格质量很差，需要事先根据所要产生的网格类型 将目标区域分割为一系列可映射的子区域，这一步通常需要手工完成，因而自动化程度低，不适 合网格全自动生成。提供结构化网格生成的软件有[55]3DMAGGS（NASA，Langley & Lockhead ， 四边形和六面体单元） 、 CAGI （ERC，Miss. State，四边形和六面体单元） 、CSCMDO（CSC/NASA LaRC，四边形单元） 、EMC2（INRIA，四边形单元） 、GENIE++（ERC，Miss. State，四边形和六
・19・
上海交通大学博士学位论文 面体单元） 、ICEM CFD（ICEM CFD Engineering，四边形和六面体单元） 、QulkGrid （Perspective Edge，四边形单元） 、TrueGrid（XYZ Scientific Applications ，Inc.，四边形和六面体单元） 、VGM （NASA，Langley & Lockhead ，四边形和六面体单元）等。
β = xξ xη + yξ yη ,
γ = xξ2 + yξ2
和代数法相比，偏微分方程法的计算较为复杂，需要求解偏微分方程，右端强迫函数的选取 也不方便，可以通过强迫函数的选取得到希望的映射网格。
・21・
上海交通大学博士学位论文
3.2.3 超单元映射法
等 参 超 单 元 网 格 生 成 法 是 由 Zienkiewicz 和 Phillips 最早提出的[31,42]，其基本思想是将划分区域分 成更简单的子块（超单元） ，子块被映射为自然坐标系 中的单位直角形（正方形） ，根据每个方向给定的分级 权系数，将自然坐标系中的直角形（正方形）离散并 反向变换回超单元。在三维情况使用超单元生成法时， 要划分的实体用 20 节点六面体单元分成数个子区域 （图3-3） 。在超单元内一点的坐标 x，y 和 z 和自然坐 标之间的关系由下式确定
计算区域 物理区域
ζ η ξ
z y x
图3-1计算区域和物理区域 Fig. 3-1 Computational doma in and physical domain 计算网格 物理网格
x = x(ξ, η, ζ) y = y(ξ, η, ζ) z = z(ξ, η, ζ)
ζ η
z
ξ
y x
图3-2计算区域网格和物理区域网格 Fig. 3-2 Grids in the computational domain and in the physical domain
(3-2)
右端的强迫函数为
N ξ − ξ m −c m ξ −ξ m e − ∑ bi ξ − ξm m=1 i =1 M N η − η m − cm η−ηm Q (ξ ,η ) = − ∑ am e − ∑ bi η − ηm m =1 i =1
P (ξ , η ) = − ∑ a m
M
ξ − ξ m − d i [(ξ −ξ i )2 +(η−ηi )2 ]2 e ξ − ξm 1 η − η m −d i [(ξ −ξi )2 +(η−ηi )2 ]2 e η −ηm
3.2.2 偏微分方程法
类似于代数法，偏微分方程法也是通过映射将曲的边界映射为网格坐标线，只不过这里使用 的映射函数满足某一类微分方程。以二维椭圆偏微分方程为例[41]，满足最大条件的椭圆型微分方 程组可以为
ξ xx + ξ yy = P(ξ , η ) η xx + η yy = Q(ξ , η )
x = ∑ N i xi , y = ∑ N i y i , z = ∑ N i z i
i =1 i=1 i =1 20 20 20
20 5 16 13 12 1
8
19 18
ζ η ξ
7
17
6 15 14 3
4 9
11 10 2
(3-5)
式中，Ni(ξ, η, ζ)是与按曲线坐标系ξ，η和ζ定义的与每 个节点相联系的形函数， ξ， η和ζ的取值范围为−1 到 1。 20 节点六面体单元形函数在各角点定义为
1
(3-3)
在变换空间中，这些函数写成为
αxξξ − 2 βxξη + γxηη = − J 2 (Pxξ + Qxη ) αyξξ − 2 β yξη + γyηη = − J 2 (Pyξ + Qyη )
(3-4)
这里 J 是变换的 Jacobi 矩阵行列式的值，J = xξxη− yξyη，而
(
)
(
)
(3-7)
(
)
如节点的坐标已知，则与特定坐标ξ，η和ζ相对应的点的坐标可由式(3-5)计算出。 沿每个局部自然坐标，给定划分的数目和划分的间距，这些子区域（超单元）被划分为更小 的单元。以 ξ为例，划分时自然坐标按下式增加
ξ i = ξ i −1 + 2 (Wξ )i WξT
(3-8)
Baidu Nhomakorabea
这里WξT = ∑ (Wξ )i ，ξ0 = −1。 超单元映射法计算简单，通过对自然坐标的不同划分，可以得到密度变化的网格，但超单元 描述的边界和划分区域的实际边界存在误差，不能适应复杂的边界。
3.3.1 映射法
当划分区域的几何特性合适时，四边形或六面体映射网格生成法可以得到相当好的网格[56]。 虽然映射法被认为是结构化网格生成法，但许多非结构网格划分软件都提供这种方法。当应用映 射法时，要求划分区域的对边具有相同的子划分。在三维情况下，将要划分区域的拓扑结构的立 方体的每个相对的面必须有相同的表面网格。这种情况在几何形状复杂是不可能达到，或者需要 人工干预将几何区域划分成可映射的子区域并分配适当的间隔。 Yang 等[57,58]使用的方法是将划分区域分成子区域，每个子区域对应于不同的模块，每个模 块中的网格已预先划分，模块中的网格映射为子区域中的实际网格。这里模块中的网格可以是结 构化的，也可以是非结构化的，但要求相邻子区域的交接面上的表面网格具有相同的拓扑结构。 该方法的特点是可以针对子区域的几何特性和物理特性选择相应的映射模块，生成的网格质量 高，缺点是通用性差。提供映射法生成网格的软件有 ANSYS（ANSYS Inc.）、CUBIT （Sandia National Laboratories ）、goemagic Wrap（Raindrop Geomagic Inc.）、Hypermesh （Altair Computing Inc.）等。
・22・
第 3 章 六面体单元生成方法及相关技术
3.3 非结构网格四边形/六面体网格生成法
由于算法的特点，自动非结构网格生成算法主要是三角形或四面体网格划分，有许多文献讨 论了三角形或四面体网格划分方法和相关的技术[43,44,45,46,47,48,49,50,51,52,53,54]， 许多软件也提供了三角 [55] 形或四面体网格划分功能 。非结构网格生成的特点是能适应复杂的边界，需要的人工干预少， 计算费时。 Steven. J. Owen[55]总结了当前学术界和工业界使用的网格划分方法和可以利用的软件。 非结构六面体网格生成方法主要有映射法、转换法、铺层算法和须段编织算法等。不象三角 形网格划分算法向四面体网格划分算法扩展那样直接，从二维四边形网格划分算法向三维六面体 网格算法的扩展并不是直接的。
・20・
第 3 章 六面体单元生成方法及相关技术
x (ξ , η , ζ ) F (ξ ,η , ζ ) = y (ξ , η , ζ ) z (ξ ,η , ζ )
(0 ≤ ξ ≤ 1,
0 ≤ η ≤ 1, 0 ≤ ζ ≤ 1)
(3-1)
多变量的坐标变换式(3-1)使用无限插值(Transfinite Interpolation)来完成。在三维情况下，使 用混合函数和与之相联系的参数（特定点的位置及偏导数）来显式地决定式(3-1)，然后通过对每 个单变量的循环完成无限插值。一般来说，这些混合函数和参数都选取为区域边界处的函数和参 数。文献[40]给出了几种无限插值函数及其参数的选取。 代数法无限插值网格生成法的特点是计算简单，可以采用中间变量的方法方便地控制网格的 密度，对边界简单的区域，可以生成质量较高的网格，但缺点是不适应复杂边界的划分，边界不 规则时生成的网格的质量很差，并可能产生奇异性。可以通过将划分区域分解为子区域的方法， 在子区域上应用，可以在一定程度上克服这些缺点，但不易实现自动划分。
扫描线
目标
源
图3-4 扫描法生成六面体单元 Fig. 3-4 Hex elements generated by sweeping
3.2.1 代数法网格生成
代数法网格生成技术[40]是通过插值函数将理想的直角坐标系表示的计算区域变换为实际物 理区域来实现的（图3-1） 。它的具体做法是将计算区域的直角坐标用均匀的间隔划分为计算区域 的网格，通过变换，计算区域上均匀分布的坐标被映射为物理区域上的坐标（图 3-2） 。网格点数 目的多寡并不影响变换的特性。从计算区域到物理区域的变换是通过一个向量函数完成的
N i (ξ , η , ζ ) =
图3-3 20 节点六面体单元 Fig. 3-3 20-node hexahedral element
1 (1 + ξξi )(1 +ηηi )(1 + ζζi ), i = 1, 2, L, 8 8
(3-6)
在各棱边的中间节点定义为
1 1 − ξ 2 (1 + ηη i )(1 + ζζi ), i = 9, 11, 17, 19 4 1 N i (ξ , η , ζ ) = (1 + ξξi ) 1 − η 2 (1 + ζζi ), i = 10 , 12, 18, 20 4 1 N i (ξ , η , ζ ) = (1 + ξξi )(1 + ηηi ) 1 − ζ 2 , i = 13, 14, 15, 16 4 N i (ξ ,η , ζ ) =
第3章 六面体单元生成方法及相关技术
3.1 引言
有限元法是求解偏微分方程描述的连续体问题的一种近似工程方法。为克服实际连续体问题 难于处理的问题，它将分析区域离散化，将偏微分方程转化为线性方程组，采用数值计算方法求 出连续体问题的近似解。用有限元法进行工程分析的主要过程包括三个阶段： （1）有限元模型的 建立和数据准备； （2）用软件分析计算； （3）分析结果的判断和评定。迅速而合理地划分有限元 网格是完成有限元分析的前提和保证。 有限元网格生成技术发展到现在，已经出现了大量的、不同的实现方法 [22,23,24,25,26,27,28,29,30,31, 32,33,34,35,36,37,38,39] 。网格生成方法一般可以分为两类：结构化网格生成方法和非结构网格生成方法。 所谓结构化网格，严格的讲，是指网格内部节点具有相同的邻接单元的网格，即网格内部节点具 有相同的拓扑结构。结构化网格生成通常使用复合的迭代平滑技术，用边界或物理区域来排列单 元，划分区域的边界不能太复杂，以便能将区域分解为拓扑结构的块。所谓非结构化网格，放松 了对节点的连接要求，允许一个节点属于任意数量的单元。虽然四边形和六面体单元也可以是非 结构网格，但通常非结构网格是指三角形和四面体网格。目前二维区域网格自动划分发展较为成 熟，许多商业软件都提供了成熟的工具。而三维区域的网格自动划分，除个别商业软件提供了四 面体网格的自动划分外，还不成熟，特别是形状较规范的六面体网格的自动生成。 一般来说，由程序自动划分的网格总有一部分的单元的形状不是很好，需要对网格的质量进 行优化。节点序号的的标注方法直接影响刚度矩阵存储容量的大小从而影响计算的效率，因而需 要对网格的节点重新编号，以得到具有较小带宽的网格系统。 本章主要讨论了三维六面体网格生成的方法及其相关技术。
3.2 结构化网格生成
结构化网格生成方法主要分为两类：代数法和 PDE（Partial Differential Equation，偏微分方 程） 方法。 代数法包括超限插值法、 等参映射法和保角映射法； PDE 方法包括椭圆系统方程、 Poisson 系统方程、抛物线系统方程和双曲线系统方程方法。应用最多的是等参映射法。结构化网格的优 点是易于实现，在每个子区域内网格可以得到很好的控制，生成规则的结构化网格，而且能生成 曲面网格。缺点是对不规则的形体有时生成的网格质量很差，需要事先根据所要产生的网格类型 将目标区域分割为一系列可映射的子区域，这一步通常需要手工完成，因而自动化程度低，不适 合网格全自动生成。提供结构化网格生成的软件有[55]3DMAGGS（NASA，Langley & Lockhead ， 四边形和六面体单元） 、 CAGI （ERC，Miss. State，四边形和六面体单元） 、CSCMDO（CSC/NASA LaRC，四边形单元） 、EMC2（INRIA，四边形单元） 、GENIE++（ERC，Miss. State，四边形和六
・19・
上海交通大学博士学位论文 面体单元） 、ICEM CFD（ICEM CFD Engineering，四边形和六面体单元） 、QulkGrid （Perspective Edge，四边形单元） 、TrueGrid（XYZ Scientific Applications ，Inc.，四边形和六面体单元） 、VGM （NASA，Langley & Lockhead ，四边形和六面体单元）等。