How mesh(x) work in Openems and Octave

[Dianebonk]

Hi Thomas, I am a beginner in using OpenEMS and I want to know how to properly use the mesh.x(number) command in OpenEMS. I am writing code to simulate a PCB trace antenna for my internship and I am having problems with some commands. This antenna must work at 2.45GHz for Bluetooth low-energy communication. I change the length of my antenna from 27mm to 30mm to see how the frequency varies. I wrote a message on the OpenEMS forum but I didn't get much response. Only one person replied and told me that my nf2ff dump box is not well defined and therefore my antenna is not centred. This is why I have negative values for gain and efficiency. Here are some images of my result:

Here is the code I am trying to correct. If you have an idea on how I could solve it, I take it:

% EXAMPLE / antennas / ) 2.4GHz

%

% This example demonstrates how to:

% - calculate the reflection coefficient of an antenna
% - Give the antenna impedance
% - calculate farfield of an antenna

%

close all

clear

clc

%% setup the simulation

physical_constants;

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

substrate.width = 10; % width of substrate

substrate.length = 10; % length of substrate

substrate.thickness = 0.8; % thickness of substrate

substrate.cells = 4; % use 4 cells for meshing substrate

ifa.h = 2.1; % height of short circuit stub

ifa.l = 8.875; % length of radiating element1

ifa.w1 = 0.9; % width of short circuit stub

ifa.w2 = 0.3; % width of radiating element

ifa.wf = 0.3; % width of feed element

ifa.fp = 4; % position of feed element relative to short

% circuit stub

ifa.e = 2.7; % distance to edge

ifa.b = 6.075; %Radiating element 3_5_7 length

ifa.a = 0.9; %Radiating element 2_4_6 length

ifa.c = 0.3; % distance between feeding elementand circuit stub

ifa.d = 0.3; %distance between feeding element and radiating element 2...

ifa.gnd = 0.3; %distance antenna to the ground plane

ifa.Ledge = 0.9; % distance to the left edge 1.1 Previously
ifa.cut = 0 ;  %Cutting element

ifa.sbt= 0.4; %reduce substrate length

length = 125 ;
% substrate setup

substrate.epsR = 4.3;

substrate.kappa = 1e-3 * 2*pi*2.45e9 * EPS0*substrate.epsR;

%setup feeding

feed.R = 50; %feed resistance

%open AppCSXCAD and show ifa

show = 1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% size of the simulation box

SimBox = [substrate.width*6 substrate.length*6 60];

%% setup FDTD parameter & excitation function

f0 = 2.4e9; % center frequency

fc = 1e9; % 20 dB corner frequency

c0 = 3e8; %Speed of light

%FDTD = InitFDTD(90000, 1e-3, 'TimeStep', 5e-13);
%FDTD = InitFDTD('OverSampling', 400) %highly discourage you to change the timestep.
%we either make the simulation slower or unstable.
%openEMS chooses the fastest possible and stable timestep. So change the oversampling

FDTD = InitFDTD('NrTS', 200000, 'EndCriteria', 1e-3);

%FDTD = InitFDTD('NrTS', 100000 );

FDTD = SetGaussExcite( FDTD, f0, fc );

BC = {'MUR' 'MUR' 'MUR' 'MUR' 'MUR' 'MUR'}; % simple absorbing boundary conditions because our antenna being destined to be in open space

FDTD = SetBoundaryCond( FDTD, BC );

%% setup CSXCAD geometry & mesh

CSX = InitCSX();

%initialize the mesh with the "air-box" dimensions

mesh.x = [-SimBox(1)/2 SimBox(1)/2];

mesh.y = [-SimBox(2)/2 SimBox(2)/2];

mesh.z = [-SimBox(3)/2 SimBox(3)/2];

%% create substrate

CSX = AddMaterial( CSX, 'substrate');

CSX = SetMaterialProperty( CSX, 'substrate', 'Epsilon',substrate.epsR, 'Kappa', substrate.kappa);

start = [-substrate.width/2 -substrate.length/2 0];

stop = [ substrate.width/2 substrate.length/2 substrate.thickness];

CSX = AddBox( CSX, 'substrate', 1, start, stop );

% add extra cells to discretize the substrate thickness

mesh.z = [linspace(0,substrate.thickness,substrate.cells+1) mesh.z];

%% create ground plane

CSX = AddMetal( CSX, 'groundplane' ); % create a perfect electric conductor (PEC)

