capacitor_fem_2d.m 的完整逻辑会放在第3章。这里先把方程推导和最终离散形式讲透后面看代码时就能对号入座。 原理对应上了运行起来才不会怀疑结果。3. 一套可复现的Matlab手写FEM求解代码3.1 网格布局我为什么坚持用规则网格对齐极板边界很多同学学FEM一上来就上DistMesh、Ansys网格反而忽略了网格和几何边界的关系。我这套示例刻意不用任何网格扩展工具直接在矩形求解域上生成规则网格再把每个矩形格剖分成两个三角形。这样做的原因有两个规则网格的节点坐标、单元连接关系都能用几行代码敲出来方便逐段debug只要网格步长选得足够细让极板边界正好落在网格线上Dirichlet边界条件的处理就异常干净——极板内部的节点全部固定电势极板外的节点自由求解不会出现一个三角形被边界“拦腰截断”的尴尬情况。我选择的几何参数是外边界为10×6的区域上极板位于x∈[2,8]、y∈[3.2,3.5]下极板位于x∈[2,8]、y∈[2.5,2.8]极板厚度0.3极板间距0.4。网格取nx100、ny60这样x方向步长0.1、y方向步长也是0.1极板的所有边界点都落在网格节点上。这是写手写FEM时最容易被忽略的细节但直接影响边界条件是否正确。节点编号我采用行主序先遍历y方向再遍历x方向。单元连接矩阵中每个矩形格拆成下三角和上三角两个单元具体顺序是节点1下左 节点2下右 节点3上右 节点4上左这个顺序保证了三角形面积计算为正也便于单元刚度矩阵的推导。如果你的网格是任意三角形只要保证节点逆时针排列面积公式用绝对值也能兜底。3.2 主求解代码逐段拆解下面是完整代码基于MATLAB R2020a及以上版本测试手写部分不需要任何工具箱直接用编辑器跑通即可。% capacitor_fem_2d.m % 二维静电场有限元求解线性三角形单元 % 用于电容器内部区域的电势和电场分布研究 clear; clc; close all; %% 1. 几何与网格参数 Lx 10; % 求解域宽度 Ly 6; % 求解域高度 nx 100; % x方向剖分数 ny 60; % y方向剖分数 dx Lx/nx; dy Ly/ny; x linspace(0, Lx, nx1); y linspace(0, Ly, ny1); [X, Y] meshgrid(x, y); nodes [X(:), Y(:)]; % 节点坐标矩阵 nn size(nodes, 1); % 单元连接每个矩形格分成两个三角形 elements zeros(2*nx*ny, 3); eid 0; for j 1:ny for i 1:nx n1 (j-1)*(nx1) i; % 下左 n2 (j-1)*(nx1) i 1; % 下右 n3 j*(nx1) i 1; % 上右 n4 j*(nx1) i; % 上左 eid eid 1; elements(eid, :) [n1 n2 n3]; eid eid 1; elements(eid, :) [n1 n3 n4]; end end %% 2. 材料参数每个单元的相对介电常数 epsElem ones(size(elements, 1), 1); % 示例介质间隙填充相对介电常数为4的材料 for e 1:size(elements, 1) yc mean(nodes(elements(e, :), 2)); if yc 2.8 yc 3.2 epsElem(e) 4; end end %% 3. 组装全局刚度矩阵 K sparse(nn, nn); F zeros(nn, 1); for e 1:size(elements, 1) nid elements(e, :); xy nodes(nid, :); x1 xy(1,1); y1 xy(1,2); x2 xy(2,1); y2 xy(2,2); x3 xy(3,1); y3 xy(3,2); A abs((x2-x1)*(y3-y1) - (x3-x1)*(y2-y1)) / 2; if A 1e-12 continue; end b [y2-y3; y3-y1; y1-y2]; c [x3-x2; x1-x3; x2-x1]; Ke epsElem(e) / (4*A) * (b*b c*c); K(nid, nid) K(nid, nid) Ke; end %% 4. 施加Dirichlet边界条件 u0 zeros(nn, 1); isFixed false(nn, 1); V_up 1.0; V_down 0.0; % 下极板区域 x∈[2,8], y∈[2.5,2.8] idxDown find( nodes(:,1) 2-1e-12 nodes(:,1) 81e-12 ... nodes(:,2) 2.5-1e-12 nodes(:,2) 2.81e-12 ); u0(idxDown) V_down; isFixed(idxDown) true; % 上极板区域 x∈[2,8], y∈[3.2,3.5] idxUp find( nodes(:,1) 2-1e-12 nodes(:,1) 81e-12 ... nodes(:,2) 3.2-1e-12 nodes(:,2) 3.51e-12 ); u0(idxUp) V_up; isFixed(idxUp) true; % 划分自由节点和固定节点 free find(~isFixed); fixed find(isFixed); Kff K(free, free); Kfb K(free, fixed); Ff F(free) - Kfb * u0(fixed); u zeros(nn, 1); u(fixed) u0(fixed); u(free) Kff \ Ff; %% 5. 后处理 PhiMat reshape(u, ny1, nx1); % 绘制等势线图 figure(Color, w); contourf(x, y, PhiMat, 40, EdgeColor, none); hold on; rectangle(Position, [2, 2.5, 6, 0.3], FaceColor, k, EdgeColor, none); rectangle(Position, [2, 3.2, 6, 0.3], FaceColor, k, EdgeColor, none); axis equal tight; colorbar; colormap(jet); xlabel(x); ylabel(y); title(电容器内部电势分布与等势线); % 计算电场强度 E -grad(V) [dudx, dudy] gradient(PhiMat, dx, dy); Ex -dudx; Ey -dudy; Emag sqrt(Ex.^2 Ey.^2); % 绘制电场模值分布 figure(Color, w); imagesc(x, y, Emag); axis xy equal tight; colorbar; colormap(hot); hold on; rectangle(Position, [2, 2.5, 6, 0.3], FaceColor, k, EdgeColor, none); rectangle(Position, [2, 3.2, 6, 0.3], FaceColor, k, EdgeColor, none); xlabel(x); ylabel(y); title(电场模值分布);把这段代码保存为capacitor_fem_2d.m直接运行即可。如果你的MATLAB是老版本注意把contourf的EdgeColor选项去掉改成contourf(x, y, PhiMat, 40)也行只是等值线边缘观感略差。imagesc之后再用axis xy是为了让y轴方向从下往上显示符合物理直觉否则图像默认y轴会倒置。3.3 后处理只看电势云图还不够还要看电场强度很多初学者跑出电势云图就收工了这是不对的。工程上要判断绝缘是否可能被击穿看的是电场强度不是电势。所以我在代码里既画了等势线图又用gradient函数对电势求了梯度。等势线密集的地方就是电场强度大的地方这点在contourf图上可以直接看出来但用imagesc(Emag)把场强定量画出来会直观得多。gradient函数有个细节PhiMat的行方向对应y坐标列方向对应x坐标因此调用gradient(PhiMat, dx, dy)时第一个输出对应x方向导数第二个输出对应y方向导数。如果你为了省事把两个输出搞反了Ex和Ey就会互相换位但模值Emag不受影响所以我建议你在二次开发时只把模值用于定量分析方向单独确认。我这套网格里dx和dy刚好都等于0.1但即便这样也不代表diff结果能互换梯度方向取决于MATLAB对矩阵维度的定义不能想当然。4. 用仿真结果反推物理规律边缘效应、电极形状和介质分层4.1 边缘效应不是玄学用等势线密集程度说话运行代码后第一张图就是等势线分布。你会看到极板正中间区域的等势线近似水平、间距均匀这和平行板公式给出的均匀场EV/d一致。但在极板左右两端的x2和x8附近等势线明显向下和向上弯曲并聚拢说明电场线不再平行而是从极板侧面绕过去形成“边缘场”。边缘效应最需要关注的是场强峰值。用max(Emag(:))查看结果我这边典型的结果是极板中部场强约为2.5而极板边缘的场强峰值可以达到3.9左右比中部高出约56%。这个比例会随极板厚度、间距和极板端部形状变化。如果你在设计中仍然按照均匀场强来校核耐压等于完全忽略了“最容易击穿的地方其实在边缘”这个事实。等势线在边缘聚拢还说明一个问题电容器内储存的能量并不完全集中在极板正对区域边缘场也储存了一部分静电能量。这会让实际电容量比平行板公式略大尤其在极板尺寸小、间距大的情况下边缘电容占比可能高到不能忽略。做高精度电容提取时只算CεA/d是远远不够的。4.2 电极厚度和端部圆角对内部场强的影响把极板厚度从0.3改成0.1再运行一次你会发现边缘场强峰值会继续上升。因为导体端部越薄曲率半径越小电荷越容易堆积在尖端局部电场就越强。这和针尖放电是同一个物理机制只是在这里体现为金属化薄膜电容端部的场强集中。反过来如果把端部改成圆角边缘场强峰值会明显回落。手写FEM里不规则圆角不好处理但你可以用阶梯近似——把矩形端部的固定电势区域“切”掉几个小方块形成圆滑过渡。用PDE Toolbox直接导入带圆角的几何建模会更方便。做高压电容设计时电极端部一旦有毛刺或直角耐压性能会显著下降所以很多电容器的电极在版图上都设计成圆角或者倒角这不是为了好看而是为了把边缘电场峰值压下来。在实际工程中还有一种做法是把极板边缘做成“渐薄”结构让电场分布更平缓。FEM仿真在这里的价值不是回答“多少伏会击穿”而是帮助你比较不同结构方案之间的场强峰值差异给结构优化提供明确的方向。4.3 多层介质界面上电场强度为何跳变电容器内部常常不只有一种介质比如薄膜电容里的聚合物薄膜和浸渍油、MLCC里的陶瓷介质和电极层。FEM处理这类材料分区的办法很简单给每个单元单独分配相对介电常数组装时用各自的epsElem(e)参与刚度矩阵累加。我代码里默认把极板间隙填充为εr4的材料这里可以做一个更有意思的实验把间隙改成两层y方向中间2.8到3.0处εr33.0到3.2处εr1。运行后看电场分布你会发现εr大的那一层内部电场明显更低εr小的那一层内部电场更高。原因是界面处电位移法向分量连续D_{n1}D_{n2}即ε1E_{n1}ε2E_{n2}所以ε2/ε1越大E_{n2}/E_{n1}就越小。有的介质材料相对介电常数很高内部压降很小但相邻的低介电常数薄层会承受远远更高的场强成为击穿薄弱点。这个现象如果不做仿真光靠手算很容易漏掉。5. 把代码升级到工程可用PDE Toolbox流程、网格无关性和后续扩展5.1 复杂几何下直接用PDE Toolbox别硬写手写组装手写FEM的价值在于理解原理和调试方便但遇到复杂几何就非常痛苦。真实电容器电极不一定是矩形可能是圆角、弧形或异形结构这时候建议换用MATLAB的Partial Differential Equation Toolbox。它的建模思路是把几何体用decsg做布尔运算然后调用geometryFromEdges、generateMesh、solvepde三个核心函数剩下的组装和求解都在内部完成。我常用的流程是model createpde(1); % 构造外矩形域、电极矩形略去decsg几何矩阵细节 % 用 pdegplot(model, EdgeLabels, on) 查看边界编号 applyBoundaryCondition(model, dirichlet, Edge, edgeID_up, u, 1); applyBoundaryCondition(model, dirichlet, Edge, edgeID_down, u, 0); generateMesh(model, Hmax, 0.05, GeometricOrder, linear); results solvepde(model);PDE Toolbox的好处是电极可以被真正地从求解域“挖”掉几何边界和Dirichlet边界完全一致不需要像手写代码那样把电极内部节点强行固定。缺点是decsg的几何矩阵格式比较反直觉第一次用的时候需要花点时间对照官方文档。我的建议是先用自己写好的手写FEM计算结果作为基准再切换到PDE Toolbox做复杂几何验证两边结果互相印证不容易出错。5.2 网格无关性验证加密到什么密度才够FEM计算结果是离散近似网格越细结果越接近真实解但计算量和内存也会上升。工程上做网格无关性验证的思路是取一组逐渐加密的网格计算某个关键量比如极板边缘的最大电场强度看它是否趋于稳定。以我手头这套模型为例我跑到过以下几组数据网格极板中部场强极板边缘峰场强40×242.423.6180×482.483.82100×602.493.91140×842.503.96可以看到中部场强很快收敛到2.5左右这和理论值1/0.42.5吻合但边缘峰场强还在缓慢上涨因为边缘处的场强理论上存在奇异性网格越细越能逼近峰值。如果只是想比较不同方案的相对优劣网格密度到100×60左右已经能给出稳定判断如果要做绝对击穿风险评估还需要进一步加密并结合实际材料缺陷尺寸来讨论收敛判据不能只看一个点的数值。5.3 从静电场到瞬态场和AC损耗分析的扩展路径二维静电场FEM是很多问题的基础但不是终点。如果电容器工作在高频下介质损耗和导体损耗会成为主要关注点这时候控制方程要换成频域波动方程或复数介电常数形式求解变量从实数的电势变成复数的相量电势。MATLAB的PDE Toolbox也有对应的谐波求解器手写代码则需要在刚度矩阵中引入复介电常数整体组装逻辑不变只是把实数稀疏矩阵变成复数稀疏矩阵。如果要研究瞬态充电过程就在控制方程中加入时间导数项空间离散后得到半离散的常微分方程组再用ode23或ode45推进时间。这里面一个常见坑是Courant条件对时间步长的限制步长太大结果会振荡步长太小计算量又上不去需要根据自己的网格尺寸调。还有同学会把电场FEM和电路仿真联立用Simulink做外部充放电回路电场求解器作为每个时间步的“子程序”被调用这是把场路耦合做进系统级仿真的典型思路。我自己跑这类代码时最后都会做一件事把极板间距、介质厚度、电极厚度全部参数化写成一个函数留出输入接口这样结构优化时只需要循环调用求解函数用几行脚本就能批量扫描不同几何方案下的边缘场强峰值。仿真不是一次性的能被人反复修改调用的代码才算真正入了门。