一维FDTD模拟平面波

1、上海大学 20 20 学年 秋 季学期电磁场理论A仿真报告题 目： 一维FDTD模拟平面波 任课教师： XXX 学 号： XXXXXXXX 姓 名： XXX 日 期： 2019 成 绩： 一维FDTD模拟平面波1、 仿真的目的和意义自1873年Maxwen发表电磁学通论以来，人们在研究电磁频谱、拓展电磁频段、开发电磁资源以及应用研究方面作出了大量的工作。其中具有代表性的有1888年Hertz首次通过实验验证了Maxwen的电磁理论，证明了无线电辐射具有波的所有特性，并记录下了观测到的无线电波波长。在微波工程中，人们常常以麦克斯韦方程及其解作为开始。然而，正是这些方程引起了数学上的复杂性，因为麦

2、克斯韦方程包含了作为空间坐标函数的矢量微分或积分的运算，而这些计算常常又很复杂，所以一个很好的计算方法尤其重要。常用的计算方法有数值积分法、有限差分法（FEM）、有限元法（FDM）、矩量法（MOM）、边界元法等方法，但是这些方法仍然是基于数值计算方法的，其内容在结合电磁场理论、计算机科学和计算数学方面显得不足。计算电磁学是20世纪90年代在电磁学、计算数学和计算机科学基础上产生的新的基于近代计算机平台的交叉科学，创造了对世界模拟和理解的新的途径，创造了新的思维和实验方法。FDTD算法是以Yee元胞为空间电磁场离散单元，将麦克斯韦旋度方程转化为差分方程，表述简明，容易理解，结合计算机技术能耐处理

3、十分复杂的电磁问题；在时间轴上逐步推进求解，有很好的稳定性和收敛性，因而在工程电磁学各个领域备受重视。2、 仿真过程方法理论介绍：时域有限差分方法 （Finite-Difference Time-Domain Method，FDTD）由K.S.Yee于1966年首次提出。他是利用这种计算方法将麦克斯韦方程的中的磁通连续方程和电磁感应定律在计算机上模拟得出，从而可以求得电场和磁场的值。FDTD方法是求麦克斯韦微分方程的直接时域的方法，在计算中将空间某一样本点的电场（或磁场）与周围格点的磁场（或电场）直接相关联，求介质参数已赋值给空间每一个元胞，因此这个方法可以处理复杂形状目标和非均匀介质物体的电

4、磁散射、辐射等问题。同时，FDTD随时间推进可以方便地给出电磁场的实践演化过程，在计算机上以伪色彩方式显示，这种电磁场可视化结果清楚地显示了物理过程，便于分析和设计。（1）（2）其中E为电场强度，单位为伏特/米（V/m）； D为电通量密度，单位为库仑/平方米（C/）； H为磁场强度，单位为安培/米（A/m）； B为磁通量密度，单位为韦伯/平方米（Wb/）； J为电流密度，单位为安培/平方米（A/）；本构关系： FDTD法对电磁场E、H分量在空间和时间上采取交替抽样的离散方式，每一个E（或H ）场分量周围有四个H（或E ）场分量环绕，应用这种离散方式将含时间变量的麦克斯韦旋度方程转化为一组差分方

5、程。（3）（4）如上述方程（3）或（4），上n所对应的为第n个时刻，而x所对应的为第x个网格。图1如图1所示，如果要得到第n个时刻第k个网格的磁场就要用第n-1个时刻第K和第k+1个格子的电场值对距离的微分来的到。然后利用第n个时刻的磁场值就可以得到第n+1个时刻第k+1个格子的电场值3、 吸收边界条件：由于电磁场是具有反射特性的，在对边界的网格计算就不能一味的继续利用式（3）、（4）进行计算，而是要补充吸收边界条件。边界条件的补充可以有不同的方法，比如利用电磁波的传播。4、FDTD仿真结果示意图：图a图b采用matlab软件进行一维FDTD的模拟平面波仿真注：蓝色为电场，红色为磁场经计算：速

