主题039：磁纳米颗粒仿真

引言

磁纳米颗粒是尺寸在1-100纳米范围内的磁性材料，具有独特的超顺磁性、高比表面积和可调控的磁学性质。这些特性使其在生物医学（磁热疗、靶向药物递送、MRI造影）、环境治理（污染物吸附）、数据存储等领域有广泛应用。本教程将系统介绍磁纳米颗粒的物理模型、弛豫机制和生物医学应用仿真。

---

1. 磁纳米颗粒基础

1.1 尺寸效应

当磁性颗粒尺寸减小到纳米尺度时，会表现出与块体材料截然不同的性质：

超顺磁性：颗粒尺寸小于临界值时，热扰动能克服磁晶各向异性能，磁矩可自由翻转，表现为超顺磁性。

高比表面积：纳米颗粒具有极高的比表面积，表面原子占比显著增加。

量子尺寸效应：电子能级离散化，影响磁学性质。

1.2 典型材料

材料	饱和磁化强度	各向异性常数	应用
Fe₃O₄	480 kA/m	10⁴ J/m³	MRI、磁热疗
γ-Fe₂O₃	380 kA/m	5×10³ J/m³	磁记录
CoFe₂O₄	400 kA/m	10⁵ J/m³	高频应用
Fe	1700 kA/m	5×10⁴ J/m³	高磁矩应用

---

2. 超顺磁性模型

2.1 Langevin函数

超顺磁性颗粒的磁化强度由Langevin函数描述：

$$M=M_{s}L(\alpha )$$

$$L(\alpha )=\coth (\alpha )-\frac{1}{\alpha }$$

其中 $\alpha =\frac{{\mu }_{0}mH}{k_{B}T}$ 为无量纲参数，$m=M_{s}V$ 为颗粒磁矩。

2.2 实现代码

class Superparamagnetism:
    def __init__(self, params):
        self.params = params
        self.V = (4/3) * np.pi * (params.diameter * 1e-9 / 2)**3
        self.m = params.Ms * self.V
        self.KV = params.K * self.V
    
    def langevin(self, alpha):
        # 避免除零
        alpha = np.where(np.abs(alpha) < 1e-10, 1e-10, alpha)
        L = np.cosh(alpha) / np.sinh(alpha) - 1 / alpha
        return L
    
    def magnetization(self, H):
        alpha = self.mu0 * self.m * H / (self.kB * self.params.temperature)
        M = self.params.Ms * self.langevin(alpha)
        return M

2.3 阻塞温度

阻塞温度是超顺磁性与铁磁性转变的临界温度：

$$T_B=\frac{KV}{25k_B}$$

当 $T>T_B$ 时，颗粒表现为超顺磁性；当 $T<T_B$ 时，磁矩被"阻塞"在易磁化方向。

仿真结果：

• 20nm Fe₃O₄颗粒：$T_B$ = 12.1 K
• 阻塞温度与体积成正比
• 室温超顺磁要求 $d<$ 25 nm

---

3. 磁弛豫机制

3.1 布朗弛豫

颗粒在液体中旋转导致的磁矩弛豫：

$${\tau }_B=\frac{3V\eta }{k_{B}T}$$

其中 $\eta$ 为液体粘度。

3.2 尼尔弛豫

颗粒内部磁矩翻转导致的弛豫：

$${\tau }_N={\tau }_0\exp \left(\frac{KV}{k_{B}T}\right)$$

其中 ${\tau }_0\approx 10^{-9}$ s 为尝试时间。

3.3 有效弛豫时间

两种机制并联：

$$\frac{1}{{\tau }_{eff}}=\frac{1}{{\tau }_B}+\frac{1}{{\tau }_N}$$

主导机制判断：

• 小颗粒/低粘度 → 尼尔弛豫主导
• 大颗粒/高粘度 → 布朗弛豫主导

3.4 交流磁化率

Debye模型描述频散特性：

$${\chi }^{\prime }(\omega )=\frac{{\chi }_0}{1+(\omega {\tau }_{eff}{)}^2}$$

$${\chi }^{\prime \prime }(\omega )=\frac{{\chi }_0\omega {\tau }_{eff}}{1+(\omega {\tau }_{eff}{)}^2}$$

---

4. 磁热疗SAR计算

4.1 比吸收率(SAR)

SAR表示单位质量颗粒吸收的功率：

$$SAR=\frac{P}{m}[W/kg]$$

4.2 线性响应理论

$$SAR=\frac{{\mu }_0\pi {\chi }^{\prime \prime }{H}^{2}f}{\rho }$$

仿真结果：

• 20nm Fe₃O₄ @ 20kA/m, 300kHz
• SAR (线性响应)：0.21 W/kg
• SAR (磁滞损耗)：80535 W/kg

4.3 尺寸优化

最优颗粒尺寸与磁场参数相关：

• 最优尺寸：17.7 nm
• 最大SAR：2.94×10⁶ W/kg

4.4 临床应用考虑

安全限制：

• 磁场强度：H × f < 5×10⁹ A·m⁻¹·s⁻¹
• 典型参数：H = 10-30 kA/m, f = 100-500 kHz

治疗温度：

• 目标温度：42-45°C
• 过热损伤：>50°C

---

5. 颗粒动力学模拟

5.1 受力分析

磁力：$\mathbf{F}=(\mathbf{m}\cdot \nabla )\mathbf{B}$

