电磁场数值分析与计算05-有限元方法介绍

文章目录
  1. 1. 目录
  2. 2. 引言
  3. 3. 伽辽金有限元法
  4. 4. 子区域展开法
  5. 5. 总结
  6. 6. 商业软件分析过程
  7. 7. 有限元法的步骤
  8. 8. 参考资料

电磁场数值分析与计算课程笔记
无PPT,按照板书整理,若有错误敬请指正。

目录

1.电磁场数值分析与计算01-场论
2.电磁场数值分析与计算02-Maxwell方程组
3.电磁场数值分析与计算03-电磁场数值分析的定解问题
4.电磁场数值分析与计算04-边界条件
5.电磁场数值分析与计算05-有限元方法介绍
6.电磁场数值分析与计算06-2D有限元分析

引言

对于实际工程问题,可以根据相关的数学物理定律建立模型,列写PDE偏微分方程求其解析解(精确解),但复杂问题的解析解的求解十分困难。
有限元法是用于获得数学物理模型近似解(数值解)的数值方法。有限元法于二十世纪四十年代首次提出,用于求解飞机涉设计相关问题。

边值问题指有边界条件的定解问题,其形式为\(\delta \Phi = f\)
其中\(\Phi\)为未知函数,\(f\)为激励,\(\delta\)为微分算子。
对于电磁场实际工程问题,其涉及的方程式为 \(\left\{ {\begin{array}{*{20}{c}} {\bf{\nabla}}^2 \Phi = -\frac{\rho}{\varepsilon} & 电场 \\ {\bf{\nabla}} \times {\bf{\nabla}} \times {\bf{A}} = \mu {\bf{J}}& 磁场 \end{array}} \right.\)

常见的有限元方法包括里兹(Ritz)变分法与伽辽金(Galerkin)有限元法,本课介绍的是伽辽金有限元法。

伽辽金有限元法

伽辽金法于1915年提出,又被称作加权余量法(Weighted residual method)。
选取有限多项式函数(基函数/形函数)将其叠加,再要求结果在求解域内及边界上的加权积分满足原方程,便可以求得一组易于求解的线性代数方程,且自然边界条件能够自动满足。
假设\(\tilde \Phi\)\(\delta \Phi = f\)的一个近似解,即\(\delta \tilde \Phi \approx f\),将其代入\(\delta \Phi = f\)的方程中,必然会产生一个非零的余量(残差,residual error)\(r=\delta \tilde \Phi - f \ne 0\)
对于\(\tilde \Phi\)的最好近似即为在求解区域的每一个点上使得\(r\)为最小值
\(R_i=\int_\Omega{w_i(\delta \tilde \Phi - f)d\Omega}\)
其中\(w_i\)为选定的权重函数(weighting function),权重函数的表达形式与近似解一致
\(R_i\)为权重余量的积分。

【举例】求解两个平行板之间的静态电位\(\Phi\)
上极板\(x=1 \quad \Phi=1V\);下极板\(x=0 \quad \Phi=0V\)
图1 平行板电容器

两个平行板电容器之间填充介质,介电常数为\(\varepsilon F/m\),电荷密度\(\rho(x)=-(x+1)\varepsilon C/m^3\)

该问题为静电场问题,故列出控制方程\({\bf{\nabla}}^2 \Phi = -\frac{\rho}{\varepsilon}\),该方程为静电场电位的泊松方程。
由于本问题只涉及\(x\)轴一个维度,故方程可写为\(\frac{\partial^2 \Phi}{\partial x^2}=-\frac{\rho}{\varepsilon}\)
代入题干条件及边界条件,得到待求解的方程组为 \[\left\{ {\begin{array}{*{20}{c}} \frac{\partial^2 \Phi}{\partial x^2} = x+1 \quad (0<x<1) \\ \Phi |_{x=0}=0 \\ \Phi |_{x=1}=1 \end{array}} \right.\] 【解法1】列方程求其精确解
\(\Phi '' (x)=x+1 \Rightarrow \Phi ' (x)=\frac{1}{2}x^2+x+C_1 \Rightarrow \Phi (x)=\frac{1}{6}x^3+\frac{1}{2}x^2+C_1 x+C_2\)
\(\Phi (0)=C_2=0 \quad \Phi (1)=\frac{2}{3}+C_1+C_2=1 \Rightarrow C_1=\frac{1}{3} \quad C_2=0\)
解得\(\Phi (x)=\frac{1}{6}x^3+\frac{1}{2}x^2+\frac{1}{3}x\)

【解法2】伽辽金有限元法
加权余量公式\(\int_0^1{w_i(\frac{d^2 \Phi}{d x^2} - x-1)dx=0}\)
注:只有\(x\)一个维度,因此可以将\(\frac{\partial^2 \Phi}{\partial x^2}\)写作\(\frac{d^2 \Phi}{d x^2}\)
\(\tilde \Phi\)的展开式为\(C_1+C_2 x+C_3 x^2+C_4 x^3\)
\(\Phi(0)=0 \Rightarrow C_1=0\)
\(\Phi(1)=1 \quad C_1=0 \Rightarrow C_2+C_3+C_4=1 \Rightarrow C_2=1-C_3-C_4\)
\(\Rightarrow \tilde \Phi = (1-C_3-C_4)x+C_3 x^2+C_4 x^3=x+C_3(x^2-x)+C_4(x^3-x)\)
\(\frac{d^2 \Phi}{d x^2} - x-1=2C_3+6C_4 x-x-1\)
假设\(w_1=x^2-x\)\(\int_0^1{(x^2-x)(2C_3+6C_4 x-x-1)dx=0} \Rightarrow \frac{1}{3}C_3+\frac{1}{2}C_4-\frac{1}{4}=0\)
假设\(w_2=x^3-x\)\(\int_0^1{(x^3-x)(2C_3+6C_4 x-x-1)dx=0} \Rightarrow \frac{1}{2}C_3+\frac{4}{3}C_4-\frac{23}{60}=0\)
\(\Rightarrow C_3=\frac{1}{2} \quad C_4=\frac{1}{6} \quad C_2=1-C_3-C_4=\frac{1}{3}\)
\(\Rightarrow \tilde \Phi = \frac{1}{3}x+\frac{1}{2}x^2+\frac{1}{6}x^3\)
本例子的数值解与精确解一样,是因为假设的多项式形式一致,区域简单。

子区域展开法

对于控制方程复杂,区域复杂的实际工程问题,伽辽金法不能很好的应对,最好将整个区域分成若干个子区域逐一分析求解,即子区域展开法。
将整个区域分成\(n\)个子区域,区间端点为\(x_1,x_2,\cdots,x_n,x_{n+1}\),对应端点上的函数值为\(\Phi_1,\Phi_2,\cdots,\Phi_n,\Phi_{n+1}\),其中\(\Phi_1\)\(\Phi_{n+1}\)值已知。
\(n\)个子区域的区间为\(x_i \leqslant x_{i+1},\quad i=1,2,3,\cdots,n\)
假设在每一个区域中内物理量\(\Phi\)呈线性变化,那么对第\(i\)个区域的两个端点\((x_i,\Phi_i)\)\(x_{i+1},\Phi_{i+1}\)进行一阶插值,得到该子区域的基函数为
\(\tilde \Phi (x)=\Phi_i\frac{x_{i+1}-x}{x_{i+1}-x_i}+\Phi_{i+1}\frac{x-x_i}{x_{i+1}-x_i}\)
其中\(\tilde \Phi (x_i)=\Phi_i \quad \tilde \Phi (x_{i+1})=\Phi_{i+1}\),两个端点的函数值已知。
该子区域的基函数图像就是两个端点\((x_i,\Phi_i)\)\((x_{i+1},\Phi_{i+1})\)构成的一条线段。
图2 子区域的基函数
该子区域的权重函数为
\(w_i=\left\{ {\begin{array}{*{20}{c}} \frac{x-x_{i-1}}{x_i-x_{i-1}} & x_{i-1}\leqslant x \leqslant x_i \\ \frac{x_{i+1}-x}{x_{i+1}-x_i} & x_i\leqslant x \leqslant x_{i+1} \end{array}} \right.\)
\(x=x_i\)时,\(w_i=1,\Phi=\Phi_i\);当\(x=x_{i-1}\)\(x=x_{i+1}\)时,\(w_i=0\)
图3 子区域的权重函数
【解法3】子区域展开法
将原有区域等举例拆分为3部分,得到四个端点\(x_1=0,x_2=\frac{1}{3},x_3=\frac{2}{3},x_4=1\),设对应的函数值为\(\Phi_1,\Phi_2,\Phi_3,\Phi_4\),其中已知\(\Phi_0=0,\Phi_4=1\)
图4 子区域的拆分

对加权余量公式进行整理:
\(\int_{x_{i-1}}^{x_{i+1}}{w_i(\frac{d^2 \tilde \Phi}{d x^2} - x-1)dx}=\int_{x_{i-1}}^{x_{i+1}}{w_i\frac{d^2 \tilde \Phi}{d x^2}dx}-\int_{x_{i-1}}^{x_{i+1}}{w_i(x+1)dx}=0\)
其中由分部积分知\(\int_{x_{i-1}}^{x_{i+1}}{w_i\frac{d^2 \tilde \Phi}{d x^2}dx}=\left. w_i\frac{d\tilde \Phi}{dx}\right|_{x=x_{i-1}}^{x=x_{i+1}}-\int_{x_{i-1}}^{x_{i+1}}{\frac{d w_i}{d x}\frac{d \tilde \Phi}{dx}dx}=-\int_{x_{i-1}}^{x_{i+1}}{\frac{d w_i}{d x}\frac{d \tilde \Phi}{dx}dx}\)
\(\Rightarrow \int_{x_{i-1}}^{x_{i+1}}{\frac{d w_i}{d x}\frac{d \tilde \Phi}{dx}dx}+\int_{x_{i-1}}^{x_{i+1}}{w_i(x+1)dx}=0\)
权重函数:\(w_2=\left\{ {\begin{array}{*{20}{c}} \frac{x-x_1}{x_2-x_1} & x_1\leqslant x \leqslant x_2 \\ \frac{x_3-x}{x_3-x_2} & x_2\leqslant x \leqslant x_3 \end{array}} \right. \quad w_3=\left\{ {\begin{array}{*{20}{c}} \frac{x-x_2}{x_3-x_2} & x_2\leqslant x \leqslant x_3 \\ \frac{x_4-x}{x_4-x_3} & x_3\leqslant x \leqslant x_4 \end{array}} \right.\)
\(\Rightarrow \frac{dw_2}{dx}=\left\{ {\begin{array}{*{20}{c}} \frac{1}{x_2-x_1} & x_1\leqslant x \leqslant x_2 \\ -\frac{1}{x_3-x_2} & x_2\leqslant x \leqslant x_3 \end{array}} \right. \quad \frac{dw_3}{dx}=\left\{ {\begin{array}{*{20}{c}} \frac{1}{x_3-x_2} & x_2\leqslant x \leqslant x_3 \\ -\frac{1}{x_4-x_3} & x_3\leqslant x \leqslant x_4 \end{array}} \right.\)
基函数:
\(x_1\leqslant x \leqslant x_2 \quad \tilde \Phi (x)=\Phi_1\frac{x_2-x}{x_2-x_1}+\Phi_2\frac{x-x_1}{x_2-x_1} \quad \Rightarrow \quad \frac{d\tilde \Phi}{dx}=-\frac{\Phi_1}{x_2-x_1}+\frac{\Phi_2}{x_2-x_1}\)
\(x_2\leqslant x \leqslant x_3 \quad \tilde \Phi (x)=\Phi_2\frac{x_3-x}{x_3-x_2}+\Phi_3\frac{x-x_2}{x_3-x_2} \quad \Rightarrow \quad \frac{d\tilde \Phi}{dx}=-\frac{\Phi_2}{x_3-x_2}+\frac{\Phi_3}{x_3-x_2}\)
\(x_3\leqslant x \leqslant x_4 \quad \tilde \Phi (x)=\Phi_3\frac{x_4-x}{x_4-x_3}+\Phi_4\frac{x-x_3}{x_4-x_3} \quad \Rightarrow \quad \frac{d\tilde \Phi}{dx}=-\frac{\Phi_3}{x_4-x_3}+\frac{\Phi_4}{x_4-x_3}\)
\(x=x_2\)时,分析\(w_2\)
\(\int_0^{\frac{1}{3}}{\frac{d w_2}{dx}\frac{d \tilde \Phi}{dx}dx}+\int_{\frac{1}{3}}^{\frac{2}{3}}{\frac{d w_2}{dx}\frac{d \tilde \Phi}{dx}dx}+\int_0^{\frac{1}{3}}{w_2(x+1)dx}+\int_{\frac{1}{3}}^{\frac{2}{3}}{w_2(x+1)dx}=0\)
\(\Rightarrow \int_0^{\frac{1}{3}}{3(-3\Phi_1+3\Phi_2)dx}+\int_{\frac{1}{3}}^{\frac{2}{3}}{-3(-3\Phi_2+3\Phi_3)dx}+\int_0^{\frac{1}{3}}{3x(x+1)dx}+\int_{\frac{1}{3}}^{\frac{2}{3}}{3(\frac{2}{3}-x)(x+1)dx}=0\)
\(\Rightarrow (-3\Phi_1+3\Phi_2)+(3\Phi_2-3\Phi_3)+\left.(x^3+\frac{3}{2}x^2)\right|_0^\frac{1}{3}+\left.(-x^3-\frac{1}{2}x^2+2x)\right|_\frac{1}{3}^\frac{2}{3}=0\)
\(\Rightarrow 6\Phi_2-3\Phi_3+\frac{4}{9}=0\)
\(x=x_3\)时,分析\(w_3\)
\(\int_{\frac{1}{3}}^{\frac{2}{3}}{\frac{d w_3}{dx}\frac{d \tilde \Phi}{dx}dx}+\int_{\frac{2}{3}}^1{\frac{d w_3}{dx}\frac{d \tilde \Phi}{dx}dx}+\int_{\frac{1}{3}}^{\frac{2}{3}}{w_3(x+1)dx}+\int_{\frac{2}{3}}^1{w_3(x+1)dx}=0\)
\(\Rightarrow \int_{\frac{1}{3}}^{\frac{2}{3}}{3(-3\Phi_2+3\Phi_3)dx}+\int_{\frac{2}{3}}^1{-3(-3\Phi_3+3\Phi_4)dx}+\int_{\frac{1}{3}}^{\frac{2}{3}}{(3x-1)(x+1)dx}+\int_{\frac{2}{3}}^1{3(1-x)(x+1)dx}=0\)
\(\Rightarrow (-3\Phi_2+3\Phi_3)+(3\Phi_3-3\Phi_4)+\left.(x^3+x^2-x)\right|_\frac{1}{3}^\frac{2}{3}+\left.(3-x^3)\right|_\frac{1}{3}^\frac{2}{3}=0\)
\(\Rightarrow 3\Phi_2-6\Phi_3+\frac{22}{9}=0\)
由①②\(\Rightarrow \Phi_2=\frac{14}{81} \quad \Phi_3=\frac{40}{81}\)
所得函数\(\tilde \Phi\)的图像为依次过端点\((0,0),(\frac{1}{3},\frac{14}{81}),(\frac{2}{3},\frac{40}{81}),(1,1)\)的线段的集合。
增加子区域的数目和选用高阶权函数可以提高所得函数图像的精确度。

总结

伽辽金法:
①测试函数是基函数\((1,x,x^2,x^3)\)的线性组合。
②可在全域表示,近似真实解,满足边界条件。

子区域展开法:
①测试函数也是基函数的线性组合。
②仅在子区域中表示,再组合成整个区域。
③在有限元中,测试函数是定义在子区域中的一系列基函数的线性组合。如果子区域足够小的话,定义在子区域的基函数可以很简单。

商业软件分析过程

①前处理(preprocess)
模型简化⇒构建模型⇒添加材料⇒施加激励⇒边界条件⇒剖分(meshing)
其中高亮部分是有限元的离散化过程
②求解设置
③后处理

有限元法的步骤

区域离散、插值函数的选择、系统方程的组装与求解
插值方法包括线性插值、二次插值、Cubic插值、Spline插值等。

参考资料