限元方法进行电磁-热耦合模拟

1. 使用有限元方法进行电磁-热耦合模拟。

2. 求解麦克斯韦方程组得到电磁场分布，然后耦合傅里叶热传导方程计算温度场。

3. 模型考虑了木材的介电特性（随温度、含水率等变化，但论文中假设为恒定以简化）。

4. 边界条件：微波腔体的金属壁为理想导体，采用阻抗边界条件；木材表面与空气的对流换热。

5. 网格划分需要满足电磁波计算的精度（通常每个波长内至少10个网格）。

import numpy as np
import matplotlib.pyplot as plt
import scipy.sparse as sp
import scipy.sparse.linalg as spla
from scipy.constants import pi, c, epsilon_0, mu_0
import pandas as pd
import time
from typing import Dict, Tuple, List, Optional
import logging
from dataclasses import dataclass
from enum import Enum
import sys
import os# 配置中文字体支持
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']  # 使用SimHei字体，如果不存在则使用DejaVu Sans
plt.rcParams['axes.unicode_minus'] = False  # 正常显示负号# 配置日志
logging.basicConfig(level=logging.INFO, format='%(asctime)s - %(levelname)s - %(message)s')
logger = logging.getLogger(__name__)class WoodSpecies(Enum):"""木材种类枚举"""PAULOWNIA = "paulownia"  # 泡桐，论文中使用的材料class MicrowaveMode(Enum):"""微波模式枚举"""TE10 = "TE10"  # 横电模式，论文中使用的模式@dataclass
class WoodProperties:"""木材材料属性（论文表2）"""density: float = 237.0  # kg/m³specific_heat: float = 1200.0  # J/(kg·K)thermal_conductivity: float = 0.105  # W/(m·K)epsilon_prime: float = 2.0  # 介电常数实部epsilon_double_prime: float = 0.1  # 介电常数虚部（损耗因子）electrical_conductivity: float = 0.01  # 电导率 S/m@dataclass
class MicrowaveCavity:"""微波腔体参数（论文表1）"""length: float = 0.5  # 腔体长度(m)width: float = 0.4  # 腔体宽度(m)height: float = 0.3  # 腔体高度(m)frequency: float = 2.45e9  # 频率2.45GHzmode: MicrowaveMode = MicrowaveMode.TE10@dataclass
class WoodSample:"""木材样品参数（论文2.1节）"""tangential: float = 0.1  # 弦向100mmradial: float = 0.05  # 径向50mmlongitudinal: float = 0.08  # 纵向80mmspecies: WoodSpecies = WoodSpecies.PAULOWNIAinitial_density: float = 237.0  # kg/m³class MicrowaveWoodHeatingSimulator:"""微波加热木材电磁-热耦合模拟器基于论文《微波加热中木材内部温度场的电磁热耦合模拟》实现"""def __init__(self, wood_properties: WoodProperties,cavity: MicrowaveCavity,wood_sample: WoodSample,grid_resolution: float = 0.005):"""初始化模拟器"""self.wood_properties = wood_propertiesself.cavity = cavityself.wood_sample = wood_sampleself.grid_resolution = grid_resolution# 边界条件参数（论文1.4节）self.environment_temperature = 20.0  # °Cself.convection_coefficient = 10.0  # W/(m²·K)self.initial_temperature = 20.0  # °Cself.setup_parameters()self.initialize_mesh()self.setup_matrices()logger.info("Microwave wood heating simulator initialized successfully")def setup_parameters(self):"""设置模拟参数"""logger.info("Setting up simulation parameters...")# 计算网格维度（论文1.6节）self.nx = int(self.wood_sample.tangential / self.grid_resolution) + 1self.ny = int(self.wood_sample.radial / self.grid_resolution) + 1self.nz = int(self.wood_sample.longitudinal / self.grid_resolution) + 1logger.info(f"Grid dimensions: {self.nx} × {self.ny} × {self.nz}")logger.info(f"Total grid points: {self.nx * self.ny * self.nz}")# 计算时间步长（CFL条件）alpha = self.wood_properties.thermal_conductivity / \(self.wood_properties.density * self.wood_properties.specific_heat)dx = self.grid_resolutiondt_max = 0.25 / (alpha * (1 / dx ** 2 + 1 / dx ** 2 + 1 / dx ** 2))self.dt = min(dt_max, 0.1)  # 最大0.1秒logger.info(f"Heat conduction time step: {self.dt:.4f} s")def initialize_mesh(self):"""初始化网格系统"""logger.info("Initializing mesh system...")# 创建网格坐标（论文1.5节）self.x = np.linspace(0, self.wood_sample.tangential, self.nx)self.y = np.linspace(0, self.wood_sample.radial, self.ny)self.z = np.linspace(0, self.wood_sample.longitudinal, self.nz)# 创建网格点矩阵self.X, self.Y, self.Z = np.meshgrid(self.x, self.y, self.z, indexing='ij')# 初始化温度场self.temperature = np.ones((self.nx, self.ny, self.nz)) * self.initial_temperaturelogger.info("Mesh system initialized successfully")def setup_matrices(self):"""设置有限差分矩阵"""self.setup_thermal_matrices()self.setup_electromagnetic_matrices()def setup_thermal_matrices(self):"""设置热传导矩阵（论文公式11）"""logger.info("Setting up thermal matrices...")dx = self.x[1] - self.x[0]dy = self.y[1] - self.y[0]dz = self.z[1] - self.z[0]# 构建三维拉普拉斯算子self.laplacian_matrix = self._build_3d_laplacian(dx, dy, dz)self.laplacian_matrix = sp.csr_matrix(self.laplacian_matrix)logger.info("Thermal matrices setup completed")def setup_electromagnetic_matrices(self):"""设置电磁场矩阵（论文第1.2节）"""logger.info("Setting up electromagnetic matrices...")self.propagation_constant = self._calculate_propagation_constant()logger.info(f"Propagation constant: {self.propagation_constant:.4f} rad/m")def _build_3d_laplacian(self, dx: float, dy: float, dz: float) -> sp.csr_matrix:"""构建三维拉普拉斯算子的稀疏矩阵"""n_total = self.nx * self.ny * self.nzdata, row_ind, col_ind = [], [], []dx2, dy2, dz2 = dx ** 2, dy ** 2, dz ** 2for i in range(self.nx):for j in range(self.ny):for k in range(self.nz):idx = self._get_1d_index(i, j, k)# 中心点系数coeff_center = -2 / dx2 - 2 / dy2 - 2 / dz2data.append(coeff_center)row_ind.append(idx)col_ind.append(idx)# 邻居点for di, dj, dk in [(-1, 0, 0), (1, 0, 0), (0, -1, 0), (0, 1, 0), (0, 0, -1), (0, 0, 1)]:ni, nj, nk = i + di, j + dj, k + dkif 0 <= ni < self.nx and 0 <= nj < self.ny and 0 <= nk < self.nz:coeff = 1 / dx2 if di != 0 else (1 / dy2 if dj != 0 else 1 / dz2)data.append(coeff)row_ind.append(idx)col_ind.append(self._get_1d_index(ni, nj, nk))return sp.csr_matrix((data, (row_ind, col_ind)), shape=(n_total, n_total))def _get_1d_index(self, i: int, j: int, k: int) -> int:"""将三维索引转换为一维索引"""return i * (self.ny * self.nz) + j * self.nz + kdef _calculate_propagation_constant(self) -> complex:"""计算传播常数（论文公式14-15）"""a, b = 0.1, 0.05  # 波导尺寸m, n = 1, 0  # TE10模式v_c = (c / 2) * np.sqrt((m / a) ** 2 + (n / b) ** 2)v = self.cavity.frequencyif v > v_c:return 2 * pi * np.sqrt((v / c) ** 2 - (v_c / c) ** 2)return 0  # 低于截止频率def calculate_electric_field_distribution(self, microwave_power: float) -> np.ndarray:"""计算电场分布（论文第3.1节）"""logger.info(f"Calculating electric field distribution, power: {microwave_power} kW")power_scaling = np.sqrt(microwave_power / 3.6)E_field = np.zeros((self.nx, self.ny, self.nz), dtype=complex)for i in range(self.nx):for j in range(self.ny):for k in range(self.nz):x_norm, y_norm, z_norm = i / (self.nx - 1), j / (self.ny - 1), k / (self.nz - 1)# 三个波导干涉叠加E1 = np.sin(2 * pi * x_norm) * np.cos(2 * pi * y_norm) * np.exp(1j * 2 * pi * z_norm)E2 = np.cos(2 * pi * x_norm) * np.sin(2 * pi * y_norm) * np.exp(1j * (2 * pi * z_norm + pi / 2))E3 = np.sin(2 * pi * x_norm + pi / 4) * np.cos(2 * pi * y_norm + pi / 4) * np.exp(1j * (2 * pi * z_norm + pi))E_field[i, j, k] = (E1 + E2 + E3) / 3.0# 功率缩放E_max_base = 8.73e3E_max = E_max_base * power_scalingcurrent_max = np.max(np.abs(E_field))if current_max > 0:E_field = (E_field / current_max) * E_maxlogger.info(f"Max electric field: {np.max(np.abs(E_field)):.2f} V/m")return E_fielddef calculate_power_density(self, E_field: np.ndarray) -> np.ndarray:"""计算功率密度（论文公式10）"""omega = 2 * pi * self.cavity.frequencyepsilon_double_prime = self.wood_properties.epsilon_double_primepower_density = 0.5 * omega * epsilon_0 * epsilon_double_prime * np.abs(E_field) ** 2logger.info(f"Max power density: {np.max(power_density):.2f} W/m³")return power_densitydef solve_heat_conduction(self, temperature: np.ndarray,power_density: np.ndarray,time_step: float) -> np.ndarray:"""求解热传导方程（论文公式11）"""alpha = self.wood_properties.thermal_conductivity / \(self.wood_properties.density * self.wood_properties.specific_heat)T_flat = temperature.flatten()Q_flat = power_density.flatten()# 应用边界条件T_flat = self.apply_boundary_conditions(T_flat)n = self.nx * self.ny * self.nzI = sp.eye(n)A = I - alpha * time_step * self.laplacian_matrixsource_term = (time_step / (self.wood_properties.density *self.wood_properties.specific_heat)) * Q_flatrhs = T_flat + source_termtry:T_new_flat = spla.spsolve(A, rhs)except:T_new_flat, info = spla.gmres(A, rhs, maxiter=1000, tol=1e-6)if info != 0:logger.warning("Iterative solver failed, using original values")T_new_flat = T_flatreturn T_new_flat.reshape((self.nx, self.ny, self.nz))def apply_boundary_conditions(self, temperature_flat: np.ndarray) -> np.ndarray:"""应用边界条件（论文1.4节）"""temperature_new = temperature_flat.copy()for i in range(self.nx):for j in range(self.ny):for k in range(self.nz):if (i == 0 or i == self.nx - 1 or j == 0 or j == self.ny - 1 or k == 0 or k == self.nz - 1):idx = self._get_1d_index(i, j, k)T_current = temperature_flat[idx]dT = T_current - self.environment_temperaturetemperature_new[idx] = T_current - 0.05 * dTreturn temperature_newdef simulate_heating_process(self, microwave_power: float,total_time: float,save_interval: float = 1.0) -> Dict:"""执行微波加热过程模拟"""logger.info(f"Starting simulation: {microwave_power} kW, {total_time} s")start_time = time.time()n_steps = int(total_time / self.dt)save_every = max(1, int(save_interval / self.dt))# 预计算电场和功率密度E_field = self.calculate_electric_field_distribution(microwave_power)power_density = self.calculate_power_density(E_field)temperature = np.ones((self.nx, self.ny, self.nz)) * self.initial_temperaturetime_points, temperature_history = [], []for step in range(n_steps + 1):current_time = step * self.dtif step % save_every == 0:time_points.append(current_time)temperature_history.append(temperature.copy())if step % (10 * save_every) == 0:max_temp = np.max(temperature)avg_temp = np.mean(temperature)logger.info(f"Time: {current_time:.1f}s, Max temp: {max_temp:.1f}°C, Avg temp: {avg_temp:.1f}°C")if step < n_steps:temperature = self.solve_heat_conduction(temperature, power_density, self.dt)temperature_history = np.array(temperature_history)metrics = self.calculate_simulation_metrics(temperature_history, np.array(time_points))results = {'microwave_power': microwave_power,'total_time': total_time,'time_points': time_points,'temperature_history': temperature_history,'electric_field': E_field,'power_density': power_density,'metrics': metrics,'simulation_time': time.time() - start_time}logger.info(f"Simulation completed in {results['simulation_time']:.2f} s")return resultsdef calculate_simulation_metrics(self, temperature_history: np.ndarray,time_points: np.ndarray) -> Dict:"""计算模拟性能指标"""final_temperature = temperature_history[-1]surface_temp = final_temperature[0, self.ny // 2, self.nz // 2]core_temp = final_temperature[self.nx // 2, self.ny // 2, self.nz // 2]sub_surface_temp = final_temperature[1, self.ny // 2, self.nz // 2]metrics = {'surface_temperature': surface_temp,'core_temperature': core_temp,'sub_surface_temperature': sub_surface_temp,'max_temperature': np.max(final_temperature),'min_temperature': np.min(final_temperature),'temperature_difference': core_temp - surface_temp,'temperature_std': np.std(final_temperature),'heating_rate': (core_temp - self.initial_temperature) / time_points[-1]}logger.info("=== Simulation Results ===")logger.info(f"Surface temperature: {surface_temp:.1f}°C")logger.info(f"Core temperature: {core_temp:.1f}°C")logger.info(f"Temperature difference: {metrics['temperature_difference']:.1f}°C")return metricsclass ResultVisualizer:"""结果可视化类"""@staticmethoddef plot_electric_field_distribution(electric_field: np.ndarray,microwave_power: float,wood_sample: WoodSample):"""绘制电场分布图"""fig, axes = plt.subplots(2, 2, figsize=(12, 10))E_magnitude = np.abs(electric_field)# XY平面im1 = axes[0, 0].imshow(E_magnitude[:, :, E_magnitude.shape[2] // 2].T,cmap='plasma', origin='lower',extent=[0, wood_sample.tangential * 1000, 0, wood_sample.radial * 1000])axes[0, 0].set_title('(a) XY Plane Electric Field')axes[0, 0].set_xlabel('Tangential (mm)')axes[0, 0].set_ylabel('Radial (mm)')plt.colorbar(im1, ax=axes[0, 0], label='Electric Field (V/m)')# XZ平面im2 = axes[0, 1].imshow(E_magnitude[:, E_magnitude.shape[1] // 2, :].T,cmap='plasma', origin='lower',extent=[0, wood_sample.tangential * 1000, 0, wood_sample.longitudinal * 1000])axes[0, 1].set_title('(b) XZ Plane Electric Field')axes[0, 1].set_xlabel('Tangential (mm)')axes[0, 1].set_ylabel('Longitudinal (mm)')plt.colorbar(im2, ax=axes[0, 1], label='Electric Field (V/m)')# YZ平面im3 = axes[1, 0].imshow(E_magnitude[E_magnitude.shape[0] // 2, :, :].T,cmap='plasma', origin='lower',extent=[0, wood_sample.radial * 1000, 0, wood_sample.longitudinal * 1000])axes[1, 0].set_title('(c) YZ Plane Electric Field')axes[1, 0].set_xlabel('Radial (mm)')axes[1, 0].set_ylabel('Longitudinal (mm)')plt.colorbar(im3, ax=axes[1, 0], label='Electric Field (V/m)')# 功率密度power_density = 0.5 * 2 * pi * 2.45e9 * epsilon_0 * 0.1 * E_magnitude ** 2im4 = axes[1, 1].imshow(power_density[:, :, power_density.shape[2] // 2].T,cmap='viridis', origin='lower',extent=[0, wood_sample.tangential * 1000, 0, wood_sample.radial * 1000])axes[1, 1].set_title('(d) Power Density Distribution')axes[1, 1].set_xlabel('Tangential (mm)')axes[1, 1].set_ylabel('Radial (mm)')plt.colorbar(im4, ax=axes[1, 1], label='Power Density (W/m³)')plt.suptitle(f'Microwave Power {microwave_power} kW - Electric Field Distribution')plt.tight_layout()plt.show()@staticmethoddef plot_temperature_distribution(results: Dict, wood_sample: WoodSample):"""绘制温度分布图"""fig, axes = plt.subplots(2, 2, figsize=(12, 10))T_final = results['temperature_history'][-1]# XY平面im1 = axes[0, 0].imshow(T_final[:, :, T_final.shape[2] // 2].T,cmap='hot', origin='lower',extent=[0, wood_sample.tangential * 1000, 0, wood_sample.radial * 1000])axes[0, 0].set_title('(a) XY Plane Temperature')axes[0, 0].set_xlabel('Tangential (mm)')axes[0, 0].set_ylabel('Radial (mm)')plt.colorbar(im1, ax=axes[0, 0], label='Temperature (°C)')# XZ平面im2 = axes[0, 1].imshow(T_final[:, T_final.shape[1] // 2, :].T,cmap='hot', origin='lower',extent=[0, wood_sample.tangential * 1000, 0, wood_sample.longitudinal * 1000])axes[0, 1].set_title('(b) XZ Plane Temperature')axes[0, 1].set_xlabel('Tangential (mm)')axes[0, 1].set_ylabel('Longitudinal (mm)')plt.colorbar(im2, ax=axes[0, 1], label='Temperature (°C)')# YZ平面im3 = axes[1, 0].imshow(T_final[T_final.shape[0] // 2, :, :].T,cmap='hot', origin='lower',extent=[0, wood_sample.radial * 1000, 0, wood_sample.longitudinal * 1000])axes[1, 0].set_title('(c) YZ Plane Temperature')axes[1, 0].set_xlabel('Radial (mm)')axes[1, 0].set_ylabel('Longitudinal (mm)')plt.colorbar(im3, ax=axes[1, 0], label='Temperature (°C)')# 等值线x = np.linspace(0, wood_sample.tangential * 1000, T_final.shape[0])y = np.linspace(0, wood_sample.radial * 1000, T_final.shape[1])X, Y = np.meshgrid(x, y)Z = T_final[:, :, T_final.shape[2] // 2].Tcontour = axes[1, 1].contourf(X, Y, Z, levels=20, cmap='hot')axes[1, 1].set_title('(d) Temperature Contours')axes[1, 1].set_xlabel('Tangential (mm)')axes[1, 1].set_ylabel('Radial (mm)')plt.colorbar(contour, ax=axes[1, 1], label='Temperature (°C)')plt.suptitle(f'Microwave Power {results["microwave_power"]} kW - Temperature Distribution after {results["total_time"]}s')plt.tight_layout()plt.show()@staticmethoddef plot_temperature_evolution(results: Dict):"""绘制温度随时间演化图"""plt.figure(figsize=(10, 6))surface_temp = [T[0, T.shape[1] // 2, T.shape[2] // 2] for T in results['temperature_history']]core_temp = [T[T.shape[0] // 2, T.shape[1] // 2, T.shape[2] // 2] for T in results['temperature_history']]sub_surface_temp = [T[1, T.shape[1] // 2, T.shape[2] // 2] for T in results['temperature_history']]plt.plot(results['time_points'], surface_temp, 'b-', label='Surface Temperature', linewidth=2)plt.plot(results['time_points'], sub_surface_temp, 'g-', label='Sub-surface Temperature', linewidth=2)plt.plot(results['time_points'], core_temp, 'r-', label='Core Temperature', linewidth=2)plt.xlabel('Time (s)')plt.ylabel('Temperature (°C)')plt.title(f'Microwave Power {results["microwave_power"]} kW - Temperature Evolution')plt.legend()plt.grid(True, alpha=0.3)metrics = results['metrics']stats_text = (f'Final Surface: {metrics["surface_temperature"]:.1f}°C\n'f'Final Core: {metrics["core_temperature"]:.1f}°C\n'f'Max Difference: {metrics["temperature_difference"]:.1f}°C\n'f'Heating Rate: {metrics["heating_rate"]:.2f}°C/s')plt.annotate(stats_text, xy=(0.02, 0.98), xycoords='axes fraction',bbox=dict(boxstyle="round,pad=0.3", fc="white", alpha=0.8),verticalalignment='top')plt.tight_layout()plt.show()@staticmethoddef plot_power_comparison(all_results: Dict):"""绘制不同功率条件下的比较图"""fig, axes = plt.subplots(2, 2, figsize=(12, 10))powers = list(all_results.keys())# 最大温度max_temps = [np.max(results['temperature_history'][-1]) for results in all_results.values()]axes[0, 0].plot(powers, max_temps, 'bo-', linewidth=2)axes[0, 0].set_xlabel('Microwave Power (kW)')axes[0, 0].set_ylabel('Max Temperature (°C)')axes[0, 0].set_title('(a) Power vs Max Temperature')axes[0, 0].grid(True, alpha=0.3)# 温差temp_diffs = [results['metrics']['temperature_difference'] for results in all_results.values()]axes[0, 1].plot(powers, temp_diffs, 'ro-', linewidth=2)axes[0, 1].set_xlabel('Microwave Power (kW)')axes[0, 1].set_ylabel('Surface-Core Difference (°C)')axes[0, 1].set_title('(b) Power vs Temperature Difference')axes[0, 1].grid(True, alpha=0.3)# 加热速率heating_rates = [results['metrics']['heating_rate'] for results in all_results.values()]axes[1, 0].plot(powers, heating_rates, 'go-', linewidth=2)axes[1, 0].set_xlabel('Microwave Power (kW)')axes[1, 0].set_ylabel('Heating Rate (°C/s)')axes[1, 0].set_title('(c) Power vs Heating Rate')axes[1, 0].grid(True, alpha=0.3)# 温度均匀性temp_stds = [results['metrics']['temperature_std'] for results in all_results.values()]axes[1, 1].plot(powers, temp_stds, 'mo-', linewidth=2)axes[1, 1].set_xlabel('Microwave Power (kW)')axes[1, 1].set_ylabel('Temperature STD (°C)')axes[1, 1].set_title('(d) Power vs Temperature Uniformity')axes[1, 1].grid(True, alpha=0.3)plt.suptitle('Microwave Heating Performance Comparison')plt.tight_layout()plt.show()def main():"""主函数"""logger.info("=== Microwave Wood Heating Simulation Started ===")# 设置参数wood_props = WoodProperties()cavity = MicrowaveCavity()wood_sample = WoodSample()# 创建模拟器simulator = MicrowaveWoodHeatingSimulator(wood_props, cavity, wood_sample)# 模拟参数power_levels = [3.6, 4.5, 5.4]heating_time = 30all_results = {}# 执行模拟for power in power_levels:logger.info(f"\n{'=' * 50}")logger.info(f"Simulating microwave power: {power} kW")logger.info(f"{'=' * 50}")results = simulator.simulate_heating_process(power, heating_time)all_results[power] = results# 可视化visualizer = ResultVisualizer()# 绘制详细结果main_power = 3.6if main_power in all_results:main_results = all_results[main_power]visualizer.plot_electric_field_distribution(main_results['electric_field'], main_power, wood_sample)visualizer.plot_temperature_distribution(main_results, wood_sample)visualizer.plot_temperature_evolution(main_results)# 功率比较visualizer.plot_power_comparison(all_results)logger.info("=== Microwave Wood Heating Simulation Completed ===")return all_resultsif __name__ == "__main__":# 执行主程序results = main()# 保存结果摘要summary_data = []for power, result in results.items():summary_data.append({'Power (kW)': power,'Surface Temp (°C)': result['metrics']['surface_temperature'],'Core Temp (°C)': result['metrics']['core_temperature'],'Max Difference (°C)': result['metrics']['temperature_difference'],'Heating Rate (°C/s)': result['metrics']['heating_rate'],'Simulation Time (s)': result['simulation_time']})df_summary = pd.DataFrame(summary_data)print("\nSimulation Summary:")print(df_summary.to_string(index=False))# 保存详细结果到CSV文件output_dir = "simulation_results"os.makedirs(output_dir, exist_ok=True)for power, result in results.items():# 保存温度数据temp_data = result['temperature_history'][-1].flatten()positions = []for i in range(result['temperature_history'].shape[1]):for j in range(result['temperature_history'].shape[2]):for k in range(result['temperature_history'].shape[3]):positions.append(f"({i},{j},{k})")df_temp = pd.DataFrame({'Position': positions,'Temperature (°C)': temp_data})filename = f"{output_dir}/temperature_distribution_{power}kW.csv"df_temp.to_csv(filename, index=False)logger.info(f"Temperature data saved to: {filename}")logger.info("All simulation results saved successfully")