做静磁场仿真尤其是变压器、电机这类带铁磁材料的模型碰到非线性迭代不收敛几乎是家常便饭。很多人以为把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_{k1} 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 时HG很好解当 τ1 时HF回到目标问题。τ 从 0 慢慢变到 1 的过程中我们每增加一点 τ就把前一时刻的解当作当前时刻的初值再去解 H(A,τ)0。这样一个参数一个参数地递推过去就是“延拓”。它相当于在山脚下找到一条缓坡路每一步都小到不会偏离路径最终走到山顶。在静磁场问题里这个思路对应非常自然的物理图景。最简单的同伦构型是负载步进H(A, τ) S(A) A - τ b 0τ0 时没有鼓励解是 A0τ 逐步增大相当于把绕组电流一点点加满。这个方案我先给结果好用直观而且绝大多数静磁场问题都可以这样解决。2.2 三种常用同伦构型的选型对比实际工程实现里我见过三种同伦构型用得最多。第一种就是上面说的负载步进。它适合电机、变压器、螺线管这类有明确电流激励的问题。物理意义清晰初始解零场求几乎不需要额外处理。缺点是如果材料本身是永磁体、或模型里存在剩磁没有外部激励时磁场也不为零τ0 的解并不是一个自然的起点。第二种是材料参数同伦。把实际非线性材料的磁阻率做一个带参数的插值ν_eff(τ) τ ν(B) (1-τ) ν_linearτ0 时全模型用线性磁导率解起来非常轻松τ 增加到 1 时材料逐渐变回真实非线性 B-H 曲线。这种做法的好处是不依赖外部激励对永磁结构一样适用。代价是中间状态没有物理意义纯粹是一种数值手段。第三种是经典的牛顿同伦也叫全局同伦H(A, τ) F(A) - (1-τ) F(A₀)其中 A₀ 是初始猜测。τ0 时H(A₀,0)0A₀ 天然就是一个解τ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 0b 节点由电流密度 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_end1.0, dt_init0.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_{k1} A_k α ΔAα ∈ (0,1]先从 α1 开始计算新残差的范数如果残差不下降就把 α 减半最多减五六次。伪代码如下def newton_at_fixed_t(A0, mesh, load, t, bh_model, max_iter30, tol_rel1e-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 逐步求解的过程记录我用上一节的代码跑了一遍记录下了每一步的收敛细节外步目标 τ本次步长 Δτ内层牛顿迭代次数铁芯最大磁密 B_max状态0.050.0530.61 T顺利0.100.0540.95 T顺利0.200.1061.42 T进入拐点区0.350.1591.63 T收敛变慢0.500.15——发散步长减半0.4250.07571.74 T重试成功0.550.12581.78 T顺利0.700.15101.81 T饱和收敛变慢0.800.1091.83 T顺利0.900.1081.84 T顺利1.000.1071.85 T最终收敛这个表很有意思。τ 从 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.01~0.05或改用材料参数同伦某步内层牛顿不下降大步跨过非线性拐点初值离开收敛域步长减半加入线搜索残差始终卡在同一量级雅可比矩阵缺 dν/dB 项检查修正项是否正确组装接近 τ1 时反复失败中间容差给太松误差沿路径累积把中间容差收紧到 1e-5末尾用 1e-6解出来磁密分布明显不对称网格剖分不对称或单元方向性检查网格质量对称问题用对称边界条件铁芯磁密超过 B-H 数据上限电流太大或材料模型太软扩展 B-H 数据或限制求解过程中磁密不超过上限这个表里的每一条我都踩过。最坑的是第一行有人看到 τ0 时方程几乎是线性就大胆把初始步长设成 0.3结果从零解跳到 30% 负载铁芯直接进入饱和区牛顿法第一步就把解带飞。所以初始步长宁可保守。第二行也很经典。某一步牛顿迭代失败不一定是算法坏了而是路径在这个区域曲率变大步长需要缩小。自适应步长正是为了这个场景设计的不要手动硬扛。5.2 两个容易被忽视的坑第一个坑只组装残差不重新组装雅可比。我在调试早期犯过一个错误为了省时间牛顿迭代里雅可比矩阵只在每个 τ 步开始时组装一次后面内层迭代一直复用。结果残差下降两个数量级后就再也不动了怎么增加迭代步数都没用。原因很简单非线性系统里雅可比矩阵也是解的函数必须每步更新。后来我看到代码里复用了旧雅可比立刻改成每次迭代重新组装收敛速度马上恢复正常。第二个坑对 B-H 曲线直接用分段线性插值。很多初学者从材料库拿到的是离散的 B-H 数据直接用线性分段连起来就丢给求解器。在强饱和区附近相邻数据点之间导数变化非常大甚至出现折线拐点。这会让牛顿法陷入一种“来回弹跳”的状态残差不再下降。我用三次样条或者 Akima 插值处理 B-H 数据同时保证导数连续这个问题就消失了。如果你必须保留原始数据点作为拟合依据也千万不要让样条出现过冲。B-H 曲线的过冲会产生伪负微分磁导率物理上完全不可接受求解结果也会特别诡异。6. 个人经验收尾同伦延拓是备胎还是主力用的次数多了我对同伦延拓的定位是它应该作为求解器里的“保底策略”而不是每个问题都主动上。对大多数常规静磁场模型一个好的阻尼牛顿法加连续光滑的 B-H 数据三四步就能收敛没必要引入延拓。但是一旦这些常规手段失败同伦延拓就是我最先调用的救场工具。我养成了一个调试习惯面对新模型时先跑一次同伦延拓拿到一个可靠的参考解然后再关掉延拓直接跑普通牛顿观察它在哪个环节出问题、收敛域边界在哪里。这个过程能帮你积累很多关于模型本身的认识。如果一个模型在同伦延拓的辅助下依然发散那我会立刻停止调算法转而检查建模问题——大概率是网格有负面积、气隙没加密、材料数据区域有空缺或者边界条件加载有误。最后分享一个小技巧把每个 τ 步的 B_max 和内层迭代次数打出来。这两个量是极好的“探路石”能清晰地告诉你整个求解路径上的难度分布。什么时候该放大步长、什么时候该收着走一看输出就知道。比盲试参数高效得多。这个方法我也用在别的非线性问题上比如超导磁体的场分布计算、磁滞损耗近似模型原理相通。非线性方程组的求解思路是通用的同伦延拓只是其中最可靠的一条路希望这次记录能帮你少走点弯路。 同伦延拓破解静磁场非线性仿真收敛难题
做静磁场仿真尤其是变压器、电机这类带铁磁材料的模型碰到非线性迭代不收敛几乎是家常便饭。很多人以为把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_{k1} 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 时HG很好解当 τ1 时HF回到目标问题。τ 从 0 慢慢变到 1 的过程中我们每增加一点 τ就把前一时刻的解当作当前时刻的初值再去解 H(A,τ)0。这样一个参数一个参数地递推过去就是“延拓”。它相当于在山脚下找到一条缓坡路每一步都小到不会偏离路径最终走到山顶。在静磁场问题里这个思路对应非常自然的物理图景。最简单的同伦构型是负载步进H(A, τ) S(A) A - τ b 0τ0 时没有鼓励解是 A0τ 逐步增大相当于把绕组电流一点点加满。这个方案我先给结果好用直观而且绝大多数静磁场问题都可以这样解决。2.2 三种常用同伦构型的选型对比实际工程实现里我见过三种同伦构型用得最多。第一种就是上面说的负载步进。它适合电机、变压器、螺线管这类有明确电流激励的问题。物理意义清晰初始解零场求几乎不需要额外处理。缺点是如果材料本身是永磁体、或模型里存在剩磁没有外部激励时磁场也不为零τ0 的解并不是一个自然的起点。第二种是材料参数同伦。把实际非线性材料的磁阻率做一个带参数的插值ν_eff(τ) τ ν(B) (1-τ) ν_linearτ0 时全模型用线性磁导率解起来非常轻松τ 增加到 1 时材料逐渐变回真实非线性 B-H 曲线。这种做法的好处是不依赖外部激励对永磁结构一样适用。代价是中间状态没有物理意义纯粹是一种数值手段。第三种是经典的牛顿同伦也叫全局同伦H(A, τ) F(A) - (1-τ) F(A₀)其中 A₀ 是初始猜测。τ0 时H(A₀,0)0A₀ 天然就是一个解τ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 0b 节点由电流密度 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_end1.0, dt_init0.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_{k1} A_k α ΔAα ∈ (0,1]先从 α1 开始计算新残差的范数如果残差不下降就把 α 减半最多减五六次。伪代码如下def newton_at_fixed_t(A0, mesh, load, t, bh_model, max_iter30, tol_rel1e-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 逐步求解的过程记录我用上一节的代码跑了一遍记录下了每一步的收敛细节外步目标 τ本次步长 Δτ内层牛顿迭代次数铁芯最大磁密 B_max状态0.050.0530.61 T顺利0.100.0540.95 T顺利0.200.1061.42 T进入拐点区0.350.1591.63 T收敛变慢0.500.15——发散步长减半0.4250.07571.74 T重试成功0.550.12581.78 T顺利0.700.15101.81 T饱和收敛变慢0.800.1091.83 T顺利0.900.1081.84 T顺利1.000.1071.85 T最终收敛这个表很有意思。τ 从 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.01~0.05或改用材料参数同伦某步内层牛顿不下降大步跨过非线性拐点初值离开收敛域步长减半加入线搜索残差始终卡在同一量级雅可比矩阵缺 dν/dB 项检查修正项是否正确组装接近 τ1 时反复失败中间容差给太松误差沿路径累积把中间容差收紧到 1e-5末尾用 1e-6解出来磁密分布明显不对称网格剖分不对称或单元方向性检查网格质量对称问题用对称边界条件铁芯磁密超过 B-H 数据上限电流太大或材料模型太软扩展 B-H 数据或限制求解过程中磁密不超过上限这个表里的每一条我都踩过。最坑的是第一行有人看到 τ0 时方程几乎是线性就大胆把初始步长设成 0.3结果从零解跳到 30% 负载铁芯直接进入饱和区牛顿法第一步就把解带飞。所以初始步长宁可保守。第二行也很经典。某一步牛顿迭代失败不一定是算法坏了而是路径在这个区域曲率变大步长需要缩小。自适应步长正是为了这个场景设计的不要手动硬扛。5.2 两个容易被忽视的坑第一个坑只组装残差不重新组装雅可比。我在调试早期犯过一个错误为了省时间牛顿迭代里雅可比矩阵只在每个 τ 步开始时组装一次后面内层迭代一直复用。结果残差下降两个数量级后就再也不动了怎么增加迭代步数都没用。原因很简单非线性系统里雅可比矩阵也是解的函数必须每步更新。后来我看到代码里复用了旧雅可比立刻改成每次迭代重新组装收敛速度马上恢复正常。第二个坑对 B-H 曲线直接用分段线性插值。很多初学者从材料库拿到的是离散的 B-H 数据直接用线性分段连起来就丢给求解器。在强饱和区附近相邻数据点之间导数变化非常大甚至出现折线拐点。这会让牛顿法陷入一种“来回弹跳”的状态残差不再下降。我用三次样条或者 Akima 插值处理 B-H 数据同时保证导数连续这个问题就消失了。如果你必须保留原始数据点作为拟合依据也千万不要让样条出现过冲。B-H 曲线的过冲会产生伪负微分磁导率物理上完全不可接受求解结果也会特别诡异。6. 个人经验收尾同伦延拓是备胎还是主力用的次数多了我对同伦延拓的定位是它应该作为求解器里的“保底策略”而不是每个问题都主动上。对大多数常规静磁场模型一个好的阻尼牛顿法加连续光滑的 B-H 数据三四步就能收敛没必要引入延拓。但是一旦这些常规手段失败同伦延拓就是我最先调用的救场工具。我养成了一个调试习惯面对新模型时先跑一次同伦延拓拿到一个可靠的参考解然后再关掉延拓直接跑普通牛顿观察它在哪个环节出问题、收敛域边界在哪里。这个过程能帮你积累很多关于模型本身的认识。如果一个模型在同伦延拓的辅助下依然发散那我会立刻停止调算法转而检查建模问题——大概率是网格有负面积、气隙没加密、材料数据区域有空缺或者边界条件加载有误。最后分享一个小技巧把每个 τ 步的 B_max 和内层迭代次数打出来。这两个量是极好的“探路石”能清晰地告诉你整个求解路径上的难度分布。什么时候该放大步长、什么时候该收着走一看输出就知道。比盲试参数高效得多。这个方法我也用在别的非线性问题上比如超导磁体的场分布计算、磁滞损耗近似模型原理相通。非线性方程组的求解思路是通用的同伦延拓只是其中最可靠的一条路希望这次记录能帮你少走点弯路。