caiyu/openEMS
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

drude type material: losses added

Browse Source
This commit is contained in:
Thorsten Liebig 2012-04-24 11:04:30 +02:00
parent cd3c8baefd
commit fb5b330d25
2 changed files with 63 additions and 23 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

37
FDTD/extensions/operator_ext_lorentzmaterial.cpp
View File

@ -52,9 +52,10 @@ bool Operator_Ext_LorentzMaterial::BuildExtension()
unsigned int numLines[3] = {m_Op->GetNumberOfLines(0),m_Op->GetNumberOfLines(1),m_Op->GetNumberOfLines(2)};
CSPropLorentzMaterial* mat = NULL;
double w_plasma;
double w_plasma,t_relax;
bool b_pos_on;
double L_D[3], C_D[3];
double R_D[3], G_D[3];
vector<double> v_int[3];
vector<double> v_ext[3];
vector<double> i_int[3];
@ -73,6 +74,7 @@ bool Operator_Ext_LorentzMaterial::BuildExtension()
for (int n=0; n<3; ++n)
{
L_D[n]=0;
R_D[n]=0;
coord[0] = m_Op->GetDiscLine(0,pos[0]);
coord[1] = m_Op->GetDiscLine(1,pos[1]);
coord[2] = m_Op->GetDiscLine(2,pos[2]);
@ -82,18 +84,24 @@ bool Operator_Ext_LorentzMaterial::BuildExtension()
if ((mat = prop->ToLorentzMaterial()))
{
w_plasma = mat->GetEpsPlasmaFreqWeighted(n,coord) * 2 * PI;
if (w_plasma)
if (w_plasma>0)
{
b_pos_on = true;
m_volt_ADE_On = true;
L_D[n] = 1/(w_plasma*w_plasma*m_Op->EC_C[n][index]);
}
t_relax = mat->GetEpsRelaxTimeWeighted(n,coord);
if ((t_relax>0) && m_volt_ADE_On)
{
R_D[n] = L_D[n]/t_relax;
}
}
}
for (int n=0; n<3; ++n)
{
C_D[n]=0;
G_D[n]=0;
coord[0] = m_Op->GetDiscLine(0,pos[0],true);
coord[1] = m_Op->GetDiscLine(1,pos[1],true);
coord[2] = m_Op->GetDiscLine(2,pos[2],true);
@ -103,12 +111,17 @@ bool Operator_Ext_LorentzMaterial::BuildExtension()
if ((mat = prop->ToLorentzMaterial()))
{
w_plasma = mat->GetMuePlasmaFreqWeighted(n,coord) * 2 * PI;
if (w_plasma)
if (w_plasma>0)
{
b_pos_on = true;
m_curr_ADE_On = true;
C_D[n] = 1/(w_plasma*w_plasma*m_Op->EC_L[n][index]);
}
t_relax = mat->GetMueRelaxTimeWeighted(n,coord);
if ((t_relax>0) && m_curr_ADE_On)
{
G_D[n] = C_D[n]/t_relax;
}
}
}
@ -117,16 +130,26 @@ bool Operator_Ext_LorentzMaterial::BuildExtension()
for (unsigned int n=0; n<3; ++n)
{
v_pos[n].push_back(pos[n]);
v_int[n].push_back(1);
i_int[n].push_back(1);
if (L_D[n]>0)
v_ext[n].push_back(dT/L_D[n] *m_Op->GetVI(n,pos[0],pos[1],pos[2]));
{
v_int[n].push_back((2*L_D[n]-dT*R_D[n])/(2*L_D[n]+dT*R_D[n]));
v_ext[n].push_back(dT/(L_D[n]+dT*R_D[n]/2)*m_Op->GetVI(n,pos[0],pos[1],pos[2]));
}
else
{
v_int[n].push_back(1);
v_ext[n].push_back(0);
}
if (C_D[n]>0)
i_ext[n].push_back(dT/C_D[n] *m_Op->GetIV(n,pos[0],pos[1],pos[2]));
{
i_int[n].push_back((2*C_D[n]-dT*G_D[n])/(2*C_D[n]+dT*G_D[n]));
i_ext[n].push_back(dT/(C_D[n]+dT*G_D[n]/2)*m_Op->GetIV(n,pos[0],pos[1],pos[2]));
}
else
{
i_int[n].push_back(1);
i_ext[n].push_back(0);
}
// cerr << v_int[n].back() << " " << v_ext[n].back() << " " << i_int[n].back() << " " << i_ext[n].back() << endl;
}
}

49
matlab/examples/other/Metamaterial_PlaneWave_Drude.m
View File

