有限元方法课件

1、精品教学课件设计| Excellent teaching plan第1讲抛物问题有限元方法1、椭圆问题有限元方法考虑椭圆问题边值问题：(1)f x ,0,(2)问题(1)的变分形式：求uH0使满足a u, vf,v, V H0a u,v的性质，广义解的正则性结果。致剖分条件。Pvh Sh区域 的剖分，矩形剖分，三角剖分，剖分规则，正则剖分条件，拟 剖分区域上分片k次多项式构成的有限元空间ShH0。Sh的逼近性质，逆性质：u Ihu m,p Chk1m u ”， 0 m k,1 k, 1 n nm l max (0,)vh l,qChp q vh m,p 1 P,q , m l,这里，Ihu S

2、h为u的插值逼近。(3)问题(2)的有限元近似： 求uh Sh使满足a uh,vhf,vh , vhSh(3)的解唯一存在，且满足|uh|1M| f| oN(3)的解uh xui i x所满足的矩阵方程(离散方程组)形式：i 1N(4)a i, j U f, j , j 1,2, N i 1K u f刚度矩阵K a i, j N N的由单元刚度矩阵组装而成。H 1模误差分析：由(2) ( 3)可得(5a(u Uh,Vh) QVhSh由(5)可首先得到2r| u Uh 1 a u Uh,u Uh a u Uh,u I hUM u Uh 1| u IhU 1则得到u uh 1 C u Ihu1 C

3、hk u3 k 1L2 一模误差分析设w H； H 2满足w u uh, in , w 0, Iw2 Cllu uh|用u uh与此方程做内积，由(5)式和插值逼近性质得到| u uh|2 A w,u uh A u uh ,wA u uh, w I hWC u uh 1 w I hW 1Ch u uh 11| w 2Ch u uh, 1| u /再利用H1模误差估计结果，得到| u uh|Ch|u uh|1 Chk1|u|k1, k 1最优阶误差估计和超收敛估计概念。当u u t与时间相关时(为抛物问题准备)，由(5)式得a(uuh)t,Vh)0, VhSh利用(7),类似分析可得|u uh

4、t|h| uuh t|1Chk1|ut|k1, k 12、抛物问题半离散有限元方法考虑抛物型方程初边值问题：(10)- u f t,x ,u x0,u 0,x u0 x ,t,x0,Tt,x 0,Tx(10)的变分形式：求 u t ：(0,Th0u(0)u0(x)使满足(11)ut,va u, vf ,v ,v h0(11)的半离散有限元近似：求 uh t :(0,TSh使满足(12)Uh,t,Vha Uh,Vhf,Vh ,Uh(0)ShVhShN令Uh t Ui t i x ,代入(12),依次取Vhj可导出常微分方程组：i 1(13) M 号) Ku(t) f(t),0 t T其中M (

5、i, j)N N为质量矩阵，K为刚度矩阵。f (f1, ,fN)T,fjf, j。求解常微分方程组(13),得到utU1, ,UN T,代回Uh t的表达式，即得半离散有限元解 uh t。定理1.问题(12)的解Uh唯一存在且满足稳定性估计:(14) |Uht| | Uh 0| 0t|f |d , t 0证明：在(12)中取vh uh得到1 d 2rUh t a Uh,Uhf,Uh2 dt整理为(注意a u,v是正定的)ddtUht对此式积分，证毕。误差分析。引进解U的椭圆投影逼近：RhU Sh满足(15) a uRhu,vh0, vh Sh根据椭圆问题的有限元结果可知(16) |d： u U

6、h | Chk 1|Dtsu|k1, s 0,1, k 1分解误差：u uh u Rhu Rhu uh的估计由(16)式给出，只须估计Sh。由(11), (12)和(15)知， 满足t,Vh a ,Vht,Vh , Vh Sh取Vh,类似稳定性论证可得(17)0 Rhuuh 0 Rhu 0uh 0uh 0可取为u 0 u0的L2投影，插值逼近等。由（17）式，三角不等式和(16),得到(18)u t uh tuh 0Chk 1(ut1d )3、抛物问题全离散有限元近似剖分时间区间:0 t0t1tMT, tn n t引进差分算子:n tu规定，当为连续函数时,u tnn tu由此得到(19)(2

7、0)定义问题11)）将uh1Nnuii 1(22)其中Utun ut(tn)ut(t1)n _ 2I tun(u；,tnUtt s dsdtnt.Ut(tn)tntn12utttt,tn2utttUtts dstun ut(tn 2)的全离散向后ntuh,Vh0 Qu hSh ,定的。可逐层求解。Eulertntn 1Uttt s ds有限元近似：求na uh,Vhn 1,2,nf ,Vhx代入（21）可导出全离散方程组MU nn Tun)误差分析。令u tnn uhSh，使满足ShtKU n MU n 1 tF n,Fn (f1, ,fN)T,fjfuh1u(tn) Rhu tnRh1,2,

