（飞行器设计专业论文）三维混合网格生成与ns方程求解.pdf

摘要 本论文对三维非结构粘性网格的生成方法进行了研究，并在此基础上发展了 一套具有一定通用性的适用于混合网格的流场求解器。为了验证方法和程序的正 确性，论文中对m 6 机翼和f 4 翼身组合体的流场进行了数值模拟，并与试验结 果进行了对比。 在网格生成方面，本文发展了用于n a v i e r s t o k e s ( n s ) 方程数值计算的非 结构粘性网格的自动生成方法。在物面附近的粘性作用区域内使用推进层方法生 成三棱柱或六面体网格，并使用一定的网格生长控制参数实现流场内网格尺寸的 平滑过渡。远离物面的流场其他区域，采用了阵面推进方法生成常规的四面体网 格。 在流场求解中，本文发展了混合网格的求解技术，流场求解程序具有一定的 通用性，并且对于把b a l d w i n l o m a x ( b l ) 湍流模型用于非结构网格做了研究。 在n s 方程的数值计算上，采用了中心差分格式和迎风格式，用格心格式的有 限体积法进行了空间离散，用四步龙格一库塔法作显式时间推进，并采用了当地 时问步长和隐式残差光顺等加速收敛措施。 本文以m 6 机翼和f 4 翼身组合体为算例，对e u l e r 方程和n s 方程的计 算结果进行了比较，从算例上可以看出，n s 方程的计算结果更接近实验值；对 中心差分格式和迎风格式的比较，可以看出，在激波的捕捉方面，迎风格式比中 心差分格式更精确。 关键词： 非结构网格粘性网格n a v i e r - s t o k e s 方程 b a l d w i n l o m a x 湍流模型 迎风格式 a b s t r a c t i nt h i st h e s i s ，as t r a t e g yf o rg e n e r a t i n gt h r e e d i m e n s i o n a lu n s t r u c t u r e dv i s c o u s g r i d s i s i n v e s t i g a t e d a n daf l o wf i e l ds o l v e r b ys o l v i n g t h r e e d i m e n s i o n a l n a v i e r - s t o k e s ( n - s ) e q u a t i o n s o n h y b r i dg r i d si sd e v e l o p e d t h ev i s c o u sf l o wa r o u n d m 6 w i n ga n df 4w i n g - b o d yc o n f i g u r a t i o n a r es i m u l a t e da n dc o m p a r e dw i t ht h e r e s u l t so fw i n dt u n n e le x p e r i m e n t i no r d e rt os o l v en - s e q u a t i o n s ，as t r a t e g y i s p r e s e n t e d f o r g e n e r n i n g t h r e e d i m e n s i o n a l u n s t r u c t u r e dv i s c o u s g r i d sa u t o m a t i c a l l y , w h i c hu s i n g t h e a d v a n c i n g l a y e r sm e t h o d ( a l m ) f o rg e n e r a t i n gp r i s m o rh e x a h e d r a lg r i d si nv i s c o u s d o m i n a t e df l o wr e g i o n sa n du s i n gt h ea d v a n c i n g - f r o n tm e t h o d ( a f m ) f o rg e n e r a t i n g c o n v e n t i o n a lt e t r a h e d r a lg r i d si nt h eo t h e rf i e l do fc o m p u t a t i o n a lr e g i o n i nt h i s p r o c e s s ，t h es t r e t c h e df u n c t i o n i su s e dt om a k e g r i d ss m o o t ha n dg r a d u a l a st h ef l o wf i e l ds o l v e ri sd e v e l o p e df o rh y b r i dg r i d s ，w h i c hs o l v i n gn - s e q u a t i o n s o nu n s t r u c t u r e dg r i d s ，t h eb a l d w i n l o m a x ( b - l ) t u r b u l e n tm o d e li s a d o p t e do nu n s t r u c t u r e dg r i d s i n t h i s p a p e r ，t h eu p w i n ds c h e m ea n dt h e c e n t r a l s c h e m ea r e p r e s e n t e df o rs o l v i n g3 - dn - se q u a t i o n su s i n g t h ec e l l - c e n t e rf i n i t e v o l u m es p a t i a ld i s c r e t i z a t i o na n df o u r s t a g er u n g e k u t t a t i m es t e p p i n gs c h e m e ，w i t h s t a n d a r dc o n v e r g e n c ea c c e l e r a t i o nt e c h n i q u e ss u c h a sl o c a lt i m es t e p p i n ga n di m p l i c i t r e s i d u a ls m o o t h i n g t a k i n gm 6w i n ga n df 4w i n g b o d yc o n f i g u r a t i o na sc o m p u t a t i o n a le x a m p l e s ， s o m er e s u l t sb ys o l v i n ge u l e ra n dn se q u a t i o n sa r ep r e s e n t e d c o m p a r e dw i t ht h e r e s u i t so fw i n dt u n n e le x p e r i m e n t ，t h er e s u l t so f n - s s o l v e ra r em u c hb e r e rt h a nt h o s e o fe u l e rs o l v e r , a n dt h eu p w i n ds c h e m ei sb e a e rt h a nt h ec e n t r a ls c h e m ei nc a t c h i n g s h o c kw a v e k e y w o r d s ： u n s t r u c t u r e dg r i d v i s c o u sg r i d n a v i e r - s t o k e s e q u a t i o n s b a l d w i n l o m a xt u r b u l e n tm o d e lu p w i n d s c h e m e 西北工业大学硕士学位论文 第一章绪论 1 i 计算流体力学的发展现状 计算流体力学( c f d ) 是现代流体力学中的一个重要分支。作为- - 多学科交 叉融合的新兴学科，它是计算机科学、流体力学和计算数学相结合的产物。二十 多年来，随着计算机性能的飞速提高以及计算方法的不断发展，计算流体力学技 术也同趋成熟。与流体力学的实验研究相比，计算流体力学所采用的数值模拟方 法具有周期短、成本低、通用性好的优点。目前，计算流体力学作为一种基本的 研究方法和设计手段，已经逐步进入航空、航天、船舶、气象、能源、建筑等许 多应用研究领域。尤其是在航空、航天领域里的飞行器气动设计中，流体力学的 数值模拟、理论推导和实验研究这三大分支相互补充、相互促进、相互验证，推 动着航空航天技术和流体力学研究的不断深入与发展。 计算流体力学经历了数值求解拉普拉斯方程、小扰动速势方程、全速势方程、 欧拉方程和n s 方程等发展阶段。2 0 世纪8 0 年代以前由于计算机水平的限制，计 算流体力学的数值模拟主要以求解拉普拉斯方程、小扰动速势方程、全速势方程 为，主其中有代表性的是基于拉普拉斯方程的面源法以及有限差分法求解小扰动 速势方程和全速势方程。在随后的二十多年中，由于计算机技术的发展，广大计 算流体力学工作者在求解欧拉方程和n s 方程以及数值模拟复杂流场方面取得 了重大突破。计算流体力学数值模拟的方法主要有：有限差分法、有限体积法和 有限元法。随着诸如t v d 、e n o 年h n n d 等高精度高分辨率差分格式的提出，计算流 体力学对激波、漩涡等复杂问题的模拟能力也有了很大的提高。目前计算流体力 学： 作者正致力于研究和发展更高精度的计算格式和方法，以适应更精细更复杂 的流动研究和设计的需要。 与此同时，计算流体力学的空间划分技术也获得了飞速的发展，结构网格方 面出现了代数生成网格法、解微分方程生成网格法、保角变换法等多种网格生成 方法、网格类型也由单一的c 型网格、o 型网格、h 型网格发展到嵌套网格和多块 对接网格等。非结构网格方面，出现了阵面推进法、d e l a u n a y 方法、八叉树方法 等多种非结构网格生成方法。近年来各种类型的结构非结构杂交网格技术更是 层出不穷，使计算流体力学解决复杂流动和精细流动问题的能力有了很大的提 两北工业大学硕士学位论文 高。 1 2 非结构网格技术的发展 计算网格的生成技术一直是计算流体力学研究的一个重要分支，是计算流 体力学走向工程实用阶段所必须面临的关键技术之一。适合工程使用的网格生成 技术应该具备以下特点： 1 、网格生成过程直观明了、简单易行、效率高、自动化程度好。 2 、通用性、普适性好，对复杂外形、复杂流动的适应能力强。 3 、网格几何灵活性好，尺度变化易于控制，网格自适应加密简便易行。 对于比较简单的几何外形的处理，结构网格生成技术已经发展的比较完善， 人们利用结构网格已经可以获得相当准确和稳定的流场结果。但对于比较复杂的 几何外形，常用的结构网格在处理这些外形时往往会遇到困难。为此计算流体力 学工作者发展了基于结构网格的多块网格分区搭接和嵌套等技术来解决这一问 题。但这些技术仍然有自己的缺点，主要表现是：分区搭接网格为了保证各区块 的网格搭接，常常会造成较大的网格扭曲和变形，使网格的质量难以保证，而且 网格搭接过程中的人工工作量较大；多块嵌套网格虽然可以保证各个区间内的网 格质量，但在重叠区内流场参数的合理传递将成为一项复杂的工作。而且，这两 类网格的生成过程都取决于特定的几何外形，一旦物体的几何拓扑结构发生变 化，整个网格的生成过程也要随之发生较大的变化，因此多块结构网格的搭接和 嵌套很难实现网格生成的通用化和自动化。 非结构网格的提出与应用，成功地克服了结构网格的上述缺点，大大推动了 应用c f d 技术的向前发展。在非结构网格的生成中不再需要考虑网格连接的结构 性和正交性，使的网格单元的大小、形状及网格点的位置更加容易控制，因此非 结构网格具有很大的几何灵活性，对复杂外形具有非常好的普适性，其随机的数 据存储结构使网格的疏密分布和自适应计算变得十分方便，有利于提高计算精 度。f 因为有这一系列的优点，十多年来非结构网格技术受到了计算流体力学界 的广泛重视，并且得到了迅速的发展和应用。目前，国外已经推出了多款具有商 业价值的非结构网格生成软件包；国内也有数家单位开发成功了非结构网格生成 和流场计算软件，并且成功的解决了许多复杂外形区域的流动问题。 4 西北工业大学硕士学位论文 常规的非结构网格在二维问题中以三角形为基本网格单元，在三维问题中 以四面体为基本的网格单元来填充流体的空间区域。常用的非结构网格生成方法 有两类：d e l a u n a y 方法和阵面推进法。d e l a u n a y 方法的基本思路是：首先在空间 布满初始的网格单元，然后按照一定的准则不断引入新点，通过节点重组形成新 的网格单元。这罩的准则通常指的是d e l a u n a y 球面准则，对二维情况该准则要求 每个网格单元的外接圆不包含其它节点；对于三维情况要求四面体单元的外接球 不包含其它网格节点。d e l a u n a y 方法的主要优点是网格生成速度快、网格质量 较高，其缺点是需要事先提供初始网格，网格的疏密变化不易控制对比较复杂的 外形，网格生成时的壁面穿透问题也很难避免。阵面推进方法( a d v a n c i n gf r o n t m e t h o d ) 是一种逐个生成网格的方法，其基本过程是首先将物面离散为一系列初 始推进面元，从中选取某一最合适的面元作为当i j 推进面元，按照一定的准则在 计算空间选取或生成一个节点，这个节点与当前推进面元组合形成网格，同时更 新推进面元系列并重复上述步骤，直到所有流体空间区域都布满网格为止。阵面 推进方法具有网格生成原理简单、自动化程度高、网格疏密容易控制的优点，相 对d e l a u n a y 方法而言，其网格生成的效率较低、网格质量稍差，但通过改进数 据结构和算法以及对生成的网格进行光顺和优化等后置处理，上述阵面推进方法 的缺点是可以克服的。 考虑到结构网格和非结构网格的优缺点，尤其是在进行n s 方程求解时，对 附面层内的流场模拟方面，结构网格有着非结构网格不可比拟的优点。所以，目 前许多研究者开始使用结构网格与非结构网格的混合网格，即在距物面较近的范 围内，生成结构网格，而在流场的其它区域使用非结构网格填充。 1 3 流场数值模拟技术的发展 流场数值模拟的实质就是通过差分或其它办法把流动的控制方程由非线性的 偏微分方程组转化为只包含简单加减乘除等运算的代数方程组，并且将这些代数 方程组和具体的流动边界条件结合起来，编制成程序后让计算机进行求解，从而 获得所需流场参数的过程。 最近二十年里以欧拉方程和n s 方程为控制方程的流场数值模拟技术获得 了巨大的成功，取得了丰硕的应用成果。基于欧拉方程的流场数值模拟技术不但 两北工业大学硕士学位论文 能够模拟和捕捉激波，而且能够有效地模拟有旋流和自动捕获由压力差引起的分 离涡，自动形成分离线和尾迹涡面，揭示许多流动现象的本质特征。相对欧拉方 程而吉，基- 于n - - s 方程的数值模拟技术发展相对缓慢一些。由于计算机容量和速 度的限制完全n s 方程的求解在短时期内还难以进入实用阶段。目前n s 方程的 求解还是建立在雷诺平均n s 方程的基础上，借助适当的湍流模型来实现的。而 湍流问题至今仍然是困扰整个流体力学界的一个难题，湍流的基本机理至今还没 有完全弄清，这就决定了目前的各种湍流模型都有各自的局限。目前尚不存在一 个通用的湍流模型，既能在相当广泛的流动情况中反映出比较正确的物理特征， 从而得到足够准确的计算结果，又能够为当前的计算资源所接受。可以这么说， 湍流模型问题在一定程度上制约了n s 方程数值模拟技术的发展。 相对结构化网格而言，基于非结构网格的n s 方程求解技术还处在起步阶 段，许多结构网格中成熟的流场解算方法尚不能简单的移植到非结构网格之中。 目前国内外许多计算流体力学工作者都在进行这方面的研究，以图使基于非结构 网格的n s 方程求解技术完善化、实用化。 当前，数值模拟欧拉方程和n s 方程的空间离散方法主要有三种：有限差 分法、有限体积法和有限元方法，在非结构网格中常用的是后面两种。有限体积 法的主要思想是把流体力学的积分型控制方程应用到每一个由网格剖分得到的 控制体微元上，从而把积分方程的求解转化为代数方程或常微分方程的求解，再 利用各种时间推进格式进行迭代，最后在每个控制体微元上得到流场的稳定解。 有限元方法的数学基础是变分原理和加权余量法，其主要思想是在每个网格单元 上选择若干合适的点作为解函数的插值节点，方程的近似解由各网格单元的近似 函数逼近，而单元中的近似函数又由单元基函数的线性组合得到。这两种方法各 有所长，有限元方法的数学概念较浓厚，对网格形状的适应性较好，有限体积法 的物理概念十分明确，方法比较简便，容易被接受和掌握。许多研究者已经证明 了分片线性插值近似的情况下的有限元方法和有限体积法是等价的。本文采用的 是有限体积法。 西北工业大学硕士学位论文 1 4 本文的主要工作 本文的主要工作： 1 生成用于n s 方程求解的混合网格。 在三维空间中，由于非结构网格舍去了网格节点的结构性限制，易于控制网 格单元的大小、形状及网格点的位置，因此比结构网格具有更大的灵活性，对复 杂外形的适应能力更强。但是，对于常规非结构四面体网格来说，由于其形状的 限制，使其很难对附面层内的流场进行描述，不适宜作粘性流场的模拟。为此， 发展一套既具有对复杂外形强适应能力，又能够用于n - s 方程求解的粘性网格， 就势在必行了。本文在前人的基础上，发展了一套用于n - s 方程求解的混合网格， h f j 内层足三棱柱网格，外层是常规四面体网格。 2 在混合网格的基础上，编写不局限于网格形状的通用求解器。 由于本文所生成的网格是由三棱柱网格和四面体网格，两种非结构网格组合 而成的，如果在程序的编制上分别处理，不仅增加程序编制的复杂性，而且考虑 到以后进一步发展自适应网格和多块网格技术，原来的求解器不具有通用性，在 使用上要做很大的调整，因此本文发展了一套不局限于网格形状，在某种程度上 具有一定通用性的n s 方程求解器。 3 分别编写二阶中心差分格式和二阶迎风格式的求解器，并就计算结果进行比 较。 中心差分格式是相对而言较为简单的一种差分格式，在程序的编写上较为容 易，同样网格的条件下，其计算时间相对来说也要短一些。但是，其人工粘性项 系数的调节，较为麻烦，并且考虑到本文工作目的是做n s 方程的求解，若人工 粘性项系数给定的不合适，会对计算结果造成较大影响，甚至有可能会覆盖真实 粘性的作用。基于上述原因，本文对迎分格式也作了一些研究。 4 加入b a l d w i n l o m a x ( b l ) 零方程湍流模型，并进行n s 方程的求解。 b l 模型是最简单的一种湍流模型，在工程计算中应用非常广泛，其准确性 和可靠性也得到了较多实验数据的验证，并且它的精度并不比更复杂的两方程模 型差，所以本文选择了b l 零方程模型作为对n s 方程计算的初步尝试。 西北工业大学硕十学位论文 第二章三维混合网格的生成技术 由于非结构网格在处理复杂外形时具有结构网格难以比拟的优点，因此它 在c f d t - t 4 算中得到了广泛而迅速的应用。国内经过十多年的发展和完善以后，在 处理复杂外形的无粘绕流问题上也已经比较成熟。 目前常见的非结构网格生成方法可以分为两类，即d e l a u n a y 三角化法和推进 阵面法。由于d e l a u n a y - 三角化法需要事先提供初始网格，且对于复杂外形容易产 生壁面穿透等缺点，因此本文选择了阵面推进法。 2 1 非结构无粘网格 在三维非结构网格的生成过程中，首先要给定一个边界即给定一个计算域， 并在给定的边界上生成二维的非结构网格即三角形网格作为初始的推进阵面( 二 维的非结构网格本文也是采用阵面推进法生成的) ，然后再用阵面推进法生成空 间的三维非结构网格。 2 1 1 初始推进阵面的形成 初始推进阵面是物面、远场以及对称面等边界上的三角形网格所组成的序 列，他们构成一套封闭的非结构曲面网格。具体生成过程如下： 1 将所求流场的各种边界划分为若干个平面和曲面。 2 在划分好的各个面上采用不同的方法生成表面三角形网格。对于一个平面而 言，可以用常规生二维网格的方法，生成表面的三角形网格；对于曲面，表 面网格的生成较为复杂，方法也很多，本文使用的是投影变换的方法，即先 将曲面投影到二维的平面，在这个平面上生成二维非结构网格，然后再把这 些非结构网格投影回曲面上，就生成了三维曲面上的表面网格。 3 对上面生成的所有网格进行统一编号( 分别对节点和面元进行编号) ，作为 初始的推进阵面。 西i t ：i 业人学硕七学位论文 2 1 2 网格背景参数的确定 在使用阵面推进法生成网格前，一般需要先给定网格生成的背景信息，以 用来控制网格的数目、疏密以及质量等。 确定网格背景参数的方法有很多种，在二维情况下较常见的有背景网格法 和势函数法。但在三维情况下，两种方法实现起来都比较困难，相对于二维来说 效果都不是很理想。本文使用的是通过计算空间点距边界的距离来确定背景参数 的方法，这种方法比较直观，而且实现起来比较容易。其具体过程如下： 1 分别计算物面边界和外边界上每个三角形网格的平均边长，并将其作为该网 格格心处的特征长度。 。 2 定义一个空间参考距离d 1 和一个相应的参考长度s ，d i 群i s 分别表示空间网格 在距物面d 1 的位置上网格的尺度应为s 。 3 网格生成过程中先确定出当前推进面元的面心点距物面的最近距离d ，如果d 大于d 1 ，则当地的背景长度在外边界最近位置的特征长度l 外擐近和参考长度s 间通过某一函数插值得到，如果d d , 于d 1 ，则当地的背景长度在物面最近处 的特征长度三自嘏近和参考长度s 间通过某一函数插值获得。 在生成嘲格的过程中，通过调节空间参考距离d 1 和相应参考长度s 来控制 网格的疏密变化，以得到较为满意的网格。 2 1 3 空间网格的“预长” 为了提高网格生成的效率，合理分布网格节点，保证边界面附近的网格质 量，本文在使用阵面推进法之前，采用了空间网格“预长”的方法。它与阵面推 进法的区别在于：“预长”只是对初始推进阵面序列中的每一个阵面进行一次网 格生成，不像阵面推进法那样对推进阵面序列按某一规则进行排序，按先后顺序 进行网格生成。其过程如下： 1 对初始阵面序列中的每一个三角形，按照网格质量最佳的原则在其前方( 即 该三角行网格的外法线方向) 计算出一个最佳的“预长”点，一般来讲此点 位于通过该网格中心的外法线上，其距面心的距离等于该点所在位置的网格 西j l t 业人学硕：【1 学位论文 背景参数。 2 逐个检查阵面与相对应的“预长”点所组成的四面体，是否与其它阵面相交， 如果不相交，记录此四面体，并更新推进阵面；如果相交，则放弃在此阵面 上的网格“预长”。 3 检查完全部的阵面，则“预长”结束，转入正常的阵面推进。 2 1 4 阵面推进法生成空间网格 阵面推进法生成空间网格的过程如下 1 按照某一规则选取推进面序列中的一个推进面作为“当前推进面”。选取规则 很多，可以自己定义，一般有按最小内角、最小面积以及最短高等方式选取 “当前推进面”的，但不论什么方式，其主要目的就是为了避免较大的网格 单元覆盖较小的网格单元。本文以最小面积作为选取规则。 2 由“当前推进面”的形心所在位置计算出当地的网格背景长度，再由背景长 度和形心坐标计算出一个理想生长点。通常这个点位于通过推进面面心的外 法向线上，其到面心的距离等于网格背景长度。 3 建立候选点序列。从推进面的点序列中找出所有与理想点的距离小于一给定 长度的点，在这些找出的点中再找出位于“当前推进面”外法向一侧的点 然后把这些点作为候选点，建立候选点序列。 4 为了保证网格生成的质量和效率，通常要按某一标准对候选点进行排序。常 见的排序标准是候选点与“当前推进面”所构成的网格质量的好坏。但是这 个标准没有考虑到相交性检查( 见5 ) ，如果候选点的排序不合理，不仅会 降低网格生成的效率，还会引起后续网格生成的困难。本文采用了一种复合 式标准，即不仅考虑所生成网格的质量，还要考虑“当前推进面”与其它推 进面有无共边关系以及候选点与“当前推进面”的距离。 5 相交性检查。逐个检查排序后的候选点与“当前推进面”所构成的四面体单 元的各个面是否与其它推进面以及已生成的网格的各个面相交，直到找到第 一个满足不相交条件的候选点。如果所有候选点都不满足条件，则应该对局 部网格和推进面进行删除和修改，然后重新生成网格。 6 记录新生成的网格并更新推进面。把通过相交性检查的候选点与“当前推进 l o 西北二r ：业大学硕士学位论文 面”构成的网格，以一定的顺序存储起来，并更新推进面序列，加入新的推 进面，同时删除被新网格覆盖的推进面。 7 检查推进面序列中是否还有推进面，如果还有，回到步骤1 。直到推进面序 列为空，网格生成即告结束。 图2 1 阵面推进法生成空间网格的主要流程图 2 1 5 非结构无粘网格生成的后置处理 三维非结构网格生成后，会发现在某些位置的网格质量不太理想，这些质 量较差的网格可能会影响流场求解，因而必须对网格进行光顺和优化处理，以提 高整体网格的质量；同时为了流场求解时的方便，还需要给出网格单元的相邻关 两北工业大学颂i ：学位论文 系以及对边界面的类型进行标定。 本文存储的数据包括所有网格点的坐标，各个网格面是由1 日1 1 ) l 点组成的， 这个i f j ! f 被哪两个网格所共用，以及网格单元由几个面构成，分别是哪几个面等。 对于单个网格而言，好的网格质量主要是指网格的形状很接近正四面体。 对于流场的求解，并不是要求每一个网格都接近正四面体，而是要求流场内的网 格整体质量要高，网格大小的过渡要光滑。为此，给出下列几种质量评估标准： 1 单元体质量系数q 。这个系数是用来评估单个四面体网格质量。其定义为 q = 3 r l r ( 2 1 ) 其中，r 是四面体的外接球半径，r 是四面体的内切球半径。虽然上式是质量 系数的准确定义，但是计算量很大，所以实际应用中一般采用下面的简化定 义 q = 卷 ( 2 z ) 其中，v 是四面体的体积，l 是四面体的最长边的边长。对于越接近正四面体 的网格，其质量系数q 越接近于l ；反之越小。由于网格的体积总是大于零， 所以q 值恒大于零。 2 平均质量系数瓯。此系数用于反映流场内全部网格的总体质量状况。定义为 绋：i 1 nq j ( 2 3 ) v l ，1 即所有网格单元体质量的平均值，q 表示第i 个网格的单元体质量系数n 是流场内的网格总数。 3 关联质量系数q ，。其定义是 上：上上 易n 智q ( 2 4 ) q r 表示第i 个网格的单元体质量系数，n 是流场内的网格总数。此系数是用 来衡量质量较差的网格在全体网格中所占的比率。尤其，当存在一些质量特 别差的网格时，此系数一般不会太高。 本文采用的网格优化方法包括i 。8 p l a c i a n 光顺法和网格变换法。 西北工业大学硕十学位论文 l a p l a c i a n 光顺法，就是将流场内部的网格节点，移动到包含这一节点的所 有网格组成的多面体的形心，一般此时可以得到最佳的平均质量系数和关联质量 系数。但是，实际上这种情况并不总是成立，有时候还会造成网格质量的下降， 甚至会有“负体积”的网格出现。为了解决这个问题，本文采用了一种试探性的 l a p l a c i a n 光顺法，具体步骤如下： i ) 找出与当前网格节点有相邻连接关系的所有节点。 2 ) 计算这些相邻点所形成的多面体的形心坐标。 3 ) 试探性的把当前网格节点移动到多面体的形心位置。 4 ) 验证平均网格质量系数和关联质量系数是否得到提高。 5 ) 检查是否有“负体积”的网格。 6 ) 如果上述条件满足，则移动成功。 7 ) 如果不满足，则利用折半法求出新的试探点( 即上一次的试探点坐标与移 动节点的中点) 。 8 ) 判断试探次数，如果小于设定次数，则回到4 ) ，对网格进行检查，否则 停止试探。 可以看出，即使节点移动失败，网格的质量也不会下降。所以只要移动了， 网格的质量就一定会提高。 网格变换法主要是通过修改空间点之间的连接关系，得到一种网格质量最佳 的组合。对于三维空间，一般是五个点，通过改变点的连接关系，可以得到两种 四面体的组合，一种是两个四面体，一种是三个四面体。 从图2 - 2 可以看出，左图是两个四面体的组合a b c d 和a b c e t 右图三个四面 体的组合a b d e 、a d e c 和c d e b 。当线段d e 和三角形a b c 相交时，这两种组合就 可以互换了。 因此，只要互换的条件成立，我们就可以检查这两种组合，选择一种网格质 量较好的作为我们的网格连接方式。 总的来说，通过上述两种网格的优化方法，无论是网格的平均质量系数还是 关联质量系数都得到了很大的提高。但是需要指出的是，无论你用什么方式优化， 网格质量的改善都是有限的。最终决定网格质量的因索还是阵面推进过程中选点 的好坏，只有好的选点标准才是提高网格质量的关键a 西北工业人学硕士学位论文 2 1 6 小结 图2 2 空间四面体的转化 本节主要介绍了用阵面推进法生成三维空间常规四面体网格的方法。由于篇 幅的限制，对于二维网格以及曲面网格的生成方法没有作详细的介绍，下面给出 一一些附图对其进行一些简单的描述。 附图2 1 和2 2 ，给出了一些用阵面推进法生成的二维平面网格，其基本思 路和三维的是样的，只是三维的表面网格在二维中是用线段表示的。 附图2 3 ，是把m 6 机翼先投影在二维平面，也就把m 6 机翼的四条边框进行投 影，然后再用阵面推进法生成的非结构三角形网格，同一般的平面网格的生成方 法是一样的。 附图2 4 ，是把附图2 3 的网格，在反投影回去得到的真正m 6 机翼的表面网 格，这是一种常见的生成曲面表面网格方法。 附图2 5 ，是某翼身组合体的表面网格，生成方法和前面是一样的。 附图2 6 ，给出了三维空间网格生成过程中的网格图。 2 2 非结构粘性网格 虽然四面体网格有着不少的优点，生成过程也比较简单，但是它在粘性绕流 的模拟方面还不是很理想。尤其是传统的四面体非结构网格，在粘性流场的模拟 方面有其不可克服的缺点，很难用来描述附面层内的流场细节。因此设计一套既 西北丁业大学t b o j ：学位论文 具有传统非结构网格对复杂外形的强适应性，又能用于粘性流场模拟的网格，就 成了问题的关键。 2 0 世纪9 0 年代中期，国外逐步开展三维非结构粘性网格的研究，大体上提 出了两种方案：一类是在物面附近( 粘性主导作用区域内) 采用适合粘性计算的 大“展弦比”扁平四面体网格，在其它区域生成常规意义上的四面体网格，但从 整个流场来看采用的都是统一的四面体非结构网格；另一类是在物面附近( 粘性 主导作用区域内) 采用结构网格或三棱柱型的非结构网格，在流场其它区域采用 常规意义上的四面体网格，整个流场是两种不同类型网格的混合。这两种方法在 具体应用中都获得了成功，都有自己的优缺点。第一类方法的优点在于因为整个 流场是统一的网格类型，所以流场求解程序的编制相对要简单一点，不足的是网 格的数目相对后一种方法要多一点，粘性区域内离散精度的保持要相对困难一 些；而后一类方法的优点在于网格数目较少，粘性区域内的离散精度较高，缺点 在于流场中的网格类型不统一，增加了流场求解程序编制的难度， 考虑到上述两类方法的优缺点，本论文选择了第二类方法，这主要是考虑 到编制的是n - s 方程的求解器，对精度的要求较高，而且n s 方程的求解，相对 于e u i ，e r 方程的求解，更耗时间，在能够达到同等精度的条件下，网格数目要少 则所需要的时间也会短一些。本文生成非结构粘性网格的主要思路是：利用推进 层方法，在物面附近生成适合粘性计算的三棱柱型的非结构网格，在流场的其它 区域则采用阵面推进法生成常规意义上的四面体网格。 2 。2 1 网格生成前的准备工作 推进层方法需要确定推进层上各个节点的生长方向。对于比较复杂的三维外 形，特别是那些位于外形变化较大位置的节点，节点生长方向的确定必须保证为 所有与该点相连的推进层面元的“可视”，即点位于推进层面元的外法线方向， 否则就有可能出现“负体积”。 一个理想的节点生长方向矢量应满足下述条件：该矢量同其对应的节点所连 接的所有三角形面元的夹角都相等。可表示为： 一v i ：c ；i = l ，2 ，(25)d n 1 5 西j t ；t , l k 大学硕十学位论文 其中，v ：是点p 的生长方向矢量，是一个待求的量，i 是与p 点相连的第i 个三角形面元的单位外法向矢量，c 是一个常数，n 是与p 点相连的三角形面元的 个数。v 。的求解，可以用下面的迭代公式： v - k p + 1 ：n w k + l n f + ( 1 一) ； ( 2 6 ) i = l 其中是松弛迭代因子，w 是进行地k + 1 次迭代计算时第i 个三角形面元上 的权重系数，它可由下面的计算公式求得： w ：- - w k ，+ 。? ( 1 一n - - w ，k “) ( 2 7 ) 坚、 口 ( 2 8 ) 其中口，是；，与第i 个三角形面元外法向单位矢量磊，的夹角，五是口。的算术平 均值。所有节点的；，值，可以在开始进行推进之前求出来。n 推进t - - n n n 后，所有节点的；，值要重新计算。 2 2 2 粘性网格的生成 确定好节点的生长方向以后，就可以开始生粘性网格了。具体过程如下： 1 计算当前推进层推进的高度占， 民= a o o + c 1 ( 1 + c 2 ) “) ”1( 2 9 ) 其中瓯是第n 层的推进高度，即当前推进层的高度，瓯是第一层的推进高度， c i 和c 2 是非负的常数，用来控制推进高度逐层变化的速度。一般c 的取值范 围在0 1 到o 2 之间，它控制着推进层内层生长高度的变化速度；c ，的取值范围 在0 0 5 n o 1 5 2 _ 间，它控制着推进层外层生长高度的变化速度。 2 根据推进高度d 和生长方向，预计算出新添节点的位置，然后检查这些节点 是否满足继续推进条件。在几何外形较为简单的情况下，一般推进停止的条 件是当前推进高度占大致等于当地网格生长的背景长度( 在尘粘性网格中并 1 6 西：l b i j k 大学烦士学位论文 不用背景网格的长度，这里的背景长度是指用阵面推进法时计算出来的) ， 但是在复杂外形中，如存在拐角和狭缝，这时还要检查新生长的节点是否很 接近或穿透了对面的推进层面元，以决定是否停止推进。 3 如果通过检查，则记录新生长节点的坐标并对这些节点进行编号和记录；否 则，该点停止推进，将包含该点的三角形面元从推进层中去除，并存储这些 三角形面元作为阵面推进法的阵面。 4 存储三棱柱网格，如果有网格的一个面位于对称面上，则将这个面加入对称 面网格，用于后面的阵面推进法生网格。同时还要更新推进层面元，即用新 生成节点所连接的面元代替旧的面元。 5 重复上述步骤，直到每个节点都满足停止生长的条件。 2 2 3 粘性网格生成后的工作 粘性网格生成结束以后，在进行阵面推进之前，还需要做一些整理工作。 对于对称面，由前面存储的对称面网格找出新的推进阵面序列，在对称面内 用阵面推进法生成二维的对称面网格，并把生成的网格加入阵面推进法的阵面序 列。 阵面推进法的阵面序列还包括推进层法结束时的最外层的面网格，以及远场 的网格。把这些网格组合在一起形成一个封闭的空间作为阵面推进法生成非结构 无粘网格的初始推进面。 2 2 4 小结 从上一节的附图中，我们可以看出来，用阵面推进法生出来的常规的非结 构网格，无论是二维的三角形网格，还是三维的四面体网格，对于粘性流场的模 拟，都是不合适的：在这些附图中，可以看到，贴近物面的第一层网格的高度就 已经很大了，最小的也只能到1 o e 4 这个量级，一般都在1 o e 3 这个量级，以m 6 机翼为例，在给定的实验数据中其雷诺数为1 1 7 e + 6 ，为了使粘性项的贡献不至 于落入截断误差的范围类，其临界网格尺度应在1 o e 4 这个量级，而其第一层网 格的高度则应在1 0 e 一5 这个量级，对于常规的非结构网格，是很难达到这个量级 西北t 业大学硕士学位论文 的。另外，常规的非结构网格的形状也不适合对粘性层内的流场进行模拟。 附图2 7 ，给出了适合于二维粘性流场进行模拟的混合网格的前缘局部图， 附图2 8 ，给出的是后缘局部图。 附图2 9 和2 1 0 ，是用推进层方法进行粘性网格生长前的表面网格上推进 方向的矢量图。 其它一些关于三维粘性网格生成的图，可以见后面的算例。 2 3 网格的整合 当用阵面推进法生成网格以后，所有的网格的生成工作就完成了，在进行 流场的计算之前还需要把两种网格整合在一起，作为计算的输入数据。由于本文 采用的是混合网格，所以数据格式不同于通常所见的数据存储格式。在本文中所 采用的数据格式如下： n n p ( 总的网格点数)n f ( 总的网格面数)n g ( 总的网格数) “n u m 、m l l ” 物面网格数( 表面网格) “n u m f a r ” 远场网格数( 表面网格) “n u m s y m ” 对称面网格数( 表面网格) x 1 y 1z l ( 网格点坐标) x y f jz 1 1 n 1n p l ln p l 2n p l 3 n p l 4n f l ln f l 2 ( 面的信息) ( n l 表示这个面是有几点构成的，n p l l ，n p l 2 ，n p l 3 ，n p l 4 表示是哪几点构成的 这个面，如果这个面是由三点构成的，则使n p l 4 = o ，如果多于4 点可以再增加 n p l 5 等，n f l l ，n f l 2 表示这个面被哪两个网格所共用，如果这个面在边界上，则 使n f l 2 等于一个负值，例如，设置物面上n f l 2 = 一1 。) 两北工业人学硕”i ：学位论文 n nn p n ln p n 2n p n 3n p n 4n f n ln f n 2 n g ln b f lr g f 2n g f 3n g f 4n 6 f 5 n g f 6 ( 网格的信息) ( n g l 表示这个网格是有几个面构成的，n g f l ，n g f 2 ，n g f 3 ，n g f 4 ，n g f 5 ，n g f 6 表示是由哪几个面构成的，如果n g l = 5 ，则使n g f 6 = 0 ，如果多于6 个而可以增 加。) 1 9 西北工业大学硕士学位论史 附录网格生成中所需的几何判据 在网格的生成过程中往往需要用到大量的几何判据，这些几何判据使用的 正确与否直接关系到网格生成的质量、速度以及网格生成的成败。下面介绍一些 常用的几何判据： 1 点与面的关系。用来判断点是位于推进面的外法线方向还是内法线方向，像在 选取候选点时就要用到这一判据。 见图2 3 所示，三角形a b c 代表当前推进面，p 是空间中的一点，三角形 a b c 的外法线方向由公式行= a b a c 计算。 图2 3 点与面的关系图 点与面的关系可由公式s = a p l t 判断，如果f 0 ，则p 点位于三角形的 外法线方向的一侧，反之位于另一侧。 2 线与面的关系。用来判断推进面是否与线段相交，在网格生成的过程中常常需 要判断所生成的网格的各个面是否和其它的面相交，就需要先检查面是否和线相 交。 山图2 4 可见三角形a b c 和线段p q ，h 是推进面的外法线方向。下面分别 计算s 。= 乃；和s ：= 面；，如果f 。和s ：同号，则表示点p 和q 位于三角形 a b c 的同侧，即线段阳与三角形所代表的面不相交，反之位于异侧，即相交。 若线段p q 与三角形所代表的面相交，则交点m 可用如下公式求得： 如= ( 14 - t ) x ，+ t x o = ( 1 + f ) 耳+ f y 。j 2 1 o ) 西北工业大学硕士学位论文 z m = ( 1 + t ) z ，+ 亿口 其中，f ：d ( d + p ) ，d = 否；，p = 一a q ；。 图2 4 线与面的关系图 3 面的交点位置。若线与面相交，则还需要判断交点是位于三角形的内部还是外 部，即线段是否与三角形相交。 如图2 4 所示，分别计算： 晶= ( a b a m ) s 2 = ( b c b m ) n ( 2 1 1 ) e 3 = ( c a c m ) i 1 如果s ，s ：，岛全部都大于零，则表示交点m 位于三角形a b c 的内部，即 三角形a b c 与线段p q 相交，否则不相交。 4 两三角形的关系。由上述一些判据可知，要判断两三角形是否相交，只 要判断其中一个三角形中的各条边是否与另一三角形相交即可，若有一条边与另 一三角形相交，则两三角形相交，若都不相交，则两三角形不相交。 西北工业大学硕士学位论文 第三章基于三维混合网格的n s 方程求解 本文采用格心格式的有限体积法对n s 方程进行空间离散，用四步龙格一 库塔方法作显式的时间推进，采用了当地时间步长、隐式残差光顺等加速收敛措 施。 3 1 控制方程 本文应用的控制方程是三维可压缩非定常n - s 方程的积分形式。在直角坐标 系下，其表达式为： 昙小缈+ 皿地) 勰= j n g ( q ) 勰 ( 3 ) 其中 q = p p u p v p w p g