《电磁场数值方法电子教案:有限元(FEM)》由会员分享,可在线阅读,更多相关《电磁场数值方法电子教案:有限元(FEM)(119页珍藏版)》请在金锄头文库上搜索。
1、有限元(FEM),概述,历史 1943 Courant 最早提出思想 20世纪50年代 用于飞机设计 1960 Clough在著作中首先提出名称 19641965年间数学家冯康独立地开创有限元方法并奠定其数学基础 1965 Winslow首次应用于电气工程问题 1969 Silvester推广应用于时谐电磁场问题 应用范围 广泛地被应用于各种结构工程 成功地用来解决其他工程领域中的问题 热传导、渗流、流体力学、空气动力学、土壤力学、机械零件强度分析、电磁工程问题等等,电磁工程应用及发展 静态场时变场,闭域开域,线性非线性,散射,波导、腔体、传输线 标量有限元发展到矢量有限元 高阶矢量有限元 单
2、一方法发展到混合方法 (快速算法) 频域求解发展到时域求解 (区域分解技术) 商用软件:比如HFSS、ANSYS,有限元思想1 有限元法是函数逼近理论、偏微分方程、变分与泛函分析的巧妙结合。从数学上分析,有限元法是Rayleigh-Ritz-Galerkin法的推广。 传统的有限元以变分原理为基础 变分问题就是求泛函极值的问题 直接解法把变分问题化为普通多元函数求极值的问题Ritz 寻找一组在全域上解析、而又要在边界上满足强加边界条件的基函数 间接解法 变分原理 变分问题与对应的边值问题等价,有限元思想2 有限元法采取了与变分问题间接解法相反的途径,把所要求的微分方程型数学模型边值问题,首先转
3、化为相应的变分问题,即泛函求极值问题;然后利用剖分插值,离散化变分问题为普通多元函数的极值问题,即最终归结为一组多元的代数方程组,解之即得待求边值问题的数值解。,有限元思想3 有限元法的核心在于剖分插值,它是将所研究的连续场分割为有限个单元,用比较简单的插值函数来表示每个单元的解,但是它并不要求每个单元的试探解都满足边界条件,而是在全部单元总体合成后再引入边界条件。这样,就有可能对内部和边界上的单元采用同样的插值函数,使方法构造极大地得到简化。,有限元思想4 由于变分原理的应用,使第二、三类及不同媒质分界面上的边界条件作为自然边界条件在总体合成时将隐含地得到满足,也就是说,自然边界条件将被包含
4、在泛函达到极值的要求之中,不必单独列出,而唯一考虑的仅是强制边界条件(第一类边界条件)的处理,这就进一步简化了方法的构造。,有限元法主要特点1 离散化过程保持了明显的物理意义。因为变分原理描述了支配物理现象的物理学中的最小作用原理(如力学中的最小势能原理、静电学中的汤姆逊定理等)。因此,基于问题固有的物理特性而予以离散化处理,列出计算公式,当可保证方法的正确性、数值解的存在与稳定性等前提要素。,有限元法主要特点2 优异的解题能力。与其他数值方法相比较,有限元法在适应场域边界几何形状以及媒质物理性质变异情况的复杂问题求解上,有突出的优点:不受几何形状和媒质分布的复杂程度限制;不同媒质分界面上的边
5、界条件是自动满足的;不必单独处理第二、三类边界条件;离散点配置比较随意,通过控制有限单元剖分密度和单元插值函数的选取,可以充分保证所需的数值计算精度。,有限元法主要特点3 可方便地编写通用计算程序,使之构成模块化的子程序集合。容易并行。 从数学理论意义上讲,有限元作为应用数学的一个分支,它使微分方程的解法与理论面目一新,推动了泛函分析与计算方法的发展。,&3.1变分原理与尤拉方程,在微积分学形成的初期,以数学物理问题为背景,与多元函数的极值问题相对应,就已经在几何、力学上提出了若干求解泛函极值的问题。 例如最速降线问题,即在于研究当质点从定点A自由下滑到定点B时,为使滑行时间最短,试求指点应延
6、着怎样形状的光滑轨道下滑。,沿曲线滑行弧线所需时间为,滑行总时间为,泛函的极值(max或min)问题就称为变分问题。,间接解法是将变分问题转化为尤拉方程 (微分方程)的定解问题,即边值问题来求解。,极值,简写为,只差一个数值因子,等同于,泛函的极值问题就称为变分问题 变分问题与边值问题等价 有限元正是间接求解变分问题过程的逆过程 泛函取极值的过程中 第二、第三类边界条件为自然边界条件 无条件变分问题 第一类边界条件为强加边界条件 条件变分,&3.2与线性问题等价的变分问题,与齐次边值问题等价的变分问题,与泊松方程齐次第三类边值问题等价的变分问题 与泊松方程齐次第二类边值问题等价的变分问题,与泊
7、松方程齐次第一类边值问题等价的变分问题 混合型边界条件,与非齐次边值问题等价的变分问题 与泊松方程非齐次第三类边值问题等价的变分问题,与泊松方程非齐次第二类边值问题等价的变分问题 与泊松方程非齐次混合型边值问题,二维问题 体积分变成面积分,面积分变成线积分 轴对称场 分层介质中的变分问题 变为问题中,由于介质分界面上的边界条件为齐次自然边界条件,所以泛函取极值时自动满足,不必另行处理。,半线系统 无遗漏、无多余的覆盖,&3.3基于变分原理的差分方程,梯形积分公式,线性Lagrange插值公式,梯形公式,当(i,j)为内点,右边点而非角点,上边点而非角点(i,N),角点(M,N),&3.4有限元
8、法求解,给出与待求边值问题相应的泛函及其等价变分问题 应用有限单元剖分场区域,并选取相应的插值函数 将变分问题离散化为一个多元函数的极值问题,导出一组联立的代数方程 单元分析 /总体合成 /强加边界条件处理 选择适当的代数解法,解有限元方程,即得待求边值问题的近似解(数值解) 检验(附加计算),&3.4.1场域剖分,区域离散的方式将影响计算机内存需求、计算时间和数值结果的精确度 一维区域 短直线段 二维区域 小三角形或矩形 三维区域 四面体、三棱柱或矩形块,划分单元的基本原则 总原则 在满足精度要求的前提下,尽量减少剖分单元数目以节省存储量和计算时间 其他 需要详尽了解的部位要切分得细小,其他
9、部位可以粗糙一些,几何形状变化剧烈的地方电磁场也变化大,单元要细小一些 所有单元应接近等边三角形 当有曲线边界和复杂形状边界时就需要把复杂形状用标准形状来逼近 节点、分割线或分割面应设置在几何形状和介质形状发生突变处,注意事项 各单元只能在顶点处相交 不同单元在边界处相连,既不能相互分离又不能相互重叠,单元、节点编号 单元编号 、全局节点编号 、局部节点编号 单元编号 用一组整数给单元编号 没有特殊要求,只要方便 一般按内部单元、第一类、第二类边界条件单元的顺序进行 节点编号 完整描述应包括它的坐标值、局部编号和全局编码 局部编号表示它在单元中的位置:逆时针方向 用两组整数编号,引入3M的整型
10、数组,&3.4.2分片插值与形状函数,插值 用一个简单函数去近似代替真实函数,二者在某些积分点上具有相同的函数值甚至直到某阶导数值 插值函数:一阶(线性)、二阶(二次)、或高阶多项式,: 待定系数,: P点的插值函数 (也称为展开函数或基函数), 通常代表了有限单元上 用来逼近带求场分布的近似规律,它们只有在单元e内才不为零,而在单元e外均为零,形状函数 基函数与单元的形状尺寸有关,线性三角形单元 单元剖分得足够小,以致可将其上的场量看作不变,三角形单元节点数与插值函数的待定系数的个数相等。,的几何意义,形状函数性质,保证了单元两侧解的连续性,&3.4.3 有限元方程建立,齐次自然边界条件下拉
11、氏方程的有限元方程,静电场的边值问题,二维直角坐标系,用三角形剖分,得到E个单元,个节点,每个单元上的能量泛函为,单元分析,单元上的泛函表达式,总体合成,迭加:交汇于同一节点的单元泛函方程,例:如何由单元矩阵组合成总体矩阵,所对应的等价变分问题为,零元素由不在一个单元、不相干的节点产生,对称,主对角线元素占优,正定,稀疏,带状,减小系数矩阵的最大半带宽,先选与其他节点联系最少的节点作为起始节点,然后将相邻节点编为紧接着的号数,使相邻节点的编号数均为相差不多的数,节点全局编号的顺序,D=2,B=3,D=3,B=4,一般沿场域的窄边(节点数少的)编号,非齐次自然边界条件下拉氏方程的有限元方程,二维
12、直角坐标系,剖分时应尽量使场域边界分段线性化,应先编不含非齐次边界条件的单元,再编含非齐次边界条件的单元,应先编不含非齐次边界条件的节点1,2,N, 再编含非齐次边界条件节点N+1,N+2,。,将正对着边界的那个三角形单元顶点编为i,jm的长度为,写为矩阵形式,的值只和与它相联系的节点有关,齐次自然边界条件下泊松方程的有限元方程,静电场的边值问题,二维直角坐标系,三角元剖分,提出积分号外,近似为常量,积分 x-y平面,平面,由二重积分的变换式雅克比式,得到面积元素变换式为,为N阶列阵,非齐次自然边界条件下泊松方程的有限元方程,二维直角坐标系,非齐次边界条件的边界单元才存在,&3.4.3有限元方
13、程的计算,迭代法,系数矩阵正定,非齐次边界条件下的拉氏方程的有限元方程式,任一节点i的第n+1次场值的GS迭代公式,非齐次自然边界条件的节点编号为,其他节点编为,SOR迭代格式,强加边界条件的处理,设第m个节点上具有强加边界条件,处理强加边界条件小结 将对角线元素 置1 该行的右端项改为强加电位值 第m行与m列的其他元素全部置零 除第m行外,其他各行右端项为原右端项减去强加电位置 与对应第m列未变换前的系数的乘积,若强加边界条件节点共有N0个,则如法处理N0次,系数矩阵的存储 利用其对称性 、稀疏性,等带宽存储,变带宽存储,具体作法,将K矩阵的下三角部分带内元素按行的顺序, 从各行的第一个非零
14、元素起, 至主对角元素(包括带内的零元素及对角元素)为止的元素 (共有MD+1),依次存入一维数组KL(n)中,说明KL(n)数组中各元素在原矩阵K中的行、列位置,需要计算的量 K的下三角阵各行的半带宽MD 第i行的半带宽MD=i-min(j) 主对角元素在KL(n)的地址JO JO=L(i)=L(i-1)+MD+1 非主对角元素在KL(n)的地址JO JO=L(i)-(i-j) 非主对角元素的地址为该行主对角元素的地址向前推移(i-j)个位置,K阵中各元素在一维数组KL(n)中的位置为,高斯消元法求解线性代数方程组 消元 将原方程的系数矩阵化为对角元素为1的上三角阵 回代 从最末一个方程(只
15、包含一个未知数)求出一个未知数,再代入上一个方程,求出另一个未知数,消元公式,第一个方程的系数 及右端项的处理,回代公式,电场强度的计算,一单元中E 的数值解,单元的节点一般与相邻几个单元相联系,节点场强必须将有关单元对场强综合贡献,场强的方向,剖分细,对称,分界面处细剖分,高精度:,例1:计算方同轴线间的电位分布,边值问题,等价变分问题,有限元方程,&3.5 应用举例,每个单元的三个顶点编号,每个顶点的xy坐标,内边界节点编号,边界,外边界节点编号,节点总数、单元总数,1、将几何信息读入代码,2、计算总体系数矩阵K,计算每个三角形单元的矩阵系数(子程序),计算三角形单元总面积,计算中间矩阵Se,单个三角形单元的矩阵系数合成总体系数矩阵,3、处理强加边界条件,找到内、外边界节点,为其电位赋值,节点编号,重新排序,4、求解有限元方程,例2:波导场的有限元解,用有限元法求解的优越 将截面逐步细分,将导致本征值向极限单调减小 可保证有一个更快的收敛速度趋于本征值 对难处理的边界形状更容易处理,而不导致非对称矩阵 波导中的奇点无需特别处理,假设 波导壁为完纯导体 波导内的介质系均匀、线性且各向同性的理想介质 波导中无自由电荷和传导电流 波导工