（流体机械及工程专业论文）工业送风系统中管道突扩出口均匀气流的优化计算与实验研究.pdf

浙江大学2 0 0 1 届硕士学位论文 摘要 摘要 本文根据工业送风系统存在的问题，采用s i m p l e 方法， 计算模拟了送风系统中不同边界形状的亟遭对管道突扩出 口流场的影响，分析研究了在一定范围内出口流场均匀度随 流道长度的变化曲线，得到了一个合适的流道长度，并且对 流场中有无导流板的情况进行了计算比较，最后通过各种参 数对比得到了一组满意的设计方案。f 在本文的计算中，湍流 模式选择了r n gf 一模型，边界上采用了壁函数技术。本 文随后对这种情况做了实验测试，通过一条生产线上六个骤 冷却单元的数据分析，表明了计算模拟得到的流道设计方案 是可行的，完全达到了工程上的要求；将实验数据与计算模 拟结果相比较发现两者基本吻合但又略有差别，引起这个结 果的主要原因是实验时流道中加入了紊流网和蜂窝板，这两 者的作用使得出口速度有所降低而均匀度变好，这个结果是 完全合理的。卜勺7 一 关键词：管道突扩，优化计算，s i m p l e 方法 塑竺! 竺! ! 笪! ! 兰! 受! 磐g 型! ! 竺! 生堡 一兰坐堡型 a b s t r a c t i nt h i sp a p e r 。“s i m p l e m e t h o dh a sb e e n u s e dt os i m u l a t et h e f l o wi nt h ea i rb l o w n i n gs y s t e mu n d e rv a r i o u sp i p es h a p e st og e t m o r eu n i f o r m v e l o c i t ya tt h eo u t l e t a n a l y z i n gt h ev a r i e t yo ft h e u n i f o r m i t yo ft h ev e l o c i t yo no u t l e t s e c t i o na l o n gw i t ht h ep i p e l e n g t h a n d c o m p a r i n g t h er e s u l tw h e nd i v e r s i o n p l a n e w a s i n s e r t e di n t ot h ef l o wf i e l d ，a no p t i m i z e dd e s i g no fa i rb l o w n i n g p i p e w a so b t a i n e d a c c o r d i n g t ot h eo p t i m i z e dd e s i g nw ed i dt h ee x p e r i m e n tt o v e r i f y o u rc o m p u t a t i o nr e s u l t c o m p a r i n gt h e e x p e r i m e n ta n d c o m p u t a t i o n d a t a ，w e f o u n dt h a tt h e e x p e r i m e n t d a t aw a s a p p r o x i m a t e l y c o n s i s t e n tw i t ht h e c o m p u t a t i o nr e s u l t a n dt h e r e q u i r e m e n t o f i n d u s t r y h a sb e e n e n t i r e l y s a t i s f i e d 浙江大学2 0 0 1 届硕士学位论文 引言 工程背景： 引 在工业生产中，很多地方都用到了骤冷却系统，一般分为水冷却和风 冷却两种方法。水冷却应用比较广泛效果也较好，但是，有很多时候，生 产要求或者生产条件不允许用水进行冷却，因而只能进行风冷却。例如， 化纤、医药、集成电路生产等许多方面，都要求使用风冷却方法。 用风进行冷却一般比较方便，但是对于精度比较高的情况下就出现了 一些问题。譬如，一些医药产品的生产线上，要求送风口出来的冷却风分 布均匀，送风不均匀的地方往往会出现不合格产品，这样一来，使得如何 设计冷却送风系统的问题变得很关键。 本文的课题是从化纤行业中产生的：一种与头发丝粗细相若的玻璃丝 在生产线上要经过骤冷却过程，达到一定的强度要求，这样，很多股玻璃 丝拧在一起不但比单根同样粗细的玻璃丝耐拉耐剪切而且透气性也好。在 骤冷却过程中，用水进行冷却的方法是不可行的，我们选择风冷却的方法。 但是用风进行冷却的过程中，如果冷却风的流速不均匀，那么在风速变化 较大的地方容易使玻璃丝产生粗细不均，甚至使玻璃丝严重变形，这样造 成的最终结果就是玻璃丝的强度不足，达不到厂家的要求。所以，本课题 的任务是设计一个合适的骤冷却系统使得出口流场均匀。 二 问题的提出与解决方法： 在本文的冷却送风系统中，空气在马达的带动下进入回风室，经过骤 冷却单元的表冷器到达管道突扩出口，我们认为通过表冷器后的气流是均 匀的，希望通过管道突扩后的气流仍然是均匀的，所以本文需要设计和改 造的部分是表冷器后面管道突扩这一部分但我们的任务首先是将工程实 际问题转化为力学问题。 本文管道突扩中的流动在工程实际中属于三维不可压定常湍流流动， 如果不进行任何简化那么文章中的优化计算工作量将是十分巨大的。本文 在第一章中，首先建立物理模型和控制方程，然后在第二章中选择合适的 计算方法，选择湍流模式并推导和离散流场的控制方程，接着第三章 对方程使用s i m p l e 算法进行计算，第四章对改造的送风装置进行实 验测试，将实验数据和计算结果相比较，发现不足，进一步改进送风 浙江大学2 0 0 1 届硕士学位论文 引言 装置。 三文献综述： 管道内的二维不可压缩粘性流动，前人已有较多的研究成果，本 文是根据比较成熟的理论和算法解决工业生产中骤冷却系统中的实 际问题，有关骤冷却送风系统这方面的文章还未见到，但是关于流体 力学在工程方面应用的基础书目却有很多【卜1 0j 。在设计流道方面，文 献【7 1 虽然没有讲到对扩张管道的设计方法，但是有关风洞设计方面的 知识，诸如收缩段和实验段的设计原理，蜂窝板和紊流网的作用等等， 对本文都有重要的启示。 计算流体力学发展至今已经有很多成熟的理论和算法【l - 3 1 ，各种新算 法也层出不穷，有关计算流体力学的国内外专著也很多，对不同的情况适 合有不同的计算方法，本文采用s v p a t a n k a r 的s i m p l e 算法【4 1 j 。 工程实际中碰到的流体力学问题大部属于湍流流动，现在关于湍流的 理论和计算研究很多，解决方法各异，但是在工程问题上一般选用较成熟 的工程湍流模式。文献叫5 1 主要是湍流流动基础方面的一些知识，对各种 湍流模型也有细致的讲解和比较，而文献【1 7 2 4 i 则主要是两种湍流模型在不 同算例下的比较。文献【1 8 ) 先后采用了r n gk f 湍流模型和标准k 一5 模型 两种不同的湍流模式，并将计算结果与文献【1 9 】的实测结果进行了比较，发 现用r n gk s 模型来模拟固壁湍流剪切流动更加符合实测数据。对于本 文的不可压缩湍流粘性流动，可由雷诺时均n s 方程和连续性方程来描述 h ”1 ，对于边界层和带回流的流动，目前较多地采用标准k s 二方程湍 流模式1 2 副及其修正模式【2 4 】进行数值模拟。但该模式不能反映湍流的各向异 性，在模拟带有回流的剪切流动时，得到的回流区长度普遍偏短1 2 0 l 。由 y o k h o t 与o r s z a g 2 应用重整化群的方法，导出了一个新型的r n gk g 湍 流模式。对于高雷诺数湍流流动，该模式具有与标准k 一模式相同的形式， 只是将雷诺应力项和耗散生成项系数扩展，即在方程中出现了一个附加的 生成项，当流动快速畸变时，这一项显著增大。该模式也适用低雷诺数区， 可直接积分至壁面，而且它的模式系数均由理论计算得到。本文采用r n g k s 湍流模式来进行计算，模式系数【2 ”2 2j 见第二章。 在工程项目中，仅仅靠计算是没有说服力的，本文在计算之后还进行 了实验工作，在测速方面，根据对各种实验仪器和方法特性2 5 2 8 1 的分析以 及实际情况的需要，本文选择使用了热敏风速仪【26 1 。 堑垩盔兰! 塑! 堡堡主堂垡笙奎 苎二童塑堡堡型垦丝型互堡 第一章物理模型及控制方程 在引言中已经讲了本文的任务，本章主要讲如何建立计算模型，采用 什么样的方法使骤冷却系统的送风口流场变得均匀，以及相关的控制方 程。 1 1 物理模型的建立 简化为二维模型 下面是骤冷却系统的装置图 y 1 图1 1 骤冷却系统装置图 如图1 1 所示，这是骤冷却系统的一个剖面，右图是流道出口的横截 面。由于x 方向是气流的主流方向，流道截面沿y 方向变化，z 方向的截 面如右图所示，本文认为气流速度在z 方向上较小相比x 方向上的流速可 以忽略，又因玻璃丝沿图中的节点线纵向运动，而且生产要求只是对运动 方向上( 纵向) 有限制，横向上速度的不均匀性对玻璃丝的冷却并无太大 影响，所以我们将物理模型简化为二维。 均匀流场的判断标准 我们要得到均匀的出口流场，那么首先必须知道如何判断出口流场是 塑坚查兰! ! ! ! 星堡主堂垡丝奎 笙= 童塑望堡型墨塑型互矍 否均匀，即判断出1 3 流场均匀的标准是什么。 本文在流场出口纵向线上设置了n 个点，将这些点上的x 方向上的速 度值u 进行平均，若平均值与各个点速度值的差值的平方和的平方根 u v a r i a n c e 在1 5 之内，则我们认为出口流场是均匀的。 2 v s u m 至于出口处各点上y 方向的速度值v ，我们将各点上的v 与u 的比值 的平方和的平方根作为一个参数v - s u m ，这个参数越小越好，从后面 的计算和实验我们可以发现v s u m 的值一般比较小，也就是说在大多 数情况下v 相比u 来说，可以忽略。 流道现状分析 要使得出口流场均匀，显而易见，增加流道长度和增大进出口面积比 是最好的办法，也是最根本的方法，但由于生产空间的限制及减少生产成 本的角度，进出口面积已经确定，而且要求我们使用尽量短的流道产生尽 量均匀的流场。 1 2 流道的设计方法 上一节分析了骤冷却系统中现存的问题，这一节里主要介绍我们对流 道设计采用的方法： 流道形状设计： 现在的流道形状不合理，那么，如何设计合适的流道形状呢，怎样判 定所设计流道的性能优劣呢? 设计流道形状的方法有很多，例如复变函数法和保角映射等等，我们 这里采用两种比较简单的方法。 ( a ) 轴向加速度法 这种方法是由b a t c h e l o r s h a w 提出的( 参见 7 ) ，根据理想不可压流 的一维流公式( 认为每个截面上的气流参数是均匀的) ，先假定一个轴向 加速度分布，然后推导而获得响应的计算公式。 流道几何断面如图1 2 ( a ) 所示，进口截面的面积和速度为f i 和v 1 ，出 塑坚茎兰! 塑! 旦堡主堂垡笙奎 兰二童塑堡堡型壁丝鱼! ! ! 笪 口截面为f 2 和v 2 ，扩张段长为l 。假定沿轴向加速度的分布为 d ：ks i n 2 竺 ( 1 1 ) 工 ( a ) 几何特性( b ) 轴向加速度分布 图1 2 扩张段 气流轴向加速度分布如图1 2 ( b ) ，加速度为一个正弦平方的分布。当 x = 0 或x = l 时，a = 0 。即气流在扩张段进i ：1 处的加速度为零，这是符合设 计要求的。k 为系数，其值取决于气流经过扩张段的速度平方增量及扩张 段长度。当x = 0 5 l ，即气流在扩张段中点处加速度达最大值。根据一维流 方程，有 d ：”_ d v( 1 2 ) 因而v d v = ksin2孚出(1-3) 从x = o 至任意x 处积分( 1 3 ) 式 f v d v = r 删n2 詈出 得 v 2 一v i = k ( x - 2 s i n 丁2 7 o c ) ( 1 - 4 ) 当x = l 时，v = v ，因而 k ：! i 二互 ( 1 5 ) 浙江大学2 0 0 1 届硕士学位论文 第一章物理模型及控制方程 考虑到不可压一维流动，因而有 v f = v l e = v 2 ( 1 6 ) 将( 1 5 ) ，( 1 6 ) 式代入( 1 - 4 ) 式，则任意x 处的截面积为 f = f 1 虿鼋再 ( 1 7 ) 这是根据一维流动，且预先假定加速度变化为正弦平方分布，而推导 出来的截面变化公式。若扩张比较小，扩张段又比较长，这个公式还是合 理的。 其实，我们可以从b a t e h e l o r s h a w 的方法中得到提示，可以遵循一定 的原则，得到更多的流道形状，当然，这样未必有很严格的数学证明，但 是我们不妨尝试我们提出的流道形状设计，应用计算机进行模拟，可望获 得多种形式的流道形状。 遵循的原则是：两端气流加速度为零； 气流加速度在中间某一点达到最大： 可以按要求设定加速度最大点的位置； 可以预先设定加速度函数类型。 下面我们将第三章中用到的另外一组曲线边界的方程在这里先行推 导出来，后面就直接应用而不进行具体阐述了。 与上面的原理相同，只是将轴向加速度a 的分布改为两条折线，如图 1 3 所示，这样可以方便的改变a 达到最大值的点k 的位置，我们将k 与 o 点的距离l x 作为一个可变参数，得到一个一般的截面积公式，这样有 助于我们对这种轴向加速度方法的整体认识。 塑堑查堂! ! ! ! 旦堡主堂垡堡苎 苎二童丝堡堡型垦丝型! 三堡 由图1 3 得 图l ，3轴线加速度分布 o x 时：口= l x z ，” l x x _ l 时一壶h 兰 根据一维方程，有 d v 口”忑 即 得 o _ x l x 时：( v 咖= r 毒础 k 时：r v 咖= c 一壶h 兰皿 0 x l x 时： ( 1 8 ) ( 1 9 ) 耶斛时一击( l 2 _ x 2 ) + 兰功( 1 - - 1 4 ) 当x = l 时，( 1 一t 3 ) + ( 1 1 4 ) 得 考虑到一维流动，因而 v f = v 0 = v i e = 1 2 将( 1 15 ) ，( 1 1 6 ) 代入( i 一1 3 ) ，( 1 1 4 ) 中，得 x o k = 2 0 v p 华 i l 七 塑婆奎堂! 塑! 旦堡圭堂垡笙壅 墓二童塑翌堡型墨塑型查堡 0 x l j 时： l z l 时： 若令。l = 1 n ，则( 1 0 x l n 时： l n 石l 时： 当x = l z 时 嘲i f o1 ；x 2 嘲藤 蓐焉 嘲厍) 2 _ 瑟1 1n ( i l - x ) 2 广t 一 7 j _ j 海- 1 1 1 + 弘 ( 1 1 7 ) ( 1 18 ) ( 1 2 0 b ) ( 1 2 1 ) ( b ) 插值函数法 这是一种先用插值函数作为流道形状函数的方法，我们可以先按经验 确定流道边界线所要求的插值条件，然后进行插值，即可获得流道的函数 曲线，然后用计算机进行模拟，观察气流的均匀性效果，再进行改进。 在第三章的算例中，我们会用到的一种曲线边界形状就属于插值函数 法，这里先讲一下制定这种曲线边界的原则。 为简单起见，我们采用多项式插值函数法。由于进出口的大小已经确 定，而我们要达到的要求是出口流场均匀，那么插值条件可以在进出口处 给出： x = 0 ，y = 0 x = ，v = h ( 1 2 2 ) ( 1 2 3 ) 塑翌盔堂! ! ! ! 旦堡主堂焦堡塞 兰二兰塑堡堡型墨塑型! ! 墨 x ；0 ，坐：0 ( 1 2 4 ) 巩 x ：，里：o( 1 2 5 ) 出 由这四个插值条件我们可以定出三次多项式： y = c l y 3 + h 2 + 甜+ d ( 1 2 6 ) 将( 1 2 2 ) ，( 1 2 3 ) ，( 1 2 4 ) ，( 1 2 5 ) 代入( 1 2 6 ) 中得： y ：一娑x3 + j 3 h x 2 ( 1 2 7 ) - 2 y 一可。+ 可。 t 1 其中： 2 五3 厅 ”一可扣可 c = 0，d = 0 另外，本文还用到了折线边界，由于流道比较短，使用曲线边界未必 会有很好得效果，所以我们将这些情况都进行了比较，详见第三章。 增设内部导流板和其他一些方法： 从后面得第三章中我们可以看出，由于进出口面积比太大，所以无论 怎么改变流道的形状，效果均不是很好，流道出口中心位置的速度都很大， 而出口边缘上的速度又较小，且v 向速度过大。所以，要从根本上来解决 这个问题，必须想办法减小截面积比，而整个装置的尺寸基本上已经被限 定，于是本文在流道内部加了三块导流板。这样，相当于将个宽大的流 道变为了一些狭长的流道，从后面得计算及实验结果看，这是解决问题最 有效的方法。 除此之外，为了导直气流，降低气流的紊流度，我们在此骤冷却系统 中加入了蜂窝器、纱网等设施。 蜂窝器由许多方形、圆形或六角形的等截面小管道并列组成，形状如 同蜂窝，故名蜂窝器。蜂窝器的作用在于导直气流，使其平行于流道轴线， 把气流中的大尺度漩涡分割成小漩涡，因而有利于加快漩涡的衰减。同时， 蜂窝管道对气流的摩擦还有利于改善气流的速度分布。 纱网的作用是降低气流的紊流度，故又名紊流网。紊流网的作用与蜂 窝器类似，将较大的漩涡分割成小漩涡，以利于衰减。但紊流网对气流没 塑翌奎兰! ! ! ! 旦堡主兰垡笙奎 篁二童望里堡型垦丝鱼! ! ! 堡 有导直作用。网越细，网的层数越多，整流的效果也就越好，但气流经过 紊流网的损失也就越大。具体的证明和推导可参见 7 】。 1 3 控制方程 连续性方程 d “ = o = o d x 动量方程： 警+ 瓦嚣= 一葶+ 爿( v ) ( 菪+ 期 式中巧平均速度，f 是平均压力，v 与v ，分别为分子粘性系数和涡粘性系数 ，i _汐“d u ”r2 q 了h 2 弘”丢云 搞功苷邑万程： 瓦ok + 虿差= 牡+ 丢吾v 朋r a i k i瓦+ 虿瓦2 只。6 + 瓦j l 吾。a 。i 湍动能耗散率方程： 享+ 巧轰：最( v + 习蠹 + 詈( c t 只一q e ) 式中最是龇项：p 。= 2 v r 粥，峨= 数器+ 矧 上面使用的是r n gk - - * 湍流模型，关于r n gk 一模型的优点及引言中提到 的诵用方棵形式，本立熔存第二童和篦= 童中详细阐诛。 浙江大学2 0 0 1 届硕士学位论文 第二章流场的计算方法 第二章流场的计算方法 计算流体力学是六十年代以来- - n 蓬勃发展的新兴学科，是伴随着高速计算 机的出现而产生的。对流体力学问题传统的研究方法是理论分析和实验研究，但 由于目前对于大量实际的流体动力学问题在数学上仍未解决，而在工程实践方面 很多时候做试验又显得耗资太大，因此计算流体力学就应运而生，它以物理基本 定律和流体力学现有理论为指导，对具体的流体动力学问题作离散化处理，借助 高速电子计算机来进行数值试验或数值模拟，然后把计算得到的结果与已有的对 自然现象的定性了解和部分实验，理论分析结果作比较，来决定数值计算模型的 好坏和数值试验结果的准确性。 计算流体力学常用的方法有有限差分法，有限元方法，边界元方法，有限分 析法等等：有限差分方法发展得较早，使用比较广泛。 本文要进行的是二维不可压缩粘性流动的数值求解。在不可压缩流动问题中 最困难的是质量守恒问题( 即v v = 0 ) ，若流动可压缩，连续性方程中有p ，v 两 项使得求解过程对于误差的敏感性由所降低，但在不可压缩问题数值求解中如果 质量守恒满足不好，则计算误差会很快在计算过程中扩散，造成整个计算失败。 二维不可压缩粘性流动的数值求解方法用的最多的是涡量流函数法、原始变量法 和压力校正法，其立足点都与v v = 0 有关。本文采用了压力校正法。 由于本文是湍流流动问题，因此在求解流场时首先要选择适当的湍流模式。 下面先来讨论这个问题，然后介绍压力修正法s i m p l e 算法。 2 1 湍流模式的选择 湍流是空间上不规则和时间上无序的一种非线性的流体运动，这种运动表现 出非常复杂的流动形态。国内外许多研究者从不同的角度对它的机理进行了研 究，诸如：混沌、分形、重整化群的方法；切变湍流的拟序结构；直接数值模拟； 大涡模拟：概率密度分布函数模拟等。这些湍流机理，有的虽然概念清晰、机理 也较为完整，但距解决工程实际问题为期尚远。总的来说，工程中最常用的方法 仍然是各种湍流模式。 浙江大学2 0 0 1 届硕士学位论文 第二章流场的计算方法 由流场中各瞬时量的n s 方程组出发，忽略密度脉动及输运系数脉动可得 r e n o l d s 平均方程。不可压缩流体的雷诺方程为： 。鲁= 一善+ 杀一毒( 。订) + 瓦 c z 叫 它与n s 方程相比，多出了一个附加项暑( 一历巧) ，雷诺把( 一历巧) 项解释为 作用在流体上的附加应力，称为雷诺应力，它是一个二阶的对称应力张量。这样 方程组就不能自行封闭，因为雷诺应力项一p v j v j ) 的规律是未知的。所谓湍流模 式就是以r e n o l d s 平均守恒方程中的湍流输运项的规律作出公设性的假定，以使 联立方程封闭。假设二阶关联量与平均参数的梯度成比例，变量中具有平均变量， 这种模式称为一阶模式。湍流应力输运方程中，对三阶关联量用模拟公式，变量 包括平均变量梯度和二阶关联量，这种模式即为二阶模式。目前采用的大多是一 阶和二阶模式，湍流模式的问题集中在如何应用模拟的方法求解未知的湍流有效 粘性系数或者各个r e n o l d s 应力分量的问题上。求得湍流输运项而无需增加偏微 分方程的模式称为零方程模式，增加一个参量湍流动能k 的偏微分方程模式为单 方程模式，增加两个湍流参量湍动能k 和湍流耗散率占的偏微分方程的模式则为 双方程模式。 一零方程模式 零方程模式是湍流模式中最早提出的，它至今仍被应用于工程实际问题。这 种近似方法分为两类，一类是由b o u s s i n e s q 提出的湍流粘性理论，另一类为由 p r a n d t l 提出的混合长度理论。 ( 1 ) 湍流粘性理论 l8 7 7 年b o u s s i n e s q 提出了湍流粘性系数的概念，即将r e n o l d s 应力写做如下 形式： 一丽哳( 善+ 割 式中的系数，称做湍流粘性系数。 斗r 。r 了 ( 2 - 2 ) 同理可得出湍流运动粘性系数的概念，即 r 2 - 3 ) 浙江大学2 0 0 1 届硕士学位论文 第二章流场的计算方法 为了求解雷诺方程，还必须对，或者v ，给出新的补充关系式。 ( 2 ) p r a n d t l 的混合长度理论 1 9 2 5 年p r a n d t l 提出了混合长度理论，他通过比较湍流扩散与分子扩散过程， 应用了类似于气体分子动力论中分子自由程的概念，引进了新参数称为混合长 度，它表示这样的一段距离，即脉动微团在经历这段距离内保持有不变的脉动速 度值，这和湍流标尺的物理意义是基本相同的，均表示湍流微团的作用范围，也 可以说是湍流微团的尺寸。针对边界层运动，有如下的关系式： 肾2 倒 c z 叫 只要知道，。的变化规律，方程组就能封闭ap r a n d t l 和t a y l o r 等很多研究者就是 根据这个代换，由实验或假设得出了f 。的规律，对各种工程问题进行了求解。 另外，t a y l o r 根据广义旋涡理论得到的湍流粘性系数是： v r = 三1 l 。2 j d 砂i ( 2 5 ) 可见，t a y l o r 旋涡理论中的混合长度为p r a n d t l 理论中的混合长度的 2 倍，即f ，= 2 z 。 v o nk a r m a n 曾于1 9 8 0 年提出另一种湍流局部相似性假设，他认为在自由湍 流中各空间点的湍流脉动具有几何相似性，即流场中各点的湍流脉动对于同一时 间尺度和长度尺度只有比例系数的差别，因而只用一个时间尺度和长度尺度就能 确定流场的湍流结构。对于二维湍流场v o n k a r m a n 求得的湍流粘性系数为： ( 2 6 ) 式中，k 为卡门常数，k = o 4 0 。 各种零方程模型是用平均速度场的局部速度梯度来表示局部雷诺应力的，使 雷诺方程式中六个雷诺应力分量都变成平均速度场的量，从而达到封闭湍流运动 微分方程组的。这种方法简单，它不增加附加的方程，故称为零方程湍流模型。 浙江大学2 0 0 1 届硕士学位论文第二章流场的计算方法 虽然零方程模型存在不少缺点，但对一些简单的工程实际i 司题，仍是可以用的。 它的最大优点是计算简单。 二单方程模型 p r a n d t l 在1 9 4 5 年建议把湍流粘性系数写成如下的形式： v ，= c k 2 f f 2 - - 7 ) 式中，c 。为实验常数，k 则为湍流脉动动能的时均值，即 = 毒( 可+ 可+ 可) ( 2 - 8 ) ，为湍流特性长度。再加上湍流脉动动能守恒方程： 等= 一l ( 了2 一v ，- 叫- - - ) 百o w k c d 竽+ 苦i v r 瓦o k ) + v 筹( z 一，) 就能使方程组封闭。但是实际上，的值仍须由实验确定，所以这种封闭是不完全 的。 三双方程模型 在单方程的基础上，研究者提出再增加一个方程来使方程完全封闭。湍流双 方程模型中的湍流脉动动能k 方程仍是单方程模型中的七方程，它是在不可压缩 各向同性湍流运动前提下推导出的，适用于高雷诺数的湍流。 。等= 砉c 尝善，+ 文善+ 刮吾+ 筹一。c 。等c z 一。， ! 当引入组合量g = 等后，上式改写成： p 等= 丢 ( 鼍刳+ “，( 毒十割善一。c o e c z 刊， 令“咿= “+ p r ，则上式成为： 。瓦d k = g 吖；, t w 瓦引j + 斗，( 善+ 割吾一p c d e ( 2 - - 1 2 ) 湍流动能耗散率s 的微分方程为： 。告= 毒( 詈剖+ 昙( g 只。s ) c z 邗， 塑望奎兰! ! ! ! 星堡主堂垡堡奎 兰三童堕堑塑盐墨互鲨 式中，a 。为脉动动能耗散翠扩散的p r a n d t l 毅，问 ， 鲁+ 刳吾( 2 - 1 4 ) 为湍流动能生成率。综上所述，经过对湍流运动引用了一s 方程后，湍流运动 微分方程就得到封闭，它们是： 璺：0 卿 。鲁：一暑+ 肛亲一昙( 。啊) + f j 9 百一瓦+ 肛丽i 一面州”j j 。鲁：丢c 等善卅 岳+ 刳善一嘛s ( 2 - - 1 5 ) 。告= 丢( 等刳+ i - ( c t p k 却e ) u t ：c 。p 生 方程组中各常数的取值见下表 q c dc lg哝呸 l0 0 91 01 4 4 1 9 21 o1 3 上面的女一s 双方程湍流模型只适用于高雷诺数的湍流运动。j o n e s 和 l a u n d e r 曾于1 9 7 3 年提出一种低雷诺数湍流运动的k o o 双方程模型，即： 。等：跏+ 讣“，( 刳2 - 9 s - 2 p - 协- 西- - ( 2 1 6 ) 。去= 跏+ - i 一- t _ z _ r 讣c ， “刳2 譬n 吣，( 割( 2 - - 1 7 ) 丽u t = c f 。p _ k z 式中， = 1 0 ， = 1 0 0 3 e x p ( 一r e ) 2 及 浙江大学2 1 3 0 1 届硕士学位论文 第二章流场的计算方法 ( 2 1 9 ) 用k s 双方程湍流模型来计算湍流射流绕钝体流动和湍流边界层类型的 问题，所得到的计算结果和实验值符合得较好。 2 2 通用微分方程 在欧拉法描写的流场中，单位控制体内单位质量的特性量中的输运公式为： 瓦c 3 ( 删+ 毒( p v f 中) = 一夏0 。、j 。) 城( 2 - - 2 0 ) ( 1 )( 2 )( 3 )( 4 ) 式中，( 1 ) 项表示单位体积内特性量的局部变化率，( 2 ) 项表示通过控制面的净通 量，( 3 ) 项表示通过控制面由分子效应7 1 起的输运项厶的散度，( 4 ) 项表示任一内 部和外部过程或源对控制体内特性量变化所作的贡献。根据扩散的f i c k s 定律， 可以把扩散项写成下列形式： 厶= 一r 。g r a d m( 2 - - 2 1 ) 把上式代入式( 2 - - 2 0 ) ，即得： 昙( p 。) + 未( p v m ) = 杀( r 彤) ，+ & ( 2 2 2 ) 上式通用微分方程中的四项分别是不稳定项、对流项、扩散项以及源项。因变量 可以代表各种不同的量( 如p ，口，k ，占) ，只要把不同的因变量和相应的扩散系数 r 。以及源项s 。代入，就可以得到描述各不同物理量的微分方程式。 在二维椭圆型流动情况下，式( 2 - - 2 1 ) 成为： 鲁( 一中) + 知m ) + 扣。) = 昙( k 割+ 参( k 期+ & c z 畅， 通用微分方程中各项的具体含义如下表： 1 4 商 ，l p xe i i ， 浙江大学2 0 0 l 届硕士学位论文第二章流场的计算方法 rs o 连续性方程 10o x 轴动量 “ 劫 ar 抛、 ar 加、 一言+ 瓦l 斗刮+ 丽l p 刮 y 轴动量 v 一蚶 一万+ 瓦l 冽+ 万l 刮 湍流动能七 咿 g k d 盯 湍流动能耗散率 占 咿 c 。只詈一c 2 9 i e ；- 表中： p 曰| 2h + h t 度为 2 3 二维问题的离散化方程推导 在上一节中，我们已经得到了通用微分方程的一般形式 昙( 。) + 知。) + 参( p 仲) = 昙( l 罚+ 号( r 。期战c ：叫， 在如图2 1 所示的控制容积内来推导它的离散化方程形式。定义总流量密 、ll加一缸 + 抛一砂 + 1，j 、 加一砂 + 、ll塑缸 生。批 q j j 兰 蜥 最 浙江大学2 0 0 1 届硕士学位论文 第二章流场的计算方法 则式( 2 2 4 ) 成为如下的形式： 珈) + 警+ 等哉 在控制容积内对方程( 2 2 6 ) 进行积分 ( 2 2 5 ) f 2 2 6 、 业掣小l + 以一纠+ s e c l e ) 缸妙( 2 - 2 7 ) 式中源项已经按常用的方法进行了线性化。对于非稳态的项，假定p ，与中，代表 整个控制容积内的值。“老”的值( 即时间步开始的值) 用p 0 ，司u 0 ，表示。按照全隐 式的做法，所有其它的值( 即不带下标的值) 都应该看成是“新”值。量也、l 、 和以是整个控制容积面上的积分总流量；即也代表整个界面p 上的f 以咖， 依此类推。 类似地，我们可以在整个控制容积内积分连续性方程： 害+ ( 鸭) = o c z 咄， 蛙型a t 竺竺+ e 只+ 一只：。( 2 - - 2 9 ) 式中f 。、f 。、f 。t f ，是通过控制容积面的质量流量。如果在点p 的y u 代表整 个界面e 上的值，我们就可以写出： c = ( p “) 。妙 ：2 燃妙 。：吲， c = ( p v ) 。缸 ” f = ( p v ) 。出 现在以中，乖方程( 2 2 9 ) ，并从方程( 2 2 7 ) 中减去所乖得的结果，我们就 揶一屠挪一彦 民 一 一 中 中 专 三 以 几 浙江大学2 0 0 1 届硕士学位论文第二章流场的计算方法 ( 。，一。；) 丛a t 血y + ( j o c 巾，) 一( j w 只巾，) + ( 1 ，。一，) 一( j 。一只由，) = ( s 。+ s ，) x a y 另外有 j 。一e 中，= d 。( m ，一o 。) j 。一凡巾，= 口，( 巾，一中，) 以一e 中，= a 。( ，一o 。) 以一f 巾，= 口。( 巾。一m ，) 二维的离散化方程成为： a p 中p = d e e 十口矿中矿+ 口中+ 口s o s + b 式中： 。= z ) 。a ( i z l ) + - z ，0 a ，= d ，a ( i p i ) + f w ，0 】 d 。= 乜爿( 1 只1 ) 十 _ e ，0 】 吩= 皿一( + f c ，0 】 b = s c 缸缈+ 口p 0 0 p 口：地 。 ， d p = d e + 口咿+ 口+ d s + 口；一s p z y 相应的传导性定义为 1 1 。妙 也2 两 。 r 。缈 乩。面f n 广。缸 以2 丽 b 2 簖 贝克列数定义为： e 2 每，只= 每，只= 每，只= 每 ( 2 3 2 ) ( 2 3 3 ) r 2 3 4 ) ( 2 3 5 ) ( 2 3 6 ) 7 浙江大学2 0 0 1 届硕士学位论文第二章流场的计算方法 者。 函数爿( 俐) 取值为爿( 例) = 0 ，( 1 - 0 1 1 p i ) 5 ，符号 4 ，b 代表取爿，b 中的大 2 4 流场计算中交错网格的采用 流场的速度场是由动量方程控制的。虽然动量方程对应于上节介绍的通用微 分方程的特殊情况( 取o = “，f = ，等等) ，但真正求解起来远非如此简单。计算 速度场的困难在于未知压力场，压力梯度构成动量方程中源项的一部分。然而还 没有一个用来求解压力场的明显方程，压力场是间接地通过连续性方程规定的。 当正确的压力场代入动量方程时，所得到的速度场满足连续性方程。但是这种间 接地规定只有当我们去求解动量方程与连续性方程所推得的整个离散化方程组 时，才是有用的。而我们并不打算这么做，我们乐意用迭代法来进行求解。 一、计算中的几点困难 1 压力梯度项的表达 对于如图2 2 所示的一维情况， 如果我们着手构成x 方向动量方程的 离散化公式时，那么唯一的新问题是 一咖出项对整个控制容积的积分表 达式。我们假设压力是分段线性分布 的，控制容积面g 和w 选在相应网格点 的中点，那么有： p 。一p 。= 尘二；j 旦上一2 ；上= 兰尘手 x 控制容积 图2 2三网格点群 ( 2 3 7 ) 塑翌盔堂! ! ! ! 旦堡主兰垡堡塞 篁三童堕塑塑堕翌查堂 p _ l q qi 日 10 05 0 01 0 05 0 0 , 1 矿_ c l _ j o ) _ o j k 图2 3锯齿形压力场 这意味着动量方程将包含有两个相间而不是相邻网格点之间的压力差。从图 2 - - 3 可看到由此产生的严重后果，由于各相间网格点上的压力值相等，在任一 网格点p 处，所对应的p ，一p 。= 0 。而动量方程对这样一个波形压力场的“感 受”竟然与对均匀压力场的“感受”一样。在二维情况下也是如此，图2 - - 4 的 压力场是由布置成棋盘形的四个任意压力值组成的，这个高度不均匀的压力场被 动量方程的这种特殊离散化形式处理成为一个均匀的压力场了。由于动量方程对 于这样一种不均匀的压力场毫无反应，所以一旦在迭代过程中形成这样的压力 场，它就会一直保留到解完全收敛为止。 52 752 7 52 7 0 ln n n n1n n n n1n n1 n n 5趣 52 75 2 7 图2 4 棋盘形压力场 2 连续性方程的表达 当我们构成连续性方程的离散化公式时，一个与上述动量方程类似的困难发 生了。对于一维定常的常密度流场，连续性方程简化为： 坐：o 黜 f 2 3 8 ) 如果我们在图2 2 的控制容积内积分这一方程，并且对“采用分段线性假 浙江大学2 0 0 1 届硕士学位论文 第二章流场的计算方法 设和控制容积面取在中点位置，有： “e 一“= 0 ( 2 3 9 ) 于是离散化方程要求相间点而不是相邻点上的速度相等。其结果是如图2 5 所示的那样，完全不符合实际的速度场却满足离散化连续性方程。 u 弛嗡蜗l u l g 吐霉之皿l 图2 5 波形的速度场 二、解决困难的方法一交错网格 把速度分量放在与其它变量不同的网格上会带来相当大的好处。采用交错网 格时，速度分量是对于位于控制 容积表面上的点进行计算的。如图 2 6 所示，“的位置用水平短箭头 表示，而网格点( 以下称为主网格 点) 用小圆圈表示，虚线表示控制 容积面。相对于主网格点而言，“ 的位置只在x 方向是错位的。换句 话说，“的位置落在x 方向连接相 邻主网格点的连线上。v 的位置同 理可得。 图2 - - 6 对u 与v 的交错位置 1 “；叶y ；o 其它的变量 采用交错网格的直接结果是： 无需对有关的速度分量进行任何内插运算，就可以计算出通过控制容积表面的质 量流量。交错网格具有两方面的优点。一方面是离散化的连续性方程将含有相邻 速度分量的差，这样就避免了诸如图2 5 所示的那样波形速度场会满足连续性 方程的情况。另一方面是两相邻网格点之间的压力差现在成了位于这两个网格点 之间的速度分量的自然驱动力，波形压力场不会再被当做均匀的压力场了，也不 会成为最后的解了。 三、动量方程在交错网格下的形式 采用交错网格后，使得离散化的动量方程与在主网格点处计算的其它变量 2 0 浙江大学2 0 0 1 届硕士学位论文第二章流场的计算方法 的离散化方程略有不同。这种不同只是细节性的，而不是本质的。 图2 7 对u 的控制体积图2 - - 8 对v 的控制体积 在如图2 - - 7 ，图2 - - 8 所示的“，v 的控制容积下的动量方程离散化形式为： a c 虬= 。6 “m + 6 + ( p p 一既) 爿。( 2 叫o ) a n = 口。6 v 。6 + 6 + ( p p p n ) 爿。( 2 4 1 ) “，v 的四个相邻点示于控制容积的外边。项6 的定义与方程( 2 3 4 ) 中的定义相同， 但是压力梯度不包括在源项量s c 与品中，压力梯度构成离散化方程的最后一项。 动量方程只有在压力场己知或是按照某种方法估计出来的情况下才能求解。 除非采用正确的压力场，所得的速度场将不会满足连续性方程。用“，v 表示这 样的一个以估计的压力场p 为基础的不完善的速度场： a e “：= 口。“二+ 6 + ( p ；一p ：) 爿。( 2 叫2 ) 。v ：= 口。v 二+ 6 + ( p ；一p ：) 爿。( 2 - - 4 3 ) 2 5 压力修正方程 为了寻找到一个改进压力估计值p + 的方法，以使所算得的带星号的速度场 逐渐地接近于满足连续性方程。让我们假定正确的压力p 由下式得到： p 2 p + p ( 2 叫4 ) 塑坚查堂! ! ! ! 星受主兰垡堡奎 笙三主堕垄堕盐兰查堂 其中p 称之为压力修正。引入相应的速度修正“，v 。，使得 “= “+ “ v = v + v + 如果我们从方程( 2 4 0 ) 减去方程( 2 - - 4 2 ) ，得到 掣如+ ( p j p ；) 爿。 从上式中摒弃项口。“- ，这样做的结果是： 玑= 也( p ；一p ：) 其中d 。；尘 a 。 对v 。同理可得。最后得到： “。= “：+ d 。( p ；一p ：) ( 2 - - 4 8 ) v 。= v ：+ d 。( p j p j ) ( 2 - - 4 9 ) 下面来把连续性方程转换为压力修正方程。在式( 2 2 9 ) 中已得到连续性方 程的离散化形式为： 竺学+ ( 州。缈一( ) ，缈+ ( ) 。缸一( ) ，缸= 。 ( 2 5 。) 把式( 2 4 8 ) 与( 2 4 9 ) 的速度修正公式代入，得到下列对p 的离散化方程： 其中 a ：pp = a e p e + q p + an p n + ns p 。+ b “e = p 。d 。砂 = p 。d 。缈 a = p 。d 。k x 日s = p ，d 。a x( 2 5 2 ) n p = qe + n 七a n + d s 。：堡生! a t 竺竺+ ( 彤) 。一p u ) 。】缈+ ( + ) ，一p v ) 。】缸 由式( 2 5 2 ) 的最后一式可见：压力修正方程中的项b 实际上是按带星号的 速度取值的离散化连续性方程左侧的负值。如果b 值为o ，那就意味着带星号的 速度值满足连续性方程，从而不再有必要对压力进一步的修正。于是项b 代表一 浙江大学2 0 0 1 届硕士学位论文 第二章流场的计算方法 个“质量源”，该质量源必须由压力修正去消除。 2 6 “s i m p l e ”算法和“s i m p l e r ”算法 一“s i m p l e ”算法 “s i m p l e ”算法的全称是s e m i - - i m p l i c i t m e t h o df o rp r e s s u r e - - l i n k e d e q u a t i o n s ，意思是解压力耦合方程的半隐式法。 主要的计算步骤是： ( 1 ) 估计压力场p + 。 ( 2 ) 求解动量方程( 如式2 4 2 与2 4 3 ) 以得到“，v + 。 ( 3 ) 解p 方程。 ( 4 ) 由方程( 2 4 4 ) ，以p 加p + 得到p 。 ( 5 ) 利用速度修正公式( 2 4 8 ) 与( 2 4 9 ) 由带星号的速度值计算“，v ( 6 ) 求解那些通过源项、流体物性等等影响流场的其它一些物理量o ( 如湍动 能，湍动能耗散率) 的离散化方程。 ( 7 ) 把经过修正的压力p 处理成一个新的估计的压力p ，返回到第二步，重 复全部过程，直到求得收敛的解为止。 二“s i m p l e r ”算法 为了力求改进计算收敛的速度，已经制定出了一个修订的版本，这个版本叫 做s i m p l e