（流体力学专业论文）基于预处理技术的微可压缩液体流场模拟方法.pdf

i ：海人学硕上学位论文 摘要 一般情况下，液体流动可近似作不可压缩流场处理。压力修正法是求解不 可压缩流场的最常用的方法；也有学者提出采用预处理技术，利用时间推进法 进行不可压流场的求解。但是真正的不可压流体是不存在的，如果将液体工质 看作可压缩流体，采用常规的时间推进法求解这类问题，一般很难得到收敛解。 其原因在于水的音速相当高，般流体机械内的流速相对其音速很小，方程组 的对流项系数矩阵具有很强的刚性。通过引入预处理矩阵，改变方程组的特征 值，使对流项系数矩阵的特征值处于同一量级，从而可以达到消除刚性与加速 计算收敛性的目的。 本文以水为工质，考虑低速流动时水的微压缩性，通过以密度和压力为原 始变量的预处理技术，采用时间推进方法，进行流场模拟。利用水和蒸汽的热 力学性质状态参数关系，由压力和密度求解流场模拟过程中所需的其他状态参 数，以封闭可压缩的流体力学方程组。对水中的n a c a 0 0 15 翼型绕流、低雷诺 数的圆柱绕流和n a c a 0 0 1 2 翼型振荡绕流问题的模拟表明：本计算方法能够较 好地模拟水中的低速流动问题。 关键词：微可压缩 预处理技术流场模拟 v 1 - i 1 1 人学颀上学化论文 a b s t r a c t l i q u i df l o wi so f t e nr e g a r da si n c o m p r e s s i b l ea p p r o x i m a t e l y t h ep r e s s u r e c o r r e c t i o nm e t h o di sw i d e l yu s e dt os i m u l a t et h i sk i n do ff l o w f ie l d h o w e v e r ， t h ep r e c o n d i t i o n i n gt e c h n o l o g yc a na l s ob eu s e dt os i m u l a t et h ei n c o m p r e s s i b l e f i o w i nf a c t ，t h ea b s o l u t e l yi n c o m p r e s s i b l ef l u i dd o e sn o te x i s t a 1t h o u g h t h e f l o w i n gli q u i dc a nb er e g a r d e da sc o m p r e s s i b l e ，t h ec o m p u t a t i o nu s u a l l y d o e sn o tc o n v e r g ew h e nt h et i m e m a r c h i n gm e t h o di su s e d b e c a u s et h ev e l o c i t y o fs o u n di nw a t e ri sm u c hh i g h e rt h a nt h ev e l o c i t yo ff l u i di nt h em a c h i n e r y ， t h ec o e f f i c i e n t m a t r i xo fc o n v e c t i o nt e r mi sv e r ys t i f f b yi n d u c i n ga p r e c o n d it i o n i n gm a t r i x ，t h em a g n i t u d eo f t h ee i g e n v a l u e so f t h ec o e f f i c i e n t m a t r i xo fc o n v e c t i o nt e r mc a nb em o d i f i e dt ob ei n t h es a m eo r d e r ，s ot h e s t i f f n e s so ft h ee q u a t i o n si se l i m i h a t e da n dt h et i m e - m a r c h i n gm e t h o dw o r k s i nt h i sd i s s e r t a t i o n ，w a t e ri su s e da sw o r k i n gf l u i d ，ap r e c o n d i t i o n i n gm e t h o d u s i n gd e n s i t y ，v e l o c i t ya n dp r e s s u r ea st h ep r i m i t i v ev a r i a b l e si sa d o p t e dt o s i m u l a t et h ef l o w i no r d e rt oe n c l o s et h ee q u a t i o ns y s t e m ，t h er e l a t i o n s h i p s b e t w e e nd e n s i t y ，p r e s s u r ea n do t h e rs t a t ef u n c t i o n sa r ei n d u c t e db yu s i n gt h e t h e r m o d y n a m i c sp r o p e r t yo fw a t e ra n ds t e a m t h r e et e s tc a s e s ，t h ef l o wa r o u n d n a c a 0 0 1 5a i r f o i1i nw a t e r ，t h el o wr e y n o l d sn u m b e rf l o wa r o u n dc y l i n d e ra n d t h ef l o wa r o u n do s c i l l a t i n gn a c a 0 0 1 2a i r f o i l ，v a li d a t et h ec a p a b i l i t yo ft h i s m e t h o dt os i m u l a t et h el i q u i df l o w k e y w o r d s ：s l i g h tc o m p r e s s i b i l i t y , p r e c o n d i t i o n i n gt e c h n o l o g y , s i m u l a t i o no ff l o w f i e l d v l i ：海人学硕l p 位论文 原创性声明 本人声明：所呈交的论文是本人在导师指导下进行的研究工作。 除了文中特别加以标注和致谢的地方外，论文中不包含其他人已发 表或撰写过的研究成果。参与同一工作的其他同志对本研究所做的 任何贡献均己在论文中作了明确的说明并表示了谢意。 签名：丝逛日期：丝壁：! 兰： 本论文使用授权说明 本人完全了解上海大学有关保留、使用学位论文的规定，即： 学校有权保留论文及送交论文复印件，允许论文被查阅和借阅；学 校可以公布论文的全部或部分内容。 ( 保密的论文在解密后应遵守此规定) 签名：趋盔翮签名：瑾、瑶日期： op g f2 乙 ：海大学硕士学位论文 1 1 课题来源 第一章绪论 本课题来源于国家自然科学基金项目，项目编号：5 0 5 7 6 0 4 9 。 1 2 课题研究的目的和意义 本文通过采用预处理技术，求解可压缩的流体力学方程组，来计算以水为 工质的低速流动问题。利用水和蒸汽热力学性质图表，用插值的方法由压力和 密度求得温度等其他状态参数，用以封闭方程组。该方法旨在统一可压与不可 压流场的计算方法，最终能够模拟液体工质发生空化现象的复杂流动，并将其 应用于流体机械的流场模拟中，课题组正在进行相关研究。 因为流体都是可压缩的，所以液体流动及气体的低马赫数流动都具有一定 的压缩性，用预处理方法对可压缩n s 方程组进行求解比求解不可压缩n s 方程 组更符合流体流动的物理特性。此外，从低速到超音速很大范围内的流动，预处 理方法只需采用同一组控制方程描述，选择相同的计算格式，从而实现了可压与 不可压流场计算方法上的统一，该方法最终可以应用于可压与不可压、高速与低 速、气体与液体同时存在的复杂流场。 1 3 国内外研究概况 1 。3 1 国外研究概况 预处理方法是由b r i l e y ，m c d o n a l d 和s h m a r o t h 最先提出来的。通过对流 体力学方程组旌加一个预处理矩阵，总体上提高了计算过程的收敛速率。t u r k e l 基于原始变莓和基本的流体力学方程组，针对应用于可压缩方程的人工可压缩 方法，设计了一个普遍性的预处理矩阵12 | 。c h o i 与m e r k l e 将预处理技术用于粘 性流场的模拟之中1 3 j ，实践表明：经过预处理后的n s 方程组，计算收敛速率 提高了两个量级。w e i s s 与s m i t h 报道了一个适用于变密度低速流场求解的预处 l 海大学硕：f ：学位论文 理方法【4 】ob r i l e y ，t a y l o r 和w h i t f i e l d 提出了一个基于全马赫数流场的预处理 矩阵，该预处理矩阵形式简单，且易于实现低速与高速流动的自动切换p j 。 对于微可压缩液体的流场模拟，国外学者一般将微可压缩液体工况流场简化 为不可压流动，采用不可压缩n s 方程组进行求解。典型的解法有人工压缩性方 法、涡量流函数法和压力修上e 法。而采用可压缩n s 方程组求解微可压缩流动的 方法暂时没有报道。 1 3 2 国内研究概况 国内学者对预处理矩阵的构造比较少，黄典贵教授在2 0 0 5 年构造了一个基于 压力与焓为原始变量的预处理矩阵，该预处理矩阵适用于求解全马赫数流场1 6 j 。 大多数国内学者的研究，都是针对国外学者构造好的预处理矩阵进行直接应用。 比如：两北工业大学韩忠华、乔志德等比较了c h o i 预处理矩阵与t u r k e l 预处理矩 阵的功能1 7 j ，得出如下结论：( 1 ) 预处理方法能十分有效地改善低马赫数粘性流动 计算的收敛性和收敛精度；( 2 ) 预处理方法同样适用于高速流动，且有一定的收 敛效果；( 3 ) c h o i 预处理和t u r k e l 预处理对于低马赫数流动计算的性能基本相当； ( 4 ) 在较高马赫数计算中c h o i 预处理的稳定性和加速性能略优于t u r k e l 预处理。北 京航空航天大学严明副教授将文献【5 】的预处理矩阵应用于叶轮机械的三维数值 模拟，对于低压压气机算例的计算结果表明，采用预处理方法后可以得到与实验 结果吻合较好的数值模拟结果【8 】。中国空气动力发展研究中一1 1 , 肖中云、牟斌等【9 】， 将文献【4 】的预处理矩阵用于计算无粘b u m p 流和绕二维圆柱的非定常流动，得出结 论：数值结果表明预处理技术显著提高了低速流场求解的效率及求解精度，对所 模拟低速的非定常流动而言，所需代价只是传统时间推进法的几分之一。中国空 气动力发展研究中心肖中云、江雄等【1 ，将w e i s s s m i t h 预处理矩阵应用于悬停 旋翼的模拟中，得到如下结论：( 1 ) 在加入低速预处理方法以后，流场的收敛速 率得到了明显提高，在同等迭代步数下得到的下洗速度场更为准确。( 2 ) 通过尾 涡耗散的大小验证了低速预处理避开了因特征值差异带来的数值耗散过大的问 题，使可压缩方法具有在计算低速问题时保持精度的特点。 对于微可压缩液体的流场模拟，国内学者一般将微可压缩液体工况流场近似 为不可压流动，采用不可压n s 方程组进行求解。而采用可压缩n s 方程求解的 2 上海人学硕上学位论文 方法暂时也没有报道。 1 4 论文的主要研究内容 本论文是以作者攻读硕士学位期间承担课题的工作为基础，在第一章中阐 述了课题研究的来源、目的、意义以及国内外研究的现状。第二章阐述了预处 理方法及本文预处理矩阵的构造，顺带介绍下求解不可压缩流动的人工压缩性 方法、涡量流函数方法及压力修正法。第三章阐述了状态方程，如何利用水和 蒸汽性质表获得水的其余状态参数。第四章讲述了本文控制方程的数值方法， 采用时间推进法进行求解，采用r o e 格式离散对流项，因为控制方程的特征值 比较小而引入熵修正，通过限制器使离散格式精度提高。第五章讲述了三个算 例的数值模拟结果与分析。利用本文方法对n a c a 0 0 1 5 翼型绕流、低速圆柱绕 流及n a c a 0 0 1 2 翼型振荡绕流进行数值模拟，然后对得到的数值模拟结果与文 献中的实验结果或数值模拟结果进行对比分析。最后第六章总结全文。 j ：海人学顾十学位论文 第二章预处理方法 2 1 不可压缩流动解法 当流体运动速度比声波传播速度小很多时，流体的密度变化不大，可近似 为常数，此时描述流体运动规律的n s 方程组可以大为简化，称其为不可压n s 方程组，方程如下： 锄o vo w 一 + + = 0 a x 却 a z 塑4 - 甜塑+ v 塑+ w - - 锄：一! 望+ 闪：搿 甜+ v + = 一一土+ d v 。搿 o to x o y 0 z p 0 x ( 2 1 ) 尘+ z ，生+ v 堡4 - w 塑：一三望+ o v ：1 ，+ z ，+ v w = 一一二+ 。1 ， 魂a x 劫 娩 p 却 一0 w + 甜坐+ v 业+ w 丝：一上望+ o v 2 w 一十甜一+ v 一十w = 一一二+ 。w a t融 两0 zp0 z 其中u = 。 从形式上看不可压n s 方程组比可压n s 方程组较为简单，但足在求解上 引入了新的困难，由于祟+ 芸+ 譬：o 使方程组不再是进化型，难以采用己发 t 强 u y o z 展比较成熟的时间推进法求解。原来的进化型连续方程变成了在每个时间步上 对速度矢晕的约束条件。此外对压力缺少边界条件也给不可压n s 方程组的求 解带来一定的困难。针对这问题，学者们发展了许多处理方法，比如人工压 缩性方法12 1 、涡鼍流函数法、压力修正法 1 3 , 1 4 等等。 在低速情况下，如果把流体仍然看作是可压缩的，不做任何简化直接求解 可压缩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 5 】。下面分别介绍上述提到的各种方法。 4 ：海人学硕1 ：学位论文 2 1 1 人工压缩性方法【i 0 2 1 s 人工压缩性方法是把不可压缩流体看作一种虚拟的可压缩流体。在求解中 仍采用非定常的动量方程和经过人工非定常化的连续方程。这样处理的结果使 得方程组变为进化型，可以通过通常的求解进化型方程的方法进行求解。这一 方法最早由c h o r i n 和y a n e n k o 提出。 经过变型的连续方程为 塑+ c 2 弋7 莎：o a t ( 2 2 ) 这里c 2 是一任意的自由参数。在未达到定常状态之前该式没有任何意义。在达 到定常状态时雩磊= 0 ，即v 莎= 0 。其中参数c 2 的选取应使求解过程收敛和 加速收敛。 2 1 2 涡量一流函数法【6 1 通过引入涡晕方程和流函数，对二维不可压n s 方程组进行变换，写成涡量 流一函数的形式。对如下二维不可压n s 方程组： 丝+ 翌：0 叙 砂 塑+ 材塑+ v 塑：一至+ 上一v ：“ 一十材+ v 一= 一一二+ 一“ a ta x a v p8 x r e 堡+ 甜鱼+ v 鱼：一三至+ 上一vz p 十甜一十v = 一一二+ 一 。l ， a ， 出 砂p 砂 r e 引入流函数妒和涡量力： 8 p8 p z f = 二，v = 一二 砂 缸 ( 2 3 ) ( 2 。4 ) ( 2 5 ) ( 2 6 ) j ：海人学硕l ：学位论文 国：鱼一一o u ( 2 7 ) 国= 一一 t ， 苏 方 ( 2 4 ) 式对y 求导，( 2 5 ) 式对x 求导，然后消去压力项得到方程 丝+ z ，塑+ v 塑：上一v2 c _ o( 2 8 )+ z ，+ v = k 二o j o t融 却 r e 将( 2 6 ) 式代入式( 2 7 ) 得 v ：妒：磐+ 宴：一国 ( 2 9 ) 钡。 咖 方程( 2 8 ) 与( 2 9 ) 构成了涡量一流函数型的n s 方程组。 可用不同的方法对( 2 8 ) 与( 2 9 ) 式进行离散和求解。可以看出，方程组中不 再出现压力项，使得方程更为简单。当妒、国以及u 和v 求出后，可利用动量 方程或对应的压力p o s s i o n 方程求出压力。 虽然涡量一流函数方法在计算二维不可压缩n s 方程组中取得了成功，但是 把该方法推广到三维时却出现了困难。其原因是三维流动的流函数无法如二二维 一样直接定义，三维问题需要引入多个流函数，而且物理意义也不如二维明确。 2 ，1 3 压力修正法【1 3 ，1 4 ”】 该方法由p a t a n k e r 与s p a l d i n g 于1 9 7 2 年提出，是一种主要用于求解不可压 缩流场的数值方法。它的核心是采用“猜测一修正 的过程。在交错网格的基 础上来计算压力场，从而达到求解n s 方程组的目的。 压力修正法的基本思想可描述如下：对于给定的压力场( 它可以是假定的 值，或是上一次迭代计算所得到的结果) ，求解离散形式的动量方程，得出速度 场。因为压力场是假定的或不精确的，这样，由此得到的速度场一般不满足连 续性方程，囚此，必须对给定的压力场加以修上e 。修正的原则是：与修正后的 压力场相对应的速度场能满足这一迭代层次上的连续方程。据此原则，把动量 方群的离散形式所规定的压力与速度耦合的关系代入连续方程的离散形式，从 6 t - 海大学硕士学位论文 而得到压力修正方程，由压力修正方程得出压力修l l 三值。接着，根据修正后的 压力场，求得新的速度场。然后检查速度场是否收敛。若不收敛，用修正后的 压力值作为给定的压力场，开始下一层次的计算。如此反复，直到获得收敛的 解。 ( 1 ) 交错网格 在上述求解过程中，如何获得压力修正值( 即如何构造压力修正方程) ，以 及如何根据压力修正值确定“正确”的速度( 即如何构造速度修正方程) ，是 s i m p l e 算法的两个关键问题。为此，下面先解决这两个问题，然后给出s i m p l e 算法的求解步骤。 在解决s i m p l e 算法的两个关键问题之前，首先介绍一下交错网格。所谓 交错网格，是指将标量( 如压力p 、温度t 和密度p 等) 在正常的网格节点上储 存和计算，而将速度的各分量分别在错位后的网格上存储和计算，错位后的网 格的中心位于原控制体积的界面上。这样，对于二维问题，就有三套不同的网 格系统，分别用于存储p 、u 和v 。对于三维问题，就有四套网格系统，分别用 于存储p 、u 、v 和w 。 图2 1 给出了二维流动计算的交错网格系统实例。其中，主控制体积为求 解压力p 的控制体积，称为标量控制体积。控制体积的节点p 称为主节点或标 鼍节点，如图2 1 ( a ) 所示。速度u 在主控制体积的东、西界面e 和w 上定义 和存储，速度v 在主控制体积的南、北界面s 和1 1 上定义和存储。u 和v 各自 的控制体积则分别以速度所在位置( 界面e 和界面n ) 为中心的，分别称为u 控制体积和v 控制体积，如2 1 ( b ) 和2 1 ( c ) 表示。可以看到，u 控制体积 和v 控制体积是与主控制体积不一致的，u 控制体积与主控制体积在x 方向有 半个网格步长的错位，而v 控制体积与主控制体积在y 方向上有半个错位。 n s - 1 1 l v n w 一彩- u e i v 。 w e n s 州 v 。 钐么 w 一 一 p 一 ， ：朔y ， i v 。 - s e n s 矿， f a v 。 w 一钐么 - - u wu e l l v 。 e f i f , 人学硕j ：学位沦文 ( a ) 主控制体积 ( b ) u 控制体秋 ( c ) v 控制体积 图2 1 交错网格示意图 现以笛# 尔坐标系下二维稳态问题为例，交错网格及其编码系统如图2 2 所示。图中实线表示计算网格线，实心小圆点表示计算节点( 即主控制体积中 心) ，虚线表示主控制体积界面。 v 控制体积 - - 划- - - ll nl ! 7 lj i 丝e l i - w山 lisl ii ll7xi - - - - s - 、v 控目 一 ii i ilr 一 - 控制体积 体积 i - i 上- i1土 i + il + l 图2 2 交错网格及其编码系统 x 方向的动鼍方程使用u 控制体积，在位置( i ，j ) 处的关于速度甜川的动量方 程的离散形式为： j u i ab u ，jn + b ( 2 1 0 ) ，j2 乙月柚+ ，- 1 j 一，j ，+ ( 2 ) 式中：a ，是u 控制体积的东界面或两界面的面积，即a ，= y 川一y _ ，： 川，一p 川n 川为u 单元控制体所受的x 方向的总压力差；b 为不包括压力在 内的源项部分。 y 方向的动量方程使用v 控制体积，其离散形式为： 1 弋，、 a t j v t j2 乞a h h v 柚+ ( p ，j - 1 p ，】归，+ b ， ( 2 11 ) ( 2 ) 速度修上e 方程 设有初始的猜测压力场p + ，动最方程的离散方程可借助该压力场得以求解， 从而求出相应的速度分量u 和v 。 口川u i + , j = a 柚z ，二+ 山p l + 扣+ 6 川 ( 2 1 2 ) 1 1 1 i l 弘 弘 j ， r 卜 ：海人学硕上学位论丈 口，=。v+njv ab v n b + 川- p ；如+ 6 ， 口， ，。乙n+ p ，一1 p ，+ 6 ， 定义压力修正值p 为正确的压力场p 与猜测的压力场p + 之差，有： p ：p + + p 同样的，我们定义速度修正值u 和v ，则有： u = u + + u v = v + + v 将方程( 2 1 0 ) 减去( 2 1 2 ) ，方程( 2 1i ) 减去( 2 1 3 ) ，得到： 吼，甜：，= a ：。+ o p j h ，口，j 甜，- ，2 己柚”肭+ p ，- | - ，p ，m ， a l 。= q n 奠南七b l 王 一p l 淞。0 1 - 、| ，j 。乙o n 心曲七p l 。p l j p l ，j 为简化求解过程，略去z ，二和z a n b v n b ，于是有： “：，= z ，0 ：吐j p ) ，) v ，。 其中， ，) 川 d j ：红 a i j 速度修正后： p j ) 一a l ? j 口，j ，= “i ，+ 珥，0 ) - 1 ，p j ) ，= v ：。，+ 4 ，) 厶，) 对于“川，与v ， j + 1 ，存在类似的表达式： 弘i + 1 i2 “i + 1 j+ 谚“ j b 。p “1 、 v ，j 卅= v ：小，+ 4 _ 。0 j 厂p j 冉。) 其中， d “j ，2 a ，+ l ，j 口f + 1 j ，d p l = ( 3 ) 压力修正方程 a l ? j a ，j + 9 ( 2 1 3 ) ( 2 1 4 ) ( 2 1 5 ) ( 2 1 6 ) ( 2 1 7 ) ( 2 1 8 ) ( 2 1 9 ) ( 2 2 0 ) ( 2 2 1 ) ( 2 2 2 ) ( 2 2 3 ) ( 2 2 4 ) ( 2 2 5 ) ( 2 2 6 ) i ：海大学硕上学位论文 二维稳态问题连续性方程为： 皇+ 旦：o ( 2 2 7 ) 苏 砂 钊。对图2 2 所示的标量控制体积，连续性方程满足如下离散形式： i ( 刎) 川厂( 刚l ，小l ( 刚) 厂( p v a ) j = 0 ( 2 2 8 ) 将正确的速度值，即式( 2 2 2 ) ( 2 2 5 ) 带入连续性方程( 2 2 8 ) ，整理后 得到p r 的离散方程： a d p ；= 口，+ l ，p j + l ，+ 口l - i , d p ) - l _ ，+ 口i , j + l p i ，+ l + 口1 d - 1 p ；，j i + 6 ；， ( 2 2 9 ) 其中， 口，+ 。、，= ( p d a ) ，+ 1 ， ( 2 3 0 ) a ，- l = ( 舢) u ( 2 3 1 ) 口“+ l = ) + i ( 2 。3 2 ) 口一= ( p d a ) ， ( 2 3 3 ) a ，2 日，+ l ， ，+ 口，一1 ，+ a ，+ l + a ，j i ( 2 3 4 ) 6 j ，= b + 么) ，厂妇彳h ，+ b 彳) ，厂b 么k 。 ( 2 3 5 ) 其后有许多学者改进了s i m p l e 算法，改进的目的都是为了提高计算收敛 速率，比如s i m p l e c 、s i m p l e r 、p o s i 算法等。 2 2 预处理方法 2 2 1 预处理方法介绍 预处理方法是由b r i l e y ，m c d o n a l d 和s h m a r o t h 最先提出来的。预处理方 法是育接对可压缩n a v i e r s t o k e s 方程组施加一个预处理矩阵，通过改变方程组 对流项系数矩阵的特征值，使特征值的量级相同或接近，从而消除方程组的刚 l o 洳人乎顾 ：。z 化论文 性，起到加速收敛的目的。 对于一般曲线坐标系下可压缩n s 方程组 8 w 8 e l8 f 8 g t8 ra s t 8 t 一+ 一+ 一+ 一= 一+ 一+ 酞8 芎 a q8 8 芎a r l8 乞 在时间导数项乘一预处理矩阵r 后的n s 方程组为 a q a w ta e t a f ta g ta r a s ja t t 8 za | a 芎8 珥8 8 芎a q8 本文采用双时间步迭代，f 为虚拟时间，t 为物理时间。 如果采用压力、速度和密度为原始变量， 令o = p 列 v w p ，9 = f ： 专e ( 2 3 6 ) ( 2 3 7 ) ( 2 3 8 ) 方程进一步化为 _ a q + r - ，婴型+ r 一一墨望+ r 一，堡型+ r ，篓型： a f a q a t坦鸳坦却弛a f ( 2 3 9 ) r t ( 堡型+ 堡型+ 竺望1 、a q | 8 芎a q ? a qa q l8 。 此时对流项系数矩阵为 彳：r - l 竺，艿：r 一- 堡，c ：一塑 a q a q 7 a q 可以验证，只要预处理矩阵r 构造合理，就可以保证对流项系数矩阵a ，b ，c 的特征值保持在同量级，从而达到消除方程刚性，起到加速收敛的效果。然 后可以对预处理后的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 a v i e r s t o k e s 方程组更符合流体流动的物理特性。此 外，低速到超音速很大范围内的流动，预处理方法只采用l 司一个控制方程描述， 求解也采用同个差分格式和时间推进格式，对于这种既有高速又有低速的流 动问题( 如收缩性很强的喷管内流场) 提供了很大的方便，因此预处理方法也 l ：海人学硕i ：学位论文 得到了许多研究者的高度重视和推广。预处理方法的关键是预处理矩阵的构造， 不同的预处理矩阵，其功能也有一定的区别，比如有专门针对低速流动的预处 理矩阵，有针对全马赫数流动的预处理矩阵等等，下面介绍本文预处理矩阵的 构造。 2 。2 2 本文预处理矩阵的构造 预处理矩阵的合理构造是预处理方法成功的前提。本文采用压力、速度和密 度为原始变量， 令曰= ( 夕，甜，1 ，比p ) 丁 守恒性变量为q = ( p ，脚，p ) 7 其中e = p h - p = p 卜1 + v 2 + w 2 ) 3 一p 假设速度u ，v ，w ，压力和密度是线性无关的五个变量，可得转换矩阵为 m ：墼： o q 令：尺= 10 一i t 1 1 ， 0 一w 0 一( 向一型2 一“lj 得到k m = l 0000 0 p 0 00 00 p 00 000 d 0 政口o o o 力尸一1 1 2 o0 oo 0o 1o 一1 ，一w1 一 厅 0 0 0 o 尸 d o o o p o o p o 伽 。 夕o o 。 玎 v w 叩 上二海人学硕 ：学何论文 令。：去一三， u ，二a = 当t 1 2 + v 2 + w 2 a v l - i u r f 2 当u 2 + v 2 + w 22a v n u m r 2 和a v r l 为控制参数，上述开关函数能够使得时间推进法在流场低速区能够 收敛。 将k m 矩阵的右上角的元素。用。= 万1 一事取代，得到： l 。= 1000 0 p 000 00 p 00 00 0p0 础弋p 00 0 础j p 一1 由矩阵k ，可以求得矩阵k 的逆矩阵为： k ： 1 0 u 1 v 0 w 0 卜叠2 笪) “ lj 将k 与l 。相乘得到本文预处理矩阵 r = 髟= 0 0 0 1 w o o o o 1 l 何+ 础r 口p up 4 , 鹏+ 力尸一1j 可以验证，经过本文预处理后的可压缩n a v i e r - s t o k e s 方程组，其对流项系数矩 阵的特征值保持在同一量级。 景 囝舾国 娜饱帕 0 o 0 p o o p o o p o 0 “ y w 海人学坝l ：学位论文 第三章状态方程 3 1 不同工况状态方程 状态方程被用来封闭流体力学方程组。对于理想气体采用p u = r t 、实际 气体可以采用范德瓦尔气体状态方程等。对于液体工质，如水，由于没有简单 准确的状态方程，因此，可以采用“水和蒸汽热力学性质表”，即i a p w s - i f 9 7 ， 来表达状态参数之间的关系。 3 2i a p w s i f 9 7 介绍f 1 8 】 在2 0 世纪6 0 年代，水和蒸汽热力学性质的工业公式称为“适用于工业使 用的19 6 7 i f c 公式( z f c 6 7 ) 9 9 0 自19 6 7 年以来，i f c 6 7 正式用于计算任何应用 场合的水和蒸汽的热力学性质，例如动力循环的性能计算。除此之外，i f c 一6 7 还广泛地用于其他工业应用中。然而，在最近几年中，i f c 6 7 显露出许多不足。 这些不足及用于发展准确状态方程的数学方法所取得的进展导致了由水和蒸汽 性质国际协会( i a p w s ) 发起和协调了该国际研究课题，以发展新的公式。 开发新工业公式这一课题是由i a p w s “工业计算”工作组发起和连续指导， 此工作组的许多成员是国际动力循环公司的代表。德国的b 卢克斯为现任主席， 该工作组其前身是“i a p w s 工业计算分委员会”。在1 9 9 0 年布宜若斯艾利召开 的i a p w s 年会上，成立了新工业公式工作组，从那时起开始了该综合课题的研 究。由w 瓦格纳任主席，来自七个国家的1 2 名成员组成了这个工作组，直接 负责i a p w s 。i f 9 7 的研究。日本的k m i y a g a w a 主持了“新工业公式评价”组， 根据所需性能和对实际动力循环计算的影响，检验了该公式，完成了整个课题 中非常重要的工作。 1 9 9 7 年，在德国埃朗根召开的水和蒸汽性质国际协会年会上，采用了新公 式作为水和蒸汽热力性质的国际工业标准。新公式i a p w s i f 9 7 取代了前i f c 一6 7 工业标准公式。与i f c 一6 7 相比，i a p w s i f 9 7 有效提高了热力学性质的计算精 1 4 l ：海人。学硕一 二学位论文 度和速度【1 8 1 。 3 3 本文状态方程 对于新的“水和蒸汽热力学性质的工业公式”，i a p w s i f 9 7 提供的公式不 能直接应用于c f d 的计算中。本文利用根据i a p w s ，i f 9 7 编写的动态链接库函 数的调用接i ：3 ，把提供的库文件( w a s p c n l i b ) 直接导入f o r t r a n 库文件中。把 提供的模块文件( w a s p c n f g o ) d h 虱j 本文的生成表格程序中。在这些库文件中， 由压力和密度可以求得其余状态参数。本文采用以( p ，p ) 为变量插值求得其余 状态参数来封闭n a v i e r - s t o k e s 方程组。 上海大学坝士学他论义 4 。l 控制方程 第四章数值方法 控制方程组由连续性方程、动量方程和能量方程组成，在三维贴体坐标系 下的具体形式为： 警+ 箦+ 筹+ 嚣= 姜- t i - 蓦- i i - 善 c 4 t ， 叶一一叶一一十一： t t i ， 执8 薹a 玎a a 鼍a 玎a ：一1 j p p u p v p w e 1 乞= j n 1 月= 一 j p u p u u + 善，p p v u + 乞，p p w u + 参p ( p e + p ) u 0 彘f x k 彘z - y k 甑f z l 专k 8 t 。 1 6 = 一 j 1 ，= 一 j 0 7 7 iz x k 刁t f y k 7 7 t f z l qk 8k p v p u v + 7 7 ，p p v v + 7 7 ，p p w v + 7 ：p ( p e + p ) v 丁：! j o t f 溅 氕f y k k 呶 。p 。 一 1 u = j p w p u w j + 乞x p p v w + f 、，p p w w ”+ ：p ( p e + p ) w 气。= 去c ( 考+ 警 + 见罄磊，一； l 反= 一g ，+ u z y o c 十v f m ，+ w f 岛= 一q y + “r 掣+ v r ，+ w f 胆 4 乙= 一g ：+ 甜f 。+ ，丁。+ w f 。 v = u i + 巧+ w k ，u = v 孝y ， v = vr l v ， w = v f 矿 其中，= ( p , o up y 夕wp ) r 为守恒型变量，夕为密度，u ，v ，w 为直 角坐标系下的速度分量，e = p 一尸= p e 1 i ( , 1 2 1 4 2 2 ) 一p ，为总内能，h ：沁人。学颂 ：学f _ ：f ：沦丈 为焓。雅可比矩阵j ：娶掣，e 、f 与g 分别为孝、7 7 和f 方向上的无粘通 o ( x ，y ，z ，f j 量，尺、s 与丁分别为善、刁和f 方向上的粘性通量。 另外，必须补充合适的状态方程使得上述方程组封闭。根据水和蒸汽的热 力学性质图表，利用插值的方法可建立温度、粘度和内能等参数与压力和密度 之间的关系。 4 2 计算方法与差分格式 4 2 1 时间推进法 1 9 2 0 2 1 , 2 2 1 如今，计算流体力学中广泛应用时间推进法数值求解e u l e r n a v i e r - s t o k e s 方程组。时间推进法是一种渐进方法，其基本思想是从非定常e u l e r 方程组或 非定常n a v i e r - s t o k e s 方程组出发，沿时间方向推进求解，由此得到的当时间趋 近于无穷大时的渐进解，即为所要求的定常解。它既能求得流动定常解又能模 拟流体运动的非定常过程。因而是应用范围极广的一般性方法。 时问推进方法可分显式和隐式两类。显式方法的优点每推进步，计算量 和存储量都较少，程序简单；缺点是时间推进步长受稳定性条件限制，计算c f l 数过小，效率比较低。目前使用较多的显式方法为多步r u n g e k u t t a 法。隐式 方法的优点是在对模型方程进行线性稳定性分析时一般是无条件稳定的，在数 值求解多维e u l e r n a v i e r - s t o k e s 方程组时，计算步长也可较大而使整体效率较 高：缺点是在每时间推进步都需要解线性方程组，因而计算量和存储量都较 大。著名的隐式方法有b e a m w a r m i n g 隐式格式、m a c c o r m a c k 隐式格式、近似 因子分解法a f 、谱半径分裂的l u s g s 方法等等。 这螳隐式时间推进法在无须考虑时间精度的定常流场计算中得到j 泛应用 并取得很大成功，但是各种高效率的近似处理都将使迭代过程丧欠时间精度。 进入2 0 世纪9 0 年代以来，国外相继发展了双时间步隐式牛顿迭代方法，用来 求解非定常的e u l e r n a v i e r - s t o k e s 方程组，促使了非定常数值模拟的巨大发展。 为了提高非定常流动的时间计算精度，同时又要具有较高的计算效率，j a m e s o n 提出了一种蚁时间步( d u a l t i m e s t e p ) 方法，在冻结的真实时刻点上引入类似牛顿 迭代的虚拟时间迭代过程，通过这种内迭代过程来提高在a f a d l 和l u 处理 过程中所损失的时间精度。 在贴体坐标系下三维非定常n a v i e r - s t o k e s 方程组进行空间离散后，可以写 成如下半离散形式f 1 9 _ = 8 u _ + r ( 【，) ：0 d f 式中万：j u 。 将( 4 2 ) 式改写为 j - ) a _ q _ u + 尺( u ) ：o 西 、7 ( 4 2 ) ( 4 3 ) 在n + l 时刻，对式( 4 - 3 ) 采用时间二阶精度的隐式三点后差离散，得到时间二阶 精度的离散方程为 j 一- 三竺：! 二! 竺：竺：1 2 f + r ( u 1 ) = 0 若对粘性通鼍采用显式处理，则式( 4 4 ) 可表示为 j t 兰竺：! 二兰竺：竺：1 2 ， ( 4 4 ) + ( a 0 面”+ 1 + 6 j ，i ”+ 1 + 0 i ；”+ 1 ) 2 ( 4 5 ) 一婶；ey + 6 ；fy + 6 # g v ) n 式中，艿是差分算子，丁表示中心差分算子。 对于式( 4 5 ) 的非定常隐式差分方程，j a m e s o n 引入内迭代求解方法，即在式( 4 5 ) 左端增加虚拟时间导数项，并采用一阶前差处理，得到 j t 竺：! 二竺：+ j 一，三竺：! 二兰望：竺：! +， 一十j一十 a f 2 a t ( 4 6 ) ( 睡一e n l + 万，芦尸+ + 8 c - 否) = 一( 夏i r + 虿7 矿+ 虿否r ) 尸 式中f ，a t 分别为虚拟时间步长和真实时间步长；p , n 分别为虚拟时间迭代步 数和真实时问推进步数。 j a m e s o n 引入方程( 4 6 ) 的目的是：方程( 4 6 ) 可以理解为是方程( 4 5 ) 的定常 上二海人学硕 j 学位论文 解，因为方程( 4 6 ) 的定常解意味着j - i 旦二二! 三里：0 ，此时方程( 4 6 ) 幂i ：l * 一f - ( 4 5 ) 2 _ 等价，由于方程( 4 5 ) 是时间二阶精度的非定常流动差分方程，此时方程( 4 6 ) 的 定常解就是所要求的二阶精度的非定常解。对于方程( 4 6 ) 的定常解，可以使用 很成熟的多种时间推进方法及加速收敛方法得到定常解。 4 2 2 l u 分解与追赶法【2 3 】 对带型矩阵进行三角分解时，三角矩阵也是带型的，因而可以使算法大大 简化。对于三对角矩阵 a = b lc l 口，b ， 口3 c 盯一j b 。 存在三角分解a = l u ，l 是下三角矩阵，u 是上三危矩阵。 l =，u = ( 4 7 ) ( 4 8 ) l 中与u 中待定元素总数恰好为3 n 一2 。由a = l u ，得出如下确定l 与u 的元素 公式： v f = c ，i _ 1 ，2 ，n - 1 ；( 4 9 ) z f - ：6 ，= - 1 ，甜，= b i - 饥l ，i _ 2 ，3 ，一n o ( 4 1o ) 实现a 的l u 分解后，即可应用与求解以a 为系数矩阵的方程组a x = d ， 等价于依次求解l y = d ，u x = y ，这里d = ( 4 ，d ：，j ，) 7 。是已知向量，这种求解方 法也叫追赶法。 1 9 一 玎 甜 f ：j f ：f 人学珂! 上学化论文 4 2 3r o e 格式【1 9 2 4 , 2 5 , 2 6 ,