（应用数学专业论文）统一坐标系下多介质流体力学计算方法研究.pdf

摘要 近年来，h u i 等提出了统一坐标系的思想，并对单介质流体力学问题进行了 数值模拟，取得了比较理想的结果。在统一坐标系下，流体动力学的各物理量看 成是时间和伪粒子( p s e u d o p a r t i c l e s ) 的某静固有特征的函数，伪粒子的运动速度 为h q ，q 是流体质点的速度。统一坐标包括e u l e r 坐标和l a g r a n g e 坐标两种特殊 情况，当h = 0 时为e u t e r 坐标，h = t 时为l a g r a n g e 坐标。在二维无粘流动中， 当自由函数h 的选择为保持网格角时，可以避免e u l e r 坐标系下接触间断处的过 度数值耗散以及l a g r a n g e 坐标系下严重的网格扭曲，同时在接触间断附近仍然 具有较高的分辨率。本文的主要目的是将统一坐标系下二维单介质的流体力学计 算扩展到多介质的流体力学计算。我们详细地讨论了统一坐标系下基于，模型的 二维扩展欧拉方程的形式以及双曲性，给出了统一坐标系下维数分裂后的一维扩 展欧拉方程r i e m a n n 闯题的精确解，并采用带有m u s c l 修正的g o d u n o v 方法 求解二维扩展欧拉方程，通过数值实验表明了统一坐标系下多介质流体力学计算 的优势。同时，我们给出了确定h 值的一种不同于网格保角的方法，这种方法是 种基于网格节点之间的吸引和排斥来重新分布节点的移动网格方法，方法简 单，易于实现，数值结果表明，不仅可以显著地提高计算效率，而且仍然具有较 高的流场分辨率。此外，我们针对一维情况给出了统一坐标系下扩展欧拉方程的 j a c o b i a n 矩阵及其特征值私特征向量，同时采用一种二阶精度的两步t v d 格式 来求解统一坐标系下的一维多介质流体力学方程组，这种格式不需要通常t v d 格式中所用的特征分鳃，数值实验表明了方法的简单和实用，也验证了统一坐标 的优势。 关键词：统一坐标，多介质漉体，扩姨e u t e r 方程，移动网格方法，g o d u n o v 型格式，t v d 格式 a b s t r a c t r e c e n t l y , h u ii n t r o d u c e d as o - c a l l e du n i f i e dc o o r d i n a t e s y s t e m f o rc o m p r e s s i b l e s i n g l e f l u i dc o m p u t a t i o n s ，i nt h i ss y s t e m ，t h ef l o wv a r i a b l e sa r ec o n s i d e r e d t ob ef u n c t i o n so f t i m ea n do fs o m ep e r m a n e n ti d e n t i f i c a t i o no fp s e u d o p a r t i c l e sw h i c hm o v ew i t hv e l o c i t y h q ，qb e i n gt h ev e l o c i t yo ff l u i dp a r t i c l e s i ti n c l u d e st h ee u l e r i a nc o o r d i n a t e sa ss p e c i a l c a s ew h e nh = 0a n dt h el a g r a n g i a nw h e nh = 1 f o rt w o d i m e n s i o n a li n v i s c i df l o w , t h e f l e ef u n c t i o nhi sc h o s e ns oa st op r e s e r v eg r i da n g l e s t h i sr e s u l t si nac o o r d i n a t es y s t e m w h i c ha v o i d se x c e s s i v en u m e r i c a 【d i f f u s i o na c r o s ss l i pi i n e s i ne u l e r i a nc o o r d i n a t e sa n d a v o i d ss e v e r eg r i dd e f o r m a t i o ni nl a g r a n g i a nc o o r d i n a t e s ，y e ti tr e t a i n ss h a r pr e s o l u t i o no f s i i pl i n e s i nt h ed i s s e r t a t i o n ，w e e x t e n dt w o d i m e n s i o n a ls i n g l e f l u i dc o m p u t a t i o n st o m u t t i m a t e r i a lf l u i dc o m p u t a t i o n si nt h eu n i f i e dc o o r d i n a t es y s t e m ，w ed i s c u s st h ef o r ma n d h y p e r b o l i c i t yo ft h e ，一b a s e de x t e n d e de u l e re q u a t i o n si nt h eu n i f i e dc o o r d i n a t e ，a n dg i v e t h es o l u t i o nt ot h er i e m a n np r o b l e mf o rt h e1 - de x t e n d e de u l e re q u a t i o n so b t a i n e df r o mt h e 2 - de x t e n d e de u l e re q u a t i o n sa f t e rd i m e n s i o n a ls p l i t t i n g t h e2 - de x t e n d e de u l e re q u a t i o n s a r en u m e r i c a l l ys o l v e du s i n gt h eg o d u n o vs c h e m ew i t hm u s c l u p d a t e t h en u m e r i c a lt e s t s s h o wt h ea d v a n t a g e so ft h eu n i f i e dc o o r d i n a t e m o r e o v e r , w eu s eam o v i n gm e s hm e t h o d o t h e rt h a np r e s e r v i n gg r i da n g l e sf o rd e t e r m i n a t i o no fh n u m e r i c a lr e s u l t ss h o wg o o d e f f i c i e n c ya n da c c u r a c yo ft h em e t h o d i nt h el a t t e rp a r to ft h ed i s s e r t a t i o n ，w eg i v e t h ej a c o b i a nm a t r i x ，e i g e n v a l u e sa n dc o r r e s p o n d i n g e i g e n v e c t o r s f o rt h e1 一d e x t e n d e de u l e re q u a t i o n si nt h eu n i f i e dc o o r d i n a t e ，a n ds o l v et h ee q u a t i o n sa p p l y i n g as e c o n d o r d e ra c c u r a t e ，t w o s t e pt v ds c h e m e t h es c h e m ed o e sn o tn e c e s s i t a t et h e c h a r a c t e r i s t i cd e c o m p o s i t i o n so ft h eu s u a lt v ds c h e m e s af e wn u m e r i c a lr e s u l t s s h o ws i m p l i c i t ya n dp r a c t i c a b i l i t yo ft h es c h e m e ，a n ds h o wt h ea d v a n t a g e so ft h e u n i f i e dc o o r d i n a t e k e yw o r d s ：u n i f i e dc o o r d i n a t e ，m u l t i m a t e r i a lf l o w ，e x t e n d e de u l e re q u a t i o n s ， m o v i n gm e s hm e t h o d ，o o d u n o v t y p es c h e m e ，t v ds c h e m e 独创性声明 本人声明所呈交的学位论文是本人在导师指导下进行的研究工作及取得的 研究成果。据我所知，除了文中特别加以标注和致谢的地方外，论文中不包含其 他人已经发表或撰写过的研究成果，也不包含为获得中国工程物理研究院或其他 教育机构的学位或证书使用过的材料。与我一同工作的同志对本研究所做的任何 贡献均已在论文中作了明确的说明并表示谢意。 学位论文作者签名：髫鹇秀 签字日期：p 修年多月拥 学位论文版权使用授权书 本学位论文作者完全了解并接受中国工程物理研究院研究生部有关保存、使 用学位论文的规定，允许论文被查阅、借阅和送交国家有关部门或机构，同时授 权中国工程物理研究院研究生部可以将学位论文全部或部分内容编入有关数据 库进行检索，可以影印、缩印或扫描等复制手段保存、汇编学位论文。 学位论文作者签名：嘴鹏秀 导师签名 签字日期：砖年z 月矽咱签字日期：力百年占月巧目 中国 一程物理研究院博士学位论文 第一章绪言 本章的主要内容是简要介绍可压缩多介质流体力学计算方法的研究状况， 第二节介绍了统一坐标的思想，第三节介绍多介质流体力学问题模型，最后一 节介绍本文研究的主要内容。 1 。1 可压缩多介质流体力学计算方法的研究概况 计算流体力学在2 0 世纪7 0 年代以来有了突飞猛进的发展，推动这一发展的 原因一方面是实际问题的需求，另一方面是计算技术的飞速发展和高速巨型计 算机的出现，使其发展成为门独立的学科。计算流体力学是多种领域的交叉 学科，它所涉及的学科有流体力学、偏微分方程的数学理论、数值方法、计算 机科学等。它的发展进一步促进了这些学科的发展。 任何流体运动的动力学特性都是由质量守恒律、动量守恒律和能量守恒德所 确定的。这些基本定律可由数学方程组( 偏微分方程组或积分方程组) 来描述， 如歇拉方程( e u l e r 方程) ，纳维一斯托克斯方程( n a v i e r s t o k e s 方程，简称为 n - s 方程) 等。计算流体力学就是利用数值方法通过计算机求解描述流体运动的 数学方程，揭示流体运动的物理规律，研究定常流体运动的空间物理特性和非 定常流体运动的时一空物理特性。 计算流体力学的兴起推动了流体力学研究工作的发展。从1 6 8 7 年牛顿发现 牛顿定律以来，到2 0 世纪s o 年代初，研究流体力学的主要方法是两种：是 实验研究，以地面实验为研究手段；另一是理论分析方法，利用简化的流动模 型假设，给出所研究闯题的解析解。理论工作者在研究流体运动基本规律的基 础上，提出了各种简化流动模型，给出了系列解析解和数值方法。这些研究 成果推动了流体力学的发展，奠定了今天计算流体力学的基础，很多方法仍是 目前解决实际问题时常采用的方法。然而，利用解析方法求解的数学问题的解 析解和近似解的范围是极其有限的，一般只能考虑一些很简单的问题。利用实 验方法来测量数据是有限的而且来之不易。仅采用这些方法研究复杂非线性流 中国工程物理研究院博二e 学位论文 2 体运动规律是不够的，它已不能满足2 0 世纪5 0 年代开始高速发展起来的近代 科学技术的需求。 今天随着实际的需要，计算流体力学己发展成为流体力学第三种研究方法。 它的兴起促进了实验研究和理论分析方法的发展，将实验研究和理论分析方法 联系起来，为简化流动模型的建立提供了更多的依据，使很多简化方法得到了 发展和完善。计算流体力学采用它独有的数值模拟方法直接求解描述流体运动 基本规律的非线性数学方程组来研究流体运动的物理特性，因此在某种意义上 比理论与实验对运动过程认识得更为深刻，更为细致，不仅可以了解运动的结 果，而且可以了解运动整体的与局部的细致过程。应该指出，要建立正确的数 学方程，特别是对于复杂的流动问题，如燃烧、多相流、湍流等，还必须与实 验研究和理论分析方法相结合。更重要的是计算流体力学所求解的非线性偏微 分( 积分) 方程组，其数值方法的现有数学理论尚不够充分，严格的稳定性分 析、误差估计和收敛性证明等理论工作的发展还跟不上数值模拟方法的进展。 所以在计算流体力学中，一方面仍必须依靠对一些简单的、线性的、与原有问 题有一定关系的数学方程进行严格的数学分析，依靠启发性的推理，分析非线 性问题，给出数值解的理论依据；另方面，依靠对线性和非线性方程的数值 实验及数值解与地面试验值或他人典型算例的计算结果的比较和物理特性分 析，验证计算结果，进步改进计算方法。所以实验研究、理论分析和数值模 拟方法是研究流体运动的三种基本方法，它们的发展是相互依赖相互促进的。 计算流体力学的发展与计算机技术的发展直接相关。这是因为采用数值方法 可阱模拟物理问题的复杂程度，解决问题的广度、深度和所能给出数值解的精 度都与计算机的速度、内存和外围设备( 如图象输出的能力) 直接相关。一般 来说，只有计算机的速度、内存和外围设备达到新的水平时，才会有计算流 体力学新阶段的出现。随着计算技术的提高，巨型计算机的出现，计算流体力 学求解问题的深度和广度不断发展，它不但可用于研究一些物理问题的机理， 解决实际流动中的各种问题，而且可用于发现新的物理现象。 描述流体运动的方程组是拟线性双曲型方程组，对于拟线性双曲型方程组， 一般来说，不管初值如何光滑，解在有限时间内可以发生间断，并且间断也可 能消失。这种间断的产生与消失反映了流体运动中冲击波间断的产生与消失， 中国一l ：程物理研究院博j 一学位论文 这种特性使求解流体力学运动方程组有它特殊的问题与困难，具有挑战性。 二维非定常可压缩理想流体力学计算方法的研究开始于五十年代中期，到了 六十年代，二维流体力学计算方法的研究进入了个新的时期，人们发表了大 量的关于计算格式的文章，并研制了许多计算程序。由于二维非定常的流体运 动是很复杂的运动，在流场上可以发生扭曲、涡流与滑移等现象。在有多种介 质时，这些机制会显出复杂性来。不同介质的界面运动的不稳定性使不同介质 混淆起来，而一维流体运动因为流体介质前后按序不能相互超越而简单得多。 这些情况使得多介质流体力学运动方程数值求解时，二维问题比一维问题复杂 很多。对于一维流体力学问题的计算，往往一种有限差分格式可以计算多种流 动问题，但是对于二维流体力学问题的计算，很难有一些统一的格式可以计算 好各种闽题，对于不周的模型往往需要采用不同的格式进行计算，有时甚至对 于一个运动中的模型的不同发展阶段，还要采用不同的格式，才能把整个运动 过程计算出来。 在二维流体力学计算中，常用的两种坐标系是e u l e r 坐标系和l a g r a n g e 坐标 系。欧拉坐标系是在固定的空间坐标系中讨论流体运动，因此很好地保持了网 格的几何性质( 均匀性、正交性等) ，适合于有大变形的流体运动场的计算，不 会出现网格相交的问题，但是由于流体微团会穿过网格，因此在对接触间断的 计算过程中，会造成较严重的数值耗散影响计算的准确性。拉格朗日坐标系是 跟踪流体质点来研究流体运动，因此可以用来计算包含多种物质的系统，而且 能很好地计算接触闽断问题，不同物质间的界面能清晰地表示出来，但在计算 大变形的流场时，会出现网格的严重变形，甚至可能出现网格相交，使得计算 不能进行下去。而且网格的扭曲又会引起计算误差的增长，使得l a g r a n g e 方法 失去精确性。人们曾经设想过各种各样的办法来克服网格相交。其中的类办 法是设想用某种人为粘性来阻止网格的扭曲和变形，以达到推迟和防止网格相 交的目的。尽管这类办法对于防止网格相交有一定效果，但是对于物理真实解 的扰动和变形，这类办法所加入流场的“人为粘性”可能改变流场的性质，歪 曲物理图象。 由于e u l e r 方法和l a g r a n g e 方法各有其优缺点，因此很自然地发展了一些将 e u l e r 方法和l a g r a n g e 方法相结合的方法( 7 ，1 6 ，1 5 ，6 7 ，5 8 ，5 9 1 ，这些方法发挥各 中国j ：程物理研究院博士学位论文 4 种方法的长处而避免其弱点。f r a n k ，l a z a r u s ( 1 9 6 4 ) ( 5 8 1 ) 提出了以一个空间坐 标取作固定的e u l e r 坐标，而另一个空间坐标耿作l a g r a n g e 坐标的混合 e u | e r l a g r a n g e 方法：n o b ( 1 9 6 4 ) ( i s 9 ) 的耦合e u l e r l a g r a n g e 方法( c e l 方 法) 则是将求解区域划分为若干个子区域，在一些子区域上用e u l e r 方法，而在 另一些子区域上用l a g r a n g e 方法。另外还有在实际应用中采用比较多的任意 l a g r a n g e e u l e r 方法( 简称a l e 方法 7 ，1 6 ，15 ，6 7 ) ，这种方法可以象普通的 l a g r a n g e 方法样，让网格跟随流体一起运动，或者象e u l e r 方法一样，让网格 固定，也可以让网格以任意方式运动，因此和纯粹的l a g r a n g e 方法相比，有更 大的优点，可以处理较大畸变的流体运动，而又比纯粹的e u l e r 方法提供更细致 的结果。但是这种方法在执行的过程中，需要进行网格的重新划分和网格之间 输运量的计算，而这些计算中采用的插值方法所带来的误差影响了结果的准确 性( 【6 ) 。 在过去的数十年里，尽管计算机的速度和内存都有很大提高，但很多问题的 数值模拟仍需要在很多简化下才可完成。也就是说，现有的计算机能力对付过 多的自由度仍然有很多困难，特别是三维问题。在这种情形下，自适应算法应 运而生，并在理论研究和实际应用上得到了广泛的重视。如果个偏微分方程 的解有足够的光滑度，则致网格就可以给出满意的解。但是，也有些很重 要的问题，其解的光滑性并非很好，比如说解间断或解有大梯度的情形。在这 种情形下，一致网格的计算将是非常昂贵的，而自适应算法则是一种很有效的 工具。很多实际问题在局部有奇性，需要很多的计算网格点，平均分布网格将 造成不必要的计算时间和数据存储上的浪费，合理分布网格对于高效和精确计 算起着极其重要地作用。 通常有两种方法可以达到自适应效果，一种是局部加密的方法( h r e f i n e m e n t ) ( 例如【3 4 ，3 5 ，3 6 ，3 7 1 ) ，一种是移动嗣格方法( m o v i n gm e s hm e t h o d ) ( 例如 【2 3 ，2 4 ，2 5 ，5 2 ，5 3 ，5 7 】) ，前者的想法是在原有网格的基础上，通过近似解的某种后 验估计和解在局部的误差表现，在有需要的地方做局部加密。 移动网格方法的基本思想是保持求解过程中网格节点数不变，但网格节点的 位最随着时间的变化而变化，并且将较多的网格点移动到解的性质较奇异的地 方。通过这种合理分布网格点，可以使得解变化较大的局部区域有较多的网格， 中国j - 程物理研究院博士学位论文 从而使整体的误差减小，使数据存储量小，计算速度加快。对于移动网格方法， 虽然得到合适的控制网格生成的方程比较困难，同时增加了对移动网格偏微分 方程的求解，但是它在数值计算中的优势也是很显然的，比如，它的原理实现 起来是简单的，并且很多基于固定网格的求解偏微分方程的方法可以很容易地 扩展到移动网格中。 对于移动网格方法般可以分为两类( 【3 3 】) ：一类是基于网格位置的方法， 这种方法可以直接确定网格的位置；另一类典型的方法是基于网格速度的方法， 这种方法先由网格控制方程确定出网格的速度，然后根据速度得到网格的位置， 并且所求解的偏微分方程通过组坐标变换，在计算平面上来计算。 在移动网格方法中有很多值得注意的问题，其中一个重要的问题就是光滑 化技术的应用。在从粗网格到纫网格移动的过程中，网格的变化应该是平稳的， 较缓慢的，过分快速的网格变化会减低计算的精度，但过分缓慢的变化又起不 到自适应的效果。如何应用一些光滑化技术是提高移动网格计算效果的一个重 要手段。一个常见的光滑化技术就是网格的平均化，即把每个区间的左右端点 ( 或邻近端点) 进行合理的平均化。另外，这一光滑化过程可能还需要进行适 当的迭代处理，也就是既，每个时间层上形成的网格需要通过几次( 比如说3 到4 次) 迭代磨合最后形成。 近年来，h u i 等( 【2 ，3 ，9 ，8 ，2 0 ) 提出了统一坐标系的思想，并对单介质流 体力学问题进行了数值模拟，取得了比较理想的结果。统一坐标方法可以说是 一种基于网格速度的移动网格方法。在统一坐标系下，流体动力学的各物理量 看成是时间和伪粒子【p s e u d o p a r t i c l e s ) 的某种固有特征的函数，伪粒子的运动 速度为h q ，q 是流体质点的速度。统一坐标包括e u l e r 坐标和l a g r a n g e 坐标两种 特殊情况，当h = 0 时为e u l e r 坐标，h = 1 时为l a g r a n g e 坐标。在统一坐标系中， 由于计算网格是跟随伪粒子变化的，丽h 值限制在0 h l 之间，伪粒子的运动 方向与流体微团的运动方向是一致的，因此在计算二维无粘流动中，可以避免 e u l e r 坐标系下接触间断处的过度数值耗散。当自由函数h 的选择为保持网格角 时，可以避免l a g r a n g e 坐标系下严重的网格扭曲，同时在接触阊甑附近仍然具 有较高的分辨率。而且，由于在压缩波区域，计算网格会跟随伪粒子的运动而 中国j ：程物理研究院博士学傍论文6 自动加密，团此，在计算激波时，统一坐标系比欧拉坐标系有所改善。 目前，统一坐标方法仅仅针对单介质，本文的主要目的就是将统坐标系 下二维单介质的流体力学计算扩展到多介质的流体力学计算，希望得到好的计 算结果，从而为二维多介质流体力学计算问题提供一种途径。 1 2 统一坐标介绍 在流体力学中，常用的两种坐标系是欧拉坐标和拉格朗日坐标。欧拉坐标系 是在固定的空间坐标系中讨论流体运动，流体力学中的各物理量看成是时间和 固定的空间坐标的函数。拉格朗日坐标系是跟踪流体质点来研究流体运动，流 体力学中的各物理量和空闻位置看成是时间和流体质点的某种固有特征的函 数，例如，可以耿质点的初始位置作为坐标，也可以取质点的某特定质量i l l 作 为坐标，等等。 现在我们将欧拉坐标系中的空间坐标( z ，y ，z ) 和时间f 通过坐标变换变为空 间坐标( ；，q ，毒) 和时阊 ，变换如下( 【2 】) ： f d t = d 2 j 出= “以+ 4 蟛+ 工期+ p 西( 1 1 ) 方= h v d 2 + b 霹+ m d r t + q d ( 1 d z = m 烈十c d f + n d r l + r d ( 其中，“，v 和w 分别表示流体速度q 在r ，y 和z 方向的分量。 我们假定一种不同于流体质点的微粒，称为伪粒子( p s e u d o p a r t i c l e s ) ，它的 运动速度为 q ，令： 生：旦+ h “旦+ h v 旦+ h w 旦 ( 1 2 ) d t a t缸卸8 z 则( 1 2 ) 表示跟随伪粒子的全微商，而且，我们很容易得到： 警= o ，等= o ，警= o d td t d t 中国r 程物理研究院博士学位论文 7 这说明，在伪粒子的流动过程中，坐标( f ，7 7 ，f ) 是不变的，因此表示了伪粒子 的固有特征。同时，计算网格随着伪粒子的流动而发生变化，不同于拉格朗日 坐标系下网格随流体微粒的变化而变化，因此当选择合适的h 值时，就可以避免 拉格朗日坐标系下严重的网格扭曲。 几点说明： ( a ) 在变换( 1 1 ) 中，h 是坐标( a ，# ，叩，亭) 的任意函数， 的不同选择就 可以获得不同的网格分布，因此可以通过确定h 的值得到优化的网格结构。 ( b ) 在变换( 1 1 ) 中，( 爿，l ，j p ，b ，m ，q ，c ，n ，r ) 是坐标( ， 川， 0 ；在介质2 的区域中，g l 时，特征值玎：的值与h l 时的值符号相反，从而 表明信号的传播方向是错误的；当h o 时，伪粒子( p s e u d o p a r t i c l e s ) 的运动 方向与流体质点的运动方向相反，因此在数值模拟过程中，计算不能继续下去。 同时，在前面的讨论中我们看到，当h = 1 时，二维扩展欧拉方程只是弱双曲 系统，因此在数值计算中，我们将h 的取值限制在o h 0 的解；( b ) 姐果h 是方程的解，则冬也是方程的鳃( c 是任意常数) 。 由于 的取值需限制在o h o 的解，可以做如下的变 量替换； 令：g = i n ( h q ) 则方程( 2 7 ) 变为： s 2 ( a s i n o - b c o s o ) 善+ 丁2 ( m c o s o - l s i n o ) 。o 鍪g ，7 ( 2 8 。) 毋c b 警一爿警m 2c m 等一上等l 其中， 口为质点流动方向角，g = “2 + v 2 ， “= q c o s o ，v = q s i n 0 为了保持网格正交，应满足下面的条件： 箜：塑 叙 础 a 善一 a 叩 印 缸 从而可以得到：4 = m ， b = 一l ，将其代入方程( 2 ，g a ) 得： 中国j ：程物理研究院障士学位论文 2 0 ( m s i n 9 + 细s 9 ) 袅+ ( 4 c o s 口+ b s i n d 要 d ；d 叩 ( 2 8 b ) = 一( 工0 c o s 0 + m m 8 s i n 0 ) 一( a0 c o s o ，+ 甘b s s i n o ) ， 、 e 8 。o q8 q 。 票+ ( s 枷+ 三c o s 曰) 翼+ ( 爿c o s 口+ 胤n 目) 祟 0 c 磊o s 0 屿o s m 。s 叱) - ( a 掣枷霉0s m o i ：一( 三+ m 塑塑+ b 、 一 、8d08 ”8 玎” 将方程( 2 9 ) 离散后，采用迭代方法进行求解，当古( g ”1 一g ”) ( s 时( 占为 当得到解g 的值后，我们也就得到了 的值： = 。么取c 为整个计算空间 中等的最大值，则 = 哆c ，并且满足o l 之间。 2 4 求解统一坐标系下二维扩展欧拉方程的方法 2 4 1 维数分裂近似法 对于多维流体力学r i e m a n n 问题的求解比较常用的种方法是采用维数分 裂近似法来得到它的近似解。也就是将一个多维问题的解分裂为几个一维问题 的顺序解。在对统一坐标系下二维扩展欧拉方程( 2 3 ) 进行求解时，我们使用 了s t r a n g 分裂近似法( 2 7 】) 。 中国j i 程物理研究院博士学位论文 2 1 首先将二维问题( 2 3 ) 分裂成f 方向的一维阔题帮刁方l 各】的一维问题，然后 分别对一维问题求解。 假定三：o 得到孝方向的一维方程组为： d 票+ 婴_ o ， d ，l d 毒 假定昙：o 得到即方向的一维方程组为： d 告 嚣+ 熹= 。 ( 2 1 0 ) ( 2 1 1 ) 令： 时i i i 步长为 ，堪。表示i - 掌平面上一维方程组的解算子，表示 五一玎平面上一维方程组的解算子，e = ( 岛b “，v ，儿a ，b ，l m ) 7 ， 则二维问题( 2 3 ) 的解为： 豆“= 己i 。q 。三二。豆 ( 2 1 2 ) 下面我们将以i - f 平面上一维方程组( 2 1 0 ) 为例，给出求解方法，2 - u 平 面上维方程组( 2 1 1 ) 的解可类似得到。 2 4 2g o d u n o v 间断分解方法 经过维数分裂后，善方向的一维方程组为； 0 e 挪 + = 0 a za f 7 这里， ( 2 13 ) 中国j 一程物理研究院博士学位论文 e =f = p ( i h ) l p ( 1 ) m 十p m p ( 1 一 、，v p l p o h ) l e + p l 二( 卜h ) l y l h u h v 0 o 我们看到，方程组( 2 1 3 ) 是双曲型的守恒律组，因此可以用现有的适用于 双曲型守恒律组的方法来求解。在本文中，我们将采用带有m u s c l ( 2 9 ) 修正 ( 即初始间断两侧的状态通过线性插值的方式得到) 的o o d u n o v 间断分解方法 ( 【2 8 】) 求解方程组( 2 1 3 ) 。 令：k 表示对闯步数，i 表示善方向的网格号，j 表示巧方向的网格号。 = “一为时间步长，毒= 善，一手。和刁，= 叩，一叩。分别表示善方向和卵 j + 一 一 ，十一，一一 方向的步长。在我们的数值计算中六和刁，均取了相同的值，即： 矢= # ，v i ；，7 ，= 7 7 ，w 将计算的物理变量q = ( p ，p ，“，v ) ”和几何变量k = ( 爿，b ，上，m ) 1 定义在网格 中心，即，q 7 和k 。k ，并且具有按体积元平均的意义。 用，表示任意一个变量，则厂的网格平均值为： 工：2 巧l 1 似，觇) 必， 厂的时间平均值为： 争击斯) 撇 刚方程组( 2 1 3 ) 的差分格式为： 班础础础出一川4日l m 中国j ：秽物理研究院博士学位论文 鲋藿，筹c 一，幺川2 ，n 坪 其中，网格边界上的数值通量帚4 由r i e m a n n 问题的精确解来计算，如何得 h i 。7 到r i e m a n n 问题的精确解在后面有详细的讨论。 2 。4 。3t s e 近似方法 当我们对一维方程组( 2 1 3 ) 整体求r i e m a n n 问题精确解时，存在很大的困 难，因此我们采用了t s e 近似方法( t i m es t e p w i s ee u l e r i a n 即p f o x i m 烈i o n ) ( ( 2 3 ) 。 这种方法的基本思想是：当由时刻计算 ( 兄坩“) 时刻的变量值时， 首先假定几何变量k = ( 爿，b ，l ，m ) 7 和h 保持不变，但仍然是孝和印的函数( 叩是 参变量) ，即 k = k ( , t k ，f ，叩) ，h = h ( a k ，孝，7 7 ) 求解物理量守恒律组，得到q = p ，p ，虬v ，力钓值，然后修正几何变量 k = ( a ，b ，厶肘) 的值和h 的值。 2 5 五善平面的r i e m a n n 问题的解 由于经过坐标变换后，统一坐标系下的物理量守恒律组比欧拽坐标系孛的耍 复杂一些，因此求解统一坐标系下的r i e m a n n 问题时，会比较麻烦。下面我们 详细讨论a 掌平面的r i e m a n n 问题的解。 根据前面的讨论，在对r i e m a n n 问题( 2 】3 ) 进行隶解时采用了t s e 近似 方法，因此我们得到下面的一维r i e m a n n 问题； 孚十要地 加0 ，一 掌 0 器；： 为了容易求解方程组( 2 1 5 ) ，我们做以下近似处理，用h l 和h ，的平均值来代 替它们，即；叠= 南，= 三国+ 自，) ，同样， a e = ar = i ( a l + a , ) r b t = b ? = 。 与= t = 圭( 厶+ l ) ，鸩= m ，= 丢( m ，+ 必，) 。 这样，方程组( 2 1 5 ) 中的系数是常数。 为了求解r i e m a n n 问题，现将方程组( 2 15 ) 进行变换： 平面亭= c o n s t a n f 的法线方向为： n - 嵩叫，棚s 其中，v 亭：( 荨，当 中国r ：程物理研究院博士学位论文 2 5 平面f = e o 嬲t a n t 切线方向为： 将流体速度q = 0 ，v ) 7 投影到法线方向n 和切线方向t 得到： j 2q m2 ( 洲一v l s ，( 2 1 6 ) 卜= q t = ( u l + v m ) s 其中，s ：压可 通过变换( 2 。1 6 ) ，方程组( 2 1 5 ) 变为： 其中 e= 并且，s 和h 是常数。 n 口 d f 丛e a 7 1 f p ( 1 一h ) c o s p ( 1 一h ) c o i s + p s p ( 1 一h ) c o r s p ( 1 一h ) e s + c o p s 二_ p ( 1 一h ) c o s ( 2 1 7 ) 从对方程组( 2 1 7 ) 的观察，我们可以看到，方程组( 2 1 7 ) 与欧拉坐标系 下经过维数分裂后得到的维方程组： a 。 新 p 硎 矾 p e p ，一1 a 一 融 删 俐2 + p 翻v p u e + u p p u ，一 = o 具有同样的形式，因此可以用类似的方法求解( 6 2 ，6 5 ，6 6 9 。 下面详细地给出r i e m a n n 问题的解( o s h 1 ) ( 2 1 8 ) 0 o 宇f p p参力力 rll 几 = 曲 啦r l l r 堡苗叩 o p堕觐 中国上程物理研究院博士学位论文 ( a ) 特征域 首先求方程组( 2 17 ) 的特征值和相应的右特征向量。 将方程组( 2 1 7 ) 写为基本变量的形式 a 箜+ b 翌：o 8 a 8 孝 兵甲，u = ( p ，p ，国，f ，) 。 ，1o00o 脚o p o o f00 p 0 小伊，击p 甜p r 一南 l 击。一寿 b = 1 ( 1 - h ) 亡1 o ，一 口 r 2 一h ) p c o , o r ( 3 _ h ) p r 0 2 + ：1p f 2 p y 一1 由d e t ( 匹a b ) = 0 得到方程组( 2 1 9 ) 的特征值为： 旷半国，( 3 重) d ：= 丢【( 1 一 ) 口1 与盯，相对应的右特征向量为： r 。= ( 1 ，0 ，0 ，0 ，0 ) 7 r2 = ( 0 ，o ，0 ，1 ，0 ) 1 r3 = ( o ，0 ，0 ，0 ，1 ) 与仃相对应的右特征向量为： p c 专上去，o ， 0 0 ( 1 - h ) p c o + 吉p ( 1 - h ) p c 0 7 0 ( 2 1 9 ) o o 0 1 一h 一湎丁印 1 一h 一而r p 曲 y 一， o，o山一。 啪矿帅“，埘，mm ”呦 中国= 程物理研究院博士学位论文2 7 显然，盯：特征域是线性退化的，而盯。特征域是真正j 线性的。因此，r i e m a n n 问题( 2 1 7 ) 的解中将包含接触间断、中心稀疏波与激波( 或激波) ，接触间断 位于中间，解的基本结构如图2 1 所示： s l i pl i n e |3 s h o c k、 j j l 、 o 图2 1r i e m a n n 问题解的基本结构 ( b )中心稀疏波 盯：特征域中含中心稀疏波的解可通过求解下面的方程组得到： ( 2 2 0 ) 假定初始状态：q 。= ( p o ，p o ，0 9 0 ，) ，波后状态：q = ( 麒p ，0 3 ，f ，) 7 ，则 其中，a = 旦， 声速日。： p o ( 2 2 1 ) 土妒 ，一矗 + 一 o n i i i l = = 塑咖塑咖生印立勿 一1 = ，2一y 矿 雕 孙 中国j 一程物理研究院博士学位论文 从方程组( 2 2 1 ) 可以看到，f 和y 穿过中心稀疏波是不变的，物理量q 的 值与几何变量k 以及h 的值是无关的，并且解的表达式与一维扩展欧拉方程的解 的表达式是一致的。 设( a ，亭) 是中心稀疏波区内的任意一点，则特征线的斜率为： a 芝芝 翥2 i 2 q 中心稀疏波区内的物理解为： 。 篇千两y 丽- i 卜蛾一铷寿 p = p o a 。， 卯= 峨粤( o ! 2 y 一1 ) ， r = f o ， y = ，o 对于间断解，方程组( 2 1 7 ) 的r a n k i n e h u g o n i o t 条件为： c 陋】_ 【( 1 一h ) p o , s l c 胁】= ( 1 一h ) p c o2 s + p s l c 沁r 】= 【( 1 一h ) m z s l c 胁e 】= 【( 1 一h ) p e o e s + c o p s 】， c 爵 b ”啪出 其中， 】表示括弧内的物理量在间断两侧的状态差，c = 萼是间断速度。 d l ( c ) 激波 用qo = ( 风，p o ，0 3 0 ，t o , ) 7 表示波前状态，q = ( n p ，r ，y ) 7 表示波后状态，则 解的表达式为： 中国丁程物理研究院博十学位论文 2 9 。篇， 。丽a 萧o ( 口- 1 蓊) ， 我们再一次看到，r 和y 穿过激波是不变的，物理量q 的值与几何变量k 以 及h 的值是无关的，并且解的表达式与一维扩展欧拉方程的解的表达式是一致 的。 ( d )接触间断 在这种情况下，间断两边的压力和法向速度是相同的，即： p = p 。， p2 而间断两边的密度、切向速度和比热比y 是任意的，它们的值可以是相同的， 也可以是不同的。同样可以看到，间断解与一维扩展欧拉方程的间断解是一致 的，并且和几何变量k 以及h 的值是无关的。 为了得到r i e m a n n 问题的整体解，我们还需要给出接触间断附近的压力以及 激波速度的计算公式。 对于压力的求解与一维欧拉方程求解接触间断附近压力的方法是一样的( 可 参见文献【6 2 ) 。将接触间断附近的压力和法向速度分别记为尸和，则： 对于左波( 向左激波及向左中心稀疏波) ，有： 国一国f + 兰善：0 ，( 2 1 2 2 ) 0 l 对于右波( 向右激波及向右中心稀疏波) ，有： 一国，一单：0 ，( 2 1 2 3 ) d 其中， 中国j ：程物理研究院博士学位论文3 0 令 口t = p d i pk dk l 一二 j 口t _ t 面 格舻j 雁2 y 。正pk 下， p p k = ，7 f ( p ；n ) ：生盟，t ：f ， 吼 则( 2 2 2 ) 、( 2 2 3 ) 式改写成： 甜一棚，= 一f ( p ；p p ，) ， 一鲫，= f ( p ；p ，所) 上两式相减，消去0 9 ，得到尸的方程： t o f t o ，= f ( p ；p ，p f ) + 厂( ，；p ，p ，) i f ( p ) ( 2 2 4 ) 可以证明，f ( p ) 是p 的单调上升的有连续导数的凸函数，因此可以用牛顿切 线法求方程( 2 2 4 ) 的根，取适当的初值，用下面的迭代计算公式： p 。“：p s 一! ! ! ：2 ：! 竺! 二竺1 2 ， f ( p 。) 迭代收敛后得到p 的值。 下面给出激波速度的计算公式： 统一坐标系中的间断速度c 与欧拉坐标系中的间断速度d 之间的关系为 c = 丢( 心) 则左激波速度： c - 要i ( 1 - 啪，一日 。霈y,+i(p_ll+i 五成 中国二r ：程物理研究院博士学位论文 右激波速度 铲永 l 。y 。r 乃+ i 、p p ，1 ) ，+ 1j yz ，l p ，j 通过上面的讨论，我们已经得到了r i e m a n n 问题( 2 1 7 ) 的解，并且看到， 穿过中心稀疏波、激波或接触间断后，物理量q = ( p ，p ，国，l ，) 7 的值与几何变量 k = ( 爿，b ，l ，m ) 以及h 的值是无关的，但同时我们也应该注意激波速度、接触间 断速度和中心稀疏波区内的流动结构( 如：波头和波尾速度等) 与k = ( a ，b ，l ，m ) 以及h 的值是有关的。此外，从r i e m a n n 问题( 2 。1 7 ) 的精确鳃来看，我们可以 说，统一坐标系下的一维扩展欧拉方程，在采用t s e 近似方法后所得到的弱解 与欧拉坐标系下的一维扩展欧拉方程的弱解是等价的，它们解的形式之间的任 何差别