第二十六篇：压电耦合分析

摘要

压电耦合是多物理场耦合仿真中的重要研究领域，涉及机械应力与电场之间的双向相互作用。本主题系统阐述压电效应的物理机制，包括正压电效应和逆压电效应的基本原理。详细推导压电材料的本构方程，建立力-电耦合的数学模型。介绍压电传感器、执行器、能量收集器等典型器件的仿真方法。通过六个工程案例的Python实现，展示压电耦合问题的数值求解技术，包括静力分析、谐响应分析、瞬态动力学仿真等。特别地，本主题将生成GIF动图展示压电振子的动态响应过程。最后探讨压电耦合在智能结构、微机电系统、能量收集等前沿领域的应用前景。

关键词

压电耦合，压电效应，压电传感器，压电执行器，能量收集，力-电耦合，谐响应分析，有限元方法

---

1. 引言

1.1 压电效应概述

压电效应是指某些晶体材料在受到机械应力时产生电极化（正压电效应），或在施加电场时产生机械应变（逆压电效应）的物理现象。这一效应于1880年由居里兄弟首次发现，现已成为现代工程应用中不可或缺的技术基础。

压电效应的物理本质在于晶体结构的不对称性。在32种晶体点群中，有20种具有压电性。常见的压电材料包括：

单晶材料：石英（SiO₂）、罗谢尔盐、磷酸二氢钾（KDP）等，具有优异的稳定性和低损耗特性。

陶瓷材料：钛酸钡（BaTiO₃）、锆钛酸铅（PZT）等，具有高压电系数和易于加工成型的优点，是应用最广泛的压电材料。

聚合物材料：聚偏二氟乙烯（PVDF）等柔性压电材料，适用于曲面结构和生物医学应用。

复合材料：压电陶瓷与聚合物复合，兼顾高压电性能和机械柔性。

1.2 压电耦合的工程应用

压电耦合在众多高科技领域具有关键应用价值：

传感器应用：压电加速度计、压力传感器、力传感器等，利用正压电效应将机械量转换为电信号，具有响应快、灵敏度高的特点。

执行器应用：精密定位平台、喷墨打印头、超声电机等，利用逆压电效应实现纳米级精度的位移控制。

能量收集：从环境振动中收集能量，为无线传感器网络供电，是可再生能源领域的重要研究方向。

超声技术：超声换能器、医学成像、无损检测等，利用压电材料的高频振动特性。

智能结构：压电层合板、主动振动控制、形状自适应结构等，实现结构的智能化。

1.3 本主题内容安排

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

第2节阐述压电效应的物理机制，包括晶体结构、本构关系等；第3节建立压电耦合的数学模型，推导力-电耦合的控制方程；第4节介绍数值求解方法，包括有限元离散、谐响应分析、瞬态求解等；第5节通过六个工程案例展示压电耦合分析的Python实现，包含GIF动图展示动态过程；第6节讨论结果分析与工程应用；第7节总结全文并展望未来发展方向。

---

2. 压电效应物理机制

2.1 晶体结构与对称性

2.1.1 压电性的晶体学基础

压电性要求晶体结构不具有对称中心。在受到机械应力时，正负电荷中心发生相对位移，产生宏观电极化。对于具有对称中心的晶体，正负电荷中心的位移相互抵消，不产生净极化。

压电张量的独立分量数目取决于晶体的对称性：

• 三斜晶系：18个独立分量
• 单斜晶系：10个独立分量
• 正交晶系：5个独立分量
• 四方晶系（4mm）：3个独立分量
• 六方晶系（6mm）：3个独立分量
• 立方晶系（43m）：1个独立分量

2.1.2 极化过程

铁电陶瓷（如PZT）需要经过极化处理才能表现出压电性。极化过程包括：

1. 升温：将材料加热到居里温度以上，使电畴随机分布
2. 施加强电场：在居里温度以下施加直流电场，使电畴沿电场方向取向
3. 冷却：在保持电场的条件下冷却到室温，冻结极化状态

极化后的材料具有剩余极化强度$P_r$，表现出宏观压电性。

2.2 压电本构方程

2.2.1 基本本构关系

压电材料的本构方程描述了应力$T$、应变$S$、电场$E$和电位移$D$之间的关系。常用的形式有：

应变-电场形式（第一类本构方程）：

$$S_{ij}=s_{ijkl}^E{T}_{kl}+d_{kij}{E}_k$$

$$D_i=d_{ikl}{T}_{kl}+{\epsilon }_{ik}^T{E}_k$$

应力-电场形式（第二类本构方程）：

$$T_{ij}=c_{ijkl}^E{S}_{kl}-e_{kij}{E}_k$$

$$D_i=e_{ikl}{S}_{kl}+{\epsilon }_{ik}^S{E}_k$$

其中，$s_{ijkl}^E$为恒电场柔顺系数，$c_{ijkl}^E$为恒电场弹性系数，$d_{kij}$为压电应变系数，$e_{kij}$为压电应力系数，${\epsilon }_{ik}^T$为恒应力介电常数，${\epsilon }_{ik}^S$为恒应变介电常数。

2.2.2 矩阵表示

对于工程应用，常采用简化的矩阵表示（Voigt记号）：

$$\left\{\begin{array}{c}S\\ D\end{array}\right\}=\left[\begin{array}{cc}{s}^E& d^t\\ d& {\epsilon }^T\end{array}\right]\left\{\begin{array}{c}T\\ E\end{array}\right\}$$

或

$$\left\{\begin{array}{c}T\\ D\end{array}\right\}=\left[\begin{array}{cc}{c}^E& -e^t\\ e& {\epsilon }^S\end{array}\right]\left\{\begin{array}{c}S\\ E\end{array}\right\}$$

对于PZT-5H材料，典型参数值为：

• $c_{11}^E=126$ GPa, $c_{12}^E=79.5$ GPa, $c_{13}^E=84.1$ GPa
• $c_{33}^E=117$ GPa, $c_{44}^E=23.0$ GPa
• $e_{31}=-6.5$ C/m², $e_{33}=23.3$ C/m², $e_{15}=17.0$ C/m²
• ${\epsilon }_{11}^S=1700{\epsilon }_0$, ${\epsilon }_{33}^S=1470{\epsilon }_0$

2.3 压电材料的机电耦合系数

机电耦合系数$k$是衡量压电材料能量转换效率的重要参数：

