（计算数学专业论文）二维问题的一种积分形式的流量重构算法.pdf

中国民航大学硕士学位论文 摘要 二阶椭圆型边界值问题一v k v p = f 中p 可解释为压力， u = 1 i p 可解释为流量。在实际应用中，此问题的数值求解往往不 仅仅限于计算压力p 的逼近解既，更重要的是寻求流量u 的逼近解 u 。标准的g a l e r k i n 有限元方法产生的线性方程组的系数矩阵是对 称正定的，因而是易于求解的，但它不能自动得到流量的有效逼近。 混合有限元方法是实现我们目标的常用方法，具有所得逼近流量满足 局部守恒性质的优越性，但是各种混合有限元方法有一个共同的缺陷 是产生的线性方程组的系数矩阵是非对称正定的，因而通常的数值代 数( 例如共轭梯度法) 方法不易处理，并且计算量较大。本文以二维 问题为例提出一种积分形式的流量重构算法，分别计算逼近压力p 和逼近流量u 。首先，采用标准的g a l e r k i n 非协调线性有限元格式 求得逼近压力p 。；然后，在每一个单元k 上，直接构造一个极其简 单的逼近流量的局部计算公式，直接得到精确流量u 的有效逼近u 。 此方法不仅有效地克服了通常混合有限元方法的上述缺陷，大大节省 了计算工作量，而且保持了混合有限元方法满足局部守恒性质的优越 性。本文严格分析了算法的收敛性，数值模拟试验显示了该算法的优 越性。 关键词有限元方法，流量重构算法，局部守恒性，收敛性 中国民航大学硕士学位论文 a b s t r a c t f o rt h es e c o n do r d e re l l i p t i ce q u a t i o n v r v p2f t h ev a r i a b l e p ，a sas t a t ev a r i a b l e ，c a nb ei n t e r p r e t e da sp r e s s u r e ，a n dt h ev a r i a b l e u = 叫c 即，a saf l u xv a r i a b l e ，c a nb ei n t e r p r e t e da sd i f f u s i v ef l u x i n a p p l i c a t i o n s ，i ti sd e s i r a b l et oh a v en o to n l ya na c c u r a t ea p p r o x i m a t e p r e s s u r ep h ，b u ta l s o a na c c u r a t e a p p r o x i m a t e f l u x a w h i l et h e s t a n d a r dg a l e r k i nf i n i t ee l e m e n tm e t h o dr e s u l t s i n “e a s y - t o - s o l v e ” s y m m e t r i cp o s i t i v ed e f i n i t ef m i t ee l e m e n ts y s t e m s ，i td o e sn o tp r o v i d e a c c u r a t ef l u xa u t o m a t i c a l l y t h em i x e dm e t h o da p p r o a c hc a np r o v i d e a c c u r a t ef l u xa n di sl o c a l l yc o n s e r v a t i v e b u ti nt h em i x e da p p r o a c h ，o n e h a st of a c es o l v i n ga ni n d e f i n i t es y m m e t r i cs y s t e mr e s u l tf r o mt h es a d d l e p o i n tf o r m u l a t i o n w ee v a l u a t et h ea p p r o x i m a t ef l u xu t o t h ee x a c tf l u x a = 却b y as i m p l eb u tp h y s i c a l l yi n t u i t i v ef o r m u l ao v e re a c hf i n i t e e l e m e n t o n c ep hi sg o t ，t h e n u c a n b ec o m p u t e dd i r e c t l yv i aav e r y s i m p l ef o r m u l a o v e re a c hf m i t ee l e m e n tw i t h o u th a v i n gt ou s ea n ym i x e d f i n i t ee l e m e n tm e t h o d t h u st h i sm e t h o ds a v e sl o t so fc o m p u t i n gw o r k b e s i d e si tm a i n t a i n sl o c a lc o n s e r v a t i o np r o p e r t ya tt h ee l e m e n tl e v e l s i n c et h ec r u xo ft h i sm e t h o di s c o m p u t i n gt h ev a l u e so fu a tt h e b o r d e rb yt h ei n t e g r a lf o r m u l a s ，s ow ec a l lt h i sm e t h o daf l u xr e c o v e r y m e t h o d t a k i n gi n t e g r a lf o r mi nt h i sp a p e r , w ew i l li n t r o d u c et h i sm e t h o d t h r o u g h a ni n t u i t i v ei d e a t h e nw ea n a l y z et h ec o n v e r g e n c eo ft h i s m e t h o d n u m e r i c a lr e s u l t sa r ep r o v i d e dt os u p p o r tt h et h e o r y k e yw o r d s ：f i n i t ee l e m e n tm e t h o d ，f l u xr e c o v e r ym e t h o d ，l o c a l c o n s e r v a t i o np r o p e r t y , c o n v e r g e n c e 1 1 中国民航大学学位论文独创性声明 本人声明所呈交的学位论文是我个人在导师指导下进行的研究工作及取得的研究成果。尽我所 知，除了文中特别加以标注和致谢的地方外，论文中不包含其他人已经发表或撰写过的研究成果， 也不包含为获得中国民航大学或其它教育机构的学位或证书而使用过的材料。与我一同工作的同志 对本研究所做的任何贡献均己在论文中作了明确的说明并表示了谢意。 研究生签名： 枣赵 日期：2 芝芝2 ：室：；7 中国民航大学学位论文使用授权声明 中国民航大学、中国科学技术信息研究所、国家图书馆有权保留本人所送交学位论文的复印件 和电子文档，可以采用影印、缩印或其他复制手段保存论文。本人电子文档的内容和纸质论文的内 容相一致。除在保密期内的保密论文外，允许论文被查阅和借阅，可以公布( 包括刊登) 论文的全 部或部分内容。论文的公布( 包括刊登) 授权中国民航大学研究生部办理。 研究叠签名：垂奁 导师签名：日期：幽：! ：刁 第一章绪论 有限元方法是数值求解偏微分方程的基本方法之一，有限元方法 最突出的优点是能够容易地处理求解复杂区域上的问题，这也正是在 过去的几十年中，有限元方法在弹性力学，流体力学以及大量的工程 问题的计算中获得极大成功的主要原因之一。它的数值分析理论已相 当完善。众多的文献资料报告了有限元方法的研究成果。例如， c i a r l e t “3 总结了椭圆型问题的有限元解法，t h o m e e 乜1 总结了抛物型问 题的有限元解法，而g i r a u l t 口3 ，c h u n g h l ，t e m a n 嵋1 等则阐述了流体动 力学中的有限元方法。 在用有限元方法数值模拟流体流动时，常常需要同时求解两个未 知函数，而且常常同时求一个标量、一个向量，例如如下二阶椭圆型 边界值问题即为此类问题的一个典型代表。 r v ( k v p )= 厂x e q 、p=0x e a q ( 1 1 ) 其中q 为r 2 中的多边形区域，铀为q 之边界，- c = k ( x ) 是一对称正定 的矩阵函数，并且满足条件：存在两个正常数峨和p ：，使得如下不等 式成立 q 号7 芎毛k ( x 圬曝飞，芎r 2 ，x q 。 ( 1 2 ) 引进流量函数u = 咪v p ，则问题( 1 1 ) 可以等价地表示为如下一阶偏 微分方程组 iv u = 厂，x q ， u + k v p = 0 ，x q ，( 1 3 ) 【p = o ，x 触 问题( 1 3 ) 可以解释为忽略重力作用的情况下，不可压缩单向流体在 介质中流动的数学模型。( 1 3 ) 中矩阵k 是介质浸透率与流体粘性系 数之比，其值决定于介质与流体的物理属性，未知函数p 可解释为压 力，而未知流量函数u = - - r v p 可解释为d a r c y 流速。在实际应用中， 数值求解问题( 1 3 ) 往往不仅仅限于计算压力p 的逼近解p ，更重要 的是寻求流量u 的逼近解u 。众所周知，标准的g a l e r k i n 有限元方法 产生的线性方程组的系数矩阵是对称正定的，因而是易于求解的，但 是它不能自动得到流量的有效逼近。这类问题常常用混合有限元方法 处理，基于( 1 3 ) 的混合有限元方法不仅给出精确的逼近流量，而且 逼近流量满足局部守恒性，这是各种混合有限元方法的主要优势。也 正是由于这个原因，混合有限元方法的理论与应用研究一直是有限元 方法中的重要研究方向之一1 。另方面，也可以采用混合有限体 积方法数值求解问题( 1 3 ) 睁“1 ，混合有限体积方法同样具有混合有 限元方法的上述优越性( 事实上，局部守恒性本来就是有限体积方法 的出发点! ) 。但是各种混合有限元方法和混合有限体积方法有一个共 同的缺陷产生的线性方程组的系数矩阵是非对称正定的，因而通常的 数值代数( 例如共轭梯度法) 方法不易处理，另外，由于同时求解两 个未知函数，混合有限元方法和混合有限体积方法计算量都较大。 为了克服混合有限元与混合有限体积方法的上述缺陷，自然提出 非常有实际价值的这样一个问题：可否分别计算逼近压力和逼近流 量? 也就是说，首先通过对问题( 1 1 ) 采用标准的协调或非协调 g a l e r k i n 有限元方法，得到压力逼近解以，其次在每一个有限元单元 上，构造简单明了的显式公式，直接计算精确流量 u h ( d i v ；t a ) = w ：w l 2 ( q ) ，v w e l 2 ( n ) 的逼近解u ，并且逼 、， 近流量是局部守恒的，从而无需求解任何大型线性方程组。换言之， 这种想法如果得以实现，新方法将保留混合有限元方法或混合有限体 积方法的优越性，但有效地克服克服了它们的局限性，这种方法就是 所谓流量重构算法。 美籍华人s o - h s i a n gc h o u 教授提出一种满足局部守恒性的流量 重构算法“”1 ，其基本思想简述如下。 假设对多边形区域q 作三角剖分，并且用足表示任意取定的一个 三角形单元。假设逼近压力p 。已经通过标准的协调或非协调g a l e r k i n 有限元方法获得，一个自然的想法是通过形如 中国民航大学硕士学位论文 u 。= 一心弧+ 修正项 ( 1 4 ) 的定义于三角形单元k 上的局部表达式直接计算出逼近流量u ，其 中石= 志正磁是心= 圪( x ) 在三角形单元置上的平均值，而f 五l 表 示三角形单元k 的面积。( 1 4 ) 中第一项一麻源于模拟u = _ c 晒的 朴素想法，另加一个修正项则是满足局部守恒性的需要。事实上，若 不加修正项，则显然在眉上恒成立v u h = 0 ，但是精确流量u 满足 驴u = f 或者在单元z 上成立，v u a x = r 弘。也就是说，逼近 v ， 流量应当通过加修正项方法，使其在任意一个三角形单元置上满 足v u h = 厶，亦即使得逼近流量u 满足局部守恒性质，这里 石= 高z 弦是在三角形单元石上的平均值。 现在，用t a y l o r 展开的观点来看( 1 4 ) 应当具有何种形式。假设 在任意一个三角形单元量上逼近流量取线性函数，我们可以把u z 在三角形单元鼻的重心口点处展开，有 ( 矽= u 矗0 矽+ z 协 ( 撇一x x 石( 1 5 ) 其中伽。( = 矽是例的j a c o b i a n 矩阵。如果我们用最低阶连续的( i n n o r m a lc o m p o n e n t ) r a v i a r t t h o m a s 空间作为流量u 的逼近空间，亦即 在三角形单元五上，逼近流量u = 似+ 魄c + b y ) 2 其中仍b 厦， x = 瓴力1 f ，那么我们容易得到关系式 1 d u h ( b ) ( x x 口) = b ( x x b ) = ：誓u & ( x x 8 ) o0 1 勘 为了保证u 。的局部守恒性，直接令 兮ua d x 2 l x 皿= | 灭| 足， 而这个等式又等价于 矿= 石。 把这些代入到( 1 5 ) 中可以得到下面的计算公式“朝 中国民航人学硕士学位论文 u d p = 一 ：v p h + 伍一如) ，x = ( x ，y ) 7 k ( 1 7 ) 么 其中，我们已经用唯代替了( 1 5 ) 中的第一项u h ( b ) 。( 1 7 ) 即为 满足局部守恒性质的一种流量重构算法。 我们知道，解释为d a r c y 流速的精确流量u 总是连续的，但是 ( 1 7 ) 不能保证逼近流量h ；在两个相邻的三角形单元公共边的法向 分量是连续的，这表明算法( 1 7 ) 还有一个不足之处，这就是算法( 1 7 ) 还不能模拟精确流量的这一属性。为此，s o h s i a n gc h o u 把( 1 7 ) 进一步改进为如下形式“6 u h d ，= 一为p + 三笋f x - x b ) + ，x = ( 马y ) 7e k ( 1 8 ) 其中c - 是一个常数向量修正量，其动因在于克服算法( 1 7 ) 的上述缺 点。显然，( 1 8 ) 仍然满足局部守恒性要求矿= 石，g 的值可以 通过强制要求相邻三角形单元公共边的法向分量的连续性来确定。 s o h m a n gc h o u 证明了求解q 问题的相容性( 注意显然有相容性问 题。事实上，若三角形单元眉在有限元剖分的内部，则眉之三条 边都是与其它三角形单元的公共边界，而( 是二维常数向量，于是 在形式上看就出现三个条件决定两个未知数的情形! ) ，得到了g 的 具体表达式，并且证明了ck = o ( h ) 。s o - h s i a n gc h o u 还给出了算 法( 1 8 ) 的误差估计，并进行了大量的数值模拟试验，数值结果显示 了该算法的优越性。 如果我们考虑d ( d = 1 ，2 或3 ) 维空间中的问题，那么按照上 述思想，可以得到卫中逼近流量更一般的表达式“6 3 一 ， uh ( x ) = - r v p l + 王詈( x - - x 8 ) + c k l 1 9 ) 其中q 是一个d 维常数向量修正量。 何松年对一维问题提出一种所谓积分形式的流量重构算法“，其 基本思想更为自然、简洁与明快，更容易为工程技术人员所接受，其 中国民航大学硕士学位论文 理论分析也更为方便。现在，我们简述这种算法的基本思想如f 。 考虑一维二阶椭圆边界值问题 严絮梨艇0 m 埘 【p ( d ) = p ( 6 ) = ”7 任给定亿缈的一个分割a = x o 而 一。 = 西，令 | i ，= m 。a 。x ( x , 一五一t ) ，q ，表示通常相应于节点的形状函数，即吼是连续 分段线性函数且满足嚷魄) = 6 持，这里6 。= 孟：i ： t 量= o ，l ，见 取压力p 的有限元逼近空间为 w = s p a n 鲰，甲2 叩。1 ) ， 取流量u 的有限元逼近空间为 v = s p a n q o ，币2 吼 假设通过标准的g a l e r k i n 有限元方法已经得到逼近压力最，何松 年采用直接模拟通过有限元单元边界的流量的思想构造逼近流量的 局部表达式。首先由( 1 1 0 ) 的两边同乘v a , p 。并在区间i x 。，1 上分部 积分可得 “( x t ) = 一p ( t ) p ( t ) = 一p p 和胁+ i 。如，d x ， ( 1 1 1 ) 吨一id t j 若在区间阮，x 。】上分部积分又可得 “( t ) = 一) p ( t ) = rp p q ；d x r 加，d x ( 1 1 2 ) 由于“( z ) ，p ( 工) 的连续性，自然( 1 1 1 ) ，( 1 1 2 ) 右端的值是相同的。 这表明( 1 1 1 ) 或( 1 1 2 ) 的右端可精确地“恢复”流量“( 砷在置点的 值u ( x ，) ，这启发我们自然地定义逼近流量在这点的值为 “ i ) = 一pp 知咎+ 。向，d x ， ( 1 1 3 ) 以及 ( x ? ) = rp p 知胁一r 加，d x ( 1 1 4 ) 中国民航大学硕士学位论文 这里非常自然而朴素的想法是：若p 。是p 的好的逼近，则( 1 1 3 ) ， ( 1 1 4 ) 的右端自然应当分别是( 1 1 1 ) ，( 1 1 2 ) 的右端的好的逼近。 何松年证明了u h ( z i ) = ( x ? ) ，这样可用( 1 1 3 ) 或者( 1 1 4 ) 右端积分 值作为( x ) 在一，i = 0 ，1 ，玎上的值( 同时也自然保证了逼近流量的 连续性! ) ，再把各点( t ，1 h ( t ) ) ，i = 0 ，1 ，2 ，盯依次线性连接得到连续 的分段线性函数即为逼近流量函数u 。( 工) ，显然圪从而得到逼近 流量蚝的局部表达式“7 1 h ( 工) = “ ( x ，一1 ) + z ( x t 1 ) ，x ( x ，i ，工。) ( 1 1 5 ) 其中 ( 钆) = 。p 幽二。出一e 。肌出，, - - x l r他，扭1，2，聍i-x 何松年证明了算法( 1 1 5 ) 的局部守恒性，严格证i _ i 明* - 了算法的收敛性以 及p = p ( x ) 为常函数时算法的超收敛性“”。数值模拟试验结果显示了 这种算法的优越性。由于算法( 1 1 5 ) 的核心是通过( 1 1 3 ) 或者( 1 1 4 ) 右端的积分确定( 石) 在一，i = 0 ，1 ，以上的值，因此称这种算法为一 种积分形式的流量重构算法。 显然积分形式的流量重构算法只限于解决一维问题有很大的局 限性，只有把这种算法推广到高维问题才能有理论意义和实际应用价 值。本文旨在系统研究高维空间积分形式的流量重构算法。我们将以 二维二阶椭圆型边值问题为例，把多边形区域q 作三角形剖分，采用 标准的非协调g a l e r k i n 线性有限元计算格式率先获得逼近压力，然后 沿用积分形式直接模拟通过三角形单元流量的思想，构造逼近流量在 每一个三角形单元上的表达式，给出流量的一种极其简单的有效逼 近。本文证明了逼近流量的局部守恒性以及在相邻三角形单元公共边 上法向分量的连续性，严格证明了所提出的算法的收敛性，给出了逼 近流量的严格的误差估计。此外，本文还给出了大量的数值模拟试验 结果，这些数值计算结果证明了算法的有效性。本文的理论分析方法 和技巧可以比较容易地推广到三维问题和发展型问题，事实上，本文 给出了积分形式流量重构算法分析的一般性框架。 本文将使用通常s o b o l e v 空间w i 9 的标准符号：分别以| | | | 。和h 。j 表示日一的范数和半范数，肼= 0 , 1 ，2 。当k = q 时，我们将省略下标k 。 我们也用j 1 1 | 和( ，) 分别表示r 一范数和r 一内积，也用h 州打埘表示 日( d 沁q ) 的半范数。 中国民航大学硕士学位论文 第二章算法构造及其合理性分析 2 1 算法构造 在这一节，我们的主要任务是提出积分形式的流量重构算法，给 出逼近流量的局部表达式。为实现此目标，首先做一些准备工作。 假设q e r 2 如前所述，是一个多边形区域，我们用x = ( x ，y ) 。表 示区域q 中的点，仍然考虑问题( 1 1 ) 或其等价形式( 1 3 ) ，假设非 负对称正定函数矩阵k = k ( x ) 满足条件( 1 2 ) 。 假设以= f 卅是通常的对区域q = 乩k 的不相重叠的有限元三 角形剖分，且这种剖分是正规的，即存在与乃无关的正的常数f ，满 足条件m i n 。巩p p ( 2 f ，其中p 是三角形彤的直径，以矽是三 角形髟的内切圆直径，五= 脚a ：【。“p 。只( k ) 表示定义于f 上的二元 线性函数全体。 我们将采用标准的线性非协调g a l e r k i n 有限元计算格式率先获 得逼近压力最。众所周知，标准的线性非协调g a l e r k i n 有限元逼近空 间为 k = f 岛曰( 矽j 露在每条边f 嵇的中点处连续，。 ( 2 1 ) 考虑到问题( 1 1 ) 的边界条件要求，则压力有限元逼近空间显然应取 k 的如下子空间 k 。= f 矿五，并且三角形单元落在e ol - 各边中点处取值为0 ，( 2 2 ) 流量口的逼近空间采用最低阶连续的( i nn o r m a lc o m p o n e n t ) r a v i a r t t h o m a s 空间“们 = f u e h ( d i v ；( 2 ) ：u l x r t o ( k ) ， v k e 以j ( 2 3 ) 其中 r r o ( x ) = u = ( “1 ，2 ) ：1 = a + b x ，2 = c + b y ，x = ，) ，) e k l 。 相应于图形2 - 1 ，设k = k 。，s = s ，e = e ，( 三角形足中与顶点s 对 r 中国民航大学硕士学位论文 应的边) 。我们用k 表示通常与顶点s 对应的线性基函数，即当k 限 制在三角形单元足上时，是在顶点s 处取值为1 ，而在其余两个顶点 处取值为o 的线性函数。令9 。( x ) = l 一2 九( x ) ，则( p ，是压力有限元 逼近空间瓦的与边e 对应的线性基函数，即当叩。限制在三角形单元k 上时，是在边e 的中点取值为1 ，而在另外两条边的中点取值为o 的 线性函数。容易验证如下等式成立 图2 - 1 三角单元的局部形状 v 似x ) = 斟n ( 2 4 ) 其中，i e i 表示边e 的长度，i k l 表示三角形k 的面积，而n 。表示边e 上 关于三角形k 的单位外法向量。于是，对于任何g ( x ) e ，我们有如 下局部表达式“6 1 g 。( x ) = g 。甲。( x ) ，x k ( 2 5 ) f e o k 其中，吼是q 。( x ) 在边e 中点的值。 对于任意三角形单元k ，我们总是如图2 1 所示( 例如令k = k 。) 那样，约定取逆时针方向为正向。我们现在确定流量有限元逼近空间 的相应于三角形单元三边的三个局部基函数如下。例如，对于若 k = ，则与边e = e a 对应的局部基函数为 气x ) = 习1 k i l f y x 一- 儿x , v ，y ) k ( 2 6 ) 中国民航大学硕十学位论文 注意到n 。是边e - - e a 的关于三角形k 的单位外法线向量，容易验 证 f 1 mv x p = 是墨 气，。( x ) 以= 0v x s 是( 2 7 ) 1 0 v x s s s 与其它两边对应的两个基函数气。( 功，f = 2 , 3 类似定义。由( 2 6 ) 立即可得 f 1 ，e = q ， i k ，。( x ) n e d j = to ，p = 乞， 【0 ，e = e s 对于任意u 。( x ) k ，我们得到它在三角形k 上的局部表达式 u 。( x ) l 。= 心最，( x ) ， ( 2 8 ) 其中心= l u n 凼是通过边e 的流量。 现在，我们来考虑原始问题( 1 1 ) 的弱形式：求p 域使得 a ( p ，g ) = ( ，g ) v g 崩 ( 2 9 ) 其中 a ( p ，日) 2l ( k 印) v q d x( 2 1 0 ) 于是，问题( 2 9 ) 的标准的线性非协调g a l e r k i n 有限元计算格式为： 寻求p h k 0 ，使得 a h ( p h ，q h ) = ( 厂，q h )v 吼k o( 2 1 1 ) 其中 a h ( p h ，q 一) = 莓l ( 唧。) w q h d x ( 2 1 2 ) 现在，我们转向讨论逼近流量的格式构造。假设已经通过求解格 式( 2 1 1 ) 得到逼近压力a ，来考虑如何构造逼近流量u h 的局部( 每一 个单元上) 计算公式。这种计算公式是通过一种积分形式模拟通过 三角形各边的流量而获得的。由( 2 8 ) 可以看出，解决流量重构问题 中国民航犬学硕士学位论文 的关键在于重构逼近流量u 。通过三角形各边的流量。 仍然以k = k ：，e = e ( 如图2 1 ) 为例来考虑。方程( 1 1 ) 的两 边同时乘以饥( x ) ，并在k 上积分得到 l - - v - 水印却。d x = 【，向。d x ( 2 1 3 ) 再对( 2 1 3 ) 左端的积分应用g r e e n 公式，有 l l - 4 c v p 御，d s = 。咖v ( p 。d x + 。如。d x ( 2 1 4 ) 其中n 。表示吒之边界的单位外法线向量。注意到u = 一k 跏连续性以 及纯( x ) 在f 上恒为1 ，而在置的另两条边的中点为零，保持( 2 1 4 ) 左 边在f 上的积分不变，而对其它两条边上的积分用矩形公式近似替代 ( 注意，这两边上矩形公式近似计算，其值均为01 ) 。则在其它两边 近似后，左边积分只剩下在边f 上的积分，于是得到如下的近似公式 i - k v p 吨面z 上，却v 妒, a k + ，鹿破 ( 2 1 6 ) 也就是 l 吨d r “【，一唧v 1 9 0a k + l 鹿么 ( 2 1 6 ) 这表明在三角形髟上通过边f 的精确流量可以用( 2 1 6 ) 式右边的积 分来近似表达。这为我们找到逼近流量的有效模拟方法提供了有益的 启发。 模拟( 2 1 6 ) ，我们规定 硭= i u h n 。d s = 上，一x v p h v c pd x + l ，厄级 ( 2 1 7 ) 即彩为我们要寻求的u 通过边f 的流量的模拟。从而，我们有逼近流 量的如下局部计算公式 u h ( x ) i k = u ：p k , e ( x ) ， ( 2 1 8 ) e e 孤 其中，彩由( 2 1 7 ) 计算得到。注意到( 2 1 7 ) 右端积分表达式中均为已 知量，并且即。和v 吼均为常数向量，因此比较容易计算。 中国民航大学硕士学位论文 2 2 算法的合理性 在这一节中，我们主要讨论算法( 2 1 7 ) 一( 2 1 8 ) 的合理性。所谓 合理性是指算法是否满足这样两个条件： 第一，由算法( 2 1 7 ) 一( 2 1 8 ) 得到的逼近流量是否具有连续性? 即逼近流量在相邻三角形单元的公共边上的法分量方向是否连续? 第二，由算法( 2 1 7 ) 一( 2 1 8 ) 得到的逼近流量是否满足局部守恒 性? 对于这两个问题的回答都是肯定的。首先，验证逼近流量在相邻 三角形单元公共边上的法分量是连续的。仍以图2 1 为例加以说明， e = e 。为两个相邻三角形墨和的公共边，n l 和分别为e 关于k 和 琢的单位外法线向量。分别关于三角形墨和应用公式( 2 1 7 ) ，我 们可分别得到 i u 。n 。幽= 上：一k v 叩。叔+ 。加，叔， 以及 f u 。n 。d s = t 哺v 吼d x + t 加。d x 。 显然，所谓逼近流量在边e 上的法分量是否连续，就是 i u h n l + l u 。n 。= o 是否成立! 我们有如下结果。 定理2 1 假设三角形单元k 和酶相邻，其公共边p l u 。n 。+ i u 。n 。= o 。 ( 2 1 9 ) 其中，n 工和分别为e 关于瓦和磁的单位外法线向量。 证明在式子( 2 1 1 ) 中，令q h = 饥。则有 k ( k 吼) v 叩f i x = 莓工加r 叙。 注意到基函数平。的支撑为k 。uk 。，因此上式又可简化为 l k ? ”v c p 8 d x + f k k k 勺p h v c pe d x = 0 婶e d x + 0 婶e d x ， 亦即 中周民航大学硕十学位论文 上。吼v q 。d x + l 。加。叙+ 上。_ k 。v ( p 。叙+ 。加。次= o 。 上式与( 2 1 7 ) 比较容易发现，上式正是i r l l + i u h n r = 0 。 我们注意到，从逼近压力的有限元计算格式本身，就可以非常简 洁地证明逼近流量的连续性，这说明本文提出的流量重构算法具有理 论分析非常方便的优越性。 其次，我们来验证这种流量逼近算法的局部守恒性质。 定理2 2 由算法( 2 1 7 ) ( 2 1 8 ) 得到的逼近流量u h 满足局部守恒 性质，即在任一三角形单元k 上成立 v u = 以， ( 2 2 0 ) 其中 为，在置上的平均。 证明由于 d i v u 一出= 上d i v ( 毒( f f u h n i d x ) p r 。( x ) ) 叙 = 上毒( i ，u h n i d x ) d i v p k x ) d x 而由式子( 2 4 ) 我们很容易得到 d i v 取一( x ) = 二i k i ，扛1 ，2 ，3 故有 f 【d ；。叔= 丽1 善3 ( u 。u ，出，【叔= 喜( 1 u q 出) = 窆【一慨v 叙+ 【r f t p , , d x 注意到在三角形眉上，恒成立窆叩。( x ) ；1 ，则有 喜【嘲呷吼出= 丘呷喜出= 【r - - k v p h v ( 1 ) d x = o ， 中国民航大学硕士学位论文 从而有 喜加。农= ，喜斑= 胁。 亦即( 2 2 0 ) 成立。 f k d i v u h d x = l x ，叔 1 4 中国民航大学硕士学位论文 第三章误差分析 在这一章中，我们主要讨论流量重构算法( 2 1 7 ) 一( 2 1 8 ) 的收敛 性，给出严格的误差估计。这一部分中出现的c 均表示与精确压力p 、 精确流量u 以及h 均无关的某个正的常数( 可以与己知函数矩阵 圪= k ( x ) 有关) ，在不同的地方出现，可以表示不同的值。 为证明方便起见，我们先把逼近流量的局部计算公式表示为另一 形式。 定理3 1 逼近流量有如下局部表示式 u 。( x ) k = 一以+ 窆( 加。叔) ，。( x l v x k ( 3 1 ) 其中，靠= 志k ( x ) 叔。 i = l 证明由( 2 1 8 ) 我们知道的局部表达式为 u 。( x ) i 茁= u 段，。( x ) ， 其中 u 。h = l u 一饥凼= l 吨。v 吼d x + l 加。d x 。 因此( 2 1 8 ) 又可表示为 u 。( x ) i 。= 艺( 却。w 。叙+ 加。d x ) p 。，。( x ) v x k ( 3 2 ) 由于蹶和v 吼在三角形单元足上均为常向量，因此利用( 2 4 ) 可得 i l c - r , v p h v , 叙一一l k i a x v p h v 叩。一以n 。斟 = 一以v p h 1 1 。i q i ， 于是 u 。( x ) i 。= ( 一a k v p h 。n qi e , i + f , r f , p j x ) p k 。( x ) v x e k ( 3 3 ) 中国民航大学硕士学位论文 为证( 3 1 ) ，我们只需证 3 - 4 v p h = ( - a k v p h n 。l 岛i ) k ，( x ) v x c k ( 3 4 ) i = 1 事实上，两端向量均属于8 t o ( g ) ，欲证它们为同一向量，只需验证两 端向量通过k 之三边流量相同，而这又只需验证两端向量在三边中点 与n 日，n 岛，n 岛之内积分别相等即可( 注意到( 3 4 ) 之两端向量在x 上 都是线性的! ) 。由( 2 7 ) 可得 厂3 i ( - a 石v p h - - 。l qi ) k ，。( m ，) l n 。= 一以le ii 最而( ) n 。 l t = l j = 一以n 。 同理 ( 一- 4 x v p h n blq 隅 ( m ，) f ，n = 一4 k v p h n ， j = l j 3 1 ( 一以n 。ie il ) 艮。( m ，) l n 。= 盘v p ， 其中皿表示边岛的中点( i = 1 ，2 ，3 ) 。因此( 3 4 ) 成立。 现在，我们分析逼近流量u 。对于精确流量u 的误差估计。定义半 范数 i q l 。= ( i ) v g 磁。瑁 ( 3 5 ) k “ 对于逼近压力，我们有如下误差估计“m - 2 2 。 i f ，一p h l l o + l p p 。k 2 | | ，i i 。 逼近流量u 。对于精确流量u 成立如下误差估计。 定理3 2 慨- u l l 蛳训州。皿+ i p l ，皿) ， 1 1 1 - - i a h 。丑。， - c h l f i 。 证明由( 3 4 ) 以及u = 叫印可知，当) 【时， u 。( x ) - - t l ( x ) = ( - - a x v p h 七衙p ) + 皇( k 鲰d ) 【) p 。( x ) 对于( 3 9 ) 右端第一部分，容易得到 ( 3 6 ) ( 3 7 ) ( 3 8 ) ( 3 9 ) 生里垦堕查兰堡主兰竺堡塞 i a k v p h 一蚓。恢一蛾b 慨一划晰 卜刽。m 。+ f 似阮一驯，对 s c 鼎譬| p | l 盖| k | l 一占+ c lp 一p | l 占 s c 镌| p , x + c p h p | i ? r s c k k p h - - p | l 量+ c h x p | l 点+ c | p i p l l 盖 s c h x p t u + c | p h p | l x o 对于( 3 9 ) 右端第二部分，首先来估计耳。酬i 。，由于三角剖分 是正则的，所以存在正常数q ，巳使得 c 鼻0 马x l - c 童， 因此 x 忆制1 ( y x - 一x 也, = 高峪：疆 2 高m 卜) 2 + ( y 一) 2 皿 高腼= 高同2 而h k 去 蚓气。( x ) k f 0 = i , z 3 ) a 再来估计l m 刮，有 i j 矿j - - u ：i i 。- l l ：l l 。圳吖 - - - i i a i 。j , v - 鬲- i 【) i 强e k 因此有 0 喜( 【见破) ( x ) 0 。c | 1 0 。缸 综合以上各式，得到 l i 一叫l 。 - c 强, a a 。+ i p , x ) + c i p 一，6 z ( 3 1 0 ) ( 3 i 0 ) 式两边平方，并按照三角形单元量求和，再应用( 3 6 ) 可知( 3 7 ) 成立。 1 7 中国民航大学硕士学位论文 最后，我们做h ( d i v ；l i ) 范数下的误差估计。显然成立 i i v u - v t l l 。，= i v a i i o 心二榷洲。， 上式两边平方并按照石求和立得( 3 8 ) 。 中国民航丈学硕士学位论文 第四章数值结果 这章我们主要对此类微分方程给出几个具体的例子，用这种流量 重构算法进行数值运算，比较其与真解之间的误差大小，用图表与图 形的形式进行对比。 4 1 区域的划分与离散范数的选择 这节我们对正方形区域o = o , 1 x o , 1 1 进行简单的类似如图4 - 1 的 三角划分。 强。j。 x 一” 图4 - 1 区域的三角划分形式 设 ，y p 是边( i ，力的中点，其中t = ( f 一1 ) h 2 ，y j = ( - ，一1 ) h 2 ，h = l n ， i , j = l 2 ，2 n + l 。再设p o 是( 五，y j ) 点的近似压力。我们定义： p e r t 一厶_ 【( p ( 薯，y j ) - p v ) 2 】“2 ( 4 1 ) , j = l 是p p h 绝对误差的r 范数。我们也可以定义流量通过各边的沿着法 线方向的，2 范数的误差估计式为： u e r r 一聊：= 【i ( u - - g h ) n d s 2 1 ，2 ( 4 2 ) 故我们可以用下面的式子计算p 和“的相对误差： 中国民航大学硕十学位论文 p e r r 一儿【( p ( ，y , ) - p u ) 2 】“2 t e a r , ，y j ) 2 】“2 ( 4 3 ) ，= l1 i i u e r r 一肛# 莓娶l ( u - - u h ) n d s “ t e 。e 。 p n d s 2 ) ”2 ( 4 4 ) 4 2 具体微分方程的数值结果 在这一节中我们主要对d i r i c h l e t 问题的具体方程通过我们构 造的算法进行数值计算，并且得出相对误差。对p 和在各个节点上 的值，我们用图形的形式，让其与精确值进行比较，效果更直观。 下面我们考虑问题( 1 1 ) 。 设精确的压力p = ( ，一x x y 2 一力，下面我们就不同的k 的选择进行 分别的计算。我们分别对正方形区域q = o ，l 】【o ，l 】进行h = l l o ， h ：1 2 0 ，h ：1 3 0 ，h ：1 4 0 进行计算。 如图4 1 的划分。其中h 表示每一单元直角边的长度。 例1k 取单位矩阵， c = d i a g ( 1 ，1 ) 。 当h ：1 1 0 时，通过构造算法计算，精确的虬与计算所得蚝h 如下图： 一 一 。 ，。，。v。 一“r 姊 图4 - 2 = l l o 精确的与计算所得对照图 当h = l 2 0 时，通过构造算法计算，精确的心与计算所得啦如下图： 秽 、。 ”4“”m 一 i t ，g l g 降1 删匿形对比 黼解 一y 00 x 一* 4 。4，一 图4 - 3h = l 2 0 精确的“。与计算所得矿对照图 当是= 1 3 0 时，通过计算，精确的与计算所得彩如下图： 6 5“ q 图4 - 4 = 1 3 0 精确的心与计算所得硭对照图 p 鼍 一 中国民航大学硕十学位论文 压力p 与流量”的相对误差如下表： 表4 1 压力p 与流量的相对误差 h = 1 1 0h = 1 2 0h = 1 3 0h = 1 4 0 o r d e r p e r r r l 0 0 0 5 60 0 0 1 46 3 9 8 8 e - 0 0 43 6 0 4 4 e 一0 0 4 2 u e r r 皿 0 0 0 3 48 9 0 7 0 e - 0 0 4 3 9 8 8 7 e - 0 0 42 2 4 9 6 e - 0 0 4 2 例2取k = d i a g ( 1 + 1 0 ：d + y 2 ，l + x 2 + l o y 2 ) 当h = 1 1 0 时，通过构造算法计算，精确的嘭与计算所得