企业官网建设流程全解析

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。网格取nx=100、ny=60，这样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, nx+1); y = linspace(0, Ly, ny+1); [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)*(nx+1) + i; % 下左 n2 = (j-1)*(nx+1) + i + 1; % 下右 n3 = j*(nx+1) + i + 1; % 上右 n4 = j*(nx+1) + 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) <= 8+1e-12 & ... nodes(:,2) >= 2.5-1e-12 & nodes(:,2) <= 2.8+1e-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) <= 8+1e-12 & ... nodes(:,2) >= 3.2-1e-12 & nodes(:,2) <= 3.5+1e-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, ny+1, nx+1); % 绘制等势线图 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 边缘效应不是玄学，用等势线密集程度说话

运行代码后，第一张图就是等势线分布。你会看到极板正中间区域的等势线近似水平、间距均匀，这和平行板公式给出的均匀场E=V/d一致。但在极板左右两端的x=2和x=8附近，等势线明显向下和向上弯曲并聚拢，说明电场线不再平行，而是从极板侧面绕过去形成“边缘场”。

边缘效应最需要关注的是场强峰值。用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)参与刚度矩阵累加。

我代码里默认把极板间隙填充为εr=4的材料，这里可以做一个更有意思的实验：把间隙改成两层，y方向中间2.8到3.0处εr=3，3.0到3.2处εr=1。运行后看电场分布，你会发现ε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计算结果是离散近似，网格越细结果越接近真实解，但计算量和内存也会上升。工程上做网格无关性验证的思路是：取一组逐渐加密的网格，计算某个关键量（比如极板边缘的最大电场强度），看它是否趋于稳定。

以我手头这套模型为例，我跑到过以下几组数据：

可以看到中部场强很快收敛到2.5左右，这和理论值1/0.4=2.5吻合，但边缘峰场强还在缓慢上涨，因为边缘处的场强理论上存在奇异性，网格越细越能逼近峰值。如果只是想比较不同方案的相对优劣，网格密度到100×60左右已经能给出稳定判断；如果要做绝对击穿风险评估，还需要进一步加密，并结合实际材料缺陷尺寸来讨论收敛判据，不能只看一个点的数值。

5.3 从静电场到瞬态场和AC损耗分析的扩展路径

二维静电场FEM是很多问题的基础，但不是终点。如果电容器工作在高频下，介质损耗和导体损耗会成为主要关注点，这时候控制方程要换成频域波动方程或复数介电常数形式，求解变量从实数的电势变成复数的相量电势。MATLAB的PDE Toolbox也有对应的谐波求解器，手写代码则需要在刚度矩阵中引入复介电常数，整体组装逻辑不变，只是把实数稀疏矩阵变成复数稀疏矩阵。

如果要研究瞬态充电过程，就在控制方程中加入时间导数项，空间离散后得到半离散的常微分方程组，再用ode23或ode45推进时间。这里面一个常见坑是Courant条件对时间步长的限制，步长太大结果会振荡，步长太小计算量又上不去，需要根据自己的网格尺寸调。还有同学会把电场FEM和电路仿真联立，用Simulink做外部充放电回路，电场求解器作为每个时间步的“子程序”被调用，这是把场路耦合做进系统级仿真的典型思路。

我自己跑这类代码时，最后都会做一件事：把极板间距、介质厚度、电极厚度全部参数化，写成一个函数，留出输入接口，这样结构优化时只需要循环调用求解函数，用几行脚本就能批量扫描不同几何方案下的边缘场强峰值。仿真不是一次性的，能被人反复修改调用的代码，才算真正入了门。