6、度（光速）：=2.9979108 m/s波阻抗：=376.75、分析与总结FDTD的提出，不得不说是对电磁界的又一个突破，使得纸面上的理论同高速、准确的计算机科学相联系起来，完成了电磁在现实中的模拟。在刚开始收集资料的过程中由于自己对于FDTD以及相关概念并不是多了解，导致浪费了很多的精力和时间，后面我发现应该采取的方法是多去实践，不要总是在那里纠结于理论的东西，要通过实践来验证理论并加深对理论的理解。同时，程序出来后也不是一开始就会很顺利的运行，出现了错误要仔细分析其中的原因，并不断去调试错误直到可以完全运行。总之，通过此次仿真自己不仅加深了FDTD对平面波传播的理解，同时也加强了对编程的掌

7、握。6、参考文献1蒋祺. FDTD算法的若干关键技术研究D.上海交通大学,2015.2蔡明娟. 时域积分方程在分析介质问题中的算法研究与应用D.国防科学技术大学,2006.附录：程序源代码clcclear% 定义初始常数eps_0 = 8.e-12; % 真空介电常数mu_0 = 4*pi*1e-7; % 自由空间磁导率c = 1/sqrt(mu_0*eps_0); % 光速% 定义问题几何和参数domain_size = 1; % 一维问题空间长度dx = 1e-3; % 长度单位为米dt = 3e-12; % 时间步长以秒为单位number_of_time_steps = 2000; %

8、迭代次数nx = round(domain_size/dx); % 一维问题空间中的单元数source_position = 0.5; % 电流源Jz的位置% 初始化Ceze = zeros(nx+1, 1);Cezhy = zeros(nx+1, 1);Cezj = zeros(nx+1, 1);Ez = zeros(nx+1, 1);eps_r_z = ones (nx+1, 1); % 自由空间sigma_e_z = zeros(nx+1, 1); % 自由空间Chyh = zeros(nx, 1);Chyez = zeros(nx, 1);Chym = zeros(nx, 1);Hy

9、= zeros(nx, 1);My = zeros(nx, 1);mu_r_y = ones (nx, 1); % 自由空间sigma_m_y = zeros(nx, 1); % 自由空间 % 计算FDTD的更新参数Ceze = (2 * eps_r_z * eps_0 - dt * sigma_e_z) . ./(2 * eps_r_z * eps_0 + dt * sigma_e_z); Cezhy = (2 * dt / dx) . ./(2 * eps_r_z * eps_0 + dt * sigma_e_z);Chyh = (2 * mu_r_y * mu_0 - dt * sigm

10、a_m_y) . ./(2 * mu_r_y * mu_0 + dt * sigma_m_y); Chyez = (2 * dt / dx) . ./(2 * mu_r_y * mu_0 + dt * sigma_m_y);% 定义高斯源波形time = dt * 0:number_of_time_steps-1.;Ez_waveform = exp(-(time-2e-10)/5e-11).2)*1e-3/dx;source_position_index = round(nx * source_position/domain_size)+1;Ez_positions = 0:nx*dx;Hy

11、_positions = (0:nx-1+0.5)*dx;% 初始化画图程序axis(0 1 -1 1);lez = line(Ez_positions, Ez*0, Ez, Color, b, linewidth, 1.5);lhy = line(Hy_positions, 377*Hy, Hy*0, Color, r, LineWidth, 1.5, LineStyle,-.); set(gca, fontsize, 12, fontweight, bold);set(gcf,Color,white);axis square;xlabel(x m);ylabel(E V/m);legend

12、(Electric field, Magnetic field, location, northeast);grid on; % FDTD 循环for time_step = 1:number_of_time_steps % 更新Ez Ez(source_position_index) = Ez_waveform(time_step); % 更新磁场 Hy(1:nx) = Chyh(1:nx) .* Hy(1:nx) . + Chyez(1:nx) .* (Ez(2:nx+1) - Ez(1:nx); % 更新电场 Ez(2:nx) = Ceze (2:nx) .* Ez(2:nx) . + Cezhy(2:nx) .* (Hy(2:nx) - Hy(1:nx-1); Ez(1) = 0; % 在x = 0 m处应用PEC边界条件 Ez(nx+1) = 0; % 在x = 1 m处应用PEC边界条件 % 绘制一维瞬态场 delete(lez); delete(lhy); lez = line(Ez_positions, Ez, Color, b,