（流体力学专业论文）BGC圆筒热对流的多稳定解现象和RM不稳定性研究.pdf

摘要 摘要 本硕士论文主要通过数值模拟的方法研究流场流动的稳定性问题。本文从 直接数值模拟和整体稳定性角度着重对三维圆柱b r i d g m a ng r o w t h c o n f i g u r a t i o n 条件下热对流的流动特性，临界g r a s h o f 数，流动分叉类型， 多模态共存等方面进行深入分析。此外，本文还直接数值模拟了激波作用下， 单模r i c h t m y e r - m e s h k o v 不稳定性交界面的发展。本文的主要工作及研究成果 包括： 1 基于柱坐标下的有限差分方法，发展了一套研究b g c 条件下圆筒热对流 全三维整体稳定性分析方法。该方法利用柱坐标系下三维不可压流动的改进投 影算法，迭代得到一定研数下三维圆柱b g c 轴对称热对流的定常流场，再用 a r n o l d i 方法求出该定常流场的优势特征值和特征向量。 2 采用分叉分析和直接数值模拟相结合的方法，对p r 数为0 ：0 1 5 ，尺度 比彳= 4 情况下三维圆柱b g c 热对流多态共存进行了细致的研究。研究表明，轴 对称流在西在约5 2 1 4 1 时产生第一次叉式分叉，得到一支非轴对称定常流动分 支。在g ，数增加到约1 1 4 9 8 2 ，轴对称定常流又变为稳定的，与此同时，稳定 的非轴对称定常流动仍然存在。在1 1 4 9 8 2 研1 4 0 2 9 8 的区域里，随着函数增 加，先后有稳定的轴对称定常流与稳定的非轴对称定常流，稳定的轴对称定常 流与非轴对称周期流，稳定的轴对称定常流与非轴对称调制振荡流，稳定的轴 对称定常流与混沌流共存的现象。当g ，数大于1 4 0 2 9 8 的时候，轴对称流出现 第二次失稳，分叉出另一支非轴对称周期流动。此后，顺序出现了周期流与混 沌流，两个不同的周期流，两个混沌流共存的现象。 3 采用高精度的多介质g h o s t f l u i d 方法，对马赫数为1 1 5 的激波分别 作用于单模大扰动a i r c o ：，a i r s f 。，a i r n ：和a i r h e 交界面后的 r i c h t m y e r m e s h k o v 不稳定现象进行了数值研究，得到了不同时刻扰动界面的 演化图像，给出了流场的密度等值线和密度纹影图，同实验结果符合得较好。 本文给出了交界面的扰动增长随时间变化情况，并同理论模型进行了对比。对 激波从轻气体进入到重气体的情况，扰动增长可采用s a d o t 模型描述线性阶段 和早期非线性阶段；对于弱激波同密度接近交界面相互作用，线性阶段时间较 长，可以用线性模型描述。 关键词：b r i d g 髓ng r o w t hc o n f i g u r a t i o n ，叫l t 卜s t a t e ，流动稳定性，非轴 对称，r ic h t m y e r m e s h k o v 不稳定性，单模，非线性扰动 + 国家自然科学基金项目( 1 0 6 0 2 0 5 6 ，1 0 3 7 6 0 3 5 ) 资助 a b s t r a ( 玎 柚da i r - h ea r ep r e s e m c d 谢t ht l l e i re v o i u 吐o n so fm ei n t e r f a c e ，i n c l u d i n gd e l l s i t ) r 锄d s t 砌i n gc o n t o u r s t h es i m u l a _ t i o n 锄p l i t u d e sa r cs h o w nt 0b ei ng o o da g ，e e m e n t w i t l lt l l ee x p e m e n t mi i a ：t aa i l d 、i t l ln l ep r e d i c t i o 螂o f 廿l e o r e t i c a lm o d e l s t h e p e 咖r b a t i o ng r 0 、t l if o rt h el i g h 讹e a v yc 勰ea g r e e sw e l l 埘lt 1 1 en o i l l i n e 盯m o d e lo f s a d o t 甜口，i l lt 1 1 ea i r n 2c a s ew i mw e a ks h o c k l ，a v e ，t h es l o we v o m o nc a nb e i l c s c r i b e db yt h el i n e a rm e o r y k e yw o r d s ：b r i d 鲫a 1 1 觚吣c o n f i 毋删i o n ，m u l t i s t a t e ，h y d r o d 跚吼i c 赋) i l 咄 n o n a ) 【i 跚姗e 仃y 慰c h 哑y c r - m e s l l l ( 0 vi n s 协b i l 啵 s i l l 西em o d e ， n o i d i a rp e r 心b a t i o n m 论文原创性和授权使用声明 本人声明所呈交的学位论文，是本人在导师指导下进行研究工 作所取得的成果。除已特别加以标注和致谢的地方外，论文中不包 含任何他人已经发表或撰写过的研究成果。与我二同工作的同志对 本研究所做的贡献均已在论文中作了明确的说明。 本人授权中国科学技术大学拥有学位论文的部分使用权，即： 学校有权按有关规定向国家有关部门或机构送交论文的复印件和电 子版，允许论文被查阅和借阅，可以将学位论文编入有关数据库进 行检索，可以采用影印、缩印或扫描等复制手段保存、汇编学位论 文。 保密的学位论文在解密后也遵守此规定。 作者签名： 年月日 第l 章绪论 1 2 圆筒热对流的研究现状 目前晶体生长热对流不稳定性理论研究的大部分都是集中在圆柱内 r a y l e i g h b e n a r d 流动条件下的不稳定性( r b i ) 。r b i 热对流加热条件是：上壁面 低温，下壁面高温，侧壁绝热，图卜1 。圆筒内r a y l e i g h b e n a r d 对流的多解现 象广泛存在于实验【5 2 ，3 2 ，5 3 ，5 4 】当中，并成为各种理论分析【5 l ，5 5 ，】和数值模 拟【5 7 ，3 3 ，5 8 ，3 4 ，3 5 ，5 9 ，6 0 】的研究对象。线性稳定性分析【3 4 ，5 7 ，5 8 】已经给出了 圆筒内热平衡静止流体的第一次失稳分岔的情况，临界r a y l e i g h 数是尺寸比的 函数，随着尺寸比的增加，逐渐趋近于1 7 0 8 这个固定值。 图卜l ：圆柱r a y l e i g h - b e 腿r d 热对流条件示意图 当流动第一次分岔出现热对流后，随着r a y l e i g l l 数的增加，流体运动变得 高度非线性，并且强烈依赖于圆筒的尺寸比和p m i l d t l 数。目前的研究表明，即 使是弱非线性理论也很难准确预测第二次分岔的位置和情形【5 8 ，3 4 ，3 5 】。h o f 等 人【6 l 】首先在实验中发现中等尺寸比( 尺寸比接近1 ) 时也存在着非常类似于小尺 寸比时的多解现象。图卜2 为文 6 1 的实验结果，在肋= 1 4 2 0 0 时，实验给出了 同时存在的8 种不同定常解：包括轴对称解( a x i s y 咖e t r i c ) ，一个三涡解( t h r e e r 0 1 1 s ) ，一对相反的两涡解( t w or o l l s ) ，一对相反的四涡解( f o u rr o l l s ) ，一 对相反的三辐条( s p o k e ) 模式解。 2 第1 章绪论 f 王g 1 s h a d o w g r a p hi m a g e so fc o n v e c t i o np a t t e m s i n( a ) 一( h )月 = 142 0 0 ( a ) t h r e cm l l s ：( b ) t w or o l j sw i t hh o tf 1 1 l i dr i s i n ga l o n gt l l cc e n t e r ： ( c ) t 、v or o l l sw i t hc o i df i u i df a l l i n ga i o n gt h ec e n t e r ：( d ) f o u rr 0 1 l sw i t l lc o l d i l u i df a l l i n ga l 锄gt i l ec e n t e r ；( e ) f o u rr o i i sw i 1h o tn u i dr i s i n ga l o n gt i l e c e n t e e ( dt h r e es p o k ep a t t e mw i t hc o l df l u i df a l i i n ga l o n gt h es p o k e s ；( g ) m r e es p o k ep a n e mw i t h + h o tf l u i dr i s i n ga l o n gt 1 1 es p o k e ：( h ) a x i s 蛐m e t r i c p a t t e mw i t l lh o tn u i df i s i n j 乏 i nt 1 1 e c e n t e r ；( i )r o t a t i n gp a t t e ma tr = 2 66 0 0 ：( i ) p u l s i n gs p o k e dp a t t e ma tr = 3 30 0 0 图卜2 ：圆筒内r a y l e i g h b e n a r d 热对流的多解现象实验 分岔理论 5 1 是研究类似r a y l e i g h b e n a r d 对流等非线性复杂系统的强有力 的方法。如果流动具有对称性的话，群论分析也是一个非重重要的辅助分析工具。 为了求得所有的稳定和不稳定解，不仅需要局部的线性稳定性分析，求解特征值 和特征向量来判断流动的稳定性和发展趋势；而且需要整体的数值延续方法，获 得整个分岔分支的发展演化情况。 而更复杂，但是更接近实际情况的加热条件为b r i d g m a ng r o w t h c o n f i g u r a t i o n 条件( b g c ) ：上下壁面为低温，对应于晶体的熔点；。侧壁为抛物 型温度分布，对应于熔体加热条件，图卜3 。热对流的稳定性同边界条件是密切 相关的，r b i 的失稳振荡机理与b g c 是不相同的：与金属熔体对应的小p r a n d t l 数的r b i 的第一次失稳是从静止分叉到轴对称或非轴对称的定常流动，第二次分 叉才可能出现周期振荡；而b g c 由于侧壁加热条件始终存在流动，在小尺度比( a = 高度h 半径r ) 时流动由于动力不稳定性直接失稳振荡，在大尺度比时才表现为 r a y l e i 曲一b e n a r d 不稳定性起主导作用，需要经过定常的对称性破缺分叉后才出 现周期流动，而且我们发现中间尺寸比时两种机制相互作用可能有复杂的迟滞现 象出现。 到目前为止，对于b g c 条件下圆柱热对流问题的研究不多。而且大多都集中 于使用基于周向模态分解进行线性稳定性分析的方法对轴对称流的问题进行研 究。对于b g c 条件下轴对称热对流的问题，文献 2 0 3 6 中发现当某些特定参数， 比如r a y l e i g h ，超过特定值的时候，对流轴对称出现分叉或者流动过度到振荡 第l 章绪论 流动。而r a 数的临界值强烈的依赖于圆柱的尺度比。 图1 _ 3 ：圆柱b g c 热对流条件示意图 1 3 交界面问题的意义 多流体界面( 或流动) 问题几十年来一直是计算流体力学领域的热点问题，特 别是考虑互不掺混的多流体界面问题的高精度计算方法。简单的来说当两种理 想，互不混合的流体共处在空间某一区域中时，它们之间会形成一个两流体的交 界面，两侧的流体具有不同的物理和热力学性质，人们感兴趣的是这个界面在 外力如激波、重力等驱动下的演化情形，这就是一个可压缩多流体流动的简单模 型。交界面问题具有重要的理论意义和广泛的应用背景。一方面，交界面的问题 涉及到高精度的数值算法的发展。早期的l a g r a n g e 方法精度都比较低，而且难于 处理交界面的复杂演化，特别是当交界面发生融合或破裂等拓扑结构剧烈变化的 情况。现代的方法一方面可以处理复杂交界面的演化，同时也更注重高精度地捕 捉交界面的位置。交界面上一般存在密度间断等数学上比较奇异的问题，如何准 确地模拟交界面间断一直是近几十年来研究的重要课题之一。另一方面，交界面 的问题也涉及到广泛的实际应用背景。比如r a y l e i g h t a y l o r 不稳定性和 r i c h t m y e r m e s h k o v 不稳定性对于湍流机理的研究有这重要的作用。一个典型的 多流体跚i 问题是激波与界面的相互作用。跚i 在天体物理中的超新星的爆炸和惯 4 第l 章绪论 性约束核聚变中都起着重要的作用，研究它不仅可以帮助人们解释一些自然现象 而且有助于人们对于核爆炸机理更深入的理解进而实现对它的控制。同时，在各 种工业实际应用中也存在着广泛的运动界面，在晶体生长中枝晶的形成，演化和 发展，等等都会遇到多流体界面问题。科学研究的动力来源于实际的需求，研 究的成果反过来又推动着实践向着更高的水平发展，可压缩多流体界面问题的研 究工作同样遵循这样一个普适的规律，人们研究它的最终目的是揭示这种现象的 规律以利于更好地处理实际问题。 1 4 交界面问题的研究现状 近几十年来，交界面的数值模拟由于其深刻的现实意义和广阔的应用前景而 被广泛的研究。主要的交界面数值模拟方法可以分为两大类：l a g r a n g e 界面跟 踪方法( f r o n t t r a c k i n g ) 和e u l e r 界面捕捉方法( f r o n t c a p t u r i n g ) 。前者能准 确跟踪交界面的演化，间断处理无数值耗散，但需要显式地在交界面上处理界面 跳跃条件，难于处理界面拓扑结构的变化( 如界面融合和破碎) ，且扩展到高维情 况时算法非常复杂。g l i 咖等在这方面做出了系统的工作 4 2 4 5 ，主要是进行 可压缩流动如r a y l e i g h t a y l o r 和r i c h t m y e r - m e s h k o v 不稳定性的数值模拟，但 是算法非常复杂。后者计算简单，格式统一，易于处理界面拓扑改变，但不能精 确定位多介质交界面，不适合对界面位置要求很高的数值模拟。界面捕捉方法由 于容易实现得到了广泛的应用，早期的界面捕捉方法是借用多组分流体模型中基 于质量分数的方法来描述不同的流体，但是这样会导致在多流体交界面上的速度 和压力振荡 4 6 。为防止在界面附近的压力振荡而提出了众多的数值方法，例如 基于质量分数模型的修正方法 4 7 ，4 8 。p u c k e t t 等 4 9 将v o f 方法应用到多流 体可压缩的计算中来，而m i l l e r 等 5 0 给出了守恒的v o f 模型计算方法，并推 广到一般状态方程的情况。 为了更细致的刻画交界面的运动，l d e r 等 3 7 首先将l e v e 卜s e t 方法引 入到多介质可压缩流动中，但是该方法会在不同介质交界面附近产生振荡限制了 其的应用，f e d k i w 等 3 8 将界面追踪的思想引入到l e v e 卜s e t 方法中，提出了 g h o s t f l u i d 方法，利用接触界面两侧压力平衡和速度一致，并辅以等熵插值来 处理多介质在交界面上的内部边界条件，该方法可以避免交界面振荡的同时将界 面宽度控制在一个网格以内。近来，r i d e r 指出了数值方法的守恒性在多介质可 压缩流动模拟中的重要性，通过将各种数值模拟结果同实验结果的对比表明，在 大尺度统计参数上守恒计算方法的模拟结果同实验结果符合得更好 3 9 。 g l i 舳 4 0 指出传统的多介质计算方法都不是守恒的，并给出了一类守恒数值方 第l 章绪论 法，该方法结合f r o n t t r a c k i n g 方法能准确处理间断和f r o n t _ c a p t u r i n g 方法 计算界面简单的优点，较传统数值方法在界面附近提高了一阶精度；l i u 3 4 采 用类似的思想结合l e v e 卜s e t 方法也给出了相近的结果。 1 5 本文的主要工作内容 本文的工作主要包含b g c 圆筒热对流和交界面r m 不稳定性两部分。对于晶 体生长b g c 条件下加热模型的热对流不稳定性进行系统研究，采用动力系统分叉 分析和直接数值模拟相结合的方法来研究流动稳定性，给出相应的稳定性曲线和 分叉路径，得到失稳机理，探讨其中可能的非线性现象。由于本文采用的是真正 全三维的稳定性分析，同现有广泛使用的基于周向模态分解进行线性稳定性分析 不同，不仅能得到轴对称基本流线性稳定性分析的结果，还可以考虑三维非轴对 称分叉问题，预测和解释迟滞现象。因此，我们的研究得到了较多的新结果。此 外，对于交界面上r i c h t m y e r m e s h k o v 不稳定性问题，本文采用直接数值模拟的 方法，得到交界面随时间的演化过程，流场流动的图像，与实验结果和理论模型 预测相一致。 本文内容安排如下： 第一章，介绍圆筒热对流和交界面问题研究意义和现状 第二章，介绍使用基于柱坐标下的高精度有限差分方法，利用柱坐标下三 维不可压流动的改进投影算法，求解定常解，用a r n 0 1 d i 方法求解 以定常解为基本流的线性小扰动n s 方程的j a c o b i a n 矩阵的特征 值和对应特征向量的整体稳定性分析方法，并进行数值验证。 第三章，用第二章叙述的方法对不同g r a s h o f 下的轴对称流和非轴对称流 进行直接数值模拟，从流场结构，速度变化，频谱分析，相图等 等讨论各种流动模式的流动特征，发现七种双模态共存的现象。 第四章，用l e v e 卜s e t 同w e n o - a c 计算格式相结合的方法数值模拟单模大 扰动r m 不稳定性问题。详细讨论界面随时间的演化，给出扰动增 长率的变化。 第五章，总结性的给出结论和展望。 6 第2 章b g c 数值方法和校核算例 第2 章b g c 数值方法和校核算例 本章主要考虑圆筒中热对流的数学模型，给出控制方程，采用在计算中引入 新变量g ，= ，甜，消除了柱坐标下三维控制方程的奇异性。并利用已有的整体稳 定性和分岔分析的数值工具，采用柱坐标下三维不可压流动的改进投影算法，使 用非均匀网格，进行直接数值模拟，最后使用a r n o l d i 方法进行稳定性分析，研 究侧壁加热的圆筒内热对流的稳定性和失稳振荡。 本章使用上述数值方法对圆筒内b ri d g 眦ng r o w t hc o n fig u r ati o n 轴对称热 对流采用两种不同的无量纲公式进行了验证性研究，考虑了不同收敛时间步长和 e p s 收敛条件对计算结果的影响，计算了其三个不同的临界g ，数，并与已有的数 值模拟计算结果进行比较，考核了计算程序的正确性。 2 1 控制方程，算法及网格 2 1 1 控制方程 由于本章主要考虑圆筒内热对流问题，在柱坐标下进行数值研究显得更为 方便。基于b o u s s 扭l e s q 假设的不可压流动的无量纲控制方程在柱坐标下可以写 成为矢量形式： v u = 0 等+ u v u = 一跏+ v 2 u + g ，铊 ( 2 _ 1 ) 西 。 、厶u 塑+ u v 秒= 土v 2 口 国 p r 对于柱坐标的三个方向( ，仍力，控制方程可以展开为连续性方程： 三掣+ 三兰+ 誓：o ( 2 - 2 ) ，务 ，a 矽 瑟 、。 三个动量方程： 第2 章b g c 数值方法和校核算例 能带来奇异性问题。v e r z i c c o o r l a n d i ( 1 9 9 6 ) 11 提出的柱坐标下三维不可压 流动的有限差分法成功地解决了这个问题。该算法建立在k i m m o i n ( 1 9 8 5 ) 1 2 和r a i m o i n ( 1 9 9 1 ) 1 3 的工作基础上，用分数步法( f r a c t i o n a 卜s t e p ) 和近似 分解技术( a p p r o x i m a t i o nf a c t o r i z a t i o nt e c h n i q u e ) 求解原始变量下的 n a v i e r s t o k e s 方程。他们采用了交错网格，并在计算中引入新变量g ，= 埘，代替 通常的轴向速度甜，( 在z 处轴向变量吼= 0 ) ，使得坐标奇异性得到了处理。 为了处理柱坐标( ，仍z ) 下n a v i e r s t o k e s 方程在对称轴，- = o 处的奇异性，引 入以下新变量代替速度分量，蚱，： = ，吼= m ，吼= “： 此时连续性方程如下： 冬+ 篓+ ，冬：o ( 2 啕 关于吼的动量方程如下： 鲁= 弓嚣+ 以导，誓) - 笋+ 吉等+ 争号器i z ) f ，a 矽l ，ia p a rj，2，2a 矽2a 2 2 ，a 矽l 鲁害+ 愕b 誓) + 专等+ 碧号等l 6 ， 历务l 毋i ，务r za 矿 瑟2 ，a 矽l 鲁= 一刮吾新誓) + 专等+ 参i + 函p 历如l ，勿勿，2a 矿瑟2 i 其中的非线性项均司以写成守恒形式： 一弛1 呱绋。l 嘲弛丝 一 -i d ta tp却r8 pa z 鲁暑鲁+ 昙+ 品( 芋) + 警一露7 ， 历研 务i ，ja 矽i ，j瑟 叩 ” 观一却：。1a 甄吼1 却伊吼。a g ； 励一街。，西厂a 伊。如 同时温度方程可以写为如下形式： 脚一a 口。1a g ，秒1a 秒。a 吼秒 lflaf ，a 口、la ；a z 秒 ( 2 8 ) 2 百l7 石l ，石j + 7 可+ 萨i p r l ，勿l 跏，2a 矿勿2l 9 第2 章b g c 数值方法和校核算例 2 1 3 分数步法 本文中的算法是采用v e r z i c c o 0 r l a n d i ( 1 9 9 6 ) 提出的柱坐标下三维不可 压流动的有限差分法，使用混合的三步r u n g e k u t t a c r a n k n i c h 0 1 s o n 。o r a n l d i 提出的这个柱坐标下的有限差分方法精度好效率高，是非常有效的直接数值模拟 方法。这个算法建立在k i m m o i n 和r a i n m o i n 的工作基础上，用分数步法 ( f r a c t i o n a 卜s t e p ) 和近似分解技术( a p ”o x i 腿t i o nf a c t o r i z a t i o nt e c h n i q u e ) 求解原始变量下的n a v ie r s t o k e s 方程。 第一步： 譬= 陋护咖c 厶蛳4 ，华 酽一钟= 一蛹1 。 ( 2 _ 9 ) d ( 酊) = 0 第二步： 譬= 陋咖蒯护咖2 州向+ 4 ，+ 4 ) 华 钟一岔= 哆f q 2 ( 2 1 0 ) d ( 钟) = o 第三步： 譬= 陋枷蒯护卿3 州如+ 4 圳华 矿1 一钟= 一妈3 ( 2 一1 1 ) d ( 1 ) = 0 其中e 包括了非线性项和非交叉导数的简单速度导数的粘性项，如，4 ，4 ，q 分 别表示粘性和压力项的离散微分关系，( ，= l ，2 ，3 ) 是与压力有关的标量。d 是散 度算子： d = 吾嘉+ 吾砉+ 砉 c 2 一- 2 ， r6 r6 r6 z 嘶，乃，岛是时间推进r u n g e k u t t a 法的系数( s p a l a r te 舌a j ( 1 9 9 l 1 4 ) ) ： l o 第2 章b g c 数值方法和校核算例 = 8 1 5 ，乃= 8 1 5 = 2 1 5 ，儿= 5 1 2 ，仍= 一1 7 ，6 0 嘞= 1 3 ，乃= 3 4 ，岛= _ 5 1 2 此方法对粘性项的时间离散采用了二阶精度的c r a n k n i c h o l s o n 格式，对流项采 用了三阶精度的r u n g e k u t t a 法。空间离散均在交错网格上采用中心差分。所以， 整个算法是时空二阶精度的。 利用( 2 9 ) ，( 2 1 0 ) 和( 2 1 1 ) 第一个式子计算中间速度场毒：时，可以把这三 个公式写成： ”觚p 一铲磐馗一乡， ( 2 _ 1 3 ) = 出【乃叫+ 岛叫一嘶q + 嘶( 如+ 4 + 以) 彰】 式中届= 嘶，2 上式左边的求逆是比较耗时的计算，对于定常问题需要大步长计算时，可以 利用周向均匀网格周期性边界条件采用f f t 离散，然后得到可分离变量的椭圆型 方程应用f i s h p a c k 软件包直接快速求解；对于非定常问题，由于时间步长比较小， 采用近似分解技术讲上式左边的大型稀疏矩阵求逆运算转换为三个三对角矩阵 的求逆，误差只有d ( f 3 ) ： ( 1 一届如： ( 1 - 届苎x l 一届，( 参一咖， ( 2 _ 1 4 ) = m 乃叫+ 岛叫一q q + 嘶( 如+ 4 ，+ 4 ) 彰1 用( 2 9 ) ，( 2 _ 1 0 ) 和( 2 1 1 ) 第二个式子求解速度场酊，口? ，g ：l “时，先结合连续 性方程得到关于的p o i s s o n 方程： 三= 毛脚7 鲁q ( 五，左，五) ( 2 一1 5 ) 口l ， 。 其中 三亭鲁七砉厂砉皤 埘 三2 7 万+ 7 石厂石+ 虿 ( 2 一1 6 ) 工，左，五分别表示( 伊，- ，z ) 方向的网格结点。计算( 2 1 5 ) 是在的7 时，对有限差分 方程可以在p 和z 方向用离散f o u 褥罚薛蓉羹鍪翼矽囊篓 2 1 7) 其中五= l ，2 ，l ，五= 1 ，2 ，2 ，五= l ，2 ，3 。将( 2 一1 7 ) 式带入( 2 _ 1 5 ) 式，得到：e 熹，熹+ 导一多以垂( ( 2-18) ( 7 石 第2 章b g c 数值方法和校核算例 方程组的直接求解还是可以接受的。 2 压力预估步：在投影算法的速度预估步之前，首先通过压力p o i s 湖方程 估一个压力值p 肿1 。根据连续性方程得到的压力p o i s s o n 方程为： 厶痧肿1 = v 卜 “) + v 2 矿+ ，州】 i nq ( 2 2 3 ) 其边界条件 等一卜华州m 1 2 4 ， 由于采用交错网格，这时求解实际上并不需要压力边界条件。 3 速度预估步：这时的计算同投影算法一致，通过动量方程樽到一个预估 速度，这个速度的散度不为零， 等+ 脚协一哕+ 1 + v 2 玉+ 1 i nq ( 2 - 2 5 ) 速度边界条件采用 番= 甜l 舳 o nq ( 2 2 6 ) 材2 甜l a 晓 o n 2 【z 一2 6 ) 4 速度校正步：将上面计算得到的中间速度投影到最终的无散速度场， 掣：- v ( p 槲叫+ 1 ) ( 2 - 2 7 ) 一= 一v i 乃一n - l7 7 ，i f 一 v 甜肿1 = 0 ( 2 2 8 ) 为了计算速度，需要定义一个中间变量：矽= 缸( p 斛1 一多肘1 ) ，根据上述两个 方程呢个可以得到一个p o i s s o n 方程 = v 疗 i nq ( 2 2 9 ) 边界条件采用n e u 腿n n 边界条件 掣：o 。nq ( 2 3 0 ) 最后根据中间变量矽计算得到最终的无散速度场： 材打+ 1 = 五一v ( 2 3 1 ) 2 1 5 非均匀网格 为了快速而准确的求解p o i s s o n 方程和h e l l i l h o l t z 方程，数值离散时在周向缈 对n a v i e r s t o k e s 方程的定常解“进行稳定性分析，其加上小扰动后的 性稳定性方程 相应的特征方程 色= ( - + 三) “ ( 虬+ 三) 甜= 知 ( 2 3 2 ( 2 3 3 这里a 为彳= ( 玑+ 三) 的特征值。可以根据实部最大的特征值丸来判断定常解 稳定性。当丸实部为负时流动稳定；当丸实部为正且虚部为0 ，流动不稳定 加上小扰动后会收敛与另一个定常解；当丸实部为正且虚部不为0 ，最终的稳 流动是一个非定常振荡流，其虚部表征初始振荡频率。考虑特征值还有其他重 意义，比如分叉点附近的正实部特征值数目可以揭示流动的多次分叉路径等等 系统( 2 3 3 ) 的所有特征值以及特征向量可以用q r 方法求解，但对于复杂流动 析中比较大的矩阵，q r 方法的计算量始非常大的。而且实际上稳定性分析并 需要所有的特征值，而是最主要的几个特征值对及其相应的特征向量部。在计 中通常要求8 1 6 个特征值，当前一般求解算法采用a r n o l d i 方法或者正交子 间迭代方法。在本文的稳定性分析计算中，求解前1 6 个最优特征值进行分析 线化扰动方程( 2 3 2 ) 的解可以写成为如下形 甜o + f ) = e x p ( f ( 虬+ 上) ) ( ， 当缸l 时，有以下线化方程近 o + 出) = “o ) + 缸( ，一圮) 一( 玑+ 三) “ = “( f ) + ，( 虬+ 三) “( f ( 2 3 5 兰矿- + 工( ， 由于这个近似将j a c o b i a n 矩阵么的最大实部特征值( 对) 变换为e x p ( 纠) 的 大模特征值( 对) ，可以采用a r n o l d i 方法来计算这些优势特征值。可以使a r n 0 1 d 方法中的每一步迭代相当于求解一次线化的n a v i e r s t o k e s 方程 基本的a r n 0 1 d i 方法是向l ( r y l o v 子空间进行正交投影的不变子空间法 1 5 该方法可简述如下，给定单位1 1 i | ：范数的初始向量嵋尺”，逐列的形成矩 第2 章b g c 数值方法和校核算例 伯格矩阵。这个过程可概括为a r n 0 1 d i 分解 彳圪= 圪皿+ 魄“j 蚝+ l 彳 ( 2 3 6 ) 圪的列构成了k r y l o v 子空间的妒册 甜i ，彳，彳2 地，4 h 甜1 ) 的一组标准正交基。 如果噍小= 0 ，则它是彳的不变子空间，风的所有特征值是么的特征值的子集， 即a ( 巩) 名似) 。若噍+ l j o ，则根据上式可得 出一触= 噍+ l j + l y ( 2 3 7 ) 式中，y 为巩的单位特征向量，即巩y = 砂，x = 圪j ，。称五为r i t z 值，x 为相 应的r i t z 向量。显然 i i 出一允x l = i 阪+ ，j + ，蠢y 0 = 魄+ 。l y i ( 2 3 8 ) 基本的a r n 0 1 d i 方法如下： 1 ) 给定j ，卜约| | 1 1 2 2 ) f o r ，= 1 ，2 ，七 国净彳吩 f o rf = 1 ，2 ，七 = l f 国；国车国一； e n d 吩“j = 9 力0 2 ：吩+ l = 缈以“_ ，； e n d 3 ) 计算矩阵风= ) 的特征对( 丑，乃) 及彳的关于丑的r i t z 向量而= 圪乃。 风中包含的特征信息主要取决于初始向量，在一般情况下，我们想要的 特征信息可能在七很大时才出现。然而，当后很大的时候，不但需要巨大的存储 空间和计算量，而且保持k 中各列向量正交也相当困难。为此，；需要对基本 a r n 0 1 d i 方法进行改进，使得在较小的固定七维子空间中包含我们想要的，( 厂 o ，虚部 为o ，所以有非轴对称基本流出现。 口 g 芦5 2 1 0 0 o 4 g f 5 2 1 1 0 e g 卢5 2 1 2 0 - g 声5 2 1 3 0 g 户5 2 1 4 0 g 芦5 2 1 5 0 0 2 ( ) g 芦5 2 2 0 0 mi ii _ - 。 ，r j 一 rv o 。 l r e o 0 0 5o - 0 0 0 5 i 2 n 4= - 。t 。 图3 - 3 ：不同g ，数下的流动特征值谱 最后，对不同函数下最优特征值的实部进行插值，如图3 4 。得到第一次 轴对称临界g ，数g ，= 5 2 1 4 1 。 二：二彘5 二么 。 r e 2 05 2 1 2 55 2 1 3 0 5 2 1 3 5 5 2 1 4 0 5 2 1 4 55 2 1 图3 - 4 ：通过对最优特征值实部进行插值得到临界g r 数 第3 章b g c 圆筒热对流中的多稳定解现象 3 2 轴对称定常流和非轴对称定常流 当g r a s h o f 数刚开始大于轴对称基本流第一次分叉的临界g ，值的时候， 轴对称的定常解是不稳定的。这个时候通过对轴对称基本流加扰动，然后进行 非定常迭代，我们可以得到非轴对称基本流的定常解。由于g ，数还比较小，所 以这个时候非轴对称的定常解是稳定的。而且，一开始的时候非轴对称基本流 的流场，温度场，速度场偏离轴对称的情况是比较小的。但是，随着g ，数的慢 慢增加，非轴对称的基本流的流场，温度场，速度场开始越来越偏离轴对称状 态，偏离的幅度越来越大。通过对其稳定性分析得到的前几个最优特征值实部 的结果也表明其基本流越来越趋于不稳定的状态。 ( = 8 5 0 0 0 由于在相当大的一段g ，数范围之内，轴对称基本流的流场，温度场，速度 场变化不大，所以在这里我们只给出当西= 8 5 0 0 0 的时候的计算结果。表3 3 给出g ，= 8 5 0 0 0 的时候，通过对其定常解求解线化小扰动方程的j a c o b i a n 矩阵 得到的前两个优势特征值。这两个优势特征值的实部大于0 ，说明此时轴对称 定常解是不稳定的。 表3 3 ：g ，= 8 5 0 0 0 情况下轴对称定常解的前2 个优势特征值 r e q ) ，增长率i m ( 力) ，线性圆频率 13 4 0 5 0 3 1 4 2 9 1 9 7 2 7 9 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 23 4 0 5 0 3 1 4 1 4 9 2 9 1 1 10 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 第3 章b g c 圆筒热对流中的多稳定解现象 r ：1r = or = l 图3 _ 5 ：g r = 8 5 0 0 0 轴对称定常流子午面的流线图( 左) 和等温线图( 右) ( b ) z = h 2 第3 章b g c 圆筒热对流中的多稳定解现象 ( c )z = 3 h 4 图3 - 6 ：g ，= 8 5 0 0 0 轴对称定常流横截面的流线图( 左) ，等温线图( 中) ，轴向速度等值 线( 右) ，虚线为负值。( a ) z - h 4 ( b ) z = h 2 ( c ) z = 3 h 4 图3 5 的等温线并不是上下对称的，而是在上半部分相对稀疏，下半部分 相对密集，我们所给外部温度场是上下对称的，所以有可能是由于重力场作用 所导致的。子午面的流线场也有相似的现象，其下半部分的流线比较密集。图 3 6 的轴向速度等值线表明，流体在中间部分沿着轴线向下运动，然后从四周 的环线向上运动。在高度z = h 4 的横截面上的流线图表明，流体从中心点向外 呈辐射状直线扩散，但是并没有到达壁面。在距离壁面相当一段距离内。在横 截面上流体几乎是静止的，这个时候内部流体并没有通过流动与壁面进行热交 换，其所得到的能量主要是通过轴向方向上的流体交换以及热传导。这样直接 导致了温度在壁面附近变化的比较大，而在圆心内部附近比较小。在圆柱的中 上部分z = h 2 ，z = 3 h 4 ，横截面内流体是从外部向内部流动，进行热交换。 g ，= 5 5 0 0 0 表3 4 给出g ，= 5 5 0 0 0 的时候，通过对其定常解求解线化小扰动方程的 j a c o b i a n 矩阵得到的前三个优势特征值。这三个优势特征值的实部小于o ，说 明此时非轴对称定常流是稳定的。 第3 章b g c 圆筒热对流中的多稳定解现象 表3 - 4 ：g ，= 5 5 0 0 0 情况下非轴对称定常解的前3 个优势特征值 r e ( 名) ，增长率i m ( 名) ，线性圆频率 11 6 4 4 2 5 0 1 7 1 7 6 8 1 0 8 6 e 0 0 4o 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 7 2 0 4 0 4 6 35 4 5 9 6 8 6 92 4 3 5 2 2 7 7 3 5 7 4 3 3 4 7 3- 7 2 0 4 0 4 6 3 5 4 5 9 6 8 6 9- 2 4 3 5 2 2 7 7 3 5 7 4 3 3 4 7 图3 7 ：g ，= 5 5 0 0 0 非轴对称定常流子午面的流线图( 左) 和等温线图( 右) ( a ) z h 4 3 3 第3 章b g c 圆筒热对流中的多稳定解现象 ( b ) z 爿i 2 图3 - 8 ：g ，= 5 5 0 0 0 非轴对称定常流横截面的流线图( 左) ，等温线图( 中) ，轴向速度等 值线( 右) ，虚线为负值。( a ) z - h 4 ( b ) z = h 2 ( c ) z = 3 h 4 从图3 7 我们可以看出，在g ，不是很大的时候，g ，= 5 5 0 0 0 ，非轴对称定 常流子午面上的等温线和流线场，相对于中心轴还是对称的，其非轴对称现象 并不是很明显。而在圆柱下部分的横截面，如图3 8 a 所示，其温度场出现明显 的非轴对称性，在中心部分开始出现“心型的等温线。此时，轴向速度等值 线的中心点开始偏移圆心的位置，向壁面靠拢。流线图也不是像轴对称流那样 呈辐射状直线扩散，而是开始出现弯曲。有的流线弯曲1 8 0 度，然后到达壁面； 有的流线弯曲9 0 度，然后沿着半径直线或向内流向圆点，或向外流向壁面。整 个流线图开始出现比较复杂的图形。而在圆柱的中上部分的横截面上，如图 3 8 b ，c 所示，其等温线呈现准圆形，但是中心点的位置开始偏离圆柱的圆心。 而轴向速度等值线图的非轴对称现象就比较明显了。整个轴向速度负等值线开 始向边上偏移。而在相反方向上，轴向速度正等值线开始出现单连通域。这时 的流动模式是一种单涡结构。在中间偏侧部分向下运动，然后从另一侧的壁面 上向上运动。且在越靠上部分，非轴对称现象越明显，其单涡的结构也越大。 与此同时，横截面上的流线图也开始复杂，并不像轴对称流时沿着半径向圆心 第3 章b g c 圆筒热对流中的多稳定解现象 流动。有的流线直接从壁面流向圆心；有的流线绕过圆心后，或者沿着半径流 向圆心，或者沿着半径流向壁面。在z = h 2 的横截面上的流线图可以看到，在 壁面的对称位置上，有一个点源和一个点汇。从点源流出的流线大部分流向圆 心，少部分流向对称位置上的点汇。而在z = 3 h 4 的横截面上的流线图可以看到， 在靠近圆心的部分出现两个对称的小涡。 g ，= 8 5 0 0 0 表3 - 5 ：g r = 8 5 0 0 0 情况下非轴对称定常解的前3 个优势特征值 i k ( 五) ，增长率h ( 名) ，线性圆频率 11 9 0 5 8 1 9 8 3 7 0 9 5 5 6 1 9 e 0 0 50 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 29 7 6 3 6 8 5 3 6 9 5 3 4 3 8 2 18 5 5 3 0 4 5 4 6 2 3 4 3 9 7 39 7 6 3 6 8 5 3 6 9 5 3 4 3 8 218 5 5 3 0 4 5 4 6 2 3 4 3 多7 图3 - 9 ：g r = 8 5 0 0 0 非轴对称定常流子午面的流线图( 左) 和等温线图( 右) 3 5 ( b ) z - h 2 ( c )z = 3 h 4 图3 一1 0 ：g ，= 8 5 0 0 0 非轴对称定常流横截面的流线图( 左) ，等温线图( 中) ，轴向速度等 值线( 右) ，虚线为负值。( a ) z - h 4 ( b ) z - h 2 ( c ) z = 3 h 4 从图3 9 我们可以看出，在g ，达到一定数值的时候，g ，= 8 5 0 0 0 ，非轴对 称定常流子午面上的等温线和流线的非轴对称现象就比较明显。子午面上的等 温线开始出现一定程度上的扭曲，像被挤压的“x ”。而流线图也不像在轴对称 的时候两边是相互独立对称的。这个时候左右两边的流线是相互相交的，左边 的流线通过中心轴流向右边，反之亦然。而且，此时流线图并不是关于中心轴 第3 章b g c 圆筒热对流中的多稳定解现象 对称的。左边子午面的流线有两个涡，在中下部，有一个比较大的涡，其中心 大约在z = h 5 处。而且，此大涡的流线与中心轴右边子午面涡的流线是相连的。 在左半边的子午面的左上角，有个小涡的存在。此时，在圆柱下部分的横截面， 如图3 1 0 a 所示，其温度场在中心部分出现两个涡，一个大，一个小，两个涡 大小相差比较大。而此横截面上轴向速度等值线出现两涡结构。流体在中间部 分向下运动，然后从两个相反方向的侧