主题026：变压器电磁场仿真——三相电力变压器电磁场分析

目录

1. 引言
2. 变压器基础理论
3. 变压器结构与工作原理
4. 变压器电磁场数学模型
5. 案例一：变压器几何建模
6. 案例二：空载磁场计算
7. 案例三：漏磁场与短路阻抗分析
8. 案例四：损耗计算与效率分析
9. 案例五：综合性能评估
10. 代码深度解析
11. 结果分析与讨论
12. 常见报错与解决方案
13. 进阶挑战与思考
14. 总结与展望

---

引言

电力变压器是电力系统中最重要的电气设备之一，承担着电压变换、电能传输和电气隔离的关键作用。从发电厂到终端用户，电能需要经过多次变压才能高效、安全地输送和使用。据统计，电力系统中变压器的总容量通常是发电机总容量的6-8倍，其性能直接影响着电网的安全稳定运行和电能质量。

变压器电磁场分析是变压器设计、优化和故障诊断的核心技术。通过电磁场仿真，工程师可以：

1. 优化铁芯设计：分析主磁通分布，避免局部饱和，降低铁损
2. 计算短路阻抗：评估变压器的抗短路能力和电压调整率
3. 评估漏磁场：预测绕组受力、局部过热和电磁兼容问题
4. 计算损耗分布：铁损、铜损的精确计算是效率评估和温升预测的基础
5. 故障分析：通过场分布异常识别绝缘缺陷、绕组变形等故障

传统的变压器设计主要依赖于经验公式和等效电路模型，这些方法虽然计算简便，但难以准确反映变压器内部复杂的电磁现象。随着有限元法等数值计算技术的发展，基于场的分析方法已成为现代变压器设计的主流手段。

本教程将带领读者深入学习三相电力变压器的电磁场仿真技术。我们将从基础理论出发，逐步构建完整的变压器电磁场分析框架，包括几何建模、磁场求解、损耗计算和性能评估等环节。通过本教程的学习，读者将能够：

1. 理解三相变压器的基本结构和工作原理
2. 掌握变压器电磁场的数学建模方法
3. 学会使用Python实现变压器磁场数值计算
4. 能够分析空载磁场、漏磁场和计算短路阻抗
5. 具备进行变压器损耗分析和效率评估的能力

本教程采用理论与实践相结合的方式，所有代码均经过实际运行验证，确保读者可以直接复制使用并获得预期结果。

---

变压器基础理论

2.1 电磁感应定律

变压器的工作原理基于法拉第电磁感应定律。当交变电流通过一次侧绕组时，产生交变磁通，该磁通通过铁芯耦合到二次侧绕组，在二次侧感应出电动势。

法拉第电磁感应定律：

$$e=-N\frac{d\Phi }{dt}$$

其中，$e$ 是感应电动势，$N$ 是绕组匝数，$\Phi$ 是磁通量。

对于正弦交流电，感应电动势的有效值为：

$$E=4.44fN{\Phi }_m=4.44fN{B}_{m}A$$

其中，$f$ 是频率，${\Phi }_m$ 是磁通幅值，$B_m$ 是磁密幅值，$A$ 是铁芯截面积。

2.2 变压器基本方程

电压方程：

一次侧：
$${\dot{U}}_1=-{\dot{E}}_1+{\dot{I}}_1{Z}_1$$

二次侧：
$${\dot{U}}_2={\dot{E}}_2-{\dot{I}}_2{Z}_2$$

其中，$Z_1$ 和 $Z_2$ 是绕组漏阻抗。

磁动势平衡：

$$N_1{\dot{I}}_1+N_2{\dot{I}}_2=N_1{\dot{I}}_0$$

其中，${\dot{I}}_0$ 是励磁电流。

等效电路参数：

• 励磁电阻 $R_m$：反映铁损的等效电阻
• 励磁电抗 $X_m$：反映主磁通的等效电抗
• 漏电阻 $R_1,R_2$：绕组电阻
• 漏电抗 $X_1,X_2$：反映漏磁通的等效电抗

2.3 磁路分析

变压器的磁路包括铁芯柱和铁轭。根据安培环路定律：

$$\oint \mathbf{H}\cdot d\mathbf{l}=NI$$

对于变压器磁路：

$$H_{leg}{l}_{leg}+H_{yoke}{l}_{yoke}=N_1{I}_0$$

其中，$H_{leg}$ 和 $H_{yoke}$ 分别是铁芯柱和铁轭的磁场强度，$l_{leg}$ 和 $l_{yoke}$ 是对应的磁路长度。

磁化特性：

铁磁材料的磁化特性通常用B-H曲线描述。在工程计算中，常用以下经验公式：

$$B=\frac{{\mu }_0{\mu }_{r}H}{1+\alpha |H|}$$

其中，$\alpha$ 是饱和系数。

2.4 漏磁场分析

漏磁通是指不链绕两个绕组的磁通，主要存在于：

1. 绕组间漏磁：高低压绕组之间的间隙
2. 绕组端部漏磁：绕组端部与铁轭之间的区域
3. 相间漏磁：三相绕组之间的漏磁

漏磁场的计算对于短路阻抗、绕组受力和杂散损耗的评估至关重要。

---

变压器结构与工作原理

3.1 三相变压器结构

三相电力变压器主要由以下部分组成：

铁芯部分：

• 铁芯柱：三个（或五个）垂直的铁芯柱，用于套装绕组
• 铁轭：上下水平部分，连接各铁芯柱，构成闭合磁路
• 夹件：固定铁芯的结构件

绕组部分：

• 一次侧绕组（高压侧）：匝数多，线径细，绝缘要求高
• 二次侧绕组（低压侧）：匝数少，线径粗，绝缘要求相对较低
• 绕组排列：同心式（低压在内、高压在外）或交叠式

其他部件：

• 油箱：容纳变压器油，提供绝缘和冷却
• 套管：引出绕组端子，提供绝缘
• 分接开关：调节变比，适应电网电压变化
• 冷却系统：散热器、风扇、油泵等

3.2 磁路结构类型

根据铁芯结构，三相变压器可分为：

三相三柱式（芯式）：

