静磁场仿真-主题004_磁路理论与等效电路-磁路理论与等效电路

kkchenjj

283人浏览 · 2026-04-06 21:32:05

kkchenjj · 2026-04-06 21:32:05 发布

第4篇：磁路理论与等效电路

摘要

磁路理论是电磁学工程应用中的核心基础理论，它将复杂的磁场分布问题简化为类似于电路分析的磁路模型，极大地降低了电机、变压器、电磁铁等电磁装置的设计难度。本教程系统讲解磁路的基本理论，包括磁路欧姆定律、磁阻计算、磁动势概念、磁路串并联分析方法，以及气隙磁阻的计算方法。通过Python数值仿真，深入分析磁芯线圈的电感计算、磁路饱和效应，并建立磁路与电路的类比关系。本教程提供完整的Python仿真代码，生成6-10张高质量仿真图，帮助读者直观理解磁路的工作原理，掌握磁路等效电路的建模与分析方法，为电磁装置的设计与优化奠定坚实基础。

关键词

磁路理论、磁路欧姆定律、磁阻、磁动势、等效电路、气隙磁阻、磁饱和、电感计算、磁路仿真

在这里插入图片描述

一、教程目标

1.1 学习目标

完成本教程的学习后，您将能够：

理解磁路基本概念：掌握磁路、磁阻、磁动势、磁通等核心概念及其物理意义
应用磁路欧姆定律：熟练运用 $\Phi R_m$ 进行磁路计算
计算磁路参数：能够计算各种磁路结构的磁阻，包括带气隙磁路
分析磁路串并联：掌握复杂磁路的等效变换与分析方法
计算线圈电感：理解电感与磁路参数的关系，能够计算带磁芯线圈的电感
分析磁饱和效应：理解磁饱和现象，掌握饱和磁路的分析方法
建立磁路-电路类比：理解磁路与电路的对应关系，建立等效电路模型

1.2 应用场景

本教程的知识可应用于以下工程领域：

电机设计：计算电机磁路的磁通分布、磁阻和电感参数
变压器设计：分析变压器铁芯的磁化特性、励磁电流
电磁铁设计：计算电磁铁的吸力、磁通密度分布
电感器设计：优化电感器的磁芯结构，提高电感量
传感器设计：分析磁路型传感器的灵敏度与线性度
电磁兼容：评估设备的磁泄漏与屏蔽效果

二、前置知识与环境

2.1 前置知识

学习本教程前，建议具备以下基础知识：

电磁学基础：
- 磁场的基本概念（磁感应强度 $B$ 、磁场强度 $H$ ）
- 安培环路定律
- 磁介质的磁化特性
- 磁导率 $μ\mu$ 的概念
电路分析基础：
- 欧姆定律
- 电阻串并联
- 基尔霍夫定律
- 等效电路变换
数学基础：
- 微积分基本运算
- 线性代数基础
- 常微分方程基础

2.2 软件环境

本教程使用以下软件环境：

Python 3.8+：编程语言
NumPy：数值计算库
Matplotlib：数据可视化库
SciPy：科学计算库（用于插值和优化）

2.3 安装依赖

pip install numpy matplotlib scipy

三、原理极简讲解

3.1 磁路的基本概念

3.1.1 什么是磁路

磁路是指磁通所经过的闭合路径。与电路中电流沿导体流动类似，磁通主要在磁导率高的材料（如铁磁材料）中通过。磁路理论的核心思想是将磁场问题类比为电路问题，从而简化分析。

磁路的典型结构包括：

磁芯：提供低磁阻通路（如硅钢片、铁氧体）
气隙：磁路中的高磁阻部分
励磁线圈：产生磁动势的源

3.1.2 磁路欧姆定律

磁路欧姆定律是磁路理论的基础，其数学表达式为：

$\Phi \cdot R_m$

其中：

$F$ ：磁动势（Magnetomotive Force, MMF），单位：安匝（At）
$Φ\Phi$ ：磁通量（Magnetic Flux），单位：韦伯（Wb）
$R_m$ ：磁阻（Reluctance），单位：安匝/韦伯（At/Wb）或亨 $^{-1}$ （ $H^{-1}$ ）

这与电路欧姆定律 $\cdot R$ 形式完全对应：

电路	磁路
电动势 $V$	磁动势 $F$
电流 $I$	磁通 $Φ\Phi$
电阻 $R$	磁阻 $R_m$
电导 $G = 1/ R$	磁导 $Λ=1/Rm\Lambda = 1/R_m$

3.2 磁动势

磁动势是产生磁通的源，由通过线圈的电流产生：

$\cdot I$

其中：

$N$ ：线圈匝数
$I$ ：线圈电流（A）

磁动势的单位是安匝（Ampere-Turns, At），表示线圈产生磁场的能力。

3.3 磁阻的计算

3.3.1 均匀磁路的磁阻

对于截面均匀、材料均匀的磁路，磁阻计算公式为：

$Rm=lμA=lμ0μrAR_m = \frac{l}{\mu A} = \frac{l}{\mu_0 \mu_r A}$

其中：

$l$ ：磁路长度（m）
$A$ ：磁路截面积（m²）
$μ\mu$ ：材料磁导率（H/m）
$μ0\mu_0$ ：真空磁导率， $μ0=4π×10−7\mu_0 = 4\pi \times 10^{-7}$ H/m
$μr\mu_r$ ：相对磁导率（无量纲）

3.3.2 气隙磁阻

气隙是磁路中的重要组成部分，其磁阻计算需要考虑边缘效应。

忽略边缘效应时：

$Rg=gμ0AR_{g} = \frac{g}{\mu_0 A}$

其中 $g$ 为气隙长度。

考虑边缘效应时：

气隙边缘的磁力线会向外扩散，等效增大了气隙面积。经验公式为：

$A_{eff} = (a + g)(b + g) = ab + g(a + b) + g^2$

对于圆形截面：

$Aeff=π(r+g2)2A_{eff} = \pi \left(r + \frac{g}{2}\right)^2$

修正后的气隙磁阻：

$Rg,eff=gμ0AeffR_{g,eff} = \frac{g}{\mu_0 A_{eff}}$

3.4 磁路的串并联

3.4.1 串联磁路

多个磁路段串联时，总磁阻等于各段磁阻之和：

$Rm,total=Rm1+Rm2+⋯+RmnR_{m,total} = R_{m1} + R_{m2} + \cdots + R_{mn}$

串联磁路中，各段的磁通相同，磁动势在各段上分配：

$\Phi(R_{m1} + R_{m2} + \cdots) = \Phi R_{m1} + \Phi R_{m2} + \cdots = F_1 + F_2 + \cdots$

3.4.2 并联磁路

多个磁路段并联时，总磁阻的倒数等于各段磁阻倒数之和：

$1Rm,total=1Rm1+1Rm2+⋯+1Rmn\frac{1}{R_{m,total}} = \frac{1}{R_{m1}} + \frac{1}{R_{m2}} + \cdots + \frac{1}{R_{mn}}$

