couple line_fdtd.m

FDTD法计算平行耦合双微带线电容电感参数的程序

% The main function
clear all;
Timestart=datestr(now)
% 定义仿真区，网格大小，电参数等
% [dx,dy,dz,dt,Nt,W,S,H,TT,Ls,Nx,Ny,Nz,t0,T,NH,eps,u,segema]=Boxsize(0.05e-3,0.05e-3,0.05e-3,0.05e-12,1600,0.3e-3,0.3e-3,0.25e-3,0.05e-3,2e-3,2.5e-3,0.85e-3,6e-3,4.5,1,0);
dx=0.1e-3; dy=0.1e-3; dz=0.1e-3;  dt=0.05e-12; % step length
Lx=4e-3;Ly=1.2e-3; Lz=12e-3; H=0.2e-3;
Nx=Lx/dx+1; Ny=round(Ly/dy)+1; Nz=Lz/dz+1;NH=H/dy+1;  Nt=1600;  % 空间时间网格数
W=0.3e-3; S=0.3e-3;  TT=0.1e-3;  Er=4.5; Ls=4e-3;  % 结构参数
t0=500*dt; T=60*dt;  % parameters of Gauss inpulse 
u0=pi*4e-7; ur=1; %  真空磁导率
u(1:Nx,1:Ny,1:Nz)=u0*ur; 
segema0=0;
segema(1:Nx,1:Ny,1:Nz)=segema0;
eps0=8.854e-12;  %真空介电常数
epsr(1:Nx,NH+1:Ny,1:Nz)=1;
epsr(1:Nx,1:NH,1:Nz)=Er;
eps=eps0*epsr;
                                                                
% 差分方程跌代系数
[Exee,Exeh,Eyee,Eyeh,Ezee,Ezeh]=Coefficient_E(eps,segema,dt,Nx,Ny,Nz);
[Hxhe,Hxhh,Hyhe,Hyhh,Hzhe,Hzhh]=Coefficient_H(u,segema,dt,Nx,Ny,Nz);
% 附初值
Ex(1:Nx,1:Ny,1:Nz)=0.0;  Ey(1:Nx,1:Ny,1:Nz)=0.0;  Ez(1:Nx,1:Ny,1:Nz)=0.0; 
Hx(1:Nx,1:Ny,1:Nz)=0.0;  Hy(1:Nx,1:Ny,1:Nz)=0.0;  Hz(1:Nx,1:Ny,1:Nz)=0.0;

% 导体边界
[Exc,Ezc]=PEC(Nx,Ny,Nz,dx,dy,dz,W,S,TT,NH);

% 激励源位置
[Sxsta,Sxend,Systa,Syend,Szsta,Szend]=Source(round((Nx-(2*W+S)/dx)/2),round((Nx-S/dx)/2),1,NH,Ls/dz,Ls/dz); 

% 观察电压的位置
[Oxsta1,Oxend1,Oysta1,Oyend1,Ozsta1,Ozend1]=Observe(round((Nx-S/dx-W/dx)/2),round((Nx-S/dx-W/dx)/2),1,NH-1,Szsta+20,Szsta+20);
[Oxsta2,Oxend2,Oysta2,Oyend2,Ozsta2,Ozend2]=Observe(round((Nx+S/dx+W/dx)/2),round((Nx+S/dx+W/dx)/2),1,NH-1,Szsta+20,Szsta+20);

for n=1:Nt
  % 保存前时刻的场值，备Mur吸收边界所用
[preExys1,preExys2,preExye1,preExye2,preExzs1,preExzs2,preExze1,preExze2]=PreviousEx(Ex,Ey,Ez,Nx,Ny,Nz);
[preEyxs1,preEyxs2,preEyxe1,preEyxe2,preEyzs1,preEyzs2,preEyze1,preEyze2]=PreviousEy(Ex,Ey,Ez,Nx,Ny,Nz);
[preEzxs1,preEzxs2,preEzxe1,preEzxe2,preEzys1,preEzys2,preEzye1,preEzye2]=PreviousEz(Ex,Ey,Ez,Nx,Ny,Nz);
 