粘性阻力：${\mathbf{F}}_d=-\gamma \mathbf{v}$，其中 $\gamma =6\pi \eta r$

布朗力：随机热涨落导致的扩散

5.2 运动方程

$$\gamma \frac{d\mathbf{r}}{dt}={\mathbf{F}}_{mag}+{\mathbf{F}}_{Brown}$$

5.3 均方位移(MSD)

$$MSD=4Dt$$

扩散系数：$D=\frac{k_{B}T}{6\pi \eta r}$

---

6. 仿真结果分析

6.1 超顺磁性

图1：磁化曲线

• 小颗粒更易饱和
• Langevin函数描述非线性响应
• 磁化率随尺寸增大

6.2 弛豫特性

图2：弛豫分析

• 布朗弛豫时间：~10⁻⁴ s
• 尼尔弛豫时间：~10⁻³ s
• 有效弛豫时间：~10⁻⁴ s
• 主导机制：混合弛豫

6.3 磁热疗

图3：SAR特性

• SAR与H²成正比
• 存在最优频率
• 尺寸优化可提高效率

6.4 动力学

图4：颗粒运动

• 初始随机分布
• 磁场梯度导致聚集
• MSD呈线性增长

---

7. 常见报错与解决方案

7.1 数值溢出

现象：尼尔弛豫时间计算溢出

原因：指数项过大

解决方案：

• 检查参数单位
• 使用对数计算
• 限制温度范围

7.2 收敛问题

现象：动力学模拟发散

解决方案：

• 减小时间步长
• 检查力计算
• 添加数值阻尼

7.3 SAR异常

现象：SAR值过大或负值

解决方案：

• 检查磁场单位
• 验证磁化率计算
• 确认频率范围

---

8. 总结

本教程介绍了磁纳米颗粒的核心仿真方法：

1. 超顺磁性：Langevin模型描述磁化行为
2. 弛豫机制：布朗与尼尔弛豫的竞争
3. 磁热疗：SAR计算与尺寸优化
4. 动力学：磁力与布朗运动的耦合

关键要点

• 尺寸决定磁学性质
• 弛豫机制影响应用选择
• 磁热疗需优化尺寸和场参数
• 动力学模拟指导靶向递送

应用建议

• MRI造影：超顺磁Fe₃O₄，d=5-10nm
• 磁热疗：d=15-20nm，优化SAR
• 靶向递送：表面功能化，磁场引导
• 分离纯化：高磁矩，磁场梯度

---

附录

完整代码见 run_simulation.py，包含：

• 超顺磁性模型
• 弛豫分析
• 磁热疗SAR计算
• 颗粒动力学模拟

运行：

python run_simulation.py

生成文件：

• fig1_superparamagnetism.png：超顺磁性分析
• fig2_relaxation.png：弛豫分析
• fig3_hyperthermia.png：磁热疗SAR
• fig4_dynamics.png：颗粒动力学
• nanoparticle_results.json：结果数据

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
主题039：磁纳米颗粒仿真
=======================
本代码实现磁纳米颗粒的关键仿真功能：
1. 超顺磁性模型 (Langevin函数)
2. 布朗弛豫与尼尔弛豫
3. 磁热疗SAR计算
4. 颗粒动力学模拟

