Excitation of a resonator in polar coordinates.

How to use openEMS. Discussion on examples, tutorials etc

Moderator: thorsten

Post Reply
Hale_812
Posts: 171
Joined: Fri 13 May 2016, 02:54

Excitation of a resonator in polar coordinates.

Post by Hale_812 » Fri 31 Jan 2020, 02:54

I am trying to simulate a resonator excited by a magnetic loop. Since I am interested more in fundamental modes, and do not know the impedance yet, I place a magnetic probe inside (see the wedges in the picture) instead of wrappint the loop around.
Taking one probe was not possible, because it creates a field of opposite direction on opposite side of the disk. So I made another loop with "[0, -1,0]" excitation.
Still, it does not work, as you see in the figure.

Probably I should take a rectangular (XYZ) grid, But it will have excessive density owing to curvature of the resonator.

Is there a way to solve the problem in cylindrical coordiantes?
キャプチャ.PNG
キャプチャ.PNG (201.66 KiB) Viewed 1356 times

Code: Select all

fclose ("all")
%% prepare simulation folder & run
Sim_Path = ['tmp_' mfilename];
Sim_CSX  = [mfilename '.xml'];
 
[status, message, messageid] = rmdir( Sim_Path, 's' ); % clear previous directory
[status, message, messageid] = mkdir( Sim_Path ); % create empty simulation folder

res=40;

physical_constants;
unit = 1e-3; % all length in mm

%%Frequency band
f0 = 60e6; % center frequency
fc = 40e6; % 20 dB corner frequency
max_res = c0 / (f0+fc) / unit / 4/ sqrt(15)/ res;

%%    Configure the size, substrate properties, and feed point
 
%substrate setup
conf.eps   = 15;
conf.kappa  = 1 * 2*pi*2.45e9 * EPS0*conf.eps;;%1e-3 * 2*pi*2.45e9 * EPS0*conf.eps; % strong loss for analysis
conf.r  = 33;
conf.h = 2;
conf.feed.r=conf.r*.95;
conf.feed.a=max_res/conf.r;
conf.feed.y=0;
conf.port.R = 50;     %port impedance
conf.ridge.l=0;%conf.r/3;

max_ang = conf.feed.a;%max_res/conf.r/4; % max_res in radians

% setup FDTD parameter & excitation function
FDTD = InitFDTD('CoordSystem', 1,'MultiGrid',conf.r/2, 'MaxTime', (60*60*10)); % init a cylindrical FDTD
FDTD = SetGaussExcite( FDTD, f0, fc );
 
%% setup CSXCAD geometry & mesh
% init a cylindrical mesh
CSX = InitCSX('CoordSystem',1);

%% create geometry
vert=[];mesh=[];
vert.r=[];vert.a=[];vert.z=[];
mesh.r=[];mesh.a=[];mesh.z=[];

%% create dielectric
CSX = AddMaterial( CSX, 'diel' );
CSX = SetMaterialProperty( CSX, 'diel', 'Epsilon', conf.eps, 'Kappa', conf.kappa );
start = [0,0,0];
stop  = [conf.r, 2*pi, conf.h ];
CSX = AddBox( CSX, 'diel', 0, start, stop);
vert = addnodes(vert, start, stop);

%% add plates
CSX = AddMetal( CSX, 'PEC' );
% bottom PEC plane
start = [0,0,-min(max_res/4,conf.h/4)];
stop  = [conf.r, 2*pi, 0 ];
CSX = AddBox( CSX, 'PEC', 999, start, stop);
vert = addnodes(vert, start, stop);
start = [0,0,conf.h];
stop  = [conf.r, 2*pi, conf.h+min(max_res/4,conf.h/4) ];
CSX = AddBox( CSX, 'PEC', 999, start, stop);
vert = addnodes(vert, start, stop);

%Add an excitation port

CSX = AddExcitation( CSX, 'DipoleP', 2, [0 1 0] );
start = [max_res/4, -conf.feed.a/6, 0];
stop  = [conf.feed.r, +conf.feed.a/6, conf.h];
CSX = AddBox( CSX, 'DipoleP', 1, start, stop );
vert = addnodes(vert, start, stop);

CSX = AddExcitation( CSX, 'DipoleN', 2, [0 -1 0] );
start = [max_res/4, deg2rad(180)-conf.feed.a/6, 0];
stop  = [conf.feed.r-conf.ridge.l, deg2rad(180)+conf.feed.a/6, conf.h];
CSX = AddBox( CSX, 'DipoleN', 1, start, stop );
vert = addnodes(vert, start, stop);