$$k^2=\frac{\text{转换的机械能}}{\text{输入的电能}}=\frac{\text{转换的电能}}{\text{输入的机械能}}$$

对于不同振动模式：

纵向振动：$k_{33}=\frac{d_{33}}{\sqrt{s_{33}^E{\epsilon }_{33}^T}}$

横向振动：$k_{31}=\frac{d_{31}}{\sqrt{s_{11}^E{\epsilon }_{33}^T}}$

剪切振动：$k_{15}=\frac{d_{15}}{\sqrt{s_{55}^E{\epsilon }_{11}^T}}$

PZT材料的典型耦合系数：$k_{33}\approx 0.7$，$k_{31}\approx 0.35$，$k_{15}\approx 0.68$。

---

3. 压电耦合数学模型

3.1 控制方程

3.1.1 力学平衡方程

$$\nabla \cdot \mathbf{T}+\mathbf{f}=\rho \frac{{\partial }^2\mathbf{u}}{\partial t^2}$$

其中，$\mathbf{T}$为应力张量，$\mathbf{f}$为体积力，$\rho$为密度，$\mathbf{u}$为位移矢量。

3.1.2 静电学方程

在准静态近似下（电磁波长远大于结构尺寸）：

$$\nabla \cdot \mathbf{D}=q_v$$

$$\mathbf{E}=-\nabla \varphi$$

其中，$q_v$为体积电荷密度，$\varphi$为电势。

3.1.3 几何方程

应变-位移关系：

$$S_{ij}=\frac{1}{2}\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)$$

3.2 耦合方程组

综合上述方程，得到压电耦合控制方程组：

