电磁场有限元分析PPT课件.ppt

工程电磁场数值分析 4 电磁场有限元法 华中科技大学电机与控制工程系陈德智Email dzhchen Tel Office Room108 电机楼2010 12 第4章电磁场有限元法 FiniteElementMethod FEM 有限元法可以基于变分原理导出 也可以基于加权余量法导出 本章以加权余量法作为有限元法的基础 以静电场问题的求解为例介绍有限元法的基本原理与实施步骤 并介绍有限元法中的一些特殊问题 第4章电磁场有限元法 FEM 有限元基本原理与实施步骤 1DFEM有限元基本原理与实施步骤 2DFEM有限元方程组的求解二维有限元工程应用三维有限元原理与工程应用矢量有限元 加权余量法回顾 对算子方程用作为该方程的近似解 试探解 代入方程得余量 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 有限元的系数阵总是对称的吗 与有限差分法 FDM 相比 有限差分法是对点的离散 得到一系列离散点上的解 而有限元 FEM 是对区域的离散 单元 尽管所求的是节点上的自由度 但它的解在场域中每一个点上都有定义 所以 即是有限元节点上的解是精确的 有限元的整个解仍然是近似的 好的数据处理技术可以从该近似解中提取更精确的分析结果 线性单元中 如果所求的自由度是电位j 单元中的电场E是场量 节点上的E取邻近单元的平均 一些补充说明 关于有限元的解 计算系数阵是有限元分析的主要工作量 所涉及到的积分 如果不是解析可积的 通常要用到数值积分 其中最常用的数值积分方法是Gauss数值积分 一些补充说明 高斯数值积分 先将积分区间变换到 1 1 上 按照固定的积分点计算若干函数值P xi 以固定权值wi累加即可 具 2n 1 阶精度 n 4x 1 0 8613d0 x 2 0 339981043584856d0w 1 0 347854845137454d0w 2 0 6526d0n 5x 1 0 9664d0 x 2 0 538469310105683d0 x 3 0 0d0w 1 0 236926885056189d0w 2 0 478628670499366d0w 3 0 568888888888889d0n 6x 1 0 932469514203152d0 x 2 0 6612d0 x 3 0 238619186083197d0w 1 0 1770d0w 2 0 3639d0w 3 0 4679d0 n 16x 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 0950125098d0w 1 0 0271524594d0w 2 0 0622535239d0w 3 0 0951585117d0w 4 0 1246289713d0w 5 0 1495959888d0w 6 0 1691565194d0w 7 0 1826034150d0w 8 0 1894506105d0 一些Gauss积分点和权值 关于x 0对称 只给出一半 为提高有限元分析精度 有两种方法 其一 增加节点 细化网格 称为h方法 其二 增加有限元的阶数 称为p方法 一些补充说明 线性单元与高阶单元 一些补充说明 二阶单元 一些补充说明 三阶单元 h方法和p方法的求解精度 ByJianmingJin TheFiniteElementMethodinElectromagnetics 2ndEd 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上 Ni 0 故积分贡献为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 稳定性 选主元 迭代法 迭代法的基本思想 等价方程组 从一组猜测的初值开始迭代 直至不再变化为止 即为方程组的解 收敛 好的迭代法应该 对初值不敏感 收敛速度快 例如 高斯 赛德尔迭代 有限差分法常用 高斯 赛德尔迭代 实际运算过程 这就是通常的高斯塞德尔迭代格式 矩阵中的零元素不参与运算 矩阵甚至可以不出现 大大减少了内存需求量和计算量 共轭梯度法 ConjugateGradientMethod CG法 共轭梯度法在原理上可以通过n步迭代得到方程的准确解 因而也称为半直接法或半迭代法 把迭代法表示为更一般的形式 称为步长 p称为搜索方向 用r表示残差 预优共轭梯度法 PreconditionedConjugateGradientMethod PCG法 当系数阵的特征值较均匀地分布在一个很长的区间上时 称矩阵具有很大的条件数 此时共轭梯度法的收敛速度可能很慢 解决的办法是选取合适的非奇异矩阵C进行处理 矩阵称为预处理矩阵 预优矩阵M应具有如下特性 稀疏性与A相近 矩阵的特征值分布集中 形如的方程组容易求解 易于寻找 目前公认有效的方法是对系数阵A做不完全Cholesky分解 以M LDLT为预优矩阵 这种方法称为不完全分解预优共轭梯度法 IncompleteCholeskydecompositionpreconditionedConjugateGradientMethod ICCG法 ICCG法 稀疏矩阵技术没有稀疏矩阵技术就没有有限元的成功 带状矩阵技术通过合适的节点编码 可使系统矩阵的非零元素集中于主对角线附近的带形区域内 在使用LU分解法求解方程组的过程中 带形区域以外的0元素无需计算 缺点 节点编码优化 带形区域内仍然有大量零元素 非零元素存储技术只存储矩阵的非零元素 一般用一个一维数组存储矩阵的非零元素 用另一个索引数组存储这些非零元素在原矩阵中的行和列值 例如 非零元素存储技术在迭代法中 系统矩阵参与的运算只有矩阵左乘以某个相量 如果采用上面的存储方法 则实现q Ap的算法如下 CComputethematrix vectorproductDOk 1 noti 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刘万勋 刘长学 华伯浩等 大型稀疏线性方程组的解法 国防工业出版社 1981RP梯华森 稀疏矩阵 科学出版社 1981 有限元分析精度的影响因素静电场问题静磁场问题涡流问题波的传输与散射 4 二维有限元的工程应用 有限元分析精度的影响因素 1 数学模型对工程问题的近似 2 材料电磁参数的不确定性 3 数学模型的有限元近似 场域拟合精度 单元大小 未知数个数与局部场的变化情况 边界拟合精度 曲线边界 系数阵计算过程中数值积分精度 方程求解精度 数字误差 计算结果的数据处理 好的处理技术可以提高分析精度 h方法 p方法 高次单元 单元阶数与计算精度的关系ByJin 曲边三角形单元 更好地拟合曲面边界 数值积分 三角元的高斯数值积分单元比较小的时候 单元内的函数可以近似认为是常数 通常可以获得满意的精度 当单元内函数变化比较快 或者采用曲边三角形单元 高次单元时时 都会用到数值积分 高斯 勒让德数值积分公式 L1 L2 L3是位置的面积坐标 权值wi如下页表格 三角形单元高斯勒让德数值积分点和权值 三角形单元高斯勒让德数值积分点和权值 二维边值问题的通用形式 在媒质交界面G12上 单元矩阵元素计算公式 单元右端项计算公式 二维边值问题一般形式 在媒质交界面G12上的条件自动满足 第一类边界条件需要强加 第三类边界条件的计算参看金建铭 电磁场有限元方法 第三类边界条件中 如果g 0则成为第二类边界条件 如果再有q 0则自动满足 无需任何处理 二维静电场 无限远边界条件 静电场分析 电力线平行与垂直边界条件 注意与下式的区别 轴对称静电场 无限远边界条件 静电场分析 如果不将r看作系数 有系数阵将不再对称 二维静磁场 无限远边界条件 静磁场分析 磁力线垂直与平行边界条件 轴对称静磁场 静磁场分析 在平行平面场中 等A线平行于磁力线 在轴对称静磁场中 等rA线平行于磁力线 注意 与轴对称静电场的表达式非常不同 如果直接以A为求解量 系数阵将不再对称 故 轴对称静磁场问题无限远边界条件 小电流环产生的磁位 一阶渐近边界条件 适合于电流回路 球心位于回路中心的球面 其中 二维静磁场分析 包括轴对称静磁场分析 一般不使用磁标位 因为正如所看到的 磁矢位A只有一个分量 而且激励电流和边界条件都很容易处理 使用磁标位时必须单独处理电流 有时候还需要规定磁障碍面 对有限元来说很不方便 但是三维的情况将非常不同 此为后话 二维涡流问题 涡流场分析 相量的含义 参数是位置的函数 直流与交流的区别 透入深度与单元尺寸 电场与涡流的计算 电压激励线圈的处理 载流块导体的处理 轴对称涡流问题 涡流场分析 在不引起混淆的情况下 可以省略表示相量的点 线圈电压与阻抗的计算 另一种形式 涡流场分析中的无限远边界条件 如果所有的电流 包括源电流和涡流 都位于有限区域内 有限元模型的外围区域为空气 那么涡流问题的无限远边界条件与恒定磁场的无限远边界条件一致 如果导体区域一直延伸到无限远处 有限元模型的外围区域包含导体 普通的无限远边界条件是无效的 需要重新推导 轴对称涡流场中线圈电压与阻抗的计算 1 2 二维波导问题 波的传输与散射 二维波散射问题 波的传输与散射 散射场一阶吸收边界条件 ABC 总场等于入射场加散射场 二维有限元的工程问题 非线性问题电磁力和力矩的计算永磁问题功率与阻抗ANSYS中的二维电磁场分析 非线性问题的求解 设求解静磁场问题 如果磁导率是常数 使用有限元法可将上述问题化为线性代数方程组求解 但如果磁导率是磁感应强度B的函数 从而也是磁矢位A的函数 当A还没获得之前 是未知的 从而矩阵 K 也未知 非线性问题需要用迭代法求解 非线性问题的求解 使用有限元法求解非线性问题的思路是 k 0 猜测一个 计算系数阵 求解有限元方程 得到各点的B k 1 k k 1 估计新的 以上是线性迭代的流程 收敛缓慢 更好的方法是牛顿 拉普逊方法 非线性问题的求解 使用有限元法求解非线性问题的思路是 k 0 猜测一个 计算系数阵 求解有限元方程 得到各点的B k 1 k k 1 估计新的 以上是线性迭代的流程 收敛缓慢 更好的方法是牛顿 拉夫逊方法 求解非线性方程 从一个初值x0出发 求得近似方程的解 依次迭代 逐渐逼近问题的真解x 牛顿 拉夫逊法 Newton Raphsonmethod 左端记为 是一个向量函数 如果 K 不是 u 的函数 则导数矩阵 称为雅克比矩阵 记做 J 方程组的解一次即可获得如果 K 是 u 的函数 则雅克比矩阵的元素为迭代求解方程组 非线性方程组的牛顿 拉夫逊算法 牛顿 拉夫逊法对初值很敏感 有一些改进的算法 但是目前没有从根本上解决问题的路子 使用ANSYS计算铁磁材料磁场时 如果磁化曲线和初始点不合适 也会碰到类似问题 牛顿 拉夫逊法不收敛的例子 B H曲线的表示 为了计算非线性有限元矩阵 必须输入铁材料的磁化曲线 磁化曲线可以通过若干点插值获得 如ANSYS 也可以采用多项式或者其他方式 如指数形式 的拟合表示 磁化曲线通常不可能很准确 计算时如果需要可以根据实际情况作适当的调整 以保障非线性迭代的收敛和计算的精度 正弦涡流问题B H曲线的等效处理 严格说来 含有非线性材料的场不可能是正弦电磁场 当只关心力 力矩 功率 损耗等随时间平均意义的量 而不怎么关心波形的时候 正弦近似有很好的精度 否则只能按照瞬态场计算 计算量很大 正弦近似B H曲线的等效处理方法 令等效磁导率得到的磁能密度与瞬时形式磁能在一个周期内的平均值相等 得到 He为磁场有效值 可用数值积分获得 作业 计算等效磁化曲线 H A m 1 B T 0090 00000000 50000000270 0000001 00000000318 2500001 10000000384 5000001 20000000479 5000001 30000000608 5620001 38750000755 4370001 45000000939 1850001 500000001188 930001 545000001407 930001 575000002077 310001 627500003117 930001 673750003969 370001 70225000 H A m 1 B T 4843 660001 727500006081 340001 758250008581 090001 8087500011066 40001 8500000014985 70001 9025000033003 30