做磁路设计的时候最常被问到的问题之一就是两块长圆柱永磁体端面相对中间那一点气隙里的磁场到底怎么算。手边没有大型有限元软件时我一般会直接用Matlab加磁荷模型顶上——把永磁体表面的磁荷离散成一个个“单极”再逐点叠加这就是单极表面电荷密度方法。这篇文章拿它和更粗糙的点磁单极近似方法做对比给出完整推导、可运行的Matlab代码以及一组“哪些地方能用近似、哪些地方必须老老实实积分”的实用结论。这个方法特别适合两类人一是做永磁体结构方案预研的工程师手头只需要快速评估气隙磁场够不够强二是学电磁场数值计算的学生想找一个比有限元更直观、能自己在Matlab里写出来的入门案例。文章不用任何商业软件所有结果单靠一个脚本就能复现。1. 两块长圆柱永磁体之间的气隙场到底长什么样1.1 建模假设同向磁化、端面相对、只算空气区先把模型说清楚。我们研究的是两根轴线都在同一直线上的圆柱永磁体半径都取R 10 mm长度都取L 60 mm端面之间留一个气隙g。默认取g 2 mm后面会专门讨论g变化时两种方法的误差走向。两个磁体的磁化方向都沿轴向向上。这样布置之后下圆柱的上端面是N极上圆柱的下端面是S极两个极尖面对面气隙里的磁力线从下往上走。这个构型其实就是磁路中最常见的“N-S极对”在磁分离、磁夹具、粒子聚焦装置里都经常遇到。我们只关心气隙空气区域里的磁感应强度B不计算磁体内部的场。空气区没有磁化电流也没有自由电流所以磁场完全可以由磁荷模型给出。这也是磁荷法最大的优势不需要像有限元那样把整个求解域都画上网格只需要在磁体表面布置源然后空气中任意一点的场都可以用叠加积分直接算出来。1.2 为什么磁荷法比有限元更适合快速估算有人可能会问现在有限元软件已经很成熟了为什么要自己写磁荷法我的体会是方案阶段做参数扫描的时候有限元建模的时间成本太高了。改一次气隙尺寸、换一个磁体半径都要重新画网格、重新求解一次几秒到几分钟不等。磁荷法则是一劳永逸只要把几何参数写进脚本气隙从0.5 mm扫到50 mm所有场点作为独立变量并行算几百个工况也就是一瞬间的事。尤其在轴对称结构里表面离散可以用解析面积公式直接做不需要什么网格剖分库。当然有限元能处理非线性材料、复杂边界、涡流等问题磁荷法做不了。但单就“两块永磁体之间的空气磁场”这种静磁问题磁荷法的精度和速度都足够而且能清晰看到每个表面磁荷对总场的贡献脑子里特别有数。2. 单极表面电荷密度法的推导从面积分到面元求和2.1 均匀磁化圆柱为什么只有端面有磁荷磁荷模型的出发点是永磁体可以看成内部有磁化强度M的介质空间磁场由体磁荷密度ρ_m -∇·M和面磁荷密度σ_m M·n共同产生n是磁体表面的外法向单位向量。如果磁化是均匀的体内每一处M都是常矢量散度为零所以体磁荷为零。再看侧面轴向磁化时M与侧面法向垂直内积为零侧面磁荷也等于零。最后只剩下两个端面磁化向量穿出的那个端面带正磁荷磁化向量穿入的哪个端面带负磁荷。所以“长且均匀磁化圆柱”这个前提条件很有价值。它把一个三维体积分问题直接退化成了二维面积分问题而且是两个端面的面积分。如果磁化不均匀或者圆柱是倾斜磁化、局部退磁的那情况就会复杂很多后面我们在扩展部分再讲。2.2 面元源强公式与库仑式叠加为了避开单位上的混乱我这里直接引入“磁通源强”的概念。对于一个微元面积dS面磁荷产生的磁通源强为dΦ B_r · dS方向由N/S极性决定。N极取正S极取负。其中B_r是永磁体剩磁例如钕铁硼N42一般取1.2 T。为什么可以用B_r而不用磁化强度M因为标准磁荷公式给出B (μ0 / 4π) ∫ σ_m (r - r) / |r - r|³ dS而μ0 σ_m μ0 (M·n) B_r所以公式可以直接改写为B (1 / 4π) ∫ dΦ (r - r) / |r - r|³这就是库仑定律的磁学版本。每个磁荷微元都像电荷一样向外发射场场强与距离平方成反比方向沿着源到场点的连线。把面积分离散化我们就得到单极表面电荷密度方法的数值形式把端面划分成很多小面元每个小面元中心放一个“元磁单极”源强是ΔΦ B_r ΔS然后将所有面元的贡献叠加起来。面元越密结果越接近连续积分解。2.3 圆盘轴线的解析解给数值结果一把尺子数值方法写出来之后必须有一把尺子来验证它对不对。对于均匀面磁荷圆盘轴线上任意一点的场其实有严格的解析解。考虑一个半径为R的面磁荷圆盘面密度为B_r我们计算轴线上距盘面距离为z处的轴向场强。取圆盘上一小环带2πρdρ该环带对轴线上点的距离是√(ρ² z²)沿轴向的分量比例是z / √(ρ² z²)于是B_z (B_r / 2) [1 - z / √(z² R²)]这个式子在后面非常有用。对于两个端面相对的气隙轴线上任意一点y处的总场就是下极面正磁荷和上极面负磁荷分别贡献的叠加B_axis(y) (B_r/2)[1 - z₁/√(z₁²R²)] (B_r/2)[1 - z₂/√(z₂²R²)]其中z₁ y g/2是到下极面的距离z₂ g/2 - y是到上极面的距离。两个方向都指向上都是正贡献。这个解析式可以用来直接检查离散代码。只要离散面元足够多数值结果应该逐步逼近这个解析值。后面我在Matlab里也把这条解析曲线画出来做对比。3. 点磁单极近似图省事的方案边界在哪里3.1 把整个极面缩成一个点的数学过程点磁单极近似的思路非常朴素既然每个面元都按库仑定律贡献场那不如把整个极面的磁荷集中到圆心当成一个“超级磁单极”来算。单个极面的总源强是Φ B_r · πR²那么下极面的点单极位于y -g/2源强为Φ上极面的点单极位于y g/2源强为-Φ。气隙中任意点的场就是这两个点源的叠加B (Φ / 4π) (r - r₁)/|r - r₁|³ - (Φ / 4π) (r - r₂)/|r - r₂|³这个方法的优势是疯狂简单手算、Excel甚至心算都能出量级所以在工程估参时很受欢迎。但它牺牲了一个关键信息——极面的尺寸。3.2 用解析几何直接估算中心场点磁单极近似在轴线上的解析表达式也很漂亮。在气隙轴线上的y点下极面点源和上极面点源对轴向场强的贡献方向一致都是正方向于是B_mono_axis(y) (Φ / 4π) [1/(y g/2)² 1/(g/2 - y)²]特别注意这里是相加不是相减。因为正点源在下方场向上负点源在上方对下方场点产生的场也向上。很多人第一次写这个公式会写成相减得到中心场为零的可笑结果那就完全错了。在气隙中心y0处代入Φ B_r πR²可得B_mono_center 2 B_r R² / g²这个式子特别直观地暴露了点磁单极近似的毛病当气隙g远小于半径R时中心场按1/g²爆炸。物理上不可能。3.3 什么时候“省事”会翻车从解析式里已经能看出来气隙越小点磁单极近似越离谱。默认参数R 10 mm, g 2 mm时仅按两个点源叠加中心场算出来是60 T而按圆盘解析公式中心场只有约0.965 T。误差超过60倍。但反过来当气隙远大于半径时极面相对场点来说确实“小得像个点”点磁单极近似就会趋近于精确解。比如g 100 mm时点单极结果为0.024 T圆盘解析结果是0.0233 T误差只有3%左右。所以点磁单极近似不是不能用关键是得会判断边界。判据很简单看g/R的比值。我的经验是g/R 5时点单极近似可以作为快速估算g/R在2到5之间可以用于量级估计g/R 1时必须老老实实用表面积分否则结果会大到你怀疑自己单位写错了。4. 离散与数值细节这些坑会让结果完全不对4.1 圆盘面元的划分方式既然极面是圆盘最自然的离散方式就是极坐标网格半径方向分成Nr份圆周方向分成Nphi份每个小方格是一个环形扇区面元。第i个半径环的边界是r1 (i-1)/Nr · R和r2 i/Nr · R第j个角度扇区边界是φ1 (j-1)/Nphi · 2π和φ2 j/Nphi · 2π。该面元的面积是ΔS 0.5 (r2² - r1²) (φ2 - φ1)这个面积公式是精确的不是数值近似因为环形扇区的面积本来就能写成两个扇形面积之差。4.2 面元代表点取几何中心还是质量中心这里有个很容易被忽略、但直接影响收敛速度的细节每个面元里的“元磁单极”到底放在哪个位置。很多入门代码会直接取面元的几何中心即半径取(r1r2)/2角度取(φ1φ2)/2。但环形扇区的几何中心和整个面源分布的质量中心并不重合尤其是最内圈、靠近圆心的地方偏差很大。正确做法是取面元的质量中心rc 2/3 · (r2³ - r1³) / (r2² - r1²)这个公式是把面元上所有点的半径做加权平均得到的。当r10时rc 2R/3而不是R/2当r1接近r2时它退化到(r1r2)/2。使用质量中心之后同样的网格密度轴线场精度能明显提升一个档次。角度方向则直接取中值φ_m (φ1φ2)/2因为扇形在角度方向是对称的。4.3 网格密度与收敛性检查网格密度怎么选我的默认做法是先用一组比较粗的网格算一遍比如Nr20, Nphi40再用加密网格Nr80, Nphi160算一遍对比轴线场。如果两个网格结果在感兴趣区域里相对偏差小于0.1%就说明收敛了。对于默认参数R10 mm, g2 mmNr40, Nphi80已经足够把中心场精确到0.1%量级。但如果气隙压到0.5 mm场点就会离面元非常近这时候需要适当加密否则面元集总误差会让曲线出现锯齿状波动。4.4 远端磁荷要不要算标题里强调“长”圆柱这里终于派上用场了。每根圆柱其实有两个端面靠近气隙的极面和远离气隙的远端端面。我们一直说均匀磁化圆柱有两个端面带磁荷那远端端面对气隙里的场到底贡献多大计算一下就知道。下圆柱的远端端面位于气隙下方L处源强为负上圆柱的远端端面位于气隙上方L处源强为正。这两个远端点源在气隙中心产生的场方向都是向下也就是抵消主磁场。用圆盘公式估算当L 60 mm, g 2 mm时远端面对中心场的总贡献约-0.0159 T相对于主场的0.965 T大约占1.6%。所以“长圆柱”的物理意义就在这里只要圆柱长度远大于气隙远端磁荷就可以安全忽略。我的代码里留了一个includeFarEnds开关默认关闭。如果你的模型里圆柱不够长、或者你追求更高精度就打开它。5. Matlab代码可直接运行的单极法与点单极对比程序5.1 主脚本与参数定义下面这份代码是完整可运行的版本我把关键函数都写进去了。参数部分集中在最前面想改几何或材料直接动这里。%% 单极表面电荷密度方法 对比 点磁单极近似 % 模型两根轴向磁化圆柱永磁体端面相对气隙磁场计算 clear; clc; close all; % ---------------- 物理与几何参数 ---------------- Br 1.2; % 剩磁 [T]钕铁硼典型值 R 10e-3; % 圆柱半径 [m] L 60e-3; % 圆柱长度 [m]长圆柱假设远端可忽略 g 2e-3; % 气隙长度 [m] includeFarEnds false; % 是否计入远端端面磁荷false按长圆柱简化 % ---------------- 极面定义 ---------------- % 下圆柱上端面: y -g/2, sigma Br (N极) % 上圆柱下端面: y g/2, sigma -Br (S极) plates(1).y0 -g/2; plates(1).sigma Br; plates(2).y0 g/2; plates(2).sigma -Br; if includeFarEnds % 下圆柱下端面: y -(Lg/2), sigma -Br % 上圆柱上端面: y (Lg/2), sigma Br plates(3).y0 -(L g/2); plates(3).sigma -Br; plates(4).y0 (L g/2); plates(4).sigma Br; end % ---------------- 轴线场点 ---------------- % 略避开极面避免场点距面元过近造成波动 yAxis linspace(-g/2 1e-5, g/2 - 1e-5, 201); fieldAxis [zeros(size(yAxis)), yAxis, zeros(size(yAxis))]; % 单极表面电荷法 [~, BySurf] surfaceChargeField(plates, R, fieldAxis); % 圆盘解析解只含两个近端面 ByAna analyticTwoDisks(Br, R, g, yAxis); % 点磁单极近似 Phi Br * pi * R^2; % 单极源强 [Wb] BmonoAxis monopoleField([0, -g/2, 0], Phi, fieldAxis) ... monopoleField([0, g/2, 0], -Phi, fieldAxis); ByMono BmonoAxis(:,2);5.2 面元求和函数 surfaceChargeField核心函数就是做扇形面元离散并叠加库仑场。注意我专门用了质量中心半径rc而不是平均半径。function [Bx, By] surfaceChargeField(plates, R, points) % 单极表面电荷密度方法 % 输入plates为结构数组含y0和sigmaR为圆盘半径points为Nx3场点 % 输出Bx, By本文模型轴向对称z方向分量理论为零不单独返回 Nr 40; % 半径离散数 Nphi 80; % 角度离散数 N size(points, 1); Btot zeros(N, 3); for k 1:numel(plates) y0 plates(k).y0; sigma plates(k).sigma; for ir 1:Nr r1 (ir-1)/Nr * R; r2 ir/Nr * R; % 扇形面元面积 dA 0.5 * (r2^2 - r1^2) * (2*pi/Nphi); % 面元质量中心半径比几何中心收敛更快 rc 2/3 * (r2^3 - r1^3) / (r2^2 - r1^2); dPhi sigma * dA; % 面元磁通源强 [Wb] for ip 1:Nphi phiM (ip - 0.5) * (2*pi/Nphi); src [rc*cos(phiM), y0, rc*sin(phiM)]; % 源点 rvec points - src; % 源指向场点 r3 sqrt(sum(rvec.^2, 2)).^3; Btot Btot (dPhi/(4*pi)) .* rvec ./ r3; end end end Bx Btot(:,1); By Btot(:,2); end这段循环在Nr40, Nphi80时单个极面有3200个面元两个极面就是6400次迭代跑轴线200个点没有任何压力。要进一步提高速度的话可以把角度循环向量化但为了直观我没有这么做。5.3 点单极函数与解析校验点单极函数和圆盘解析式都非常短function B monopoleField(src, q, points) % 点磁单极库仑场 % src: 1x3 源点坐标 % q : 磁通源强 [Wb] rvec points - src; % Nx3 r3 sqrt(sum(rvec.^2, 2)).^3; B (q/(4*pi)) .* rvec ./ r3; endfunction By analyticTwoDisks(Br, R, g, y) % 两个均匀面磁荷圆盘在气隙轴线上产生的轴向场解析解 % 仅考虑两个近端面用于校验离散代码 z1 y g/2; % 到下极面距离 z2 g/2 - y; % 到上极面距离 By (Br/2) .* (1 - z1./sqrt(z1.^2 R^2)) ... (Br/2) .* (1 - z2./sqrt(z2.^2 R^2)); end主脚本里可以顺手做一次校验% 数值解与解析解对比 disp(数值解与圆盘解析解的最大绝对误差:); disp(max(abs(BySurf - ByAna)));在默认网格下这个误差应该在1e-5 T量级甚至更小。如果数量级对不上多半是面元中心取错或者单位混了。5.4 二维场分布绘图只算轴线太不过瘾我把气隙区域的二维轴对称场也画出来。这里需要注意二维图里横轴代表轴向y纵轴代表径向r我们的场点取在xr, z0这个截面上。%% 二维场分布径向-轴向切片 rAxis linspace(0, 2*R, 41); yAxis2 linspace(-1.5*g, 1.5*g, 51); [YY, RR] meshgrid(yAxis2, rAxis); pts2d [RR(:), YY(:), zeros(numel(RR),1)]; [Bx2, By2] surfaceChargeField(plates, R, pts2d); Bx2 reshape(Bx2, size(YY)); By2 reshape(By2, size(YY)); Bmag sqrt(Bx2.^2 By2.^2); figure; contourf(YY*1e3, RR*1e3, Bmag, 30, LineStyle, none); colorbar; hold on; quiver(YY*1e3, RR*1e3, By2, Bx2, k, AutoScaleFactor, 0.8); xlabel(轴向 y (mm)); ylabel(径向 r (mm)); title(气隙区域的磁感应强度幅值与方向); axis equal; xlim([-1.5*g, 1.5*g]*1e3); ylim([0, 2*R]*1e3);画出来之后你会在气隙中心看到一片颜色较亮、箭头密集且方向几乎一致的区域这就是有效磁场区。越靠近极面边缘箭头会出现明显弯曲那就是边缘效应。6. 结果对比轴线分布、误差曲线、气隙尺寸影响6.1 小气隙下两种方法的巨大分歧先用默认参数R 10 mm, g 2 mm跑一遍。轴线上从下极面到上极面取201个点单极表面电荷法的结果和圆盘解析解几乎重合这验证了代码本身是可靠的。点磁单极近似在气隙中心给出的场强是60 T而单极表面电荷法只有0.965 T。这个60 T不是物理真实值而是“把每个极面压成一个点”这个近似在近场区域彻底失效的体现。物理上两块钕铁硼永磁体之间不可能出现60 T的磁场这个数字本身就是警报g远小于R时不能用点单极近似。6.2 误差随位置的变化轴线上的误差不是均匀分布的。把两种方法做相对误差relErr abs(ByMono - BySurf) ./ abs(BySurf) * 100;你会发现一个“两头高、中间相对低”的形状。气隙中心处点单极近似误差虽然也大到百分之几千但真正爆炸的是极面附近——距离极面不到0.1 mm的位置点单极近似的场随1/d²增长而真实的面磁荷分布的场增长率要慢得多两者能差几个数量级。这也是为什么我在第4节反复强调“场点不能离面元太近”。即使对单极表面电荷法自己来说当面元尺寸相对气隙不够小时离散误差也会在极面附近被放大。6.3 G/R 比值决定能否用点单极用中心场作为评价指标把两种方法的解析值做成表格趋势一目了然气隙 g (mm)g/R面电荷法中心场 (T)点单极法中心场 (T)相对误差0.50.051.170960820倍20.20.9656062倍1010.6632.4262%3030.2020.26732%100100.02330.02403%表格里“面电荷法”一列用的是圆盘解析公式也就是单极表面电荷法网格无限密时的极限值。可以看到当g/R从0.2增长到10点单极近似的误差从62倍快速收敛到3%左右。工程上的经验阈值我一般取g/R 5。小于这个值就老老实实跑表面积分代码反正脚本也就几十毫秒的事没必要为省这一点计算量引入不可控误差。6.4 二维场的形态二维场图能直观看出气隙磁场的分布特征。中心区域磁力线几乎平行轴向场占主导这正是设计磁路时最希望看到的“均匀工作区”。在极面边缘处磁力线向外拱起径向分量变大形成典型的边缘磁场。如果气隙进一步缩小比如g0.5 mm你会发现中心区域的均匀性反而更好边界效应集中在更窄的环带里。这就是为什么很多磁路设计会把气隙尽量做小来换取更高的中心场强和更好的均匀度。代价是机械安装精度要求更高。7. 实操经验单位、极性校验、扩展思路7.1 最容易错的符号与单位这个模型我前前后后帮人调试过很多次出现最多的错误有两个。第一个错误是正负号。下极面N为正磁荷上极面S为负磁荷。在写叠加循环时一旦把sigma的符号搞反轴线中心的场会变成两个贡献相减结果小得离谱。最稳妥的自检方法就是验算中心场气隙中心处磁场方向必须一致向上而且数值接近Br/2的量级。如果方向或量级不对先查符号。第二个错误是单位。B_r用特斯拉长度用米面元源强ΔΦ的单位就是韦伯Wb公式比例系数是1/(4π)。有人喜欢保留μ0和磁化强度M的写法那也行但很容易把M B_r/μ0代进去又乘以μ0结果偏大几个数量级。我的建议是代码里统一用B_r作为面磁荷源强不要出现μ0省心。7.2 场点在极面附近怎么处理数值计算时场点不要恰好落在极面上也不要落在极面内部因为那里是磁荷分布所在的位置库仑公式不适用。轴线场点通常从-g/2 δ取到g/2 - δδ取气隙的1/100到1/1000都行。如果确实需要计算极面表面或极面外很近处的场需要更精细的处理。简单做法是把网格细分到场点距离的十分之一以下复杂做法是提取奇异积分的解析部分。对大多数气隙应用来说避开极面是明智的选择毕竟你真正关心的是空气隙里的场。7.3 向锥形极尖、Halbach 结构扩展这套方法最大的价值是可以快速迁移到其他轴对称结构。比如你把圆柱端面加工成锥形极尖只需要改plates的几何描述把每个面元所在半径和轴向位置都作为半径的函数来离散原来的一维y0变成随半径变化的数组叠加逻辑完全不用动。又比如你想算Halbach阵列外侧的气隙场那就需要把每根圆柱的磁化方向改成非轴向这时候侧面磁荷不再为零需要把侧面也加入离散。原理仍然是同一个积分公式只是需要按实际几何把表面分成若干段。我后来做磁分离器方案的时候就是在这个代码基础上加了锥形极尖的几何参数半天时间就把不同极尖角度的气隙场扫了一遍出图、出曲线、出结论速度比开有限元软件快得多。这个方法本身不复杂但把推导过程吃透之后改起来非常顺手。