弱可压缩流体流动模型及数值求解

第１　８卷第１期　抽ｎ“０ｆ　长Ｙ￣ｇｉ江ｚｅ　Ｒｉ科ｖｅｒ　学知Ｒｅ院　院拒　Ｉｎｔａｉｔｕｔｅ　．１８　Ｎｏ　１　２　０　０１车２月　Ｆｅｂ　２　０　０　１　文章编号：ｌ００１—５４８５｛２００１　ｌ０１—０００６—０４　弱可压缩流体流动模型及数值求解　于波　（武汉水利电力大学动力工程系，湖北武汉４３００７２）　摘要：介绍了小马赫数下弱可压缩流体流动的模型，对该模型的数值求解方法也作了较为详尽的叙述。讨论了紊　流模型模拟涡流的方法。简要介绍了水力机械内部流场的计算实例。　关键词：有限体积法；可压缩；数值模拟　中固分类号：Ｏ２，ｔ２　１　文献标识码：Ａ　由于水流中声速和对流速度在数值上相差较　ｐ　Ⅱ　·“　＝０；　（４）　大，因此小马赫数下的不稳定可压缩流动很难计算。　假定流动为不可压缩，将方程简化，那么，推导出来　Ｓ　＋（“　·　）“　＋　＋　ａ￡　０　的方程不再适用于描述与时间紧密相关的流动（例　如水力机械的水力过渡过程）。而且不可压缩流动　Ｍ　南　０，（５）　方程组中的压力项Ｐ在连续方程中没有出现，数值　方程中：Ｍ＝ｕｏ／ａ０，是马赫数；Ｓ：ｌｏ／（ｔｏ“０），是斯　求解时处理比较麻烦。基于上述原因，美国Ｃｈａｒｌｅｓ　特若哈数；＊表示无量纲化的量。如果Ｍ很小，方　Ｃ　Ｓ．Ｓｏｎｇ教授提出了弱可压缩流体流动模型。该　程（４）中的第２项和方程（５）中的最后１项可以忽略　模型介于可压缩和不可压缩之间，既保留了真实流　不计；如果Ｓ较大，方程（４）中的第１项不能忽略，　体可压缩的特性，又可使方程组简化，易于求解。　比如计算水力过渡过程或水声传播的情况时就是这　样。因此，对于弱可压缩流体无粘性流动，其控制方　１基本方程的推导　程可简化为　＋ＫＶ·“：０；　（６）　对于正压流动非粘性流体的连续方程和运动方　程可写为　譬＋（ｕ　“．　）“＋　卫＝０，　（７）　警＋（ｎ·　）　十ｐａ　Ａ·ｎ＝ｏ；　（１）　式中：Ｋ＝　是容积弹性模量，Ｐ和。为已知常数。　．△）　詈＋寿Ａｐ：０＇（２）　对无粘性的弱可压缩流体流动而言，如果Ｍ　较小，方程（６）和（７）可达到Ｍ的一阶精度，如果流　式中：Ｐ为压力；“为速度；ｎ为声速，以下式表示　动是稳定的，则连续方程退化成不可压缩流动，即　ｎ２＝　０　（３）　０。　（８）　只要Ｋ不为零或无穷大且为有限值，弱可压缩　如果流体不可压缩，则ｎ一一，方程（１）中的前２项　流动方程就可能被用于计算稳定的不可压缩流动，　和方程（２）的最后１项没有。弱可压缩流体流动被　而不管Ｋ或ｎ假定为何值。用原始变量法进行数　定义为：ｐ和ｎ都是变量且改变都非常小，以至于有　值求解时，压力梯度是动量方程中源项的重要组成　时可以把它们当作已知常数对待。　部分，但没有压力方程用于计算，求解压力场有较大　方程（１）和方程（２）中的变量可以用ｐＤ　５，“０，　ｐ０，　０，ｆ０等参数无量纲化，其结果为下列无量纲方　困难。而弱可压缩流动可对压力项导出显式的求解　程：　方程，这是与其它方法的显著区别。　ＳＭ２＋告＋Ｍ　（“　·　）　＋　收稿日期：ｌ９９９．１１—２４　作者简介：于渡（１９５４一）．男，河北威县＾．武汉水利电力大学动力工程系教授，主要从事流体机械动力特性和水电机组ＣＡＤ硬ＣＡＴ技术　研究．　维普资讯 http://www.cqvip.com

