同伦延拓破解静磁场非线性仿真收敛难题

做静磁场仿真，尤其是变压器、电机这类带铁磁材料的模型，碰到非线性迭代不收敛几乎是家常便饭。很多人以为把B-H曲线和网格都给对了，剩下的就交给求解器，结果牛顿迭代在半路发疯，要么震荡要么直接溢出。这次分享的是我基于“静磁场仿真-主题070_同伦延拓求解非线性问题”的完整技术记录，核心就是用同伦延拓（Homotopy Continuation）把强非线性、大励磁工况下的静磁场方程一步步“扶上马”。文章会从非线性问题的物理来源入手，把同伦延拓的数学思路拆开讲透，再给出一套我自己常用的落地实现流程和调试策略。如果你正在做电磁有限元仿真，又被收敛性问题卡得焦头烂额，这篇文章应该能帮你打开一条新路。

1. 先把非线性静磁场为什么难算说清楚

1.1 非线性是哪里来的

说同伦延拓之前，必须先把静磁场的非线性本质讲清楚。电磁场仿真里用的方程组并不复杂，核心就是：

∇ × H = J

∇ · B = 0

本构关系 B = μ(H) H

前两个方程是麦克斯韦方程组的静磁版本，第三个方程才是一切麻烦的源头。空气的磁导率是常数 μ0，但铁磁材料不一样，它的 μ 会随着磁场强度 H 的增大而明显下降。低场强时，硅钢的相对磁导率能到两三千，一旦磁场加大接近饱和区，μr 可能掉到几十甚至十几。B-H 曲线因此不是一条直线，而是一条先陡后缓、最终接近水平饱和段的曲线。

在有限元实现里，更习惯用磁阻率 ν 代替磁导率 μ，ν = 1/μ。引入磁矢量位 A，令 B = ∇ × A，静磁场方程就能改写成：

∇ × [ν(B) ∇ × A] = J

这一步非常关键，因为 ν 本身依赖 B，而 B 又由 A 的空间导数决定。把问题离散到有限元网格上之后，最终会得到这样一个代数系统：

S(A) A = b

这里的 S(A) 是刚度矩阵，b 是激励向量。注意，矩阵 S 是解向量 A 的函数。普通线性问题是一次方程，一次性求解就行；而这里未知量 A 出现在等号左边，同时又藏在矩阵系数里，所以我们只能通过迭代逼近真实解。

很多人第一次接触这个方程会觉得奇怪：“矩阵和未知数捆在一起，怎么解？”这就是非线性有限元的核心困难。求解器每一步都在做同一个事情：给一个当前的 A，算一下矩阵系数，再求一个线性系统，得到新的 A，然后重复。问题在于，这个迭代不保证收敛，尤其当初始猜测离真实解很远的时候。

1.2 常规牛顿法在哪些地方翻车

目前绝大多数静磁场有限元求解器都用牛顿法或牛顿法变体。牛顿法的思路很清晰：在某个近似解 A_k 附近，把 F(A)=0 做 Taylor 展开，忽略二阶以上项，得到一个线性方程组：

F(A_k) + J(A_k) ΔA = 0

其中 J 是 Jacobian 矩阵，ΔA 是修正步。解出 ΔA 之后更新：

A_{k+1} = A_k + ΔA

牛顿法有一个很漂亮的特性：如果初始猜测足够靠近真实解，它会以二次速度收敛，通常四五步就能把残差压到机器精度。但前提是“初始猜测足够靠近真实解”，这句话在强非线性静磁场问题里往往不成立。

我实际调试中总结出四个最容易翻车的场景。

第一，初值给得太差。很多人习惯从零场开始迭代，这在线性不太强的模型中没问题，但如果励磁电流很大，第一次计算就可能把铁芯磁密推到好几特斯拉，远超实际饱和点，整个 B-H 曲线外推部分完全失控，后续计算自然越走越偏。

第二，饱和区 Jacobian 矩阵条件数恶化。饱和以后 ν 迅速增大，刚度矩阵对角元素和非对角元素的比值发生变化，矩阵接近奇异。这时候即使是牛顿法的线性代数求解环节本身，也会因为条件数太差而给出不可靠的修正步。

