Python和Matplotlib电磁一维二维三维模拟

电磁一维二维三维模拟"/>

Python和Matplotlib电磁一维二维三维模拟

自由空间中传播的脉冲一维模拟

自由空间的时间相关麦克斯韦旋度方程为

∂ E ∂ t = 1 ε 0 ∇ × H ∂ H ∂ t = − 1 μ 0 ∇ × E . \begin{aligned}\frac{\partial \boldsymbol{E}}{\partial t} &=\frac{1}{\varepsilon_{0}} \nabla \times \boldsymbol{H} \\\frac{\partial \boldsymbol{H}}{\partial t} &=-\frac{1}{\mu_{0}} \nabla \times \boldsymbol{E} .\end{aligned} ∂t∂E​∂t∂H​​=ε0​1​∇×H=−μ0​1​∇×E.​

E 和 H 是三个维度的向量，因此，一般来说，等式 (1.1a) 和 (1.1b) 分别代表三个方程。 我们将从仅使用 Ex 和 Hy 的简单一维案例开始，因此 等式(1.1a) 和 (1.1b) 变为

∂ E x ∂ t = − 1 ε 0 ∂ H y ∂ z , ∂ H y ∂ t = − 1 μ 0 ∂ E x ∂ z . \begin{aligned}&\frac{\partial E_{x}}{\partial t}=-\frac{1}{\varepsilon_{0}} \frac{\partial H_{y}}{\partial z}, \\&\frac{\partial H_{y}}{\partial t}=-\frac{1}{\mu_{0}} \frac{\partial E_{x}}{\partial z} .\end{aligned} ​∂t∂Ex​​=−ε0​1​∂z∂Hy​​,∂t∂Hy​​=−μ0​1​∂z∂Ex​​.​

这些是沿z方向行驶的平面波的方程，其电场在X方向上定向，并以y方向定向磁场。

对时间和空间导数进行中心差分近似，给出

E x n + 1 / 2 ( k ) − E x n − 1 / 2 ( k ) Δ t = − 1 ε 0 H y n ( k + 1 2 ) − H y n ( k − 1 2 ) Δ x H y n + 1 ( k + 1 2 ) − H y n ( k + 1 2 ) Δ t = − 1 μ 0 E x n + 1 / 2 ( k + 1 ) − E x n + 1 / 2 ( k ) Δ x \begin{aligned}\frac{E_{x}^{n+1 / 2}(k)-E_{x}^{n-1 / 2}(k)}{\Delta t} &=-\frac{1}{\varepsilon_{0}} \frac{H_{y}^{n}\left(k+\frac{1}{2}\right)-H_{y}^{n}\left(k-\frac{1}{2}\right)}{\Delta x} \\\frac{H_{y}^{n+1}\left(k+\frac{1}{2}\right)-H_{y}^{n}\left(k+\frac{1}{2}\right)}{\Delta t} &=-\frac{1}{\mu_{0}} \frac{E_{x}^{n+1 / 2}(k+1)-E_{x}^{n+1 / 2}(k)}{\Delta x}\end{aligned} ΔtExn+1/2​(k)−Exn−1/2​(k)​ΔtHyn+1​(k+21​)−Hyn​(k+21​)​​=−ε0​1​ΔxHyn​(k+21​)−Hyn​(k−21​)​=−μ0​1​ΔxExn+1/2​(k+1)−Exn+1/2​(k)​​

代码模拟

import numpy as np
from math import exp
from matplotlib import pyplot as pltke = 200
ex = np.zeros(ke)
hy = np.zeros(ke)# Pulse parameters
kc = int(ke / 2)
t0 = 40
spread = 12nsteps = 100# Main FDTD Loop
for time_step in range(1, nsteps + 1):# Calculate the Ex fieldfor k in range(1, ke):ex[k] = ex[k] + 0.5 * (hy[k - 1] - hy[k])# Put a Gaussian pulse in the middlepulse = exp(-0.5 * ((t0 - time_step) / spread) ** 2)ex[kc] = pulse# Calculate the Hy fieldfor k in range(ke - 1):hy[k] = hy[k] + 0.5 * (ex[k] - ex[k + 1])# Plot the outputs as shown in Fig. 1.2
plt.rcParams['font.size'] = 12
plt.figure(figsize=(8, 3.5))plt.subplot(212)
plt.plot(hy, color='k', linewidth=1)
plt.ylabel('H$_y$', fontsize='14')
plt.xlabel('FDTD cells')
plt.xticks(np.arange(0, 201, step=20))
plt.xlim(0, 200)
plt.yticks(np.arange(-1, 1.2, step=1))
plt.ylim(-1.2, 1.2)plt.subplots_adjust(bottom=0.2, hspace=0.45)
plt.show()

通量密度的一维模拟

横向电磁二维模拟

矢量场电磁三维模拟

环形相控阵热疗模拟

源代码

参阅 - 亚图跨际