start = [-substrate.width/2 -substrate.length/2 substrate.thickness];

stop = [ substrate.width/2 substrate.length/2-ifa.e-ifa.sbt substrate.thickness];

CSX = AddBox(CSX, 'groundplane', 10, start,stop);

%% create ifa
CSX = AddMetal( CSX, 'ifa' ); % create a perfect electric conductor (PEC)
tl = [ifa.Ledge,substrate.length/2-ifa.e,substrate.thickness];   % translate

start = [(-substrate.width/2)+ifa.w1+ifa.c (ifa.h-3*ifa.a+2*ifa.w2-0.3) 0] + tl;
stop = start + [ifa.wf ifa.h 0];
CSX = AddBox( CSX, 'ifa', 10,  start, stop);  % feed element

start = [-substrate.width/2  (ifa.h-3*ifa.a+2*ifa.w2)-ifa.sbt 0] + tl;
stop =  start + [ifa.w1 ifa.h+ifa.sbt 0];
CSX = AddBox( CSX, 'ifa', 10,  start, stop);  % short circuit stub

start = [(-substrate.width/2) ifa.h 0] + tl;
stop = start + [ifa.l -ifa.w2 0];
CSX = AddBox( CSX, 'ifa', 10, start, stop);   % radiating element1

start = [(-substrate.width/2+ifa.l) ifa.h 0] + tl;
stop = start + [-ifa.w2 -ifa.a 0];
CSX = AddBox( CSX, 'ifa', 10, start, stop);   % radiating element2

start = [(-substrate.width/2+ifa.l) ifa.h-ifa.a 0] + tl;
stop = start + [-ifa.b ifa.w2 0];
CSX = AddBox( CSX, 'ifa', 10, start, stop);   % radiating element3

start = [(-substrate.width/2+ifa.l-ifa.b+ifa.w2) ifa.h-ifa.a+ifa.w2 0] + tl;
stop = start + [-ifa.w2 -ifa.a 0];
CSX = AddBox( CSX, 'ifa', 10, start, stop);   % radiating element4

start = [(-substrate.width/2+ifa.l) (ifa.h-2*ifa.a+ifa.w2) 0] + tl;
stop = start + [-ifa.b ifa.w2 0];
CSX = AddBox( CSX, 'ifa', 10, start, stop);   % radiating element5

start = [(-substrate.width/2+ifa.l) ifa.h-2*ifa.a+2*ifa.w2 0] + tl;
stop = start + [-ifa.w2 -ifa.a 0];
CSX = AddBox( CSX, 'ifa', 10, start, stop);   % radiating element6

start = [(-substrate.width/2+ifa.l) (ifa.h-3*ifa.a+2*ifa.w2) 0] + tl;
stop = start + [-ifa.b+ifa.cut ifa.w2 0];
CSX = AddBox( CSX, 'ifa', 10, start, stop);   % radiating element7

%% apply the excitation & resist as a current source
start = [ (-substrate.width/2)+ifa.w1+ifa.c (ifa.h-3*ifa.a+2*ifa.w2-0.3) 0] + tl;
stop  = start + [ifa.wf -ifa.sbt+0.3 0];
[CSX port] = AddLumpedPort(CSX, 5 ,1 ,feed.R, start, stop, [0 1 0], true); %

% MESH left
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
mesh.x = [mesh.x linspace(-substrate.width*3,(-substrate.width/2)+ifa.Ledge-1,40.0)]; %tab 15 elements step of 0.4

mesh.y = [mesh.y linspace(-substrate.length*3,(substrate.length/2)-ifa.e-ifa.sbt-1,40.0)];%mesh.z = [mesh.z linspace(-10.0,5,40.0) + 0];
CSX = DefineRectGrid(CSX, 0.001, mesh);  %define the mesh with a drawing unit of 1mm (1e-3)

% MESH left2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
mesh.x = [mesh.x linspace((-substrate.width/2)+ifa.Ledge-1,(-substrate.width/2)+ifa.Ledge,4.0)]; %tab 15 elements step of 0.4

mesh.y = [mesh.y linspace((substrate.length/2)-ifa.e-ifa.sbt-1,(substrate.length/2)-ifa.e-ifa.sbt,5.0)];%mesh.z = [mesh.z linspace(-10.0,5,40.0) + 0];
CSX = DefineRectGrid(CSX, 0.001, mesh);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MESH right
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
mesh.x = [mesh.x linspace((-substrate.width/2)+ifa.l+ifa.Ledge+1,substrate.width*3,35.0)]; % a=4.97
mesh.y = [mesh.y linspace(substrate.width/2-ifa.e+ifa.h+1,substrate.length*3,35.0) ]; % a=5.4

