企业官网建设流程全解析

前一阵调试一台自制高压试验装置，耐压测试做到一半，能听到腔体里有细微的嘶嘶声。把设备拆开，放电的痕迹不在极板正对区域，而是沿着上下极板的外侧边缘，一条发亮的细线。这个现象在设想里其实出现过：平行板电容器内部的电场，教科书上写 E=U/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 - y3，b2 = y3 - y1，b3 = y1 - y2 c1 = x3 - x2，c2 = x1 - x3，c3 = 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, Nx+1); y = linspace(0, Ly, Ny+1); 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(i+1,j); n3 = nodeIdx(i+1,j+1); n4 = nodeIdx(i,j+1); 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 = yU，x在[xL, xR]内 upId = find(abs(Yg - yU) < tol & Xg > xL - tol & Xg < xR + tol); % 下极板节点：y = yL，x在[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 - (xL+xR)/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 按求解域长宽比自动生成：

可以看到两个重要现象。第一，边缘电容增量是收敛的，网格加密后稳定在 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，跑一趟大约一两分钟，配合这几十组几何参数，足够把边缘电容增量和极板宽度的关系摸清楚了。