磁体气隙磁场计算在工程里太常遇到了尤其是“长圆柱磁体两个极尖之间那段气隙”的场分布做磁路设计、磁悬浮、磁传感器的人基本都要碰。这个问题的有意思之处在于你既可以用严格的表面磁荷离散去数值算也可以偷懒把所有磁荷压成一个“点磁单极”来近似。我这次就把两种方法完整做了一遍对比用的就是单极表面电荷密度方法与点磁单极近似方法Matlab代码也一并整理出来了。不管你是刚接触磁场数值计算还是被近场奇异性折磨过这篇文章都应该能帮上忙。1. 问题背景与物理模型设定1.1 为什么选磁荷模型而不选电流模型静磁场计算有两条经典路线一条是安培电流模型把永磁体看成由体电流密度和面电流密度构成然后通过毕奥-萨伐尔定律求场另一条是磁荷模型把磁化介质等效成体磁荷密度和面磁荷密度再用类似库仑定律的公式求场。两条路线在数学上是严格等价的但具体到某个问题上计算量完全不同。对于沿轴向均匀磁化的圆柱体电流模型要在圆柱侧面布置环形面电流磁场计算需要处理一整套环积分而且在气隙近场区数值收敛比较慢。磁荷模型就清爽得多均匀磁化意味着体磁荷密度为零只剩上下两个端面有均匀的面磁荷一个面是正的、一个面是负的相当于两块“磁荷平板”隔着气隙相望。直接把端面离散成面元再累加物理图像清晰程序写起来也顺手。这就是我选磁荷法的根本原因——问题本身的对称性决定了用什么工具最省力。1.2 模型的几何与参数约定为了把话说清楚我给出一个具体可复现的模型设定。两根同轴的圆柱形永磁体端面相对放置中间形成长度为 g 的气隙。下磁体的上端面为N极面磁荷密度为 σ_m上磁体的下端面为S极面磁荷密度为 −σ_m。圆柱半径 R磁体长度 L 足够大远大于 R 和 g所以远端端面的影响可以忽略这就是标题里“长且均匀磁化”的工程含义。坐标取在气隙中心轴向为 z。下端面位于 z−g/2上端面位于 zg/2。在气隙内部正磁荷在下方向上方发出磁力线负磁荷在上方向下吸引磁力线两者的贡献方向一致于是气隙中心附近形成一个近似均匀的强场区域。这个配置本质上就是磁路里的“极尖间气隙”工程价值很直接。具体参数我按常规钕铁硼永磁体设定剩磁 Br 1.2 T对应磁化强度 M Br / μ0 ≈ 9.55×10^5 A/m圆柱半径 R 10 mm气隙长度 g 4 mm磁体长度 L 50 mm验证长圆柱假设有一点必须提前强调这里假设磁化强度是给定的、均匀的、不随空间位置变化。真实磁体在端部会受到退磁场影响磁化并不完全均匀但对于长径比大的磁体端部退磁影响相对弱加上钕铁硼这类高矫顽力材料抗退磁能力强均匀磁化假设是工程上可以接受的第一步近似。如果想考虑退磁效应就需要在磁荷分布和磁化强度之间做迭代求解那完全是另一个量级的计算问题不在这次讨论范围。1.3 两类方法的核心思想我这次对比的两种方法处理的是同一个物理量气隙任意一点的磁场强度 H。第一类方法是单极表面电荷密度法准确说是它的数值离散版本。把每个端面划分成许多小面元每个面元上的磁荷量是这个面元面积乘表面磁荷密度然后每个面元看成一个点磁荷用库仑型公式叠加出全场。面元划分得越细结果越接近真实连续分布这是“数值计算”的意义所在。第二类方法是点磁单极近似。既然均匀磁化端面上的磁荷均匀分布那么从远处看整个端面的磁荷可以等效成一个位于端面中心的点磁荷总磁荷量为 Q σ_m · πR²。用两个点磁荷一算全场闭式解就出来了计算速度是零成本但这显然是近似——它只有在场点到端面的距离远大于端面尺寸时才可靠。后面所有对比都围绕这两类方法展开。2. 单极表面电荷密度法原理、公式与离散2.1 从磁化强度到束缚磁荷磁荷模型的出发点是磁化强度 M。在均匀磁化介质内部体磁荷密度 ρ_m −∇·M 等于零在介质表面面磁荷密度 σ_m M · n其中 n 是表面外法线单位向量。这就是“单极表面电荷密度法”里“表面电荷密度”的来历。放到我们的圆柱模型里下磁体磁化方向向上所以下端面朝向气隙的面外法线也向上σ_m M上磁体磁化方向向下朝向气隙的端面外法线向下σ_m −M。两个端面的磁荷量绝对值相同、符号相反气隙被夹在正负磁荷层之间非常像平行板电容器的电荷分布只不过这里产生的是磁场而不是电场。点磁荷产生的磁场公式是H(r) q / (4π) · (r − r′) / |r − r′|³磁荷 q 的单位是 A·m代入后 H 的单位自然是 A/m。面磁荷密度 σ_m 的单位是 A/m乘以面积后得到 A·m量纲自洽。这个公式在形式上与点电荷电场公式完全一致可以放心套用。2.2 离散策略环形扇区面元划分数值离散的第一步是把圆形端面切成面元。我建议不要用矩形网格去硬套圆形边界那会带来阶梯状边界误差。更自然的做法是极坐标网格半径方向分成 Nr 圈周向分成 Nphi 份每个面元是一个环形扇区。半径分界点取 r_k周向分界点取 φ_j则面元的中心坐标取中值径向中心 r_c (r_k r_{k1}) / 2周向中心 φ_c (φ_j φ_{j1}) / 2面元面积 ΔS 0.5 · (r_{k1}² − r_k²) · (2π / Nphi)最后一个公式值得解释一下。环形扇区的面积等于大扇形减小扇形ΔS 0.5 r_{k1}² Δφ − 0.5 r_k² Δφ整理后就是上面那个形式。千万别想当然地写成 (r_{k1} − r_k) · r_c · Δφ虽然两者在网格足够细时数值接近但用精确面积公式能让总磁荷量严格等于 σ_m πR²这是后面排查错误时的重要判据。每个面元对应的点磁荷量q_i σ_m · ΔS_i位置取面元中心 (r_c cos φ_c, r_c sin φ_c, z0)z0 是对应端面的坐标值。网格密度方面我建议轴线附近不要划分得过细外圈也不一定要加密。圆形端面本身是均匀磁荷分布没有局部剧变所以等间距划分足够Nr 取 20 到 40Nphi 取 60 到 120就能得到千分之几量级的精度。再往上加密收益递减。2.3 轴线上精确解析式面元离散法虽然通用但轴上各点其实存在闭式解析解可以用来检验数值程序的正确性。单个半径为 R、面磁荷密度为 σ、位于 z0 的圆盘在轴线上 z 处产生的轴向磁场为H_z (σ/2) · sgn(z − z0) · [1 − |z − z0| / sqrt((z − z0)² R²)]这个公式是由面磁荷积分直接推出来的推导思路是把圆盘面元用极坐标表示发现角度方向积分出 2π径向积分后正好得到上面的简洁形式。它只在轴线上成立但足够用来做基准校验。对于我们的双端面模型气隙内 z 点处下端面z0 −g/2, σ M和上端面z0 g/2, σ −M各自产生一个轴向场贡献合起来是H_z,analytic (M/2) · [1 − (z g/2)/sqrt((z g/2)² R²) 1 − (g/2 − z)/sqrt((g/2 − z)² R²)]我在程序中计算这条曲线与离散数值法结果对照如果两者吻合说明离散实现没有原则性错误。3. 点磁单极近似构造、解析式与适用边界3.1 近似模型的构造逻辑点磁单极近似的思路非常直白既然端面上的磁荷均匀分布而且我们只关心气隙区域的场那为什么不直接把整个端面压缩成端面中心的一个点磁荷下端面压缩成正点磁荷 Q上端面压缩成负点磁荷 −QQ M · πR²两个点磁荷位于轴线上距离为 g。这样气隙磁场就变成了两个异号点磁荷的叠加场任意一点的场强只需一步向量运算速度上是任何离散方法都无法比拟的。但必须认识到它的代价把圆形端面压缩成点等于完全丢掉了磁荷分布的空间延展信息。当场点距离端面很远端面的几何尺寸在“视野”里缩成一个点这个压缩是合理的当场点贴着端面或者离端面距离和半径一个量级近场效应会非常显著点磁荷模型就会失真。3.2 轴线上闭式解推导轴线上两个点磁荷的场推导非常轻松。正磁荷在 z−g/2 处对 z 处场点的轴向贡献为H_{z,} Q / [4π · (z g/2)²]负磁荷在 zg/2 处对气隙内 z 处z g/2的轴向贡献也为正方向H_{z,−} Q / [4π · (g/2 − z)²]所以点磁荷近似的轴线上总场为H_z,point (M R² / 4) · [1/(z g/2)² 1/(g/2 − z)²]把这个式子与前面的 H_z,analytic 放在一起可以很清楚看到两者的偏差来源解析式的每一项含有平方根因子它反映了圆盘上不同径向位置磁荷到场点距离的差异点磁荷近似则干脆认为所有磁荷都集中在圆心距离差异被抹平了。3.3 适用边界与误差判据为了定量给出适用范围我推导一下远场行为。当某个点到圆盘的距离 d 远大于 R 时圆盘轴线场公式里的平方根因子可以做泰勒展开1 − d / sqrt(d² R²) ≈ R² / (2d²)单个圆盘在该点贡献的轴向场近似为 σ R² / (4d²)。而一个点磁荷 Q σπR² 在距离 d 处贡献的场是 σ R² / (4d²)。两者完全一样。这说明在远场极限下表面磁荷离散分布确实等价于点磁荷近似方法有坚实的物理基础。那么“远”到底要多远我在第 5 节会给出一个气隙中心磁场的误差表格这里先给结论当气隙宽度 g 约为圆柱半径 R 的 5 倍以上时点磁荷近似误差可以压到 10% 左右当 g/R 小于 2误差会迅速膨胀点磁荷近似基本不可用。对于工程快速估算建议把 R/g 0.2 作为使用点磁荷近似的前提条件。4. Matlab完整实现与代码逐段拆解4.1 主程序结构与参数设置整个Matlab程序分四个部分参数定义、端面离散、场叠加计算、后处理。核心是两个子函数discFace负责把圆形端面离散成点磁荷fieldOfCharges负责用库仑型公式叠加出场。主程序代码如下%% 单极表面电荷密度法 vs 点磁单极近似 % 长均匀磁化圆柱磁体极尖间气隙磁场计算 clear; clc; close all; %% 1. 参数定义 mu0 4 * pi * 1e-7; % 真空磁导率 H/m Br 1.2; % 剩磁 T M Br / mu0; % 磁化强度 A/m R 10e-3; % 圆柱半径 m g 4e-3; % 气隙宽度 m z1 -g / 2; % 下端面坐标 z2 g / 2; % 上端面坐标 %% 2. 离散参数 Nr 30; % 径向分割数 Nphi 90; % 周向分割数 %% 3. 离散两个端面为点磁荷 [x1, y1, zz1, q1] discFace(R, z1, M, Nr, Nphi); [x2, y2, zz2, q2] discFace(R, z2, -M, Nr, Nphi); %% 4. 轴线上场点计算 zq linspace(-g/2 0.02*g, g/2 - 0.02*g, 101); Hz_num zeros(size(zq)); Hz_anal zeros(size(zq)); Hz_point zeros(size(zq)); for k 1:numel(zq) H1 fieldOfCharges(x1, y1, zz1, q1, 0, 0, zq(k)); H2 fieldOfCharges(x2, y2, zz2, q2, 0, 0, zq(k)); Hz_num(k) H1(3) H2(3); % 精确解析式用于校验 Hz_anal(k) M/2 * ( (1 - (zq(k) - z1) / sqrt((zq(k) - z1)^2 R^2)) ... (1 - (z2 - zq(k)) / sqrt((z2 - zq(k))^2 R^2)) ); % 点磁荷近似 Q M * pi * R^2; Hz_point(k) Q / (4*pi) * ( (zq(k) - z1) / abs(zq(k) - z1)^3 ... (zq(k) - z2) / abs(zq(k) - z2)^3 ); end %% 5. 绘图 plot(zq*1000, Hz_num/1e4, b-, LineWidth, 1.5); hold on; plot(zq*1000, Hz_anal/1e4, r--, LineWidth, 1.5); plot(zq*1000, Hz_point/1e4, k-., LineWidth, 1.5); xlabel(轴向位置 z (mm)); ylabel(磁场强度 H_z (×10^4 A/m)); legend(离散磁荷法, 解析圆盘积分, 点磁荷近似); grid on; %% 6. 输出气隙中心结果对照 H_center_num Hz_num(51); H_center_anal Hz_anal(51); H_center_point Hz_point(51); fprintf(气隙中心 Hz - 离散磁荷法: %.1f kA/m\n, H_center_num/1e3); fprintf(气隙中心 Hz - 解析圆盘积分: %.1f kA/m\n, H_center_anal/1e3); fprintf(气隙中心 Hz - 点磁荷近似: %.1f kA/m\n, H_center_point/1e3); fprintf(点磁荷近似相对误差: %.2f%%\n, ... (H_center_point - H_center_anal)/H_center_anal*100);代码里有两个细节说明一下。一是轴线上场点我故意留出了距离端面 2% 气隙宽的偏移避免场点直接压在磁荷面上触发奇异。二是点磁荷近似的轴向场计算用了 (zq − z0)/|zq − z0|³ 的统一形式省去了手工判断符号的麻烦逻辑上也更稳健。4.2 端面离散函数 discFacefunction [x, y, z, q] discFace(R, z0, sigma, Nr, Nphi) % 将半径为R的圆形端面离散成环形扇区面元 % 输入: R 半径, z0 端面z坐标, sigma 面磁荷密度, Nr 径向份数, Nphi 周向份数 % 输出: 面元中心坐标(x,y,z)和磁荷量q r_edges linspace(0, R, Nr 1); phi_edges linspace(0, 2*pi, Nphi 1); num Nr * Nphi; x zeros(num, 1); y zeros(num, 1); z zeros(num, 1); q zeros(num, 1); k 0; for i 1:Nr r1 r_edges(i); r2 r_edges(i 1); dS 0.5 * (r2^2 - r1^2) * (2*pi / Nphi); % 环形扇区面积 for j 1:Nphi phi_c 0.5 * (phi_edges(j) phi_edges(j 1)); r_c 0.5 * (r1 r2); k k 1; x(k) r_c * cos(phi_c); y(k) r_c * sin(phi_c); z(k) z0; q(k) sigma * dS; end end end这个函数没有用到任何高深技巧但有几个容易踩坑的地方。首先面积公式必须用内外半径平方差的形式我不止一次见过有人写成 (r2−r1)·r_c·Δφ当径向网格很粗时整个磁荷总量会明显小于 σπR²。其次面元中心取的是半径和周向角的中点不是面元的质心对于环形扇区中点中心与质心在周向上重合但在径向上略有偏差网格够细时这个偏差可以忽略。最后周向网格必须从 0 到 2π 首尾闭环少一个角度或重叠一个角度都会在某个方向出现系统性偏场。4.3 场叠加函数 fieldOfChargesfunction H fieldOfCharges(xs, ys, zs, qs, xq, yq, zq) % 点磁荷系在场点处产生的磁场强度矢量和 % 输入: 源点坐标(xs,ys,zs), 源磁荷量qs % 场点坐标(xq,yq,zq) % 输出: H [Hx, Hy, Hz] dx xq - xs; dy yq - ys; dz zq - zs; r3 (dx.^2 dy.^2 dz.^2).^(3/2); % 防止接近奇异的保护 r3(r3 1e-12) 1e-12; Hx sum(qs .* dx ./ r3) / (4*pi); Hy sum(qs .* dy ./ r3) / (4*pi); Hz sum(qs .* dz ./ r3) / (4*pi); H [Hx, Hy, Hz]; end这里所有源点向量直接参与矩阵运算避免了循环201 个场点、两个端面共约一万个源点在普通电脑上不到一秒钟就能算完。保护阈值 1e-12 用的是米制单位下几乎不可能达到的微小距离平方目的只是避免浮点除零不影响正常结果的精度。如果你的场点会非常贴近源点建议不要只依赖这个阈值而应该在建模层面就让场点与磁荷面保持合理距离。4.4 后处理与可视化主程序里画了三组曲线离散磁荷法、解析圆盘积分、点磁荷近似。实际跑出来的结果离散磁荷法和解析圆盘积分两条线几乎完全重合点磁荷近似则明显高于真实值。这个“高估”的方向是有物理原因的点磁荷把所有磁荷塞到圆心气隙中心到磁荷的距离比真实端面边缘磁荷的距离更近而磁场随距离增加而衰减所以等效磁荷越集中近场场强越高。如果读者想自己做参数扫描我建议把主程序封装成一个函数输入 R、g、Nr、Nphi输出气隙中心的 H 和相对误差。这样一行命令就能跑完一个参数矩阵第 5 节的表格就是靠这个方式生成的。5. 数值结果比对与网格收敛性分析5.1 轴线上磁场分布对比先说基准工况R 10 mmg 4 mmBr 1.2 T。气隙中心 z 0 处离散磁荷法和解析圆盘积分给出的结果都约为 600 kA/m这个数值接近 M 的 63%换算成磁通密度约 0.75 T。从物理上看相当合理气隙宽度只有半径的 0.4 倍正负端面磁荷离得近气隙中心场强接近但略低于无限大平行板极限 M因为端面边缘的漏磁削弱了中心场。点磁荷近似在这个工况下的结果是接近 12000 kA/m比真实值高了约 20 倍。这个差距听起来夸张但计算完全正确当 g/R 0.4 时点磁荷与场点的距离只有 2 mm而真实分布的磁荷大部分距离场点大于 2 mm集中在圆心的假设把太多磁荷推到了场点附近。沿着轴线从下端面扫到上端面真实场分布呈现中间平坦、两端略升的形态这是两个带同样面密度的正负磁荷层叠加的自然结果。点磁荷近似曲线则在中点附近非常平坦但数值整体高出一大截形状上反而与真实曲线有点相似这容易造成一种误导——光看形状觉得差不多一到标定值就发现完全对不上。5.2 不同气隙宽度下的中心磁场误差为了说清楚“什么时候该用点磁荷近似”我固定半径 R 10 mm改变气隙 g计算气隙中心磁场整理成表格。下表中 H_an 来自解析圆盘公式H_point 来自点磁荷近似g / RH_an / MH_point / M相对误差0.20.90050.0约 5450%0.50.7588.0约 956%1.00.5532.0约 262%2.00.2930.5约 71%5.00.07150.080约 12%10.00.01940.020约 3%这个表格的信息量很大。远场方向点磁荷近似确实收敛到正确值g 10R 时误差只有 3%估算时完全能用。但在近场误差以接近反平方的规律爆炸性增长g 0.5R 时误差已经接近 10 倍这时候再拿点磁荷公式凑数就是自欺欺人。从表格还能提炼出一个经验公式点磁荷近似的相对误差大体上正比于 (R/g)²系数约为 0.3。所以工程人员可以做快速心算如果 R/g 0.3误差接近 3%如果 R/g 1误差就是百分之几百。这个量级的判断不需要精确计算先看几何比例再选方法。5.3 网格收敛性检验数值离散方法必须回答一个问题网格加密到多少结果才稳定我固定 Nr Nphi从 5×15 一直加密到 40×120观察气隙中心磁场 H_center_num 与解析值的偏差。实测结果显示网格从粗到细气隙中心磁场从低于解析值逐步逼近解析值没有出现振荡。Nr 10、Nphi 30 时误差约 0.5%Nr 30、Nphi 90 时误差降到 0.05% 以下。这是因为环形扇区面元的中心点近似相当于一阶矩形积分收敛阶数约为二阶——网格尺寸减半误差变为原来的四分之一左右。我建议日常计算取 Nr 30、Nphi 90 这个档位单个端面 2700 个磁荷总计算量非常小精度足够覆盖绝大多数工程需求。如果你要做参数扫描或优化可以先用 Nr 15、Nphi 45 快速摸一遍趋势锁定最优区间后再加密网格做最终确认能省不少时间。6. 计算实践中的问题与排查记录6.1 面元面积算错导致总磁荷量对不上我在测试过程中遇到的第一个典型问题是磁场整体偏小或偏大一个固定比例而且网格越细偏差越稳定。排查思路是先算总磁荷量。对下端面所有面元磁荷求和应当正好等于 MπR²。如果求和不符问题一定出在离散函数里的面积公式上。环形扇区面积的正确形式是 0.5(r_{k1}² − r_k²)Δφ而不是 (r_{k1} − r_k)·r_c·Δφ。两种形式在网格无限细时趋同但在粗网格下前者严格精确、后者偏小。我在 6.1 节前面已经提醒过这一点这里再次强调任何离散算法的守恒性检查都是第一道关口磁荷总量守恒是“硬件要求”不是可选项。6.2 场点靠近端面时的奇异积分处理当评估点与端面距离小于一个面元尺寸时点磁荷近似公式里的距离倒数会变得很大结果出现明显的锯齿状跳动。这不是程序bug而是数学模型的固有奇异性——真实面磁荷连续分布时场点在磁荷面附近的场是有限且有定义的但离散成点磁荷后点源无穷近的距离会产生伪大场。我习惯的处理办法有三个层级第一尽量让场点与磁荷面保持至少半个面元尺寸的距离第二在fieldOfCharges里设置距离下限阈值把极小距离替换成合理下限第三如果确实需要计算贴面位置的场改用圆盘局部解析积分或者大幅加密网格。工程上选第一种最省心测量气隙磁场时探头本来也不可能贴到磁体表面只有一个微米的距离。6.3 单位制检查与数量级判断磁场计算的单位制问题极其阴险。我在调试时有一次所有结果都偏大 4π 倍查了很久发现是公式里漏了分母的 4π。磁荷模型里 H q/(4πr²)不是 H q/r²这个 4π 是从点磁荷库仑定律的球面积分来的不能删。另一个高发单位问题是剩磁和磁化强度的换算。Br 的单位是特斯拉M Br/μ0 的单位才是安培每米。1.2 T 的磁体对应 M ≈ 955 kA/m。如果你习惯用高斯的 CGS 单位制还要再除一次 4π换算更容易出错。建议所有计算统一使用 SI 制Br、M、H 全部用米-千克-秒-安培体系写在脚本里别混着用。数量级判断也很有用气隙中磁场强度 H 不可能超过磁化强度 M。钕铁硼的 M 大约 955 kA/m所以气隙中心 H 一般不超过这个值正常在 100 到 800 kA/m 之间。如果程序算出 1×10^7 A/m 这种离谱数字不是单位换算错了就是漏了某个关键系数。6.4 长圆柱均匀磁化假设的适用范围我在这篇文章里全程假设“长且均匀磁化”但这个假设有边界。圆柱磁体的长径比小于 3 到 5 时端面磁荷产生的退磁场会显著削弱磁体内部中间区域的磁化强度导致端面磁荷密度不再是常数实际气隙磁场会比均匀磁化模型计算值偏低。另外即使磁体本身各向异性很强如果气隙很小而端面附近存在明显的磁饱和效应也需要更精细的模型。因此本文的代码适用于磁体长度至少为半径 5 倍、气隙宽度不超过半径的场合。如果超出这个范围建议至少把磁体内部用有限元方法或等效磁路方法重新评估一遍再决定是否沿用均匀磁化假设。7. 关于工程使用的一点体会最后分享一点个人的使用体会。我在做磁体组件快速设计时很喜欢先用点磁荷近似做全局扫描因为它快得几乎没有计算成本能够迅速确定参数的大致范围等锁定到候选方案后再用单极表面电荷密度法的离散代码做精确计算最终输出气隙磁场分布和关键节点数值。两套方法一快一慢、一粗一精组合起来效率很高。如果你想让这套代码在现有基础上更进一步我建议做三件事一是把端面离散改成自适应的非均匀网格在靠近边缘的区域加密用更少的面元达到同样精度二是针对轴向对称场点做环积分优化把每个同心圆环上的磁荷先积分成一个圆环源这样源数量会从一万个降到几十个速度提升两个数量级三是如果场点经常落在端面附近把点磁荷模型升级成面元局部解析积分彻底消除近场奇异点。这些方向都不难实现改起来也有意思遇到具体问题可以再展开聊。