多物理场耦合仿真-主题019-磁-流耦合分析

kkchenkx

91人浏览 · 2026-03-25 21:30:47

kkchenkx · 2026-03-25 21:30:47 发布

第十九篇：磁-流耦合分析

摘要

磁-流耦合分析研究磁场与流体流动之间的相互作用，是磁流体力学（MHD）的核心内容。这种耦合现象广泛存在于等离子体物理、冶金工业、核聚变反应堆、磁流体发电和电磁泵送等领域。本主题系统阐述磁-流耦合的基本原理，包括洛伦兹力、感应电流、磁压效应和阿尔芬波等物理机制，建立磁-流耦合的数学模型和控制方程。通过六个典型案例——哈特曼流动分析、磁流体发电通道、电磁泵送系统、磁约束等离子体、磁流体稳定性分析和磁-流耦合有限元实现，深入探讨磁-流耦合问题的数值求解方法。所有案例均配有完整的Python实现代码和可视化结果，帮助读者掌握磁-流耦合仿真的核心技术和工程应用方法。

关键词

磁-流耦合，磁流体力学，洛伦兹力，哈特曼数，磁流体发电，电磁泵送，等离子体约束，阿尔芬波，有限元法，感应电流

在这里插入图片描述

1. 磁-流耦合理论基础

1.1 磁-流耦合现象概述

磁-流耦合是指磁场与导电流体之间相互影响、相互作用的物理现象。当导电流体在磁场中运动时，流体中的自由电荷受到洛伦兹力作用；同时，流体运动切割磁力线产生感应电场和感应电流，这些电流又产生附加磁场，从而改变原有的磁场分布。这种双向耦合构成了磁流体力学（Magnetohydrodynamics, MHD）的研究对象。

洛伦兹力是磁-流耦合的核心物理机制。导电流体中的电荷在电磁场中受到的洛伦兹力密度为：

$f=ρeE+J×B\mathbf{f} = \rho_e \mathbf{E} + \mathbf{J} \times \mathbf{B}$

其中， $ρe\rho_e$ 是电荷密度（C/m³）， $E\mathbf{E}$ 是电场强度（V/m）， $J\mathbf{J}$ 是电流密度（A/m²）， $B\mathbf{B}$ 是磁感应强度（T）。在大多数MHD问题中，电荷密度可以忽略，洛伦兹力简化为：

$f=J×B\mathbf{f} = \mathbf{J} \times \mathbf{B}$

感应电流由法拉第电磁感应定律描述。当导电流体以速度 $u\mathbf{u}$ 在磁场 $B\mathbf{B}$ 中运动时，产生的感应电场为：

$Eind=u×B\mathbf{E}_{ind} = \mathbf{u} \times \mathbf{B}$

根据欧姆定律，流体中的电流密度为：

$J=σ(E+u×B)\mathbf{J} = \sigma (\mathbf{E} + \mathbf{u} \times \mathbf{B})$

其中， $σ\sigma$ 是流体的电导率（S/m）， $E\mathbf{E}$ 是外加电场。

磁压效应是磁-流耦合的重要表现。磁场对导电流体产生等效压力，称为磁压：

$pmag=B22μ0p_{mag} = \frac{B^2}{2\mu_0}$

其中， $μ0\mu_0$ 是真空磁导率（ $4π×10−74\pi \times 10^{-7}$ H/m）。磁压与流体动压的比值定义了磁-流相互作用的重要无量纲参数。

阿尔芬波是磁-流耦合特有的波动现象。当导电流体中存在磁场时，磁力线像弹性弦一样具有张力，可以支持沿磁力线传播的横波。阿尔芬波速为：

$vA=Bμ0ρv_A = \frac{B}{\sqrt{\mu_0 \rho}}$

其中， $ρ\rho$ 是流体密度（kg/m³）。阿尔芬波在等离子体物理和天体物理中有重要应用。

1.2 磁-流耦合控制方程

磁-流耦合问题的控制方程包括流体力学方程和电磁学方程的耦合系统。

连续性方程（质量守恒）：

$∂ρ∂t+∇⋅(ρu)=0\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{u}) = 0$

动量方程（纳维-斯托克斯方程+洛伦兹力）：

$ρ∂u∂t+ρ(u⋅∇)u=−∇p+μ∇2u+J×B+ρg\rho \frac{\partial \mathbf{u}}{\partial t} + \rho (\mathbf{u} \cdot \nabla) \mathbf{u} = -\nabla p + \mu \nabla^2 \mathbf{u} + \mathbf{J} \times \mathbf{B} + \rho \mathbf{g}$

其中， $p$ 是压力（Pa）， $μ\mu$ 是动力粘度（Pa·s）， $g\mathbf{g}$ 是重力加速度（m/s²）。

能量方程：

$ρcp∂T∂t+ρcpu⋅∇T=k∇2T+J2σ+Φ\rho c_p \frac{\partial T}{\partial t} + \rho c_p \mathbf{u} \cdot \nabla T = k \nabla^2 T + \frac{\mathbf{J}^2}{\sigma} + \Phi$

其中， $c_p$ 是比热容（J/(kg·K)）， $k$ 是热导率（W/(m·K)）， $Φ\Phi$ 是粘性耗散函数。

麦克斯韦方程组：

$∇×E=−∂B∂t\nabla \times \mathbf{E} = -\frac{\partial \mathbf{B}}{\partial t}$

$∇×H=J\nabla \times \mathbf{H} = \mathbf{J}$

$∇⋅B=0\nabla \cdot \mathbf{B} = 0$

$J=σ(E+u×B)\mathbf{J} = \sigma (\mathbf{E} + \mathbf{u} \times \mathbf{B})$

其中， $H=B/μ0\mathbf{H} = \mathbf{B}/\mu_0$ 是磁场强度（A/m）。

感应方程（由麦克斯韦方程和欧姆定律推导）：