第１期　于　渡　弱可压缩流体流动模型及数值束解　定了解的唯一性。求解方程（１Ｏ）采用Ｍａｃ　Ｃｏｒｍａｃｋ　２弱可压缩流体三维非定常流动欧拉　方程的数值解　描述流体一般情况的运动方程常用Ｎａｖｉｅｒ－　Ｓｔｏｋｅｓ方程。该方程能适用于实际中任何流动问　题，但求解较为困难。欧拉方程提供了求解比较困　难问题的有效途径，如果处理适当，欧拉方程可以用　来描述任意马赫数下三维非定常流动，且不需作无　旋（ｖ·口＝０）条件假设。　２．１计算网格　网格中的一个有代表性的控制体如图１所示。　３对正反表面很　明显且由下标　ｌ，２，３标识。８　个顶点中的１个　称为正角，相交　图１控制体示意图　Ｌ　的显式预估．校正法。根据方程（１０），在已知当前时　间Ｆ值的情况下求解下一时间的Ｇ。求解方法的　显式特性表明，在给定时刻，所有Ｆ值为已知并且　等可由１个方程求出（隐式方法需要求解联立方　程）。因为计算的时间步长极小，而且只知道某些点　上的Ｆ值（不是连续分布），所以方程（１０）的时间导　数和表面积分只能是近似值。预估．校正法的本质　表明方程（１Ｏ）被使用了２次，第１次预测出Ｇ和Ｆ　在后一时刻的值，然后通过预测值校正Ｇ和Ｆ在该　时刻的值。这些校正值对时间和空间步长来说大致　为二阶精度。这对大多数应用来说已足够了。　２．３控制方程的离散及数值求解　数直离散方法用ＭａｃＣｏｒｍａｃｋ两步显式迭代　法。时间步长用Ｃｏｕｒａｎｔ准则确定。由基本方程　于这个正角的３　（９），式中Ｇ代表流动变量，Ｆ代表通量矢量。为得　到第一时刻整个计算区域内Ｇ的值，阐述求解Ｇ　的过程如下：　Ｆｉｇ　１　Ｓｋｅｔｃｈ　ｅｏｎｔｍｌｌｉｎｇ　ｂｏｄｙ　个表面的面积矢　量称为ｓ　，ｓ　，ｓ　，另外３个表面的面积矢量称为　ｓ　，ｓ　，ｓ　。当２个有限体有公共面时，一个有限　体的面积矢量称为ｓ　，另一个有限体的面积矢量　称为ｓ　。这种偶联形成３个“伪方向”而贯穿计算　网格，并能够不依赖于关于坐标轴的网格方向选择　另外３个下标对流动变量进行分类。　２　２控制方程　弱可压缩流体基本方程（６）、（７）可统一写为　假设在某些时刻ｔ时已知每一网格点Ｇ的值，　后一时刻Ｇ的值可通过Ｔａｙｌｏｒ级数得到，即　Ｇ（ｔ＋Ａｔ）＝ａ（ｔ）＋（　）　Ａｔ十　专（等）　ｚ＋Ｉ＿．，　对于有限时间步长，方程（１１）可近似地写为　Ｇ（　＋　）＝Ｇ（　）＋（　）　￡，…）　（１２）　式中导数项的下标ｆ代表计算时刻ｔ时的导数值。　＋Ｖ·Ｆ＝０，　（９）　式中：Ｇ是矢量，代表各独立流动变量　，“，　，Ｗ；　Ｆ为通量度量。如果方程（９）被应用于一个很小的　固定的有限体，比如Ｖ的网格，则可写为　式中（等）　是从ｆ到￡＋Ａｔ时刻导数的有效值或　平均值。寻找新的Ｇ值的问题变为寻找（等）　，　这可通过方程（１１）来求解。方程（１１）中的矢量Ｆ　ｕ　ｄ　ｒ　Ｉ　．Ｆｄ　（１０）　和流动变量Ｇ只是简单的代数关系。如果已知时　刻ｔ时的Ｇ值，那末同一时亥４各处的Ｆ值也就知道　了。为了在离散点上应用方程（１１），有必要做一些　式中：ｎ是包围体积ｖ的表面Ｓ的单位外法线；Ｇ　取有限体中的平均值，该平均值用中心点的值来近　似代替　方程（１０）是欧拉流动的控制方程，它描述了三　维非定常流动，但常常被用于计算定常流动。对Ｇ　和Ｆ需要给出初始条件，使用方程（１０）可确定Ｇ的　时间变化率，而后可求得Ｇ的值。允许时间一直向　前推进直到某一时刻Ｇ的值与预期的收敛值十分　处理。将方程（１１）在一个小的、固定的网格有限体　上积分，则有　Ｊ　警ｄＶ—Ｊ　Ｖ·ＦｄＶ。　（１３）　考虑等式左边，因为假定有限体是固定的，所以时间　导数项可以移出来。应用积分中值定理，可得　接近为止。通常假设上游为均匀自由流动并且在径　向上远离物体，给出物体表面一个特定的法向速度　分量（通常为零）和一些下游条件。这些边界条件决　』ｖ　３。Ｇ　ｄｖ：　Ｖ＝　（　』ｖＧｄＶ）＝　（１４　（　Ｖ）＝Ｖ等，　维普资讯 http://www.cqvip.com

