简介：这是一份面向电磁仿真初学者与科研人员的三维时域有限差分（3D FDTD）MATLAB实现代码包，聚焦于DNG（可能指双负材料或特定算法缩写）相关电磁响应建模，适用于天线设计、微波器件分析及计算电磁学教学实践。压缩包仅含2个文件（1个MATLAB主程序DNG.m + 1个license.txt），总大小仅2KB，轻量简洁，便于快速部署与二次开发；其中DNG.m封装了完整的3D FDTD核心流程——包括Yee网格初始化、Courant稳定条件约束下的电/磁场交替更新、基础边界处理及源激励设置，可直接运行并支持参数化修改。已有249人学习下载，适合具备电磁场理论基础与MATLAB编程能力的学习者，用于理解FDTD离散原理、验证数值稳定性、开展简单介质建模实验或作为课程设计脚本基础。
1. 项目背景与核心目标：从“DNG.zip”到三维FDTD仿真
最近在整理一个老项目的资料时，翻出了一个名为“DNG.zip”的压缩包。看到这个名字，很多从事电磁场仿真或者超材料研究的朋友可能会心一笑。DNG，即“双负材料”，指的是介电常数和磁导率同时为负的人工复合材料，它能够实现许多奇特的电磁现象，比如负折射、完美透镜等。这个压缩包里，存放的正是一套用MATLAB编写的、用于模拟三维DNG材料电磁特性的时域有限差分法代码。
FDTD，全称Finite-Difference Time-Domain，时域有限差分法，是计算电磁学领域的基石方法之一。它的核心思想非常直观：将麦克斯韦方程组中的微分算子用差分近似代替，从而在离散的时空网格上，一步步地“推进”计算电磁场的演化。这种方法特别适合处理复杂的几何结构、非线性材料以及宽频带响应问题。而“3D”则意味着这是一个全三维的仿真，电场和磁场在空间的三个方向上都有分量，模拟的复杂度和计算量远高于二维情况。
这个项目的目标很明确：构建一个能够稳定、准确运行的三维FDTD仿真环境，专门用于分析DNG材料构成的简单结构（如平板、立方体）在电磁波照射下的散射、透射和场分布特性。对于研究者或学生而言，拥有一套自己理解透彻、可以随意修改的底层FDTD代码，远比单纯使用商业软件的黑箱来得有价值。它不仅能让你验证理论，更能让你深入理解数值计算中的每一个细节，比如稳定性条件、吸收边界设置、激励源引入以及材料参数的处理。
2. 三维FDTD算法的核心原理与Yee网格
要读懂甚至修改“DNG.zip”里的代码，必须对三维FDTD的底层原理有清晰的把握。这一切都始于1966年Kane S. Yee提出的Yee网格。这个网格的巧妙之处在于，它将电场E和磁场H的各个分量在空间和时间上进行了交错放置。
2.1 Yee网格的空间与时间交错
在一个三维直角坐标系中，我们用一个立方体网格来离散空间。Yee规定：
- 电场E的三个分量（Ex, Ey, Ez）分别位于立方体棱边的中心。
- 磁场H的三个分量（Hx, Hy, Hz）分别位于立方体表面的中心。
更具体地说，假设网格索引为(i, j, k)，空间步长为Δx, Δy, Δz，那么：
- Ex分量位于
(i+1/2, j, k)。 - Ey分量位于
(i, j+1/2, k)。 - Ez分量位于
(i, j, k+1/2)。 - Hx分量位于
(i, j+1/2, k+1/2)。 - Hy分量位于
(i+1/2, j, k+1/2)。 - Hz分量位于
(i+1/2, j+1/2, k)。
这种交错放置使得每个磁场分量被四个环绕的电场分量所“包围”，反之亦然。这完美契合了麦克斯韦方程组中安培环路定律和法拉第电磁感应定律的积分形式：磁场的变化由电场的旋度决定，电场的变化由磁场的旋度决定。在差分计算时，我们可以直接用相邻点的场值来计算旋度，无需复杂的插值。
在时间上，Yee网格也采用了交错采样。通常，我们假设电场E在整数时间步n*Δt被采样，而磁场H则在半整数时间步(n+1/2)*Δt被采样。这样，计算过程就形成了一个“蛙跳”式的迭代：
- 在已知第n步电场E^n和第n-1/2步磁场H^{n-1/2}的情况下，利用法拉第定律更新得到第n+1/2步的磁场H^{n+1/2}。
- 接着，利用安培环路定律，由已知的E^n和刚求出的H^{n+1/2}，更新得到第n+1步的电场E^{n+1}。 如此循环往复。
2.2 从麦克斯韦方程组到更新方程
我们来看各向同性、无耗散介质中的麦克斯韦旋度方程： ∇ × E = -μ ∂H/∂t ∇ × H = ε ∂E/∂t
将旋度算子展开，并把微分用中心差分近似，我们就可以得到六个标量更新方程。以Ex和Hz的更新为例（假设介质参数ε和μ在网格点处定义）：
Hz分量的更新（位于(i+1/2, j+1/2, k)）：H_z^{n+1/2}(i+1/2,j+1/2,k) = H_z^{n-1/2}(i+1/2,j+1/2,k) + (Δt/μ) * [ (E_x^n(i+1/2, j+1, k) - E_x^n(i+1/2, j, k)) / Δy
- (E_y^n(i+1, j+1/2, k) - E_y^n(i, j+1/2, k)) / Δx ]
Ex分量的更新（位于(i+1/2, j, k)）：E_x^{n+1}(i+1/2,j,k) = E_x^{n}(i+1/2,j,k) + (Δt/ε) * [ (H_z^{n+1/2}(i+1/2, j+1/2, k) - H_z^{n+1/2}(i+1/2, j-1/2, k)) / Δy
- (H_y^{n+1/2}(i+1/2, j, k+1/2) - H_y^{n+1/2}(i+1/2, j, k-1/2)) / Δz ]
其他四个分量的更新方程形式类似，只是旋度计算所依赖的相邻场分量不同。在MATLAB实现中，我们通常会为Ex, Ey, Ez, Hx, Hy, Hz分别创建三维数组，然后利用数组切片操作来高效地实现这些差分计算，避免使用低效的多重循环。
2.3 稳定性条件与数值色散
FDTD方法不是无条件稳定的。著名的Courant-Friedrichs-Lewy (CFL)稳定性条件在三维各向同性网格中表现为： Δt ≤ 1 / (c * sqrt( (1/Δx)^2 + (1/Δy)^2 + (1/Δz)^2 )) 其中c是介质中的光速。对于均匀立方体网格（Δx=Δy=Δz=Δ），条件简化为 Δt ≤ Δ / (c * sqrt(3))。在实际编程中，我们通常会取一个安全系数，比如0.99倍的极限值，以确保长期迭代的稳定性。
另一个重要概念是数值色散。在FDTD中，不同频率的波以略微不同的相速度在网格中传播，这会导致脉冲波形在传播过程中发生畸变。网格分辨率越高（Δ越小），数值色散效应越弱。一个经验法则是，为了较好地模拟一个波长为λ的电磁波，网格步长Δ应小于λ/10，对于高精度要求，可能需要λ/20甚至更小。这对于模拟DNG材料尤其重要，因为其特性往往在特定频段附近变化剧烈，需要更精细的网格来捕捉。
3. DNG材料的FDTD实现：Drude与Lorentz模型
常规材料的介电常数ε和磁导率μ是常数或随频率缓变的实数。但DNG材料要求ε(ω) < 0 且 μ(ω) < 0，这通常只在谐振频率附近实现。在FDTD这种时域方法中，我们不能直接使用频域的参数，必须将频率相关的ε(ω)和μ(ω)转换成时域的本构关系，通常引入极化强度P和磁化强度M作为辅助变量。
3.1 色散材料的时域处理
对于色散材料，电通量密度D与电场E的关系不再是简单的D=εE，而是通过卷积或微分方程关联。最常用的方法是辅助微分方程法。以Drude模型为例，它常用来描述等离子体或DNG材料的介电行为： ε(ω) = ε∞ - ω_p^2 / (ω^2 - iωγ) 其中ε∞是高频极限介电常数，ω_p是等离子体频率，γ是碰撞频率。对应的时域微分方程为： d^2P/dt^2 + γ dP/dt = ε0 ω_p^2 E 这里P是极化强度，满足 D = ε0 ε∞ E + P。
在FDTD中，我们需要将这个二阶微分方程差分化，与主更新方程联立求解。具体步骤是：
- 在更新电场E之前，先利用当前的E和过去的P值，通过差分形式的Drude模型方程更新P。
- 然后利用更新后的P和E来计算D（如果使用D-H形式的更新方程）或直接用于修正电场更新方程中的系数。
3.2 在MATLAB中实现Drude模型
在“DNG.zip”的代码中，你可能会看到类似以下结构的片段（概念性代码）：
% 参数定义 epsilon_inf = 1.0; omega_p = 2*pi*1e15; % 等离子体频率 gamma = 0.01*omega_p; % 阻尼系数 dt = ...; % 时间步长 % 计算与Drude模型相关的更新系数 C_E = (2 - gamma*dt) / (2 + gamma*dt); % P的更新系数1 C_P = (2*epsilon0*omega_p^2*dt^2) / (2 + gamma*dt); % P的更新系数2 % 在时间循环内 for n = 1: Nsteps % 1. 更新磁场H (与传统FDTD相同) Hx = Hx + (dt/mu0) * ( (Ey(:,:,2:end)-Ey(:,:,1:end-1))/dz - (Ez(:,2:end,:)-Ez(:,1:end-1,:))/dy ); % ... 更新Hy, Hz % 2. 更新辅助变量极化强度P (在DNG材料区域) % 假设Px, Py, Pz是与Ex, Ey, Ez同位置的三维数组 Px(DNG_region) = C_E * Px(DNG_region) + C_P * Ex(DNG_region); % ... 同样更新Py, Pz % 3. 更新电场E (考虑P的贡献) % 电通量密度 D = epsilon0*epsilon_inf*E + P % 我们需要从更新后的H计算D的旋度，然后反解出E % 通常做法是重新整理更新方程，将P作为已知源项 curl_Hx = (Hz(:,2:end,:)-Hz(:,1:end-1,:))/dy - (Hy(:,:,2:end)-Hy(:,:,1:end-1))/dz; Ex_new = (Ex + (dt/(epsilon0*epsilon_inf)) * curl_Hx - (dt/(epsilon0*epsilon_inf)) * (Px - Px_old)/dt ) / (1 + (dt*gamma/(2*epsilon_inf))); % 这是一个简化示意，实际公式更复杂 % ... 同样更新Ey_new, Ez_new % 4. 场量迭代 Ex = Ex_new; Px_old = Px; % 为下一步保存旧的P值 end这段代码展示了核心思想：在DNG区域，电场的更新不再仅仅依赖于磁场的旋度，还强烈依赖于由辅助微分方程描述的极化强度P的演化历史。磁导率μ的负值特性通常通过类似的模型（如磁性的Drude或Lorentz模型）实现，引入磁化强度M作为辅助变量。
3.3 模型选择与参数设置心得
对于DNG材料，单Drude模型可能不足以在所需频段同时实现ε和μ为负，且满足一定的带宽。更常用的是Lorentz模型或其变体（如Drude-Lorentz组合模型）。Lorentz模型能产生谐振型的色散曲线，更容易在谐振点附近设计出双负特性。
在MATLAB中实现Lorentz模型会稍微复杂一点，因为它涉及二阶极点。其频域形式为： ε(ω) = ε∞ + (ω_p^2) / (ω_0^2 - ω^2 - iωγ) 对应的时域微分方程也是二阶的。编程时，需要引入两个辅助变量（例如极化强度P和它的时间导数dP/dt）来进行差分更新。
实操心得一：参数归一化。在编写和调试这类代码时，强烈建议使用归一化单位。例如，设定光速c=1，空间步长Δx=1，那么时间步长Δt就由CFL条件决定（如0.99/sqrt(3)）。频率、波长等参数也都相对于一个特征长度或时间进行归一化。这能极大简化公式，减少因量纲和数值过大过小带来的错误，也便于在不同尺度的问题间复用代码。
实操心得二：先验证常规材料，再挑战DNG。不要一开始就运行完整的DNG仿真。首先，用你的FDTD代码模拟真空或常规介质（ε, μ为正常数）中的波传播，例如计算平面波在自由空间的传播速度，或者模拟一个点源的辐射场，并与解析解对比。确保基础框架绝对正确后，再将材料模型替换为DNG模型。这样，当结果异常时，你可以迅速定位问题是出在材料模型实现上，还是出在基础的FDTD算法上。
4. 三维FDTD仿真环境的关键组件构建
一套完整、可用的三维FDTD仿真程序，远不止核心更新循环那么简单。它就像一个精密的仪器，需要多个关键组件协同工作。
4.1 激励源：如何让电磁波“产生”
在FDTD中，我们必须通过某种方式在网格中“注入”电磁波。最常用的激励源是软源和硬源。
- 硬源：直接强制设定某个网格点或区域的电场（或磁场）值为一个给定的时间函数，如E_z(source_point) = sin(2pif*t)。这种方法简单粗暴，但会在源点产生强烈的非物理反射，因为强制设定的场值破坏了该点的麦克斯韦方程组。
- 软源：作为电流密度项J或磁流密度项M添加到更新方程中。例如，在电场更新方程中加入一项：E_new = E_old + ... + (Δt/ε) * J。其中J = J0 * sin(2pif*t)。软源更物理，产生的反射较小。
对于宽带脉冲仿真，常用高斯脉冲调制正弦波作为源函数：
t0 = 20*dt; % 脉冲中心时间 spread = 8*dt; % 脉冲宽度 % 在时间循环中计算源值 n = current_time_step; source_time = (n*dt - t0) / spread; gaussian_envelope = exp(-0.5 * source_time^2); carrier = sin(2*pi*f0 * n*dt); Jz(source_position) = J0 * gaussian_envelope * carrier;这样产生的波既有中心频率f0，又有一定的带宽，一次仿真就能通过傅里叶变换得到一定频带内的响应。
实操心得三：源的放置与总场/散射场分离技术。直接把源放在计算区域中间，会同时产生向四周传播的波。为了只研究一个方向入射的平面波，需要使用总场/散射场分离技术。该技术用一个虚拟的连接边界将计算区域分为总场区（包含入射波和散射波）和散射场区（只包含散射波）。源被加在这个连接边界上，精确地注入一个理想的平面波。这是FDTD中模拟平面波入射的标准方法，实现起来需要仔细处理边界上场值的加减。在“DNG.zip”的代码中，如果涉及平面波照射DNG平板，很可能用到了此技术。
4.2 吸收边界条件：让波“消失”而非反射
计算区域总是有限的，当波传播到边界时，如果不做处理就会反射回计算区域，污染内部结果。因此，我们需要在边界处设置吸收边界条件，模拟波无反射地传播到无限远处。最经典的是完全匹配层。
PML的基本思想是在计算区域外围包裹一层特殊的有耗介质层，其波阻抗与相邻内部区域完全匹配，因此无论入射波以何种角度、何种频率入射，都不会在交界处发生反射。波进入PML层后，其振幅会按指数规律迅速衰减，在到达PML外边界前就已衰减到可忽略不计。
在实现上，PML通过在更新方程中引入分场和复坐标拉伸来实现。以电场Ex为例，在PML层内，我们将其分裂为两个子分量Exy和Exz，分别对应由磁场Hy和Hz的旋度产生的部分。然后对这两个子分量引入电导率σ和磁损耗σ*，其更新方程变为： ∂Exy/∂t + σy/ε * Exy = (1/ε) * ∂Hz/∂y ∂Exz/∂t + σz/ε * Exz = -(1/ε) * ∂Hy/∂z 其中σy, σz在PML层内从0逐渐增大到某个最大值。这种分裂场的方法使得PML的实现非常模块化，但代价是增加了需要存储和更新的变量数量（在三维情况下，每个场分量最多分裂成两个子分量）。
实操心得四：PML参数调优。PML的性能取决于其层数N、理论反射系数R0以及电导率剖面的形状（多项式或几何级数增长）。层数越多、R0越小，吸收效果越好，但计算开销也越大。一个常见的起点是设置8-16层PML，R0=1e-6，使用多项式（通常阶数为3或4）剖面。在MATLAB中，你需要为计算区域的六个外表面分别设置PML参数。调试时，可以运行一个点源仿真，观察波传播到PML边界后是否被干净吸收，没有明显的反射波纹回到中心区域。
4.3 数据采集与后处理：从时域场到频域信息
FDTD给出的是时域场值，而我们通常关心频域响应，如散射截面、透射率、反射率、场分布等。
近场监视器：在仿真过程中，我们在关心的位置（如DNG材料前后、内部特定平面）记录电场或磁场随时间变化的数据。这些就是原始的时域信号。
傅里叶变换：仿真结束后，对记录的时域信号进行离散傅里叶变换，得到其频谱。例如，在透射率计算中，我们在DNG材料后方某个平面记录透射场的时域波形E_trans(t)，在入射波路径上（无样品时）相同位置记录入射场的时域波形E_inc(t)。然后分别做FFT得到E_trans(f)和E_inc(f)。透射系数T(f) = |E_trans(f)| / |E_inc(f)|。同理可得反射系数。
场分布可视化：这是MATLAB的强项。使用slice,isosurface,contourslice,quiver3等函数，可以生动地展示三维空间中某个时刻的电场强度、磁场强度或坡印廷矢量的分布。对于DNG材料，一个经典的验证是观察电磁波通过DNG平板时的负折射现象：波前在界面处不是向法线另一侧偏折，而是向同侧偏折。
% 示例：绘制电场强度在x-y中心截面的分布 figure; imagesc(squeeze(abs(Ex(:, Ny/2, :)))'); % 假设Ex是三维数组，取y方向中心切片 axis equal; axis tight; xlabel('x index'); ylabel('z index'); title('|Ex| distribution at y-center'); colorbar;实操心得五：运行时间与信号长度。FDTD仿真需要运行足够长的时间，以确保所有瞬态过程（如从源激发到稳定）结束，并且感兴趣的频率分量有足够的采样。一个经验法则是，仿真时间至少应持续到最慢的传播路径（如经过多次反射）上的波穿过计算区域数个来回。对于频域结果，时域信号的长度决定了频率分辨率。如果关心低频特性，就需要更长的仿真时间。在MATLAB中，可以使用trapz或累加的方式实时计算频域响应，避免存储所有时间步的庞大场数据，但会牺牲一些灵活性。
5. MATLAB代码实现框架与性能优化技巧
一个结构清晰的三维FDTD MATLAB代码，通常包含以下几个模块：
- 初始化模块：定义物理常数、计算区域大小、网格步长、时间步长、总时间步数。初始化场量数组（Ex, Ey, Ez, Hx, Hy, Hz）和材料参数数组（epsilon, mu）。定义DNG材料区域并设置其Drude/Lorentz模型参数。设置PML参数和系数。
- 更新系数预计算模块：根据CFL条件计算Δt。为常规区域和PML区域分别预计算电场和磁场的更新系数。这些系数通常是空间位置的函数，在时间循环外一次性算好，可以大幅提升循环内速度。
- 时间循环主模块：
- 更新磁场H（包括PML层）。
- 在源位置添加激励源（软源）。
- 更新辅助变量（P, M等，针对DNG区域）。
- 更新电场E（包括PML层，并考虑辅助变量的贡献）。
- 在监视器位置记录场值。
- 可选：实时可视化（每N步绘制一次场图，但会显著降低速度）。
- 后处理模块：对记录的时域数据进行FFT，计算S参数、透射/反射率、场分布图等。
性能优化是三维FDTD在MATLAB中面临的巨大挑战。网格点动辄数十万甚至上百万，时间步数上万，直接使用多重循环会慢得无法忍受。
核心技巧：向量化操作。尽可能避免对i, j, k的三重循环。利用MATLAB的数组切片和矩阵运算。例如，更新整个Hz分量的核心差分可以写成：
% 假设 Ex, Ey 是三维数组 curl_E_x = (Ey(:,:,2:end) - Ey(:,:,1:end-1)) / dz; curl_E_y = (Ex(:,2:end,:) - Ex(:,1:end-1,:)) / dy; % 注意：由于Yee网格交错，Hz数组比Ex, Ey在x,y方向上各多一个元素，需要调整索引匹配 Hz(2:end-1, 2:end-1, :) = Chzh * Hz(2:end-1, 2:end-1, :) ... + Chze * ( (Ex(2:end-1, 3:end, :) - Ex(2:end-1, 2:end-1, :))/dy ... - (Ey(3:end, 2:end-1, :) - Ey(2:end-1, 2:end-1, :))/dx );这里Chzh和Chze是预先计算好的更新系数数组（考虑了材料参数和PML）。通过这种向量化操作，更新一个场分量只需要几行代码，速度比循环快几十到上百倍。
进阶技巧：使用MEX函数。对于无法完全向量化的复杂操作（如某些复杂的材料模型更新），可以将这部分核心循环用C或C++编写，编译成MEX文件供MATLAB调用。这能带来数量级的性能提升。但对于学习和研究目的，纯向量化的MATLAB代码通常已足够，且更易于调试和修改。
内存管理心得：三维数组非常消耗内存。在初始化时，使用zeros(Nx, Ny, Nz, 'single')指定单精度浮点数，可以将内存占用减半，对于大多数FDTD仿真，单精度已足够。及时清除不再需要的大变量。如果网格太大，可以考虑使用磁盘存储或分块计算，但这会极大增加复杂度。
6. 典型仿真案例：DNG平板的负折射验证
让我们以一个具体的案例，串联起上述所有内容：仿真一束高斯波束斜入射到一个DNG材料平板上，观察其负折射现象。
步骤1：问题定义与参数设置
- 计算区域：设为 200Δx × 200Δy × 300Δz。Δx=Δy=Δz=10 nm（归一化单位下可设为1）。
- DNG平板：在z方向居中放置，厚度为40个网格。平板材料采用Drude模型，参数设计为在归一化频率0.3附近实现ε和μ同时为负。
- 激励源：使用总场/散射场技术，在z方向某一位置设置连接边界，注入沿z方向传播、在x-y面内具有一定宽度的 Gaussian-Hermite 模式波束（模拟准平面波），入射角设为30度。
- 边界条件：除注入面外，其余五个面均设置10层PML。
- 监视器：
- 在平板前后各放置一个二维场监视器（记录Ez），用于计算反射和透射系数。
- 在多个时间步，记录整个计算区域或某个截面的场分布，用于制作动画。
步骤2：仿真执行与现象观察运行FDTD程序。在时域动画中，你可以看到波束从左侧入射，在DNG平板前表面发生反射和折射。关键现象：折射波束在平板内不是偏向法线另一侧，而是偏向同侧，这就是负折射。波束从平板后表面射出时，再次发生负折射，最终出射波束与入射波束位于法线同侧。
步骤3：数据处理与验证
- 透射/反射谱：对平板前后监视器记录的时域信号做FFT，计算频域的透射率T(ω)和反射率R(ω)。在DNG频段，由于阻抗匹配可能不完美，R可能不为零，但T应显示出该频段有能量通过。
- 折射角测量：在稳态场分布图中，画出平板内外波前的等相位面。测量入射角θ_i和折射角θ_t。根据斯涅尔定律 n1 sinθ_i = n2 sinθ_t。由于n2为负，sinθ_t将为负，即θ_t为负角，证实了负折射。
- 与解析模型对比：将仿真得到的透射谱与基于Drude模型、利用传输矩阵法计算的解析透射谱进行对比。两者在趋势和共振峰位置上应基本一致，差异主要来自数值色散和边界反射。
常见问题与排查：
- 仿真发散：首先检查CFL条件是否满足。然后检查DNG材料参数是否设置得过于极端（如等离子体频率过高、损耗过小），导致更新方程不稳定。可以尝试增加损耗因子γ。
- 没有观察到负折射：检查入射波频率是否落在设计的DNG频段内。检查材料参数实现是否正确，特别是辅助微分方程的差分格式是否稳定。确保网格分辨率足够高，能分辨DNG材料中的波（在DNG材料中，波长可能非常短）。
- 边界反射严重：检查PML设置，确保层数足够，电导率剖面光滑。检查源是否离PML太近，一般源与PML之间至少应留有几个波长的距离。
通过这样一个完整的案例，你不仅能验证代码的正确性，更能深刻理解DNG材料的奇异电磁特性以及FDTD方法在模拟复杂电磁问题中的强大能力。这个“DNG.zip”里的代码，正是通往这一理解过程的宝贵工具。
本文还有配套的精品资源，点击获取
