（计算数学专业论文）利用多介质ppm方法数值模拟斜激波与物质交界面的相互作用.pdf

摘要 发展适合于多介质流动的数值方法及其对复杂流动现象的模拟是近几年来的研究 热点之一，比较常用的方法是利用质量分数模型方程。为防止在界面附近的压力振荡而 提出了众多的数值方法。a b g r a l l 通过接触界面两侧的压力平衡和速度一致假设导出流 体体积分数等参数的演化方程，结台e u l e r 方程组，得出了基于准守恒型控制方程组的 计算模型。s h y u e 将该计算方法推广到基于一般刚性气体状态方程的，模型和体积分数 模型的准守恒型e u l e r 方程组。 本文结合a b g r a l l s h y u c 的计算模型和经典的高精度p p m 方法数值模拟了多介质可 压缩流动的情况。一维和二维问题的算例表明，本方法可以较好的模拟高密度比水气两 相介质、带强激波的可压缩流动，对于激波和接触间断分辨率较高。文中重点数值模拟 了斜激波与不同物质交界面的相互作用的演化过程，并给出了交界面上由于切向流效应 产生的涡列的演化过程。 关键字：多介质p p m 方法、斜激波、涡列演化、切向流效应、界面 a b s t r a c t i ti so n eo ft h eh o t s p o t st od e v e l o pt h eh u m e r i c a lm e t h o d sf o rm u l t i - f l u i d sa n dt o s i m u l a t et h ec o m p l e xf l o w s t h ec o m m o n m e t h o di su s i n gt h em a s sf r a c t i o nm o d e l e q u a t i o n s i no r d e rt oa v o i dt 1 1 ep r e s s u r eo s c i l l a t i o nn e a rt h ei n t e r f a c e ，m a n yn u i n e r i c a lm e t h o d sa r ep u t f o r w a t d a b g r a l lo b t a i n e d t h ee v o l u t i o ne q u a t i o n sf o rv o l u m e f r a c t i o no n 也ea s s u m p t i o nt h a t t h ep r e s s u r er e t a i n si ne q u i l i b r i u ma n dt h ev e l o c i t yi sc o n s i s t e n t t h e nh ec o m b i n e dt h e e v o l u t i o ne q u a t i o n sw i 伍t h ee u l e re q u a t i o n s 、a n df i n a l l yg o t 也em o d e le q u a t i o n sb a s e do nt h e q u a s i c o n s e r v a t i v ee q u a t i o n s s h y u ed e v e l o p e dt h en l e t h o dt ot h eq u a s i - c o n s e r v a t i v ee u l e r e q u a t i o n sg r o u n d e do nt h ey m o d e la n dt h ev o l u m e f r a c t i o nm o d e lf o rt h em u l t i f l u i d si n s t i f f e n e dg a se q u a t i o n so fs t a t e i nt h i sp a d e r , as i m p l ef l u l d m i x t u r et y p ep p m a l g o r i t h mf o rm u l t i c o m p o n e n tf l o w , w h i c h c o m b i n e st h e 丑u i d - m i x t u r em o d e lo fa b g r a l l s h y u ea n dt h ec l a s s i cp i e c e w i s e p a r a b o l i c - m e m o d 口p m ) ，i sd e v e l o p e dt os i m u l a t et h ec o m p r e s s i b l em u l t i f l u i df l o w s 1 da n d2 d n t l m e r i c a lr e s u l t ss h o wt h a tt h em e t h o dc a l ls o l y et h em u l t i f l u i d si n t e r f a c i a lf l o w sc o r r e c t l y a n db er o b u s tf o rh i 曲d e n s i t y - r a t i of l o wa n dl a r g ed i f f e r e n c eo fe o sa tb o t hs i d e so ft h e i n t e r f a c e i nt h i sp a d e r , w ea l s op r e s e n tt h er e s u l t so fn u m e r i c a lc o m p u t a t i o n so ft h er e f r a c t i o n o fap l a n es h o c kw a v ea tt i l ed i f i e r e n ti n t e r f a c e su s i n gt h em u l t i c o m p o n e n t sp p m m e t h o d w e p r e s e n t t h ee v o l v e m e n to ft h ei n t e r f a c e sb e c a u s eo ft h eb a r o e l i n i ce f f e c t so i lt h ei n t e r f a c e e s p e c i a l l yt h ec a s eo f t h es 拄o n gs h o c ks t r i k i n ga tt h ei n t e r f a c e sa sw e l l ， k e yw o r d s ：m u l t i c o m p o n e n t sp p m ，i n c l i n e ds h o c k , e v o l u t i o no ft h ev o n e x b a r o c l i n i ce f f e c t s i n t e r f a c e 独创性声明 y 6 6 2 9 8 9 本人声明所呈交的学位论文是本人在导师指导下进行的研究工作及取得的 研究成果。据我所知，除了文中特别加以标注和致谢的地方外，论文中不包含其 他人已经发表或撰写过的研究成果，电不包含为获得中国工程物理研究院或其他 教育机构的学位或证书使用过的材料。与我一同工作的同志对本研究所做的任何 贡献均已在论文中作了明确的说明并表示谢意。 学位论文作者签名 五景告 签字日期：御u 年r 月i 日 学位论文版权使用授权书 本学位论文作者完全了解并接受中国工程物理研究院研究生部有关保存、使 用学位论文的规定，允许论文被查阅、借阅和送交匡家有关部门或机构，同时授 权中国工程物理研究院研究生部可以将学位论文全部或部分内容编入有关数据 库进行检索，可以影印、缩印或扫描等复制手段保存、汇编学位论文。 学位论文作者签名 至景争 导师签名 住啉军 签字日期：聊毕年r 月日 签字日期：如牛年岁月r 1 日 利用多介质p p m 方法数值模拟斜激波与物质交界面的相互佳旦 第1 章绪论 计算流体力擎的诞生和发展是和电子计算机的诞生和发展密切相关的。计算流体力 学的开创者v o nn e u m a n n 也是近代计算机的创始人之一绝不是偶然的。近十年来计算流 体力学经历了一个蓬勃发展的时期，在计算机科学和超级计算机研制，计算数学和数值 方法以及它在计算流体和工程中应用都取得了实质性进展。目前，计算流体力学已无可 争辩地确立了自己应有的地位。理论、实验和数值模拟已被普遍认为是科学研究和应用 开发的三大手段1 1 。 1 1 界面失稳研究的意义 界面失稳是指密度不同或速度不同的两股流体相遇后由于小扰动所引起的 复杂流动现象，它包括k e l v i n f i e i m h o l t z 失稳、r a y l e i g h t a y l o r 失稳和 r ic h t m y e r m e s h k o v ( 简称r m ) 失稳。本文的研究对象是r m 失稳，它是指当激 波以一定的速度扫过密度分界面后，加速界面所导致的不稳定性。激波与密度分 界面的作用有着深厚的实用背景。激波扫过密度界面后导致界面两边流体混合增 加的机理可直接应用来增加航空发动机燃料与空气的混合从而提高燃烧效率；而 对激波导致界面失稳的研究则有助于对激光控制热核聚变的了解。这一切证实了 对界面失稳的研究具有重要的意义。 对激波与密度分界面相互作用的问题，数值模拟和实验研究是两个主要的研 究手段。实验研究具有客观性和可靠性等优势，而数值模拟的优点则体现在可重 复性以及细致性。由于数值模拟条件的限制，使得人们对算法提出了较高的要求。 由于流场中出现强间断，为了防止在流场中产生非物理振荡，我们必须在算法上 采取相应的措施。 1 2 模型方程提出的背景 发展适合于多介质流动的数值方法及其对复杂流动现象的模拟是近几年来 的研究热点之一。这些问题通常是所谓的接触间断问题，对它们的数值模拟要比 对强间断的激波问题更为困难。特别是要有精细而锐利的高分辨率图像实为不 易，因而对它的数值研究和分析受到目益强烈的关注。一个比较常用的方法是利 用质量分数模型方程，即在g u l e r 方程组的基础上加上描述各流体组分运动的守 恒输运方程进行求解。但是在处理多组分的流体混合时，这种方法特别是利用激 波捕捉的方法，在多介质界面附近的网格上不能保持压力平衡，从而会产生一系 型岂童坌重! 型查鲨塑堕蔓型塑堂婆量塑堕奎墨亘笪查坚型兰塑一 列的问题 2 。为了防止在界面附近的压力产生振荡而提出了众多的数值方法， 例如基于y 模型和1 e v e l s e t 模型的数值方法 3 ，以及基于质量分数模型的修正 方法 4 ，5 。 a b g r a l “6 ，s a u r e l 7 和s h y u e 8 基于界面上的物理运动规律和能量方程， 提出了一种多介质可压缩流动的简单方法，并可以处理一般的可压缩材料 9 。 a b g r a l l 6 通过接触界面两侧压力平衡和速度一致假设导出流体体积分数等参 数的演化方程，结合e u l e r 方程组，得出了基于准守恒型控制方程的计算模型。 该模型保证了界面附近能量方程的致近似并避免压力振荡，从而成功解决了多 介质流动界面模拟中的问题。s h y u e 8 推广该计算方法到基于一般刚性气体状态 方程的，模型和体积分数模型的准守恒型e u l e r 方程组，形式如下： 娑+ 昙( 删) = o 昙( p z ，) + 昙( p “! + p ) = o d f0 咒 暑( 加+ 杀n 洲- o 鲁( 击卜嘉 吉卜 鲁( 告卜旦a 。f 逸r - u = 。 _ o p + 昙( p f ) ：0 昙( ) + e 。2 ( p “2 + p ) = o 导( p e ) + - i ( p e + p ) “ ：o 甜u o 型+ 。型：u ni 一1 2 m 一十“一2 f = 上面两式分别对应于y 模型和体积分数模型的模型方程组，具体的推导过程在后 文中会有叙述。从原则上讲，由于该模型方程组比较简单而且适用范围广，可以 充分利用目前发展成熟的高精度的求解双曲守恒律的激波捕捉方法进行求解，例 如t v d 格式，m u s c l 格式等等。在文章中，我们使用了p p m 方法进行求解。 1 3 高分辨率格式的发展概况 一维非定常流体运动方程组是拟线性双曲型方程组，因而不管初值如何光滑 解都可能会出现间断。长期以来，这种特性使求解流体力学方程组有其特殊的困 难。自1 9 5 0 年v o nn e u m a n n 和r i c h t m y e r 提出人工粘性法 1 0 的差分方法计算 流体力学的激波问题以后，人们一直在高精度高分辨率格式的研究上进行大量的 工作。高阶精度是指在连续可微区域中差分解与精确解相比所达到的精度；高分 辨率是指差分解在表示解的间断面时所表现的明显程度。 v o nn e u m a n n 和r i c h t m y e r 提出的人工粘性法是最早的激波捕捉法，他们的 2 型旦童坌堕! ! 塑查鲨墼笪燮型塑鲎茎兰塑星窒墨耍塑塑呈笪塑l _ 一 工作也标志着激波捕捉法研究的开始。其思想是在描述流体运动的拟线性双曲型 方程组中加上一项人工粘性项，其作用是将运动中出现的激波间断抹光滑，使它 变成有一定宽度的连续的过渡区，这样就可以用有限差分方法进行数值模拟计 算。数值结果中物理量梯度变化剧烈的地方就表示激波所在之处。用这种人工粘 性法进行计算，激波过渡区一般占据五个空间网格步长，在激波波阵面后会出现 非物理的数值振荡。 g o d u n o v 于1 9 5 9 年提出的间断分解方法是有效处理激波、接触间断等间断 问题的有效方法之一，它把偏微分方程的解当成是分段光滑的，然后在每个区域 的边界解一个一维r i e m a n n 问题，并利用迎风格式得到下一步的解。具体地说， g o d u o v 方法的基本思想是：将。时刻已知的离散分布量 蠕，： ， 店【2 ， 躁t ，： 分别 看成是每个小网格区间x ， z 0 s ， 。 l甘_ 垒僖押。 ”7 以保证略，：落在町和畋，之间。 兰给出a l , y 和口一。后，在大多数情况下，由( 3 3 ) 得到的a l , j 和d 蹦，就是二次多项 式( 3 t 2 ) 中的两端值。但在以下情况下，由上式确定的抛物线分布可能在某些区域的取 值不在吼，和a r , y 之间，所以还要对口幻和“蹦重新赋值进行调整： ( i ) 当白材一a ，n ) i ，n - a s , g ) 蛆“， w j2 1m l n 【p ，p h ) 7 0其他情况 j = m a x ( o ，国1 ( c o “一e ) ) 从而系数，的计算 其中西：! 兰二呈! 二! ， 。p ，+ 2 一p ，一2 o o p 乃 一 一 鼽缈 l 一 。，、l = ， s 利用多介质p p m 方法数值模拟斜激波与物质交界面的相互作用 f = m a x ( z + z ，z + 。)其峨= m i n ( w ，q ，) 蜥，= 半孵f 些害芋型 在实际计算中，可以使用文献 1 2 中给出常数的值： 盯( 1 ) = 2 ，芷( 2 ) = 0 0 1 ，( 1 ) = 0 7 5 ，彩( 2 ) = 1 0 ，s = 0 3 3 3 1 2 l a g r a n g e 坐标的流体力学方程组 考虑l a g r a n g e 坐标的流体力学方程组 r ：土为比容。 口 塑一她：n 8 ta m 尝嘉= g s ， 丝3 t + 型o 塑m = 耀 这里，r = 去为比容，崩密度，“为速度，p 为压力，e = p + 昙甜：为总能，。状态方程取 口 1 “ p 2 ( y 一1 ) p e ，g 为体积力，r 为空间坐标，r f l 为l a g r a n g e 质量坐标，它们之间的关系为 = f p ( r ) r 。d r + ，t 为时间，对t 的求导为l a g r a n g e 导数a d o t ：“，扛：o ，1 ，2 分别对应于 平面对称、柱对称、球对称的一维情况。 - 令d x m ，为- ，网格的质量，j 网格两端的质量坐标为： ”1 ) - 1 1 2 ，“2 ，z 眈。w 2 一_ 2 这些量不随时间变化。 现设给了，“时刻的网格均值分布： u ；72 i i ：o 州勺砌 箸 oxa盯 i i 盯 利用多介质p p m 方法数值模拟斜激波与物质交界面的相皇塑 其中u ；( ，。e ) ，m 川j ：m 。，以及网格端点的空间坐标以。 。要计算 r s ， f ”1 = f “+ a t 时刻的网格平均值及端点位置。 p p m 方法的计算步骤如下： ( i ) 由网格平均值f ? ，“? ，e ? ，用状态方程计算压力的网格平均值p ? 伍) 应用线性方程所描述的构造分段抛物线函数的方法，分别对所给的网格均值分布 钟) ，蟛 ， 矿) ，构造相应的分段抛物线函数，确定有关参数。这只涉及标量的插值算 法，与描述这些量的方程式无关的。且在l a g r a n g e 坐标中，接触间断是自动保持的， 所以也不必进行接触间断的处理。 ( i i i ) 计算网格端点m = m 川，2 两侧物理量( f ，“，p ) 的时间平均值( i m 亿l ，i m 2 ，多2 c ) 及 ( a j 具甲( 7 川亿l ，“，2 ，l ，p 川悒lj 表不在j 删格的c + 族特征线d d t = r 。肛影响f 的时 间平均值。取其特征速度为 旷h l = n 罴+ 1 n “+ 1 聊 于是， “m 厂，一i 1 画a a t 卜v 飞厂筹卜， “2 _ - “n ，一i 面卜一”u li 。了意卜i ， 。，芦。，的计算公式形式与它完全相同。 ( b ) ( 川忆w ，破川亿n ，p 川亿。) 表示在j + 1 网格的c 一族特征线驯m = 一，4 p c 影响下的时 j 司平均值。取其特征速度为一q 于是 扎m ，一一；筹 d r , l + - - x l , , + l - - ，2a 耐+ 1 a = u l j + i b l r , + l - - 1 x i , , ! + 1 1 - 町+ ， ， 叶崛8 一j 蔷。3 锄。j zj i 川缇。，p 川亿。的计算公式形式与它完全相同。 ( i v ) 通过求解两侧状态为( t “：i 川忆。，p 川缇。) 和( 乇。：用川偶。，p j + i 2 , r ) 的r i e m a r m 问 题来计算在肌= t q l j r 1 , , 2 处的u ，p 的时间平均值t + ：，万。， ( v 讲算r ”1 时刻的分布。 利用多介质p p m 方法数值模拟斜激波与物质交界面的相互作用 _ 誊2 = ”+ i 2 + m j “2 了? ：广一。广“ “”“2 ( 口+ 1 ) ，瓦2 ，。：盥坐二丝坐： 。 f 口+ 1 ) a m j “y = “；+ ；( 巧。：+ 刁。：石a 万t ( p _ h ，：一瓦。：) + 了a t ( g ，n + g ，n + 1 ) e = e j + s 篇- - - ，( a j “：乃w z 乃w ：一万“：乃w ：乃。：) + 丁a t ( “，n g ，n + “，n + l 毋n + 1 ) 3 2 一维多介质p p m - l r 格式 考虑将一维准守恒型多介质e u l e r 方程组( 2 1 2 ) 写成如下的l a g r a n g e 形式，这 时附加的输运方程可以简化为相应变量的随体导数为0 。 塑一旦! 塑：o o to m 坐+ ，a 里：0 o td ，卵 o 西e + o ( r a 。m u p ) ：。 ( 3 _ 9 ) d y 西 = 0 方程组中f ，“，e ，1 分别为比容、速度、l - l 总能、i 组分的体积分数，d 在平面、轴 对黧球对称时铡帅l ，2 。巨夕，为空间蛳满足象叫蚍对 刚性气体状态方程，密度、压力、忱容、比内能为 多= 1 f ，p = ( y 一1 ) p e 一7 p 。，p = e = u2 2 。 空间坐标r ，体积坐标v 和质量坐标m 之间的关系为 d m = l - - i x 。d x = f d v 。 设f ”时刻第j 个网格流动变量u = ( f ，“，e ，f ) 的质量平均为 叼2 麟u ( ) 咖。 一维l a g r a n g i a n - r e m a p p i n g 型p p m 方法分为四步： 利用多介质p p m 方法数值模拟斜激波与物质交界面的相互作旦 ( 1 ) 物理量的分段抛物插值； ( 2 ) 近似r i e m a n n 问题求解网格边界时间平均的压力和速度； ( 3 ) l a g r a n g e 方程组推进求解； ( 4 ) 最后将物理量变回静止的e u l e r 网格上。 物理量的插值与前面所描述的过程一样，这时需要插值的网格平均量为 p ? ，“? ，p ? ，( j ，r ；而且这时无需利用接触间断的检测，n n 捌, 2 格朗日步的计算中，接 触间断是自动保持的。 为获得网格边界上时间平均的压力和速度，p p m 方法需要在网格边界上求解近似 r i e m a n n