%% define voltage calc boxes %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
CSX = AddProbe(CSX,'ut1',0);
start = [ conf.r/2, deg2rad(120), 0 ];
stop  = [ conf.r/2+max_res/4, deg2rad(120), conf.h ];
CSX = AddBox(CSX,'ut1', 0 ,start,stop);
vert = addnodes(vert, start, stop);
CSX = AddProbe(CSX,'ut2',0);
start = [ (conf.r-conf.ridge.l)/2, deg2rad(180), 0 ];
stop  = [ (conf.r-conf.ridge.l)/2+max_res/4, deg2rad(180), conf.h ];
CSX = AddBox(CSX,'ut2', 0 ,start,stop);
vert = addnodes(vert, start, stop);




mesh=vert;
mesh.z=[mesh.z,0:conf.h/ceil(conf.h/max_res+3.5):conf.h];
mesh.a=mesh.a-floor(mesh.a./(2*pi)).*2*pi;
mesh.a(end+1)=min(mesh.a)+2*pi;
mesh.r=unique(mesh.r);
mesh.a=unique(mesh.a);
mesh.z=unique(mesh.z);
%mesh = SmoothMesh(mesh, [max_res, max_ang, max_res], 1.4);
mesh.r = AutoSmoothMeshLines(mesh.r, max_res,1.4);
mesh.a = AutoSmoothMeshLines(mesh.a, max_ang);
if not(rem(size(mesh.a,2),2))
  mesh.a=sort(mesh.a);
  temp=diff(mesh.a);
  temp=find( (temp-max(temp))==0 ,1);
  mesh.a=[mesh.a(1:temp),(mesh.a(temp)+mesh.a(temp+1))/2,mesh.a((temp+1):end)];
  temp=[];  clear temp;
end

%% define dump boxes... %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
CSX = AddDump(CSX,'Et_','DumpType',0, 'DumpMode',0);

ipr=findi(mesh.a,deg2rad(120));
start = [mesh.r(1) , mesh.a(ipr) , mesh.z(1)];
stop = [mesh.r(end) , mesh.a(ipr) , mesh.z(end)];
CSX = AddBox(CSX,'Et_',0 , start,stop);

CSX = AddDump(CSX,'Ht_','DumpType',1,'DumpMode',0);
CSX = AddBox(CSX,'Ht_',0,start,stop);

%current calc, for each position there are two currents, which will get
%averaged to match the voltage position in between (!Yee grid!)
CSX = AddProbe(CSX,'it1a',1);

ih2=findi(mesh.z,conf.h/2);
start = [mesh.r(1),   mesh.a(ipr-1),  0 ];
stop  = [mesh.r(end), mesh.a(ipr-1), mesh.z(ih2) ];
CSX = AddBox(CSX,'it1a', 0 ,start,stop);
CSX = AddProbe(CSX,'it1b',1);
start = [mesh.r(1),   mesh.a(ipr),   0 ];
stop  = [mesh.r(end), mesh.a(ipr),   mesh.z(ih2) ];
CSX = AddBox(CSX,'it1b', 0 ,start,stop);


% start = [ (conf.r-conf.ridge.l)/2, deg2rad(180), 0 ];
% stop  = [ (conf.r-conf.ridge.l)/2+max_res, deg2rad(180), conf.h ];
ipr=findi(mesh.a,deg2rad(180));
CSX = AddProbe(CSX,'it2a',1);
start = [ 0,  mesh.a(ipr-1),   0 ];
stop  = [ (conf.r-conf.ridge.l), mesh.a(ipr-1), mesh.z(ih2) ];
CSX = AddBox(CSX,'it2a', 0 ,start,stop);
CSX = AddProbe(CSX,'it2b',1);
start = [ 0,   mesh.a(ipr),   0 ];
stop  = [ (conf.r-conf.ridge.l), mesh.a(ipr), mesh.z(ih2) ];
CSX = AddBox(CSX,'it2b', 0 ,start,stop);


%Add a dump box for the surface current
CSX = AddDump(CSX, 'Hv','DumpType',1,'FileType',0);
start = [mesh.r(1), mesh.a(1), mesh.z(1)];
stop  = [mesh.r(end), mesh.r(end),  mesh.z(end)];
CSX = AddBox( CSX, 'Hv', 0, start, stop );


mesh = AddPML(mesh, [0 8 0 0 0 0],1); % add 8 lines in both z-directions
BC = {'MUR' 'PML' 'MUR' 'MUR' 'PEC' 'PEC'}; % boundary conditions
FDTD = SetBoundaryCond( FDTD, BC );

disp(['Num of cells: ' num2str(numel(mesh.r)*numel(mesh.a)*numel(mesh.z))]);
CSX = DefineRectGrid( CSX, unit, mesh );