$∂B∂t=∇×(u×B)+1μ0σ∇2B\frac{\partial \mathbf{B}}{\partial t} = \nabla \times (\mathbf{u} \times \mathbf{B}) + \frac{1}{\mu_0 \sigma} \nabla^2 \mathbf{B}$

这个方程描述了磁场如何被流体运动对流（第一项）和磁扩散（第二项）所改变。

1.3 磁-流耦合无量纲参数

磁-流耦合问题涉及多个重要的无量纲参数：

哈特曼数（Hartmann Number, Ha）：

$\sqrt{\frac{\sigma}{\mu}}$

哈特曼数表征电磁力与粘性力的比值。当Ha>>1时，电磁力主导流动；当Ha<<1时，粘性力主导。

磁雷诺数（Magnetic Reynolds Number, Rm）：

$\mu_0 \sigma u L$

磁雷诺数表征对流输运与磁扩散的比值。当Rm>>1时，磁场被流体强烈扭曲；当Rm<<1时，磁场几乎不受流体影响。

相互作用参数（Interaction Parameter, N）：

$\frac{\sigma B^2 L}{\rho u} = \frac{Ha^2}{Re}$

相互作用参数表征电磁力与惯性力的比值。

阿尔芬马赫数（Alfvén Mach Number, Ma_A）：

$MaA=uvA=uμ0ρBMa_A = \frac{u}{v_A} = u \sqrt{\frac{\mu_0 \rho}{B}}$

阿尔芬马赫数表征流动速度与阿尔芬波速的比值。

斯图亚特数（Stuart Number, St）：

$\frac{\sigma B^2 L}{\rho u} = N$

斯图亚特数与相互作用参数等价，在磁流体发电中常用。

1.4 磁-流耦合数值方法

磁-流耦合问题的数值求解面临以下挑战：

多场耦合：需要同时求解流场、磁场和电场，计算量大。

时间尺度分离：流体对流时间尺度和磁扩散时间尺度可能相差很大，导致刚性问题。

守恒性要求：需要保持磁场无散度（ $∇⋅B=0\nabla \cdot \mathbf{B} = 0$ ）的约束条件。

常用的数值方法包括：

有限体积法（FVM）：在计算流体力学中广泛应用，可以很好地处理对流-扩散问题。

有限元法（FEM）：适用于复杂几何，可以处理多场耦合问题。

谱方法：高精度方法，适用于简单几何和光滑解。

投影法：用于保持磁场无散度约束，通过求解泊松方程修正磁场。

算子分裂法：将耦合问题分解为流体步和磁步，分别求解后再耦合。

2. 案例分析

2.1 案例1：哈特曼流动分析

哈特曼流动是磁-流耦合的经典问题，描述在横向磁场作用下，导电流体在平行平板间的层流流动。

问题描述：

两平行平板间距2a=0.02m
流体：液态汞，电导率σ=1.0×10⁶ S/m，粘度μ=1.5×10⁻³ Pa·s，密度ρ=13500 kg/m³
外加磁场B₀=0.1 T，垂直于平板
压力梯度dp/dx=-100 Pa/m驱动流动
分析速度分布和哈特曼层特性

控制方程：

对于充分发展的稳态流动，x方向动量方程为：

$μd2udy2−σB02u=dpdx\mu \frac{d^2 u}{dy^2} - \sigma B_0^2 u = \frac{dp}{dx}$

解析解：

$\frac{a^2}{\mu Ha^2} \left(-\frac{dp}{dx}\right) \left(1 - \frac{\cosh(Ha \cdot y/a)}{\cosh(Ha)}\right)$

其中，哈特曼数Ha = B₀a√(σ/μ)。

Python实现：

import numpy as np
import matplotlib.pyplot as plt

print("=" * 60)
print("案例1：哈特曼流动分析")
print("=" * 60)

# 物理参数
a = 0.01  # 半间距 [m]
sigma = 1.0e6  # 电导率 [S/m]
mu = 1.5e-3  # 动力粘度 [Pa·s]
rho = 13500  # 密度 [kg/m³]
B0 = 0.1  # 外加磁场 [T]
dp_dx = -100  # 压力梯度 [Pa/m]

# 计算哈特曼数
Ha = B0 * a * np.sqrt(sigma / mu)
print(f"哈特曼数 Ha = {Ha:.2f}")

# 数值求解
ny = 100
y = np.linspace(-a, a, ny)

# 解析解
u_analytical = (a**2 / (mu * Ha**2)) * (-dp_dx) * (1 - np.cosh(Ha * y / a) / np.cosh(Ha))

# 无磁场时的泊肃叶流动（参考）
u_poiseuille = (-dp_dx) * (a**2 - y**2) / (2 * mu)

# 计算边界层厚度（哈特曼层）
delta_H = a / Ha
print(f"哈特曼层厚度 δ_H = {delta_H*1000:.3f} mm")

# 计算壁面剪切应力
tau_w = mu * np.abs(np.gradient(u_analytical, y)[0])
print(f"壁面剪切应力 τ_w = {tau_w:.4f} Pa")

# 计算流量
Q = np.trapz(u_analytical, y)
print(f"单位宽度流量 Q = {Q*1000:.4f} L/(m·s)")

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

# 图1：速度分布
ax1 = axes[0, 0]
ax1.plot(u_analytical, y*1000, 'b-', linewidth=2, label='Hartmann Flow')
ax1.plot(u_poiseuille, y*1000, 'r--', linewidth=2, label='Poiseuille Flow (B=0)')
ax1.set_xlabel('Velocity u (m/s)', fontsize=11)
ax1.set_ylabel('Position y (mm)', fontsize=11)
ax1.set_title(f'Velocity Profile (Ha={Ha:.1f})', fontsize=12, fontweight='bold')
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3)

