1、第3章 有限元法基础3.6 有限元方程组的求解有限元方程组的求解l 用有限元方法求解电磁场问题，在导出有限元方程组后，即归结为求解用有限元方法求解电磁场问题，在导出有限元方程组后，即归结为求解一个方程组的问题，方程组阶数等于未知节点参数的节点数。一个方程组的问题，方程组阶数等于未知节点参数的节点数。dAdKRAKII121I11常数时：线性方程组常数时：线性方程组常数时：非线性方程组常数时：非线性方程组3.6.1 有限元线性方程组解法有限元线性方程组解法l 高斯高斯-若尔当解法若尔当解法l 注意有限元方程组系数矩阵的特点：注意有限元方程组系数矩阵的特点：(1)对称性对称性(2)稀疏性与正定性稀
2、疏性与正定性第3章 有限元法基础下面主要介绍牛顿下面主要介绍牛顿-拉夫逊迭代法。拉夫逊迭代法。3.6.2 有限元非线性方程组解法有限元非线性方程组解法l 线性化方法：如逐次线性化方法，牛顿线性化方法：如逐次线性化方法，牛顿-拉夫逊迭代法及改进型的牛顿拉夫逊迭代法及改进型的牛顿-拉夫逊迭代法等。拉夫逊迭代法等。l 通过求解目标函数的极小值点来获得非线性方程组的解：如最速下降法通过求解目标函数的极小值点来获得非线性方程组的解：如最速下降法，共轭梯度法等。，共轭梯度法等。逐次线性化方法逐次线性化方法特点：比较简便，有良好的收敛性，但收敛速度慢，适特点：比较简便，有良好的收敛性，但收敛速度慢，适宜在求
3、解方程的阶数不是特别高时采用；宜在求解方程的阶数不是特别高时采用；牛顿牛顿-拉夫逊迭代法拉夫逊迭代法特点：收敛速度快，但在形成方程组时需要很大的计特点：收敛速度快，但在形成方程组时需要很大的计算量，并要求有很好的初值，否则迭代过程可能不收敛。为了克服这个缺算量，并要求有很好的初值，否则迭代过程可能不收敛。为了克服这个缺点，但又能发挥其特点，产生了多种点，但又能发挥其特点，产生了多种改进型的牛顿改进型的牛顿-拉夫逊迭代法拉夫逊迭代法。1.牛顿牛顿-拉夫逊迭代法拉夫逊迭代法第3章 有限元法基础对矢量磁位对矢量磁位 A，使用直角坐标系，使用直角坐标系，RKA F AKAR 000-AFF AF AA
4、 AA 000F AF AJ AAA 00100F AJ AAA 11000AAJ AF A一般认为一般认为 A 是近似解，对于真解必有一个不为零的余矢量是近似解，对于真解必有一个不为零的余矢量在此方程中，在此方程中，K 和和 F 都是都是 A 的函数，如果的函数，如果 A(0)是是 的近似解，将的近似解，将F(A)用泰勒级数表示，取一次项，则有用泰勒级数表示，取一次项，则有若用若用 J(A(0)表示表示 F(A)的雅可比矩阵在的雅可比矩阵在 A(0)处的值，则有处的值，则有RKA 令此余矢量对于令此余矢量对于 A(1)解时为解时为0，有，有第3章 有限元法基础11KKKKAAJ AF A式中
5、式中11112222()1212KKKnKKKKnKKKnnnnFFFAAAFFFAAAFFFAAAJ A求余矢量的雅可比矩阵，需要从一个单元的余矢量函数讨论起。求余矢量的雅可比矩阵，需要从一个单元的余矢量函数讨论起。一个单元的余矢量为：一个单元的余矢量为：eeeeeiiiijjimmiieeeeeeeeejiijjjjmmjjeeeeemmiimjjmmmmKAKAKARfKAKAKARffKAKAKAR FKAR第3章 有限元法基础单元雅可比矩阵为：单元雅可比矩阵为：eeeiiiijmeeeeejjiijmeeemmmijmfffAAAfffAAAfffAAAFJAA以以 为例，看雅可比矩
6、阵第为例，看雅可比矩阵第 ii 项元素对单元项元素对单元 e 的贡献。的贡献。K 是是 A 的函数的函数，对对A的偏导数均不为零，由的偏导数均不为零，由(3.66)式可知式可知 ieiAf eijK eimK ststestccbbk41 4esttstskbbc c1因此有：因此有：第3章 有限元法基础 11144eeststststststiiiiKb bc cb bc cKAAAA mieimjieijiieiieiiieiAAKAAKAAKKAf 1eeeeeiiiiiiijjimmiifKKAKAKAA A 1eeststiiKKAA下面求下面求 。由以下各式：。由以下各式：iAii
7、ABBA12iijjmmiiiijjjjmmmmAN AN AN Aab xc y Aab xc y Aab xc y A22yAxAB第3章 有限元法基础可以求出：可以求出：22122211122221 141iiiiiiiiijjmmiijjmmiiijijjimimmBAAAAAB ABBxAxyAybcb Ab Ab Ac Ac Ac AB BbcAbbccAbbccAB BB 11 1eeeiiiijjimmeiKAKAKABgB B eeeeiiiiijjimmgKAKAKA式中式中所以所以 2e2()11 11eeeeeiiiiiiiiifgKggKAB BBB同理可求出：同理可
8、求出：2eeeeeijjiijjiggffKA BBA st21,eeeeststggJKs ti j mBB kNeeekAJAJ1 11kkkkAAJ AF A单元雅可比矩阵元素统一表示为：单元雅可比矩阵元素统一表示为：（、B 对各单元均不相同）对各单元均不相同）求出每个单元的雅可比矩阵求出每个单元的雅可比矩阵 J(e)和余矢量和余矢量 F(A(k)后，将所有单元的分析结后，将所有单元的分析结果进行总体合成，得到：果进行总体合成，得到：则牛顿则牛顿-拉夫逊迭代公式可以归结为下式：拉夫逊迭代公式可以归结为下式：第3章 有限元法基础第3章 有限元法基础KKKKKKBBBBHHHH1112112
9、1KKKKKKHBH BHHBHBBBBBBBB 11kkkkAAJ AF A（3.185)在上式中，对于在上式中，对于J(e)中各元素，必须在求出中各元素，必须在求出 后才能进行计算，因此，由矢后才能进行计算，因此，由矢量磁位量磁位 A 求解时，由每次迭代解求解时，由每次迭代解 A(K)计算出单元磁通密度计算出单元磁通密度 B，由，由 B 值查磁化曲值查磁化曲线得对应线得对应 H 值，再计算值，再计算 ，磁化曲线的查取归结为数值计算。，磁化曲线的查取归结为数值计算。对对 B=f(H)的数学处理，一般包括插值法和拟合法，比较简单而有效的办法的数学处理，一般包括插值法和拟合法，比较简单而有效的办
10、法是线性插值和拉格朗日插值等。是线性插值和拉格朗日插值等。BB当磁场磁化曲线采用逐段线性插值函数逼近，当磁场磁化曲线采用逐段线性插值函数逼近，H 值在值在 HK+1与与 HK之间时，有之间时，有第3章 有限元法基础2.改进型的牛顿改进型的牛顿-拉夫逊迭代法拉夫逊迭代法牛顿牛顿-拉夫逊迭代法的优点：拉夫逊迭代法的优点：（1）收敛速度快，按平方律收敛，每经一次迭代有效位数基本上增加一倍；）收敛速度快，按平方律收敛，每经一次迭代有效位数基本上增加一倍；（2）自校正功能，即）自校正功能，即 A(K+1)仅依赖于仅依赖于 F(A)与与 A(K)，前面迭代的舍入误差不，前面迭代的舍入误差不会一步步传递下去
11、。会一步步传递下去。其缺点：其缺点：（1）每次迭代都要形成一次）每次迭代都要形成一次 J(K)，而，而 J(K)的计算时间往往比增加迭代的次数的计算时间往往比增加迭代的次数所需要的时间多；所需要的时间多；（2）对初值要求较高，选择不当，常会引起振荡。）对初值要求较高，选择不当，常会引起振荡。为了克服所存在的缺点，现已形成多种改进型牛顿为了克服所存在的缺点，现已形成多种改进型牛顿-拉夫逊迭代法。拉夫逊迭代法。第3章 有限元法基础(1)修正的牛顿修正的牛顿拉夫逊迭代法拉夫逊迭代法 为了克服每次迭代都需要计算为了克服每次迭代都需要计算 J(K)，采取全部求解过程都用，采取全部求解过程都用J(0)，即
12、每次迭，即每次迭代不需要形成新的系数矩阵，保持其斜率不变。这将使迭代次数增加而影响代不需要形成新的系数矩阵，保持其斜率不变。这将使迭代次数增加而影响收敛速度，但如初值选取较好，总体上将节省运算时间。收敛速度，但如初值选取较好，总体上将节省运算时间。(2)采用欠松弛因子的方法采用欠松弛因子的方法欠松弛迭代的方法为欠松弛迭代的方法为式中式中为收敛因子，为收敛因子，01；A(K)为每次迭代近似解的误差。为每次迭代近似解的误差。此外，还有此外，还有牛顿牛顿-拉夫逊迭代法与最速下降法的组合求解法拉夫逊迭代法与最速下降法的组合求解法，加速的全牛顿法加速的全牛顿法，采用阻尼因子采用阻尼因子等，具体方法请参阅
13、文献等，具体方法请参阅文献2。KKKAAA1否否读入所需系数及数组读入所需系数及数组建立右端矢量建立右端矢量 R建立雅可比矩阵及右端剩余矢量建立雅可比矩阵及右端剩余矢量 R计算迭代误差计算迭代误差计算计算的初值的初值计算各单元计算各单元 b，c，是是开始开始牛顿拉夫逊迭代牛顿拉夫逊迭代迭代误差是否迭代误差是否小于控制值？小于控制值？B结束结束l 非线性方程组非线性方程组计算框图：计算框图：根据边界条件修正方程组根据边界条件修正方程组计算各单元计算各单元 A，B解方程组解方程组得到各节得到各节点点 A计算各单计算各单元元及及3.8 有限元素的自动剖分有限元素的自动剖分l 采取自动剖分的必要性：采
14、取自动剖分的必要性：第3章 有限元法基础l 在单元剖分中，为了尽可能的压缩存储，减小计算量，提高精度，必须在单元剖分中，为了尽可能的压缩存储，减小计算量，提高精度，必须注意以下问题：注意以下问题：(1)三角形各边不要相差太悬殊，避免出现尖锐的三角形；三角形各边不要相差太悬殊，避免出现尖锐的三角形；(2)一个节点周围，不宜集中过多的三角形单元，为压缩存储创造条件；一个节点周围，不宜集中过多的三角形单元，为压缩存储创造条件；(3)三角形单元内物理参数三角形单元内物理参数(或或)变化连续，即媒质交界面应与单元的边变化连续，即媒质交界面应与单元的边界重合；界重合；(4)精度要求不同的区域，元素的密度应
15、不同；精度要求不同的区域，元素的密度应不同；(5)元素节点编排应规格化。元素节点编排应规格化。本节讨论平面内的几种自动剖分方法。本节讨论平面内的几种自动剖分方法。第3章 有限元法基础3.8.1 直线内插法直线内插法l 适宜于对以直线段为边界的的场域进行自动剖分。只需给出适宜于对以直线段为边界的的场域进行自动剖分。只需给出 x,y方向上两方向上两端点的坐标，就可算出所有节点的坐标。端点的坐标，就可算出所有节点的坐标。1.确定确定 x 方向、方向、y方向节点数及总方向节点数及总节点数节点数x方向节点数方向节点数1 nNx1 mNy0(1)(1)Nnmy方向节点数方向节点数节点总数节点总数2.确定各
16、节点的坐标确定各节点的坐标x 和和 y 方向节点坐标最小值为方向节点坐标最小值为x1、y1，最，最大值为大值为xm、ym，则节点坐标增量分别为，则节点坐标增量分别为nxxm1myym1，第3章 有限元法基础可推得第可推得第 Ni 列列 Nj 行的第行的第 Nk 个节点的编号为个节点的编号为kiyj(1)NNNN式中式中yjxiN,NN,N2121111111jmkimkNmyyyNyNnxxxNx第第 Nk 个节点的坐标值为个节点的坐标值为例例3.1 设设 x1=0,y1=0,xm=7cm,ym=5cm，取，取n=7,m=5,可以求出可以求出x方向节点数方向节点数 Nx=7+1=8，y方向节点
17、数方向节点数Ny=5+1=6，总节点数，总节点数 N0=86=48。求第。求第2列第列第3行行节点数编号及坐标。节点数编号及坐标。节点编号节点编号(Ni=2,Nj=3)，按，按(3.192)式，有式，有 Nk=(2-1)6+3=9 其坐标按其坐标按(3.193)式为式为（3.193)7-0 x(9)=0+(2-1)cm=1cm75-0y(9)=0+(3-1)cm=2cm53.确定三角形单元编号及各三角形单元确定三角形单元编号及各三角形单元3顶点顶点(节点节点)按规定顺序的的总体编号按规定顺序的的总体编号第3章 有限元法基础对于图对于图(a)中类型的三角形单元，其单中类型的三角形单元，其单元编号
18、元编号 E 与行、列的关系为与行、列的关系为 E2(Ni-1)(Ny-1)+Nj 式中式中 Ni1，2，3，(Nx-1);Nj 1，2，3，(Ny-1)。1)()()()(12)1(2)(EJEMNEIEJNN/NNEIyjyyi图图(a)中第中第 E 个三角形单元的个三角形单元的3节点，节点，i,j,m的总体编号与行、列的关系分别为的总体编号与行、列的关系分别为同理，对图同理，对图(b)类型的三角形单元有类型的三角形单元有第3章 有限元法基础对于图对于图(b)中类型的三角形单元，其单元编号中类型的三角形单元，其单元编号 E 与行、列的关系为与行、列的关系为 E2(Ni-1)(Ny-1)+Nj
19、 式中式中 Ni1，2，3，(Nx-1);Nj 1，2，3，(Ny-1)。3节点节点 i,j,m的总体编号与行、列的关系分别为的总体编号与行、列的关系分别为1)()()()(2)2(2)(EJEMNEIEJN/NNEIyjyi例例3.2 (1)求上例中求上例中1列列3行（行（b）类型单元的编号）类型单元的编号 E 和和3个节点个节点 i,j,m 的总体节的总体节点编号。点编号。(2)求出求出1列列3行（行（a）类型单元的编号）类型单元的编号 E 的的3个节点个节点 i,j,m的总体节点的总体节点编号。编号。l 当场域形状复杂，要求场域内单元密度分布不均匀，又希望密度是平滑过当场域形状复杂，要求
20、场域内单元密度分布不均匀，又希望密度是平滑过渡时，这种简单的直线内插法就不适用了，广泛应用的是渡时，这种简单的直线内插法就不适用了，广泛应用的是等势剖分等势剖分的方法。的方法。第3章 有限元法基础3.8.2 等势剖分等势剖分1962年，年，W.P.Crowley1.微分方程的推导微分方程的推导等势剖分是把网格线看成平面上的两个势函数的等势线，定义两个函数等势剖分是把网格线看成平面上的两个势函数的等势线，定义两个函数 (x,y)及及(x,y)，它们在整个计算区域内部满足拉普拉斯方程：，它们在整个计算区域内部满足拉普拉斯方程：000022222222yyxxyyxxyxyx(3.198)第3章 有
21、限元法基础当边界当边界 和和s 给定后，可以解出场域内给定后，可以解出场域内 和和的分布，求出等势线。的分布，求出等势线。s,1,2,3.,ix yin,1,2,3.,jx yjml 两族曲线的交点就是剖分所要确定的节点，但是得到的节点坐标是两族曲线的交点就是剖分所要确定的节点，但是得到的节点坐标是 和和 值，我们需要的是等势线交点在值，我们需要的是等势线交点在(x,y)平面上的坐标！平面上的坐标！若将若将 x,y看做看做(,)平面上的函数平面上的函数 x(,)及及 y(,)，通过函数的变换，方，通过函数的变换，方程式程式(3.198)可化为所要求解的形式，简单步骤如下：可化为所要求解的形式，
22、简单步骤如下：设设(,)(,)xxyy 两端对两端对 x 求导求导dd1ddddddddxxxxyyyxxx第3章 有限元法基础可解得可解得yxyxyyxxxdd11ddJxxxyxyx yx yyy 令令可得可得ydxdyxJx11ddxxxddd11ddyyxJx(3.204)第3章 有限元法基础(,)(,)x yx y又因为又因为两端对两端对 x 求导求导ddddddddyxxyxyxxyxd1yJxd1xJyd1yJxd1xJy将式将式(3.204)两式展开，与上式比较可得两式展开，与上式比较可得(3.208)(3.209)(3.210)(3.211)第3章 有限元法基础将式将式(3.
23、208)、(3.209)代入代入 ，则有，则有22220 xy11x0 xdyxxyyJyJ22111dxyJyyxJJxJ 其中其中由式由式Jxyxyx yx y可得可得yxyyxxxxyxxxJxyxyxxxxxxyyxx(3.213)(3.212)(3.214)将式将式(3.213)中中 展开为展开为第3章 有限元法基础2xy xyxyyxxy202212222222223xxxyyyyJ2x02212222222223xxxyyyyJ2x将式将式(3.203)、(3.208)、(3.209)、(3.210)、(3.211)、(3.213)、(3.214)及及 表达式代入式表达式代入式(
24、3.212)，化简可得，化简可得(3.215)Jy同理可以导出同理可以导出(3.216)(3.217)第3章 有限元法基础其中其中22222222yxyxyyxxyyxxyxyx齐次方程式齐次方程式(3.216)、式、式(3.217)有唯一解的条件为有唯一解的条件为 J0，因此有，因此有0202yxyxxxl 若在若在(,)平面得到式平面得到式(3.219)的两组数值解的两组数值解x(,)和和y(,)，则对于定域，则对于定域内每一个确定的内每一个确定的 、i，就有一对就有一对 xi 和和 yi，这就是要求的某两条网格线，这就是要求的某两条网格线相交处的节点坐标。相交处的节点坐标。(3.218)
25、(3.219)il 在在(,)平面上网格与节点数与平面上网格与节点数与(x,y)平面上是相等的，但在平面上是相等的，但在(,)平面上平面上是等距的网格，在是等距的网格，在(x,y)平面上却是用相应边界确定的，各单元密度过渡是平平面上却是用相应边界确定的，各单元密度过渡是平滑的。但并不是等距网格。滑的。但并不是等距网格。第3章 有限元法基础2.等势剖分计算公式的推导等势剖分计算公式的推导图图3.19 网格与边界的划分网格与边界的划分x y平面区域平面区域逻辑网格逻辑网格（1）边界值的确定）边界值的确定第3章 有限元法基础（2）三角形网格的逻辑坐标）三角形网格的逻辑坐标在空间坐标中用三角形网格剖分
26、时，在空间坐标中用三角形网格剖分时，在在(,)坐标中，等坐标中，等 线和等线和等 线有线有多种取法，常用的有如图多种取法，常用的有如图 3.20 所示所示3种。种。图图3.20 不同的逻辑网格方式不同的逻辑网格方式(a)(b)(c)（3）对等边三角形逻辑网格求）对等边三角形逻辑网格求(3.219)式的差分方程式的差分方程0202yxyxxx第3章 有限元法基础令等边三角形令等边三角形 ,，推导的关键是对三角形网格求出，推导的关键是对三角形网格求出 ，。lmxxxyyy利用积分中值定理有利用积分中值定理有设函数设函数(,)ff SSffddddSSffddddd ddd ddSLSLff ffd
27、dLLLLLddddLfff 格林公式格林公式同理，有同理，有第3章 有限元法基础对于二次积分，同样应用积分中值定理和格林公式，可得对于二次积分，同样应用积分中值定理和格林公式，可得ffddLLLLLddddLfff LLLddddLfff现把现把 替换成替换成 x(,)或或 y(,)，可以计算出每个节点，可以计算出每个节点(l,m)的的 ，与与 。(,)f xxyy第3章 有限元法基础图图3.21 等边三角形节点编号等边三角形节点编号如图示，节点如图示，节点(l,m)由由6个等边三角个等边三角形包围，求形包围，求 的回路积分是沿的回路积分是沿6个个三角形的外边界积分，即三角形的外边界积分，即
28、x 6111611161112163456i111d2d21226iiiiiliiiiiiiiiiiixxmmx mxxlll mmmx mmxxxxxxl mm第3章 有限元法基础同理可得同理可得 6111652341112266LmiiiiLxdlxxx llxxxxxxldm ，与上面求法完全相同，即将与上面求法完全相同，即将 xi 换成换成 yi 以上式子对以上式子对 y 全部成立。由此全部成立。由此可进一步求出可进一步求出，。yy为求出为求出 ，要利用图要利用图3.21中虚线网格，即对质心构中虚线网格，即对质心构成包围成包围(l,m)节点的六边形进行积分。节点的六边形进行积分。xxx
29、yyyddddddddddddlLLLLLLLLLLLLxx mxl mxxxxxx 第3章 有限元法基础首先求出各三角形内平均的首先求出各三角形内平均的 和和 。xx3111l3111xjjjjjjjjxmmxlmm利用上式利用上式，求出图，求出图3.21式中对应三角形的式中对应三角形的 。lx11213234454665106102203304405506llllllxxxxxxxxxxxxxxxxxx第3章 有限元法基础再沿图中虚线做环路积分，求出再沿图中虚线做环路积分，求出142llxxxxx对应图示各点对应图示各点(m,l)：1(m+1/3),(l+1/3)I m,(l+1/2)2(
30、m-1/3),(l+1/3)II (m-1/2),(l+1/2)3(m-2/3),(l+1/3)III (m-1/2),l 4(m-1/3),(l-1/3)IV m,(l-1/2)5(m+1/3),(l-2/3)V (m+1/2),(l-1/2)6(m+2/3),(l-1/3)VI (m+1/2),l xxxxxxxxxxmllm221524361362xxxxxmm同理可得同理可得第3章 有限元法基础yyy同理可得到同理可得到 ，类似的关系。类似的关系。0202yxyxxx将它们代入式将它们代入式(3.219)：整理后可得整理后可得 14253614253622202220 x-xx xxx
31、xxxyyy yyyyyy上式整理后可得上式整理后可得第3章 有限元法基础 12345612345600 xx xxxxx-x xxxxyy yyyyy-y yyyy归纳成统一形式为：归纳成统一形式为：616100iiiiiiyywxxw式中式中 wi 是是，和，和 的函数。上面得到的值是对等边三角形的情况。对于的函数。上面得到的值是对等边三角形的情况。对于不同的逻辑网格，它们取值将不同。不同的逻辑网格，它们取值将不同。(3.236)3.差分方程组的求解差分方程组的求解第3章 有限元法基础方程式方程式(3.236)是非线性代数方程，是非线性代数方程，n个内部节点有个内部节点有2n个坐标，有个坐
32、标，有2n个方程，个方程，一般采用一般采用高斯高斯-赛德尔赛德尔选代。计算的步骤为：选代。计算的步骤为：(1)边界顶点坐标已知，其它节点赋初值。边界顶点坐标已知，其它节点赋初值。(2)从左下角开始计算从左下角开始计算 J，和，和。(3)计算计算(3.236)式，求式，求 x,y，经过多次迭代，最后得到各节点坐标值。，经过多次迭代，最后得到各节点坐标值。第3章 有限元法基础3.8.3 自适应剖分技术自适应剖分技术第3章 有限元法基础l 自适应技术是靠网格细分与场量计算的循环过程来实现，具体步骤为：自适应技术是靠网格细分与场量计算的循环过程来实现，具体步骤为：(1)生成开端网格，即采用网格自动细分
33、软件将场域剖成很粗的网格；生成开端网格，即采用网格自动细分软件将场域剖成很粗的网格；(2)求解场量，即形成有限元方程与求解方程；求解场量，即形成有限元方程与求解方程；(3)分析场结果计算误差分析场结果计算误差(对每个单元进行误差分析对每个单元进行误差分析)；(4)根据误差分析，确定需要细分的网格单元；根据误差分析，确定需要细分的网格单元；(5)细分局部网格；细分局部网格；(6)返回到返回到(2)，求解场量；，求解场量；(7)计算最后结果。计算最后结果。第3章 有限元法基础否否场域信息输入场域信息输入是是精度满足精度满足要求？要求？开端网格生成开端网格生成场量计算场量计算误差分析误差分析网格自动细分网格自动细分后处理后处理图图3.22 自适应剖分框图自适应剖分框图网格初始剖分网格初始剖分 采用自适应技术后网格一次细分采用自适应技术后网格一次细分 网格二次细分网格二次细分 第3章 有限元法基础
