（计算数学专业论文）运动界面的捕捉与追踪.pdf

摘要在流体力学和工程计算中有这样一类问题，求解区域中物理量有大的间断或求解区域有活动的边界( 包括激波、自由面、物质界面) 如气动力学的激波，流体流动与空气的交界面，结晶、凝固、融化等问题的活动边界，几种不相溶的流体在同一流场中物质界面等等眩些区域不同于一般的流场中物理量的大梯度区，在这些地方物理量往往是绝对间断的；它们也不同于一般流场的边界，它们是随时间运动的，而且两侧的流体会相互作用，如气动力学的激波两侧要满足r a n k i n e h u g o n i o t 间断条件、二相流( 两种不相溶的流体在同一流场中运动) 中两种流体界面两侧的压强和法向速度连续、s t e f a n问题中不同相间的界面上温度是给定的( g i b b s t h o m s o n 关系) 、流体自由面上的压强等关键是它们随时间运动，所以我们统称之为运动界面含有运动界面的问题称为运动界面问题这类问题自数值计算产生以来就引起了广大计算数学和计算力学工作者极大的兴趣和关注对于这类问题的求解方法，也已经有许许多多其中包括早期提出的求解自由面问题的格子类方法( 如p i c 、m a c ) ，后来的质点链接法，v o f 方法，l e v e ls e t 方法等等这些方法对于具体问题有着各自的优势，而且对于每种方法的细节上的处理也只能是具体问题具体解决如何利用现有的数值方法解决好方法的每一个步骤，找到更好的方法组合来求解复杂的运动界面问题，是重要的1 广7 本文较系统地介绍了求解运动界面问题的几种重要的方法，并通过大量的数值试验，分析数值结果，希望找到适合求解运动界面追踪问题的方法组合最后，我们用间断g a l e r k i n 有限元方法求解物理量的控制方程，用l e v e ls e t 方法追踪界面，对各种一维、二维激波问题进行了计算，得到了较理想的结果、，v关键词：运动界面：捕捉和追踪，t v dr u n g e k u t t s l 聃g a l e r k i n 有限元，方法，v o f 方法jl e v e ls e t 方法p 7a b s t r a c ti nf l u i dd y n m n i c sa n de n g i n e e r i n gc o m p u t a t i o n ，t h e r ea r eak i n do fp r o b l e m s ：i nt i l ec o m p u t a t i o n a la r e a ，t h ( 、i ea r eg r e a td i s c o n t i n u i t i e so rm o v i n g1 ) o u n d a r i e s t h o s ei n c l u d es h o c kw a v e ，f r e es n r f a c e ，m a t e r i a li n t e r f a c e t h e s ei n t e r f a c e sa r ed i f f e r e n tf r o mg e n e r a la r e aw i t hl g r e a tg r a d s a tt h e s ei n t e r f a c e s ，t h ep h y s i c a lv a r i a b l e sa r ea b s o l u t e l yd i s c o n t i n u o u s a n dt h e ya r en o tl i k eb o u n d a r i e so fg e n e r a lf l u i dr e g i o n ，t h e ya r em o v i n g a n dt h ef l u i d sn e i g h b o u r i n ga ti n t e r f a c ew i l la c to ne a c ho t h e r s u c ha st h ea i rs i d e b yt h es h o c kw i l ly i e l dt h er a n k i n e i t u g o n i o tj u m p i n gc o n d i t i o n t h et w ok i n d so fi n c o m p r e s s i b l ef l u i dn e a r b yt h e i l 。i n t e r f a c ew i l lh a v ec o n t i n u o u sp r e s s u r ea n dn o r m a lv e l o c i t y t h et e m p e r a t u r eo nt h em o v i n gb o u n d a r yi ns t e f a np r o b l e m sy i e l d sg i b b s t h o m s o nr e l a t i o n t h es a m ek e yc h a r a c t e ro ft h e s ep r o b l e m si st h a tt h ei n t e r f a c e sa r em o v i n g s ow ec a l lt h e ma l lm o v i n gi n t e r f a c e t h ep r o b l e m si n c l u d i n gm o v i n gi n t e r f a c e sa r ec a l l e dm o v i n gi n t e r f a c ep r o b l e m s f r o mt h eb e g i n n i n go fn u m e r i c a lc o m p u t a t i o n ，c o m p u t a t i o n a lm a t h e m a t i c a na n dc o m p u t a t i o n a ld y n a m i c i s ts h o w e dg r e a ti n t e r e s ti nt h i sk i n do fp r o b l e m s ，a n dd e v o t eg r e a t l yi n t oi t t h e r ea r eal o to fm e t h o d sa b o u tt h i sk i n do fp r o b e m s ，i n c l u d i n ge a r l yc e l l t y p em e t h o d sa n dm a r k e r - l i n k l n gm e t h o d ，v o l u m eo ff l u i dm e l ，h o d ，l e v e ls e tm e t h o d ，e t c t h e s em e t h o d sh a v et h e i rr e s p e c t i v ea d v a n t a g ef o rd i f i e r e n tp r o b l e m s m o r e o v e r ，t h ed e t a i l so fe a c hm e t h o d1 l a v em a n yt h i n g st od o i nt h i sp a p e r ，s e v e r a li m p o r t a n tm e t h o d sa r ed e s c r i b e d t h eg o a li st ol o o kf o ra na p p r o p r i a t em e t h o d s c o m b i n a t i o nt os o l v em o v i n gi n t e r f a c ep r o b -l e m st h r o u g hal o to fn u m e r i c a le x p e r i m e n t sa n da n a l y s e s f i n a l l y ，u s i n gd i s c o n t i n u o u sg a l e r k i nf e mt os o v ee u l e re q u a t i o na n du s i n gl e v e ls e tm e t b o dt ot r a c km o v i n gi n t e r f a c e ，w es o l v e dt h e1 - da n d2 - ds h o c kw a v ep r o b l e m s ，a n dg o tg o o dr e s u l t s k e y w o r d s ：m o v i n gi n t e r f a c ec a p t u r i n ga n dt ia c k i n gt v dr u n g e - k u t t ad i s c o n t i n u 0 1 1 sg a l e r l d nf e mv o fm e t h o dl e v e ls e tm e t h o d致谢本文是在刘儒勋教授的悉心指导和直接参与下完成的，在这三年里，他不仅在学业上对我严格要求，勤加督促，而且在生活上给予我无微不至的关怀刘老师渊博的知识，严谨的治学和谦逊的为人都深深地影响着我师母隋老师，在生活上更是给予我热情的关怀和细致入微的照顾，在此向他们致以最诚挚的谢意衷心感谢“长江计划”讲座教授，美国布朗大学应用数学系主任舒其望博士对我的热情指导，他在科大开的“e n o 和w e n o ”、“d i s c o n t i n u o u sg a l e r k i nf e m ”等课程，深入浅出、极易接受，是我听到的最好的讲座之一冯玉瑜、陈祖墀、陈发来等老师开的基础课使我受益非浅，得以顺利完成学业在此一并深表谢意我还要感谢师姐李宏、师兄汪继文博士以及张瑞、张磊、融华等同学给予我学习上的大力帮助感谢数学系科学计算与计算机图形学实验室全体成员，本文的大量工作都是在该实验室完成并得到他们的极大帮助另外，我在科大数学系七年一直得到了各位老师和同学的关心、帮助在此谨向他( 她) 们表示深深的谢意。我更要感谢我的父母，作为沂蒙山的农民供我读书近二十年，他们比别的父母付出了更多的艰辛，至今他们仍然在辛苦的劳作着希望我能以成功带给他们欣慰和喜悦第一章运动界面问题概述1 1 问题的提出在流体力学和工程计算中有这样一类问题，求解区域中物理量有大的间断或求解区域有活动的边界，包括激波、自由面、物质界面如气动力学的激波，流体流动与空气的交界面，结晶、凝固、融化等问题的活动边界，几种不相溶的流体在同一流场中交界面等等这些区域不同于一般的流场中物理量的大梯度区，在这些地方物理量往往是绝对间断的；它们也不同于般流场的边界，它们是随时间运动的，而且两侧的流体会相互作用，如气动力学的激波两侧要满足r a n k i n e h u g o n i o t 间断条件、二相流( 两种不相溶的流体在同一流场中运动) 中两种流体界面两侧的压强和法向速度连续、s t e f a n问题中不同相间的界面上温度是给定的( g i b b s t h o m s o n 关系) 、流体自由面上的压强等关键是它们随时间运动，所以我们统称之为运动界面含有运动界面的问题称为运动界面问题对于整个流场的求解来说，运动界面是重要的，而且往往是关键的不仅是因为运动界面往往对整个流场中流体的运动情况以及物理量影响很大；而且有时候我们最关心的就是运动界面的位置及其附近的流场由于某些物理量在界面上是间断的，在求解时我们希望得到高分辨率的解( 如果解本身是间断的就希望得到间断的解) ，并且希望计算出在任意时刻的界面的位置、法向、曲率等几何性质，以便在求解物理量的控制方程时，在界面附近做特殊处理，以提高分辨率或者要求界面上满足某种边界条件( 如s t e f a n 问题) 而一般的数值方法，包括高分辨率方法，都是将运动界面作为大梯度区进行计算由于数值耗散，模糊了间断，得不到界面的确切位置和形状所以要引进界面追踪，通过某种方法得到任意时刻的界面位置和形状，然后再根据界面来求解整个流场区域上的物理量这样问题就归结为：1 追踪运动界面由于界面随时问运动，界面上的点或者满足某个控制方程，或者以某种速度和方式运动我们希望追踪它任意时刻的位置和形状；2 求解整个流场上物理量的控制方程我们可以用已有的方法( 如高分辨率方法) 求解流场，但要在界面上或界面附近做特殊处理比如让界面2 0 0 1 年中国科学技术大学硕士学位论文第2 页第一章运动界面问题概述5 1 2 各种追踪方法概述上满足某种条件，在界面附近做偏心差分或插值计算，以提高分辨率因为实际上，在间断处产生的耗散主要是由于在计算界面一侧的流场时，用到了另一侧的流场的物理量函数值所以我们在已经追踪到的界面附近做精细地插值、差分处理，希望能得到绝对的间断这也就是如何把界面追踪与已有数值方法耦合起来求解整个流场的问题由上所述，有效的追踪运动界面( 包括激波、物质界面、活动边界、自由面) 是问题的关键而且这类问题在工程计算和流体力学中极为常见，十分重要，因此从数值计算产生以来，计算数学和计算力学家对此投入了极大的关注，提出了各种各样的追踪方法1 2 各种追踪方法概述最早出现的有效处理自由面问题的方法是h a r l o w 和w e l c h 等人提出的格子类( c e l l t y p e ) 方法 1 这类方法在e n l e r 网格的基础上对控制方程进行差分离散，在网格中布置若干质点来标记流体，通过这些标记点( m a r k e rp a r t i c l e s ) 来模拟自由面的大体位置和形状网格固定不动，是e u l e r 型的，而标记点随流体运动，是l a g r a n g e 型的格子类( c e l l t y p e ) 方法中著名的有p i c 、m a c 、f l i c 等其中p i c ( p a r t i c l e si nc e l l ) 方法是e v a n s 和h a r l o w 于1 9 5 7 年提出的标记点是有质量、动量、能量的流体质点，先在e u l e r 网格上对n s 方程进行无输运计算，然后通过面积加权平均，估计出标记点的速度，再移动标记点，并对物理量进行输运修正此后，g e n t r y ，m a r t i n 和d a l y 对p i c 方法进行简化和创新，提出了f l i c ( f l u i di nc e l l ) 方法p i c 方法需要存储每个质点的物理量，存储量很大，而1 9 6 5 年h a r l o w和w e l c h 提出的m a c ( m a r k e ri nc e l l ) 方法采用无质量、动量、能量，只有坐标位置的标记点，来表示流体所有物理量均归属网格，方程单独进行e u l e r 型差分离散求解，然后标记点做l a g r a n g e 型移动这类方法追踪运动界面的精度取决于网格内质点的多少，要得到较高的精度，必须在每个网格内布置大量的点这就要求很大的存储量，同时我们要追踪这些质点，时间复杂度也相当大在此基础上，提出了质点链接法，就是只在运动界面上布置质点，在求解过程中，这些质点以当地速度( 界面的运动速度) 运动在任意时刻，我们只要求得所有质点的位置，然后将它们顺序链接就得到该时刻的界面了这种方法在运动界面不是很复杂的时候，大2 0 0 1 年中国科学技术大学硕士学位论文第3 页第一章运动界面问题概述5 13 格子类方法大降低了存储量而且可以给出显式的界面但是，这种方法难于处理复杂界面和运动界面在运动中发生拓扑形变( 如翻转、破碎、相交等) 的情况，而且难于向三维推广七十年代，h i r t 4 1 等人提出的高度函数法只能处理运动界面是某个方向( 三维时是某两个方向) 的单值函数时的情形之后，他们又提出了v o f ( v o l u m eo ff l u i d ) 方法【4 在这种方法中，定义一个所谓v o f 函数e ( 叠，t ) 把自由面问题中的流体或多相流场中的一种流体称为“目标流体”，每个网格上的v o f 函数g 定义为该网格内目标流体的体积与网格体积( 二维时是面积，我们统称为体积) 的比值这样，含有运动界面的网格就是满足0 c 1 的网格，而在充满目标流体的网格上c = l ，在不含目标流体的网格上c = 0 v o f 函数满足v o f 方程筹+ “鬻+ 筹= 0 其中( “，u ) 是流体的速度v o f 方法就是通过求解v o f 方程追踪运动界面的v o f 方法基本上解决了上述困难，易于向三维推广，易于处理复杂界茴的情形但是，v o f 函数是一个有明确物理意义的量，在求解过程中要保持它的这一性质是困难的，而且一般只能求解不可压流体，且每种流体的密度是常数的情形九十年代左右，o s h e r 和s e t h i a n 提出了l e v e ls e t 方法1 7 ，1 8 ，1 9 1 这种方法也是定义一个函数妒( 量，t ) ( l e v e ls e t 函数) 来追踪界面，但是l e v e ls e t 函数只与运动界面有关，它的零等值面( 等值线)就是运动界面比如我们可以定义妒( 毒，t ) 是到t 时刻运动界面的符号距离在求解过程中，我们只要保持妒和界面的这种关系关系后来o s h e r 等人又提出了重新初始化的概念，使算法在求解过程中保持妒( 童，t ) 是童到t时刻运动界面的符号距离这种方法在追踪二相流运动界面、结晶问题的活动边界、激波等问题取得了很好的结果本文的目的是通过大量的数值试验，理论、结果分析，希望找到适合求解运动界面问题的方法组合我们用间断g m e r k i n 方法求解物理量控制方程，用l e v e ls e t 方法追踪界面，得到了较理想的结果1 3 格子类方法本节较详细地介绍1 9 5 7 年由e v a l l 和l l a r l o w 、提出的p i c 方法和1 9 6 52 0 0 1 年中国科学技术大学硕士学位论文第4 页第一章运动界面问题概述5 i3 格子类方法年由h a r l o w 和w e l c h 提出的m a c 方法一、p i c ( p a r t i c l ei nc e l l ) 方法p i c 方法是格子类( c e l l t y p e ) 方法中最重要的一种方法之一最早是由l o sa l a m o s 实验室的e v a n s 和h a r l o w 于1 9 5 7 1 2 1 年提出的，而后经多人作出推广和完善p i c 方法把一个连续的流场看成有限个质量集中的流体质点系，它们是分布在e u l e r 网格内的，具有l a g r a n g e 特征的离散点通过对它们的计算、追踪和调整，实现流场的数值模拟和显示因而它既可以用于可压缩流，又可以用于不可压流，以及多相、多介质流p i c 方法采用的质点，是具有质量、动量和能量的网格内流场的质量、动量和能量等特性完全是通过对该网格内的流体质点的统计累加给出的所采用的差分网格是e u l e r 型的，它不随流体运动；而网格内的质点则是流动的，在一个时间段内，有的流进格子，有的流出格子，是l a g r a n g e 型的p i c 方法总过程分为三步：第一步为e u l e r 步，不考虑格子之间质点的输运，对网格的速度进行预估，具有一般差分法的特性；第二步为l a g r a n g e 步，考虑质点的运动和网格间的质点输运，根据第一步估计出的网格速度，计算网格内质点的速度，并以此速度计算质点的新的位置；第三步对网格上各物理量进行修正，通过对网格内质点的统计累加计算出网格的密度、速度、能量、压力等具体地讲，用p i c 方法求解自由面问题的步骤如下：l 、初始状态质点的布置和物理量的赋值根据问题的实际力学背景和初条件，在整个流场中布置m 个质点r ( p =1 ，2 ，m ) ，在每个网格内质点的数目由该网格内流体量的多少和要求达到的精度决定的并对每个质点赋以确定的质量、动量、能量和初始位置然后，每个网格内的密度、动量、能量可以由网格内所含质点的统计累加得到具体地讲( 以均匀矩形网格上求解二维e u l e r 方程为例) ：p i c 方法中的质点具有位置( 。，) 、质量m 、速度( “，”) ，可表示为尸。x 。y 。m 。“。，) ；网格具有密度p 、动量( m 。，m 叼、能量e ，可表示为c i j ( 阳，m 舀，m 嚣，e 。j ) 他们之间满足如下关系：2 0 0 1 年中国科学技术大学硕士学位论文第5 页第一章这动界面l u q , n 概述5 l3 格子类方法记6 i j = 6 ( 斟一i ) a ( 岗一j )其中，a ( o ) = 1 ；j ( 。) = o ，x 0 ， x 1 表示不大于x 的最大整数定义舰，为网格( i ，j ) 内质点的质量总和，m 0 ，”z 嚣分别为网格( i ，j ) 内质点的z 方向和方向动量总和，则有而网格密度为网格内质点质量总和除以网格体积，网格速度为网格内质点动量总和除以质量总和舶= 舰j 啦j = m i m i jv i i = r n ；j m o( 1 3 2 )2 、计算的第一个阶段一一不计输运的差分计算假设z = 。时刻的状态已知露( z ：，：，m ：，“：，”：) c 焉( p ；，m ，n ，x ，m ，e 3 )现在不考虑质点的输运来对网格的速度、能量进行预估由于p i c 方法采用带有质量的点，我们不必考虑质量守恒，只考虑加入人工粘性的动量和能量方程，而使控制方程化为誊黼p 替耻0。p 酱+ ( + q ) ( 券+ 器) =其中q 是引入的v o nn e u m a n n 人工粘性项。一j ( z m 爱) 是 0采用f 1 1 c s 格式计算矗件l 2 ，j2 “1 1 2 j 一面；等尝面( p 苒一p 乙+ q 耳一q 乙)祝一22v ? , j - k l 2 一疆篇万( 吃+ l p 乙+ 联n ，抖1 一吼)j = e 0 一雨1t 硒a t p 。n + q 4 1 2 ，j ) 面i + 1 2 ，一( 踢+ q 卫1 2 ，) 面i 一1 2 ，j +石i 岛 ( 踢+ q ；n ，+ - 2 ) 晚，j + 1 e 一( p ”n + q w n 一- 2 ) 矾，j t 。 一面”基 略1 2 ，j q 2 l 2 ， 一i ”五a 百t q 己+ l 2 一q 。g 一1 2 】)( 1 3 4 )30pm豇m 脚=，um“盯4m 脚j jr um“m矗m 一=慨2 0 0 1 年中国科学技术大学硕士学位论文第6 页第一章运动界面问题概述5 l3 格子类方法其中祝u = j ( “0 + i ”) ，o i j = ；( 口马+ i u ) ，人工粘性项的离散这种半隐式的处理可以提高稳定性和精度3 、追踪质点已知质点r 在t = f 。时刻的位置，根据上一步得到的速度计算t = t 。+ t的质点位置分为以下几步：( i ) 决定质点r 所在的虚框网格( 如图1 所示) 的左下角下标m = 怎一纠1z = 若一纠1( 2 ) 利用面积加权平均计算质点巳的速度“- - 1 。1 = w 。( u ：+ i 。) 2 - - 1 2 = ”( 蠊+ ) 2( 1 3 5 )其中w l = g ( h + 1 2 ) ，w 2 = ( h + i 2 ) ( 1 一g )w 3 = ( 1 9 ) ( 1 2 一 ) ， 4 = g hg = ( x k + l 2 一) a xh = ( y t + l 2 一 a ) a yu l = u k 一1 1 2 ，i 札2 = i l k + i ，2 ，l1 1 3 = u k + 1 1 2 ，f + 11 1 4 = u 一1 ，2 ，2 + 1w ( m = 1 ，2 ，3 ，4 ) 和 。( m = 1 ，2 ，3 ，4 ) 做类似定义( 3 ) 计算t = i n + 时刻质点b 的位置z ：+ 1 = z ：+ z - z 。7 l z 鳐+ 1 = 蝠+ i ：4 、对网格内质点进行统计累加，计算网格上物理量的值= 。+ 。时刻的网格函数值由网格内质点的累加得到p 矿1 = m 矿1 u 矿1 = m 矿“m 矿1”矿1 = m 寸岫j j l 锈+ 1m科12p 矿1 i 玎+ 量1 n “+ 1 ( 。“2 + i 触”( 1 舢)r 1 = e 。+ 1 p 矿1 一j 1 ) 2 + ( t ，矿1 ) 2 ) 坨1 = ，( p 寸1 ，e 矿1 )n u，n 冲m2，p+，n 件p，一2u，i，、【| l胆诉2 0 0 1 年中国科学技术大学硕士学位论文第7 页第一章运动界面问题概述5 l3 格子类方法实际计算中，对于质点的累加，我们并不采用上述公式里表示的方法，而是先把网格上的物理量( 包括总质量、总动量等) 清零，再对质点p 。( =1 ，2 ，m ) 进行扫描，对每个质点，判断它所在的网格，然后把它累加到该网格的物理量上这样，时间复杂度会大大降低该方法的最初的网格内的质点的布置和赋值，质点的追踪，以及由此对于网格内的物理量的计算和各种异常情况( 例如，质点的过分凝聚、稀疏，边界质点、等等) 的处理是非常有特色除了经典的p i c 方法外，为了适应多介质流、稀疏气体流等等，又有g i l a 方法和s p i c 法等等，以及由此发展起来的具有示踪点与质点的其它格子类方法，如h e l p 、h e m p ，用于计算大变形问题西五j r 。“k + 1 门，1 “i1 l 量缸v t , , , t + l e ：h 兮”k + l ，i + i 2il 一1：p p j“k - 1 n 1lm k + l 2 ，1- v k l 一1 nv k +1 1 q 尼( 1 - 1 ，j + 1 )( i ，j + 1 )叫呜矗产夏( r 1 ，】)( 1 ，j )图i 利用网格速度的面积加权平均计算质点速度图2 面积加权近似m a r k e r 点的速度二、m a c ( m a r k e rc e l l ) 方法m a c 方法是由h a r l o w 和w e l c h 等人在1 9 6 5 年【3 发展起来的它是一种处理不可压缩、粘性流体、自由面流动问题有效的方法它是h a r l o w等人早期发展的p i c 方法的改进它把p i c 方法中带有质量、动量和能量的流体质点改变为只具有坐标位置的标记点，这些标记点只起着标记流体的作用，所有物理量均归属于网格，在固定的e u l e r 网格上求解而标记点以其所在的当地速度运动，自由面形状可以通过标记点的位置来确定这样，m a c 方法比p i c 方法所需的内存大大减少，更具有实用价值同p i c 方法一样，m a c 方法是一种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 型的运动从而能反映每一时刻自2 0 0 1 年中国科学技术大学硕士学位论文第8 页第一章。运动界面问题概述5 1 3 格子类方法由面的变形和运动由于计算网格是固定的，因此不存在通常纯l a g r a n g e方法会出现的网格翻转和畸变等问题此外它又克服了通常e u l e r 方法中难以了解流动细节的缺点这样，对于自由面的大变形、波浪翻卷乃至破碎等现象，m a c 方法都能较好地描述m a c 方法的重大创造在于采用了“示踪点”活动边界、自由面、以至于运动界面的追踪等等问题，都可以通过类似的方式进行特别是今天许多的运动界面重构方法的提出和应用论文，常常可以看到其中利用m a r k e r 点的思想这里，我们简单地就m a c 方法中“m a r k e r ”点的追踪技术作一介绍以二维不可压粘性流自由面问题为例，采用矩形交错网格类似于p i c方法，首先在i 两格内布置若干标记点r ( p = 1 ，2 ，m ) ，这些标记点与流体质点类似，只是没有动量、能量等，只有位置属性，可表示为巳( z 。，钆) 假设在。时刻标记点巳( z ：，鳐) 的位置已知，网格c i i 的速度( “，”) 已知示踪点r 的运动是随流而动的，它在下一个时间步，即t 。+ 时刻的新位置，应当由微分方程决定即掣= “( z ，y ，z ) 掣= ”( z m )( 1 3 7 )z ：+ 1 = t ：+ 1 + ! + 1 “( z ( ) ，( t ) ，t ) d tf 1 3 8 1f 矿1 = ：“+ 口”( z ( f ) ，( ) ，) m、给出可见，计算的关键是速度场矿= ( “( z ) ，”( t ) ) 在时间段【t 。，t 。+ 。】内，在该示踪点附近的变化情况根据数值积分的梯形公式，我们自然可以采用近似公式( 1 3 9 )进行计算其中“：，”：分别表示质点巳在t = t 。时刻所在位置的z ，方向上的流体速度先求解物理量的控制方程( n s 方程) ，得到t = t n + 时刻的网格上的速度m a c 方法也采用了同p i c 方法的“面积加权平均”的方案，来计算质点所在位置的流体速度由于采用交错网格，速度“，”分别定义在网格左右边界和上下边界如图2 所示的情形，对于x 方向的速度i t 。”咿+ 十吣嘧mo出一o+ +蛇| |一一十+鲅2 0 0 1 年中国科学技术大学硕士学位论文第9 页第一章运动界面问题概述5 l3 格子类方法有加权公式q 。= 坐业止坐等学( 1 3 1 0 )其中a k ( k = l ，2 ，3 ，4 ) 如图2 所示，是标记点巴所在的虚线网格的剖分在任意时刻，计算出了每个标记点位置，也就是全部标记点的分布情况，我们就可以把网格分为三类：一类是空网格，其中没有标记点；一类是表面网格，其中有至少一个标记点，且与空阿格相邻；一类是含有标记点且不与任何空网格相邻从表面网格的分布，我们可以看出自由面的大体形状对于表面网格的物理量的计算( 边界条件的处理) ，请参见文献11 第二章间断g e l e r k i n 有限元方法8 0 年代末期，c o c k b u r n ，s h u 和他们的合作者提出了解决非线性双曲方程的间断g a k l e r k i n 方法，它自然地融入了高分辨率有限差分方法和有限体积方法中如数值流通量和r i e m a n n 解、高阶t v dr u n g e k u t t a 时间离散和间断解问题的非线性限制器等思想最初的间断g m e r k i n 方法是为解决双曲问题提出的而近年来，这种方法已经成功地推广到求解其它很多类型的偏微分方程，诸如抛物和椭圆问题19 9 8 年，首届间断g a l e r k i n 方法的国际会议在美国罗德艾兰州n e w p o r t 召开，会议论文集由s p r i n g e r 出版社出版 1 3 】，其内容包括算法的发展、分析与应用，以及该方法的并行计算的实现，同时会议组织者c o c k b u r n ，k a r n i a d a k i s 和s h u 给出了一篇综述性文章这本论文集可以作为详细了解该方法的一般性参考书目本章我们以守恒型方程( 如e u l e r 方程) 为例具体介绍这一方法2 1 一维情形一维守恒律方程及初值条件( 。，t ) n ，b 【o ，t r 。、。1j空间离散，将区间【a ，b 剖分为a = z o 。l x n = b ，网格定义为= 【q 一1 2 ，x j + l 2 】取有限元空间坛= v 三v l 1 ( o ，1 ) ： i l p k ( 厶) ，i = 1 ，2 ，)其中，p ( ，) 表示区间，上的次数不超过e 的多项式空间取任意的光滑函数乘以方程两边，并且在上积分，然后分步积分得到e 掣”( z ) 如一 m ( 。，f ) ) 掣如+ f ( u ( x j + l 2 ，z ) ) ”( 。再。) 一f o t ( 。川。，t ) ) ”( 吐_ l 2 ) = 0( 2 1 2 )et c ( 。，o ) ”( 。) 如= “。( z ) ”( z ) 如1 02 0 0 1 年中国科学技术大学硕士学位论文第l l 页第二章间断g e l e r k i n 有限元方法5 21 一维情形现在以v hx ) u 代替v ( x ) ，以数值流通量代替网格边界点上的函数值，得到离散形式l 毪笋v h ( 。) 如一e 巾c n ( z ，f ) ) 掣如+ 扎z ，1 ，2l t 。4 - ，。v h ( z 再。) 一凡c 。- u ，+ 1 t h ( x ，9 ) ( 。) 出= 正，l $ 0 ( z ) ( z ) 如其中，数值流通量，( u 川- u ，+ + 1 1 2 ) 要满足l i p s c h i t z 连续之外，还要求相容性即，( u ，“) = ，( “)单调性即对第一个变量非减，对第二个变量非增非常著名的满足上述条件的数值流通量有，g 。4 ( t c 一，“+ ) = m i nf ( u ) , i i f ，“u - 一 ooe玎叮出2 0 0 1 年中国科学技术大学硕士学位论文第1 4 页第二章间断g e l e r k i n 有限元方法5 2 2 二维情形2 2 二维情形二维守恒律方程1t c f + ，( “) 。+ 9 ( “) 。= 0( 2 2 1 )设计算区域是【a ，b 【c ，棚，采用矩形网格剖分n = x o x l x n 。=6 ，c = y o y l e ，e 是某个小量，则置“。= 乱。，t ”= 矗”，t l 。= 0 ，“。”= 0 ，“”。= 02 0 0 1 年中国科学技术大学硕士学位论文第1 7 页第二章间断g e l e r k i n 有限元方法5 2 3 时间的离散2 如果是方程组，问题要复杂些 1 4 】我们采用一种简化的方法对u 。，z l y进行限制，t ，u y 一旦真正被限制，则高阶项的系数全部置为零首先我们要用左特征向量，把向量投影到它的特征空间，进行限制，然后再返回来设a = 磊o f ，其特征值和左特征向量分别是sa =则有如下关系s a s = a通过”= s “，把+ “= “0 + - ，一t 。，一“= “屯一“巴，u 0 变换到其特征空间+ ”，一”，”嚣对蝽，v 。k 进行限制 嚣= m i n r n o d ( v i ，+ ，一 )然后把它返回到原空间中，即矗。= s 一1 ”。对u y 做同样的限制，如果( 一t c 0 ) 2 + ( 略一“嚣) 2 e ，e 是某个小量，则置。= 也。，u ”= 矗。，“。= 0 ，“。”= 0 ，“”。= 02 3 时间的离散时间的离散可以采用二阶或三阶的t v d r u n g e k u t t a 多步方法方法采用显式格式，空间离散全部在t 。时间层上进行设空间离散已经完成，记右端项为l n ( “) 常微分方程鲁= l j 。( u ) ，具体离散步骤如下：二阶t v d r u n g e k u t t a 方法t0o ，叶一，也“姒k 州托l，u心、j0、二卜“f，+“广“0儿+扩+瓠。矿= 矿矿2 0 0 1 年中国科学技术大学硕士学位论文第1 8 页第二章。间断g e l e r k i n 有限元方法2 4 数值算例三阶t v d r u n g e - k u t t a 方法“1 = “”+ a t l h ( u “1“2 = ；“+ ；( “1 + a t l h ( u 1 ) )( 2 3 2 )i t 1 = ；“”+ ；( “2 + a t l h ( u 2 ) )本文的算例全部用三阶的t v d r u n g e k u t t a 方法离散时间导数2 4 数值算例算例1 、m a c h3 前台阶问题 1 3 ，1 5 m a c h 数为3 的流动流过一个含有台阶的风洞管道，管道长3 宽1 ，台阶在距管道左端0 6 处，高为0 2 初始时刻管道内充满理想气体，7 = 1 4 ，密度p = 1 4 。压力p = 1 0 ，速度i t = 3 0 ，”= 0 左、右边界分别用入流和出流条件，上、下边界( 包括台阶) 用反射边界条件在2 4 0 8 0的网格上，用三阶问断g a l e r k i n 方法( 多项式p 2 ，时间离散用三阶t v dr u n g e - k u t t a 方法) 计算到时间t = 1 0 ，t = 2 0 ，t = 3 0 ，t = 4 0 ，密度等值线如图1 所示图中画出3 0 条等值线算例2 、双m a c h 反射问题1 3 ，1 5 1m a c h 数为1 0 的激波，以6 0 。角在墙上反射，未扰流体的密度p = 1 4 ，压力p = 1 0 计算区域 0 ，4 x 0 ，1 】，反射墙位于区域底边1 6 。 4 左右边界分别以入流出流条件，在底边上0 z 1 6 的部分赋以准确的波后条件，其余部分是反射墙，上边界为m a c h 为1 0 的激波精确解，即以激波速度移动，波头左边赋以波后条件，右边赋以波前条件在3 2 0 8 0 的网格上，计算到t = 0 2 ，图2 中只画出了区域【0 ，3 0 ，l 】上的3 0 条密度等值线算例3 、稳态斜激波 1 6 如图3 所示的区域f 0 ，4 】 0 ，1 】，其中充满理想气体，y = 1 4 ，三个子区域a ，b ，c 中，初值分布如下：。p = 1 0 ，i t = 2 9 ， = 0 ，p = 1 o 1 4 ，在区域ap = 1 7 ，“= 2 6 1 9 3 ， = 一0 5 0 6 3 2 ，p = 1 5 2 8 2 ，在区域bp = 2 6 8 7 2 ，i t = 2 4 0 1 5 ，v = 0 ，p = 2 9 3 4 0 ，在区域c2 0 0 1 年中国科学技术大学硕士学位论文第1 9 页第二章间断g e l e r k i n 有限元方法5 2 4 数值算例左边界和上边界用固定边界条件，分别设置为上述区域a 和区域b 的初值底边为反射条件，右边为出流条件采用1 2 0 3 0 的网格，稳定状态解的压力等值线如图4 所示在y = 0 5 处的压力系数c p = 2 ( 老一1 ) ( 7 m 。2 )( m o o = 2 9 ，p 。= 1 0 1 4 ) 如图5 所示图lm a c h3 前台阶问题，问断g a l e r k i n 方法求解至t = l0 ，2 0 ，3 0 ，4 02 0 0 1 年中国科学技术大学硕士学位论文第二章问断g e l e r k i n 有限元方法第2 0 页5 2 4 数值算例图2 双m a c h 反射问题，计算到时间t = 0 2 的密度等值线图3 稳态斜激波问题的计算区域和物理量分布图4 斜激波问题稳定解压力等值线圈5 斜激波问题稳定解在”= 0 5 的压力系数第三章v o f 方法3 1v o f 方法概述在解决自由面问题或不可压二相流问题中，格子类方法( 如m a c 、p i c等) 1 】模拟界面比较粗略，由于所用的网格内的质点1 数目有限，容易出现虚假的空网格和错误的密度而且每个网格内都布置若干个质点，要求的存储量十分巨大在此基础上发展的一种“质点链结方法”是，只在流体界面上布置质点，这些质点随时问以当地速度运动，在任意时刻，将这些点顺序链结( 比如样条插值) 就构成运动界面的近似这种方法降低了存储量，但质点的分布会随时间变化，特别是速度很不均匀的流体，需要定时在质点过于密集的地方删除一些，而在稀疏的地方插入一些，这会增加时间复杂度，而且带来较大的误差更重要的是这种方法难于向三维推广，难于处理界面的翻卷和破碎等拓扑上的逻辑困难对于自由面问题，n i c h o l s 和h h t ( 1 9 7 3 ) 4 1 年提出了一种所谓“高度函数法”，把运动界面和某个参考面的距离定义为高度函数h ( x ，t ) 在自由面问题中，高度函数满足方程警蜘瓦o h = ”( 3 1 1 )丽十“瓦2 ”b 1 1j推广为三维情形，高度函数为h = h ( x ，y ，f ) 这种方法降低了存储量，只需要一个一维数组( 三维时是二维数组) 来存放高度函数，方程( 3 1 1 ) 也比较容易求解。但是只能处理界面是。或y 的单值函数的情况七十年代末，h i r t 和n i c h o l s 等人提出了v o f ( v o l u m eo ff l u i d ) 方法4 1 v o f 方法是在整个流场中定义一个函数c ，在每个网格中，这个函数定义为一种流体( 我们称之为目标流体) 的体积与网格体积的比值不包含这种流体的网格称为空网格( c = 0 ) ，充满这种流体的网格称为满网格( c = 1 ) ，包含界面的网格称为半网格( 0 c 1 ) 在任意时刻，知道了这个函数在每个网格上的值，也就可以通过某种途径构造出运动界面然后求解物理方程时可以在界面附近做11 特殊的精细处理，以提高分辨率和精度设计算区域是q ，流体a 所在的区域记为n 1 ，而流体b 所在的区域2 12 0 0 1 年中国科学技术大学硕士学位论文第2 2 页第三章v o f 方法5 3 1v o f 方法概述记为f 2 2 首先定义这样一个函数乖，归f 濮未由于两流体是不相溶的，a 满足如下方程等= 等+ 加a = o现在可以在每个网格岛( 以矩形网格为例) 上定义函数c 0 为n 在上的积分文心，t 、我们称之为v o f 函数同样c 满足方程瓦o c + 矗vc r = o( 3 1 2 )f 3 1 3 1网格，( 3 1 4 )( 3 1 5 )这个方程称为v o f 方程可以看出c i j 其实是网格 ，中流体a 所占的体积比其中充满流体a 的网格上v o f 函数c = 1 ，称为满网格；不含流体a 的网格上的v o f 函数c = 0 ，称为空网格；而包含流体界面的网格上口( 0 ，1 ) 反之，只要求得了v o f 函数在每个网格上的值，就可以确定哪些网格充满a 流体，哪些网格中充满b 流体，哪些网格中含有界面( 即含a 流体又含b 流体) v o f 方法就是通过求解v o f 函数，实现对运动界面的追踪由于v o f函数是有物理意义的，它是网格中一种流体所占整个网格的体积比，我们虽然知道它满足的方程，但在求解方程( 3i 5 ) 时，要保持它的这一物理意义是困难的我们的求解格式首先应该是守恒的，而且每一时间步得到的函数值跨过界面从1 变到0 ，