# 图2：不同哈特曼数下的速度分布
ax2 = axes[0, 1]
Ha_range = [0, 1, 5, 10, 20]
colors = plt.cm.viridis(np.linspace(0, 1, len(Ha_range)))
for i, Ha_test in enumerate(Ha_range):
    if Ha_test == 0:
        u_test = (-dp_dx) * (a**2 - y**2) / (2 * mu)
    else:
        u_test = (a**2 / (mu * Ha_test**2)) * (-dp_dx) * (1 - np.cosh(Ha_test * y / a) / np.cosh(Ha_test))
    ax2.plot(u_test, y*1000, color=colors[i], linewidth=2, label=f'Ha={Ha_test}')
ax2.set_xlabel('Velocity u (m/s)', fontsize=11)
ax2.set_ylabel('Position y (mm)', fontsize=11)
ax2.set_title('Velocity Profiles for Different Ha', fontsize=12, fontweight='bold')
ax2.legend(fontsize=10)
ax2.grid(True, alpha=0.3)

# 图3：哈特曼层结构（近壁区放大）
ax3 = axes[1, 0]
y_wall = y[y >= 0]
u_wall = u_analytical[y >= 0]
ax3.plot(u_wall, y_wall*1000, 'b-', linewidth=2)
ax3.axhline(y=delta_H*1000, color='r', linestyle='--', label=f'δ_H = {delta_H*1000:.2f} mm')
ax3.set_xlabel('Velocity u (m/s)', fontsize=11)
ax3.set_ylabel('Position y (mm)', fontsize=11)
ax3.set_title('Hartmann Layer Structure', fontsize=12, fontweight='bold')
ax3.legend(fontsize=10)
ax3.grid(True, alpha=0.3)

# 图4：流量随哈特曼数变化
ax4 = axes[1, 1]
Ha_range_plot = np.linspace(0.1, 50, 100)
Q_range = []
for Ha_test in Ha_range_plot:
    u_test = (a**2 / (mu * Ha_test**2)) * (-dp_dx) * (1 - np.cosh(Ha_test) / np.cosh(Ha_test))
    # 简化计算
    Q_test = 2 * a * (-dp_dx) * a**2 / (mu * Ha_test**3) * (Ha_test - np.tanh(Ha_test))
    Q_range.append(Q_test)

Q_poiseuille = 2 * a * (-dp_dx) * a**2 / (3 * mu)
ax4.plot(Ha_range_plot, np.array(Q_range)/Q_poiseuille, 'g-', linewidth=2)
ax4.axhline(y=1, color='r', linestyle='--', label='Poiseuille Flow')
ax4.set_xlabel('Hartmann Number Ha', fontsize=11)
ax4.set_ylabel('Normalized Flow Rate Q/Q₀', fontsize=11)
ax4.set_title('Flow Rate vs Hartmann Number', fontsize=12, fontweight='bold')
ax4.legend(fontsize=10)
ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('case1_hartmann_flow.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 案例1完成，图片已保存")

结果分析：

哈特曼流动的主要特征：

速度分布：在哈特曼数较小时，速度分布接近抛物线（泊肃叶流动）；随着哈特曼数增大，速度分布变得平坦，在中心区域近似均匀，在壁面附近形成薄的哈特曼层。
哈特曼层：在壁面附近存在一个特征厚度为δ_H = a/Ha的边界层，在此层内速度迅速降为零。哈特曼层是电磁力与粘性力平衡的结果。
流量抑制：磁场对流动有抑制作用，流量随哈特曼数增加而减小。当Ha>>1时，流量与Ha⁻²成正比。

2.2 案例2：磁流体发电通道

磁流体发电（MHD发电）利用高温导电流体（等离子体或液态金属）通过磁场时产生的感应电动势直接发电。

问题描述：

发电通道尺寸：长度L=1m，宽度w=0.1m，高度h=0.05m
工作介质：燃烧气体，电导率σ=100 S/m，速度u=1000 m/s
外加磁场B₀=2 T，垂直于速度方向
分析发电功率、效率和负载特性

基本原理：

感应电动势： $E = u B_0 w$

输出功率： $Pout=(uB0)2wh(1+RloadRint)2RloadRintP_{out} = \frac{(u B_0)^2 w h}{\left(1 + \frac{R_{load}}{R_{int}}\right)^2} \frac{R_{load}}{R_{int}}$

Python实现：

import numpy as np
import matplotlib.pyplot as plt

print("\n" + "=" * 60)
print("案例2：磁流体发电通道")
print("=" * 60)

# 发电通道参数
L = 1.0  # 长度 [m]
w = 0.1  # 宽度 [m]
h = 0.05  # 高度 [m]
sigma = 100  # 电导率 [S/m]
u = 1000  # 流速 [m/s]
B0 = 2.0  # 磁场强度 [T]

# 计算参数
A_cross = w * h  # 截面积
R_internal = L / (sigma * A_cross)  # 内部电阻
V_open = u * B0 * w  # 开路电压

print(f"通道截面积: {A_cross*1e4:.1f} cm²")
print(f"内部电阻: {R_internal*1e3:.3f} mΩ")
print(f"开路电压: {V_open:.2f} V")

# 负载特性分析
R_load_range = np.linspace(0.1, 10, 100) * R_internal

P_out_range = []
I_range = []
V_out_range = []
eta_range = []

for R_load in R_load_range:
    I = V_open / (R_internal + R_load)
    V_out = I * R_load
    P_out = I**2 * R_load
    eta = R_load / (R_internal + R_load)  # 简化效率
    
    P_out_range.append(P_out)
    I_range.append(I)
    V_out_range.append(V_out)
    eta_range.append(eta)

P_out_range = np.array(P_out_range)
I_range = np.array(I_range)
V_out_range = np.array(V_out_range)
eta_range = np.array(eta_range)

# 最大功率点
max_P_idx = np.argmax(P_out_range)
print(f"\n最大功率输出: {P_out_range[max_P_idx]/1e6:.3f} MW @ R_load={R_load_range[max_P_idx]/R_internal:.2f}R_int")
print(f"对应电流: {I_range[max_P_idx]:.2f} A")
print(f"对应电压: {V_out_range[max_P_idx]:.2f} V")

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

# 图1：功率-负载曲线
ax1 = axes[0, 0]
ax1.plot(R_load_range/R_internal, P_out_range/1e6, 'b-', linewidth=2)
ax1.axvline(x=1, color='r', linestyle='--', label='R_load = R_int')
ax1.scatter(R_load_range[max_P_idx]/R_internal, P_out_range[max_P_idx]/1e6, 
            c='red', s=100, zorder=5, label=f'Max Power = {P_out_range[max_P_idx]/1e6:.2f} MW')
ax1.set_xlabel('Load Resistance Ratio R_load/R_int', fontsize=11)
ax1.set_ylabel('Output Power (MW)', fontsize=11)
ax1.set_title('Output Power vs Load Resistance', fontsize=12, fontweight='bold')
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3)

