第五章有限元素方法 51有限元素方法的基本思想 特点: (1)能处理复杂区域和复杂边界条件的求解问题。 (2)有限元素法是一套求解微分方程的系统化数值计算方法。 它比传统解法具有理论完整可靠,物理意义直观明确,解题效能强等优点。 (3)对连续体的问题采用有限元素法,是将连续问题离散化的数值求解方法。 应用范围: 由于这种方法适应性强,形式单纯、规范,因而自50年代以来,在计算机的配合下, 有限元素法在物理和工程设计计算的许多领域得到了广泛的应用。该方法不仅适用于 电磁场问题的求解,也是对其它具有复杂边值问题的数学物理方程求解时的高效能的 方法
第五章 有限元素方法 5.1 有限元素方法的基本思想 特点： （1） 能处理复杂区域和复杂边界条件的求解问题。 （2） 有限元素法是一套求解微分方程的系统化数值计算方法。 它比传统解法具有理论完整可靠，物理意义直观明确，解题效能强等优点。 （3） 对连续体的问题采用有限元素法, 是将连续问题离散化的数值求解方法。 应用范围： 由于这种方法适应性强，形式单纯、规范，因而自 50 年代以来, 在计算机的配合下， 有限元素法在物理和工程设计计算的许多领域得到了广泛的应用。该方法不仅适用于 电磁场问题的求解, 也是对其它具有复杂边值问题的数学物理方程求解时的高效能的 方法
基本思想 有限元素方法是基于变分原理,既通过求解一个泛函取极小值的变分问题。 有限元素法是在变分原理的基础上吸收差分格式的思想发展起来的,是变分问题中欧拉法的 进一步发展。它是人们在尝试求解具有复杂区域,复杂边界条件下的数学物理方程的过程中, 找到的一种比较完美的离散化方法。 它比有限差分法的矩形网格划分方法在布局上更为合理,在处理复杂区域和复杂边界条件时 更方便和适当。采用有限元素法还能使物理特性基本上被保持,计算精度和收敛性进一步得 到保证。正是由于有限元素法有这样一些优点,尽管其计算格式比较复杂,但仍然在很多场合 代替了差分法而受到计算物理工作者的偏爱。 不过需要指出的是:并不是所有有限差分法可以处理的问题都可以采用有限元素法。 泛函: 数学上,通常变量与变量间的关系称为函数,而泛函则是函数集合的函数,也就是函数的函数。 例如,静电场的势函数o(是定义在坐标空间的函数集,系统电场总能量U(o)则是定义在该函 数集中的一个泛函,可记为Ⅳ(o
基本思想： 有限元素方法是基于变分原理，既通过求解一个泛函取极小值的变分问题。 有限元素法是在变分原理的基础上吸收差分格式的思想发展起来的, 是变分问题中欧拉法的 进一步发展。它是人们在尝试求解具有复杂区域, 复杂边界条件下的数学物理方程的过程中, 找到的一种比较完美的离散化方法。 它比有限差分法的矩形网格划分方法在布局上更为合理，在处理复杂区域和复杂边界条件时 更方便和适当。采用有限元素法还能使物理特性基本上被保持, 计算精度和收敛性进一步得 到保证。正是由于有限元素法有这样一些优点, 尽管其计算格式比较复杂, 但仍然在很多场合 代替了差分法而受到计算物理工作者的偏爱。 不过需要指出的是：并不是所有有限差分法可以处理的问题都可以采用有限元素法。 泛函： 数学上，通常变量与变量间的关系称为函数，而泛函则是函数集合的函数，也就是函数的函数。 例如，静电场的势函数ϕ(r)是定义在坐标空间的函数集，系统电场总能量 (ϕ(rU ))则是定义在该函 数集中的一个泛函，可记为 (ϕ(rI ))
例子: 为了进一步说明有限元素方法的基本思想,我们考虑一个确定静电势的问题,该场域的 介质中放置了一个球形金属导体,球形金属导体的半径为,球外距离球中心r处的电位为o(r)。 当这个系统处在电荷平衡的状态下时,金属导体上的电荷分布应当是均匀的,导体表面是等电 位的。我们按照通常的作法,把从导体表面到无穷远处的球面之间的空间,作为导体外的全空间。 假定在这个导体外的空间中的体电荷密度到处为零。 则在此空间中的能量为 r (51.1) 同时该系统的能量应当取最小值,即该系统的能量变分应当满足 观()4=42列- (5.1.2) o dr ar ar ar ar 这里E为介质的相对介电常数,积分是对导体外的空间进行的。因为导体边界上的电位为常数9, 无穷远处的电位为零
例子： 为了进一步说明有限元素方法的基本思想，我们考虑一个确定静电势的问题，该场域的 介质中放置了一个球形金属导体，球形金属导体的半径为r0 , 球外距离球中心 r 处的电位为ϕ(r)。 当这个系统处在电荷平衡的状态下时，金属导体上的电荷分布应当是均匀的, 导体表面是等电 位的。我们按照通常的作法，把从导体表面到无穷远处的球面之间的空间, 作为导体外的全空间。 假定在这个导体外的空间中的体电荷密度到处为零。 则在此空间中的能量为 drr r drr r U dVE r r r 2 2 2 2 2 0 0 0 24 2 2 )( ∫∫ ∫ +∞ +∞ +∞ ⎟⎠⎞ ⎜⎝⎛ ∂∂ ⎟ = ⎠⎞ ⎜⎝⎛ ∂∂ = = ϕ εππ ε ε ϕ ϕ . (5.1.1) 同时该系统的能量应当取最小值，即该系统的能量变分应当满足 ( ) ( ) 4)( 4 0 0 0 0 2 2 2 =⎟⎟⎠⎞ ⎜⎜⎝⎛ ⎭⎬⎫ ⎩⎨⎧ ⎟⎠⎞ ⎜⎝⎛ ∂∂ ∂∂ − ∂∂ = ∂ ∂ ∂∂ = ∫ ∫ ∞+ +∞ ∞+ r r r dr r r r r rdr rr rrU δϕ ϕ δϕ ϕ επ δϕϕ επϕδ . (5.1.2) 这里ε 为介质的相对介电常数，积分是对导体外的空间进行的。因为导体边界上的电位为常数ϕ 0 , 无穷远处的电位为零
则从公式(5.12)可以得到将能量U()取最小值的势函数o必须满足特定的边界条件和如下球坐标下径 向的微分方程 9= 0\=0 ar、or (51.3) 因此,求此微分方程解的问题,可以在数学上等价于找到一个势函数卯,使得积分U(l)取极小值的 问题。实际上采用方程(512)求泛函的极值和解欧拉方程,在数学上都可以代表同一个物理问题。对 两者求得近似解都具有同样的效果。但是在实际计算中,对后者求解往往是困难的,而对前者求近似 解则常常并不太困难。 上面的势函数o门)是定义在坐标空间的函数集,系统电场总能量U(o()则是定义在该函数集中的 个泛函,记为1or)。类似于普通函数取极值的条件,若泛函1()在处取的极值,那末泛函在该 处的变分应当为零。用数学公式表示为 (()=0 上面(5)式中所表示的总电场能量U()则是势函数的一个泛函。对该泛函的变分得到的微分 方程(513)和边值条件=9和-∞-=0的第一类边值问题与场能量U变分的极值问题是等价的
则从公式(5.1.2)可以得到将能量U(ϕ)取最小值的势函数ϕ必须满足特定的边界条件和如下球坐标下径 向的微分方程： 0 1 2 2 2 ⎟ =⎠⎞ ⎜⎝⎛ ∂∂ ∂∂ =∇ r r r r ϕ ϕ . (5.1.3) 因此，求此微分方程解的问题，可以在数学上等价于找到一个势函数ϕ ，使得积分U(ϕ)取极小值的 问题。实际上采用方程(5.1.2)求泛函的极值和解欧拉方程，在数学上都可以代表同一个物理问题。对 两者求得近似解都具有同样的效果。但是在实际计算中，对后者求解往往是困难的，而对前者求近似 解则常常并不太困难。 上面的势函数ϕ(r)是定义在坐标空间的函数集，系统电场总能量 (ϕ(rU ))则是定义在该函数集中的 一个泛函，记为 (ϕ( )rI )。类似于普通函数取极值的条件，若泛函 I(ϕ)在ϕ 0 处取的极值，那末泛函在该 处的变分应当为零。用数学公式表示为 δ (ϕ 0 ( )rI ) = 0 . (5.1.4) 上面(5.1.1)式中所表示的总电场能量U(ϕ)则是势函数ϕ 的一个泛函。对该泛函的变分得到的微分 方程(5.1.3)和边值条件 0 0 = ϕϕ =rr 和 = 0 r +∞= ϕ 的第一类边值问题与场能量 变分的极值问题是等价的。 U
若介质空间存在电荷分布p,则这个静电问题的电场总能量U()为如下积分表示。该积分是未知 势函数的函数,也成为泛函。 (5.1.5) 般来讲,精确估计出此泛函在极值情况下￠函数的形式是不可能的。但是原则上我们可以采用 猜测出的函数近似表示叭=,,其中6对应于N个未知的参数日,再计算泛函(o),然后用取最小值 的条件 O(0) .2=0.(=12,M) (516) 得到N个方程,这个方程组可能用来求出参数θ的解。如同在有限差分法中一样,这个解φ仍然是场微 分方程的近似,但是,该近似方法在参数很少的时候,近似程度还是很好的。 有限元素法是将网络节点上的函数的离散值作为参数,而网络元素内的该势函数值则采用多项 式插值从周围临近节点上的这些参数值求出。例如,我们选择用三角形元素将求解区域划分为子区间 的网络,对泛函()求极小值,就得到节点上未知的势函数的值,然后采用线性插值法,则可以求出 在一个三角形元素内的任意一点(xy)上的势函数值。有限元素法的最后解是势函数在这些节点上的 估计值。由于用来求泛函极小值的函数是近似的线性迭代函数,因而所得到的节点上的势函数值并不 是精确解。该截断误差可以通过减小元素的尺寸或提高迭代函数的阶数来降低
若介质空间存在电荷分布 ρ , 则这个静电问题的电场总能量U(ϕ)为如下积分表示。该积分是未知 势函数ϕ 的函数，也成为泛函。 I( ) dV V ∫ ⎥⎦⎤ ⎢⎣⎡ −∇= ρϕϕ ε ϕ 2 2 . (5.1.5) 一般来讲, 精确估计出此泛函在极值情况下ϕ 函数的形式是不可能的。但是原则上我们可以采用 猜测出的函数近似表示 ( θφ ) r zyx ,,, , 其中θr对应于N 个未知的参数θ i , 再计算泛函I(ϕ)，然后用取最小值 的条件 ( ) ( ) Ni I i == ,...,2,1,0 ∂∂ θφ . (5.1.6) 得到 个方程 N , 这个方程组可能用来求出参数θ i的解。如同在有限差分法中一样, 这个解φ 仍然是场微 分方程的近似, 但是, 该近似方法在参数很少的时候, 近似程度还是很好的。 有限元素法是将网络节点上的函数ϕ 的离散值作为参数, 而网络元素内的该势函数值则采用多项 式插值从周围临近节点上的这些参数值求出。例如, 我们选择用三角形元素将求解区域划分为子区间 的网络，对泛函 I( ) ϕ 求极小值, 就得到节点上未知的势函数的值，然后采用线性插值法，则可以求出 在一个三角形元素内的任意一点( )上的势函数值。有限元素法的最后解是势函数在这些节点上的 估计值。由于用来求泛函极小值的函数是近似的线性迭代函数, 因而所得到的节点上的势函数值并不 是精确解。该截断误差可以通过减小元素的尺寸或提高迭代函数的阶数来降低。 , yx
许多物理问题的分析结果在数学上都可以归结为下面形式的重要微分方程 pVo)+go=p (5.17) 它在边界r上至少有部分的边界条件是个狄利克利问题,即 (51.8) 而其余的边界则满足纽曼或者混合边界条件,它们可以写为: on +g(s)=b(s) (5.1.9) 对应于上面的微分方程517和边界条件(518,(5.19)的泛函应当是 (o)=(plvof+gg2-2pp)dW+ (2-2bp)ds (5.1.10 其中r()为以r为边界的体积(对三维问题)或面积区域(对二维问题);S为边界上的一部分边 界,在S上势函数满足混合边界条件(51.9)。 在二维情况下,如果p=E,q=a,b=aB,8=0,S为整个r边界的情况下,微分方程(5.1.7)及 边界条件(5.9)可以写为 平面场域为D,L为D的边界,s为边界上的点 (51.11) +a(r,y)o=B(s)
许多物理问题的分析结果在数学上都可以归结为下面形式的重要微分方程 −∇ ∇ + = ( p g ϕ ϕρ ) . (5.1.7) 它在边界 Γ上至少有部分的边界条件是个狄利克利问题，即 ϕ = F s( ) . (5.1.8) 而其余的边界则满足纽曼或者混合边界条件，它们可以写为： sbsq )()( n =+ ∂ ∂ ϕ ϕ r . (5.1.9) 对应于上面的微分方程(5.1.7)和边界条件(5.1.8), (5.1.9)的泛函应当是 ( ) ( ) 2 2 2 () ( 2 ) ( 2 ) V S I ϕ ϕ ϕ ρϕ ϕ p g dV q b dS ϕ . (5.1.10) Γ Γ′ = ∇+ − + − ∫ ∫ 其中 为以 V( ) Γ Γ 为边界的体积（对三维问题）或面积区域（对二维问题）； S′为 Γ 边界上的一部分边 界，在 S′上势函数满足混合边界条件(5.1.9)。 在二维情况下，如果 p = ε ，q = εα , b = εβ , g = 0 , S’为整个 Γ边界的情况下，微分方程（5.1.7）及 边界条件（5.1.9）可以写为 ⎪ ⎪ ⎩ ⎪ ⎪ ⎨ ⎧ ⎟ = ⎠ ⎞ ⎜ ⎝ ⎛ + ∂ ∂ −= ∂ ∂ + ∂ ∂ )(),( 2 2 2 2 yx s n xx L βϕα ϕ ε ρϕϕ r 平面场域为 D，L 为 D 的边界， s 为边界上的点。 (5.1.11)
根据公式(5110),此时的泛函可以取为 1()(2o+g.m(x02 (51.12) 证明:求泛函(51.12)的极值与满足上述边界条件下的微分方程511)的求解是等价的。 做泛函(5)的变分 0(7)+2(-1 (51.13) D(L) 利用格林公式 uV v+ Vv. Vuydxdy=p, uVvnd (n为D区域的边界L上的外法线方向的单位矢量) 公式(5113)变为 8/=-2 evp+pyopdrd+2cnf dg ESpdl-2sh. l sodl (51.15) 由(510)式的泛函取极值的条件和>的任意性就得到了公式(51)中的偏微分方程
根据公式(5.1.10)，此时的泛函可以取为 ( ) ( ) ( ) 2 2 ( ) ( ) 2 2 ϕ dl . (5.1.12) D L S I d ϕ ε ϕ ρϕ ε αϕ β V ′ Γ = ∇− + − ∫ ∫ 证明：求泛函(5.1.12)的极值与满足上述边界条件下的微分方程(5.1.11)的求解是等价的。 做泛函(5.1.12)的变分 { ( ) } ( ) ( ) 2 2 D L L δ I d = ∇ ⋅∇ − + − ε ϕ δϕ ρδϕ ε αϕδϕ βδϕ xdy dl ∫ ∫ { } ( ) 2 2 D L L dxdy dl nϕ ε ϕ δϕ ρδϕ ε δϕ ∂ = ∇ ⋅∇ − − ∂ ∫ ∫ r . (5.1.13) 利用格林公式 { } u v v u dxdy u v ndl D L ∇ +∇ ⋅∇ = ∇ ⋅ ∫ ∫ ∫ 2 $ . (5.1.14) ( 为 D 区域的边界 L 上的外法线方向的单位矢量) n$ 公式(5.1.13)变为 2 ( ) 2{ } 2 2 D L L L I dxdy dl dl n n ϕ ϕ δ ε ϕ ρ δϕ ε δϕ ε δϕ ∂ ∂ =− ∇ + + − ∂ ∫ ∫∫ r ∂r . (5.1.15) 由(5.1.10)式的泛函取极值的条件和δϕ 的任意性就得到了公式(5.1.11)中的偏微分方程
当偏微分方程满足第一类边界条件时,即=0,由于0的任意性,公式(53)中的第二项的 变分为零,所以和第一类边值问题等价的变分泛函为: (5116 对第二类边值问题,即a=0时,等价的泛函为 5.1.17) 特别是当边界为导体面时,由于导体面是等电位的,则在边界上电位为常数g 此时(5117)式可以化为 )J0-c+-9 公式(5118)中的q为导体表面上的电荷量。 由变分原理可以知道对上述平面泊松方程的第一、二、三类边值问题都可以等价地化为求泛函极 值(或称为变分问题)来处理。我们从上面的分析可以看到,对泛函()求极值会自动保证满足边界 条件。同有限差分法中的边界问题比较,特别是节点不在边界上时会带来很大麻烦相比较,这是有限 元素法的最大的优点。若在此基础上再进行离散化,这就导致了有限元素方法。这种离散方法是通过 网格离散化的处理,用构造分片光滑的基函数{来以变分法求得近似解的
当偏微分方程满足第一类边界条件时，即 = 0 ∂ ∂ L n r ϕ ，由于δϕ 的任意性，公式(5.1.13)中的第二项的 变分为零，所以和第一类边值问题等价的变分泛函为： ( ) ( ) 2 ( ) 2 D L I ϕ ε ϕ ρϕ = ∇− ∫ dx dy . (5.1.16) 对第二类边值问题，即 α = 0时，等价的泛函为 ( ) ( ) 2 ( ) 2 2 D L L I d = ∇− − ε ϕ ρϕ ε βϕ x dy dl ϕ ∫ ∫ . (5.1.17) 特别是当边界为导体面时，由于导体面是等电位的，则在边界上电位 ϕ 为常数 ϕ 0 。 此时(5.1.17)式可以化为 ( ) 2 0 ( ) 1 D L 2 I ϕ ε ϕ ρϕ ϕ dxdy q ⎛ ⎞ = ∇− ⎜ ⎟ − . (5.1.18) ⎝ ⎠ ∫ 公式(5.1.18)中的 q 为导体表面上的电荷量。 由变分原理可以知道对上述平面泊松方程的第一、二、三类边值问题都可以等价地化为求泛函极 值（或称为变分问题）来处理。我们从上面的分析可以看到，对泛函 I(ϕ)求极值会自动保证满足边界 条件。同有限差分法中的边界问题比较，特别是节点不在边界上时会带来很大麻烦相比较，这是有限 元素法的最大的优点。若在此基础上再进行离散化，这就导致了有限元素方法。这种离散方法是通过 网格离散化的处理，用构造分片光滑的基函数{ϕ k }来以变分法求得近似解的
5.2二维场的有限元素法 在这节中以满足第一类边界条件的二维平面场泊松方程为例具体地来讨论如何构造有限元素法的 计算格式。 假定该问题的求解场域为D的区域。对该问题的变分法处理可以归结为求解满足边界条件叫=F 的ox,),使得对于任意的,公式516)所示的泛函变分为零。即 (o)=0 在找出与边值问题相对应的泛函及其变分问题以后,就需要对待求解区域进行划分,将其离散为 有限个元素的集合,然后进行分片插值建立计算格式。 1.场域划分的约定 采用有限元素法对平面场域D的分割时,常用的办法是用一些分割直线将D划分为许多三角形单 元(如图521所示) 因为三角形子区间的计算格式最为简单的。采用三角形元素划分场域时,我们允许场域内各三角 形元素的大小及形状可以不一样。三角形元素越小,场域的分割就越细,计算的精度就会越高。因而 在实际应用中是按精度的要求来决定场域内各处三角形元素的大小
5．2 二维场的有限元素法 在这节中以满足第一类边界条件的二维平面场泊松方程为例具体地来讨论如何构造有限元素法的 计算格式。 假定该问题的求解场域为 D 的区域。对该问题的变分法处理可以归结为求解满足边界条件ϕ ( ) L = F s 的ϕ(x y, ) , 使得对于任意的δϕ ，公式(5.1.16)所示的泛函变分为零。即 δI( ) ϕ = 0 . 在找出与边值问题相对应的泛函及其变分问题以后，就需要对待求解区域进行划分，将其离散为 有限个元素的集合，然后进行分片插值建立计算格式。 1． 场域划分的约定 采用有限元素法对平面场域 D 的分割时，常用的办法是用一些分割直线将 D 划分为许多三角形单 元(如图 5.2.1 所示)。 因为三角形子区间的计算格式最为简单的。采用三角形元素划分场域时，我们允许场域内各三角 形元素的大小及形状可以不一样。三角形元素越小，场域的分割就越细，计算的精度就会越高。因而 在实际应用中是按精度的要求来决定场域内各处三角形元素的大小
划分要求: 般规定每个三角形元素的三个边的边长尽量地接近,尽量避免三角形元素具有大的钝 角,一般最长的一条边不得大于最短边的三倍。 (2)在分割场域时要求各三角形元素之间只能以顶点相交,即两相邻的三角形元素有两个公 共的顶点及一条等长的公共边。但是不能把一个三角形的顶点取在另一个三角形的边上。 (3)在边界上,我们可以将三角形元素的两个顶点放在边界曲线上,近似地用这两个顶点间 的三角形边来代替边界上这段曲线。 (4)划分时还应当注意要尽量地使由相邻边界节点之间的线段所近似构成的曲线足够光滑 如果在场域D内有不同的介质,则需要将介质的交面线选为分割线。 按照上述三角形单元分割原则,我们可以看出这样的分割是适应于各种复杂几何形状的场域的。 图521允许的三角形元素的划分 图522不允许的三角形元素的划分
划分要求： （1） 一般规定每个三角形元素的三个边的边长尽量地接近，尽量避免三角形元素具有大的钝 角，一般最长的一条边不得大于最短边的三倍。 （2） 在分割场域时要求各三角形元素之间只能以顶点相交，即两相邻的三角形元素有两个公 共的顶点及一条等长的公共边。但是不能把一个三角形的顶点取在另一个三角形的边上。 （3） 在边界上，我们可以将三角形元素的两个顶点放在边界曲线上，近似地用这两个顶点间 的三角形边来代替边界上这段曲线。 （4） 划分时还应当注意要尽量地使由相邻边界节点之间的线段所近似构成的曲线足够光滑。 如果在场域 D 内有不同的介质，则需要将介质的交面线选为分割线。 按照上述三角形单元分割原则，我们可以看出这样的分割是适应于各种复杂几何形状的场域的。 0 x y m j i 图 5.2.1 允许的三角形元素的划分 图 5.2.2 不允许的三角形元素的划分