８　长江科学院院报　２００１丘　式中石是Ｇ在有限体ｖ中的平均值。如果有限体　Ｖ很小，则　可近似地认为和有限体中心的值一　样，因此有　式中通量积分由某一合适表面后面的通量矢量逼　近。方程（１７）中的　（ｔ十　）现在可用来求通量矢　』　ａＶ＝Ｖ誓。　）　Ａｔ）：　量的预估值，即　（ｔ＋Ａｔ）。这些新的通量矢量被　用于预估导数值，而该导数值又用来修正　（￡＋　方程（１３）和（１４）是通过严格的数学方法由方程（９）　推导而来，方程（１１）亦是如此（Ｇ是可微分、积分　的）。方程（１２）和（１５）引人了近似条件。方程（１３）　的等号右边，应用散度定量，将体积积分变为表面积　分，为　（著）一专（　均值问题。　ｓ：＋　Ｆ　）　２）　这里通量积分由某一合适表面后面的通量矢量逼　近，方程（２２）是校正步，并解决了前面提出的通量平　ｆ　．ＦｄＶ：ｆ　ｓ，　（１６）　至此可总结求解的数值算法如下：基本的数值　算法是Ｍａｃ　Ｃｏｒｍａｅｋ的显式预估　校正方法。将方　程（１Ｏ）对通量积分，用面积矢量的标量积和恰当地　式中ｎ是Ｓ的单位法向矢量。将源于方程（１３）的　方程（１５），（１６）除以体积ｖ，得到　ａｆ　Ｊ』　…。　．Ｆｄｓ。　Ｖ（…　１７）　估计通量矢量来近似逼近。方程的解从时间　到　＋１向前推进，即　方程（１７）提供了一种用于近似求解时间ｔ时有限体　中心处Ｇ的导数方法。　若已知Ｇ（ｆ），求Ｇ（￡＋Ａｔ），可由方程（１２）求　Ｇ”　＝　‘　“＝Ｇ　’　（∑Ｆ－ｓ；＋∑Ｆ一－ｓｉ）；　＾　＋（１一　）　一　Ｖ（　●　ｓ　＋　（２３）　得。但该方程需要求解（等）　，而该项此时还未被　求出。不过，特定的时间导数可通过方程（１７）求得。　有理由认为，为了达到某一特定阶数的精度，以Ｍａｃ　Ｃｏｒｍａｃｋ中步骤的有效性和精确性来说，其基本上　有二阶精度。　Ｆ－ｓｉ），　式中：Ｆ是有限体中心处流动通量；Ｆ　是和ｓ　共　面的有限体中心处流动通量；Ｆ类似地和ｓ　对　应。　公式（２３）中引入了参数　帮助使数值方法稳　定。　＝１时，公式退化为Ｍａｃ　Ｃｏｎｎａｃｋ数值方法　方程（１７）中积分决定了定解问题的通过有限体　的净通量　因为每个有限体有３对相对的面．所以　上述积分可认为是６个面上的积分之和，即　的有限体通量形式。增加　值就会增加截断误差，　起到了人工粘性的作用。　时间步长　受Ｃｏｕｒａｎｔ－Ｆｒｉｅｄｒｉｃｈ—ｌｅｗｙ稳定判　Ｊ　Ｓ　唰ｓ　唰ｓ＋　唰ｓ，　（１８）　据限制，所允许的最大时间步长计算式为　１，　这样的分解是精确的。但是，每个面上的Ｆ值并不　知道，只知道有限体中心处的Ｆ值，因此各项积分　只能采用数值逼近。Ｍａｃ　Ｃｏｒｍａｃｋ法基本上采用　Ｆ和Ｆ的平均值。　预估一校正是为了求Ｇ（ｆ＋Ａｔ），首先预测，然　后该预测被修正。有如下定义：　０Ｇ）　Ａｔ　是速度矢量。　Ｔ　可ｃ　（２４）　全体网格中，对每一个有限体的第一个表面，ｑ　２．４边界条件和初始条件　边界条件是通过使用计算网格外部的映像点来　加上去的。例如自由流动给定时，其值出现在影像　壁或物体表面边界，一般也使用镜像条件。当边界　吉［（著）　＋（　）］。（１９）　点上。在轴对称平面条件下，可用镜像条件。在边　任意几何形状的表面目前还没有找到一种有救的方　上式（等）　预估一个值，然后（等）修正这一结果。　不是平面时，这将导致计算结果的不精确。然而，对　Ｇ（ｔ＋Ａｔ）的预估值为　Ｇ＿（ｆ＋△　）＝Ｇ（　）＋（等）　由方程（１７），（１８）可求得　（２０）　法解决常规或参数型运动方程的边界条件问题。　（等）　一专（　Ｆ·ｓｌ＋　ｓ　（２１）　３关于粘性流动的讨论　若考虑粘性流动，动量方程可写为　维普资讯 http://www.cqvip.com