%% write openEMS compatible xml-file
WriteOpenEMS( [Sim_Path '/' Sim_CSX], FDTD, CSX );

% show the structure
CSXGeomPlot( [Sim_Path '/' Sim_CSX] );

RunOpenEMS(Sim_Path, Sim_CSX);

%%
freq = linspace( f0-fc, f0+fc, 501 );
U = ReadUI({'ut1','ut2'},[Sim_Path '/'],freq);
I = ReadUI({'it1a','it1b','it2a','it2b'},[Sim_Path '/'],freq);
Exc = ReadUI('et',Sim_Path,freq);

%% plot voltages
figure(1)
plot(U.TD{1}.t, U.TD{1}.val,'Linewidth',2);
hold on;
grid on;
plot(U.TD{2}.t, U.TD{2}.val,'r--','Linewidth',2);
xlabel('time (s)')
ylabel('voltage (V)')
legend('u_1(t)','u_2(t)')

%% calculate incoming and reflected voltages & currents
ZL_a = ones(size(freq))*Z0/2/pi/sqrt(conf.eps); %analytic line-impedance of a coax

uf1 = U.FD{1}.val./Exc.FD{1}.val;
uf2 = U.FD{2}.val./Exc.FD{1}.val;
if1 = 0.5*(I.FD{1}.val+I.FD{2}.val)./Exc.FD{1}.val;
if2 = 0.5*(I.FD{3}.val+I.FD{4}.val)./Exc.FD{1}.val;

uf1_inc = 0.5 * ( uf1 + if1 .* ZL_a );
if1_inc = 0.5 * ( if1 + uf1 ./ ZL_a );
uf2_inc = 0.5 * ( uf2 + if2 .* ZL_a );
if2_inc = 0.5 * ( if2 + uf2 ./ ZL_a );

uf1_ref = uf1 - uf1_inc;
if1_ref = if1 - if1_inc;
uf2_ref = uf2 - uf2_inc;
if2_ref = if2 - if2_inc;

% plot s-parameter
figure(2)
s11 = uf1_ref./uf1_inc;
s21 = uf2_inc./uf1_inc;
plot(freq,20*log10(abs(s11)),'Linewidth',2);
xlim([freq(1) freq(end)]);
xlabel('frequency (Hz)')
ylabel('s-para (dB)');
% ylim([-40 5]);
grid on;
hold on;
plot(freq,20*log10(abs(s21)),'r','Linewidth',2);
legend('s11','s21','Location','SouthEast');

% plot line-impedance comparison
figure(3)
ZL = uf1./if1;
plot(freq,real(ZL),'Linewidth',2);
xlim([freq(1) freq(end)]);
xlabel('frequency (Hz)')
ylabel('line-impedance (\Omega)');
grid on;
hold on;
plot(freq,imag(ZL),'r--','Linewidth',2);
plot(freq,ZL_a,'g-.','Linewidth',2);
legend('\Re\{ZL\}','\Im\{ZL\}','ZL-analytic','Location','Best');


%% transform functions
paren = @(x, varargin) x(varargin{:});
curly = @(x, varargin) x{varargin{:}};

function ix=findi(meshx,x)
    diff=abs(meshx-x);
    ix=find((diff-min(diff))==0);
end

function a=aw(w,r)
%Angular Width
% transforms tangential width w to angular dimension a at radius r
    a=2.*atan( (w./2)./r);
end

function nout=addnodes(nin,varargin)
% Add Nodes
% Add coordinates of vertices vx to a list of nodes n
nout=nin;
    for i=1:size(varargin,2)
            temp=varargin{i};
        if isstruct(temp)
            nout.r=[nout.r,temp.r];
            nout.a=[nout.a,temp.a];
            nout.z=[nout.z,temp.z];
        else
            nout.r=[nout.r,temp(1)];
            nout.a=[nout.a,temp(2)];
            nout.z=[nout.z,temp(3)];            
        end
    end
end

thorsten
Posts: 1393
Joined: Mon 27 Jun 2011, 12:26

Re: Excitation of a resonator in polar coordinates.

Post by thorsten » Fri 31 Jan 2020, 08:48

I'm not sure what you mean or what is not working?? The wedges should not look like wedges? Or what do you mean?

Hale_812
Posts: 171
Joined: Fri 13 May 2016, 02:54

Re: Excitation of a resonator in polar coordinates.

Post by Hale_812 » Mon 03 Feb 2020, 02:23

the field. The sourceit is broken into two dipoles. it leaks in the center, giving preference to wrong modes. I was wondering, if there is a right way of exciting the field uniformly along the diameter line? But looks like no.
at the moment I am changing the mesh generator script to rectangular coords, but the mesh becomes menacingly dense.

Post Reply