8、系数矩阵tK是对称正tnn uhno Rhu t为u的有限元椭圆投影,只须估计n Sh 。由方程（11）知Un U tn满足(23)ntU ,Vh_ na u ,Vh fnn,Vh R ,Vh ,VhShRntU tntun。则利用 aun,Vha RhUn,Vh，从（23）和（21）得到n满足n,Vhn ,VhRnn,Vh得到n)ntaRn或写为Rntn对上式求和且利用19）式得tnt0Utttn利用椭圆投影的逼近性质得到00Uh UChk 1tnUtk1dtnUtt |d再利用三角不等式即得全离散误差估计(24)U(tn) Uhn00U Uhtnt0UttChk1tnUtk1d全离散向后E

9、uler格式关于时间方向只有一阶精度。二阶精度的Crank-Nicolson格式:求unSh使满足(25)ntUh,Vh-na Uh,Vh12,Vh ,Sh0,1,其中un方程2（25）的矩阵方程形式为Un Un 1M 1,2,Un 1n 2tF 21,2,t t 1处，从方程(11)知，精确解满足 n .21 n _ nn2tU ,Vha U ,Vh(f 2,Vh)nR1n,vha(un u12,Vh)n nR tU tu(t 1)。则椭圆投影满足n _21n _"RhUn,Vha RhUn,Vh(f 2t n RnM) a(u1 n _U 2,Vh)分解误差：U(tn) Uhn

10、n nRhU(tn) Uh”。从(25)式得到1 n _t n Rn，Vha(Unu 2,Vh)(26)"t n,Vha-n,Vh在(26)中取Vh- n,注意六 n2 n12(27)1 n _ n2u us1 Utt s - t2d dst Utttn 1s ds11nn -a(Un u 2,Vh) C|un u 2112 vh则在(26)中取vh n得到1tt n|R1n|C|Un u |2求和，且利用(20), (27)和椭圆投影的逼近性质，得到1nl | 0| Chk1：|ut|dC t2otn HuttL |uttt|d再利用三角不等式，得Crank-Nicolson格式的

11、误差估计：(28)|u(tn) Uhnl Ct2 hk 1 n0,1,第2讲有限体积元方法椭圆问题考虑椭圆边值问题：u f, inu 0, on设R2为凸多边形区域。当 f L2 时，u H0H 2唯一存在，且满1剖分与对偶剖分是内节设Jh是 的一个正则三角剖分， Nh是所有剖分节点集合，N： Nh点集合，对每一个 pi Nh做包含Pi的有限体积Vi （见图1）.,1其中q - Pi Pj Pk是单3元K Pi, Pj, Pk的重心，1mj 2 Pi Pj 是单兀边 p?pj的中点，所有有限体积单元Vi构成区域 的一个剖分，称为Jh的对偶剖分，记为* 一Jh ,则K VK JhV Jh*图12

