（计算数学专业论文）退化抛物方程的局部间断有限元方法.pdf

南开大学学位论文版权使用授权书 本人完全了解南开大学关于收集、保存、使用学位论文的规定， 同意如下各项内容：按照学校要求提交学位论文的印刷本和电子版 本；学校有权保存学位论文的印刷本和电子版，并采用影印、缩印、 扫描、数字化或其它手段保存论文：学校有权提供目录检索以及提供 本学位论文全文或者部分的阅览服务；学校有权按有关规定向国家有 关部门或者机构送交论文的复印件和电子版；在巧i 以赢利为目的的前 提下，学校可以适当复制论文的部分或全部内容奠j 于学术活动。 学位论文作者签名： 怖矗 矽；年箩月? 。日 经指导教师同意，本学位论文属于保密，在年解密后适用 本授权书。 l i 指导教师签名：学位论文作静签名： i 勰密时间：年月日 备密级的最长保密年限及书写格式规定如下 南开大学学位论文原创性声明 本人郑重声明：所呈交的学位论文，是本人在导师指导下，进行 研究工作所取得的成果。除文中已经注明引用的内容外，本学位论文 的研究成果不包含任何他人创作的、已公开发表或者没有公开发表的 作品的内容。对本论文所涉及的研究工作做出贡献的其他个人和集 体，均已在文中以明确方式标明。本学位论文原创性声明的法律责任 由本人承担。 学位论文作者签名：坝1 书 妒( 年r 月；口日 摘要 诸多物理现象和工程问题可归结某种退化抛物型方程，譬如多孔介质中 的毛细吸附、金属焊接面分析、以及生物种群的演变等问题。这类方程兼具抛 物与双曲方程的特点，真解常具有随时间迁移的突变界面等“复杂结构”。许 多传统的数值方法不能高效准确地模拟退化抛物型方程，故而寻求能够准确 追踪突变界面的位置、宽度及其性态的数值方法是十分必要的。这是科学工程 计算中一个较富有挑战性的课题，具有广泛的应用意义和科学意义。 本文利用局部间断g a l e r k i n 有限元( l o c 出d i s c o n t i n u o u sg a l e r k i n ，简 称l d g ) 方法数值模拟退化型抛物型方程。l d g 方法是求解双曲守恒律方 程的间断有限元( d i s c o n t i n u o u sg a l e r b n ，简称d g ) 方法的成功延续与推 广，具有强健的数值稳定性和高阶精度。l d g 方法与d g 方法具有统一的数 值实现框架，可以简单地应用于退化抛物型方程的数值模拟 以退化多孔介质方程为主要模型，本文利用l d g 方法进行了大量的数值 模拟。数值结果表明，l d g 方法可以有效地应用于多孔介质方程的数值求 解，可以准确地模拟真解的“复杂结构”最后，本文还给出l d g 方法的自 适应实现策略及其应用。 关键字：多孔介质方程，退化抛物方程，l d g 方法，自适应策略 a b s t r a c t m a i l yp h y s i c a lp h e n o m e n aa n de n 百n e e r i n gp m b l e m s ，f o re ) a m p l e ，c a p i l l a r ys o r p t i o ni np o r o u sm e d i a ，t h ew e l di nc o h e r ej o i n t i n ga n dt h ee v o l u t i o o fs p e c i e s ，c a nb ed e s c r i b e db yd e g e n e r a t ep a r a b 0 1 i ce q u a t i o n s t h e s ek i n do f e q u a t i o nh a st h ec h a r a c t e r so fb o t hp a r a b o 【i ca n dh y p e r b o u e ，a n dt h es o l u t i o nm a yh a eac o m p l e xs t r u c t u r e 】f o re x a m p l e ，i n n e rl a 弘r ，b o u n d a 叫l a y e r a n dd i s c o n t i n u o u si n t e r f a c e i ng e n e r a l ，t h ec l a s s i c a ln u m e r i c a lm e t h o d sc a n n o ts i m u l a t et h e s ee q u a t i o n si nap r o p e rw a mt h e r e f o r ei ti sn e c e s s a r ya n d i m p o r t a n tt o 缸dag o o dn u m e r i c a lm e t h o df o rt h i sp r o b l e m s i nt h i sp 叩e rw es h a uc o n s i d e rt h es i m u 】a t i o nf o rt h ed e g e n e r a t ep a r a b 0 1 i c e q u a t i o n s ，b yt h el o c a ld i s c o n t i n u o u sg a l e r k i n ( l d gf o rs h o r tb e l o w ) 丘i t ee l e m e n tm e t h o d ，w h i c hi st h ee x t e n s i o no fd i s c o n n n u o u sg a l e r h nf d gf o rs h o r t b e i o w ) 丘n i t ee l e m e n tm e t h o df o rt h ec o n s e r v a t i o nl a w s t h el d gm e 七1 1 0 dn o t o n l yh a st h er o b u s tn u m e r i e a ls t a b l i l i t y jb u ta l s oh a sah i g h - o r d e ra c c u r a c y s i n c el d gm e t h o dh a st h es a m e 矗a m e w d r ka sd gm e t h o d ，j ti se a s yt oi m p i e m e l l tf o rt h ed e g e n e r a t ep a r a b o l i ce q u a t i o n s i nt h i sp a p e rw ew o u l dl i k et ot a k ea sm o d e lo n e _ d i m e n s 沁n a ld e g e n e r a t ep o r o u sm e d i ae q u a t i o n n u m i r a c a le x p e r i m e l l t ss h o wt h el 】) gm e t h o d c a ns i m l a t ee 伍c i e n t l yt h es o l u t i o no fd e g e n e r a t ep a r a b o l i ce q u a t i o nw i t h o u t n u m e r i c a lo s c 订a t i o ni nt h es h a r pl a y e r t h i si m p l i e st h a tt 1 1 i sm e t h o di so n e o fg o o dn u m e r i c a lm e t h o df b rt h ed e g e n e r a t ep a r a b o l i ce q u a t i o n f i n a l l 弘w e a l s oc o n s i d e rt h ea d a p t i v ev e r s l o no fi d gm e t h o da n dp r e s e n ts o m en u m e r i c a l r e s l l h | s k e y w o r d s ：p o r o u sm e d i ae q u a t i o n ，d e g e n e r a t ep a r a b o l i ce q u a t i o n ，l d g m e t h o d ，a d a p t i v es t r a t e g y 2 o 1 引言 诸多物理现象和科学工程( 多孔介质中的毛细吸附、金属焊接面分析、生 物种群演变等) 问题在数学上可归结为某种非线性退化抛物型方程。相应的方 程通常兼具抛物与双曲的特点：当退化程度越来越强时，真解常常具有随时间 迁移的突变界面等“复杂结构”，甚至出现间断的激波面。即使方程缺少对流 项，其真解也会出现突变界面。 退化抛物型方程的数值模拟非常具有挑战性：传统的数值方法往往不能 准确刻画出真解的具体变化，特别是突变界面的准确位置及其附近的具体性 态。因此，退化抛物型方程的数值模拟一直是数值分析的相当活跃的领域。 在有限元领域，近2 0 年来的解决途径主要集中在算法构造和网格技术两个方 面： 1 采用新型的非标准有限元算法，建立与“复杂结构”相适应的高分辨数值 格式，消除传统数值方法的数值振荡。注意到退化抛物型方程特有的混合 特征，理想的数值方法应在统一框架下，准确模拟抛物型方程以及双曲型 方程。许多非标准算法从不同的角度出发，较为有效地解决了由于真解的 “复杂结构”所导致的数值计算困难，并取得了不同程度的成功，如s u p g 方法或流线扩散方法f 1 5 1 、迎风有限元方法、稳定有限元方法、时空有限 元方法 1 6 1 ，以及各种版本的间断有限元方法 6 】。 2 建立自适应有限元算法，自动地调整计算网格，在保证数值效果的前提下 尽量降低数值计算的工作量。网格的改善可以有效提升数值逼近效果， 尤其是突变界面的动态捕捉。网格调整通常依赖于某种后验误差指标而进 行，具体方式主要有如下两种途径：其一是跟踪峰面而移动网格结点，其 二是网格的局部加密和粗化。 局部间断ga l e r k i n 有限元( l o c a ld i s c o n t i n u o u sg a l e r 虹n ，简称l d g ) 方 法是传统间断有限元( 简称d g ) 方法的延续，由双曲守恒律方程f 7 ，8 ，9 ，1 0 ，1 1 1 成功地推广到抛物型方程降，6 1 。两者具有统一的数值实现和理论分析框架， 均具有强健的数值稳定性和高阶精度，对突变界面的模拟非常尖锐且没有数 值振荡。l d g 方法允许检验函数和试探函数在跨越单元边界时间断，使得数 值格式对各单元之间的关联要求很弱，甚至允许网格剖分节点出现所谓的“悬 挂”现象而形成非结构网格因此，l d g 方法非常适宜自适应算法，其网格 局部加密和粗化处理非常容易。 4 本文将以多孔介质方程为主要模型，探讨l d g 方法求解退化抛物型方程 的数值模拟效果。针对问题本身的特点，本文对l d g 方法的数值实现进行少 许细微的调整，并进行大量的数值实验。在此基础上，我们也尝试l d g 方法 的自适应实现和数值验算。数值结果表明，l d g 方法可以简单而高效地用于 多孔介质方程的数值计算，尤其是突变界面的模拟效果令人满意。 本文整体结构如下：2 给出多孔介质方程的简单理论和数值困难，并列 举出某些常见的特解；3 详细介绍l d g 方法的具体实现过程，包括空间的 离散，时间推进，边界处理以及斜率限制器的实现；4 具体给出l d g 方法 求解多孔介质方程的数值模拟结果；5 讨论l d g 方法的自适应策略及其数 值模拟。最后，在6 给出本文的结论。 2 退化多孔介质方程 2 1理论背景 退化多孔介质方程是典型的非线性方程，它可以描述土壤中的吸附线性， 多孔介质中气体的流动等现象。下面给出其简短的推导过程，详见f 2 1 。 当单一气体通过均匀多孔介质时，其将满足所谓的质量守恒定律和d a r c y 定律： ( 1 ) ( 2 ) 其中v 是流速向量，7 表示气体的密度，芦是气体的粘性，而厂和奄分别是 多孔介质的孔隙率和渗透率。假设气体的密度7 与气体压力p 满足状态方程 1 = 询p o 其中 b 是已知常数。如果是等温过程，q = 1 ；如果是绝热过程，0 1 众所周知，问题( 5 ) 通常不具有古典 解，而其唯一弱解在豫x o ，+ 。) 上连续非负。 若初值( z ，t ) 的支集 z l 蛳( z ，t ) o ) 有界，则任何时刻的u ( z ，t ) 均具 有界支集( 0 ，q ) 。由文献 1 】可知存在等待时间虹o ，使得对于t 妊 时，q 保持定值，而当t 妊时，q 以有限的速度单调变化。 例如，多孔介质方程( 4 ) 具有著名的b 盯e n b l a t t 型解 晰玎2n 掣警) + r 1 就具有上述性质，其中u + = m a x ( u ，o ) ，k = ( m + 1 ) 显然，b 。( z ，t ) 具有 随时间流逝而增大的有界支集( 一。( ) ，。( t ) ) ，其中血。( ) = 、i 考栌。 图1 给出b a r e n b l a t t 型解在t = 1 0 和t = 2 0 的图形，其中m = 2 ，3 ，5 ，8 。 随m 的增加，突变界面变得越来越陡峭。 若初值“o ( z ) 的支集无界，多孔介质方程( 4 ) 的解将不再具有上述性质。 本文将考虑如下两类著名的特解： 1 i 型特解；核心多项式次数为1 ，它是一个行波解 其中c 是任意常数图2 给出i 型特解( 7 ) 对应m = 2 ，3 ，5 ，8 的图形，其 中c = 1 ；左图时间为t = 0 ，右图时间为t = 0 5 。随m 的增加，突变 界面也变得越来越陡峭。 2i i 型特解：核心多项式次数为2 ，在有限时间内属于古典解。 ( 8 ) 其中= 端。图3 给出i i 型特解( 8 ) 对应m = 3 ，5 ，8 的图形；左 列的时间为= o ，右列的时间为= 2 。 6 土d “ 等 = 一 、 j 一 ，i、 f f z 呱 图1b a r e n b l a t t 型解： m = 2 ，3 ，5 ，8 图2i 型特解：左图t = 0 ，右图t = 0 5 。 2 2 数值困难 模拟退化抛物型方程时，传统数值方法的效果常常令人难以满意，譬如标 准有限元方法( 简称s f e m ) 。下面给出s f e m 数值模拟多孔介质方程( 4 ) 的 一些数值结果，其中空间网格剖分采用均匀网格，时间变量采用预测校正进行 离散。 例2 1 b a r e n b l a t t 型解( 6 ) 图4 给出s f e m 求解b a r e n b l a t t 型解在t = 2 的数值结果，其中初始时 刻均为t = 1 。图4 ( a ) 中m = 2 ，而图4 ( b ) 一图4 ( e ) 中m = 8 。 数值结果表明：当m 取值较低( 参见图4 ( a ) ) 时，在支集边界处未发现 明显的数值振荡；但当m 取值变大( 参见图4 ( b ) 一图4 ( e ) ) 时，数值振荡已 7 图3i i 型特解：左列t = 0 ，右列t = 码2 ；m 从上往下依次为3 ，5 ，8 。 经清晰可见。我们指出：有限元空间次数的增高和网格尺度的缩减，并不能有 效改善支集边界处的数值振荡。图4 ( b ) 至4 ( d ) 给出有限元空间分别为蚂，p z 和p 3 的数值结果，其中单元个数锁定为e = 8 0 ；图4 ( e ) 给出采用如有限 元求解时，空间网格被加密( 单元个数e = 1 6 0 ) 后的数值结果。 例2 2 i 型特解( 7 ) 图5 ( a ) 给出s f e m 求解i 型特解在t = o 5 时刻的数值结果，其中 m = 8 。在数值模拟中，初始时刻为t = 0 ，有限元空间为如元，单元个数 为e = 4 0 。对比问题的真解( 参见图2 ) ，可以看出支集边界处出现虚假的 数值振荡。 例2 3 i i 型特解( 8 ) 图5 ( b ) 给出s f e m 求解i i 型特解在t = t o 2 = o 0 2 4 3 0 5 5 5 5 的数值结 果，其中m = 8 。在数值模拟中，初始时刻为= 0 ，有限元空间为如元， 单元个数为= 4 0 。此时数值解未发生振荡，但与问题真解( 参见图2 ) 相比较，不难注意到其最低点( 尖点) 的位置被提升。 3l d g 方法 l d g 方法是基于求解双曲守恒律的r u n g e - k u t t a 间断有限元方法 7 ，8 ， 9 ，l o ，1 1 】而发展起来，两者具有统一的数值实现和理论框架。其思想首现于 8 ( a ) m = 2 ，p 2 ，札= 8 0 1 广二= 二二二二二= = = 二= _ l 一、l ”f厂1 1 = l 卜寺j 爿 ( b ) m = 8 ，p 1 ，e = 8 0 1 厂了二二二= = = = = = 二 “5 f 。l ：j j 。爿 ( c ) m = 8 ，p 2 ，。= 8 0 1r 1 1 。1 j 一、 乱s f1 ff o l = = = 。- = 一 - 64 2 0246 ( d ) m = 8 ，p 3 ，。= 8 0 ( e ) m = 8 ，p 2 ，肌= 1 6 0 图4s f e m 求解b a r e n b l a t t 型解 ( a ) i 型特解：m = 8 ，p 2 ，e = 4 0( b ) i i 型特解：m = 8 ，p 2 ，。= 4 0 图5s f e m 求解i 型特解和i i 型特解 9 】9 9 7 年，b a s s i 和r e b a y 3 】成功地将间断有限元( d g ) 方法应用于可压缩 n a 订e r _ s t o k e s 方程的数值模拟在1 9 9 8 年， c o c k b u r n 和s h u 【6 率先给出 l d g 方法的基本数值实现框架和理论分析。 l d g 方法和d g 方法均允许试探函数和检验函数在单元边界处间断，可 以很好地模拟数值解在突变区域的“坠落”现象。具体地，两者具有很多共同 之处： l 局域并行：l d g 的空间离散是局部的，形成的质量矩阵具有块对角结构。 当l d g 空间离散配以显式r k 方法时间离散时，实际求解具有内在的并 行特征，其程序实现非常便捷。 2 高阶精度：在光滑解区域，l d g 方法至少具有阶精度；而在突变界面 附近，l d g 方法将保持解的单调性，而放弃高阶精度的追求。 3 强健的稳定性：l d g 方法本身就具有良好的l 2 稳定性，而斜率限制器的 使用更加增强了格式的稳定性。通过确保数值解的均值总变差不增，l d g 方法可以有效地控制突变区域的数值振荡。 下面给出一维l d g 方法的具体实现过程，包括l d g 空间离散、边界处 理、时间推进及其斜率限制器等。高维l d g 方法的实现可类似处理，详见 6 ，1 2 。考虑，= c ，d 】上的守恒型对流扩散方程 u + ( ，( z ) 一n ( 札) u 。】。= o ，z ，t o 乱( z ，o ) ：u o ( z ) ， z ， ( 9 a ) ( 9 b ) 其中o ( u ) 2o 是系统的扩散系数，( u ) 是系统的流动流通量 引进中间变量口= 、，五而u 。，将方程( 9 ) 改写为等价的一阶偏微分方程 组： 毗十 ，( 乱) 一再砑q i 。= o ， z ， o g 一啦( 珏) = 0 ，z ， o ( z ，o ) = u o ( z ) ， z ， ( 1 0 a ) ( 1 0 b ) ( 1 0 c ) 其中9 ( “) = ，“、，石而d u 。求解( 9 ) 的l d g 方法就是借鉴混合元方法的思 想，利用间断有限元方法求解方程组( 1 0 ) 。 1 0 3 1l d g 空间离散 设丁= l2h 一 ，弓+ ：j = l ，m ) 是求解区间f c ，明的某个网格 刹分，定义相应的间断有限元空间为 k = u 厶2 ( o ，6 | ：u f 如p k ( ) ，v 弓丁) ， ( 1 1 ) 其中巩( 弓) 是单元乃上的所有次数不超过k 的多项式全体。需强调指出， 的函数允许在跨越单元节点时间断。在l d g 方法经常使用局部l 2 投影，其 定义方式如下：对给定的函数u ( z ) l 2 ( j ) ，其局部l 2 投影p h u 是中满 足 五“( z ) ( z ) 出2z 耽珏( z ) ( z ) 如，协撬( 。) 致( ) ，j = 】，2 ，e 的唯一元素。 设w h = ( u h ，q ) 是l d g 方法所给出的数值解，其中u h ( z ，o ) = p h “o ( z ) ；当t o 时，w h ( t ) 由如下变分约束 五( 删z 2 z ( 伽沪隔以疵 h u ( w h ) j + 1 2 u h ( 。再1 2 ) + k ( w h ) j 一1 2 u h ( 。二1 2 ) ，h ； ( 1 2 a ) z 础= 六( 刮m 也血 k ( w h h ，2 强( 咯1 2 ) + k h 2 ( 尊l ，2 ) ， 机k ； f 1 2 b ) 所惟一确定，其中数值流通量矗( w ) = ( 矗。( w ) ；矗。( w h ) ) 是l d g 方法取得 成功的关键技术。通常所选取的数值流通量依赖于边界两侧的函数值，即 h ( w n ) ，+ 严h ( w 一( z j ，w 一( 嚷) 为保持d g 与l d g 方法的承继性，( 1 2 ) 中的数值流通量通常包含对流 和扩散两部分之和。换言之，其常定义为 矗( w ，w + ) = 矗。( w 一，w + ) + 矗戒，( w 一，w + ) ，( 1 3 ) 其中各部分的定义如下： 1 对流数值流通量：沿袭求解双曲守恒律时d g 方法的处理方式，定义 矗一。( w 一，w + ) = ( ，( 牡，矿) ，o ) 。， ( 1 4 ) 这里，( u 一，札+ ) 是与，( u ) 相容的满足局部“p t s c h i t z 条件的b 数值流通 量，即 ，( ) 一，( u 一，u + ) 】( u + 一一) o ，v 乱【m i n ( 札+ ，一) ，m a x ( “+ ，“一) 2 扩散数值流通量：记乒= ( 矿+ p 一) 2 ，纠= 矿一p 一是函数p 在单元边 界处的两侧平均和跳跃，其中p 土= p ( 勖士o ) 是函数在该点的右( 左) 极 限。定义 喝加一，w 卜( 一样驴丽) - 【w 】 ( 1 5 ) 其中 ，也鬈 ， 任选的参数c 1 2 = c 1 2 ( w 一，w + ) 应是局部l i p s c h i t z 连续，且当( u ) 三。 时c 1 2 兰o 注3 1 按如上方式选取的数值流通量关于各个分量是局部l i p s c h i t z 连 续的，并且相容于问题的真实流通量( ，( 札) 、，坛( u ) g ，g ( u ) ) 。，保证l d g 方 法的驴稳定性。此时的吼可以由u 局部解出，这也是l d g 方法称之为“局 部”的原因之一。当n ( 让) ；o 时，数值流通量退化到双曲守恒律的d 数值流 通量。 注3 2 对于简单的纯抛物线性问题( 即方程( 9 ) 中，；0 ，o ( ) = ) ， 文献( 6 】给出了l 2 模误差估计。当取定参数c 。三。时，若七为偶数则误差 为k + 1 阶，若为奇数则误差为阶。 3 2r u n g e k u t t a 时间离散 既然吼可以用“ 局部地显式表示，取定有限元空间的基函数后，l d g 格式可以简单地形成自由度( 为简便仍记做“h ) 关于时间的常微分方程组： 象让 = a f 一1 r h s ( 乜 ) = c ( n ) ， ( o ，丁) ， ( 1 7 ) 1 2 其中r h s ( 仙h ) 对应( 1 2 ) 等号右侧的有限元离散，质量矩阵m 具有简单的块 对角结构。特别地，若选取的局部基函数彼此正交，质量矩阵是对角阵。 求解常微分方程组( 1 7 ) 可以用各种熟知的方法求解；本文采用常用的显 式t o t “v a r i a t i o nd i m i n i s h i n gr n n 铲k u t t a ( 简称t v d r k ) 方法【1 3 。它 是e u l e r 向前积分方法的某种凸组合，具有非常强的数值稳定性。具体结构如 下： 其中r 1 是t v d r k 格式的阶数，护为时间步长。为保证格式是t v d 的，系数蚴和风i 非负且苫：蚴= 1 ，其中i = 1 ，2 ，r 。例如，求解 ( 1 7 ) 的三阶显式t v d r k 方法具体定义为 u ( 1 ) = 矿+ t c h ( 乱“) ， 札担，= ；u ”+ ；“n ，+ j t n ( “q ，) “”+ 1 = ；“”+ ；u 。，十；t c n ( 札犯，) ( 1 9 a ) ( 1 9 b ) ( 1 9 c ) 这是本文数值模拟主要采用的时间离散方式。 为保证计算的数值稳定性，t v d r k 格式要求时间步长t 必须满足c f l 条件 竺掣+ 堡粤些曼c f l m ( 2 0 ) 其中的最大值取遍此时刻的数值解覆盖区间。文献f 1 2 1 的表格列出用巩有限 元求解双曲守恒律方程的最大c f l 数c f l 。若时间离散采用k + 1 次显 式t v d r k 方法，最大的c f l 数c f l 。约等于i 南。 注3 3 t v d r k 方法较普通r 儿n 妒k u t t a 方法具有更高的数值稳定性， 特别适用于高精度格式的数值模拟，是高精度数值模拟时问离散的首选方法 之一。l4 1 称这种方法为强稳定性保持方法。 3 3 斜率限制器 整体的l 2 稳定性还不足以保证数值解不产生震荡，还需考虑更好的数值 稳定性。例如，一维的全变差有界可以更好地描述数值稳定性。为此， l d g 1 3 墩 f f 哪h uc扩 风一蚴 + u 1 | h 坩 h 血 ，l弋、0， 喀 【 | | | | 拶嘏 方法继续使用d g 中的斜率限制器，对数值解进行事后调整。它可以有效 地改善数值解的性态，消除可能出现的虚假数值振荡。 首先构造分片线性函数的斜率限制器a ：。设在每个单元l 上，有分片 线性函数 其中巧是在单元内的均值。斜率限制器a ：构造如下： 峨圳如= 州z 刮伉( ，觜，爷) ，( 2 2 ) 这里的佩是所谓的t v b 型m i n m o d 函数，具体定义是 f 梳( a 1 ，。2 ，a 3 ) = 【 f 。1 f 肛 2 ， l 1 i 肛h 2 ，s = s g n 。1 = s g n 0 2 = s g n 3 其它情形 ( 2 3 ) 其中p 是适当取定的正常数，本文统一取其为肛= 1 。相应地，单元节点处 的数值被修正为 a ：嘞( 嘻1 2 ) = 巧+ 虎( 咯l 2 一巧，巧一巧一l ，巧+ 1 一吗) a ：( _ _ 1 2 ) = 吗一伉( 吗一? 二，2 ，吗一吗一，哆+ - 一吗) 而每个单元的均值吗保持不变。 f 2 4 a ) ( 2 4 b ) 高次元( 七2 ) 的斜率限制器 8 ，1 2 】构造要利用线性元的处理技巧，其 斜率限制器a h 的具体定义如下： l 利用( 2 4 ) 计算a 珥( 嘻。) ，其中吗表示在单元乃内的均值，嚷； 表示在单元j ，内部方向的取值； 2 如果a ：( 嘻1 2 ) = ( 嘻1 2 ) ，数值解无需调整，定义a l j = l o ；否则，令a n n b = a ：畦，其中以是到线性有限元空间 的局部l 2 投影。 3 4 边界条件的处理 l d g 方法的边界条件处理非常简便：将边界条件引入到数值流通量中， 而计算格式在形式上保持不变。若边界条件为 u ( z 1 ) = 血( t ) ，u ( z 2 ) = p ( t ) ，( 2 5 ) 】4 、) 3 0 2 00n，r 1 l i 0 s o 只需重新定义边界z ，z 。处的数值流通量为 危。l 。，= 一再面两一 。l 。= 一9 ( a ( ) ) h 。1 。= 撕：丽两q 一，k 1 。= 一g ( p ( t ) ) ( 2 6 a ) ( 2 6 b ) 注3 4 边界条件对斜率限制器也有影响。此时，我们引入两端的虚拟单 元如和，+ 1 ，并定义相应的均值面为( 1 2 】 面o = 2 q ( t ) 一面1 ， 面v 。+ 1 = 2 卢 ) 一面。，( 2 7 ) 然后利用3 3 的斜率限制器对每个单元的数值解进行调整。 3 5数值实验 下面数值验证l d g 方法的精度。考虑具周期边界条件的纯线性抛物问题 初值问题 毗。z z ， 钍( z ，0 ) = s i n o + 5 ， z ( 一7 r ，7 r ) ，t o z ( 一7 r ，7 r ) ， ( 2 8 a ) ( 2 8 b ) 其真解为乱( z ，) = 5 + e 。s i n z 表1 给出l d g 方法在t = o 9 8 6 9 6 0 时刻 的数值误差和误差阶，其中e 表示均匀网格所含的单元个数，有限元空间为 p l j 啦和峨元，时间变量采用3 阶显式t v d r k 格式( 1 9 ) 离散。 数值分析表明，l d g 方法的误差阶与有限元空间的次数七密切相关。 求解纯抛物问题时，若为奇数，则l d g 方法仅有阶精度；而若为偶 数，则l d g 方法具有丰满的+ 1 阶误差这与f 6 1 的理论分析结果是一致 的。仔细分析= 3 与七= 2 的数值结果，我们不难发现：尽管两者在误差阶 数上没有明显区别，但前者的误差还是较后者有了明显的改善，基本达到前者 的1 1 0 换言之，高次有限元的计算精度还是较低次有限元要更为准确些。 4 多孔介质方程的l d g 数值模拟 本节给出利用l d g 方法求解多孔介质方程的数值模拟算例。针对退化多 孔介质方程的特点，本文在l d g 算法基础上提出“正值修正器”，以确保数 值计算中扩散系数保持非负，更为符合实际问题的物理背景。 1 5 1 03 0 6 9 4 1 4 e - 0 25 1 7 3 6 1 8 e - 0 463 0 6 5 4 6 e _ 0 5 2 0 14 9 2 6 4 9 争0 2 1 _ 0 4 0 l6 4 8 0 6 9 0 e - 0 52 9 9 7 07 ，9 4 1 1 1 7 e - 0 62 9 8 9 4 4 074 0 7 9 8 7 e _ 0 31 0 1 0 78 1 0 4 7 4 8 e - 0 6 2 9 9 9 3 9 ，9 4 4 0 了1 e _ 0 729 9 7 4 8 03 6 9 6 9 9 6 e _ 0 31 0 0 2 7l ，0 1 3 2 1 2 e - 0 62 9 9 9 81 2 4 3 5 5 8 e - 0 72 9 9 9 4 1 6 01 8 4 7 6 2 1 e 0 31 0 0 0 71 2 6 6 5 5 1 e - 0 73 0 0 0 01 5 5 4 6 1 8 e _ 0 82 9 9 9 8 表1 误差阶：札为单元个数 4 1 正值修正器 在求解退化多孔介质方程时，数值解保持非负可以有效地改善数值效果； 这也是问题本身的特殊要求，否则扩散系数可能为负，严重违反了实际的物理 意义。故本文认为：合理的计算格式应该尽可能地避免负值的出现。 由于l d g 方法使用间断有限元空间，每个单元可以便捷地进行事后的调 整。因此，本文借用斜率限制器的思想，引进正值修正器| p h u n ，事后修正数 值解在每个单元上的负值。其基本思想是检验每个单元内是否有负值存在， 然后在保持均值不变的前提下，将其修正为正值函数。具体做法因有限元次数 的不同而稍有区别。 1 线性有限元p 1 ：设蛳l = 吗+ ( z q ) 。由端点取值的符号，我们 可以直接看出单元内是否存在负值。若均值奶o 且端点取值非负，即 可断言在单元厶上u 2o 。因此，线性元的正值修正器p h u h 定义为 ( a ) 若吗 o ，令p h “ 1 如= o ； ( b ) 若吗o 且u h ( q 士 ) o ，则p 胁= ( 1 千掣 吗 ( c ) 其他情形不作任何调整，即p h 让 i ，j = 扎 如。 2 高阶有限元巩( 2 ) ：此时的负值检验和修正将比较困难，我们将其简 单处理为相应的线性元再进行讨论。具体做法如下： ( a ) 若吗 o ，本文统一取做e = 0 0 0 0 l 。首先，由左向右扫描初值分布，找到均值超过e 的首个单元，并断定 支集边界点位于此单元内，不妨取单元左端点作为支集边界位置在随后的时 间推进中，不断地重复上述扫描过程。当首个单元发生变更时，支集边界已转 移到相邻单元，取两单元公共节点为支集边界位置。依次记录每个时刻的支集 边界位置，我们便可得到支集边界的位置曲线。 1 8 ( a ) m = 2 f b l m ：3 ( c ) m = 5 ( d ) m = 8 图8l d g 方法求解i i 型特解：单元个数为他：3 2 0 ( a ) b a r e n b l a t t 型解：m= ( b ) i 型特解；m ：8 2 ，3 ，5 8 图9 有界支集边界点的数值追踪 1 9 4 ，2 2 算例i i 既然l d g 方法可以有效地应用于退化多孔介质方程( 5 ) 的数值模拟，我 们利用l d g 方法进行了更多的数值计算以下算例均采用p 。元和3 阶显式 t v d r k 时间离散，c f l 条件数为c f l = 0 0 8 。 算例4 1 等高的两波融合：取参数m = 5 ，相应的初值为 “牡r 裂以7 ，。0 7 - 7 3 7 ( 2 9 ) 图1 0 给出不同时刻的数值结果，其中单元个数为e = 2 2 0 。 算例4 2 不等高的两波融合：取参数m = 8 ，相应的初值为 f 1 ， z ( 一4 ，1 ) “o ( z ) = 15 ，z ( o ，3 ) ( 3 0 ) 【o ，其他 图1 1 给出不同时刻的数值结果，其中单元个数为e = 2 4 0 。 算例4 3 具有等待时间的边界扩张：取参数m = 8 ，相应的初值为 嘣垆：美哪胆，衫2 ( 3 1 ) 图1 2 给出不同时刻的数值结果，不难看到：在等待时间之前，支集边界位置 不动，参见图1 2 ( a ) 一1 2 ( d ) ；而在这个等待时间之后，支集的边界以有限速度 开始向外扩张，参见图1 2 ( e ) ， 算例4 4 带有吸收项的多孔介质方程： 鬻= 坩一砒 ( 3 2 ) 其中的p 】e 为给定的正常数。文献1 18 首先从数值计算中发现了问题真解的支 集分离现象，而后不久f 17 1 从理论上证明这种现象的存在，并给出当m + p = 2 ，0 p 蜘：。，则两侧单元被标 注为1 ，准备进行局部加密。 为避免相邻网格出现大小悬殊的情况，算法规定相邻单元的加密 等级差距不超过1 。为此，当网格标志为1 时，调用子程序c 1 1 e c k 检查 相邻单元是否符合上述规则。当相邻单元的加密等级小于此单元的加 密等级，则相邻单元需标志为1 ，然后，相邻单元也调用子程序c h e c k 检查，重复上述过程。 ( b ) 再搜索全部单元节点。若节点两侧单元是同级同源单元，且未标注为 1 ，则需判断数值跳跃是否满足1 m 2 i m s 。“。若成立，则同 级同源的两个单元将被标注为一l ，准备进行单元局部合并。 4 根据每个单元的标注，对网格进行局部加密和合并单元合并时需将两个 单元的u h 局部l 2 投影到父单元上；单元加密时只需将u h 限制到两个新 的子单元上。具体过程详见文献1 2 0 。 5 2 自适应数值算例 本节采用上述l d g 自适应策略数值求解多孔介质方程的b a r e n b l a t t 型 解( 6 ) ，i 型特解( 7 ) 和i i 型特解( 8 ) 。在所有算例中，空间均采用赐有限元 空间，时间采用3 阶显式t v d r k 方法离散，c f l 条件数为c f l = 0 0 8 ， 自适应参数均取为且玷面= 4o 和a 始。洲= 0 5 。本节所有图例，沿用以前的 表示方式，用口表示数值解在各单元的均值，用实线表示真解。数值结果表 明：上述自适应策略能够准确地捕捉突变界面的移动，而在平坦区域使用较少 的计算网格，比4 的均匀网格计算工作量有明显降低 例5 ，1 b a r e n b l a t t 型解f 6 ) 图1 5 一图1 5 ( d ) 给出自适应l d g 方法求解b a r e n b l a t t 型解在= 2 的 数值结果，其中m = 2 ，3 ，5 ，8 ，初始时刻为t = 1 。初始网格是单元个数为 e = 4 0 的均匀网格，网格加密最高次数为3 。图1 6 显示m = 8 情况下网格 随时间的变化，其中纵轴表示时间。 ( a ) m = 2 ，= 2 ( b ) m = 3 ，t = 2 ( c ) m = 5 ，t = 2 ( d ) 仇= 8 ，= 2 图1 5 自适应l d g 方法求解b 盯e n b l a t t 型解 图1 6m = 8 时的网格变化：其中纵轴表示时间 ， ：宝 。一 ， 0 1 ( a ) m = 2 ( c ) m = 5( d ) m = 8 图1 7 自适应l d g 方法求解i 型特解 例5 2 i 型特解( 7 ) 图1 7 给出自适应l d g 方法求解i 型特解在t = 0 5 的数值结果，其中 m = 2 ，3 ，5 ，8 。初始网格是单元个数为e = 2 0 的均匀网格，网格加密最高次 数为3 。 例5 3 i i 型特解f 8 ) 图1 8 给出自适应l d g 方法求解i i 型特解在= 蜀2 的数值结果，其 中m = 2 ，3 ，5 ，8 。初始网格是单元个数为e = 2 0 的均匀网格，网格加密最 高次数为4 注意，在m = 2 时自适应计算没有启动( 见图1 8 f a ) ) ，因为此 时的真解非常光滑，自适应条件没有满足。 6 结论 本文以典型的退化多孔介质方程为主要计算模型，探讨l d g 方法求解退 化抛物型方程的数值效果。大量的数值算例表明，l d g 方法可以高效地模拟 此类问题的真解，准确尖锐地追踪到真解的突变界面。同时，l d g 自适应策 略可以有效地降低计算工作量 ( a ) m = 2 f c ) m = 5 7 感谢 ( b ) m = 3 ( d ) m = 8 图1 8 自适应l d g 方法求解i i 型特解 日月如梭，光阴似箭，三年的研究生学习生涯即将结束。特别要感谢我的 导师张强老师，在我研究生三年内对我因材施教，在学业上给与悉心的指导， 并且指导我的毕业论文得以顺利完成；在生活上，也同样给与我诸多帮助。作 为良师益友，我对您的尊重感激无以言表。我还要感谢尊敬的胡健伟老师、 杨庆之老师、田春松老师，朱少红老师，由同顺老师等计算数学专业的其他老 师，感谢他们对我的指导。老师们正直的为人和一丝不苟的科研作风，值得我 尊重和学习。他们对我的教诲之恩，使我终生受益。 在学习期间，我的很多想法得益于和同学之间的交流，从他们的身上我学 到了许多东西，也得到了许多帮助，在此谨向他们表示诚挚的谢意。 “人生也有涯，学业也无涯”在南开大学就读本科和研究生这七年来， 在“允公允能，日新月异”校训的指导和南开严谨治学风气的熏陶下，我的学 业有了较大的进步，但是水平还很有限，论文还有不足之处，还请各位老师和 同学不吝赐教，惠予斧正。 参考文献 1 1s a n g e n e n t ，a n a l y t i c i t yo ft h ei n t e l f a c eo ft h ep o r o u sm e d i ae q u a t i o na 托e r w 出t i n gt i m e ，p r o c e e d z 佗驴吖t h ea m e n 竹m n t e m o 托c 以s o c i e t 玑1 9 8 8 ，1 0 2 ， 3 2 9 3 3 6 2 】d ，g a r o n s o ，r e g u l a r 姆p r o p e r t i e so f 丑o w st h r o u 曲p o r o