• 三个铁芯柱在同一平面内
• 中间柱（B相）磁路短，磁阻小
• 两侧柱（A、C相）磁路长，磁阻大
• 适用于容量不太大的变压器

三相五柱式：

• 增加两个旁轭，为三相磁通提供独立回路
• 三相磁路对称，零序磁通有通路
• 适用于大容量变压器和需要承受不平衡负载的场合

3.3 工作原理

空载运行：

当一次侧施加额定电压，二次侧开路时：

• 一次侧流过励磁电流（约为额定电流的2-5%）
• 产生主磁通，在铁芯中闭合
• 主磁通在一次、二次侧感应电动势
• 感应电动势与匝数成正比：$\frac{E_1}{E_2}=\frac{N_1}{N_2}$

负载运行：

当二次侧接负载时：

• 二次侧产生电流，建立去磁磁动势
• 一次侧电流自动增加，维持磁动势平衡
• 能量从一次侧传递到二次侧
• 变比：$\frac{U_1}{U_2}\approx \frac{N_1}{N_2}=k$

3.4 关键性能指标

1. 额定容量 $S_N$：变压器能持续输出的视在功率（kVA）
2. 额定电压比 $k$：一次与二次额定电压之比
3. 短路阻抗 $u_k$：反映变压器内阻抗压降的参数（%）
4. 空载损耗 $P_0$：铁损，与负载无关
5. 负载损耗 $P_k$：铜损，与负载电流平方成正比
6. 效率 $\eta$：输出功率与输入功率之比
7. 电压调整率：负载变化引起的二次电压变化

---

变压器电磁场数学模型

4.1 问题描述与假设

变压器电磁场分析通常采用以下假设：

1. 二维场假设：忽略轴向变化，假设场量在横截面内分布
2. 正弦稳态：电压、电流为正弦波形，采用相量法分析
3. 准静态近似：忽略位移电流，满足 $\nabla \times \mathbf{H}=\mathbf{J}$
4. 材料线性化：铁芯磁导率分段线性或采用等效值

求解区域：

典型的二维求解区域包括：

• 铁芯区域（铁芯柱和铁轭）
• 绕组区域（一次侧和二次侧绕组）
• 绝缘和油隙区域
• 油箱壁区域（边界条件）

4.2 控制方程

在二维直角坐标系下，矢量磁位 $A_z$ 满足泊松方程：

$$\frac{\partial }{\partial x}\left(\frac{1}{\mu }\frac{\partial A_z}{\partial x}\right)+\frac{\partial }{\partial y}\left(\frac{1}{\mu }\frac{\partial A_z}{\partial y}\right)=-J_z$$

磁感应强度分量：

$$B_x=\frac{\partial A_z}{\partial y},B_y=-\frac{\partial A_z}{\partial x}$$

4.3 边界条件

齐次狄利克雷边界：

在远场边界或对称边界上：

$$A_z=0$$

周期性边界条件：

对于三相变压器，可利用对称性简化计算。

不同材料交界面：

• 磁感应强度法向分量连续：$B_{n1}=B_{n2}$
• 磁场强度切向分量连续：$H_{t1}=H_{t2}$

4.4 短路阻抗计算

短路阻抗是变压器的重要参数，可通过漏磁场能量计算：

$$Z_k=\sqrt{R_k^2+X_k^2}$$

其中，漏电抗：

$$X_k=\omega L_k=\omega \frac{2W_{leak}}{I_k^2}$$

漏磁场能量：

$$W_{leak}=\frac{1}{2}{\int }_V\frac{B^2}{{\mu }_0}dV$$

---

案例一：变压器几何建模

5.1 案例目标

本案例将学习如何建立三相电力变压器的几何模型。准确的几何模型是电磁场仿真的基础，直接影响计算结果的可靠性。

5.2 变压器参数定义

首先定义变压器的基本参数：

@dataclass
class TransformerParameters:
    """三相电力变压器参数"""
    # 额定参数
    rated_power: float = 1000.0  # 额定容量 [kVA]
    primary_voltage: float = 10.0  # 一次侧额定电压 [kV]
    secondary_voltage: float = 0.4  # 二次侧额定电压 [kV]
    frequency: float = 50.0  # 频率 [Hz]
    
    # 铁芯参数 (单位: mm)
    core_width: float = 200.0  # 铁芯柱宽度
    core_depth: float = 150.0  # 铁芯叠厚
    core_height: float = 400.0  # 铁芯窗高
    window_width: float = 120.0  # 铁芯窗宽
    yoke_height: float = 80.0  # 铁轭高度
    
    # 绕组参数
    primary_turns: int = 1000  # 一次侧匝数
    secondary_turns: int = 40  # 二次侧匝数
    winding_thickness: float = 30.0  # 绕组厚度
    winding_height: float = 350.0  # 绕组高度
    
    # 材料属性
    core_mu_r: float = 5000.0  # 铁芯相对磁导率
    core_conductivity: float = 2.0e6  # 铁芯电导率 [S/m]
    core_density: float = 7650.0  # 铁芯密度 [kg/m³]
    
    # 损耗系数
    core_loss_coeff: float = 2.5  # 铁损系数 [W/kg]
    copper_resistivity: float = 1.72e-8  # 铜电阻率 [Ω·m]

参数说明：

1. 额定容量1000kVA：典型的配电变压器容量
2. 电压比10/0.4kV：标准的中压配电变压器
3. 匝数比1000:40=25:1：与电压比匹配
4. 铁芯尺寸：根据容量和经验公式确定
5. 材料参数：硅钢片相对磁导率约5000，铜电阻率1.72×10⁻⁸Ω·m

5.3 几何建模类

创建变压器几何建模类：

class TransformerGeometry:
    """变压器几何建模"""
    
    def __init__(self, params: TransformerParameters):
        self.params = params
        
        # 转换为米
        self.core_w = params.core_width / 1000
        self.core_d = params.core_depth / 1000
        self.core_h = params.core_height / 1000
        self.window_w = params.window_width / 1000
        self.yoke_h = params.yoke_height / 1000
        
        # 计算总尺寸
        self.total_width = 3 * self.core_w + 2 * self.window_w
        self.total_height = self.core_h + 2 * self.yoke_h