作者：AI仿真专家
日期：2026-03-23
"""

import numpy as np
import matplotlib
matplotlib.use('Agg')  # 无头模式，不显示GUI
import matplotlib.pyplot as plt
from matplotlib.patches import Circle, FancyBboxPatch
from matplotlib.collections import PatchCollection
import json
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


class RelaxationMechanism(Enum):
    """弛豫机制"""
    BROWN = "布朗弛豫"
    NEEL = "尼尔弛豫"
    BOTH = "混合弛豫"


@dataclass
class ParticleParams:
    """纳米颗粒参数"""
    name: str
    diameter: float       # 颗粒直径 [nm]
    Ms: float            # 饱和磁化强度 [A/m]
    K: float             # 各向异性常数 [J/m^3]
    viscosity: float     # 液体粘度 [Pa·s]
    temperature: float   # 温度 [K]
    density: float = 5000  # 密度 [kg/m^3]


class Superparamagnetism:
    """超顺磁性模型"""
    
    def __init__(self, params: ParticleParams):
        self.params = params
        
        # 物理常数
        self.kB = 1.38e-23  # 玻尔兹曼常数
        self.mu0 = 4 * np.pi * 1e-7
        
        # 计算颗粒体积
        self.V = (4/3) * np.pi * (params.diameter * 1e-9 / 2)**3
        
        # 磁矩
        self.m = params.Ms * self.V
        
        # 各向异性能垒
        self.KV = params.K * self.V
    
    def langevin(self, alpha: np.ndarray) -> np.ndarray:
        """
        Langevin函数
        
        L(alpha) = coth(alpha) - 1/alpha
        
        参数:
            alpha: 无量纲参数 m*H/(kB*T)
        返回:
            L: Langevin函数值
        """
        # 避免除零
        alpha = np.where(np.abs(alpha) < 1e-10, 1e-10, alpha)
        
        L = np.cosh(alpha) / np.sinh(alpha) - 1 / alpha
        
        # 小值近似
        small_alpha = np.abs(alpha) < 0.1
        L[small_alpha] = alpha[small_alpha] / 3 - alpha[small_alpha]**3 / 45
        
        return L
    
    def magnetization(self, H: np.ndarray) -> np.ndarray:
        """
        计算磁化强度 (超顺磁性)
        
        M = Ms * L(alpha)
        alpha = mu0 * m * H / (kB * T)
        
        参数:
            H: 磁场强度 [A/m]
        返回:
            M: 磁化强度 [A/m]
        """
        alpha = self.mu0 * self.m * H / (self.kB * self.params.temperature)
        M = self.params.Ms * self.langevin(alpha)
        return M
    
    def susceptibility(self, H: float) -> float:
        """
        计算磁化率
        
        参数:
            H: 磁场强度 [A/m]
        返回:
            chi: 磁化率
        """
        alpha = self.mu0 * self.m * H / (self.kB * self.params.temperature)
        
        if abs(alpha) < 0.01:
            # 低场极限
            chi = self.mu0 * self.m**2 / (3 * self.kB * self.params.temperature * self.V)
        else:
            # 数值微分
            dH = H * 0.01
            dM = self.magnetization(np.array([H + dH]))[0] - self.magnetization(np.array([H]))[0]
            chi = dM / dH
        
        return chi
    
    def blocking_temperature(self) -> float:
        """
        计算阻塞温度
        
        TB = KV / (25 * kB)
        
        返回:
            TB: 阻塞温度 [K]
        """
        TB = self.KV / (25 * self.kB)
        return TB


class MagneticRelaxation:
    """磁弛豫分析"""
    
    def __init__(self, params: ParticleParams):
        self.params = params
        self.sp = Superparamagnetism(params)
        
        # 物理常数
        self.kB = 1.38e-23
        self.mu0 = 4 * np.pi * 1e-7
    
    def brown_relaxation_time(self) -> float:
        """
        布朗弛豫时间
        
        tau_B = 3 * V * eta / (kB * T)
        
        返回:
            tau_B: 布朗弛豫时间 [s]
        """
        tau_B = 3 * self.sp.V * self.params.viscosity / (self.kB * self.params.temperature)
        return tau_B
    
    def neel_relaxation_time(self) -> float:
        """
        尼尔弛豫时间
        
        tau_N = tau_0 * exp(KV / (kB * T))
        
        返回:
            tau_N: 尼尔弛豫时间 [s]
        """
        tau_0 = 1e-9  # 尝试时间 [s]
        tau_N = tau_0 * np.exp(self.sp.KV / (self.kB * self.params.temperature))
        return tau_N
    
    def effective_relaxation_time(self) -> float:
        """
        有效弛豫时间
        
        1/tau_eff = 1/tau_B + 1/tau_N
        
        返回:
            tau_eff: 有效弛豫时间 [s]
        """
        tau_B = self.brown_relaxation_time()
        tau_N = self.neel_relaxation_time()
        
        tau_eff = 1 / (1/tau_B + 1/tau_N)
        return tau_eff
    
    def dominant_mechanism(self) -> RelaxationMechanism:
        """
        判断主导弛豫机制
        
        返回:
            主导弛豫机制
        """
        tau_B = self.brown_relaxation_time()
        tau_N = self.neel_relaxation_time()
        
        if tau_B < tau_N * 0.1:
            return RelaxationMechanism.BROWN
        elif tau_N < tau_B * 0.1:
            return RelaxationMechanism.NEEL
        else:
            return RelaxationMechanism.BOTH
    
    def ac_susceptibility(self, frequencies: np.ndarray) -> Tuple[np.ndarray, np.ndarray]:
        """
        计算交流磁化率
        
        参数:
            frequencies: 频率数组 [Hz]
        返回:
            chi_real: 实部
            chi_imag: 虚部
        """
        tau_eff = self.effective_relaxation_time()
        
        # 静态磁化率
        chi_0 = self.sp.susceptibility(1e3)  # 1 kA/m
        
        omega = 2 * np.pi * frequencies
        
        # Debye模型
        chi_real = chi_0 / (1 + (omega * tau_eff)**2)
        chi_imag = chi_0 * omega * tau_eff / (1 + (omega * tau_eff)**2)
        
        return chi_real, chi_imag


class MagneticHyperthermia:
    """磁热疗SAR计算"""
    
    def __init__(self, params: ParticleParams, concentration: float = 1e15):
        """
        参数:
            params: 颗粒参数
            concentration: 颗粒浓度 [particles/m^3]
        """
        self.params = params
        self.sp = Superparamagnetism(params)
        self.relaxation = MagneticRelaxation(params)
        self.concentration = concentration
        
        self.mu0 = 4 * np.pi * 1e-7
    
    def sar_linear_response(self, H: float, f: float) -> float:
        """
        线性响应理论计算SAR
        
        SAR = (mu0 * pi * chi'' * H^2 * f) / rho
        
        参数:
            H: 磁场幅值 [A/m]
            f: 频率 [Hz]
        返回:
            SAR: 比吸收率 [W/kg]
        """
        chi_real, chi_imag = self.relaxation.ac_susceptibility(np.array([f]))
        
        # 单位质量SAR
        mass_per_particle = self.sp.V * self.params.density
        
        # 每个颗粒的功率
        P_per_particle = np.pi * self.mu0 * chi_imag[0] * H**2 * f * self.sp.V
        
        # SAR
        SAR = P_per_particle / mass_per_particle
        
        return SAR
    
    def sar_hysteresis(self, H_max: float, f: float, n_points: int = 100) -> float:
        """
        磁滞损耗计算SAR
        
        参数:
            H_max: 最大磁场 [A/m]
            f: 频率 [Hz]
            n_points: 采样点数
        返回:
            SAR: 比吸收率 [W/kg]
        """
        # 一个周期的磁滞回线
        H_cycle = H_max * np.sin(np.linspace(0, 2*np.pi, n_points))
        M_cycle = self.sp.magnetization(H_cycle)
        
        # 计算面积 (磁滞损耗)
        dH = np.diff(H_cycle)
        area = np.sum(M_cycle[:-1] * dH)
        
        # 每周期能量损耗 [J/m^3]
        W = self.mu0 * area
        
        # SAR [W/kg]
        mass_density = self.params.density
        SAR = W * f / mass_density
        
        return abs(SAR)
    
    def optimal_size(self, H: float, f: float, 
                    d_min: float = 5, d_max: float = 50) -> Dict:
        """
        寻找最优颗粒尺寸
        
        参数:
            H: 磁场 [A/m]
            f: 频率 [Hz]
            d_min, d_max: 尺寸范围 [nm]
        返回:
            最优尺寸分析
        """
        diameters = np.linspace(d_min, d_max, 100)
        sars = []
        
        for d in diameters:
            params_d = ParticleParams(
                name=self.params.name,
                diameter=d,
                Ms=self.params.Ms,
                K=self.params.K,
                viscosity=self.params.viscosity,
                temperature=self.params.temperature,
                density=self.params.density
            )
            
            mh = MagneticHyperthermia(params_d, self.concentration)
            
            try:
                sar = mh.sar_linear_response(H, f)
                sars.append(sar)
            except:
                sars.append(0)
        
        sars = np.array(sars)
        idx_opt = np.argmax(sars)
        
        return {
            'diameters_nm': diameters,
            'sar_W_per_kg': sars,
            'optimal_diameter_nm': diameters[idx_opt],
            'max_sar_W_per_kg': sars[idx_opt]
        }


class ParticleDynamics:
    """颗粒动力学模拟"""
    
    def __init__(self, 
                 n_particles: int = 100,
                 box_size: float = 1e-6,  # 1 um
                 params: ParticleParams = None):
        self.n = n_particles
        self.L = box_size
        self.params = params if params else ParticleParams(
            name="Fe3O4", diameter=20, Ms=480e3, 
            K=1e4, viscosity=1e-3, temperature=300
        )
        
        # 初始化位置
        self.positions = np.random.rand(n_particles, 2) * box_size
        self.velocities = np.zeros((n_particles, 2))
        
        # 磁矩方向
        self.moments = np.random.randn(n_particles, 2)
        self.moments /= np.linalg.norm(self.moments, axis=1)[:, np.newaxis]
        
        self.kB = 1.38e-23
        self.mu0 = 4 * np.pi * 1e-7
    
    def magnetic_force(self, B_field: callable) -> np.ndarray:
        """
        计算磁力
        
        F = (m · ∇)B
        
        参数:
            B_field: 磁场函数 B(x,y) -> (Bx, By)
        返回:
            F: 力 [N]
        """
        F = np.zeros_like(self.positions)
        
        dx = self.L * 0.001
        
        for i in range(self.n):
            x, y = self.positions[i]
            
            # 数值计算梯度
            Bx_plus, By_plus = B_field(x + dx, y)
            Bx_minus, By_minus = B_field(x - dx, y)
            dBx_dx = (Bx_plus - Bx_minus) / (2 * dx)
            dBy_dx = (By_plus - By_minus) / (2 * dx)
            
            Bx_plus, By_plus = B_field(x, y + dx)
            Bx_minus, By_minus = B_field(x, y - dx)
            dBx_dy = (Bx_plus - Bx_minus) / (2 * dx)
            dBy_dy = (By_plus - By_minus) / (2 * dx)
            
            mx, my = self.moments[i]
            
            # F = (m · ∇)B
            F[i, 0] = mx * dBx_dx + my * dBx_dy
            F[i, 1] = mx * dBy_dx + my * dBy_dy
        
        # 磁矩大小
        m = self.params.Ms * (4/3) * np.pi * (self.params.diameter * 1e-9 / 2)**3
        F *= m * self.mu0
        
        return F
    
    def brownian_motion(self, dt: float) -> np.ndarray:
        """
        布朗运动
        
        参数:
            dt: 时间步长 [s]
        返回:
            displacement: 位移 [m]
        """
        # 扩散系数
        kB = 1.38e-23
        T = self.params.temperature
        eta = self.params.viscosity
        r = self.params.diameter * 1e-9 / 2
        
        D = kB * T / (6 * np.pi * eta * r)
        
        # 随机位移
        sigma = np.sqrt(2 * D * dt)
        displacement = np.random.randn(self.n, 2) * sigma
        
        return displacement
    
    def step(self, B_field: callable, dt: float):
        """
        演化一步
        
        参数:
            B_field: 磁场函数
            dt: 时间步长 [s]
        """
        # 磁力
        F = self.magnetic_force(B_field)
        
        # 粘性阻力 (Stokes定律)
        eta = self.params.viscosity
        r = self.params.diameter * 1e-9 / 2
        gamma = 6 * np.pi * eta * r
        
        # 速度更新
        self.velocities = F / gamma
        
        # 位置更新
        self.positions += self.velocities * dt
        
        # 布朗运动
        self.positions += self.brownian_motion(dt)
        
        # 边界条件
        self.positions = np.mod(self.positions, self.L)
    
    def simulate(self, B_field: callable, 
                dt: float = 1e-6, 
                n_steps: int = 1000) -> Dict:
        """
        模拟颗粒运动
        
        参数:
            B_field: 磁场函数
            dt: 时间步长
            n_steps: 步数
        返回:
            轨迹数据
        """
        trajectory = np.zeros((n_steps, self.n, 2))
        
        for i in range(n_steps):
            self.step(B_field, dt)
            trajectory[i] = self.positions.copy()
        
        return {
            'trajectory': trajectory,
            'dt': dt,
            'n_steps': n_steps,
            'box_size': self.L
        }


def create_visualizations():
    """创建可视化图表"""
    
    print("\n" + "="*70)
    print("生成可视化")
    print("="*70)
    
    # 图1: 超顺磁性
    fig1, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    # 不同尺寸的颗粒
    diameters = [10, 20, 30]  # nm
    colors = ['blue', 'green', 'red']
    
    ax = axes[0, 0]
    for d, color in zip(diameters, colors):
        params = ParticleParams(
            name="Fe3O4", diameter=d, Ms=480e3,
            K=1e4, viscosity=1e-3, temperature=300
        )
        sp = Superparamagnetism(params)
        
        H = np.linspace(0, 500e3, 100)  # A/m
        M = sp.magnetization(H)
        
        ax.plot(H / 1e3, M / 1e3, color=color, linewidth=2, 
                label=f'd = {d} nm')
    
    ax.set_xlabel('磁场 H [kA/m]', fontsize=12)
    ax.set_ylabel('磁化强度 M [kA/m]', fontsize=12)
    ax.set_title('超顺磁性磁化曲线', fontsize=14)
    ax.legend()
    ax.grid(True, alpha=0.3)
    
    # Langevin函数
    ax = axes[0, 1]
    alpha = np.linspace(-10, 10, 100)
    L = sp.langevin(alpha)
    ax.plot(alpha, L, 'b-', linewidth=2)
    ax.axhline(y=1, color='r', linestyle='--', alpha=0.5)
    ax.axhline(y=-1, color='r', linestyle='--', alpha=0.5)
    ax.set_xlabel('α = μ₀mH/(kBT)', fontsize=12)
    ax.set_ylabel('L(α)', fontsize=12)
    ax.set_title('Langevin函数', fontsize=14)
    ax.grid(True, alpha=0.3)
    
    # 阻塞温度
    ax = axes[1, 0]
    d_range = np.linspace(5, 50, 100)
    TBs = []
    for d in d_range:
        params = ParticleParams(
            name="Fe3O4", diameter=d, Ms=480e3,
            K=1e4, viscosity=1e-3, temperature=300
        )
        sp = Superparamagnetism(params)
        TBs.append(sp.blocking_temperature())
    
    ax.semilogy(d_range, TBs, 'g-', linewidth=2)
    ax.axhline(y=300, color='r', linestyle='--', label='室温 (300K)')
    ax.set_xlabel('颗粒直径 [nm]', fontsize=12)
    ax.set_ylabel('阻塞温度 TB [K]', fontsize=12)
    ax.set_title('阻塞温度与颗粒尺寸', fontsize=14)
    ax.legend()
    ax.grid(True, alpha=0.3)
    
    # 磁化率
    ax = axes[1, 1]
    chi_values = []
    for d in d_range:
        params = ParticleParams(
            name="Fe3O4", diameter=d, Ms=480e3,
            K=1e4, viscosity=1e-3, temperature=300
        )
        sp = Superparamagnetism(params)
        chi_values.append(sp.susceptibility(1e3))
    
    ax.semilogy(d_range, chi_values, 'm-', linewidth=2)
    ax.set_xlabel('颗粒直径 [nm]', fontsize=12)
    ax.set_ylabel('磁化率 χ', fontsize=12)
    ax.set_title('初始磁化率与颗粒尺寸', fontsize=14)
    ax.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('fig1_superparamagnetism.png', dpi=150, bbox_inches='tight')
    print("  fig1_superparamagnetism.png 已保存")
    plt.close()
    
    # 图2: 弛豫分析
    fig2, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    params = ParticleParams(
        name="Fe3O4", diameter=20, Ms=480e3,
        K=1e4, viscosity=1e-3, temperature=300
    )
    relax = MagneticRelaxation(params)
    
    # 弛豫时间vs尺寸
    ax = axes[0, 0]
    tau_Bs = []
    tau_Ns = []
    for d in d_range:
        params_d = ParticleParams(
            name="Fe3O4", diameter=d, Ms=480e3,
            K=1e4, viscosity=1e-3, temperature=300
        )
        r = MagneticRelaxation(params_d)
        tau_Bs.append(r.brown_relaxation_time())
        tau_Ns.append(r.neel_relaxation_time())
    
    ax.loglog(d_range, tau_Bs, 'b-', linewidth=2, label='布朗弛豫')
    ax.loglog(d_range, tau_Ns, 'r-', linewidth=2, label='尼尔弛豫')
    ax.set_xlabel('颗粒直径 [nm]', fontsize=12)
    ax.set_ylabel('弛豫时间 τ [s]', fontsize=12)
    ax.set_title('弛豫时间与颗粒尺寸', fontsize=14)
    ax.legend()
    ax.grid(True, alpha=0.3)
    
    # 交流磁化率
    ax = axes[0, 1]
    frequencies = np.logspace(2, 8, 100)
    chi_real, chi_imag = relax.ac_susceptibility(frequencies)
    
    ax.loglog(frequencies, chi_real, 'b-', linewidth=2, label='实部 χ\'')  # 使用普通单引号
    ax.loglog(frequencies, chi_imag, 'r-', linewidth=2, label='虚部 χ\'\'')  # 使用普通单引号
    ax.set_xlabel('频率 f [Hz]', fontsize=12)
    ax.set_ylabel('磁化率 χ', fontsize=12)
    ax.set_title('交流磁化率频谱 (d=20nm)', fontsize=14)
    ax.legend()
    ax.grid(True, alpha=0.3)
    
    # 主导机制相图
    ax = axes[1, 0]
    mechanisms = []
    for d in d_range:
        params_d = ParticleParams(
            name="Fe3O4", diameter=d, Ms=480e3,
            K=1e4, viscosity=1e-3, temperature=300
        )
        r = MagneticRelaxation(params_d)
        mech = r.dominant_mechanism()
        mechanisms.append(0 if mech == RelaxationMechanism.BROWN else 
                         1 if mech == RelaxationMechanism.NEEL else 0.5)
    
    ax.fill_between(d_range, 0, mechanisms, alpha=0.3)
    ax.plot(d_range, mechanisms, 'g-', linewidth=2)
    ax.set_xlabel('颗粒直径 [nm]', fontsize=12)
    ax.set_ylabel('弛豫机制', fontsize=12)
    ax.set_title('主导弛豫机制 (0=布朗, 1=尼尔)', fontsize=14)
    ax.set_ylim(-0.1, 1.1)
    ax.grid(True, alpha=0.3)
    
    # 温度影响
    ax = axes[1, 1]
    T_range = np.linspace(250, 400, 50)
    tau_B_T = []
    tau_N_T = []
    
    for T in T_range:
        params_T = ParticleParams(
            name="Fe3O4", diameter=20, Ms=480e3,
            K=1e4, viscosity=1e-3, temperature=T
        )
        r = MagneticRelaxation(params_T)
        tau_B_T.append(r.brown_relaxation_time())
        tau_N_T.append(r.neel_relaxation_time())
    
    ax.semilogy(T_range, tau_B_T, 'b-', linewidth=2, label='布朗弛豫')
    ax.semilogy(T_range, tau_N_T, 'r-', linewidth=2, label='尼尔弛豫')
    ax.set_xlabel('温度 [K]', fontsize=12)
    ax.set_ylabel('弛豫时间 τ [s]', fontsize=12)
    ax.set_title('弛豫时间与温度 (d=20nm)', fontsize=14)
    ax.legend()
    ax.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('fig2_relaxation.png', dpi=150, bbox_inches='tight')
    print("  fig2_relaxation.png 已保存")
    plt.close()
    
    # 图3: 磁热疗
    fig3, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    mh = MagneticHyperthermia(params, concentration=1e15)
    
    # SAR vs 磁场
    ax = axes[0, 0]
    H_range = np.linspace(1e3, 50e3, 50)
    frequencies = [100e3, 300e3, 500e3]
    colors_f = ['blue', 'green', 'red']
    
    for f, color in zip(frequencies, colors_f):
        sars = [mh.sar_linear_response(H, f) for H in H_range]
        ax.plot(H_range / 1e3, sars, color=color, linewidth=2, 
                label=f'f = {f/1e3:.0f} kHz')
    
    ax.set_xlabel('磁场 H [kA/m]', fontsize=12)
    ax.set_ylabel('SAR [W/kg]', fontsize=12)
    ax.set_title('SAR与磁场强度', fontsize=14)
    ax.legend()
    ax.grid(True, alpha=0.3)
    
    # SAR vs 频率
    ax = axes[0, 1]
    f_range = np.logspace(4, 6, 50)
    H_values = [5e3, 10e3, 20e3]
    
    for H, color in zip(H_values, colors_f):
        sars = [mh.sar_linear_response(H, f) for f in f_range]
        ax.loglog(f_range / 1e3, sars, color=color, linewidth=2,
                 label=f'H = {H/1e3:.0f} kA/m')
    
    ax.set_xlabel('频率 f [kHz]', fontsize=12)
    ax.set_ylabel('SAR [W/kg]', fontsize=12)
    ax.set_title('SAR与频率', fontsize=14)
    ax.legend()
    ax.grid(True, alpha=0.3)
    
    # 最优尺寸
    ax = axes[1, 0]
    opt_result = mh.optimal_size(20e3, 300e3)
    ax.plot(opt_result['diameters_nm'], opt_result['sar_W_per_kg'], 
            'g-', linewidth=2)
    ax.axvline(x=opt_result['optimal_diameter_nm'], color='r', 
               linestyle='--', label=f'最优尺寸 = {opt_result["optimal_diameter_nm"]:.1f} nm')
    ax.set_xlabel('颗粒直径 [nm]', fontsize=12)
    ax.set_ylabel('SAR [W/kg]', fontsize=12)
    ax.set_title(f'尺寸优化 (H=20kA/m, f=300kHz)', fontsize=14)
    ax.legend()
    ax.grid(True, alpha=0.3)
    
    # 磁滞损耗
    ax = axes[1, 1]
    H_max_range = np.linspace(1e3, 50e3, 30)
    sar_hyst = [mh.sar_hysteresis(H, 300e3) for H in H_max_range]
    sar_linear = [mh.sar_linear_response(H, 300e3) for H in H_max_range]
    
    ax.plot(H_max_range / 1e3, sar_hyst, 'b-', linewidth=2, label='磁滞损耗')
    ax.plot(H_max_range / 1e3, sar_linear, 'r--', linewidth=2, label='线性响应')
    ax.set_xlabel('最大磁场 H_max [kA/m]', fontsize=12)
    ax.set_ylabel('SAR [W/kg]', fontsize=12)
    ax.set_title('SAR计算模型对比', fontsize=14)
    ax.legend()
    ax.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('fig3_hyperthermia.png', dpi=150, bbox_inches='tight')
    print("  fig3_hyperthermia.png 已保存")
    plt.close()
    
    # 图4: 颗粒动力学
    fig4, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    # 定义磁场
    def gradient_field(x, y):
        B0 = 0.1  # T
        L = 1e-6
        Bx = B0 * (x / L)
        By = B0 * (y / L)
        return Bx, By
    
    # 模拟
    dynamics = ParticleDynamics(n_particles=50, box_size=1e-6, params=params)
    result = dynamics.simulate(gradient_field, dt=1e-7, n_steps=500)
    
    # 初始分布
    ax = axes[0, 0]
    ax.scatter(result['trajectory'][0, :, 0] * 1e6, 
              result['trajectory'][0, :, 1] * 1e6,
              c='blue', s=50, alpha=0.6)
    ax.set_xlabel('x [μm]', fontsize=12)
    ax.set_ylabel('y [μm]', fontsize=12)
    ax.set_title('初始分布', fontsize=14)
    ax.set_aspect('equal')
    ax.grid(True, alpha=0.3)
    
    # 最终分布
    ax = axes[0, 1]
    ax.scatter(result['trajectory'][-1, :, 0] * 1e6,
              result['trajectory'][-1, :, 1] * 1e6,
              c='red', s=50, alpha=0.6)
    ax.set_xlabel('x [μm]', fontsize=12)
    ax.set_ylabel('y [μm]', fontsize=12)
    ax.set_title('最终分布', fontsize=14)
    ax.set_aspect('equal')
    ax.grid(True, alpha=0.3)
    
    # 轨迹
    ax = axes[1, 0]
    for i in range(min(10, dynamics.n)):
        traj = result['trajectory'][:, i, :] * 1e6
        ax.plot(traj[:, 0], traj[:, 1], alpha=0.5, linewidth=1)
    ax.set_xlabel('x [μm]', fontsize=12)
    ax.set_ylabel('y [μm]', fontsize=12)
    ax.set_title('颗粒轨迹 (前10个)', fontsize=14)
    ax.set_aspect('equal')
    ax.grid(True, alpha=0.3)
    
    # 均方位移
    ax = axes[1, 1]
    msd = np.zeros(result['n_steps'])
    for i in range(result['n_steps']):
        displacements = result['trajectory'][i] - result['trajectory'][0]
        msd[i] = np.mean(np.sum(displacements**2, axis=1))
    
    time = np.arange(result['n_steps']) * result['dt'] * 1e3  # ms
    ax.plot(time, msd * 1e12, 'b-', linewidth=2)  # nm^2
    
    # 理论扩散 (MSD = 4*D*t for 2D)
    D = 1.38e-23 * 300 / (6 * np.pi * 1e-3 * 10e-9)
    msd_theory = 4 * D * np.arange(result['n_steps']) * result['dt']
    ax.plot(time, msd_theory * 1e12, 'r--', linewidth=2, label='理论扩散')
    
    ax.set_xlabel('时间 [ms]', fontsize=12)
    ax.set_ylabel('MSD [nm²]', fontsize=12)
    ax.set_title('均方位移', fontsize=14)
    ax.legend()
    ax.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('fig4_dynamics.png', dpi=150, bbox_inches='tight')
    print("  fig4_dynamics.png 已保存")
    plt.close()


def save_results(results: Dict, filename: str = 'nanoparticle_results.json'):
    """保存结果到JSON文件"""
    
    def convert_to_serializable(obj):
        if isinstance(obj, np.ndarray):
            return obj.tolist()
        elif isinstance(obj, dict):
            return {k: convert_to_serializable(v) for k, v in obj.items()}
        elif isinstance(obj, list):
            return [convert_to_serializable(item) for item in obj]
        elif isinstance(obj, (np.integer, np.floating)):
            return float(obj)
        elif isinstance(obj, np.bool_):
            return bool(obj)
        elif isinstance(obj, bool):
            return obj
        elif isinstance(obj, Enum):
            return obj.value
        return obj
    
    serializable_results = convert_to_serializable(results)
    
    with open(filename, 'w', encoding='utf-8') as f:
        json.dump(serializable_results, f, ensure_ascii=False, indent=2)
    
    print(f"\n  结果已保存到: {filename}")


# =============================================================================
# 主程序
# =============================================================================

if __name__ == "__main__":
    print("="*70)
    print("主题039: 磁纳米颗粒仿真")
    print("="*70)
    
    results = {}
    
    # 1. 超顺磁性分析
    print("\n[1/4] 超顺磁性分析...")
    
    params = ParticleParams(
        name="Fe3O4",
        diameter=20,  # nm
        Ms=480e3,     # A/m
        K=1e4,        # J/m^3
        viscosity=1e-3,  # Pa·s
        temperature=300  # K
    )
    
    sp = Superparamagnetism(params)
    
    print(f"  颗粒: {params.name}")
    print(f"  直径: {params.diameter} nm")
    print(f"  饱和磁化强度: {params.Ms/1e3:.0f} kA/m")
    print(f"  颗粒体积: {sp.V*1e27:.3f} nm³")
    print(f"  磁矩: {sp.m:.3e} A·m²")
    print(f"  阻塞温度: {sp.blocking_temperature():.1f} K")
    
    # 计算磁化率
    chi = sp.susceptibility(1e3)
    print(f"  初始磁化率: {chi:.3e}")
    
    results['superparamagnetism'] = {
        'material': params.name,
        'diameter_nm': params.diameter,
        'Ms_kA_per_m': params.Ms / 1e3,
        'volume_nm3': float(sp.V * 1e27),
        'moment_A_m2': float(sp.m),
        'blocking_temperature_K': float(sp.blocking_temperature()),
        'susceptibility': float(chi)
    }
    
    # 2. 弛豫分析
    print("\n[2/4] 弛豫分析...")
    
    relax = MagneticRelaxation(params)
    
    tau_B = relax.brown_relaxation_time()
    tau_N = relax.neel_relaxation_time()
    tau_eff = relax.effective_relaxation_time()
    mechanism = relax.dominant_mechanism()
    
    print(f"  布朗弛豫时间: {tau_B:.3e} s")
    print(f"  尼尔弛豫时间: {tau_N:.3e} s")
    print(f"  有效弛豫时间: {tau_eff:.3e} s")
    print(f"  主导机制: {mechanism.value}")
    
    results['relaxation'] = {
        'brown_tau_s': float(tau_B),
        'neel_tau_s': float(tau_N),
        'effective_tau_s': float(tau_eff),
        'dominant_mechanism': mechanism.value
    }
    
    # 3. 磁热疗分析
    print("\n[3/4] 磁热疗分析...")
    
    mh = MagneticHyperthermia(params, concentration=1e15)
    
    # 计算SAR
    H_applied = 20e3  # A/m
    f_applied = 300e3  # Hz
    
    sar_linear = mh.sar_linear_response(H_applied, f_applied)
    sar_hyst = mh.sar_hysteresis(H_applied, f_applied)
    
    print(f"  磁场: {H_applied/1e3:.0f} kA/m")
    print(f"  频率: {f_applied/1e3:.0f} kHz")
    print(f"  SAR (线性响应): {sar_linear:.2f} W/kg")
    print(f"  SAR (磁滞损耗): {sar_hyst:.2f} W/kg")
    
    # 最优尺寸
    opt = mh.optimal_size(H_applied, f_applied)
    print(f"  最优颗粒尺寸: {opt['optimal_diameter_nm']:.1f} nm")
    print(f"  最大SAR: {opt['max_sar_W_per_kg']:.2f} W/kg")
    
    results['hyperthermia'] = {
        'H_kA_per_m': H_applied / 1e3,
        'f_kHz': f_applied / 1e3,
        'sar_linear_W_per_kg': float(sar_linear),
        'sar_hysteresis_W_per_kg': float(sar_hyst),
        'optimal_diameter_nm': float(opt['optimal_diameter_nm']),
        'max_sar_W_per_kg': float(opt['max_sar_W_per_kg'])
    }
    
    # 4. 颗粒动力学
    print("\n[4/4] 颗粒动力学模拟...")
    
    dynamics = ParticleDynamics(n_particles=50, box_size=1e-6, params=params)
    
    def uniform_field(x, y):
        return 0.01, 0.01  # 均匀磁场 0.01 T
    
    result = dynamics.simulate(uniform_field, dt=1e-7, n_steps=500)
    
    print(f"  颗粒数量: {dynamics.n}")
    print(f"  模拟步数: {result['n_steps']}")
    print(f"  时间步长: {result['dt']*1e9:.1f} ns")
    print(f"  总模拟时间: {result['n_steps'] * result['dt'] * 1e6:.2f} μs")
    
    results['dynamics'] = {
        'n_particles': dynamics.n,
        'n_steps': result['n_steps'],
        'dt_ns': float(result['dt'] * 1e9),
        'total_time_us': float(result['n_steps'] * result['dt'] * 1e6)
    }
    
    # 生成可视化
    create_visualizations()
    
    # 保存结果
    save_results(results)
    
    print("\n" + "="*70)
    print("主题039仿真完成！")
    print("="*70)
    print("\n生成的文件:")
    print("  - fig1_superparamagnetism.png: 超顺磁性分析")
    print("  - fig2_relaxation.png: 弛豫分析")
    print("  - fig3_hyperthermia.png: 磁热疗SAR")
    print("  - fig4_dynamics.png: 颗粒动力学")
    print("  - nanoparticle_results.json: 结果数据")
    print("="*70)