（计算数学专业论文）守恒型跟踪法在切向间断处的算法改进.pdf

摘要 茅德康近几年发展了一种守恒型的间断跟踪法7 ，8 ，9 1 该跟踪法是以解的守 恒性质作为跟踪的机制，而不是像传统的跟踪法利用r a n k i n e - h u g o n i o t 条件来跟 踪间断。刘妍将该算法应用到一维e u l e r 方程组上，设计了一个完整而强健的算法， 编制了相应的程序，取得了很好的计算结果1 0 1 。本论文对刘妍所做算法做进一步 改进和完善，将切向间断两侧的质量分开来算，改进后的算法完全保证了切向间断 两侧流体的质量守恒，最后用数值算例来验证该算法的可行性。 关键词：守恒型方程组，间断跟踪法，切向间断，网格平均 i a b s t r a c t i nr e c e n ty e a r s ，m a od e v e l o p e dac o n s e r v a t i v ef r o n t t r a c k i n gm e t h o d ，w h i c h w a sb a s e do nt h ec o n s e r v a t i o np r o p e r t i e so ft h es o l u t i o n s ，r a t h e rt h a no i lt h eh u g o - n i o tc o n d i t i o n sa st h et r a d i t i o n a lf r o n t t r a c k i n gm e t h o d sw e r e l i ua p p l i e dt h i s m e t h o dt ot h ee u l e re q u a t i o n so ff l u i dd y n a m i c si no n es p a c ed i m e n s i o n ，a n dd e v e l - o p e dar o b u s ta n dg e n e r a l - p u r p o s e da l g o r i t h ma n di t sf o r t r a nc o d e ，a n do b t a i n e d g o o dn u m e r i c a lr e s u l t s 1 0 i nt h i sd i s s e r t a t i o n ，w ei m p r o v el i u sa l g o r i t h mo nt h e c o m p u t a t i o na tc o n t a c td i s c o n t i n u i t i e s t h ei m p r o v e da l g o r i t h mm a i n t a i n st h ec o n - s e r v a t i o no fm a s so nt h et w os i d e so ft r a c k e dc o n t a c td i s c o n t i n u i t i e s n u m e r i c a l e x a m p l e sa r ed i s p l a y e dt os h o wt h ee f f i c i e n c yo ft h em o d i f i e dm e t h o d k e yw o r d s ：c o n s e r v a t i v et r a c k i n gm e t h o d ，c o n t a c td i s c o n t i n u i t i e s ，d i s c o n t i n u i t i e sp o s i t i o n ，c e l l - a v e r a g e i i 原创性声明 本人声明：所呈交的论文是本人在导师指导下进行的研究工作除了文中特别加 以标注和致谢的地方外，论文中不包含其他人已发表和撰写过的研究成果参与同一 工作的其他同志对本研究所做的任何贡献均已在论文中作了明确的说明并表示了 谢意 期：2 1 选5 ：9 本论文使用授权说明 本人完全了解上海大学有关保留、使用学位论文的规定，即：学校有权保留论文 及送交论文复印件，允许论文被查阅和借阅，学校可以公布论文的全部或部分内容 ( 保密的论文在解密后应遵守此规定) 签名：j 蠡址导师签名：獭 第一章绪论及背景知识 1 1守恒型方程组和e u l e r 方程组 具有如下形式的一维空间的偏微分方程组 饥- t - 厂( 让) z = 0 ，( 1 1 1 ) 称为守恒型方程组，其中u = ( u 1 ，u 2 ，u m ) 是关于x 和t 的m 维的矢量函数， 称为守恒量，或状态变量，如流体动力学中的质量，动量和能量。更精确点就是u f 是第j 个状态变量的密度函数。2u j ( x ，t ) d x 表示的是该状态变量在区间 x l ，x 2 】 中在时刻t 的总量。我们称这些状态变量是守恒的指的是熙嘶( z ，t ) d x 关于t 是 不变的。，( ) = ( ( u ) ，如( u ) ，厶( 乱) ) 称为流函数，当m = 1 时，( 1 1 1 ) 即为单 个守恒律。该守恒型方程( 组) 是由物理定律在任意两点z ，和x 2 之间的如下积分 形式 罢上，u ( z ，t ) d x = m ( 亡) ) 一f ( u ( x 2 亡) ) ( 1 1 2 ) 得到的。表示在区间陋1 ，z 2 】中总的流体的量( 如质量，动量，能量等) 的变化仅仅与 两端点处的流有关，这就是守恒的基础，其中，( 乱( z l ，亡) ) 和厂( “( z 2 ，+ 亡) ) 分别表示在 x 1 和x 2 点的流入流出的量。 守恒型方程( 组) ( 1 1 1 ) 是拟线性双曲型方程( 组) ，对于它，不论所给的初值 及边值条件如何光滑，它的解都会产生间断，而且间断既可以产生( 自发生激波) ， 又可以消失。另外解还有连续部分，它可以是简单波，如压缩波，稀疏波等。间断在 物理和力学中就对应着激波和切向间断。 其中最著名的守恒型方程组是e u l e r 方程组，它是描述流体运动规律的基本方 程组，其在航空，航天，大气科学等领域有着很广泛的应用，下面我们简单介绍一 维e u l e r 方程组是怎样由物理定律得到的。假设在一长的管道中充满了气体，并且气 体的密度或者速度穿过管道的每一个横截面都是常数。令z 表示管道的距离，并且 p ( x ，t ) 是气体在x 点处t 时刻的密度。质量用以下方式来定义，即气体在任意部分 x ，和z 2 之间的总的质量由以下的密度的积分给定： t 时刻在 z 1 ，z 2 中气体的质量= p ( z ，t ) d t ( 1 1 3 ) ，x l 如果我们假设管壁是不可渗透并且质量既没有产生也没有减少，则气体的质量在某 一部分中的变化量仅仅是由流过端点x 1 和z 2 的气体的流量来决定的。 1 2 0 0 8 上海大学硕士学位论文 2 令v ( x ，t ) 是气体在点z 和t 时刻的速度，则气体流过这一点的速率或者流量 由下式决定： 在( z ，亡) 的质量流= p ( x ，t ) v ( x ，亡) ( 1 1 4 ) 气体的质量在【x l ，x 2 】中的变化率由z 1 处和x 2 处的流量的差来决定 爰z 叫) 如= 出“) 巾“) 一出2 巾2 ， ( 1 1 5 ) 这是质量守恒定律的积分形式。 假设p ( x ，t ) 和v ( x ，t ) 都是可微的函数，则可得 o 珈+ 跏叫m 州) ) 】如- 0 ( 1 1 6 ) 由于( 1 1 6 ) 对任意一部分x l ，x 2 】和任何时刻都成立，因此被( 1 1 6 ) 中的被积函 数必须等于零，即 p t + ( ) z = 0 ( 1 1 7 ) 这就是质量守恒定律的微分形式，也称为连续性方程。 按同样的方式可以推导出动量守恒方程 ( ) t + ( p v 2 + p ) z = 0 ，( 1 1 8 ) 及能量守恒方程 易+ ( v ( e + p ) ) 茁= 0 ，( 1 1 9 ) 这里p 是气体的压力，e = 互1 2 + p e 是总能量密度，其中丢2 是总能量中的动能 部分，而e 为内能密度。在上面方程的推导中我们考虑了压力作用对动量变化的影 响及压力做功对能量变化的影响。 把以上方程结合起来就得至 j e u l e r 方程组 p p v e + p v 2 + p v ( e + p ) = 0 ( 1 1 1 0 ) e u l e r 方程组是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 方程组。 2 0 0 8 上海大学硕士学位论文 3 在e u l e r 方程组中我们假设气体是化学和热力学平衡的，则内能是压力和密度 的已知函数 e = e ( v ，p ) ( 1 1 1 1 ) 这称为气体的状态方程，是涉及气体具体性质的热力学量之间的关系式，是由气体 本身的性质决定的。e u l e r 方程组必须与状态方程( 2 2 2 4 ) 结合起来才得到封闭的 方程组。对理想气体( 粒子之间的相互作用很微弱以致该作用可以忽略不计的气体) 来说，内能只是温度的函数，即e = e ( t ) ，而由理想气体定律知温度t 与p 和p 有 关，即p = r p t ，r 为常数。一个很好的近似形式是 e = c v l( 1 1 1 2 ) 其中c t ，是一常数，表示单位体积上的比热。此式表明，内能与温度成正比，这样的 气体称为多方气体，对多方气体来说，状态方程只和气体的比热比有关，即 7 = c v c v ， ( 1 1 1 3 ) 称为比热比，其中c p 为单位压力上的比热，这样的气体也称为“7 - l a w 气体”。多方 气体的状态方程是 p ：( ，y 一1 ) ( e 一昙2 ) ， ( 1 1 “) p = a ( s ) 矿， ( 1 1 1 5 ) 其中a ( s ) 是熵s 的函数。熵s 定义为 s = a v l o g ( p p 7 ) + 常数 ( 1 1 1 6 ) 易验证熵s 是u 的凸函数。则状态方程变为 = p ( p ) = 1 ，(1117pk p) = j2，l 1 ， 其中忌：k e s c t ，k 是一常数，当s 三常数时，则e u l e r 方程组就简化为只有两个方 程的方程组，称为等熵方程组 ”k p vp 卜 8 ， 此时状态方程( 1 1 1 7 ) 中的忌是常数。 2 0 0 8 上海大学硕士学位论文 4 1 2 有限体积法和守恒型格式 我们的算法将以有限体积法的差分格式作为基本格式，并且整个算法也是在有 限体积法的框架下描述的，因此需要介绍一下这一方法及其相关的一些基本概念和 理论。详细情况请参考f 6 1 0 有限体积法和有限差分法密切相关，但实际上有限体积法是在守恒律的基础上 得到的。下面简单介绍一下它的基本思想。 首先将空间区域划分成一些小的区间，称为网格，又称为有限体积。然后再 将时间区域进行剖分，即在( z ，t ) 平面上，做两组平行线，将整个( z ，t ) 平面划 分成长方形网格，其中网格线的交点称为节点。在空间方向上，第j 个网格记为 b 一1 2 ，巧+ 1 2 】，分别以巧一1 2 和+ i 2 为左右网格边界，网格中心为，网格长 度记为b = 巧+ 1 2 一巧一1 2 。在本文中我们采用的是正规网格划分，即网格长度就 是空间步长，h = = z ，则x j = l = l 2 = 4 - 宝= ( j4 - 1 2 ) h ，而时间步长7 _ = a t 是在计算中根据c f l 条件求得的。c f l 条件会在下面给出简单介绍。 对有限差分法来说，所求得的数值解让? 舻通常是对原方程( 组) 在节点 ( x j ，t n ) 处的精确解u ( x s ，t n ) 的逼近。而在有限体积法中，数值解叼是在时刻t n 对 精确解在网格( 有限体积) x s i 2 ，+ 1 2 】上的积分的平均( 网格平均) 晌逼近，即 n 型去 u ( z ，如) 出 ( 1 2 1 9 ) ”j q 一1 2 采用网格平均不但可以很容易地运用守恒型方程( 组) 的一些重要性质，并可以保 证得到的数值方法是守恒的。 由守恒律的积分形式( 1 1 2 ) 得： _ d x j + i 2u ( z ，亡) 如= m ( 舢卅弛( x j + i 2 , t ) ) ( 1 2 2 0 ) 将( 1 2 2 0 ) 从k 到t 竹+ l 积分可得 丢x j + l 27 比( x ，t 州) 如= 丢璧劣u ( z ，t n ) 如 一引口+ 1f ( u ( x j + i 2 ，t ) ) d t 一它+ 1 ，( u ( 巧_ 1 2 ，亡) ) 叫 ( 1 2 2 1 ) 一般来说我们不能精确的对( 1 2 2 1 ) 式的右端的时间积分进行赋值，因为u ( x i + 1 2 ，t ) 在每一网格边界上随着时间变化。但我们可以研究具有以下形式的数值方法 叼+ 1 = 叼一入( 兢1 2 一魂。：) ， ( 1 2 2 2 ) 2 0 0 8 上海大学硕士学位论文 5 此处入= 丁 是网格步长比，屠1 2 称为数值流函数，是在网格边界x j + i 2 处对精 确解的流函数f ( u ) 的网格边界平均的逼近 露+ 1 2 型 ( x j l - 1 - 1 2 , 让) ：；1 厂k + 1m ( x j l + l 2 , t ) ) d t ( 1 2 2 3 ) ，t n ( 1 2 2 2 ) 称为守恒型格式，因为它模拟了精确解的性质( 1 2 2 1 ) 。如果我们在所有网 格上求九哼+ 1 的加和，可得 矿1 = 呼一7 i ( 露+ 1 2 一乜。2 ) ( 1 2 2 24 ) 由此可知，除了在边界处的流外在整个区域上我们得到了精确的守恒型。 上述数值流函数是一个2 k 个变量的函数，即魂。2 = ，( 呼南+ l ，吼七) ，此 时有限体积格式是2 七+ 1 点的，并且与方程组( 1 1 1 ) 中的流函数，( u ) 是相容的， 即有： ，( 让，钍，u ) = ，( u ) ( 1 2 2 5 ) 再简单介绍一下c f l 条件：对任意的有限体积法或是有限差分法来说，c f l 条件是 保证算法稳定和收敛到偏微分方程的解的一个必要条件。若方程组( 1 1 1 ) 为非线性 的单个守恒律，贝j j c f l 条件是 一 u 2 瓦- 严l ，7 ( 咛) i 1 ， ( 1 2 2 6 ) 比率v 称为c f l 数或者更多的被称为c o u r a n t 数。若方程坌f i ( 1 1 1 ) 为双曲型方程组， 贝i j c f l 条件为 u = m a x 。m a xi ”( 钆) i 1 (1227)hu e u ”j = l ，2 ，3 、 其中，( u ) 是方程组( 1 1 1 ) 关于状态u 的第歹个特征值，并且矿n 是在k 时刻的 所有的网格匕的网格平均的集合。 1 3间断捕捉法和间断跟踪法 求解方程组( 1 1 1 ) 的数值方法可分为两大类：间断捕捉法和间断跟踪法。间断 捕捉法的特点是：在计算的过程中不考虑间断的存在而在整个流体的任何地方都用 几乎同样的数值格式，借助方法所固有的数值耗散性效应( 或者数值粘性效应) ，自 动捕获到所要计算的间断。它期望间断在数值解中表现为很窄的过渡层。该类方法 的思想比较简单，便于编写程序。古典的人工粘性法以及各种具有数值耗散性的有 2 0 0 8 上海大学硕士学位论文 6 限差分法都属于此类。近十几年来，人们研究设计了一系列求解方程组( 1 1 1 ) 的高 分辨率间断捕捉法，如t v d 格式，t v b 格式，e n o 格式，p p m 格式和带子网格技 巧的e n o 格式和w e n o 格式等，见 2 5 ，1 ，2 ，3 ，6 】。用这些差分格式计算得到的 数值结果精度比较高，基本上消除了在间断附近的非物理的震荡现象，同时对间断 的磨损程度也比较低。但是间断捕捉法不能分辨出间断在网格内的确切位置和解的 结构，并且在间断附近数值解的精度是0 ( 1 ) 阶的。 间断跟踪法是把间断作为移动的内边界来处理。这样整个求解区域就被间断线 分为若干区域，在每个区域中，解被假设为是光滑的，用计算光滑解有效的数值方 法来求解这些区域中的数值解，而间断的移动要满足r a n k i n e - h u g o n i o t 条件。这可 通过求解r i e m a n n 问题计算间断的移动速度，从而就可以确定间断在每一时间层上 的间断位置。但这样处理后，方程组( 1 1 - 1 ) 的守恒性一般不能保证。 较早的跟踪法思想见文f 2 7 】。在过去的三十多年中，人们研究设计了很多这种 类型的间断跟踪法f 1 5 ，1 6 ，1 7 ，1 8 ，1 9 ，2 0 1 ，我们称这类方法为传统的间断跟踪法。而 关于跟踪法如何在一维空间的实现可见f 2 9 ，2 2 ，2 6 ，2 4 ，4 1 。传统的间断跟踪法可对 解的光滑部分和间断位置同时获得高精度。但是算法很复杂，并且用该类方法求解 方程组( 1 1 1 ) 时会遇到很多问题。 近几年来，茅德康研究设计了一种基于解的守恒性质的间断跟踪法8 ，9 ，1 1 ，1 2 ， 1 3 ，1 0 1 ，它是利用守恒性来确定间断位置而不是像传统的跟踪法那样利用r a n k i n e - h u g o n i o t 间断条件。其基本思想来源于对间断捕捉法的如下观察： 1 对间断捕捉法，当解光滑时，差分格式与方程组( 1 1 1 ) 是相容的，但当解含有间 断和计算中产生很大误差时，守恒性能够帮助控制误差，从而使得即使相容性 不再成立时，也能够求得正确的间断位置 2 1 】。只是这些方法在计算中用到了间 断两侧的数据，从而有振荡和磨损的问题。 2 h a r t e n 的子网格技巧1 1 表明守恒性能够帮助确定间断在网格内的确切位置。 以上的两种观察是这一种守恒型间断跟踪法的思想基础。目前将该方法运用到 一维空间中已经得到了很好的结果 9 ，1 1 ，1 2 】，与传统的间断跟踪法相比较，该方法 具有如下优点： 如果把此方法看作是一种数值技巧，则该技巧可用到任意的间断捕捉法上去， 其作用相当于为基本格式在间断处提供了最优的人工粘性和人工压力，因此比 较适于按间断捕捉法的风格对该方法编制程序； 2 0 0 8 上海大学硕士学位论文 7 由于该跟踪法利用解的守恒性作为跟踪的机制并采用了“幽灵”技巧，从而该 方法的整个计算始终在正规网格上进行，因此克服了一直困扰大多数间断跟踪 法的“小网格 问题，在很大程度上简化了传统的间断跟踪法。编程简洁，稳定 性也很好，并且使得建立一个一般的能处理任意情况的间断跟踪法成为可能； 由于保证了守恒性，因此对间断的相互作用和自发生激波的处理简单、精确和 稳定。 该方法的整个计算始终在正规网格上进行，在很大程度上简化了传统的间断跟 踪法，从而编程简洁，稳定性也很好。另外解的守恒性自动满足。并且在处理间断 相互作用时，该方法显得非常简单，从而使得建立一个一般的强健的间断跟踪法成 为可能。茅德康最近一直在致力于发展二维的e u l e r 方程组的一个一般的强健的守 恒型间断跟踪法，目前已取得相当大的进展，见f 7 ，8 1 。 1 4 本文的工作 本文考虑的是一维e u l e r 方程组，在文 1 0 】中，切向间断处间断两端质量总体守 恒，但是不能保证间断两端分别守恒。本文在切向间断处间断网格上引入了两个普 通网格平均，并在以后的每步计算中做了相应的修正，做到了间断两端质量分别守 恒。整个算法的具体过程简要介绍如下： 通过解r i e m a n n 问题得到初始时间层的数值解，并且根据基本格式得出下一时 间层的数值解； 根据守恒性计算间断位置，判断间断位置是否移出原间断网格，如果还在原网 格，则进行下一步的计算，如果间断位置移到相邻网格，则新间断网格上的数 值解做相应修正； 给出了两个间断移动到同一间断网格后的处理方法，并用数值算例来验证我们 算法的有效性。 第二章守恒型跟踪法在切向间断处的算法改进 考虑一维e u l e r 方程组( 1 1 1 0 ) ，e u l e r 方程组是流体力学的重要方程组，其解含 有三类波，因此有三类间断，左向激波，右向激波和切向间断。本文我们要对算法 的切向间断部分进行修改，为此首先简单回顾一下文 1 0 】中关于切向间断部分的算 法。由于切向间断本质上来说是一个质量密度的间断，因此我们将主要回顾算法中 质量密度的计算，本文对算法的修改也主要是在质量密度的计算格式上。 2 1 算法回顾 文 1 0 】中的计算采用正规网格剖分，网格的空间步长为h ，而时间步长7 - 是在 计算中根据c f l 条件求得的，入= 7 - 是网格步长比。假设我们考虑的间断在如 时刻位于第歹1 个网格中。文 1 0 】中数值解是对精确解的网格平均的逼近，记为： 吆= 瞧 哝 吆 1 r 件1 1 2 竺_ ， h j 哟一1 2 p ( x ，t n ) ( ) ( z ，t n ) e ( x ，t n ) 跟踪法需要一个基本格式，它可取为任意一个有限体积格式： 如 ( 2 1 1 ) 孵1 = 略一a ( 露+ 。2 一露一。2 ) ， ( 2 1 2 ) 其中露+ 。2 为上面引言部分介绍的数值流函数，为了表述方便，我们记 层+ 1 2 = ( 瑶+ 1 21 ( 露+ 1 22 ( 霸+ 1 2 ) 3 ( 2 1 3 ) 其中( 露+ 1 2 ) t 表示数值流函数的第i 个分量，i = 1 ，2 ，3 。 假设我们考虑的是一个含有切向间断的网格，在间断网格中，定义了解的左， 右网格平均和普通网格平均，它们分别是对精确解的左，右和普通网格平均的逼近。 这里以质量密度为例，它所对应的左，右和普通网格平均分别记做瞻一，磅+ 和 8 2 0 0 8 上海大学硕士学位论文 9 鳐，它们是对精确解的质量密度的左，右网格和普通网格平均的逼近： 磅一一1 ，f z j l i ，i 1 2 - - 3 1 - - 1 1 。p 一( z ，t n ) 如， 磅+ 皇去z 二二2 p + ( z ，亡n ) d z ， ( 2 1 4 ) 鳍竺私飞 + ，胆九叫胤 这里p 一( z ，t n ) 和矿( z ，t n ) 是时间层上精确解在间断两侧的质量密度函数及其 在间断另一侧的光滑延拓，厶是该时间层上的间断位置。类似可以定义动量和能 量函数的左，右和普通网格平均。 在算法中左，右和普通网格平均的计算如下：首先通过插值分别“光滑地 延 拓间断两侧的数值解到间断的另一侧，用咳；1 ，u j 忆l 一+ 2 ，咳；知表示左侧数值 解在间断右侧的延拓，用咯：七，吆：知+ 1 ，u n ，, - 一 - 1 表示右侧数值解在间断左侧 的延拓。延拓数据用一种修正的l a g r a j l g e 插值得到，它可保证插值所得到的状态物 理上是合理的，具体情况请参考 1 0 】。 间断两侧的数值解及间断网格上的左，右网格平均分别用两侧的数值解及其在 另一侧的延拓值来计算，例如，对质量密度来说，其具体的数学公式如下： 而间断网格上的普通网格平均的计算在左右网格边界上用的是对应的左右两侧的 流，其解计算如下， 在( 2 1 5 ) ，( 2 1 6 ) 式中 材1 = 瞻一入( ( 露矗。) - 一( 靠 二1 2 ) ) ， ( 2 1 6 ) 间断位置n + 1 由以下方程来求得： 丢( 。群1 一如+ 坍材1 + 删= 材1 ， ( 2 工8 ) 事实上切向间断本质上来说是一个质量密度的间断，其间断位置可由质量密度的守 恒性计算得到。 勋q h h 仡 乜一汁o 瑶瑶 一 一 h h 屈 绝一，“ “ 露露 n n 一 一 一 + 磅磅 = = 一 十 噼噼 力 j d 咯咳盯 一 + 嗡嗡 + + 咳嗡 ， | i = h h 您 您 一，汁 瑶瑶 2 0 0 8 上海大学硕士学位论文 1 0 若计算得到的间断位置p + 1 仍在原间断网格中，则如+ - 层的数值解的计算就 已完毕，但是间断可能在传播，在c f l 条件的限制下，间断位置p + 1 可能会移动到 原间断网格的左、右相邻网格中。此时旧的间断网格应该被删除，新的间断网格应 该被建立。 x j l - - 1 2px j l + l 2 n + 1 x j i - i 2px j l + 1 2巧1 + 3 2 ( b ) 图2 1 ( a ) 间断n + 1 传到了左相邻网格中； ( b ) 间断p + 1 传到了右相邻网格中 不妨以p + 1 传到左边相邻网格的情况为例，即有： l 一3 2 卅1 1 1 2 ， ( 2 i 9 ) 成立，此时b 。一3 2 ，。一1 2 】成为新的间断网格，而b 。一1 2 ，。+ i 2 】变成普通网格。 在这一过程中， 巧i - 3 2 ，吻，一l 2 】上的原普通网格平均吆品变成了新的间断网格的 左网格平均，而b 如一1 2 ，z ，+ 1 2 】上的原右网格平均变成该网格上的普通网格平均。 新旧间断网格上数值解的的修正按以下方式保持守恒性：新、旧间断网格上的普通 网格平均之和在修正前后相等。而对于新的间断网格上的右网格平均则由其右侧各 相邻网格中普通网格平均向左延拓而得到。我们仍以质量密度为例，其数学表达如 2 0 0 8 上海大学硕士学位论文 1 1 卜： ， i 稿一：= 约p j 。+ _ 1 l 壁譬壁+ 群1 一醛1 + ，( 2 1 1 0 ) l 材1 ：= 群1 + ， 7 i 刀0 + ：= 来自右侧的插值值 若间断位置移到原间断网格的右相邻网格中，修正可以类似得到。显然如此设计的 数值方法保证了间断两侧总质量的守恒。 以上讨论的是跟踪单个间断的情况，但在实际问题中，往往会同时存在多个不 同类型的间断，且它们之间会相互作用，并由此产生新的间断。和传统的跟踪法不 同的是，文【1 0 1 中的算法在处理间断的相互作用时，大量的采用了“幽灵”技巧。通 过引进虚拟的状态来克服因空间不够所造成的困难。其中最为重要的是所谓的“网 格堆积 技巧。下面先从如何处理两个间断的相互作用说起，而处理多个间断的相 互作用可通过两个两个的处理间断相互作用来完成。设t = t n 时，两间断的位置 分别为6 和矗，并且6 靠，邑 x j ，一1 2 ，巧，+ 1 2 ，靠 x j 。一1 2 ，巧。+ 1 2 】，这里 一1 2 ，。+ ，2 和【。一1 2 ，巧。+ 1 1 2 是两间断分别对应的间断网格。 2 0 0 8 上海大学硕士学位论文 1 2 & 靠 7 舌n + 1 t n 1 1 2x i + 1 2 x j l + 1 2 6靠 亡n + 1 t n 巧1 一i 2x j i + i 2x j l + 1 2 6靠 亡n + 1 瓮n 一1 2x j l + 1 2 z n 一1 2 x j l + 1 2 6矗 亡n + 1 t 竹 z h i 2z h + 1 2x j l + l 2 图2 2 两间断移入同一个网格的几种形式 2 0 0 8 上海大学硕士学位论文 1 3 当两间断的间断位置6 和矗移入同一网格b 。一1 2 ，巧，+ 1 2 】之后，此时往往 仍然成立 矗，为此我们设计了一种“网格堆积”的技巧，使得在此情况下网 格内解的结构能够被分辨出来。“网格堆积”的技巧认为，当已 矗时，在网格 b 。一1 2 ，。+ 1 2 】中仍然有两个独立的间断网格，分别对应着左间断6 和右间断6 ， 它们各自有自己的左、右网格平均和普通网格平均，且有一个隐含的中间状态u + 将 它们相连，以质量密度为例来说，p 。既是左间断网格上的右网格平均的质量密度， 又是右间断网格上的左网格平均的质量密度。而b ，一1 2 ，。+ 1 2 】上的实际的普通 网格平均的质量密度为 p j - i + 1 = 甜+ 约p j ，+ ，r a p 。 ( 2 1 1 1 ) 其中吩p j “+ l 和稼分别为左、右两间断网格的普通网格平均的质量密度。上式这样 的定义保证了守恒性。 下面来介绍网格堆积的形成，以右边间断传到了左间断网格中为例：在t n 时 刻，间断6 和间断矗所对应的间断网格为相邻网格，分别称之为左，右间断网格。 而在t n + l 时刻，间断靠传到了左间断网格 ，一1 2 ，。+ 1 2 】中，此时左间断网格没 有发生变化，其上的左、右网格平均和普通网格平均无需修正。而右间断网格却发 生了变化，在t n 时刻，右间断矗所对应的间断网格为 x j t + 1 2 ，巧2 + 1 2 】，而在t n + l 时刻，右间断靠所对应的间断网格变为 巧。一1 2 ，巧，+ 1 2 】此时隐含的中间状态应 取为p ，= 甜+ ( 材+ 为左间断的右网格平均的质量密度) 。此时由于右间断网 格发生了移动，因此其质量密度的相关量进行修正如下： i 彬一：= p + ， p j + r l ：一- - - p 一稿一p i n + l , + ，( 2 1 1 2 ) l 乃n ，+ l = 群+ ， 、 【约p j ，+ ，r l + ：= 来自右侧的插值值 同理，如果左间断传到右间断网格中、或者左右两个间断同时传出自己的原间 断网格而传到原两间断网格之间的网格中、或者左右两个间断分别传到各自原来所 对应的间断网格的右相邻和左相邻网格中，相应的修正可以类似得到。对于堆积后 的间断网格上数值解的计算以及间断的相互作用，详细情况请参考3 1 0 1 。 2 2 改进后算法 为了保证切向间断两侧气体分别守恒，我们必须把两侧质量密度单独计算，但 是动量和能量与以前 1 0 】文中一样计算。为此我们在切向间断处间断网格上放弃原 2 0 0 8 上海大学硕士学位论文 1 4 来的普通网格平均，引进了左右普通网格平均磅一，哪! ”，它们定义为： 张一竺去e m p ， ( 2 2 ) 磅+ 竺j 1 。x j l + i 2 矿( z ，t n ) 如， ( 2 2 1 4 ) 其中p 一( z ，t n ) ，矿( z ，t n ) 意义如( 2 1 4 ) 式下面所述，同样是如时间层上精确解在间 断两侧的质量密度函数及其在间断另一侧的光滑延拓，p 为精确解的间断位置。 由于改进后的算法在切向间断处间断网格上放弃了原来的普通网格平均，因此 在计算过程中切向间断网格上质量密度普通网格平均的计算将由左右普通网格平 均的计算来代替。左右普通网格平均的计算公式推导如下：假设厶时刻精确解的间 断位置为p ，t n + l 时刻间断位置移动到了p + 1 。设流体的运动速度为钉( z ，t ) ，间 断线的方程为z = z ( 亡) ，由于我们考虑的是切向间断处的数值解，间断的传播速度 就是间断处流体的移动速度，因此在间断线上成立： i d x ：u ( 。( 亡) ，亡) ( 2 2 1 5 ) 面。u 艺川 瞄z l b j 图2 3 网格 x j ，一1 1 2 ，x j ，+ 1 2 】 t n ，亡叶1 】和其中的间断线。 假设我们考虑的切向间断位于网格 。一1 2 ，乃，+ i 2 t n ，亡1 】中( 如上图) ，对e u l e r 方 程组( 1 1 1 ) 的第一个方程来说， ( ( ) ( z ，亡) ) z + ( p ( x ，亡) ) t = 0 ( 2 2 1 6 ) 设间断左侧区域为q ，那么 z ( ( 训叫m + ( 出) 廊出) = 。， ( 2 2 1 7 ) 2 0 0 8 上海大学硕士学位论文1 5 由格林公式，可以得到： j r 厶( ( ) ( z ( 亡) ，亡) ) 茁+ ( p ( z ( 亡) ，亡) ) t d x d t ) = 正m ( 一p ( z ( 亡) ，t ) d x + ( ) ( z ( 亡) ，t ) ) d t ( 2 2 1 8 ) = 位( ) + 局( ) + 局( ) + 局( ) ， 根据( 2 2 1 7 ) ，有： 厶( ) + 厶( ) + 厶( ) + 厶( ) _ o - ( 2 2 1 9 ) 首先看式( 2 2 1 9 ) 中左侧部分第一项： 厶( ) 2 厶( _ 以叫d 一) 如+ ( 舢) ( 叫幻一) 鳓2 厶1 - 1 ：一从如) 如 ( 2 2 2 0 ) ，厂，( n 而对于左侧部分第二项来说： 厶( 。) 2 厶( 一p ( z ( 州如+ ( ) ( z ( 州出) ( 2 2 2 1 ) 将( 2 2 1 5 ) 式代入上式，可以得出式( 2 2 2 1 ) 的值为零，该式从数学角度证明了从k 时 刻到t n + 1 时刻没有气体流过切向间断，问断两侧气体质量没有交换，这也是我们算 法改进的思想基础。 类似的，我们得到 苫五( ) = j 苫右( 一p ( z ( 亡) ，t ) d x + ( ) ( z ( t ) ，亡) 班) = j f p x j + l ，- 1 7 2 ( 一p ( z ，t 竹+ 1 ) 如) ( 2 2 2 2 ) = 鹱二0 p ( z ，亡州) 如) ， 局( ) = 局( 一j d ( z ( 亡) ，t ) d x + ( ) ( z ( 亡) ，t ) d t ) = 盘。( ( ) ( 巧。一1 2 ，t ) d t ) ( 2 2 2 3 ) = 眨“( ( 一) ( 巧。一1 2 ，t ) d t ) 将( 2 2 2 0 ) ，( 2 2 2 1 ) ，( 2 2 2 2 ) ，( 2 2 2 3 ) 代x ( 2 2 1 9 ) ： o 2 鹱- l ：一p ( 酬蚺群：( p ( x , t n + 1 ) d x )( 2 2 2 4 ) + 露“( ( 一) ( z 沪1 2 ，t ) ) d t ， 、7 把上式右端第二项移到等号另一侧，然后等式两边同时乘以一丢，得到： 丢：施 “如= 丢出如+ 丢r n + l p c x j l _ i 2 , t m x a _ l 2 , t 皿 ( 2 2 2 5 ) 2 0 0 8 上海大学硕士学位论文 1 6 此式表明切向间断左侧气体质量的变化只和左边界的流入量有关。根据左普通网格 平均的定义( 2 2 1 3 ) 及( 2 1 7 ) ，为此我们的计算格式可以表示为： 乃- n ，+ 1 一= 磅一+ 入( f 2 2 1 1 2 1 ( 2 2 2 6 ) 同理可以得到右普通网格平均的计算公式： 乃- n 。+ 1 + = 磅+ 一入( 寇矗2 ) 1 ( 2 2 2 7 ) 得出了左右普通网格平均的计算公式后，下面再简要介绍改进后的算法： 我们所作的修改只是针对切向间断的质量密度，因此，只有切向间断所对应间 断网格上的质量密度的数值解需要改变，激波以及切向间断处动量和能量的计算与 文【1 0 1 中一样。以下叙述的重点是切向间断处的质量密度的改进。首先考虑只含有 一个切向间断的情况，t 竹时刻此间断位于第j 1 个网格中。在该间断网格中，质量密 度对应着左、右网格平均，普通网格平均由左、右普通网格平均代替。上述各网格 平均的定义见( 2 1 4 ) ，( 2 2 1 3 ) ，( 2 2 1 4 ) 。算法中左、右网格平均和普通左、右网格 平均的计算由( 2 1 5 ) ，( 2 2 2 6 ) ，( 2 2 2 7 ) 得出。 接下来考虑间断位置的计算。切向间断处的间断位置仍由质量密度的守恒性计 算得到，由于我们在间断网格上引入了两个普通质量网格平均，根据左、右普通网 格平均可以得出两个间断位置。例如对左普通网格平均来说，根据其定义( 2 2 1 3 ) 及 采用对数值解的分片常数重构可计算得数值的间断位置r 1 铲1 = 匆，_ 1 2 + 簪 ( 2 2 2 28 ) 同理，利用右普通网格平均可以得到间断位置器“： g + 1 = x j l + l 2 一舞尼 ( 2 2 2 29 ) 最终间断位置p + 1 取为r 1 和等+ 1 的平均， p 1 书“埘1 ) = x j l - - - 互1 ( 转一耢m ( 2 2 s o ) 若计算得到的间断位置p + 1 仍在原间断网格中，则t n + 1 层的数值解的计算就已完 毕。 同样的，也可能会出现p + 1 移到左边或者右边相邻网格的情况，具体的情况分 析与上节算法回顾的相应部分类似，只是数值解的修正会有一定的改变。以间断位 2 0 0 8 上海大学硕士学位论文 1 7 置移到原间断网格的左相邻网格为例，修正后的数值解用数学式子表示为： p i n + l , 一：= 稿， 材1 ：= 群1 + ， 露3 + ：= 来自右侧的插值值， ( 2 2 3 1 ) 稿一：= p i n + l , 一+ 群1 一， 稿+ ：= 群1 + 一材1 + 如果p + 1 移到右边相邻网格，相应的修正类似可得。 以上讨论的是跟踪单个间断的情况，实际情况往往会复杂的多，比如，有可能 所考虑的问题中含有不止一个间断，而且两个间断移到同一个间断网格中，对于这 种情况我们采取以下方式进行处理： 假若在时刻t n 两个间断位置6 和矗移入同一网格b 。一1 2 ，巧，+ 1 2 ，由于 此时往往仍然成立矗 矗，因此，我们仍然采用“网格堆积”技巧，认为网格 【，一1 2 ，+ 1 2 】中有两个独立的间断网格，分别对应着左间断自和右间断靠。我们 算法的改进是在切向间断处的间断网格上添加了左、右普通网格平均，由于两个量 的代数和即为文f 1 0 1 中间断网格上质量密度的普通网格平均，于是，在堆积后的间 断网格上，切向间断处质量密度的普通网格平均的修正也由左普通网格平均或者右 普通网格平均的修正来代替。例如，如果靠为一切向间断，那么其右普通网格平均 没有变化，左普通网格平均则要做一定的修正。仍然用p 。表示中间状态的质量密度， 数值解的修正用数学式子来表示： 如果左间断6 为一切向间断，相应的修正可以类似得到。 间断之间相互作用的处理同文 1 0 1 ，在此不再赘述。 ( 2 2 3 2 ) 甜 一 的稿，蒯+ 。自风风卅卅来 j l j i 叫 j | 秽秽稼秽 第三章数值算例 在本章所有数值算例中，基本格式是文 5 】中的两步r 肌n g 争k u t t a 法的二阶e n o 格 式，c o u r a n t 数选为0 5 。 算例1 ：单介质流的激波管问题 1 4 】 初值为 钆。c z ，= ：二三三三： c 3 。1 ， 其中 p t = 1 0 ，胁= 0 1 2 5 ， 口l = 钐= 0 0 ， ( 3 0 2 ) p t = 1 0 p r = 0 1 此例中间断两侧均为1 = 1 4 的气体。图3 1 是用1 0 0 个点算得的t = 0 4 时刻的数 值结果，其中( a ) 是密度，( b ) 是速度，( c ) 是压力，精确解是通过解r i e m a n n 问题 得到的。 一 算例2 ：多介质流s o d 激波管f 司是亟 2 8 ，2 3 】 间断两侧分别是7 为1 4 和，y 为1 7 的气体。初值为 “。c z ，= ：2 三三兰： c 3 。3 ， 其中 p t = 1 0 ，胁= 0 1 2 5 ， q t = 钐= 0 0 ， ( 3 0 4 ) p t = 1 0 ，肼= 0 1 图3 2 是在t = 0 4 时刻用1 0 0 个点算得的数值结果，其中“实线 表示通过解 r i e m a n n 问题得到的精确解，而“点”代表我们算得的数值结果。从图中可以看出， 算得的数值结果是理想的。 算例3 - 单介质流的b l a s t w a v e f 驱 2 5 】 1 8 2 0 0 8 上海大学硕士学位论文 1 9 此例的初始条件是这样的：在距离为1 的两个反射固壁之间，有一包含三个常数 状态的7 一l a w ( 7 = 1 4 ) 气体，并且气体处于静止状态( 即速度为o ) ，密度处处为1 ， 从左到右三个常数状态的压力分别为1 0 0 0 ，0 0 1 和1 0 0 即e u l e r 方程组的初值为 其中 l u l ，0 0 z 0 1 ； u o ( x ) = 钆m ，o 1 z o 9 ； ( 3 0 5 ) iu r ，0 9 z 1 0 p t = 胁= p r = 1 ， g l = q m = q r = 0 ， ( 3 0 6 ) p l = 1 0 3 ，p m = 1 0 一，p r = 1 0 2 在边界x = 0 和z = 1 0 处是反