（计算数学专业论文）多介质流体计算的laxfriedrichs格式.pdf

摘要 多介质流体的计算是近年来计算流体力学的一个热点相较于单介质的流体 它的主要难点在于介质界面的处理上很多对于单介质流体十分成功的计算格式应 用到多介质流体上时，在介质界面处都产生伪振荡如何消除这些不必要的振荡也 就成了算法研究的重点 在这篇论文中，我们设计了一个基于非平衡态的l a x f r i e d r i c h s 格式用以对多介 质流体进行数值模拟我们构造的格式保证了多个守恒律：总质量和分质量守恒， 总动量和总能量守恒它还可以保证分质量均非负，更重要的是它消除了速度和压 强在介质界面处的伪振荡数值例子表明这一算法是有效的 关键词： 多介质流，伪振荡，非平衡态 a b s t r a c t n u m e r i c a ls i m u l a t i o no fm u l t i f l u i di sah o ts p o ti nc o m p u t a t i o n a lf l u i dd y n a m i c s c o m p a r e dw i t hs i n g l ef u i dc o m p u t a t i o n ，t h em a i nd i f f i c u l t yo fm u l t i f l u i ds i m u l a t i o ni s t ot r e a tt h em a t e r i a li n t e r f a c e s s c h e m e st h a ta l es u c c e s s f u li ns o l v i n gs i n g l ef l u i df l o w s w i l lp r o d u c es p u r i o u so s c i l l a t i o n si l e a l m a t e r i a li n t e r f a c ei nt h el a t e r h o wt oe l i m i n a t e t h e s es p u r i o u so s c i l l a t i o n st h e nb e c o m e st h em a i nc o n c e r ni nd e s i g n i n gg o o dn u m e r i c a l m e t h o d s i nt h i sp a p e rw eb u i l dl a x - f r i e d r i c h sas c h e m eb a s e do i ln o n e q u i l i b r u mt h e r m a l s t a t e sf o rn u m e r i c a ls i m u l a t i o no fm u l t i f l u i df l o w s t h es c h e m ed e s i g n e db yt h i sa p p r o a c h i sf u l l yc o n s e r v a t i v ei nt h es e n s et h a ti tc o n s e r v e st h et o t a la sw e l la si n d i v i d u a lm a s s e s 。 m o m e n t u ma n dt o t a le n e r g y , i tm a i n t a i n st h ep o s i t i v i t yo ft h et o t a la n di n d i v i d u a ld e n s i t i e s ， a n dm o r ei m p o r t a n ti tm a i n t a i n st h eu n i f o r m n e s so ft h ev e l o c i t ya n dp r e s s u r e a c r o s sm a t e r i a li n t e r f a c e n u m e r i c a le x a m p l e sa l ed i s p l a y e dt os h o wt h ee f f i e n c yo ft h e s c h e m ei ns i m u l a t i n gt w o - c o m p o n e n tf l o w k e yw o r d s ：m u l t i f l u i df l o w ，s p u r i o u so s c i l l a t i o n ，n o n e q u i l i b r i u m 原创性声明 本人声明：所呈交的论文是本人在导师指导下进行的研究工作除了文中特别 加以标注和致谢的地方外，论文中不包含其他人已发表和撰写过的研究成果参与 同一工作的其他同志对本研究所做的任何贡献均已在论文中作了明确的说明并表示 了谢意 本论文使用授权说明 本人完全了解上海大学有关保留、使用学位论文的规定，即：学校有权保留论 文及送交论文复印件，允许论文被查阅和借阅，学校可以公布论文的全部或部分内 容 ( 保密的论文在解密后应遵守此规定) 签名：庭乒冠因导师签名：趁日期： 第一章绪论 数值模拟可压缩多介质流体流动的复杂现象是近年来计算流体力学中的一大研 究热点在这一章中，我们将主要介绍多介质流的问题描述，困难所在以及这一领 域的发展现状 1 1 多介质流体动力学问题的介绍以及它的发展现状 可压缩多介质流体计算的主要难点在于数值模拟可能造成速度和压强在介质界 面附近的伪振荡因为即使是最好的具有一阶精度的格式，例如l a x f r i e d r i c h s 格式 或者是g o d n o v 格式，在计算多介质流体介质界面时已经产生了伪振荡，而从这些 格式扩展而来的高阶的格式则更难消除这些振荡关于这一点在文献1 ，【2 】， 【3 】， 8 】，f 9 】9 ， 1 0 中有较详细的论述在接触面以外的位置，流体实际上是单流 体，对于单介质流体已有很多有效的数值计算方法而如何处理多介质流介质边界 的间断问题，则成了主要的研究对象处理间断的数值方法大致可以分为两类：跟 踪法和捕捉法跟踪法可以精确确定介质界面的位置，但是我们很难把它推广到高 维空间中去；捕捉法可以直接推广到高维空间上去，但是利用这种方法得到的介质 界面不如跟踪法精确它们各有利弊 1 2 问题的数学描述 首先我们给出本文所研究问题的数学描述一般的无粘可压缩流体运动都可以 用e u l e r 方程组来描述，它由流体的三个守恒律所组成，即质量守恒，动量守恒以 及能量守恒具体方程如下： 肌+ ( p t ) z = 0 ， ( 肚) t + ( p u 2 + p ) 。= 0 ， ( 1 2 1 ) 西+ ( u ( e + p ) ) 。= 0 ， 其中p 是密度，u 是速度，p 是压强，e 是总能量，e = ! p u 2 + e ，这里e 是内 能如果再加上一个状态方程，那么整个方程组就是封闭的，并且在初始条件和边 界条件给定的情况下可以求得唯一解 对于单介质理想流体，我们有它的状态方程： p = ( 7 一a ) e ( 1 2 2 ) 对于单介质s t i f f 气体，我们有它的状态方程： p = ( 7 1 ) e 一7 b ( 1 2 3 ) 1 2 0 0 8 上海大学硕士学位论文 2 对于单介质b a r o t r o p i c 气体，我们有它的状态方程： p ：( 7 1 ) ( e + 譬) 一7 b ( 1 2 4 ) 对于单介质v a n - d e r - w a a l s 气体，我们有它的状态方程： p = ( ，y - 1 ) 竿鲁一n 矿 ( 1 2 5 ) 其中，y ，b ，a ，b ，p o 为一些状态常量前三个状态方程可以写成一个一般的 形式： p = n 一1 ) ( e + p a ) 一 t b ( 1 2 6 ) 在本文中我们将考虑仅含上述两种状态方程相互作用的流体，而模型的多种流体介 质可以简单地予以推广 如果流体是由多种介质构成的，并且在开始状态各种介质是互相分开的，没有 发生混合那么由于它们没有粘性以及互相绝热的特性，在以后的时间里它们将仍 然相互分开，不发生混合在不同的介质处我们将有不同的状态方程，例如所有的 理想气体都有( 1 2 2 ) 形式的状态方程，唯一的区别是不同的理想气体会有不同的热 容比，y 当然所有的流体介质均在介质界面处与别的介质相连接，而事实上这个连 接处也正是大多数应用问题主要关注的 然而在对多介质流的大多数有限体积法的数值模拟中，迄今为止所有的数值格 式都将引入一定量的数值耗散，不同介质的流体在介质边界处会发生一定量的数值 混合因此人们提出了一些改进的模型来代替方程组( 1 2 1 ) 来描述多介质流体，以 应付必然要发生的数值耗散至今为止所有的改进模型都是基于混合区域达到热平 衡这一假设之上的由于在单介质流体的时候，所有区域都处于热平衡状态，而且 所有的热力学量如压强、温度、内能、熵和焓等只有在热平衡的时候才有定义，因 此这么做似乎无可厚非于是基于这个假设，混合区域的状态方程也能够被推导出 来例如如果两种介质均为理想气体，那么混合气体也是理想气体，因而有它自己 的热容比，y ，它是组成它的单种介质的热容比的某种意义上的加权平均： 7 = 笔笔 7 2 特1 群p 2 ( 7 1 1 ( 1 2 7 ) p l (一) +一) ” 其中，y 1 ，仇，p 1 ，比分别是两种介质的热容比和密度，第二章中我们将对它作具体推 导 一个常用的描述两种介质流体的改进模型叫做质量分配模型( 见【l 】) ，它的具 体形式如下： ( p 1 ) t + ( p l u ) 霉= 0 ， ( p 2 ) t + ( 化冀) z = o ，( 1 2 8 ) ( 肚) t + ( p u 2 + p ) 2 = 0 ， 、7 毋+ ( u ( e + p ) ) z = 0 2 0 0 8 上海大学硕士学位论文3 其中p 1 ，, 0 2 是两种介质各自的密度，它们满足p = p l + 化，u 是速度，p 是压强， e = 肚2 + e 是总能量，e 是内能当然还有其他很多改进模型，在【1 】中都有具体 的介绍当气体没有粘性，不发生热传导的时候这些模型都相互等价但是如果 有粘性或者热传导的时候，即使很小，它们将不再相互等价 一个好的数值模拟应该是守恒的，对于总质量以及各种介质的分质量、总动量 以及总能量均守恒但是所有那些对于单介质流体的守恒型格式推广到多介质流体 的时候都产生了伪振荡这些振荡首先出现在混合区域流体的压强上，然后再影响 到该区域内所有的其它物理量在第二章中我们以理想气体为例说明这些振荡的确 是由于数值耗散所引起，特别是其中两种不同气体介质之间的伪热传导所造成的 但是由e u l e r 方程组( 1 2 1 ) 所描述的流体是无粘且绝热的，因此不会有任何混合 发生，在介质界面处也就不应该有任何振荡即使是由n a v i e r - s t o k e s 方程组所描述 的实际气体，虽然不同气体介质间存在粘性和热传导，但是混合区域的宽度十分狭 小，当粘性和热传导很弱的时候压强在这个区域的振荡是不可见的问题是对于数 值格式，粘性和热传导都被放大到了网格宽度的程度，这远远大于实际，因此流体 压强上的振荡就变得可见了 为了解决这个问题，很多物理现象的实质在数值模拟的时候不得不被牺牲掉 近年来，个普遍的做法是牺牲解的守恒性，很多非守恒以及半守恒的方法被发展 出来，我们可以在【l 】 f 2 】 8 】，【9 】，【1 0 】等参考文献中找到这些方法但是在很多问题 中，解的守恒性其实是不可取代的，它直接关系到间断位置是否正确所以对于那 些对间断位置十分敏感的实际应用，现有的方法显然难以满足 王长春和茅德康 2 9 】针对多介质流体中的两种理想气体设计了一种基于非平衡 态的l a x f r i e d r i c h s 格式，该格式保证了多个守恒律：总质量和分质量守恒，总动量 和总能量守恒它还可以保证分质量非负更重要的是它消除了速度和压强在介质 界面的伪振荡在对不同理想气体的数值模拟中取得了良好的计算效果 本文是上述工作的继续我们所做的主要工作就是将这种基于非平衡态设计格 式的思想推广到更为一般的流体中去我们提出了一种计算多流体的l a x f r e i d r i c h s 格式，它是针对一般流体的而言的我们对所设计的格式讨论了四种不同流体两两 相互作用的情况： ( 1 ) p o l y t t o p i c 气体；( 2 ) b a r o t r o p i c 气体；( 3 ) s t i f f 气体；( 4 ) v a nd e rw a a l s 气体，并作了数值实验在设计中我们对流体作了几个必要的假设 为了跟踪介质界面的运动，我们引入了跟踪介质界面的运动方程： k + u 圪= 0 ， ( 1 2 9 ) 其中y 为混合区域中某一种介质所占的质量份额我们利用我们得出的格式在对 上述几种流体相互作用的数值模拟中取得了良好的计算效果 2 0 d 8 上海大学硕士学位论文4 本文的结构如下：第一章为绪论，第二章为基于平衡态的l a x - f r i e d r i c h s 格式， 第三章为基于非平衡态的物理模型及l a x - f r i e d r i c h s 格式的建立第四章为基于非平 衡态的l a x - f r i e d r i c h s 格式的理论分析，第五章为数值算例，第六章为结论和展望 第二章基于平衡态的l a x - f r i e d r i c h s 格式 在这一章我们将介绍一些相关的热力学背景，推导出基于平衡态的l a x - f r i e d r i c h s 格式并以理想气体为例推导出速度和压强产生振荡的原因 2 1 热力学知识背景 为了更深入地理解混合流体问题，有必要先回顾一些相关的热力学知识在热 力学中，我们的研究对象为闭合的热力学系统，它可以具备任意量的物质，并且性 质能唯一且完全地被一定的宏观量所描述这些物质被物理墙壁与周围分开，系统 可以与周围交换能量，但不能和周围交换质量热力学中最重要也是最基本的概念 是热力学平衡，它是一个系统经过足够长的时间后自动到达的一种状态，当系统达 到热平衡后，它的热力学特征，比如温度、压强、熵等描述系统状态的参数，将不 再随着时间而改变，同时整个系统的热力学状态量也只有在平衡状态下才有定义 相对的，个封闭系统在达到热力学平衡前所有的状态都不是平衡态，它的热力学 状态量无法定义，因为系统将随着时间不断改变从热力学平衡这个概念伴随而来 的个重要结论是：如果两个系统相互之间处于热平衡状态，那么他们将有相同的 温度，这也是为什么温度计可以测量温度的原因一个系统若是处于热力学平衡状 态，那么它的任意子系统都将处于热力学平衡状态而对于一个处于非平衡状态的 系统，它有可能可以被分成若干个处于热平衡状态的子系统例如一个温度计和被 测量组成的系统，在温度计停止变动前系统处于非平衡态，但是每一时刻温度计和 被测量物体本身构成的子系统都处于热平衡态，他们有自己的温度从上面这个例 子我们也可以发现，个系统达到热平衡的主要手段是自身以及和周围环境进行热 交换更多的概念可以在热力学教科书中找到 2 2混合理想气体的热容比 处于平衡态的理想气体必须有相同的温度假设我们有一由两种理想气体构成 的混合气体，分别用下标1 和2 表示，并且有不同的等温比热和等压比热g ，g ，以 及不同的热容比7 两种气体的密度分别为p 1 ，p 2 假设混合气体处于热平衡并且 也是一种理想气体，我们计算混合气体新的热容比7 首先，我们知道，当混合气体处于热平衡的时候，混合气体中的每一种介质作 为个子系统也处于热乎衡状态并且它们热力学性质由它们自身的状态方程所控 制因此，各种介质的分压强可以由理想气体状态方程得到： 当= e l 兰：e 2 ( 2 2 1 2 ) 百了2e 1 孑i 2 眈( j j 5 2 0 0 8 上海大学硕士学位论文6 并且两种介质相同的温度可以分别用下面两个公式计算， 而6 1 = 正磊e 2 = t ( 2 2 2 ) 石而纠劢i 2 卫) 由道尔顿( d a l t o n ) 定律，我们得到 p = p l + p 2 ( 2 2 3 ) 同时由能量守恒定律，我们得到 e = e l + e 2 ( 2 2 4 ) 注意到混合气体处于热平衡状态并且被看作是一种理想气体，因此( 1 2 2 ) 式对它同 样成立将( 2 2 3 ) ，( 2 2 4 ) 代入其中，我们得到 雩警：e l + e 2 ( 2 2 5 ) 1 1 。 、。7 然后将( 2 2 1 ) ，( 2 2 2 ) 带入( 2 2 3 ) ，可得到 p l + p 2 = ( 7 l 一1 ) p l g ，1 t + ( 1 ，2 1 ) m ( 乃，2 t ( 2 2 6 ) 由( 2 2 4 ) 可得到 6 1 + e 2 = p l g ，i t + m c ，2 t ( 2 2 7 ) 最后把( 2 2 6 ) ，( 2 2 7 ) 带入( 2 2 5 ) 得到 7 一- = 虹絮器甍铲 仫2 固 由热力学的相关知识对理想气体我们可以得到 q ，1 一g ，1 = r ， c p ，2 一c v ，2 = r ， 1 1 ：孕， ( 2 2 9 ) 他：罄， 其中r 为气体常数经过整理我们可以得到 c v 户击，= 击( 2 2 1 0 ) 将( 2 2 1 0 ) 代入( 2 2 8 ) ，经过整理我们可以最终可以得到 7=型芝苷黼(2211p1(72 1 m ( 7 1 1 ) 一) +一) 7 这一混合气体的热容比的形式被很多计算流体力学的相关研究论文所引用，我 们在这里介绍它的推导只是为了保证我们叙述的完整性接下来我们将要看到，正 是由于使用这一热容比，y 的形式，最终导致数值方法中在介质接触面处压强的不一 致 2 0 0 8 上海大学硕士学位论文7 2 3 格式的导出 l a x - f r i e d r i c h s 格式是一个经典的守恒形格式，如 1 1 】和【1 2 】中所描述的那样， l a x f r i e d r i c h s 格式可以被看成是g o d u n o v 格式的一种变换两种方法都通过将数 值解作为网格平均的逼近，然后解各个网格交界处的r i e m a n n 问题，并在下一时 刻对这些解取网格平均作为这一层上的数值解唯一的区别在于l a x f r i e d r i c h s 格 式定义在偏移的网格上，而g o d u n o v 格式定义在普通的正交网格上精确地说， l a x - f r i e d r i c h s 格式在一个时间层上定义在以z f 为中一5 - 的网格上，然后下一个时间 层上定义在以z 件l 2 为中心的网格上，其中x j = j h ，x j + 1 2 = ( 歹+ 1 2 ) h 下面我们 将直接写出方程组( 1 2 8 ) 所对应的l a x - f r i e d r i c h s 格式的离散形式，这里气体在混 合区域被认为是热平衡的，并且使用理想气体的状态方程，热容比，y 由混合理想气 体热容比所决定我们还将分析通过介质边界的时候压强的一致性是怎样被数值耗 散，特别是两种流体成份间的热传导所破坏的 基于平衡态的l a x f r i e d r i c h s 格式的离散形式如下： ( p 1 ) 麓2 ( ( p 1 ) + ( p 1 ) 苒1 ) 一害( ( p 1 ) 苒1 略一( p i ) 7 u ? ) ， 象 、穗n + l ：三l ( 。( 露p , 2 札) n ；+ 蠢。n i 。一a 害。废n ：赢u nj 。+ 毒n 。? n 二。叼，2 一刃，c 2 s ，( 肚) 麓= ( 露札 + 癣l 略1 ) 一宝( 饵1 ( u 髯1 ) 2 + 略1 一力( 叼) 2 一刃) ， 、 f 麓2 烈1q n + e h l ) 一含( u 苒1 曙1 一u 曰) 一砉( u 孙1 p ”i + l 一钆 n n n ) 在( 2 3 1 ) 中佗+ 1 层上的压强p 凄是按下面的状态方程来计算的： 碟= ( 域1 ) ( 嘣一知“s n + l 讹) 划n + l 缸蔫) 2 ) ( 2 3 2 ) 其中飞) i n + + l 是混合流体的热容比，它由( 2 2 1 1 ) 式来计算，也即在该式中以( p 1 ) 葛 和( p 2 ) 蔫来替代p l 和化 众所周知，l a x f r i e d r i c h s 格式相较于g o d u n o v 格式有一个很大的缺点，就是 有太多数值耗散但它为格式的理论分析提供了足够的信息，而且高阶的格式也可 在此基础上简单地发展出来，见【2 8 】 2 4 基于平衡态的l a x - f r i e d r i c h s 产生振荡的理论分析 接下来我们将分析在个纯粹的介质交界面压强的一致性是如何被基于平衡态 的l a x f r i e d r i c h s 格式所破坏的换句话说，就是它如何被计算格式产生的数值耗 散特别是两种介质间的热传导所破坏的事实上，同样的论述在很多与多介质流体 计算相关的文章中都有所提到，在本文的绪论中我们列举了相关论文但它们绝大 2 0 0 8 上海大学硕士学位论文 8 多数都是以g o d u n o v 格式的形式出现，这里我们将用l a x - f r i e d r i c h s 的代数形式描 述不管怎样，他们的核心思想是一致的 p i 0 m i e 1 0 ，p r m r e r x j 。x j o + 1 2 t i m e - - - t n x j o + l 2x j o + 1 t i m e 2 t n + 1 囝2 4 1 一个介质边界以速度雹 0 从左向右移动，相应地，在它两面流体的状态为渤，0 ，m l ，e 1 ) ta n d ( o ，p r ，m ，b ) t ，其中m l = p j 豇a n dm ，= p r 豇分别是两边流体的动量边界在t 。时位于x j o + 1 2 并且在一 个时间间隔内向右移动a = ；忱 我们考虑一个纯粹的介质界面问题，并且流体在介质边界处的传输速度为u 0 如图( 2 4 1 ) 所示，一开始它位于网格边界+ l 2 处，在接触面两侧，流体是两 个常状态， 叼= 豫p l , m o , p l u 舭, e 蚓1 ) t , r ，雾 j 如o ( 2 4 1 ) 其中 局= 三( 胁) u 2 - k - 告，e r = 圭( 肼) “去， ( 2 4 2 ) p 为接触面两边共同的压强，而7 1 和 7 2 是两种气体各自的热容比 我们以带上划杠的符号来表示下一时刻在j o + 1 2 网格中的数值因为接触面 以速度u 运动，我们由格式( 2 3 1 ) 可得到 卢1 = 言( 1 + 入乱) p z ，庇= 去( 1 一入u ) j d ， ( 2 4 3 ) 2 0 0 8 上海大学硕士学位论文 9 以及 历= 言( 1 + a u ) p l u + 言( 1 一a t i ) 肼u ( 2 4 4 ) 注意到卢= 卢1 + 庇，从( 2 4 3 ) 和( 2 4 4 ) 可得到，当气体在偏移网格内充分混合后有 面：一p u ：u ， ( 2 4 5 ) p 这意味着速度在计算中保持一致 然而，压强却无法保持这性质下面由能量入手分析流体的压强由( 2 3 1 ) 中最后一式得到 e = 丢( 1 + 入u ) 局+ 去( 1 一入u ) e ， ( 2 4 6 ) 注意到 豆= 去芦面2 + 百，芦= 芦l + 历，面= 让，( 2 4 7 ) 我们可以得到 百= 三( 等+ 岩协 ( 2 4 8 ) e2 互( 了可+ 彳i ) p 苎j 现在在这个偏斜的网格内混合气体在t n + l 时刻达到热平衡，它拥有( 1 2 2 ) 式的状 态方程，并且彳可以像( 2 2 1 1 ) 中那样计算最终我们得到 乒一1丛脊滁高等凳黔掣c岩+d害u12 111111m 1j=一一l一 i t ， ( + a u ) 见( 一y 2 一) + (一a u ) p ，( 一y 1 一)、，y 1 一7 2 7 ( 2 4 9 ) 在上式中，我们可以看到，p ，f 几乎不可能相等，除非是在单介质流的时候- y 1 = 他， 这时压强才能保持一致 我们从上面的讨论中可以看到在通过介质边界时压强的一致性被混合气体在 混合区域达到热平衡这一假设所破坏这也最终导致了数值解在介质边界附近的剧 烈振荡因为混合气体的热容比7 与两种介质的热容比1 1 ，1 2 均不相同，所以由 ( 2 4 9 ) 所计算的压强在这个混合区域也和两边的压强均不相同由于( 2 4 9 ) 中压强 的振荡在单介质流体的情况时并不存在，我们可以断定振荡的产生正是由于两种介 质间的热传导所造成的 在以e u l e r 方程组所刻划的流体内，流体被认为是无粘性且绝热的，因此在介 质边界附近并不发生混合即使在由n a v i e r - s t o k e s 方程所刻划的实际气体中，虽然 由于粘性和热传导，存在一个气体混合并局部达到热平衡的区域，但在r e y n o l d s 数 很大，热传导很弱或气体几乎绝热的情形，这个混合区域将非常之小因此在这些 流体中压强的振荡应该是肉眼所不可见的然而在基于平衡态的l a x - f r i e d r i c h 格式 中，在每个小网格内气体充分混合且达到热力学平衡，这使得气体的粘性和热传导 都被扩大到网格长度a x 的程度，从而使得压强的振荡变得可见 第三章基于非平衡态的物理模型及l a x - f r i e d r i c h s 格式的建 立 这一章是本文的核心部分在这一章中将介绍我们所设计的用于计算多介质流 体基于非平衡态的l a x - f r i e d r i c h s 格式就像我们在绪论中所说的那样，我们将放弃 混合气体处于热平衡这一假设，而不是牺牲解的守恒性我们设计的格式在对不同 流体的计算过程中得到良好的计算效果 3 1 基于非平衡态方法的前提假设 在前一章中，我们看到压强的振荡是由于数值耗散所引起的，最主要的就是其 中两种介质间的数值热传导为了保证格式的稳定性我们不能完全消除格式的数值 粘性同样的原因，我们也不能取消单种介质内的数值热传导但是我们能够通过 使两种介质之间绝热来阻止两种介质问的数值热传导现在，两种气体介质在热力 学意义上被分割开来，从而他们之间没有任何热交换这样在有混合发生的区域内 的气体不再处于平衡态，为此我们作出如下的假设： ( 1 ) 两种流体在空间和运动学意义上充分混合，也即流体各处的总密度均等于 两种流体的分密度之和，两种流体的速度一致 ( 2 ) 每种单一的流体自身达到热平衡，然而两种介质之间绝热，从而他们之间 没有统一的温度，因此总流体处于非热平衡态 ( 3 ) 混合气体的压力按分压方式计算，也即总压力等于各流体产生的分压力之 和 ( 4 ) 各流体的分压力是这样来计算的：将总流体视为同一种流体，计算其压力， 则单流体的分压力即为所计算的压力乘上了该单流体在总流体中所占的密度份额 ( 5 ) 压力对流体做功所传递的能量按流体的密度份额分配到各流体上去 我们对上述假设做如下四点解释： 一，在混合区域中，流体在空间意义上充分混合，即流体各处的密度等于两种 流体的分密度之和我们可以设想在一个封闭的容器中有两种不同流体，经过充分 长的时间后两种流体都均匀充满整个容器那么在容器中流体的质量应为两种流体 的质量之和利用密度和质量的关系： p = 号，( 3 1 1 ) ， 我们可以得到 =百m：ml矿+m2：p1+p2p ( 3 1 2 )2 百2 i 21 十l 5 l z j 】0 2 0 0 8 上海大学硕士学位论文 1 1 其中m ，m ，m 2 分别为装置中混合流体的总质量，第一种流体的质量和第二种流体 的质量它们满足关系式； m = m 1 + m 2 ， ( 3 1 3 ) 即混合流体各处的密度等于两种流体的分密度之和 两种流体充分混合后，它们的速度相等也是很好理解的因为这里的速度是表 征流体运动的宏观速度，当经过充分长的时间后，流体在宏观上达到了一个运动学 上的平衡态，其平均速度为一恒定量这时候由于流体充分混合，我们测得的速度 应为两种流体共同的速度 二，在混合区域中，我们对流体做了一个非物理的假设，我们假设两种流体相 互间绝热，因此即使经过充分长的时间，混合流体也不处于热平衡态，从宏观上讲 也就是说两种流体不具有相同的温度之所以会作出这个非物理的假设是基于我们 上一章对平衡态的分析，在前面我们看到正是由于流体在混合区域达到热平衡态导 致了压强产生伪振荡由于流体在混合区域中没有相同的温度，因而两种流体的能 量必须分开计算为此，我们必须引入第一种流体的总能量e 1 和第二种流体的总 能量场由能量的守恒和转化定律我们可以得到： e = 马+ 场，( 3 1 4 ) 其中e 为混合流体的总能量 三，为闭合方程组，我们必须要计算混合区域的总压强由于此时混合流体处 于非平衡态，因此我们压强的计算不可能像平衡态一样有一个统一的计算公式混 合区域中总压强的计算是根据分压方式进行的：即混合区域中的总压强等于两种流 体产生的分压强之和然而，在混合流体之中各流体的分压强是很难计算的，为此 我们有了前提假设( 3 ) 和( 4 ) ，它们事实上是和如下的两个假设等价的 ( a ) 在混合流体之中任何一种流体都没有意识到另一流体的存在，任何一种流 体的热力学行为就如同它处于同类流体的包围之中一样这显然是一种g h o s t 的思 想，见【1 1 】和【1 2 】，然而这儿我们将之用于物理的建模之中，而不是在算法的设计 之中 ( b ) 在单一流体之中，部分流体的分压强对总压强的比例等于该部分流体的密 度对总密度的比显然，对理想气体这一结论是正确的，见【1 9 ，对一般的流体我 们沿用了这一结论 四，最后我们解释前提假设( 5 ) ，即由于压力做功而传递的能量的分配问题 注意，对于单流体运动并不存在这一问题，在单流体中一部分流体对另一部分流体 作功所传递的能量为另一部分流体所完全吸收问题存在于多流体的运动之中 2 0 0 8 上海大学硕士学位论文 1 2 首先，在? 昆合流体之中一种流体作功所传递的能量并不只传递给同种流体例 如在一管道之中，如果用一薄膜将两侧的流体隔开，且薄膜随流体运动此时如果 左面的流体对右面的流体作功，右面的流体得到能量，那么右面的流体是无法区分 哪一部分能量是第一种流体所传递的，哪一部分能量是第二部分流体传递的 其次，如前所述，在混合区域中两种流体在空间上是充分混合的，因此当一部 分? 昆合流体因其它流体对它作功而得到能量时，这些能量应该按照各流体所占的质 量份额分配给各流体 下面我们将根据这些前提假设推导出我们基于非平衡态的l a x - f r i e d r i c h s 格式 3 2 格式的导出 按照上述假设我们得到了如下的六个方程的流体力学方程组： ( p 1 ) t + ( p l u ) 茁= 0 ， ( 化) t + ( p 2 u ) 2 = 0 ， ( 叫) t + ( 旷+ p ) z = o ， ( 3 2 1 1 ( 目) t + ( u e l ) z + y ( u p ) 。= 0 ， 、。 ( 岛) t + ( 心岛) 。+ ( 1 一y ) ( 印) 。= 0 ， k + 让圪= 0 ， 其中p l ，p 2 ，e 1 ，e 2 分别表示混合流体中第一种流体和第二种流体的密度和能 量，“和p 分别表示流体的速度和总压强，y 为第一种流体的密度份额，而1 一y 为第二种流体的密度份额这儿前三个方程是根据假设( 1 ) 得出来的；而第四个和 第五个方程是根据假设( 5 ) 得出来的最后一个方程摸拟的是介质界面的运动 为闭合方程组，我们还需要个计算总压力的状态方程由假设( 2 ) ，( 3 ) 和( 4 ) ， 压力的计算原则上是按下述方式来进行的， p = p 1 + p 2 = y p l ( p ，e ) + ( 1 一y ) 岛( p ，e ) ，( 3 2 2 ) 其中p = 只( 肛e ) 和p = p 2 ( p ，e ) 是两种流体的各自的状态方程，而e 是内能密度， e = l z p u 2 + e ( 3 2 3 ) 状态方程( 3 2 2 ) 右端出现的p ，e 应是总质量密度和总内能密度，然而如果在( 3 2 2 ) 右端的计算中出现y p ，( 1 一y ) p ，y e ，( 1 一y ) e ，我们分别代之以分质量密度 和内能p l ，p 2 ，e l 和e 2 例如，我们将在本文中考虑p o l y t r o p i c 气体( 见f 1 】) ， b a r o t r o p i c 气体( 见【1 0 ) ，s t i f f 气体( 见【1 】) 和v a i ld e rw a a l s 气体( 见【9 】) 根据前 2 0 0 8 上海大学硕士学位论文1 3 面讨论知，前三种流体的状态方程可统一写成( 1 2 6 ) 此时( 3 2 2 ) 式中的第一种流 体的分压力p 1 为 fh 一1 ) ( e l + a 1 p 1 ) 一y 7 i b l ，状态方程( 1 2 6 ) ， p 12y p l p e ) 21 ( 7 l 一1 ) 鱼掣一n l p l p ，状态方程 1 2 5 孓 ( 3 七4 ) 同理可得( 3 2 2 ) 式中的第二种流体的分压力p 2 的计算 i ( 7 1 1 ) ( e 2 + a 2 p 2 ) 一( 1 一y ) y 2 8 2 ，状态方程( 1 2 6 ) ， 沈2 ( 1 一y ) 岛 e ) 。1m 1 ) 鱼掣口2 化p 状态方程；1 2 5 i 、 一u 2 p ( 3 2 5 ) 注3 1 方程组( 3 2 1 ) 中的密度份额y 是对p l 化和e l e 的逼近，在数值计算 中它将独立计算，因此在数值上y 和p l p 并不一致我们之所以需要一个独立的 y 是为了在计算纯密度流时能保证速度和压力不变，见下文定理4 2 如果两种流 体均是p o l y t r o p i c 气体，则l ，是不需要的，见f 2 9 1 注3 2 如此建立的物理模型是非现实的，它是对真实物理现象的近似需要指 出，我们所建立的模型和【3 】中的模型本质上是完全不同的 注3 3 对混合流体压强的计算，我们是以前面提到的四种流体为例来进行说明 的但并不是说我们设计的格式只对上述四种流体成立实际上我们设计的格式对 一般的流体方程也是成立的只是这四种流体方程比较典型，别人也计算过它们两 两相互作用的情况，所以我们以这四种流体方程为例来介绍我们的算法 3 3 基于非平衡态的l a x f r i e d r i c h s 格式的具体离散形式 籀毒燕一- u i p 岬i ) , a k + 2 主( m + m + 1 ) 一三! 生之竺型( k + 1 一m ) ， 2 0 0 8 上海大学硕士学位论文 1 4 厩+ 2y i + p l ( z i + ，毛+ ) + ( 1 一k + 言) 马( 反+ ；，咆+ ) ， ( 3 3 2 ) 其中p l ，b 状态方程如( 3 2 4 ) 和( 3 2 5 ) 中定义， 反+ 吾和毛十吾分别为第n + l 层上 的总密度和总内能 须要注意的两点是： 一，在我们的计算格式中气体流体介质互相绝热，因而我们也将不用涉及到取 两种气体介质热容比饥，忱某中间值的混合气体热容比，y 在下一节的理论分析中 我们将看到正是由于这一特性使得我们的格式在介质接触面附近不会产生任何振 荡 二，( 3 3 1 ) 中最后两个方程并不是以守恒形式出现的，但是这并不说明我们的 格式是非守恒的，将它们相加我们会得到总能量仍旧是守恒的，而事实上在上一章 所描述的最初始的l a x f r i e d r i c h s 格式中两种气体介质的各自能量也是不守恒的， 它们通过作功和热交换的途径交换能量，而我们所做的就是把能量交换限制在做功 上，两种流体之间没有热交换 第四章基于非平衡态的l a x - f r i e d r i c h s 格式的理论分析 关于格式( 3 3 1 ) 我们有如下的理论结果： 定理4 1 当c f l 条件满足时，所设计的格式( 3 3 1 ) 满足： ( 1 ) 总质量及分质量守恒，总动量和总能量守恒， ( 2 ) 分质量非负， ( 3 ) 质量密度份额y 在【0 ，1 】之间取值 ( 4 ) 当y = 0 ( y = 1 ) 时，或者当p 1 ( p ，e ) = b ( p ，e ) 时，格式( 3 3 1 ) 退化为单流 体的l a x - f r i e d r i c h s 格式 证明：( 1 ) 由格式( 3 3 1 ) 第一式和第二式相加，我们可以得到总质量的计算公 式 履+ 丢= 吉( 风+ p i + 1 ) 一号( 店+ 1 让件l p i u i ) ( 4 1 ) 将( 3 3 1 ) 的第四式和第五式相加，我们可以得到总能量的计算公式 毫+ = 三( 晟+ 晟+ ) 一害( 日+ l u i + t 一厩地) 一害( u 件p i + 一“t p i ) ( 4 2 ) 格式( 3 3 1 ) 的第三式即为总动量的计算公式由于对总质量以及分质量，总动 量和总能量的计算格式都是守恒型的，因此数值解将仍旧满足守恒律 ( 2 ) 由( 3 3 1 ) 的第一式我们不难得到 卢。，；+ ；= 丢( 1 + a u t ) p - ，t + 壶( 1 一a “t + - ) p ，t + - ( 4 3 ) 当c f l 条件成立时 a m a x i u i i 1 ， ( 4 4 ) 因此( 4 3 ) 中的两系数( 1 + a u i ) 和( 1 一a u 件1 ) 均大于0 ，故所得的关于第一种流体 的分质量非负 由( 3 3 1 ) 的第二式不难得到 忍，件2 考( 1 + a t ) 化，i + 主( 1 一a u 件1 ) 比一l ( 4 5 ) 11 由( 4 4 ) 我们可以得到( 4 5 ) 中的两系数( 1 + a u i ) 和( 1 一a u i + 1 ) 均大于0 ，故所得的 关于第二种流体的分质量非负 ( 3 ) 由( 3 3 1 ) 的最后一式可得 = 三( 1 + 垫掣) k + 互1 ( 1 一掣) ( 4 6 ) 1 5 2 0 0 8 上海大学硕士学位论文 1 6 同样当c f l 条件( 4 4 ) 满足时，( 4 6 ) 式中的两系数均非负，且它们的和为1 ， 从而我们可以得到质量密度份额y 在( 0 , i 】之间取值 ( 4 ) 将格式( 3 3 1 ) 中的第一，二两式相加，第四，五两式相加就得到了计算总密 度和总能量的公式如果再加上计算总动量的第三式，那它们与单流体计算的l a x - f r i e d r i c h s 格式形式上是一致的，只是格式( 3 3 1 ) 计算压力的方式与单介质流体不 一样我们计算的压力是通过分压方式得到的，而单流体的l a x - f r i e d r i c h s 是通过单 流的状态方程得到的然而由( 3 2 2 ) ( 3 2 6 ) 式的讨论不难得出，当y = 0 ( y = 1 ) 时，此时对应的流体即为单流体，压强的计算是按照单流体的状态方程进行的 当p 1 ( p ，e ) = p 2 ( p ，e ) = p ( p ，e ) 时，虽然此时对应的是流体是有混合发生的情 况，但两种流体具有相同的状态方程，由基于非平衡态压强的计算公式( 3 2 2 ) ，我 们得到混合流体的总压强p 满足： p = p l + p 2 = y 只( p ，e ) + ( 1 一y ) 最( 店e ) = p ( p ，e ) ( 4 7 ) 此时流体的压强也退化为单流体计算的l a x - f r i e d r i c h s 格式压强的计算公式 综上所述，当y = 0 ( y = 1 ) 时，或者当p 1 ( p ，e ) = 马( p ，e ) 时，格式( 3 3 1 ) 退化 为单流体的l a x - f r i e d r i c h s 格式 定理4 2 当两种流体的状态方程都具有( 1 2 6 ) 的形式时，那么对纯密度流，也 即速度和压强是常数的流，格式( 3 3 1 ) 可保持速度和压强不变 证明：根据条件可假设在第n 层上的数值解如下； w , o ： ( 儿t ，, 0 2 两p i u , e 蛳，e 2 ，t ) ， ( 4 8 ) 【( p l ，i + 1 ，成，i + 1 ，p i + l