第三，B-H 曲线处理不光滑。材料库提供的 B-H 数据往往是几十个离散点，如果直接用分段线性插值，ν 的导数在数据点之间是跳变的。而雅可比矩阵里需要 dν/dB 这一项，导数一抖，牛顿迭代就容易进入振荡，永不收敛。

第四，接近磁通突变或分支点。某些几何结构在特定励磁下会出现非线性响应转折，雅可比矩阵奇异，普通牛顿法直接丢失解。

遇到这些情况，常规反应是加密网格、调松弛因子、换初值，这些手段有时候管用，但碰到强非线性问题往往只是拖延时间。我们需要的是一个从原理上更稳妥的思路，这就是同伦延拓登场的理由。

2. 同伦延拓到底在做什么

2.1 用“变形”代替“硬解”

同伦延拓的思想用一句话概括就是：不要直接解难问题，把难问题和简单问题用参数连起来，再慢慢变形过去。

这个思想在拓扑里叫同伦，本身是指两个映射之间可以连续变化。拿来做数值求解时，我们要解的方程是：

F(A) = 0

这个方程太难，找不到合适的迭代初值。那我们就造一个容易解的方程 G(A)=0，然后把二者组合：

H(A, τ) = (1-τ) G(A) + τ F(A)，τ ∈ [0,1]

当 τ=0 时，H=G，很好解；当 τ=1 时，H=F，回到目标问题。τ 从 0 慢慢变到 1 的过程中，我们每增加一点 τ，就把前一时刻的解当作当前时刻的初值，再去解 H(A,τ)=0。

这样一个参数一个参数地递推过去，就是“延拓”。它相当于在山脚下找到一条缓坡路，每一步都小到不会偏离路径，最终走到山顶。

在静磁场问题里，这个思路对应非常自然的物理图景。最简单的同伦构型是负载步进：

H(A, τ) = S(A) A - τ b = 0

τ=0 时没有鼓励，解是 A=0；τ 逐步增大，相当于把绕组电流一点点加满。这个方案我先给结果：好用，直观，而且绝大多数静磁场问题都可以这样解决。

2.2 三种常用同伦构型的选型对比

实际工程实现里，我见过三种同伦构型用得最多。

第一种就是上面说的负载步进。它适合电机、变压器、螺线管这类有明确电流激励的问题。物理意义清晰，初始解零场求，几乎不需要额外处理。缺点是如果材料本身是永磁体、或模型里存在剩磁，没有外部激励时磁场也不为零，τ=0 的解并不是一个自然的起点。

第二种是材料参数同伦。把实际非线性材料的磁阻率做一个带参数的插值：

ν_eff(τ) = τ ν(B) + (1-τ) ν_linear

τ=0 时全模型用线性磁导率，解起来非常轻松；τ 增加到 1 时材料逐渐变回真实非线性 B-H 曲线。这种做法的好处是不依赖外部激励，对永磁结构一样适用。代价是中间状态没有物理意义，纯粹是一种数值手段。

第三种是经典的牛顿同伦，也叫全局同伦：

H(A, τ) = F(A) - (1-τ) F(A₀)

其中 A₀ 是初始猜测。τ=0 时，H(A₀,0)=0，A₀ 天然就是一个解；τ=1 时回到目标方程。这种构型理论上是通用的，但实际用起来我发现中间路径经常绕弯，收敛稳定性反而比物理同伦差。所以我个人只在无法构造前两种构型时才会用它。

三选一的经验法则：有明确励磁先选负载步进；有永磁或零载磁场再考虑材料参数同伦；前两种都构建困难时，才轮到牛顿同伦兜底。

2.3 为什么多绕一步反而更稳

有人会问：每一步都走在牛顿迭代里，同伦只是换个方式工作，它不是多此一举吗？

关键区别在于每步的起点位置。直接求解 F(A)=0 时，我们从某种猜测出发，但这个猜测可能距离真实解十万八千里，牛顿法跨出的一大步可能把解带向完全错误的区域。同伦延拓通过控制 τ 的步幅，保证了每一步开始时的近似解和这一步的目标解足够接近——这个“足够接近”满足牛顿法的局部收敛条件。

再强调一下本质：牛顿法是寻找局部线性化后的解，它天然是局部的；而同伦延拓把局部线性化的范围通过参数连续覆盖到了全域。所以延拓不是替代牛顿法，而是给牛顿法不断提供合格的初值。

