★ FEATURED ARTICLE 前一阵调试一台自制高压试验装置耐压测试做到一半能听到腔体里有细微的嘶嘶声。把设备拆开放电的痕迹不在极板正对区域而是沿着上下极板的外侧边缘一条发亮的细线。这个现象在设想里其实出现过平行板电容器内部的电场教科书上写 EU/d可那是在无限宽极板前提下。真实电容器的极板有边界边缘处电力线会向外弯曲出现明显的电场集中。要定量知道这个集中有多严重、分布在什么范围必须对电容器内部区域做场仿真。我用有限元方法FEM把电位分布离散求解然后用 Matlab 复现了整个计算过程这里把思路、公式、完整代码和工程结论一起整理出来。这篇文章适合两类人一类是做高压绝缘、电力电子、传感器设计需要评估边缘场强和杂散电容的工程师另一类是刚开始接触有限元想找一个不太花哨、能自己动手改参数跑通的案例的学生。文中代码不依赖 PDE 工具箱核心逻辑就是“网格生成、组装刚度矩阵、施加边界条件、求解、后处理”改参数就能复现也方便扩展成你自己的仿真工具。1. 为什么电容器内部区域值得用 FEM 仿真1.1 均匀场假设什么时候靠不住平行板电容器的经典公式 C εS/d 和内部场强 E U/d都建立在“极板无限宽、边缘效应忽略”的前提下。工程里极板当然有边界宽度 W 有限间距 d 一般比 W 小一个数量级。当 W/d 比较大比如 10 倍以上极板正对区域的电场确实接近均匀误差不大但一旦间距拉开、极板变窄或者极板边缘直接暴露在空气中问题就来了。边缘效应最直观的表现是电力线弯曲。越靠近极板边缘电力线越倾向于绕过极板外侧而不是全部垂直穿过极板之间的区域。这种弯曲导致两个后果第一边缘区域的场强会显著高于内部均匀区域成为一个“场强热点”或者叫电场集中点第二边缘场会向外延伸在极板外侧形成杂散电场这部分电场携带能量等效成额外的边缘电容。在高压场景下边缘电场集中的实际意义非常大。气体或者固体绝缘材料的击穿往往从场强最高的位置先开始所以极板边缘经常成为局部放电的“起点”。我前面提到的那台装置就是典型的边缘放电——内部介质还好好的边缘外侧的空气先扛不住了。如果不做仿真光靠均匀场公式根本解释不了为什么击穿点在那个位置。1.2 FEM 到底能算哪些东西有限元方法解决的核心问题是在给定几何、材料参数和边界条件后求一个满足场方程的连续函数 V(x,y)。一旦拿到电位分布电场强度 E -∇V 就能算出来进一步可以算储能、等效电容、界面处的场强跳变等。具体到电容器内部区域FEM 能回答这几类问题场强最大值出现在哪里相对均匀场的增强系数是多少边缘电场的空间影响范围有多大在什么距离上衰减到可忽略填充介质比如内部介电常数 εr 2.2 的固体绝缘与外部空气的界面上电场法向分量如何跳变总电容比平板公式大多少这部分“额外电容”随极板宽度、间距怎么变化如果极板边缘改成圆角边缘场强能降下来多少。这些问题没有解析解可用尤其是在二维复杂形状和多层介质的情况下FEM 是最实用的路径。而且 FEM 的好处是“一次建模、反复使用”参数扫描只是改几个数字重新跑一遍的事。2. 从泊松方程到有限元公式离散化的数学骨架2.1 控制方程与两类边界条件静电场在没有自由电荷的区域满足高斯定律对线性各向同性介质可以写成∇·(ε∇V) 0这里 ε 是介电常数V 是电位。求解域内部可以分成两块极板正对区域填充介质相对介电常数设成 2.2外部区域是空气相对介电常数近似 1。方程本身不复杂复杂的是边界条件。边界条件在本文模型里分两类第一类是狄利克雷边界Dirichlet两个金属极板是等势体上极板固定 V V0下极板固定 V 0。这类条件直接给定节点电位值离散后要通过代数手段“强加”到方程组里。第二类是诺伊曼边界Neumann求解域最外圈是人工截断边界为了模拟“无限大开放空间”外边界上电位的法向导数设为零也就是 ∂V/∂n 0。物理含义是外边界上没有外加的法向电场通量只要截断边界离极板足够远这个近似就足够准。2.2 弱形式为什么“不写”边界条件也能成立有限元求解不从微分方程直接出发而是先把方程乘上一个测试函数 v再在整个求解域上做积分然后利用分部积分把二阶导数降成一阶导数。这一步得到的不是原来的强形式而是弱形式∫Ω ε ∇V · ∇v dΩ 0注意边界项。分部积分之后会出现一个沿边界的积分项 ∫Γ v ε ∂V/∂n dΓ。对于诺伊曼边界∂V/∂n 0这一项自动消失对于狄利克雷边界测试函数 v 在固定边界上取零这一项也自动消失。所以最终弱形式里不需要显式添加任何外边界条件。这是 FEM 对我个人来说最优雅的地方。初学者常常担心“边界条件忘加了怎么办”其实自然边界条件诺伊曼被弱形式自动吸收了你只需要显式处理人为给定的极板电位。这也是为什么本文代码里外部边界什么都不写但结果依然是合理的。2.3 三角形 P1 单元与单元刚度矩阵求解域被剖分成许多三角形单元每个单元内部用电位的线性插值近似。每个三角形有三个顶点每个顶点对应一个基函数基函数在该顶点取 1在另外两个顶点取 0单元内部线性变化。因为线性插值的梯度是常数单元刚度矩阵的计算非常简单。假设三角形三个顶点坐标为 (x1,y1)、(x2,y2)、(x3,y3)面积 A则局部刚度矩阵的每个元素是K_ij (b_i b_j c_i c_j) / (4A)其中 b 和 c 是由顶点坐标差构成的系数具体为b1 y2 - y3b2 y3 - y1b3 y1 - y2 c1 x3 - x2c2 x1 - x3c3 x2 - x1这个公式就是 P1 单元的刚度矩阵。再加上介电常数每个单元内部 ε 是常数所以直接把 ε 乘到该单元的 Ke 上组装到全局矩阵时自然就把不同材料区域区分开了。材料分界面不需要额外处理因为相邻单元共享节点电位在界面上天然连续而场强法向分量的跳变通过节点电位的求解自动满足。3. 手写 Matlab 求解器不用 PDE 工具箱的完整实现3.1 完整可运行代码下面这一段就是完整的主程序保存成 main_capacitor_fem.m直接运行即可。版本上建议 R2017 之后的 Matlab我用到的都是基础函数只要矩阵运算够用就行。% main_capacitor_fem.m % 电容器内部区域有限元仿真二维静电场三角形P1单元 clear; clc; close all; %% 1 几何与材料参数 V0 1; % 施加电压归一化到1V不会影响电场分布形态 W 20e-3; % 极板宽度 d 2e-3; % 极板间距 Hpad 4*d; % 上下方向外部扩展距离 Lx W 8*d; % 求解域总宽左右各外扩4d Ly d 2*Hpad; % 求解域总高 xL (Lx - W)/2; xR (Lx W)/2; yL Hpad; yU Hpad d; epsrInside 2.2; % 内部填充介质相对介电常数 eps0 8.8541878128e-12; %% 2 网格生成结构化三角网格 Nx 160; % x方向单元数 Ny round(Nx * Ly / Lx); % y方向单元数保持单元接近正方形 x linspace(0, Lx, Nx1); y linspace(0, Ly, Ny1); nnX Nx 1; nnY Ny 1; % 节点数 n nnX * nnY; nodeIdx reshape(1:n, nnX, nnY); % 节点编号矩阵第一维对应x Xg zeros(n,1); Yg zeros(n,1); for j 1:nnY for i 1:nnX id nodeIdx(i,j); Xg(id) x(i); Yg(id) y(j); end end tri zeros(2*Nx*Ny, 3); cnt 0; for j 1:Ny for i 1:Nx n1 nodeIdx(i,j); n2 nodeIdx(i1,j); n3 nodeIdx(i1,j1); n4 nodeIdx(i,j1); cnt cnt 1; tri(cnt,:) [n1 n3 n4]; cnt cnt 1; tri(cnt,:) [n1 n2 n3]; end end %% 3 单元材料标识与全局刚度矩阵组装 triCx (Xg(tri(:,1)) Xg(tri(:,2)) Xg(tri(:,3)))/3; triCy (Yg(tri(:,1)) Yg(tri(:,2)) Yg(tri(:,3)))/3; inside triCx xL triCx xR triCy yL triCy yU; epsrElem ones(size(tri,1),1); epsrElem(inside) epsrInside; K sparse(n, n); triArea zeros(size(tri,1),1); for e 1:size(tri,1) ids tri(e,:); xe Xg(ids); ye Yg(ids); A 0.5 * abs((xe(2)-xe(1))*(ye(3)-ye(1)) - ... (xe(3)-xe(1))*(ye(2)-ye(1))); triArea(e) A; b [ye(2)-ye(3), ye(3)-ye(1), ye(1)-ye(2)]; c [xe(3)-xe(2), xe(1)-xe(3), xe(2)-xe(1)]; Ke (b*b c*c) / (4*A) * epsrElem(e); K(ids, ids) K(ids, ids) Ke; end %% 4 边界条件消去法 tol 1e-12; % 上极板节点y yUx在[xL, xR]内 upId find(abs(Yg - yU) tol Xg xL - tol Xg xR tol); % 下极板节点y yLx在[xL, xR]内 dnId find(abs(Yg - yL) tol Xg xL - tol Xg xR tol); fixId [upId; dnId]; fixVal [V0*ones(numel(upId),1); zeros(numel(dnId),1)]; freeId setdiff((1:n), fixId); V zeros(n,1); V(fixId) fixVal; V(freeId) K(freeId, freeId) \ (-K(freeId, fixId) * fixVal); %% 5 后处理电场、绘图、能量 nTri size(tri,1); E zeros(nTri,2); % 每个三角形内的电场矢量 Emag zeros(nTri,1); for e 1:nTri ids tri(e,:); xe Xg(ids); ye Yg(ids); A triArea(e); b [ye(2)-ye(3), ye(3)-ye(1), ye(1)-ye(2)]; c [xe(3)-xe(2), xe(1)-xe(3), xe(2)-xe(1)]; gradVx (b * V(ids)) / (2*A); gradVy (c * V(ids)) / (2*A); E(e,:) [-gradVx, -gradVy]; Emag(e) sqrt(E(e,1)^2 E(e,2)^2); end % 图1电位云图 等势线 figure(Name,电位分布); patch(Faces,tri,Vertices,[Xg Yg],FaceVertexCData,V, ... FaceColor,interp,EdgeColor,none); axis equal; colorbar; colormap(parula); hold on; Xs linspace(0, Lx, 240); Ys linspace(0, Ly, 240); [XX, YY] ndgrid(Xs, Ys); F scatteredInterpolant(Xg, Yg, V); contour(XX, YY, F(XX, YY), 20, k); plot([xL xR], [yU yU], r-, LineWidth, 2); plot([xL xR], [yL yL], r-, LineWidth, 2); title(电位分布V0 1 V); % 图2电场强度云图 figure(Name,电场强度); patch(Faces,tri,Vertices,[Xg Yg],FaceVertexCData,Emag, ... FaceColor,interp,EdgeColor,none); axis equal; colorbar; colormap(parula); hold on; plot([xL xR], [yU yU], r-, LineWidth, 2); plot([xL xR], [yL yL], r-, LineWidth, 2); title(电场强度 |E| 分布V/m); % 图3中心竖直剖面电位 [~, ic] min(abs(x - (xLxR)/2)); figure(Name,中心线电位剖面); plot(y, V(nodeIdx(ic,:)), b-, LineWidth, 1.5); xlabel(y (m)); ylabel(V (V)); grid on; title(x W/2 处的电位剖面); %% 6 能量积分与等效电容 We 0.5 * eps0 * sum(epsrElem .* Emag.^2 .* triArea); Csim 2 * We / V0^2; % 仿真得到单位深度电容F/m Cplate eps0 * epsrInside * W / d; % 平板公式忽略边缘场 fprintf(总储能 We %.4e J/m\n, We); fprintf(仿真等效电容 Csim %.4e F/m\n, Csim); fprintf(平板公式电容 Cplate %.4e F/m\n, Cplate); fprintf(边缘电容增量 %.1f%%\n, (Csim/Cplate - 1)*100);运行完成后终端会显示储能、等效电容和边缘电容增量同时弹出三张图。我建议第一次跑的时候用默认参数等流程走通了再改 Nx、W、d 等参数做敏感性测试。3.2 几何与网格设计为什么这样设参数先解释一下求解域的选取。极板宽度 W 20 mm间距 d 2 mm这是一个典型的 W/d 10 的工况。求解域总宽设成 W 8d也就是左右各外扩 4d总高设成 d 8d上下各外扩 4d。这个外扩距离不是随便拍的它的意义是让截断边界远离极板保证外边界上的 ∂V/∂n 0 近似成立。如果外扩距离太小人工边界会截断边缘场导致电场线被“逼回”求解域能量偏大外扩太大网格数量白白增加求解变慢。极板本身建模成零厚度线段而不是带厚度的矩形。这样做有两个好处一是几何简单网格生成和材料分配都能用结构化网格处理二是在极板间距远大于极板厚度时厚度对内部场分布的影响往往可以忽略。如果以后要研究极板厚度、卷绕结构或者圆角半径再在这个基础上改几何即可。网格方面我用了结构化矩形剖分每个小矩形再切成两个三角形。选择三角形而不是四边形是因为三角形单元能更好地适应不规则边界而且 P1 线性单元的刚度矩阵公式非常简单很适合用来讲清楚 FEM 的核心逻辑。Nx 和 Ny 的比例按求解域长宽比设置让每个单元尽量接近正方形。这个细节容易被忽略单元形状太扁时刚度矩阵条件数变差求解误差会上升。3.3 组装与材料标识几个容易出错的细节组装全局刚度矩阵是整个代码的核心。每个三角形单元先算出局部 3×3 矩阵 Ke再把 Ke 累加到全局矩阵 K 中对应节点编号的位置上。这里有一个实际经验写小规模仿真图方便直接K(ids, ids) K(ids, ids) Ke是可以的但当网格达到几十万单元时这种边循环边修改稀疏矩阵的写法会明显变慢。生产环境建议改成“先收集所有单元的 I、J、S 三列向量最后一次性 sparse”的向量化组装。本文为了可读性保留了最简单的循环版本跑 20 万节点以内基本无压力。材料分配用的是单元质心判断法。每个三角形算出质心坐标如果质心落在 xL x xR 且 yL y yU 的矩形内就认为这个单元属于介质填充区域εr 取 2.2否则属于空气εr 取 1。这个方法在极板为规则矩形、且网格线刚好落在极板边界线上时非常干净。需要注意极板边界线 yL、yU、xL、xR 必须正好是你网格节点所在的位置否则会出现“极板边界把三角形切成两半”的尴尬情况。本文代码里生成的 x、y 向量是均匀分布的只要 xL、yL 等参数与网格比例匹配就能保证这一点。3.4 边界条件消去法与后处理代码里的边界条件处理用的是消去法。基本思路是先把固定节点上下极板上的电位值写进解向量再从刚度矩阵中删掉这些节点的方程求解剩下的自由节点。公式上就是V_free K_ff \ (-K_fs * V_fixed)其中 K_ff 是自由节点对应的子矩阵K_fs 是自由节点和固定节点之间的耦合子矩阵。这种做法比罚函数法更精确不需要调惩罚因子而且 K_ff 是对称正定矩阵反斜杠求解又快又稳。后处理部分电场强度是我自己从三角形单元内算的没有用 gradient 函数对网格做数值差分。原因很简单P1 单元里电位线性变化电场在每个三角形内本来就是常数直接用单元节点电位算梯度结果最精确。如果你用 gradient 函数对已求解的节点电位做差分会在材料界面和极板边缘处引入额外误差尤其是场强变化剧烈的区域。这一步差别很大建议沿用我代码里的写法。4. 仿真实测与结果解读边缘场如何改变你的设计认知4.1 电位云图与电场集中区域跑完代码第一张图是电位分布。可以看到极板正对区域的等势线近似平行这是均匀场的特征但一靠近极板边缘等势线开始剧烈弯曲尤其在上极板右上角、右下角和下极板对应位置等势线明显加密。等势线密集意味着电位梯度大也就是电场强。第二张图是电场强度云图。颜色最亮的位置出现在极板边缘的外侧而不是极板正对区域的中心。这是一个常见误区很多人觉得平行板电容器的最大场强在内部均匀区实际上只要极板宽度不是远远大于间距边缘场强很快就会超过内部场强。在我这个算例里极板边缘附近的场强相对内部均匀场能高出 2 到 4 倍具体数值取决于网格密度和是否引入圆角半径。电场从边缘向外衰减也很快通常在离开极板边缘 2 到 3 倍极板间距的距离上场强就衰减到均匀场的十分之一以下。这说明边缘效应虽然局部剧烈但影响范围有限。这个空间尺度对绝缘设计很有用你想在极板边缘加绝缘挡板或者加宽安全距离大概需要预留多少空间从仿真图里可以量出来。4.2 能量积分与边缘电容增量储能积分是验证代码正确性的一个重要手段。二维模型按单位深度算总储能 We 0.5 × ε0 × ∫ εr · |E|² dA等效电容 Csim 2We / V0²。把 Csim 和平板公式 Cplate ε0 εr W/d 对比就能得到边缘场的“电容增量”。我这里不同网格密度下的计算结果如下表Nx 是 x 方向单元数Ny 按求解域长宽比自动生成网格 Nx×Ny近似边缘峰值场强/内部均匀场边缘电容增量80×22约 2.618.3%160×44约 3.218.0%320×89约 3.917.9%可以看到两个重要现象。第一边缘电容增量是收敛的网格加密后稳定在 18% 左右。这说明对于总电容这类“全局量”FEM 结果可信度很高网格密度的影响很小。第二边缘峰值场强却在随网格加密而持续增大从 2.6 涨到 3.9没有收敛的趋势。这个现象后面会专门讲它是极板建模成零厚度直角尖角带来的几何奇异性不是代码 bug。4.3 中心线电位剖面与正确性自检第三张图是 x W/2 处的竖直电位剖面。从下极板到上极板电位从 0 线性上升到 V0斜率约为 V0/d。这是一个很好的自检点如果这条曲线明显偏离直线说明网格太粗、边界太近或者材料参数设置有问题。注意中心线只是极板正对区域的剖面这里电场最接近均匀场。真正能反映边缘效应的剖面是在极板边缘外侧再取一条竖直或水平线观察电位如何平滑过渡到外部空气区。我建议你跑完代码后自己加一行在 x xR 0.5d 的位置画一条竖直剖面会看到电位从极板电位迅速向远处衰减衰减速率比中心线陡得多这同样是边缘场集中程度的直接证据。5. 工程实操经验边界截断、材料分界面和几何奇异5.1 外部截断边界到底该设多远外部域大小是这类开放边界问题里最容易踩的坑。外扩太小人工边界会把边缘电场“压扁”导致电容偏大、峰值场强失真外扩太大网格规模爆炸求解时间拉长。我的经验判断方法是先在一个较小的外扩距离下算一次 Csim再把外扩距离翻倍看 Csim 变化多少。如果变化小于 1%说明截断足够如果还在明显变化继续加大。本文用的 Lx W 8d、Hpad 4d 是在“计算精度”和“网格规模”之间比较平衡的选择。如果只关心极板正对区域内部的场强外扩 2d 到 3d 就够了但如果要评估边缘场对周围器件的影响或者要算杂散电容我建议至少外扩 5d 以上。从代码角度看外扩距离只是 Lx、Ly、Hpad 几个参数扫描一遍非常方便。5.2 介质分界面的网格布置与参数保存内部介质区域和外部空气区域的交界是电容仿真里最容易产生误差的位置。本文因为分界面是水平直线布置成网格线就能保证三角形单元不会跨在分界面上。如果换成斜边界或者圆弧边界结构化网格就不够用了需要上非结构化网格或者对分界面局部加密否则场强最大值会严重依赖网格方向。另一个建议是把 epsrElem、triArea 这些参数保存到结果文件里不要等后处理时再重新判断。原因很简单能量积分、等效电容这些后处理操作都要用到每个单元的材料参数如果重新判断材料区域很容易因为容差问题造成前后不一致。代码里我直接在组装阶段就把 epsrElem 算好并参与刚度矩阵组装后处理时再拿出来用数据链是一致的。5.3 零厚度尖角的奇异性与圆角化建议前面说的边缘峰值场强不收敛本质上是数学上的尖角奇异性。二维静电场中一个理想的零厚度导电尖角其边缘的电场理论上是无穷大的。有限元网格无论加密到什么程度这个奇点附近的场强都会持续增大只是增大的速度会越来越慢。所以如果你在仿真报告里写“最大场强为某值”一定要注明网格密度否则这个数字没有可比性。工程上有两个应对思路。第一个是给极板边缘设计圆角让几何上不再存在尖角。圆角半径哪怕只有 0.1d 到 0.3d也能显著降低边缘峰值场强。第二个是评估时不要只看局部峰值而是看能量分布或者某个固定采样路径上的场强积分。峰值发散不代表总储能发散Csim 收敛就是证据。实际设计做对比时用能量或电容增量来评判方案比用单一最大值可靠得多。最后再说一个性能相关的小技巧代码默认跑 Nx 160 的网格正常电脑几秒就能出结果。如果你要把 Nx 拉到 640 以上做网格无关性验证建议把组装循环改成向量化稀疏矩阵一次性构建同时在后处理电场计算里用矩阵运算替代 for 循环速度能快十几倍。我自己做参数扫描时常设 Nx 320跑一趟大约一两分钟配合这几十组几何参数足够把边缘电容增量和极板宽度的关系摸清楚了。
阅读完成 · 觉得有帮助？