几何计算方法：

def get_core_regions(self) -> List[Dict]:
    """获取铁芯区域定义"""
    regions = []
    
    # 三个铁芯柱 (A相、B相、C相)
    for i in range(3):
        x_center = self.core_w/2 + i * (self.core_w + self.window_w)
        regions.append({
            'type': 'core_leg',
            'phase': ['A', 'B', 'C'][i],
            'x': x_center - self.core_w/2,
            'y': self.yoke_h,
            'width': self.core_w,
            'height': self.core_h
        })
    
    # 上下铁轭
    # ...
    
    return regions

5.4 区域判断方法

def is_inside_core(self, x: float, y: float) -> bool:
    """判断点是否在铁芯内"""
    # 检查是否在铁芯柱内
    for i in range(3):
        x_center = self.core_w/2 + i * (self.core_w + self.window_w)
        if (x_center - self.core_w/2 <= x <= x_center + self.core_w/2 and
            self.yoke_h <= y <= self.yoke_h + self.core_h):
            return True
    
    # 检查是否在铁轭内
    # ...
    
    return False

---

案例二：空载磁场计算

6.1 案例目标

本案例将学习如何计算变压器的空载磁场分布。空载磁场主要由励磁电流产生，反映了变压器的主磁通分布特性。

6.2 励磁电流计算

空载时，一次侧流过励磁电流：

def solve_no_load_field(self):
    """求解空载磁场"""
    # 励磁电流 (约为额定电流的2-5%)
    I_excitation = self.params.primary_current * 0.03
    
    # 计算主磁通密度
    mu0 = 4 * np.pi * 1e-7
    mu_core = mu0 * self.params.core_mu_r
    core_length = 2 * self.geom.core_h + 2 * self.geom.yoke_h
    B_main = mu_core * self.params.primary_turns * I_excitation / core_length

物理意义：

• 励磁电流产生磁动势 $F=N_1{I}_0$
• 磁动势在磁路中产生磁通
• 磁通密度 $B=\mu H=\mu F/l$

6.3 磁场分布计算

# 填充磁场分布
for i in range(self.ny):
    for j in range(self.nx):
        x, y = self.X[i, j], self.Y[i, j]
        
        if self.geom.is_inside_core(x, y):
            # 铁芯内磁场
            self.Bx[i, j] = B_main * 0.8
            self.By[i, j] = B_main * 0.6
        else:
            # 铁芯外磁场 (漏磁)
            dist = self._distance_to_core(x, y)
            decay = np.exp(-dist / 0.1)
            self.Bx[i, j] = B_main * 0.1 * decay
            self.By[i, j] = B_main * 0.1 * decay

关键要点：

1. 铁芯内磁场：幅值高，方向沿磁路
2. 漏磁场：在铁芯外衰减，反映不完全耦合的磁通
3. 衰减模型：采用指数衰减模拟漏磁分布

---

案例三：漏磁场与短路阻抗分析

7.1 案例目标

本案例将学习如何计算变压器的漏磁场分布，并基于漏磁场能量计算短路阻抗。短路阻抗是变压器设计和运行中的关键参数。

7.2 漏磁场计算

短路工况下，绕组流过很大的短路电流，产生强烈的漏磁场：

def solve_short_circuit_field(self):
    """求解短路工况磁场"""
    # 短路电流
    Z_sc = self.params.short_circuit_impedance / 100
    I_sc = self.params.primary_current / Z_sc
    
    mu0 = 4 * np.pi * 1e-7
    
    # 漏磁场计算
    for i in range(self.ny):
        for j in range(self.nx):
            x, y = self.X[i, j], self.Y[i, j]
            
            winding_info = self._get_nearest_winding(x, y)
            
            if winding_info['in_winding']:
                # 绕组内部磁场
                self.Bx[i, j] = mu0 * self.params.primary_turns * I_sc / self.geom.winding_h * 0.5
            else:
                # 绕组外部漏磁场
                dist = winding_info['distance']
                B_leak = mu0 * self.params.primary_turns * I_sc / (2 * np.pi * max(dist, 0.01))
                self.Bx[i, j] = B_leak * 0.4

漏磁场特征：

1. 绕组内部：磁场强度与电流成正比
2. 绕组外部：随距离衰减，衰减规律与载流导体类似
3. 主要分布：集中在高低压绕组之间的间隙

7.3 短路阻抗计算

def calculate_short_circuit_impedance(self) -> Dict:
    """计算短路阻抗"""
    mu0 = 4 * np.pi * 1e-7
    
    # 漏感计算
    h = self.geom.winding_h
    delta = self.geom.winding_t * 2 + 0.01
    a = self.geom.core_w
    
    L_leak = mu0 * self.params.primary_turns**2 * h * delta / a
    
    # 短路阻抗
    omega = 2 * np.pi * self.params.frequency
    X_sc = omega * L_leak
    R_sc = 0.1  # 简化值
    Z_sc = np.sqrt(R_sc**2 + X_sc**2)
    
    # 标幺值
    Z_base = (self.params.primary_voltage * 1000)**2 / (self.params.rated_power * 1000)
    Z_sc_pu = Z_sc / Z_base * 100
    
    return {
        'R_sc': R_sc,
        'X_sc': X_sc,
        'Z_sc': Z_sc,
        'Z_sc_percent': Z_sc_pu,
        'L_leak': L_leak
    }

工程意义：

• 短路阻抗决定短路电流大小
• 影响电压调整率和并联运行
• 标准值通常为4%-10%

---

案例四：损耗计算与效率分析

8.1 案例目标

本案例将学习如何计算变压器的铁损和铜损，并评估变压器效率。损耗分析是变压器热设计和经济性评估的基础。

8.2 铁损计算

铁损（空载损耗）包括磁滞损耗和涡流损耗：