用个不恰当的类比，爬陡坡时你要么选择从坡脚一步步小碎步走上去，每一步都稳稳踏在坡面上；要么试图三步并作两步冲上去，很可能在半路失去平衡。同伦延拓就是前者。

代价当然也有：延拓过程要解若干次牛顿迭代，计算量会变大。但相比牛顿法直接发散、报错、整个仿真白跑，多算几次数值迭代的代价完全值得。

3. 从公式到程序：同伦延拓的落地细节

3.1 有限元离散后的非线性系统到底长什么样

我们要把上文的方程真正写到程序里。以二维静磁场问题为例，磁矢量位只剩 z 方向分量，设为 A_z(x,y)。单元采用三角单元，对每个单元 e，单元刚度矩阵的元素是：

S^e_ij = ∫_Ω^e ν(B^e) (∂N_i/∂x ∂N_j/∂x + ∂N_i/∂y ∂N_j/∂y) dxdy

其中 N_i 是形函数，B^e 是这个单元内的磁感应强度模值。对二维线性三角单元来说，B 在单元内是常数，所以 ν(B^e) 可以提到积分号外面。

组装完整系统后得到残差方程：

F(A) = S(A) A - b = 0

b 节点由电流密度 J 积分而来。

求 Newton 迭代所需的 Jacobian 时，除了保留上面的 S(A)，还必须加入因为 ν 随 A 变化引起的修正项：

∂F_i/∂A_j = S_ij + ∫_Ω^e (∂ν/∂B) (∂B/∂A_j) (∇N_i · ∇A) dxdy

在二维单元里，B² = (∂A/∂y)² + (∂A/∂x)²，对 A_j 求导后可以得到：

∂B/∂A_j = (∇A · ∇N_j) / B

于是雅可比矩阵在每个单元上的贡献是：

J^e_ij = ν ∫ ∇N_i · ∇N_j dΩ + (∂ν/∂B) / B ^e ∫ (∇A · ∇N_i)(∇A · ∇N_j) dΩ

这个第二项千万不能漏。在饱和区，∂ν/∂B 数值可以很大，它主导了整个雅可比矩阵的行为。我见过不少人在实现时只用了第一项近似，结果牛顿迭代收敛极慢，而且一到强饱和就发散。真正的收敛域和收敛速度差距明显。

3.2 外层步长怎么控制

同伦延拓的实现框架大概是这样的：

import numpy as np def solve_magnetostatic_homotopy(mesh, current_load, bh_model, t_end=1.0, dt_init=0.1): A = np.zeros(mesh.n_node) t = 0.0 dt = dt_init while t < t_end - 1e-12: t_target = min(t + dt, t_end) converged, A, n_iter = newton_at_fixed_t( A, mesh, current_load, t_target, bh_model ) if not converged: dt *= 0.5 if dt < 1e-6: raise RuntimeError("延拓步长已无法再缩小，问题可能已经病态") continue if n_iter <= 3: dt = min(dt * 1.5, 0.4) else: dt = max(dt * 0.7, dt_init * 0.01) t = t_target return A

外层循环控制 τ 的推进。每次推进的目标是某个固定的 t_target，然后用当前解 A 作为初始猜测去解 H(A, t_target)=0。如果成功，根据内层迭代次数决定下一步走大步还是小步；如果失败，就把目标步长减半，原地重试。

这个自适应策略很粗糙，但非常有效。它的逻辑是：收敛快说明路径平坦，可以走快点；收敛慢甚至失败说明遇到了非线性强的区域，要放慢脚步。

几个经验值供参考：初始步长 dt_init 通常取 0.05~0.1；成功且迭代不超过 3 次时，最大步长我限制到 0.4，避免一次性推进太多被打回原形；失败时直接把步长砍半，最多允许退到 1e-6，再小就真的没法算了。

3.3 内层牛顿与线搜索的组合细节

内层牛顿决不能是裸的牛顿法。我的建议是在每个内层步里加一条简单的线搜索，公式是：

A_{k+1} = A_k + α ΔA，α ∈ (0,1]

先从 α=1 开始，计算新残差的范数，如果残差不下降，就把 α 减半，最多减五六次。伪代码如下：