并联磁路中，各支路的磁动势相同，磁通在各支路中分配。

3.5 电感的磁路计算

3.5.1 电感与磁路的关系

线圈的电感定义为磁链与电流之比：

$\frac{\Psi}{I} = \frac{N\Phi}{I}$

根据磁路欧姆定律 $Φ=F/Rm=NI/Rm\Phi = F/R_m = NI/R_m$ ，可得：

$\frac{N^2}{R_m} = N^2 \Lambda$

其中 $Λ=1/Rm\Lambda = 1/R_m$ 为磁导。

3.5.2 带气隙磁芯的电感

对于带气隙的磁芯，总磁阻为铁芯磁阻与气隙磁阻之和：

$Rm=Rcore+Rg=lcoreμ0μrA+gμ0AeffR_m = R_{core} + R_g = \frac{l_{core}}{\mu_0 \mu_r A} + \frac{g}{\mu_0 A_{eff}}$

电感为：

$\frac{N^2}{R_{core} + R_g}$

当 $μr≫1\mu_r \gg 1$ 时， $Rcore≪RgR_{core} \ll R_g$ ，电感主要由气隙决定：

$\approx \frac{N^2 \mu_0 A_{eff}}{g}$

3.6 磁饱和效应

3.6.1 B-H曲线

铁磁材料的磁化特性用B-H曲线描述：

$\mu_0 \mu_r(H) \cdot H$

其中 $μr(H)\mu_r(H)$ 是磁场强度的函数，随H增大而减小。

3.6.2 饱和对磁路的影响

当磁路饱和时：

磁导率 $μ\mu$ 下降
磁阻 $R_m$ 增大
电感 $L$ 减小
磁通 $Φ\Phi$ 增长变缓

饱和磁路的分析需要采用迭代法或图解法。

3.7 磁路-电路类比总结

电路量	符号	单位	磁路量	符号	单位
电动势	$V$	V	磁动势	$F$	At
电流	$I$	A	磁通	$Φ\Phi$	Wb
电阻	$R$	Ω	磁阻	$R_m$	At/Wb
电导	$G$	S	磁导	$Λ\Lambda$	Wb/At
电流密度	$J$	A/m²	磁通密度	$B$	T
电场强度	$E$	V/m	磁场强度	$H$	At/m
欧姆定律	$V = I R$	-	磁路定律	$F=ΦRmF=\Phi R_m$	-
基尔霍夫电流定律	$∑I=0\sum I = 0$	-	磁通连续性	$∑Φ=0\sum \Phi = 0$	-
基尔霍夫电压定律	$∑V=0\sum V = 0$	-	安培环路定律	$∑F=∑NI\sum F = \sum NI$	-

四、完整步骤（Step by Step）

4.1 仿真案例设计

本教程设计以下仿真案例：

案例1：简单磁路分析 - 验证磁路欧姆定律
案例2：磁路串联分析 - 含气隙磁路的计算
案例3：磁路并联分析 - 分支磁路的磁通分配
案例4：电感计算 - 带磁芯线圈的电感特性
案例5：磁饱和效应 - 非线性磁路的分析
案例6：磁路-电路类比可视化 - 对比展示

4.2 仿真步骤

步骤1：定义磁路基本参数

首先定义磁路的几何参数和材料参数：

# 磁路几何参数
l_core = 0.2      # 磁芯长度 (m)
A_core = 1e-4     # 磁芯截面积 (m²)
g = 1e-3          # 气隙长度 (m)
mu_r = 5000       # 相对磁导率
N = 100           # 线圈匝数
I = 1.0           # 电流 (A)

步骤2：计算磁阻

根据磁阻公式计算各段磁阻：

# 真空磁导率
mu_0 = 4 * np.pi * 1e-7  # H/m

# 磁芯磁阻
R_core = l_core / (mu_0 * mu_r * A_core)

# 气隙磁阻（考虑边缘效应）
# 假设磁芯截面为正方形，边长为a
a = np.sqrt(A_core)
A_gap_eff = (a + g)**2  # 考虑边缘效应的等效面积
R_gap = g / (mu_0 * A_gap_eff)

步骤3：计算磁动势和磁通

# 磁动势
F = N * I

# 总磁阻
R_total = R_core + R_gap

# 磁通
Phi = F / R_total

步骤4：计算磁通密度和磁场强度

# 磁通密度
B = Phi / A_core

# 磁芯中的磁场强度
H_core = B / (mu_0 * mu_r)

# 气隙中的磁场强度
H_gap = B / mu_0

步骤5：计算电感

# 电感
L = N**2 / R_total

步骤6：可视化结果

使用Matplotlib绘制磁路参数分布图、B-H曲线等。

五、Python仿真代码

以下是完整的Python仿真代码，包含详细的注释说明：

# -*- coding: utf-8 -*-
"""
Topic 004: Magnetic Circuit Theory and Equivalent Circuit Simulation
====================================================================
This program implements basic magnetic circuit theory simulation:
1. Magnetic circuit Ohm's law verification
2. Magnetic circuit series/parallel analysis
3. Air gap reluctance calculation
4. Inductance calculation
5. Magnetic saturation effect analysis
6. Magnetic-electric circuit analogy visualization

Author: Simulation Tutorial
Date: 2026
"""

import matplotlib
matplotlib.use('Agg')  # Use non-interactive backend
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import FancyBboxPatch, Circle, Rectangle, FancyArrowPatch
import warnings
warnings.filterwarnings('ignore')
import os

# Set font support
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans', 'Arial Unicode MS']
plt.rcParams['axes.unicode_minus'] = False

# Create output directory
output_dir = r'd:\文档\500仿真领域\工程仿真\静磁场仿真\主题004_磁路理论与等效电路'
os.makedirs(output_dir, exist_ok=True)

print("=" * 60)
print("Topic 004: Magnetic Circuit Theory and Equivalent Circuit Simulation")
print("=" * 60)

# ============================================================================
# Part 1: Basic Parameter Definition
# ============================================================================
print("\n[1] Defining basic physical constants and material parameters...")

# Vacuum permeability
mu_0 = 4 * np.pi * 1e-7  # H/m

# Material parameters
mu_r_iron = 5000         # Relative permeability of iron core
mu_r_ferrite = 2000      # Relative permeability of ferrite
mu_r_air = 1             # Relative permeability of air

# Geometric parameters (typical EE core dimensions)
l_core = 0.1             # Core magnetic path length (m)
A_core = 2e-4            # Core cross-sectional area (m^2), 20mm x 10mm
g = 1e-3                 # Air gap length (m), 1mm

# Coil parameters
N = 100                  # Number of turns
I_range = np.linspace(0.01, 5, 100)  # Current range (A)

print(f"  Vacuum permeability mu_0 = {mu_0:.4e} H/m")
print(f"  Iron core relative permeability mu_r = {mu_r_iron}")
print(f"  Core length l = {l_core*1000:.1f} mm")
print(f"  Core cross-sectional area A = {A_core*1e6:.1f} mm^2")
print(f"  Air gap length g = {g*1000:.1f} mm")
print(f"  Number of turns N = {N}")

