第二十四篇：电-热耦合分析

摘要

电-热耦合是多物理场耦合仿真中的重要研究领域，涉及电能向热能的转换过程及其热效应对电学性能的影响。本主题系统阐述电-热耦合的基本物理机制，包括焦耳热效应、电阻温度特性、热电效应等核心概念。详细推导电-热耦合的控制方程组，建立电场与温度场的双向耦合数学模型。介绍顺序耦合与直接耦合两种求解策略，讨论电阻加热、电子器件热管理、热电转换等典型应用场景。通过六个工程案例的Python仿真实现，展示电-热耦合问题的数值求解方法，包括稳态与瞬态分析、非线性材料特性处理、多尺度建模等技术。最后探讨电-热耦合在新能源、微电子、电力系统等领域的应用前景与发展趋势。

关键词

电-热耦合，焦耳热，电阻加热，热电效应，温度依赖性，非线性分析，电子散热，有限元方法

---

1. 引言

1.1 电-热耦合现象概述

电-热耦合是指电场与温度场之间的相互作用现象，是自然界和工程应用中广泛存在的多物理场耦合问题。当电流通过导体时，由于电阻的存在，电能会转化为热能，这就是著名的焦耳热效应。反过来，温度变化会影响材料的电学性质，如电阻率、介电常数等，从而改变电场分布。这种双向耦合机制在电力电子、微电子器件、新能源系统等领域具有重要影响。

电-热耦合问题的复杂性主要体现在以下几个方面：

非线性特性：材料的电阻率通常随温度变化，导致电场方程呈现非线性特征。例如，金属的电阻率随温度升高而增大，而半导体的电阻率则随温度升高而减小。这种温度依赖性使得电-热耦合问题本质上是非线性的，需要迭代求解。

多时间尺度：电场传播的时间尺度极短（电磁波传播），而热传导的时间尺度相对较长（热扩散）。在大多数工程应用中，可以采用准静态近似，忽略电磁波的传播效应，仅考虑焦耳热的产生和热传导过程。

多空间尺度：从纳米级的微电子器件到千米级的输电线路，电-热耦合问题涉及多个数量级的空间尺度。不同尺度下，主导物理机制可能不同，需要采用不同的建模方法。

材料多样性：导体、半导体、绝缘体、复合材料等不同类型材料的电-热特性差异显著。超导材料在临界温度以下电阻为零，不产生焦耳热；而正温度系数（PTC）材料在特定温度范围内电阻急剧增大，具有自限温特性。

1.2 工程应用背景

电-热耦合分析在众多工程领域具有重要应用价值：

电力电子器件：IGBT、MOSFET等功率半导体器件在工作时产生大量热量，需要有效的散热设计。电-热耦合分析可以预测器件结温，评估热可靠性，优化散热结构。

集成电路：芯片集成度不断提高，功率密度持续增大，热管理成为制约性能提升的关键因素。电-热耦合仿真可以分析芯片温度分布，识别热点位置，指导布局优化。

电阻加热设备：电热水器、电炉、电阻焊接等设备利用焦耳热效应实现加热功能。电-热耦合分析可以优化加热元件设计，提高能量转换效率，实现温度精确控制。

热电转换系统：热电发电机利用塞贝克效应将热能直接转换为电能，热电制冷器利用珀尔帖效应实现固态制冷。电-热耦合分析是热电系统性能评估和优化的基础。

电力系统：输电线路、变压器、电缆等电力设备在运行中因电阻损耗而发热。电-热耦合分析可以评估设备载流能力，预测温升，确保运行安全。

新能源领域：锂离子电池、燃料电池等电化学器件存在显著的电-热耦合效应。热失控是锂离子电池安全性的主要问题，电-热耦合分析对于电池热管理设计至关重要。

1.3 本主题内容安排

本主题将系统介绍电-热耦合分析的理论基础和数值方法，具体安排如下：

第2节阐述电-热耦合的物理机制，包括焦耳热效应、电阻温度特性、热电效应等；第3节建立电-热耦合的数学模型，推导控制方程组；第4节介绍数值求解方法，包括有限元离散、耦合策略、非线性处理等；第5节通过六个工程案例展示电-热耦合分析的Python实现；第6节讨论结果分析与工程应用；第7节总结全文并展望未来发展方向。

---

2. 电-热耦合物理机制

2.1 焦耳热效应

焦耳热效应是电-热耦合最基本的物理机制，描述了电流通过导体时电能向热能的转换过程。这一现象由英国物理学家詹姆斯·普雷斯科特·焦耳于1841年通过实验发现。

2.1.1 焦耳定律

焦耳定律指出，电流通过导体产生的热量与电流的平方、导体电阻以及通电时间成正比：

$$Q=I^{2}Rt$$

其中，$Q$为产生的热量（J），$I$为电流（A），$R$为电阻（Ω），$t$为时间（s）。

从功率角度，焦耳热功率密度可以表示为：

$$q^{\prime \prime \prime }=\mathbf{J}\cdot \mathbf{E}=\sigma |\mathbf{E}{|}^2=\frac{|\mathbf{J}{|}^2}{\sigma }$$

其中，$q^{\prime \prime \prime }$为体积热生成率（W/m³），$\mathbf{J}$为电流密度（A/m²），$\mathbf{E}$为电场强度（V/m），$\sigma$为电导率（S/m）。

这个公式表明，热功率密度与电场强度的平方成正比，与电导率成反比。对于给定电流密度的情况，电阻率越大的材料产生的热量越多。

2.1.2 微观机制

从微观角度看，焦耳热源于载流子（电子或离子）与晶格原子的碰撞。在电场作用下，载流子获得定向漂移速度，同时通过与晶格振动（声子）的碰撞将动能传递给晶格，表现为热能的产生。

在金属导体中，自由电子在电场作用下运动，与晶格中的离子实发生碰撞，将能量传递给晶格，导致温度升高。在半导体中，电子和空穴都参与导电，两者都对焦耳热有贡献。在电解质溶液中，离子的定向迁移也会产生焦耳热。

2.1.3 能量守恒视角

从能量守恒的角度，焦耳热效应体现了电能向热能的转换。在稳态情况下，输入的电功率等于输出的热功率：

$$P_{elec}={\int }_V\mathbf{J}\cdot \mathbf{E}dV=P_{heat}$$

这一能量平衡关系是电-热耦合分析的基础，确保了数值计算的能量守恒。

2.2 电阻温度特性

材料的电阻率通常随温度变化，这种温度依赖性是电-热双向耦合的重要机制。温度变化改变电阻率，进而影响电流分布和热生成率。

2.2.1 金属导体的温度特性

对于大多数金属，电阻率随温度升高而线性增大：

$$\rho (T)={\rho }_0[1+\alpha (T-T_0)]$$

