有限元方法课件
- 格式:doc
- 大小:1.58 MB
- 文档页数:31
有限元⽅法(课件)第⼀章有限元概貌与发展有限元⽅法是近似求解数理边值问题的⼀种数值技术。
这种⽅法⼤约有60年的历史。
它⾸先在本世纪40年代被提出,在50年代开始⽤于飞机设计。
后来,该⽅法得到了发展并被⾮常⼴泛地⽤于结构分析问题中。
⽬前,作为⼴泛应⽤于⼯程和数学问题的⼀种通⽤⽅法,有限元已相当著名。
有限元法应⽤于电磁场中,最先是⽤结点上的插值基函数来表征该结点上的⽮量电场或磁场分量的,称为结点有限元。
但是,在使⽤结点有限元进⾏电磁仿真时,会有⼏个严重的问题。
⾸先,⾮物理的或所谓伪解可能会出现。
其次,在材料界⾯和导体表⾯强加边界条件很不⽅便。
再次,处理导体和介质边缘及⾓也很困难,这是由与这些结构相关的场的奇异性造成的。
在这些问题中,最后⼀个问题⽐其它两个问题更严重,因为它缺少通⽤的处理⽅法。
即使对前两个问题,⽬前的处理状况也不能完全令⼈满意。
因此,有必要探讨其它的可能性或其它⽅法,⽽不仅仅是改进,从⽽将电磁场有限元分析引⼊⼀个新的时代。
幸运的是,⼀种崭新的⽅法已经被发现。
这种⽅法使⽤所谓⽮量基或⽮量元,它将⾃由度(未知量)赋予棱边⽽不是单元结点。
因为这个原因,它也叫棱边元(edge element )。
虽然Whitney 早在35年前就描述过这些类型的单元,但它们在电磁学中的应⽤及其重要性直到前⼏年才被认识到。
在80年代初,Nedelec 讨论了四⾯体和矩形块棱边元的构造。
Bossavit 和Verite 将四⾯体棱边元应⽤于三维涡流问题。
Hano 独⽴地导出了矩形棱边元,并⽤于介质加载波导的分析。
Mur 和de Hoop 考虑了⾮均匀媒质中的电磁场问题。
Van Welij 和Kameari 应⽤六⾯体棱边元进⼀步考虑了棱边元在涡流计算中的应⽤。
Barton 和Cendes 将四⾯体棱边元应⽤于三维磁场计算,同时,Crowley 提出了⼀种更复杂的单元类型,即所谓的协变(covariant )投影单元,它允许单元带有弯曲的棱边。
第1讲 抛物问题有限元方法1、椭圆问题有限元方法考虑椭圆问题边值问题:(1) ()⎩⎨⎧Ω∂∈=Ω∈=∆-x u x x f u ,0,问题(1)的变分形式:求()Ω∈10H u 使满足(2) ()()()Ω∈∀=1,,,H v v f v u a ()v u a ,的性质,广义解的正则性结果。
区域Ω的剖分,矩形剖分,三角剖分,剖分规则,正则剖分条件,拟一致剖分条件。
剖分区域上分片k 次多项式构成的有限元空间()Ω⊂10H S h 。
h S 的逼近性质,逆性质:∞≤≤≤≤≤≤-+-+p k k m uCh uI u pk m k pm h 1,1,0,,11,h h pm hqn p n l m ql hS v l m q p v Chv ∈∀≤∞≤≤≤---,,,1,),0(max ,这里,h h S u I ∈为u 的插值逼近。
问题(2)的有限元近似:求h h S u ∈使满足 (3) ()()h h h h h S v v f v u a ∈∀=,,,(3)的解唯一存在,且满足f M u h ≤1。
(3)的解()()∑==Ni i i h x u x u 1φ所满足的矩阵方程(离散方程组)形式:()()N j f u a jNi ijiΛ,2,1,,,1==∑=φφφ(4) f u K ϖϖ=刚度矩阵()()NN ji a K ⨯=φφ,的由单元刚度矩阵组装而成。
-1H 模误差分析:由(2)-(3)可得(5) h h h h S v v u u a ∈∀=-,0),(由(5)可首先得到()()1121,,u I u u u M u I u u u a u u u u a u u r h hh h h h h--≤--=--≤-则得到(6) 1,111≥≤-≤-+k uCh u I u C u u k k h h2L -模误差分析设210H H w I ∈ 满足h h u u C w win u u w -≤=Ω-=∆-Ω∂2,0,,用h u u -与此方程做内积,由(5)式和插值逼近性质得到()()w u u A u u w A u u h h h,,2-=-=-()hhhh h h h u u u u Ch wu u Ch w I w u u C w I w u u A --≤-≤--≤--=12111,再利用-1H 模误差估计结果,得到 (7) 1,111≥≤-≤-++k uCh u u Ch u u k k hh最优阶误差估计和超收敛估计概念。
当()t u u =与时间相关时(为抛物问题准备),由(5)式得 (8) h h h t h S v v u u a ∈∀=-,0),)((利用(7),类似分析可得 (9)()()1,111≥≤-+-++k u Ch u u h u u k tk th t h2、抛物问题半离散有限元方法考虑抛物型方程初边值问题:(10) ()()()()()()()()⎪⎪⎩⎪⎪⎨⎧Ω∈=Ω∂⨯∈=Ω⨯∈=∆-∂∂x x u x u T x t x u T x t x t f u t u,,0,0,,0,0,,,0(10)的变分形式:求 ()())()0(,],0(:010x u u H T t u =Ω→ 使满足(11) ()()()()Ω∈∀=+1,,,,H v v f v u a v u t(11)的半离散有限元近似:求 ()h h S T t u →],0(:使满足 (12) ()()()⎩⎨⎧∈∈∀=+hh hh h h h h t h S u S v v f v u a v u )0(,,,,,令()()()∑==Ni ii h x t u t u 1φ,代入(12),依次取j hv φ=可导出常微分方程组:(13) T t t f t u K dtt u d M ≤<=+0),()()(ρρρ其中N N j i M ⨯=)),((φφ为质量矩阵,K 为刚度矩阵。
()j j T N f f f f f φ,,),,(1==Λρ。
求解常微分方程组(13),得到()(),,,1TN u u t u Λρ=代回()t u h 的表达式,即得半离散有限元解()t u h 。
定理1. 问题(12)的解h u 唯一存在且满足稳定性估计:(14) ()()()0,00>+≤⎰t d f u t u th h ττ证明:在(12)中取h h u v =得到()()()h h h h u f u u a t u dtd ,,212=+ 整理为(注意()v u a ,是正定的)()f t u dtdh ≤ 对此式积分,证毕。
误差分析。
引进解u 的椭圆投影逼近:h h S u R ∈满足 (15) ()h h h h S v v u R u a ∈∀=-,0,根据椭圆问题的有限元结果可知 (16) ()1,1,0,11≥=≤-++k s uD Chu u D k s t k h st分解误差:θη+=-+-=-h h h h u u R u R u u uη的估计由(16)式给出,只须估计h S ∈θ。
由(11),(12)和(15)知,θ满足()()()h h h t h h t S v v v a v ∈∀-=+,,,,ηθθ取θ=h v ,类似稳定性论证可得(17) ()()()⎰+≤tt d t 00ττηθθ()()()()()()()0000000h h h h u u u u R u u R -+-=-=θ()0h u 可取为()()x u u 00=的2L 投影,插值逼近等。
由(17)式,三角不等式和(16),得到(18) ()()()()()())0(000111⎰+++++-≤-tk t k k h h d u u Chu u t u t u ττ3、抛物问题全离散有限元近似剖分时间区间:M T t t n t T t t t n M /,,010=∆∆==<<<=Λ。
引进差分算子:tu u u n n nt ∆-=∂-1规定,当()t u 为连续函数时,()n nt u u =,则有()⎰⎰-∆-=-∂n n nt t t tt n t nt dsd s u t t u u 11)(ττ()()()()⎥⎥⎦⎤⎢⎢⎣⎡-+-∆=-∂⎰⎰-----n n n n t t ttt n t t ttt n n t nt d u t d u t t t u u 212112212121)(ττττττ 由此得到(19) ()()⎰-∆=≤-∂nn t t tt n t nt t O ds s u t u u 1)((20) ()()()22112)(t O ds s u t tu u n n t t ttt n t nt ∆=∆≤-∂⎰--定义问题(11)的全离散向后Euler 有限元近似:求h nh S u ∈,使满足(21) ()()()⎩⎨⎧=∈∈∀=+∂Λ,2,1,,,,,0n S u S v v f v u a v u h h hh h n h nh h n h t将 ()∑==Ni i n i nhx u u 1φ代入(21)可导出全离散方程组(22) Λ,2,1,1=∆+=∆+-n tF MU tKUMU n n nn其中()j nj TN nTnN nnf f f f F u u U φ,,),,(,),(11===ΛΛ。
系数矩阵tK M ∆+是对称正定的。
可逐层求解。
误差分析。
令()()()nn n h n h n h n n h n u t R t u R t u u t u θη-=-+-=-)(。
()t u R h 为u 的有限元椭圆投影,只须估计h n S ∈θ。
由方程(11)知()n nt u u = 满足(23) ()()()()h h h n h n h n h n t S v v R v f v u a v u ∈∀+=+∂,,,,,()nt n t nu t u R ∂-∂=ρ。
则利用()()h n h h n v u R a v u a ,,=,从(23)和(21)得到n θ满足()()()h n t n h n h ntv R v a v ,,,ηθθ∂-=+∂取nh v θ=得到()()nnt n nn nn n R t ta θηθθθθθ∂+∆≤∆+--ρ,),(12或写为()n t n n n R t ηθθ∂+∆+≤-1写 ()⎰-∆=∂nn t t t nt d t 11ττηη 对上式求和且利用(19)式得()()⎰⎰+∆+≤nn t t t tt n d d u t 00ττηττθθ利用椭圆投影的逼近性质得到()()⎰⎰∆+⎪⎭⎫ ⎝⎛++-≤+++nnt tt t k t k k h n d u t d u uCh u u 00110100ττττθ 再利用三角不等式即得全离散误差估计(24)()()Λ,1,0,)(0110100=⎪⎭⎫ ⎝⎛++∆+-≤-⎰⎰+++n d u uCh d u t u u u t u nnt k t k k t tt h nhn ττττ全离散向后Euler 格式关于时间方向只有一阶精度。
二阶精度的Crank-Nicolson 格式:求h nh S u ∈使满足(25) ()()⎪⎩⎪⎨⎧=∈∈⎪⎪⎭⎫ ⎝⎛=+∂-Λ,2,1,,,,,021n S u S v v f v u a v u h hhh h n h nh h n h t其中 ()121-+=n nnu u u 。
方程(25)的矩阵方程形式为21112---=⎪⎪⎭⎫⎝⎛++⎪⎪⎭⎫ ⎝⎛∆-n n n n n F U U K t U U M K ,2,1,22211=∆+⎪⎭⎫ ⎝⎛∆-=⎪⎭⎫ ⎝⎛∆+--n tF U K t M U K t M n n n21-=n tt 处,从方程(11)知,精确解满足()()()),(,),(,,21121h n nh n h n h nh ntv uu a v R v fv u a v u---++=+∂)(211-∂-∂=n t n t n tu u R 。
则椭圆投影满足()()),(),(,,21121h n nh n nt n h nh h nhtv uu a v R fv u R a v uR ---++∂-=+∂η分解误差:n h n h n n n nh n u t u R u t u -=+=-)(,)(θθη。
从(25)式得到(26) ()()()),(,,,211h n n h n n t h n h n t v u u a v R v a v --++∂-=+∂ηθθ在(26)中取nh v θ=,注意()()()112122121,---+-∆=⎪⎭⎫ ⎝⎛-∆=∂n nn n n n nnttt θθθθθθθθ(27) ()()⎰⎰⎰--∆≤=-∆--nn nn t t tt s t s tt t tn nds s u t ds d u uu 1212121ττh n nh n nv u u C v uu a 22121||||),(---≤-则在(26)中取nh v θ= 得到⎪⎪⎭⎫ ⎝⎛-++∂∆≤---22111||||n n n n t n n u u C R t ηθθ求和,且利用(20),(27)和椭圆投影的逼近性质,得到()⎰⎰+∆++≤+nn t ttt tt t t k n d u ut C d u Ch 022010ττθθ再利用三角不等式,得Crank-Nicolson 格式的误差估计: (28) ()Λ,1,0)(12=+∆≤-+n h t C u t u k nhn第2讲 有限体积元方法一、 椭圆问题考虑椭圆边值问题:⎩⎨⎧Ω∂=Ω=∆-on u in f u ,0,设2R ⊂Ω为凸多边形区域。