# 图2：电流-电压特性
ax2 = axes[0, 1]
ax2.plot(I_range, V_out_range, 'g-', linewidth=2)
ax2.scatter(I_range[max_P_idx], V_out_range[max_P_idx], 
            c='red', s=100, zorder=5, label='Max Power Point')
ax2.set_xlabel('Current (A)', fontsize=11)
ax2.set_ylabel('Voltage (V)', fontsize=11)
ax2.set_title('I-V Characteristic', fontsize=12, fontweight='bold')
ax2.legend(fontsize=10)
ax2.grid(True, alpha=0.3)

# 图3：效率-负载曲线
ax3 = axes[1, 0]
ax3.plot(R_load_range/R_internal, eta_range*100, 'm-', linewidth=2)
ax3.set_xlabel('Load Resistance Ratio R_load/R_int', fontsize=11)
ax3.set_ylabel('Efficiency (%)', fontsize=11)
ax3.set_title('Efficiency vs Load Resistance', fontsize=12, fontweight='bold')
ax3.grid(True, alpha=0.3)

# 图4：不同磁场强度下的功率
ax4 = axes[1, 1]
B_range = np.linspace(0.5, 5, 50)
P_max_range = []
for B_test in B_range:
    V_open_test = u * B_test * w
    P_max_test = V_open_test**2 / (4 * R_internal)
    P_max_range.append(P_max_test)

ax4.plot(B_range, np.array(P_max_range)/1e6, 'c-', linewidth=2)
ax4.scatter(B0, P_out_range[max_P_idx]/1e6, c='red', s=100, zorder=5)
ax4.set_xlabel('Magnetic Field B (T)', fontsize=11)
ax4.set_ylabel('Max Output Power (MW)', fontsize=11)
ax4.set_title('Max Power vs Magnetic Field', fontsize=12, fontweight='bold')
ax4.grid(True, alpha=0.3)

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

结果分析：

磁流体发电的关键特性：

负载匹配：最大功率输出发生在负载电阻等于内部电阻时（R_load = R_int），这与电路理论一致。
功率密度：输出功率与磁场强度的平方成正比，与流速的平方成正比。提高磁场强度和流速是提高功率密度的关键。
效率限制：磁流体发电的理论效率受卡诺效率限制，实际效率还受限于电导率、热损失等因素。

2.3 案例3：电磁泵送系统

电磁泵利用洛伦兹力驱动导电流体流动，广泛应用于液态金属冷却反应堆、冶金工业等领域。

问题描述：

泵送通道：矩形截面，宽度w=0.05m，高度h=0.02m，长度L=0.3m
流体：液态钠，电导率σ=1.0×10⁷ S/m，密度ρ=850 kg/m³
外加磁场B₀=0.5 T，外加电场E₀=100 V/m
分析泵送压力、流量和效率

工作原理：

洛伦兹力驱动： $\times B = \sigma (E + u \times B) \times B$

泵送压力： $Δp=JBL=σ(E0−uB0)B0L\Delta p = J B L = \sigma (E_0 - u B_0) B_0 L$

Python实现：

import numpy as np
import matplotlib.pyplot as plt

print("\n" + "=" * 60)
print("案例3：电磁泵送系统")
print("=" * 60)

# 电磁泵参数
w = 0.05  # 宽度 [m]
h = 0.02  # 高度 [m]
L = 0.3  # 长度 [m]
sigma = 1.0e7  # 电导率 [S/m]
rho = 850  # 密度 [kg/m³]
B0 = 0.5  # 磁场强度 [T]
E0 = 100  # 外加电场 [V/m]

# 计算参数
A_cross = w * h
print(f"通道截面积: {A_cross*1e4:.1f} cm²")

# 分析不同流速下的性能
u_range = np.linspace(0, 20, 100)

delta_p_range = []
J_range = []
P_mech_range = []
P_elec_range = []
eta_range = []

for u in u_range:
    # 感应电场
    E_ind = u * B0
    # 净电场
    E_net = E0 - E_ind
    # 电流密度
    J = sigma * E_net
    # 泵送压力
    delta_p = J * B0 * L
    # 机械功率
    P_mech = delta_p * A_cross * u
    # 电功率
    P_elec = J * E0 * A_cross * L
    # 效率
    eta = P_mech / P_elec if P_elec > 0 else 0
    
    delta_p_range.append(delta_p)
    J_range.append(J)
    P_mech_range.append(P_mech)
    P_elec_range.append(P_elec)
    eta_range.append(eta)

delta_p_range = np.array(delta_p_range)
J_range = np.array(J_range)
P_mech_range = np.array(P_mech_range)
P_elec_range = np.array(P_elec_range)
eta_range = np.array(eta_range)

# 最大压力点（零流速）
print(f"最大泵送压力: {delta_p_range[0]/1e5:.3f} bar")
print(f"最大电流密度: {J_range[0]/1e4:.2f} A/cm²")

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