def calculate_core_loss(self) -> float:
    """计算铁芯损耗"""
    # 铁芯体积
    core_volume = (3 * self.geom.core_w * self.geom.core_h * self.geom.core_d +
                  2 * self.geom.total_width * self.geom.yoke_h * self.geom.core_d)
    
    # 铁芯质量
    core_mass = core_volume * self.params.core_density
    
    # 铁损 (Steinmetz方程简化)
    B_avg = np.mean(self.B_mag[self.B_mag > 0.1])
    P_core = self.params.core_loss_coeff * core_mass * (B_avg / 1.0)**2 * \
             (self.params.frequency / 50)**1.3
    
    return P_core

铁损组成：

1. 磁滞损耗：$P_h=k_{h}f{B}_m^n$，与频率和磁密幅值有关
2. 涡流损耗：$P_e=k_e{f}^2{B}_m^2$，与频率平方成正比
3. 附加损耗：由于局部饱和和异常涡流产生

8.3 铜损计算

铜损（负载损耗）是电流流过绕组电阻产生的损耗：

def calculate_copper_loss(self) -> float:
    """计算铜损"""
    # 绕组平均匝长
    mean_turn_length = 2 * (self.geom.core_w + self.geom.core_d) * np.pi
    
    # 一次侧绕组电阻
    R_primary = self.params.copper_resistivity * \
                self.params.primary_turns * mean_turn_length / \
                (np.pi * (0.005)**2)
    
    # 二次侧绕组电阻
    R_secondary = self.params.copper_resistivity * \
                  self.params.secondary_turns * mean_turn_length / \
                  (np.pi * (0.025)**2)
    
    # 铜损
    P_cu_primary = R_primary * self.params.primary_current**2
    P_cu_secondary = R_secondary * self.params.secondary_current**2
    
    return P_cu_primary + P_cu_secondary

8.4 效率计算

total_loss = core_loss + copper_loss
efficiency = params.rated_power * 1000 / (params.rated_power * 1000 + total_loss) * 100

效率特性：

• 配电变压器效率通常在95%-99%
• 最高效率点通常在50%-75%负载
• 大型变压器效率可达99%以上

---

案例五：综合性能评估

9.1 案例目标

本案例将综合前面所学内容，对变压器进行全面的性能评估，包括电气性能、损耗特性和运行特性分析。

9.2 性能参数汇总

基于前面的计算，可以得到以下性能参数：

results = {
    'transformer_params': {
        'rated_power_kva': params.rated_power,
        'primary_voltage_kv': params.primary_voltage,
        'secondary_voltage_kv': params.secondary_voltage,
        'turns_ratio': params.turns_ratio
    },
    'performance': {
        'core_loss_w': core_loss,
        'copper_loss_w': copper_loss,
        'total_loss_w': total_loss,
        'efficiency_percent': efficiency,
        'short_circuit_impedance_percent': sc_results['Z_sc_percent']
    }
}

9.3 运行特性分析

电压调整率：

$$\Delta U=\frac{U_{20}-U_2}{U_{20}}\times 100\%$$

与短路阻抗和负载功率因数有关。

温升估算：

基于损耗和散热条件，可以估算绕组和铁芯的温升。

---

代码深度解析

10.1 整体架构

代码采用模块化设计：

1. TransformerParameters：数据类，存储所有参数
2. TransformerGeometry：几何建模类
3. TransformerFieldSolver：磁场求解器类
4. 可视化函数：创建各种图表

10.2 关键算法

磁场叠加：

空载磁场和漏磁场分别计算，反映了不同的物理机制。

损耗模型：

采用简化工程模型，便于快速估算。

---

结果分析与讨论

11.1 磁场分布特征

从仿真结果可以观察到：

1. 空载磁场：主要集中在铁芯内，漏磁较小
2. 漏磁场：分布在绕组周围，影响短路阻抗
3. 饱和效应：本模型未考虑，实际需关注

11.2 性能评估

仿真得到的性能参数可用于：

1. 验证设计是否满足要求
2. 与其他方案比较
3. 指导优化设计

---

常见报错与解决方案

12.1 数值不稳定

**现象：**磁场计算结果异常大或异常小

**解决：**检查单位换算，确保使用国际单位制

12.2 内存不足

**现象：**程序崩溃或运行缓慢

**解决：**减小网格数量，或采用稀疏矩阵存储

---

进阶挑战与思考

13.1 思考题

1. 改变铁芯材料，观察对铁损的影响
2. 优化绕组排列，降低漏磁场
3. 考虑非线性磁化特性

13.2 扩展应用

1. 暂态过程分析（励磁涌流）
2. 多物理场耦合（电磁-热）
3. 优化设计算法

---

总结与展望

本教程系统介绍了变压器电磁场仿真的基本原理和实现方法。通过五个案例，读者可以掌握：

1. 变压器几何建模方法
2. 空载磁场和漏磁场计算
3. 短路阻抗和损耗计算
4. 性能评估方法

变压器仿真技术正在向多物理场耦合、智能化设计方向发展，为电力系统的高效、可靠运行提供支撑。

---

附录：完整代码

"""
主题026: 变压器电磁场仿真
========================
三相电力变压器电磁场分析

本代码介绍：
1. 变压器几何建模
2. 铁芯磁场计算
3. 漏磁场分析
4. 损耗计算
5. 短路阻抗计算
"""

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle, FancyBboxPatch, Circle
from matplotlib.collections import PatchCollection
import json
import time
from typing import Dict, Tuple, List, Optional
from dataclasses import dataclass
from enum import Enum
import warnings
warnings.filterwarnings('ignore')

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False


# ============================================================
# 第一部分：变压器参数定义
# ============================================================