12、有限体积元空间1*设Sh H0是Jh上分片线性多项式构成的有限元空间，在对偶剖分Jh上，定义分片常数的有限体积元空间:*.Sh vL2vV constant, V Jh设i是分片线性基函数,i是Vi的特征函数，则Shspan i x : pi0Nh0*Nh ; Sh span i x : Pi插值算子定义插值算子:Ih:C-*Ih :CSh, IhuPi*Sh , Ih uuN0PiPiuN0Pi i x显然I hU PiPi ,PiPi ,0PiNh因此,*对任息vh*Sh存在唯一 VhSh*IhvhVh,即 Sh* 一h Sh双线性形式由 u, vK Jhu,V Ku,v则定义与u,v相应

13、双线性形式:u,v A u, vVivdu,vh0*JhVids,u, v*Sh Sh问题的有限体积元近似Auh,Vhf,Vh,*Vh*Sh(1)由于Sh IhSh ,则(1)等价于A u h, I h vhf , I hvh ,则有误差方程*A u uh, I hVh插值算子Ih的性质i).*u Ihu0, PCh u.p,ii)*Vh IhVhdK 0,Kiii)IhvhI；Vh ,0,Vh守恒性质：在(1)中取vh 得到uhdsVif,Vi Jh7解uh的存在唯一性vhVhShSh * .一IhVh dl0,VhSh(2)Jh,lK*引理1 设 Uh,VhSh,则 AUh,VhAUh,I

14、hVh证明：利用Green公式和（3）式，注意uh为常数，得到 n*,*Uh , 1 hvh KUh , 1h vh Vj K*JK Jh Vj JhK Jhh *一IhvhdsnVj JhVjUh . *Ihvhds n-h-vhds AUh,Ihvhk Jh K n*uh, vh k AUhJhvhAuh,vhAUhJhvhK Jh定理1有限体积元解uhSh唯一存在，且满足|uh' C| f|证明：由引理1,方程(1)和(4)式知uhSh满足C|uh|：| uh|2Auh,uhA(uh,I；uh) (f,I；uh)IlflllIhuhl2|f|uh|2 flML8 H1模误差分析定

15、理2 u和uh满足误差估计u uh 1 Ch f证明：利用引理1得到uhC u uh 12 A u uh ,u uh A u uh ,u Ihu A u uh, I huA uuh,uI huf,IhuuhA uh, Ihuuh*Auuh,uIhuf,IhuuhA(uh,Ih(Ihuuh)*Auuh,uIhu(f, (Ihuuh)Ihguuh)C uuh 1 u Ihuh 1 ChfIhuuh 1Ch u uh 1 u 2 Ch f|( Ihu u 1 u uh 1)Ch|u uhlJIfll Ch2|f|2 Ch f Illi u uhl利用柯西不等式，定理 2得证。9 L2 模误差分析定理

16、3 u和uh满足误差估计：|u uhl Ch2| f|1w 0, IW 2 C| u Uh|A u uh,w证明：设w H 0 H 2满足w u Uh, in ,则利用引理1和(3)式得A u uh,w Ihw A u uh,w Ihw A u uh,w Ihw A u uh,w Ihw|u Uh|2 Aw,u UhA u uh ,I hw*f ,I hwA(uh, I hI hw)*.f,Ihw IhIhw*.f fK,Ihw IhIhw其中fK1 f, x K,K Jh,是f的分片常数逼近。则利用逼近性质得到 K KhII u uh|2 I u uh|1|w %w111f fK|Ihw I

17、hIhwCh2llfllliwi2时旧1|岛1 Ch2llfll1lwl2 Ch2ll fllJu uh|证毕。抛物问题u r(5)-uft,x,t,x0,Tu0 ,t,x0,T0u 0, xu x ,x2*_ *其中 R为凸多边形区域。设区域剖分Jh和对偶剖分Jh ,有限体积元空间Sh和Sh如 椭圆问题情形。抛物问题(5)的有限体积元近似为：求 uh t : 0,TSh使满足*(tuh,IhVh) A(uh,IhVh)(f,IhVh),Vh Sh”、(6)uh(0) Sh其中双线性形式Au,v如椭圆问题。引理2下述结论成立*i) (uh, IhVh) (Ihuh,Vh), uh,Vh Shi