%   以下是电场的迭代 
 Ex(1:Nx,2:Ny-1,2:Nz-1)=Exee(1:Nx,2:Ny-1,2:Nz-1).* Ex(1:Nx,2:Ny-1,2:Nz-1)+Exeh(1:Nx,2:Ny-1,2:Nz-1).*((Hz(1:Nx,2:Ny-1,2:Nz-1)-Hz(1:Nx,1:Ny-2,2:Nz-1))/dy-(Hy(1:Nx,2:Ny-1,2:Nz-1)-Hy(1:Nx,2:Ny-1,1:Nz-2))/dz);
 Ey(2:Nx-1,1:Ny,2:Nz-1)=Eyee(2:Nx-1,1:Ny,2:Nz-1).* Ey(2:Nx-1,1:Ny,2:Nz-1)+Eyeh(2:Nx-1,1:Ny,2:Nz-1).*((Hx(2:Nx-1,1:Ny,2:Nz-1)-Hx(2:Nx-1,1:Ny,1:Nz-2))/dz-(Hz(2:Nx-1,1:Ny,2:Nz-1)-Hz(1:Nx-2,1:Ny,2:Nz-1))/dx);   
 Ez(2:Nx-1,2:Ny-1,1:Nz)=Ezee(2:Nx-1,2:Ny-1,1:Nz).* Ez(2:Nx-1,2:Ny-1,1:Nz)+Ezeh(2:Nx-1,2:Ny-1,1:Nz).*((Hy(2:Nx-1,2:Ny-1,1:Nz)-Hy(1:Nx-2,2:Ny-1,1:Nz))/dx-(Hx(2:Nx-1,2:Ny-1,1:Nz)-Hx(2:Nx-1,1:Ny-2,1:Nz))/dy); 
 % 激励信号
 Ey(Sxsta:Sxend,Systa:Syend,Szsta:Szend)=-exp(-((n-1)*dt-t0)^2/T^2)/((Syend-Systa)*dy); 
%  Ey(Sxsta:Sxend,Systa:Syend,Szsta:Szend)=exp(-((n-1)*dt-t0)^2/T^2); 

%*****  first order  Mur absorb boundary condition   
c=3e8;
%% for x=0
Ey(1,:,:)=preEyxs2(1,:,:)+(c*dt-dx)/(c*dt+dx)*(Ey(2,:,:)-preEyxs1(1,:,:));
Ez(1,:,:)=preEzxs2(1,:,:)+(c*dt-dx)/(c*dt+dx)*(Ez(2,:,:)-preEzxs1(1,:,:));
%% for x=Lx
Ey(Nx,:,:)=preEyxe2(1,:,:)+(c*dt-dx)/(c*dt+dx)*(Ey(Nx-1,:,:)-preEyxe1(1,:,:));
Ez(Nx,:,:)=preEzxe2(1,:,:)+(c*dt-dx)/(c*dt+dx)*(Ez(Nx-1,:,:)-preEzxe1(1,:,:));
%% for y=0
Ex(:,1,:)=preExys2(:,1,:)+(c*dt-dy)/(c*dt+dy)*(Ex(:,2,:)-preExys1(:,1,:));
Ez(:,1,:)=preEzys2(:,1,:)+(c*dt-dy)/(c*dt+dy)*(Ez(:,2,:)-preEzys1(:,1,:));
%% for y=Ly
Ex(:,Ny,:)=preExye2(:,1,:)+(c*dt-dy)/(c*dt+dy)*(Ex(:,Ny-1,:)-preExye1(:,1,:));
Ez(:,Ny,:)=preEzye2(:,1,:)+(c*dt-dy)/(c*dt+dy)*(Ez(:,Ny-1,:)-preEzye1(:,1,:));
%% for z=0
Ex(:,:,1)=preExzs2(:,:,1)+(c*dt-dz)/(c*dt+dz)*(Ex(:,:,2)-preExzs1(:,:,1));
Ey(:,:,1)=preEyzs2(:,:,1)+(c*dt-dz)/(c*dt+dz)*(Ey(:,:,2)-preEyzs1(:,:,1));
%% for z=Lz
Ex(:,:,Nz)=preExze2(:,:,1)+(c*dt-dz)/(c*dt+dz)*(Ex(:,:,Nz-1)-preExze1(:,:,1));
Ey(:,:,Nz)=preEyze2(:,:,1)+(c*dt-dz)/(c*dt+dz)*(Ey(:,:,Nz-1)-preEyze1(:,:,1));

 %  导体边界条件
