工程电磁场数值分析4(有限元法1).ppt

工程电磁场数值分析（4） （电磁场有限元法） 华中科技大学电机与控制工程系 陈德智 Email: Tel:Office: Room 108, 电机楼 2010.12 第4章 电磁场有限元法 （Finite Element Method, FEM） 有限元法可以基于变分原理导出，也可以基于加权 余量法导出，本章以加权余量法作为有限元法的基础， 以静电场问题的求解为例介绍有限元法的基本原理与实 施步骤。并介绍有限元法中的一些特殊问题。 第4章 电磁场有限元法（FEM） 有限元基本原理与实施步骤：1D FEM 有限元基本原理与实施步骤：2D FEM 有限元方程组的求解 二维有限元工程应用 三维有限元原理与工程应用 矢量有限元 加权余量法回顾： 对算子方程 用 作为该方程的近似解（试探解）： 代入方程得余量： 1. 有限元法基本原理与实施步骤：一维问题 在有限元法中，基函数一般用 表示 。采用Galerkin方案，取权函数与基函数相同。使与余量正 交化： 设L为线性算子，代入 ，得 或 记 得代数方程组： 加权余量法回顾（续） 利用有限元法求解一维边值问题： （1）单元剖分 如图5个单元，6个节点 （2）选取基函数 （3）方程离散 （计算系数阵 K 和右端项 b） 基函数 Ni 只是一阶可导 的，不能严格满足微分方 程，称为“弱解”。 （3）方程离散 第一项在 xj 处为0，在 xi 处的值 被来自 (i-1) 单元的贡献抵消，故只剩下第二项。 由于基函数 Ni 局域支撑，显见只有 不为0。 使用分步积分： （3）方程离散 故 类似，当 j = i 时 右端项： 总体方程 强加边界条件：u1 = 0, u6 = 0 （4）求解方程 思考：（1）有限元的解跟有限差分法的解有何根本不同？ （2）有限元的系数阵总是对称的吗？ x 000 0.20.03610.0360 0.40.06280.0625 0.60.07100.0708 0.80.05250.0523 1.000 与有限差分法（FDM）相比，有限差分法是对点的 离散，得到一系列离散点上的解；而有限元（FEM ）是对区域的离散（单元），尽管所求的是节点上 的自由度，但它的解在场域中每一个点上都有定义 。 所以，即是有限元节点上的解是精确的，有限元的 整个解仍然是近似的。好的数据处理技术可以从该 近似解中提取更精确的分析结果。 线性单元中，如果所求的自由度是电位j，单元中 的电场 E是场量；节点上的 E 取邻近单元的平均。 一些补充说明： 关于有限元的解 计算系数阵是有限元分析的主要工作量。所涉及到的积 分，如果不是解析可积的，通常要用到数值积分。其中 最常用的数值积分方法是Gauss数值积分。 一些补充说明： 高斯数值积分 先将积分区间变换到 -1,1上；按照固定的 积分点计算若干函数 值 P(xi), 以固定权值 wi 累加即可。具 (2n+1)阶精度。 n=4 x(1)= 0.861136311594053d0 x(2)= 0.339981043584856d0 w(1)= 0.347854845137454d0 w(2)= 0.652145154862546d0 n=5 x(1)= 0.906179845938664d0 x(2)= 0.538469310105683d0 x(3)= 0.0d0 w(1)= 0.236926885056189d0 w(2)= 0.478628670499366d0 w(3)= 0.568888888888889d0 n=6 x(1)= 0.932469514203152d0 x(2)= 0.661209386466265d0 x(3)= 0.238619186083197d0 w(1)= 0.171324492379170d0 w(2)= 0.360761573048139d0 w(3)= 0.467913934572691d0 n=16 x(1) = 0.9894003948d0 x(2) = 0.9445750231d0 x(3) = 0.8656312024d0 x(4) = 0.7554044084d0 x(5) = 0.6178762444d0 x(6) = 0.4580167777d0 x(7) = 0.2816035508d0 x(8) = 0.0950125098d0 w(1) = 0.0271524594d0 w(2) = 0.0622535239d0 w(3) = 0.0951585117d0 w(4) = 0.1246289713d0 w(5) = 0.1495959888d0 w(6) = 0.1691565194d0 w(7) = 0.1826034150d0 w(8) = 0.1894506105d0 一些Gauss积分点和权值： （关于x=0对称，只给出一半） 为提高有限元分析精度，有两种方法： 其一：增加节点，细化网格称为h方法。 其二：增加有限元的阶数称为p方法。 一些补充说明： 线性单元与高阶单元 一些补充说明： 二阶单元 一些补充说明： 三阶单元 h方法和p方法的求解精度 By Jianming Jin. The Finite Element Method in Electromagnetic s, 2nd Ed., 2002 作业：要独立完成，凡雷同者没分！ 编写有限元程序，计算一维边值问题。改变剖分单元 数目，观察解的精度变化。（建议也同时做一个有限 差分法的程序，比较二者的精度差别） 以二维静电场泊松方程的求解为例 。 2. 有限元法基本原理与实施步骤： 二维问题 目标：依据加权余量法，利用分域基，建立离散的代数 方程组，即确定系数Kij 和bi。 场域离散 二维问题常使用三角形单元离散，便于处理复杂的场域 形状，容易实现。 单元：互不重叠，覆盖全部场域；每个单元内介质 是 单一、均匀的。 节点：网格的交点，待求变量的设置点。 该步骤需要记录的信息： 节点编号、节点坐标 节点属性(激励源、是否边界等) 单元编号 单元节点编号 单元介质 基函数 有限元采用分片逼近的思想 ，类似于一维情况下使用折线 逼近一条任意曲线。 使用分域基Ni，基函数的个数 等于节点的个数；每个基函数 Ni的作用区域是与该节点i相 关联的所有单元。 三角形单元内的基函数 设三角形三个顶点处待求函数值 分别为u1, u2, u3。如果单元足够 小，可以采用线性近似，将单元 内任意p点的u(x,y)表示为 代入三个顶点的坐标和函数值， 可以解出a、b、c。得到 单元节点的编号按 逆时针方向排列！ 其中， 记住我们的任务 寻找基函数 对比 可得 基函数Ni常被称为插值函数或者形状函数，具有以下性质： （1）是插值的； （2） （3）在相邻单元的公共边界上， Ni是连续的，从而通过Ni构造的逼近函数也是连续的 。 在积分 中，对于确定的 i，j的有效取值为i 本身以及与节点i相联的周围节点，积分的有效区域为以i、 j 为公共节点的所有三角形单元 ，在这些单元中Ni、Nj才有 交叠。 计算系数阵 这些积分可以分单元进行。例如对 右图所示的局部编码，K01、K00以 及b0的计算公式为： 计算系数阵 以下把单元e的贡献记为 这样，就有 每个 或 的计算都在具体的单元内单独考虑（ 称为单元分析）。 单元分析：计算单元内积分对系数阵和右端项元素的贡献 。 系数阵元素： 当L为拉普拉斯算子时，由于Ni在 单元内是(x, y)的线性函数，经 Laplace算子作用后值为0。 但是，在相邻单元的边界上， Ni是 连续但是不光滑的，因此对积分的 贡献主要来自边界。 为考虑单元边界的影响，需要借助 于格林公式： 故 ， 格林公式： 因： 写成一般形式，若一个三角形三个顶 点编号为i, j, m（逆时针顺序），则 从而 再看边界部分： （1）在节点 i 的对边Gjm上，Ni0，故积分贡献为0； 结论：单元边界对积分的贡献为0。 所以单元e对系数阵元素的贡献为： （2）在节点 i 的邻边Gij上，由于计 算Kij时需要把具有公共邻边的单元 的积分累加，此二单元的Ni是连续 的；对于单一均匀媒质，要求相邻 单元满足 ，故积分 的贡献相互抵消。 由于单元很小，做单元分析时通常可以取 f (e) 为常数值 （可以认为等于三个顶点上的平均值）。因此 右端项元素： 公式： 上述以节点为序的分析过程对于有限元原理的说明是易 于理解的。而在实际编程中，更有效率的是以单元为序 ，逐个计算单元系数阵K(e)，然后合成整体系数阵K。 单元系数阵K(e)定义为 设 i, j, m 是节点的整体编号，元素Kij在整体矩阵中的实 际位置是第i行、j列；因此 必须合成到整体矩阵的第 i行、j列元素上。 单元矩阵： 整体矩阵合成： 通过上述过程，对于一个“正常”的内部节点就建 立起了一个代数方程。“非正常”的节点包括：媒 质交界面衔接条件和场域边界条件。 对于静电场问题，媒质分界面衔接条件为 媒质交界面衔接条件 第一个条件是自动满足的（Why？），无须格外处理。 对于第二个条件，前面计算单元边界上积分 时，默认两边 u 的法向导数相等，使内边界上的积分结 果抵消。因此只要把泊松方程写成 或 满足的条件将是 , 从而也无需另行处理 。 由于有限元方法能够自动满足媒质交界面条件， 因此有限元法特别适合于处理多层复杂媒质问题 。这是其它方法无可比拟的。 媒质交界面衔接条件 第一类边界条件（强加边界条件） 第一类边界节点是指边界上函数值 已知。因此处理 方法是，合成整体系数阵之后，将该节点所在行的主元素置 1，其它元素均置零，同时将右端项中对应元素设为已知函 数值。 要保持对 称性；有 更简便的 做法 第二类边界条件（自然边界条件） 第二类边界节点是指边界上函数法向导数 已知 。对于内部单元，相邻单元边界的积分 相互抵 消。但是对于场域边界，如果给定第二类边界条件不为0， 则积分结果要计入右端项中。但是若给定的是齐次第二类边 界条件，则积分结果为0，无需另行处理，非常方便。 这是ANSYS中自动满足的边界条件。 有限元方法的推导过程虽然看起来有些复杂，但是最终 结果是非常简单而且优美的。因为边界条件的处理和媒 质交界面条件的处理都非常方便，使得有限元方法在处 理复杂媒质问题和复杂场域问题时得心应手，获得了广 泛的应用，成为最重要的数值分析手段，广泛应用于各 个领域。有人用“功盖四方”来形容有限元，实不为过。 中国人在有限元的发明中有自己独特的贡献。 作业： (2) 对于研究方向为数值计算的同学： 编写一个二维静电场有限元程序， 计算右图所示问题，或其它自己找 一个问题。 （1）推导三角形单元的2次和3次插值函数。 3. 有限元方程组的求解 代数方程组求解方法概述 所有的数值方法最终都归结为求解一个代数方程组： 系数阵 A也称系统矩阵或刚度矩阵。不同离散方法得到的 系统矩阵具有不同的特点，方程组的解法也就不同。 基于微分方程（如FEM、FDM等）得到的系统矩阵是稀 疏的，有时还是对称的； 而基于积分方程得到的系统矩阵则是稠密的，如BEM、 模拟电荷法等。 代数方程组的求解是数值计算（计算数学）研究的核心 内容。求解代数方程组的方法归纳起来有两类：直接法 和迭代法。 3. 有限元方程组的求解 直接法：直接法都是基于高斯消去法，经过确定次数的 运算，理论上可以得到方程组的精确解。适用于小型、 稠密方程组的计算。 迭代法：是一种间接方法，从某个预定的初值出发，按 照一定的迭代步骤，逐渐逼近方程组的真解。得到一个 满足给定精度要求的近似解。适用于大型、稀疏方程组 的计算。 3. 有限元方程组的求解 直接法（LU分解算法） LU分解算法： 回带： 消元： 计算量：需要的 乘除法次数： O(n3) 稳定性：选主元 迭代法 迭代法的基本思想： （等价方程组） 从一组猜测的初值 开始迭代 直至 不再变化为止，即为方程组的解（收敛）。 好的迭代法应该：对初值不敏感；收敛速度快。 例如：高斯赛德尔迭代（有限差分法常用） 高斯赛德尔迭代 实际运算过程： 。 这就是通常的高斯塞德尔迭代格式，矩阵中的零元素不参与运 算，矩阵甚至可以不出现。大大减少了内存需求量和计算量。 共轭梯度法（Conjugate Gradient Method, CG法） 共轭梯度法在原理上可以通过n步迭代得到方程的准确解 ，因而也称为半直接法或半迭代法。 把迭代法表示为更一般的形式： 称为步长，p 称为搜索方向，用r 表示残差。 预优共轭梯度法（Preconditioned Conjugate Gradient Method, PCG法） 当系数阵的特征值较均匀地分布在一个很长的区间上时，称 矩阵具有很大的条件数；此时共轭梯度法的收敛速度可能很 慢。解决的办法是选取合适的非奇异矩阵C进行处理： 矩阵 称为预处理矩阵。 预优矩阵M应具有如下特性：稀疏性与A相近；矩阵 的特征值分布集中；形如 的方程组容易求解；易 于寻找。目前公认有效的方法是对系数阵A做不完全 Cholesky分解，以M = LDLT 为预优矩阵。这种方法称为不 完全分解预优共轭梯度法（Incomplete Cholesky decomposition preconditioned Conjugate Gradient Method， ICCG法）。 ICCG法 稀疏矩阵技术 没有稀疏矩阵技术就没有有限元的成功。 带状矩阵技术 通过合适的节点编码 ，可使系统矩阵的非零 元素集中于主对角线附 近的带形区域内。在使 用LU分解法求解方程组 的过程中，带形区域以 外的0元素无需计算。 缺点：节点编码优化； 带形区域内仍然有大量 零元素 非零元素存储技术 只存储矩阵的非零元素。一般用一个一维数组 存储矩阵的非零元素，用另一个索引数组存储这 些非零元素在原矩阵中的行和列值。例如： 非零元素存储技术 在迭代法中，系统矩阵参与的运算只有矩阵左乘 以某个相量。如果采用上面的存储方法，则实现 q = A p 的算法如下： C Compute the matrixvector product DO k=1, not i=IA(k) j=JA(k) q(i)=q(i)+A(k)*p(j) ENDDO 非零元素存储技术 下列的二维存储方案是我的发明：直接将矩阵 挤扁。 数组A有n行m列，m是A各行非零元素个数的最大值。对角 元素放在该行的第一个位置上，且不管是否为0，都要存贮 。用一个整型数组JA存放各元素的列号。JA中第一列可以 用来记录每行的非零元素个数。 本节更多的参考文献： 金建铭. 电磁场有限元方法，西安电子科技大学出版社， 1998 徐树方，矩阵计算的理论与方法，北京大学出版社， 1995 杨绍祺，谈根林，稀疏矩阵算法及其程序实 现 ,高 等教育出版社, 1985 刘万勋, 刘长学, 华伯浩等，大型稀疏线性方程组的解法. 国防工业出版社, 1981 R P 梯华森, 稀疏矩阵，科学出版社，1981 有限元分析精度的影响因素 静电场问题 静磁场问题 涡流问题 波的传输与散射 4. 二维有限元的工程应用 有限元分析精度的影响因素 （1）数学模型对工程问题的近似； （2）材料电磁参数的不确定性； （3）数学模型的有限元近似： 场域拟合精度单元大小、未知数个数与局部 场的变化情况； 边界拟合精度曲线边界； 系数阵计算过程中数值积分精度； 方程求解精度、数字误差； 计算结果的数据处理；好的处理技术可以提高分 析精度。 h方法、p方法；高次单元 单元阶数与计算精度的关系 By Jin。 曲边三角形单元：更好地拟合曲面边界 数值积分 三角元的高斯数值积分 单元比较小的时候，单元内的函数可以近似认为是常数 ，通常可以获得满意的精度。当单元内函数变化比较快 ，或者采用曲边三角形单元、高次单元时时，都会用到 数值积分。 高斯-勒让德数值积分公式： L1、L2、L3 是位置的面积坐标。 权值 wi 如下页表格。 三