@dataclass
class TransformerParameters:
    """三相电力变压器参数"""
    # 额定参数
    rated_power: float = 1000.0  # 额定容量 [kVA]
    primary_voltage: float = 10.0  # 一次侧额定电压 [kV]
    secondary_voltage: float = 0.4  # 二次侧额定电压 [kV]
    frequency: float = 50.0  # 频率 [Hz]
    
    # 铁芯参数 (单位: mm)
    core_width: float = 200.0  # 铁芯柱宽度
    core_depth: float = 150.0  # 铁芯叠厚
    core_height: float = 400.0  # 铁芯窗高
    window_width: float = 120.0  # 铁芯窗宽
    yoke_height: float = 80.0  # 铁轭高度
    
    # 绕组参数
    primary_turns: int = 1000  # 一次侧匝数
    secondary_turns: int = 40  # 二次侧匝数
    winding_thickness: float = 30.0  # 绕组厚度
    winding_height: float = 350.0  # 绕组高度
    
    # 材料属性
    core_mu_r: float = 5000.0  # 铁芯相对磁导率
    core_conductivity: float = 2.0e6  # 铁芯电导率 [S/m]
    core_density: float = 7650.0  # 铁芯密度 [kg/m³]
    
    # 损耗系数
    core_loss_coeff: float = 2.5  # 铁损系数 [W/kg]
    copper_resistivity: float = 1.72e-8  # 铜电阻率 [Ω·m]
    
    # 短路阻抗 [%]
    short_circuit_impedance: float = 4.0
    
    @property
    def turns_ratio(self) -> float:
        """匝数比"""
        return self.primary_turns / self.secondary_turns
    
    @property
    def primary_current(self) -> float:
        """一次侧额定电流 [A]"""
        return self.rated_power * 1000 / (np.sqrt(3) * self.primary_voltage * 1000)
    
    @property
    def secondary_current(self) -> float:
        """二次侧额定电流 [A]"""
        return self.rated_power * 1000 / (np.sqrt(3) * self.secondary_voltage * 1000)


# ============================================================
# 第二部分：变压器几何建模
# ============================================================

class TransformerGeometry:
    """变压器几何建模"""
    
    def __init__(self, params: TransformerParameters):
        self.params = params
        
        # 转换为米
        self.core_w = params.core_width / 1000
        self.core_d = params.core_depth / 1000
        self.core_h = params.core_height / 1000
        self.window_w = params.window_width / 1000
        self.yoke_h = params.yoke_height / 1000
        self.winding_t = params.winding_thickness / 1000
        self.winding_h = params.winding_height / 1000
        
        # 计算总尺寸
        self.total_width = 3 * self.core_w + 2 * self.window_w
        self.total_height = self.core_h + 2 * self.yoke_h
        
        print(f"变压器几何参数:")
        print(f"  额定容量: {params.rated_power} kVA")
        print(f"  电压比: {params.primary_voltage}/{params.secondary_voltage} kV")
        print(f"  铁芯柱宽: {params.core_width} mm")
        print(f"  窗高: {params.core_height} mm")
        print(f"  窗宽: {params.window_width} mm")
        print(f"  总宽度: {self.total_width*1000:.1f} mm")
        print(f"  总高度: {self.total_height*1000:.1f} mm")
    
    def get_core_regions(self) -> List[Dict]:
        """获取铁芯区域定义"""
        regions = []
        
        # 三个铁芯柱 (A相、B相、C相)
        for i in range(3):
            x_center = self.core_w/2 + i * (self.core_w + self.window_w)
            regions.append({
                'type': 'core_leg',
                'phase': ['A', 'B', 'C'][i],
                'x': x_center - self.core_w/2,
                'y': self.yoke_h,
                'width': self.core_w,
                'height': self.core_h
            })
        
        # 上铁轭
        regions.append({
            'type': 'core_yoke',
            'position': 'top',
            'x': 0,
            'y': self.yoke_h + self.core_h,
            'width': self.total_width,
            'height': self.yoke_h
        })
        
        # 下铁轭
        regions.append({
            'type': 'core_yoke',
            'position': 'bottom',
            'x': 0,
            'y': 0,
            'width': self.total_width,
            'height': self.yoke_h
        })
        
        return regions
    
    def get_winding_regions(self) -> List[Dict]:
        """获取绕组区域定义"""
        regions = []
        
        for i in range(3):
            x_center = self.core_w/2 + i * (self.core_w + self.window_w)
            
            # 一次侧绕组 (高压侧，靠近铁芯)
            regions.append({
                'type': 'primary_winding',
                'phase': ['A', 'B', 'C'][i],
                'x': x_center + self.core_w/2 + 0.005,  # 5mm间隙
                'y': self.yoke_h + (self.core_h - self.winding_h)/2,
                'width': self.winding_t,
                'height': self.winding_h
            })
            
            # 二次侧绕组 (低压侧，外侧)
            regions.append({
                'type': 'secondary_winding',
                'phase': ['A', 'B', 'C'][i],
                'x': x_center + self.core_w/2 + 0.005 + self.winding_t + 0.005,
                'y': self.yoke_h + (self.core_h - self.winding_h)/2,
                'width': self.winding_t,
                'height': self.winding_h
            })
        
        return regions
    
    def is_inside_core(self, x: float, y: float) -> bool:
        """判断点是否在铁芯内"""
        # 检查是否在铁芯柱内
        for i in range(3):
            x_center = self.core_w/2 + i * (self.core_w + self.window_w)
            if (x_center - self.core_w/2 <= x <= x_center + self.core_w/2 and
                self.yoke_h <= y <= self.yoke_h + self.core_h):
                return True
        
        # 检查是否在上铁轭内
        if (0 <= x <= self.total_width and
            self.yoke_h + self.core_h <= y <= self.total_height):
            return True
        
        # 检查是否在下铁轭内
        if (0 <= x <= self.total_width and
            0 <= y <= self.yoke_h):
            return True
        
        return False
    
    def is_inside_winding(self, x: float, y: float) -> Tuple[bool, str]:
        """判断点是否在绕组内，返回(是否在绕组内, 绕组类型)"""
        for region in self.get_winding_regions():
            if (region['x'] <= x <= region['x'] + region['width'] and
                region['y'] <= y <= region['y'] + region['height']):
                return True, region['type']
        return False, ''


# ============================================================
# 第三部分：磁场求解器
# ============================================================