18、i) Nh|hJ(uh,I；uh)是Sh上与等价的范数。由引理2可知，tuh,IhVh导出的质量矩阵是对称正定的，则常微分方程组(6)*唯一可解。由引理1也知道由A uh, IhVhA uh,Vh导出的刚度矩阵也是对称正定的，这可保证全离散格式的唯一可解。误差分析 设u为问题（6）的解,引进u的有限体积元投影：Rhu Sh满足*A(u Rhu,I hvh)0,vhSh对（7）关于t求导得到*A(utRhut, I hvh) 0,VhSh(8)则利用椭圆问题有限体积元的结论可得（注意u f ut)|uRhU|Ch2(|f| |ut|i)(9)uRhU tCh2( ut 2ftUtt i)(10)

19、分解误差：u uh u Rhu Rhu uh的估计已知，只须估计Sh o由方程（5） （ 7）可知,满足取vh由（4）,利用引理1 2得到2dt式和引理2知vhSht IhIh 2 2Ch,则从上式得到0 h C 0tdtt C( 00 tdt)则利用三角不等式和（10）式得到u t uh t I |u Rhu Rhu uhC u 0 uhCh2 f 0 1 ut 0 1 f 1 ut(ut2ft1utt|L)d(1):（1, a）T ,代表流体时刻，任一初值点x 0 x0都对应一条特征线x t atx0。第3讲、一阶线性双曲方程的差分方法1、一阶双曲问题a 0,0 tt xu(0,x) x一

20、阶双曲方程的一个重要概念是特征。方程（1）的特征方向是dx t流动的万向。牛寸征万程te： a,对应的特征曲线为：x t atdt函数u t, x沿特征线的微分为:du udt tu dx u也可理解为沿特征方向的导数:dt uu axua x由方程（1）得知，沿特征线du * * t, x tdtu t, x t C利用初值条件u0, x 0 u 0, x0xo可知，沿从xO出发的特征线x t at x0, (1)的解u恒为常数u t,x txo ,从而得到（1）的解。2、间断性质（激波）当 x有间断时，随着t的发展，u保持间断，若u函数值间断，称为强间断，若u连续，但导数间断，称为 弱间断

21、。例1强间断，给定初值0,x0x1,x0解见图1。例2 弱间断，给定初值1, x 1x - 1 x ,1 x 120,1 x解见图2。了 二缁十I图2 弱间断解3、差分格式在x方向和时间t方向进行剖分，剖分节点:tnn t, n 0,1,xjjh,j 0,1,t,h为剖分步长。可对方程（1）构造如下三种差分格式n 1 nUj UjtnnUj 1Ujah0,(2)n 1 nUj Ujtn nUj Uj 1ah0,(3)n 1 nUj Ujtn nUj 1 Uj 1a2h0,h2(4)稳定性分析 采用Fourier方法，令nuj为参数。记网比r将U；表示代入（2） （4）,可分别导出增长因子:Gi

22、i harearar 4arsin2G2G3GiG2G3则在ari harear ,G222 har 4ar sin 21,iar sinh,G399.9a r sin h0,0,a,r1条件下，当是不稳定的。（2）和（3）arar0时，采用格式（称为迎风差分格式n 1 nujujan nr uj uj 1 a其中 a max( a,0), a网比条件依赖域，见下图。2）,当a 0时，采用格式（3）。格式（4）,可统一写为:nnr uj 1 uj 0min（a,0）。迎风差分格式（2） -（3）是依据特征线走向取值的,（特征线的斜率）是为了保证差分解的依赖域包含精确解的t“一特征线/ / /如