第１期　于波　弱可压绾流体流动模型及数值求解　９　等＋　一！　ｖｐ＝　中７３＝　紊流模型。　（２５）　４水力机械内部流场计算实例简介　文献［１］应用弱可压缩流动理论，采用大涡模拟　方程的离散及求解类似于非粘性流动的情况，　但对于粘性项需要专门考虑。可令ｆ＝　ｖ　，式　，　为层流运动粘性系数，　为紊流　运动粘性系数。确定紊流运动粘性系数，需要讨论　研究表明，湍流流动由不同尺度的旋涡组成。　大尺度的旋涡对湍流能量和雷诺应力的产生以及对　各种量的湍流扩散起主要作用，大涡的行为强烈地　依赖于边界条件，它随流动类型而异。小涡主要对　紊流模型，对射流泵内部流动（二维）进行了计算。　计算过程是用非定常解逐步逼近定常解。网格划分　为（６ｏ×１２）＋（６０×１２），马赫数％＝０．２。计算得　出了射流泵流场流速分布和轴向压力及速度变化。　文献　推导了基于弱可压缩流动的Ｎａｖｉｅｒ—Ｓｔｏｋｅｓ方　程，同样也采用大涡模拟紊流模型，数值模拟了混流　式水轮机从转轮到尾水管的整个流动过程，得到了　耗散起作用，在高雷诺数下小涡近似于均匀各向同　性，受边界条件的影响较小。应该说，雷诺平均的处　理方法并不能真实反映湍流流动的上述基本特点，　因为由此出现的雷诺应力项包含了全部大涡引起的　脉动，因而从根本上说难于用一种通用的湍流模型　去描述强烈依赖于边界条件和流动类型的大涡的行　转轮叶片正、反面的压力分布，尾水管各截面的流速　矢量图，尾水管压力恢复系数曲线等重要数据。对　水轮机的主要特性而言，预估值（计算值）与试验数　据比较吻合。　基于大涡模拟的弱可压缩流动模型对流动的模　拟功能更强，可以记录速度、压力、涡量等与流动有　为。更合理的方法是用计算机计算大涡而用模型　（即所谓亚格子ＳＧＳ模型）模拟小涡，这就是大涡模　拟的基本思想。大涡模拟的主要优点是大涡部分可　以精确求解，整个计算受人工模型的影响比较小。　计算模型为　关的参数的瞬时值和时问历程，能较为精确地预测　过流部件的流速场、压力场，转动部件和固定部件之　间的相互影响以及流动过程中产生的水力损失。这　些对揭示水力机械内部的流动机理具有重要意义。　（ｅ△）　（２ｓ　Ｓ　）　；　（２６）　参考文献　陆宏圻．流体机械与流体动力工程　Ｍ］．武汉：武汉水　利电力大学出版社．１９９８　ｓ　＝吉（誉＋　），　（２７）　…　式中：Ｃ为子涡扩散系数，对不同的问题，ｃ　＝０　１　０．４；△为能直接计算的大涡流最小尺度，与网格　大涡模拟需要网格足够细，否则不足以反映紊　流脉动详情，这样必然导致计算工作量巨大，一般微　的尺度有关。　粘性流体动力学基础［Ｍ］北京：高等教育出　［２］　陈矛章．版社，１９９３　Ｓ，ＹＵＡＮ　ＭＩＮＧＳＨＵＮ　Ａ　ｗｅａｋｌｙ　ｔｏｎｉ—　［３］　ＣＨＡＲＬＥＳ　Ｃ　Ｓ　ｐｒｅｓｓｉｂｌｅ　ｌｆｏｗ　ｍｏｄｅｌ　ａｎｄ　ｒａｐｉｄ　ｃｏｎｖｅｒｇｅｎｃｅ　Ｍｅｔｈｏｄｓ［Ｊ］．　ｄ￣ｕｒｎａ［ｏｆ　Ｆｌｕｉｄｓ　Ｅｎｇｉｎｅｅｒｉｎｇ，１９９８，１１０：４４１—４４５．　机和小型机都难以满足要求。　（编辑：聂文）　Ｂｒｉｅｆ　ｐｒｅｓｅｎｔａｔｉｏｎ　ｏｎ　ｗｅａｋｌｙ　ｃｏｍｐｒｅｓｓｉｂｌｅ　ｆｌｏｗ　ｍｏｄｅｌ　ａｎｄ　ｉｔｓ　ｎｕｍｅｒｉｃａｌ　ｓｏｌｕｔｉｏｎ　ＹＵ１３ｏ　（Ｗｕｂａｎ　Ｕｎｉｖｅｒｓｉｔｙ　ｏｆ　Ｈｙｄｒａｕｌｉｃ　ａｎｄ　Ｅｌｅｃｔｒｉｃ　Ｅｎｇｉｎｅｅｒｉｎｇ，Ｗｕｂａｎ　４３００７２，Ｃｈｉｎａ）　Ａｂｓｔｒａｃｔ：Ｉｎ　ｔｈｅ　ｐａｐｅｒ，ａ　ｗｅａｋｌｙ　ｃｏｍｐｒｅｓｓｉｂｌｅ　ｆｌｏｗ　ｍｏｄｅｌ　ｆｏｒ　ｓｍａｌｌ　Ｍａｃｈ　ｎｕｍｂｅｒ　ｉｓ　ｉｎｔｒｏｄｕｃｅｄ；ｔｈｅ　ｍｅｔｈｏｄ　ｏｆ　ｍｕｍｅｒｉｃａｌ　ｓｏｌｕｔｉｏｎ　ｉｓ　ｇｉｖｅｎ　ｉｎ　ｄｅｔａｉｌ；ｔｅ　ｈｍｅｔｈｏｄ　ｆｏｒ　ｓｉｍｕｌａｔｉｇ　ｖｏｒｔｎｅｘ　ｌｆｏｗ　ｉｓ　ｄｉｓｃｕｓｓｅｄ；ａｎｄ　ｃａｌｃｕｌａｔｉｏｎ　ｉｌｌｕｓｔｒａ－　ｔｉｏｎｓ　ｏｎ　ｆｌｏｗ　ｆｉｅｌｄｓ　ｉｎｓｉｄｅ　ｈｙｄｒａｕｌｉｃ　ｍａｃｈｉｎｅｒｉｅｓ　ａｒｅ　ｂｒｉｅｎＹ　ｒｅｐｒｅｓｅｎｔｅｄ．　Ｋｅｙ　ｗｏｒｄｓ：ｆｉｎｉｔｅ　ｖｏｌｕｍｅ　ｍｅｔｈｏｄ；ｍｍｐｒｅｓｓｉｈｉｌｉｗ；ｎｕｍｅｒｉｃａｌ　ｓｉｍｕｌａｔｉｏｎ　①ＣｈａｄｍＣＳＳ．Ｈｅ　Ｊ】ａｎ　，ＣｔｅｌｌＸｉａｎｇｙｉｎｇ　Ｕａｋｕｌａｔｉｏｎ　ｏｆｔｕｒｂｕｌｅｎｔ　ｎ　ｔｈｒｏｕｇｈ　ａｆｒａｎｃｉｓｔｕｒｂｉｎｅ　ｒｔｍｎｅ￣ａｎｄｄｂｏｗ　ｄｒａｆｔｅｒ　ｅａｔｅｒ　ＡＳＭＥ　ｐａｐｅｒ　Ｉｍｅｍａｆｉｏｎ￣ｌＰｏｗｅｒＧｅ￣ｅｒａｔｋｍ　Ｃａｆｆｅｒｅｎｃｅ．ＵＳＡ：ＡＳＭＥ，１９９１