第五章 有限元素方法 5.1 有限元素方法的基本思想 特点： （1 ） 能处理复杂区域和复杂边界条件的求解问题。 （2 ） 有限元素法是一套求解微分方程的系统化数值计算方法。 它比传统解法具有理论完整可靠，物理意义直观明确，解题效能强等优点。 （3 ） 对连续体的问题采用有限元素法 , 是将连续问题离散化的数值求解方法。 应用范围： 由于这种方法适应性强，形式单纯、规范，因而自 50 年代以来, 在计算机的配合下， 有限元素法在物理和工程设计计算的许多领域得 到了广泛的应用。该方法不仅适用于 电磁场问题的求解, 也是对其它具有复杂边值问题的数学物理方程求解时的高效能的 方法。 基本思想： 有限元素方法是基于变分原理 ，既通过求解一个泛函取极小值的变分问题。 有限元素法是在变分原理的基础上吸收差分格式的思想发展起来的 , 是变分问题中欧拉法的 进一步发展。它是人们在尝试求解具有复杂区域 , 复杂边界条件下的数学物理方程的过程中 , 找到的一种比较完美的离散化方法。 它比有限差分法的矩形网格划分方法在布局上更 为合理，在处理复杂区域和复杂边界条件时 更方便和适当。采用有限元素法还 能使物理特性基本上被保持, 计算精度和收敛性进一步得 到保证。正是由于有限元素法有这样一些优点 , 尽管其计算格式比较复杂 , 但仍然在很多场合 代替了差分法而受到计算物理工作者的偏爱。 不过需要指出的是： 并不是所有有限差分法可以处理的 问题都可以采用有限元素法。 泛函： 数学上，通常变量与变量间的关系称为函数，而泛函则是函数集合的函数，也就是函数的函数。 ( )( )rU ? 则是定义在该函 例如，静电场的势函数 ?(r)是定义在坐标空间的函数集，系统电场总能量 ( )( )rI ? 。 数集中的一个泛函，可记为 例子： 为了进一步说明有限元素方法的基本思想，我们考虑一个确定静电势的问题，该场域的 介质中放置了一个球形金属导体，球形金属导体的半径为 , 球外距离球中心 r 处的电位为 ?(r)。 0 r 当这个系统处在电荷平衡的状态下时，金属导体上的电荷分布应当是均匀的, 导体表面是等电 位的。我们按照通常的作法，把从导体表面到无穷远处的球面之间的空间 , 作为导体外的全空间。 假定在这个导体外的空间中的体电荷密度到处为零。 则在此空间中的能量为 drr r drr r dVEU rrr 2 2 2 2 2 000 24 22 )( ∫∫∫ +∞+∞+∞ ? ? ? ? ? ? ? ? =? ? ? ? ? ? ? ? == ? εππ ?εε ? . (5.1.1) 同时该系统的能量应当取最小值， 即该系统的能量变分应当满足 () () 044)( 000 222 = ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? = ? ? ? ? = ∫∫ ∞+ +∞ ∞+ rrr dr r r rr rdr rr rrU δ? ? δ? ? επ δ?? επ?δ . (5.1.2) ε 为介质的相对介电常数，积分是对导体外的空间进 行的。因为导体边界上的电位为常数 0 ? , 这里 无穷远处的电位为零。 则从公式 (5.1.2)可以得到将能量 ( )?U 取最小值的势函数? 必须满足特定的边界条件和如下球坐标下径 向的微分方程： 0 1 2 2 2 = ? ? ? ? ? ? ? ? ? ? =? r r rr ? ? . (5.1.3) 因此， 求此微分方程解的问题，可以在数 学上等价于找到一个势函数 ? ，使得积分 ( )?U 取极小值的 问题。 实际上采用方程 (5.1.2)求泛函的极值和解欧拉方程，在数学上都可以代表同一个物理问题。对 两者求得近似解都具有同样的效果。但是在实际计算中 ，对后者求解往往是困难的，而对前者求近似 解则常常并不太困难。 上面的势函数?(r) 是定义在坐标空间的函数集，系统电场总能量 ( )( )rU ? 则是定义在该函数集中的 一个泛函，记为 ()()rI ? 。类似于普通函数取极值的条件，若泛函 ( )?I 在 ? 0 处取的极值，那末泛函在该 处的变分应当为零。用数学公式表示为 ()( ) 0 0 =rI ?δ . (5.1.4) 上面(5.1.1) 式中所表示的总电场能量 ( )?U 则是势函数 ? 的一个泛函。 对该泛函的变分得到的微分 方程(5.1.3) 和边值条件 0 0 ?? = =rr 和 0= +∞=r ? 的第一类边值问题与场能量 变分的极值问题是等价的。 U 若介质空间存在电荷分布 ρ , 则这个静电问题的电场总能量 ( )?U 为如下积分表示。该积分是未知 势函数 ? 的函数，也成为泛函。 () dVI 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, 许多物理问题的分析结果在数学上都可以归结为下面形式的重要微分方程 ( ) ?? ? + =pg??ρ . (5.1.7) 它在边界 Γ 上至少有部分的边界条件是个狄利克利问题，即 ( ) ? = Fs . (5.1.8) 而其余的边界则满足纽曼或者混合 边界条件，它们可以写为： )()( sbsq n =+ ? ? ? ? r . (5.1.9) 对应于上面的微分方程 (5.1.7)和边界条件(5.1.8), (5.1.9) 的泛函应当是 ()() 2 22 () ( 2 ) ( 2 ) VS I p g dV q b dS???ρ??? . (5.1.10) ′ΓΓ =?+? + ? ∫∫ 其中 为以 () V Γ Γ 为边界的体积（对三维问题）或 面积区域（对二维问题）； S′为 Γ 边界上的一部分边 界，在 S ′ 上势函数满足混合边界条件(5.1.9) 。 在二维情况下 ，如果 p ε= ， q εα= , b εβ= , 0g = , S’为整个 Γ边界的情况下，微分方程（ 5.1.7）及 边界条件（5.1.9 ）可以写为 ? ? ? ? ? ? ? =? ? ? ? ? ? + ? ? ?= ? ? + ? ? )(),( 2 2 2 2 syx n xx L β?α ? ε ρ?? r 平面场域为 D，L 为 D 的边界， s 为边界上的点。 (5.1.11) 根据公式 (5.1.10)，此时的泛函可以取为 () ( ) ( ) 2 2 () () dl? . (5.1.12) DL S IdV?ε?ρ?εα?β ′ Γ =?? + ? ∫∫ 证明： 求泛函(5.1.12)的极值与满足上述边界条件下的微分方程 (5.1.11)的求解是等价的。 做泛函(5.1.12)的变分 ( ){ } ( ) () 22 DL L Idxdydlδ ε ? δ? ρδ? ε α?δ? βδ?=???? + ? ∫∫null {} () 22 DL L dxdy dl n ? ε ?δ?ρδ? ε δ? ? =???? ? ? ∫ r null∫ . (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 DL L L Idxdydldl nn ? ? δ ε?ρδ? ε δ? ε δ? ? ? =? ? + + ? ∫∫∫ r nullnull r. (5.1.15) 由(5.1.10) 式的泛函取极值的条件和 δ? 的任意性就得到了公式 (5.1.11)中的偏微分方程。 当偏微分方程满足第一类边界条件时，即 0= ? ? L n r ? ，由于 δ? 的任意性，公式 (5.1.13)中的第二项的 变分为零，所以和第一类边值问题等价的变分泛函为： () ( ) 2 () 2 DL Idxy . (5.1.16) ?ε?ρ?=?? ∫ 对第二类边值问题，即 α = 0时，等价的泛函为 () ( ) 2 () 22 DL L Idxdydlε? ρ? εβ?=?? ? ∫null∫ . (5.1.17) ? 特别是当边界为导体面时，由于导体面是等电位的，则在边界上电位 ? 为常数 0 ? 。 此时(5.1.17) 式可以化为 () 2 0 () 1 2 DL I dxdy q? ε? ρ? ? ?? =??? . (5.1.18) ?? ?? ∫ 公式(5.1.18) 中的 q 为导体表面上的电荷量。 由变分原理可以知道对上述平面泊松方程的第一、二、三类边值问题都可以等价地化为求泛函极 值（或称为变分问题）来处理。我们从上面的分析可以看到，对泛函 ( )?I 求极值会自动保证满足边界 条件。同有限差分法中的边界问题比较，特别是节点不 在边界上时会带来很大麻烦相比较，这是有限 元素法的最大的优点。若在此基础上再进行离散化，这就导致了有限元素方法。这种离散方法是通过 网格离散化的处理，用构造分片光滑的基函数 { }? k 来以变分法求得近似解的。 5． 2 二维场的有限元素法 在这节中以满足第一类边界条件的二维平面场泊松方程为例具体地来讨论如何构造有限元素法的 计算格式。 假定该问题的求解场域为 D 的区域。对该问题的变分法处 理可以归结为求解满足边界条件 ()? L Fs= 的 ), 使得对于任意的( ? xy, δ? ，公式(5.1.16) 所示的泛函变分为零。即 () 0=?δI . 在找出与边值问题相对应的泛函及其变分问题以后， 就需要对待求解区域进行划分，将其离散为 有限个元素的集合，然后进行分片插值建立计算格式。 1． 场域划分的约定 采用有限元素法对平面场域 D 的分割时， 常用的办法是用一些分割直线将 D 划分为许多 三角形单 元(如图 5.2.1 所示 )。 因为三角形子区间的计算格式最为简单的。采用三角形元素划分场域时，我们允许场域内各三角 形元素的大小及形状可以不一样。三角形元素越小，场 域的分割就越细，计算的精度就会越高。因而 在实际应用中是按精度的要求来决定场域内各处三角形元素的大小。 划分要求： （1 ） 一般规定每个三角形元素的三个边的边长 尽量地接近，尽量避免三角形元素具有大的钝 角，一般最长的一条边不得大于最短边的三倍。 （2 ） 在分割场域时要求各三角形元素之间只能 以顶点相交，即两相邻的三角形元素有两个公 共的顶点及一条等长的公共边。 但是不能把一个三角形的顶点取在另一个三角形的边上。 （3 ） 在边界上，我们可以将三角形元素的两个 顶点放在边界曲线上，近似地用这两个顶点间 的三角形边来代替边界上这段曲线。 （4 ） 划分时还应当注意要尽量地使由相邻边界 节点之间的线段所近似构成的曲线足够光滑。 如果在场域 D 内有不同的介质，则需要将介质的交面线选为分割线。 按照上述三角形单元分割原则，我们可以看出这样的分割是适应于各种复杂几何形状的场域的。 0 x y m j i 图 5.2.1 允许的三角形元素的划分 图 5.2.2 不允许的三角形元素的划分 对元素和节点分别进行统一的编号： 节点： 三角形单元的顶点。 节点参数：函数 ),( yx? 在节点处 ),,( mjil = 的值 ),( ll yx? 。 （1 ） 为了计算的简便和格式的统一，各个元素 节点的编号顺序要一致，一般规定各个元素的 节点编号顺序选取逆时针的方向，如图 5.2.1 中某元素（ e）的三个顶点的编号顺序取为 mji ,, 。 （2 ） 在对所有节点进行总体编号时，还应当考 虑尽量使得在此编号下，各元素中最大编号值 和最小编号值之差的最大值要尽可能小。这样可以使得最后形成的代数方程组的系数矩 阵的带宽较小，从而减少计算机的存贮量。 （3 ） 为了统一计算格式，还规定每个边界元素 只应当有一条边落在边界曲线上，这一边相对 的顶点取编号为 。 i 完成场域的划分之后，等效的泛函就可以是各个三角形单元泛函的和 . (5.2.1) () () () ∑ = = 0 1 e e e e II ?? 2． 计算格式的建立 对于任一个元素（e ），我们采用符号 )来标记三角形元素三个顶点上的函 数值（节点参数）。如果（e ）元素很小，函数 ),( )()( ll ee l yx?? ≡ (mjil ,,= ),( yx? 在（ e）内近似认为是随 yx, 线性变化的。这相当 于在这个局域范围内，场可以看成是近似均匀的。这样我们可以用 线性插值法 来构造在元素（e ）内 部任一点上的势函数值 ()yx,? (以下我们略去元素(e) 的上标) ，即 ygxggyx 321 ),( ++==?? . (5.2.2) 其中 和 由元素（ e）上的三个节点的函数值来决定。由上式可以得到方程组 21 ,gg 3 g ( ) () () ? ? ? ? ? ==++ ==++ ==++ mmmmm jjjjj iiiii yxygxgg yxygxgg yxygxgg ?? ?? ?? , , , 321 321 321 . 由此很容易解出： ( ) () () ? ? ? ? ? Δ++= Δ++= Δ++= 2 2 2 3 2 1 mmjjii mmjjii mmjjii cccg bbbg aaag ??? ??? ??? . (5.2.3) 其中 为(e) 元素的三角形面积。 Δ () Δ= = ? 1 2 1 1 1 1 2 xy xy xy bc bc ii jj mm ij ji ， . (5.2.4) axxxx byy cxx ijmmj ijm imj =? =? =? ? ? ? ? ? 其余的 j 及 m 则可以由公式 (5.2.5)按下标 abc jj ,, ,,abc mm mji ,, 的顺序轮换得到。 引入三角形型函数(the shape functions for the triangle) ，其定义为 ( ) ( ) . (5.2.5) ),,( mjil =Δ++≡ 2, ycxbayxN llll , 利用上式，并将 (5.2.3)式代入(5.2.1) 式中就得到了 ( )e 元素内任意一点 ( )yx, 的势函数的插值 . (5.2.6) () ∑ = = mj il ll Nyx , , ?? 即在三角形元素 内任意一点的函数值 ()e ( )yx,? 是由该元素的三个节点参数 ( ) ),,(, mjilyx ll =? 唯一确定下 来的，函数 ( )yx,? 在每个元素中是局域线性的。 并且通过这种线性插值法构造的插值函数显然可以保 证在两相邻元素的公共边上函数值连续。 求解该场微分方程问题时所对应的泛函的有限元素法表达式推导： 首先，我们将此第一类边值的两维平面场的变分问题的提法重新列于下面公式中： ( ) () ? ? ? ? ? ? ? = ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? + ? ? ? ? ? ? ? ? = = ∫∫ 0 2 2 | 2 0 ?? ρ? ??ε ? ?δ L D dxdy yx I I . (5.2.7) 很明显，如同公式(5.2.1) ，这时 总的泛函应为各三角形元素上泛函的代数和，即 () ∑ ∫∫ ∫∫ = ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? + ? ? ? ? ? ? ? ? = 0 1 2 2 2 e e ee dxdydxdy yx I ρ? ??ε ? . (5.2.8) () ()[]() ()???? 21 1 21 0 IIII e e ee ?=?= ∑ = 用向量记法表示公式(5.2.8) 中的 ( )? 1 I ：我们设三角形单元 的节点编号为 ()e mji ,, ，并定义列向量 ， . () ? ? ? ? ? ? ? ? ? ? =Φ m j i e ? ? ? () ? ? ? ? ? ? ? ? ? ? = m j i e N N N N 则公式（ 5.2.6）可以写为 ( )( )( ) ( ) e T ee T e NNyx Φ=Φ= )(,? . (5.2.9) 由公式(5.2.6) 可以求得在元素(e)内对函数 ?的偏微商为 ∑ = Δ = ? ? mj il ll b x , 2 1 ? ? ， ∑ = Δ = ? ? mj il ll c y , 2 1 ? ? (5.2.10) 若记列向量 () , ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? =? y x e ? ? ? 并定义 () ? ? ? ? ? ? ? ? Δ = mji mji e ccc bbb B 2 1 . 则公式(5.2.10)可以改写为 ( )()() eee B Φ=?? . (5.2.11) 于是有 ()() ()()[]()()[] ∫∫∫ ΦΦ=??= e ee T ee e e T ee dxdyBBdxdyI 22 1 ε ?? ε () ()() () ()()() ee T ee e e T e T e KdxdyBB ΦΦ=Φ ? ? ? ? ? ? Φ= ∫∫ 2 1 2 1 ε . (5.2.12) 其中 不是坐标的函数，因而我们可以将它移出积分号外。 () e Φ 在(5.2.12) 式中我们定义了 () ()() . (5.2.13) ∫∫ = e e T ee dxdyBBK ε 同样， 也不是坐标的函数，也可以移出积分号外。 () e B 导出 的表达式 () e K () ? ? ? ? ? ? ? ? ? ? = e mm e mj e mi e jm e jj e ji e im e ij e ii e KKK KKK KKK K ? ? ? ? ? ? ? ? ? ? +++ +++ +++ Δ = 22 22 22 4 mmjmjmimim mjmjjjijij mimijijiii cbccbbccbb ccbbcbccbb ccbbccbbcb ε . (5.2.14) 由此可见 ( ) e K 是一个三阶正定对称方阵，它的一般形式可以写为 () srsr e sr e rs ccbbKK + Δ == 4 ε ， (r, s= i, j, m) . (5.2.15) 最后我们得到 ()? 1 I 的向量表示为 () () ()()() ∑∑ == ΦΦ== 00 11 11 2 1 e e ee T e e e e KII ?? . (5.2.16) 考虑 ( )? 2 W 的向量记法： 假定三角形元素足够小， ρ 值可以取等于 ji ρρ , 和 m ρ 的平均值 e ρ 。 将公式(5.2.9) 代入（ 5.2.8)式中 的表示，则可以得到 () e I 2 ( ) ( ) ( ) ∫∫∫∫ ∫∫ Φ=Φ== e ee T e ee e T eee dxdyNdxdyNdxdyI ρρρ? 2 . (5.2.17) 定义 ( ) ( ) ∫∫ ≡ e eee dxdyNP ρ . (5.2.18) 于是(5.2.17) 式可以写为 ()() () ? ? ? ? ? ? ? ? ? ? Δ Φ=Φ= 1 1 1 3 2 e T ee T ee PI ρ . (5.2.19) 在上式的推导中，我们用到如下三角形型函数的积分公式(参见附录 D) 3 )( Δ == ∫∫ e k e k dxdyNt ， ),,( mjik = . (5.2.20) 从(5.2.19) 式中可以看出 () ? ? ? ? ? ? ? ? ? ? Δ = ? ? ? ? ? ? ? ? ? ? = 1 1 1 3 )( )( )( e e m e j e i e p p p P ρ . (5.2.21) 最后我们得到 ()? 2 I 的向量记法为： () () ()() ∑∑ == Φ== 00 11 22 e e e T e e e e PII ?? . (5.2.22) 综合(5.2.16) 和(5.2.22) 式，我们就得到泛函表达式 ( )()( )??? 21 III ?= ()()() ()() ∑∑ == Φ?ΦΦ= 00 11 2 1 e e e T e e e ee T e PK . (5.2.23) 为将元素 (e)上的表示用总体向量来表示 , 引入一个 n×3 阶的辅助矩阵 ( ) e R ， n 为总的节点数。 , () ? ? ? ? ? ? ? ? ? ? = 0......1...0...0...00 0......0...1...0...00 0......0...0...1...00 e R | | | i j (5.2.24) m 则有 ( )()()Φ=Φ ee R . (5.2.25) 其中 是场域内所有节点上的函数值向量 , 它表示为 ()Φ ( )( ) T n ??? ...,, 21 =Φ . (5.2.26) 重新表示公式 (5.2.23)的泛函为 () () ()()()()() ()() ? ? ? ? ? ? Φ?Φ ? ? ? ? ? ? Φ= ∑∑ == 00 11 2 1 e e e T e T e e ee T e T PRRKRI ? ()()()()()PK TT Φ?ΦΦ= 2 1 . (5.2.27) 其中 是场域内所有节点上与 ()P ρ 的值相关的向量 (参见式 (5.2.21)), 它表示为 ( )( ) T n PPPP ,...,, 21 = (5.2.28) 当对公式 (5.2.27)所示泛函取极值 , 就需要满足条件 ()()0=? ? I d d i , ),...,2,1( ni = . (5.2.29) 由微分方程 (5.2.29)可以得到必须满足的线性代数方程组 ( )( ) ( )PK =Φ . (5.2.30) 0=ρ 时对应的方程为拉普拉斯方程。公式 (5.2.30)中的向量 ( )P 为 ( )0 零向量， 即拉普拉斯方程对 应的有限元素方程为齐次线性代数方程组。 ( )( ) ( )0=ΦK . (5.2.31) 由上面的讨论可以看出： ( )K 的矩阵元素是由所相关的三角形元素 对该矩阵元的贡献之和。 具体地讲，其对角线上的某矩阵元 是以 为节点的各三角形元素对该矩阵元的贡献和，即 ∑ ll K l . (5.2.32) = = 为节点的三角形元素以 le e llll KK 矩阵元素 是以 边为邻边的某两个三角形元素的贡献 之和。因此，如果和 节点相邻的节点有 ，那末 )的第 行中除了对角矩阵元 和与第 列相交处的矩阵元非零外，其它的均为 零。所以 lm K lm e lm K l i mmm ,...,, 21 (K l ll K i mmm ,...,, 21 ( )K 是大型稀疏矩阵。又由于 ( ) e K 是正定对称的，因此 ( )K 也应当是正定对称的。在电磁场问题 中，其对应的泛函描述了电磁场的能量，因而在离散化处理后矩阵 ( )K 仍然是正定的。同样我们可以 知道： 的各分量也是各相关的三角形元素贡献之和。 ()P 3． 边界条件处理 由于我们处理的问题本身还要满足第一类边界条件 0 ?? = L ，因此必须把这一要求强制性地综合到 有限元素方程中去。通过在对节点编号时，使 个总节点中的前 个为内部节点，从n 0 n 1 0 +n 到 为边界节 点。即 n ,0 0 ?? = +in ),...,2,1( 0 nni ?= . (5.2.33) 公式 (5.2.33)可以改写为向量形式。为此定义 ( ) () T nnn ??? ,...,, 212 00 ++ ≡Φ ， ( ) T nn )(0,02010 0 ...,, ? ≡Φ ??? . (5.2.34) 上面的 (5.2.33)式就可写为 ( )( ) 02 Φ=Φ . (5.2.35) 我们进一步再定义 ( ) () T n 0 ,...,, 211 ???≡Φ . (5.2.36) 把 都写成相应的分块形式，则线性代数方程组 (5.2.30)变为 ()()PK , ( )( ) ()() ( ) () ( ) () ? ? ? ? ? ? ? ? = ? ? ? ? ? ? ? ? Φ Φ ? ? ? ? ? ? ? ? 2 1 2 1 2221 1211 P P KK KK . (5.2.37) 它的第一个方程为： ( )( ) ( ) ( )( ) 2121111 Φ?=Φ KPK . (5.2.38) 根据边界条件，我们可以强制性地命令上式中 ( ) ( ) 02 Φ=Φ ，这样就得到了强加边界条件处理后的有限 元方程： ( )( ) ( ) ( )( ) ()() ? ? ? Φ=Φ Φ?=Φ 02 2121111 KPK ， (5.2.39) 其中 是 () 11 K 00 nn × 阶的对称方阵， ( ) )1( P 是 维列向量。 0 n 显式地写出公式 (5.2.39)的第一个方程为 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 00000 0 0 2 1 21 22221 11211 ... ....................... ....................... ... ... nnnnn n n KKK KKK KKK ? ? ? M M = . (5.2.40) ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ???? ???? ???? ?++ ?++ ?++ )(002)2(01)1()( )(0202)2(201)1(2)2( )(0102)2(101)1(1)1( 0000000 000 000 ... ............... ............... ... ... nnnnnnnnn nnnnn nnnnn KKKP KKKP KKKP ??? ??? ??? 公式 (5.2.40)还可以简单地记为 ( )( )() 111 PK ′ =Φ . (5.2.41) 4.有限元方程的求解 对有限元方程 (5.2.41)求解。在采用划分场域三角形元素的约定后，矩阵 ( )K 应当是正定对称的大型 稀疏矩阵，可以采用直接法求出有限元的线性方程组的解 , 但是通常仍使用在有限差分法中讲述过的 迭代法来求解。 对高斯 -赛德尔迭代法有如下公式 ii i j n ij i m jij m jij m i KpKK ? ? ? ? ? ? ? ? ?+?= ∑∑ ? =+= ++ 1 11 )()1()1( 0 ??? ， ),...,2,1( 0 ni = .. (5.2.42) 采用超松驰迭代法时，有公式 )()()1( m i m i m i Rω?? += + ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? +??+= ∑∑ ? =+= + )( 1 11 )()1()( / 0 m iii i j n ij i m jij m jij m i KpKK ???ω? . (5.2.43) () ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? +??+?= ∑∑ ? =+= + ii i j n ij i m jij m jij m i KpKK /1 1 11 )()1()( 0 ??ω?ω 如同在有限差分法中求解差分方程组时的情况，求解方 程组(5.2.41)采用超松弛迭代法更为有效。 有限元素法处理复杂边界条件时具有很好的灵活性，并且在划分三角形元素时人们还可以增加在 函数变化剧烈的区域内节点的密度，以得到较高精度的数值结果， 因而这种方法的优点是十分显著的。 5.有限元素法的一般步骤 首先，推导出与给定边界条件的偏 微分方程等价的泛函表示； 第二，把求解的区域用三角形元素划分为小的单元。然后对每个节点和三角形元素按照约定的规 则分别进行编号。 第三，利用公式 (5.2.14-15)和 (5.2.18-21)，计算出各个三角形元素 的系数矩阵 e ( ) e K 和 ( ) e P 。 第四，将各个三角形单元的系数矩阵 ( ) e K 和 ( ) e P 装配成总矩阵 ( )K 和 ( )P ，形成有限元方程组，然 后利用强加边界条件法对有限元方程组进行修正。 最后，利用超松弛迭代法求解有限元方程组，则得到域内各个节点上的函数 ? 值。 5.3 有限元素法与有限差分法的比较 数学方法的不同： 有限元素法是基于数学上的变分原理，将所要求解的物理问题化为对泛函求极值的一个变分方 程；再利用差分法中的区域划分的离散化方法，并通过 元素划分所构造的插值函数，把求解连续的 变分方程问题离散化为求解线性方程组。按照这样的有 限元素方法来处理物理问题，就不再需要通 过建立偏微分方程这一道步骤，并且其物理问题在离散化的整个过程中就始终具有明确的物理意 义。 有限差分法来求解物理问题的数值解时，必须首先从物理模型出发，列出相应的偏微分方程及 定解条件，然后通过网格划分将偏微分方程的求解问题离散化为对差分方程组的数值求解。 对区域的离散化方法不同： 有限元素法采用的一般是三角形划分的方法。这样的划分对节点在区域内的配置方式比较任 意，其配置方式可以根据边界条件的情况来选择。这 样就可以在边界形状比较复杂时，仍然可以选 择边界节点完全处在区域的边界上，从而在边界上可以做到较好的逼近。特别是在由不同介电常数 的介质构成的静电场域内求解时，我们可以将节点取在 不同介质区域的交界面上，并在电位梯度较 大的区域，节点还可配置密一些，以实现较好的计算精 度。一般在有限元素法中采用三角划分时， 如果三角单元无限缩小，有限元近似解收敛于精确解， 近似解的平均导数也收敛于精确解的导数。 在有限差分法中，通常采用的是矩形网络区域划 分，因而很难实现网络节点在区域中的配置 与边界（不同介质界面）的良好逼近。 计算精度比较： 在有限元素法中，实践证明选择三角形不要太狭长，三角形元素边长比越接近 1，计算得到 的数值质量越好。有限元素法的计算精度于三角形划分时最大三角形边长 有关，若精确解有二阶 导数，则函数误差与 同阶，导数误差与 同阶（当 趋于无穷小时）。用有限元素法求解物理问题 时，它是用统一的观点对区域内的节点和边界节点列出 计算格式。这就使得各节点的计算精度总体 上比较协调。此外，有限元素法的计算格式中的矩阵 h 2 h h h ( )K 具有比较好的性质，即它是一个对称正定的 大型稀疏矩阵。这就给求解有限元方程组带来方便。 有限差分法在采用直交网络时其计算精度与矩形最大边长 有关，此时列出的计算格式比有限元 素法简单方便。而有限差分法则是孤立地对微分方程及定解条件分别列差分方程，因而各节点精度 总体上不够一致。 h 应用范围比较： 在对边界形状规则的求解区域，自然采用有限差分法就比较合适。事实上，有限差分法的适 用范围要比有限元素法广泛得多。有很多物理问题目前 还不能用有限元素法求解，但是人们总是可 以采用有限差分法。特别是在边界形状比较规则时，采用有限差分法是最合适的 数学形式以及度计算机的要求： 有限元素法的节点配置比较任意，计算格式要列出来就 要复杂得多。但是这些计算格式都可 以在电子计算机上自动形成，也容易将程序标准化，因而这并不会影响它的实际应用。但是有限元 素法要求的计算机内存量比较大，需要准备输入的数据量也比较大，这是它的缺点之一。