# 图1：压力-流量曲线
ax1 = axes[0, 0]
Q_range = u_range * A_cross * 1000  # L/s
ax1.plot(Q_range, delta_p_range/1e5, 'b-', linewidth=2)
ax1.set_xlabel('Flow Rate Q (L/s)', fontsize=11)
ax1.set_ylabel('Pressure Difference Δp (bar)', fontsize=11)
ax1.set_title('Pressure-Flow Characteristic', fontsize=12, fontweight='bold')
ax1.grid(True, alpha=0.3)

# 图2：电流密度-流速关系
ax2 = axes[0, 1]
ax2.plot(u_range, J_range/1e4, 'r-', linewidth=2)
ax2.axhline(y=0, color='k', linestyle='-', linewidth=0.5)
ax2.set_xlabel('Flow Velocity u (m/s)', fontsize=11)
ax2.set_ylabel('Current Density J (A/cm²)', fontsize=11)
ax2.set_title('Current Density vs Flow Velocity', fontsize=12, fontweight='bold')
ax2.grid(True, alpha=0.3)

# 图3：功率-流速曲线
ax3 = axes[1, 0]
ax3.plot(u_range, P_mech_range/1e3, 'g-', linewidth=2, label='Mechanical Power')
ax3.plot(u_range, P_elec_range/1e3, 'm--', linewidth=2, label='Electrical Power')
ax3.set_xlabel('Flow Velocity u (m/s)', fontsize=11)
ax3.set_ylabel('Power (kW)', fontsize=11)
ax3.set_title('Power vs Flow Velocity', fontsize=12, fontweight='bold')
ax3.legend(fontsize=10)
ax3.grid(True, alpha=0.3)

# 图4：效率-流速曲线
ax4 = axes[1, 1]
valid_eta = eta_range >= 0
ax4.plot(u_range[valid_eta], eta_range[valid_eta]*100, 'c-', linewidth=2)
ax4.set_xlabel('Flow Velocity u (m/s)', fontsize=11)
ax4.set_ylabel('Efficiency (%)', fontsize=11)
ax4.set_title('Efficiency vs Flow Velocity', fontsize=12, fontweight='bold')
ax4.grid(True, alpha=0.3)

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

结果分析：

电磁泵的性能特点：

压力-流量特性：电磁泵的压力随流量增加而线性下降，类似于直流电机的转矩-转速特性。
电流特性：当流速增加时，感应电动势抵消外加电场，电流密度减小。当u = E₀/B₀时，电流为零，此时达到最大流速。
效率优化：电磁泵的效率与流速有关，存在最优工作点。设计时需要平衡压力、流量和效率。

2.4 案例4：磁约束等离子体

磁约束是控制高温等离子体的主要方法，应用于托卡马克等核聚变装置。

问题描述：

圆柱形等离子体，半径a=0.5m，长度L=2m
等离子体密度n=1.0×10²⁰ m⁻³，温度T=10 keV
外加纵向磁场B_z=5 T
分析磁镜效应、粒子轨道和约束时间

物理模型：

磁矩守恒： $μ=mv⊥22B=const\mu = \frac{m v_\perp^2}{2B} = const$

磁镜比： $Rm=BmaxBminR_m = \frac{B_{max}}{B_{min}}$

损失锥角度： $sin⁡2θc=1Rm\sin^2 \theta_c = \frac{1}{R_m}$

Python实现：

import numpy as np
import matplotlib.pyplot as plt

print("\n" + "=" * 60)
print("案例4：磁约束等离子体")
print("=" * 60)

# 等离子体参数
a = 0.5  # 半径 [m]
L = 2.0  # 长度 [m]
n = 1.0e20  # 密度 [m^-3]
T_keV = 10  # 温度 [keV]
T = T_keV * 1.602e-16 / 1.38e-23  # 转换为K [约1.16e8 K]
B_z = 5.0  # 纵向磁场 [T]

# 基本物理常数
m_e = 9.109e-31  # 电子质量 [kg]
m_i = 1.673e-27  # 质子质量 [kg]
e = 1.602e-19  # 元电荷 [C]
k_B = 1.38e-23  # 玻尔兹曼常数 [J/K]
mu_0 = 4 * np.pi * 1e-7  # 真空磁导率

# 计算等离子体参数
v_th = np.sqrt(2 * k_B * T / m_i)  # 热速度
omega_c = e * B_z / m_i  # 回旋频率
r_L = v_th / omega_c  # 拉莫尔半径

print(f"等离子体温度: {T_keV} keV ({T:.2e} K)")
print(f"离子热速度: {v_th/1e6:.2f} × 10^6 m/s")
print(f"回旋频率: {omega_c/1e8:.2f} × 10^8 rad/s")
print(f"拉莫尔半径: {r_L*1000:.3f} mm")

# 磁镜场分布（简化模型）
z = np.linspace(-L/2, L/2, 200)
B_min = B_z
B_max = B_z * 2.5  # 磁镜比2.5
B_field = B_min + (B_max - B_min) * (z / (L/2))**2

# 磁镜比
R_m = B_max / B_min
print(f"磁镜比: {R_m:.2f}")

# 损失锥角度
theta_c = np.arcsin(np.sqrt(1 / R_m))
print(f"损失锥角度: {np.degrees(theta_c):.2f}°")

# 粒子轨道模拟（简化）
theta_particle = np.linspace(0, np.pi/2, 100)
v_parallel = v_th * np.cos(theta_particle)
v_perp = v_th * np.sin(theta_particle)

# 判断约束条件
mu = 0.5 * m_i * v_perp**2 / B_min  # 磁矩
B_reflection = 0.5 * m_i * v_perp**2 / mu  # 反射点磁场
confinement = B_reflection <= B_max  # 是否被约束

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