def newton_at_fixed_t(A0, mesh, load, t, bh_model, max_iter=30, tol_rel=1e-5): A = A0.copy() norm0 = None for it in range(max_iter): r, J = assemble_residual_and_jacobian(A, mesh, t * load, bh_model) norm_r = np.linalg.norm(r) if it == 0: norm0 = norm_r if norm0 < 1e-14: return True, A, it if norm_r <= tol_rel * norm0: return True, A, it delta = np.linalg.solve(J, -r) alpha = 1.0 A_trial = A + alpha * delta r_trial, _ = assemble_residual_and_jacobian( A_trial, mesh, t * load, bh_model ) count = 0 while (np.linalg.norm(r_trial) >= norm_r) and (count < 6): alpha *= 0.5 A_trial = A + alpha * delta r_trial, _ = assemble_residual_and_jacobian( A_trial, mesh, t * load, bh_model ) count += 1 A = A + alpha * delta return False, A, max_iter

别小看这个 α 折半过程。饱和区附近的牛顿修正步经常偏大，直接加上去残差不降反升，如果没有线搜索，整个迭代序列马上发散。加了线搜索以后，我们至少保证每一步都在向减小残差的方向走。

另外内层迭代的收敛判据，我不用固定的绝对残差阈值，而是用相对下降量：

‖r_k‖ ≤ tol_rel × ‖r_0‖

因为不同网格单元数量、不同激励大小，残差绝对值可以差好几个数量级。相对下降更容易跨网格适用。中间延拓阶段 tol_rel 设 1e-3 或 1e-4 就够了，最后 τ=1 时必须收紧到 1e-6，否则最终解精度不够，后处理出来的磁密可能带着明显误差。

4. 一个二维电磁铁实例走一遍

4.1 模型参数与 B-H 曲线设定

理论讲得再多，不如跑一个真实案例。我随手搭了一个 C 形电磁铁：

• 铁芯截面 20 mm × 20 mm
• 气隙 2 mm
• 磁路长度约 160 mm
• 绕组安匝数 3000 A
• 铁芯材料为电工钢，B-H 曲线用三次样条平滑，数据覆盖到 2.2 T

网格采用三角形单元，约 2400 个节点。气隙局部加密，因为整个模型的磁压降主要集中在气隙区域，网格太少会把磁密分布算得过于平滑，影响饱和位置的判断。

两根关键设定说明一下。第一，安匝数 3000 A 是故意选到让铁芯进入明显饱和的水平，这样同伦延拓才有用武之地。第二，B-H 曲线我用连续导数光滑的三次样条插值，而不是最简单的折线连接。前面说过，雅可比矩阵里需要 dν/dB，如果导数不连续，牛顿迭代靠近解的时候会异常振荡。这是一个必须提前处理的问题，不要等求解报了发散再去查。

4.2 逐步求解的过程记录

我用上一节的代码跑了一遍，记录下了每一步的收敛细节：

这个表很有意思。τ 从 0 到 0.2 时，磁路基本处于线性段，磁密从 0 升到 1.4 T，每步牛顿迭代都在 6 次以内，走得很轻松。到 0.35 这一步，9 次迭代才收敛，明显感觉到拐点区的非线性变强。

到了尝试 τ=0.5 的时候，我故意让固定步长去冲一个大步，结果牛顿迭代第一轮残差反而上升，最后直接判定失败。这一步失败很正常，因为我们从 0.35 跳到 0.5，落点正好处在 B-H 曲线最陡的饱和拐点之后，雅可比矩阵变化太剧烈，初值已经不在收敛域里了。

自动步长控制生效后，步长折半到 0.075，从 0.35 原地重试到 0.425，成功。之后系统重新提速，虽然中间的迭代次数始终偏多，但再没有发散。

如果不用同伦，直接对满载 τ=1 做一次牛顿法，我实测的结果非常难看：第一次迭代算出的 B_max 超过 7 T，残差不降反升，后面连续五六步都在振荡，最终报错。这个对比足以说明延拓的价值。

4.3 最终解与物理合理性检查

τ=1 处的解出来后，不能直接收工。我习惯做两个后处理检查。