23、一1TXa >0 ./-I格式（4）是不稳定的，可以改造为:n 11 nn U j 2(U j 1 U j 1 )tnnUj 1 Uj 1 八a0,2hR 0(h21)(5)称为Lax-Friedrichs 格式，稳定性条件为上述格式都是一阶的，一个二阶格式是Lax-Wendroff格式：n 1 nnnUj Uj Uj 1 Uj 1at2hnUj 1 ,O t2h2(6)可以看作是对原方程的粘性近似方程:u ax的逼近，称为粘性差分格式,右端项称为人工粘性项。稳定性条件是事实上，迎风格式和LF格式都可以看作为粘性差分格式，它们可分别写为n 1 nUj Ujtn nUj 1 Uj 1 a2

24、ha2hnnUj 1 2UjnUj 1迎风格式n 1 nUj UjtnnUj 1 Uj 1a2hnnUj 1 2UjnUj 1LF 格式人工粘性项分别为:1 ah21 h22 t x2对变系数at,x情形，差分格式同样改造，稳定性分析时，先冻结系数n aja,可导出同样的稳定性条件:ht 1对有界区域的边值问题, 左边界给定边值。要依特征走向给定边界条件。例如，当a 1时，只能在区域第4讲 非线性守恒律方程的差分方法1、非线性守恒律方程0 t,(1)u f u c 0,t xu(Q x) U0 x其中f （u）是u的非线性函数，称为 通量函数，问题（1）的光滑解可能不存在，但可以存在弱解：u