# 图1：磁场分布
ax1 = axes[0, 0]
ax1.plot(z, B_field, 'b-', linewidth=2)
ax1.axhline(y=B_min, color='g', linestyle='--', label=f'B_min = {B_min} T')
ax1.axhline(y=B_max, color='r', linestyle='--', label=f'B_max = {B_max} T')
ax1.fill_between(z, B_min, B_field, alpha=0.3, color='blue')
ax1.set_xlabel('Position z (m)', fontsize=11)
ax1.set_ylabel('Magnetic Field B (T)', fontsize=11)
ax1.set_title('Magnetic Mirror Field Distribution', fontsize=12, fontweight='bold')
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3)

# 图2：损失锥
ax2 = axes[0, 1]
theta_deg = np.degrees(theta_particle)
ax2.plot(theta_deg, confinement.astype(int), 'b-', linewidth=2)
ax2.axvline(x=np.degrees(theta_c), color='r', linestyle='--', 
            label=f'Loss Cone θ_c = {np.degrees(theta_c):.1f}°')
ax2.fill_between(theta_deg, 0, confinement.astype(int), 
                 where=(theta_deg >= np.degrees(theta_c)), alpha=0.3, color='red')
ax2.set_xlabel('Pitch Angle θ (degrees)', fontsize=11)
ax2.set_ylabel('Confinement (1=confined, 0=loss)', fontsize=11)
ax2.set_title('Loss Cone and Particle Confinement', fontsize=12, fontweight='bold')
ax2.legend(fontsize=10)
ax2.grid(True, alpha=0.3)

# 图3：不同磁镜比下的损失锥
ax3 = axes[1, 0]
R_m_range = np.linspace(1.1, 10, 100)
theta_c_range = np.degrees(np.arcsin(np.sqrt(1 / R_m_range)))
ax3.plot(R_m_range, theta_c_range, 'g-', linewidth=2)
ax3.scatter(R_m, np.degrees(theta_c), c='red', s=100, zorder=5, 
            label=f'Current R_m = {R_m:.1f}')
ax3.set_xlabel('Mirror Ratio R_m', fontsize=11)
ax3.set_ylabel('Loss Cone Angle θ_c (degrees)', fontsize=11)
ax3.set_title('Loss Cone Angle vs Mirror Ratio', fontsize=12, fontweight='bold')
ax3.legend(fontsize=10)
ax3.grid(True, alpha=0.3)

# 图4：粒子轨道（简化示意）
ax4 = axes[1, 1]
theta_orbit = np.linspace(0, 2*np.pi, 100)
r_orbit = r_L
x_orbit = r_orbit * np.cos(theta_orbit)
y_orbit = r_orbit * np.sin(theta_orbit)

ax4.plot(x_orbit*1000, y_orbit*1000, 'b-', linewidth=2, label='Larmor Orbit')
ax4.scatter(0, 0, c='red', s=50, zorder=5, label='Guiding Center')
ax4.set_xlabel('x (mm)', fontsize=11)
ax4.set_ylabel('y (mm)', fontsize=11)
ax4.set_title(f'Particle Larmor Orbit (r_L = {r_L*1000:.2f} mm)', fontsize=12, fontweight='bold')
ax4.legend(fontsize=10)
ax4.grid(True, alpha=0.3)
ax4.axis('equal')

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

结果分析：

磁约束等离子体的关键特性：

磁镜效应：在磁场强度较大的区域，粒子受到反射力，形成磁镜。磁镜比越大，损失锥越小，约束效果越好。
损失锥：只有速度方向在损失锥内的粒子才能逃逸。提高磁镜比可以减小损失锥，改善约束。
拉莫尔轨道：带电粒子在均匀磁场中做圆周运动，轨道半径与粒子能量成正比，与磁场强度成反比。

2.5 案例5：磁流体稳定性分析

磁流体稳定性是磁约束聚变和天体物理中的重要问题。本案例分析瑞利-泰勒不稳定性和开尔文-亥姆霍兹不稳定性。

问题描述：

两层不可压缩导电流体界面
上层流体密度ρ₁=1000 kg/m³，下层ρ₂=2000 kg/m³
界面处存在磁场B₀=0.1 T，平行于界面
分析磁场对界面稳定性的影响

稳定性判据：

瑞利-泰勒增长率： $γ2=gkρ2−ρ1ρ2+ρ1−(k⋅B0)2μ0(ρ2+ρ1)\gamma^2 = g k \frac{\rho_2 - \rho_1}{\rho_2 + \rho_1} - \frac{(k \cdot B_0)^2}{\mu_0 (\rho_2 + \rho_1)}$

Python实现：

import numpy as np
import matplotlib.pyplot as plt

print("\n" + "=" * 60)
print("案例5：磁流体稳定性分析")
print("=" * 60)

# 流体参数
rho_1 = 1000  # 上层密度 [kg/m³]
rho_2 = 2000  # 下层密度 [kg/m³]
g = 9.8  # 重力加速度 [m/s²]
B0 = 0.1  # 磁场强度 [T]
mu_0 = 4 * np.pi * 1e-7  # 真空磁导率

# 波数范围
k = np.linspace(0.1, 50, 200)

# 瑞利-泰勒不稳定性增长率
# 无磁场情况
gamma_RT_no_B = np.sqrt(g * k * (rho_2 - rho_1) / (rho_2 + rho_1))

# 有磁场情况（磁场平行于界面）
gamma_RT_with_B = np.sqrt(np.maximum(0, g * k * (rho_2 - rho_1) / (rho_2 + rho_1) - 
                                     (k * B0)**2 / (mu_0 * (rho_2 + rho_1))))

# 临界波数
k_critical = np.sqrt(g * (rho_2 - rho_1) * mu_0) / B0
print(f"临界波数 k_c = {k_critical:.3f} m^-1")
print(f"临界波长 λ_c = {2*np.pi/k_critical:.3f} m")

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

