Finite Difference Time Domain with multi materials

5 visualizzazioni (ultimi 30 giorni)
zina shadidi
zina shadidi il 19 Set 2021
Modificato: darova il 19 Set 2021
This is the finite element time domain for one barrier
I want to add a barrier , the loop was doubled, But I can not see the effect of the second barrier on figure.
I have no error but the results was not resonable
Any advice please
clear;
clc;
%%%%% for the first material
%epsr is the dielectric constant
apsilonr=41.04;
%epsz0 is the permitivity of space
apsilonz0=8.854223e-12;
%sigma is the conductivity
sigma=0.86674;
%KE is the number of cells to be used
KE=400;
dx=zeros(KE,1);
Ex=zeros(KE,1);
Hy=zeros(KE,1);
ix=zeros(KE,1);
gax=ones(KE,1);% amplitude of wave in +ve y axis direction
gbx=zeros(KE,1); %amplitude of wave in -ve y axis direction
%kc is the number of cell in which asinusoidal wave originates
p = zeros(KE,1);
sar = zeros(KE,1);
kc=20;
%f0 the corresponding frequancy of the wave
f0=1800e6;
%c0 the speed of light in avacuum
c0=3e8;
lambdaM1=c0/(sqrt(apsilonr)*f0);
ddx=lambdaM1/10;
dt=ddx./(2*c0);
Ex_low_m2=0;
Ex_low_m1=0;
T =0;
c=0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for the second material
%epsr is the dielectric constant
apsilonr1=52.725;
%epsz0 is the permitivity of space
apsilonz0=8.854223e-12;
%sigma is the conductivity
sigma1=0.942;
%KE is the number of cells to be used
KE=400;
dx1=zeros(KE,1);
Ex1=zeros(KE,1);
Hy1=zeros(KE,1);
ix1=zeros(KE,1);
gax1=ones(KE,1);% amplitude of wave in +ve y axis direction
gbx1=zeros(KE,1); %amplitude of wave in -ve y axis direction
%kc is the number of cell in which asinusoidal wave originates
p1 = zeros(KE,1);
sar1 = zeros(KE,1);
kc=20;
%f0 the corresponding frequancy of the wave
f0=1800e6;
%c0 the speed of light in avacuum
c0=3e8;
lambdaM1=c0/(sqrt(apsilonr1)*f0);
ddx1=lambdaM1/10;
dt=ddx./(2*c0);
Ex1_low_m2=0;
Ex1_low_m1=0;
T =0;
c=0;
%mu= % permeability
apsilonrr1= sigma1/(2*pi*f0) ; %imaginery factor of apsilonr
% attenuation factor alpha(Np/m), phase factor beta(rad/m)
%alpha= (2*pi* f0)*sqr((mu*apsilonr/2)*(sqr(1+(apsilonr/apsilonrr)^2)-1))
%beta= (2*pi* f0)*sqr((mu*apsilonr/2)*(sqr(1+(apsilonr/apsilonrr)^2)+1))
%FDTD loop
%%iteration = 3000; %%% here is the division to different layers
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%mu= % permeability
apsilonrr1= sigma1/(2*pi*f0) ; %imaginery factor of apsilonr
% attenuation factor alpha(Np/m), phase factor beta(rad/m)
%alpha= (2*pi* f0)*sqr((mu*apsilonr/2)*(sqr(1+(apsilonr/apsilonrr)^2)-1))
%beta= (2*pi* f0)*sqr((mu*apsilonr/2)*(sqr(1+(apsilonr/apsilonrr)^2)+1))
%FDTD loop
%%iteration = 3000; %%% here is the division to different layers
for k=120:KE
gax(k,1)=1/(apsilonr1+(sigma1*dt/apsilonz0));
gbx(k,1)=sigma1*dt/apsilonz0;
end
for k=100:KE
gax(k,1)=1/(apsilonr+(sigma*dt/apsilonz0));
gbx(k,1)=sigma*dt/apsilonz0;
end
for m=100:500:1000
c=c + 1;
for n=1:m
T = T + 1;
for k=2:KE
dx(k)=dx(k)+.5*(Hy(k-1)-Hy(k));
end
source=sin(2*pi*f0*n*dt);
dx(kc)=dx(kc)+source;
%calculate the E field of the FDTD
for k=2:KE
Ex(k)=gax(k,1)*(dx(k)-ix(k));
ix(k)=ix(k)+gbx(k,1)*Ex(k);
end
%absorbing boundary condition added
Ex1(1)=Ex1_low_m2;
Ex1_low_m2=Ex1_low_m1;
Ex1_low_m1=Ex1(2);
%calulate the H field of the FDTD
for k=1:(KE-1)
Hy(k)=Hy(k)+.5*(Ex(k)-Ex(k+1));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
for k=2:KE
dx1(k)=dx1(k)+.5*(Hy1(k-1)-Hy(k));
end
source1=dx(kc);
dx1(kc)=dx(kc)+source1;
%calculate the E field of the FDTD
for k=2:KE
Ex1(k)=gax1(k,1)*(dx1(k)-ix(k));
ix1(k)=ix(k)+gbx(k,1)*Ex1(k);
end
%absorbing boundary condition added
Ex(1)=Ex_low_m2;
Ex_low_m2=Ex_low_m1;
Ex_low_m1=Ex(2);
%calulate the H field of the FDTD
for k=1:(KE-1)
Hy1(k)=Hy1(k)+.5*(Ex1(k)-Ex1(k+1));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(c); plot(Ex,'b'); hold on;
xa=[120 120];
ya=[-2 2];
xb=[100 100];
yb=[-2 2];
plot(xa,ya,xb,yb,'r','LineWidth',2,'LineStyle','--')
title(['T = ',num2str(n)]);
ylabel('Propagation of E')
xlabel('FDTD ')
axis([1 KE -2 2])
text(250,-1.15,'Skin')
text(250,-1.3,'in head')
text(30,-1.3,'Epsilon=1','color','k')
hold off
end
end
for k = 1:19
p (k,1) = (sigma * (Ex(k,1)))^2 / 2;
end
for k = 1:20
sar(k,1) = p(k,1) / k;
end
figure;
plot(Ex);
ylabel(' Ex')
xlabel('x ')
figure;
plot(Hy);
ylabel('Hx')
xlabel('y ')
figure;
plot(p);
ylabel('P')
xlabel('x ')
figure;
plot(sar);
ylabel('SAR')
xlabel('x ')

Risposte (0)

Categorie

Scopri di più su Programming in Help Center e File Exchange

Prodotti

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by