# ============================================================================
# Part 2: Magnetic Circuit Ohm's Law Simulation
# ============================================================================
print("\n[2] Magnetic Circuit Ohm's Law simulation...")

# Calculate reluctance
R_core = l_core / (mu_0 * mu_r_iron * A_core)  # Core reluctance
R_gap = g / (mu_0 * A_core)                      # Air gap reluctance (no fringing)
R_total = R_core + R_gap                         # Total reluctance

# Permeance
Lambda_core = 1 / R_core
Lambda_gap = 1 / R_gap
Lambda_total = 1 / R_total

print(f"  Core reluctance R_core = {R_core:.4e} At/Wb")
print(f"  Air gap reluctance R_gap = {R_gap:.4e} At/Wb")
print(f"  Total reluctance R_total = {R_total:.4e} At/Wb")
print(f"  Core permeance Lambda_core = {Lambda_core:.4e} Wb/At")
print(f"  Air gap permeance Lambda_gap = {Lambda_gap:.4e} Wb/At")

# Calculate flux for different currents
F_range = N * I_range  # MMF range
Phi_range = F_range / R_total  # Magnetic flux
B_range = Phi_range / A_core   # Flux density

print(f"  At I = 1A:")
print(f"    MMF F = {N*1:.1f} At")
print(f"    Flux Phi = {N*1/R_total:.4e} Wb")
print(f"    Flux density B = {N*1/(R_total*A_core):.4f} T")

# ============================================================================
# Figure 1: Magnetic Circuit Ohm's Law Verification
# ============================================================================
fig, axes = plt.subplots(2, 2, figsize=(12, 10))
fig.suptitle('Magnetic Circuit Ohm\'s Law Verification', fontsize=14, fontweight='bold')

# 子图1：磁动势-磁通关系
ax1 = axes[0, 0]
ax1.plot(F_range, Phi_range * 1e4, 'b-', linewidth=2)
ax1.set_xlabel('Magnetomotive Force F (At)', fontsize=11)
ax1.set_ylabel('Magnetic Flux Φ (10⁻⁴ Wb)', fontsize=11)
ax1.set_title('F-Φ Characteristic (Ohm\'s Law)', fontsize=12)
ax1.grid(True, alpha=0.3)
ax1.set_xlim(0, max(F_range))
ax1.set_ylim(0, max(Phi_range) * 1e4 * 1.1)

# 子图2：电流-磁通密度关系
ax2 = axes[0, 1]
ax2.plot(I_range, B_range, 'r-', linewidth=2)
ax2.set_xlabel('Current I (A)', fontsize=11)
ax2.set_ylabel('Flux Density B (T)', fontsize=11)
ax2.set_title('I-B Characteristic', fontsize=12)
ax2.grid(True, alpha=0.3)
ax2.set_xlim(0, max(I_range))
ax2.set_ylim(0, max(B_range) * 1.1)

# 子图3：磁阻分布（饼图）
ax3 = axes[1, 0]
labels = ['Core\nReluctance', 'Gap\nReluctance']
sizes = [R_core/R_total*100, R_gap/R_total*100]
colors = ['#ff9999', '#66b3ff']
explode = (0.05, 0.05)
ax3.pie(sizes, explode=explode, labels=labels, colors=colors, autopct='%1.1f%%',
        shadow=True, startangle=90)
ax3.set_title('Reluctance Distribution', fontsize=12)

# 子图4：磁导对比（柱状图）
ax4 = axes[1, 1]
categories = ['Core\nPermeance', 'Gap\nPermeance', 'Total\nPermeance']
values = [Lambda_core*1e6, Lambda_gap*1e6, Lambda_total*1e6]
bars = ax4.bar(categories, values, color=['#ff9999', '#66b3ff', '#99ff99'], edgecolor='black')
ax4.set_ylabel('Permeance (μWb/At)', fontsize=11)
ax4.set_title('Permeance Comparison', fontsize=12)
ax4.grid(True, alpha=0.3, axis='y')
# 添加数值标签
for bar, val in zip(bars, values):
    height = bar.get_height()
    ax4.text(bar.get_x() + bar.get_width()/2., height,
             f'{val:.2f}', ha='center', va='bottom', fontsize=10)

plt.tight_layout()
plt.savefig(f'{output_dir}/fig1_ohms_law.png', dpi=150, bbox_inches='tight')
plt.close()
print("  [OK] Figure 1 saved: fig1_ohms_law.png")

# ============================================================================
# Part 3: Air Gap Fringing Effect Analysis
# ============================================================================
print("\n[3] Air gap fringing effect analysis...")

# Air gap length range
g_range = np.linspace(0.1e-3, 5e-3, 100)  # 0.1mm to 5mm

# Assume square cross-section for core
a = np.sqrt(A_core)  # Side length

# Air gap reluctance without fringing
R_gap_simple = g_range / (mu_0 * A_core)

# Air gap reluctance with fringing (effective area method)
A_gap_eff = (a + g_range)**2  # Effective area
R_gap_edge = g_range / (mu_0 * A_gap_eff)

# Fringing effect coefficient
edge_factor = R_gap_simple / R_gap_edge

print(f"  Air gap length range: {g_range[0]*1000:.1f} mm ~ {g_range[-1]*1000:.1f} mm")
print(f"  Fringing effect factor range: {edge_factor[0]:.3f} ~ {edge_factor[-1]:.3f}")

# ============================================================================
# 图2：气隙边缘效应
# ============================================================================
fig, axes = plt.subplots(1, 3, figsize=(15, 5))
fig.suptitle('Air Gap Fringing Effect Analysis', fontsize=14, fontweight='bold')

# 子图1：气隙磁阻对比
ax1 = axes[0]
ax1.plot(g_range*1000, R_gap_simple, 'b-', linewidth=2, label='Without Fringing')
ax1.plot(g_range*1000, R_gap_edge, 'r--', linewidth=2, label='With Fringing')
ax1.set_xlabel('Gap Length g (mm)', fontsize=11)
ax1.set_ylabel('Reluctance (At/Wb)', fontsize=11)
ax1.set_title('Gap Reluctance vs Length', fontsize=12)
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3)

# 子图2：边缘效应系数
ax2 = axes[1]
ax2.plot(g_range*1000, edge_factor, 'g-', linewidth=2)
ax2.set_xlabel('Gap Length g (mm)', fontsize=11)
ax2.set_ylabel('Edge Factor', fontsize=11)
ax2.set_title('Fringing Effect Factor', fontsize=12)
ax2.grid(True, alpha=0.3)
ax2.axhline(y=1, color='k', linestyle='--', alpha=0.5)

# 子图3：等效面积增加比例
ax3 = axes[2]
area_ratio = A_gap_eff / A_core
ax3.plot(g_range*1000, area_ratio, 'm-', linewidth=2)
ax3.set_xlabel('Gap Length g (mm)', fontsize=11)
ax3.set_ylabel('A_eff / A_core', fontsize=11)
ax3.set_title('Effective Area Increase', fontsize=12)
ax3.grid(True, alpha=0.3)
ax3.axhline(y=1, color='k', linestyle='--', alpha=0.5)

