一维电磁波的有限差分时域 (FTDT) 方法答案

【问题标题】：Finite Difference Time Domain (FTDT) method for 1D EM Wave一维电磁波的有限差分时域 (FTDT) 方法
【发布时间】：2014-06-25 12:45:09
【问题描述】：

我尝试编写代码来求解以下耦合偏微分电磁波方程：

代码采用有限差分时域，使用 Yee 算法，可以在以下两个在线文档中阅读：

http://www.eecs.wsu.edu/~schneidj/ufdtd/ufdtd.pdf

http://www.eecs.wsu.edu/~schneidj/ufdtd/chap3.pdf

我希望左侧边界处的源是一个正弦波，其参数如下：

Ex(0,t) = E0 sin(2πft)      for 0 ≤ t ≤ 1/f

代码在每个循环中更新波的电场和磁场属性。

我的初始代码如下：

%FDTD Yee algorithm to solve coupled EM wave equations
clear
clc


G=50;           %Specify size of the grid
f=10^3;         %choose initial frequency of wave in hertz
e=1;            %specify permitivity and permeability (normalised condition)
u=1;            
Nt=150;         %time steps
E0=1;           %Electric Field initial amplitude

%Specify step sizes using corruant condition
c=3*10^8;
dx=0.01;
dt=dx/2*c;

%make constant terms
c1=-dt/(dx*e);
c2=-dt/(dx*u);

%create vgector place holders
Ex=zeros(1,G);
Hy=zeros(1,G);

%create updating loop
M=moviein(Nt);
for t=1:Nt

    % Spatial Ex

    for k=2:G-1
        Ex(k)=Ex(k)+c1*(Hy(k)-Hy(k-1));
    end
    Ex(G)=0; %PEC boundary condition
    %E Source at LHS boundary 
    Ex(1)=E0*sin(2*pi*f*t);
    %Spatial Hy
    for n=1:G-1
        Hy(n)=Hy(n)+c2*(Ex(n)-Ex(n+1));
    end
    Hy(G)=0; %PMC boundary condition

plot(Ex);
M(:,t) = getframe;
end
movie(M,1);

基本上我想创建一个更新电影，显示传播到右侧边界的正弦波，编码为完美的电导体；因此反射波，然后传播回左侧边界，编码为完美绝缘体；吸收波。

我遇到的问题如下：

1) 我不确定如何正确实现我想要的源。它似乎不是纯粹的正弦曲线。

2) 我编码的波开始传播，但在模拟的大部分时间里它很快就消失了。我不知道为什么会发生这种情况

3) 我对运行电影模拟知之甚少，并且在运行解决方案时情节会摇摆不定。我怎样才能阻止这种情况？

【问题讨论】：

标签： matlab numerical-methods

【解决方案1】：

你的波衰减是因为微分方程没有正确实现；而是：

Ex(k)=Ex(k)+c1*(Hy(k)-Hy(k-1));

你应该使用

Ex1(k)=Ex(k)+c1*(Hy(k)-Hy(k-1));

而不是：

Hy(n)=Hy(n)+c2*(Ex(n)-Ex(n+1));

你应该使用：

Hy1(n)=Hy(n)+c2*(Ex(n)-Ex(n+1));

并且，在循环结束时更新当前的“数据框”：

Hy = Hy1;
Ey = Ey1;

（您还应该注意边界条件）。

如果您想要一个在数据更改时不会更改的固定图框，请首先创建一个axis，您可以在其中进行绘图，并定义 xmin/max 和 ymin/max，请参阅http://www.mathworks.com/help/matlab/ref/axis.html

【讨论】：

【解决方案2】：

您应该将 Courant 数设置为更接近 1，例如 0.995。因此 delta_t = 0.995*delta_x/c。还假设 delta_x 是公制单位，那么 e 和 u 应该是公制单位。我不知道使用的具体编码语言，但在 c 或 c++ 中不需要中间变量 Ey1 等。此外，每个波长至少应有 10 个样本以确保准确性（最好是 60 个）。因此波长 = 60*delta_x，因此频率大致等于 10 的 9 次方。另外，我认为正弦源应该是 E0 * sin(2* pi * f * t * delta_t)。你需要调整你的常数，然后再试一次

【讨论】：