其中，$\rho (T)$为温度$T$时的电阻率，${\rho }_0$为参考温度$T_0$时的电阻率，$\alpha$为电阻温度系数（1/K）。

典型金属的电阻温度系数：

• 铜：$\alpha \approx 0.0039$ K⁻¹
• 铝：$\alpha \approx 0.0043$ K⁻¹
• 铁：$\alpha \approx 0.0050$ K⁻¹
• 钨：$\alpha \approx 0.0045$ K⁻¹

这种正温度系数特性意味着金属导体在温度升高时电阻增大，在恒定电压条件下电流减小，热生成率降低，形成一定的自稳定效应。

2.2.2 半导体的温度特性

半导体的电阻率随温度升高而减小，表现出负温度系数特性：

$$\rho (T)={\rho }_0\exp \left(\frac{E_g}{2k_{B}T}\right)$$

其中，$E_g$为禁带宽度，$k_B$为玻尔兹曼常数。

在较低温度下，本征激发较弱，电阻率主要由掺杂浓度决定。随着温度升高，本征载流子浓度指数增加，电阻率迅速下降。这种特性使得半导体器件在高温下更容易产生热失控。

2.2.3 正温度系数材料

正温度系数（PTC）材料在特定温度范围内电阻率急剧增大，常用于自限温加热器和过流保护器件。典型的PTC材料包括掺杂的钛酸钡陶瓷。

PTC材料的电阻率-温度关系可以表示为：

$$\rho (T)={\rho }_0\exp \left[A(T-T_c)\right],T>T_c$$

其中，$T_c$为居里温度，$A$为材料常数。在居里温度以上，电阻率随温度指数增长，实现自限温功能。

2.3 热电效应

热电效应描述了温度梯度与电势之间的相互转换关系，包括塞贝克效应、珀尔帖效应和汤姆逊效应。

2.3.1 塞贝克效应

塞贝克效应（1821年发现）指出，当两种不同导体构成闭合回路且两个接点处于不同温度时，回路中会产生电动势。塞贝克系数定义为：

$$S=\frac{\Delta V}{\Delta T}$$

其中，$S$为塞贝克系数（V/K），$\Delta V$为产生的电压，$\Delta T$为温差。

塞贝克效应是热电发电的基础。热电发电机的效率由优值系数ZT决定：

$$ZT=\frac{S^2\sigma T}{k}$$

其中，$k$为热导率。ZT值越大，热电转换效率越高。目前商用热电材料的ZT值约为1-2，实验室材料可达2-3。

2.3.2 珀尔帖效应

珀尔帖效应（1834年发现）是塞贝克效应的逆效应。当电流通过两种不同导体的接点时，接点处会吸热或放热，热量与电流成正比：

$$Q_{Peltier}={\Pi }_{AB}It=(S_B-S_A)TIt$$

其中，${\Pi }_{AB}$为珀尔帖系数。

珀尔帖效应用于固态制冷（热电制冷器），具有无运动部件、无制冷剂、精确控温等优点，但效率通常低于压缩式制冷。

2.3.3 汤姆逊效应

汤姆逊效应（1851年发现）描述了单一导体中存在温度梯度且有电流通过时的吸热或放热现象。汤姆逊热功率为：

$$q_{Thomson}^{\prime \prime \prime }=\tau \mathbf{J}\cdot \nabla T$$

其中，$\tau$为汤姆逊系数。汤姆逊效应通常比焦耳热和珀尔帖效应小得多，在一般工程分析中可以忽略。

2.4 其他电-热耦合机制

2.4.1 介电损耗

在交变电场中，电介质材料会因极化弛豫而产生热量，称为介电损耗。介电损耗功率密度为：

$$q_{dielectric}^{\prime \prime \prime }=\omega {\epsilon }_0{\epsilon }^{\prime \prime }|\mathbf{E}{|}^2$$

其中，$\omega$为角频率，${\epsilon }_0$为真空介电常数，${\epsilon }^{\prime \prime }$为损耗因子。

介电加热广泛应用于微波加热、射频加热、高频焊接等工艺。

2.4.2 磁滞损耗

铁磁材料在交变磁场中会因磁畴壁移动和磁矩转动而产生热量，称为磁滞损耗。磁滞损耗功率与磁滞回线面积成正比：

$$P_{hysteresis}=f\oint HdB$$

其中，$f$为磁场频率。磁滞损耗是变压器、电机等电磁设备的主要热源之一。

2.4.3 涡流损耗

导体在交变磁场中会感应出涡流，产生焦耳热。涡流损耗功率为：

$$P_{eddy}=\frac{{\pi }^2{f}^2{B}_{max}^2{d}^2}{6\rho }V$$

其中，$d$为导体厚度，$B_{max}$为最大磁感应强度。为减小涡流损耗，电机铁芯通常采用叠片结构。

---

3. 电-热耦合数学模型

3.1 基本控制方程

电-热耦合问题由电场方程和温度场方程组成，两者通过焦耳热项和电阻温度特性相互耦合。

3.1.1 电场方程

在准静态近似下，电场满足电流连续性方程和欧姆定律：

$$\nabla \cdot \mathbf{J}=0$$

$$\mathbf{J}=\sigma (T)\mathbf{E}=-\sigma (T)\nabla \varphi$$

其中，$\varphi$为电势。合并得到电势方程：

$$\nabla \cdot [\sigma (T)\nabla \varphi ]=0$$

这是一个椭圆型偏微分方程，当电导率随温度变化时呈现非线性特征。

边界条件通常包括：

• 狄利克雷边界：指定电势$\varphi ={\varphi }_0$
• 诺伊曼边界：指定电流密度$\mathbf{J}\cdot \mathbf{n}=J_n$
• 混合边界：$\alpha \varphi +\beta \mathbf{J}\cdot \mathbf{n}=\gamma$

3.1.2 温度场方程

温度场满足能量守恒方程（热传导方程）：

$$\rho c_p\frac{\partial T}{\partial t}=\nabla \cdot (k\nabla T)+q^{\prime \prime \prime }$$

其中，$\rho$为密度，$c_p$为比热容，$k$为热导率，$q^{\prime \prime \prime }$为体积热生成率。

对于电-热耦合问题，热源项包括焦耳热：

$$q^{\prime \prime \prime }=\mathbf{J}\cdot \mathbf{E}=\sigma (T)|\nabla \varphi {|}^2$$

边界条件包括：

• 指定温度：$T=T_0$
• 指定热流：$-k\nabla T\cdot \mathbf{n}=q_n$
• 对流换热：$-k\nabla T\cdot \mathbf{n}=h(T-T_{\infty })$
• 辐射换热：$-k\nabla T\cdot \mathbf{n}=\epsilon {\sigma }_{SB}(T^4-T_{\infty }^4)$