plt.tight_layout()
plt.savefig(f'{output_dir}/fig2_fringing_effect.png', dpi=150, bbox_inches='tight')
plt.close()
print("  [OK] Figure 2 saved: fig2_fringing_effect.png")

# ============================================================================
# Part 4: Magnetic Circuit Series and Parallel Analysis
# ============================================================================
print("\n[4] Magnetic circuit series and parallel analysis...")

# Define circuit parameters
# Series circuit: core + gap
R_series = R_core + R_gap

# Parallel circuit: two identical cores in parallel
R_parallel = R_core / 2

# Mixed circuit: two cores in parallel, then in series with gap
R_mixed = R_gap + R_core / 2

print(f"  Series reluctance: {R_series:.4e} At/Wb")
print(f"  Parallel reluctance: {R_parallel:.4e} At/Wb")
print(f"  Mixed reluctance: {R_mixed:.4e} At/Wb")

# 计算不同磁路结构的磁通
I_test = np.linspace(0.1, 3, 50)
F_test = N * I_test

Phi_series = F_test / R_series
Phi_parallel = F_test / R_parallel
Phi_mixed = F_test / R_mixed

# ============================================================================
# 图3：磁路串并联对比
# ============================================================================
fig, axes = plt.subplots(2, 2, figsize=(12, 10))
fig.suptitle('Magnetic Circuit Series and Parallel Analysis', fontsize=14, fontweight='bold')

# 子图1：不同磁路结构的磁通-电流特性
ax1 = axes[0, 0]
ax1.plot(I_test, Phi_series*1e4, 'b-', linewidth=2, label='Series (Core+Gap)')
ax1.plot(I_test, Phi_parallel*1e4, 'r-', linewidth=2, label='Parallel (2 Cores)')
ax1.plot(I_test, Phi_mixed*1e4, 'g-', linewidth=2, label='Mixed (Gap+2 Cores)')
ax1.set_xlabel('Current I (A)', fontsize=11)
ax1.set_ylabel('Magnetic Flux Φ (10⁻⁴ Wb)', fontsize=11)
ax1.set_title('Flux vs Current for Different Circuits', fontsize=12)
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3)

# 子图2：磁阻对比（柱状图）
ax2 = axes[0, 1]
circuits = ['Series\n(Core+Gap)', 'Parallel\n(2 Cores)', 'Mixed\n(Gap+2 Cores)']
reluctances = [R_series*1e-6, R_parallel*1e-6, R_mixed*1e-6]
colors = ['#ff9999', '#66b3ff', '#99ff99']
bars = ax2.bar(circuits, reluctances, color=colors, edgecolor='black')
ax2.set_ylabel('Reluctance (10⁶ At/Wb)', fontsize=11)
ax2.set_title('Reluctance Comparison', fontsize=12)
ax2.grid(True, alpha=0.3, axis='y')
for bar, val in zip(bars, reluctances):
    height = bar.get_height()
    ax2.text(bar.get_x() + bar.get_width()/2., height,
             f'{val:.2f}', ha='center', va='bottom', fontsize=10)

# 子图3：磁导对比
ax3 = axes[1, 0]
permeances = [1/R_series*1e6, 1/R_parallel*1e6, 1/R_mixed*1e6]
bars = ax3.bar(circuits, permeances, color=colors, edgecolor='black')
ax3.set_ylabel('Permeance (μWb/At)', fontsize=11)
ax3.set_title('Permeance Comparison', fontsize=12)
ax3.grid(True, alpha=0.3, axis='y')
for bar, val in zip(bars, permeances):
    height = bar.get_height()
    ax3.text(bar.get_x() + bar.get_width()/2., height,
             f'{val:.2f}', ha='center', va='bottom', fontsize=10)

# 子图4：磁通密度对比（在I=1A时）
ax4 = axes[1, 1]
B_series = (N*1/R_series) / A_core
B_parallel = (N*1/R_parallel) / A_core
B_mixed = (N*1/R_mixed) / A_core
B_values = [B_series, B_parallel, B_mixed]
bars = ax4.bar(circuits, B_values, color=colors, edgecolor='black')
ax4.set_ylabel('Flux Density B (T)', fontsize=11)
ax4.set_title(f'Flux Density at I=1A', fontsize=12)
ax4.grid(True, alpha=0.3, axis='y')
for bar, val in zip(bars, B_values):
    height = bar.get_height()
    ax4.text(bar.get_x() + bar.get_width()/2., height,
             f'{val:.3f}', ha='center', va='bottom', fontsize=10)

plt.tight_layout()
plt.savefig(f'{output_dir}/fig3_series_parallel.png', dpi=150, bbox_inches='tight')
plt.close()
print("  [OK] Figure 3 saved: fig3_series_parallel.png")

# ============================================================================
# Part 5: Inductance Calculation and Characteristics
# ============================================================================
print("\n[5] Inductance calculation and characteristics...")

# Calculate inductance for different turns
N_range = np.arange(10, 510, 10)  # 10 to 500 turns

# Different air gap lengths
g_values = [0.5e-3, 1e-3, 2e-3, 5e-3]  # 0.5mm, 1mm, 2mm, 5mm
g_labels = ['0.5mm', '1mm', '2mm', '5mm']
colors_gap = ['#ff6b6b', '#4ecdc4', '#45b7d1', '#96ceb4']

# Store inductance values
L_results = {}

for g_val, label in zip(g_values, g_labels):
    R_gap_val = g_val / (mu_0 * A_core)
    R_total_val = R_core + R_gap_val
    L_val = N_range**2 / R_total_val
    L_results[label] = L_val
    print(f"  Gap {label}: L(N=100) = {100**2/R_total_val*1e3:.2f} mH")

# ============================================================================
# 图4：电感特性分析
# ============================================================================
fig, axes = plt.subplots(2, 2, figsize=(12, 10))
fig.suptitle('Inductance Calculation and Characteristics', fontsize=14, fontweight='bold')

# 子图1：电感与匝数的关系
ax1 = axes[0, 0]
for label, color in zip(g_labels, colors_gap):
    ax1.plot(N_range, L_results[label]*1e3, '-', linewidth=2, 
             color=color, label=f'Gap={label}')
ax1.set_xlabel('Number of Turns N', fontsize=11)
ax1.set_ylabel('Inductance L (mH)', fontsize=11)
ax1.set_title('Inductance vs Turns', fontsize=12)
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3)

# 子图2：电感与气隙的关系（N=100）
ax2 = axes[0, 1]
g_range_L = np.linspace(0.1e-3, 10e-3, 100)
R_gap_range = g_range_L / (mu_0 * A_core)
R_total_range = R_core + R_gap_range
L_range = N**2 / R_total_range
ax2.plot(g_range_L*1000, L_range*1e3, 'b-', linewidth=2)
ax2.set_xlabel('Gap Length g (mm)', fontsize=11)
ax2.set_ylabel('Inductance L (mH)', fontsize=11)
ax2.set_title(f'Inductance vs Gap (N={N})', fontsize=12)
ax2.grid(True, alpha=0.3)

