微波加热内部温度场的电磁−热耦合模拟

一个三维微波木材加热模拟器，采用有限差分法求解热传导方程，结合电磁场分布计算木材加热过程。主要功能包括：1) 建立三维网格系统并初始化参数；2) 计算电磁场分布和功率密度；3) 求解热传导方程；4) 可视化3D温度场、电场分布和功率密度。模拟器可分析不同微波功率(3.6-5.4kW)下的加热效果，输出温度分布、加热速率等关键指标。通过稀疏矩阵求解器和边界条件处理，实现了高效的三维热传导模拟，并提供了多种3D可视化方法展示模拟结果。

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import scipy.sparse as sp
import scipy.sparse.linalg as spla
from scipy.constants import pi, c, epsilon_0
import pandas as pd
import time
from typing import Dict
import logging
from dataclasses import dataclass
from enum import Enum# 配置日志
logging.basicConfig(level=logging.INFO, format='%(asctime)s - %(levelname)s - %(message)s')
logger = logging.getLogger(__name__)class WoodSpecies(Enum):PAULOWNIA = "paulownia"OAK = "oak"PINE = "pine"class MicrowaveMode(Enum):TE10 = "TE10"TM01 = "TM01"@dataclass
class WoodProperties:density: float = 237.0specific_heat: float = 1200.0thermal_conductivity: float = 0.105epsilon_prime: float = 2.0epsilon_double_prime: float = 0.1electrical_conductivity: float = 0.01@dataclass
class MicrowaveCavity:length: float = 0.5width: float = 0.4height: float = 0.3frequency: float = 2.45e9mode: MicrowaveMode = MicrowaveMode.TE10@dataclass
class WoodSample:tangential: float = 0.1radial: float = 0.05longitudinal: float = 0.08species: WoodSpecies = WoodSpecies.PAULOWNIAinitial_temperature: float = 20.0class 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_resolutionself.environment_temperature = 20.0self.convection_coefficient = 10.0self.initial_temperature = wood_sample.initial_temperatureself.setup_parameters()self.initialize_mesh()self.setup_matrices()logger.info("Microwave heating simulator initialized")def setup_parameters(self):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) + 1alpha = 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)logger.info(f"Grid dimensions: {self.nx} × {self.ny} × {self.nz}")logger.info(f"Time step: {self.dt:.4f} s")def initialize_mesh(self):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_temperaturedef setup_matrices(self):self.setup_thermal_matrices()self.setup_electromagnetic_matrices()def setup_thermal_matrices(self):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)def setup_electromagnetic_matrices(self):self.propagation_constant = self._calculate_propagation_constant()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:a, b = 0.1, 0.05m, n = 1, 0v_c = (c / 2) * np.sqrt((m/a)**2 + (n/b)**2)v = self.cavity.frequencyif v > v_c:beta = 2 * pi * np.sqrt((v/c)**2 - (v_c/c)**2)alpha_d = (pi * v * np.sqrt(self.wood_properties.epsilon_prime * self.wood_properties.epsilon_double_prime)) / creturn beta + 1j * alpha_dreturn 0 + 1j * 10def calculate_electric_field_distribution(self, microwave_power: float) -> np.ndarray:power_scaling = np.sqrt(microwave_power / 3.6)E_field = np.zeros((self.nx, self.ny, self.nz), dtype=complex)waveguide_positions = [(0.1, 0.2, 0.15),(0.3, 0.1, 0.15), (0.4, 0.3, 0.15)]for i, pos in enumerate(waveguide_positions):phase_shift = i * 2 * pi / 3for i_idx in range(self.nx):for j_idx in range(self.ny):for k_idx in range(self.nz):dx = self.x[i_idx] - pos[0]dy = self.y[j_idx] - pos[1]dz = self.z[k_idx] - pos[2]distance = np.sqrt(dx**2 + dy**2 + dz**2)if distance > 0:attenuation = 1.0 / (1.0 + distance)E_waveguide = (np.sin(pi * self.x[i_idx] / self.wood_sample.tangential) *np.exp(1j * (np.real(self.propagation_constant) * distance + phase_shift)))reflection_x = np.sin(2 * pi * self.x[i_idx] / self.wood_sample.tangential)reflection_y = np.cos(2 * pi * self.y[j_idx] / self.wood_sample.radial)reflection_z = np.sin(2 * pi * self.z[k_idx] / self.wood_sample.longitudinal)E_field[i_idx, j_idx, k_idx] += (E_waveguide * attenuation * reflection_x * reflection_y * reflection_z)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: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: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: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,'mesh_info': {'nx': self.nx, 'ny': self.ny, 'nz': self.nz,'x': self.x, 'y': self.y, 'z': self.z}}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]metrics = {'surface_temperature': surface_temp,'core_temperature': core_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] if time_points[-1] > 0 else 0}return metricsclass ThreeDVisualizer:@staticmethoddef plot_3d_temperature_distribution(results: Dict, time_index: int = -1):fig = plt.figure(figsize=(16, 12))T = results['temperature_history'][time_index]x, y, z = results['mesh_info']['x'], results['mesh_info']['y'], results['mesh_info']['z']X, Y, Z = np.meshgrid(x, y, z, indexing='ij')# 1. 三维散点温度分布ax1 = fig.add_subplot(231, projection='3d')scatter = ax1.scatter(X[::2], Y[::2], Z[::2], c=T[::2], cmap='hot', alpha=0.6, s=10)ax1.set_title('(a) 3D Temperature Distribution')ax1.set_xlabel('Tangential (m)')ax1.set_ylabel('Radial (m)')ax1.set_zlabel('Longitudinal (m)')plt.colorbar(scatter, ax=ax1, label='Temperature (C)')# 2. 修复的切片视图 - 确保形状匹配ax2 = fig.add_subplot(232, projection='3d')slice_idx = T.shape[2] // 2# 正确的网格生成方式X_slice, Y_slice = np.meshgrid(x, y, indexing='ij')Z_slice = T[:, :, slice_idx]# 修复形状匹配问题surf = ax2.plot_surface(X_slice, Y_slice, Z_slice, cmap='viridis', alpha=0.8)ax2.set_title('(b) Longitudinal Slice Temperature')ax2.set_xlabel('Tangential (m)')ax2.set_ylabel('Radial (m)')ax2.set_zlabel('Temperature (C)')plt.colorbar(surf, ax=ax2, label='Temperature (C)')# 3. 温度梯度矢量场ax3 = fig.add_subplot(233, projection='3d')grad_x, grad_y, grad_z = np.gradient(T)skip = max(1, T.shape[0]//8)ax3.quiver(X[::skip, ::skip, ::skip], Y[::skip, ::skip, ::skip], Z[::skip, ::skip, ::skip],grad_x[::skip, ::skip, ::skip], grad_y[::skip, ::skip, ::skip], grad_z[::skip, ::skip, ::skip], length=0.01, normalize=True, color='blue', alpha=0.6)ax3.set_title('(c) Temperature Gradient Field')ax3.set_xlabel('Tangential (m)')ax3.set_ylabel('Radial (m)')ax3.set_zlabel('Longitudinal (m)')# 4. 简化的等温面可视化ax4 = fig.add_subplot(234, projection='3d')levels = np.linspace(np.min(T), np.max(T), 5)for i, level in enumerate(levels[1:-1]):mask = T >= levelpoints = np.argwhere(mask)if len(points) > 1000:points = points[::len(points)//1000]if len(points) > 0:ax4.scatter(x[points[:,0]], y[points[:,1]], z[points[:,2]], c=[level]*len(points), cmap='viridis', alpha=0.3, s=1)ax4.set_title('(d) Isothermal Surface')ax4.set_xlabel('Tangential (m)')ax4.set_ylabel('Radial (m)')ax4.set_zlabel('Longitudinal (m)')# 5. 电场分布ax5 = fig.add_subplot(235, projection='3d')E_field = np.abs(results['electric_field'])scatter_e = ax5.scatter(X[::2], Y[::2], Z[::2], c=E_field[::2], cmap='plasma', alpha=0.6, s=10)ax5.set_title('(e) 3D Electric Field')ax5.set_xlabel('Tangential (m)')ax5.set_ylabel('Radial (m)')ax5.set_zlabel('Longitudinal (m)')plt.colorbar(scatter_e, ax=ax5, label='Electric Field (V/m)')# 6. 功率密度分布ax6 = fig.add_subplot(236, projection='3d')power_density = results['power_density']scatter_p = ax6.scatter(X[::2], Y[::2], Z[::2], c=power_density[::2], cmap='inferno', alpha=0.6, s=10)ax6.set_title('(f) 3D Power Density')ax6.set_xlabel('Tangential (m)')ax6.set_ylabel('Radial (m)')ax6.set_zlabel('Longitudinal (m)')plt.colorbar(scatter_p, ax=ax6, label='Power Density (W/m³)')plt.suptitle(f'Microwave Power {results["microwave_power"]} kW - 3D Analysis')plt.tight_layout()plt.show()@staticmethoddef plot_temperature_evolution_3d(results: Dict):fig = plt.figure(figsize=(12, 8))ax = fig.add_subplot(111, projection='3d')T_history = results['temperature_history']time_points = results['time_points']x, y, z = results['mesh_info']['x'], results['mesh_info']['y'], results['mesh_info']['z']# 选择核心点跟踪温度演化core_i, core_j, core_k = len(x)//2, len(y)//2, len(z)//2core_temps = [T[core_i, core_j, core_k] for T in T_history]# 创建时间-空间三维图 - 修复形状匹配time_grid, x_grid = np.meshgrid(time_points, x, indexing='ij')temp_surface = np.array([T[:, core_j, core_k] for T in T_history])surf = ax.plot_surface(time_grid, x_grid, temp_surface, cmap='coolwarm', alpha=0.8)ax.set_xlabel('Time (s)')ax.set_ylabel('Tangential Position (m)')ax.set_zlabel('Temperature (C)')ax.set_title('Temperature Evolution Over Time and Space')plt.colorbar(surf, ax=ax, label='Temperature (C)')plt.tight_layout()plt.show()@staticmethoddef plot_comparative_analysis(all_results: Dict):fig = plt.figure(figsize=(15, 10))powers = list(all_results.keys())# 1. 最大温度对比ax1 = fig.add_subplot(231, projection='3d')max_temps = [np.max(results['temperature_history'][-1]) for results in all_results.values()]x_pos = range(len(powers))y_pos = [0] * len(powers)z_pos = [0] * len(powers)dx = [0.8] * len(powers)dy = [0.8] * len(powers)dz = max_tempsax1.bar3d(x_pos, y_pos, z_pos, dx, dy, dz, alpha=0.8)ax1.set_xticks([i + 0.4 for i in range(len(powers))])ax1.set_xticklabels([f'{p}kW' for p in powers])ax1.set_title('(a) Max Temperature Comparison')ax1.set_xlabel('Microwave Power')ax1.set_ylabel('')ax1.set_zlabel('Temperature (C)')# 2. 温差对比ax2 = fig.add_subplot(232, projection='3d')temp_diffs = [results['metrics']['temperature_difference'] for results in all_results.values()]ax2.bar3d(x_pos, y_pos, z_pos, dx, dy, temp_diffs, alpha=0.8)ax2.set_xticks([i + 0.4 for i in range(len(powers))])ax2.set_xticklabels([f'{p}kW' for p in powers])ax2.set_title('(b) Core-Surface Difference')ax2.set_xlabel('Microwave Power')ax2.set_ylabel('')ax2.set_zlabel('Temperature Diff (C)')# 3. 加热速率对比ax3 = fig.add_subplot(233, projection='3d')heating_rates = [results['metrics']['heating_rate'] for results in all_results.values()]ax3.bar3d(x_pos, y_pos, z_pos, dx, dy, heating_rates, alpha=0.8)ax3.set_xticks([i + 0.4 for i in range(len(powers))])ax3.set_xticklabels([f'{p}kW' for p in powers])ax3.set_title('(c) Heating Rate Comparison')ax3.set_xlabel('Microwave Power')ax3.set_ylabel('')ax3.set_zlabel('Heating Rate (C/s)')# 4. 温度均匀性ax4 = fig.add_subplot(234, projection='3d')temp_stds = [results['metrics']['temperature_std'] for results in all_results.values()]ax4.bar3d(x_pos, y_pos, z_pos, dx, dy, temp_stds, alpha=0.8)ax4.set_xticks([i + 0.4 for i in range(len(powers))])ax4.set_xticklabels([f'{p}kW' for p in powers])ax4.set_title('(d) Temperature Uniformity')ax4.set_xlabel('Microwave Power')ax4.set_ylabel('')ax4.set_zlabel('Temperature STD (C)')# 5. 三维散点对比ax5 = fig.add_subplot(235, projection='3d')for i, (power, results) in enumerate(all_results.items()):T_final = results['temperature_history'][-1]x, y, z = results['mesh_info']['x'], results['mesh_info']['y'], results['mesh_info']['z']X, Y, Z = np.meshgrid(x, y, z, indexing='ij')ax5.scatter(X[::3], Y[::3], Z[::3], c=T_final[::3], cmap='hot', alpha=0.3, s=5, label=f'{power}kW')ax5.set_title('(e) Temperature Distribution Comparison')ax5.set_xlabel('Tangential (m)')ax5.set_ylabel('Radial (m)')ax5.set_zlabel('Longitudinal (m)')ax5.legend()# 6. 温度分布直方图ax6 = fig.add_subplot(236, projection='3d')for i, (power, results) in enumerate(all_results.items()):T_final = results['temperature_history'][-1].flatten()hist, bins = np.histogram(T_final, bins=20)xs = [i] * len(hist)ys = range(len(hist))ax6.bar3d(xs, ys, [0]*len(hist), [0.3]*len(hist), [0.8]*len(hist), hist, alpha=0.6)ax6.set_title('(f) Temperature Histogram')ax6.set_xlabel('Power')ax6.set_ylabel('Temperature Range')ax6.set_zlabel('Frequency')ax6.set_xticks([i + 0.15 for i in range(len(powers))])ax6.set_xticklabels([f'{p}kW' for p in powers])plt.suptitle('Microwave Heating Performance 3D Comparison')plt.tight_layout()plt.show()@staticmethoddef plot_simple_2d_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('Max Temperature Comparison')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('Core-Surface Diff (C)')axes[0, 1].set_title('Temperature Uniformity')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('Heating Efficiency')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('Temperature Distribution')axes[1, 1].grid(True, alpha=0.3)plt.suptitle('Microwave Heating Performance 2D Comparison')plt.tight_layout()plt.show()def main():logger.info("=== Microwave Wood Heating 3D 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"Simulating microwave power: {power} kW")results = simulator.simulate_heating_process(power, heating_time)all_results[power] = resultsvisualizer = ThreeDVisualizer()main_power = 3.6if main_power in all_results:main_results = all_results[main_power]try:visualizer.plot_3d_temperature_distribution(main_results)visualizer.plot_temperature_evolution_3d(main_results)except Exception as e:logger.error(f"3D visualization error: {e}")logger.info("Using 2D comparison as alternative")visualizer.plot_simple_2d_comparison(all_results)try:visualizer.plot_comparative_analysis(all_results)except Exception as e:logger.error(f"Comparative analysis error: {e}")visualizer.plot_simple_2d_comparison(all_results)logger.info("=== Simulation Completed ===")summary_data = []for power, result in all_results.items():summary_data.append({'Power (kW)': power,'Surface Temp (C)': result['metrics']['surface_temperature'],'Core Temp (C)': result['metrics']['core_temperature'],'Max Temp (C)': result['metrics']['max_temperature'],'Temp 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))return all_resultsif __name__ == "__main__":results = main()