Ex=Ex.*Exc;
Ez=Ez.*Ezc; 
 
%   以下为磁场迭代
Hx(1:Nx,1:Ny-1,1:Nz-1)=Hxhh(1:Nx,1:Ny-1,1:Nz-1).*Hx(1:Nx,1:Ny-1,1:Nz-1)-Hxhe(1:Nx,1:Ny-1,1:Nz-1).*((Ez(1:Nx,2:Ny,1:Nz-1)-Ez(1:Nx,1:Ny-1,1:Nz-1))/dy-(Ey(1:Nx,1:Ny-1,2:Nz)-Ey(1:Nx,1:Ny-1,1:Nz-1))/dz);
Hy(1:Nx-1,1:Ny,1:Nz-1)=Hyhh(1:Nx-1,1:Ny,1:Nz-1).*Hy(1:Nx-1,1:Ny,1:Nz-1)-Hyhe(1:Nx-1,1:Ny,1:Nz-1).*((Ex(1:Nx-1,1:Ny,2:Nz)-Ex(1:Nx-1,1:Ny,1:Nz-1))/dz-(Ez(2:Nx,1:Ny,1:Nz-1)-Ez(1:Nx-1,1:Ny,1:Nz-1))/dx);
Hz(1:Nx-1,1:Ny-1,1:Nz)=Hzhh(1:Nx-1,1:Ny-1,1:Nz).*Hz(1:Nx-1,1:Ny-1,1:Nz)-Hzhe(1:Nx-1,1:Ny-1,1:Nz).*((Ey(2:Nx,1:Ny-1,1:Nz)-Ey(1:Nx-1,1:Ny-1,1:Nz))/dx-(Ex(1:Nx-1,2:Ny,1:Nz)-Ex(1:Nx-1,1:Ny-1,1:Nz))/dy);

% 记录电压
 Vs(n)=-sum(Ey(Sxsta,Systa:Syend-1,Szsta)*dy);
V11(n)=-sum(Ey(Oxsta1,Oysta1:Oyend1,Ozsta1)*dy);
V12(n)=-sum(Ey(Oxsta1,Oysta1:Oyend1,Ozsta1+10)*dy);
V13(n)=-sum(Ey(Oxsta1,Oysta1:Oyend1,Ozsta1+20)*dy);
V14(n)=-sum(Ey(Oxsta1,Oysta1:Oyend1,Ozsta1+30)*dy);
V15(n)=-sum(Ey(Oxsta1,Oysta1:Oyend1,Ozsta1+40)*dy);
V21(n)=-sum(Ey(Oxsta2,Oysta2:Oyend2,Ozsta2)*dy);
V22(n)=-sum(Ey(Oxsta2,Oysta2:Oyend2,Ozsta2+10)*dy);
V23(n)=-sum(Ey(Oxsta2,Oysta2:Oyend2,Ozsta2+20)*dy);
V24(n)=-sum(Ey(Oxsta2,Oysta2:Oyend2,Ozsta2+30)*dy);
V25(n)=-sum(Ey(Oxsta2,Oysta2:Oyend2,Ozsta2+40)*dy);
n
Ey(Sxsta,Systa,Szsta)
Ey(Oxsta1,round((Oysta1+Oyend1)/2),Ozsta1)
end

clear n;
n=200:1:1600;
figure(1);
plot(n,Vs(200:1600));
figure(2);
plot(n,V11(200:1600));
hold on;
plot(n,V12(200:1600));
plot(n,V13(200:1600));
plot(n,V14(200:1600));
plot(n,V15(200:1600));
plot(n,V21(200:1600),'-k');
plot(n,V22(200:1600),'-k');
plot(n,V23(200:1600),'-k');
plot(n,V24(200:1600),'-k');
plot(n,V25(200:1600),'-k');
Timeend=datestr(now)
Timecosted=Timeend-Timestart;
Timecosted