⭐ 欢迎来到虫虫下载站！ | 📦 资源下载 📁 资源专辑 ℹ️ 关于我们
⭐ 虫虫下载站
📤 上传资源

📄 coupleline1stmur.m

📁 FDTD法计算平行耦合双微带线电容电感参数的程序
💻 M
字号:
🔍 
% 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.05e-3; dy=0.05e-3; dz=0.05e-3;  dt=0.05e-12; % step length
Lx=4e-3;Ly=1.15e-3; Lz=12e-3; H=0.25e-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.05e-3;  Er=4.5; Ls=4e-3;  % 结构参数
t0=500*dt; T=100*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+40,Szsta+40);
[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+40,Szsta+40);
% 观察电流的位置
[OIxsta1,OIxend1,OIysta1,OIyend1,OIzsta1,OIzend1]=Observe(round((Nx-S/dx)/2-W/dx-2),round((Nx-S/dx)/2+2),NH-2,NH+2,Szsta+40,Szsta+40);
[OIxsta2,OIxend2,OIysta2,OIyend2,OIzsta2,OIzend2]=Observe(round((Nx+S/dx)/2-2),round((Nx+S/dx)/2+W/dx+2),NH-2,NH+2,Szsta+40,Szsta+40);

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+20)*dy);
V13(n)=-sum(Ey(Oxsta1,Oysta1:Oyend1,Ozsta1+40)*dy);
V14(n)=-sum(Ey(Oxsta1,Oysta1:Oyend1,Ozsta1+60)*dy);
V15(n)=-sum(Ey(Oxsta1,Oysta1:Oyend1,Ozsta1+80)*dy);
V21(n)=-sum(Ey(Oxsta2,Oysta2:Oyend2,Ozsta2)*dy);
V22(n)=-sum(Ey(Oxsta2,Oysta2:Oyend2,Ozsta2+20)*dy);
V23(n)=-sum(Ey(Oxsta2,Oysta2:Oyend2,Ozsta2+40)*dy);
V24(n)=-sum(Ey(Oxsta2,Oysta2:Oyend2,Ozsta2+60)*dy);
V25(n)=-sum(Ey(Oxsta2,Oysta2:Oyend2,Ozsta2+80)*dy);

% 记录电流
% 环路积分
I11(n)=sum(Hx(OIxsta1:OIxend1,Oysta1,Ozsta1)*dy);
I11(n)=I11(n)+sum(Hx(OIxsta1:OIxend1,Oyend1,Ozsta1)*dy);
I11(n)=I11(n)+sum(Hy(OIxsta1,Oysta1:Oyend1,Ozsta1)*dy);
I11(n)=I11(n)+sum(Hy(OIxend1,Oysta1:Oyend1,Ozsta1)*dy);

I12(n)=sum(Hx(OIxsta1:OIxend1,Oysta1,Ozsta1+20)*dy);
I12(n)=I12(n)+sum(Hx(OIxsta1:OIxend1,Oyend1,Ozsta1+20)*dy);
I12(n)=I12(n)+sum(Hy(OIxsta1,Oysta1:Oyend1,Ozsta1+20)*dy);
I12(n)=I12(n)+sum(Hy(OIxend1,Oysta1:Oyend1,Ozsta1+20)*dy);

I13(n)=sum(Hx(OIxsta1:OIxend1,Oysta1,Ozsta1+40)*dy);
I13(n)=I13(n)+sum(Hx(OIxsta1:OIxend1,Oyend1,Ozsta1+40)*dy);
I13(n)=I13(n)+sum(Hy(OIxsta1,Oysta1:Oyend1,Ozsta1+40)*dy);
I13(n)=I13(n)+sum(Hy(OIxend1,Oysta1:Oyend1,Ozsta1+40)*dy);

I14(n)=sum(Hx(OIxsta1:OIxend1,Oysta1,Ozsta1+60)*dy);
I14(n)=I14(n)+sum(Hx(OIxsta1:OIxend1,Oyend1,Ozsta1+60)*dy);
I14(n)=I14(n)+sum(Hy(OIxsta1,Oysta1:Oyend1,Ozsta1+60)*dy);
I14(n)=I14(n)+sum(Hy(OIxend1,Oysta1:Oyend1,Ozsta1+60)*dy);

I15(n)=sum(Hx(OIxsta1:OIxend1,Oysta1,Ozsta1+80)*dy);
I15(n)=I15(n)+sum(Hx(OIxsta1:OIxend1,Oyend1,Ozsta1+80)*dy);
I15(n)=I15(n)+sum(Hy(OIxsta1,Oysta1:Oyend1,Ozsta1+80)*dy);
I15(n)=I15(n)+sum(Hy(OIxend1,Oysta1:Oyend1,Ozsta1+80)*dy);

