（流体力学专业论文）方腔热对流和方柱地效绕流的有限元分析.pdf

摘要 摘要 本文采用亳隰垂鎏数值求解了b o u s s i n e s q 近似下的原始变量不可压缩 n a v i e r s t o k e s 方程，研冗了赢堕堡凼自筮型薅和考虑地面效应的燮的流动。 通过求解二维n a v i e r s t o k e s 方程， 狭缝内的自然对流，狭缝长细比a = i ：2 0 本文首先模拟了一侧铅直板加热的二维 流动参数取p r = 0 7 1 ( 空气) ，r a = 1 0 3 1 0 6 。 随着r a 数的升高，狭缝内涡的个数由一个变到多个再回到一个，流动由稳定 发展为不稳定。 在三维底板加热的方腔自然对流中，流动参数取p r = o7 1 ，r a = l5 x1 0 ，腔 体底板为正方形，计算了高度分别为底面边长l ，2 ，l 3 ，i 4 ，1 5 的方腔中的 流动现象和热现象。饼。算发现，高度为1 2 和1 4 时，流场得到定常解，h - 腔内 分别产生2 个和4 个长直涡；而高度取为l 3 和l 5 时，得到的却是周期解。在 高度为l 3 的方腔内，存在一个绕方腔中心的主涡环，而高度为l 5 时则得到内 外两个同轴而旋向相反的涡环。；、。 绕方柱的流动，考虑了地面边界层的影响，流动参数r e = 1 0 0 0 ，方柱高与底 边之比为2 ：l 。在二维流动中，从柱体顶端生成的顺时针涡与从底部生成的逆时 针涡周期性交替脱泻。三维方柱的绕流中，得到了流场中迎风面上的马蹄涡， 顶部和侧壁的分离涡以及背风涡等复杂流动现象。 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 ，t h en a t u r a lc o n v e c t i o nf l o wi nac a v i t yi n d u c e db yt h eh e a t e d b o s o ma n dt h es e p a r a t e df l o w sa r o u n dav e r t i c a ls q u a r ec y l i n d e ri ng r o u n de f f e c t sa r e i n v e s t i g a t e db ym e a n s o ft h ef i n i t ee l e m e n ts o l u t i o nf o rt h e i n c o m p r e s s i b l ef l o w e q u a t i o n s t h en a t u r a lc o n v e c t i o ni nat w od i m e n s i o n a lv e r t i c a ls l o ti ss t u d i e df i r s t l yt h e a s p e c tr a t i o i s2 0a n dt h ep r a n d t ln u m b e ri so71 a st h er a y l e i g hn u m b e ri n c r e a s e s f r o m10 3t o10 6 ，t h et r a n s i t i o nf r o mu n i c e l l u l a rf l o wt oam u l t i c e l l u l a rf l o w , a sw e l la s f r o ms t e a d yf l o wt ou n s t e a d yf l o wo c c u r sa st h er a y l e i g hn u m b e rf u r t h e ri n c r e a s e s t h ei n v e r s et r a n s i t i o nf r o mm u l i t c e l l u l a rt ou n i c e l l u l a rf l o w a p p e a r t h et h r e ed i m e n s i o n a ln a t u r a lc o n v e c t i o nf l o w si na c a v i t yi n d u c e db yt h eh e a t e d b o t t o m ，w i t hf o u rd i f f e r e n ta s p e c tr a t i o s ( t h eb o t t o mi sas q u a r ew i t ht h eu n i tl e n g t h s i d ew a l l sa n dt h eh e i g h t so f t h e c a v i t ya r e05 ，03 3 3 ，02 50 2r e s p e c t i v e l y ) h a v ea l s o b e e ns t u d i e di nt h ec a v i t yw i t ha s p e c tr a t i o05a n do2 5 t h ea e a d ys o l u t i o n sa r e o b t a i n e d i nt h ec a v i t yw i t ha s p e c tr a t i o03 3 3a n d0 2 ，t h es o l u t i o n sp e r i o d i ci nt i m e a r eg a i n e di na d d i t i o n ，t h ea x e so ft h ev o r t i c e si nt h ec a v i t yf l o ww i t ha s p e c tr a t i o05 a n d02 5a r es t r a i g h ti i n e sa n dp a r a l l e lt ot h eh o r i z o n t a la x i s h o w e v e r i , t h e v o r t i c e sn t h ec a v i t yw i t ha s p e c tr a t i o0 3 3 3a n d02a r ev o r t e xr i n g sw i t ht h e i rc e n t r a ll i n e s p a r a l l e lt ot h eh o r i z o n t a lp l a n e u p o nt h et w od i m e n s i o n a l f l o wa r o u n dav e r t i c a l s q u a r ec y l i n d e r i n g r o u n d e f f e c t s ，ap e r i o d i c a ls o l u t i o ni so b t a i n e d t h ev o r t e xg e n e r a t e do nt h et o pa n do nt h e l e es i d eo ft h es q u a r ec y l i n d e rs h e df r o mt h ec y l i n d e ra l t e r n a t i v e l y i nt h en u m e r i c a ls i m u l a t i o no f t h ef l o wp a s s e dat h r e ed i m e n s i o n a lv e r t i c a ls q u a r e c y l i n d e ri ng r o u n de f f e c t s ，t h ev o r t i c e sa r eg e n e r a t e db yt h es h e a rl a y e r ss e p a r a t e d f r o mt h et o p e d g ea n ds i d ee d g eo ft h ec y l i n d e rr e s p e c t i v e l y a l lt h e s ev o r t i c e s i n t e r a c te a c ho t h e ra n df o r mt h ec o m p l e xf l o ws t r u c t u r e s 第一章绪论 第一章绪论 1 1 引言 流体力学是研究流体的机械运动和力的作用的科学，它广泛应用于航空， 航天，丈，海洋和生物科学等领域。早先的流体力学研究大部是通过理论分 析的方法进行的，即先建立数学模型( 控制方程和边界条件) ，然后将问题转 化为对微分方程组的求解，此时解决的问题大多局限于线性问题。但是大多数 的流体力学问题在理论上表述为非线性的数学问题，求其精确解或理论解是 | 分困难的。 而且传统的理论方法和各类解析解都是在各种简化假设条件下得到 的， 这些方法和结果无法描述大量复杂的流动现象。随着科学技术的发展和人 类对物理世界研究和思考不断深入，人们迫切需要对更多更复杂的流动现象有 深入细致的了解。幸运的是现代的实验手段和电子计算机在力学领域的应用为 流体力学的研究注入了新的活力。特别是流动的数值模拟以其成本低，应用范 围广和灵活性高等优点倍受青睐而得到广泛的研究和应用，从而形成和发展了 计算流体力学这一新的学科，并逐渐成为研究流体力学的有力工具。 计算流体力学( c o m p u t a t i o n a lf l u i dd y n a m i c s ，c f d ) 是将流体力学方程， 定解条件一边界条件和初始条件和求解域进行离散，得到离散的代数方程组 和离散的计算区域，再用代数方程组的解求得原来微分方程组的近似解。它给 出的是某种特定的流体运动的数值解，而非能描述一类流体运动的解析解，故 从一次数值解中无法看出它随流动参数变化的趋势。但它可以以比实验少得多 的花费给出流场内细节的定量描述，且若数值模拟的数学提法适当，则可在较 广的流动参数范围内较快地给出流场的定量结果，而不受实验中固有的约束冬 件的影响。此外从实际当中抽象出来的数学提法往往是十分复杂的多维非线性 偏微分方程组，其数值解的数学理论尚研究得不够充分，如严格的稳定性分析， 收敛性，唯陛和误差估计等理论的发展还跟不上数值模拟的需要，因而对于 数值计算的结果，也往往要与实验结果、物理分析相结合来验证数值解的可靠 性。 自然界中的流动主要有层流和湍流两大类，对层流的数值模拟方法主要为 有限差分法( f i n i t ed i f f e r e n c em e t h o d ) 和有限单元法仃i n i t ee l e m e n tm e t h o d l ： 对湍流，还要设计适当的湍流模型。有限差分法在c f d 中是主流方法，发展电 较为完善，有一整套的定性分析理论。针对有限差分法，人们已经发展了各种 高精度，高分辨率格式，如总变差减小格式( t v d ) 1 ，无自由参数的耗散格 式( n n d ) 2 ，基本无振荡格式( e n o ) 3 】3 以及矢通量分裂格式 4 】等。这些格 第一章绪论 式可以较好的模拟包含激波，旋涡等现象的非光滑流场。有限差分法的本质是 用a f l a x 这样的差分代替微分方程中的导数项可瓠。显然， 这种近似易于 在正交的直线网格上进行。因此有限差分的计算常常要将物理空间中非正交的 曲线网格变换到计算空f 司中的直线正交网格，同时也要将方程变换到计算空间 中去。结果，这样的变换给计算带来一定的复杂度。 如果说有限差分法用的是“点”近似，只考虑差分网格点上的参数值而并 不给出节点附近参数的变化，有限元法则用的是分段( 块) 近似，每一段( 块) 用某种多项式来逼近。有限单元法的思想始于四十年代，首次明确提出有限单 元法的概念是c l o u g h 5 。日前有限元法已广泛应用于固体和结构力学领域。它 的数学基础是变分原理和加权余量法，由此导出的变分表达式或加权余量积分 表达式是有限冗解法的出发点。固体和结构力学中的有限元法大多足从变分原 理出发来构成的，通常称为变分有限元法。流体力学中由于其控制方程一一 n a v i e r - s t o k e s 方程中对流项的非对称性，可能不存在相应的变分原理，即不存 在这样一个泛函，使该泛函的欧拉拉各朗日方程为n a v i e r s t o k e s 方程。这一点 已被f i n l a y s o n 对定常的n s 方程所证明 6 。因此，对于n s 方程，应用有限 元法多数是从加权余量法出发。加权余量法在流体力学中的引入源于两位固体 力学家z i e n k i e w i c z 和c h e u n g 7 的开创性工作，他们论证了有限元法解决势流 问题的可能性。g a l e r k i n 方法是最为简洁直接的加权余量法【8 ，它是将流场的 求解区域剖分为互不重叠的子域( 单元) ，用单元基函数的线性组合来逼近微分 力程的解。有限元法应用于流体力学有以下几个好处：【有限单元法对网格的 正交性要求不象差分法那么严格，因为它是针对每一个网格单元进行计算，只 要能保证彤成刚度矩阵后对角元素占优，网格有一定程度的斜交仍是可以接受 的。2 有限元对所考虑的计算区域的形状没有什么特殊的要求，易于处理复杂 区域形状的流动问题。3 有限元在选择近似函数时不涉及边界条件，这使得每 个单元具有相同形式的单元方程，它的求解步骤几乎是统一的，易于编制成通 用程序。这类程序，如国外n a s t r a n ，a d i n a ，s a p d 等商用软件。几十年 来，有限元法在流体力学中的应用有了长足的发展，但也存在很多问题，有待 进步解决。例如在求解对流占优的对流扩散问题时，由于对流效应大于粘性 扩散效应 9 ，流场中往往包含薄的边界剪切层和内伏剪切层，且方程组非线性， 其流场内察的复杂度增大，计算容易发生不稳定。尤其是在求解不可压缩流时， 由于连续方程形式的特殊性，反而增加了数值求解的难度：加之每一个时闻步 鄙要经过总刚矩阵的求解和合成，计算量很大，计算缓慢。 第一童绪论 1 2 流体力学的有限元方法 有限元法是在六十年代开始应用于流体力学的，从最初的只能求解位势流已 发展到现在的能够求解不可压的n s 方程，高速的可压缩粘性流 1 0 和气动 热力学 1 2 1 3 等复杂问题。数值方法也由单一的g a l e r k i n 方法，发展到 p e t r o v g a l e r k i n ，s u p g ( s t r e a m l i n e u p w i n d p e t r o v g a l e r k i n ) 1 4 1 5 1 1 1 6 】， g l s ( g a l e r k i n l e a s t s q u a r e s ) i7 】 i 吼c a u ( c o n s i s t e n ta p p r o x i m a t eu p w i n d ) i9 等 方法。 g a l e r k i n 有限元方法被最早引入流体力学的有限元计算 2 0 2 1 。研究和数 酒实验表明，当计算高r e 数流动时，由于方程中的非线性项占主导地位，数 值解在节点问产生伪振荡，一般认为用对称的g a l e r k i n 格式处理方程中的非对 称项是产生振荡的原因之一。g a l e r k i n 格式类似于有限差分法中的中心差分 2 2 1 ，差分法中处理这类数值振荡的办法是将微分算子离散为非对称的迎风算 子。1 9 7 6 年，z i e n k i w i c z 等 2 3 】指出了上述g a l e r k i n 方法中的缺点并首次提出 了一维格式的修正。随后，h e i n r i c h 等 2 4 提出了针对无源项一维的定常常系 数对流扩散方程的p e t r o v g a l e r k i n 格式，该格式在节点上得到方程的精确解。 有0 于g a l e r k i n 格式的是，p e t r o v g a l e r k i n 权函数不同于基函数，其方法是将 权函数取为基函数加一个高阶扰动函数，高阶函数的加入，使单元的权重更偏 向于来流方向的节点，起到了迎风的作用。不久，h u g h e s 2 5 通过修改对流项 的数值积分，将数值积分中值点由原来的单元中点取在靠近上游的某处，从而 成功的求解了一维定常对流扩散方程。到此，不难看出，上述几种方法仅仅限 1 二定常无源的对流扩散方程，仍难以向多维流和非定常流推广，因而学术界对 它们的评价也是毁誉参半。进入八十年代后，b r o o k s & h u g h e s 2 6 1 和j o h n s o p & n a v e r t 2 7 分别独立发现，为消除数值振荡，只需在流线方向加入人工粘性即 可。他们在原有o a l e r k i n 有限元的基础上，在流线方向上加一个迎风的扰动项。 这便是s u p g 方法或称s d 方法( s t r e a m l i n e d i f l u s i o nm e t h o d ) 。后来，h u g h e s 又 在s u p g 中加入了间断面捕捉项 2 8 1 ，可以较好的模拟边界层和剪切层。为了 加速g a l e r k i n 解的收敛，h u g h e s 将最小二乘项加入到g a l e r k i n 方法中去，从 而形成了g l s 方法 1 8 】。最小二乘项的加入同时提高了解的稳定性和精度。 第章绪论 1 3 不可压缩流n a v i e r - s t o k e s 方程的求解 非定常不可压缩流n a v i e r s t o k e s 方程的求解，受到了计算流体力学界的 广泛关注，其重要特性在于，不可压缩非定常n s 方程组是混合型的，即含有 抛物型和椭圆型两种性质的微分方程 2 9 】。连续性方程中缺少动量方程所含有 的时间导数项，因而不能对连续性方程用时间推进求解。 = 维n s 方程的有限元数值解法大体可分为以下三种类型3 0 ：即涡一流 函数形式，流函数形式和原始变量形式。“涡一流函数法”就是在二维流情况 卜，通过引入流函数来满足连续性方程。本来，n s 方程组要求解三个方程， 汁算量较大，而涡一流函数方程组是由一个对流扩散型的涡量方程和一个 p o i s s o n 型的流函数方程组成，使求解变得较为容易并减少了计算机内存要求。 这是不少文献在平面和轴对称流动的情形下采用涡一流函数方法求解不可压缩 n s 方程的原因【3 1 ，3 2 。其缺陷在于，围壁边界上的涡量边界条件难以规定得 很合理。“流函数”方法在数值求解时有很多缺点，目前很少有人采用。求解 原始变量形式的n s 方程的工作一直受到人们的重视，一是因为求解的变量是 物理的，非常直观，避免了壁涡边器条件处理的困难，且适合于强调求解压力 场的情况。二是此法便于由二维流向三维流的推广。原始变量形式的n s 方程 的有限元数值解浊大体可分为两种类型，即s i m p l e ( s e m ii m p l i c i tp r e s s u r e l i n k e de q u a t i o n s ) 方法和分裂步数法 3 3 。这两种方法都是用由不可压条件得 出的压力来修正速度，在这一点上，它们是一致的。但在求解压力方程的过 霹中s i m g p l e 方法的每一个计算时间步都要先进行压力方程系数矩阵的计 算，因此，对计算量本来就很庞大的有限元法采用s i m p l e 方法不啻是雪j ： 加霜 3 4 ，3 5 ，3 6 。分裂步数法则是通过对动量方程取散度得到压力p o i s s o n 方程， 然后分别离散动量方程和压力p o i s s o n 方程，交替求解速度，v 和压力p 的值。 丰h 比而言，分裂步数法的压力系数矩阵可在计算开始时一次求出，可谓一劳永 逸 3 3 ，3 7 。 用有限元法求解原始变量n s 方程，有速度和压力的插值函数等阶和不 等阶两种处理方法 3 0 。阶数选择不当也可能导致数值振荡。在早期用原始 变量求解n s 方程的工作中，t a y l o r 和h o o d 利用等阶的速度和压力插值函 数解得的速度场较好， 但压力场的分布没有规律 3 s 。s a n ie ta l 【3 9 】指出等 阶插值函数往往容易产生奇异的总刚矩阵。为了克服上述缺陷，在选择速度和 压力插值函数时，须满足b b ( b a b u s k a b r e z z i ) 相容性条件 4 0 1 ，即对速度“， 第一章绪论 v 和压力p 要取不同阶的插值函数。一般，压力的插值函数比速度的插值函数 低阶。这种条件的引入，无疑增加了编程的复杂程度，降低了计算效率，并 目使边界条件的引入也带来了困难。为了解决这种矛盾，近年来，许多学者又 重新选用等阶的速度和压力插值函数。例如，h u g h e s 等利用等阶插值函数并 结合p e t r o v - g a t e r k i n 方法成功地模拟了s t o k e s 流动问题 4 0 。d es a m p a i o 4 i 1 等 证明当压力p o i s s o n 方程由最小二乘法给出时，用等阶插值函数将不受b b 条 件的限制。王演兴等 4 2 年1 j 用等阶原始变量的分裂步数格式，成功求解了二维 大攻角翼形的非定场绕流问题。 对于n s 方程的有限元求解，许多学者部曾参与讨论过如何处理压力方程 的边界条件这棘手问题。早期由于提压力物理边界条件十分困难，只好给出 人为的压力边界条件。后来，g r e s h o s a n i ( 1 9 8 7 ) 4 3 】通过对原始变量n s 方 程的分析，讨论了多种压力边界条件的提法，使压力方程的边值提法初步得以 澄清。 1 4 本文主要工作 本文以g a l e r k i n 有限元法数值求解三维不可压缩流n a v i e r s t o k e s 方程为手 段t 模拟j 方腔体内自然对流和带地面摩擦效应的绕方柱的三维流动。 一、，维和三维方腔体内的自然对流 二维晴况，研究了p r = 07 l ( 空气) ，r a y l e i g h 数在1 0 3 15 1 0 5 范围内边长比 为a 3 l ：2 0 的方腔内的对流运动。再现了随着r a 数的升高涡由一个发展到多个 再到一个以及流动由稳定到不稳定再到稳定的变化过程。 在乏维腔体底板加热的自然对流中，模拟了p r = 0 7 1 ，r a = i 5 1 0 s 的自然对 流在几种不同长高比( 底面为边长为l 的正方形不变，对高度进行改变，分别 取为底面边长的1 2 ，1 3 ，i 4 ，1 5 ) 的方腔体中的流动现象和温度分布。数值 计算结果表明，方腔内流场解的性态和涡结构受高度的影响很大，而且似乎有 规律可循。 二、方柱地效绕流的数值模拟 对于：维方柱的绕流，从方柱顶端产生的顺时针涡和从方柱背风面底部产 生的逆时针涡周期性交替脱出，形成涡街。流场得到周期解。 对于三维方柱的绕流，模拟得到了流场中的复杂涡结构，在方柱迎风面靠 近地面的地方有马蹄涡的生成，马蹄涡的形成较方柱顶部和侧壁的分离涡以及 背风涡晚；方柱的顶端和侧壁均有分离涡的产生和脱泻。 第二章控制方程和数值方法 第二章控制方程和数值方法 2 1 1 控制方程 2 1 数学模型 本文采用原始变量形式的不可压缩流n a v i e r s t o k e s 方程，其表达式为 掣+ ( 口v ) 口： d ， 可u = o l v p + u v2 0 + 7 p ( 21 ) ( 22 ) 其中，u 一= u i + v j 一+ w 云，为流场的速度；p 为压力ip 是流体的密度，在 不可压流动中，p 为常数；u 是流体的运动粘性系数，为流体受到的体积力。 般假设，体积力歹有势，即夕= 一v q b ，代入上式与压力项合并后有p + 肿 记为p ，这样体积力项不再在方程中显示出现。为书写方便压力项直接用p 表 示，不再加一撇。因此n a v i e r - s t o k e s 方程写为： 掣+ ( 口v ) 巧：一! 印+ 闪：口 o ，口 为计算方便，x 寸) b - n 进行无量纲化， 俨老班仁老 其中上为特征长度，u 。为参考速度，r e 为雷诺数。将各无量纲量代入方程 ( i ) 一( 2 ) 得到笛卡尔坐标下无量纲的不可压缩流n a v i e r s t o k e s 方程： 0 u ：o o x 亟+ 。丑：一至+ 上垫。 a i 瓠i 瓠j r e 瓠i 瓠i 这里为书写方便，略去了无量纲量上的撇。 ( 23 ) ( 24 ) 第二章控制方程和数值方法 2 1 2 边界条件和初始条件 1 边界条件： 速度边界条件的表达式如下： 蚓= 一，r ，刀， r = ，( 25 ) 其n = i 1c 等等， f 和i - 2 表示分段光滑的全场边界r 上的两段子边界，满足r l u l = r 和 一n l = 0 。6 。是d i r i c h l e t 边界r 1 上的已知速度分量；，是n e u m a n n 边界上 的切应力，n 。是边界r 2 上的单位法矢量。 2 初始条件： t = 0 时，取厅g ，o ) = 瓦g ) ，且满足条件v - 玩= 0 。( 26 ) 2 2 数值方法 本文前面曾经提到，当计算的网格r e y n o l d s 数r e 。 l 时( 其中 r e 。= l i i 1h v ) ，g a l e r k i n 有限元会导致数值解的振荡和失真。这是因为g a l e r k i n 有限元法等价于中心差分，可以证明，在均匀网格下，g a l e r k i n 有限元离散将 引起负的额外耗散，从而会导致数值不稳定。后来，h u g h e s 1 6 ，1 7 ，1 8 ，2 2 ，2 5 ，2 8 ，4 0 等人发展的s u p g 方法，等价于在流动方向加了一个人2 1 2 粘性u7 。在线性的一 维标量对流一扩散方程中，可以证明，当人工粘性系数取为 u = 半( r e 。) 其中 肥e 沪洲孚卜击 时，在均匀网格的节点上，中心差分法将给出精确解1 ”：“。这种方法在有限元 第二章控制方程和数值方法 计算中的实现只需将权函数取成 w = n + p 其中，为插值函数，尸称为扰动函数，且 7 尸= r a n 。，r = ( r e 。) “ 就会得到同样的超收敛解。这种方法就称为s u p g 方法。这里的萨常数，是方 程中的对流速度。 从一维线性对流扩散方程的分析中得出的这种处理方法，对构造求解 n a v i e r - s t o k e s 方程的稳定高精度有限元格式在理论上给出了一个重要提示。后 来，有限元学者围绕这一方法的推广和改进做了大量工作 1 7 ，1 8 ，1 9 ，4 4 】。遗憾的 是，这一方法在高维和非线性方程的求解中仍遇到不少困难，最突出的两点是： i 它在求解高维非线性的n s 方程时，仍然难以完全避免数值振荡和解的 失真： 2 在时间推进过程中，每一个时间步都要重新计算扰动函数p ，并改变方 程的系数矩阵，这导致了庞大的计算量。 注意到当网格充分细密以致r e f 1 时，函数l f ，( r e 。) - - - 0 ，u 7 斗0 ；因此， 采用细密的网格应同样可以获得避免人工振荡的高精度g a l e r k i n 有限元解。于 是，我们用普通g a l e r k i n 方法和s u p g 方法进行了数值试验，发现对于求解n s 方程，为实现无振荡的高精度有限元解 用s u p g 法更费时或精度更低。所以， 的g a l e r k i n 方法，而未用s u p g 方法。 采用加密网格的g a l e r k i n 方法未必比 本文在正式计算中采用了适当加密网格 换句话说，提高s u p g 方法的计算效率 和消除s u p g 方法中的经验系数仍然是两个有待解决的问题。不解决这两个问 题，s u p g 方法的应用将有很大的局限陛。 不可压缩流动的n a v i e r s t o k e s 方程的解法大体可分为两种，一是同时求觯 连续性方程和动量方程，t a y l e r & h o o d 曾指出 4 5 ，这种解法要求速度的插值 函数比压力的插值函数高一阶。另一种解法是分裂步数法f 3 3 ，3 7 ，即在一个时 间步内将求解过程分为几步进行，先求解出中间速度( 这一速度一般不能保证 不可压缩条件) ，然后利用由动量方程和连续性方程导出的压力p o i s s o n 方程解 出压力，最后用压力来矫正速度使得散度为零的条件得以满足。 本文采用四步分裂格式对原始变量的不可压缩n s 方程进行g a l e r k i n 有限 元法数值求解。下面对这一格式详细阐述。 第二章控制方程和数值方法 2 2 1 四步分裂格式的一般形式 第一步，先将压力从n a v i e r s t o k e s 方程中解耦，采用半隐式求出第一个中间速 度： 华+ 喇= 叫+ 去j ( 27 ) 第二步，利用上个时间步的压力求出第二个中间速度： 掣叫 ( 2 s ) 第三步，根据上述两步计算，得到了两个中间速度“ 1 ) 和”j ”，但无论是“还是 “p 都不满足连续性方程，即不可压条件。为此，对式垡之 兰= 一“求散度， 并令v 驴”1 = 0 ， 则得到如下形式的压力p o i s s o n 方程： n n ，+ l ：坚a t 第四步，用压力对速度进行矫正，得第n + l 时间步的速度”+ 掣叫 2 2 2 压力方程的有限元格式 ( 29 ) ( 21 0 ) 压力的求解通常有两种方法。一种是直接解压力p o i s s o n 方程( 2 9 ) 式，在( 29 ) 式两边乘以权函数并进行分步积分，将压力降阶。这就需要人为地加入一个压 力边界条件，n g i r 梯度为零的条件，即印& t = o 。另种解法是避开求解压 力p o i s s o n 方程，从连续方程和压力p o i s s o n 方程推导出求解压力的替代表达式。 过程如下： 连续方程的g a l e r k i n 有限元形式为：l p ，“i , r 擒= 0 ( 2 1 1 ) 考虑到b ；妒a = 妒。，+ p ，和散度为0 的条件，( 2 t 】) 式可写成： 扣，“? “m = j ，哆m ，d f( 21 2 ) 第二蕈控制方程和数值方法 珂方程( 28 ) ，( 21 0 ) 得? “= “? 一o ( a t 2 ) ( 21 3 ) 将方程( 21 0 ) 代入( 21 2 ) ，并结合式( 21 3 ) 得到下式： 出l 吼p ：“艘：f 0 妒。叫2 ) d f ! 一r 妒，“d i - ( 21 4 ) 在计算中，用方程( 21 4 ) 代替方程( 29 ) 在四步分裂法中的位置。使用该表 达式求解压力，不需要压力边界条件。但随之而来的是需要给出边界上的速度。 对于类似腔体等内部绕流，壁面上速度为零，速度边界条件易于给出，方程 ( 21 4 ) 中右端第二项直接变成零。但对于机翼等绕流问题，除非远场给得足 够远，可以提d i r i c h l e t 条件，否则，出口的速度值难以直接给定。所以两种办 法求压力，各有利弊，应根据具体问题进行选择。 2 2 3 四步分裂的g a l e r k i n 有限元形式 本文采用g a l e r k i n 有限兀格式对分裂得到的各方程进行离散求解，抽值函 数祁权函数均取线性函数，( i = l ，8 ) ，破的形式将在本节的第四部分介绍。将分 裂格式中各步方程( 27 ) ( 2 1 0 ) ，以及( 21 4 ) 的两端乘以庐，并在单元内积分得到 g a l e r k i n 形式的有限元离散方程： 第一步， f ：卜( 型枷。+ + 聂1 饥， 艘一去：吼正= 。c z - s a ， 第二步， 如。瓮坐犯= n 妒m s n ) 第三步， 型f 。，p ：“锄= 一i a g o j 。1 2 d u z + f f r 妒，p 和。订 或 ，l p 。p ：“锄= l 妒。“ 2 抛一妒，“，肝。订 ( 2 15 c ) 第四步， 如笋艘= 一胁1 艘 皿例， 自然边界条件已经隐含在以上的方程组中，计算中只需处理本征边界条件 第二章控制方程和数值方法 2 2 4 有限元总刚系数矩阵和总刚方程组 在得到每个单元的系数矩阵后，不能孤立地求解一个个单独的小单元，而 ，必须考虑到单元之间的关系。因此，为了建立单元之问的关系，需将单元的系 数矩阵合并成总体系数矩阵，即总刚矩阵。用总刚系数矩阵表达的总刚方程组 形式如下： 甜p 1 坩j 1 ：2 v w ；2 + a tb = a u 据；1 v j l 3 1 少 v 尹 w 尹 f 一 十聂 + a t 几口 f 2 口 f 3 u 甜；1 ) v ；1 ) 1 = a 。 描? v ? ，1 f 一叫，2 4 f 3 。 ( 2 1 6 b ) 可= 古( 0 1 。“；2 ) + q 2 , jv p + 0 3 0 w ；：) ) 一古厶( 2 1 6 c ) “ v ? “ w 芦 f 1 口 ，2 。 f 3 f 譬 共o p ，音系烈旭眸；3 - 义利彬a 划r ； 爿。“= d 。，中，中，i 1 d 影，彬 b f = - “女+ c 水8 k + d 肚9 w ，a 巾， ”k k ，中，巾t ( q ，- v 。坞看妈，倒彬f c 咖也巾舭。：等坞。鲁坞：掣彬f 列州= ，姒杀鸲；鲁坞，等，掣梢 硝) 屯啪，等坞。杀坞、争d 倒彬f 门。= l 洲等坞：鲁坞：杀) 蝴 ( 21 6 d ) ( 21 7 a ) f 2 1 7 b ) ( 21 7 c ) ( 21 7 c 1 ) ( 2 1 7 e ) ( 2 1 7 t 3 佗1 7 9 ) 蜘，蚋j吣， ，、，、，l 凡 | i 第二童控制方程和数值方法 o ：等鹄，等坞，- 警- ) a 4 a 4 时1 = 舢篑鹄。等坞，等m 倒嘁 既= “等鹄：杀坞：等心d 彰悄 ( 21 7 h ) ( 21 7 0 ( 21 7 j ) 鲶= “等坞等坞，詈心蜊悄 犯t ， = 南l 扣罢坞嚣坞詈慨詈坞鼍坞鲁删栲+ 南“q 二詈鸲二等岷詈：詈岷鲁峨争d 彰彬+ 高“器坞，普坶，詈知，等坞，鲁鸲，等删孵 蛤鼻。中，u 。d r ( 21 7 0 f 21 7 r n ) a 。】为时间系数矩阵；b v k * 1 , c i j k 。j , d 一分别是x ，y ，= 三个方向上的对流系 数矩阵：e 叫是总的对流系数矩阵：毛是粘性系数矩阵；f 1 。，f 2 。，f 3 。是压力 系数矩阵：q l 口，q 2 。，q 3 ，l ，是求解压力时的中间矩阵。其中，中，为插值函数，本 文采用线性插值，中，的具体表达式如下： = 妻( 1 + f ，舌) ( 1 + 仇，7 ) ( 】十f ，f ) ，i = 1 ，8 ，单元节点编号如图所示 第二章控制方程和数值方法 8 个单元结点上的，r t ，取值见下表 l2345678 i 一111一li111 兀 一l，l111一l11 e 11111l1i 一 在得到总刚系数矩阵后，将本征边界条件嵌入到总刚矩阵中然后求解总刚 方程组。总刚方程组实际上是一个高阶的代数方程组。本文采用g a u s s s e d e i 迭 代对方程进行求解。采用g a u s s s e d e l 求解代数方程组，要求左端的系数矩阵对 角占优。只要单元的网格畸变不太厉害，一般情况下，有限元的总刚系数矩阵 能自然满足对角占优的条件。 第三章方腔中自然对流的有限元分析 第三章方腔中自然对流的有限元分析 3 1 热对流运动的数学模型 处于重力场中的流体，当发生热力不平衡或不稳定的温度分布时，会产生 热对流运动。在流场中的温度差不太大时，我们可以采用不可压缩流假设，重 力方向的动量方程也可以取b o u s s i n e s q 近似，这时控制方程为： 甲u = 0r 31 1 _ 0 u + ( 盯v ) 口= 一上一v p 十u v ：口+ 亭口。( 32 ) c t p o 掣十( 口v ) o ：胛：o( 33 ) f 其中，0 = t 一瓦，p 。= 常数为基准密度，季为重力加速度，a 是热膨胀系数， r 是热扩散系数。在自然对流中，不存在特征速度，因此不能采用r e 作为无 量纲参数。然而，这里存在温度差7 1 = 7 ；一兀，i 与瓦分别是相对两板上的设 定温度。正是温差的大小和重力加速度决定了流体运动的强弱。无量纲化的特 征长度取为腔体的三个方向的最大长度，记为，取特征速度为k i l ，特征时 m i 为上：f ，方程无量纲化后，可以得到两个无量纲参数： g r = g a a y 才v 2 ，p r = u k 分别称作g r a s h o f 数和p r a n d t l 数。代替国，也可以用另外一个无量纲数 r a = g a a t l 3 m ，r a = g rp r ，称为r a y l e i g h 数。 无量纲的控制方程为： _ a u + ( 口v ) 疗= 一即+ p r v ! 口+ j r d p r 。 ( 34 ) “ 掣+ ( 口v ) o ：v z 。 a ， 、( 35 ) 计算所取的初场为静止的均匀流体加上一个微小的初始扰动( 露，“，o ：) ，它 4 第三章方腔中自然对流的有限元分折 们为模数s 1 0 _ 5 的随机量，且v 瓦= 0 。试算表明，当r a 充分小时 ( r a 1 7 0 0 ) ，这些小扰动随时f 司衰减，只有在大r a 数下，才会发展成为有 序的自然对流。 方柱的壁面边界上速度均取无滑移条件，温度取底板o = 1 ，顶板 = 0 ， 其它壁面均为绝热面( 在二维狭缝流动中，取左壁石= 0 的温度0 为l ，右壁 x = l 的温度为0 ，上下两个壁面为绝热壁，仅此例外) 。 3 2 二维方腔中的自然对流 所n - - 维方腔是指个截面为方形的直管道，管道沿轴向两端无限伸长。 横截面取为驯平面，z 轴平行于直管道的方向。 3 2 1 二维正方形腔体内的自然对流 本文先进行了二维正方形腔体的一竖直板和一底板分别加热时两种自然对 流的数值模拟。前者旨在测试程序，后者的目的在于模拟二维正方形腔体底板 加热时内部流场的演化过程。竖直板加热的流场和温度场如图3 1 ( a ) ，( b ) ，图 3 2 ( a ) ，( b ) 为h g c h o i 等的计算结果，可见二者基本一致。 ，蓟。翘 图3 1 本文结果 画圉 图3 2 i - i g c h o i 等的结果 第三章方腔中自然对流的有限元分析 对于底板加热产生的自然对流，参数取r a = l5 1 0 5 ，p r = 07 1 ，时间步长 d t ：10 1 0 s 。流场的演化过程如图33 所示。由于流场的左右边界条件对称于 x ：o5 ，所以起始阶段( t = o0 2 ) 出现两个对称的大涡，左侧涡顺时针旋转，右侧涡 逆时针旋转，各占据半个流场。这样的两个大涡不稳定，很快，腔体顶部分离 出两个小涡，分布在两个角上。此时，流场的对称性己受到破坏，右侧稍强一 些。t = 00 5 时，角上的两个小涡得到发展而变大，并且右侧的小涡发展更快一 些。原来的两个大涡受到压缩而变小。t = 0 0 7 时，上边两个小涡得到充分发展， 流场被四个大小大体相等的涡所控制，腔体中心出现一个鞍点。t = 00 8 时，右 卜和左下两个涡合二为一，形成一个顺时针方向的大涡。此后，大涡不断加强， 压迫着右下和左上两个涡，使之变扁变小，如t = 00 9 时刻的流场图。t = 01 0 时大涡继续加强，右下和左上两个扁涡破裂，分布于四个角上。t = 01 1 时，右 上和左下角的两个小涡耗散掉，大涡占据了流场主体。t = 01 2 时，顺时针大涡 的结构基本形成，右下和左上有两个很弱的小涡，这以后，流场基本保持不变。 流场内一点速度随时间的变化曲线如图33 ，由图可见，t = 0 1 5 时，流场已无大 的变化。t = 02 后，流场已经达到稳定。 第三章方腔中自然对流的有限元分析 图33 底板加热的腔体内流场演化图 t r n ( a ) 图34 中点速度u 随时间的变化曲线 3 2 2 细长狭缝内的自然对流 本文还研究了铅直方向细长的矩形狭缝内自然对流涡结构随r a 数的变化。 此时，温差加在狭缝的两侧壁上，上下底为绝热边界条件。狭缝的长细比为 x ：y = l ：2 0 ，流体p r = 07 i ，r a 数的变化范围为1 0 3 1 0 6 。具体数值实验所选取的 r a 数为1 0 3 ，4 2 1 0 3 ，1 0 4 ，2 0 1 0 4 ，1 0 5 ，1 5 1 0 5 ，1 0 6 。网格数1 6 3 2 0 ，时 间步长出为5 l o 。5 第三章方腔中自然对流的有限元分析 边界条件：h = v = 0 ，当x ，y = o 或l 时 。i 。吐。i 。- o ，乱= 乱= 。 数值计算得到的流场如图35 示，流场中心点速度“随时间的变化曲线如图 3 6 ，涡的个数与r a 的关系见下表： a = 2