1. 为什么静磁场仿真绕不开解泊松方程这道坎拿到静磁场仿真-主题005_有限差分法求解泊松方程这个题目时我先说一个可能让新人意外的结论静磁场仿真里真正折磨人的部分往往不是模型怎么画、网格怎么剖而是求解器在后台默默算的那套偏微分方程。换句话说懂不懂泊松方程的数值解法直接决定了你是会点仿真软件还是真的懂仿真。静磁场的控制方程从麦克斯韦方程组出发在忽略位移电流、仅考虑恒定磁场的条件下会简化成两组旋度方程[ abla \times \mathbf{H} \mathbf{J}, \quad abla \cdot \mathbf{B} 0 ]其中 (\mathbf{B} \mu \mathbf{H})。把这两组方程联立引入矢量磁势 (\mathbf{A})满足 (\mathbf{B} abla \times \mathbf{A})并选定库仑规范 ( abla \cdot \mathbf{A} 0)最终会得到一个矢量形式的泊松方程[ abla^2 \mathbf{A} -\mu \mathbf{J} ]如果是二维问题——比如电机横截面、变压器窗口、电磁铁轴向剖面这类典型场景——矢量磁势只剩一个 z 方向分量 (A_z)方程就进一步退化成标量形式[ \frac{\partial^2 A_z}{\partial x^2} \frac{\partial^2 A_z}{\partial y^2} -\mu J_z ]这就是标准形式的二维泊松方程。换句话说静磁场仿真的核心任务本质上就是在给定电流源 (J_z)、给定材料磁导率 (\mu)、给定边界条件的前提下把 (A_z) 的分布求出来。有了 (A_z)磁通密度 (\mathbf{B}) 就能通过一阶导数直接得到。而有限差分法FDM是目前求解这类方程最经典、最容易理解的数值方法之一。它的核心思路可以用八个字概括化微分为差分化连续为离散——把求解区域划分成规则的网格用网格节点上的未知量近似连续函数用差分公式替代偏导数最终把偏微分方程变成一个大型线性代数方程组。解方程组就得到每个节点上的 (A_z) 值。这篇文章我会完整过一遍从二维静磁场控制方程的推导逻辑到五点差分格式的构造过程再到 Dirichlet 边界和 Neumann 边界的离散实现、超松弛迭代SOR加速收敛的做法以及我自己在实操中积累的收敛判定、松弛因子选取、失效排查方面的经验。适合刚接触电磁场数值计算的研究生、做电机或电磁器件设计的工程师也想把这个方法吃透的仿真爱好者。这不是一篇软件操作教程而是一篇把仿真器后台到底在算什么讲清楚的原理加实战笔记。2. 二维静磁场位函数方程从 nabla 算子到可离散的差分算子2.1 麦克斯韦方程组如何一步步退化成泊松方程很多教材一上来就丢出泊松方程导致初学者总觉得它是从天而降的。实际上整个过程是可以一步步推下来的而且每一步都有明确的物理含义。先从安培环路定律的微分形式说起。在静磁场条件下电流密度不随时间变化电场旋度为零因此磁场强度 (\mathbf{H}) 与电流密度 (\mathbf{J}) 满足[ abla \times \mathbf{H} \mathbf{J} ]同时磁通密度 (\mathbf{B}) 永远是无散场磁力线闭合没有磁单极子所以有[ abla \cdot \mathbf{B} 0 ]对于一个无散矢量场总可以把它写成某个矢量场的旋度。于是引入矢量磁势 (\mathbf{A})令[ \mathbf{B} abla \times \mathbf{A} ]接着利用本构关系 (\mathbf{B} \mu \mathbf{H})代入安培环路定律的微分形式[ abla \times \left( \frac{1}{\mu} abla \times \mathbf{A} \right) \mathbf{J} ]如果材料是线性的、各向同性的且磁导率 (\mu) 在求解区域内为常数或者分片常数那么上式可以简化。利用矢量恒等式[ abla \times ( abla \times \mathbf{A}) abla ( abla \cdot \mathbf{A}) - abla^2 \mathbf{A} ]再选定库仑规范 ( abla \cdot \mathbf{A} 0) 消掉第一项就得到[ abla^2 \mathbf{A} -\mu \mathbf{J} ]这就是静磁场的矢量泊松方程。在二维问题中电流方向只有 z 方向即 (\mathbf{J} (0, 0, J_z))因此矢量磁势也只剩下 z 分量 (\mathbf{A} (0, 0, A_z))拉普拉斯算子作用在矢量上就退化为作用在标量上[ \frac{\partial^2 A_z}{\partial x^2} \frac{\partial^2 A_z}{\partial y^2} -\mu J_z ]这才到了我们真正要解的标量泊松方程。2.2 为什么在无电流区域方程会变成拉普拉斯方程这里有一个非常关键且容易忽略的细节并不是整个求解区域都有电流。比如电机气隙、铁芯、永磁体以外的空气区域这些地方 (J_z 0)方程自动退化成了拉普拉斯方程[ \frac{\partial^2 A_z}{\partial x^2} \frac{\partial^2 A_z}{\partial y^2} 0 ]这个细节在编程实现时非常重要因为你在组装系数矩阵的时候需要根据每个网格节点是否位于电流区域来决定源项是否为 0。很多新手跑出来的结果一团糟排查半天发现哦漏了判断电流区域这个坑我后面还会专门提。泊松方程统一 源项置零的处理方式比分别实现两套方程要稳得多。你在写代码时只需要维护一个 (J_z) 数组有电流的节点填真实值无电流的节点填零然后用统一的离散格式求解即可。2.3 从偏微分方程到代数方程那一步关键的离散化连续形式下(A_z(x, y)) 是定义在求解区域内每一个点上的二元函数有无限多个自由度没法直接丢给计算机。有限差分法的做法是在求解区域上画一张规则的网格比如均匀矩形网格把连续函数的值只保留在网格的节点上。如果沿着 x 方向以步长 (h_x) 划分沿着 y 方向以步长 (h_y) 划分那么节点坐标就是 ((x_i, y_j) (x_0 i h_x, y_0 j h_y))。把 (A_z) 在节点 ((x_i, y_j)) 上的值记作 (u_{i,j})。接下来是关键一步用相邻节点的函数值来近似二阶偏导数。以 (\partial^2 A_z / \partial x^2) 为例利用泰勒展开[ u_{i1,j} u_{i,j} h_x \frac{\partial u}{\partial x}\bigg|{i,j} \frac{h_x^2}{2} \frac{\partial^2 u}{\partial x^2}\bigg|{i,j} \frac{h_x^3}{6} \frac{\partial^3 u}{\partial x^3}\bigg|_{i,j} O(h_x^4) ][ u_{i-1,j} u_{i,j} - h_x \frac{\partial u}{\partial x}\bigg|{i,j} \frac{h_x^2}{2} \frac{\partial^2 u}{\partial x^2}\bigg|{i,j} - \frac{h_x^3}{6} \frac{\partial^3 u}{\partial x^3}\bigg|_{i,j} O(h_x^4) ]两式相加消去一阶导数项和三阶导数项整理得到[ \frac{\partial^2 u}{\partial x^2}\bigg|{i,j} \frac{u{i1,j} - 2u_{i,j} u_{i-1,j}}{h_x^2} O(h_x^2) ]这就是二阶中心差分格式截断误差是 (O(h_x^2))。同理对 y 方向有[ \frac{\partial^2 u}{\partial y^2}\bigg|{i,j} \frac{u{i,j1} - 2u_{i,j} u_{i,j-1}}{h_y^2} O(h_y^2) ]把这两个式子代入泊松方程忽略截断误差项就得到了离散后的代数方程[ \frac{u_{i1,j} - 2u_{i,j} u_{i-1,j}}{h_x^2} \frac{u_{i,j1} - 2u_{i,j} u_{i,j-1}}{h_y^2} -\mu J_z ]如果采用均匀正方形网格即 (h_x h_y h)方程可以进一步写成非常漂亮的形式[ 4u_{i,j} - u_{i-1,j} - u_{i1,j} - u_{i,j-1} - u_{i,j1} \mu J_z h^2 ]这个格式就是俗称的五点差分格式——中心节点 (u_{i,j}) 加上上下左右四个邻居节点。它的物理意义也很直观某一点的磁势值近似等于其四个邻居的平均值再叠加一个由电流产生的源项贡献。在没有电流的区域(u_{i,j}) 就严格等于四个邻居的平均值这正是调和函数的离散表达形式。到这里我们已经把一个偏微分方程转化成了一个线性代数方程。对求解区域内的每一个内部节点都可以列出这样一条方程把所有方程联立起来就组成了一个以 (u_{i,j}) 为未知量的大型线性方程组。方程组有多少个未知量取决于网格规模一个 (100 \times 100) 的网格就有 1 万个未知量也就是 1 万个方程。3. 五点差分格式的完整构造从公式到程序实现的映射3.1 网格编号与未知量排列在真正动手写代码之前必须先解决一个工程问题怎么把二维的网格节点存进一维的数组里。这是个看似简单但异常重要的问题因为它直接决定了系数矩阵的结构和求解效率。常见的做法是按行优先row-major编号。假设网格在 x 方向有 (N_x) 个节点编号 (i 0, 1, \dots, N_x - 1)y 方向有 (N_y) 个节点编号 (j 0, 1, \dots, N_y - 1)。那么节点 ((i, j)) 的全局序号可以定义为[ k j \times N_x i ]总未知量个数为 (N N_x \times N_y)。这样排列的好处是相邻的 x 方向邻居在数组里也是相邻的序号相差 1而 y 方向邻居序号相差 (N_x)。在组装稀疏矩阵时这种规律性能显著降低索引计算的复杂度。对于均匀正方形网格内部节点 ((i, j)) 对应的方程可以写成[ 4u_k - u_{k-1} - u_{k1} - u_{k-N_x} - u_{kN_x} \mu J_z h^2 ]注意这里 (k-1) 和 (k1) 对应左右邻居(k - N_x) 和 (k N_x) 对应上下邻居。3.2 Dirichlet 边界条件的离散实现接下来是最关键的环节之一边界条件怎么处理。静磁场仿真中最常见的边界条件有两种。第一种是 Dirichlet 边界条件也就是给边界节点直接指定磁势值 (u u_0)。典型场景是在电机分析中如果把求解区域取得足够大可以近似认为远处磁势为零或者在具有对称性的模型中对称轴上的磁力线方向已知对应某条等磁势线。处理方式非常直接把边界节点的方程直接替换为 (u_k u_0)。在系数矩阵中这相当于把该行主对角元设为 1其余元素置 0右端项设为 (u_0)。这里有一个新手常犯的错误觉得反正是边界值固定没错但内部求解时边界值不影响系数矩阵。这是不对的。当边界节点作为某个内部节点的邻居出现时内部节点的方程里包含了边界值而边界值是已知数应该挪到右端项。比如节点 ((i, j)) 的左邻居恰好是边界节点那么在它的差分方程中(u_{i-1,j}) 是已知的 (u_0)这个值不能留在未知量一侧而是要移到右端项去。如果忘了这一步你求出来的内部磁势完全是错的。3.3 Neumann 边界条件与虚节点技巧第二种是 Neumann 边界条件给的是磁势的法向导数 (\partial u / \partial n g)。在静磁场中这对应于磁场在边界上的切向分量已知。典型的场景是模型的对称边界——如果磁场关于某条线对称那么沿对称线方向磁势的法向导数为零。处理 Neumann 边界比 Dirichlet 要麻烦一些。最经典的方法是虚节点法ghost point method。以左边界(i 0)为例假设边界条件为[ \frac{\partial u}{\partial x}\bigg|_{0,j} g_j ]在 (x x_0) 处对 (u) 做泰勒展开用中心差分近似这个一阶导数[ \frac{u_{1,j} - u_{-1,j}}{2h} g_j ]其中 (u_{-1,j}) 是求解区域之外一个虚拟节点上的值。由这个式子可以得到[ u_{-1,j} u_{1,j} - 2h g_j ]然后回到五点差分格式。对于左边界节点 ((0, j))它的差分方程原本需要用到左邻居 (u_{-1,j})现在可以用上式把 (u_{-1,j}) 替换掉最终得到只含内部节点值的方程。我直接给出最终的离散结果以无电流区域、均匀网格为例[ 4u_{0,j} - 2u_{1,j} - u_{0,j-1} - u_{0,j1} -2h g_j ]注意系数从 4 变成了 4主对角元不变但右端项多了一项 (-2h g_j)。如果是零法向导数齐次 Neumann 条件右端项就为零方程简化为[ 4u_{0,j} - 2u_{1,j} - u_{0,j-1} - u_{0,j1} 0 ]这个疏远邻居的系数变化是代码里最容易排查半天的问题为什么边界附近的结果看起来不对劲先检查 Neumann 边界条件有没有把系数从 1 改成 2。3.4 内部介质分界面的处理磁导率突变时怎么做实际工程问题中求解区域往往包含多种材料——铁芯、铜线、空气、永磁体它们的磁导率差别可以达到几千倍。在介质分界面上磁通密度法向连续、磁场强度切向连续但对磁势方程来说等价于一个系数跳变的系数型偏微分方程。这时候前面推导的均匀介质五点差分格式就不能直接用了。更一般的形式是求解[ abla \cdot (\nu abla A_z) -J_z ]其中 (\nu 1/\mu) 是磁阻率。在两种介质交界处(\nu) 发生跳变。离散时需要采用调和平均来处理界面处的等效磁阻率。以 x 方向为例节点 ((i, j)) 和 ((i1, j)) 之间的界面磁阻率 (\nu_{i1/2, j})如果用算术平均会严重高估等效磁阻因为磁通优先走低磁阻路径正确做法是取调和平均[ \nu_{i1/2, j} \frac{2 \nu_{i,j} \nu_{i1,j}}{\nu_{i,j} \nu_{i1,j}} ]如果网格正好落在界面上那么 (\nu_{i1/2, j}) 就等于紧邻界面两侧的两个磁阻率的调和平均。y 方向同理。采用调和平均后离散方程变为[ \frac{\nu_{i1/2,j}(u_{i1,j} - u_{i,j}) - \nu_{i-1/2,j}(u_{i,j} - u_{i-1,j})}{h^2} \frac{\nu_{i,j1/2}(u_{i,j1} - u_{i,j}) - \nu_{i,j-1/2}(u_{i,j} - u_{i,j-1})}{h^2} -J_z ]整理成五点格式的通用形式[ a_P u_{i,j} a_E u_{i1,j} a_W u_{i-1,j} a_N u_{i,j1} a_S u_{i,j-1} b ]其中系数为[ a_E \frac{\nu_{i1/2,j}}{h^2}, \quad a_W \frac{\nu_{i-1/2,j}}{h^2}, \quad a_N \frac{\nu_{i,j1/2}}{h^2}, \quad a_S \frac{\nu_{i,j-1/2}}{h^2} ][ a_P a_E a_W a_N a_S, \quad b J_z ]这套系数写法是很多商业软件的底层实现逻辑。如果你使用 COMSOL 或 ANSYS Maxwell 这类软件界面里看到的材料分界处理本质上就是在这里调和平均的细节。4. 大型稀疏线性方程组的求解从高斯消元到 SOR 超松弛迭代4.1 为什么不能直接高斯消元把五点差分格式应用到所有内部节点会得到一个大型稀疏线性方程组 (\mathbf{A} \mathbf{u} \mathbf{f})。对于 (100 \times 100) 的网格未知量是 1 万个系数矩阵是 (10000 \times 10000)。如果是 (500 \times 500) 的网格就是 25 万个未知量矩阵规模达到 625 亿个元素。然而虽然矩阵规模惊人但矩阵中每一行最多只有 5 个非零元素。也就是说非零元素占比只有百万分之二十左右。用一个稠密矩阵来存储它是对内存的极大浪费(500 \times 500) 网格的稠密矩阵需要约 5 GB 内存双精度浮点数而稀疏存储只需要大约 (N \times 5 \times 8) 字节即 10 MB。两者相差 500 倍。更致命的是计算复杂度。高斯消元法的计算量是 (O(N^3))对于 (N 250000) 的情况(N^3 \approx 1.56 \times 10^{16})哪怕是每秒能算 1 万亿次浮点运算的计算机也要跑 4 个多小时。虽然高斯消元可以利用稀疏矩阵的填充元优化如 Cholesky 分解但在规则网格上迭代法几乎是毫无疑问的首选。4.2 雅可比迭代、高斯-赛德尔迭代与逐步更新思想迭代法的核心思想是不给一个精确解而是从一个初始猜测出发反复更新每个节点的值直到变化足够小。它的优势是每轮迭代的代价极小只需遍历所有节点做一次算术运算。雅可比迭代Jacobi iteration是最朴素的形式。根据五点差分格式[ u_{i,j}^{(k1)} \frac{1}{4} \left( u_{i-1,j}^{(k)} u_{i1,j}^{(k)} u_{i,j-1}^{(k)} u_{i,j1}^{(k)} - \mu J_z h^2 \right) ]注意右边所有值全部使用的是第 k 轮的旧值。这样做的优点是天然适合并行化但缺点也很明显——收敛速度太慢。对于 (N \times N) 的网格雅可比迭代需要大约 (O(N^2)) 次迭代才能收敛意味着网格每加密一倍迭代次数就要增加到四倍。高斯-赛德尔迭代Gauss-Seidel iteration做了一个很简单的改进在计算 (u_{i,j}^{(k1)}) 时如果它的左邻居和下方邻居在当前轮次已经被更新过了就直接用新值而不是还等着用旧值[ u_{i,j}^{(k1)} \frac{1}{4} \left( u_{i-1,j}^{(k1)} u_{i1,j}^{(k)} u_{i,j-1}^{(k1)} u_{i,j1}^{(k)} - \mu J_z h^2 \right) ]这个改动看起来微不足道但收敛速度几乎翻倍。更关键的是它不需要额外存储一份旧值副本原地更新即可内存占用更小。唯一的代价是失去了并行性因为你更新某个节点时依赖的邻居可能刚被更新过。4.3 超松弛迭代松弛因子的选取与实战建议在 Gauss-Seidel 的基础上再加一把火就得到了逐次超松弛迭代Successive Over-Relaxation, SOR。SOR 的核心思想是在每次 Gauss-Seidel 更新之后不直接接受新值而是在旧值和新值之间做一个加权平均[ u_{i,j}^{(k1)} u_{i,j}^{(k)} \omega \left( u_{i,j}^{(GS)} - u_{i,j}^{(k)} \right) ]其中 (u_{i,j}^{(GS)}) 是 Gauss-Seidel 迭代得到的中间值(\omega) 是松弛因子。当 (\omega 1) 时SOR 退化为标准 Gauss-Seidel当 (\omega 1) 时称为超松弛可以加速收敛当 (\omega 1) 时称为欠松弛通常用于迭代发散或收敛不稳定时强制压低更新幅度。对于矩形区域上的泊松方程理论上存在一个最优松弛因子[ \omega_{opt} \frac{2}{1 \sqrt{1 - \rho_{GS}^2}} ]其中 (\rho_{GS}) 是 Gauss-Seidel 迭代矩阵的谱半径。对于均匀正方形网格、矩形区域加 Dirichlet 边界谱半径可以解析表达为[ \rho_{GS} \frac{\cos\left(\frac{\pi}{N_x}\right) \cos\left(\frac{\pi}{N_y}\right)}{2} ]我以实际计算来说明问题。假设 (N_x N_y 50)则[ \rho_{GS} \cos\left(\frac{\pi}{50}\right) \approx 0.998 ]最优松弛因子为[ \omega_{opt} \frac{2}{1 \sqrt{1 - 0.998^2}} \approx \frac{2}{1 \sqrt{0.00398}} \approx \frac{2}{1 0.0631} \approx 1.881 ]这个值非常接近 2。我们再看收敛表现。Gauss-Seidel 的话每轮迭代的误差衰减因子约为 0.998也就是说想降三个数量级误差缩小到千分之一需要约[ k \approx \frac{\ln(0.001)}{\ln(0.998)} \approx \frac{-6.91}{-0.00202} \approx 3420 \text{ 次迭代} ]如果采用 SOR 且松弛因子取到最优值 1.881误差衰减因子可以降到约 (\omega_{opt} - 1 0.881)。同样的精度要求需要约[ k \approx \frac{\ln(0.001)}{\ln(0.881)} \approx \frac{-6.91}{-0.1267} \approx 55 \text{ 次迭代} ]迭代次数从 3420 次降到 55 次差了 60 多倍。这就是超松弛的价值所在——在网格规模较大的时候SOR 几乎是规则网格上的最优解。当然这个解析公式只适用于均匀网格加矩形求解区域的理想情况。对于复杂几何、多介质区域最优松弛因子无法直接解析计算。工程上的做法通常是先跑一个粗网格试算观察不同 (\omega) 值下的收敛速度找到拐点再套用到细网格上。我自己经常用 (\omega 1.7) 作为默认起点——在多数工程问题上它表现相当稳健不会有发散风险。4.4 收敛判据与停止准则不能只看迭代次数写迭代法有一个常见的纠结到底迭代多少轮才算收敛迭代次数不能一上来就拍脑袋定死而是要动态地根据残差来判断。残差的定义是[ r_{i,j} \frac{u_{i1,j} u_{i-1,j} u_{i,j1} u_{i,j-1} - 4u_{i,j}}{h^2} \mu J_z ]理论上如果 (u_{i,j}) 是精确解所有节点上的 (r_{i,j}) 都应该为零。实际迭代中会逐渐趋近于零。常用的停止准则是计算所有节点残差的某种范数比如二范数并与初始残差比较[ \frac{| \mathbf{r}^{(k)} |_2}{| \mathbf{r}^{(0)} |_2} \epsilon ]工程上我习惯取 (\epsilon 10^{-6}) 或更严格到 (10^{-8})。注意这里比较的是残差范数的相对下降而不是相邻两次迭代解的变化量。两者虽然大多数时候差别不大但在某些病态问题上解的变化量可能已经很小了但残差还很大——这时候如果你因为解没变化就停止而退出迭代算出来的磁场是错的。我踩过这个坑后面会细说。5. 完整的静磁场仿真流程从建模到后处理的落地实现5.1 仿真域的设定与边界条件选取策略一个完整的静磁场有限差分仿真流程可以归纳为六个步骤每一步都有需要注意的操作细节。首先是定义求解区域。用一个矩形区域框住你关心的电磁装置区域要取得足够大让边界影响降到可接受程度。我在做电机槽内磁场分析时通常会把边界取到装置外轮廓的 3 到 5 倍。边界外磁势衰减很快取到 5 倍之后Dirichlet 零边界引入的误差通常可以忽略。然后是确定边界条件类型。四条外边界每一条都要明确是 Dirichlet 还是 Neumann。这不能随意定必须依据对称性。以电机的一个磁极周期为例如果模型关于 x 轴对称那么在对称轴上可以设置齐次 Neumann 条件法向导数为零而在远场边界上设置 Dirichlet 零边界。边界条件选取不合理算出来的磁场分布会有明显的畸变。第三步是划分网格。均匀网格是有限差分法的天然选择但步长的选择需要权衡精度与计算量。一个粗略的指导原则是在磁场变化剧烈的区域如气隙附近、铁芯拐角处需要更细的网格。均匀网格做不了局部加密所以通常的做法是整体用细网格或者在不同区域用不同步长的网格拼接。5.2 程序设计框架与伪代码我自己常用 Python 写原型验证用 NumPy 和 SciPy 做矩阵操作和迭代求解。下面给出一个紧凑但完整的实现框架可以直接套用。import numpy as np def solve_poisson_fdm(Nx, Ny, h, Jz, mu, bc_type, bc_value): # 初始化磁势场 u np.zeros((Ny, Nx)) # 初始化磁阻率场假设均匀介质 nu np.ones((Ny, Nx)) / mu # 初始化源项 f Jz.copy() # 设置边界条件 # bc_type: dirichlet 或 neumann # bc_value: 边界值或边界法向导数值 max_iter 10000 tol 1e-8 omega 1.7 # 松弛因子 for it in range(max_iter): u_old u.copy() max_res 0.0 # 内部节点更新不含边界 for j in range(1, Ny - 1): for i in range(1, Nx - 1): res (u[j, i1] u[j, i-1] u[j1, i] u[j-1, i] - 4 * u[j, i]) / (h * h) mu * Jz[j, i] u[j, i] u_old[j, i] omega * (res * (h * h / 4)) max_res max(max_res, abs(res)) # 边界节点更新 # 左边界 j np.arange(Ny) if bc_type[0] dirichlet: u[:, 0] bc_value[0] else: # neumann u[:, 0] u[:, 1] - h * bc_value[0] # 右边界 if bc_type[1] dirichlet: u[:, -1] bc_value[1] else: u[:, -1] u[:, -2] h * bc_value[1] # 下边界 i np.arange(Nx) if bc_type[2] dirichlet: u[0, :] bc_value[2] else: u[0, :] u[1, :] - h * bc_value[2] # 上边界 if bc_type[3] dirichlet: u[-1, :] bc_value[3] else: u[-1, :] u[-2, :] h * bc_value[3] if max_res tol: print(fConverged after {it1} iterations) break return u这里的内部节点更新公式简化处理为均匀磁导率的情况。如果涉及多介质需要替换成 3.4 节给出的系数型离散公式把每个节点的邻居系数按照调和平均的磁阻率分别计算。还要特别说明的是上面的代码把 Neumann 边界条件处理成了单向差分的形式(u_{0,j} u_{1,j} - h g_j)。在一阶精度下这个处理是可用的但如果要追求更高的精度建议还是用 3.3 节介绍的虚节点中心差分法——虽然实现复杂度高一些但误差从 (O(h)) 降到了 (O(h^2))。5.3 后处理从磁势到磁通密度和磁力线求解得到磁势分布 (A_z(x, y)) 之后磁场信息还需要经过一步后处理才能提取出来。磁通密度的两个分量就是磁势的一阶偏导[ B_x \frac{\partial A_z}{\partial y}, \quad B_y -\frac{\partial A_z}{\partial x} ]数值上可以用中心差分近似。对内部节点[ (B_x){i,j} \frac{u{i,j1} - u_{i,j-1}}{2h}, \quad (B_y){i,j} -\frac{u{i1,j} - u_{i-1,j}}{2h} ]边界节点则用单侧差分近似。得到 (\mathbf{B}) 之后可以算出磁通密度幅值 (B \sqrt{B_x^2 B_y^2})用颜色图或者等值线图展示。磁力线的绘制就直接用 (A_z) 的等值线——因为在二维静磁场中磁力线就是等磁势线。还有一个很方便的指标是电感或磁链的计算。对于线圈区域磁链可以写成[ \Psi \frac{L_z}{S_c} \iint_{S_c} A_z(x, y) , dx , dy ]其中 (L_z) 是线圈在 z 方向的长度(S_c) 是线圈截面积。这个积分在离散网格上就是简单的求和。通过磁链除以电流就能得到电感值。这个量在电机和变压器设计中是很有参考价值的。6. 我踩过的坑收敛失败、网格不匹配与结果验证6.1 迭代发散排查顺序从松弛因子到边界条件我最早用 SOR 做静磁场求解时最常遇到的就是迭代发散——残差一轮比一轮大最后直接爆掉。第一次遇到这个问题的反应是调小松弛因子但这只对了一半。实际上发散的原因可能有三种排查顺序很关键。先查边界条件。如果你用的是 Dirichlet 边界但边界值设置得和内部物理过程严重不匹配比如边界全设为零但内部有强电流源那么泊松方程的解在边界附近就会被迫产生急剧变化迭代过程中这些变化会在网格上反复振荡、放大最终导致发散。解决方法是把求解区域扩大让边界远离电流源。再查松弛因子。如果边界条件没问题但依然发散把松弛因子先调回 1.0等效于 Gauss-Seidel测试。如果 Gauss-Seidel 能收敛但很慢再逐步增大 (\omega)比如 1.1、1.3、1.5观察残差曲线。如果 Gauss-Seidel 也不收敛那就不是松弛因子的问题了需要检查离散格式本身是否写错了。最后查网格均匀性。有限差分法对网格的质量要求比有限元高得多。如果网格在某个区域过度扭曲长宽比过大差分格式的局部截断误差会急剧增大导致不满足离散极大值原理迭代矩阵的主对角占优性质被破坏。出现这种情况不要犹豫重新设计网格剖分策略。6.2 解的精度验证解析解对照与能量守恒检验任何数值结果不加验证就下结论都是不负责任的。我的习惯是至少做两道验证。第一道是解析解对照。对于简单的模型比如矩形区域中心有一根无限长直导线这个问题存在解析解。设导线位于坐标原点求解区域是 ([-a, a] \times [-b, b])边界上磁势为零。此时磁势的解析解可以用镜像法或者级数展开表示[ A_z(x, y) \frac{\mu I}{4\pi} \ln \left( \frac{r_0^2}{x^2 y^2} \right) ]其中 (r_0) 是某个参考半径由边界条件确定。把它和数值解画在同一张图上看重合度这是最快的数据正确性检查。第二道是能量守恒检验。计算磁场总能量有两种途径一是从场量出发[ W_1 \frac{1}{2\mu} \iint_S |\mathbf{B}|^2 , dS ]二是从源和磁势出发[ W_2 \frac{1}{2} \iint_S J_z A_z , dS ]两者理论上相等。数值解不会严格相等但如果它们的相对偏差超过几个百分点说明网格太粗或者收敛精度不够。这个检验还有一个好处能量是一个整体量对局部小误差不敏感如果能量都对不上那问题一定很严重。我记得有一次算出来的磁场分布图看起来完全正常磁力线形状、最大磁密位置都和预期一致但能量偏差在 15% 以上。排查了很久最后发现问题出在磁导率赋值上——铁芯区域用的相对磁导率是 1000但界面处的调和平均处理写错了导致漏磁通的路径算得偏多。这再次说明图像好看不代表数值对。6.3 网格独立性分析什么时候该加密什么时候没必要最后一个非常重要的实操要点是网格独立性分析。很多人上来就用很细的网格结果算了一天一夜精度提升却微乎其微。正确的做法是从粗网格出发逐步加密观察关键量比如最大磁通密度、电感值随网格步长的变化趋势。以均匀正方形网格为例如果采用二阶中心差分误差理论上应该以 (O(h^2)) 的速度下降。也就是说步长减半误差应该缩小到原来的四分之一。用这个规律可以快速判断你的解是否已经收敛到网格无关。假设在 (h 0.1) 时算得最大磁密 (B_{max} 1.250) T加密到 (h 0.05) 时算得 (B_{max} 1.238) T再加密到 (h 0.025) 时算得 (B_{max} 1.234) T。从 0.05 到 0.025变化量只有 0.3% 左右这时候继续加密意义已经不大了。但如果从 0.1 到 0.05 变化了 1%说明粗网格的结果还不太可靠至少需要加密到 0.05 才有工程参考价值。我见过有的新手一上来就追求特别精细的网格把每根导线都单独建模、剖分得极细结果计算时间从几分钟涨到几小时精度只提升了不到 0.5%。仿真这件事先粗后细先验证后加密才是真正高效的工作流。7. 有限的网格无限的逼近关于迭代效率与扩展方向的经验之谈还是要回到那个被问过无数次的初学问题有限差分法求解泊松方程到底够不够用我的回答是在二维规则几何、均匀或分片均匀介质的静磁场问题中它的精度和效率都足够好但在复杂几何、多尺度结构面前它的劣势也很明显——规则网格在曲面边界和局部小特征面前会显得笨拙。这也是为什么真正复杂的工程仿真大多用有限元法FEM。但 FDM 的价值从来不只是能算而是它把数值方法的核心逻辑暴露得最彻底——从微分方程到差分方程从连续到离散从大型方程组到迭代收敛每一步都可以被亲手验证。那些在有限元软件里被黑盒封装掉的细节在 FDM 里全靠自己掌控。把 FDM 真正吃透之后再学习有限元、边界元或其他数值方法会有一种原来如此的豁然。我自己现在做电磁场数值计算时遇到规则几何的验证性分析仍然会先用 FDM 快速跑一版结果用来交叉验证有限元软件的输出。两个结果对上了才能放心地把数据拿去做设计决策对不上一定是某一个环节出了问题——这个排查过程往往能发现很多更深层的物理问题。最后分享一个扩展方向如果求解区域不是矩形而是圆形或扇形比如旋转电机的气隙可以考虑在极坐标系下建立差分格式。极坐标下泊松方程的离散形式会多一个 (1/r) 项边界条件处理也略有不同但核心思想完全一样。从直角坐标系扩展到极坐标系是理解有限差分法普适性的最好下一课。 有限差分法求解静磁场泊松方程：从原理到工程实践
1. 为什么静磁场仿真绕不开解泊松方程这道坎拿到静磁场仿真-主题005_有限差分法求解泊松方程这个题目时我先说一个可能让新人意外的结论静磁场仿真里真正折磨人的部分往往不是模型怎么画、网格怎么剖而是求解器在后台默默算的那套偏微分方程。换句话说懂不懂泊松方程的数值解法直接决定了你是会点仿真软件还是真的懂仿真。静磁场的控制方程从麦克斯韦方程组出发在忽略位移电流、仅考虑恒定磁场的条件下会简化成两组旋度方程[ abla \times \mathbf{H} \mathbf{J}, \quad abla \cdot \mathbf{B} 0 ]其中 (\mathbf{B} \mu \mathbf{H})。把这两组方程联立引入矢量磁势 (\mathbf{A})满足 (\mathbf{B} abla \times \mathbf{A})并选定库仑规范 ( abla \cdot \mathbf{A} 0)最终会得到一个矢量形式的泊松方程[ abla^2 \mathbf{A} -\mu \mathbf{J} ]如果是二维问题——比如电机横截面、变压器窗口、电磁铁轴向剖面这类典型场景——矢量磁势只剩一个 z 方向分量 (A_z)方程就进一步退化成标量形式[ \frac{\partial^2 A_z}{\partial x^2} \frac{\partial^2 A_z}{\partial y^2} -\mu J_z ]这就是标准形式的二维泊松方程。换句话说静磁场仿真的核心任务本质上就是在给定电流源 (J_z)、给定材料磁导率 (\mu)、给定边界条件的前提下把 (A_z) 的分布求出来。有了 (A_z)磁通密度 (\mathbf{B}) 就能通过一阶导数直接得到。而有限差分法FDM是目前求解这类方程最经典、最容易理解的数值方法之一。它的核心思路可以用八个字概括化微分为差分化连续为离散——把求解区域划分成规则的网格用网格节点上的未知量近似连续函数用差分公式替代偏导数最终把偏微分方程变成一个大型线性代数方程组。解方程组就得到每个节点上的 (A_z) 值。这篇文章我会完整过一遍从二维静磁场控制方程的推导逻辑到五点差分格式的构造过程再到 Dirichlet 边界和 Neumann 边界的离散实现、超松弛迭代SOR加速收敛的做法以及我自己在实操中积累的收敛判定、松弛因子选取、失效排查方面的经验。适合刚接触电磁场数值计算的研究生、做电机或电磁器件设计的工程师也想把这个方法吃透的仿真爱好者。这不是一篇软件操作教程而是一篇把仿真器后台到底在算什么讲清楚的原理加实战笔记。2. 二维静磁场位函数方程从 nabla 算子到可离散的差分算子2.1 麦克斯韦方程组如何一步步退化成泊松方程很多教材一上来就丢出泊松方程导致初学者总觉得它是从天而降的。实际上整个过程是可以一步步推下来的而且每一步都有明确的物理含义。先从安培环路定律的微分形式说起。在静磁场条件下电流密度不随时间变化电场旋度为零因此磁场强度 (\mathbf{H}) 与电流密度 (\mathbf{J}) 满足[ abla \times \mathbf{H} \mathbf{J} ]同时磁通密度 (\mathbf{B}) 永远是无散场磁力线闭合没有磁单极子所以有[ abla \cdot \mathbf{B} 0 ]对于一个无散矢量场总可以把它写成某个矢量场的旋度。于是引入矢量磁势 (\mathbf{A})令[ \mathbf{B} abla \times \mathbf{A} ]接着利用本构关系 (\mathbf{B} \mu \mathbf{H})代入安培环路定律的微分形式[ abla \times \left( \frac{1}{\mu} abla \times \mathbf{A} \right) \mathbf{J} ]如果材料是线性的、各向同性的且磁导率 (\mu) 在求解区域内为常数或者分片常数那么上式可以简化。利用矢量恒等式[ abla \times ( abla \times \mathbf{A}) abla ( abla \cdot \mathbf{A}) - abla^2 \mathbf{A} ]再选定库仑规范 ( abla \cdot \mathbf{A} 0) 消掉第一项就得到[ abla^2 \mathbf{A} -\mu \mathbf{J} ]这就是静磁场的矢量泊松方程。在二维问题中电流方向只有 z 方向即 (\mathbf{J} (0, 0, J_z))因此矢量磁势也只剩下 z 分量 (\mathbf{A} (0, 0, A_z))拉普拉斯算子作用在矢量上就退化为作用在标量上[ \frac{\partial^2 A_z}{\partial x^2} \frac{\partial^2 A_z}{\partial y^2} -\mu J_z ]这才到了我们真正要解的标量泊松方程。2.2 为什么在无电流区域方程会变成拉普拉斯方程这里有一个非常关键且容易忽略的细节并不是整个求解区域都有电流。比如电机气隙、铁芯、永磁体以外的空气区域这些地方 (J_z 0)方程自动退化成了拉普拉斯方程[ \frac{\partial^2 A_z}{\partial x^2} \frac{\partial^2 A_z}{\partial y^2} 0 ]这个细节在编程实现时非常重要因为你在组装系数矩阵的时候需要根据每个网格节点是否位于电流区域来决定源项是否为 0。很多新手跑出来的结果一团糟排查半天发现哦漏了判断电流区域这个坑我后面还会专门提。泊松方程统一 源项置零的处理方式比分别实现两套方程要稳得多。你在写代码时只需要维护一个 (J_z) 数组有电流的节点填真实值无电流的节点填零然后用统一的离散格式求解即可。2.3 从偏微分方程到代数方程那一步关键的离散化连续形式下(A_z(x, y)) 是定义在求解区域内每一个点上的二元函数有无限多个自由度没法直接丢给计算机。有限差分法的做法是在求解区域上画一张规则的网格比如均匀矩形网格把连续函数的值只保留在网格的节点上。如果沿着 x 方向以步长 (h_x) 划分沿着 y 方向以步长 (h_y) 划分那么节点坐标就是 ((x_i, y_j) (x_0 i h_x, y_0 j h_y))。把 (A_z) 在节点 ((x_i, y_j)) 上的值记作 (u_{i,j})。接下来是关键一步用相邻节点的函数值来近似二阶偏导数。以 (\partial^2 A_z / \partial x^2) 为例利用泰勒展开[ u_{i1,j} u_{i,j} h_x \frac{\partial u}{\partial x}\bigg|{i,j} \frac{h_x^2}{2} \frac{\partial^2 u}{\partial x^2}\bigg|{i,j} \frac{h_x^3}{6} \frac{\partial^3 u}{\partial x^3}\bigg|_{i,j} O(h_x^4) ][ u_{i-1,j} u_{i,j} - h_x \frac{\partial u}{\partial x}\bigg|{i,j} \frac{h_x^2}{2} \frac{\partial^2 u}{\partial x^2}\bigg|{i,j} - \frac{h_x^3}{6} \frac{\partial^3 u}{\partial x^3}\bigg|_{i,j} O(h_x^4) ]两式相加消去一阶导数项和三阶导数项整理得到[ \frac{\partial^2 u}{\partial x^2}\bigg|{i,j} \frac{u{i1,j} - 2u_{i,j} u_{i-1,j}}{h_x^2} O(h_x^2) ]这就是二阶中心差分格式截断误差是 (O(h_x^2))。同理对 y 方向有[ \frac{\partial^2 u}{\partial y^2}\bigg|{i,j} \frac{u{i,j1} - 2u_{i,j} u_{i,j-1}}{h_y^2} O(h_y^2) ]把这两个式子代入泊松方程忽略截断误差项就得到了离散后的代数方程[ \frac{u_{i1,j} - 2u_{i,j} u_{i-1,j}}{h_x^2} \frac{u_{i,j1} - 2u_{i,j} u_{i,j-1}}{h_y^2} -\mu J_z ]如果采用均匀正方形网格即 (h_x h_y h)方程可以进一步写成非常漂亮的形式[ 4u_{i,j} - u_{i-1,j} - u_{i1,j} - u_{i,j-1} - u_{i,j1} \mu J_z h^2 ]这个格式就是俗称的五点差分格式——中心节点 (u_{i,j}) 加上上下左右四个邻居节点。它的物理意义也很直观某一点的磁势值近似等于其四个邻居的平均值再叠加一个由电流产生的源项贡献。在没有电流的区域(u_{i,j}) 就严格等于四个邻居的平均值这正是调和函数的离散表达形式。到这里我们已经把一个偏微分方程转化成了一个线性代数方程。对求解区域内的每一个内部节点都可以列出这样一条方程把所有方程联立起来就组成了一个以 (u_{i,j}) 为未知量的大型线性方程组。方程组有多少个未知量取决于网格规模一个 (100 \times 100) 的网格就有 1 万个未知量也就是 1 万个方程。3. 五点差分格式的完整构造从公式到程序实现的映射3.1 网格编号与未知量排列在真正动手写代码之前必须先解决一个工程问题怎么把二维的网格节点存进一维的数组里。这是个看似简单但异常重要的问题因为它直接决定了系数矩阵的结构和求解效率。常见的做法是按行优先row-major编号。假设网格在 x 方向有 (N_x) 个节点编号 (i 0, 1, \dots, N_x - 1)y 方向有 (N_y) 个节点编号 (j 0, 1, \dots, N_y - 1)。那么节点 ((i, j)) 的全局序号可以定义为[ k j \times N_x i ]总未知量个数为 (N N_x \times N_y)。这样排列的好处是相邻的 x 方向邻居在数组里也是相邻的序号相差 1而 y 方向邻居序号相差 (N_x)。在组装稀疏矩阵时这种规律性能显著降低索引计算的复杂度。对于均匀正方形网格内部节点 ((i, j)) 对应的方程可以写成[ 4u_k - u_{k-1} - u_{k1} - u_{k-N_x} - u_{kN_x} \mu J_z h^2 ]注意这里 (k-1) 和 (k1) 对应左右邻居(k - N_x) 和 (k N_x) 对应上下邻居。3.2 Dirichlet 边界条件的离散实现接下来是最关键的环节之一边界条件怎么处理。静磁场仿真中最常见的边界条件有两种。第一种是 Dirichlet 边界条件也就是给边界节点直接指定磁势值 (u u_0)。典型场景是在电机分析中如果把求解区域取得足够大可以近似认为远处磁势为零或者在具有对称性的模型中对称轴上的磁力线方向已知对应某条等磁势线。处理方式非常直接把边界节点的方程直接替换为 (u_k u_0)。在系数矩阵中这相当于把该行主对角元设为 1其余元素置 0右端项设为 (u_0)。这里有一个新手常犯的错误觉得反正是边界值固定没错但内部求解时边界值不影响系数矩阵。这是不对的。当边界节点作为某个内部节点的邻居出现时内部节点的方程里包含了边界值而边界值是已知数应该挪到右端项。比如节点 ((i, j)) 的左邻居恰好是边界节点那么在它的差分方程中(u_{i-1,j}) 是已知的 (u_0)这个值不能留在未知量一侧而是要移到右端项去。如果忘了这一步你求出来的内部磁势完全是错的。3.3 Neumann 边界条件与虚节点技巧第二种是 Neumann 边界条件给的是磁势的法向导数 (\partial u / \partial n g)。在静磁场中这对应于磁场在边界上的切向分量已知。典型的场景是模型的对称边界——如果磁场关于某条线对称那么沿对称线方向磁势的法向导数为零。处理 Neumann 边界比 Dirichlet 要麻烦一些。最经典的方法是虚节点法ghost point method。以左边界(i 0)为例假设边界条件为[ \frac{\partial u}{\partial x}\bigg|_{0,j} g_j ]在 (x x_0) 处对 (u) 做泰勒展开用中心差分近似这个一阶导数[ \frac{u_{1,j} - u_{-1,j}}{2h} g_j ]其中 (u_{-1,j}) 是求解区域之外一个虚拟节点上的值。由这个式子可以得到[ u_{-1,j} u_{1,j} - 2h g_j ]然后回到五点差分格式。对于左边界节点 ((0, j))它的差分方程原本需要用到左邻居 (u_{-1,j})现在可以用上式把 (u_{-1,j}) 替换掉最终得到只含内部节点值的方程。我直接给出最终的离散结果以无电流区域、均匀网格为例[ 4u_{0,j} - 2u_{1,j} - u_{0,j-1} - u_{0,j1} -2h g_j ]注意系数从 4 变成了 4主对角元不变但右端项多了一项 (-2h g_j)。如果是零法向导数齐次 Neumann 条件右端项就为零方程简化为[ 4u_{0,j} - 2u_{1,j} - u_{0,j-1} - u_{0,j1} 0 ]这个疏远邻居的系数变化是代码里最容易排查半天的问题为什么边界附近的结果看起来不对劲先检查 Neumann 边界条件有没有把系数从 1 改成 2。3.4 内部介质分界面的处理磁导率突变时怎么做实际工程问题中求解区域往往包含多种材料——铁芯、铜线、空气、永磁体它们的磁导率差别可以达到几千倍。在介质分界面上磁通密度法向连续、磁场强度切向连续但对磁势方程来说等价于一个系数跳变的系数型偏微分方程。这时候前面推导的均匀介质五点差分格式就不能直接用了。更一般的形式是求解[ abla \cdot (\nu abla A_z) -J_z ]其中 (\nu 1/\mu) 是磁阻率。在两种介质交界处(\nu) 发生跳变。离散时需要采用调和平均来处理界面处的等效磁阻率。以 x 方向为例节点 ((i, j)) 和 ((i1, j)) 之间的界面磁阻率 (\nu_{i1/2, j})如果用算术平均会严重高估等效磁阻因为磁通优先走低磁阻路径正确做法是取调和平均[ \nu_{i1/2, j} \frac{2 \nu_{i,j} \nu_{i1,j}}{\nu_{i,j} \nu_{i1,j}} ]如果网格正好落在界面上那么 (\nu_{i1/2, j}) 就等于紧邻界面两侧的两个磁阻率的调和平均。y 方向同理。采用调和平均后离散方程变为[ \frac{\nu_{i1/2,j}(u_{i1,j} - u_{i,j}) - \nu_{i-1/2,j}(u_{i,j} - u_{i-1,j})}{h^2} \frac{\nu_{i,j1/2}(u_{i,j1} - u_{i,j}) - \nu_{i,j-1/2}(u_{i,j} - u_{i,j-1})}{h^2} -J_z ]整理成五点格式的通用形式[ a_P u_{i,j} a_E u_{i1,j} a_W u_{i-1,j} a_N u_{i,j1} a_S u_{i,j-1} b ]其中系数为[ a_E \frac{\nu_{i1/2,j}}{h^2}, \quad a_W \frac{\nu_{i-1/2,j}}{h^2}, \quad a_N \frac{\nu_{i,j1/2}}{h^2}, \quad a_S \frac{\nu_{i,j-1/2}}{h^2} ][ a_P a_E a_W a_N a_S, \quad b J_z ]这套系数写法是很多商业软件的底层实现逻辑。如果你使用 COMSOL 或 ANSYS Maxwell 这类软件界面里看到的材料分界处理本质上就是在这里调和平均的细节。4. 大型稀疏线性方程组的求解从高斯消元到 SOR 超松弛迭代4.1 为什么不能直接高斯消元把五点差分格式应用到所有内部节点会得到一个大型稀疏线性方程组 (\mathbf{A} \mathbf{u} \mathbf{f})。对于 (100 \times 100) 的网格未知量是 1 万个系数矩阵是 (10000 \times 10000)。如果是 (500 \times 500) 的网格就是 25 万个未知量矩阵规模达到 625 亿个元素。然而虽然矩阵规模惊人但矩阵中每一行最多只有 5 个非零元素。也就是说非零元素占比只有百万分之二十左右。用一个稠密矩阵来存储它是对内存的极大浪费(500 \times 500) 网格的稠密矩阵需要约 5 GB 内存双精度浮点数而稀疏存储只需要大约 (N \times 5 \times 8) 字节即 10 MB。两者相差 500 倍。更致命的是计算复杂度。高斯消元法的计算量是 (O(N^3))对于 (N 250000) 的情况(N^3 \approx 1.56 \times 10^{16})哪怕是每秒能算 1 万亿次浮点运算的计算机也要跑 4 个多小时。虽然高斯消元可以利用稀疏矩阵的填充元优化如 Cholesky 分解但在规则网格上迭代法几乎是毫无疑问的首选。4.2 雅可比迭代、高斯-赛德尔迭代与逐步更新思想迭代法的核心思想是不给一个精确解而是从一个初始猜测出发反复更新每个节点的值直到变化足够小。它的优势是每轮迭代的代价极小只需遍历所有节点做一次算术运算。雅可比迭代Jacobi iteration是最朴素的形式。根据五点差分格式[ u_{i,j}^{(k1)} \frac{1}{4} \left( u_{i-1,j}^{(k)} u_{i1,j}^{(k)} u_{i,j-1}^{(k)} u_{i,j1}^{(k)} - \mu J_z h^2 \right) ]注意右边所有值全部使用的是第 k 轮的旧值。这样做的优点是天然适合并行化但缺点也很明显——收敛速度太慢。对于 (N \times N) 的网格雅可比迭代需要大约 (O(N^2)) 次迭代才能收敛意味着网格每加密一倍迭代次数就要增加到四倍。高斯-赛德尔迭代Gauss-Seidel iteration做了一个很简单的改进在计算 (u_{i,j}^{(k1)}) 时如果它的左邻居和下方邻居在当前轮次已经被更新过了就直接用新值而不是还等着用旧值[ u_{i,j}^{(k1)} \frac{1}{4} \left( u_{i-1,j}^{(k1)} u_{i1,j}^{(k)} u_{i,j-1}^{(k1)} u_{i,j1}^{(k)} - \mu J_z h^2 \right) ]这个改动看起来微不足道但收敛速度几乎翻倍。更关键的是它不需要额外存储一份旧值副本原地更新即可内存占用更小。唯一的代价是失去了并行性因为你更新某个节点时依赖的邻居可能刚被更新过。4.3 超松弛迭代松弛因子的选取与实战建议在 Gauss-Seidel 的基础上再加一把火就得到了逐次超松弛迭代Successive Over-Relaxation, SOR。SOR 的核心思想是在每次 Gauss-Seidel 更新之后不直接接受新值而是在旧值和新值之间做一个加权平均[ u_{i,j}^{(k1)} u_{i,j}^{(k)} \omega \left( u_{i,j}^{(GS)} - u_{i,j}^{(k)} \right) ]其中 (u_{i,j}^{(GS)}) 是 Gauss-Seidel 迭代得到的中间值(\omega) 是松弛因子。当 (\omega 1) 时SOR 退化为标准 Gauss-Seidel当 (\omega 1) 时称为超松弛可以加速收敛当 (\omega 1) 时称为欠松弛通常用于迭代发散或收敛不稳定时强制压低更新幅度。对于矩形区域上的泊松方程理论上存在一个最优松弛因子[ \omega_{opt} \frac{2}{1 \sqrt{1 - \rho_{GS}^2}} ]其中 (\rho_{GS}) 是 Gauss-Seidel 迭代矩阵的谱半径。对于均匀正方形网格、矩形区域加 Dirichlet 边界谱半径可以解析表达为[ \rho_{GS} \frac{\cos\left(\frac{\pi}{N_x}\right) \cos\left(\frac{\pi}{N_y}\right)}{2} ]我以实际计算来说明问题。假设 (N_x N_y 50)则[ \rho_{GS} \cos\left(\frac{\pi}{50}\right) \approx 0.998 ]最优松弛因子为[ \omega_{opt} \frac{2}{1 \sqrt{1 - 0.998^2}} \approx \frac{2}{1 \sqrt{0.00398}} \approx \frac{2}{1 0.0631} \approx 1.881 ]这个值非常接近 2。我们再看收敛表现。Gauss-Seidel 的话每轮迭代的误差衰减因子约为 0.998也就是说想降三个数量级误差缩小到千分之一需要约[ k \approx \frac{\ln(0.001)}{\ln(0.998)} \approx \frac{-6.91}{-0.00202} \approx 3420 \text{ 次迭代} ]如果采用 SOR 且松弛因子取到最优值 1.881误差衰减因子可以降到约 (\omega_{opt} - 1 0.881)。同样的精度要求需要约[ k \approx \frac{\ln(0.001)}{\ln(0.881)} \approx \frac{-6.91}{-0.1267} \approx 55 \text{ 次迭代} ]迭代次数从 3420 次降到 55 次差了 60 多倍。这就是超松弛的价值所在——在网格规模较大的时候SOR 几乎是规则网格上的最优解。当然这个解析公式只适用于均匀网格加矩形求解区域的理想情况。对于复杂几何、多介质区域最优松弛因子无法直接解析计算。工程上的做法通常是先跑一个粗网格试算观察不同 (\omega) 值下的收敛速度找到拐点再套用到细网格上。我自己经常用 (\omega 1.7) 作为默认起点——在多数工程问题上它表现相当稳健不会有发散风险。4.4 收敛判据与停止准则不能只看迭代次数写迭代法有一个常见的纠结到底迭代多少轮才算收敛迭代次数不能一上来就拍脑袋定死而是要动态地根据残差来判断。残差的定义是[ r_{i,j} \frac{u_{i1,j} u_{i-1,j} u_{i,j1} u_{i,j-1} - 4u_{i,j}}{h^2} \mu J_z ]理论上如果 (u_{i,j}) 是精确解所有节点上的 (r_{i,j}) 都应该为零。实际迭代中会逐渐趋近于零。常用的停止准则是计算所有节点残差的某种范数比如二范数并与初始残差比较[ \frac{| \mathbf{r}^{(k)} |_2}{| \mathbf{r}^{(0)} |_2} \epsilon ]工程上我习惯取 (\epsilon 10^{-6}) 或更严格到 (10^{-8})。注意这里比较的是残差范数的相对下降而不是相邻两次迭代解的变化量。两者虽然大多数时候差别不大但在某些病态问题上解的变化量可能已经很小了但残差还很大——这时候如果你因为解没变化就停止而退出迭代算出来的磁场是错的。我踩过这个坑后面会细说。5. 完整的静磁场仿真流程从建模到后处理的落地实现5.1 仿真域的设定与边界条件选取策略一个完整的静磁场有限差分仿真流程可以归纳为六个步骤每一步都有需要注意的操作细节。首先是定义求解区域。用一个矩形区域框住你关心的电磁装置区域要取得足够大让边界影响降到可接受程度。我在做电机槽内磁场分析时通常会把边界取到装置外轮廓的 3 到 5 倍。边界外磁势衰减很快取到 5 倍之后Dirichlet 零边界引入的误差通常可以忽略。然后是确定边界条件类型。四条外边界每一条都要明确是 Dirichlet 还是 Neumann。这不能随意定必须依据对称性。以电机的一个磁极周期为例如果模型关于 x 轴对称那么在对称轴上可以设置齐次 Neumann 条件法向导数为零而在远场边界上设置 Dirichlet 零边界。边界条件选取不合理算出来的磁场分布会有明显的畸变。第三步是划分网格。均匀网格是有限差分法的天然选择但步长的选择需要权衡精度与计算量。一个粗略的指导原则是在磁场变化剧烈的区域如气隙附近、铁芯拐角处需要更细的网格。均匀网格做不了局部加密所以通常的做法是整体用细网格或者在不同区域用不同步长的网格拼接。5.2 程序设计框架与伪代码我自己常用 Python 写原型验证用 NumPy 和 SciPy 做矩阵操作和迭代求解。下面给出一个紧凑但完整的实现框架可以直接套用。import numpy as np def solve_poisson_fdm(Nx, Ny, h, Jz, mu, bc_type, bc_value): # 初始化磁势场 u np.zeros((Ny, Nx)) # 初始化磁阻率场假设均匀介质 nu np.ones((Ny, Nx)) / mu # 初始化源项 f Jz.copy() # 设置边界条件 # bc_type: dirichlet 或 neumann # bc_value: 边界值或边界法向导数值 max_iter 10000 tol 1e-8 omega 1.7 # 松弛因子 for it in range(max_iter): u_old u.copy() max_res 0.0 # 内部节点更新不含边界 for j in range(1, Ny - 1): for i in range(1, Nx - 1): res (u[j, i1] u[j, i-1] u[j1, i] u[j-1, i] - 4 * u[j, i]) / (h * h) mu * Jz[j, i] u[j, i] u_old[j, i] omega * (res * (h * h / 4)) max_res max(max_res, abs(res)) # 边界节点更新 # 左边界 j np.arange(Ny) if bc_type[0] dirichlet: u[:, 0] bc_value[0] else: # neumann u[:, 0] u[:, 1] - h * bc_value[0] # 右边界 if bc_type[1] dirichlet: u[:, -1] bc_value[1] else: u[:, -1] u[:, -2] h * bc_value[1] # 下边界 i np.arange(Nx) if bc_type[2] dirichlet: u[0, :] bc_value[2] else: u[0, :] u[1, :] - h * bc_value[2] # 上边界 if bc_type[3] dirichlet: u[-1, :] bc_value[3] else: u[-1, :] u[-2, :] h * bc_value[3] if max_res tol: print(fConverged after {it1} iterations) break return u这里的内部节点更新公式简化处理为均匀磁导率的情况。如果涉及多介质需要替换成 3.4 节给出的系数型离散公式把每个节点的邻居系数按照调和平均的磁阻率分别计算。还要特别说明的是上面的代码把 Neumann 边界条件处理成了单向差分的形式(u_{0,j} u_{1,j} - h g_j)。在一阶精度下这个处理是可用的但如果要追求更高的精度建议还是用 3.3 节介绍的虚节点中心差分法——虽然实现复杂度高一些但误差从 (O(h)) 降到了 (O(h^2))。5.3 后处理从磁势到磁通密度和磁力线求解得到磁势分布 (A_z(x, y)) 之后磁场信息还需要经过一步后处理才能提取出来。磁通密度的两个分量就是磁势的一阶偏导[ B_x \frac{\partial A_z}{\partial y}, \quad B_y -\frac{\partial A_z}{\partial x} ]数值上可以用中心差分近似。对内部节点[ (B_x){i,j} \frac{u{i,j1} - u_{i,j-1}}{2h}, \quad (B_y){i,j} -\frac{u{i1,j} - u_{i-1,j}}{2h} ]边界节点则用单侧差分近似。得到 (\mathbf{B}) 之后可以算出磁通密度幅值 (B \sqrt{B_x^2 B_y^2})用颜色图或者等值线图展示。磁力线的绘制就直接用 (A_z) 的等值线——因为在二维静磁场中磁力线就是等磁势线。还有一个很方便的指标是电感或磁链的计算。对于线圈区域磁链可以写成[ \Psi \frac{L_z}{S_c} \iint_{S_c} A_z(x, y) , dx , dy ]其中 (L_z) 是线圈在 z 方向的长度(S_c) 是线圈截面积。这个积分在离散网格上就是简单的求和。通过磁链除以电流就能得到电感值。这个量在电机和变压器设计中是很有参考价值的。6. 我踩过的坑收敛失败、网格不匹配与结果验证6.1 迭代发散排查顺序从松弛因子到边界条件我最早用 SOR 做静磁场求解时最常遇到的就是迭代发散——残差一轮比一轮大最后直接爆掉。第一次遇到这个问题的反应是调小松弛因子但这只对了一半。实际上发散的原因可能有三种排查顺序很关键。先查边界条件。如果你用的是 Dirichlet 边界但边界值设置得和内部物理过程严重不匹配比如边界全设为零但内部有强电流源那么泊松方程的解在边界附近就会被迫产生急剧变化迭代过程中这些变化会在网格上反复振荡、放大最终导致发散。解决方法是把求解区域扩大让边界远离电流源。再查松弛因子。如果边界条件没问题但依然发散把松弛因子先调回 1.0等效于 Gauss-Seidel测试。如果 Gauss-Seidel 能收敛但很慢再逐步增大 (\omega)比如 1.1、1.3、1.5观察残差曲线。如果 Gauss-Seidel 也不收敛那就不是松弛因子的问题了需要检查离散格式本身是否写错了。最后查网格均匀性。有限差分法对网格的质量要求比有限元高得多。如果网格在某个区域过度扭曲长宽比过大差分格式的局部截断误差会急剧增大导致不满足离散极大值原理迭代矩阵的主对角占优性质被破坏。出现这种情况不要犹豫重新设计网格剖分策略。6.2 解的精度验证解析解对照与能量守恒检验任何数值结果不加验证就下结论都是不负责任的。我的习惯是至少做两道验证。第一道是解析解对照。对于简单的模型比如矩形区域中心有一根无限长直导线这个问题存在解析解。设导线位于坐标原点求解区域是 ([-a, a] \times [-b, b])边界上磁势为零。此时磁势的解析解可以用镜像法或者级数展开表示[ A_z(x, y) \frac{\mu I}{4\pi} \ln \left( \frac{r_0^2}{x^2 y^2} \right) ]其中 (r_0) 是某个参考半径由边界条件确定。把它和数值解画在同一张图上看重合度这是最快的数据正确性检查。第二道是能量守恒检验。计算磁场总能量有两种途径一是从场量出发[ W_1 \frac{1}{2\mu} \iint_S |\mathbf{B}|^2 , dS ]二是从源和磁势出发[ W_2 \frac{1}{2} \iint_S J_z A_z , dS ]两者理论上相等。数值解不会严格相等但如果它们的相对偏差超过几个百分点说明网格太粗或者收敛精度不够。这个检验还有一个好处能量是一个整体量对局部小误差不敏感如果能量都对不上那问题一定很严重。我记得有一次算出来的磁场分布图看起来完全正常磁力线形状、最大磁密位置都和预期一致但能量偏差在 15% 以上。排查了很久最后发现问题出在磁导率赋值上——铁芯区域用的相对磁导率是 1000但界面处的调和平均处理写错了导致漏磁通的路径算得偏多。这再次说明图像好看不代表数值对。6.3 网格独立性分析什么时候该加密什么时候没必要最后一个非常重要的实操要点是网格独立性分析。很多人上来就用很细的网格结果算了一天一夜精度提升却微乎其微。正确的做法是从粗网格出发逐步加密观察关键量比如最大磁通密度、电感值随网格步长的变化趋势。以均匀正方形网格为例如果采用二阶中心差分误差理论上应该以 (O(h^2)) 的速度下降。也就是说步长减半误差应该缩小到原来的四分之一。用这个规律可以快速判断你的解是否已经收敛到网格无关。假设在 (h 0.1) 时算得最大磁密 (B_{max} 1.250) T加密到 (h 0.05) 时算得 (B_{max} 1.238) T再加密到 (h 0.025) 时算得 (B_{max} 1.234) T。从 0.05 到 0.025变化量只有 0.3% 左右这时候继续加密意义已经不大了。但如果从 0.1 到 0.05 变化了 1%说明粗网格的结果还不太可靠至少需要加密到 0.05 才有工程参考价值。我见过有的新手一上来就追求特别精细的网格把每根导线都单独建模、剖分得极细结果计算时间从几分钟涨到几小时精度只提升了不到 0.5%。仿真这件事先粗后细先验证后加密才是真正高效的工作流。7. 有限的网格无限的逼近关于迭代效率与扩展方向的经验之谈还是要回到那个被问过无数次的初学问题有限差分法求解泊松方程到底够不够用我的回答是在二维规则几何、均匀或分片均匀介质的静磁场问题中它的精度和效率都足够好但在复杂几何、多尺度结构面前它的劣势也很明显——规则网格在曲面边界和局部小特征面前会显得笨拙。这也是为什么真正复杂的工程仿真大多用有限元法FEM。但 FDM 的价值从来不只是能算而是它把数值方法的核心逻辑暴露得最彻底——从微分方程到差分方程从连续到离散从大型方程组到迭代收敛每一步都可以被亲手验证。那些在有限元软件里被黑盒封装掉的细节在 FDM 里全靠自己掌控。把 FDM 真正吃透之后再学习有限元、边界元或其他数值方法会有一种原来如此的豁然。我自己现在做电磁场数值计算时遇到规则几何的验证性分析仍然会先用 FDM 快速跑一版结果用来交叉验证有限元软件的输出。两个结果对上了才能放心地把数据拿去做设计决策对不上一定是某一个环节出了问题——这个排查过程往往能发现很多更深层的物理问题。最后分享一个扩展方向如果求解区域不是矩形而是圆形或扇形比如旋转电机的气隙可以考虑在极坐标系下建立差分格式。极坐标下泊松方程的离散形式会多一个 (1/r) 项边界条件处理也略有不同但核心思想完全一样。从直角坐标系扩展到极坐标系是理解有限差分法普适性的最好下一课。