# 图1：增长率-波数曲线
ax1 = axes[0, 0]
ax1.plot(k, gamma_RT_no_B, 'b-', linewidth=2, label='Without Magnetic Field')
ax1.plot(k, gamma_RT_with_B, 'r--', linewidth=2, label='With Magnetic Field')
ax1.axvline(x=k_critical, color='g', linestyle=':', label=f'k_critical = {k_critical:.2f}')
ax1.set_xlabel('Wavenumber k (m^-1)', fontsize=11)
ax1.set_ylabel('Growth Rate γ (s^-1)', fontsize=11)
ax1.set_title('Rayleigh-Taylor Instability Growth Rate', fontsize=12, fontweight='bold')
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3)

# 图2：不同磁场强度下的稳定性
ax2 = axes[0, 1]
B_range = [0, 0.05, 0.1, 0.2, 0.5]
colors = plt.cm.viridis(np.linspace(0, 1, len(B_range)))
for i, B_test in enumerate(B_range):
    if B_test == 0:
        gamma_test = np.sqrt(g * k * (rho_2 - rho_1) / (rho_2 + rho_1))
    else:
        gamma_test = np.sqrt(np.maximum(0, g * k * (rho_2 - rho_1) / (rho_2 + rho_1) - 
                                       (k * B_test)**2 / (mu_0 * (rho_2 + rho_1))))
    ax2.plot(k, gamma_test, color=colors[i], linewidth=2, label=f'B₀={B_test} T')
ax2.set_xlabel('Wavenumber k (m^-1)', fontsize=11)
ax2.set_ylabel('Growth Rate γ (s^-1)', fontsize=11)
ax2.set_title('Growth Rate for Different B₀', fontsize=12, fontweight='bold')
ax2.legend(fontsize=10)
ax2.grid(True, alpha=0.3)

# 图3：界面扰动发展（简化）
ax3 = axes[1, 0]
x_interface = np.linspace(0, 2*np.pi, 100)
times = [0, 0.5, 1.0, 1.5]
for i, t in enumerate(times):
    amplitude = 0.1 * np.exp(gamma_RT_with_B[10] * t)
    y_interface = amplitude * np.sin(x_interface)
    ax3.plot(x_interface, y_interface, linewidth=2, label=f't = {t} s')
ax3.set_xlabel('x (rad)', fontsize=11)
ax3.set_ylabel('Interface Displacement (m)', fontsize=11)
ax3.set_title('Interface Evolution (Simplified)', fontsize=12, fontweight='bold')
ax3.legend(fontsize=10)
ax3.grid(True, alpha=0.3)

# 图4：稳定性区域图
ax4 = axes[1, 1]
k_range = np.linspace(0.1, 100, 100)
B_range_plot = np.linspace(0.01, 0.5, 100)
K, B = np.meshgrid(k_range, B_range_plot)
Stability = g * (rho_2 - rho_1) - (K * B)**2 / mu_0

im = ax4.contourf(K, B, Stability, levels=20, cmap='RdYlGn')
ax4.contour(K, B, Stability, levels=[0], colors='black', linewidths=2)
ax4.set_xlabel('Wavenumber k (m^-1)', fontsize=11)
ax4.set_ylabel('Magnetic Field B₀ (T)', fontsize=11)
ax4.set_title('Stability Diagram', fontsize=12, fontweight='bold')
plt.colorbar(im, ax=ax4, label='Stability Parameter')

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

结果分析：

磁流体稳定性的主要结论：

磁场稳定作用：磁场对瑞利-泰勒不稳定性有抑制作用。当磁场强度足够大时，可以稳定短波扰动。
临界波数：存在一个临界波数，小于该波数的扰动被稳定，大于该波数的扰动仍然不稳定。
稳定性区域：在k-B参数空间中，存在稳定区域和不稳定区域的分界线。

2.6 案例6：磁-流耦合有限元实现

本案例实现磁-流耦合问题的有限元求解，展示全耦合数值方法。

问题描述：

二维方腔，边长L=0.1m
导电流体，电导率σ=100 S/m，粘度μ=0.01 Pa·s，密度ρ=1000 kg/m³
外加磁场B₀=0.5 T，垂直方向
上壁面以速度U=0.01 m/s移动（顶盖驱动流）
求解稳态速度场和感应磁场

Python实现：

import numpy as np
import matplotlib.pyplot as plt
from scipy.sparse import csr_matrix
from scipy.sparse.linalg import spsolve

print("\n" + "=" * 60)
print("案例6：磁-流耦合有限元实现")
print("=" * 60)

# 物理参数
L = 0.1  # 方腔边长 [m]
sigma = 100  # 电导率 [S/m]
mu = 0.01  # 粘度 [Pa·s]
rho = 1000  # 密度 [kg/m³]
B0 = 0.5  # 外加磁场 [T]
U = 0.01  # 顶盖速度 [m/s]

# 计算哈特曼数
Ha = B0 * L * np.sqrt(sigma / mu)
print(f"哈特曼数 Ha = {Ha:.2f}")

# 计算磁雷诺数（假设特征速度为U）
mu_0 = 4 * np.pi * 1e-7
Rm = mu_0 * sigma * U * L
print(f"磁雷诺数 Rm = {Rm:.2e}")

# 简化求解：使用有限差分法求解耦合问题
nx = 50
ny = 50
dx = L / (nx - 1)
dy = L / (ny - 1)

x = np.linspace(0, L, nx)
y = np.linspace(0, L, ny)
X, Y = np.meshgrid(x, y)

# 初始化场
u = np.zeros((ny, nx))  # x方向速度
v = np.zeros((ny, nx))  # y方向速度
B_x = np.zeros((ny, nx))  # x方向感应磁场
B_y = np.zeros((ny, nx))  # y方向感应磁场