class TransformerFieldSolver:
    """变压器磁场求解器"""
    
    def __init__(self, geometry: TransformerGeometry, params: TransformerParameters):
        self.geom = geometry
        self.params = params
        
        # 创建计算网格
        self.nx = 80
        self.ny = 60
        
        self.x = np.linspace(-0.05, self.geom.total_width + 0.05, self.nx)
        self.y = np.linspace(-0.05, self.geom.total_height + 0.05, self.ny)
        
        self.X, self.Y = np.meshgrid(self.x, self.y)
        
        # 磁场分量
        self.Bx = np.zeros_like(self.X)
        self.By = np.zeros_like(self.X)
        self.B_mag = np.zeros_like(self.X)
        
        # 磁位
        self.Az = np.zeros_like(self.X)
        
        print(f"\n磁场求解器初始化:")
        print(f"  水平网格: {self.nx}")
        print(f"  垂直网格: {self.ny}")
        print(f"  总网格数: {self.nx * self.ny}")
    
    def calculate_magnetomotive_force(self, phase: str, time: float = 0.0) -> float:
        """
        计算指定相的磁动势
        
        Parameters:
        -----------
        phase : 'A', 'B', 或 'C'
        time : 时间 [s]
        """
        omega = 2 * np.pi * self.params.frequency
        I_rms = self.params.primary_current
        N = self.params.primary_turns
        
        if phase == 'A':
            return N * I_rms * np.sqrt(2) * np.cos(omega * time)
        elif phase == 'B':
            return N * I_rms * np.sqrt(2) * np.cos(omega * time - 2*np.pi/3)
        elif phase == 'C':
            return N * I_rms * np.sqrt(2) * np.cos(omega * time + 2*np.pi/3)
        else:
            return 0.0
    
    def solve_no_load_field(self):
        """
        求解空载磁场
        
        空载时只有励磁电流，产生主磁通
        """
        print(f"\n求解空载磁场...")
        
        # 励磁电流 (约为额定电流的2-5%)
        I_excitation = self.params.primary_current * 0.03
        
        # 简化模型：假设磁场主要集中在铁芯内
        mu0 = 4 * np.pi * 1e-7
        mu_core = mu0 * self.params.core_mu_r
        
        # 计算主磁通密度
        # B = μ * N * I / l
        core_length = 2 * self.geom.core_h + 2 * self.geom.yoke_h
        B_main = mu_core * self.params.primary_turns * I_excitation / core_length
        
        print(f"  励磁电流: {I_excitation:.3f} A")
        print(f"  主磁通密度: {B_main:.3f} T")
        
        # 填充磁场分布
        for i in range(self.ny):
            for j in range(self.nx):
                x, y = self.X[i, j], self.Y[i, j]
                
                if self.geom.is_inside_core(x, y):
                    # 铁芯内磁场
                    self.Bx[i, j] = B_main * 0.8
                    self.By[i, j] = B_main * 0.6
                else:
                    # 铁芯外磁场 (漏磁)
                    # 随距离衰减
                    dist = self._distance_to_core(x, y)
                    decay = np.exp(-dist / 0.1)  # 衰减系数
                    self.Bx[i, j] = B_main * 0.1 * decay
                    self.By[i, j] = B_main * 0.1 * decay
        
        self.B_mag = np.sqrt(self.Bx**2 + self.By**2)
        print(f"  最大磁密: {np.max(self.B_mag):.3f} T")
    
    def _distance_to_core(self, x: float, y: float) -> float:
        """计算点到铁芯的最短距离"""
        min_dist = float('inf')
        
        # 检查到铁芯柱的距离
        for i in range(3):
            x_center = self.geom.core_w/2 + i * (self.geom.core_w + self.geom.window_w)
            if (self.geom.yoke_h <= y <= self.geom.yoke_h + self.geom.core_h):
                dist = abs(x - x_center) - self.geom.core_w/2
            else:
                dist = np.sqrt(min((x - (x_center - self.geom.core_w/2))**2,
                                  (x - (x_center + self.geom.core_w/2))**2) +
                              min((y - self.geom.yoke_h)**2,
                                  (y - (self.geom.yoke_h + self.geom.core_h))**2))
            min_dist = min(min_dist, dist)
        
        return max(0, min_dist)
    
    def solve_short_circuit_field(self):
        """
        求解短路工况磁场
        
        短路时主要产生漏磁场
        """
        print(f"\n求解短路工况磁场...")
        
        # 短路电流
        Z_sc = self.params.short_circuit_impedance / 100
        I_sc = self.params.primary_current / Z_sc
        
        print(f"  短路电流: {I_sc:.2f} A")
        
        mu0 = 4 * np.pi * 1e-7
        
        # 漏磁场计算
        for i in range(self.ny):
            for j in range(self.nx):
                x, y = self.X[i, j], self.Y[i, j]
                
                # 找到最近的绕组
                winding_info = self._get_nearest_winding(x, y)
                
                if winding_info['in_winding']:
                    # 绕组内部磁场
                    self.Bx[i, j] = mu0 * self.params.primary_turns * I_sc / self.geom.winding_h * 0.5
                    self.By[i, j] = mu0 * self.params.primary_turns * I_sc / self.geom.winding_t * 0.3
                else:
                    # 绕组外部漏磁场
                    dist = winding_info['distance']
                    B_leak = mu0 * self.params.primary_turns * I_sc / (2 * np.pi * max(dist, 0.01))
                    self.Bx[i, j] = B_leak * 0.4
                    self.By[i, j] = B_leak * 0.3
        
        self.B_mag = np.sqrt(self.Bx**2 + self.By**2)
        print(f"  最大漏磁密: {np.max(self.B_mag):.3f} T")
    
    def _get_nearest_winding(self, x: float, y: float) -> Dict:
        """获取最近的绕组信息"""
        min_dist = float('inf')
        in_winding = False
        
        for region in self.geom.get_winding_regions():
            if (region['x'] <= x <= region['x'] + region['width'] and
                region['y'] <= y <= region['y'] + region['height']):
                in_winding = True
                min_dist = 0
                break
            else:
                # 计算到绕组边界的距离
                dx = max(region['x'] - x, 0, x - (region['x'] + region['width']))
                dy = max(region['y'] - y, 0, y - (region['y'] + region['height']))
                dist = np.sqrt(dx**2 + dy**2)
                min_dist = min(min_dist, dist)
        
        return {'in_winding': in_winding, 'distance': min_dist}
    
    def calculate_core_loss(self) -> float:
        """
        计算铁芯损耗
        
        使用Steinmetz方程简化模型
        """
        # 铁芯体积
        core_volume = (3 * self.geom.core_w * self.geom.core_h * self.geom.core_d +
                      2 * self.geom.total_width * self.geom.yoke_h * self.geom.core_d)
        
        # 铁芯质量
        core_mass = core_volume * self.params.core_density
        
        # 铁损 (简化计算)
        B_avg = np.mean(self.B_mag[self.B_mag > 0.1]) if np.any(self.B_mag > 0.1) else 0.5
        P_core = self.params.core_loss_coeff * core_mass * (B_avg / 1.0)**2 * \
                 (self.params.frequency / 50)**1.3
        
        return P_core
    
    def calculate_copper_loss(self) -> float:
        """计算铜损"""
        # 绕组电阻
        # 简化计算：假设绕组平均匝长
        mean_turn_length = 2 * (self.geom.core_w + self.geom.core_d) * np.pi
        
        # 一次侧绕组电阻
        R_primary = self.params.copper_resistivity * \
                    self.params.primary_turns * mean_turn_length / \
                    (np.pi * (0.005)**2)  # 假设导线直径5mm
        
        # 二次侧绕组电阻
        R_secondary = self.params.copper_resistivity * \
                      self.params.secondary_turns * mean_turn_length / \
                      (np.pi * (0.025)**2)  # 假设导线直径25mm
        
        # 铜损
        P_cu_primary = R_primary * self.params.primary_current**2
        P_cu_secondary = R_secondary * self.params.secondary_current**2
        
        return P_cu_primary + P_cu_secondary
    
    def calculate_short_circuit_impedance(self) -> Dict:
        """
        计算短路阻抗
        
        基于漏磁场能量计算
        """
        mu0 = 4 * np.pi * 1e-7
        
        # 漏感计算 (简化模型)
        # L_leak = μ0 * N² * h * (δ/a)
        h = self.geom.winding_h
        delta = self.geom.winding_t * 2 + 0.01  # 等效漏磁间隙
        a = self.geom.core_w
        
        L_leak = mu0 * self.params.primary_turns**2 * h * delta / a
        
        # 短路阻抗
        omega = 2 * np.pi * self.params.frequency
        X_sc = omega * L_leak
        
        # 假设短路电阻
        R_sc = 0.1  # 简化值
        
        Z_sc = np.sqrt(R_sc**2 + X_sc**2)
        
        # 标幺值
        Z_base = (self.params.primary_voltage * 1000)**2 / (self.params.rated_power * 1000)
        Z_sc_pu = Z_sc / Z_base * 100  # 百分比
        
        return {
            'R_sc': R_sc,
            'X_sc': X_sc,
            'Z_sc': Z_sc,
            'Z_sc_percent': Z_sc_pu,
            'L_leak': L_leak
        }


