1856阅读
3回复

[转载]X波段标准波导BJ100主模式FDTD仿真求解程序

离线tingo

临渊羡鱼，不如退而结网

UID ：28673

注册：2009-03-29
登录：2019-09-21
发帖：126
等级：仿真三级

0楼发表于: 2011-10-25 15:33:03

最近初学FDTD的人挺多的，经常遇到问问题的朋友，在此转载一个自己在微网分享的二维FDTD仿真求解X波段标准波导BJ100主模程序，希望对大家有所帮助。
其中包含的内容有：FDTD仿真TE波，调制的高斯信号，离散傅里叶变换。

[code]% FDTD 仿真求解X波段标准波导BJ-100主模式
%编写者：何国强；
%Email:guoqianghe22@gmail.com;
%编写时间：2011年10月25日；
clear;
clc;
%% 计算中用到的一些电磁常量；
mu=4*pi*(1e-7); %磁导率；
epsilon=(1/36/pi)*(1e-9); %介电常数；
Vvacuum=1/sqrt(mu*epsilon); % 真空中的波传播速度；
%% 网格数据；
a=22.86e-3;  %波导宽边；
b=10.16e-3;  %波导窄边；
nx=100;
ny=50;
dx=a/nx;
dy=b/ny;

dt=1/Vvacuum/sqrt(1/dx/dx+1/dy/dy);

io=round(nx/4); %观察点位置
jo=round(ny/2); %观察点位置
%% 激励源；
fc=6.5e9; % GHz;
fmax=8e9; %  GHz;
T=1/2/fmax; % 激励源周期；
tau=3*T;
ntau=tau/dt;
ndelay=3*ntau;
omega=fc*2*pi; %角频率
itmax=round(20*tau/dt); %最大时间步进数目；
t=(1:itmax)*dt; %时域时间采样点；
source=zeros(1,itmax);
%调制的高斯信号；
for n=1:round(6*ntau)+1;
    source(n)=sin(omega*(n-ndelay)*dt)*exp(-((n-ndelay)^2/ntau^2));
end
is=round(nx/2);
%
%% 内存分配；
% EM场值
ex=zeros(nx,ny+1);
ey=zeros(nx+1,ny);
hz=zeros(nx,ny);
%系数变量；
caex=zeros(nx,ny+1);
cbex=zeros(nx,ny+1);
caey=zeros(nx+1,ny);
cbey=zeros(nx+1,ny);
cahz=zeros(nx,ny);
cbhz=zeros(nx,ny);
eyobsv=zeros(1,itmax);
%% 波导中的介质参数；
sigma=0;
rou=0;
eaf=sigma*dt/2/epsilon;
ca=(1-eaf)/(1+eaf);
cb=dt/epsilon/(1+eaf);
haf=rou*dt/2/mu;
cp=(1-haf)/(1+haf);
cq=dt/mu/(1+haf);
caex=ca*ones(size(caex));
cbex=cb*ones(size(cbex));
caey=ca*ones(size(caey));
cbey=cb*ones(size(cbey));
cahz=cp*ones(size(cahz));
cbhz=cq*ones(size(cbhz));
% PEC边界处的ex,ey计算系数；
caex(:,1)=1;
cbex(:,1)=0;
caex(:,ny+1)=1;
cbex(:,ny+1)=0;
caey(1,:)=1;
cbey(1,:)=0;
caey(nx+1,:)=1;
cbey(nx+1,:)=0;
%% FDTD Looping
for it=1:itmax
    %% Update Hz；
    hz(1:nx,1:ny)=cahz(1:nx,1:ny).*hz(1:nx,1:ny)-...
        cbhz(1:nx,1:ny).*( ( ey(2:nx+1,:)-ey(1:nx,:) )/dx - ( ex(:,2:ny+1)-ex(:,1:ny) )/dy );
    %% Update Ex , Ey;
    ex(:,2:ny)=caex(:,2:ny).*ex(:,2:ny)+...
        cbex(:,2:ny).*( hz(:,2:ny)-hz(:,1:ny-1) )/dy;
    ey(2:nx,:)=caey(2:nx,:).*ey(2:nx,:)-...
        cbey(2:nx,:).*( hz(2:nx,:)-hz(1:nx-1,:) )/dx;
    ey(is,:)=ey(is,:)+source(it);  % 加载激励源;
    eyobsv(it)=ey(io,jo);    % 记录观察点的场值；
end
%% 观察点场值做离散傅里叶变换（时域到频域）
nf=1000;  %频率采样点数目；
fmin=5e9; % 5GHz;
df=(fmax-fmin)/(nf-1); %频率采样步长；
Ef=zeros(1,nf); %为频域场量分配内存；
J=sqrt(-1); %虚数i;
for n=1:nf
    coe=exp(-J*2*pi*( fmin+(n-1)*df ).*t );
    Ef(n)=dt*sum(eyobsv(1,:).*coe );
end
freq=fmin:df:fmax;
Ef=abs(Ef);
[Efmax,idmax]=max(Ef);
%% 输出计算结果
% figure1 Ey在观察点的瞬态响应
figure(1)
plot(t,eyobsv);
xlabel('时间(s)')
ylabel('电场(V/m)')
title('观察点Ey瞬态响应')

% figure2 Ey的频域响应
figure(2)
plot(freq,Ef,'c');
hold on;
plot(freq(idmax),Ef(idmax),'r^');
xlabel('频率(Hz)');
ylabel('电场（V/m）')
title('电场的频域值')
legend('Ey频域值','主模')
% 输出主模式TE10模截止频率的理论值和仿真值；
Theorynamdac=4.572e-2;
Theoryfc=Vvacuum./Theorynamdac;
Simfc=freq(idmax);
disp('=*******************************************************=');
disp('** ..

未注册仅能浏览部分内容，查看全部内容及附件请先登录或注册

离线tingo

临渊羡鱼，不如退而结网

UID ：28673

注册：2009-03-29
登录：2019-09-21
发帖：126
等级：仿真三级

1楼发表于: 2011-10-25 20:17:49

这比较适合于初学FDTD者做参考，有问题大家可以在此一起讨论！

仿真新人！

离线lx84

呵呵

UID ：30566

注册：2009-04-22
登录：2020-09-16
发帖：866
等级：积极交流六级

2楼发表于: 2011-10-26 23:20:23

前排支持楼主

继续！

离线tingo

临渊羡鱼，不如退而结网

UID ：28673

注册：2009-03-29
登录：2019-09-21
发帖：126
等级：仿真三级

3楼发表于: 2011-10-27 22:22:27

回 2楼(lx84) 的帖子

谢谢支持！

仿真新人！

发帖回复