（计算数学专业论文）双曲守恒律方程组高精度weno有限体积格式研究.pdf

南京航空航天大学硕士学位论文 摘要 本文在结构网格w e n o 格式的基础上，构造了一类非结构网格w e n o 有限体积 格式。这种格式不是选择最光滑的节点模板，而是对所有的多项式进行一次凸组合， 并将w e n o 格式推广到多维欧拉方程组，在光滑区域达到高阶精度，在不连续区域 权近似于零。通过减小n e w t o n 插值多项式模值的幅度，达到减小震荡的效果，又克 服了t v d 格式精度不高的缺点。 本文采用非结构网格中较为成熟的d e l a u n a y 三角化方法和阵面推进法 ( a d v a n c i n g f r o n tm e t h o d ) ，同时采用较为合理的数据结构以提高效率，能够迅速简 单地生成具有复杂外形的c f d 非结构网格。，为提高双曲守恒律方程组解的精度和激 波捕捉能力，以流场密度解的梯度为激波探浏器，发展了非结构网格的自适应技术。 深入研究了非结构网格上高阶w e n o 格式的构造，先在不同的三角形单元上构造插 值多项式，然后进行加权。对可能出现的病态方程组，分析了产生的原因，并提出了 解决的方法。给出了二元函数光滑因子的计算方法，采用了三角形作为计算网格，避 免了在复杂多边形上进行计算。对双曲守恒律方程组进行了数值实验，结果表明本文 的方法能很好的模拟复杂流场，准确的捕捉到激波，体现了格式的高精度和无震荡性 质o ，卜一，一 关键词：双曲守恒律方程组、w e n o 格式、有限体积、非结构网格、自动加密技术 _ - 一一一一一一一一一 双曲守恒律方程组高精度w e n o 有限体积格式研究 a b s t r a c t i nt h i sp a p e r ，w es t u d yt h ef i n i t ev o l u m es c h e m e so f h i g ho r d e ra c c u r a t ee s s e n t i a l l y n o n o s c i l l a t o r ys h o c k - c a p t u r i n gt ot w od i m e n s i o n a li n i t i a lb o u n d a r yv a l u ep r o b l e m s a c l a s so fe s s e n t i a l l yn o n - o s c i l l a t o r ys c h e m e sf o rt h en u m e r i c a ls i m u l a t i o no fh y p e r b o l i c e q u a t i o n sa n ds y s t e m sh a s b e e nc o n s t r u c t e d af e we x t e n s i o n sh a v eb e e na p p l i e dt o m u l t i d i m e n s i o n a ls i m u l a t i o n so fc o m p r e s s i b l ef l o w so nr e g u l a rs t r u c t u r e dm e s h e sa n d u n s t r u c t u r e dm e s h e s t h ee n os c h e m e i s c o m p o s e d o fa p i e c e w i s ep o l y n o m i a l r e c o n s t r u c t i o nf r o mt h ep i e c e w i s ec o n s t a n tv a l u e so ft h el a s tt i m es t e p i ti sb a s e do na n a d a p t i v e s e l e c t i o no fas u i t a b l es t e n c i lf o re a c hc e l ls u c ht h a ts t r o n go s c i l l a t i o n sn e a r d i s c o n t i n u i t i e sa r ea v o i d e d i ns m o o t h p a r t so f t h es o l u t i o nt h es c h e m e i so f h i g h e ro r d e r i nf i r s t c h a p t e r , o nr e c t a n g u l a rp a r t i t i o no f t h ex - yp l a n e ，t h ef i n i t ev o l u m ee n o s c h e m e sa r ed e v e l o p e di ns e m i - d i s c r e t ef o r m ，e m p l o y i n gh i g ho r d e r r u n g e - k u t t a m e t h o d s f o rt e m p o r a la c c u r a c y w es t u d yt h ee r r o r ，a c c u r a t ea n a l y s i s ，t h et w od i m e n s i o n a l r e c o n s t r u c t i o nb a s e do nt h ec e l l a v e r a g eo ft h e s o l u t i o na tt i m eta n dy i e l dag l o b a l ， p i e c e v f i s ep o l y n o m i a l o nt h eb a s eo fe n os c h e m e s ，w es t u d yt h ew e n os c h e m e s ， e s p e c i a l l yu s i n gt h en e ws m o o t h n e s sm e a s u r e m e n ti nr e s o l v i n gc o m p l i c a t e ds h o c ka n d f l o ws t m c t n r e s i nt h ec h a p t e r3 , w eg e n e r a t eu n s t r u c t u r e dm e s h e s ，m a i n l yt h et r i a n g u l a rm e s h e sb y a d v a n c i n g f r o n tm e t h o d sa n dd e l a u n a ym e t h o d i nc h a p t e r4 ，w ec o n s t r u c t h i g h o r d e rw e i g h t e d e s s e n t i a l l yn o n - o s c i l l a t o r y s c h e m e so n t w od i m e n s i o n a lu n s t r u c t u r e dm e s h e si nt h ef m i t ev o l u m ef o r m u l a t i o no fl i n e a r p o l y n o m i a l sa n dq u a d r a t i cp o l y n o m i a l s ，t h ec o m p u t a t i o no fs m o o t h n e s si n d i c a t o ra n d n o n - n e g a t i v ew e i g h tc o e f f i c i e n t so fp o l y n o m i a l s i nt h ep a p e r , w ec h o i c et h et r i a n g u l a r m e s h e sa sc o m p u t a t i o nm e s h a v o i d i n gt h ef o r m u l a t i o no f p o l y g o n i nt h ec h a p t e r5 ，w ee n f o r c en u m e r i c a le x p e r i m e n t so f g a sd y n a m i c sp r o b l e m s ，s u c h a sas t e a d ys t a t er e g u l a rs h o c kr e f l e c t i o na n dam a c h3w i n dt u n n e lw i t has t e pw h i c h d e m o n s t r a t ea c c u r a c i e sa n dr o b u s t n e s so f t h em e t h o d sf o rs h o c kc a l c u l a t i o n s k e yw o r d s ：h y p e r b o l i cs y s t e m s ，w e i g h t e de s s e n t i a l l yn o n - o s c i l l a t o r y , f i n i t ev o l u m e ， u n s t r u c t u r e dm e s h ，a d a p t i v em e t h o d i i 南京航空航天大学硕士学位论文 1 1 问题的背景 第一章绪论 作为现代流体力学的新兴学科分支一计算流体力学( c f d ) 在近二三十年得到了 迅速的发展。它涉及到计算机科学、流体力学、偏微分方程、计算几何、数值分析、 泛函分析等学科。流体力学涉及的物理现象是多方面的，如：激波、湍流、旋涡、非 定常运动等。 1 9 8 7 年，h a r t e n 和o s h e r 等人提出了三阶或更高精度的基本无振荡格式 ( e s s e n t i a l l y n o n - o s c i l l a t o r y , 简称e n o ) 1 l 。其特点是守恒的，基本上具有t v d 性质， 具有一致的高阶精度的差分格式。e n o 与t v d 相比，除了精度提高，总变差由递减 变为有界外，主要是在极点处保持了空间的高精度。l i u 等提出了w e n o 格式【”， 即在多个插值多项式中，不是选择其中的一个，而是对这些多项式进行凸组合。而后， s h u 等人又不断的进行推广和完善【5 】。 有限体积法不同与有限差分法，主要采用单元的离散，取单元的中心值( 通常也 取单元的平均值) ，而不是单元节点值：它又不同于有限元法，避免了有限元法插值 函数的构造。它既具有有限元方法网格剖分的灵活性，能逼近几何形状的复杂区域， 又具有有限差分格式构造上的多样性。由于非结构网格处理复杂外形问题特别灵活， 随着处理的气动外形日趋复杂，采用非结构网格更有生命力。网格生成是c f d 作为 空气动力学工程应用的有效工具所面临的关键技术之一。网格质量的好坏直接影响数 值结果的精度。县前，数值计算所采用的网格大体上可分为结构网格和非结构网格。 在差分格式中，基本采用结构网格，这种网格所覆盖的区域一般要求在拓扑上与一矩 形区域等价。当计算区域比较复杂时，生成高质量的结构性网格是比较困难的。非结 构网格的生成不受区域复杂性的限制，特别适用于有限体积法和有限元法的计算。目 前，非结构网格的生成方兴未艾，出现了分区对接网格、分区重叠网格。d e l a u n a y 三角化方法和阵面推进法 1 0 , 1 i , 1 2 , i 7 , i 8 , 2 0 , 3 8 , 3 9 , 5 1 1 是较为成熟的非结构网格生成方法。近年 来，在非结构网格上开展高精度格式的研究工作也有了进展。a b g r a l l l 6 】h u c h a n g q i n g 3 h ，h a r t e n 1 在非结构网格上进行了研究。自此，在非结构网格上开展高精 度格式的研究已进入了一个新的阶段。 双曲守恒律方程组高精度w e n o 有限体积格式研究 1 2 国内外研究的现状 n a r t e l l 和o s h e r 等人提出的e n o 格式【1 】，重点论述了一维有限体积的构造， 自此，人们对于高精度格式的研究进入了一个新的阶段。e n o 格式已被成功的用来 解决双曲型守恒律问题和其它的对流扩散问题瞰“。后来，c a s p e r t ”1 和h a r t e n i 6 5 把 e n o 格式有限体积格式的构造推广n - 维守恒律方程。在间断解区域，e n o 格式能 有效的捕捉到激波，但此格式也有缺陷，就是在光滑解区域不很理想。“u 等提出了 w e n o 格式”，可以克服这种不足，这种方法不是选择最光滑的网格模板，而是取它 们的一个凸组合，对组合中的每组网格模板上构造的插值多项式设置一个系数，以此 决定对数值通量逼近函数的贡献，在光滑区域可以达到( r + 1 ) 阶精度( 在e n o 格式 达到，阶精度的情况下) 。j i a n g 对w e n o 格式进行了改进和提高【4 】，提出了一种新的 光滑措施，并获得一致的高精度和无震荡，能使这种格式在光滑的区域达到高阶精度 ( 在空间上e n o 格式达到，阶精度的情况下，w e n o 格式可以达到( 2 ，一1 ) 阶精度) ， 在间断解区域权设置系数近似于零，可达到基本无震荡的效果。w e n o 格式的这种 特性在间断解区域有e n o 格式的特性，而在光滑区域类似于逆风中心格式。 上述进行的研究，大多局限于结构网格上。与结构网格相比，非结构网格处理 复杂外形问题特别灵活，随着处理的气动外形日趋复杂，采用非结构网格更有生命办。 近年来，在非结构网格上开展这方面的工作也有了一些进展，如a b g r a l l t “，h a t t e n t 4 5 】 对非结构网格上的e n o 格式进行了研究h uc h a n g q i n g ”1 】等人对w e n o 格式进行了 研究。但总起来说开展这方面的研究工作不多，还不很成熟，没有形成全面系统的理 论体系。 1 3 本文的内容和特点 首先，研究了结构网格上的高精度w e n o 格式的构造，对一维问题网格模板的 选取、插值多项式的构造、精度进行了详细的分析和论证，然后推广n - 维双曲守恒 型问题。在构造二元插值多项式时，给出了网格模板的选取方法、插值多项式光滑因 2 南京航空航天大学硕士学位论文 子的选取方法、权系数的计算，并对二维方程组w e n o 格式构造进行了论述。 其次，对二维非结构网格的生成，采用了比较成熟的d e l a u n a y 三角化方法和阵 面推进法( a d v a n c i n g f r o n tm e t h o d ) 生成非结构网格，并讨论了网格的加密技术。在 此基础上深入研究了非结构网格上高阶w e n o 格式的构造，先在不同的三角形单元 上构造插值多项式，然后进行加权。对可能出现的病态方程组，分析了产生的原因， 并提出了解决的方法。给出了二元函数光滑因子的计算方法，采用了三角形作为计算 网格，避免了在复杂多边形上进行计算。时间上采用r u n g e k u t t a 时间离散。 最后，利用生成的非结构网格对双曲守恒律方程组问题进行了数值实验。 1 4 本文的难点和处理方法 本文的主要工作是研究了结构网格上的高精度w e n o 格式的构造，对精度在理 论上进行了详细的分析和论证，特别是在构造二元插值多项式时，给出了网格模板的 选取方法，并对双曲守恒型方程组问题的插值多项式的构造进行了研究。 论述非结构网格上的w e n o 格式的构造，对二元插值多项式提出一种新的光滑 因予的计算方法。由此，得到了多项式组合的权系数计算公式。对可能出现的病态方 程组，分析了其产生的原因，并提出了解决问题的办法。构造二次插值多项式时，三 角形网格模板的选取与h uc h a n g q i n g ”l 的略有不同，本文尽可能选择靠近控制体网 格点周围的三角形网格。通过数值实验可以发现，本文的格式能更准确地捕捉到激波。 本文采用了三角形作为计算网格，避免了在三角形网格的基础上构造多边形计算 网格，从而简化了计算程序，节省了计算内存和计算时间。 3 双曲守恒律方程组高精度w e n o 有限体积格式研究 第二章结构网格上w e n o 格式 在c f d 中，提高数值算法的精度对于准确的模拟流场的物理特性具有重要的意 义，是计算数学和计算流体力学追求的目标。针对双曲型守恒律数值方法，以b r a i n v a i l l e e r 提出的m u s c l 格式为先导1 2 2 1 ，出现了全新的高分辨率守恒型差分格式，特别是 h a r t e n 提出的t v d 格式1 ，使流体力学的发展进入了个新的阶段。s w c b y t “】、 c h a k r a v a r t h y l 2 4 1 等人相继提出了一系列高分辨率t v d 差分格式。t v d 格式提高了数 值解的精度，又能有效的捕捉到激波，但是必须保证数值解的总变差不增这一特性。 h a r t e n 和o s h e r 等人在1 9 8 7 年提出了三阶或更高精度的基本无振荡格式0 1 ( e s s e n t i a l l yn o n o s c i l l a t o r y , 简称e n o ) 。后来s h u 在此基础上进行了扩展1 5 , ? 1 。这种 格式的特点是守恒的，基本上具有t v d 性质，具有一致的高阶精度的差分格式。e n o 与t v d 相比，除了精度提高，总变差由递减变为有界外，主要是在极点处保持了空 间的高精度。e n o 方法不是通过加密节点，而是逐次增加节点数的节点模板，从而 达到提高分辨率和实现高精度的目的，在每次增加节点扩充模板时，要进行n e w t o n 差商的绝对值比较。取其极小值的扩充作为新的模板。 在e n o 的基础上，l i u 等提出了w e n o 格式i2 1 ，即在多个插值多项式中，不是 选择其中的一个，而是对这些多项式进行凸组合。适当的选取振荡因子，能使这种格 式在间断解附近基本无振荡( e n o ) ，在光滑的区域有尽可能高的精度。j i a n g 等人 对其进行了改进和提高【4 】。用网格平均值构造的格式，需要从网格平均到点值的恢复 过程。本章重点对二维问题进行论述。 本章构造了一类w e n o 格式，运用数值通量，通过自适应节点模板构造差分， 分别对一维和二维双曲型守恒问题，节点模板的选取、插值多项式的构造等问题进行 了详细的论述。在时间上采用t v d r u n g e k u t t a 离散方法。 2 1 一维守恒方程的离散 考虑一维守恒方程 “= ( “。，“：，) 7， ) = ( z ( ”) ， ( “) ，工( ”) ) 7 ( 2 1 ) 4 南京航空航天大学硕士学位论文 写成半离散的形式 其中 知卜去c o q 砸) = 去鹾增) 嘶 ( 2 2 ) ( 2 3 ) 由于构造的插值多项式函数h ( x ) 不能保证在网格边界上连续，因而不能由上式直接积 分。选择一数值通量函数逼近在区间边界上的值。本文选择l a x - f r i e d r i c h s 通量夕， 满足夕是l i p s c h i t z 连续的；夕与厂是相容的，即夕( 叫) ：厂 ) 夕= 八卅八加圭( m ) + 伽) + 丢( 他) 一俐) 其中口= m a ) , l f ) i 2 2e n o 插值多项式的构造 众所周知，因只用到了两个节点，逆风差分格式只能达到一阶精度。在一维情况 下，t v d 格式最高只能达到二阶精度。为了得到三阶或更高精度的格式，惟有通过扩 展节点模板的方法。e n o 和w e n o 格式的构造则是利用了这种思想，从而能达到高阶 精度。 2 2 1 节点模板的选取 节点模板通常有固定模板和优化模板两种选取方法。本文采用优化模板，设选择 的节点模板为 s r ( f ) = x ，t t ) ，= o ，1 ，k 共有k + 1 组节点模板，每组有k + 1 个节点。例如对于k = 2 时，节点模板选择如下图 所示 共有三组节点模板分别为 5 捧 双曲守恒律方程组高精度w e n o 有限体积格式研究 构造e n o 插值多项式时，选择的是k + 1 组模板中的一组。 2 2 2 插值多项式的构造 看下面的两个定理1 6 定理1 设s 是r 2 中的一组相容网格，s = l ，a 2 ，a 。 ，j 是，的最大跨度 i = 1 , 2 ，n 。p = s u p r i ，设甜( x ) 在s 上有 + 1 阶导数，即d ( “) 成立，并设 m = s u p l l o u ( x ) b s o ，则对于逼近酌数p ( 砷，它的某 一系数q 满足l q | c 害，其中c 是不依赖于s 的常数。 这个定理表明，在间断曲线附近，如果网格的跨度p 趋向零，插值多项式函数 p ( 工) 的某一系数将趋向无穷。 根据特征线的走向选取初始点，如果f ( ) 0 ，选择s 2 ( 1 ) = 饥，x 。) ；否则选择 j 2 ( f ) = n 十z ，) 。不妨设，( “) 0 ，则j 2 ( f ) = “，x f + 1 ) 。在s 2 ( f ) 上的n e w t o n 线性多项 式可表示成 尸1 ( x ) = v x ，】+ v x ，x ，“ ( x 一工，) 下一步，我们将增加一个点到节点模板中。或者加入左边相邻的点x 。，或者加 入右边相邻的点x 。构成插值多项式。为达到减小震荡的目的，我们选择均差较小的 一个。即，如果 南京航空航天大学硕士学位论文 v x ，x 。，x 。+ l 】 f 矿 x ，x ，+ 1 ，z ，+ 2 】i 取屯( f ) = “+ x ，z 。) ，构成下面的多项式 s ( x ) = p 1 ( x ) + 矿【x i 一1 ，x ，x 。“l ( x x ，) ( x x ，+ 1 ) 否则，取5 ，( f ) = 融，x 。，x 。) ，构造相应的多项式。于是得到网格边界值 夕+ = p ( 一) 夕一= p ( z 。) 由上述优化模板的选取过程来看，上述模板必为七+ 1 组中的一组。 具体算法可描述如下 1 0 计算均差 日+ 【一，x m 】= 厂+ 阻。】= j 1 ( ，阻( t ) 1 + 。以扛。d 日一【t ，x j + 1 1 = 厂一【蚝1 = 圭( 厂阻( 毛) 卜。训x 】) 日【x ，j ，+ 。】= ，+ m ，”，。1 = l ( f u ( x ，) ， ( x ，+ 。) 】黜口，x 。1 ) = = 1 ，r 2 0 对于h + ，初始点取为l a a = i 3 。p o ) ( 功= h + 【地，地+ 。 ( x 一地) 4 0 根据归纳法，假设，盘1 和p “已知。则令 口”= h + 【。躲”3 删，+ 】 b ”= h + b f 掣- 1 ，x j 灿】 如果l l a 5 ，取c ”= 6 ，z 盘= z 盘1 - l 否贝t l c 忙1 = a 耻，r a t a = 蓝1 p 肚= p 耻。+ c 仲n o x ，) 5 0 p + ( x ) = p ( 功 6 0 同理，对于一，初始点取为勰= i 十1 ，重复3 0 一5 0 的步骤，, - i p a 寻至u p 。( 曲 ，。最后分忑d 州x 牝+ l 双曲守恒律方程组高精度w e n o 有限体积格式研究 2 3w e n o 插值多项式的构造 w e n o 格式的构造是建立在e n o 的基础之上，它不是选择最光滑的节点模板构 造插值多项式，而是对所有的多项式进行一次凸组合。比较早的研究w e n o 格式的 有l i u 等人。后来j i a n g 和s h u 等人对w e n o 格式进行扩展和提高。e n o 格式 的构造通过选择最光滑的节点模板达到减小振荡的目的。振荡因子，的选择是关键。 要求在间断解附近，要无振荡( e n o ) ；在光滑的区域，有尽可能高的精度；权系数 之和为1 ，即q = 1 。 设节点模板的选取如2 2 1 中所示。定义区间平均值为 = 廖( f ) 蹭 令玩= h 挑) ( | = i - 2 ，i 一1 ，i ，i + 1 ，f + 2 ) ，构造的w e n o 插值多项式应具有如下性质： 1 0 如果厅( x ) 在xl 附近是光滑的，应有hi = h ( x1 ) + o ( 缸2 1 ) 什2什2什2 2 0 如果厅( x ) 在x 。附近是间断的，应无振荡( e n o ) 。 什i 考虑守恒型方程( 2 1 ) ，正 ) 在节点z = 一处由守恒通量近似表示为 ，：( ”) l z 。- - m + ；一乒一；) 定义如下的三阶数值通量 分扣。孤。，+ 芸m j z ：；= 一吉厂( “一一- ) + 詈，( ) + ；厂( 甜“。) ( 2 4 ) h ：ooj 分扣，+ 沁。一扣m ， 如果以均差作为尺度构造e n o 格式，其通量为( 2 4 ) 中一个。五阶w e n o 数 值通量为 z + ! = ，乒：! + z 乒：! + ，z i 2222 8 南京航空航天大学硕士学位论文 对于权系数，的计算，采用j i a i l g 和s h u 的方法”1 ，一股的 铲轰一2 南 腰。为插值多项式的光滑因子，s 一般取为1 0 。“1 0 ，例如 ，：2 时，c g ：= 1 ，口：； jj ，= s 时，c ；= 去，c ? = ；，c ；= 素 2 4 光滑因子四。的计算 光滑因子是表示插值多项式函数光滑程度的度量， 度如何。四。越大，振荡越剧烈。 耻喜守妒) ) 2 出 即表明插值多项式函数振荡程 例如当r = 3 时 玛= ( m m ) 一2 f ( u l - 1 ) + ，) 2 + 去( ，昧：) 一4 f ( u i - 1 ) + 3 f ( u 1 ) ) 2 尽：= ( 厂( 札。) 一2 f ( 蚴) + ( 。) ) 2 + 去( ，( ，) 一，( 。) ) 2 玛= ( 厂( “) 一2 ，( “。) + ( “。) ) 2 + i 1 ( 3 ，( ) 一4 ，( “m ) + ，( 甜。) ) 2 2 5 一维守恒方程的w e n o 格式 ( 2 5 ) 本文选择l a ) 【f r i e 谢c l l s 通量 ，= ，+ ) + 厂 ) = 圭( 厂 ) + 删) + 圭( ， ) 一倒) 其中 口= 峄叭“) l ，未八啦o ，盖，- ( 哪。 对于夕。+ ，如果选择2 2 1 中的网格模板。取一作为初始点，对于七= 2 时， h j 相应的节点模板为。x 。，_ ，x 。，x 。) ；对于夕1 一，取x 。作为初始点，相应的节 ” 一2 9 双曲守恒律方程组高精度w e n o 有限体积格式研究 点模板为“+ _ ，x 。，x 。，x 。 ，这样近似夕。需要六个网格点的值，分别为 件j 运用上节构造的w e n o 插值多项式得到两个插值多项式函数z ：! ，所以 乒+ ：= 乒：+ 声：2 去厂( 札：) 一嚣，( 札。) + 器厂( + 杀，( 。) 一去，( ：) + 击，( 虬，) 一丽1 3 ，( + 丽4 7 ，( 。) + 云厂( 札：) 一去厂( ，) = 去厂( 雌：) 一丽1 1 厂( 。) + 丽3 4 厂( + 丽7 4 ，( 。) + 詈，( ：) 一i f ( ，) 将声+ ! 代入到( 2 2 ) 中得到五阶w e n o 格式的空间离散。 高阶r u n g e k u t t a 时间离散在2 7 中详尽描述。l a x w e n d r o f f理表明，( 2 1 ) 的解如果收敛。则收敛到它的弱解。 2 6 误差和精度分析 在间断解区域，e n o 格式能有效的捕捉到激波，但在光滑解区域不很理想。 w e n o 格式可以克服这种不足，并具有一致的高精度和无振荡。这种方法不是选择 最光滑的网格模板，而是取它们的一个凸组合。对组合中的每组网格模板上构造的插 值多项式设置一个系数，以此决定对数值通量逼近函数的贡献，能使这种格式在光滑 的区域达到高阶精度( 在空间上e n o 格式达到r 阶精度的情况下，w e n o 格式可以达 到( 2 r 一1 ) 阶精度) ，在间断解区域权系数设置近似于零。w e n o 格式的这种特性在间 断解区域有e n o 格式的特性，而在光滑区域类似于逆风中心格式。例如对于，= 3 时， 可以达到五阶精度。 设网格模板为 s k = 扛。t i ，x ，+ 一，+ 2 ，x ，+ i ) k = 0 , 1 ，r 一1 一 r l 厂+ ! = “( z + t - ，+ t ) = q z + 。一，+ ，+ 。 ( 2 6 ) 2 f t o 例如，= 2 时 七= o g ；= 一j 1 ，一。+ 三一 i = 1 g ? = 三，+ j 1 山 1 0 南京航空航天大学硕士学位论文 r = 3 时 t = 。3 = j 1 ，：一j 工一。+ 芸， 七= 1 g ? = 一i 1z 一。+ 昙z + ；+ 七= 2 q ；= j i ，+ i 5 ，。一吉，+ ：joo 为了克服e n o 格式自身的缺陷，我们不是选择最光滑的网格模板，而是利用所 有的，组网格模板的信息。取( 2 r 1 ) 个网格节点扛。，x 。一1 ，记 乒+ ! = g ：i 1 ( ，l ，一，z 一- ) ( 2 7 ) 在光滑区域可达到( 2 r 一1 ) 阶精度，而在间断解区域采用r 阶e n o 格式。我们对 每组网格模s k ( 七= 0 , i ，r 一1 ) 上的插值函数“( ，= 0 , 1 ，一1 ) ，给出相应的权 系数峨，通过它决定对逼近函数贡献的大小。对“( ，= 0 , 1 ，一1 ) 进行凸组合即 得w e n o 插值逼近多项式 42 荟喇t r ( 2 8 ) 为了达到高精度无震荡的目的，吼应满足 i o 在间断区域，纨近似于零。 2 。在光滑区域，调整魄，( 2 8 ) 中的数值通量尽可能逼近( 2 7 ) 中的夕。，即达到 什2 ( 2 r 一1 ) 阶精度。 3 。嘶= 1 。 k = o 令 比较( 2 8 ) 与( 2 9 ) 有 只有 厂+ 1 2 g 野( z + 。， r - 1 ，：+ 。) = c ；靠( + 。、，+ ，z + 。)( 2 9 ) t = 0 r - i ，；一。) + ( c o 。一c ：) “( ，+ 。，z + 。) ( 2 1 0 ) k = o 双曲守恒律方程组高精度w e n o 有限体积格式研究 0 9 i = q + o ( h 1 ) ( 2 1 1 ) 成立，才能满足( 2 7 ) ，即可达到( 2 r 一1 ) 阶精度。要满足( 2 1 1 ) ，要有下式成立 i s k = c ( 1 + o ( h 1 ) ) ，| = 0 , 1 ，一，r - 1 ，c 为常数。( 2 1 2 ) i s 。、c o 女、c ：的定义见2 3 和2 4 。 如果i s 。按照l i u 2 1 的定义方法，即 耻薯善皿掣 厂【，】表示均差，定义如下 f j ，0 】= 一 ( 2 1 3 ) f j ，刀= f j + 1 ，一1 卜f j ，一1 】 当r = 3 时 i s t = 去( ( ，【f + k 一2 , 1 1 ) 2 + ( 厂【i + 女- 1 , 1 1 ) 2 ) + ( a i + k - 2 ，2 1 ) 2 k = 0 , 1 ，2 在光滑区域，即厂0 ，将( 2 1 3 ) t a y l o r 展开 1 8 k = ( 厂知) 2 ( 1 + o ( ) ) ，七= 0 , 1 ，一1 。 o ( ) 不再是( 2 1 1 ) 式要求o ( h 1 ) ，这样c o i = q + o ( _ 1 1 ) ( 七= 0 , 1 ，r 一1 ) 。 由( 2 1 0 ) ，在光滑区域可以达到( ，+ 1 ) 阶精度。对于r = 2 时，达到三阶精度，这与 ( 2 7 ) 式一致。但是当r = 3 时，只能达到r + l = 4 阶精度。而不是2 r 一1 = 5 阶精度。 在临界点处，即厂= 0 ，当，= 3 时，对i s i ( 七= 0 , 1 ，2 ) 进行t a y l o r 展开 i s 0 = 三( ，h 一三厂锄2 ) 2 + ( ， 一圭厂锄2 ) 2 + ( 厂2 ) 2 + o ( h 5 ) i s l 丢( 厂h 一丢，锄2 ) 2 + ( 厂知+ 圭厂 1 2 ) 2 + ( 厂2 ) 2 + o ( h 5 ) 路：= 圭( + 圭) 2w h + 吾) 2 + ( ，2 + o ( h 5 ) 由于二阶求导系数随着网格模板选择的不同有变化，按这种光滑因子选取方法， 在间断区域只能达n - - 阶精度。 1 2 南京航空航天大学硕士学位论文 由2 4 中的光滑因子s k 的选取方法，在光滑区域，即f 0 ，t a y l o r 展开有 s - 2 ( ，锄2 ) 2 + 去( 2 ，知一詈厂3 ) 2 + o ( h 6 ) 11 s ：= ( ，2 ) 2 + 三( 2 ，知+ l f h 3 ) 2 + o ( h 6 ) s ，= ( 厂2 ) 2 + 圭( 2 ，知一f h 3 ) 2 + o ( h 6 ) 由上面三个展开式知8 。= ( 厂h ) 2 ( 1 + o ( 2 ) ) k = 1 , 2 ，3 这意味着，可得五阶精度的w e n o 格式。 s i = 驴知) 2 ( 1 + o ( 2 ) ) k = 0 , 1 ，一1 在临界点处，即当f 7 = 0 ，有 俗i = ( 厂2 ) 2 ( 1 + o ( h 2 ) k = 1 , 2 ，3 这意味着由( 2 5 ) 的方法，即使在临界点处也能达到高阶精度，满足( 2 1 1 ) 式。 总之，在r 阶e n o 的基础上，由( 2 1 3 ) 选取光滑因予的方法得到的w e n o 格 式，在光滑区域内有( ，+ 1 ) 阶精度，但不能达至l j ( 2 r - 1 ) 阶精度；按2 4 中的光滑因 子s 。的选取方法，在光滑区域磁l j ( 2 r - 1 ) 阶精度；在不连续区域内与e n o 格式 一致，达到，阶精度。本文采用2 4 中的光滑因子琏的选取方法 下面我们对一维守恒方程”，+ ”，= 0 ，i g o ( x ) = s i n ( 2 瓜) ，用三阶e n o 和五阶 w e n o 格式分别在网格4 0 、1 6 0 、3 2 0 上计算到时间t = 2 0 。误差和精度分别如下表 所示 三阶e n o五阶w 默0 网格 4 08 01 6 04 08 01 6 0 误差 2 9 3 e 33 7 2 e 44 7 3 e 48 8 4 e 59 6 7 e 墙3 1 1 e - 9 精度 2 9 82 9 93 0 14 9 24 9 64 9 9 1 3 双曲守恒律方程组高精度w e n o 有限体积格式研究 2 7r u n g e - k u t t a 时间离散 在时间上采用t v dr t m g e k u t t a 时间离散方法，可以写成如下形式 卜l “= ( 口m “肚+ a t f l 肫l ( u ( 七) ) i = 1 ，肌“= “”，“伽= 甜即+ 1 k ，0 c f l 系数c = r 皆凳，t 0 ，屈，t 。下是t v 。的。三阶离散如下 “1 = “”+ a t l ( u “1 “( 2 ) ：寻n 圮( 44 、 甜扣+ ) ：三“一+ 三m ( z ) + 三圮( 一) 333 、 c f l 系数取1 。四阶和五阶离散可相应给出1 $ 1 。 2 8 二维w e n o 格式的推广 2 8 1 插值多项式的构造 设控制体积为q r 2 ，二维欧拉方程可写成 “，+ f “) ，+ 厂( “) ，= 0( 2 1 4 ) 通过前面构造一维插值多项式的方法，将一维问题推广到二维守恒型方程问题。设 ( x ，y ，) 为矩形网格区域o = x 。，x ，】抄，。，】的形心，分段函数u ( x ，y ) 在网格区 一i+ i一22 域上的平均值为面，。 瓦= 上i 1 m 队( x ，y ) 桫 其中l o 口i 表示网格区域o f 的面积，节点模板的选取如下 1 4 南京航空航天大学硕士学位论文 可能用到的网格模板 对于固定的hy ，一； y y ，t ，以瓦表示在y 方向上的平均值，乃2 士z x r 知( x ，孝) 蟛 ，一一，+ 一 。 是x 的分段连续函数。对于每一乃，设巧是它的高阶近似，由前面一维j y 插 值- 7 多项式的 构造方法得到 r ( x ；f f ) = 巧( x ) = 瓦( x ) + o ( h r ) ( 2 1 5 ) 其中h x = m a x ， 。z ，这样得到一族多项式函数e ) 。 类似上面的方法，进行第二次重构，有 r ( y ；可( x ) ) = v ( x ，y ) + o ( h ，r ) ，害中杉= m a x ， ，) ，) ( 2 1 6 ) 这样可以得到二元分段光滑函数u ( x ，y ) 的高阶近似，即 r 2 ( x ，y ，石) = r ( _ y ；r ( x ；百) ) = u ( x ，y ) + o ( h 7 ) 2 8 2 二维e u l e r 方程组插值多项式的构造 在这一节中，将上节二维标量重构方法推广到双曲型守恒方程组，设控制体积为 q r 2 ，二维双曲型守恒方程组可写成 u ，+ ，( ”) ，+ ，( “) ，= 0( 2 1 7 ) “2 ( l ，“2 ，一，“。) f ( u ) = f ( u l ，1 , 1 2 ，。一，u 。) g ( u ) = g ( u 1 ，”2 ，“。) 一，b 分别为j a c o b i 矩阵，相应特征值分别为刀( “) 碟( 甜) 和矸( “) 砖( “) 。 a 的左右特征向量分别为，? ( “) ，1 2 和8 ( 甜) ，b 的左右特征向量分别为 ls 双曲守恒律方程组高精度w e n o 有限体积格式研究 宁( “) ，砖和r l b ( “) ，且满足 翘= 器暑即m = 器暑 ( 2 - 1 8 ) 通过点积，定义 ? = “，吒b = 譬“ k = 1 ，m ( 2 1 9 ) 由( 2 1 8 ) 可知 “= 卅= 国? 矿 ( 2 2 0 ) k = lk t l 方程组( 2 1 7 ) 在常系数情况下，采用( 2 1 9 ) 式的变形，可化为一族独立的标 量方程来考虑。因而，方程组( 2 1 7 ) 的重构可以用下面的方法。 给定向量“的网格平均值瓴) ，x 方向的网格平均值：对固定的，回= 巧乃。 由( 2 1 5 ) 标量方程插值多项式的构造方法，用这些均值，并利用( 2 2 0 ) 进行第一 次重构 r ( x ；订) = r ( e 国? ) = 巧( x ) ，对所有的，。 k - i 与上面类似，对固定的x ，在y 方向上有 0 , - 7 = z ：i i 由( 2 1 6 ) ，进行第二次重构，可得关于y 的重构多项式 r ( _ y ；哥( x ”= 胄( y ；群) k - i 这样，我们完成7 - 维线性方程组的插值多项式的重构，写成如下的形式 2 9 二维w e n o 格式 方程( 2 1 4 ) 和方程组( 2 1 7 ) 的半离散形式为 妄哪卜毒h 。聃0c 硼埘。夸聃宫。孚酬 2 ， 1 6 南京航空航天大学硕士学位论文 其中口f 是网格区域。”的面积，乃2 ：l 。f f 一：“( x ，y ) 咖出 夕。(百，f)-f厂(“。，鸳，雪。(石，f)=f霄(“li，)嘴i y + 一，一j”i，+ j“一2”j 令 厂( 沪圭( m ) 倒) 朋= m a x m 脚a x 2 ：( “) j g + ( “) = j ( g ( ”) 艄，= r 警麟k ( “) l 利用2 2 2 中插值多项式