$$\{\begin{array}{ll}{c}_{ijkl}^E\frac{{\partial }^2{u}_k}{\partial x_l\partial x_j}+e_{kij}\frac{{\partial }^2\varphi }{\partial x_k\partial x_j}+f_i=\rho \frac{{\partial }^2{u}_i}{\partial t^2}& \text{(力学方程)}\\ e_{ikl}\frac{{\partial }^2{u}_k}{\partial x_l\partial x_i}-{\epsilon }_{ik}^S\frac{{\partial }^2\varphi }{\partial x_k\partial x_i}=q_v& \text{(电学方程)}\end{array}$$

这是一组耦合的偏微分方程，需要同时求解位移场和电势场。

3.3 边界条件

3.3.1 力学边界条件

位移边界：$u_i={\bar{u}}_i$ on ${\Gamma }_u$

应力边界：$T_{ij}{n}_j={\bar{t}}_i$ on ${\Gamma }_t$

3.3.2 电学边界条件

电势边界：$\varphi =\bar{\varphi }$ on ${\Gamma }_{\varphi }$

电位移边界：$D_i{n}_i=-\bar{\sigma }$ on ${\Gamma }_D$

开路条件：${\int }_{{\Gamma }_e}{D}_i{n}_{i}d\Gamma =0$（电极上的总电荷为零）

---

4. 数值求解方法

4.1 有限元离散

4.1.1 弱形式

力学方程的弱形式：

$${\int }_{\Omega }\delta S_{ij}{c}_{ijkl}^E{S}_{kl}d\Omega +{\int }_{\Omega }\delta S_{ij}{e}_{kij}{E}_{k}d\Omega ={\int }_{\Omega }\delta u_i{f}_{i}d\Omega +{\int }_{{\Gamma }_t}\delta u_i{\bar{t}}_{i}d\Gamma -{\int }_{\Omega }\delta u_i\rho {\ddot{u}}_{i}d\Omega$$

电学方程的弱形式：

$${\int }_{\Omega }\delta E_i{e}_{ikl}{S}_{kl}d\Omega -{\int }_{\Omega }\delta E_i{\epsilon }_{ik}^S{E}_{k}d\Omega =-{\int }_{\Omega }\delta \varphi q_{v}d\Omega +{\int }_{{\Gamma }_D}\delta \varphi \bar{\sigma }d\Gamma$$

4.1.2 有限元矩阵方程

离散后得到耦合的矩阵方程：

$$\left[\begin{array}{cc}{\mathbf{M}}_{uu}& 0\\ 0& 0\end{array}\right]\left\{\begin{array}{c}\ddot{\mathbf{u}}\\ \ddot{\mathbf{\varphi }}\end{array}\right\}+\left[\begin{array}{cc}{\mathbf{K}}_{uu}& {\mathbf{K}}_{u\varphi }\\ {\mathbf{K}}_{\varphi u}& {\mathbf{K}}_{\varphi \varphi }\end{array}\right]\left\{\begin{array}{c}\mathbf{u}\\ \mathbf{\varphi }\end{array}\right\}=\left\{\begin{array}{c}{\mathbf{F}}_u\\ {\mathbf{F}}_{\varphi }\end{array}\right\}$$

其中：

• ${\mathbf{K}}_{uu}$：刚度矩阵
• ${\mathbf{K}}_{u\varphi }={\mathbf{K}}_{\varphi u}^T$：压电耦合矩阵
• ${\mathbf{K}}_{\varphi \varphi }$：介电刚度矩阵
• ${\mathbf{M}}_{uu}$：质量矩阵

4.2 谐响应分析

对于简谐激励$\mathbf{F}=\hat{\mathbf{F}}{e}^{j\omega t}$，假设响应形式为$\mathbf{u}=\hat{\mathbf{u}}{e}^{j\omega t}$，$\mathbf{\varphi }=\hat{\mathbf{\varphi }}{e}^{j\omega t}$，得到频域方程：

$$\left[\begin{array}{cc}{\mathbf{K}}_{uu}-{\omega }^2{\mathbf{M}}_{uu}& {\mathbf{K}}_{u\varphi }\\ {\mathbf{K}}_{\varphi u}& {\mathbf{K}}_{\varphi \varphi }\end{array}\right]\left\{\begin{array}{c}\hat{\mathbf{u}}\\ \hat{\mathbf{\varphi }}\end{array}\right\}=\left\{\begin{array}{c}{\hat{\mathbf{F}}}_u\\ {\hat{\mathbf{F}}}_{\varphi }\end{array}\right\}$$

4.3 模态分析

自由振动时，求解特征值问题：

$$\left[\begin{array}{cc}{\mathbf{K}}_{uu}& {\mathbf{K}}_{u\varphi }\\ {\mathbf{K}}_{\varphi u}& {\mathbf{K}}_{\varphi \varphi }\end{array}\right]\left\{\begin{array}{c}{\mathbf{\psi }}_u\\ {\mathbf{\psi }}_{\varphi }\end{array}\right\}={\omega }^2\left[\begin{array}{cc}{\mathbf{M}}_{uu}& 0\\ 0& 0\end{array}\right]\left\{\begin{array}{c}{\mathbf{\psi }}_u\\ {\mathbf{\varphi }}_{\varphi }\end{array}\right\}$$

由于电学方程没有质量项，这是一个特殊的特征值问题，需要采用特殊方法求解。

---

5. Python仿真实现

5.1 案例1：压电双晶片梁弯曲分析

问题描述：分析压电双晶片悬臂梁在电压作用下的弯曲变形。

# -*- coding: utf-8 -*-
"""
主题026：压电耦合分析 - Python仿真代码
包含6个案例：
1. 压电双晶片梁弯曲分析
2. 压电传感器响应分析
3. 压电能量收集器仿真
4. 压电换能器谐响应分析
5. 压电层合板振动模态
6. 压电微泵流固耦合

包含GIF动图展示动态响应
"""
import matplotlib
matplotlib.use('Agg')
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation, PillowWriter
from scipy.linalg import eigh, solve
from scipy.integrate import odeint
import warnings
warnings.filterwarnings('ignore')
import os

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

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

print("="*60)
print("主题026：压电耦合分析")
print("="*60)

# ============================================================
# 案例1：压电双晶片梁弯曲分析
# ============================================================
print("\n" + "="*60)
print("案例1：压电双晶片梁弯曲分析")
print("="*60)

# 几何参数
L = 0.1  # 梁长度 (m)
b = 0.01  # 梁宽度 (m)
h_p = 0.0005  # 单层压电片厚度 (m)
h_sub = 0.001  # 基板厚度 (m)
h_total = 2*h_p + h_sub  # 总厚度

# 材料参数 - PZT-5H
E_p = 60e9  # 压电片弹性模量 (Pa)
d31 = -274e-12  # 压电系数 (C/N)
eps33 = 3400 * 8.854e-12  # 介电常数 (F/m)

# 基板材料（铝）
E_sub = 70e9  # 弹性模量 (Pa)

# 等效弯曲刚度
I_total = b * h_total**3 / 12  # 惯性矩
E_eff = (2*E_p*h_p + E_sub*h_sub) / h_total  # 等效弹性模量
D = E_eff * I_total  # 弯曲刚度

# 等效压电耦合系数
e_eff = d31 * E_p * b * (h_p + h_sub/2)  # 等效压电系数

# 电压-位移关系
V_range = np.linspace(0, 100, 50)  # 电压范围 (V)

# 悬臂梁端部挠度（解析解）
# w_max = (3 * e_eff * V * L^2) / (2 * D)
w_max = (3 * e_eff * V_range * L**2) / (2 * D) * 1000  # mm

print(f"等效弯曲刚度 D: {D:.4f} N·m²")
print(f"等效压电系数 e_eff: {e_eff:.6e} C/m")
print(f"100V电压下最大挠度: {w_max[-1]:.4f} mm")

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

ax1 = axes[0]
ax1.plot(V_range, w_max, 'b-', linewidth=2)
ax1.set_xlabel('Voltage (V)', fontsize=11)
ax1.set_ylabel('Tip Deflection (mm)', fontsize=11)
ax1.set_title('Bimorph Actuator: Voltage vs Deflection', fontsize=12, fontweight='bold')
ax1.grid(True, alpha=0.3)

# 梁变形形状
ax2 = axes[1]
x_beam = np.linspace(0, L, 100)
V_demo = 50  # 演示电压
w_shape = (e_eff * V_demo / (2*D)) * (L**2 * x_beam - x_beam**3/3) * 1000  # mm

ax2.plot(x_beam*1000, w_shape, 'r-', linewidth=2, label='Deformed')
ax2.axhline(y=0, color='k', linestyle='--', alpha=0.5, label='Original')
ax2.fill_between(x_beam*1000, 0, w_shape, alpha=0.3)
ax2.set_xlabel('Position (mm)', fontsize=11)
ax2.set_ylabel('Deflection (mm)', fontsize=11)
ax2.set_title(f'Bimorph Beam Deflection (V={V_demo}V)', fontsize=12, fontweight='bold')
ax2.legend()
ax2.grid(True, alpha=0.3)
ax2.set_aspect('equal', adjustable='box')

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

5.2 案例2：压电传感器响应分析

问题描述：分析压电加速度传感器的频率响应特性。

print("\n" + "="*60)
print("案例2：压电传感器响应分析")
print("="*60)

# 传感器参数
m = 0.001  # 质量块 (kg)
k = 1e6  # 等效刚度 (N/m)
c = 0.1  # 阻尼系数 (N·s/m)
d33 = 593e-12  # 压电系数 (C/N)
C_p = 100e-12  # 压电片电容 (F)

# 固有频率
f_n = np.sqrt(k/m) / (2*np.pi)
print(f"固有频率: {f_n:.2f} Hz")

# 频率响应
f_range = np.logspace(1, 4, 500)
omega = 2 * np.pi * f_range

# 机械传递函数
H_mech = 1 / (k - m*omega**2 + 1j*c*omega)

# 电压灵敏度 (V/(m/s²))
S_v = d33 * k * np.abs(H_mech) / C_p

# 电荷灵敏度 (C/(m/s²))
S_q = d33 * k * np.abs(H_mech)

print(f"低频电压灵敏度: {S_v[0]:.4f} V/(m/s²)")
print(f"峰值灵敏度: {np.max(S_v):.4f} V/(m/s²)")

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

# 幅频特性
ax1 = axes[0, 0]
ax1.semilogx(f_range, S_v*1000, 'b-', linewidth=2)
ax1.axvline(x=f_n, color='r', linestyle='--', label=f'f_n={f_n:.1f}Hz')
ax1.set_xlabel('Frequency (Hz)', fontsize=11)
ax1.set_ylabel('Sensitivity (mV/(m/s²))', fontsize=11)
ax1.set_title('Voltage Sensitivity vs Frequency', fontsize=12, fontweight='bold')
ax1.legend()
ax1.grid(True, alpha=0.3)

# 相频特性
phase = np.angle(H_mech) * 180 / np.pi
ax2 = axes[0, 1]
ax2.semilogx(f_range, phase, 'g-', linewidth=2)
ax2.axvline(x=f_n, color='r', linestyle='--', label=f'f_n={f_n:.1f}Hz')
ax2.set_xlabel('Frequency (Hz)', fontsize=11)
ax2.set_ylabel('Phase (deg)', fontsize=11)
ax2.set_title('Phase Response', fontsize=12, fontweight='bold')
ax2.legend()
ax2.grid(True, alpha=0.3)

# 阻尼比影响
ax3 = axes[1, 0]
zeta_range = [0.01, 0.05, 0.1, 0.3, 0.7]
colors = plt.cm.viridis(np.linspace(0, 1, len(zeta_range)))
for i, zeta in enumerate(zeta_range):
    c_zeta = 2 * zeta * np.sqrt(k*m)
    H_zeta = 1 / (k - m*omega**2 + 1j*c_zeta*omega)
    S_zeta = d33 * k * np.abs(H_zeta) / C_p
    ax3.semilogx(f_range, S_zeta*1000, color=colors[i], linewidth=2, label=f'ζ={zeta}')
ax3.set_xlabel('Frequency (Hz)', fontsize=11)
ax3.set_ylabel('Sensitivity (mV/(m/s²))', fontsize=11)
ax3.set_title('Effect of Damping Ratio', fontsize=12, fontweight='bold')
ax3.legend()
ax3.grid(True, alpha=0.3)

# 时域脉冲响应
ax4 = axes[1, 1]
t_resp = np.linspace(0, 0.1, 1000)
zeta = 0.1
omega_d = f_n * 2 * np.pi * np.sqrt(1 - zeta**2)
h_t = (1/(m*omega_d)) * np.exp(-zeta*2*np.pi*f_n*t_resp) * np.sin(omega_d*t_resp)
ax4.plot(t_resp*1000, h_t*1e6, 'm-', linewidth=2)
ax4.set_xlabel('Time (ms)', fontsize=11)
ax4.set_ylabel('Impulse Response (μm·s)', fontsize=11)
ax4.set_title('Impulse Response (ζ=0.1)', fontsize=12, fontweight='bold')
ax4.grid(True, alpha=0.3)

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

5.3 案例3：压电能量收集器仿真（含GIF动图）

问题描述：分析压电悬臂梁能量收集器在振动环境中的发电性能，并生成GIF动图展示动态响应。

print("\n" + "="*60)
print("案例3：压电能量收集器仿真（含GIF动图）")
print("="*60)

# 能量收集器参数
L_h = 0.05  # 梁长度 (m)
b_h = 0.01  # 梁宽度 (m)
h_pzt = 0.0002  # PZT厚度 (m)
h_shim = 0.0001  # 金属垫片厚度 (m)

# 材料参数
E_pzt = 66e9  # PZT弹性模量 (Pa)
rho_pzt = 7800  # PZT密度 (kg/m³)
d31_h = -190e-12  # 压电系数 (C/N)
eps33_h = 1800 * 8.854e-12  # 介电常数 (F/m)

E_shim = 200e9  # 钢垫片弹性模量
rho_shim = 7850

# 等效参数
h_total_h = h_pzt + h_shim
E_eff_h = (E_pzt*h_pzt + E_shim*h_shim) / h_total_h
rho_eff = (rho_pzt*h_pzt + rho_shim*h_shim) / h_total_h
I_h = b_h * h_total_h**3 / 12

# 等效刚度
k_eff = 3 * E_eff_h * I_h / L_h**3

# 等效质量
m_eff = 0.236 * rho_eff * b_h * h_total_h * L_h

# 固有频率
f_1 = (1.875**2 / (2*np.pi)) * np.sqrt(E_eff_h*I_h / (rho_eff*b_h*h_total_h*L_h**4))
print(f"一阶固有频率: {f_1:.2f} Hz")

# 机电耦合系数
theta = -d31_h * E_pzt * b_h * (h_pzt + h_shim) / (2*L_h)
C_pzt = eps33_h * b_h * L_h / h_pzt
k_em = theta / np.sqrt(k_eff * C_pzt)
print(f"机电耦合系数 k: {k_em:.4f}")

# 振动激励
f_exc = f_1  # 共振激励
omega_exc = 2 * np.pi * f_exc
Y0 = 0.001  # 基础位移幅值 (m)

# 负载电阻扫描
R_load = np.logspace(2, 8, 100)

# 输出功率计算
omega_n = 2 * np.pi * f_1
zeta_m = 0.02  # 机械阻尼比

P_out = np.zeros(len(R_load))
for i, R in enumerate(R_load):
    # 等效阻尼
    zeta_e = (k_em**2 * omega_n * R * C_pzt) / (2 * (1 + (omega_n*R*C_pzt)**2))
    zeta_total = zeta_m + zeta_e
    
    # 位移响应
    r = omega_exc / omega_n
    X = Y0 * r**2 / np.sqrt((1-r**2)**2 + (2*zeta_total*r)**2)
    
    # 输出电压
    V_out = theta * omega_exc * X * R / np.sqrt(1 + (omega_exc*R*C_pzt)**2)
    
    # 输出功率
    P_out[i] = V_out**2 / R

P_max = np.max(P_out)
R_opt = R_load[np.argmax(P_out)]
print(f"最大功率输出: {P_max*1e6:.2f} μW")
print(f"最优负载电阻: {R_opt/1e3:.2f} kΩ")

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

# 功率-电阻曲线
ax1 = axes[0, 0]
ax1.semilogx(R_load/1e3, P_out*1e6, 'b-', linewidth=2)
ax1.axvline(x=R_opt/1e3, color='r', linestyle='--', label=f'R_opt={R_opt/1e3:.1f}kΩ')
ax1.set_xlabel('Load Resistance (kΩ)', fontsize=11)
ax1.set_ylabel('Power (μW)', fontsize=11)
ax1.set_title('Power Output vs Load Resistance', fontsize=12, fontweight='bold')
ax1.legend()
ax1.grid(True, alpha=0.3)

# 频率响应
ax2 = axes[0, 1]
f_sweep = np.linspace(f_1*0.5, f_1*1.5, 200)
omega_sweep = 2 * np.pi * f_sweep
R_fixed = R_opt

P_freq = np.zeros(len(f_sweep))
for i, w in enumerate(omega_sweep):
    r = w / omega_n
    zeta_e = (k_em**2 * omega_n * R_fixed * C_pzt) / (2 * (1 + (omega_n*R_fixed*C_pzt)**2))
    zeta_total = zeta_m + zeta_e
    X = Y0 * r**2 / np.sqrt((1-r**2)**2 + (2*zeta_total*r)**2)
    V_out = theta * w * X * R_fixed / np.sqrt(1 + (w*R_fixed*C_pzt)**2)
    P_freq[i] = V_out**2 / R_fixed

ax2.plot(f_sweep, P_freq*1e6, 'g-', linewidth=2)
ax2.axvline(x=f_1, color='r', linestyle='--', label=f'f_1={f_1:.1f}Hz')
ax2.set_xlabel('Frequency (Hz)', fontsize=11)
ax2.set_ylabel('Power (μW)', fontsize=11)
ax2.set_title('Power vs Frequency (R=R_opt)', fontsize=12, fontweight='bold')
ax2.legend()
ax2.grid(True, alpha=0.3)

# 耦合系数影响
ax3 = axes[1, 0]
k_em_range = np.linspace(0.05, 0.3, 50)
P_max_k = np.zeros(len(k_em_range))
for i, k_em_i in enumerate(k_em_range):
    zeta_e_opt = k_em_i**2 / 4
    zeta_total = zeta_m + zeta_e_opt
    X_opt = Y0 / (2 * zeta_total)
    P_max_k[i] = (theta**2 * omega_n**2 * X_opt**2) / (4 * C_pzt) * (k_em_i / k_em)**2

ax3.plot(k_em_range, P_max_k*1e6, 'm-', linewidth=2)
ax3.axvline(x=k_em, color='r', linestyle='--', label=f'Current k={k_em:.3f}')
ax3.set_xlabel('Electromechanical Coupling k', fontsize=11)
ax3.set_ylabel('Max Power (μW)', fontsize=11)
ax3.set_title('Effect of Coupling Coefficient', fontsize=12, fontweight='bold')
ax3.legend()
ax3.grid(True, alpha=0.3)

# 效率分析
ax4 = axes[1, 1]
zeta_e_range = np.linspace(0, 0.2, 100)
eta = 4 * zeta_e_range * zeta_m / (zeta_m + zeta_e_range)**2
ax4.plot(zeta_e_range, eta*100, 'c-', linewidth=2)
ax4.axvline(x=zeta_m, color='r', linestyle='--', label=f'ζ_m={zeta_m}')
ax4.set_xlabel('Electrical Damping ζ_e', fontsize=11)
ax4.set_ylabel('Efficiency (%)', fontsize=11)
ax4.set_title('Energy Conversion Efficiency', fontsize=12, fontweight='bold')
ax4.legend()
ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig(f'{output_dir}/case3_energy_harvester.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 能量收集器静态分析完成")

# 生成GIF动图 - 压电梁振动响应
print("\n正在生成GIF动图...")

fig, ax = plt.subplots(figsize=(10, 6))

# 时间和空间网格
x_anim = np.linspace(0, L_h, 50)
t_anim = np.linspace(0, 2/f_1, 100)  # 2个周期

# 计算振动响应（基频模态）
def beam_mode_shape(x, L):
    """悬臂梁一阶模态形状"""
    beta = 1.875 / L
    return np.cosh(beta*x) - np.cos(beta*x) - 0.734*(np.sinh(beta*x) - np.sin(beta*x))

phi_1 = beam_mode_shape(x_anim, L_h)
phi_1 = phi_1 / np.max(phi_1)  # 归一化

# 时间响应（包含机电耦合的衰减）
w_n = 2 * np.pi * f_1
zeta_total = zeta_m + k_em**2/4  # 包含机电阻尼
A0 = 0.005  # 初始振幅 (m)

line, = ax.plot([], [], 'b-', linewidth=3, label='Beam')
base_line, = ax.plot([], [], 'k-', linewidth=5, label='Base')
voltage_text = ax.text(0.02, 0.95, '', transform=ax.transAxes, fontsize=12)

ax.set_xlim(-0.01, L_h*1.2)
ax.set_ylim(-0.015, 0.015)
ax.set_xlabel('Position (m)', fontsize=11)
ax.set_ylabel('Displacement (m)', fontsize=11)
ax.set_title('Piezoelectric Energy Harvester Vibration', fontsize=12, fontweight='bold')
ax.legend(loc='lower right')
ax.grid(True, alpha=0.3)
ax.axhline(y=0, color='k', linewidth=0.5)

def init():
    line.set_data([], [])
    base_line.set_data([-0.005, 0.005], [0, 0])
    voltage_text.set_text('')
    return line, base_line, voltage_text

def update(frame):
    t = t_anim[frame]
    # 衰减振动
    decay = np.exp(-zeta_total * w_n * t)
    w_t = A0 * decay * np.cos(w_n * t)
    
    # 梁变形
    y_beam = w_t * phi_1
    
    # 基础激励
    y_base = Y0 * np.cos(2*np.pi*f_exc*t)
    
    line.set_data(x_anim, y_beam + y_base)
    base_line.set_data([-0.005, 0.005], [y_base, y_base])
    
    # 计算输出电压
    V_out = theta * w_n * w_t / np.sqrt(1 + (w_n*R_opt*C_pzt)**2)
    voltage_text.set_text(f'V_out = {V_out*1000:.2f} mV\nTime = {t*1000:.1f} ms')
    
    return line, base_line, voltage_text

anim = FuncAnimation(fig, update, frames=len(t_anim), init_func=init, 
                     blit=True, interval=50, repeat=True)

# 保存GIF
gif_path = f'{output_dir}/case3_harvester_vibration.gif'
anim.save(gif_path, writer=PillowWriter(fps=20))
plt.close()
print(f"✓ GIF动图已保存: {gif_path}")

5.4 案例4：压电换能器谐响应分析

问题描述：分析压电超声换能器的谐响应特性。

print("\n" + "="*60)
print("案例4：压电换能器谐响应分析")
print("="*60)

# 换能器参数（Langevin型）
# 由前后质量块、压电陶瓷堆和预应力螺栓组成

# 前质量块（铝）
L_front = 0.02  # m
A_front = 0.0005  # m²
E_front = 70e9
rho_front = 2700

# 压电陶瓷堆（PZT-8）
n_pzt = 2  # 陶瓷片数
L_pzt_stack = n_pzt * 0.005  # m
A_pzt = 0.0005  # m²
E_pzt_8 = 100e9
rho_pzt_8 = 7600
d33_8 = 225e-12
k33_8 = 0.64

# 后质量块（钢）
L_back = 0.03  # m
A_back = 0.0005  # m²
E_back = 210e9
rho_back = 7850

# 计算等效参数
# 纵向振动等效刚度
k_front = E_front * A_front / L_front
k_pzt = E_pzt_8 * A_pzt / L_pzt_stack
k_back = E_back * A_back / L_back

# 串联刚度
k_total = 1 / (1/k_front + 1/k_pzt + 1/k_back)

# 等效质量
m_front = rho_front * A_front * L_front
m_pzt = rho_pzt_8 * A_pzt * L_pzt_stack
m_back = rho_back * A_back * L_back
m_total = m_front + m_pzt + m_back

# 谐振频率
f_res = np.sqrt(k_total / m_total) / (2*np.pi)
print(f"谐振频率: {f_res:.2f} Hz")

# 反谐振频率
C_0 = eps33_h * A_pzt / L_pzt_stack
f_anti = f_res / np.sqrt(1 - k33_8**2)
print(f"反谐振频率: {f_anti:.2f} Hz")

# 导纳圆分析
f_admit = np.linspace(f_res*0.9, f_anti*1.1, 500)
omega_admit = 2 * np.pi * f_admit

# 机械品质因数
Q_m = 500

# 导纳计算
Y = np.zeros(len(f_admit), dtype=complex)
for i, w in enumerate(omega_admit):
    Z_mech = 1j*w*m_total + k_total/(1j*w) + w*m_total/Q_m
    Z_elec = 1/(1j*w*C_0)
    n = d33_8 * A_pzt / L_pzt_stack
    Z_coupled = Z_elec + n**2 / Z_mech
    Y[i] = 1 / Z_coupled

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

# 导纳幅频特性
ax1 = axes[0, 0]
ax1.semilogy(f_admit/1e3, np.abs(Y)*1e3, 'b-', linewidth=2)
ax1.axvline(x=f_res/1e3, color='r', linestyle='--', label=f'f_r={f_res/1e3:.2f}kHz')
ax1.axvline(x=f_anti/1e3, color='g', linestyle='--', label=f'f_a={f_anti/1e3:.2f}kHz')
ax1.set_xlabel('Frequency (kHz)', fontsize=11)
ax1.set_ylabel('|Y| (mS)', fontsize=11)
ax1.set_title('Admittance Magnitude', fontsize=12, fontweight='bold')
ax1.legend()
ax1.grid(True, alpha=0.3)

# 导纳圆（Nyquist图）
ax2 = axes[0, 1]
ax2.plot(Y.real*1e3, Y.imag*1e3, 'r-', linewidth=2)
ax2.set_xlabel('Re(Y) (mS)', fontsize=11)
ax2.set_ylabel('Im(Y) (mS)', fontsize=11)
ax2.set_title('Admittance Circle (Nyquist Plot)', fontsize=12, fontweight='bold')
ax2.grid(True, alpha=0.3)
ax2.set_aspect('equal')

# 阻抗特性
Z = 1 / Y
ax3 = axes[1, 0]
ax3.semilogy(f_admit/1e3, np.abs(Z), 'g-', linewidth=2)
ax3.axvline(x=f_res/1e3, color='r', linestyle='--', label=f'f_r={f_res/1e3:.2f}kHz')
ax3.axvline(x=f_anti/1e3, color='g', linestyle='--', label=f'f_a={f_anti/1e3:.2f}kHz')
ax3.set_xlabel('Frequency (kHz)', fontsize=11)
ax3.set_ylabel('|Z| (Ω)', fontsize=11)
ax3.set_title('Impedance Magnitude', fontsize=12, fontweight='bold')
ax3.legend()
ax3.grid(True, alpha=0.3)

# 相位特性
ax4 = axes[1, 1]
phase_Z = np.angle(Z) * 180 / np.pi
ax4.plot(f_admit/1e3, phase_Z, 'm-', linewidth=2)
ax4.axvline(x=f_res/1e3, color='r', linestyle='--', alpha=0.5)
ax4.axvline(x=f_anti/1e3, color='g', linestyle='--', alpha=0.5)
ax4.axhline(y=0, color='k', linewidth=0.5)
ax4.set_xlabel('Frequency (kHz)', fontsize=11)
ax4.set_ylabel('Phase (deg)', fontsize=11)
ax4.set_title('Impedance Phase', fontsize=12, fontweight='bold')
ax4.grid(True, alpha=0.3)

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

5.5 案例5：压电层合板振动模态（含GIF动图）

问题描述：分析压电层合板的振动模态，并生成GIF动图展示模态振型。

print("\n" + "="*60)
print("案例5：压电层合板振动模态（含GIF动图）")
print("="*60)

# 层合板参数
a = 0.2  # 板长 (m)
b = 0.15  # 板宽 (m)
h_base = 0.002  # 基板厚度 (m)
h_piezo = 0.0005  # 压电层厚度 (m)

# 材料参数
E_base = 70e9  # 铝基板
rho_base = 2700
nu_base = 0.33

E_pzt_plate = 60e9
rho_pzt = 7800
nu_pzt = 0.31
d31_plate = -274e-12

# 等效弯曲刚度（简化的层合板理论）
D_base = E_base * h_base**3 / (12 * (1 - nu_base**2))
D_piezo = E_pzt_plate * ((h_base + 2*h_piezo)**3 - h_base**3) / (12 * (1 - nu_pzt**2))
D_eff = D_base + D_piezo

# 等效密度
rho_eff = (rho_base * h_base + 2 * rho_pzt * h_piezo) / (h_base + 2*h_piezo)

# 简支矩形板固有频率
# f_mn = (π/2) * sqrt(D_eff/(rho_eff*h_total)) * (m²/a² + n²/b²)
h_total = h_base + 2*h_piezo

modes = []
freqs = []
for m in range(1, 4):
    for n in range(1, 4):
        f_mn = (np.pi/2) * np.sqrt(D_eff / (rho_eff * h_total)) * np.sqrt((m/a)**2 + (n/b)**2)
        modes.append((m, n))
        freqs.append(f_mn)

# 排序
sorted_indices = np.argsort(freqs)
modes = [modes[i] for i in sorted_indices]
freqs = [freqs[i] for i in sorted_indices]

print("前9阶固有频率:")
for i, (mode, freq) in enumerate(zip(modes[:9], freqs[:9])):
    print(f"  Mode {i+1} ({mode[0]},{mode[1]}): {freq:.2f} Hz")

# 可视化模态
fig, axes = plt.subplots(3, 3, figsize=(14, 12))
axes = axes.flatten()

x_plate = np.linspace(0, a, 50)
y_plate = np.linspace(0, b, 40)
X_p, Y_p = np.meshgrid(x_plate, y_plate)

for idx in range(9):
    m, n = modes[idx]
    f = freqs[idx]
    
    # 模态振型
    Z = np.sin(m * np.pi * X_p / a) * np.sin(n * np.pi * Y_p / b)
    
    ax = axes[idx]
    im = ax.contourf(X_p*1000, Y_p*1000, Z, levels=20, cmap='RdBu_r')
    ax.set_title(f'Mode {idx+1} ({m},{n}): {f:.1f}Hz', fontsize=10, fontweight='bold')
    ax.set_xlabel('x (mm)', fontsize=9)
    ax.set_ylabel('y (mm)', fontsize=9)
    ax.set_aspect('equal')
    plt.colorbar(im, ax=ax, shrink=0.6)

plt.tight_layout()
plt.savefig(f'{output_dir}/case5_plate_modes.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 模态分析静态图完成")

# 生成GIF动图 - 第一阶模态振动
print("\n正在生成模态振动GIF动图...")

fig, ax = plt.subplots(figsize=(10, 8))

# 选择第一阶模态
m, n = modes[0]
f_mode = freqs[0]
Z_mode = np.sin(m * np.pi * X_p / a) * np.sin(n * np.pi * Y_p / b)

# 时间参数
t_mode = np.linspace(0, 2/f_mode, 60)  # 2个周期

im = ax.contourf(X_p*1000, Y_p*1000, Z_mode, levels=20, cmap='RdBu_r', vmin=-1, vmax=1)
ax.set_xlabel('x (mm)', fontsize=11)
ax.set_ylabel('y (mm)', fontsize=11)
ax.set_title(f'Plate Mode ({m},{n}) Vibration - f={f_mode:.1f}Hz', fontsize=12, fontweight='bold')
ax.set_aspect('equal')
cbar = plt.colorbar(im, ax=ax)
cbar.set_label('Normalized Displacement', fontsize=10)

def update_mode(frame):
    t = t_mode[frame]
    Z_t = Z_mode * np.cos(2 * np.pi * f_mode * t)
    ax.clear()
    im = ax.contourf(X_p*1000, Y_p*1000, Z_t, levels=20, cmap='RdBu_r', vmin=-1, vmax=1)
    ax.set_xlabel('x (mm)', fontsize=11)
    ax.set_ylabel('y (mm)', fontsize=11)
    ax.set_title(f'Plate Mode ({m},{n}) - t={t*1000:.1f}ms', fontsize=12, fontweight='bold')
    ax.set_aspect('equal')
    return im,

anim_mode = FuncAnimation(fig, update_mode, frames=len(t_mode), 
                          interval=50, blit=False, repeat=True)

gif_mode_path = f'{output_dir}/case5_plate_mode_vibration.gif'
anim_mode.save(gif_mode_path, writer=PillowWriter(fps=20))
plt.close()
print(f"✓ 模态振动GIF已保存: {gif_mode_path}")

5.6 案例6：压电微泵流固耦合（含GIF动图）

问题描述：分析压电驱动微泵的工作原理，并生成GIF动图展示泵送过程。

print("\n" + "="*60)
print("案例6：压电微泵流固耦合（含GIF动图）")
print("="*60)

# 微泵参数
L_pump = 0.02  # 泵腔长度 (m)
W_pump = 0.01  # 泵腔宽度 (m)
H_cavity = 0.0005  # 泵腔高度 (m)
h_diaphragm = 0.0002  # 压电薄膜厚度 (m)

# 材料参数
E_diaphragm = 60e9  # PZT
rho_diaphragm = 7800
d31_pump = -190e-12

# 流体参数（水）
rho_water = 1000
mu_water = 0.001

# 驱动参数
f_drive = 1000  # 驱动频率 (Hz)
V_drive = 50  # 驱动电压 (V)

# 薄膜等效刚度
D_diaphragm = E_diaphragm * h_diaphragm**3 / (12 * (1 - 0.31**2))
k_diaphragm = 64 * D_diaphragm / L_pump**2  # 简支板近似

# 压电驱动力
F_piezo = d31_pump * E_diaphragm * V_drive * W_pump

# 薄膜位移幅值
w_diaphragm = F_piezo / k_diaphragm
print(f"薄膜位移幅值: {w_diaphragm*1e6:.2f} μm")

# 泵腔体积变化
Delta_V = w_diaphragm * L_pump * W_pump / 2
print(f"体积变化: {Delta_V*1e9:.4f} nL")

# 流量估算
Q_peak = 2 * np.pi * f_drive * Delta_V
print(f"峰值流量: {Q_peak*1e6:.2f} μL/s")

# 压力生成
# 考虑流体惯性
omega_pump = 2 * np.pi * f_drive
fluid_mass = rho_water * L_pump * W_pump * H_cavity
fluid_stiffness = rho_water * 340**2 * W_pump * H_cavity / L_pump  # 声顺

# 简化模型：薄膜-流体耦合
m_eff_pump = rho_diaphragm * h_diaphragm * L_pump * W_pump + fluid_mass
k_eff_pump = k_diaphragm + fluid_stiffness

f_pump_res = np.sqrt(k_eff_pump / m_eff_pump) / (2*np.pi)
print(f"泵共振频率: {f_pump_res:.2f} Hz")

# 频率响应分析
f_pump_range = np.linspace(100, 5000, 200)
omega_pump_range = 2 * np.pi * f_pump_range

Q_pump = np.zeros(len(f_pump_range))
P_pump = np.zeros(len(f_pump_range))

zeta_pump = 0.1  # 阻尼比

for i, w in enumerate(omega_pump_range):
    r = w / (2*np.pi*f_pump_res)
    H_response = 1 / np.sqrt((1-r**2)**2 + (2*zeta_pump*r)**2)
    w_disp = w_diaphragm * H_response
    Q_pump[i] = w * w_disp * L_pump * W_pump / 2
    P_pump[i] = Q_pump[i] * 8 * mu_water * L_pump / (np.pi * (H_cavity/2)**4)

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

# 流量-频率特性
ax1 = axes[0, 0]
ax1.plot(f_pump_range, Q_pump*1e6, 'b-', linewidth=2)
ax1.axvline(x=f_pump_res, color='r', linestyle='--', label=f'f_res={f_pump_res:.1f}Hz')
ax1.set_xlabel('Frequency (Hz)', fontsize=11)
ax1.set_ylabel('Flow Rate (μL/s)', fontsize=11)
ax1.set_title('Pump Flow Rate vs Frequency', fontsize=12, fontweight='bold')
ax1.legend()
ax1.grid(True, alpha=0.3)

# 压力-频率特性
ax2 = axes[0, 1]
ax2.plot(f_pump_range, P_pump/1e3, 'g-', linewidth=2)
ax2.axvline(x=f_pump_res, color='r', linestyle='--', label=f'f_res={f_pump_res:.1f}Hz')
ax2.set_xlabel('Frequency (Hz)', fontsize=11)
ax2.set_ylabel('Pressure (kPa)', fontsize=11)
ax2.set_title('Pump Pressure vs Frequency', fontsize=12, fontweight='bold')
ax2.legend()
ax2.grid(True, alpha=0.3)

# 电压-位移特性
ax3 = axes[1, 0]
V_range_pump = np.linspace(0, 100, 50)
w_range = d31_pump * E_diaphragm * V_range_pump * W_pump / k_diaphragm * 1e6
ax3.plot(V_range_pump, w_range, 'm-', linewidth=2)
ax3.set_xlabel('Voltage (V)', fontsize=11)
ax3.set_ylabel('Displacement (μm)', fontsize=11)
ax3.set_title('Diaphragm Displacement vs Voltage', fontsize=12, fontweight='bold')
ax3.grid(True, alpha=0.3)

# 泵效率分析
ax4 = axes[1, 1]
# 机械功率
P_mech = 0.5 * k_diaphragm * (w_diaphragm * np.ones(len(f_pump_range)))**2 * f_pump_range
# 流体功率（简化）
P_fluid = Q_pump * P_pump
eta_pump = np.abs(P_fluid / (P_mech + 1e-10)) * 100
ax4.plot(f_pump_range, eta_pump, 'c-', linewidth=2)
ax4.axvline(x=f_pump_res, color='r', linestyle='--', label=f'f_res={f_pump_res:.1f}Hz')
ax4.set_xlabel('Frequency (Hz)', fontsize=11)
ax4.set_ylabel('Efficiency (%)', fontsize=11)
ax4.set_title('Pump Efficiency vs Frequency', fontsize=12, fontweight='bold')
ax4.legend()
ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig(f'{output_dir}/case6_micropump_analysis.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 微泵静态分析完成")

# 生成GIF动图 - 微泵工作过程
print("\n正在生成微泵工作GIF动图...")

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6))

# 时间参数
t_pump = np.linspace(0, 2/f_drive, 80)  # 2个周期

# 空间网格
x_pump = np.linspace(0, L_pump, 30)
y_pump = np.linspace(0, H_cavity, 20)
X_pump, Y_pump = np.meshgrid(x_pump, y_pump)

def update_pump(frame):
    t = t_pump[frame]
    
    # 薄膜位移
    w_t = w_diaphragm * np.sin(2 * np.pi * f_drive * t)
    
    # 薄膜形状（简支梁近似）
    diaphragm_shape = w_t * np.sin(np.pi * x_pump / L_pump)
    
    # 速度场（简化模型）
    u_vel = np.zeros_like(X_pump)
    v_vel = np.zeros_like(Y_pump)
    
    for i in range(len(y_pump)):
        for j in range(len(x_pump)):
            # 简化的速度分布
            u_vel[i, j] = -2 * np.pi * f_drive * w_t * np.cos(np.pi * x_pump[j] / L_pump) * (y_pump[i] / H_cavity)
            v_vel[i, j] = np.pi * f_drive * w_t * np.sin(np.pi * x_pump[j] / L_pump) * (1 - y_pump[i] / H_cavity)
    
    # 左图：薄膜变形
    ax1.clear()
    ax1.fill_between(x_pump*1000, 0, diaphragm_shape*1e6, alpha=0.5, color='blue', label='Diaphragm')
    ax1.axhline(y=0, color='k', linewidth=1)
    ax1.axhline(y=-H_cavity*1e6, color='brown', linewidth=2, label='Cavity bottom')
    ax1.set_xlim(0, L_pump*1000)
    ax1.set_ylim(-H_cavity*1.5*1e6, w_diaphragm*1.5*1e6)
    ax1.set_xlabel('x (mm)', fontsize=11)
    ax1.set_ylabel('y (μm)', fontsize=11)
    ax1.set_title(f'Diaphragm Deflection - t={t*1000:.1f}ms', fontsize=12, fontweight='bold')
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    
    # 右图：速度场
    ax2.clear()
    speed = np.sqrt(u_vel**2 + v_vel**2)
    im = ax2.contourf(X_pump*1000, Y_pump*1e6, speed*1e3, levels=20, cmap='jet')
    ax2.quiver(X_pump[::2, ::3]*1000, Y_pump[::2, ::3]*1e6, 
               u_vel[::2, ::3]*1e3, v_vel[::2, ::3]*1e3, 
               scale=50, color='white', alpha=0.7)
    ax2.set_xlabel('x (mm)', fontsize=11)
    ax2.set_ylabel('y (μm)', fontsize=11)
    ax2.set_title(f'Flow Field - t={t*1000:.1f}ms', fontsize=12, fontweight='bold')
    ax2.set_aspect('equal')
    
    return im,

anim_pump = FuncAnimation(fig, update_pump, frames=len(t_pump), 
                          interval=40, blit=False, repeat=True)

gif_pump_path = f'{output_dir}/case6_micropump_operation.gif'
anim_pump.save(gif_pump_path, writer=PillowWriter(fps=25))
plt.close()
print(f"✓ 微泵工作GIF已保存: {gif_pump_path}")

print("\n" + "="*60)
print("所有案例仿真完成！")
print("="*60)
print("\n生成的图片文件:")
print("  - case1_bimorph_beam.png")
print("  - case2_sensor_response.png")
print("  - case3_energy_harvester.png")
print("  - case3_harvester_vibration.gif (GIF动图)")
print("  - case4_transducer_harmonic.png")
print("  - case5_plate_modes.png")
print("  - case5_plate_mode_vibration.gif (GIF动图)")
print("  - case6_micropump_analysis.png")
print("  - case6_micropump_operation.gif (GIF动图)")
print("="*60)