CSX = DefineRectGrid(CSX, 0.001, mesh);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MESH right2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
mesh.x = [mesh.x linspace((-substrate.width/2)+ifa.l+ifa.Ledge,(-substrate.width/2)+ifa.l+ifa.Ledge+1,4.0)];
mesh.y = [mesh.y linspace(substrate.width/2-ifa.e+ifa.h,substrate.width/2-ifa.e+ifa.h+1,5.0) ];

CSX = DefineRectGrid(CSX, 0.001, mesh);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%MESH

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%mesh.x = [mesh.x linspace(-4,4,1.0)];
mesh.x = [mesh.x linspace(-substrate.width/2+ifa.Ledge,-substrate.width/2+ifa.Ledge+ifa.l,35.0) ];  % a-b=4.4; Step=

mesh.y = [mesh.y linspace((substrate.length/2)-ifa.e-ifa.sbt,substrate.width/2-ifa.e+ifa.h,14.0) ]; % a-b=2.5; Step=

mesh.z = [mesh.z linspace(-substrate.width*3,-1,45.0) ];
mesh.z = [mesh.z linspace(-1,0,4) ];
mesh.z = [mesh.z linspace(substrate.thickness,substrate.thickness+1,4.0) ];
mesh.z = [mesh.z linspace(substrate.thickness+1,substrate.width*3,45.0) ];

CSX = DefineRectGrid(CSX, 0.001, mesh);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% When third rule does't work, use high resolution meshing

%% finalize the mesh with detec edge DOESN'T WORK Hand meshing is better.

%ifa_mesh = DetectEdges(CSX, [], 'SetProperty','ifa');

%mesh.x = [mesh.x SmoothMeshLines(ifa_mesh.x, 0.5)];

%mesh.y = [mesh.y SmoothMeshLines(ifa_mesh.y, 10)];

% generate a smooth mesh with max. cell size: lambda_min / 20

%mesh = DetectEdges(CSX, mesh);

%mesh = SmoothMesh(mesh, c0 / (f0+fc) / unit / 20);

% Always problem to put smoothmeshing problem error: max_recursion_depth exceeded

%CSX = DefineRectGrid(CSX, unit , mesh); % define the mesh with a drawing unit of 1mm (1e-3)

%% add a nf2ff calc box; size is 3 cells away from MUR boundary condition
%%%%%%%It's implements a near field to far field transformation for near field antenna measurement setups
start = [mesh.x(4)     mesh.y(4)     mesh.z(4)];
stop  = [mesh.x(end-3) mesh.y(end-3) mesh.z(end-3)];
[CSX nf2ff] = CreateNF2FFBox(CSX, 'nf2ff', start, stop);

%% prepare simulation folder

Sim_Path = 'tmp_LIGHT';

Sim_CSX = 'IFA.xml';

try confirm_recursive_rmdir(false,'local'); end

[status, message, messageid] = rmdir( Sim_Path, 's' ); % clear previous directory

[status, message, messageid] = mkdir( Sim_Path ); % create empty simulation folder

%% write openEMS compatible xml-file

WriteOpenEMS( [Sim_Path '/' Sim_CSX], FDTD, CSX );

%% show the structure

if (show == 1)

CSXGeomPlot( [Sim_Path '/' Sim_CSX] );

end

%%%TO START SIMULATION AND FIND S11, IMPEDANCE, GAIN

%% run openEMS
RunOpenEMS( Sim_Path, Sim_CSX);  %RunOpenEMS( Sim_Path, Sim_CSX, '--debug-PEC -v');

%% postprocessing & do the plots
freq = linspace( max([1e9,f0-fc]), f0+fc, 501 );
port = calcPort(port, Sim_Path, freq);

Zin = port.uf.tot ./ port.if.tot;
s11 = port.uf.ref ./ port.uf.inc;
P_in = real(0.5 * port.uf.tot .* conj( port.if.tot )); % antenna feed power

%% NFFF contour plots %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%find resonance frequncy from s11
f_res_ind = find(s11==min(s11));
f_res = freq(f_res_ind);

% plot feed point impedance