# 迭代求解
max_iter = 5000
tolerance = 1e-6
for iteration in range(max_iter):
    u_old = u.copy()
    v_old = v.copy()
    
    # 内部节点更新（简化处理）
    for i in range(1, ny-1):
        for j in range(1, nx-1):
            # 速度更新（包含洛伦兹力项）
            # 洛伦兹力: f_x = J_y * B0, J_y = sigma * (E_y + u * B0)
            # 简化: 假设感应电场较小，主要考虑u*B0项
            J_y = sigma * u_old[i, j] * B0
            f_x = J_y * B0
            
            # 简化的速度更新
            u[i, j] = 0.25 * (u_old[i+1, j] + u_old[i-1, j] + 
                              u_old[i, j+1] + u_old[i, j-1]) - \
                      0.25 * dx**2 * f_x / mu
            
            v[i, j] = 0.25 * (v_old[i+1, j] + v_old[i-1, j] + 
                              v_old[i, j+1] + v_old[i, j-1])
    
    # 边界条件
    # 顶盖驱动
    u[-1, :] = U
    v[-1, :] = 0
    # 其他壁面无滑移
    u[0, :] = 0  # 底面
    u[:, 0] = 0  # 左面
    u[:, -1] = 0  # 右面
    v[0, :] = 0
    v[:, 0] = 0
    v[:, -1] = 0
    
    # 检查收敛
    if np.max(np.abs(u - u_old)) < tolerance and np.max(np.abs(v - v_old)) < tolerance:
        print(f"收敛于迭代次数: {iteration}")
        break

# 计算感应电流密度
J_z = sigma * u * B0

print(f"最大x速度: {np.max(u):.6f} m/s")
print(f"最大感应电流: {np.max(J_z):.4f} A/m²")

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

# 图1：速度场云图
ax1 = axes[0, 0]
im1 = ax1.contourf(X*1000, Y*1000, u, levels=20, cmap='RdBu_r')
ax1.set_xlabel('x (mm)', fontsize=11)
ax1.set_ylabel('y (mm)', fontsize=11)
ax1.set_title('Velocity u Distribution (m/s)', fontsize=12, fontweight='bold')
plt.colorbar(im1, ax=ax1, label='u (m/s)')

# 图2：速度矢量图
ax2 = axes[0, 1]
skip = 3
ax2.quiver(X[::skip, ::skip]*1000, Y[::skip, ::skip]*1000, 
           u[::skip, ::skip], v[::skip, ::skip], scale=0.5)
ax2.set_xlabel('x (mm)', fontsize=11)
ax2.set_ylabel('y (mm)', fontsize=11)
ax2.set_title('Velocity Vector Field', fontsize=12, fontweight='bold')
ax2.set_aspect('equal')

# 图3：感应电流密度
ax3 = axes[1, 0]
im3 = ax3.contourf(X*1000, Y*1000, J_z, levels=20, cmap='hot')
ax3.set_xlabel('x (mm)', fontsize=11)
ax3.set_ylabel('y (mm)', fontsize=11)
ax3.set_title('Induced Current Density J_z (A/m²)', fontsize=12, fontweight='bold')
plt.colorbar(im3, ax=ax3, label='J_z (A/m²)')

# 图4：沿中心线的速度分布
ax4 = axes[1, 1]
center_idx = nx // 2
ax4.plot(u[:, center_idx], y*1000, 'b-', linewidth=2, label='With B field')
# 无磁场参考（简化）
u_no_B = U * np.sin(np.pi * y / L) * np.sin(np.pi * y / L)
ax4.plot(u_no_B*0.5, y*1000, 'r--', linewidth=2, label='Without B field (ref)')
ax4.set_xlabel('Velocity u (m/s)', fontsize=11)
ax4.set_ylabel('y (mm)', fontsize=11)
ax4.set_title('Velocity Profile at Center', fontsize=12, fontweight='bold')
ax4.legend(fontsize=10)
ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('case6_mhd_fem.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 案例6完成，图片已保存")

结果分析：

磁-流耦合有限元实现的关键要点：

耦合机制：流体运动产生感应电流，感应电流与外加磁场相互作用产生洛伦兹力，洛伦兹力改变流体运动。
哈特曼层：在壁面附近形成哈特曼层，速度梯度较大。
数值挑战：磁-流耦合问题需要同时求解流场和电磁场，计算量大，需要采用适当的迭代策略。

AtomGit开源社区

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

更多推荐

Harness Engineering：继 Context Engineering 之后，AI Agent 时代的新工程范式

摘要： 2026年2月提出的Harness Engineering标志着AI交互的第三阶段，从Prompt Engineering（优化提问）到Context Engineering（优化上下文），再到Harness Engineering（优化系统环境）。其核心是通过外部基础设施（如约束、验证、工具集成）确保Agent的长期稳定运行，而非依赖模型改进。OpenAI和CLI-Anything等项目

AtomGit开源社区

AI编程助手FittenCode：高效开发新利器

这两年 AI 发展迅猛，作为开发人员，我们总是追求更快、更高效的工作方式，AI 的出现可以说改变了很多人的编程方式。AI 对我们来说就是一个可靠的编程助手，给我们提供了实时的建议和解决方案，无论是快速修复错误、提升代码质量，或者查找关键文档和资源，AI 作为编程助手都能让你事半功倍。

AtomGit开源社区

CLI-Anything 只能用于桌面端软件吗？如何为自己的软件生成 Agent 可用的 CLI

CLI-Anything 工具支持为多种软件形态生成命令行界面，包括桌面应用、本地/云端服务、AI推理引擎等。通过分析23个现有实例，文章归纳出6种后端集成范式：子进程调用、REST API交互（本地/云端）、MCP协议、封装已有CLI及本地文件解析。使用决策树可判断软件适用哪种范式，只需具备可编程接口即可接入。实操部分介绍了从代码分析到发布的全流程，包括单次生成、迭代优化和测试验证。该工具突破了