# 子图3：电感与匝数的平方关系验证
ax3 = axes[1, 0]
L_1mm = L_results['1mm']
ax3.plot(N_range**2, L_1mm*1e3, 'bo', markersize=4, label='Calculated')
# 理论直线
ax3.plot(N_range**2, (N_range**2/R_total)*1e3, 'r--', linewidth=2, label='Theory: L=N²/R')
ax3.set_xlabel('N² (turns²)', fontsize=11)
ax3.set_ylabel('Inductance L (mH)', fontsize=11)
ax3.set_title('L vs N² (Linear Relationship)', fontsize=12)
ax3.legend(fontsize=10)
ax3.grid(True, alpha=0.3)

# 子图4：储能对比
ax4 = axes[1, 1]
I_fixed = 1.0  # 固定电流
energy = {}
for label in g_labels:
    energy[label] = 0.5 * L_results[label] * I_fixed**2
    
x_pos = np.arange(len(g_labels))
bars = ax4.bar(x_pos, [energy[label][9]*1e3 for label in g_labels], 
               color=colors_gap, edgecolor='black')
ax4.set_xticks(x_pos)
ax4.set_xticklabels([f'g={l}' for l in g_labels])
ax4.set_ylabel('Energy W (mJ)', fontsize=11)
ax4.set_title(f'Energy Storage at I={I_fixed}A (N=100)', fontsize=12)
ax4.grid(True, alpha=0.3, axis='y')
for bar, label in zip(bars, g_labels):
    height = bar.get_height()
    ax4.text(bar.get_x() + bar.get_width()/2., height,
             f'{height:.2f}', ha='center', va='bottom', fontsize=10)

plt.tight_layout()
plt.savefig(f'{output_dir}/fig4_inductance.png', dpi=150, bbox_inches='tight')
plt.close()
print("  [OK] Figure 4 saved: fig4_inductance.png")

# ============================================================================
# Part 6: Magnetic Saturation Effect Analysis
# ============================================================================
print("\n[6] Magnetic saturation effect analysis...")

# Define B-H curve (with saturation)
# Use hyperbolic tangent to simulate saturation
def B_H_curve(H, B_sat=1.5, mu_i=5000*mu_0):
    """
    Simulate B-H curve for ferromagnetic material
    H: Magnetic field strength (At/m)
    B_sat: Saturation flux density (T)
    mu_i: Initial permeability (H/m)
    """
    # Use modified Langevin function for saturation
    a = B_sat / (mu_i * H + 1e-10)  # Avoid division by zero
    B = B_sat * np.tanh(H / (B_sat/mu_i))
    return B

# Generate H range
H_range = np.linspace(0, 10000, 500)  # 0 to 10000 At/m

# Calculate B values
B_nonlinear = B_H_curve(H_range, B_sat=1.5, mu_i=mu_r_iron*mu_0)

# Linear approximation (no saturation)
B_linear = mu_0 * mu_r_iron * H_range

# Calculate effective permeability
mu_eff = B_nonlinear / (H_range + 1e-10)

# Calculate flux for different currents (with saturation)
I_sat = np.linspace(0.01, 10, 200)
F_sat = N * I_sat

# Average magnetic field strength in core
H_core_avg = F_sat / l_core

# Nonlinear flux density
B_core_sat = B_H_curve(H_core_avg, B_sat=1.5, mu_i=mu_r_iron*mu_0)
Phi_sat = B_core_sat * A_core

# Linear approximation flux
Phi_linear = F_sat / R_total

print(f"  Saturation flux density B_sat = 1.5 T")
print(f"  Initial relative permeability mu_r = {mu_r_iron}")

# ============================================================================
# 图5：磁饱和效应
# ============================================================================
fig, axes = plt.subplots(2, 2, figsize=(12, 10))
fig.suptitle('Magnetic Saturation Effect Analysis', fontsize=14, fontweight='bold')

# 子图1：B-H曲线
ax1 = axes[0, 0]
ax1.plot(H_range, B_nonlinear, 'b-', linewidth=2, label='Nonlinear (Saturation)')
ax1.plot(H_range, B_linear, 'r--', linewidth=1.5, alpha=0.7, label='Linear (No Saturation)')
ax1.axhline(y=1.5, color='g', linestyle='-.', alpha=0.7, label='B_sat = 1.5T')
ax1.set_xlabel('Magnetic Field H (At/m)', fontsize=11)
ax1.set_ylabel('Flux Density B (T)', fontsize=11)
ax1.set_title('B-H Curve with Saturation', fontsize=12)
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3)
ax1.set_xlim(0, max(H_range))
ax1.set_ylim(0, 2)

# 子图2：有效磁导率
ax2 = axes[0, 1]
ax2.semilogy(H_range[1:], mu_eff[1:]/mu_0, 'm-', linewidth=2)
ax2.axhline(y=mu_r_iron, color='r', linestyle='--', alpha=0.7, label=f'Initial μr={mu_r_iron}')
ax2.set_xlabel('Magnetic Field H (At/m)', fontsize=11)
ax2.set_ylabel('Relative Permeability μr', fontsize=11)
ax2.set_title('Effective Permeability vs H', fontsize=12)
ax2.legend(fontsize=10)
ax2.grid(True, alpha=0.3)
ax2.set_xlim(0, max(H_range))

# 子图3：电流-磁通特性（线性与非线性对比）
ax3 = axes[1, 0]
ax3.plot(I_sat, Phi_sat*1e4, 'b-', linewidth=2, label='With Saturation')
ax3.plot(I_sat, Phi_linear*1e4, 'r--', linewidth=1.5, alpha=0.7, label='Linear (No Saturation)')
ax3.set_xlabel('Current I (A)', fontsize=11)
ax3.set_ylabel('Magnetic Flux Φ (10⁻⁴ Wb)', fontsize=11)
ax3.set_title('Flux vs Current (Saturation Effect)', fontsize=12)
ax3.legend(fontsize=10)
ax3.grid(True, alpha=0.3)

# 子图4：电感随电流变化（饱和导致电感下降）
ax4 = axes[1, 1]
L_sat = N * Phi_sat / I_sat  # 动态电感
L_const = N * Phi_linear / I_sat  # 线性电感
ax4.plot(I_sat, L_sat*1e3, 'b-', linewidth=2, label='With Saturation')
ax4.plot(I_sat, L_const*1e3, 'r--', linewidth=1.5, alpha=0.7, label='Linear')
ax4.set_xlabel('Current I (A)', fontsize=11)
ax4.set_ylabel('Inductance L (mH)', fontsize=11)
ax4.set_title('Inductance vs Current (Saturation)', fontsize=12)
ax4.legend(fontsize=10)
ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig(f'{output_dir}/fig5_saturation.png', dpi=150, bbox_inches='tight')
plt.close()
print("  [OK] Figure 5 saved: fig5_saturation.png")