25、f dtdx 0, 0 t xC。 0,满足（2）式的弱解存在也不一定唯一。如果在间断线x （t）上，弱解还满足 嫡条件：f u f u f u f u f u f uu uu uu u(3)则称其为可容许解或物理解，物理解是唯一的。（1）的解通常存在间断（激波），可以是强 间断或弱间断。2、差分离散在网格节点tn,xjn t, jh处，方程（1）可离散为n 1ujnuj tf(un1)f(un1)j jr2- R , R是截断误差tun1 un r(f(u? f" tR, rTt舍去截断误差项，且对 f（un 1）采用某种形式的网格节点近似，可得到守恒型差分格式: j2其中?n1j

26、 2un1 un r(Pn1?n1)j 2 j 2?un,un 1称为数值通量，要求 ？满足(4)f（u,u） f （u）（5）称（5）为相容性条件。3、相容性与守恒性当差分方程（4）满足相容性条件（5）时，它与微分方程是相容。事实上，对任何光滑函数u和f（u）,利用泰勒展开得到（记 ? ？a,b ）? Uj,Uj 1f?Uj ,UjUj 1 Uj ?b Uj,Uj.2O Uj 1 Ujf ujhuj ?)uj ,uj O h2? uj 1 ,ujf uj huj f uj ,ujO h2则（4）式为n 1 n .0 n nujuj rhuj fa uj ,uj? n n2fb u j ,u

27、jrO hnd ? n nuj t ? uj ,uj O h t jdx或者写为nujnuj令t,h 0 ,即知差分方程逼近于微分方程。称差分方程为守恒型，如果nUjnUj(6)对相容的差分格式（4）,如果un 0,当j容性条件知时，则它为守恒型的。事实上，利用相P n n? n nf-uj ,uj 1? uj 1,djj 0Jim ? un, un 1 f?un 1, un 0对（4）式求和，即知（6）式成立。4、几个典型的差分格式例 1 Lax Friedrichs 格式n 1Ujnnuj 1 uj 12nnf 5 1 f 5 12h可改写为守恒型:n 1 n ?、ujuj r( ? 1f

28、- 1)j 2 j 2? 1f-uj ,uj 11 f uj 1f Uj12r uj 1 uj例2迎风格式（注意f u f u ）n 1 nUj Ujnf Ujnf Ujnnf Uj f Uj 1hnnf Uj 1 f Ujhnnf Ujf Uj 1n nUjUj 1n nf Uj 1 f UjnnUj 1Uj可改写为守恒型格式Un 1Un r( ? 1 f? 1 )j j 221 -fj 1 f?Uj ,Uj 1- f Uj 1j 22f Uj 1 f Uja. 1j 2 Uj 1 Uj例 3 EngqUist-Osher 格式Uj21a1 I Uj2Ujn 1 n . . ? UjUjr(

29、 f. 1 f 1 )j - j 22f? 1 f?Uj ,Uj 1 f Uj f Uj 1j2Umin f0s ,0 dsf u o max f s ,0 ds f(0), f u5、稳定性方程(1)可写为：那么线性方程稳定性准则:ar 1 ,可转化为max fur 1 CFL 条件U但对非线性方程，CFL条件只是稳定性的必要条件，并不充分。一般还得对差分格式附加其 它要求才能保证非线性稳定性。CFL条件的一个离散形式为f Uj 1 f Ujr IIUj 1 Uj1,定义 Uj Uj 1 Uj, Uj Uj Uj 1 ,变差若差分解满足TV(u) Ujj 0TV(un1) TV(un),n(

30、8)则称差分解是TVD的（总变差衰减）TVD格式通常保持差分解的单调性，且为高分辨格式。上述例中的三个格式在 CFL条件下都是守恒型 TVD格式。但都为一阶精度，如何构造高精度的TVD是一个难点。高精度是指在解 U的光滑区域上。将守恒型差分格式写为nUj1n n n z nH Uj i,Uj,Uj i ( Ujr(?j12? i)j -(9)如果H Uj i，Uj,Uj 1对三个参数都是单调递增白1则称差分格式（9）是单调差分格单调差分格式是L1-TVD的,且满足极值原理:inf0 ujn uj0supujj将一般的差分格式改写为n 1uj 1 ajbjnujajnuj 1nuj uj定理1如

31、果差分格式（10）的系数满足aj0,bj0,1 ajbj 0,则差分格式为单调的；如果满足aj 0,bj 0,1ajbj 1 Q(10)(11)(12)则差分格式是TVD的。证明：第一个结论是显然的,10）写为则可有（注意 ujn 1ujnujnujajnujbjnujuj)n ujajnuj 1ajbjajn ujn ujajj 1nlj 1nuj 1bjbjn ujn uj利用（12）式得到nujajbj 1nujajnujbjnuj求和得到n 1TV uajbjn ujaj0n ujbj0n ujajbjn ujaj j 0n ujbj 1n ujn ujTV上述三种差分格式在 CFL条

32、件下都是单调的 TVD格式。第5讲一阶双曲方程间断有限元方法维线性问题u a0,(1)特征方向：1,a，特征方程 一 a,特征曲线x at c。对问题（1）可以考虑纯初值 dt问题，也可以考虑初边值问题，但边值条件要根据特征走向给定。这里考虑初边值问题： 0 x l。假设a 0,在x 0左边界给定边值条件：u t,0 g t 。空间剖分：取定剖分步长 h ,剖分节点Xj jh, j0,1,N,h l/N ,半节点.1.x 1(j -)h ,单兀 I j x 1, xj 22j j 2 j1, j21,2,N 1,I0O，Xi ,In2x1,xN。2(2)Vh中函数可以是不连续的。用 Vh Vh

33、与（1）j上的内积，得至kIjVhdx au(x 1)Vh(x 1) u(xtj 2)Vh(x 1) a u j 2IjVh .dx x(3)对于间断点x 1 j2极限）vh （x 1 ） j 2上的取值，Vh（xj1）恒取为在单元2I j内侧的极限值,即 vh (x征走向取值。因u .VhdxIj t1 （右极限）。对于方程解2由连续性知u 1j 2u 1 ,则一般根据特 j 20,则取迎风值u(x 1)j 2u 1 o离散方程（3） j 2现在变为a u (x)Vhj2(Xjl)2(x. 1)Vhj2(x 1) aj 2uIjVh dx 0x这是精确解满足的变分方程，在其中用uhVh代替则

34、得到方程1)的间断有限元近定义分段r次多项式构成的间断有限元空间:VhVh L2 0,l : Vh |IjPr Ij ,似：求uh t Vh使满足Ij?Vhdx auh(x 1)Vh (x 1)tj 2j 2uh(x. 1)Vhj 2(xjl)2IjuhVhdx 0(4)VhVh,1, ,N引进单元I j上的基函数空间:Vh Ijj1 , j2 , jrjiPr Ij , ji |Ik0, j(5)则Vh为VhIjN 1 ,、一 一，j 1的直交和。可在（4）中令 Vh Vh I j , j 1,2, ,N 1。在单元 I j 上,令Uht,x1Uji1t ji (x),x I j,U j1,

35、 ,Ujr1 T ,将Uh代入（4）并依次取Vhj1,可得到jkdU dxIjt a dt1Uji t1ji(xjjk(xjU j 1i t j1i(xj1 )-2jk(xj1)-2jk dxU ji t0 ,k1,2,r(6)可将离散方程（6）整理为矩阵方程形式:dU ja AU j BU j 1 dtaCU j0,1,2,N其中M , A, B,C为r 1阶矩阵，元素分别为mik I ji jk dx,aikji(xj1)2jk(xj1)2bik j 1i (x 1 ) jk (x 1), jj -22cikIjji jk dx,i,k1,2,r例设Vh为分片常数空间，每个单元I j只有个

36、基函数j(7)dujh a u j u j 1出j j0, j 1,2, ,N(8)对t离散后，这是迎风差分格式，u0由边值确定。对方程（7）可进一步进行时间离散，例如，采用一阶向后差商，得到Un U n1 j jnnM a(A C)U j aBUjj 1,2, ,N, n 0,1,(9)这是全离散间断有限元格式，可从左至右逐单元求解Un （对固定n ）,时间方向逐层求解。二、非线性守恒律问题u f u c0,t x对方程（10）积分得到(10)7Vhf(u(x)Vh(xIj tj 2jf(u(xj)Vh(x. 1) j 2f u vh 0Ijh(11)类似线性情形，?(Uh(x 1), Uh

37、 (xj 2jUhIj t Vhf Uh Vh ?(Uh (x 1)Ijj 2,Uh (x 1 )Vhj 2(x. 1)j 2(13)?(Uh (x 1 ), Uh(x 1 )Vhj -j -(x. 1) 0,j 2VhVh由相容性条件知，问题的连续解为了使格式保持相容、稳定,U也满足(13)。TVD和高分辨性质,通常借助差分方法的思想构造数值通量函数P。例 1 Lax Friedrichs 型?LF a,b -fa2c max f0? s t?s , a?min a,b , bmax a,b例 2EngqUist Qsher 型bb? a, b min f10as ,0 ds max10f

38、s ,0 ds对分段常数间断有限元，即Uh Uj, xIj13)转化为dU j h dt?(Uj,Uj 1)?(Uj 1,Uj)再对时间离散，得到守恒型差分格式:n 1nt ?Uj Uj f(Uj,Uj1) fUj1,Uj), h为了保持格式的稳定性，还需附加CFL条件：max f u 1现考虑半离散方程(13)的求解。利用基函数可将(j 1,2,(14)13)写成常微分方程组形式：dU tdtLh U t(15)在间断点x 1处，需特殊定义f (u(x 1 )的值，记作 j 2j 21),称？为数值通量函数，要求 ？是相容的，即2? U, U f U(12)在(11)式中用Uh Vh代替u

39、, ?代替f ,则得到间断有限元格式：求uh Vh使满足初值U 0由u 0确定，一般由u 0 U° x的L2投影确定。对常微分方程组(15)可采用k阶显式Roung-Kutta方法离散，由时间方向逐层求解。例如，可采用如下 R-K算法：1 U (0) U n2 .对i 1,2, ,k计算i1i1U(i)ilU(l) t il Lh(U(l)l0l03 U n 1 U (k)为了保持非线性稳定性，避免非物理振荡，通常还要采用坡度限制器( slope limiting )算子进行处理，即在上述步2 中，令i1Uh h ii Uh(l) t ii Lh(Uh(l) , h为坡度限制器算子

40、l0对多维问题或方程组情形，间断有限元方法更复杂，最优收敛阶，超收敛，后处理等研究都还不完善。第六讲 对流占优问题的特征差分和有限元方法1 .对流占优问题uv u b u f t,x , t,x（0,Tut,x0, t,x（0,T（1）u 0, x u0 x , xv 0为扩散系数，bn，b2 ，当b v时，称（1）为对流占优问题，此时方程（1）具有双曲性质，一阶项占主要地位。这样方程的解在边界附近急速变化，称为边界层现象。 用通常的数值方法求解，数值解在边界层区域将出现剧烈振荡，从而失真。2 .方程的特征形式方程（1）双曲部分的特征方程和特征曲线是（假设4,b2为常数）dx1dx2b1，b2

41、, Xi tdtdt特征方向是b1t x1 0 , x2 tb2tx2 01 1,b1,b2 T ,|1,b12,b22 2则沿的方向导数为 t, Xi, X2- b1b2 tx1x2那么方程(1)沿特征方向可写为ufv u f(3)3.沿特征方向的离散剖分 0,T :0 t0 t1tM T,tjj t, t T/M 。由特征方程可知，在t tn时刻经过点x （x；, x2n）的特征线为:Xi t“t (X1n b1tn),X2 tb2t(X2nb2tn)(x? b1 t,xn b2 t)。那么在点则这两条特征线与t tn1平面的交点为 xx1,x2tn,x（tn,x1n,x2n）沿特征方向的

42、导数可近似为这里利用-tn,xU tn,x U tn 1, x1 U tn,x U tn 1,x O tt22tb1b21q。则方程(3)在特征方向的离散为u tn ,xu tn 1, x tV U tn,Xtn,xO t需进一步在空间区域离散方程(4)。A一维示意图4.特征差分方法设 为矩形区域，对 做矩形网格剖分，剖分节点记为(xi, yj ) ,在每个节点(xyj )处，对方程(4)进行空间离散，得到nn 1nn nUijU ijUi 1j 2uij Ui 1jv 2t,2nnUij 1 2uijh2nUij 1nfij(5)其中 Ujn1 是 tni,x U tn 1, xi bi t

43、, yj b2t的某种近似,hi,h2分别为x和y方向剖分步长。由于un1一般不是节点值，需特殊计算。取t 使 b1th1,b2 th2,这样可使点xi bi t,yj b2 t处于与节点 xyj相邻的网格区域内,见图Af(3)则可用X所在的网格单元四个顶点函数值的算术平均或双线性插值来近似ujn 1。这相当一种迎风格式，例如，当b1 0,b2 0时，是用 为，yj点左下侧节点值近似该点值。在一维情形，当b1t h时，Xxib1 t落在xi1, xib10或xi1, xb10区间内，可用x相邻的两个节点值近似计算 Uin 1。5.特征有限元方法对区域 剖分，定义有限元空间 Sh,则方程（4）的特征有限元离散为：求 un Sh使 满足 nnn 1nUh ,Vhva Uh ,VhUh x , Vht f ,Vh , Vh Sh”、0（6）u0 Shn 1,2,在计算积分un 1 x,vh时，要注意xx1b1t,x2b2t。由于un1 x是分段或分片定义的（基函数也是如此），需要知道x随x变化时所在的单元， 才能知道基函数 i x 的表达式。因此需要特殊方法计算 un1 x ,Vh