简介一套基于 MATLAB 的模拟电荷法验证工具面向电力系统电磁场计算相关的研究人员与工程师用于精准检查架空输电线路中离散模拟电荷模型的计算精度与可靠性。该程序主要通过脚本实现电荷计算、校验点电场对比与误差分析可帮助用户快速评估模拟电荷法的适用性并改进模型参数适合已掌握电磁场基础、希望用数值方法辅助分析的读者。压缩包仅包含一个MATLAB脚本文件体积约2KB轻量实用已有170人浏览学习。借助 MATLAB 的矩阵运算与绘图功能用户可直观观察电场分布、电荷偏差及校验点误差不仅能验证算法正确性还能以此为模板调整模拟电荷的数量、位置及校验点布局服务于输电线路电磁环境预测等典型场景。1. 模拟电荷法配 check_order边界点顺序错了电场解就全乱把 CAD 导出的电极边界存成 txt再用 MATLAB 读进来第一件要做的事不是画图而是检查边界点的顺序。模拟电荷法的整个流程——匹配点布置、模拟电荷定位、电位系数矩阵组装——都默认点沿闭合边界按固定方向依次排列。顺序一旦错位方程组照样解得出电荷量把解代回边界验算时电位对不上局部还会出现虚假电场尖峰残差不大画出来却明显变形。check_order 就是为这个前置步骤准备的一段代码去重、最近邻重排、统一旋转方向。适合刚用 MATLAB 实现模拟电荷法的人也适合调试已有代码但结果一直对不上的工程师。文章先把两者的配合关系说清再给出一个能对解析解验证的最小算例。2. 模拟电荷法与 check_order 的配合把边值问题写成一堆线性方程组2.1 模拟电荷法数学基础虚拟电荷替代连续表面电荷静电场计算里导体内电场为零电荷只分布在导体表面表面电荷密度连续且未知。常规做法是划分有限元网格把整个场域剖分一遍模拟电荷法换了个思路把电极表面以下的区域看成电荷允许出现的地方用一组离散的虚拟电荷去近似表面上那层真实连续电荷的对外效果。整个场域内不再需要网格只保留两类离散点匹配点和电荷点。在场域内任意一点电位由全部虚拟电荷叠加而成。既然虚拟电荷是等效源就要求它们在电极表面的匹配点上满足同样的边界条件第 i 个匹配点上叠加出来的电位必须等于已知电极电位 V_i。每个电荷在匹配点产生的电位写成q_j * p_ij于是匹配点 i 上的方程是φ_i Σ p_ij * q_j。把所有匹配点的方程写在一起就是标准的稠密线性方程组A q b。这就是为什么模拟电荷法能绕开场域网格求解对象从偏微分方程变成了一个代数系统。匹配点取在电极表面上电荷点取在电极内部两者数量相同一一对应。问题规模通常只有几十到几百个未知量MATLAB 反斜杠直接解稠密矩阵完全够用。模拟电荷法的精度高度依赖电荷点布置而电荷点布置依赖匹配点顺序check_order 在这个链条上的位置就在最前面。2.2 电位系数常用的两种形式电位系数p_ij的表达式由虚拟电荷的几何形态决定MATLAB 实现里最常用的是二维无限长线电荷和三维点电荷两种。问题类型电位系数 p_ij电场 x 分量累加式二维无限长线电荷单位 C/m-ln(r_ij) / (2π ε0)Σ q_j (x - x_j) / (2π ε0 r_ij²)三维点电荷单位 C1 / (4π ε0 r_ij)Σ q_j (x - x_j) / (4π ε0 r_ij³)r_ij是第 i 个匹配点到第 j 个电荷的距离。二维公式里的对数项来自二维格林函数解特别适合平行导体剖面和轴对称旋转体的轴向剖面三维点电荷形式则适合任意空间排布代价是同等精度下需要的电荷数更多。两类系数在r_ij → 0时都奇异所以电荷点必须放在电极表面内侧匹配点放在电极表面外侧两者不能重合。这也正是后面调整 lambda 参数的物理来源。2.3 check_order 具体检查什么重复点、闭合顺序、旋转方向比如拿到一个类似check.rar的共享算例包里面那份边界点文件数据可能来自 CAD 导出、坐标测量机甚至上一版脚本的中间结果点的顺序未必可控。顺序乱了等价于把电极表面拆成一段错位折线匹配点之间出现交叉连接相邻关系失真。这件事不会让矩阵报错A依然是方阵反斜杠也能给出一组q但q的分布正负交替、数值偏大逼近效果只在个别匹配点成立离开匹配点就迅速失效。check_order 要做的事分三步去重、恢复闭合顺序、统一方向。去重解决首尾点重复或相邻重复恢复顺序一般用基于距离的贪心最近邻统一方向用有向面积的符号来判。有向面积的计算式为A2 Σ (x_i * y_{i1} - x_{i1} * y_i)下标N1回到1。A2 0对应逆时针A2 0对应顺时针。这个判据对凹边界同样成立因为它度量的是整体绕向不要求边界严格凸。3. 用 MATLAB 跑通 check_order 与模拟电荷法同轴圆柱最小算例3.1 check_order 完整实现去重、贪心最近邻、方向统一最小算例选同轴圆柱因为它的解析解可以用初等函数写出来任何参数调整都能立即看出对错。先给出 check_order 的完整实现再给主程序。function pts check_order(pts, mode) % CHECK_ORDER 将二维离散边界点重排为有序闭合序列 % pts : n×2 矩阵每行 (x, y) % mode : ccw 强制逆时针cw 强制顺时针默认 ccw if nargin 2 mode ccw; end % 1. 去重首尾重复或同点被采多次先全部合并 [~, idx] unique(round(pts * 1e10), rows, stable); pts pts(idx, :); % 2. 最近邻重排从质心最近点出发贪心找下一个点 n size(pts, 1); if n 3 cx mean(pts(:, 1)); cy mean(pts(:, 2)); [~, start] min((pts(:,1)-cx).^2 (pts(:,2)-cy).^2); seq zeros(n, 1); seq(1) start; visited false(n, 1); visited(start) true; for k 2:n d2 sum((pts(~visited, :) - pts(seq(k-1), :)).^2, 2); [~, j] min(d2); tmp find(~visited); seq(k) tmp(j); visited(seq(k)) true; end pts pts(seq, :); end % 3. 统一旋转方向用有向面积判断 x pts(:, 1); y pts(:, 2); area2 sum(x .* circshift(y, -1) - circshift(x, -1) .* y); if strcmpi(mode, ccw) area2 0 pts flipud(pts); elseif strcmpi(mode, cw) area2 0 pts flipud(pts); end end代码里三个要点。去重用round(pts * 1e10)而不是直接对浮点数做unique因为边界点坐标带小数直接判等会因为截断误差漏掉重合点重排用最近邻贪心起点固定在质心最近的边界点上降低首尾接错的风险定向那一步里circshift(y, -1)把 y 整体上移一格下标N1自然回到 1正好对应闭合求和。提示贪心最近邻在凹边界上偶尔会“抄近路”跨过凹陷把闭合曲线连成自交叉。凸边界可以完全交给 check_order处理凹边界时配合按弧长参数化重采样更稳。3.2 主程序边界点、模拟电荷与电位系数矩阵组装主程序把同轴圆柱的内外电极各取 N 个匹配点模拟电荷放在电极内部角度错开半格再用二维无限长线电荷的电位系数组装矩阵。% CSM_DEMO.M 同轴圆柱静电场模拟电荷法最小算例 clear; clc; % 几何与边界条件单位m / V R1 0.01; % 内电极半径 R2 0.05; % 外电极半径 V1 100; % 内电极电位 V2 0; % 外电极电位 N 24; % 每圈匹配点数 % 内外边界匹配点先统一方向 th linspace(0, 2*pi, N1); th(end) []; P_in [R1*cos(th), R1*sin(th)]; P_out [R2*cos(th), R2*sin(th)]; P_in check_order(P_in, ccw); P_out check_order(P_out, ccw); % 模拟电荷角度错开半格半径按 lambda 回缩/外推 lambda 0.3; thC th pi/N; Q_in [R1*(1-lambda)*cos(thC), R1*(1-lambda)*sin(thC)]; Q_out [R2*(1lambda)*cos(thC), R2*(1lambda)*sin(thC)]; Q [Q_in; Q_out]; % 电位系数矩阵 Pm [P_in; P_out]; M 2 * N; A zeros(M, M); for i 1:M for j 1:M r2 sum((Pm(i, :) - Q(j, :)).^2); A(i, j) -0.5 * log(r2) / (2 * pi * 8.8541878128e-12); end end % 解出电荷量 Vb [V1*ones(N,1); V2*ones(N,1)]; q A \ Vb; fprintf(cond(A) %.3e\n, cond(A));主程序里几个参数值得说清。N是每圈匹配点数量总未知量是2Nlambda控制电荷离电极边界的远近越小越贴近表面拟合越细腻但矩阵越病态越大越远离边界矩阵稳定但拟合变钝0.2 到 0.5 是常规区间。thC故意错开半格让电荷点和匹配点不在同一角度上避免对称排列引起矩阵秩退化。电位系数直接实现二维无限长线电荷的-ln(r) / (2π ε0)用r2的平方形式省一次开方。还要提一句量纲ln(r2)里的r2带m²这里为了让代码短小直接由ε0吸收量纲工程代码建议先把坐标归一化到参考长度再组装矩阵避免不同量级坐标让条件数进一步恶化。3.3 用解析解校核检查点上的电场对比同轴圆柱的严格解是E(r) V1 / (r * ln(R2/R1))。把检查点放在中位半径圆上用叠加电荷计算电场再和解析解比较。% 检查点中位半径圆上取 100 个点 rMid (R1 R2) / 2; thc linspace(0, 2*pi, 100); Xc rMid * cos(thc); Yc rMid * sin(thc); % 用模拟电荷叠加电场 Ex zeros(100, 1); Ey zeros(100, 1); for j 1:M dx Xc - Q(j, 1); dy Yc - Q(j, 2); r2j dx.^2 dy.^2; coef q(j) ./ r2j / (2 * pi * 8.8541878128e-12); Ex Ex coef .* dx; Ey Ey coef .* dy; end En sqrt(Ex.^2 Ey.^2); % 解析解 Ea V1 / (rMid * log(R2 / R1)); fprintf(最大相对误差: %.3f%%\n, max(abs(En - Ea)) / Ea * 100);电场叠加写成显式循环比封装函数更直白每个检查点把 M 个电荷的贡献逐项累加。检查点放在中位半径处是有意的离匹配点一段距离能反映模型对外场的整体逼近而不是只看边界附近被方程组强制拉平的效果。N24、lambda0.3时最大相对误差通常能压到 0.5% 以内具体数值随 MATLAB 浮点环境略有浮动。4. 模拟电荷法参数怎么调电荷数量、位置比例和边界分布的取舍4.1 用检查点量化误差调整参数的前提是先把误差量化。判断一组模拟电荷解好不好不能只盯着匹配点上的电位残差因为匹配点电位被方程组强制满足残差天然很小。合理的做法是在场域内选一圈与电极边界不相交的检查点用模拟电荷解算出电位或电场再与参考解对比。没有解析解时把检查点加密一倍看结果是否稳定或者把检查点放在场域中轴附近观察场强是否光滑都能暴露局部拟合缺陷。检查点数量取 100 到 200 个足够。最大相对误差和平均相对误差分开记录最大误差反映最坏点平均误差反映整体水平调参时优先看最大误差因为工程上电场峰值往往出现在最容易失真的地方。4.2 电荷数与位置比例的参考取值N和lambda是两个最需要反复调的量而且互相制约。N决定未知量规模lambda决定矩阵病态程度。下面的经验量级来自类似结构的大量算例具体结构会有偏移但趋势稳定。每圈匹配点数 N总未知量 2N常见用途误差经验量级816结构筛选、快速估算百分之几1632常规电场计算1% 左右3264关注电场峰值的算例0.1% 附近64128高精度验证万分之一量级lambda的影响更直接它决定电荷点与边界的贴合程度。lambda行为特征建议0.05 ~ 0.1条件数可能到 1e15 以上电荷量正负交替不要用0.2 ~ 0.5精度与稳定性平衡最好首选区间0.6 ~ 0.8电荷远离电极边界拟合变钝误差回升少用开始一个新模型时先用N32、lambda0.3跑一遍再把两个参数分别向两边扫描观察最大误差的单调性。如果误差随N增大不降反升先查 check_order 有没有重排错再查检查点是否离电荷太近。4.3 三个常见坑与对应判断第一个坑是lambda取得太小。电荷点贴近匹配点时矩阵相邻两行几乎线性相关反斜杠仍能给出解但电荷量数值异常、正负交替cond(A)会显示 1e15 以上。判断依据就是条件数和电荷量分布正常解的电荷量同号且变化平滑异常解会出现大量反号。第二个坑是边界点分布不均匀。等角采样在圆形边界上没问题实际电极在曲率大的角点处点很密在直段很稀。电荷点跟随匹配点位置稀疏区域局部拟合误差被放大。解决方式是按等弧长重采样后再交给 check_order 重排方向。第三个坑是误以为 check_order 能纠正任意乱序。重复点能去、方向能统一但随机打乱的几千个点交给最近邻贪心只保证局部最优可能拼出与原边界相差很大的折线。最可靠的顺序信息来自 CAD 导出时沿轮廓的采样方向而不是事后让算法去猜。5. 从最小算例到实际电极三个值得固化的调试技巧5.1 多边界模型用 check_order 统一方向实际结构不止两个电极。每条边界先独立执行 check_order 并统一成同一个旋转方向再合并矩阵。判断内外边界的方法很简单取该条边界点的平均坐标作为形心计算边界切线与形心指向边界点的径向向量之间的叉积符号逆时针内边界的符号与逆时针外边界的符号正好相反根据目标方向逐条翻转。多边界合并后矩阵维数等于所有边界的匹配点总数分段组装时每段对应一个电极的电位值别在拼接处把行顺序弄混。5.2 用 fminbnd 自动搜索电荷位置比例lambda的搜索可以自动化。把主程序封装成一个输入lambda、返回检查区域最大相对误差的函数再用 fminbnd 扫描。function e csmError(lam) % 内部复用主程序逻辑组装 A、解 q、算检查点误差 % 需要把 R1, R2, V1, V2, N, 检查点全部做成可访问变量 % ... end lamBest fminbnd(csmError, 0.05, 0.8);fminbnd 基于黄金分割搜索和抛物线插值适合一维连续目标函数。这里它比优化工具箱里的 fmincon 更合适因为未知量只有一个函数调用次数少几秒内就能收敛。实际使用中建议先手动画一条lambda从 0.05 到 0.8 的误差曲线确认目标函数没有多个局部极小点再设搜索区间否则可能收敛到次优解。5.3 检查点布局决定误差结论的可信度检查点至少放两处一处电极之间的中场区反映整体趋势一处靠近高场强但不在匹配点上的位置暴露局部拟合缺陷。只在中场放检查点会把边界附近的失真全部漏掉只在近场放检查点又会把局部噪声误判成整体误差。两组检查点的最大相对误差分开记录调参时两者都要降下来才算真正收敛。把代码里的解析解替换成实验测量值或有限元参考解之后这套 check_order 加模拟电荷法的工作流就能直接迁移到任意二维剖面遇到边界顺序乱掉的算例第一步跑 check_order第二步再动参数。本文还有配套的精品资源点击获取