new lumped elements using new CSPropLumpedElement CSXCAD-class
Thorsten Liebig
parent
62acf5f1b3
commit
7b34a8706a
|
|
@ -824,6 +824,8 @@ int Operator::CalcECOperator( DebugFlags debugFlags )
|
|||
|
||||
CalcPEC();
|
||||
|
||||
Calc_LumpedElements();
|
||||
|
||||
bool PMC[6];
|
||||
for (int n=0; n<6; ++n)
|
||||
PMC[n] = m_BC[n]==1;
|
||||
|
|
@ -1161,6 +1163,165 @@ bool Operator::Calc_EffMatPos(int ny, const unsigned int* pos, double* EffMat) c
|
|||
return true;
|
||||
}
|
||||
|
||||
bool Operator::Calc_LumpedElements()
|
||||
{
|
||||
vector<CSProperties*> props = CSX->GetPropertyByType(CSProperties::LUMPED_ELEMENT);
|
||||
for (size_t i=0;i<props.size();++i)
|
||||
{
|
||||
CSPropLumpedElement* PLE = dynamic_cast<CSPropLumpedElement*>(props.at(i));
|
||||
if (PLE==NULL)
|
||||
return false; //sanity check: this should never happen!
|
||||
vector<CSPrimitives*> prims = PLE->GetAllPrimitives();
|
||||
for (size_t bn=0;bn<prims.size();++bn)
|
||||
{
|
||||
CSPrimBox* box = dynamic_cast<CSPrimBox*>(prims.at(bn));
|
||||
if (box)
|
||||
{ //calculate lumped element parameter
|
||||
|
||||
double C = PLE->GetCapacity();
|
||||
if (C<=0)
|
||||
C = NAN;
|
||||
double R = PLE->GetResistance();
|
||||
if (R<0)
|
||||
R = NAN;
|
||||
|
||||
if ((isnan(R)) && (isnan(C)))
|
||||
{
|
||||
cerr << "Operator::Calc_LumpedElements(): Warning: Lumped Element R or C not specified! skipping. "
|
||||
<< " ID: " << prims.at(bn)->GetID() << " @ Property: " << PLE->GetName() << endl;
|
||||
break;
|
||||
}
|
||||
|
||||
int ny = PLE->GetDirection();
|
||||
if ((ny<0) || (ny>2))
|
||||
{
|
||||
cerr << "Operator::Calc_LumpedElements(): Warning: Lumped Element direction is invalid! skipping. "
|
||||
<< " ID: " << prims.at(bn)->GetID() << " @ Property: " << PLE->GetName() << endl;
|
||||
break;
|
||||
}
|
||||
int nyP = (ny+1)%3;
|
||||
int nyPP = (ny+2)%3;
|
||||
|
||||
unsigned int uiStart[3];
|
||||
unsigned int uiStop[3];
|
||||
// snap to the native coordinate system
|
||||
SnapToMesh(box->GetStartCoord()->GetNativeCoords(),uiStart);
|
||||
SnapToMesh(box->GetStopCoord()->GetNativeCoords(),uiStop);
|
||||
|
||||
for (int n=0;n<3;++n)
|
||||
{
|
||||
if (uiStop[n]<uiStart[n])
|
||||
{
|
||||
unsigned int help = uiStart[n];
|
||||
uiStart[n]=uiStop[n];
|
||||
uiStop[n]=help;
|
||||
}
|
||||
}
|
||||
|
||||
if (uiStart[ny]==uiStop[ny])
|
||||
{
|
||||
cerr << "Operator::Calc_LumpedElements(): Warning: Lumped Element with zero (snapped) length is invalid! skipping. "
|
||||
<< " ID: " << prims.at(bn)->GetID() << " @ Property: " << PLE->GetName() << endl;
|
||||
break;
|
||||
}
|
||||
|
||||
//calculate geometric property for this lumped element
|
||||
unsigned int pos[3];
|
||||
double unitGC=0;
|
||||
int ipos=0;
|
||||
for (pos[ny]=uiStart[ny];pos[ny]<uiStop[ny];++pos[ny])
|
||||
{
|
||||
double unitGC_Plane=0;
|
||||
for (pos[nyP]=uiStart[nyP];pos[nyP]<=uiStop[nyP];++pos[nyP])
|
||||
{
|
||||
for (pos[nyPP]=uiStart[nyPP];pos[nyPP]<=uiStop[nyPP];++pos[nyPP])
|
||||
{
|
||||
// capacity/conductivity in parallel: add values
|
||||
unitGC_Plane += GetEdgeArea(ny,pos)/GetEdgeLength(ny,pos);
|
||||
}
|
||||
}
|
||||
|
||||
//capacity/conductivity in series: add reciprocal values
|
||||
unitGC += 1/unitGC_Plane;
|
||||
}
|
||||
unitGC = 1/unitGC;
|
||||
|
||||
bool caps = PLE->GetCaps();
|
||||
double kappa = 0;
|
||||
double epsilon = 0;
|
||||
if (R>0)
|
||||
kappa = 1 / R / unitGC;
|
||||
if (C>0)
|
||||
epsilon = C / unitGC;
|
||||
|
||||
if (epsilon< __EPS0__)
|
||||
{
|
||||
cerr << "Operator::Calc_LumpedElements(): Warning: Lumped Element capacity is too small for its size! skipping. "
|
||||
<< " ID: " << prims.at(bn)->GetID() << " @ Property: " << PLE->GetName() << endl;
|
||||
C = 0;
|
||||
if (isnan(R))
|
||||
break;
|
||||
}
|
||||
|
||||
for (pos[ny]=uiStart[ny];pos[ny]<uiStop[ny];++pos[ny])
|
||||
{
|
||||
for (pos[nyP]=uiStart[nyP];pos[nyP]<=uiStop[nyP];++pos[nyP])
|
||||
{
|
||||
for (pos[nyPP]=uiStart[nyPP];pos[nyPP]<=uiStop[nyPP];++pos[nyPP])
|
||||
{
|
||||
ipos = MainOp->SetPos(pos[0],pos[1],pos[2]);
|
||||
if (C>0)
|
||||
EC_C[ny][ipos] = epsilon * GetEdgeArea(ny,pos)/GetEdgeLength(ny,pos);
|
||||
if (R>0)
|
||||
EC_G[ny][ipos] = kappa * GetEdgeArea(ny,pos)/GetEdgeLength(ny,pos);
|
||||
if (R==0) //make lumped element a PEC if resistance is zero
|
||||
{
|
||||
SetVV(ny,pos[0],pos[1],pos[2], 0 );
|
||||
SetVI(ny,pos[0],pos[1],pos[2], 0 );
|
||||
}
|
||||
else //recalculate operator inside the lumped element
|
||||
Calc_ECOperatorPos(ny,pos);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// setup metal caps
|
||||
if (caps)
|
||||
{
|
||||
for (pos[nyP]=uiStart[nyP];pos[nyP]<=uiStop[nyP];++pos[nyP])
|
||||
{
|
||||
for (pos[nyPP]=uiStart[nyPP];pos[nyPP]<=uiStop[nyPP];++pos[nyPP])
|
||||
{
|
||||
pos[ny]=uiStart[ny];
|
||||
SetVV(nyP,pos[0],pos[1],pos[2], 0 );
|
||||
SetVI(nyP,pos[0],pos[1],pos[2], 0 );
|
||||
++m_Nr_PEC[nyP];
|
||||
|
||||
SetVV(nyPP,pos[0],pos[1],pos[2], 0 );
|
||||
SetVI(nyPP,pos[0],pos[1],pos[2], 0 );
|
||||
++m_Nr_PEC[nyPP];
|
||||
|
||||
pos[ny]=uiStop[ny];
|
||||
SetVV(nyP,pos[0],pos[1],pos[2], 0 );
|
||||
SetVI(nyP,pos[0],pos[1],pos[2], 0 );
|
||||
++m_Nr_PEC[nyP];
|
||||
|
||||
SetVV(nyPP,pos[0],pos[1],pos[2], 0 );
|
||||
SetVI(nyPP,pos[0],pos[1],pos[2], 0 );
|
||||
++m_Nr_PEC[nyPP];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
else
|
||||
cerr << "Operator::Calc_LumpedElements(): Warning: Primitves other than boxes are not supported for lumped elements! skipping "
|
||||
<< prims.at(bn)->GetTypeName() << " ID: " << prims.at(bn)->GetID() << " @ Property: " << PLE->GetName() << endl;
|
||||
}
|
||||
}
|
||||
return true;
|
||||
}
|
||||
|
||||
void Operator::DumpExciationSignals()
|
||||
{
|
||||
Exc->DumpVoltageExcite("et");
|
||||
|
|
|
|||
|
|
@ -188,6 +188,9 @@ protected:
|
|||
//! Calc operator at certain \a pos
|
||||
virtual void Calc_ECOperatorPos(int n, unsigned int* pos);
|
||||
|
||||
//! Calculate and setup lumped elements
|
||||
virtual bool Calc_LumpedElements();
|
||||
|
||||
//store material properties for post-processing
|
||||
float**** m_epsR;
|
||||
float**** m_kappa;
|
||||
|
|
|
|||
|
|
@ -0,0 +1,158 @@
|
|||
%
|
||||
% EXAMPLE / other / lumped elements
|
||||
%
|
||||
% This example demonstrates how to:
|
||||
% - use lumped elements
|
||||
%
|
||||
%
|
||||
% Tested with
|
||||
% - Matlab 2009b
|
||||
% - openEMS v0.0.21-3
|
||||
%
|
||||
% (C) 2010 Thorsten Liebig <thorsten.liebig@uni-due.de>
|
||||
|
||||
close all
|
||||
clear
|
||||
clc
|
||||
|
||||
%% setup the simulation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
f_max = 100e6;
|
||||
f_excite = 300e6;
|
||||
SimBox = 100;
|
||||
mesh_size = 2;
|
||||
|
||||
Lumped.R = 1000;
|
||||
Lumped.C = 10e-12;
|
||||
|
||||
% the parasitice inductance of the feeding has to be deduced with a R=0
|
||||
% simulation
|
||||
parasitic_L = 63e-9;
|
||||
|
||||
%% define openEMS options %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
openEMS_opts = '';
|
||||
% openEMS_opts = [openEMS_opts ' --debug-material'];
|
||||
% openEMS_opts = [openEMS_opts ' --debug-boxes'];
|
||||
% openEMS_opts = [openEMS_opts ' --debug-operator'];
|
||||
|
||||
Sim_Path = 'tmp';
|
||||
Sim_CSX = 'lumped.xml';
|
||||
|
||||
[status, message, messageid] = rmdir(Sim_Path,'s');
|
||||
[status,message,messageid] = mkdir(Sim_Path);
|
||||
|
||||
%% setup FDTD parameter & excitation function %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
FDTD = InitFDTD(30000,1e-6);
|
||||
FDTD = SetGaussExcite(FDTD,f_excite/2,f_excite/2);
|
||||
BC = [1 1 1 1 1 1];
|
||||
FDTD = SetBoundaryCond(FDTD,BC);
|
||||
|
||||
%% setup CSXCAD geometry & mesh %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
CSX = InitCSX();
|
||||
mesh.x = SmoothMeshLines([-SimBox/2,+SimBox/2],mesh_size);
|
||||
mesh.y = SmoothMeshLines([-SimBox/2,+SimBox/2],mesh_size);
|
||||
mesh.z = SmoothMeshLines([-SimBox/2,+SimBox/2],mesh_size);
|
||||
CSX = DefineRectGrid(CSX, 1e-3,mesh);
|
||||
|
||||
|
||||
%% create structure
|
||||
% insert curve port
|
||||
start = [ 10 -10 0];
|
||||
stop = [ 10 10 0];
|
||||
CSX = AddCurvePort(CSX,0,1,100,start,stop,'excite');
|
||||
|
||||
% insert lumped element
|
||||
CSX = AddLumpedElement( CSX, 'Capacitor', 1, 'C', Lumped.C, 'R', Lumped.R);
|
||||
start = [ -14 -4 -4];
|
||||
stop = [ -6 4 4];
|
||||
CSX = AddBox( CSX, 'Capacitor', 0, start, stop );
|
||||
|
||||
% insert feeding wire
|
||||
CSX = AddMetal(CSX,'metal');
|
||||
%first point
|
||||
points(1,1) = -10;
|
||||
points(2,1) = 4;
|
||||
points(3,1) = 0;
|
||||
%second point
|
||||
points(1,2) = -10;
|
||||
points(2,2) = 15;
|
||||
points(3,2) = 0;
|
||||
%3 point
|
||||
points(1,end+1) = 10;
|
||||
points(2,end) = 15;
|
||||
points(3,end) = 0;
|
||||
%4 point
|
||||
points(1,end+1) = 10;
|
||||
points(2,end) = 10;
|
||||
points(3,end) = 0;
|
||||
CSX = AddCurve(CSX,'metal', 10, points);
|
||||
|
||||
points(2,:) = -1*points(2,:);
|
||||
CSX = AddCurve(CSX,'metal', 10, points);
|
||||
|
||||
%% Write openEMS compatoble xml-file %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
WriteOpenEMS([Sim_Path '/' Sim_CSX],FDTD,CSX);
|
||||
|
||||
% CSXGeomPlot([Sim_Path '/' Sim_CSX]);
|
||||
|
||||
%% run openEMS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
RunOpenEMS(Sim_Path, Sim_CSX, openEMS_opts);
|
||||
|
||||
%% postproc & do the plots %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
f = linspace(1e6,f_max,1001);
|
||||
w = 2*pi*f;
|
||||
% read currents and voltages
|
||||
U = ReadUI('port_ut1','tmp/',f);
|
||||
I = ReadUI('port_it1','tmp/',f);
|
||||
|
||||
% calculate analytic impedance
|
||||
if (Lumped.R>=0)
|
||||
Z_a = Lumped.R*(1-1i*w*Lumped.C*Lumped.R)./(1+(w*Lumped.C*Lumped.R).^2);
|
||||
else
|
||||
Z_a = -1i./(w*Lumped.C);
|
||||
end
|
||||
|
||||
% calculate numerical impedance
|
||||
Z = U.FD{1}.val./I.FD{1}.val;
|
||||
|
||||
% remove parasitic feeding effects
|
||||
Z = Z - 1i*w*parasitic_L;
|
||||
|
||||
L = imag(Z)./w;
|
||||
C = -1./(w.*imag(Z));
|
||||
C(find(C<0)) = nan;
|
||||
L(find(L<0)) = nan;
|
||||
R = real(Z);
|
||||
|
||||
subplot(2,1,1);
|
||||
plot(f*1e-6,C*1e12,'Linewidth',2);
|
||||
xlabel('frequency (MHz)');
|
||||
ylabel('capacitance (pF)');
|
||||
grid on;
|
||||
subplot(2,1,2);
|
||||
plot(f*1e-6,L*1e9,'Linewidth',2);
|
||||
xlabel('frequency (MHz)');
|
||||
ylabel('inductance (nH)');
|
||||
grid on;
|
||||
|
||||
figure();
|
||||
plot(f*1e-6,R,'Linewidth',2);
|
||||
hold on
|
||||
plot(f*1e-6,imag(Z),'r--','Linewidth',2);
|
||||
|
||||
plot(f*1e-6,real(Z_a),'g-.','Linewidth',1);
|
||||
plot(f*1e-6,imag(Z_a),'m--','Linewidth',1);
|
||||
|
||||
xlabel('frequency (MHz)');
|
||||
ylabel('resistance (Ohm)');
|
||||
grid on;
|
||||
legend( '\Re\{Z\}','\Im\{Z\}','\Re\{Z_{analytisch}\}','\Im\{Z_{analytisch}\}', 'location', 'northeast' )
|
||||
|
||||
figure();
|
||||
errorR = (R-real(Z_a))./R*100;
|
||||
errorX = (imag(Z)-imag(Z_a))./imag(Z)*100;
|
||||
plot(f*1e-6,errorR,'Linewidth',2);
|
||||
hold on
|
||||
grid on;
|
||||
plot(f*1e-6,errorX,'r--','Linewidth',2);
|
||||
xlabel('frequency (MHz)');
|
||||
ylabel('error (%)');
|
||||
Loading…
Reference in New Issue