第一个看数值指标：气隙磁场约为 1.72 T，铁芯内部平均磁密约 1.79 T，最大磁密 1.85 T，出现在铁芯转角处。这个结果和磁路估算接近。因为气隙只有 2 mm，磁压降主要集中在气隙上，铁芯截面面积又与气隙面积相同，所以气隙磁密略低于铁芯内部，符合预期。

第二个做全局守恒检查：沿着铁芯中心线做一次安培环路积分，即把 H 沿闭合路径的点积分求和，应该约等于 3000 A。实际积分结果是 2985 A，误差不到 1%，说明方程组被正确求解。如果这一步误差过大，我会先怀疑边界条件施加有误，或者某些单元的 B-H 外推超出了有效范围。

还有一个细节值得注意：当 τ 接近 1，材料已经深度饱和，磁阻率上升，矩阵的条件数比线性阶段差得多。这种情况下，我建议物理解算使用直接求解器。迭代求解器在强饱和时收敛速度会显著下降，甚至不收敛，直接求解器虽然单步贵一点，但胜在稳。

5. 常见问题与排查实录

5.1 收敛问题速查表

同伦延拓也不是万灵丹，它在运行中同样会碰到各种奇怪现象。这里把我调试中常见的现象、原因和策略整理成一个速查表：

这个表里的每一条我都踩过。最坑的是第一行：有人看到 τ=0 时方程几乎是线性，就大胆把初始步长设成 0.3，结果从零解跳到 30% 负载，铁芯直接进入饱和区，牛顿法第一步就把解带飞。所以初始步长宁可保守。

第二行也很经典。某一步牛顿迭代失败，不一定是算法坏了，而是路径在这个区域曲率变大，步长需要缩小。自适应步长正是为了这个场景设计的，不要手动硬扛。

5.2 两个容易被忽视的坑

第一个坑：只组装残差，不重新组装雅可比。

我在调试早期犯过一个错误：为了省时间，牛顿迭代里雅可比矩阵只在每个 τ 步开始时组装一次，后面内层迭代一直复用。结果残差下降两个数量级后就再也不动了，怎么增加迭代步数都没用。原因很简单，非线性系统里雅可比矩阵也是解的函数，必须每步更新。后来我看到代码里复用了旧雅可比，立刻改成每次迭代重新组装，收敛速度马上恢复正常。

第二个坑：对 B-H 曲线直接用分段线性插值。

很多初学者从材料库拿到的是离散的 B-H 数据，直接用线性分段连起来就丢给求解器。在强饱和区附近，相邻数据点之间导数变化非常大，甚至出现折线拐点。这会让牛顿法陷入一种“来回弹跳”的状态，残差不再下降。我用三次样条或者 Akima 插值处理 B-H 数据，同时保证导数连续，这个问题就消失了。

如果你必须保留原始数据点作为拟合依据，也千万不要让样条出现过冲。B-H 曲线的过冲会产生伪负微分磁导率，物理上完全不可接受，求解结果也会特别诡异。

6. 个人经验收尾：同伦延拓是备胎还是主力

用的次数多了，我对同伦延拓的定位是：它应该作为求解器里的“保底策略”，而不是每个问题都主动上。对大多数常规静磁场模型，一个好的阻尼牛顿法加连续光滑的 B-H 数据，三四步就能收敛，没必要引入延拓。但是一旦这些常规手段失败，同伦延拓就是我最先调用的救场工具。

我养成了一个调试习惯：面对新模型时，先跑一次同伦延拓，拿到一个可靠的参考解，然后再关掉延拓，直接跑普通牛顿，观察它在哪个环节出问题、收敛域边界在哪里。这个过程能帮你积累很多关于模型本身的认识。如果一个模型在同伦延拓的辅助下依然发散，那我会立刻停止调算法，转而检查建模问题——大概率是网格有负面积、气隙没加密、材料数据区域有空缺，或者边界条件加载有误。

最后分享一个小技巧：把每个 τ 步的 B_max 和内层迭代次数打出来。这两个量是极好的“探路石”，能清晰地告诉你整个求解路径上的难度分布。什么时候该放大步长、什么时候该收着走，一看输出就知道。比盲试参数高效得多。

这个方法我也用在别的非线性问题上，比如超导磁体的场分布计算、磁滞损耗近似模型，原理相通。非线性方程组的求解思路是通用的，同伦延拓只是其中最可靠的一条路，希望这次记录能帮你少走点弯路。