# ============================================================================
# Part 7: Magnetic Circuit - Electric Circuit Analogy Visualization
# ============================================================================
print("\n[7] Magnetic circuit - electric circuit analogy visualization...")

# ============================================================================
# 图6：磁路-电路类比示意图
# ============================================================================
fig, axes = plt.subplots(1, 2, figsize=(14, 6))
fig.suptitle('Magnetic Circuit - Electric Circuit Analogy', fontsize=14, fontweight='bold')

# 左侧：电路图
ax1 = axes[0]
ax1.set_xlim(0, 10)
ax1.set_ylim(0, 10)
ax1.set_aspect('equal')
ax1.axis('off')
ax1.set_title('Electric Circuit', fontsize=13, fontweight='bold', pad=20)

# 绘制电池（电动势源）
battery = FancyBboxPatch((1, 4), 1, 2, boxstyle="round,pad=0.1", 
                          facecolor='#ffeb3b', edgecolor='black', linewidth=2)
ax1.add_patch(battery)
ax1.text(1.5, 5, 'V', ha='center', va='center', fontsize=14, fontweight='bold')
ax1.text(1.5, 3.3, 'EMF', ha='center', va='center', fontsize=10)

# 绘制电阻
resistor_x = [4, 4.5, 4.5, 5.5, 5.5, 6]
resistor_y = [8, 8, 7, 7, 8, 8]
ax1.plot(resistor_x, resistor_y, 'k-', linewidth=2)
ax1.plot([4.5, 4.5], [7, 6.5], 'k-', linewidth=2)
ax1.plot([5.5, 5.5], [7, 6.5], 'k-', linewidth=2)
ax1.text(5, 8.5, 'R', ha='center', va='center', fontsize=12, fontweight='bold')

# 绘制导线
ax1.plot([2, 4], [5, 5], 'k-', linewidth=2)  # 电池到电阻
ax1.plot([4, 4], [5, 6.5], 'k-', linewidth=2)
ax1.plot([6, 8], [8, 8], 'k-', linewidth=2)
ax1.plot([8, 8], [8, 2], 'k-', linewidth=2)
ax1.plot([8, 2], [2, 2], 'k-', linewidth=2)
ax1.plot([2, 2], [2, 4], 'k-', linewidth=2)

# 绘制电流箭头
arrow = FancyArrowPatch((3, 5.5), (3.5, 5.5), arrowstyle='->', 
                        mutation_scale=20, linewidth=2, color='red')
ax1.add_patch(arrow)
ax1.text(3.25, 5.9, 'I', ha='center', va='center', fontsize=12, color='red', fontweight='bold')

# 添加公式
ax1.text(5, 1, 'Ohm\'s Law: V = I × R', ha='center', va='center', 
         fontsize=12, style='italic', bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5))

# 右侧：磁路图
ax2 = axes[1]
ax2.set_xlim(0, 10)
ax2.set_ylim(0, 10)
ax2.set_aspect('equal')
ax2.axis('off')
ax2.set_title('Magnetic Circuit', fontsize=13, fontweight='bold', pad=20)

# 绘制线圈（磁动势源）
for i in range(5):
    circle = Circle((1.5, 3.5 + i*0.4), 0.3, fill=False, edgecolor='blue', linewidth=2)
    ax2.add_patch(circle)
ax2.text(1.5, 2.5, 'N turns', ha='center', va='center', fontsize=9)
ax2.text(1.5, 2, 'I', ha='center', va='center', fontsize=11, fontweight='bold', color='blue')
ax2.text(0.5, 5, 'F=NI', ha='center', va='center', fontsize=11, color='blue', fontweight='bold')

# 绘制磁芯
# 左侧
ax2.plot([2, 2], [3, 7], 'k-', linewidth=8, solid_capstyle='round')
# 顶部
ax2.plot([2, 8], [7, 7], 'k-', linewidth=8, solid_capstyle='round')
# 右侧
ax2.plot([8, 8], [7, 3], 'k-', linewidth=8, solid_capstyle='round')
# 底部（带气隙）
ax2.plot([2, 4.5], [3, 3], 'k-', linewidth=8, solid_capstyle='round')
ax2.plot([5.5, 8], [3, 3], 'k-', linewidth=8, solid_capstyle='round')

# 气隙标记
ax2.plot([4.5, 4.5], [2.5, 3.5], 'r--', linewidth=1)
ax2.plot([5.5, 5.5], [2.5, 3.5], 'r--', linewidth=1)
ax2.text(5, 2, 'g', ha='center', va='center', fontsize=11, color='red', fontweight='bold')

# 磁通箭头
arrow2 = FancyArrowPatch((3, 7.5), (4, 7.5), arrowstyle='->', 
                         mutation_scale=20, linewidth=2, color='green')
ax2.add_patch(arrow2)
ax2.text(3.5, 8, 'Φ', ha='center', va='center', fontsize=12, color='green', fontweight='bold')

# 添加公式
ax2.text(5, 0.8, 'Ohm\'s Law: F = Φ × Rm', ha='center', va='center', 
         fontsize=12, style='italic', bbox=dict(boxstyle='round', facecolor='lightblue', alpha=0.5))

plt.tight_layout()
plt.savefig(f'{output_dir}/fig6_analogy.png', dpi=150, bbox_inches='tight')
plt.close()
print("  [OK] Figure 6 saved: fig6_analogy.png")

# ============================================================================
# Figure 7: Magnetic Circuit - Electric Circuit Analogy Table
# ============================================================================
fig, ax = plt.subplots(figsize=(12, 8))
ax.axis('off')
ax.set_title('Magnetic Circuit vs Electric Circuit: Complete Analogy', 
             fontsize=14, fontweight='bold', pad=20)

# Create table data
table_data = [
    ['Quantity', 'Electric Circuit', 'Symbol', 'Unit', 'Magnetic Circuit', 'Symbol', 'Unit'],
    ['Source', 'Electromotive Force', 'V', 'Volt (V)', 'Magnetomotive Force', 'F', 'Ampere-turn (At)'],
    ['Flow', 'Current', 'I', 'Ampere (A)', 'Magnetic Flux', 'Phi', 'Weber (Wb)'],
    ['Opposition', 'Resistance', 'R', 'Ohm', 'Reluctance', 'Rm', 'At/Wb or H-1'],
    ['Conductance', 'Conductance', 'G', 'Siemens (S)', 'Permeance', 'Lambda', 'Wb/At'],
    ['Field Intensity', 'Electric Field', 'E', 'V/m', 'Magnetic Field', 'H', 'At/m'],
    ['Flux Density', 'Current Density', 'J', 'A/m2', 'Flux Density', 'B', 'Tesla (T)'],
    ['Material Property', 'Conductivity', 'sigma', 'S/m', 'Permeability', 'mu', 'H/m'],
    ['Law', 'Ohm Law', 'V=IR', '-', 'Ohm Law', 'F=Phi*Rm', '-'],
    ['KCL', 'Sum I=0', '-', '-', 'Sum Phi=0', '-', '-'],
    ['KVL', 'Sum V=0', '-', '-', 'Sum F=Sum NI', '-', '-'],
]