@ -1,7 +1,9 @@
%%%%%%%%%%%%%%%%%%%%%%%
% example demonstrating double drude meta-material
%
% author: Thorsten Liebig @ 2010
% tested with openEMS v0.0.28
%
% author: Thorsten Liebig @ 2010,2012
%%%%%%%%%%%%%%%%%%%%%%%
close all
@ -17,27 +19,27 @@ Settings.LogFile = 'openEMS.log';
pic_size = round([1400 1400/4]); %define the animation picture size
%simulation domain setup (in mm)
length = 4000;
width = 100;
mesh_res = 5; % mesh resolution
length = 500;
width = 10;
mesh_res = 0.5; % mesh resolution
height = 3*mesh_res; % hight is ony 3 lines with PEC (top/bottom) --> quasi 2D
%FDTD setup
f0 = 0.5e9; %center frequency
f0 = 5e9; %center frequency
f_BW = f0/sqrt(2); %bandwidth
MTM.eps_R = 1;
MTM.mue_R = 1;
MTM.f0 = f0; %plasma frequency of the drude material
MTM.length = 1000; %length of the metamaterial
MTM.relaxTime = 5e-9; %relaxation time (smaller number results in greater losses, set to 0 to disable)
MTM.length = 250; %length of the metamaterial
N_TS = 5e4; %number of timesteps
endCriteria = 1e-5; %stop simulation if signal is at -50dB
%constants
EPS0 = 8.85418781762e-12;
MUE0 = 1.256637062e-6;
physical_constants
%% define openEMS options %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
openEMS_opts = '';
% openEMS_opts = [openEMS_opts ' --disable-dumps'];
openEMS_opts = [openEMS_opts ' --engine=fastest'];
openEMS_opts = '-vvv';
Sim_Path = 'MTM_PW_Drude';
Sim_CSX = 'MTM_PW_Drude.xml';
@ -51,8 +53,7 @@ if (postproc_only==0)
%% setup FDTD parameter & excitation function %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
FDTD = InitFDTD(N_TS,endCriteria,'OverSampling',10);
FDTD = SetGaussExcite(FDTD,f0,f0/sqrt(2));
% FDTD = SetSinusExcite(FDTD,f0*sqrt(2));
FDTD = SetGaussExcite(FDTD,0,2*f0);
BC = [1 1 0 0 2 2];
FDTD = SetBoundaryCond(FDTD,BC);
@ -70,9 +71,9 @@ if (postproc_only==0)
CSX = AddBox(CSX,'excite',0 ,start,stop);
%% apply drude material %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
CSX = AddProperty(CSX,'LorentzMaterial','drude');
CSX = SetPropertyArgs(CSX,'LorentzMaterial','drude','PlasmaFrequency','Epsilon',MTM.f0);
CSX = SetPropertyArgs(CSX,'LorentzMaterial','drude','PlasmaFrequency','Mue' ,MTM.f0);
CSX = AddLorentzMaterial(CSX,'drude');
CSX = SetMaterialProperty(CSX,'drude','Epsilon',MTM.eps_R,'EpsilonPlasmaFrequency',MTM.f0,'EpsilonRelaxTime',MTM.relaxTime);
CSX = SetMaterialProperty(CSX,'drude','Mue',MTM.mue_R,'MuePlasmaFrequency',MTM.f0,'MueRelaxTime',MTM.relaxTime);
start=[mesh.x(1) mesh.y(1) -MTM.length/2];
stop =[mesh.x(end) mesh.y(end) MTM.length/2];
CSX = AddBox(CSX,'drude', 10 ,start,stop);
@ -91,6 +92,22 @@ if (postproc_only==0)
end
%% plot the drude type material dependency
f = linspace(0.1*f0,2*f0,501);
w = 2*pi*f;
epsr = MTM.eps_R * (1 - (2*pi*MTM.f0)^2./( w.^2 - 1j*w./MTM.relaxTime ));
muer = MTM.mue_R * (1 - (2*pi*MTM.f0)^2./( w.^2 - 1j*w./MTM.relaxTime ));
plot(f,real(epsr),'Linewidth',2);
hold on
grid on
plot(f,imag(epsr),'r--','Linewidth',2);
plot(f,real(muer),'c-.','Linewidth',2);
plot(f,imag(muer),'m-.','Linewidth',2);
ylim([-10 MTM.eps_R])
% l=legend('\Re \epsilon_r','\Im \epsilon_r','\Re \mue_r','\Im \mue_r');
l=legend('$\Re\{\varepsilon_r\}$','$\Im\{\varepsilon_r\}$','$\Re\{\mu_r\}$','$\Im\{\mu_r\}$');
set(l,'Interpreter','latex','Fontsize',12)
%% plot E-fields
freq = [f0/sqrt(2) f0 f0*sqrt(2)];
field = ReadHDF5FieldData([Sim_Path '/Et.h5']);