figure
plot( freq/1e6, abs(Zin), 'k-', 'Linewidth', 2 );
hold on
plot( freq/1e6, real(Zin), 'g--', 'Linewidth', 2 );
hold on
grid on
plot( freq/1e6, imag(Zin), 'm--', 'Linewidth', 2 );
title( 'feed point impedance' );
xlabel( 'frequency f / MHz' );
ylabel( 'impedance Z_{in} / Ohm' );
legend( 'real', 'imag' );
disp( ['-> Zin @ ', num2str(f_res/1e9), 'GHz = ', num2str(Zin(f_res_ind)), ]);
disp( ['-> |Zin| @ ', num2str(f_res/1e9), 'GHz = ', num2str(abs(Zin(f_res_ind))), ' Ohm' ]);

% plot reflection coefficient S11
figure
plot( freq/1e6, 20*log10(abs(s11)), 'k-', 'Linewidth', 2 );
grid on
title( 'reflection coefficient S_{11}' );
xlabel( 'frequency f / MHz' );
ylabel( 'reflection coefficient |S_{11}|' );
disp( ['-> S11 @ ', num2str(f_res/1e9), 'GHz = ', num2str(20*log10(abs(min(s11)))), ' dB']);

drawnow

%% %%%% 3B dB
disp( 'calculating 3D far field pattern and dumping to vtk (use Paraview to visualize)...' );
thetaRange = (0:2:180);
phiRange = (0:2:360) - 180;
nf2ff = CalcNF2FF(nf2ff, Sim_Path, f_res, thetaRange*pi/180, phiRange*pi/180,'Verbose',1,'Outfile','3D_Pattern.h5');
%nf2ff = CalcNF2FF(nf2ff, Sim_Path, f_res, [-180:2:180]*pi/180, [0 90]*pi/180, 'Mode', 1, 'Center', [0 0 length])
%nf2ff = CalcNF2FF(nf2ff, Sim_Path, f0, [-180:2:180]*pi/180, [0 90]*pi/180, 'Mode', 1, [0 0 0], [0 0 length*unit]);

figure
plotFF3D(nf2ff)

%%%%%% display power and directivity
disp( ['radiated power: Prad = ' num2str(nf2ff.Prad) ' Watt']);
disp( ['directivity: Dmax = ' num2str(nf2ff.Dmax) ' (' num2str(10*log10(nf2ff.Dmax)) ' dBi)'] );
disp( ['efficiency: nu_rad = ' num2str(100*nf2ff.Prad./real(P_in(f_res_ind))) ' %']);

E_far_normalized = nf2ff.E_norm{1} / max(nf2ff.E_norm{1}(:)) * nf2ff.Dmax;
DumpFF2VTK([Sim_Path '/3D_Pattern.vtk'],E_far_normalized,thetaRange,phiRange,1e-3);

[thomaslepoix]

Hi Diane!

Here are the official guidelines to mesh a structure.

From what I see :

• Your mesh does not comply to the thirds rule, nor to the preconized smooth factor.
• The ports require to be traversed by a meshline in each dimension, if possible at its center in X and Y axis.
• The resolution seems to me a little loose but I think applying thirds rule correctly will lead to a tighter mesh.

I have to say OpenEMS meshing is a pain since it has to be done manually, there are tools out there that can help but never automagically.

The most out of the box solution I can propose you works well with microstrip and lookalikes :

• Install Qucs (they are struggling to port it to a new version of a lib so the last version is quite outdated. If you are on Windows the last binary should be fine, if you are on Linux, use this AppImage)
• Draw your antenna using microstrip components and setup an S parameters simulation.
• Use Qucs-RFlayout to convert that to an OpenEMS script.
• Use this script as a base for your work. You will have to resize the substrate, remove a part of the groundplane under the antenna, etc., and to refine the mesh as explained here.

The benefits of going this way are :

• The special meshlines (thirds rule, port's lines, etc.) are automatically placed, you will only need to resolve conflicting lines manually but your structure does not seems too conflictual to me (what I call conflicts is when special lines are really near, but not exactly at the same place, it prevents the OpenEMS smoothmeshing functions from producing a correct smoothmesh).
• The script will be already structured and embedding a command line interface to be easily tweakable, debuggable.
• It will also embed postprocessing features, including far field calculations.
• It will save you from learning OpenEMS by falling in all its pitfalls.

The produced script is based on the octave-openems-hll Octave/Matlab library, take a look to check its capabilities.