简介基于MATLAB的电磁场数值分析.pdf是一份面向电气工程、电子工程及相关专业学生的技术文献聚焦传统解析法在复杂边界条件下难以求解的问题系统介绍有限差分法FDM与有限单元法FEM在电磁场数值计算中的实现思路。资源以接地矩形导体槽电位分布为例展示如何利用MATLAB将拉普拉斯方程离散化、迭代求解网格节点电位并通过三维曲面图、等位线和电场线可视化场域分布同时对比解析解与数值解误差帮助读者快速掌握二维静态电磁场边值问题的完整数值求解流程。资源为单份PDF文档大小约181KB内容精炼适合作为电磁场理论课程的参考材料或MATLAB数值分析入门训练。目前已有184人学习使用对需要开展电磁场仿真实验或理解数值方法工程应用的读者具有实用价值。1. 从解析解到数值解MATLAB凭什么接管电磁场边值问题一块三面接地的矩形导体槽顶盖电位被固定在100V槽内电位分布是什么用分离变量法能写出解析解但那是一串带sinh的无穷级数看式子就知道实际边界稍微歪一点就解不出来。工程上需要的不是公式本身而是能直接算到每个点、画成曲面和等位线的结果。这正是数值方法的用武之地。有限差分法FDM把拉普拉斯方程在网格上离散成线性方程组有限单元法FEM从能量极值的角度逼近同一个解而MATLAB把这两条路都做成了十几行代码或一次GUI点击。这篇博客用《基于MATLAB的电磁场数值分析》中的接地矩形槽算例从差分格式、PDE Toolbox操作到误差分析和脚本化固化完整过一遍。2. 有限差分法求解拉普拉斯方程的原理与MATLAB实现2.1 二维静态场的拉普拉斯方程与矩形槽模型矩形导体槽内部无自由电荷电位φ满足拉普拉斯方程。在直角坐标系下二维静态场的边值问题写作∂²φ/∂x² ∂²φ/∂y² 0边界条件是第一类Dirichlet条件三条边接地φ0顶部盖板电位φ100V。这类问题理论上有唯一解但解析解只在规则区域好写。矩形槽是少数能用分离变量法求解的形状之一所以非常适合用来验证数值法的精度。需要注意原文在尺寸与网格参数上没有保持严格一致复现时我采用归一化网格坐标x方向41个网孔、y方向21个网孔节点坐标取整数物理尺寸通过等比例映射换算。这样做的好处是差分方程里步长h1代码更干净也方便后面做网格加密试验。2.2 五点差分格式与正方形网格离散2.2.1 从偏微分方程到代数方程用中心差分近似二阶导数对拉普拉斯方程进行离散。对节点(i,j)x方向二阶导约等于(φ(i,j1)-2φ(i,j)φ(i,j-1))/h²y方向同理。代入拉普拉斯方程并令h1整理后得到φ(i,j) (φ(i1,j) φ(i-1,j) φ(i,j1) φ(i,j-1)) / 4这个式子说明每个内部节点的电位是上下左右四个相邻节点电位的算术平均。迭代法的思路就从这来给定边界后用这个公式反复更新所有内点直到相邻两次迭代的差足够小。网格越细这个平均过程越逼近连续场的真实光滑解但节点数也相应变大。2.2.2 网格参数设计原文采用正方形网格网格参数如下表。这些数字直接决定了方程组的规模和迭代收敛时间。参数设置值x方向网孔数41y方向网孔数21节点总数42×22 924内部待求节点40×20 800边界已知节点124迭代精度阈值1e-5最大迭代次数5000内部节点数约占总节点的87%这意味着求解过程几乎全部集中在迭代内点上。边界节点只起到约束作用不参与更新。2.3 高斯-赛德尔迭代的MATLAB代码与收敛控制2.3.1 迭代主循环常见做法是用高斯-赛德尔迭代而非雅可比迭代前者在更新过程中就地使用新值收敛速度更快。核心代码如下% 接地矩形导体槽电位 - 有限差分法 Mx 41; My 21; % x/y方向网孔数 phi zeros(My1, Mx1); % 节点电位矩阵行对y、列对x phi(end, :) 100; % 顶部边界电位 100V tol 1e-5; % 迭代精度阈值 maxIter 5000; % 防止死循环的上限 for iter 1:maxIter phiOld phi; % 保存本轮迭代前的值 % 内点四点平均更新 for i 2:My for j 2:Mx phi(i,j) 0.25 * (phi(i1,j) phi(i-1,j) ... phi(i,j1) phi(i,j-1)); end end % 强制恢复第一类边界条件 phi(1,:) 0; phi(end,:) 100; phi(:,1) 0; phi(:,end) 0; % 收敛判据全矩阵最大变化量 if max(abs(phi(:) - phiOld(:))) tol fprintf(迭代次数: %d\n, iter); break; end end % 等位线与三维曲面绘图 figure; contourf(phi, 20); colorbar; % 转置使x为横轴 figure; surf(phi); shading interp; % 三维电位曲面逻辑说明循环里的两条边界恢复语句不能省略。数值迭代会产生浮点误差如果不每轮强制重置顶部的100V会慢慢向四周扩散最终把边界条件“吃掉”。收敛判据用max取全矩阵最大绝对变化量比用平均绝对误差更严格能保证每个节点都进入稳态。参数说明tol1e-5对应原文的迭代精度要求。实际工程里如果想快速看趋势可以放宽到1e-3但做误差对比至少要到1e-6。maxIter设5000是给高密度网格留余量41×21网格在常见配置下100多次迭代就能收敛后面81×41网格大概需要几百次。2.3.2 与解析解的对比验证解析解由分离变量法得到对a1、b2的矩形槽电位函数为φ(x,y) (4×100/π) × Σ_{n1,3,5,…} sin(nπx) sinh(nπy) / (n sinh(2nπ))实际计算时取前50项即可。下面抽取原文表2中的一部分点做对比这些点的解析解、数值解和相对误差如下表网格点(x,y)解析解/V数值解/V相对误差/%(0.4,0.1)5.48315.48640.0608(0.8,0.1)8.08078.07690.0475(0.2,0.3)9.66179.67890.1784(0.6,0.3)22.237722.23280.0221(0.7,0.5)41.911841.89090.0499(1.0,0.5)18.641618.68290.2219(0.5,0.7)57.504957.46740.0651(0.9,0.7)65.198365.17320.0386从误差列看最大相对误差约为0.22%最小不足0.02%。这个量级对绝大多数工程计算都够用。注意误差并没有在边界附近单调放大而是与局部电位梯度相关梯度大的地方差分近似的截断误差更明显。3. 有限单元法在PDE Toolbox中的操作路径与边界配置3.1 有限元法与有限差分法的本质差异有限差分法直接在网格上近似导数对网格的规则性有要求正方形网格最省事但遇到倾斜边界或圆形导体时阶梯状网格会带来额外误差。有限单元法走的是另一条路先把求解区域剖分成为互不重叠的三角形单元在每个单元上构造插值函数再把拉普拉斯方程的等效变分形式——也就是能量泛函取极值——转化为代数方程组。它不需要规则的网格排列三角形单元能贴合任意复杂边界。MATLAB的PDE Toolbox自动划分三角形网格用户只需要画区域、定边界、设参数。3.2 用pdetool搭建接地矩形导体槽模型3.2.1 打开GUI并设置网格在MATLAB命令窗口输入以下命令启动PDE图形用户界面% 打开PDE Toolbox图形界面 pdetool;界面打开后先通过Options菜单下的Grid和Grid Spacing设置显示网格X-axis linear Spacings填写[-1.5:0.2:1.5]Y-axis取Auto。这个网格只是绘图辅助线不是最终的有限元剖分它帮助你在Draw模式下对矩形的位置有个参考。3.2.2 画区域与选择应用模式单击Draw下的Rectangle/Square按钮在坐标区内画一个矩形覆盖目标槽区域。默认情况下矩形四边就是后续的边界。然后在工具条中间的Generic Scalar下拉框中选择Electrostatics静电学应用模式。这一步很关键它告诉工具箱我们要求解的方程是-div(ε∇φ)ρ后面所有参数都围绕这个方程配置。3.2.3 输入Dirichlet边界条件进入Boundary Mode工具箱会高亮显示各条边界。对每个边界段选择Dirichlet条件并填写h和r。矩形槽四条边界的具体参数如下表边界位置条件类型hr左边界Dirichlet10右边界Dirichlet10下边界Dirichlet10上边界盖板Dirichlet1100h和r对应边界方程h*u ru是电位。h1表示系数归一化r直接就是边界电位数值。三条接地边r0盖板边r100。如果今后要做非零边界只需改r不需要动h。3.3 方程参数、网格剖分与求解点击PDE按钮打开PDE Specification对话框把介电常数epsilon设为1体电荷密度rho设为0。这里epsilon对应方程系数crho对应右侧源项f都为1和0时就是一个标准的拉普拉斯方程。随后选择Initialize Mesh生成初始三角形网格再单击加密按钮带三角形的图标细化一次。网格越密有限元解的精度越高但求解自由度也成倍增加对矩形槽问题加密一次足矣。求解直接单击工具栏上的等号按钮或者选择Solve菜单下的Solve PDE。求解完成后Plot菜单下的Parameter可以设置显示内容勾选Color、Height3D plot、Plot in x-y grid和Show mesh得到带网格的三维电位曲面勾选Color、Contour和Arrows得到二维等位线叠电场方向箭头图。等位线条数和Colormap都可以在对话框中调节。这就是原文图3对应的输出。4. 误差从哪里来网格密度、迭代次数与源项的处理4.1 相对误差的分布规律第2章表里的相对误差有一个特点最大误差出现在(1.0,0.5)这种靠近角落、电位梯度变化较快的位置而不是电位值最大的顶盖附近。原因是五点差分格式的截断误差正比于二阶导数的差分余项电位沿y方向按sinh增长曲率大的地方误差自然大。另一个容易被忽略的点解析解本身是无穷级数截断得到的如果只取前10项解析解自身就带有误差拿它当“真值”比较会误导。我一般取前50项并检查相邻截断项的贡献小于1e-6才停。下面的代码用来计算数值解与解析解的逐点相对误差% 假设 phiNum 来自差分法phiAna 来自解析级数 err abs(phiAna - phiNum) ./ phiAna * 100; fprintf(最大相对误差: %.4f%%\n, max(err(:))); fprintf(平均相对误差: %.4f%%\n, mean(err(:)));这段代码逻辑简单但隐藏着一个陷阱如果phiAna中有零电位点例如接地的边角直接除会得到无穷大。实际处理时通常只统计内部节点或者给分母加上一个极小量eps。对比实验要统一坐标网格否则差分的整数网格和解析解的物理坐标对不上。4.2 网格加密对收敛结果的影响有限差分法的误差按O(h²)下降h是网格步长。把网孔数从41×21加倍到81×41步长减半理论上相对误差应该缩小到原来的四分之一左右。实际因为边界和奇点的影响可能达不到理论速度但趋势不变。迭代次数也会显著增加因为相邻节点间信息传递的距离变长了。参考数据如下网格规模内部节点数相对误差上限参考迭代次数参考41×218000.22%10181×4132000.06%约400161×81128000.015%约1600这里有个工程上常见的折中不需要盲目加密。如果只关心电位的大致分布41×21就够如果要做电场强度电位梯度的仿真由于梯度是数值微分误差会被放大建议至少81×41并用有限元做交叉验证。4.3 常见误用与排错我见过不少人在这个算例上翻车归纳起来三个原因。第一个是把边界节点也放进迭代更新导致四角电位从0和100之间漂移。解决办法是每轮迭代后强制把边界行和列写回指定值顺序不能反先更新内点再恢复边界最后算收敛差。第二个是收敛阈值设得太松比如用norm(phi-phiOld)0.01看起来很快但还没真正收敛。建议用inf范数或逐点最大绝对差。第三个是搞混泊松方程与拉普拉斯方程。如果槽内有体电荷密度ρ内点公式必须改为φ(i,j) (φ(i1,j) φ(i-1,j) φ(i,j1) φ(i,j-1) ρ·h²/ε) / 4注意ρ的符号由方程形式决定MATLAB的PDE Toolbox里静电场方程是-div(ε∇φ)ρ所以这里加号是负号挪过去的。如果ρ为0公式自动退回拉普拉斯版本。5. 把pdetool操作固化成脚本命令行式有限元求解5.1 用新版PDE Toolbox编程接口重写GUI流程新版PDE Toolbox提供了更完整的命令行接口可以替代纯GUI操作。这样每次修改边界电位或材料参数不用重新点二十多次鼠标。接地矩形槽的完整脚本如下model createpde(1); % 1个因变量电位 R [3,4,0,1,1,0,0,0,2,2]; % 矩形 [x0,x1,x2,x3, y0,y1,y2,y3] gd R; ns char(R1); sf R1; [dl, bt, ml, ge] decsg(gd, sf, ns); geometryFromEdges(model, dl); % 四条边编号1-下2-右3-上4-左以实际输出为准 applyBoundaryCondition(model, dirichlet, Edge, 1, u, 0); applyBoundaryCondition(model, dirichlet, Edge, 2, u, 0); applyBoundaryCondition(model, dirichlet, Edge, 3, u, 100); applyBoundaryCondition(model, dirichlet, Edge, 4, u, 0); % 静电场方程-div(epsilon*grad(u)) rho specifyCoefficients(model, m, 0, d, 0, c, 1, a, 0, f, 0); generateMesh(model, Hmax, 0.05); results solvepde(model); u results.NodalSolution; pdeplot(model, XYData, u, ZData, u, ColorBar, on);这段脚本把第3章GUI里每一步都对应起来了decsg负责区域描述applyBoundaryCondition对应Boundary Mode里的四条边设置specifyCoefficients里的c和f对应epsilon1、rho0generateMesh的Hmax控制网格密度。重点说一下Edge编号decsg生成的边界顺序不是按你画图的习惯来的运行后先用一行pdeplot(model,EdgeLabels,on)查看编号再写边界条件否则很容易把100V加到侧面。5.2 一个加速参数扫描的技巧把上面脚本包成一个函数输入是盖板电位U0和网格参数Hmax输出是电位矩阵。用一个循环扫描多组U0就能快速画出不同激励下的电位变化趋势。例如for U0 [50 100 150] u solveChannel(U0, 0.05); % 自定义函数 pdeplot(model, XYData, u); drawnow; end这个技巧在做工程参数分析时非常实用。GUI适合单次查看脚本适合批量回归。差分法和有限元法在同一个算例上结果几乎重叠说明两条路线都没走偏。剩下的就是根据你的实际边界形状和网格需求二选一或者两者互证。本文还有配套的精品资源点击获取