北航计算流体力学大作业.docx
《北航计算流体力学大作业.docx》由会员分享，可在线阅读，更多相关《北航计算流体力学大作业.docx（23页珍藏版）》请在冰点文库上搜索。
北航计算流体力学大作业
网格生成方法及网格质量控制
（文献综述）
院系：
能源与动力工程学院
姓名：
学号：
指导老师：
宁方飞
一、前言
有限元网格生成是工程科学与计算科学相交叉的一个重要研究领域，在经历了30多年发展后的今天依然十分活跃一方面，有限元法己成为一种能够有效地求解各类工程和科学计算问题的通用数值分析方法：
另一方面，计算机硬件运算能力的不断提高也容许人们对工程和科学计算的规模、复杂度、效率、精度等方面提出更高的要求。
作为有限元走向工程应用的桥梁的有限元网格生成由此获得了源源不断的外在动力。
同时，有限元网格生成算法研究中的某些难点问题始终未能获得真正意义上的解决，它们的研究解决对计算几何与计算数学都具有重要的理论价值。
有限元网格生成方法研究领域己取得许多重要成果，形成了独特的方法论体系，提出了许多有效的算法并研制出一些成功的工程化软件产品。
近10年来，有限元网格生成方法研究不断地深入、完善和发展，各国科研人员不断尝试得到适应性强、应用范围广泛的网格生成方法。
研究重点由二维平面问题转移到三维曲面和三维实体问题，从三角形/四面体网格自动生成转移到四边形/六面体网格自动生成，在并行网格生成、自适应网格生成、贴体坐标网格生成、各向异性网格生成等方面亦取得许多重要进展[1]。
另一方面，不同的网格在有限元计算中表现各异。
网格质量对数值求解效率、收敛性和精度的巨大影响也在逐渐被人们认识到。
因此，网格生成后的质量分析、后处理成为新的研究课题。
尤其针对复杂计算区域，或者不易获得实验数据校核的计算区域，更需要获得高质量的计算网格。
二、网格生成方法
对不规则区域中的流动与传热问题进行数值计算，首先要解决如何进行区域离散化问题。
现在有多种对不规则区域进行离散生成计算网格的方法，统称为网格生成技术。
本章主要详细介绍结构与非结构网格生成技术。
2.1概述
积分区域的网格划分直接影响到方程离散的难易，数值计算的快慢和所需计算机内存的大小，也影响到数值解的收敛性和准确性。
数值计算中所用网格按网格节点排列是否有序，可分为结构化网格和非结构化网格。
在一个区域中，网格的形式可以是单一的，也可以是几种形式的组合。
网格生成的详细分类见图2.1。
对于结构化网格，常用的方法主要有：
正交曲线坐标系中的常规网格、适体坐标法和对角直角坐标法。
由于计算域的不规则性，适体坐标法的三种方法比较常用，下文会做重点介绍。
而对于非结构化网格，常用的方法主要有：
分解和映射法、前沿推进法、Delaunay三角化法和其他方法。
下文也会一一介绍[2]。
图2.1网格生成技术分类
2.2结构化网格生成方法
所谓结构化网格是指，它的网格所有节点的连接关系是规则的，例如，网格中的所有网格点能够通过索引方便的查找和确定。
生成边界固定的网格的最常用方法是需要产生一个连续的，并在边界上完全符合的网格。
也就是把一组连续毗邻的矩形可计算区域通过曲线边界映射到一个实际的物理区域，如图2.2。
图2.2复杂区域网格的映射
本节介绍和分析基于简单区域上的网格生成方法。
应用复变函数法构造的网格光滑性较好，在二维翼型计算等方面曾有过广泛的应用，但由于该方法仅限于解决二维问题，且适用的范围较狭小，已逐渐被新的网格生成方法所取代。
这里主要讨论的是代数法和微分方程法。
2.2.1代数法
这里介绍超限插值法（TFI）。
它是空间映射方法的一种，在这种方法中我们把物理空间上的坐标作为可计算空间坐标的一个函数，首先在区域的边界开始通过插值方法映射边界上的坐标。
由于计算的数据是由区域中无数个点所给出的，所以我们把这种插值算法叫做超限插值法[4]。
但是这种函数映射必须满足两个要求：
（1）这种映射必须是一一映射；
（2）计算区域的边界必须映射到物理区域的边界上，并且还存在。
二维情况下，TFI是把物理坐标x，y看成是计算坐标
，
的函数，通过物理区域边界上的离散点求该区域内部的点。
TFI就是通过两个单向插值的和与它们的积来实现的。
两个单向插值函数分别为：
这两个映射的积映射为：
由布尔和运算定理，最终的TFI函数映射关系为：
公式中，矢量
表示计算区域边界上的点对应的物理区域边界上的点，两种区域边界上的点的对应关系可预先假定。
也就是说，它们的对应关系己知，那么通过公式可以确定计算区域内部任意一点
处的物理区域的对应点的位置矢量
，这样，物理区域的网格便可以计算得到。
图2.3就是采用上面介绍的二维超限插值方法在叶轮机两个直叶片之间所生成的结构化网格示意图。
图2.3采用二维TFI方法在叶轮机两个直叶片之间生成的结构化网格
2.2.2微分方程法
下面介绍的这些网格生成方法都是基于一个思想，根据可计算空间坐标系统，为物理空间坐标系统定义一组偏微分方程，然后在可计算空间网格上求解这些方程组，从而得到在物理空间上的对应网格。
最常用的偏微分方程是椭圆形偏微分方程，这种方程适用于边界封闭的区域（对于边界不封闭的区域，我们会假设一个更大的边界来作为区域的边界进行处理）。
另外还有两种方程就是：
双曲线方程和抛物线方程。
这两种方程主要用于边界不封闭的区域。
（1）椭圆形方程网格生成
这些方法首先对可计算空间坐标
和
定义一个到物理空间坐标x1，和x2的偏微分方程（这样的方法可以推广到三维空间）。
最简单的例子就是拉普拉斯方程：
在二维情况下，极值理论保证了这样的映射关系是一一对应的，拉普拉斯算子是平滑的，所以边界上的溢出导致的不连续性不会传递到区域内部。
为了求解这个形式的方程，我们需要一个在物理区域上的网格，这样方程就可以被转变为其中的xi是依赖变量，言
是独立变量，对于
夸无就有下面的等式关系：
其中，
，等式中的
（度量张量）有这样的关系：
这个准线性方程就可以在矩形区域中，通过把物理边界上点的位置作为可计算空间的边界上点的数据来得以求解。
任何标准的求解方法都可以使用，也可以用代数方法生成的网格作为初始猜测来坚强收敛性。
基于拉普拉斯方程生成的网格不是非常的灵活。
坐标线趋向于等分，除了邻接的凸或者凹形边界，而且这种方法无法控制在边界上网格线的斜度（因为这是一个二阶方程）。
如果要控制网格的密度和斜度，需要加入一个源项（控制函数）。
（2）双曲线方程网格生成
代替椭圆形偏微分方程生成网格的方法，我们还可以用双曲线方程，例如在二维情况下：
这里第一个方程保证了网格的正交，而第二个方程控制网格单元的面积。
这些方法同样也适用于三维的情况。
他们具有一些特殊的性质：
①它们生成的网格都是正交的；②它们不对边界上的不连续点作平滑处理；③它们可以被应用到外部的区域；④网格坐标的生成过程可以作用到边界以外；⑤这一方法在处理凹边界表面时会发生错误（随着这一方法的进行，网格线会重叠，网格单元的体积会变成负值）；⑥这一方法速度快，几乎与一次迭代椭圆形网格生成时间相同。
2.3非结构化网格生成方法
非结构化网格主要是为了有限元方法的应用而发展起来的。
对于有限元方法，可以有很多适用于计算的单元体：
四面体单元，棱柱体单元，块状体单元，它们之间的连接关系是不定的，这种特性也就是非结构化网格的本质。
本节将讨论对于非结构化的网格生成的以下几种方法：
①分解和映射法，这种方法本质上是一种基于多分块区域和TFI的网格生成方法；②前沿推进法（Advancingfront）；③Delaunay三角化方法；④一些不常用的方法，这里是气泡堆积法。
目前最流行的方法是Delauney三角化法和前沿推进法。
2.3.1分解和映射法
这种方法是最早发展起来的方法，而且在一些商业FEM软件包中还在使用。
方法是首先把区域分解为一系列的子区域。
然后每一个单独的区域用映射方法进行网格化。
这种方法可以是无限插值法，但是一种更简便的方法得到了很多的应用，这种方法基于一组图形函数来定义等参的有限元单元，这种等参插值方法是点插值，这里我们给出这种插值在二维空间上应用二次映射（抛物线四边形映射）的方法。
这种情况下，被插值计算的点是在角上的点，和在曲线四边形的边线上的中点（如图2.4），插值公式如下：
其中：
其他的插值公式（如立方体和四次插值等）也可以被使用，而且这种方法可以推广到三维空间。
四边形可以退化为三角形，边界可以看作是C和Q型的网格。
图2.4等参网格
2.3.2前沿推进法
这种方法原理上基于一种很简单的思想（图2.5）对边界进行离散化（例如，在二维空间里用一组多边形来近似），这就是最初的前沿。
加入三角形或四边形到区域内，并且加入的三角形或四边形中至少有一条边与前沿重合。
在每一步中需要更新前沿。
当不再有新的前沿留下时，（例如在二维空间中没有多边形的边留下时）网格的生成也就完成了。
这种方法要求整个区域的边界是封闭的，但是对于边界不封闭的区域，前沿也可以被推进，直到前沿和区域的距离相等为止。
图2.5前沿推进法
在物理空间中，参数化就是有效的曲线弧长。
首先就是计算每个线段的长度（通常CAD软件包可以提供这些信息），然后是插值计算出每个线段上大量样点的网格参数，接着计算出每一条线段上需要的网格边的数目Ns，公式为：
在这里，我们假定没有考虑延展。
网格点的位置在弧长l：
的地方，这里
是一个整型值。
一旦边界被离散化了以后，我们就可以得到初始的前沿，从而可以开始前沿推进法的进行。
2.3.3Delaunay三角化法
Delaunay三角剖分（简称DT）是目前最流行的通用的全自动网格生成方法之一。
DT有两个重要特性：
最大-最小角特性和空外接圆特性。
DT的最大-最小角特性使它在二维情况下自动地避免了生成小内角的长薄单元，因此特别适用于有限元网格生成。
所谓空外接圆特性，就是DT中的每个三角形单元或四面体单元的外接圆（二维）或外接球（三维）都不包含其它节点，Bowyer-Watson算法正是利用了这一特性[6]。
三维Bowyer-Watson算法的基本步骤：
首先定义一个包含所有节点的初始网格，最简单的情形是单个四面体；向网格中插入一个节点，如图2.6，找出其外接球包含此节点的所有四面体单元，删除这些单元形成一个包含插入节点的空腔；将该插入节点与空腔的每个表面相连，形成新的四面体网格；重复进行上述的节点插入过程，直到全部节点插入完毕。
Delaunay三角剖分算法的计算效率与具体的实现相关。
大体上可将DT算法分为三大类：
分治算法、逐点插入法和三角网生长法。
Bowyer-Watson算法是一种典型的逐点插入法，其时间复杂度为O（N3/2），采用四（八）叉树数据结构的Bowyer-Wat-son算法可达到O（NlogN），分治算法的时间复杂度为O（NlogN），三角网生长法的时间复杂度O（N3/2），其中N为生成单元总数。
经典DT技术己经相当成熟，近年来的研究重点是约束DT的边界恢复算法，以及如何克服Bowyer-Watson算法退化现象所产生的薄元问题。
图2.6二维Bowyer-Watson增量插点过程
2.3.4气泡堆积法（bubblepackingmethod）
为了有效地提高网格的质量，基于物理思想发展了一种基于气泡堆积的非结构化网格生成算法。
该算法不同于其他典型的非结构化网格生成算法，其主要思想可以简要概括为3步：
①按一定的规则向给定区域内添加具有虚拟质量的气泡，使其紧密堆积以充满整个计算区域；②引人气泡之间的相互作用力得到各气泡的运动控制方程，通过求解该方程来调整各气泡的位置，反复迭代使所有气泡受力最小，最终达到整个系统的受力平衡，从而有效避免相邻气泡的重叠或分离，实现调整节点位置优化网格的目的；③通过Delaunay角形化方法有效地联接这些气泡的中心，生成非结构化二角形网格。
对于二维情形，其实施流程如图2.7所示。
图2.7BPM算法的流程图
如上所述，初始化完成后的气泡位置不可能非常理想，气泡之间存在部分重叠或间隙。
为了使气泡的位置尽量规则，可为每个气泡定义一个相互作用力，当相邻两个气泡彼此相切时为平衡距离，此时相互作用力为0；当相邻气泡中心之间的距离小于两半径之和时表现为斥力，相反大于两半径之和时为引力。
在力的作用下，系统中每个气泡会自动移动，最终整个系统的受力最小达到平衡，这样可以有效地消除气泡之间的重叠或间隙。
该力可仿照分子间的范德瓦尔力来定义，由于无须关注它的精确形式，Shimada和Bossen都对该力做了近似[5]。
其中，
。
式中：
l0为平衡时的中心距离；l为气泡中心距离；k0为在平衡距离l0处相应的线性弹性系数。
得到每个气泡的相互作用力后，可为每个气泡假定一个虚拟质量m，以通过牛顿第二定律得到其运动方程。
为了加快系统的收敛速度，假定每个气泡在运动过程中还要受到一个阻尼力h（t），简化起见假定该力和气泡运动速度成正比，引人阻尼系数c，得到
由牛顿第二定律得到每个气泡的运动方程为
进一步可得第i个气泡的运动控制方程为
（i=1,2,……,n）
fi为所有与其相邻的气泡作用在第i个气泡上的合力，n为气泡的总个数。
通过求解上述的运动方程，得到每个气泡的移动距离。
在求解运动方程时，为了减少截断误差，采用了四次龙格库塔方法进行求解，并且在求解过程中无须关心气泡轨迹随时间的变化，得到最终的气泡平衡位置即可。
如图2.8所示，在上述算法的基础上，针对矩形区域使用气泡堆积方法与Delaunay划分结合Laplace光顺化方法生成了不同的非结构化网格，并用网格单元规则度p（即三角形内切圆半径与外接圆半径值之比，对正三角此值最大为0.5）和网格几何不规则度（所有网格单元的规则度与最大值0.5之差的绝对值求和再比上总网格单元数）作为指标对两种网格的质量进行了对比。
很明显采用气泡堆积方法生成的三角形网格质量得到了很大的提高。
图2.8矩形区域两种方法对比（网格数2276）
3、网格质量控制
网格质量是影响数值求解效率、收敛性和精度的关键因素。
网格生成后通还需进行后处理以提升网格的质量。
网格质量问题的诱因有很多，如网格生成方法本身、边界约束的限制以及不合理的匹配等。
针对复杂计算区域获得高质量的计算网格是我们目标。
3.1网格类型对模拟结果准确性的影响
哈尔滨工业大学的梅毅、曲建俊、许明伟基于非定常Navier-Stokes方程和滑移网格技术，分别采用三角形和四边形两种网格单元类型，模拟了10m/s风速下垂直轴风力机处于最优尖速比时的风轮流场，通过比较计算最大功率系数值和风洞实验数据，分析风轮旋转过程中两种模拟方法得到压力系数分布，研究了网格单元类型对垂直轴风力机气动性能数值计算精度的影响。
风力机回转直径2.5m，风轮转轴直径0.1m，叶片高3m，弦长0.4m，采用NACA0015翼型。
选用10m/s下风洞实验获得的最大功率系数值与CFD计算数据做对比研究。
采用Gambit软件进行几何建模和网格划分。
建立图3.1所示的风轮叶片和转轴二维模型以及计算域ABCD。
风轮回转直径记为Φ，AC=BD=10Φ，AB=CD=20Φ，风轮回转中心距AC为5Φ，距BD为15Φ，风轮旋转区域直径为2Φ。
分别用四边形网格和三角形网格划分计算域，网格总数分别为345206和365548，叶片表面平均壁面函数y+≈1。
为能得到随时间变化的数值解，采用滑移网格处理旋转区域和非旋转区域之间的耦合问题。
图3.1计算域和计算网格
两种网格单元类型计算最大功率系数的误差如表1所示。
由表可知，四边形网格的计算精度比三角形网格高8.26%。
下面通过叶片表面压力系数分布来分析计算精度差异产生的原因。
表1不同网格单元类型的计算误差
网格单元
类型
最大功率系数计算误差
四边形
4.62%
三角形
12.88%
风轮上游区域转角为30o时，如图3.2所示，两种网格单元类型对叶片表面气流分离位置的预测基本一致。
转角为150o时，叶片出现了两处气流分离，两种网格单元类型对分离位置的预测差别较大。
四边形网格模拟得到的气流分离位于叶片弦长30%和81%处，而三角形网格模拟得到的气流分离位于叶片弦长63%和92%处，与四边形网格相比，三角形网格计算的负压面压力峰值更小，导致叶片气动力偏小。
转角为270o时，两种网格单元类型的压力系数分布相差不大，但气流分离前四边形网格计算的压力峰值较大。
图3.2叶片转角分别在30o、150o、270o时表面压力系数分布
结果表明，两种模拟计算值与实验值之间都存在差别，四边形网格计算最大功率系数的精度比三角形网格高8.26%。
与四边形网格相比，采用三角形网格模拟时叶片气流分离的位置更靠近前缘，负压面压力峰值更小。
在垂直轴风力机流场数值模拟中，采用四边形网格可以得到更为精确的结果。
四边形网格计算精度比三角形更高的原因在于，四边形网格每个单元都有两条平行边，叶片处于某些转角时，气流流动方向会和网格的两条边平行，有利于减小求解N-S方程时的离散误差，而三角形网格则不然，气流流动方向在任何时候都与网格的三条边成一定角度，使之与四边形网格相比，会产生更大的离散误差，而这种计算精度差别在流动复杂时更加明显。
H型垂直轴风力机流场在动态失速相对频繁的上游区域更复杂，而风力机获得能量主要在风轮上游，因此三角形网格计算得到的风力机功率系数比四边形网格小，与实验数据相比误差更大。
在垂直轴风力机流场数值模拟中，采用四边形网格可以得到更为精确的结果[7]。
3.2网格质量对模拟结果准确性的影响
北京航空航天大学的马宏伟、张辉、蒋浩康三人对一带叶尖间隙的低速轴流压气机转子三维流场进行了数值模拟，通过与实验结果对比，分析了网格质量、叶尖网格匹配对压气机转子三维流场计算的影响，发现网格正交性和网格疏密过渡会影响二阶计算精度的收敛，研究还显示压升特性线不能准确的反映网格数目对计算结果的影响[11]。
数值研究的对象为一低速大尺寸单转子轴流压气机，图3.3是该压气机的剖面图。
该压气机由17片C4叶型的转子叶片组成，外径1m，轮毂比0.6，转速1200r/min，叶尖间隙1mm，叶中弦长0.180m。
利用某商用CFD软件进行计算：
计算采用非耦合隐式解法求解雷诺平均N-S方程，采用网格中心型有限体积方法，压力修正采用SIM-PLE算法，连续方程、动量方程和能量方程的空间离散格式为二阶迎风格式，湍流模型采用S-A模型，基于转子弦长的雷诺数为750000。
图3.3低速压气机实验台图3.4带间隙转子计算域网格
计算域为转子单通道，如图3.4所示，计算域进口位于前缘上游1.5倍弦长位置，给定总压、总温和进气角边界条件；计算域出口位于尾缘下游1.5倍弦长位置，给定轮毂径向位置的静压边界条件，出口静压沿径向分布由简单径向平衡方程算出，通过改变出口边界的背压调整转子的流量状态。
轮毂、机匣以及叶片等固壁上给定绝热无滑移边界条件。
网格正交性和网格疏密过渡是网格质量的主要内容。
图3.5显示了两种不同网格正交性的转子网格，一个正交性较差，一个正交性较好。
图3.6显示了两种不同网格疏密过渡的转子前缘附近网格，一个网格过渡较平缓，一个网格过渡不平缓。
使用不同质量的网格在相同边界条件和计算设置的基础上进行二阶精度数值计算，收敛情况常会出现图3.7所示的差别。
图3.7显示了计算域出日截面流量系数的收敛情况，可以看出使用质量较差的网格经过一定步数迭代后出口流量系数围绕一定值呈现较大幅度波动，而使用质量较好的网格收敛正常。
这显示出二阶计算精度对网格正交性和网格疏密过渡情况的敏感。
图3.5两种不同正交性网格
图3.6不同网格疏密过渡情况
图3.7计算域出口流量系数收敛情况
3.3其他网格问题对模拟结果准确性的影响
其他网格问题包括网格的匹配、划分、数量等各方面的问题。
同样来自上一节的文献，由于叶尖间隙和转子通道分别属于不同的网格块，交界处网格有匹配和不匹配两种情况，图3.8显示了尾缘附近间隙和流道网格沿流向不匹配情况，有些网格生成方法还可能造成径向网格不匹配，这些网格不匹配会给计算结果带来影响。
图3.9显示了尾缘附近间隙和流道网格沿流向相匹配情况。
图3.10显示了间隙与流道不同网格匹配情况对转子总压升计算的影响，可以看出径向和流向网格都不匹配对计算结果影响最大，仅流向网格不匹配对计算结果影响较小，网格完全匹配的情况获得了与实验值较接近的结果[11]。
图3.8网格不匹配情况图3.9网格匹配情况
图3.10不同网格匹配情况对总压升计算的影响图3.11转子某轴向截面间隙内静压分布
图3.11显示了转子某轴向截面使用流向匹配和不匹配网格计算出的间隙内静压分布，可以看出沿流向不匹配的网格在交界处不但导致流场不连续而且间隙内部流动细节也发生很大变化。
造成差异的主要原因是在交界面保证流场参数通量平衡的前提下引人的插值误差，在流动情况比较复杂的区域差异更大，这说明间隙内部流动细节的准确模拟要求间隙与流道网格的完全匹配。
目前大部分涉及到叶尖间隙或处理机匣的数值计算没有考虑到叶片尖部网格匹配的问题。
逐渐增加网格数目是减少网格对计算结果影响的重要步骤。
对带间隙转子三维流场计算采用3种不同疏密程度的网格，具体情况列在表2中。
表2本文使用的计算网格
A
B
C
网格总数
194976
393406
654565
流向
92
125
137
径向
53
65
73
周向
37
45
61
间隙径向
9
13
17
间隙流向
25
35
37
y+
3～8
1～6
1～5
计算时间（单流量点）
P4，1G内存
≈8小时
≈12小时
≈20小时
图3.12显示了在相同边界条件下使用不同网格对总压升和流量系数的计算结果，可以看出随着网格数目的增加，计算结果的差异变得越来越小了，越来越稳定。
图3.12总压升和流量系数计算结果比较
四、总结
一般计算网格包括结构化网格和非结构化网格两大类。
结构化网格更有利于计算机存储数据和加快计算速度。
它的生成速度快，具有较好的正交性，网格质量较好，数据结构简单，从而计算效率高，更容易收敛，精度高。
但是它也有一些缺点：
①只适用于形状规则的图形，对于具有复杂外形的飞行器，构造时需要分块处理，非常耗时，难度大；②在结构网格上很难进行网格的自适应。
而非结构化网格可以适应十分复杂的结构，构造简单。
但是，由于网格单元不规则，网格质量受到影响，精度较低。
现代网格生成技术中包含多种网格生成方法。
结构化网格生成方法中，适应性最好的是适体坐标法下的三种方法（复变函数法、代数法、微分方程法）。
非结构化网格的生成方法中，具有代表性的是分解和映射法、前沿推进法和Delaunay三角化法。
随着复杂工程的需要和计算方法提高，诸如气泡堆积法的新兴方法也不断出现。
这些方法极大地促进了网格适用性的增强和网格质量的提高。
通过算例说明，网格类型和质量对数值计算的求解效率、收敛性和精度具有重大影响：
一般而言，四边形网格的计算精度比三角形的更高；网格正交性和网格疏密过渡的改善会消除异常的收敛情况。
网格匹配情况对计算结果有不同程度的影响：
为了减小计算偏差，内部流动细节的计算需要间隙与流道网格完全匹配。
此外，网格数量也会改变计算结果的精度以及计算成本：
网格数过少，会使精度降低；网格数过多，会使计算成本加大。
五、参考文献
[1]关振群，宋超，顾元宪，隋晓峰.有限元网格生成方法研究的新进展[J].计算机辅助设计与图形学
学报，2003，15
（1）：
1-14
[2]刘晶.结构与非结构网格的生成、转化及应用[D].南京：
南京理工大学，2006
[3]金隽.网格生成算法研究和软件实现[D].上海：
复旦大学，2008
[4]张永华.叶轮机CFD网格生成[D].南京
升级会员 