MATLAB实现模拟电荷法：输电线路电场计算与校验

简介针对电力系统架空输电线路电磁场分析中连续电荷积分求解困难、离散近似精度待校验的问题这份基于MATLAB的模拟电荷法验证程序为电气工程专业学生、科研人员及电磁场数值计算入门者提供了可直接运行的测试工具。资源包共1个文件即核心MATLAB脚本文件整体压缩包仅2KB代码精简、便于逐行阅读与二次扩展。该脚本先假设导线中心存在一个或多个模拟电荷通过高斯定律计算空间电场再在导线表面选取覆盖多个横截面的校验点对比理论值或实测值并统计误差从而全面评估模型准确性。同时融入可视化展示环节绘制电场、电荷及误差分布图帮助使用者直观掌握算法实现流程和排错思路亦可作为课程设计或数值算法对比的参考样例。目前已有170人学习下载适合需要借助轻量级代码快速上手模拟电荷法编程的读者。1. 模拟电荷法为什么值得用 MATLAB 再验证一遍做输电线路电磁环境计算的人大多绕不开这样一个矛盾用有限元软件跑一个档距的导线电场网格剖分和收敛调试要占掉大半天而换用模拟电荷法CSM把导线表面连续电荷离散成一组等效线电荷半小时就能拿到量级正确的电场分布。问题是离散电荷放在哪、放多少、校验点怎么排直接决定结果可不可信。这个 check.rar 里的 check.m就是干这件事的给定导线半径、离地高度和相电压自动完成电荷求解、表面校验点误差计算与结果可视化。它适合正在学电磁场数值方法的同学对照公式验证也适合做输变电工程电磁环境预评估的工程师快速建立基线模型。2. 离散电荷模型与 check.m 的求解流程2.1 为什么要把面电荷换成内部线电荷电荷实际上分布在高电压导线的表面但直接求解表面连续分布需要解边界积分方程计算量大且边界复杂时难以收敛。模拟电荷法换了一个思路把表面连续电荷用一组放在导体内部的虚拟电荷去近似。理论上只要虚拟电荷严格位于求解域之外且真实导体表面电位为常数Laplace 方程解的唯一性就保证了这种替代可行。对架空导线这种细长结构二维截面下等效源就是一组无限长线电荷。线电荷在空间任意一点产生的电位与到该线电荷的距离取对数关系多根线电荷叠加后就得到整个截面的电位分布。之所以要把电荷放在导线内部而不是表面是为了避免电位公式在表面出现奇异性同时保证源所在位置不属于需要满足边界条件的求解域。2.2 电位系数矩阵构造与镜像法处理架空导线下方是大地工程上把大地近似为零电位平面。处理办法是镜像法在导线关于地面的对称位置放置等量异号电荷这样地面 y0 处电场严格为零。实际计算里所有电位和电场都是“真实源 镜像源”叠加后的结果。设真实线电荷位置为 (xs,ys)镜像电荷位置为 (xs,-ys)校验点位置为 (xp,yp)。在二维场中单位长度线电荷 Q_j 在目标点产生的电位系数为P_ij (1 / 2π ε0) · ln( ρ_ij′ / ρ_ij )其中 ρ_ij 是真实源到目标点的距离ρ_ij′ 是镜像源到目标点的距离。对数比值天然满足地面零电位当地面点 y0 时两个距离相等比值为 1电位为零。构造电位系数矩阵的 MATLAB 片段如下% check.m 电位系数矩阵构造核心循环 % 输入: xp,yp 校验点坐标; xs,ys 模拟电荷坐标; eps0 真空介电常数 function P build_P(xp, yp, xs, ys, eps0) m length(xp); % 校验点数量 nq length(xs); % 模拟电荷数量 P zeros(m, nq); for i 1:m for j 1:nq rho_ij sqrt((xp(i)-xs(j))^2 (yp(i)-ys(j))^2); % 镜像电荷在 (xs(j), -ys(j))注意距离是 ypys 而非 yp-ys rho_ij_p sqrt((xp(i)-xs(j))^2 (yp(i)ys(j))^2); P(i,j) log(rho_ij_p / rho_ij) / (2*pi*eps0); end end end这里的循环写法在性能上不是最优但胜在直观电磁场数值计算里矩阵规模通常在几十乘几十两层循环完全够用。要注意镜像距离是关键写错符号会导致地面电位不为零校验时误差会整体偏大。另一种常见做法是直接构造二维坐标网格后向量化计算但可读性差不推荐在初版验证脚本里用。2.3 主流程与关键变量check.m 的主流程分四步先给物理参数再布置模拟电荷然后构造匹配点的电位系数矩阵并解线性方程组最后加校验和可视化。下面给出主脚本的可运行骨架% check.m 主脚本可独立运行 clc; clear; close all; % 1. 物理参数 R 0.013; % 导线半径单位 m约 220kV 线路常用导线 H 12.0; % 导线离地高度单位 m V0 220e3 / sqrt(3); % 相电压有效值单位 V eps0 8.854187817e-12; % 真空介电常数F/m % 2. 模拟电荷布置 nq 8; % 模拟电荷数量 th_s (0:nq-1) / nq * 2 * pi; % 周向均匀分布 xs 0.5 * R * cos(th_s); % 源电荷水平坐标 ys H 0.5 * R * sin(th_s); % 源电荷竖直坐标 % 3. 匹配点与线性求解 nm nq; % 匹配点数与电荷数一致 th_m (0:nm-1) / nm * 2 * pi; % 匹配点角度 xm R * cos(th_m); ym H R * sin(th_m); % 匹配点坐标导线表面 Pm build_P(xm, ym, xs, ys, eps0); % 电位系数矩阵 Vb ones(nm, 1) * V0; % 表面电位边界条件 Q Pm \ Vb; % 解出单位长度电荷量C/m % 4. 结果输出 fprintf(单位长度电荷量: %.3e C/m\n, sum(Q));匹配点取在导线表面、模拟电荷取在半半径圆上这个几何关系决定了源到表面最近距离是 0.5R电位系数不会出现奇异值。Pm \ Vb用的是 MATLAB 内置的 LU 分解求解对几十阶矩阵不需要关心性能。唯一要注意的是 Vb 的单位这里用伏特对应介电常数要用国际单位最后电荷量单位是 C/m。变量含义单位典型取值R导线半径m0.005~0.03H导线离地高度m6~30V0导线对地电位V相电压有效值nq模拟电荷数个4~16Q单位长度电荷量C/m由方程解出这个脚本不依赖任何扩展工具箱只要有基础 MATLAB 环境就能跑。实际工程里如果把导线视为无限长直导线二维截面下的单位长度电荷量乘以线路长度就是该段导线的总等效电荷可以直接用于后续工频电场计算。3. 校验点策略与表面电场误差分析3.1 匹配点和校验点要分开初学者最容易犯的错是把求解线性方程组用过的匹配点再一次用来评估误差这样得出的残差几乎为零没有任何参考价值。check.m 里正确的做法是把求解和校验分成两套点集匹配点用于构建边界条件方程数量与模拟电荷数相等校验点用于独立评估精度数量可以是电荷数的 2 到 4 倍。校验点均匀分布在导线表面圆周上角度间隔要小于模拟电荷角度间隔的一半否则容易漏掉电荷密度变化最剧烈的方向。如果手头有上一版计算数据一个可用的快速判断准则是把校验点角度旋转半个间隔后重新算一遍如果最大误差变化超过一个数量级说明校验点排得太疏或者模拟电荷数偏少。3.2 电位残差与表面电场计算除了电位模拟电荷法还需要验证表面电场是否满足导体边界条件。导体的电动力学边界条件要求表面切向电场为零法向电场与表面电荷密度成正比。由于电场是电位的空间导数计算出来的电场误差通常比电位误差大一个量级所以表面电场校验更能暴露模型的缺陷。校验部分代码如下% 独立校验点电场计算 mc 48; % 校验点数量 th_c (0:mc-1) / mc * 2 * pi; xc R * cos(th_c); yc H R * sin(th_c); Pc build_P(xc, yc, xs, ys, eps0); % 校验点电位系数 Vc Pc * Q; % 校验点电位 dV abs(Vc - V0) / V0 * 100; % 相对电位误差(%) % 表面电场解析计算真实源与镜像源叠加 Ex zeros(mc,1); Ey zeros(mc,1); for k 1:mc for j 1:nq % 真实源贡献正电荷向外 rx xc(k) - xs(j); ry yc(k) - ys(j); r2 rx^2 ry^2; Ex(k) Ex(k) Q(j) / (2*pi*eps0) * rx / r2; Ey(k) Ey(k) Q(j) / (2*pi*eps0) * ry / r2; % 镜像电荷贡献负电荷 rxm xc(k) - xs(j); rym yc(k) ys(j); r2m rxm^2 rym^2; Ex(k) Ex(k) - Q(j) / (2*pi*eps0) * rxm / r2m; Ey(k) Ey(k) - Q(j) / (2*pi*eps0) * rym / r2m; end end % 切向电场圆周切向量为 (-sin, cos) Et -Ex .* sin(th_c) Ey .* cos(th_c); En Ex .* cos(th_c) Ey .* sin(th_c); fprintf(电位最大相对误差: %.4e %%\n, max(dV)); fprintf(表面切向电场峰值: %.4e V/m\n, max(abs(Et))); fprintf(表面法向电场均值: %.4e V/m\n, mean(abs(En)));电场叠加时要注意镜像电荷符号为负计算点对真实源和镜像源分别用各自距离求向量场再叠加。这里不能直接对电位梯度做数值差分因为离散电荷在靠近导线表面时电位梯度变化剧烈数值差分误差会掩盖真实误差。判断结果是否可信经验阈值大致如下指标可接受范围说明电位最大相对误差 0.1%越大说明等效源逼近越差切向电场峰值 / 法向电场均值 1%超过后需增加电荷数或调整位置校验点间最大电位差 0.5%表面应近似等电位偏差大说明源布置不当3.3 校验失败的常见原因如果电位残差始终降不下来先检查模拟电荷是否离表面太近。模拟电荷越靠近表面等效精度越高但电位系数矩阵条件数会急剧恶化求解结果出现高频振荡。此时残差不降反升而且电荷量符号正负交替、数值忽大忽小。最直接的办法就是在 MATLAB 命令行里看一眼cond(Pm)超过 1e8 就要把源环半径从 0.5R 往下调或减少模拟电荷数。另一个高频原因是对称性破坏。单根导线场景下模拟电荷和校验点都应为圆周对称分布如果角度序列从 0 开始而电荷数不是 4 的倍数电场在某个特定方向会出现系统性偏差极坐标误差图上表现为单峰或双峰结构。遇到这种情况把模拟电荷的角度整体旋转半个间隔通常能明显改善。4. MATLAB 可视化云图、矢量图与误差诊断4.1 电位云图与电场矢量图MATLAB 可视化的价值不只是把结果画出来而是把“电荷布置是否合理”这件事变得肉眼可见。电位云图可以直观看到等位线的疏密电场矢量图可以看出导线附近电场集中程度。实现代码% 场域电位云图与电场矢量图 xr linspace(-2, 2, 150); % 横向范围单位 m yr linspace(0.5, 2*H, 150); % 纵向范围不含地表面 [Xg, Yg] meshgrid(xr, yr); Vg zeros(size(Xg)); % 对网格点逐点叠加真实源与镜像源 for i 1:numel(Xg) for j 1:nq rho sqrt((Xg(i)-xs(j))^2 (Yg(i)-ys(j))^2); rho_p sqrt((Xg(i)-xs(j))^2 (Yg(i)ys(j))^2); Vg(i) Vg(i) Q(j) * log(rho_p / rho) / (2*pi*eps0); end end figure(Color,w); contourf(Xg, Yg, Vg, 60); colorbar; axis equal; hold on; [Exg, Eyg] gradient(-Vg); % 网格差分粗略估计电场 quiver(Xg(1:12:end,1:12:end), Yg(1:12:end,1:12:end), ... Exg(1:12:end,1:12:end), Eyg(1:12:end,1:12:end), k); plot(0, H, ko, MarkerFaceColor,k); title(电位云图与电场矢量);逐点双重循环在 150×150 网格上耗时约一两秒对验证脚本完全够用。gradient对电位做中心差分得到的是粗略电场估计只用于定性显示定量分析仍然要用第 3 章的解析公式。另一个细节是坐标范围横轴 2m 在高电压等级线路里已经覆盖了导线附近的主要场区纵轴从 0.5m 到 2H 足够看到等位线向下半空间的延伸。4.2 表面误差极坐标图表面电位残差随角度的变化用极坐标图呈现比直角坐标图直观得多。如果模拟电荷布置对称误差曲线也对称如果某几个角度误差明显突出说明电荷在这些方向上的等效能力不足。绘制代码% 表面电位相对误差极坐标图 figure(Color,w); polarplot(th_c, dV, o-, LineWidth, 1.2); rlim([0, max(dV)*1.2]); title(表面电位相对误差随角度变化 (%));运行这段脚本前要注意th_c是弧度制的校验点角度dV是百分比数值两者长度必须一致。极坐标图不适合显示量级差异过大的数据如果最大误差是最小误差的十倍以上先检查数据而不是直接画图。误差曲线呈周期振荡通常意味着模拟电荷数不足需要回到第 2 章把 nq 调大重算。4.3 误差曲线用于参数对比把 nq 做成扫描参数每一组电荷数跑一遍校验记录电位最大相对误差和切向电场峰值最后用semilogy画成两条曲线。一般会看到误差随电荷数增加先下降到某个值后平台期再继续加电荷误差反而反弹。反弹点就是矩阵病态的起点对应cond(Pm)开始急剧上升的位置。这个扫描脚本每次只需重跑几十次线性方程求解总耗时在几秒内却能把“模拟电荷最优数量”这个经验问题变成一个可量化的选择。5. 从单导线到分裂导线的模型扩展技巧实际输电线路上分裂导线非常常见四分裂导线由四根子导线组成等间距四边形子导线半径约 0.02m间距约 0.4m。check.m 的单线模型扩展过去核心是把模拟电荷的布置逻辑从“围绕导线中心”改成“围绕每根子导线中心”其余代码几乎不用动。修改思路% 四分裂导线模拟电荷布置示意 r_sub 0.0135; % 子导线半径 spacing 0.4; % 分裂间距 cx spacing/2 * [ 1, 1, -1, -1]; % 子导线水平位置 cy spacing/2 * [ 1, -1, 1, -1]; % 子导线竖直位置 n_sub 4; n_per 4; % 子导线数量与每根子导线内电荷数 xs zeros(n_sub * n_per, 1); ys zeros(n_sub * n_per, 1); for s 1:n_sub idx (s-1)*n_per 1 : s*n_per; th (0:n_per-1) / n_per * 2 * pi; xs(idx) cx(s) 0.5 * r_sub * cos(th); ys(idx) H cy(s) 0.5 * r_sub * sin(th); end采集电位系数矩阵、求解、校验这段逻辑保持原样即可。要注意的是镜像法的前提是地面为零电位面分裂导线整体抬高后镜像位置仍然按每根模拟电荷的 y 坐标取负不需要额外处理。工程调试中如果遇到校验误差不达标还可以做一个快速检查对同一模型分别用 8、12、16 个模拟电荷计算看着cond(Pm)和电位残差的变化方向。条件数增大但残差也增大说明模型本身有问题条件数增大但残差继续下降说明正在接近最优解区间。最后的落点建议做一个电荷数扫描图把max(dV)和cond(Pm)画在同一张对数坐标图里最优电荷数就是两条曲线交汇点附近的值。本文还有配套的精品资源点击获取