# 绘制表格
table = ax.table(cellText=table_data[1:], colLabels=table_data[0],
                 cellLoc='center', loc='center',
                 colWidths=[0.12, 0.15, 0.08, 0.12, 0.18, 0.08, 0.12])

table.auto_set_font_size(False)
table.set_fontsize(10)
table.scale(1, 2.5)

# 设置表头样式
for i in range(7):
    table[(0, i)].set_facecolor('#4CAF50')
    table[(0, i)].set_text_props(weight='bold', color='white')

# 设置交替行颜色
for i in range(1, len(table_data)):
    for j in range(7):
        if i % 2 == 0:
            table[(i, j)].set_facecolor('#f0f0f0')
        else:
            table[(i, j)].set_facecolor('#ffffff')

plt.tight_layout()
plt.savefig(f'{output_dir}/fig7_analogy_table.png', dpi=150, bbox_inches='tight')
plt.close()
print("  [OK] Figure 7 saved: fig7_analogy_table.png")

# ============================================================================
# Part 8: Complex Magnetic Circuit Analysis Example
# ============================================================================
print("\n[8] Complex magnetic circuit analysis example...")

# Define a complex circuit: E-core with center leg gap
# Geometric parameters
l_side = 0.08      # Side leg magnetic path length (m)
l_center = 0.04    # Center leg magnetic path length (m)
l_yoke = 0.06      # Yoke magnetic path length (m)
A_side = 1e-4      # Side leg cross-sectional area (m^2)
A_center = 2e-4    # Center leg cross-sectional area (m^2)
g_center = 0.5e-3  # Center leg air gap (m)

# Calculate reluctance for each section
R_side = l_side / (mu_0 * mu_r_iron * A_side)
R_center_core = l_center / (mu_0 * mu_r_iron * A_center)
R_center_gap = g_center / (mu_0 * A_center)
R_yoke = l_yoke / (mu_0 * mu_r_iron * A_side)

# Equivalent circuit: two side legs in parallel, then in series with center leg (with gap)
R_side_parallel = R_side / 2
R_center_total = R_center_core + R_center_gap
R_total_complex = R_yoke + R_side_parallel + R_center_total

print(f"  E-core magnetic circuit analysis:")
print(f"    Side leg reluctance: {R_side:.4e} At/Wb")
print(f"    Center leg core reluctance: {R_center_core:.4e} At/Wb")
print(f"    Center leg gap reluctance: {R_center_gap:.4e} At/Wb")
print(f"    Yoke reluctance: {R_yoke:.4e} At/Wb")
print(f"    Total reluctance: {R_total_complex:.4e} At/Wb")

# 计算不同电流下的磁通分布
I_complex = np.linspace(0.1, 5, 100)
F_complex = N * I_complex
Phi_total = F_complex / R_total_complex

# 磁通分配（并联支路）
Phi_side = Phi_total / 2  # 每个侧柱的磁通
Phi_center = Phi_total    # 中心柱的磁通

# 磁通密度
B_side = Phi_side / A_side
B_center = Phi_center / A_center

# ============================================================================
# 图8：复杂磁路分析
# ============================================================================
fig, axes = plt.subplots(2, 2, figsize=(12, 10))
fig.suptitle('Complex Magnetic Circuit Analysis (E-Core)', fontsize=14, fontweight='bold')

# 子图1：磁通-电流特性
ax1 = axes[0, 0]
ax1.plot(I_complex, Phi_total*1e4, 'b-', linewidth=2, label='Total Flux')
ax1.plot(I_complex, Phi_side*1e4, 'r--', linewidth=1.5, label='Side Leg Flux (each)')
ax1.plot(I_complex, Phi_center*1e4, 'g-.', linewidth=1.5, label='Center Leg Flux')
ax1.set_xlabel('Current I (A)', fontsize=11)
ax1.set_ylabel('Magnetic Flux Φ (10⁻⁴ Wb)', fontsize=11)
ax1.set_title('Flux Distribution', fontsize=12)
ax1.legend(fontsize=9)
ax1.grid(True, alpha=0.3)

# 子图2：磁通密度对比
ax2 = axes[0, 1]
ax2.plot(I_complex, B_side, 'r-', linewidth=2, label='Side Leg B')
ax2.plot(I_complex, B_center, 'g-', linewidth=2, label='Center Leg B')
ax2.set_xlabel('Current I (A)', fontsize=11)
ax2.set_ylabel('Flux Density B (T)', fontsize=11)
ax2.set_title('Flux Density Comparison', fontsize=12)
ax2.legend(fontsize=10)
ax2.grid(True, alpha=0.3)

# 子图3：磁阻分布（饼图）
ax3 = axes[1, 0]
labels = ['Yoke', 'Side Legs\n(Parallel)', 'Center Core', 'Center Gap']
sizes = [R_yoke/R_total_complex*100, R_side_parallel/R_total_complex*100,
         R_center_core/R_total_complex*100, R_center_gap/R_total_complex*100]
colors = ['#ff9999', '#66b3ff', '#99ff99', '#ffcc99']
explode = (0, 0, 0, 0.1)
ax3.pie(sizes, explode=explode, labels=labels, colors=colors, autopct='%1.1f%%',
        shadow=True, startangle=90)
ax3.set_title('Reluctance Distribution', fontsize=12)

# 子图4：电感特性
ax4 = axes[1, 1]
L_complex = N**2 / R_total_complex
L_values = [L_complex*1e3] * len(I_complex)
ax4.plot(I_complex, L_values, 'b-', linewidth=2, label='Inductance')
ax4.set_xlabel('Current I (A)', fontsize=11)
ax4.set_ylabel('Inductance L (mH)', fontsize=11)
ax4.set_title(f'Inductance = {L_complex*1e3:.2f} mH', fontsize=12)
ax4.grid(True, alpha=0.3)
ax4.legend(fontsize=10)

plt.tight_layout()
plt.savefig(f'{output_dir}/fig8_complex_circuit.png', dpi=150, bbox_inches='tight')
plt.close()
print("  [OK] Figure 8 saved: fig8_complex_circuit.png")

# ============================================================================
# Part 9: Parameter Sensitivity Analysis
# ============================================================================
print("\n[9] Parameter sensitivity analysis...")

# Analyze effect of air gap length on inductance
N_sens = 100
I_sens = 1.0
g_range_sens = np.linspace(0.1e-3, 10e-3, 100)

L_vs_g = []
for g_val in g_range_sens:
    R_g = g_val / (mu_0 * A_core)
    R_c = l_core / (mu_0 * mu_r_iron * A_core)
    L_vs_g.append(N_sens**2 / (R_c + R_g))
L_vs_g = np.array(L_vs_g)

# Analyze effect of turns on inductance
g_sens = 1e-3
N_range_sens = np.arange(10, 510, 10)
R_g_sens = g_sens / (mu_0 * A_core)
R_c_sens = l_core / (mu_0 * mu_r_iron * A_core)
L_vs_N = N_range_sens**2 / (R_c_sens + R_g_sens)

