（流体力学专业论文）圆筒中RayleighBenard对流的直接数值模拟研究.pdf

摘要 摘要 本文采用直接数值模拟方法研究了底面加热、顶部冷却的竖直圆筒内 r a y e i 曲b e n a r d 湍流对流问题，圆筒上底面和下底面温度恒定，侧壁为绝热壁， 所有壁面均满足速度无滑移条件。通过对柱坐标系下的n a v i c r - s t o k e s 方程直接 数值求解，并引入径向通量代替径向速度来消除坐标奇异。得到了三维空间的速 度场，涡量场和温度场数据，分析研究了不同r 口数和形状参数下的 r a y l e i g h b e n a r d 对流。 研究了不同黝数和形状参数下r a y l e i g h b e n a r d 湍流对流问题。计算表明， 流场中存在的两种大尺度结构，即大尺度环流和轴对称涡结构，讨论了它们所呈 现的不同的流场结构及其在不同足4 数和形状参数下的变化。着重分析探讨了r a 数以及形状参数对于流场中心区域平均速度场和平均温度场、速度和温度脉动 量、上下壁面的速度和温度边界层、涡结构强度以及与湍流传热效率相关的统计 量n u s s e l t 数的影响特性。研究揭示了鼢数以及形状参数对于r a y l e i g h - b e n a r d 湍流对流的基本特性。 关键词：r a y l e i g h - b e n a r d 对流；直接数值模拟；大尺度环流结构：轴对称涡； 形状参数 a b s t r a c t a b s t r a c t i nt h i st h e s i s d i r “tn u m e r i c a ls i m u l a t i o nb ys o l v i n gt h et h r e e - d i m e n s i o n a l i n c o m p r e s s i b l en a v i e r - s t o k e se q u a t i o nu n d e rb o u s s i n e s qa s s u m p t i o ni sp e r f o r m e dt o s t u d yr a y l e i g h b e n a r dc o n v e c t i v et u r b u l e n c ei nac y l i n d r i c a lc e l lh e a t e df r o mb e l o w a n dc o o l e df r o ma b o v ea n df i l l e dw i t hg a s e o u sh e l i u m ( 丹= o 7 ) t h eb o t t o ma n dt o p p l a t e sa l et a k e nt ob ei s o t h e r m a la n dh e l da tc o n s t a n tt e m p e r a t u r c s 。a n da d i a b a t i c s i d e w a l l sa n dn o - s l i pb o u n d a r yc o n d i t i o na r ee m p l o y e d i nt h ep r e s e n tc a l c u l a t i o n ，a e o o r d i n a t es i n g u l a r i t ya tt h ea x i si se l i m i n a t e db a s e do nt h et e c h n i q u es u g g e s t e db y v e r z i c c o o d a n d ib yi n s t e a dt h er a d i a lf l u xi n s t e e do ft h er a d i a iv e l o c i t y b a s e do n o u rc a l c u l a t i o n s ，t h r e e d i m e n s i o n a lv e l o c i t y , v o r t i c i t ya n dt e m p e r a t u r ef i e i d sc a r lb e o b t a i n e da n da n a l y z e d t h eo b j e c t i v eo f t h i ss t u d yi st 0d e a lw i t ht h ei n f l u e n c eo f t h er a y l e i g hn u m b e r a n dt h ea s p e c t - r a t i oo ft h ec y l i n d e ro nt h ef l o wc h a r a c t e r i s t i c sa n dh e a tt r a n s f e ro f r a y l o i g h - b e r m r dt u r b u l e n tc o n v e c t i o n b a s e do nt h ep r e s e n tc a l c u l a t e dr e s u l t s ，t w o t y p i c a ll a r g e s c a l e s u u e t u r e so ft h ei i l c a nf l o wa r ei d e n t i f i e d o n ei sal a r g e r e i r e u l a t i o ne e l la n dt h eo t h e ra x i s y m m e t r i ct e r o i d a lv o t i c e s t h eo v e r a l lf e a t u r c so f t h em e a nf l o wa r ed e t e r m i n e db yt h er e l a t i v ei m p o r t a n c eo f t h et w os t r u c t u r e st h r o u g h t h ef o r m a t i o no f t h eb o u n d a r yl a y e r s ，s o m et y p i e a ls t a t i s t i e a lq u a n t i t i e s i n c l u d i n gt h e m e a nv e l o c i t ya n dt e m p e r a t u r e , t h e i rf l u c t u a t i o n s ，t u r b u l e n th e a tf l u x e s ，a n d t u r b u l e n c es t r u c t u t a s a r ea n a l y z e di nd e t a i l k e y w o r d s ： r a y l e i g h b e n a r dc o n v e c t i o n ，d i r e c tn u m e r i c a ls i m u l a t i o n ， l a r g e - s c a l er e c i r c u l a t i o n , a x i s y m m e t r i ct o r o i d a lv o r t i c e s , a s p e c t - r a t i o 中国科学技术太学硕士学位论文 第一章绪论 1 1 r a y l e i g h b e n a r d 对流的研究背景和意义 湍流运动是自然界和工程实际问题中存在的最为普遍的一种流体运动形式，也 是迄今为止仍属于物理及力学中其根本机理问题尚未获得解决的几个难题之一， 对它的研究长期以来一直都是流体力学学科的中心课题之一白1 8 8 3 年r e y n o l d s 的著名湍流实验以来，人们对湍流研究已有了上百年的历史，经过睦时间的研究 探索，并逐渐将其应用于解决工程等各种实际问题中。对湍流的研究不但有重大 的理论意义，而且有巨大的实用价值。飞行器的减阻、散热、复杂的化学反应流 动、河流，海洋以及大气流动等无不与湍流密切相关。 热湍流是湍流研究的重要问题之一，它既广泛存在于自然界，例如大气湍流边 界层、环境与海洋流动等，也大量出现在工程问题中，如核燃烧反应器、化工设 备等。这时热湍流和剪切揣流往往共存于同一个流动问题中通过湍流输运和热 量输运来实现其中动量、能量和物质的混合和交换( d a v i e s ，1 9 7 2 ) 。热湍流问题不 仅是流体力学的基础研究澡题，同时也是地球物理学、气象学，环境科学和海洋 科学的研究重点( b l a c k a d a r , 1 9 8 7 ) ，它与天气预报、气候预测等有着非常密切的关 系。由于人类的活动和工程应用几乎都是发生在大气层圈内，所以这类热湍流问 题的研究对于工业、农业、国防建设和环境保护等具有十分重要的现实意义。 在热湍流中，r a y l e i g h - b e n a r d 对流被公认为是研究热湍流对流的一个典型问 题，即所谓r a y l e i g h b e n a r d 对流，其基本的流动特征表现为：在一个封闭容器内， 其下表面加热而上表面冷却形成温差，从而导致容器内流体运动。r a y l e i g h b e n a r d 对流问题长期以来一直是研究工作者们热心关注的研究课题，早期对于该问题的 研究，是因为它可提供人们对纯半导体生产中对流现象的认识。 1 8 8 0 年c u r i c 兄弟首先发现了晶体的压电现象，1 9 1 7 年l a n g c v i n 利用石英板 研究了水中声波的产生和检测。此后，随着微电子t 业迅速发展，石英晶体在工 业，军事以及民用等方面的应用越来越广泛。因为天然晶体数量毕竟有限，难以 第一章绪论 满足大量的需求，人工合成就成为了解决这一矛盾的必由之路。二战以后在 n a c k z n 、b r o w n 、s a w y e r 等人的努力下，很多晶体生长技术应运而生，并逐渐走 向成熟自从1 9 世纪5 0 年代以来，各发达国家以及部分发展中国家纷纷建立自 己的晶体生产工业。由于晶体的性质在很大程度上影响着由它所生产出的半导体 性能，特剐是晶体的各项同性与否对半导体的优劣有菲常大的影响。晶体的很多 性质又与其制备过程育很大关系，因此对于晶体生长过程的研究就成为了必然。 经过大量的实验及计算研究，人们发现晶体生长在融化状态下的对流性质在其固 化过程中起着决定物质交换和热交换的巨大作用，而物质交换和热交换又对晶体 的性能起着重要的作用。所以，对于晶体生长过程中流动状态，尤其是其中的不 稳定对流的深刻坪解就变得日益重要。 另外，在很多实际问题中如大气物理问题，普遍存在不稳定分层的热输运现 象。例如大气边界层这样的湍流场，由于白天地面大量吸收太阳辐射后迅速增温， 边界层顶与地面之闻出现较大的负温度差，而这种温度逆差超过一定限度后，将 会产生不稳定分层现象。当热分层为不稳定时，会产生大尺度的热对流流动，这 也是较典型的r a y l e i 曲b e n a r d 对流。热对流的发生能够湿著改变近壁区湍流的动 力学将性，热输运机制以及传热效率。 随着对r a y l e i g h - b e n a r d 问题的研究深入，它可提供人们对星际、海洋、大气 对流现象的认识，井对耗散系统中非线性现象的精细研究提供理想的摸型。 r a y l e i 曲，b e n a r d 流动己经成为热揣流最为广泛研究的问题之一。 1 2r a y l e i g h b e n a r d 对流的研究概况 对于这一问题的研究始于本世纪初，1 9 0 1 年b e n a r d 通过实验首次观察到对流 现象，受到这一现象的徼发，1 9 1 6 年r a y l e i g h 在理论上进行了描述，建立了针对 这一现象的对流不稳定理论。自此以后，对于r a y l e i g h - b e n a r d 现象的研究就沿着 实验和理论分析两条途径发展。随着计算科学的发展，直接数值求解相应的控制 方程就成了另一研究途径，尤其是现代测试技术的发展和计算能力的提高为这 2 中国科学技术大学硕士学位论文 一问题的研究提供了强有力的手段。 与壁面边界层的发展相类似，r a y l e i g h - b e n a r d 对流亦存在层流、流动稳定性 问题及层流向湍流的转捩，湍流对流等状态。判断流动状态的参数为r a y l e i g h 数 ( 船) ，定义为 月n ：g a a t h 趔 v i f 式中g 为重力加速度，a 为热膨胀系数，为流体上下底面的温差，h 为对流容器 容器的高度，v 为流体运动粘性系数，k 为流体热扩散系数。 在r a y l e i g h b e n a r d 湍流对流的研究中m a l k u s 从理论上建立了这一问题的基 础，他假定在热传递过程中小尺度运动起主导作用，得到了表征热量特征的无量 纲参数n u s s d t 数与如数的关系式为n u r a ”。 厂 n u 2 旦a a tr 瓯。h v a 这里。代表容器内每单位面积传递的热量，q 容器内每单位体积传递的热量， 而且为热传导率。 以b o u s s i n e s q 方程为基础早期学者针对r a y l e i g h - b e n a r d 对漉的稳定性问题 进行了大量的理论分析工作。n e u m a n n 通过对稳定状态的线性化b o u s s i n e s q 方程 的求解得到了从稳定流动向非稳定流动转化的临界r a y l e i g h 数，并通过分离变量 法得到了针对非固壁边界条件以及固壁与非同壁耦台边界条件的解折解。d r a z i n 利用摄动方法求解了i 临界r a y l e i g h 数附近的流动。b e h r i n g e r 详细总结并完善了 r a y l e i g h - b e n a r d 对流稳定性问题及向湍流转捩的研究。 对于r a y l e i g h b e n a r d 对流的实验一直都是研究的基础。k r a i c h n a n 通过p e e l e r 数和r 口数确定温度边界层和粘性边界层的尺寸，得到了小p r 数和大p r 数情况下 热量、脉动速度、温度的表达式。h o w a r d 受当叫已有的实验结果的启发，建议温 度边界层的尺寸由相对于对流边界层边界稳定时所需要的尺寸来确定，得到了经 典的关系式n u r a “3 。 t o w n s e n d 也较早地对r a y t e i g h - b e n a r d 对流问题进行了实验研究，他对平均温 度和脉动温度均方根值进行了精细的测量，得到了与m a l k u s 理论分析相一致的结 果。从那时起，研究工作者选用不同形状的对流容器( 圆柱体、立方体、长方体 3 第一章绪论 等) 和不同的流体介质( 水、空气、氨气、水银等) 进行了大量的实验和数值模 拟研究。有的工作除给出温度场的特性外，还对速度场的统计特性进行了测量， 关于n u r a 关系式，c a s t i n g 等很多实验艘现公式中的r a 的指数略小于1 3 。 h e s l o t 等人在1 9 8 7 年对以低温氮气为介质，直径和高度均为8 7 c m 圆柱内的 r a y i e i g h - b e n a r d 对流问题进行了细致的研究，实验r a y l e i g h 数高达1 0 “，依据瞬 时温度信号特征和n u r a 关系曲线，他们将r a y l e i g h b e n a r d 对流分为对流开始 ( r a = 5 8 1 0 3 ) 、振荡( r a = 9 x 1 0 ) 、混沌( r a = 1 5 1 05 2 5 x t 0 5 ) 、转捩 ( r a = 2 5 1 0 一5 1 0 ) 、软湍流( r a = 5 x 1 0 5 4 x 1 0 7 ) 及硬湍流( r a 4 1 0 7 ) 等六个区域，特别是关于软湍流和硬湍流的划分，是对r a y l e i g h - b e n a r d 对流问题 研究的重要贡献。 由于粘性边界层在r a y l c i g h b e n a r d 湍流对流中很薄，并存在很大的温度梯度， 给精确测量速度分布带来了困难。x i a 等人提出了一种测量湍流对流边界层内流速 分布非常有效的方法一一双光束非相干频谱分析法，该方法具有空间分辨率高、 对流场无干扰等优点，已取得了有关1 l a y l e i g h b e n a r d 湍流对流速度场的大量信息。 此外，t i l g n e r 通过测量流场中相近的两个温度探头之间的互相关及特殊的流动显 示技术等，也得到有关流场的一些信息。 s a n o 等对半高处大尺度环流向上、向下两侧温差进行了测量。结果表明，这 一参数随r a 的变化可作为对r a y l e i g h b e n a r d 湍流对流分类的判断标准，在硬湍流 区，温差随r a 减小。通过对环流向下运动一侧温度瞬时信号的测量，得到热卷流 是温度脉动和水平温差的主要来源。b e l m o n t e 等对温度历时的分析表明，温度脉 动是由浮力引起的，并指出描述温度脉动动力学特征的时间尺度与垂向距离的关 系。 在r a y l e i g h b c n a r d 湍流对流的控制方面，c i l i b e r t o ( 1 9 9 6 ) 等的研究表明，当实 验甩矩形对流盒沿一边倾斜不超过i o 。时，大尺度流动与对流盒的相对位簧发生了 变化但并未改变对流及热量的传递特性，大尺度环流运动的周期也未发生改变。 当在对流盒内上下板沿垂向布置4 块厚1 0 m m 、高02 5 l 的网时，流动显示表明离 开边界的卷流垂直向上运动，大尺度环流被完全控制，温度t 沿垂向的分布与无 网时的结果相一致，正好说明网对热量的输送无影响。对流盒倾斜的极限状况是 4 中国科学技术大学硕士学位论文 侧板加热，b e l m o n t e ( t 9 9 5 ) 等通过温度测量和逆光摄影显示技术研究了侧壁加热对 方柱内气体对流的作用。大尺度环流在r a 3 x 1 0 7 时出现，湍流脉动限于冷热两 板附近的区域，并得到了与通常的r a y l e i 曲b e n a r d 对流相一致的边界层厚度关系 式。对流盒中流动整体上为稳定的分层流在中央处存在内波，这项工作提供了 一种在实验室研究由分层气体中昧动产生的重力波的途径。 和实验研究比，r a y l e i g h b e n a r d 湍流对流的数值模拟研究具有如下优点：( 1 ) 控制参数如置口、p f 可以线性变化，而不用考虑实验实现难度；( 2 ) 总是满足 b o u s s i n e s q 假设；( 3 ) 速度无滑移、绝热壁面等边界条件可以事先给定；( 4 ) 可以 方便的得到速度场等数据。因此，随着近些年来算法的改进和计算机运算能力的 提高，r a y l e i g h b e n a r d 湍流对流的数值模拟也取得了长足的发展。以无量纲形式 的b o u s s i n e s q 方程为基础，通过引入r a 和p r 以及采用不用的近似方法，人们已经 对一些物理机理进行了研究，例如对流不稳定的发生条件、边界条件和初始条件 对流动的影响等等，并得到了与理论预测较为吻合的结果。近年来研究的热点主 要集中在侧面的固壁对稳定性及对流状态的影响，以及各种参数对丁f 行管道、 柱状容器内的对流的影响。j o n e s 等研究了轴对称流动的不稳定行以及存在的条件， m u l l e r 和v c z i c , c o 等人分别研究了不同r 口下流动的特性。此外，r u d i g e r 等人还使 用谱方法对r a y l e i g h b e n a r d 流动的稳定性进行了研究。 1 3 本文主要研究工作 本文采用有限差分方法求解了柱坐标系下n a v i e r - s t o k e s 方程，直接数值模拟 研究了底板加热、顶扳挣却并且侧壁为绝热壁的圆柱容器内r a y l e i g h b e n a r d 湍流 封流问题。主要的研究工作有； 1 研究了不同r 口数下( r a ：2 1 0 6 2 1 0 8 ) 的网筒内r a y l e i g h b e n a r d 湍流对流问题，这里取固定p r 数( p r = 0 7 ) 和形状参数( f = l 2 ) 。 基于直接数值模拟研究，分析了在流场中存在两种大尺度结构，即充满 流场的大尺度环流和关于对称轴对称、上下反对称的涡结构，以及它们 对平均流场以及速度、温度边界层的影响。同时，着重分析探讨了r a 数 5 第一章绪论 对于流场中心区域平均速度场和平均温度场、速度和温度脉动量、上下 壁面的速度和温度边界层、涡结构强度以及与湍流传热效率相关的统计 量n u s s e l t 数的影响特性。 2 研究了不同形状参数r ( f = 2 ，l ，1 2 ，i 3 ，1 4 ) f 的圆筒内 r a y l e i g h - b e n a r d 湍流对流的问题，这里驭固定p r 数( o 7 ) 和r a 数( 1 0 7 ) 。 通过直接数值模拟所获得的流场、涡量场以及温度场数据，分析了不同 形状参数下平均流场的大尺度结构，讨论了形状参数对大尺度结构的影 响。此外，分析了形状参数对平均流场和温度场、速度和温度脉动量、 粘性和温度边界层的影响，尤其是小形状参数下，圆筒内流场的复杂化 导致湍流热传导效率降低，文中也给予了相应讨论，并从相关实验测量 数据中得到了证实。目前关于数值模拟研究形状参数对r a y l e i g h - b e n a r d 湍流对流影响的工作还不多见，也是本文工作的重点。 6 中国科学技术大学硕士学位论文 第二章数学模型和数值计算方法 2 1 湍流的直接数值模拟简介 湍流是其有多尺度结构的不规则流体运动，其不同尺度的运动之间的动力学关 系是十分复杂的。湍流的多尺度运动属于宏观流体运动，既使是最小的湍流耗教 尺度也远大于分子运动的平均自由程，因此无论多么复杂的湍流行为实际上均遵 循n a v i e r - s t o k e s 方程( 在通常的温度和压力条件下) 。然而由于该方程是强非线性的 偏微分动力系统，任何最简单的湍流运动问题的解析解都是不存在的。近代非线 性动力学系统理论证实：在某些初值和方程的参数范围内，非线性微分方程确定 性问题可以有不规则的渐近解。牛顿流体的揣流是在大r e y n o l d s 数下n a v i e r - s t o k c s 方程的渐近不规则解。因此，要得到湍流流动问题的解只能通过数值求解 n a v i e r - s t o k e s 方程，这种数值求解方法就是湍流的直接数值模拟( d n s ) 。d n s 没有 采用任何的模型近似，是一种精细程度最高的数值模拟方法，它应能给出所有尺 度的湍流脉动，而所有平均量都可以通过统计平均得到。直接数值模拟作为研究 湍流的重要手段，与实验技术、理论分析以及其它数值方法互相补充，在湍流研 究中发挥着重要的作用，并且随着计算方法的发展和计算机技术的提高，其应用 越来越广泛。 从2 0 世纪7 0 年代开始，湍流直接数值模拟成为研究湍流流动机制的主要工具 之一。最早的直接数值模拟诞生于美国的n c a r ( t h en a t i o n a lc e n t e rf o r a t m o s p h e r i cr e s e a r c h ) ，o r s z a g 和p a t t e r s o n ( 1 9 7 2 ) 完成了基于t a y l o r 微尺度定义的 r e y n o l d s 数( r e ，) 为3 5 的均匀各向同性湍流的数值模拟。他们对端流直接数值模 拟中最常用的数值方法谱方法的发展做出了重要贡献。 湍流的直接数值模拟方法要求在一个宽广的长度和时间尺度上准确地再现流 动发展过程。由于随着r e y n o l d s 数的增大，湍流运动的最小尺度越来越小，最大 和最小尺度之间的差距越来越大，湍流的直接模拟对汁算机条件和计算方法两个 方面提出了挑战。计算机不断发展的强大的计算能力促进了湍流的直接模拟= 作， 7 第二章数学模型和数值汁算方法 然而其高昂的计算成本使其目前还不可能在工程应用中占有一席之地。尽管如此， 并不能说直接数值模拟对工程应用问题毫无帮助。事实上，它对湍流机理的认识 起到了重要的作用。我们知道，在工程问题上目前使用的模拟方法主要是雷诺平 均方法( r a n s ) ，也有一些情况下采用大涡模拟方法( l e s ) ，这时湍流模型的适应性 和准确性就成为最为关键的问题。湍流模型的发展离不开直接数值模拟数据库提 供的大量信息的验证。因此，喜接数值模拟方法最切实的贡献是在工程应用的唯 象分析建模上( m o i n 和m a h c s h1 9 9 8 ) ，它能够对各种新的理论和统计模式进行检 验。雷诺平均应力方程构成封闭雷诺平均流方程的基础，它包括几个必须模型化 处理又难以通过实验测量获得的主要项，而在直接数值模拟方法中所有雷诺平均 应力方程中的分项都可以基于直接计算得到，尤其可以得到湍流扩散和耗散项的 详细信息，这是用来构造湍流模型的基础。同样直接数值模拟还可以对已有的统 计模型做出评估，进而研究模型的改进途径。d n s 不仅对湍流模型的构建具有指 导意义，也能使我们对湍流的物理规律和微观结构的认识更为深刻。 壹接数值模拟的另一重要贡献是能够获得实时的流动演化过程，提供湍流场的 高分辨的样本，有助于我们洞察形形色色的湍流结构，尤其是边界层的细微结构 的变化。比如，一个典型湍流问题所做的d n s 能够取得湍流物理的各种信息和细 节，建立起丰富的数据库，包括三维空间内的瞬时速度、压力和温度脉动等物理 量，这些数据库信息对于进行湍流机理的研究具有重要的价值。这其中有些信息 是在物理实验中难咀i 嚏| i 得或者根本无法测量的量，例如湍流场的高阶多点统计矩 就很难用实验测量，流场的脉动压强至今没有合适的测量方法，而d n s 可以提供 这些量的计算结果。直接数值模拟在湍流研究中是一种强有力的且不可替代的研 究工具，在某些情况下，如对于低r e y n o l d s 数下剪切流动中的近壁拟序结构 ( c o b t r e n ts t r u c t u r e ) 的研究，它所起的作用甚至己超过了实验。 除此之外直接数值模拟在湍流的控制方面将大有所为。目前在这方而开展了 不少研究工作，例如c h o i 等( 1 9 9 4 ) ，l u m l e y 和b l o s s e y 0 9 9 8 ) ，c h o i 等( 2 0 0 2 ) 等研 究报道。 追溯直接数值模拟的发展过程，可以注意到数值研究所涉及的湍流问题的复杂 性在明显提高。然而，虽然计算机的运算速度和容器在以几何级数的水f 迅猛发 8 中国科学技术大学硕士学位论文 展，但是由于在d n s 计算中，网格数要求达到r e y n o l d s 数的9 4 次幂的量级，所 模拟的这些流动仍然限制在低r e y n o l d s 数范围内。以槽道湍流为例，目前能够实 现直接数值模拟的流动r e y n o l d s 数在1 0 4 ( r e = u 。日y ) 左右。考虑到流动问题的 几何形状及复杂性等因素，目前的数值方法是多样的，例如有限差分格式、谱方 法和有限差分的混合格式，以及采用非结构网格的有限元方法等。 2 2 1 引言 2 2 不可压缩湍流模拟的计算方法 由于湍流运动固有的非定常特性，其数值方法上的时问推进是一个研究重点。 自从c h o r i n ( 1 9 6 9 ) 首次基于向量的h o d g e 分解发展了时间分裂步法( f r a c t i o n a ls t e p m e t h o d l 以来，这种方法已被人们广泛地使用并加以改进，但时间项的处理上都长 期沿用了一阶精度的差分。 之后，k i m & m o i n ( 1 9 8 5 ) 、l e m o l n ( 1 9 9 1 ) 以及r a i m o i n ( 1 9 9 1 ) 等人的工 作取得重要的进展，尤其是k i m m o i n ( 1 9 8 5 ) 对数值求解非定常原始变量的 n a v i c r - s t o k e s 方程做出了重要的工作。他们将基于近似分解技术 ( a p p r o x i m a t e f a c t o r i z a t l o nt e c h n i q u e ) 的分裂步格式，发展成为求解原始变量下的非 定常不可压缩n a v i e r - s t o k e s 方程的高精度差分格式，采用了a d a m s b a s h f o r t h 格式， 将时间推进提高到了二阶精度，同时进一步改进了边界条件的处理。由于计算采 用了交错网格，并通过分裂步的中间速度场来求解压力p o i s s o n 方程，使得在每个 时间步连续性方程都能准确地满足。l e m o i n ( 1 9 9 1 ) 对以上的分裂步法进一步做 了改进，将r u n g e k u t t a 预估校正技术( p r e d i e w r - e o r r e e t o rt e c h n i q u e ) 引入到三个分 裂子步的时间离散中，使得该格式在不损失计算精度的条件下具有眺下优点：1 ) 允许相对较大的时间步长( c f l = 以) ：2 ) 需要的计算内存较小。原有格式的空 间导数离散都是利用动能守恒的( k i n e t i c e n e r g y - c o n s e r v i n g ) 中心差分格式，r a i m o i n ( 1 9 9 1 ) 则将种具有高阶精度的迎风差分格式( u p w i n dd i f f e r e n c es c h e m e ) 掺合 到中心差分格式中，对空间导数进行离散。迎风格式的耗散特性可以衰减掉小幅 9 第二章数学模型和鼗值计算方法 扰动( s m a l l a m p l i t u d ed i s t u r b a n c e ) 中的高频部分，从而使得改进后的格式能够控制 计算中的混淆误差( a l i a s i n 8 e r r o r ) 。最后，r a i m o i n ( 1 9 9 1 ) 将改进后的分袭步差分 格式应用到槽道湍流的直接模拟，计算得到的结果和e e k e l m a n n ( 1 9 7 4 、1 9 7 9 ) 的实 验数据以及k i me t a l ( 1 9 8 7 ) 用谱方法得到的结果符合得非常好。 由于柱坐标下n a v i e r - s t o k e s 方程在r = 0 处具有奇异性，使不可压缩流动在 ，= 0 处的连续性方程很难精确满足。而用宜角坐标下的数值方法求解某些轴对称 流动，在数值上又有很多的不便，甚至会出现非物理的周向数值扰动。为了解决 柱坐标奇异性问题，v e r z i c c o o r l a n d i ( 1 9 9 6 ) 改进了梓坐标f - - 维不可压流动的有 限差分法。该算法是建立在k i m m o i n ( 1 9 8 5 ) 和r a i m o i n ( 1 9 9 1 ) 的工作基础上， 同样也采用了交锚网格，并在计算中引入新变量q ，= w ，代替径向速度v ，( 在r = 0 处径向变量q ，= 0 ) ，从而消除了，= 0 处的坐标奇异性。同时在求解中间速度场 的时候利用了上一步的压力值，大大简化了中间速度场的边界条件计算。o r l a n d i f a t i c a ( 1 9 9 7 ) 将这个算法应用到轴向旋转圆管辅流的直接数值模拟中，取得了很大 的成功。计算得到的结果非常准确地描述了旋转状态下管壁附近涡结构的变化规 律，湍流的统计量也和r e i c h b e e r ( 1 9 8 9 ) 的实验数据以及e g g e l sc ta 1 ( 1 9 9 4 ) 的 数值模拟结果符合得很好。 图2 1 几何模型示意图 1 0 中国科学技术大学硕士学位论文 2 2 2 控制方程 本文采用原始变量形式的n a v i e r - s t o k e s 方程来数值模拟r a y l e i g h b e n a r d 流动。 首先在直角坐标系下得到方程的无量纲形式，再将其转化到柱坐标系下。问题的 物理模型如图2 1 ，封闭圆柱单元内充满粘性不可压缩流体，侧壁以及上下底面均 为无滑移阻及绝热条件，下底面加热而上底面冷却并保持恒温。 在直角坐标下的有量纲形式的n a v i e r - s t o k e s 方程表达式如下： 塑+ 丝：o 西 a x , ，鲁= 谢，+ 导 p 尝= 毒卜詈 - p 詈柚p 百2 面【2 瓦厂p i 椰 其中瓦d = 扣丢 ( 2 1 a ) ( 2 1 b ) ( 2 i t ) q p 毛州罄+ 鲁一；詈m 罄 。= ；c 等+ 善) 2 + c a 一扣c 鼍，2 这里p 为流体的密度，p 表示压力，吒为应力张量，为流体动力学粘性系数，a 为第二粘性系数，e 为单位质量流体的内能，r 是温度， r 为热扩散系数，g 是 重力加速度，o 为粘性耗散，瓦为k r o n e c k e r 符号。对于理想气体和液体，分别 有e = c 。丁及e = c t ，q 和c 分别为气体的定容比热和液体比热。 根据b o u s s i n e s q 假设，如果在浮力驱动的流动中温度变化8 t i 4 , ，其密度变化也 很小，从而可以假设出了浮力项以外，密度的变化均可以忽略不计并假定在流 体的两个表面温度差较小时，户= 风n a ( r 一五) 】，岛为温度瓦的底板处流体的密 度，d 为体积扩散系数，并假定a 为常量。 对于大部分的真实流体在b o u s s i n e s q 近似中可以将p 、c 。、c 作为常量 看待，同时忽略第二粘性系数z 。也就是说，在热力学和动力学参量中，除了压 第二章数学模型和数值计算方法 力、密度以及浮力项以外，均可视为常数，并且在浮力项中假定密度p 随温度线 性变化。但这一近似仅适用于温差较小的流动，对温差较太的流动不再适用。 如果以u 代表特征速度，d 代表特征长度，阮一五i 代表温度差，那么由于内 摩擦与由于热交换引起的传热比为 屯“u 3 d 一。 a u -一：-二-一 p d c o ) d tp o f l r 0 一五p d “p c 阢一瓦p 在通常的情况下，由于为量级比u 1 弓- r , l a , j 、很多的量t 因此可以忽略 内摩擦所引起的传热。对于不可压缩流体，可以忽略由于可压缩性所 起的传热， 而对于理想气体，由于p = ( c 。一c ，) p t 以及d = i t ，所以 d eau，dt p 面”t “p p 面 综合上述各种假定以及流体不可压缩性质，并记p p o + g z 为新的p ，由方程 ( 2 1 a ) ( 2 1 c ) 可以简化得到b o u s s i n e s q 近似下的n a v i e r - s t o k e s 方程 o u _ z l = 0( 2 2 a ) 鲁+ 毒c 叩j ) - 一言w 若苦一唱c i 。 c z 2 b ， 塑+ 墨掣；。j 车 ( 2 2 。) 8 t 蠡，玉：玉， 其中y = u p ，对于理想气体热扩散系数k = k ( p 一。) ，对于液体e = k ( p o c ) 。这 里采用圆柱的高度 作为特征长度，上下底面温度差a t = 阢一五i 作为特征温度， 自由落体速度u = 胛砀作为特征速度，用它们对方程( 2 2 a 卜( 2 2 c ) 进行无量纲 化，并定义r a y l e i g h 致为r a 盅a g a t h 3 ( v r ) ，p r a n d t l 数为p r = y 盯，得到的无量 绸控制方程为 堕：o 觑 鲁+ 毒c 叩沪一詈+ ( 是) ”2 袭+ 矾 ，1 2 ， 中国科学技术大学颐士学位论文 塑型：! 翌生 o t 缸，f f r r a ) ”2 缸 ( 2 3 c ) 关于所研究具体问题的初始条件和无量纲控制参数等，在第三、四章将做详细 介绍。 在柱坐标下，当半径趋于零时，由于 ，并不同时趋于霉导致知，a r 出现奇 异性。为克服这一困难，采用v c r z i c c o & o r l a n d i 提出的一种柱坐标下三维不可压缩 流动的有限差分方法。他们的算法建立在k i m & m o i n 和r a i m o i n 的工作基础 上，采用分裂步法、近似分勰技术和交错网格来求解原始变量的n a v i e r - s t o k e $ 方 程，通过引入轴向通量o = 佛来代替轴向速度咋，以消除奇异性。 设柱坐标下，周向为目，径向为，轴向为：，为处理奇异性引入如下新变量： q ，= 一，q 口= p 日，吼= v ：。将方程( 2 3 a ) 一( 2 3 c ) 改写到柱坐标系下得到如下形式的 方程组，可以得到在惯性系中拄坐标下的关于吼的动量方程，其表达式为： 盟+ 盟+ ，丝：0 a ra e 恐 ( 2 4 a ) 西口i8 r q p q ，1a ( q f q f ) a ( q 口吼) d t，2 d r r 2 d 口出 ；一等+ ( 罢) “2 c c 昙r 警卜笋+ 砉鲁+ 挚吾等，4 b ) 陀 ；，耋藩8 , 爨1a q 嘉麓竹 卅 = - ，筹+ ( 甜p 石尹+ 专貉+ 争一静， ”“ 盟+ 三旦熊盟+ ! 旦堕d 盟+ 垫! 型 卉rd r ra o如 = 一篆+ ( 是厂哆昙p 等，+ 专势+ 等， 1 3 - 第二章数学模型和数值汁算方法 竺+ 三曼鲤+ 三曼盟+ 皇绁 8 tr o r r0 0a = 面矗e 蚤c ，争专害+ 窘， 2 2 3 数值方法 ( 1 ) 时间离散 ( 24 0 本文采用的时闻离散方法为混合三步r u n g e - k u g a c r a n k - n i c h o l s o n 格式，这是 由v e r z i c c o & o r l a n d i ( 1 9 9 6 ) 在k m & m o i n ( 1 9 8 5 ) 的分裂步法盼基础上改进得 到的其不同之处在于在计算中间速度时引入了上一时问子步的压力项，这样可 以简化中间速度场边界条件的处理。这种格式的优点是占用计算机内存较少。同 时，出于计算稳定性的要求，对流项的时间推进用三阶精度的显式r u n g e k u t t a 法， 而粘性项采用二阶精度的隐式c r a n k - n i c h o l s o n 格式，整个算法在时间上是二阶精 度的。具体做法描述如下，在一个时间步内的推进分成三个时问子步进行，同时 结合了近似分解技术进行求解。格式的三个子步分别为： 第一子步： g j = g ? + a t r i ( g ? ) + p h ( q 7 一) 一口i g , p “+ 口( a m + 4 ，+ a ) ( 口? + 2 7 ) j ( 2 5 a ) q ? 一g ；一g i f g ，矿( 2 5 b ) d ( q t ) = 0( 2 5 c ) 第二子步： g ；= g ? + a t e ：h ，( g ? ) + p ：日，( 口? ) 一口：g ，p 。+ 口：( 爿m + 爿。+ 也) ( 口? + g ? ) j ( 2 6 a ) q ? 一q ? = 一口2 心矿2c 6 b ) d ( 计) = 0( 2 6 c ) 第三子步 - 1 4 中国科学技术大学硕士学位论文 q ? = g ? + 口b ，n , ( q b + p , h ，( g ? ) 一口，g ，p 5 + 口，( 一。+ 爿，+ 以) ( q ? + q ? ) j 群”一牙= 一口3 f g ， d ( q 7 “1 = 0 这里，是一个有关压力p 的变量，它们二者之间的关系为 一一譬( 黟v 弘p e 和d 分别表示离散的梯度算子和散度算子，置表示非线性对流项 h ，( 菇) 还包含了浮力项z e ，。4 。、4 ：和如表示粘性项的离散微分关系， g ，2 妄，只 ) 一吼善，。= c 毒”毒 毒啪 如= ；( 罢厂蓦，以= j 1 ( 面p r ) “2 筹，以= ( 劫“2 蔷 纠瓯和j 2 掰分别表示一阶和二阶空间导数的差分算子。 n ，n ，口，为三步r u n g e k u t t a 法中的系数( s p a l a f t 等1 9 9 1 ) ，分别为： = s 1 5 ，n = o ，= 8 1 5 扎= 5 1 2 ，岛= 一7 6 0 ， = 2 5 ，= 3 4 岛= 一5 1 2 ，= i 3 由式( 2 5 a ) 、( 2 ，6 a ) ，( 2 7 a ) 在计算中间速度场g _ q 7 g ? 时可把这三个公式统一 写为， o ：- 蜥a , 4 l h , ? - + a , p 4 7 h , - _ p 一, a 。, , 叫y q ：- tq 。：。) l a | e ，a ? ，a 澜 旺、 = 啦l