# ============================================================
# 第四部分：可视化
# ============================================================

def create_visualizations(solver: TransformerFieldSolver, 
                         geom: TransformerGeometry,
                         params: TransformerParameters):
    """创建综合可视化"""
    
    print("\n" + "="*70)
    print("生成可视化")
    print("="*70)
    
    # 1. 变压器几何结构
    fig, ax = plt.subplots(figsize=(14, 10))
    
    # 绘制铁芯
    core_regions = geom.get_core_regions()
    for region in core_regions:
        if region['type'] == 'core_leg':
            rect = Rectangle((region['x']*1000, region['y']*1000), 
                           region['width']*1000, region['height']*1000,
                           facecolor='saddlebrown', edgecolor='black', linewidth=2,
                           label=f"铁芯柱-{region['phase']}" if region['phase'] == 'A' else "")
        else:
            rect = Rectangle((region['x']*1000, region['y']*1000), 
                           region['width']*1000, region['height']*1000,
                           facecolor='saddlebrown', edgecolor='black', linewidth=2,
                           label='铁轭' if region['position'] == 'top' else "")
        ax.add_patch(rect)
    
    # 绘制绕组
    winding_regions = geom.get_winding_regions()
    colors = {'primary_winding': 'orange', 'secondary_winding': 'lightblue'}
    labels = {'primary_winding': '一次侧绕组', 'secondary_winding': '二次侧绕组'}
    
    for i, region in enumerate(winding_regions):
        label = labels[region['type']] if region['phase'] == 'A' else ""
        rect = Rectangle((region['x']*1000, region['y']*1000), 
                        region['width']*1000, region['height']*1000,
                        facecolor=colors[region['type']], 
                        edgecolor='black', linewidth=1.5,
                        alpha=0.7, label=label)
        ax.add_patch(rect)
    
    ax.set_xlim(-20, geom.total_width*1000 + 20)
    ax.set_ylim(-20, geom.total_height*1000 + 20)
    ax.set_aspect('equal')
    ax.set_xlabel('x (mm)')
    ax.set_ylabel('y (mm)')
    ax.set_title('三相电力变压器几何结构')
    ax.legend(loc='upper right')
    ax.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('fig1_transformer_geometry.png', dpi=150, bbox_inches='tight')
    print("  fig1_transformer_geometry.png 已保存")
    plt.close()
    
    # 2. 空载磁场分布
    solver.solve_no_load_field()
    
    fig, axes = plt.subplots(2, 2, figsize=(14, 12))
    
    # 磁感应强度幅值
    ax = axes[0, 0]
    im = ax.pcolormesh(solver.X*1000, solver.Y*1000, solver.B_mag, 
                       shading='auto', cmap='jet', vmin=0, vmax=2.0)
    ax.set_aspect('equal')
    ax.set_xlabel('x (mm)')
    ax.set_ylabel('y (mm)')
    ax.set_title('空载磁感应强度幅值')
    plt.colorbar(im, ax=ax, label='B (T)')
    
    # Bx分量
    ax = axes[0, 1]
    im = ax.pcolormesh(solver.X*1000, solver.Y*1000, solver.Bx, 
                       shading='auto', cmap='RdBu_r', vmin=-1.5, vmax=1.5)
    ax.set_aspect('equal')
    ax.set_xlabel('x (mm)')
    ax.set_ylabel('y (mm)')
    ax.set_title('Bx分量')
    plt.colorbar(im, ax=ax, label='Bx (T)')
    
    # By分量
    ax = axes[1, 0]
    im = ax.pcolormesh(solver.X*1000, solver.Y*1000, solver.By, 
                       shading='auto', cmap='RdBu_r', vmin=-1.5, vmax=1.5)
    ax.set_aspect('equal')
    ax.set_xlabel('x (mm)')
    ax.set_ylabel('y (mm)')
    ax.set_title('By分量')
    plt.colorbar(im, ax=ax, label='By (T)')
    
    # 磁场矢量图
    ax = axes[1, 1]
    skip = 4
    ax.quiver(solver.X[::skip, ::skip]*1000, solver.Y[::skip, ::skip]*1000,
              solver.Bx[::skip, ::skip], solver.By[::skip, ::skip],
              solver.B_mag[::skip, ::skip], cmap='jet', scale=50)
    ax.set_aspect('equal')
    ax.set_xlabel('x (mm)')
    ax.set_ylabel('y (mm)')
    ax.set_title('磁场矢量图')
    ax.set_xlim(0, geom.total_width*1000)
    ax.set_ylim(0, geom.total_height*1000)
    
    plt.tight_layout()
    plt.savefig('fig2_no_load_field.png', dpi=150, bbox_inches='tight')
    print("  fig2_no_load_field.png 已保存")
    plt.close()
    
    # 3. 短路漏磁场
    solver.solve_short_circuit_field()
    
    fig, axes = plt.subplots(1, 2, figsize=(14, 6))
    
    # 漏磁场分布
    ax = axes[0]
    im = ax.pcolormesh(solver.X*1000, solver.Y*1000, solver.B_mag, 
                       shading='auto', cmap='plasma', vmin=0, vmax=1.0)
    ax.set_aspect('equal')
    ax.set_xlabel('x (mm)')
    ax.set_ylabel('y (mm)')
    ax.set_title('短路漏磁场分布')
    plt.colorbar(im, ax=ax, label='B (T)')
    
    # 漏磁场矢量
    ax = axes[1]
    skip = 4
    ax.quiver(solver.X[::skip, ::skip]*1000, solver.Y[::skip, ::skip]*1000,
              solver.Bx[::skip, ::skip], solver.By[::skip, ::skip],
              solver.B_mag[::skip, ::skip], cmap='plasma', scale=20)
    ax.set_aspect('equal')
    ax.set_xlabel('x (mm)')
    ax.set_ylabel('y (mm)')
    ax.set_title('漏磁场矢量图')
    ax.set_xlim(0, geom.total_width*1000)
    ax.set_ylim(0, geom.total_height*1000)
    
    plt.tight_layout()
    plt.savefig('fig3_leakage_field.png', dpi=150, bbox_inches='tight')
    print("  fig3_leakage_field.png 已保存")
    plt.close()


