工程电磁场数值计算4(有限元法1)

工程电磁场数值计算（4）（电磁场有限元法）第4章电磁场有限元法（FiniteElementMethod,FEM）有限元法可以基于变分原理导出，也可以基于加权余量法导出，本章以加权余量法作为有限元法的基础，以静电场问题的求解为例介绍有限元法的基本原理与实施步骤。并介绍有限元法中的一些特殊问题。第4章电磁场有限元法（FEM）1.有限元基本原理与实施步骤：1DFEM2.有限元基本原理与实施步骤：2DFEM3.有限元方程组的求解4.二维有限元工程应用5.三维有限元原理与工程应用6.矢量有限元加权余量法回顾：对算子方程用作为该方程的近似解（试探解）：代入方程得余量：1.有限元法基本原理与实施步骤：一维问题()Luf1niiiuu()RLuf在有限元法中，基函数一般用表示。采用Galerkin方案，取权函数与基函数相同。使与余量正交化：(,)[()]d0iiNRNLuf(1,2,,)in{,1,2,,}iNin设L为线性算子，代入，得11[()]d[()]d0nnijjijjjjNLNfNLNf1niiiuN或(1,2,,)in1()ddnjijijNLNNf记diibNf,()dijijKNLN得代数方程组：Kαb加权余量法回顾（续）(,)[()]d0iiNRNLuf利用有限元法求解一维边值问题：（1）单元剖分如图5个单元，6个节点（2）选取基函数22d()01d(0)(1)0uLuuxxxuu111111(,)(,)iiiiiiiiiiixxxxxxxNxxxxxxx1niiiuuN（3）方程离散（计算系数阵[K]和右端项[b]）基函数Ni只是一阶可导的，不能严格满足微分方程，称为“弱解”。diibNf,()dijijKNLNKαb,2222()dd(+)ddddddijijjijjiijKNLNNNNxNNNNx2222dddddddddddddjijjiijixjixxxjjiixxNNxNNxxNNNNxxxx（3）方程离散第一项在xj处为0，在xi处的值被来自(i-1)单元的贡献抵消，故只剩下第二项。由于基函数Ni局域支撑，显见只有不为0。使用分步积分：(1)ji,1,,1,,iiiiiiKKKdiibNf,()dijijKNLN（3）方程离散故diibNf,()dijijKNLN2,2ddddddddddjjiijijiijxxjiijxxNKNNNxNNxNNxxx(1)ji类似，当j=i时1111,ddddddiiiixxiiiiiixxNNKxNNxxx11ddiixiiixbNfNfx右端项：总体方程11121121222322323334334344454454555655656666KKubKKKubKKKubKKKubKKKubKKub强加边界条件：u1=0,u6=01221222323323334344344454554555656100010uuKKKbuKKKbuKKKbuKKKbu（4）求解方程思考：（1）有限元的解跟有限差分法的解有何根本不同？（2）有限元的系数阵总是对称的吗？ux0000.20.03610.03600.40.06280.06250.60.07100.07080.80.05250.05231.000*u22d()01d(0)(1)0uLuuxxxuu与有限差分法（FDM）相比，有限差分法是对点的离散，得到一系列离散点上的解；而有限元（FEM）是对区域的离散（单元），尽管所求的是节点上的自由度，但它的解在场域中每一个点上都有定义。所以，即是有限元节点上的解是精确的，有限元的整个解仍然是近似的。好的数据处理技术可以从该近似解中提取更精确的分析结果。线性单元中，如果所求的自由度是电位，单元中的电场E是场量；节点上的E取邻近单元的平均。一些补充说明：关于有限元的解计算系数阵是有限元分析的主要工作量。所涉及到的积分，如果不是解析可积的，通常要用到数值积分。其中最常用的数值积分方法是Gauss数值积分。一些补充说明：高斯数值积分diibNf,()dijijKNLN111()d()niiiPxxwPx先将积分区间变换到[-1,1]上；按照固定的积分点计算若干函数值P(xi),以固定权值wi累加即可。具(2n+1)阶精度。n=4x(1)=0.861136311594053d0x(2)=0.339981043584856d0w(1)=0.347854845137454d0w(2)=0.652145154862546d0n=5x(1)=0.906179845938664d0x(2)=0.538469310105683d0x(3)=0.0d0w(1)=0.236926885056189d0w(2)=0.478628670499366d0w(3)=0.568888888888889d0n=6x(1)=0.932469514203152d0x(2)=0.661209386466265d0x(3)=0.238619186083197d0w(1)=0.171324492379170d0w(2)=0.360761573048139d0w(3)=0.467913934572691d0n=16x(1)=0.9894003948d0x(2)=0.9445750231d0x(3)=0.8656312024d0x(4)=0.7554044084d0x(5)=0.6178762444d0x(6)=0.4580167777d0x(7)=0.2816035508d0x(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对称，只给出一半）22d()01d(0)(1)0uLuuxxxuu为提高有限元分析精度，有两种方法：其一：增加节点，细化网格——称为h方法。其二：增加有限元的阶数——称为p方法。一些补充说明：线性单元与高阶单元22d()01d(0)(1)0uLuuxxxuu一些补充说明：二阶单元22d()01d(0)(1)0uLuuxxxuu一些补充说明：三阶单元h方法和p方法的求解精度ByJianmingJin.TheFiniteElementMethodinElectromagnetics,2ndEd.,2002作业：要独立完成，凡雷同者没分！！22d()01d(0)(1)0uLuuxxxuu编写有限元程序，计算一维边值问题。改变剖分单元数目，观察解的精度变化。（建议也同时做一个有限差分法的程序，比较二者的精度差别）以二维静电场泊松方程的求解为例。2.有限元法基本原理与实施步骤：二维问题目标：依据加权余量法，利用分域基，建立离散的代数方程组，即确定系数{Kij}和{bi}。diibNf,()dijijKNLN122222()uuLufxyuguhn场域离散二维问题常使用三角形单元离散，便于处理复杂的场域形状，容易实现。单元：互不重叠，覆盖全部场域；每个单元内介质是单一、均匀的。节点：网格的交点，待求变量的设置点。该步骤需要记录的信息：节点编号、节点坐标节点属性(激励源、是否边界等)单元编号单元节点编号单元介质基函数•有限元采用分片逼近的思想，类似于一维情况下使用折线逼近一条任意曲线。•使用分域基Ni，基函数的个数等于节点的个数；每个基函数Ni的作用区域是与该节点i相关联的所有单元。三角形单元内的基函数设三角形三个顶点处待求函数值分别为u1,u2,u3。如果单元足够小，可以采用线性近似，将单元内任意p点的u(x,y)表示为(,)uxyabxcy312123(,)(,)(,)(,)xyxyxyuxyuuu代入三个顶点的坐标和函数值，可以解出a、b、c。得到312123(,)(,)(,)(,)xyxyxyuxyuuu12312311112xxxyyy1232311112xxxyyy2131311112xxxyyy3121211112xxxyyy单元节点的编号按逆时针方向排列！其中，112233(,)uxyNNN记住我们的任务—寻找基函数对比312123(,)(,)(,)(,)xyxyxyuxyuuu(1,2,3)i可得(,)iixyN基函数Ni常被称为插值函数或者形状函数，具有以下性质：（1）是插值的；（2）（3）在相邻单元的公共边界上，Ni是连续的，从而通过Ni构造的逼近函数也是连续的。1()(,)0()ijjijNxyij在积分中，对于确定的i，j的有效取值为i本身以及与节点i相联的周围节点，积分的有效区域为以i、j为公共节点的所有三角形单元，在这些单元中Ni、Nj才有交叠。()dijijKNLN计算系数阵diibNf,()dijijKNLN这些积分可以分单元进行。例如对右图所示的局部编码，K01、K00以及b0的计算公式为：1234560000()dKNLN160101()dKNLN12345600dbNf计算系数阵diibNf,()dijijKNLN以下把单元e的贡献记为()()()()deeeeijijKNLN()()()deeeeiibNf这样，就有(1)(2)(3)(4)(5)(6)00000000000000KKKKKKK(1)(6)010101KKK(1)(2)(3)(4)(5)(6)0000000bbbbbbb每个或的计算都在具体的单元内单独考虑（称为单元分析）。()eijK()eib单元分析：计算单元内积分对系数阵和右端项元素的贡献。()()()()deeeeijijKNLN系数阵元素：•当L为拉普拉斯算子时，由于Ni在单元内是(x,y)的线性函数，经Laplace算子作用后值为0。•但是，在相邻单元的边界上，Ni是连续但是不光滑的，因此对积分的贡献主要来自边界。•为考虑单元边界的影响，需要借助于格林公式：故，()2()dddddeeejeijijijiNKNNxyNNxyNn2()ddVSVS格林公式：1233223321()()()2xyxyyyxxxy23321()()()2yyxxNij31132()()()2yyxxNij23313213122()()()()()()4yyyyxxxxNN(,)iixyN因：写成一般形式，若一个三角形三个顶点编号为i,j,m（逆时针顺序），则13231323122()()()()()()4y