3.1.3 耦合方程组

综合电场和温度场方程，得到电-热耦合控制方程组：

$$\{\begin{array}{ll}\nabla \cdot [\sigma (T)\nabla \varphi ]=0& \text{(电场方程)}\\ \rho c_p\frac{\partial T}{\partial t}=\nabla \cdot (k\nabla T)+\sigma (T)|\nabla \varphi {|}^2& \text{(温度场方程)}\end{array}$$

这是一个非线性耦合方程组，电场方程通过$\sigma (T)$依赖于温度，温度场方程通过焦耳热源依赖于电场。

3.2 材料本构关系

3.2.1 电导率模型

电导率与温度的关系是电-热耦合的核心。常用模型包括：

线性模型（金属）：
$$\sigma (T)=\frac{{\sigma }_0}{1+\alpha (T-T_0)}$$

指数模型（半导体）：
$$\sigma (T)={\sigma }_0\exp \left[-\frac{E_g}{2k_B}\left(\frac{1}{T}-\frac{1}{T_0}\right)\right]$$

多项式模型（经验拟合）：
$$\sigma (T)={\sigma }_0[1+a_1(T-T_0)+a_2(T-T_0{)}^2+\cdots ]$$

3.2.2 热物性参数

热导率、比热容等参数也可能随温度变化：

$$k(T)=k_0[1+{\beta }_k(T-T_0)]$$

$$c_p(T)=c_{p0}[1+{\beta }_c(T-T_0)]$$

对于固体材料，这些参数的温度依赖性通常比电导率弱，在初步分析中可以视为常数。

3.3 耦合类型与求解策略

3.3.1 弱耦合（顺序耦合）

弱耦合策略交替求解电场和温度场方程：

1. 假设初始温度分布$T^{(0)}$
2. 求解电场方程得到${\varphi }^{(n)}$，计算焦耳热$q^{\prime \prime \prime (n)}=\sigma (T^{(n-1)})|\nabla {\varphi }^{(n)}{|}^2$
3. 求解温度场方程得到$T^{(n)}$
4. 检查收敛性，若不收敛返回步骤2

收敛判据：
$$\Vert T^{(n)}-T^{(n-1)}\Vert <{\epsilon }_T,\Vert {\varphi }^{(n)}-{\varphi }^{(n-1)}\Vert <{\epsilon }_{\varphi }$$

弱耦合适用于耦合强度较弱的问题，计算成本较低，但收敛速度可能较慢。

3.3.2 强耦合（直接耦合）

强耦合策略同时求解电场和温度场方程，形成统一的代数方程组：

$$\left[\begin{array}{cc}{\mathbf{K}}_{\varphi }& 0\\ \mathbf{Q}& {\mathbf{K}}_T\end{array}\right]\left[\begin{array}{c}\mathbf{\varphi }\\ \mathbf{T}\end{array}\right]=\left[\begin{array}{c}{\mathbf{f}}_{\varphi }\\ {\mathbf{f}}_T\end{array}\right]$$

其中，${\mathbf{K}}_{\varphi }$和${\mathbf{K}}_T$分别为电场和温度场的刚度矩阵，$\mathbf{Q}$为耦合矩阵。

强耦合需要处理非线性，通常采用牛顿-拉夫森迭代：

$${\mathbf{J}}^{(k)}\Delta {\mathbf{u}}^{(k)}=-{\mathbf{R}}^{(k)}$$

其中，$\mathbf{J}$为雅可比矩阵，$\mathbf{R}$为残差向量，$\mathbf{u}=[\mathbf{\varphi },\mathbf{T}{]}^T$。

强耦合收敛速度快，适用于强耦合问题，但计算成本和内存需求较大。

---

4. 数值求解方法

4.1 有限元离散

4.1.1 电场方程的弱形式

将电势方程乘以测试函数$w$并积分：

$${\int }_{\Omega }w\nabla \cdot (\sigma \nabla \varphi )d\Omega =0$$

分部积分得到弱形式：

$${\int }_{\Omega }\sigma \nabla w\cdot \nabla \varphi d\Omega ={\int }_{\Gamma }w\sigma \frac{\partial \varphi }{\partial n}d\Gamma$$

4.1.2 温度场方程的弱形式

类似地，温度场方程的弱形式为：

$${\int }_{\Omega }w\rho c_p\frac{\partial T}{\partial t}d\Omega +{\int }_{\Omega }k\nabla w\cdot \nabla Td\Omega ={\int }_{\Omega }w{q}^{\prime \prime \prime }d\Omega +{\int }_{\Gamma }wk\frac{\partial T}{\partial n}d\Gamma$$

4.1.3 有限元离散

采用伽辽金方法，设近似解：

$$\varphi \approx \sum _{j=1}^N{N}_j{\varphi }_j,T\approx \sum _{j=1}^N{N}_j{T}_j$$

其中，$N_j$为形函数，${\varphi }_j$和$T_j$为节点值。代入弱形式得到离散方程：

$${\mathbf{K}}_{\varphi }\mathbf{\varphi }={\mathbf{f}}_{\varphi }$$

$$\mathbf{M}\frac{d\mathbf{T}}{dt}+{\mathbf{K}}_T\mathbf{T}={\mathbf{f}}_T$$

其中：

$$K_{\varphi ,ij}={\int }_{\Omega }\sigma \nabla N_i\cdot \nabla N_{j}d\Omega$$

$$M_{ij}={\int }_{\Omega }\rho c_p{N}_i{N}_{j}d\Omega$$

$$K_{T,ij}={\int }_{\Omega }k\nabla N_i\cdot \nabla N_{j}d\Omega$$

$$f_{T,i}={\int }_{\Omega }{N}_i{q}^{\prime \prime \prime }d\Omega$$

4.2 时间积分

对于瞬态问题，采用时间离散方法：

向后欧拉法（一阶精度，无条件稳定）：
$$\mathbf{M}\frac{{\mathbf{T}}^{n+1}-{\mathbf{T}}^n}{\Delta t}+{\mathbf{K}}_T{\mathbf{T}}^{n+1}={\mathbf{f}}_T^{n+1}$$

Crank-Nicolson法（二阶精度，无条件稳定）：
$$\mathbf{M}\frac{{\mathbf{T}}^{n+1}-{\mathbf{T}}^n}{\Delta t}+{\mathbf{K}}_T\frac{{\mathbf{T}}^{n+1}+{\mathbf{T}}^n}{2}=\frac{{\mathbf{f}}_T^{n+1}+{\mathbf{f}}_T^n}{2}$$

4.3 非线性求解

由于电导率的温度依赖性，电-热耦合问题通常是非线性的，需要迭代求解。

4.3.1 不动点迭代

$${\mathbf{K}}_{\varphi }(T^{(k)}){\mathbf{\varphi }}^{(k+1)}={\mathbf{f}}_{\varphi }$$

$$q^{\prime \prime \prime (k+1)}=\sigma (T^{(k)})|\nabla {\varphi }^{(k+1)}{|}^2$$

$$\mathbf{M}\frac{d{\mathbf{T}}^{(k+1)}}{dt}+{\mathbf{K}}_T{\mathbf{T}}^{(k+1)}={\mathbf{f}}_T(q^{\prime \prime \prime (k+1)})$$

简单但收敛慢，适用于弱非线性问题。

4.3.2 牛顿-拉夫森迭代

构造残差向量：

$${\mathbf{R}}_{\varphi }={\mathbf{K}}_{\varphi }(T)\mathbf{\varphi }-{\mathbf{f}}_{\varphi }$$

$${\mathbf{R}}_T=\mathbf{M}\frac{d\mathbf{T}}{dt}+{\mathbf{K}}_T\mathbf{T}-{\mathbf{f}}_T(\mathbf{\varphi },T)$$

雅可比矩阵：

$$\mathbf{J}=\left[\begin{array}{cc}\frac{\partial {\mathbf{R}}_{\varphi }}{\partial \mathbf{\varphi }}& \frac{\partial {\mathbf{R}}_{\varphi }}{\partial \mathbf{T}}\\ \frac{\partial {\mathbf{R}}_T}{\partial \mathbf{\varphi }}& \frac{\partial {\mathbf{R}}_T}{\partial \mathbf{T}}\end{array}\right]$$

迭代求解：

$${\mathbf{J}}^{(k)}{\left[\begin{array}{c}\Delta \mathbf{\varphi }\\ \Delta \mathbf{T}\end{array}\right]}^{(k)}=-{\left[\begin{array}{c}{\mathbf{R}}_{\varphi }\\ {\mathbf{R}}_T\end{array}\right]}^{(k)}$$

收敛快，但需要计算雅可比矩阵，实现较复杂。

---

5. Python仿真实现

5.1 案例1：电阻丝的稳态电-热分析

问题描述：分析通电电阻丝的温度分布。电阻丝长度为$L$，直径为$d$，两端施加电压$V$，环境温度为$T_{\infty }$，表面与空气进行对流换热。

控制方程：

$$\frac{d}{dx}\left(\sigma A\frac{d\varphi }{dx}\right)=0$$

$$\frac{d}{dx}\left(kA\frac{dT}{dx}\right)+\sigma A{\left(\frac{d\varphi }{dx}\right)}^2-hP(T-T_{\infty })=0$$

其中，$A$为截面积，$P$为周长，$h$为对流换热系数。

# -*- coding: utf-8 -*-
"""
主题024：电-热耦合分析 - 案例1：电阻丝的稳态电-热分析
"""
import matplotlib
matplotlib.use('Agg')
import numpy as np
import matplotlib.pyplot as plt
from scipy.sparse import diags
from scipy.sparse.linalg import spsolve
import warnings
warnings.filterwarnings('ignore')
import os

output_dir = r'd:\文档\500仿真领域\工程仿真\多物理场耦合仿真\主题024'
os.makedirs(output_dir, exist_ok=True)

plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

print("="*60)
print("主题024：电-热耦合分析")
print("="*60)

print("\n" + "="*60)
print("案例1：电阻丝的稳态电-热分析")
print("="*60)

# 物理参数
L = 0.1  # 电阻丝长度 (m)
d = 0.001  # 直径 (m)
V_applied = 10.0  # 施加电压 (V)
T_ambient = 300.0  # 环境温度 (K)

# 材料参数（镍铬合金）
rho_0 = 1.10e-6  # 电阻率 (Ohm·m)
alpha = 0.0004  # 电阻温度系数 (1/K)
k = 15.0  # 热导率 (W/m·K)
sigma_SB = 5.67e-8  # 斯特藩-玻尔兹曼常数

# 几何参数
A = np.pi * (d/2)**2  # 截面积
P = np.pi * d  # 周长

# 对流换热系数
h_conv = 25.0  # W/m²·K

# 网格参数
nx = 100
x = np.linspace(0, L, nx)
dx = x[1] - x[0]

# 初始化
T = np.ones(nx) * T_ambient  # 初始温度分布
phi = np.zeros(nx)  # 电势分布

# 迭代参数
max_iter = 50
tol = 1e-6

print("\n开始电-热耦合迭代...")
for iteration in range(max_iter):
    T_old = T.copy()
    
    # 步骤1: 求解电场（考虑温度依赖的电阻率）
    # 电阻率随温度变化
    rho = rho_0 * (1 + alpha * (T - T_ambient))
    sigma = 1.0 / rho
    
    # 构建电势方程的系数矩阵
    # d/dx(sigma * A * dphi/dx) = 0
    a_center = np.zeros(nx)
    a_east = np.zeros(nx)
    a_west = np.zeros(nx)
    b_phi = np.zeros(nx)
    
    for i in range(1, nx-1):
        sigma_e = 0.5 * (sigma[i] + sigma[i+1])
        sigma_w = 0.5 * (sigma[i] + sigma[i-1])
        
        a_east[i] = sigma_e * A / dx
        a_west[i] = sigma_w * A / dx
        a_center[i] = -(a_east[i] + a_west[i])
    
    # 边界条件
    # x=0: phi = 0 (接地)
    a_center[0] = 1.0
    b_phi[0] = 0.0
    
    # x=L: phi = V_applied
    a_center[-1] = 1.0
    b_phi[-1] = V_applied
    
    # 构建稀疏矩阵
    diagonals = [a_west[1:], a_center, a_east[:-1]]
    offsets = [-1, 0, 1]
    K_phi = diags(diagonals, offsets, format='csr')
    
    # 求解电势
    phi = spsolve(K_phi, b_phi)
    
    # 计算电场强度和电流密度
    E_field = np.zeros(nx)
    for i in range(1, nx-1):
        E_field[i] = -(phi[i+1] - phi[i-1]) / (2*dx)
    E_field[0] = -(phi[1] - phi[0]) / dx
    E_field[-1] = -(phi[-1] - phi[-2]) / dx
    
    J_current = sigma * E_field
    I_total = J_current * A
    
    # 计算焦耳热
    q_joule = sigma * E_field**2  # W/m³
    Q_joule = q_joule * A  # W/m
    
    # 步骤2: 求解温度场
    # d/dx(k*A*dT/dx) + Q_joule - h*P*(T-T_ambient) = 0
    aT_center = np.zeros(nx)
    aT_east = np.zeros(nx)
    aT_west = np.zeros(nx)
    b_T = np.zeros(nx)
    
    for i in range(1, nx-1):
        aT_east[i] = k * A / dx
        aT_west[i] = k * A / dx
        # 添加对流项（隐式处理）
        h_term = h_conv * P * dx
        aT_center[i] = -(aT_east[i] + aT_west[i] + h_term)
        b_T[i] = -Q_joule[i] * dx - h_conv * P * dx * T_ambient
    
    # 边界条件（绝热端部，仅考虑对流）
    aT_center[0] = 1.0
    aT_east[0] = -1.0
    b_T[0] = 0.0  # 对称边界
    
    aT_center[-1] = 1.0
    aT_west[-1] = -1.0
    b_T[-1] = 0.0
    
    # 构建并求解温度场
    diagonals_T = [aT_west[1:], aT_center, aT_east[:-1]]
    K_T = diags(diagonals_T, offsets, format='csr')
    T = spsolve(K_T, b_T)
    
    # 检查收敛
    error = np.max(np.abs(T - T_old))
    if iteration % 10 == 0:
        print(f"  迭代 {iteration}: 最大温度变化 = {error:.6e} K")
    
    if error < tol:
        print(f"  收敛于迭代 {iteration}")
        break

print(f"\n计算结果:")
print(f"  电流: {np.mean(I_total):.4f} A")
print(f"  电阻: {V_applied/np.mean(I_total):.4f} Ohm")
print(f"  最高温度: {np.max(T):.2f} K ({np.max(T)-T_ambient:.2f} K 温升)")
print(f"  平均温度: {np.mean(T):.2f} K")
print(f"  总焦耳热功率: {np.sum(Q_joule)*dx:.4f} W")

# 可视化
fig, axes = plt.subplots(2, 2, figsize=(12, 10))

# 电势分布
ax1 = axes[0, 0]
ax1.plot(x*1000, phi, 'b-', linewidth=2)
ax1.set_xlabel('Position (mm)', fontsize=11)
ax1.set_ylabel('Electric Potential (V)', fontsize=11)
ax1.set_title('Electric Potential Distribution', fontsize=12, fontweight='bold')
ax1.grid(True, alpha=0.3)

# 温度分布
ax2 = axes[0, 1]
ax2.plot(x*1000, T, 'r-', linewidth=2, label='Temperature')
ax2.axhline(y=T_ambient, color='k', linestyle='--', label='Ambient')
ax2.set_xlabel('Position (mm)', fontsize=11)
ax2.set_ylabel('Temperature (K)', fontsize=11)
ax2.set_title('Temperature Distribution', fontsize=12, fontweight='bold')
ax2.legend()
ax2.grid(True, alpha=0.3)

# 电流密度
ax3 = axes[1, 0]
ax3.plot(x*1000, np.abs(J_current)/1e6, 'g-', linewidth=2)
ax3.set_xlabel('Position (mm)', fontsize=11)
ax3.set_ylabel('Current Density (MA/m²)', fontsize=11)
ax3.set_title('Current Density Distribution', fontsize=12, fontweight='bold')
ax3.grid(True, alpha=0.3)

# 焦耳热生成
ax4 = axes[1, 1]
ax4.plot(x*1000, q_joule/1e6, 'm-', linewidth=2)
ax4.set_xlabel('Position (mm)', fontsize=11)
ax4.set_ylabel('Joule Heat Generation (MW/m³)', fontsize=11)
ax4.set_title('Volumetric Heat Generation', fontsize=12, fontweight='bold')
ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig(f'{output_dir}/case1_resistor_heating.png', dpi=150, bbox_inches='tight')
plt.close()
print("\n✓ 案例1完成，图片已保存")

5.2 案例2：芯片散热分析

问题描述：分析集成电路芯片在工作时的温度分布。芯片尺寸为$10\times 10$ mm²，内部有多个热源区域，底部通过导热硅脂连接到散热器。

控制方程：二维稳态热传导方程：

$$\frac{\partial }{\partial x}\left(k\frac{\partial T}{\partial x}\right)+\frac{\partial }{\partial y}\left(k\frac{\partial T}{\partial y}\right)+q^{\prime \prime \prime }(x,y)=0$$

print("\n" + "="*60)
print("案例2：芯片散热分析")
print("="*60)

# 芯片几何参数
chip_size = 10e-3  # 芯片边长 (m)
H_chip = 0.5e-3  # 芯片厚度 (m)

# 网格
nx, ny = 50, 50
x = np.linspace(0, chip_size, nx)
y = np.linspace(0, chip_size, ny)
dx = x[1] - x[0]
dy = y[1] - y[0]
X, Y = np.meshgrid(x, y)

# 材料参数（硅）
k_chip = 150.0  # W/m·K

# 定义热源分布（模拟CPU核心布局）
q_chip = np.zeros((ny, nx))

# 核心1 - 左下角
core1_mask = ((X > 1e-3) & (X < 4e-3) & (Y > 1e-3) & (Y < 4e-3))
q_chip[core1_mask] = 5e8  # W/m³

# 核心2 - 右下角
core2_mask = ((X > 6e-3) & (X < 9e-3) & (Y > 1e-3) & (Y < 4e-3))
q_chip[core2_mask] = 5e8

# 核心3 - 左上角
core3_mask = ((X > 1e-3) & (X < 4e-3) & (Y > 6e-3) & (Y < 9e-3))
q_chip[core3_mask] = 5e8

# 核心4 - 右上角
core4_mask = ((X > 6e-3) & (X < 9e-3) & (Y > 6e-3) & (Y < 9e-3))
q_chip[core4_mask] = 5e8

# GPU核心 - 中间
gpu_mask = ((X > 3.5e-3) & (X < 6.5e-3) & (Y > 3.5e-3) & (Y < 6.5e-3))
q_chip[gpu_mask] = 8e8

# 边界条件
h_top = 1000.0  # 散热器换热系数
T_ambient = 300.0  # K

# 有限差分求解
T_chip = np.ones((ny, nx)) * T_ambient

max_iter = 5000
tol = 1e-5

print("开始求解芯片温度场...")
for iteration in range(max_iter):
    T_old = T_chip.copy()
    
    # 内部节点
    for i in range(1, ny-1):
        for j in range(1, nx-1):
            T_chip[i, j] = (k_chip * (T_chip[i+1, j] + T_chip[i-1, j]) / dy**2 +
                           k_chip * (T_chip[i, j+1] + T_chip[i, j-1]) / dx**2 +
                           q_chip[i, j]) / (2 * k_chip * (1/dx**2 + 1/dy**2))
    
    # 边界条件
    # 顶部：对流换热
    T_chip[0, :] = (k_chip * T_chip[1, :] / dy + h_top * T_ambient) / (k_chip/dy + h_top)
    # 底部：绝热（对称）
    T_chip[-1, :] = T_chip[-2, :]
    # 左右：绝热
    T_chip[:, 0] = T_chip[:, 1]
    T_chip[:, -1] = T_chip[:, -2]
    
    # 检查收敛
    if iteration % 500 == 0:
        error = np.max(np.abs(T_chip - T_old))
        print(f"  迭代 {iteration}: 最大变化 = {error:.6e}")
    
    if np.max(np.abs(T_chip - T_old)) < tol:
        print(f"  收敛于迭代 {iteration}")
        break

print(f"\n芯片温度分析结果:")
print(f"  最高温度: {np.max(T_chip):.2f} K ({np.max(T_chip)-T_ambient:.2f} K 温升)")
print(f"  平均温度: {np.mean(T_chip):.2f} K")
print(f"  总热功率: {np.sum(q_chip)*dx*dy*H_chip:.4f} W")

# 可视化
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# 温度分布
ax1 = axes[0]
im1 = ax1.contourf(X*1000, Y*1000, T_chip, levels=20, cmap='hot')
ax1.set_xlabel('x (mm)', fontsize=11)
ax1.set_ylabel('y (mm)', fontsize=11)
ax1.set_title('Chip Temperature Distribution', fontsize=12, fontweight='bold')
ax1.set_aspect('equal')
cbar1 = plt.colorbar(im1, ax=ax1)
cbar1.set_label('Temperature (K)', fontsize=10)

# 热源分布
ax2 = axes[1]
im2 = ax2.contourf(X*1000, Y*1000, q_chip/1e8, levels=20, cmap='YlOrRd')
ax1.set_xlabel('x (mm)', fontsize=11)
ax2.set_ylabel('y (mm)', fontsize=11)
ax2.set_title('Heat Source Distribution', fontsize=12, fontweight='bold')
ax2.set_aspect('equal')
cbar2 = plt.colorbar(im2, ax=ax2)
cbar2.set_label('Heat Generation (×10⁸ W/m³)', fontsize=10)

plt.tight_layout()
plt.savefig(f'{output_dir}/case2_chip_cooling.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 案例2完成，图片已保存")

5.3 案例3：热电制冷器分析

问题描述：分析热电制冷器（TEC）的性能。TEC由多个热电偶对组成，通过珀尔帖效应实现制冷。

控制方程：

$$\nabla \cdot (k\nabla T)+\rho J^2-\tau J\cdot \nabla T=0$$

其中，$\tau$为汤姆逊系数。

print("\n" + "="*60)
print("案例3：热电制冷器分析")
print("="*60)

# TEC参数
n_couples = 10  # 热电偶对数
L_leg = 2e-3  # 热电臂长度 (m)
A_leg = 1e-6  # 热电臂截面积 (m²)

# 材料参数（碲化铋基材料）
# P型
alpha_p = 200e-6  # 塞贝克系数 (V/K)
sigma_p = 1e5  # 电导率 (S/m)
k_p = 1.5  # 热导率 (W/m·K)

# N型
alpha_n = -200e-6
sigma_n = 1e5
k_n = 1.5

# 工作条件
I_current = 0.5  # 电流 (A)
T_cold = 280.0  # 冷端温度 (K)
T_hot = 320.0  # 热端温度 (K)

# 一维分析
nz = 50
z = np.linspace(0, L_leg, nz)
dz = z[1] - z[0]

# P型臂温度分布
T_p = np.linspace(T_cold, T_hot, nz)

# 能量平衡分析
# 珀尔帖热
Q_peltier_cold = n_couples * (alpha_p - alpha_n) * I_current * T_cold
Q_peltier_hot = n_couples * (alpha_p - alpha_n) * I_current * T_hot

# 焦耳热
R_total = n_couples * L_leg * (1/sigma_p + 1/sigma_n) / A_leg
Q_joule = I_current**2 * R_total

# 热传导
K_total = n_couples * A_leg * (k_p + k_n) / L_leg
Q_conduction = K_total * (T_hot - T_cold)

# 制冷量
Q_cooling = Q_peltier_cold - 0.5 * Q_joule - Q_conduction

# 散热量
Q_heating = Q_peltier_hot + 0.5 * Q_joule - Q_conduction

# 性能系数
cop = Q_cooling / (I_current * abs(alpha_p - alpha_n) * (T_hot - T_cold) + I_current**2 * R_total)

print(f"\n热电制冷器性能分析:")
print(f"  珀尔帖制冷量: {Q_peltier_cold:.4f} W")
print(f"  焦耳热损失: {Q_joule:.4f} W")
print(f"  热传导损失: {Q_conduction:.4f} W")
print(f"  净制冷量: {Q_cooling:.4f} W")
print(f"  散热量: {Q_heating:.4f} W")
print(f"  性能系数 COP: {cop:.4f}")

# 计算温度分布（考虑焦耳热）
# 简化模型：线性温度分布叠加焦耳热引起的抛物线分布
T_profile = T_cold + (T_hot - T_cold) * z / L_leg
# 焦耳热引起的温升（抛物线分布）
dT_joule = 0.5 * Q_joule * z * (L_leg - z) / (K_total * L_leg**2)
T_profile = T_profile + dT_joule

# 可视化
fig, axes = plt.subplots(1, 2, figsize=(12, 5))

# 温度分布
ax1 = axes[0]
ax1.plot(z*1000, T_profile, 'b-', linewidth=2, label='Temperature')
ax1.axhline(y=T_cold, color='b', linestyle='--', alpha=0.5, label='Cold side')
ax1.axhline(y=T_hot, color='r', linestyle='--', alpha=0.5, label='Hot side')
ax1.set_xlabel('Position (mm)', fontsize=11)
ax1.set_ylabel('Temperature (K)', fontsize=11)
ax1.set_title('TEC Temperature Profile', fontsize=12, fontweight='bold')
ax1.legend()
ax1.grid(True, alpha=0.3)

# 能量流分析
ax2 = axes[1]
categories = ['Peltier\n(Cold)', 'Joule\nHeat', 'Conduction', 'Net\nCooling', 'Peltier\n(Hot)']
values = [Q_peltier_cold*1000, -Q_joule*1000, -Q_conduction*1000, Q_cooling*1000, Q_heating*1000]
colors = ['green', 'red', 'orange', 'blue', 'purple']
bars = ax2.bar(categories, values, color=colors, alpha=0.7)
ax2.set_ylabel('Heat Flow (mW)', fontsize=11)
ax2.set_title('TEC Energy Balance', fontsize=12, fontweight='bold')
ax2.axhline(y=0, color='k', linewidth=0.5)
ax2.grid(True, alpha=0.3, axis='y')

# 添加数值标签
for bar, val in zip(bars, values):
    height = bar.get_height()
    ax2.text(bar.get_x() + bar.get_width()/2., height,
             f'{val:.2f}',
             ha='center', va='bottom' if height > 0 else 'top', fontsize=9)

plt.tight_layout()
plt.savefig(f'{output_dir}/case3_thermoelectric_cooler.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 案例3完成，图片已保存")

5.4 案例4：瞬态电-热分析

问题描述：分析电阻加热器的瞬态响应，研究通电后温度随时间的变化过程。

控制方程：

$$\rho c_p\frac{\partial T}{\partial t}=k\frac{{\partial }^{2}T}{\partial x^2}+q^{\prime \prime \prime }$$

print("\n" + "="*60)
print("案例4：瞬态电-热分析")
print("="*60)

# 参数
L = 0.05  # m
nx = 50
x = np.linspace(0, L, nx)
dx = x[1] - x[0]

# 材料参数（铜）
rho_0 = 1.68e-8  # Ohm·m
alpha = 0.0039  # 1/K
k = 400.0  # W/m·K
cp = 385.0  # J/kg·K
rho_density = 8960.0  # kg/m³

A = 1e-6  # m²
V = 0.1  # V

# 时间参数
t_max = 10.0  # s
dt = 0.01
nt = int(t_max / dt)

# 初始条件
T = np.ones(nx) * 300.0
T_ambient = 300.0
h_conv = 50.0  # W/m²·K
P = 4e-3  # 周长 (m)

# 存储结果
T_history = []
t_history = []
I_history = []

print("开始瞬态计算...")
for n in range(nt):
    t = n * dt
    
    # 计算当前电阻率
    rho = rho_0 * (1 + alpha * (T - T_ambient))
    sigma = 1.0 / rho
    
    # 计算电阻
    R = np.sum(rho * dx / A)
    I = V / R
    
    # 计算焦耳热
    q_joule = sigma * (V/L)**2
    
    # 时间推进（Crank-Nicolson）
    T_new = T.copy()
    
    # 内部节点
    for i in range(1, nx-1):
        d2T = (T[i+1] - 2*T[i] + T[i-1]) / dx**2
        # 显式欧拉（简化）
        T_new[i] = T[i] + dt * (k * d2T + q_joule[i]) / (rho_density * cp)
        # 对流冷却
        T_new[i] -= dt * h_conv * P * (T[i] - T_ambient) / (rho_density * cp * A)
    
    # 边界条件
    T_new[0] = T_new[1]  # 绝热
    T_new[-1] = T_new[-2]
    
    T = T_new
    
    # 记录结果
    if n % 100 == 0:
        T_history.append(T.copy())
        t_history.append(t)
        I_history.append(I)
        if n % 500 == 0:
            print(f"  t = {t:.2f} s: T_max = {np.max(T):.2f} K, I = {I:.4f} A")

print(f"\n瞬态分析结果:")
print(f"  最终最高温度: {np.max(T):.2f} K")
print(f"  最终电流: {I:.6f} A")
print(f"  稳态电阻: {V/I:.6f} Ohm")

# 可视化
fig, axes = plt.subplots(1, 2, figsize=(12, 5))

# 温度随时间变化
ax1 = axes[0]
for i, (T_prof, t_val) in enumerate(zip(T_history[::5], t_history[::5])):
    ax1.plot(x*1000, T_prof, label=f't={t_val:.1f}s' if i % 2 == 0 else '')
ax1.set_xlabel('Position (mm)', fontsize=11)
ax1.set_ylabel('Temperature (K)', fontsize=11)
ax1.set_title('Transient Temperature Profiles', fontsize=12, fontweight='bold')
ax1.legend(fontsize=8)
ax1.grid(True, alpha=0.3)

# 电流随时间变化
ax2 = axes[1]
ax2.plot(t_history, I_history, 'b-', linewidth=2)
ax2.set_xlabel('Time (s)', fontsize=11)
ax2.set_ylabel('Current (A)', fontsize=11)
ax2.set_title('Current vs Time', fontsize=12, fontweight='bold')
ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig(f'{output_dir}/case4_transient_analysis.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 案例4完成，图片已保存")

5.5 案例5：非线性电阻加热

问题描述：分析具有强温度依赖性电阻率的加热元件，展示非线性效应。

print("\n" + "="*60)
print("案例5：非线性电阻加热分析")
print("="*60)

# PTC材料参数
T_c = 350.0  # 居里温度 (K)
rho_0 = 1e-4  # 参考电阻率
A_ptc = 0.1  # PTC系数

# 几何参数
L = 0.02  # m
nx = 100
x = np.linspace(0, L, nx)
dx = x[1] - x[0]

# 材料参数
k = 2.0  # W/m·K

# 电压扫描
V_range = np.linspace(1, 50, 20)
T_max_list = []
I_list = []
P_list = []

print("进行电压扫描分析...")
for V_applied in V_range:
    T = np.ones(nx) * 300.0
    
    # 迭代求解
    for iteration in range(100):
        T_old = T.copy()
        
        # PTC电阻率
        rho = rho_0 * np.ones(nx)
        mask = T > T_c
        rho[mask] = rho_0 * np.exp(A_ptc * (T[mask] - T_c))
        sigma = 1.0 / rho
        
        # 求解电场
        R_total = np.sum(rho * dx / (np.pi * (0.5e-3)**2))
        I = V_applied / R_total
        
        # 焦耳热
        q_joule = I**2 * rho / (np.pi * (0.5e-3)**2)**2
        
        # 求解温度
        for i in range(1, nx-1):
            d2T = (T[i+1] - 2*T[i] + T[i-1]) / dx**2
            T[i] = (-q_joule[i] * dx**2 / k + T[i+1] + T[i-1]) / 2
        
        T[0] = 300.0  # 固定温度
        T[-1] = 300.0
        
        if np.max(np.abs(T - T_old)) < 1e-4:
            break
    
    T_max_list.append(np.max(T))
    I_list.append(I)
    P_list.append(V_applied * I)

T_max_list = np.array(T_max_list)
I_list = np.array(I_list)
P_list = np.array(P_list)

print(f"\n非线性分析结果:")
print(f"  最大温度: {np.max(T_max_list):.2f} K")
print(f"  最大电流: {np.max(I_list):.4f} A")
print(f"  最大功率: {np.max(P_list):.4f} W")

# 可视化
fig, axes = plt.subplots(1, 3, figsize=(15, 4))

# I-V特性
ax1 = axes[0]
ax1.plot(V_range, I_list, 'b-o', markersize=4, linewidth=2)
ax1.set_xlabel('Voltage (V)', fontsize=11)
ax1.set_ylabel('Current (A)', fontsize=11)
ax1.set_title('I-V Characteristic (PTC Effect)', fontsize=12, fontweight='bold')
ax1.grid(True, alpha=0.3)

# 温度-电压
ax2 = axes[1]
ax2.plot(V_range, T_max_list, 'r-s', markersize=4, linewidth=2)
ax2.axhline(y=T_c, color='k', linestyle='--', label='Curie Temp')
ax2.set_xlabel('Voltage (V)', fontsize=11)
ax2.set_ylabel('Max Temperature (K)', fontsize=11)
ax2.set_title('Temperature vs Voltage', fontsize=12, fontweight='bold')
ax2.legend()
ax2.grid(True, alpha=0.3)

# 功率-电压
ax3 = axes[2]
ax3.plot(V_range, P_list, 'g-^', markersize=4, linewidth=2)
ax3.set_xlabel('Voltage (V)', fontsize=11)
ax3.set_ylabel('Power (W)', fontsize=11)
ax3.set_title('Power vs Voltage', fontsize=12, fontweight='bold')
ax3.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig(f'{output_dir}/case5_nonlinear_heating.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 案例5完成，图片已保存")

5.6 案例6：电-热耦合优化设计

问题描述：优化电阻加热器的设计参数，实现均匀的温度分布。

print("\n" + "="*60)
print("案例6：电-热耦合优化设计")
print("="*60)

# 优化目标：实现均匀的温度分布
L = 0.1  # m
nx = 100
x = np.linspace(0, L, nx)
dx = x[1] - x[0]

# 基础参数
k = 50.0  # W/m·K
h = 20.0  # W/m²·K
T_ambient = 300.0

# 不同设计方案比较
designs = {
    'Uniform': lambda x: np.ones_like(x),
    'Linear': lambda x: 1 + 0.5 * x / L,
    'Parabolic': lambda x: 1 + 2 * (x/L - 0.5)**2,
    'Exponential': lambda x: np.exp(2 * (x/L - 0.5))
}

results = {}

for name, sigma_func in designs.items():
    print(f"\n分析设计方案: {name}")
    
    # 归一化电导率分布
    sigma_base = 1e6
    sigma_dist = sigma_func(x) * sigma_base
    
    # 施加电压
    V = 10.0
    
    # 迭代求解
    T = np.ones(nx) * T_ambient
    for iteration in range(100):
        T_old = T.copy()
        
        # 计算电流分布
        R_local = 1.0 / (sigma_dist * np.pi * (1e-3)**2)
        # 简化：假设电流均匀
        I_total = V / np.mean(R_local) * nx
        
        # 焦耳热
        q_joule = I_total**2 * R_local / (np.pi * (1e-3)**2) / nx**2
        
        # 求解温度
        for i in range(1, nx-1):
            d2T = (T[i+1] - 2*T[i] + T[i-1]) / dx**2
            # 稳态平衡
            T[i] = (-q_joule[i] / k + (T[i+1] + T[i-1]) / dx**2) * dx**2 / 2
            # 对流损失
            T[i] -= h * 2 * np.pi * 1e-3 * (T[i] - T_ambient) * dx**2 / (k * np.pi * (1e-3)**2)
        
        T[0] = T[1]
        T[-1] = T[-2]
        
        if np.max(np.abs(T - T_old)) < 1e-5:
            break
    
    # 评估均匀性
    T_mean = np.mean(T)
    T_std = np.std(T)
    uniformity = 1 - T_std / T_mean
    
    results[name] = {
        'T': T.copy(),
        'T_mean': T_mean,
        'T_max': np.max(T),
        'T_min': np.min(T),
        'uniformity': uniformity,
        'sigma': sigma_dist
    }
    
    print(f"  平均温度: {T_mean:.2f} K")
    print(f"  温度范围: {np.min(T):.2f} - {np.max(T):.2f} K")
    print(f"  均匀性指标: {uniformity:.4f}")

# 可视化
fig, axes = plt.subplots(2, 2, figsize=(12, 10))

# 温度分布比较
ax1 = axes[0, 0]
colors = ['blue', 'red', 'green', 'purple']
for (name, result), color in zip(results.items(), colors):
    ax1.plot(x*1000, result['T'], color=color, linewidth=2, label=f"{name} (U={result['uniformity']:.3f})")
ax1.set_xlabel('Position (mm)', fontsize=11)
ax1.set_ylabel('Temperature (K)', fontsize=11)
ax1.set_title('Temperature Distribution Comparison', fontsize=12, fontweight='bold')
ax1.legend()
ax1.grid(True, alpha=0.3)

# 电导率分布
ax2 = axes[0, 1]
for (name, result), color in zip(results.items(), colors):
    ax2.plot(x*1000, result['sigma']/1e6, color=color, linewidth=2, label=name)
ax2.set_xlabel('Position (mm)', fontsize=11)
ax2.set_ylabel('Conductivity (MS/m)', fontsize=11)
ax2.set_title('Conductivity Distribution', fontsize=12, fontweight='bold')
ax2.legend()
ax2.grid(True, alpha=0.3)

# 均匀性比较
ax3 = axes[1, 0]
names = list(results.keys())
uniformities = [results[n]['uniformity'] for n in names]
bars = ax3.bar(names, uniformities, color=colors, alpha=0.7)
ax3.set_ylabel('Uniformity Index', fontsize=11)
ax3.set_title('Temperature Uniformity Comparison', fontsize=12, fontweight='bold')
ax3.set_ylim([0, 1])
for bar, val in zip(bars, uniformities):
    ax3.text(bar.get_x() + bar.get_width()/2., bar.get_height() + 0.01,
             f'{val:.3f}', ha='center', fontsize=10)

# 温度统计
ax4 = axes[1, 1]
x_pos = np.arange(len(names))
T_means = [results[n]['T_mean'] - T_ambient for n in names]
T_ranges = [(results[n]['T_max'] - results[n]['T_min'])/2 for n in names]
ax4.errorbar(x_pos, T_means, yerr=T_ranges, fmt='o', capsize=5, capthick=2, markersize=8)
ax4.set_xticks(x_pos)
ax4.set_xticklabels(names)
ax4.set_ylabel('Temperature Rise (K)', fontsize=11)
ax4.set_title('Temperature Rise with Variation', fontsize=12, fontweight='bold')
ax4.grid(True, alpha=0.3, axis='y')

plt.tight_layout()
plt.savefig(f'{output_dir}/case6_optimization_design.png', dpi=150, bbox_inches='tight')
plt.close()
print("\n✓ 案例6完成，图片已保存")

print("\n" + "="*60)
print("所有案例仿真完成！")
print("="*60)
print("\n生成的图片文件:")
print("  - case1_resistor_heating.png")
print("  - case2_chip_cooling.png")
print("  - case3_thermoelectric_cooler.png")
print("  - case4_transient_analysis.png")
print("  - case5_nonlinear_heating.png")
print("  - case6_optimization_design.png")
print("="*60)