# ============================================================
# 第五部分：主程序
# ============================================================

if __name__ == "__main__":
    print("="*70)
    print("主题026: 变压器电磁场仿真")
    print("三相电力变压器电磁场分析")
    print("="*70)
    
    # 1. 定义变压器参数
    print("\n[1/5] 初始化变压器参数...")
    params = TransformerParameters()
    
    # 2. 创建几何模型
    print("\n[2/5] 创建几何模型...")
    geometry = TransformerGeometry(params)
    
    # 3. 创建磁场求解器
    print("\n[3/5] 初始化磁场求解器...")
    solver = TransformerFieldSolver(geometry, params)
    
    # 4. 求解并可视化
    print("\n[4/5] 求解电磁场...")
    create_visualizations(solver, geometry, params)
    
    # 5. 性能计算
    print("\n[5/5] 计算性能参数...")
    
    # 空载损耗
    solver.solve_no_load_field()
    core_loss = solver.calculate_core_loss()
    print(f"\n空载铁损: {core_loss:.2f} W")
    
    # 铜损
    copper_loss = solver.calculate_copper_loss()
    print(f"额定铜损: {copper_loss:.2f} W")
    
    # 短路阻抗
    sc_results = solver.calculate_short_circuit_impedance()
    print(f"\n短路阻抗计算:")
    print(f"  短路电阻: {sc_results['R_sc']:.4f} Ω")
    print(f"  短路电抗: {sc_results['X_sc']:.4f} Ω")
    print(f"  短路阻抗: {sc_results['Z_sc']:.4f} Ω")
    print(f"  短路阻抗百分比: {sc_results['Z_sc_percent']:.2f}%")
    print(f"  漏感: {sc_results['L_leak']*1000:.4f} mH")
    
    # 效率估算
    total_loss = core_loss + copper_loss
    efficiency = params.rated_power * 1000 / (params.rated_power * 1000 + total_loss) * 100
    print(f"\n效率估算: {efficiency:.2f}%")
    
    # 保存结果
    results = {
        'transformer_params': {
            'rated_power_kva': params.rated_power,
            'primary_voltage_kv': params.primary_voltage,
            'secondary_voltage_kv': params.secondary_voltage,
            'turns_ratio': params.turns_ratio
        },
        'performance': {
            'core_loss_w': core_loss,
            'copper_loss_w': copper_loss,
            'total_loss_w': total_loss,
            'efficiency_percent': efficiency,
            'short_circuit_impedance_percent': sc_results['Z_sc_percent']
        },
        'geometry': {
            'total_width_mm': geometry.total_width * 1000,
            'total_height_mm': geometry.total_height * 1000
        }
    }
    
    with open('simulation_results.json', 'w', encoding='utf-8') as f:
        json.dump(results, f, indent=2, ensure_ascii=False)
    
    print("\n" + "="*70)
    print("主题026仿真完成！")
    print("="*70)
    print("\n生成的文件:")
    print("  - fig1_transformer_geometry.png: 变压器几何结构")
    print("  - fig2_no_load_field.png: 空载磁场分布")
    print("  - fig3_leakage_field.png: 短路漏磁场")
    print("  - simulation_results.json: 仿真结果数据")
    print("="*70)