基于预处理技术的微可压缩液体流场模拟方法（可编辑）

基于预处理技术的微可压缩液体流场模拟方法（可编辑） 基于预处理技术的微可压缩液体流场模拟方法 上海大学 硕士学位论文 基于预处理技术的微可压缩液体流场模拟方法 姓名:孙建 申请学位级别:硕士 专业:流体力学 指导教师:黄典贵 20081001 I:海人学硕上学位论文 摘 要 一般情况下，液体流动可近似作不可压缩流场处理。压力修正法是求解不 可压缩流场的最常用的方法;也有学者提出采用预处理技术，利用时间推进法 进行不可压流场的求解。但是真正的不可压流体是不存在的，如果将液体工质 看作可压缩流体，采用常规的时间推进法求解这类问题，一般很难得到收敛解。 其原因在于水的音速相当高，,般流体机械内的流速相对其音速很小，’方程组 的对流项系数矩阵具有很强的刚性。通过引入预处理矩阵，改变方程组的特征 值，使对流项系数矩阵的特征值处于同一量级，从而可以达到消除刚性与加速 计算收敛性的目的。 本文以水为工质，考虑低速流动时水的微压缩性，通过以密度和压力为原 始变量的预处理技术，采用时间推进方法，进行流场模拟。利用水和蒸汽的热 力学性质状态参数关系，由压力和密度求解流场模拟过程中所需的其他状态参 数，以封闭可压缩的流体力学方程组。对水中的NACA0015翼型绕流、低雷诺 数的圆柱绕流和NACA0012翼型振荡绕流问题的模拟表明:本计算方法能够较 好地模拟水中的低速流动问题。 关键词:微可压缩 预处理技术 流场模拟 V (1-i11人学颀上学化论文 ABSTRACT flowisoften as Liquid incompressibleapproximately(Thepressure regard correctionmethodis of usedtosimulatethiskindflowfield(However， widely the canalsobeusedtosimulatethe preconditioningtechnology incompressible fiow(Infact，the fluiddoesnotexist(A1 absolutelyincompressible though the l canbe as flowingiquid regarded computationusually compressible，the doesnot whenthe methodisused(Becausethe converge time―marching velocity ofsoundinwaterismuch thanthe offluidinthe higher velocity machinery， thecoefficient matrixofconvectiontermis a veryStiff(Byinducing precondi tioning ofthe ofthecoefficient matrix，themagnitudeeigenvalues matrix ofconvectiontermcan bemodifiedtobeinthesame the order，SO stiffness ofthe is equationselimihatedandthe methodworks( time-marching Inthis isusedas method dissertation，water workingfluid， apreconditioning using and asthe variablesis to density，velocitypressureprimitive adopted simulatetheflow(In to orderenclosethe equationsystem， therelationships between density， pressureandotherstatefunctionsareinducted the byusing thermodynamicsofwaterandsteam(Threetest flowaround property cases， the NACA0015airfoi1 in low numberflowaround and water，theReynolds cylinder theflow around NACA0012 idatethe ofthis oscillatingairfoil，valcapability methodtosimulatethe flow( liquid Keywords:Slight technology, Simulationof flowfield Vl I:海人学硕L"P位论文 原创性声明 本人声明:所呈交的论文是本人在导师指导下进行的研究工作。 除了文中特别加以标注和致谢的地方外，论文中不包含其他人 已发 表或撰写过的研究成果。参与同一工作的其他同志对本研究所做的 任何贡献均己在论文中作了明确的说明并表示了谢意。 签名:丝逛日期:丝壁:!兰:? 本论文使用授权说明 本人完全了解上海大学有关保留、使用学位论文的规定，即: 学校有权保留论文及送交论文复印件，允许论文被查阅和借阅;学 校可以公布论文的全部或部分内容。 保密的论文在解密后应遵守此规定 2(乙 签名:趋盔‘翮签名:瑾、瑶日期: ?opg(f :海大学硕士学位论文 第一章 绪论 1(1课题来源 本课题来源于国家自然科学基金项目，项目编号:50576049。 1(2课题研究的目的和意义 本文通过采用预处理技术，求解可压缩的流体力学方程组，来计算以水为 工质的低速流动问题。利用水和蒸汽热力学性质图表，用插值的方法由压力和 密度求得温度等其他状态参数，用以封闭方程组。该方法旨在统 一可压与不可 压流场的计算方法，最终能够模拟液体工质发生空化现象的复杂流动，并将其 应用于流体机械的流场模拟中，课题组正在进行相关研究。 因为流体都是可压缩的，所以液体流动及气体的低马赫数流动都具有一定 的压缩性，用预处理方法对可压缩N(S方程组进行求解比求解不可压缩N―S方程 组更符合流体流动的物理特性。此外，从低速到超音速很大范围内的流动，预处 理方法只需采用同一组控制方程描述，选择相同的计算格式，从而实现了可压与 不可压流场计算方法上的统一，该方法最终可以应用于可压与不可压、高速与低 速、气体与液体同时存在的复杂流场。 1(3国内外研究概况 1。3(1国外研究概况 体力学方程组旌加一个预处理矩阵，总体上提高了计算过程的收敛速率。Turkel 基于原始变莓和基本的流体力学方程组，针对应用于可压缩方程的人工可压缩 方法，了一个普遍性的预处理矩阵12|。Choi与Merkle将预处理技术用于粘 性流场的模拟之中13J，实践表明:经过预处理后的N(S方程组，计算收敛速率 提高了两个量级。Weiss与Smith报道了一个适用于变密度低速流场求解的预处 L海大学硕:f:学位论文 理方法【4】o 矩阵，该预处理矩阵形式简单，且易于实现低速与高速流动的自动切换pj。 对于微可压缩液体的流场模拟，国外学者一般将微可压缩液体工况流场简化 为不可压流动，采用不可压缩N(S方程组进行求解。典型的解法有人工压缩性方 法、涡量(流函数法和压力修上E法。而采用可压缩N(S方程组求解微可压缩流动的 方法暂时没有报道。 1(3(2国内研究概况 国内学者对预处理矩阵的构造比较少，黄典贵教授在2005年构造了一个基于 压力与焓为原始变量的预处理矩阵，该预处理矩阵适用于求解全马赫数流场16j。 大多数国内学者的研究，都是针对国外学者构造好的预处理矩阵进行直接应用。 比如:两北工业大学韩忠华、乔志德等比较了Choi预处理矩阵与Turkel预处理矩 阵的功能17J，得出如下结论: 1 预处理方法能十分有效地改善低马赫数粘性流动 计算的收敛性和收敛精度; 2 预处理方法同样适用于高速流动，且有一定的收 敛效果; 3 Choi预处理和Turkel预处理对于低马赫数流动计算的性能基本相当; 4 在较高马赫数计算中Choi预处理的稳定性和加速性能略优于Turkel预处理。北 京航空航天大学严明副教授将文献【5】的预处理矩阵应用于叶 轮机械的三维数值 模拟，对于低压压气机算例的计算结果表明，采用预处理方法后可以得到与实验 结果吻合较好的数值模拟结果【8】。中国空气动力发展研究中一11,肖中云、牟斌等【9】， 将文献【4】的预处理矩阵用于计算无粘bump流和绕二维圆柱的非定常流动，得出结 论:数值结果表明预处理技术显著提高了低速流场求解的效率及求解精度，对所 模拟低速的非定常流动而言，所需代价只是传统时间推进法的几分之一。中国空 气动力发展研究中心肖中云、江雄等【?1，将weiss―Smith预处理矩阵应用于悬停 旋翼的模拟中，得到如下结论: 1 在加入低速预处理方法以后，流场的收敛速 率得到了明显提高，在同等迭代步数下得到的下洗速度场更为准确。 2 通过尾 涡耗散的大小验证了低速预处理避开了因特征值差异带来的数值耗散过大的问 题，使可压缩方法具有在计算低速问题时保持精度的特点。 对于微可压缩液体的流场模拟，国内学者一般将微可压缩液体工况流场近似 为不可压流动，采用不可压N(S方程组进行求解。而采用可压缩N(S方程求解的 2 上海人学硕上(学位论文 方法暂时也没有报道。 1(4论文的主要研究内容 本论文是以作者攻读硕士学位期间承担课题的工作为基础，在第一章中阐 述了课题研究的来源、目的、意义以及国内外研究的现状。第二章阐述了预处 理方法及本文预处理矩阵的构造，顺带介绍下求解不可压缩流动的人工压缩性 方法、涡量(流函数方法及压力修正法。第三章阐述了状态方程，如何利用水和 蒸汽性质表获得水的其余状态参数。第四章讲述了本文控制方程的数值方法， 采用时间推进法进行求解，采用Roe格式离散对流项，因为控制方程的特征值 比较小而引入熵修正，通过限制器使离散格式精度提高。第五章讲述了三个算 例的数值模拟结果与。利用本文方法对NACA0015翼型绕流、低速圆柱绕 流及NACA0012翼型振荡绕流进行数值模拟，然后对得到的数值模拟结果与文 献中的实验结果或数值模拟结果进行对比分析。最后第六章全文。 J:海人学顾十学位论文 第二章预处理方法 2(1不可压缩流动解法 当流体运动速度比声波传播速度小很多时，流体的密度变化不大，可近似 为常数，此时描述流体运动规律的N(S方程组可以大为简化，称其为不可压N―S 方程组，方程如下: 锄 OvOw一 ――+――+―― 0 ax 却 aZ ―― 塑4-甜塑+v塑+W--锄:一!望+闪:搿甜――+V――+ 一一土+DV。搿 Ot Ox Oy 0z p0x 2(1 。1， ――+z，――+V――尘+z，生+v堡4-w塑:一三望+O―V:1， W―― 一一二+ 魂 ax 劫 娩 p却 2w。W 一0w+甜坐+v业+w丝:一上望+O―V一十甜一+V一十W―― 一一二+ at 融 两0zp0z 其中u ,。 从形式上看不可压N(S方程组比可压N(S方程组较为简单，但足在求解上 引入了新的困难，由于祟+芸+譬:o使方程组不再是进化型，难以采用己发 t强 Uy Oz 展比较成熟的时间推进法求解。原来的进化型连续方程变成了在 每个时间步上 对速度矢晕的约束条件。此外对压力缺少边界条件也给不可压N―S方程组的求 解带来一定的困难。针对这,问题，学者们发展了许多处理方法，比如人工压 缩性方法…??121、涡鼍(流函数法、压力修正法[13,14]等等。 在低速情况下，如果把流体仍然看作是可压缩的，不做任何简化直接求解 可压缩Navier-Stokes方程组，会出现的主要困难是数值计算过程的收敛速率较 慢，甚至可能不收敛。原因在于，低马赫数条件下流动速度与声速相差很大， 导致可压缩Navier―Stokes方程组的对流项具有很强的刚性。针对这一问题，学 者们发展了预处理方法fI(5】。下面分别介绍上述提到的各种方法。 4 :海人学硕1:学位论文 2(1(1 人工压缩性方法【i02(1s] 人工压缩性方法是把不可压缩流体看作一种虚拟的可压缩流体。在求解中 仍采用非定常的动量方程和经过人工非定常化的连续方程。这样处理的结果使 得方程组变为进化型，可以通过通常的求解进化型方程的方法进行求解。这一 方法最早由Chorin和Yanenko提出。 经过变型的连续方程为 塑+c2弋7(莎:o 2(2 at 这里c2是一任意的自由参数。在未达到定常状态之前该式没有任何意义。在达 到定常状态时雩磊 0，即V(莎 0。其中参数c2的选取应使求解过程收敛和 加速收敛。 2(1(2涡量一流函数法【’61 通过引入涡晕方程和流函数，对二维不可压N(S方程组进行变换，写成涡量 流一函数的形式。对如下二维不可压N(S方程组: 丝+翌:0 2(3 叙 砂 一“ 塑+材塑+v塑:一』至+上一V:“ 2。4 一十材――+V一 一一二+―― at ax av P8x Re 。l，zp 堡+甜鱼+v鱼:一三至+上一V 2(5 ――十甜一十V―― 一一二+一 a， 出 砂 p砂 Re 引入流函数妒和涡量力: 8p 8p zf 二，V 一二 2(6 砂 缸 J:海人学硕‘l:学位论文 国:鱼一一Ou国 ――一一 t‘??，， 2(7 苏 方 2(4 式对y求导， 2(5 式对X求导，然后消去压力项得到 方程 ‘ ――+z，――+V―― ――丝+z，塑+v塑:上一V2C_O k二??oJ 2(8 ot 融 却 Re 将 2(6 式代入式 2(7 得 v:妒:磐+宴:一国 2(9 钡。 咖’ 方程 2(8 与 2(9 构成了涡量一流函数型的N(S方程组。 可用不同的方法对 2(8 与 2(9 式进行离散和求解。可以看出，方程组中不 再出现压力项，使得方程更为简单。当妒、国以及u和V求出后，可利用动量 方程或对应的压力Possion方程求出压力。 虽然涡量一流函数方法在计算二维不可压缩N―S方程组中取得了成功，但是 把该方法推广到三维时却出现了困难。其原因是三维流动的流函数无法如二二维 一样直接定义，三维问题需要引入多个流函数，而且物理意义也不如二维明确。 2，1(3压力修正法【13，14??”】 缩流场的数值方法。它的核心是采用“猜测一修正"的过程。在交错网格的基 础上来计算压力场，从而达到求解N(S方程组的目的。 压力修正法的基本思想可描述如下:对于给定的压力场 它可以是假定的 值，或是上一次迭代计算所得到的结果 ，求解离散形式的动量方程，得出速度 场。因为压力场是假定的或不精确的，这样，由此得到的速度场一般不满足连 续性方程，囚此，必须对给定的压力场加以修上E。修正的原则是:与修正后的 压力场相对应的速度场能满足这一迭代层次上的连续方程。据此原则，把动量 方群的离散形式所规定的压力与速度耦合的关系代入连续方程的离散形式，从 6 t-(海大学硕士学位论文 而得到压力修正方程，由压力修正方程得出压力修ll三值。接着，根据修正后的 压力场，求得新的速度场。然后检查速度场是否收敛。若不收敛，用修正后的 压力值作为给定的压力场，开始下一层次的计算。如此反复，直到获得收敛的 解。 1 交错网格 在上述求解过程中，如何获得压力修正值 即如何构造压力修正方程 ，以 及如何根据压力修正值确定“正确”的速度 即如何构造速度修正方程 ，是 SIMPLE算法的两个关键问题。为此，下面先解决这两个问题，然后给出SIMPLE 算法的求解步骤。 在解决SIMPLE算法的两个关键问题之前，首先介绍一下交错网格。所谓 交错网格，是指将标量 如压力P、温度T和密度P等 在正常的网格节点上储 存和计算，而将速度的各分量分别在错位后的网格上存储和计算，错位后的网 格的中心位于原控制体积的界面上。这样，对于二维问题，就有三套不同的网 格系统，分别用于存储P、u和v。对于三维问题，就有四套网格系统，分别用 于存储P、u、v和W。 图2(1给出了二维流动计算的交错网格系统实例。其中，主控制体积为求 解压力P的控制体积，称为标量控制体积。控制体积的节点P称为主节点或标 鼍节点，如图2(1 a 所示。速度U在主控制体积的东、西界面e和w上定义 和存储，速度v在主控制体积的南、北界面s和11上定义和存储。u和v各自 的控制体积则分别以速度所在位置 界面e和界面n 为中心的，分别称为U 控制体积和V控制体积，如2(1 b 和2(1 c 表示。可以看到，U控制体积 和V控制体积是与主控制体积不一致的，u控制体积与主控制体积在X方向有 半个网格步长的错位，而v控制体积与主控制体积在y方向上有半个错位。 -11 ?州 矿，，‘，， lVn V。 f,Av。,, n ―-? n ??一 钐么 n 钐么―--? ??W一 E ??W一 ??P,一 E ??W一 E Ue ‰?,， Uw Ue 彩 y,，:朔 l S S S IV。 IV。 lV。 ?? , -S ??, W f(if,,人学硕j:学位沦文 a 主控制体积 b u控制体秋 c v控制体积 图2(1交错网格示意图 现以笛#尔坐标系下二维稳态问题为例，交错网格及其编码系统如图2(2 所示。图中实线表示计算网格线，实心小圆点表示计算节点 即主控制体积中 心 ，虚线表示主控制体积界面。 v控制体积 -?-’ 划 ---?’ 弘1 控制体积 弘 l l n l!,7‖ l J I l I-w 山丝E J l I S l i I ， l l 7 I X ― 一 ----S(-???,、v控 目 体积 1^1i r 卜l I I ? ? I I l r 一 ???? - ― i-I上-I1 土 l+l i+I 图2(2交错网格及其编码系统 x方向的动鼍方程使用u控制体积，在位置 i，J 处的关于速度 甜川的动量方 程的离散形式为: ( jUi a 2( 2(10 ，J bu ，’J 2乙月柚+?，-1(J一，’J‖，((，+?n?+b? 式中:A，，(，是u控制体积的东界面或两界面的面积，即A，，(， Y川一Y_，: ?川(，一P川n川为u单元控制体所受的x方向的总压力差;b为不包括压力在 内的源项部分。 Y方向的动量方程使用v控制体积，其离散形式为: 1(弋 ， 、 at(jvt(J 2(11 2乞aHhV柚+ p，(J-1‘P，(】归，(，+b，(， 2 速度修上E方程 设有初始的猜测压力场P+，动最方程的离散方程可借助该压力场得以求解， 从而求出相应的速度分量u’和v‘。 口川Ui+,J ?a柚z，二+?山(Pl+?扣?+6川 2(12 :海人学硕上学位论丈 a 2(13 口，( 口，((((， ?。V+njV，((，。乙nbVnb+?川-p;?如?+6，(+p， ?一1‘p，?‖，(』+6，( 定义压力修正值p’为正确的压力场P与猜测的压力场p+之差， 有: 2(14 p:p++p’ 同样的，我们定义速度修正值u’和v，，则有: U U++U’ 2(15 V V++V’ 2(16 将方程 2(10 减去 2(12 ，方程 2(1I 减去 2(13 ，得到: 吼，甜:((， ?a‖:。+?o―pj?h((，口，(J甜，(-，2 己柚”肭+p，-|’-，‘p，?，M((， 2(17 01(-、|，j。乙on心曲七,pl。‘‘plJpl，jal’?。 ? Qn奠南七bl王^一pl淞。 2(18 为简化求解过程，略去?‰z，二和ZanbV'nb，于是有: “:，(， z((，0:吐J(p ，(， 2(19 ’ V，，，。 ，? 川(pj? 2(20 其中， dJ:红 一Al?j 2(21 ai(J 口，，J 速度修正后: ， “i，+珥((，0 -1((，(pj? 2(22 ， V:。，+4，，? ^((厶，(， 2(23 对于“川，，与V，'j+1，存在类似的表达式: 弘i+1(I 2(24 2“i+1(j+谚“(Jb。(p“(1、 V，(j卅 V:小，+4_。0j厂pj冉。 2(25 其中， A，+l，J Al?j d“J((，2 ，d^pl 2(26 口f+1(j a，(j+ 3 压力修正方程 9 I:海大学硕上学位论文 二维稳态问题连续性方程为: 皇四+旦四:o 2(27 苏 砂 钊。对图2(2所示的标量控制体积，连续性方程满足如下离散形式: I 刎 川厂 刚l，小l 刚 ?厂 pvA ?j 0 2(28 将正确的速度值，即式 2(22 , 2(25 带入连续性方程 2(28 ，整理后 得到pr的离散方程: 2(29 a,(dp;? 口，+l’(，pj+l，(，+口l-i,dp -l'_，+口 I,J+lP'I((，+l+口1(d-1p;，J―I+6;((， 其中， 口，+。、， pdA ，+1((， 2(30 a，-l? 舢 u 2(31 口“+l 九 ?+I 2。32 口?一 pdA ，(， 2(33 2(34 a，，(，2日，+l，'，+口，一1((，+a，((，+l+a，，J―I 6j((， b+么 ，厂妇‘彳h(，+b’彳 ，厂b’么k。 2(35 其后有许多学者改进了SIMPLE算法，改进的目的都是为了提高计算收敛 速率，比如SIMPLEC、SIMPLER、POSI算法等。 2(2预处理方法 2(2(1预处理方法介绍 法是育接对可压缩Navier(Stokes方程组施加一个预处理矩阵，通过改变方程组 对流项系数矩阵的特征值，使特征值的量级相同或接近，从而消除方程组的刚 lO (洳人’’乎顾 :。Z化论文 性，起到加速收敛的目的。 对于一般曲线坐标系下可压缩N(S方程组 8W’8El8F’8Gt8Rast8T’ 2(36 一+一+一+一 一+一+―― 酞8芎 aq 8‘8芎arl8乞 在时间导数项乘一预处理矩阵r后的N(S方程组为 (aQ’ awtaEtaFtaGtaR’ aSjaTt 2(37 8z a| a芎8珥8?8芎aq8? 本文采用双时间步迭代，f为虚拟时间，t为物理时间。 如果采用压力、速度和密度为原始变量， p 列 f: V 2(38 令O ，9’ W p 专E 方程进一步化为 _aQ'+r-，婴型+r一一墨望+r一，堡型+r，篓型: af aQ’at 坦’鸳 坦’却 弛’af 2(39 rt 堡型+堡型+竺望1 、aQ|8芎aQ?aqaQl8 。 此时对流项系数矩阵为 彳:r-l竺，艿:r一-堡，c:一塑 7 aQ’ aQ aQ’ 可以验证，只要预处理矩阵r构造合理，就可以保证对流项系数矩阵A，B，C 的特征值保持在同,量级，从而达到消除方程刚性，起到加速收敛的效果。然 后可以对预处理后的Navier-Stokes方程组采用很成熟的时间 推进法进行求解。 因为流体都足可压缩的，所以用预处理方法对可压缩Navier(Stokes方程组 进行求解比求解不可压缩Navier(Stokes方程组更符合流体流动的物理特性。此 外，低速到超音速很大范围内的流动，预处理方法只采用l司一个控制方程描述， 求解也采用同,个差分格式和时间推进格式，对于这种既有高速又有低速的流 动问题 如收缩性很强的喷管内流场 提供了很大的方便，因此预处理方法也 l:海人学硕i:学位论文 得到了许多研究者的高度重视和推广。预处理方法的关键是预处理矩阵的构造， 不同的预处理矩阵，其功能也有一定的区别，比如有专门针对低速流动的预处 理矩阵，有针对全马赫数流动的预处理矩阵等等，下面介绍本文预处理矩阵的 构造。 2。2(2本文预处理矩阵的构造 预处理矩阵的合理构造是预处理方法成功的前提。本文采用压力、速度和密 度为原始变量， 令曰 夕，甜，1，，比p 丁 守恒性变量为Q p，脚，‖，‖，P 7’ 其中e pH-P p卜1’+v2+w2 3一P 假设速度u，v，W，压力和密度是线性无关的五个变量，可得转 换矩阵为 。 。 o o 0 玎 o o 0 M:墼: v 夕o p o 0 oq o O w叩 ‖ o伽 p‖ D 一 , 厅(尸 1 0 O 0 一It 1 O O 令:尺 ―1， 0 0 O 一w 0 1 O 一 向一型2]一“L J 一 1， 一w 1 l 0 0 0 0 0 0 0 0 p 得到K(M 0 0 p 0 0 0 0 0 D 0 政口o o o力尸一1 12 上二海人学硕 :学何论文 令。:去一三， 当t12+v2+w2 avl-i??URf2 U，二a 景当U2+v2+W22avn??Umr2 ,和avrl为控制参数，上述开关函数能够使得时间推进法在流场低速区能够 收敛。 将K^M矩阵的右上角的元素。用。 万1一事取代，得到: 1 0 0 0 @ 0 0 0 0 p L。 0 0 p 0 0 0 0 0 p 0 础弋p 0 0 0础Jp一1 由矩阵K，可以求得矩阵K的逆矩阵为: 1 0 0 O U 1 0 O K,:v 0 0 O W 0 1 O ? W 1 L j 卜叠2笪 “ 将K,与L。相乘得到本文预处理矩阵 ? O O 0 “ p O O r 髟,??‘ y O p 0 w 0 O p @娜饱帕@囝舾国 ‖p4,鹏+力尸一1 L何+础R口pu J 可以验证，经过本文预处理后的可压缩Navier-Stokes方程组， 其对流项系数矩 阵的特征值保持在同一量级。 海人学坝l:学位论文 第三章状态方程 3(1不同工况状态方程 状态方程被用来封闭流体力学方程组。对于理想气体采用pu RT、实际 气体可以采用范德瓦尔气体状态方程等。对于液体工质，如水，由于没有简单 准确的状态方程，因此，可以采用“水和蒸汽热力学性质表”，即IAPWS-IF97， 来表达状态参数之间的关系。 3(2 IAPWS(IF97介绍f18】 在20世纪60年代，水和蒸汽热力学性质的工业公式称为“适用于工业使 用的1 967年以来，IFC(67正式用于计算任何应用 967IFC公式 ZFC(67 990自1 场合的水和蒸汽的热力学性质，例如动力循环的性能计算。除此之外，IFC一67 还广泛地用于其他工业应用中。然而，在最近几年中，IFC(67显露出许多不足。 这些不足及用于发展准确状态方程的数学方法所取得的进展导致了由水和蒸汽 性质国际协会 IAPWS 发起和协调了该国际研究课题，以发展新的公式。 开发新工业公式这一课题是由IAPWS“工业计算”工作组发起和连续指导， 此工作组的许多成员是国际动力循环公司的代表。德国的B(卢克斯为现任主席， 该工作组其前身是“IAPWS工业计算分委员会”。在1990年布宜若斯艾利召开 的IAPWS年会上，成立了新工业公式工作组，从那时起开始了该综合课题的研 究。由W(瓦格纳任主席，来自七个国家的12名成员组成了这个工作组，直接 根据所需性能和对实际动力循环计算的影响，检验了该公式，完成了整个课题 中非常重要的工作。 1997年，在德国埃朗根召开的水和蒸汽性质国际协会年会上，采用了新公 式作为水和蒸汽热力性质的国际工业标准。新公式IAPWS(IF97取代了前IFC一67 工业标准公式。与IFC一67相比，IAPWS(IF97有效提高了热力学性质的计算精 14 l:海人。学硕一 二学位论文 度和速度【181。 3(3本文状态方程 对于新的“水和蒸汽热力学性质的工业公式”，IAPWS―IF97提供的公式不 能直接应用于CFD的计算中。本文利用根据IAPWS，IF97编写的动态链接库函 提供的模块文件 WASPCN(fgo DH虱J本文的生成程序中。在这些库文件中， 由压力和密度可以求得其余状态参数。本文采用以 P，P 为变量插值求得其余 状态参数来封闭Navier-Stokes方程组。 上海大学坝士学他论义 第四章数值方法 4。l控制方程 控制方程组由连续性方程、动量方程和能量方程组成，在 三维贴体坐标系 下的具体形式为: t(t(I， ,叶一一叶一一十一:… c4??t， 执8薹 a玎 a a鼍 a玎 a 警+箦+筹+嚣 姜-tI'-蓦-II-善 P pU’ pV’ pW’ pu puU’+善，P puV’+77，P puWj+乞xP 1 1 1 ( ( 一 ‖:一1 pv 乞 ,pvU’+乞，P ， 一pvV’+77，P U ―― pvW’+f、，P J j J j pw pwU’+参P pwV’+7:P pwW”+?:P e pe+p u’ pe+p v’ pe+p W’ 0 0 O 彘fxk 77Iz(xk ‘tf溅 1 1 n 。 丁:! 月 一 6 一 彘z-yk 刁tfyk 氕fyk j j j 甑fzl‘ 77tfzl‘ ‘k 呶 专k8t q ‘。p。 k8k 气。 去c‖ 考+警]+见罄磊，一;‖ l反 一g，+Uz'yoc十Vfm，+wf《 岛 一qy+“r掣+Vr"，+wf胆 4乙 一g:+甜f。+’，丁。+wf。 V’ V V ui+巧+wk，U’ V孝??y，rl??V， W’ Vf??矿 py夕wP r为守恒型变量，夕为密度，u，v，W为直 其中，‖ ? p,ou 1i ,,12 角坐标系下的速度分量，e p?一尸 pE (1422 ]一P，为总内能，h :沁人。学颂 :学f_::f:沦丈 为焓。雅可比矩阵J:娶掣，E、F与G分别为孝、77和f方向上的无粘通 o x，Y，z，fJ 量，尺、S与丁分别为善、刁和f方向上的粘性通量。 另外，必须补充合适的状态方程使得上述方程组封闭。根据水和蒸汽的热 力学性质图表，利用插值的方法可建立温度、粘度和内能等参数与压力和密度 之间的关系。 4(2计算方法与差分格式 4(2(1时间推进法[19(20(21,221 方程组。时间推进法是一种渐进方法，其基本思想是从非定常Euler方程组或 非定常Navier-Stokes方程组出发，沿时间方向推进求解，由此得到的当时间趋 近于无穷大时的渐进解，即为所要求的定常解。它既能求得流动定常解又能模 拟流体运动的非定常过程。因而是应用范围极广的一般性方法。 时问推进方法可分显式和隐式两类。显式方法的优点每推进,步，计算量 和存储量都较少，程序简单;缺点是时间推进步长受稳定性条件限制，计算CFL 数过小，效率比较低。目前使用较多的显式方法为多步Runge(Kutta法。隐式 方法的优点是在对模型方程进行线性稳定性分析时一般是无条件稳定的，在数 高:缺点是在每,时间推进步都需要解线性方程组，因而计算量和存储量都较 因子分解法AF、谱半径分裂的LU(SGS方法等等。 这螳隐式时间推进法在无须考虑时间精度的定常流场计算中得到J’(泛应用 并取得很大成功，但是各种高效率的近似处理都将使迭代过程丧欠时间精度。 进入20世纪90年代以来，国外相继发展了双时间步隐式牛顿迭代方法，用来 为了提高非定常流动的时间计算精度，同时又要具有较高的计算 效率，Jameson 提出了一种蚁时间步 dual(time(step 方法，在冻结的真实时刻点上引入类似牛顿 迭代的虚拟时间迭代过程，通过这种内迭代过程来提高在AF―ADL和LU处理 过程中所损失的时间精度。 在贴体坐标系下三维非定常Navier-Stokes方程组进行空间离散后，可以写 成如下半离散形式f19’ _ 8U_+R 【， :0 4(2 df 式中万:J―U。 将 4(2 式改写为 4(3 、7 j- a_(q_u+尺 u :o 西 在n+l时刻，对式 4-3 采用时间二阶精度的隐式三点后差离散，得到时间二阶 精度的离散方程为 J一-三竺::!二!竺:?竺::1 4(4 +R U?1 0 2?f 若对粘性通鼍采用显式处理，则式 4(4 可表示为 J―t兰竺::!二兰竺:?竺::1 + a0面”+1+6j，i”+1+,0i;”+1 2 2?， 4(5 一婶;Ey+6;Fy+6#Gv n 式中，艿是差分算子，丁表示中心差分算子。 对于式 4(5 的非定常隐式差分方程，Jameson引入内迭代求解方法，即在式 4(5 左端增加虚拟时间导数项，并采用一阶前差处理，得到 J―t竺:!二竺:+J一，三竺::!二兰望:?竺::!+‘， Af 2At 4(6 一十J 睡一Enl+万，芦尸+’+8c-否,’ 一 夏ir+虿7矿+虿否r 尸 一十 式中?f，At分别为虚拟时间步长和真实时间步长;p,n分别为虚拟时间迭代步 数和真实时问推进步数。 上二海人学硕’ j学位论文 ?2_ 等价，由于方程 4(5 是时间二阶精度的非定常流动差分方程，此时方程 4(6 的 定常解就是所要求的二阶精度的非定常解。对于方程 4(6 的定常解，可以使用 很成熟的多种时间推进方法及加速收敛方法得到定常解。 4(2(2LU分解与追赶法【23】 对带型矩阵进行三角分解时，三角矩阵也是带型的，因而可以使算法大大 简化。对于三对角矩阵 bl cl ? 口， b， 口3 A 4(7 C盯一j b。 存在三角分解A LU，L是下三角矩阵，U是上三危矩阵。 L ，U ? 4(8 ? 一 ,甜 玎 L中与U中待定元素总数恰好为3n一2。由A LU，得出如下确 定L与U的元素 公式: Vf c，，i_1，2，(((，n-1; 4(9 o zf-:6‖， ,-1，甜， bi-饥l，i_2，3，((一no 4(1 实现A的LU分解后，即可应用与求解以A为系数矩阵的方程组Ax d， 等价于依次求解Ly d，Ux y，这里d 4，d:，…，J，， 7。是已知向量，这种求解方 法也叫追赶法。 19 f:jf:f人学珂!上学化论文 4(2(3Roe格式【19(24,25,26,271 1959年，苏联科学家S(K(Godunov提出了一种流体力学 数值计算的新方法， Euler方程，并分片合成整个流场的解。在有限差分和有限体积法里，一般假设 在离散的网格里每个网格内的流动变量是常数。在网格交接处，这些流动变量 就是间断的，而这刚好是Riemann问题，也就是说，网格交界处都是Riemann 问题。 因此，可以这样来进行计算，在t----n时间步，每个网格点的流动变量已知， 为分段常数。然后在每个网格交界处求解Riemann问题;接下来将得到的 Riemann问题的流动变量解在每一个网格内求平均，得到t n+1时间步的流动 变量:再对每个网格交界处求解Riemann问题，如此循环，最后收敛就得到了 每个网格的流动变量解，也就是全流场的数值解。这就是Oodunov方法或 Godunov格式的基本原理。 计算量非常大，从而阻碍了这种方法的直接应用:另外，它是一阶精度格式， 数值耗散太大。但是Godunov方法的思想却启发了后来CFD工作者。发展了 一些Riemann问题的近似解法以减小计算量、提高计算效率，比如著名的Roe 格式。 考虑Riemann间断分解问题 【，，+， U ， 0 uc工，。， 三，:‘Jx 。0’ 4(11 己， @，，…，Uk 7’，尸 己， Z 材 ，…，以 甜 r 对这个问题进行数值求解较困难，主要是因为由非线性引起的间断分解的 复杂性，而关键又是Jacobi矩阵A U 的非线性。 Roe方法的思路是:利用左右函数的常态数U一口U，(构造一个合理的常矩阵 A u，，U，( ，用来近似代替原来的A U ，这样就将复杂的非线性问题转化为线 性问题，也就是将原来的拟线性方程 20 :海人学硕i:学俯论文 掣+A u 掣:0、1 、 4(12 魂 8x 化为一个线性方程 掣+弛，u， iOU:o 4(13 df dx 其守恒形式为 型+娑: o 、 4(14 Ot 苏 式中户 彳 【，，，U， u。 为了构造这样的常矩阵j U，，U， ，Roe提出了构造这种线性化近似矩阵的 原则，要求A, U，，U， 具有如下所谓的“U’’性质: ?当Uf，U，一U时，A, u，，U， 一彳 U ，即具有相容性; ? 彳 u，，U， 己，，一U， E一巧，即具有守恒性; ? 彳 u，，U，( 的特征向量吼 膏 1，2，(((，聊 是线性无关的，即具有双曲性。 式中，U,(-Gf由下列表达式得到，即 r 一 ,Plll,七,Pru “ l―_ ――二; ，一 r ,pI?,p r 采一,p|vl+,p1(v r 心pl七,p 4(15 疗:型鱼丝?型些: Np|+4pr 万2 7一1 [疗一吉 万2+秽2 ] 上式称为Roe平均。 给出了彳 ?，U， 后，Roe格式的数值通星可以写为 2I 上海入学硕j 学位论文 7L L U，一u， 4??16 t+, 专[F u『 +F U， 卜L，Z,-1 这样就得到了Roe格式中数值通量的表达式。Roe格式一般被称为通量差 分分裂格式，由于其优秀的间断分辨率，因而成为目前应用最广、评价最高的 CFD格式之(―-。 4(2(4熵修正f‘9，27l 通量差分分裂格式 简称FDS 是--种综合性能优秀的计算方法，它的分辨 率很高，是目前CFD界最受欢迎的格式。但是它也具有大于一阶Godunov类 求解器的特点:高的间断分辨率和非物理解出现的可能。当通量Jacobi矩阵的 特征值很小时，会违反熵条件，产生非物理解，故必须引入熵修正。 本文引入了Harten(Yee型熵修正方法，这是一种传统的、经典的熵修正方 法，具体表达式为 f l彳| f五l?万 4(17 jtl 国… 肚j掣 1 【 26 由于万取值的不同，从而形成了一些不同的熵修正公式。最简单的就是直 接取万为常数，一般取为0?万?0(125。但是大多数情况下，直接取艿为常数是 不合适的。对于包含简单非稳定激波的情况时，艿 O，这是因为违反熵条件的 现象仅发生在稳定或近稳定时。对于一般情况，尤其是包含有强激波时，万的 选取十分重要。在效果是，万表现为数值耗散，当万 O时耗散最小。 在笛管尔坐标下，有 万一 1 8 u，+,I+1 式中，万为常数，一般取为O(005到0(4，在高马赫数取大些。 在一般坐标系下，6成为 u y 9 万 8一[II+ll十(导 ?丽+厮 c],+蟛 4(1 式中，U，V为逆变速度分量，万一般也取为O(005到0(4。 淞人学颂一lj学f讧论丈 j 4(2(5限制器【|9'28 Q『 Qi+-z-1仍Ag; P,… 1 。2(4 g:纨,i饵?Q+lf 和左边 nR 通过限制器插值得来。仍和饵是限制器函数，其一般形式为 么 n …(、 佴:如一尼 矽 ‘ ,。+ 1+七 矽 蚓蔗1 黑矧 @2??，P,叫 式中，。 会曼荸，_ 五AQ万i+]。计算格式的限制器精度决定于k的取值，即 k 一1 上风格式 k 0 from,(格式 k ,三阶偏上风格式 血 1 三点中心差分格式 限制器的种类很多，目前最常用的限制器有: ?Van Leer限制器 矽 厂 :娑世 j+， ?Van Albada限制器 (1:海人学彤!Ij学位沦丈 如 等 ?Minmod限制器 如， 一^D黜 ?Superbee限制器 矽 ，( [min 2r，1 ，min r，2 】 ?VanLeer平均限制器 附 :』min 2啦，o+, 厂 o 【0 ，(?0 ?混合限制器 Superbee限制器分辨率好，但容易产生振荡，而一般的限制器对接触间 断的分辨率较差，因此，为了提高接触间断的分辨率，充分利用各限制 器的优势，构造了所谓的混合限制器，对密度使用Superbee限制器，其 余变量采用VanLeer平均限制器。 对上述限制器的性能作一简述:Superbee限制器的耗散最小，Minmod限 Albada限制器， 制器的耗散最大，VanLeer限制器居中;可微性最好的是Van 而Van Leer平均限制器改进了Minmod限制器的耗散特性，混合限制器集中了 Van Leer平均限制器和Superbee限制器的特点。 4(3本文控制方程的格式离散 对于可压缩N―S方程组 4(1 ，采用有限体积法进行全隐式离散，为保证非 定常流动的时间计算精度和计算效率，计算采用双时间步法，其中f为伪时间， t为物理时间。为使得低速流动在时间推进法能够收敛，在伪时间项前面乘上以 密度、速度和压力为原始变量的预处理矩阵r。得到如下方程: r望+M望+丝+堡+箜:塑+箜+翌 4(22 8t 魂 a a芎 rl 8‘8芎aq8‘ 24 w p 为原始变量，转换矩阵 4(22 式中Q 号 专 户”y 。 O o o O 甜 o o 0 p aW ，， o O 4(23 M ――2 v O p a9 O o 一 ， p w叩,