（热能工程专业论文）基于多边形网格的流动与传热数值模拟程序开发.pdf

基于多边形网格的流动与传热数值模拟程序开发 王振兴( 热能工程) 指导教师：徐明海( 教授) 摘要 随着计算流体力学的发展，对不规则区域内流动换热问题的研究越 来越普遍，传统的结构化网格已经不能适应这一发展趋势，非结构化网 格由此得到了广泛应用。但目前广泛使用的三角形网格单元在计算时经 济性并不好，为解决这一问题，本文以基本的三角形网格单元为基础， 构造了与其对应的多边形网格单元。三角形网格单元向多边形网格单元 的转化可通过连接三角形单元的形心或者外心来实现。这样既避开了非 最简形体网格单元生成的困难，又保留了非结构化网格的灵活性。首先 根据这两种不同的转化方法研究了三角形网格向多边形网格的转化，而 后基于边的数据结构离散控制方程，依据完全压力修正算法进行编程计 算，使用了方腔自然对流、顶盖驱等经典算例对程序进行校核验证。通 过对不同网格系统的计算结果进行比较发现：多边形网格系统要比对应 的三角形网格系统收敛更快，所需计算时间更少，而计算精度并没有明 显变化；外心法转化的多边形网格计算精度要比形心法转化的多边形网 格计算精度高，但由于外心法可能造成节点位于当前单元外部的情况， 其通用性不如形心法。最后，本文使用多边形网格对常见的各种离散格 式的计算精度进行了比较，结果发现，在离散求解界面变量时插值情况 的出现，往往会导致计算精度的降低。 关键词：多边形网格，数值模拟，离散格式，完全压力修正算法 d e v e l o p m e n to f n u m e r i c a ls i m u l a t i o np r o g r a mo ff l u i d f l o wa n dh e a tt r a n s f e rb a s e do nt h ep o l y g o n a lm e s h w a n gz h e n - x i n g ( t h e r m a le n e r g ye n g i n e e r i n g ) d i r e c t e db yp r o f e s s o rx um , n g - h a i a b s t r a c t w i t ht h ed e v e l o p m e n to f t h ec o m p u t a t i o n a lf l u i dd y n a m i c s ，t h es t u d yo f t h ef l o wa n dh e a tt r a n s f e ri na b n o r m i t yd o m a i ni sm o l la n dm o l ep o p u l a r t h et r a d i t i o n a ls t r u c t u r e dg r i d sh a v en o ta d a p t e dt h i sd e v e l o p m e n tt e n d e n c y ， s o ，t h eu n s t r u c t u r e d 鲥d so b t a i n e dw i d e s p r e a da p p l i c a t i o n b u t t h e c o m p u t a t i o n a lc a p a b i l i t yo f t h et r i a n g u l a rm e s h e s ，w h i c h a r cw i d e l yu s e d , a r e n o t9 0 0 d i nt h i sp a p e r , i no r d e rt os o l v et h i sp r o b l e m , t h ec o r r e s p o n d i n g p o l y g o n a lm e s hw a sc o n s t n l c t e db a s e do nt h eb a s i ct r i a n g u l a rm e s h t h e c o n v e r s i o nf r o mt r a n g u l a rm a s ht op o l y g o n a lm e s hc a nb e d o n eb y c o n n e c t i n gt h ec e n t r o i d so rt h ec i r e u m e e n t e r so f t h et r i a n g u l a rm e s h e s i n t h i s w a y , n o to n l yc a l lt h ed i f f i c u l t yo f t h eg e n e r a t i o no f t h en o n - s i m p l e s ts h a p e s g r i du n i tb ea v o i d e d , b u ta l s ot h ef l e x i b i l i t yo ft h eu n s t r u c t u r e dg r i dg a l lh e r e t a i n e d f i r s t l y ，t w od i f f e r e n tm e t h o d sw e r eu s e dt oc o n v 溉t h e 鲥du n i t f r o mt r i a n g u l a rm e s ht op o l y g o n a lm e s hi nt h i sp a p e r t h e nt h eg o v e r n i n g e q u a t i o n w a sd i s e r e t i z e d b yu s i n gt h ee d g e - b a s e dm e t h o d a n dt h e c a l c u l a t i n gp r o g r a mw a sd e v e l o p e da c c o r d i n g t ot h ec o m p l e t ep r e s s u r e c o r r e c t i o na l g o r i t h m t h e n $ o i i l ec l a s s i cp r o b l e m sf o rn u m e r i c a lh e a tw a n s f e r s u c ha sl i d - d r i v e nc a v i t yf l o w , n a t l l r a lc o n v e c t i o nf l o wo fs q u a r ec a v i t ya n d s o0 nw e l ls i m u l a t e dt ot e s t i f yt h ep r o g r m - n f r o mt h ec o m p a r a t i v er e s u l t s w h i c hw g l eg a i n e db yu s i n gd i f f e r e n tg r i ds y s t e m s ，i tw a ss h o w nt h a tt h e p o l y g o n a lm e s hc o n v e r g e df a s t e r , a n dt h ec o m p u t a t i o np r e c i s i o nd i d n th a v e o b v i o u sv a r i a t i o n s t h ep o l y g o n a lm e s hw h i c hw a sc o n v e r t e d 姆t h e c i r c u m c e n t e rm e t h o dh a dh i g h e rc o m p u t a t i o np r e c i s i o nt h a nt h a to ft h e c e n t r o i dm e t h o d b u tt h ec i r c u m c e u t e rm e t h o dm a yc a u s et h es i t u a t i o nt h a t t h ef l e wn o d ei so u to f t h ec u r r e n tu n i t , s oi t su n i v e r s a l i t yw a si n f e r i o rt ot h e e e n t r o i dm e t h o d i nt h ee n d , t h ec o m p u t a t i o np r e c i s i o n so ft h ed i f f e r e n t d i s c r e t i z a t i o ns c h e m e sw e r ec o m p a r e db yu s i n gt h ep o l y g o n a lm e s h i tw a s s h o w nt h a tt h ec o m p u t a t i o np r e c i s i o nw o u l dr e d u c ew h e nt h ei n t e r p o l a t i o n a l g o r i t h mw a gu s e dt oc o m p u t et h ev a l u eo f t h ei n t e r f a c ei nt h ed i c r e t i z a t i o n k e y w o r d s ：p o l y g o n a lm e s h , n u m e r i c a ls i m u l a t i o n , d i s e r e t i z a t i o ns c h e m e ， c o m p l e t ep r e s s u r ec o r r e c t i o na l g o r i t h m 独创性声明 本人声明所呈交的论文是我个人在导师指导下进行的研究工作及取 得的研究成果。尽我所知，除了文中特别加以标注和致谢的地方外，论 文中不包含其他人已经发表或撰写过的研究成果，也不包含为获得中国 石油大学或其它教育机构的学位或证书而使用过的材料。与我一同工作 的同志对本研究所做的任何贡献均已在论文中作了明确的说明并表示了 谢意。 签名：墨遮丞 妒7 年歹月矽日 关于论文使用授权的说明 本人完全了解中国石油大学有关保留、使用学位论文的规定，即： 学校有权保留送交论文的复印件及电子版，允许论文被查阅和借阅；学 校可以公布论文的全部或部分内容，可以采用影印、缩印或其他复制手 段保存论文。 ( 保密论文在解密后应遵守此规定) 学生签名：至筮墨 川年歹月知日 导师签名：缮! 蜩渔妊7 年歹月弦日 中国石油大学( 华东) 硕士论文第l 章前言 1 1 课题背景 第1 章前言 对于计算流体动力学和计算传热学的研究，在某种程度上讲，其实 是对n s 方程求解方法的研究。它们所面对的首要任务就是计算区域的 离散化，即网格的生成。这种对求解区域进行离散以生成网格的方法， 我们称之为网格生成技术。其具体过程是对空间上连续的计算区域进行 剖分，把它划分为许多子区域，并确定每个区域中的节点，为方程求解 做好准备。由此可见，网格生成技术在计算流体力学和计算传热学的研 究中起到了非常重要的作用，而对这一方面的研究也一直未曾停止。从 1 9 8 6 年召开第一届国际计算流体力学网格生成会议开始，该会议每隔 2 3 年召开一次，一直延续至今，由此就可以看出网格生成技术在计算 流体动力学与计算传热学中的地位及这一方面研究的活跃程度 两格生成技术的发展经历了一个很长的过程，先后产生了结构化网 格、非结构化网格等多种网格生成技术l ”。在计算传热学的发展初期， 由于研究的传热问题大都是在规则而简单的区域中发生的，因此，结构 化网格足以满足要求，但是，随着研究的深入，不规则区域的流动换热 问题越来越普遍，结构化网格已经不能适应这一趋势。至此以后，不规 则区域网格生成技术的研究成为了一个重要的研究领域。 最初对不规则区域进行离散的方法是采用梯形边界逼近真实边界。 这种方法曾经被广泛使用，但是由于其缺点同样突出，即计算边界是带 有9 0 。角尖峰的锯齿状粗糙表面，所以后来逐渐被其他方法所取代。 1 9 7 4 年t h o m p s o n 等三人提出了生成适体坐标的方法，由此产生了 适体坐标系。所谓适体坐标系是指采用计算的方法来生成一种坐标系， 使得各坐标轴恰好与被计算物体的边界相适应。这种坐标系的生成过程 可以看成是一种变换，即把物理平面上的不规则区域变换成计算平面上 中国石油大学( 华东) 硕士论文第1 章前言 的规则区域。生成适体坐标的方法主要有代数法，微分方程法等等。但 是，采用该方法往往会使控制方程变得相当繁琐，给求解带来很大困难。 块结构化网格是另外一种重要的网格划分方法，它主要有拼接式网 格和搭接式网格。两者的区别在于块与块之间的交接面是否有重叠。块 结构化网格的主要思想是将求解区域划分为几个小的区域，在每一个区 域都采用常见的结构化网格法进行离散，方程一般也分别求解。这种网 格大大减轻了网格生成的难度。还照顾到了不同地区中空闯和对阅尺度 的不同。但是该方法的问题在于：块与块之间信息的传递成为一个难题， 这需要花费大量的功夫。 综上所述，上文提到的各种结构化网格方法，虽然在一定程度上解 决了不规则区域流动与传热问题的求解，但是在越来越复杂的求解区域 面前，正交曲线坐标系下的结构化网格显然已经无法适应实际的要求， 而非正交曲线坐标系下的结构化网格使控制方程变得十分繁琐，同时对 边界的适应性也不十分强壮。因此，网格生成技术的进一步发展成为必 然。 2 0 世纪8 0 年代，研究学者将非结构化网格从有限元方法中引入到 了有限容积法中。该方法对不规则区域显示出了特别强的适应往，而且 在物理概念上清楚、直观，正是这种结构化网格无法比拟的优势，使得 非结构化网格成为了一种处理复杂区域的有效方法，也由此得到了广泛 的应用。 1 2 国内外研究现状 2 0 世纪8 0 年代，菲结构化网格从有限元法中引入至有限体积法中 使用，由于它对不规则区域良好的适应性，使得该方法在计算流体力学 和计算传热学界得到了迅速的发展。到了2 0 世纪9 0 年代，各商业软件 s t a r - c d 、f l u e n t 等都将原有的基于结构化网格方法推广至非结构化网格 上。我国学者近年来也已经对非结构化网格方法进行了大量研究，但相 2 中国石油大学( 华东) 硕士论文第1 章前言 要有两方面的内容：一是网格生成技术的研究；二是基于网格上的数值 算法，即方程的离散与求解。 1 2 1 网格生成技术 在网格生成技术方面，三角形网格的生成技术已经相当成熟。并且 早就应用于各个商业软件中，g a m b i t 、t g r i d 就是其中的代表。目前，在 计算流体力学的数值模拟领域多采用三角形网格，然而三角形网格系统 单元数约是顶点数的2 倍。因此，从计算效率的角度讲，使用三角形网 格进行计算并不是最经济的。为了克服这一缺点，提出了生成非结构化 四边形网格单元1 3 , 4 1 ，但是四边形网格单元的自动生成只是一些经验方 法，在数学上并没有根本解决，而且它们的剖分效率并不高。近年来， 国外不少研究人员开始尝试从另一个角度解决这一问题，即将三角形网 格转化为任意的多边形网格哆一。 1 2 2 方程离散 与结构化网格相比，应用非结构化网格进行计算的主要困难在于如 何求解描述流体运动的控制方程。这是因为在非结构化网格中，对于控 靠4 方程的离散和压力与速度耦合关系的处理要相对匿难。 在具体的离散过程中，瞬态项和源项的离散方法同结构化网格几乎 完全一致，不存在很大的困难，但对流项和扩散项的离散则面临很多问 题，例如精度和稳定性。在对流项的离散上i 2 1 ，主要有以下几种处理方 案：( 1 ) 一阶迎风，( 2 ) 中心差分，( 3 ) 混合方案，( 4 ) - - 阶迎风，( 5 ) 重构梯度， 此外还有采用最小二乘方格式和把n v d 格式应用于非结构化网格提出 的g a m m a 格式。 对于上述各种格式，传统的一阶迎风格式数值耗散比较严重。已经 逐渐被人们淘汰，而中心差分格式在数值上是不稳定的。最小二乘方格 式计算量较大，而以二阶迎风或者重构梯度为基础的延迟修正理论则是 3 中国石油大学( 华东) 硕士论文第l 章前言 一种有效的方法。该方法把低阶格式作为隐式部分，把高阶格式作为显 式部分，取得了准确的计算结果。 扩散项的离散，相对于对流项要复杂很多。它要需要计算出单元界 面上的各个变量值，而各种不同离散方法的产生也都是基于这一要求而 提出的。l i e n 依据格林公式提出了在界面处构造辅助面积来求界面梯度 的思想川，c h o w 等人提出了采用方向导数离散扩散项的方法 g l ， k o b a y a s h i 提出了采用构造法向导数离散扩散项的方法1 9 1 ，w a n g 则采用 了梯度重构的方法【1 0 1 ，d e m i r d z i c 根据动量插值的思想提出了仿动量插值 算法 1 l j 离散扩散项，m u z a f e r i j a 依据最b - - 乘法思想提出了采用最d x - - 乘法1 1 lj 3 离散扩散项。此外，还有采用坐标变换【14 】的方法来离散扩散项。 目前比较流行的方法是格林函数法和最小二乘法。 在离散的过程中，还必须解决压力速度的失耦问题。为了解决这一 问题，多采用基于压力的求解方法1 l 1 5 1 8 l ，即根据连续性方程导出联系 相邻各节点压力的修正方程。但是，由于界面流速的线性插值无法达到 这一目的，r h i e 和c h o w 提出了动量插值的思想。采用这种方法，在考 虑按线性插值方式得出的界面流速的同时还引入了两项附加项，一项是 相邻单元平均压力梯度作用引起的速度变化，另外一项是在假定压力线 性变化时由两单元形心间的压力梯度引起的速度变化。这样就保证了在 任何情况下相邻两点间的压力差都会出现在界面流速的计算式中，即使 流场中出现了锯齿形的压力分布，也可在界面流速的计算中得到反映， 同时使得不合理分量得到衰减。当迭代计算接近收敛且压力呈线性分布 时，两个附加项的差是很小的，因此不会给界面速度带来什么影响。但 如果压力分布不是线性的，两个附加项间的误差就会给界面速度的计算 带来较大影响，此时多采用s i m p l e m 算法。采用这种思想，很容易得 出界面速度修正值的计算公式，同时得出压力修正方程，这样就可以很 好地解决压力速度失耦问题。 1 9 9 6 年，d a t e 分析了动量插值算法的优缺点后，提出了完全压力修 正算法( c p c ) 1 9 1o 这种方法不需要计算动量插值，所用到的界面速度 4 中国石油大学( 华东) 硕士论文第l 章前言 正算法( c p c ) 0 9 1 。这种方法不需要计算动量插值，所用到的界面速度 都采用算术平均值，这样就解决了动量插值算法中动量方程和质量方程 界面流速计算公式不统一的问题。这种方法是从离散方程的微分化形式 构造压力修正方程，从而推导出求解压力修正值的拉普拉斯方程。通过 求解该方程得到压力修正值再结合单元压力平均值，得到用于更新速度 和压力场的修正压力。一般情况下这个方程被理解为导热方程，可进行 离散，因此将该算法推广至非结构化网格后同样适用。为了能够解决压 力与速度耦合的问题，d a t e 对压力进行了修正，他通过引入各邻点压力 对计算点韵影响来达到解决压力与速度失耦的问题。正是由于完全修正 压力算法的优异性，很多学者对此展开了研究【2 0 呻”。 1 3 研究内容及意义 非结构化网格技术的迅速发展得益于其独特的适应性1 2 2 8 1 ，它的网络 拓扑结构无序复杂，表现出一种不规则、无固定结构的特点。因此，非 结构化网格对复杂外形物体适应性更好，不必受限于结构化的矩形拓扑 结构。但是另一方面，正是由于非结构化网格中节点间的关系不是固定 不变的，所以除了把各单元的信息存储外，还必须把与该单元相邻单元 的编号作为联结关系的信息存储起来。这使得非结构化网格的存储信息 量比较大，故需要较大的内存开销。 根据这些非结构化网格的特点，本文以基本三角形单元网格为基础， 构造与基本单元网格对偶的多边形网格单元为控制容积。这样既可以避 开非最简形体网格单元生成的困难，又可以提高非结构化网格的计算效 率，同时还保留了非结构化网格的灵活性。a z i z 等人研究表明，决定有 限容积法数值计算结果精度的不是单元形状而是网格步长。如果采用多 边形两格迸行计算，那么在同样网格步长的情况下，多边形单元数目将 远远少于三角形单元数目，最终使得所求解的线性代数方程组的阶数大 大降低，由此节省大量计算时间。综上所述，采用多边形网格计算流动 中国石油大学( 华东) 硕士论文第1 章前言 与传热问题是一个非常值得研究的领域。 在多边形网格系统构造成功后，本文采用边的数据结构对控制方程 进行了离散，然后分别使用三角形网格和不同的多边形网格进行了计算， 比较了它们的计算结果；最后又比较了不同离散格式带来的差异。 6 中国石油大学( 华东) 硕士论文第2 章三角形单元向多边形单元的转化 2 1 引言 第2 章三角形单元向多边形单元的转化 前文中介绍了采用多边形网格进行计算的现实意义，为了证实其正 确性，本章首先研究三角形网格单元向多边形网格单元的转化的可行性。 由于基础三角形网格单元的生成已经相当成熟，这里不再赘述基础三角 形网格系统的生成理论，转而利用商业软件生成三角形单元。在详细介 绍三角形网格向多边形网格的转化后，采用具体的网格转化实例来验证 该网格转化方法的可行性。 2 2 网格转化 2 2 1 网格转化概述 目前，三角形单元向多边形单元的转化主要有两种转化方式：将 各三角形单元的形心相连，并将连线的交点作为新多边形单元的顶点【卵： 对三角形单元的每条边做中垂线，各中垂线的交点作为新多边形单元 的顶点嘲。这两种方法生成多边形单元各有利弊，形心相连的多边形单 元容易生成，但是在离散扩散项时往往产生交叉扩散项，对于方程的离 散和最后求解带来很多不便，而且对计算精度也有较大影响。第二种方 法虽然可以避免交叉扩散项的产生，但是在网格生成的时候有可能出现 交点在三角形外部的情况。本章中将分别采用两种方法生成多边形网格 单元。 采用第一种方法转化多边形网格，对于区域的内部来讲，相对简单 明了：只需将各三角形单元的形心相连即可构造出多边形网格单元。然 而，对于临近边界的单元来说，由于边界上没有网格单元，因此多边形 网格并不封闭。本文解决这一问题的方法是：取边界边的中点，与邻边 7 中国石油大学( 华东) 硕士论文第2 章三角形单元向多边形单元的转化 界单元的形心相连，由此构造出封闭的多边形单元。 2 2 2 网格转化流程 本节选用方法1 即三角形单元形心相连的方式来描述三角形网格向 多边形网格转化的详细流程。 读入三角形网格数据 首先由商业软件生成三角形网格，导出网格数据文件，编程读入三 角形网格的顶点和单元的相关信息。 三角形网格单元顶点排序 对三角形单元的顶点进行排序，使其顶点按逆时针顺序排列。进行 该步操作是为了在后面求解界面面积矢量时能够一致，而不用另加判断。 记录边的总数 以单元为循环标准，从中选出两个顶点组成一条边，然后与已知的 边进行比较，如没有重复则添加，同时记录下边的两个顶点的编号。 记录边两侧的单元 以上文找到的边的总数为循环截止数，从中选出一条边，即可确定 边的两个顶点，然后对每个单元循环，可依次从单元中选出三组顶点( 每 组两个) ，与边的顶点进行比较，如果重复就记录下这个单元。 添加边界上的虚拟单元 添加虚拟单元的原则是在边界边上多加一点，由此可以构造出一个 面积为0 的虚拟三角形，同时记录下新加虚拟单元的编号及该单元的顶 点编号。这里添加的点选用边界边的中点。 计算单元形心坐标 采用近似的方法计算三角形单元的形心坐标，即求该单元三个顶点 坐标的平均值。此时求得的三角形形心坐标即是多边形单元顶点的坐标。 8 旦互垫盔兰! 兰垄! 堡主笙窒 兰! 童三鱼垄苎歪塑垒望型璺重箜壁些 寻找三角形顶点周围的相邻单元 以上文中记录好的边为基准进行循环，对边的两个点分别进行操作， 边两侧的单元肯定也位于该边端点的周围，但需要注意的是目前记录的 单元有没有跟已有的单元重复，如果重复，则忽略，否则就要记录到当 前顶点的周围单元中。 按逆时针的原则对顶点周围的单元排序 利用中记录好的信息对多边形网格的各个顶点进行排序，使它们 按照逆时针的顺序进行排列。选择相关的点根据右手法则来做具体判断。 输出多边形网格单元的顶点和单元信息 根据上述已经计算好的每个多边形网格单元的顶点数、顶点坐标以 及排列次序输出多边形网格的具体信息。 2 2 3 网格转化实例 1 、正方形区域网格转化 顶盖驱流动是数值传热学界验证算法正确与否的一个经典算例，其 计算区域为正方形，这里采用前面提出的网格转化方法对其进行转化。 ( a ) 三角形网格 ( b ) 多边形网格( 方法1 ) 9 中国石油大学( 华东) 硕士论文第2 章三角形单元向多边形单元的转化 ( c ) 多边形网格( 方法2 ) 图2 1正方形区域网格转化 表2 1 正方形区域网格数据对照 l 三角形单元多边形单元( 方法1 )多边形单元( 方法2 ) 【顶点数 1 3 53 0 83 0 8 f 单元数 2 2 81 3 51 3 5 2 、有凹陷区域网格转化 ( a ) 三角形网格( b ) 多边形网格( 方法1 ) 中国石油大学( 华东) 硕士论文第2 章三角形单元向多边形单元的转化 f c ) 多边形网格( 方法2 ) 图2 - 2 有凹陷区域网格转化 表2 - 2 有凹陷区域网格数据对照 三角形单元多边形单元( 方法i )多边形单元( 方法2 ) 顶点数 3 9 38 8 48 8 4 。+ 单元数 6 8 43 9 53 9 5 3 、换热器区域网格转化 ( a ) 三角形网格( b ) 多边形网格( 方法1 ) 埘 中国石油大学( 华东) 硕士论文第2 章三角形单元向多边形单元的转化 c o ) 多边形网格( 方法2 ) 图2 - 3 换热器区域网格转化 表2 - 3 换热器网格数据对照 l 三角形单元多边形单元( 方法1 )多边形单元( 方法2 ) l 顶点数 1 0 0 92 1 5 0 2 1 5 0 l 单元数 1 8 8 21 0 1 51 0 1 5 4 、区域内部有方形空腔的网格转化 ( a ) 三角形网格( b ) 多边形网格( 方法1 ) 中国石油大学( 华东) 硕士论文第2 章三角形单元向多边形单元的转化 c o ) 多边形网格( 方法2 ) 图2 - 4 内部有方形空腔区域网格转化 表2 - 4 区域内部有方形空腔的网格数据对照 三角形单元多边形单元( 方法1 ) 多边形单元( 方法坌 顶点数 4 3 59 9 09 9 0 单元数 7 5 04 3 94 3 9 一r 5 、区域内部有圆形空腔的网格转化 ( a ) 三角形网格( b ) 多边形网格( 方法1 ) 中国石油大学( 华东) 硕士论文第2 章三角形单元向多边形单元的转化 ( c ) 多边形网格( 方法2 ) 图2 5 内部有圆形空腔区域网格转化 表2 - 5内部有圆形空腔的网格数据对照 l 三角形单元多边形单元( 方法1 )多边形单元( 方法2 ) l 顶点数 5 0 01 1 2 01 1 2 0 l 单元数 8 8 05 0 05 0 0 6 、区域内部有不规则空腔的网格转化 c a ) 三角形网格 1 4 中国石油大学( 华东) 硕士论文第2 章三角形单元向多边形单元的转化 ( b ) 多边形网格( 方法1 ) ( c ) 多边形网格( 方法2 ) 图2 - 6 内部有不规则区域的网格转化 表2 - 6 内部有不规则空腔的网格数据对照 l 三角形单元多边形单元( 方法1 )多边形单元( t y 法2 ) i 顶点数 3 2 2 06 7 2 06 7 2 0 i 单元致 6 1 6 03 2 2 83 2 2 8 l 1 5 中国石油大学( 华东) 硕士论文第2 章三角形单元向多边形单元的转化 7 、喷管区域网格转化 c a ) 三角形网格 ( b ) 多边形网格( 方法1 ) ( c ) 多边形网格( 方法2 ) 图2 7 喷管区域网格转化 表2 - 7 喷管区域网格的数据对照 l 三角形单元多边形单元( 方法1 )多边形单元( 方法2 ) i 顶点数 1 5 0 93 2 8 63 2 8 6 l 单元数 2 7 4 61 5 0 91 5 0 9 8 、不规则区域网格转化 ( a ) 三角形网格 1 6 ( b ) 多边形网格( 方法1 ) 中国石油大学( 华东) 硕士论文第2 章三角形单元向多边形单元的转化 c o ) 多边形网格( 方法2 ) 表2 - 8 不规则区域的网格数据对照a i 三角形单元多边形单元( 方法i )多边形单元( 方法2 ) l 顶点数 1 3 5 l2 8 9 4 2 8 9 4 霉 l 单元数2 5 1 41 3 5 11 3 5 1 图2 - 8 不规则区域网格转化a ( a ) 三角形网格 ( b ) 多边形网格( 方法1 ) 中国石油大学( 华东) 硕士论文第2 章三角形单元向多边形单元的转化 c o ) 多边形网格( 方法2 ) 图2 - 9 不规则区域网格转化b 表2 - 9 不规则区域网格数据对照b i 三角形单元多边形单元( 方法i )多边形单元( 方法2 ) l 顶点数 2 1 76 1 06 1 0 l 单元数 2 5 42 2 92 2 9 ( a ) 三角形网格 ( b ) 多边形网格( 方法i ) 1 8 中国石油大学( 华东) 硕士论文第2 章三角形单元向多边形单元的转化 ( c ) 多边形网格( 方法2 ) 图2 一l o 不规则区域网格转化c 表2 - 1 0 不规则区域网格数据对照c 三角形单元多边形单元( 方法1 )多边形单元( 方法2 ) 顶点数8 5 4 1 9 7 8 1 9 7 8 单元数1 4 3 8 8 7 58 7 5 2 3 小结 由上述网格转化的实例可以看出，无论采用哪种方法对网格进行转 化，最终生成的顶点数和单元数都是相等的。但是，如果将转化后的网 格单元与原三角形单元相比，其单元数大大减少。三角形单元的网格单 元的数日大体上约为多边形网格单元的2 倍。因此，根据多边形网格离 散的方程形成的线性方程组阶数也将大大减少，相应的计算量也会大为 降低。特别是在网格数目较多的情况下，在同样的计算条件下，多边形 网格的计算时间将大大缩短，所以我们认为，使用多边形网格进行数值 模拟具有一定的优越性，可进行进一步的研究。 1 9 中国石油大学( 华东) 硕士论文第3 章控制方程的离散 第3 章控制方程的离散 流动传热问题的数值计算，在对计算区域进行离散后，需要在离散 区域内对控制方程进行求解。本章将简要介绍计算过程中所需的控制方 程”l 。从数学描述上讲，流动传热问题用到的控制方程有质量守恒方程、 动量守恒方程、能量守恒方程。 3 1 描写流动与传热问题的控制方程 3 1 1 质量守恒方程 望+ 型+ 型+ 必：o ( 3 - 0 西咖 砂 勿 上式中的第2 、3 、4 项是质量流密度( 单位时间内通过单位面积的流体 质量) 的散度，可以用矢量符号表示为： o 。p + d i v ( p u ) = 0 ( 3 2 ) a 对于不可压缩流体，其流体密度为常数，连续性方程简化为： d i v 缈) = 0( 3 - 3 ) 3 1 2 动量守恒方程 曼鱼型+ d i vl)：div旧1d甜)+瓯一挈ot ” 一 7 。缸 重!趔+div()：div汹dv)+s，一妻ot ”77 却 旦g 型+ d i v ( 删，) ：d i v 汹d w ) + & 一挈 a t ”77 。龙 其中、鼠、s ，为广义源项，具体如下： s 。= 丢( 玎罢 + 昙( 玎罢) + 鲁( 叩芸 + 丢) ( 3 4 a ) ( 3 - 4 b ) ( 3 - 4 c ) ( 3 5 a ) 中国石油大学( 华东) 硕士论文 第3 章控制方程的离散 鼠= 昙p 考) + 参( 叩考 + 尝( 野考 + 昙汹t v ， o - s ” & = 昙( 玎罢) + 参( 哇) + 未( 可笔 + 丢；v ) c s s c ， 对于粘性为常数的不可压缩流体，瓯= s ，= = 0 ，此时动量方程可以 简化为： 掣+ d i v ( u ) 划v d 甜) 一吉罢 ( 3 6 a ) 掣+ d i v ( v ) = d i v 旧d v ) 一石l 丙a p ( 3 6 b ) 掣砌v ( w ) 础v d w ) 一吉象0 ，- 6 c ) 3 1 3 能量守恒方程 掣+ m v h ) = 一p m v u + d i v ( 2 9 r a e t ) + 饯1 ( 3 - 7 ) 如果取_ j i = c v t ： 掣+ d i 吨, u r ) = 吐毒鲥t 心 p s 的 式中： 品- - s h + 妒 ( 3 s b ) 3 1 4 控制方程的通用形式 在流动与传热问题中，所需求解的主要变量的控制方程都可以表示 成以下通用形式： 掣+ d i v ) = d i v n g r i d 矿) + 邑 ( 3 9 ) 式中尹为通用变量，可以代表甜，v ，w ，r 等求解变量，l 为广义扩散 2 l 中国石油大学( 华东) 硕士论文第3 章控制方程的离散 系数，咒为广义源项a 3 2 控制方程的离散 对万栏( 3 - 9 ) j 珏仃积分，日j 得： n 考争y + v d i v 脚) 州 = 工d i 电v 咄矿+ f s d 矿 ( 3 1 0 ) 根据奥高定理，将体积分转化为面积分 - 挈- - - d v + t 瑚d a j = 扛季d j + 驰。d y ( 3 - 1 1 ) 式中，为通用变量，l 为广义扩散系数，为广义源项，v 为控制容 积的体积，a 为控制容积的表面积，以为控制容积界面的面积矢量，其 正向与外法线单位矢量一致。由此可得： 晦肌姜胪一w ) 鸡= 眵d 矿 ( 3 - 1 2 ) 这里，i 为控制容积界面的角标，n 为控制容积的界面数。 y 图3 - 1 面积矢量 由此可见，要获得关于当前节点变量矛的代数方程，关键是如何将 界面上的对流项与扩散项的积分结果表示成该节点与相邻节点上的关 中国石油大学( 华东) 硕士论文第3 章控制方程的离散 系式。在非结构化网格上不同的离散方法正是根据这一思想产生的。此 外，在离散过程中需要用到面积矢量的定义1 2 】，以图3 - 1 为例： s t = s 。i + s y 。j = i ) ，。一y 、 i + 0 x m + x o j q 1 3 a ) s 。= s 。+ s y ，j = 0 y 。+ y w 、i + b 。一x w 、j o 1 3 b ) 3 2 1 对流项的离散 设界面j 上的对流项为c ，则有 c j = j p 印必；p 蟛) ，以= 办 ( 3 - 1 4 ) 4 式中，下角标j 表示该矢量取值的界面位置，f 为流出j 界面的质量流 量，为界面上庐的平均值，f 则需要相关节点上的速度插值而得0 ，4 - 界面上的平均值，的计算方法取决于所采用的对流项的格式。与 结构化网格中的一些处理方法相对应，在非结构化网格中大致有这么几 种处理方案1 2 1 ：一阶迎风，中心差分，混合方案，二阶迎风， 采用重构梯度( r e c o n s t r u c t i o ng r a d i e n t ) 形成高阶插值，其他高阶格 式：如k 0 b a y 弱l l i 【9 疑出的在非结构化网格上实施的迎风最小二乘方格式 ( u p w i n dl e a s ts q u a r es c h e m e ，u l s s ) ，并且采用局部的一维函数重构的 方法实施了m i n m o d 格式。j a s a k l 2 9 1 将n v d 推广到非结构化网格并提 出了一种g a m m a 格式。 本文采用二阶迎风格式的延迟修正理论求解，具体离散公式如下： 办= 伤葛 式中，为流过当前单元界面的速度( 见图3 - 2 ) 。 。 劈= 如+ ( v 晚b 一乇) 酊= 如+ ( v 比b 一乇) ( 3 1 5 ) ( 3 - 1 价 ( 3 - 1 7 ) 中国石油大学( 华东) 硕士论文第3 章控制方程的离散 y x 图3 - 2 界面上平均值的确定 其中，吃，乃分别表示弓， 昂及j 界面的矢径( j 位于边界 的中点) 。咒点梯度( v k 则根 据只点及其邻点的值来确定。 在实际应用过程中，常采用延迟 修正的方法求解1 站s 】。该方法拔 地阶格式作为隐式部分，把高阶 格式作为显式部分，加入到源项 中去，而系数矩阵的形成依然采用一阶迎风格式。当收敛的时候，一阶 迎风格式的贡献被抵消，只留下二阶中心差分的结果，所以计算具有二 阶精度。这种方法一般和纯粹的一阶迎风格式的收敛速度基本相同。 3 1 2 2 扩散项的离散 在非结构化网格中，扩散项的离散要比结构化网格中复杂得多。对 网格单元做扩散项离散的目的在于，获得能表示该单元受其相邻单元影 响的代数关系式，具体可见图3 3 。 j ，一r p 。、二 c i 图3 3 扩散项离散示意图 中国石油大学( 华东) 硕士论文 第3 章控制方程的离散 假设图3 - 3 中j 界面上的梯度为) ，则它可以分解成为刀，方向( 垂 直界面) 的分量，称之为法向扩散( n o r m a ld i f f u s i o n ) 。及刀，方向( 垂直 于昂、p 的连线) 的扩散，称之为交叉扩散( c r o s s - d i f f u s i o n ) 。为了下 文讨论的方便，把所研究的单元称为昂，其相邻单元记为只( ，= 1 - - 栉) ， 界面记为j 。不同处理方法的区别主要体现在以下两个方面：( 1 ) 如何计 算各单元的v 妒，如果知道了k 及勺b ，则可采用线性插值获得 妒) ，；( 2 ) 如何表示交叉扩散项。 在实际离散过程中，扩散的主要离散方式有两种，分别对应式 ( 3 - 1 8 a ) 和( 3 1 8 b ) 。 i q = 一上鸩= 一r y e , 一( 3 - 。1 s a ) i i d = j 0 v 妒r 刀谢 ，( 3 1 8 b ) a ， 在方法i 中，需要计算单元内的梯度，主要的计算方法有格林函数法和 最小二乘法两种。该式是界面上扩散项引起的总通量，l 为界面上的当 量扩散系数。在j 界面上扩散的总通量可以分解成沿着刀，方向分量d ? 及 沿也方向的分量彤，即： d i = 磷+ 或q 1 9 ) 方法i i 则可以根据数学原理直接进行积分。此外，还有根据式( 3 1 9 ) 直接构造法向扩散和交叉扩散。 根据上述思想，在扩散项的离散上，主要有以下几种方法；依据 格林公式7 1 计算单元梯度 罢= 三抛= 三莓以足( 3 - 2 0 a ) 考= 一去扣= 三军蝇( 3 - 2 0 b ) 式中，a 为当前计算单元的面积，疋，s ，分别是界面上x ，y 方向的面 中国石油大学( 华东) 硕士论文第3 章控制方程的离散 积矢量。 采用最小二乘法的思想1 习计算单元梯度 赢姜南 瞥书k 甜一。 凹， a 勺妒k 鲁乩i 【ir 抿i j 。 一 式中，妒屹表示昂单元梯度在第i 个坐标轴上的分量，办为j a p o 到的有向线段，n 为与昂单元相邻的单元个数。该式要求所确定的 痧k 的三个分量应是按( v 矿) 昂计算褥出的如之值( 上式大括号内的第 二项) 与按线性分布得到的值的差的平方和最小。上式中已知的是当前 层次的如及如之值及各几何要素办，= 1 ，2 ，力，要求取的是勺k 在三个方向的投影，分别由式( 3 2 1 ) 的三个方向( i - l 、2 、3 ) 规定。 由上式所规定的k 的三个分量可用下歹l j 3 x 3 矩阵方程示之： k - - g 1 h( 3 - 2 2 ) 其中，列矢量h 的三个分量及矩阵g 的分量g 。为： 式中，j 为与只单元相邻单元的下标，k 和，则是赢角坐标系中三个 分量的角标。 采用方法i i 对扩散项进行积分，可利用方向导数分解法喁1 进行。如果 记界面两侧单元形心点的连线为f 方向，那么 丝：丝刀 f 3 2 4 ) 砌8 舯棚觯解酣雠继熙等熊用持醐髂一 一2 矽一引 纸一i 衫 一一乃一引瞥并 掣月。互芦 = = 中国石油大学( 华东) 硕士论文第3 章控制方程的离散 函和痧是笛卡尔坐标系中节点矗和尸，之问对应的距离矢量。 在求得单元梯度的情况下，w a n g l l 0 1 提出了一种新的方法来求界 面梯度，即复合梯度法。以图3 - 3 为例，假设已知界面两侧单元的梯度 分别为v 如，v 如，那么界面上的梯度可以近似认为： v 办= 丢勺如+ v 虹) ( 3 - 2 5 ) 结合方向导数可得： 黧+ 黧2 乞 ( v 硝鸭+ 勺硝鸭2 署 式中1 表示蛹方向，皿表示西旒等用警近似代晦、 式中， 表示昂乃方向，皿表示上，方向；竺 用2 学近似代畿、 法向导数法【9 删是模仿结构化网格中界面梯度的计算方法( 见图 3 - 4 ) 。其基本思想为： 型：虹垒( 3 - 2 7 ) 砌 碣2 氟= 石丸+ ( 1 一z 耽( 3 - 2 8 ) 办= 厶如+ ( 1 一五虮( 3 - 2 9 ) 式中，z ， 为界面处的插值系数。在实际应用中，由于在 、办 图3 4 法向导数构造示意图 的构造中用到远邻点信息