I21(n)=sum(Hx(OIxsta2:OIxend2,Oysta2,Ozsta2)*dy);
I21(n)=I21(n)+sum(Hx(OIxsta2:OIxend2,Oyend2,Ozsta2)*dy);
I21(n)=I21(n)+sum(Hy(OIxsta2,Oysta2:Oyend2,Ozsta2)*dy);
I21(n)=I21(n)+sum(Hy(OIxend2,Oysta2:Oyend2,Ozsta2)*dy);

I22(n)=sum(Hx(OIxsta2:OIxend2,Oysta2,Ozsta2+20)*dy);
I22(n)=I22(n)+sum(Hx(OIxsta2:OIxend2,Oyend2,Ozsta2+20)*dy);
I22(n)=I22(n)+sum(Hy(OIxsta2,Oysta2:Oyend2,Ozsta2+20)*dy);
I22(n)=I22(n)+sum(Hy(OIxend2,Oysta2:Oyend2,Ozsta2+20)*dy);

I23(n)=sum(Hx(OIxsta2:OIxend2,Oysta2,Ozsta2+40)*dy);
I23(n)=I23(n)+sum(Hx(OIxsta2:OIxend2,Oyend2,Ozsta2+40)*dy);
I23(n)=I23(n)+sum(Hy(OIxsta2,Oysta2:Oyend2,Ozsta2+40)*dy);
I23(n)=I23(n)+sum(Hy(OIxend2,Oysta2:Oyend2,Ozsta2+40)*dy);

I24(n)=sum(Hx(OIxsta2:OIxend2,Oysta2,Ozsta2+60)*dy);
I24(n)=I24(n)+sum(Hx(OIxsta2:OIxend2,Oyend2,Ozsta2+60)*dy);
I24(n)=I24(n)+sum(Hy(OIxsta2,Oysta2:Oyend2,Ozsta2+60)*dy);
I24(n)=I24(n)+sum(Hy(OIxend2,Oysta2:Oyend2,Ozsta2+60)*dy);

I25(n)=sum(Hx(OIxsta2:OIxend2,Oysta2,Ozsta2+80)*dy);
I25(n)=I25(n)+sum(Hx(OIxsta2:OIxend2,Oyend2,Ozsta2+80)*dy);
I25(n)=I25(n)+sum(Hy(OIxsta2,Oysta2:Oyend2,Ozsta2+80)*dy);
I25(n)=I25(n)+sum(Hy(OIxend2,Oysta2:Oyend2,Ozsta2+80)*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');

figure(3);
plot(n,I11(200:1600));
hold on;
plot(n,I12(200:1600));
plot(n,I13(200:1600));
plot(n,I14(200:1600));
plot(n,I15(200:1600));
plot(n,I21(200:1600),'-k');
plot(n,I22(200:1600),'-k');
plot(n,I23(200:1600),'-k');
plot(n,I24(200:1600),'-k');
plot(n,I25(200:1600),'-k');

% 下面进行傅立叶变换计算分布电容和电感参数
Fmax=25e9;
n=1:1*Nt;
k=1;
  Nf=15;
  df=Fmax/Nf;
  while k<=Nf
      FI11(k)=sum(I11.*exp(-i*2*pi*(k*df)*(dt*n)))*dt;
      FI15(k)=sum(I15.*exp(-i*2*pi*(k*df)*(dt*n)))*dt;
      FV13(k)=sum(V13.*exp(-i*2*pi*(k*df)*(dt*n)))*dt;
%       FV15(k)=sum(V15.*exp(-i*2*pi*(k*df)*(dt*n)))*dt;
      FI21(k)=sum(I21.*exp(-i*2*pi*(k*df)*(dt*n)))*dt;
      FI25(k)=sum(I25.*exp(-i*2*pi*(k*df)*(dt*n)))*dt;
      FV23(k)=sum(V23.*exp(-i*2*pi*(k*df)*(dt*n)))*dt;
%       FV25(k)=sum(V25.*exp(-i*2*pi*(k*df)*(dt*n)))*dt;
      k=k+1;
  end
  f=df*(1:Nf);
C11=abs(-j*(FI15-FI11)/(80*dz)./FV13./(2*pi*f));
C12=abs(-j*(FI15-FI11)/(80*dz)./FV23./(2*pi*f));
C21=abs(-j*(FI25-FI21)/(80*dz)./FV13./(2*pi*f));
C22=abs(-j*(FI25-FI21)/(80*dz)./FV23./(2*pi*f));
f=f*1e-9;  % 将频率由Hz换成GHz为单位
figure(4);
plot(f,C11,'-k');
hold on;
plot(f,C12,'-b');
plot(f,C21,'--r');
plot(f,C22,'.-g');
Timeend=datestr(now)
Timecosted=Timeend-Timestart;
Timecosted

⌨️ 快捷键说明

复制代码 Ctrl + C
搜索代码 Ctrl + F
全屏模式 F11
切换主题 Ctrl + Shift + D
显示快捷键 ?
增大字号 Ctrl + =
减小字号 Ctrl + -