（计算数学专业论文）数值模拟多介质可压缩流的rkdg有限元方法.pdf

摘要 本文利用r k d g 有限元方法、l e v e ls e t 方法和g h o s tf l u i d 方法数值模拟 多介质可压缩流，并且采用r k d g 有限元方法实现运动界面的追踪，对二维流 体中常见的常数流场、旋转流场和剪切流场做了追踪模拟。为了使g h o s tf l u i d 方法能够用于多介质高密度比流体力学问题的数值模拟，适当的改进了一维的 g h o s tf l u i d 方法，给出了新的带“i s e n t r o p i c ”修正的g h o s tf l u i d 方法。对一维和 二维的气气和气一液两相流进行了数值计算，得到了较高分辨率的计算结果。 关键词：r k d g 有限元方法；多介质可压缩流；g h o s tf l u i d 方法；改进的带 “i s e n t r o p i c ”修正；l e v e ls e t 方法 a b s t r a c t 。i h er k d gf i n i t ee l e m e n tm e t h o d ，l e v e ls e tm e m o da n dg h o s tf l u i dm e t h o d a r ea p p l i e dt os i m u l a t ec o m p r e s s i b l em u l t i m e d i af l u i di nt h i sp a p e r i ti sr e a l i z e db y r k d gf i n i t ee l e m e n tm e t h o dt ot r a c et h em o v i n gi n t e r f a c e ，t h et w od i m e n s i o n a l c o n s t a n tf l o w , r o t a t i n gf l o wa n ds h e a rf l o wh a v eb e e nt r a c e d t h eo n e d i m e n s i o n g h o s tf l u i dm e t h o di sm o d i f i e d an e wg h o s tf l u i dm e t h o dw i t hi s e n t r o p i cf i xi s g i v e ns ot h a ti t c a r ls i m u l a t et h eh y d r o d y n a m i c sp r o b l e mw i t hm u l t i m e d i ah i g h d e n s i t yr a t i o s o m en u m e r i c a lt e s t ss h o wt h eh i g hr e s o l u t i o nr a t i or e s u l t so no n e d i m e n s i o n a la n dt w od i m e n s i o n a lc o m p r e s s i b l et w o f l u i df l o w ss u c ha sa i r a i r ， a i r - l i q u i d k e yw o r d s ：r k d gf i n i t ee l e m e n tm e t h o d ；c o m p r e s s i b l em u l t i m e d i af l u i d ；g h o s t f l u i dm e t h o d ；i s e n t r o p i cf i x ；l e v e ls e tm e t h o d 独创性声明 本人声明所呈交的学位论文是本人在导师指导下进行的研究工作及取得的 研究成果。据我所知，除了文中特别加以标注和致谢的地方外，论文中不包含其 他人已经发表或撰写过的研究成果，也不包含为获得中国工程物理研究院或其他 教育机构的学位或证书使用过的材料。与我一同工作的同志对本研究所做的任何 贡献均已在论文中作了明确的说明并表示谢意。 学位论文作者签名：陈荣三 签字目期：l 0 0 5 年5 2 目 学位论文版权使用授权书 本学位论文作者完全了解并接受中国工程物理研究院研究生部有关保存、使 用学位论文的规定，允许论文被查阅、借阅和送交国家有关部门或机构，同时授 权中国工程物理研究院研究生部可以将学位论文全部或部分内容编入有关数据 库进行检索，可以影印、缩印或扫描等复制手段保存、汇编学位论文。 学位论文作者签名：i 臻荣三导师签名：酋喜军 签字f i 期：2 0 0 5 年5 月j 2 日签字日期：2 - y 年岁月 蝈 第一誊引青 第一章引言 1 9 7 3 年，r e e d 和h i l i 1 8 1 提出了间断有限元方法( d i s c o n t i n u o u sf e m ) ，用于 求鳃中子辕运方程： 洲+ d i v ( a u 、= f ， 这里仃是实数，百是一个常向量。1 9 7 4 年，l e s a i n t 和r a v i a r t 2 1 给出了这年中方 法麴牧敛性证萌。随藩，一些讦算工作者将它摇广到求解菲线往双瑟守恒律方 程和方程缀。1 9 8 2 年，c h a v e n t 和s a l z a n o 2 2 鼹渊龌有限元方法解决了一线标 量守恒律方程，具体做法是在空问方向采用d g ( d i s c o n t i n u o u sg a l e r k i n ) 离散， 时颡方淘莱用e u l e r 商箭差分。1 9 8 9 年，c o c k b u m 稻c h i w a n gs h u 6 构造( k 十1 ) 阶的r k d g ( r u n g e k u t t ad i s c o n t i n u o u sg a l e r k i n 舂限元方法求勰一维标量守溪 律方程。1 9 9 0 年，c o c k b u m 和c h i w a n gs h u 8 1 将r k d g 有限元方法成功地应 用翔多维褥量守沤棒方程。1 9 9 8 年，c o c k b u m 和c h i - w a n gs h u 9 1 完成了r k d g 有限元方法嚏多维守懂律方程组躲推广。男乡 ，潮叛奄羧元方法被应用于求瓣 椭圆方程口3 】，h a m i l t o n j a c o b i 方程 2 0 3 ，对流扩散方程b 4 1 等e 闯断裔限元方法允许空阉离散时产生间断，与一般的有限元相比，具有灵 活处理阙凝农易于处理复杂数区域边爨和边售竭题疑能力。阉粝套骧元鳇基函 数一般取为线性多项式，提高有限元的精度可以通过提高多项式的次数来实现。 在多介质问题的数值计算中，正确地遗踪运动界面豹位置相当重蒙。运动 赛藤鲍逡踩套缀多方法，磐橇予类方法、v o f 方法【2 6 l 、l e v e ls e t 方法 2 7 1 窝 波前追踪方法 4 4 1 等。1 9 6 5 年，h a r l o w 和w e l c h 2 5 1 等人提出t ( c e l l t y p e ) 格子 粪方法。该方法是e u l e r 网格上进行差分离散，在网格中布鬣若干个标记点 f m a r k e r ) 采标记滚傣，遴过这黧掭记患寒追踪运动器瑟。丽捂楚e u t e r 型豹，拣 记点是随流体一起运动，是属于l a g r a n g i a n 型的方法。该方法的缺点是，对于 一魉复杂的界面，需要的标记点太多，对般人米说狠难实现，局部采用还是 司+ 纷懿。t 9 8 1 年，h i r t 器n i c h o l s 2 6 提出了v o f ( v o l u m eo f f l u i d ) 方法。v o f 方 法怒在整个漉场中定义一个函数c ，在每个网格中，这个函数定义为一种流体( 目 标流体) 的体积和网格体积的比值。这样可以大大地降低存储量。v o f 方法是通 遘求簿v o f 墨数实现对运动器嚣懿遗踩。v o f 方法中熬v o f 凝数一觳是比较 第一章引蔷 难精确求解的。1 9 8 6 年，m a r s h a l l 4 4 对运动界面提出了波前追踪方法。该方法 是一种l a g r a n g e 方法，界面的运动韪通过界面上质点的运动来液征，这些质点 煞遴瘦是透过隶界器法方囱r i e m a n n 滴题簿 | i 褥到麴。对予毫维蠢题，缡稀比 较复杂而且计算爨非常大。1 9 8 8 年，o s h e r 斧口s e t h i a n 2 7 提出了l e v e ls e t 方法。 l e v e ls e t 方法用簿值面函数o ( x ，t ) 代替v o f 方法中的流体体积函数。让中( x ，t ) 以邋当的逸度运动，并认为零等值面就是物质界稍。在任意时刻，只瑟知道了 e ( x ，f ) 霹葵零等篷嚣，靛知遴了噩乏辩秘运动赛蚕。l e v e ls e t 方法辐对子揍予类 方法、v o f 方法和渡魏追踪方法煮缀多优点，蕾先该方法不嚣要显式她追踪物 质界面，从而比较容易处理复杂的物质界面，以及拓扑结构发生变化的问题； 其次很容易通过l e v e ls e t 函数求出运动赛谣的法向、曲率等几何性质，而且方 法缡程德擎，所簧要的诗算量嬲存储攫不大。惫了保证在任意时刻l e v e ls e t 蘧 数( i ，f ) 怒i 到物质界蕊的符号距离，必须在每个时间步后对l e v e ls e t 函数进 行羹新初始化。1 9 9 9 年，f e d k i w 2 筹在l e v e ls e t 方法的基础上引入了g h o s t f l u i d 方法。g h o s tf l u i d 方法豹主要懑想是秘蔫l e v e ls e t 函数将含有疆种流体 的浚场分成只含蠢一釉流体的骶个流场，每个流场国真实漩场瑚虚拟流场组成， 每个流场可以单独进行计算，从而避开了界面两边状态方程不同的问题。 许多流动滴题的研究，粥激渡与不同介质界甏磺襁闯断的裙互佟瘸、爨攀埠 与气体混合燃烧、液体中的气泡交形与倒塌、蕊秘滤体之闻黪 r i c h t m y e r w m e s h k o v 不稳定性、r a y l e i g h t a y l o r 不稳定性锋，都涉及到可压缩雾 介餍流羽数值模银。国内辨许多研究工作者对多奔质可嚣缩流逡行了数值禳叛。 国岁 ，如f e d k i w 2 、s h y u e 4 6 1 等。f e d k i w 用l e v e ls e t 方法追踪运动器匿，憋 e u l e r 方裰组和l e v e ls e t 方程分开求解，对e u l e r 方程组的空间离散采用 e n o l l f ( e s s e n t i a t l yn o n o s i l l a t o r y l o c a ll a x 。f r i e d r i c h s ) 方法，对闻离散采用t v d r k ( t o t a lv a r i a t i o nd i m i n i s h i n gr u n g e k u t t a ) 方法。s h y u e 攥蹬了，一m o d e l 方法。 该方法是将推导出来的比热魄方程刹e u l e r 方程缒组成一个薪舱非守瞧的方程 组，利用l e v e q u e 提出的w a v e p r o p a g a t i o n 4 7 方法求解方程组。在国内，张 镭髑g h o s tf l u i d 方法配合l e v e ls e t 方法计簿了可压缩双介震可压缩流f 2 9 】，他 和f 2 1 区别在于应用w e n o 方法对e u l e r 方摆组避褥空闽离散。刘德黝髑r k d g 第一章q 言 离散空间，应用g h o s t f l u i d 方法配合l e v e ls e t 方法计算了同种理想气体的问题 f 3 2 ，但是并没有应用到不同介质的气体、气体和液体的问题。 本文的主要工作是应用r k d g 有限元方法、l e v e ls e t 方法和g h o s tf l u i d 方法对多介质可压缩流进行数值模拟。具体的讲，就是利用这些方法数值模拟 了一维和二维的气一气，气液的双介质可压缩流；并且用改进的g h o s tf l u i d 方 法，数值模拟了一维大密度比问题：对二维流体中常见的常数流场、旋转流场 和剪切流场做了运动界面的追踪模拟。 第一章流体力学方程组和介质的状态方程 第二章流体力学方程组和介质的状态方程 在描述各种介质运动过程的数学模型中，流体力学方程组往往是基本的组 成部分。流体力学方程组是根据流体运动的三个守衡定律得到的，即质量守衡、 动量守衡与能量守衡。气体和高温高压下的液体都可以看作可压缩的流体，不 同的流体介质有不同的状态方程。在本章中我们将介绍流体力学方程组和相关 的状态方程。 2 1 流体力学方程组 二维可压缩流的基本方程是e u l e r 方程， p p u , o v e + p u p u 。+ p p u v r 三+ p ) u + 删v 册2 + p f e + p ) v = 0 ( 2 1 ) 这里p ，p 和e 分别表示密度，压力和单位质量的内能，“和v 分别表示x 和y 方 向的速度，r 为时间。e ：胪+ 旦攀表示单位体积的总能量。如果令。：o 就得到一维可压缩流的e u l e r 方程。 记u = ( p ，p u ，, o r ，e ) 7 ， f ( u ) = ( p u ，p u2 + p ，p u v ，u ( e + p ) ) ， g ( = ( p v ，p u v ，2 + p ，v ( e + p ) ) 7 ， 则方程组( 2 1 ) 可以写成 u ，+ f ( u ) 。+ g ( u ) 。= 0 ( 2 2 ) 一般情况下，压力可以写成密度和内能的函数p = p ( p ，p ) 或者写成密度和 温度的函数p = p ( p ，丁) ，同时可以导出p = e ( p ，p ) 和e = e ( p ，丁) 。声速c 表达式： 当p = p ( b g ) 时， c 一扛+ 等 塑三苎塑堡塑兰查堡丝型坌垦堕鲨查立篓 当p = p ( p ，t ) 时 。t ( p 7 ) 2 1 乃+ 了厂 f ( u ) 的雅可比矩阵的特征值和特征向量可以从下述统一公式中设4 ：1 和 b = 0 得到。类似的，g ( u ) 的雅可比矩阵的特征值和特征向量也可以设a = 0 和 b = 1 得到。f ( u ) 和g ( u ) 特征值和特征向量可以统一写成，特征值为 左特征向量为 右特征向量为 这卑 丑= 女一c ，丑= 如= i ，丑= d + c 卜隆丢，一了b t u 旦一生一旦蔓1 2 c 2 2 c 2j l 2 = ( 1 一b 2 ，6 。“，b v ，一b 。) ， l ，= ( 口，b ，一a ，0 ) ， 小陪丢，一竽谤a 一了b l v b “ h 一二 d ， r 3 = o 口 爿 v q 2 = “2 + v 2 ，五= 4 “+ b v ，口= a v b u r ：旦，c ： p h ：塑里 d 6 ，= 了f ，b 2 = 1 + 6 l q 2 - 6 ；日 令v = 0 得到一维e u l e r 方程的特征系统。 v + b c h + h c 5 。竺 ，l = 足 筻- 二章流体力学方程纵和奔质曲状态方程 2 。2 介质鲢状态方程 f n 理想气体的状态方程 p = ( y 一1 ) p e ( 2 3 ) ( 2 ) 剐毪气体的状态方程 p = ( y i ) p e y x( 2 。4 ) 浚状悫方程可以近似部分裹滏毫压下斡霆体署爨渡锩。濒万= 0 时就怒理想气体 的状态方程。 第三章多介质界面的描述 第三章多介质界面的描述 对于多介质的描述，我们采用l e v e ls e t 方法。在本章中，我们将具体介绍 l e v e ls e t 方法，重新初始化方法，g h o s tf l u i d 方法。 3 1l e v e l s e t 方法 假设问题的区域q = q ：u r u q ：，孑区域q ，中是a 物质，q ：中是b 物质， 工1 为q ，和q ：的分界面。界面r ( t ) 是随时间运动变化的。l e v e ls e t 方法是把随 时间运动的物质界面看作是某个函数中( 牙，r ) 的零等值面，巾( 王，f ) 满足一定的方 程。在任意时刻我们只要知道了巾( 哥，r ) 的函数值，就可以知道零等值面的位置， 即物质界面的位置。这也就是说我们要构造函数m ( i ，t ) 使其满足 r ( ，) = 扛q ：中( 豆，r ) = 0 。初值o ( i ，o ) 定义为t = 0 时刻牙到初始界面r ( o ) 的符 号距离 f d ( 孟，r ( o ) ) ，i q ， m ( i ，0 ) = 0 ，i r ( o ) ，( 3 1 ) i d ( 孟，r ( o ) ) ，i q 2 其中d ( i ，1 1 ( 0 ) ) 表示i 到r ( 0 ) 的距离。等值面函数法的思想就是o ( x ，r ) 以恰当的 速度运动。由于在活动界面f ( t ) 上的任意点王有 西( 王，f ) = 0( 3 2 ) 由链式法则，对( 3 2 ) 两边对t 求导可得 中，( i ，r ) 十v 巾( i ，r ) i ( r ) = 0( 3 3 ) 在一维两相流问题中，i ( f ) = “，( 3 3 ) 就是 o 。( 牙，t ) + “o ，( 哥，r ) = 0 ( 3 4 ) 在二维两相流运动中，冤( r ) = ( “，v ) ，( 3 3 ) 就是 ，( i ，r ) + u q b ；( i ，f ) + 椰。( i ，t ) = 0 ( 3 5 ) 辅三章辫介质界瑚的描述 因为l e v e ls e t 函数的单位外法方向为 目= 呙 b 6 ) 不妨设f 怒质点沿夕 法向的速度，那么有 i ( f ) = f( 3 7 ) 将( 3 6 ) 代入( 3 7 ) ，然后将( 3 7 ) 代入( 3 3 ) n n 得到由( 膏，f ) 的控制方程 蚤，( i ，f ) 十f v o 辑，0 = 0 ，牙爸q ( 3 8 ) ( 3 3 ) n i ( 3 8 ) 就是l e v e ls 娃丞数掰满蹩瓣秀种l e v e ls e t 方程。 3 2 重瓤初始化 一般由于数值方法的内在效应，进行几个时间步后，l e v e ls e t 函数不再是 到界面的符号距离。为了保持l e v e ls e t 函数始终为到界面的符号距离，我们必 矮褒每个时润步屡对其遴雩亍耋凝视始伲使其满足 i v 中 = 1 ( 3 9 ) 该方程为静态的h a m i l t o n j a c o b i 方程。中的初值为该时刻所计算出的l e v e l s e t 函数值。该初值阎题有很多遮代解法，如 3 5 和【3 6 】等，缎是这姥迭代方法需要 的遮代瞳阕一般魄较长，面曼编理魄较复杂。最来，m a r ks u s s m a n t 4 将该边 值问题转化为了初值问题 龋；三黑秽一o b 旧 通过求该初值问题的稳定解米实现重新初始化。该方法的计算简单而且一般只 需要几步迭代，肖时甚至迭代一步就可以满足要求。此外，际耀松【3 7 等人在 计算叁由嚣时，逶过求鳃双璃农方壤 1 7 = 0 孕k 1 ， o f t r ir = 1 ( 3 1 1 ) 来对l e v e ls e t 蕊数进行重囊镪始纯。罴用祝家麟提密懿边赛元方法 4 9 1 求鳞该 边德问题。该方法具有收敛快，精度高等优点。本文采用m a r ks u s s m a n 1 4 q b 笫三章多介质界面的描述 的方法进行重新初始化。 3 3 g h o s tf l u i d 方法 假设计算区域是q ，流体a 所在的区域是q ，流体b 所在的区域是q ， a 和b 充满整个区域q 。在q 内的网格称为是a 真实网格，在q ，中定义流体a 的所谓的g h o s t 网格。同理在q ：内的网格称为是b 真实网格，在q ，中定义流体 b 所谓的g h o s t 网格。根据间断关系和提高分辨率的要求给这些g h o s t 网格赋 以函数值。这样整个区域都有流体a 和流体b 的物理量的函数值，可以对两种 流体分别在整个计算区域上进行计算。 3 3 1 一维情形 一维e u l e r 方程有三个分量，密度p ，动量倒，能量e 。在f 。时刻，已知 区域q ，中流体a 的各个物理量的函数徨，在求下一时刻r 。时刻流体a 的各个 物理量的函数值，我们必须知道流体a 在q ，中各物理量的函数值，即要给出 “g h o s t 网格的各物理量的函数值。我们知道，由于物理原因，在跨过运动界面， 压力p 和速度“是连续的，所以“g h o s t ”网格上的压力和遮度的函数值应该是 取该区域的真实流体的压力和速度。这样在g h o s t 网格上e u l e r 方程的三个物理 量我们已经定义了两个，还有一个自由度需要确定。由于熵s 在接触间断中普 遍是间断的，所以我们选它作为第三个需要确定的自由度。 熵可以通过外插值来确定它在“g h o s t ”网格上的值，由于高阶插值会产生 所谓的“o v e r h e a t i n g ”误差，所以一般我们选择常数外插值。我们将用图形来 详细说明“g h o s t ”方法。假设运动界面落在节点i 和i + 】之间。这时流体a 的 真实值在i 点及其左边有定义，而流体b 的真实值在f + 1 点及其右侧有定义。我 们在计算时，需要知道流体a 在i + 1 点及其右侧的g h o s t 值，流体b 在i 点及 其左边的g h o s t 值。下面以计算流体a 的函数值为例，说明g h o s tf l u i d 方法。 流体a 在i + 1 点及其右侧的压力和速度值可以直接从流体b 的真实网格上复制 过来。而对于熵我们采用常数插值，即 9 燮墨尘垦堂里些堂堕 i - 2 文l = s t ，s = s j ，s 2 3 = s j ， 界面 i - i 梳体a i + 3 i ( 翊3 1 ) 为t 削弱“o v e r h e a t i n g ”现象，我们采用所谓的“i s o b a r i c 技术。即将f 点的熵 修敬为i 一1 点的熵，即 s 尸= s 。，s 三，。s 。，s ；! ! ：= s 。， 嚣蔼 i 一2 i + 3 ( 吲32 ) 对予一维情形，我们可落致避g h o s tf l u i d 方法。t gl i u 4 提高了带 “i s e n t r o p i c ”修正的g h o s t f l u i d 方法，鼹予强激波与运动雾嚣熬相互l 乍爱。该 方法是利用特征关系和激波跳跃条件，通过选代的方法求出界面的压力和速度， 避蔼求蠢弊面左边和右边酌密度和滴，再将熵延强割g h o s t 区域。本文改进的 g h o s tf l u i d 方法与tg - l i u 方法不嗣的是邋过求织器嚣黔近螅髑部r i e m a n n 翊 题得到界颟的压力、速度以及界面左边和界颟右边的密度和熵。现在我们利用f n 时刻已知的各物理量值预测，”1 时刻界面附近各物理量的值。设p ，表示界面的 压力，“，界面的速度，p j 界面左边的密度，p ? 界颇右边的密度，皿界面左边 第二章多介质界面的描述 的熵，艘界面右边的熵，则由如下方法可以求得，”1 时刻界面附近的各物理量 的值，求解r i e m a n n 问题 1 0 ：，p 三，p 品) ，界面左边， p 2 1 1 2 ，p 卫2 ，p 0 2 ) ，界面右边 这里“品，p 二l ，p 二分别表示节点( i 一1 ) 在f ”时刻的速度、密度和压力，“二：，以：， 苁。分别表示节点( i + 2 ) 在 ”时刻的速度，密度和压力。根据文献【3 3 ，上述 r i e m a n n 问题的解为五种类型解中的一种，可以求出接触间断附近的压力p 和 速度u 、接触间断左波的波后密度d l 以及接触间断右波的波后密度d r 。那么 p ，和“，分别用p 和u 代替，硝和分别用d l 和d r 来代替。以计算流体a 的各个物理量的函数值为例来说明改进的g h o s tf l u i d 方法。流体a 在( i + 1 ) 点的 速度和压力的值用界面的压力和速度代替，而( 件2 ) 点及其右边各点的压力和速 度的值可以直接从流体b 的真实网格上复制过来。而f i + 1 ) 的熵由界面左边的熵 s l 延拓过来。如图3 3 所示。为了削弱“o v e r h e a t i n g ”现象，我们采用t g l i u 2 提出的“i s e n t r o p i c ”修正。就是将f 点的熵修改为s l ，如图3 4 所示。 i + li + 2 u i ip ru l p 00 s l( ：)o i 一1 i、k = ! 一r 流体ag h o s t 罔格 ( 图3 3 ) 界面 1 十i1 十z u l l p iu lp j s loo i i：1 、y 一，r 流体ag h o s t 随格 ( 图3 4 ) 第二章多介质界面的描述 3 3 2 多维情形 一维问题的“g h o s tf l u i d ”方法和“i s o b a r i c ”技术实现比较容易，对于多 维问题，由于速度有多个分量，所以插值的方向需要选择。对于多维情况，压 力和速度的处理用通常的办法，“g h o s t ”流体的压力和速度和真实流体的值相 等。但是为了完成“g h o s t ”区域的赋值，我们还是需要把需要插值的变量延拓 到“g h o s t ”区域，我们可以通过求解下面的方程来实现： ，h v i = 0 ，( 3 1 2 ) 这里，是插值变量( 如熵s ) ，亓是每个网格点的单位法向量。由前面有关l e v e l s e t 函数的介绍，我们知道，每个网格点的单位法向量可以如下定义： 一 v 中 舻而 从i 的定义可知，它总是由o 0 的区域。所以，如果我们要 将变量从o 0 的区域延拓，( 3 1 2 ) r 扣用+ ，反之用一。“i s o b a r i c ” 技术：在求解( 3 1 2 ) 时，当用+ ，即变量从0 0 的区域延拓时， 保持巾 区域的函数值不变。s 一 般可以取为l5 a x 。 在多维情况，变量有密度p ，压力p ，速度i 。由于速度有多个方向，将 其分解到界面的法线方向亓和界面的切线方向开1 上，分别记为v ”和v 。 v ”= 厅i ，v 。= 开h 1 。我们认为法向速度是连续的，处理方法同压力p ；而切 向速度v 可能存在大的间断，需要做类似熵的处理，同样采用“i s o b a r i c ”技术， 用切向速度v 代替方程( 3 1 2 ) 中的迭代至稳定解。这样在“g h o s t ”区域的每 一个网格点我们有了两个速度场，一个是真实流体的速度，另外一个是“g h o s t ” 流体的速度。我们将真实流体网格上的法向速度v ”和“g h o s t ”网格上切线速度 v7 合成“g h o s t ”网格上的速度。 第瑙章流体力学方程组的离散 第网章流体力学方程组的离散 流体力学方程组的窳间离散有差分方法、有限体积方法和d g 有限元方法， 在差分方法中毒e n o 3 9 、w e n o 4 0 等。而时间离散一般采用t v dr u n g e k u t t a 方法。e n o 方法是遭i 毫麓商静太小蠢逶应缝选择模扳，为了褥到k 阶精度酌格 式，需要k 个单元，总体覆盖2 k 一1 个单元，包含k 个参选模板。如果覆盖2 k ，1 个单元盼攥投都被利用，可以在光潺酝域褥到2 k 。1 阶精度。w e n o 方法不蔻选 辑其中的一个模袄，甜楚剩瘸k 个横援的凸组合，扶两达蓟撬满精度豹墨的。 r k d g 有限元方法最大的优点是具有局部性，不同的剖分单元可以采用不同形 式和不固次数的多项式。两且对网格的正则性要求不离，不鬻瑟考虑一般有限 元方法中遥续牲条件就可以对弼格进行蕊密和稀蔬处理，有拳i 予自适应阚格的 生成。对于双曲型问题，利用s h u 4 8 提出的t v dr u n g e k u r a 时间离散，不仅 可以保证谯每个离散时闯中闽步格式是t v d 的，两且可啦蠢效媳糖制振端。 在本牵中，我们将掰r k d g 有陵元方法【7 】【8 】对流体力学方程组进行离散。 4 1一维流体力学方程组的d g 有限元空间离散 考虑一缳滚薅力学方疆缱载露篷囊麓 o _ 2 _ v 8 t 十塑o x _ 01 )( 4 ) u ( x ，o ) = u 。( j ) v x 若r 这墨u = ( p ，倒，切7 ，f ( u ) = ( 倒，膨2 + p ，u ( e 斗p ) ) 7 ，r 为求解区域。 厂、 将区阕r 割分为单元集合，r = u j i s , 其中2 l x ，一i 一，专j 歹= 毛2 ，v ， 单元长度觚，2 x 。! 一x ，一1 。在区间，上用分片插值多项式向堂遇近u ，记摘值 函数空潮为 磅= 访= ( p i ，p 2 ，p 3 ) 7 ：p ，最b v l 7 l 1 ，茸p ，| ，p 。( j ，) ，f = l ，2 ，3 ； 这里，p k ( ，) 表示定义在，上的次数不大于k 的多项式空间，我们选取 ( ，) = 舻8 行；文z ) 墨。，v ；乃( z ) 为，上的l e g e n d r e 多项式： 第必章溅体力学方程组的蕊数 诏b h ，译b ，= 竽硝k ，球竿h 。， 方程鳃( 4 1 ) 的近似解u 。圪定义为 u 。( 鼍o = u 知v 地x v x i j k 0 娃，3 隽了穗定g ( 0 ，耀试搽疆数致瞻去黍秘。1 ) ，菸在? ，上积分，糕鼹劳蘩积 分得 旦d tf ，u ( 列) ( z ) 威一f ( e ，( ) ) 丢圪( x 渺 + f ( v ( x1 ，0 ) 吒( x 】) 一f ( u ( xl ，r ) ) h ( x1 ) = o ( 4 4 ) 在( 4 4 ) 中，用氓代替u ，并用某一种单调数值通量 h ( u ) 1 ( t ) = h ( u ( x 一，f ) ，u ( x + 1 ， ) ) 代蘩f ( u ( x ，) ) ，缮到 o + i 鲁瓯( 黟) 玖( x ) 斑一f f 。( 邵) ) 芸巧( x ) 蠡 + h ( u ) t 9 ) ) 圪( xt ) 一h ( u s ) ；) ) 以( xt ) = o 。( 4 5 ) 姆阮( x ，吩的表达式( 4 + 3 ) 代入( 4 5 ) 薅 善备石dv ；一f ( 氓( 搿) ) 姜k ) 溉 h ( u 。) l8 ) ) 一一1 ) 。h ( u ) 擘) ) = 0 = 0 , 1 ，k( 4 + 国 通过数值积分，4 。6 ) 即怒我们空阍离激化最褥到鲍常微分方程。 上述数值通量的选取与有限元空间的形式没有关系。当k * 0 时，有限元空 闯为分片常数，t 主1 ( 4 6 ) 有 掣十竺燮掣观 ， d t呶 由此可见，当= 0 时，即为有限体积法，若网格削分为一致网格时，就可以看 1 4 第叫章流体力学方程组的离散 作是有限差分法。这里，单调的数值通量必须满足下面三个条件： ( 1 ) ( d ，6 ) 是局部l i p s c h i t z 连续，且h 和f 相容，朗h ( u ，= f ( ( 2 ) h ( a ，b ) 是第一个变元的不减函数。 ( 3 ) h ( a ，b ) 是第二个变元的不增函数。 数值通量用局部l a x f r i e d r i c h 通量： 爿( u ) ，+ ；= 圭 f ( u - ) + f ( u _ j ：_ ；) 一盘，+ ；( u j ；一u 二；) 这= 搿恢旷t ：；，川z ，是乱趣的实特弛 ( 4 6 ) 写成- g 微分方程简洁的形式为： 象乩“) ( 4 8 ) 4 2 二维流体力学方程组的d g 有限元空间离散 考虑二维无粘可压缩流体力学方程组 f 型+ 堡盟+ 里盟：0 o t 缸 劫 ( 4 9 ) i u ( x ，y ，o ) = u o ( x ，y ) ，v x ，yeq ，t ( o ，t ) 这里u = ( p ，p u ，p v ，e ) ， ，( = ( 用，p u 2 + p ，p u v ，u ( e + p ) ) 7 ， g ( u ) = ( p v ，p u v ，p v 2 + p ，v ( e 十p ) ) 设l 为区间q 的一个有限剖分，如三角形剖分或者四边形剖分，单元 k l ，y ( k ) 是k 上的局部有限元空间，一般我们取它为m 阶多项式0 ) 。 定义有限元空间 圪= 帆r ( q ) ：i 。y ( 世) ，v k f h 在( 3 9 ) 的两边同时乘以圪嘭，并且在k 上积分用u 。代替u ，得 旦d t ( x ，f ) v a x ) 出+ 击谚( 乩( 置r ) ) ( r ) 出= o ，v 瞄 ( 4 1 0 ) 第叫章流体力学方程组的离散 这里f ( u ) = ( f ( ，g ( ) tx = ( z ，y ) 7 。分部积分得 鲁( 列) ( z ) 出一夕( ( 列) ) - g r a d f ；, ( x ) 出 + 于( ( 列) ) n 眯圪( x ) d r = o , v v he 嗡 ( 4 1 1 ) 这里n 。是边界e 的外单位法向量。类似一维情况，我们用某一数值通量 h 眯( 己， ( x i n t ( k ) f ) ，u 0 “，r ) ) 来代替于( 乩( x ，f ) ) kk ，这里 u h ( x ( i n tk ) f ) = ，+ l 。i r a ，。ru ( _ y ，r ) ， 力= 隧冀苏它 其中h ( x ，f ) 是迹函数。我们得到 旦d t ( 蹦) 圪( x ) 出一于( 氓( 剐) ) d d r y ( x ) 出 + 。l h 啦( z ，r ) 圪o ) d f = 0 ，v 吒嗡 ( 4 1 2 ) 我们利用下面的数值积分 夕( u ，f ) ) 吒o ) d r z ：；玎，啦( x dr ) ( x 。) h ， ( 4 1 3 ) 于( u 。( x ，r ) ) - g r a d v h ( x ) c “：。珂。于( u 。( x 日，) ) g r d d 圪 目) l k l ( 4 1 4 ) 将( 4 1 3 ) 和( 4 1 4 ) 代入( 4 1 2 ) 得 夏d * u 一( t ，) 巧( z ) 出+ 。：疗皿，。( z “，) 矿( t f ) 。1 一：口，7 ( u 。( x 崎，r ) ) g r a d v n ( x x 。) 世l = o ，v k 嘭，v k l ( 4 1 5 ) 数值通量日。应该满足下面三个条件： ( 1 ) 日畎( u ，u ) 是n n l i p s c h i t z 连续，且与于( u ) n 。相容，即 日。( u ，u ) = 夕( u ) - n 。，。 ( 2 ) h 。( u ，u ) 是第一个变元的不减函数。 ( 3 ) h 。( u ，是第二个变元的不增函数。 ( 4 ) h 。( u ，u ) 也满足守恒性，即 第”q 章流体力学方程组的离散 h 蹦( ( x i n t ( k ) , r ) ，u ( xe x t ( k ) f ) ) + h “( ( x i n t ( k ) f ) ，u ( x e x t ( k ) ，r ) ) = 0 ( 4 ，1 6 ) 数值通量用局部l a x f r i e d r i c h s 通量： 1 r 一 一 h 啦( 口，6 ) = 去扩( 口) 心，f + 于( 6 ) 心。一t z e , k ( 6 这里口嘣是j a c 。b i 矩阵望塾鍪型塑n 。在边界e 邻近的最大特征值的估计。 o u 上述半离散格式( 4 1 5 ) 写成常微分方程为 三乩乩a ( 4 1 7 ) 4 3 r u n g e k u t t a 时间离散 在数值求解( 4 8 ) 和( 4 1 7 ) 时，采用r u n g e k u t t a 方法。我们知道在解充分光 滑的条件下，当用k 阶有限元空间时，空间离散可以达到k + 1 阶精度，因此为 了时空精度一致，我们在时间方向上用k + 1 阶的r u n g e k u t t a 方法。为了提高 非线性的稳定性，在r k d g 方法中，每一时刻的数值解也要进行限制。下面重 点介绍一下r u n g e k u t t a 方法的主要步骤。首先将时间 o ，丁 剖分为 “) ， a t ”= i f ”，竹= 0 ，n i 。我们利用k + 1 阶的t v dr u n g e k u t t a 方法来离 散( 4 8 ) 和( 4 1 7 ) ，基本步骤如下： u ? = u 对 = 0 ，n 一1 ，由u ：计算u ：“的过程如下： ( 1 ) 令u ：= u ：； ( 2 ) 对i = 1 ，k + 1 计算 u ：芝b 。，u ，) + 屈，f 一厶( u n ，f n + 4 f 一) ； ( 4 1 8 ) ( 3 ) 令u = u 。 加上限制算子a f i 。后，上述算法变为 u ：= a 1 7 u o ； 对n = 0 ，n 一1 ，由w 计算u ：“的过程如下 第蚪章流体力学方糕组的离散 ( i ) 令v ：。) = n ； ( 2 ) 对i = 1 ，k + 1 计算 “= a f i 。兰b 。u i ，) + 最，f 一三 ( ，f n + 每a f ”) 】； ( 4 1 9 ) ( 3 ) 令u ；玎。 二阶和三阶的r u n g e k u t t a 方法的参数取值如下： ( 1 ) 二阶( k 十1 = 2 ) ：口1 0 = 屈o = 1 ；甜2 0 = 甜2 1 = “2 2 叠 ； 反o = 0 ；d o = 0 ：d 1 = 1 ， ( 2 ) 三阶( + 1 = 3 ) ：口i o = 屈。= 1 ；搿2 0 = ：反o = 0 ； d ：，= 岛，= i 1 ；搿。= j 1 ；氏= 口，= 岛。= o ； 9 3 2 = 羁2 = 寺：d o = o ：d = 1 ；d 2 = 去 弓l 入耀剑嚣，在不彩蛹毙潺区域鹃精度懿祭 警下，季罄囊# 貔理薮荡。隈裁 器既可以遐t v d 的也可以是t v b 的。在一维的情况下，适当的选择限制器， 可以构造t v d 格式。由于离维双曲方程组没有t v d 性质，所以高维时采用 w 8 验疆裁嚣。维2 凌况，不戆薹；| 设蹬多颈式鹃寿隈元空阕 u f , o ，如= u ? ( ) v 。( x ) + u ：( ) v ，( x ) ，在f “时刻限制u ：( f ) 在该时刻的谯，我们甩 e ，：= 删( u ；( f ”) ，u ；+ 。( r ”) 一u ；o 。7 ) ，u ；) 一u a ( f ”) ) 来代替u ；( f ”) 。二维情况， 假设区域剖分为矩形单元和一阶多项式的有限元空间，矩形单元 k ；i x ；，x 】【y 。yl 】，阮( x ，y ，f ) = ( 砖+ u ：( f ) o 。( 并) + u ，( ) 昏， 0 的区域惩拓，一 表示中 0 的区域向 0 ，我们取。如果肝， 0 ，我们取c 。如果 ，= 0 ，我们不做选择。 同样的过程适用于。a 对于 i ，一i v ，= 0 ， 我们用一n 。代替上面的q ，一胛，代替上面的n ，重复上面的过程即可。 第六章数值实验 第六章数值实验 本章分为两部分，第一部分是一维多介质可压缩流的数值模拟。第二部分 是二维运动界面的追踪和二维多介质可压缩流的数值模拟。对于一维和二维多 介质可压缩流的数值模拟，e u l e r 方程组、l e v e ls e t 方程和重新初始化方程都采 用三阶的r k d g 有限元方法求解，延拓方程用一阶的迎风格式求解。二维运动 界面的追踪应用二阶的r k d g 方法。 6 1 一维多介质可压缩流的数值模拟 用改进的带“i s e n t r o p i c ”修正的g h o s tf l u i d 方法进行了数值实验。数值算 例的计算解用圆圈表示，精确解用实线表示。图1 图1 0 给出了各个算例的密度、 速度和压力的图像。 算例1 一维气气两相r i e m a n n 问题