# Analyze effect of permeability on inductance
mu_r_range = np.linspace(100, 10000, 100)
L_vs_mu = []
for mu_r_val in mu_r_range:
    R_c_val = l_core / (mu_0 * mu_r_val * A_core)
    R_g_val = g_sens / (mu_0 * A_core)
    L_vs_mu.append(N_sens**2 / (R_c_val + R_g_val))
L_vs_mu = np.array(L_vs_mu)

# ============================================================================
# 图9：参数敏感性分析
# ============================================================================
fig, axes = plt.subplots(1, 3, figsize=(15, 5))
fig.suptitle('Parameter Sensitivity Analysis', fontsize=14, fontweight='bold')

# 子图1：电感 vs 气隙长度
ax1 = axes[0]
ax1.semilogy(g_range_sens*1000, L_vs_g*1e3, 'b-', linewidth=2)
ax1.set_xlabel('Gap Length g (mm)', fontsize=11)
ax1.set_ylabel('Inductance L (mH)', fontsize=11)
ax1.set_title(f'Sensitivity to Gap (N={N_sens})', fontsize=12)
ax1.grid(True, alpha=0.3)

# 子图2：电感 vs 匝数
ax2 = axes[1]
ax2.plot(N_range_sens, L_vs_N*1e3, 'r-', linewidth=2)
ax2.set_xlabel('Number of Turns N', fontsize=11)
ax2.set_ylabel('Inductance L (mH)', fontsize=11)
ax2.set_title(f'Sensitivity to Turns (g={g_sens*1000:.1f}mm)', fontsize=12)
ax2.grid(True, alpha=0.3)

# 子图3：电感 vs 磁导率
ax3 = axes[2]
ax3.semilogx(mu_r_range, L_vs_mu*1e3, 'g-', linewidth=2)
ax3.axvline(x=mu_r_iron, color='r', linestyle='--', alpha=0.7, label=f'μr={mu_r_iron}')
ax3.set_xlabel('Relative Permeability μr', fontsize=11)
ax3.set_ylabel('Inductance L (mH)', fontsize=11)
ax3.set_title(f'Sensitivity to Permeability', fontsize=12)
ax3.legend(fontsize=10)
ax3.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig(f'{output_dir}/fig9_sensitivity.png', dpi=150, bbox_inches='tight')
plt.close()
print("  [OK] Figure 9 saved: fig9_sensitivity.png")

# ============================================================================
# Part 10: Summary of Results
# ============================================================================
print("\n[10] Summary of simulation results...")

# Calculate key parameters
print("\n" + "=" * 60)
print("Key Parameters Summary")
print("=" * 60)

print(f"\n[Basic Parameters]")
print(f"  Vacuum permeability mu_0 = {mu_0:.4e} H/m")
print(f"  Iron core relative permeability mu_r = {mu_r_iron}")
print(f"  Core length l = {l_core*1000:.1f} mm")
print(f"  Core cross-sectional area A = {A_core*1e6:.1f} mm^2")
print(f"  Air gap length g = {g*1000:.1f} mm")
print(f"  Number of turns N = {N}")

print(f"\n[Magnetic Circuit Parameters]")
print(f"  Core reluctance R_core = {R_core:.4e} At/Wb")
print(f"  Air gap reluctance R_gap = {R_gap:.4e} At/Wb")
print(f"  Total reluctance R_total = {R_total:.4e} At/Wb")
print(f"  Air gap reluctance ratio = {R_gap/R_total*100:.1f}%")

print(f"\n[Operating Point (I=1A)]")
F_1A = N * 1
Phi_1A = F_1A / R_total
B_1A = Phi_1A / A_core
L_calc = N**2 / R_total
print(f"  MMF F = {F_1A:.1f} At")
print(f"  Flux Phi = {Phi_1A:.4e} Wb")
print(f"  Flux density B = {B_1A:.4f} T")
print(f"  Inductance L = {L_calc*1e3:.2f} mH")

print(f"\n[Complex Circuit (E-Core)]")
print(f"  Total reluctance = {R_total_complex:.4e} At/Wb")
print(f"  Inductance = {N**2/R_total_complex*1e3:.2f} mH")

print("\n" + "=" * 60)
print("Simulation completed! Generated 9 figures")
print("=" * 60)

# Generate file list
print("\nGenerated files:")
files = [
    'fig1_ohms_law.png - Magnetic Circuit Ohm Law Verification',
    'fig2_fringing_effect.png - Air Gap Fringing Effect Analysis',
    'fig3_series_parallel.png - Magnetic Circuit Series/Parallel Analysis',
    'fig4_inductance.png - Inductance Calculation and Characteristics',
    'fig5_saturation.png - Magnetic Saturation Effect Analysis',
    'fig6_analogy.png - Magnetic-Electric Circuit Analogy Diagram',
    'fig7_analogy_table.png - Magnetic-Electric Circuit Analogy Table',
    'fig8_complex_circuit.png - Complex Magnetic Circuit Analysis',
    'fig9_sensitivity.png - Parameter Sensitivity Analysis'
]
for f in files:
    print(f"  [OK] {f}")

print(f"\nAll files saved to: {output_dir}")
print("\nSimulation completed!")

AtomGit开源社区

AtomGit 是由开放原子开源基金会联合 CSDN 等生态伙伴共同推出的新一代开源与人工智能协作平台。平台坚持“开放、中立、公益”的理念，把代码托管、模型共享、数据集托管、智能体开发体验和算力服务整合在一起，为开发者提供从开发、训练到部署的一站式体验。

更多推荐

金融科技：利用AI风控模型重构信贷审批流程

AtomGit开源社区

【继电保护】小电流接地系统故障仿真-中性点不接地与经消弧线圈接地仿真模型（Simulink仿真实现）

小电流接地系统故障仿真是电力系统中非常重要的研究领域，特别是针对中性点不接地和经消弧线圈接地的情况。这两种故障情况在电力系统中都可能发生，因此对其进行仿真模型研究具有重要意义。中性点不接地故障是指变压器或发电机中性点没有接地，这种情况下，如果出现了单相接地故障，会导致系统中产生零序电流，可能对设备和系统造成严重损坏。因此，针对中性点不接地故障，需要建立相应的仿真模型，研究其对电力系统的影响，以及采

AtomGit开源社区

金融行业 Multi-Agent 落地实践：智能投顾与风险监控的协同架构设计

凌晨2点，纽约曼哈顿下城，高盛大厦的交易大厅依旧灯火通明。交易员Alex盯着屏幕上闪烁的数字，额头上渗出细密的汗珠——市场正在经历一场罕见的"闪电崩盘"，而他负责的投资组合正面临巨大风险。"AI系统怎么还没反应？"Alex对着电话怒吼。电话那头，数据科学家Sarah正在 frantically 调试系统：“我们的风险监控模型发现了异常，但智能投顾系统还在按照原定策略执行交易……两个系统’吵架’了！