Simulating a Meandered Inverted-F Antenna, and Polar Plot Questions

[FolkSong]

I've been working on a simulation of this antenna design by TI, which includes a good description of the design with measured results:
https://www.ti.com/lit/an/swra228c/swra228c.pdf

My octave script and results are below. The script was based on the inverted-F Tutorial code, modified for this design and frequency, and with some additional outputs.

Overall the simulation seems not bad. The resonance frequency is pretty much bang on (I set the length of the last antenna segment in between their values for 868 and 915 MHz, and my result was 883 MHz). The return loss is a bit low @ 5dB, when their measured result is 25+ dB. But this could be easily matched in the real world with discrete components, and there are various potential sources of error like the FR4 material properties and the other components on the PCB which I didn't model.

One thing I'm more concerned about is the accuracy of my polar plots. Some questions:

• I understand spherical coordinates but I extended theta to 360 degrees in order to get full circular plots when phi is the parameter. Is this valid? Is there any better way?
• TI separates their plots into horizonal and vertical polarization, is there a way to produce this in the simulator? Otherwise it's hard to compare against real measurements.
• My phi=0 plot is perfectly circular which seems unlikely, any ideas of what's happening there?

Any other feedback on the simulation would be welcome, I really don't know what I'm doing as far as meshes so I'm hoping the setup from the tutorial file is still reasonable. I did make the simulation box larger since I'm working at a lower frequency (900 MHz vs 2.4 GHz).

%
% EXAMPLE / antennas / inverted-f antenna (ifa)
%
% This example demonstrates how to:
%  - calculate the reflection coefficient of an ifa
%  - calculate farfield of an ifa
%
% Tested with
%  - Octave 3.7.5
%  - openEMS v0.0.30+ (git 10.07.2013)
%
% (C) 2013 Stefan Mahr <dac922@gmx.de>

close all
clear
clc
addpath('C:\openEMS\CTB');    %Circuit Toolbox (touchstone files etc)

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  ________________________________________________    __ substrate.
% | A                                              |\  __    thickness
% | |ifa.e         __________________________      | |
% | |             |   ___   _________________| .w  | |
% | |       ifa.h |  |   | |                       | |
% |_V_____________|__|___|_|_____________________ _| |
% |                .w1   .w2 \                      | |
% |                   |.fp|  \                     | |
% |                       |    feed point          | |
% |                       |                        | | substrate.length
% |<-substrate.leftwidth->|<-substrate.rightwidth->| |
% |                                                | |
% |________________________________________________| |
%  \________________________________________________\|
%
% See TI Design Note DN023 for details on meandered antenna
%
%
substrate.leftwidth  = 7;        % width of substrate left of feed
substrate.rightwidth  = 38;       % width of substrate right of feed
substrate.length = 52;             % length of substrate
substrate.thickness = 0.8;         % thickness of substrate
substrate.cells = 4;               % use 4 cells for meshing substrate

ifa.h  = 20;           % outside height of short circuit stub (L1)
ifa.w1 = 1;            % width of short circuit stub
ifa.w = 1;             % width of radiating element
ifa.w2 = 2;            % width of feed element
ifa.fp = 5;            % inside distance from feed element to stub (L2)
ifa.e  = 21;           % vertical distance to substrate edge

ifa.L3 = 4;         # inside length of horizontal segments
ifa.L4 = 10;        # inside length of vertical segments
ifa.L6 = 6;         # outside length of final vertical segment

% substrate setup (FR4)
substrate.epsR   = 4.4;
substrate.tand   = 0.018;
substrate.kappa  = substrate.tand * 2*pi*900e6 * EPS0*substrate.epsR;

%setup feeding
feed.R = 50;     %feed resistance

%open AppCSXCAD and show ifa
show = 1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% size of the simulation box
SimBox = [300 300 300]; %aim for ~1/2 wavelength clearance in each direction

%% setup FDTD parameter & excitation function
f0 = 900e6; % center frequency
fc = 300e6; % 20 dB corner frequency

FDTD = InitFDTD('NrTS',  60000 );
FDTD = SetGaussExcite( FDTD, f0, fc );
BC = {'MUR' 'MUR' 'MUR' 'MUR' 'MUR' 'MUR'}; % boundary conditions
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
% Define origin (x,y) as geometric center of board. Origin z as bottom of substrate.
CSX = AddMaterial( CSX, 'substrate');
CSX = SetMaterialProperty( CSX, 'substrate', 'Epsilon',substrate.epsR, 'Kappa', substrate.kappa);
start = [-(substrate.leftwidth+substrate.rightwidth)/2  -substrate.length/2   0                  ];
stop  = [ (substrate.leftwidth+substrate.rightwidth)/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.leftwidth+substrate.rightwidth)/2   -substrate.length/2         substrate.thickness];
stop  = [ (substrate.leftwidth+substrate.rightwidth)/2   substrate.length/2-ifa.e    substrate.thickness];
CSX = AddBox(CSX, 'groundplane', 10, start,stop);

%% create ifa
CSX = AddMetal( CSX, 'ifa' ); % create a perfect electric conductor (PEC)
tl = [substrate.leftwidth/2-substrate.rightwidth/2, substrate.length/2-ifa.e, substrate.thickness];  % translate origin to feed element
start = [0 0.5 0] + tl;
stop = start + [ifa.w2 ifa.h-0.5 0];
CSX = AddBox( CSX, 'ifa', 10,  start, stop);  % feed element
start = [-ifa.fp 0 0] + tl;
stop =  start + [-ifa.w1 ifa.h 0];
CSX = AddBox( CSX, 'ifa', 10,  start, stop);  % short circuit stub
start = [(-ifa.fp-ifa.w1) ifa.h 0] + tl;
stop = start + [ifa.fp+ifa.w1 -ifa.w 0];
CSX = AddBox( CSX, 'ifa', 10, start, stop);   % connect feed to stub

start = [0 ifa.h 0] + tl;
stop = start + [(ifa.w2+ifa.L3) -ifa.w 0];
CSX = AddBox( CSX, 'ifa', 10, start, stop);   % initial horizontal segment

#create k meandering segments
for k=1:3
  start = stop + [0 ifa.w 0];
  stop = start + [ifa.w (-ifa.L4-2*ifa.w) 0];
  CSX = AddBox( CSX, 'ifa', 10, start, stop);   % vertical segment going down

  start = stop;
  stop = start + [(ifa.L3+ifa.w) ifa.w 0];
  CSX = AddBox( CSX, 'ifa', 10, start, stop);   % bottom horizontal segment

  start = stop;
  stop = start + [-ifa.w (ifa.L4+ifa.w) 0];
  CSX = AddBox( CSX, 'ifa', 10, start, stop);   % vertical segment going up

  start = stop;
  stop = start + [(ifa.L3+ifa.w) -ifa.w 0];
  CSX = AddBox( CSX, 'ifa', 10, start, stop);   % top horizontal segment
end

start = stop + [0 ifa.w 0];
stop = start + [ifa.w (-ifa.L6) 0];
CSX = AddBox( CSX, 'ifa', 10, start, stop);  % final vertical segment (tuneable)

ifa_mesh = DetectEdges(CSX, [], 'SetProperty','ifa');
mesh.x = [mesh.x SmoothMeshLines(ifa_mesh.x, 0.5)];
mesh.y = [mesh.y SmoothMeshLines(ifa_mesh.y, 0.5)];

%% apply the excitation & resist as a current source
start = [0 0 0] + tl;
stop  = start + [ifa.w2 0.5 0];
[CSX port] = AddLumpedPort(CSX, 5 ,1 ,feed.R, start, stop, [0 1 0], true);

%% finalize the mesh
% generate a smooth mesh with max. cell size: lambda_min / 20
mesh = DetectEdges(CSX, mesh);
mesh = SmoothMesh(mesh, c0 / (f0+fc) / unit / 20);
CSX = DefineRectGrid(CSX, unit, mesh);

%% add a nf2ff calc box; size is 3 cells away from MUR boundary condition
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_IFA';
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


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

%% postprocessing & do the plots
freq = linspace( 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

% plot feed point impedance
figure
plot( freq/1e6, real(Zin), 'k-', 'Linewidth', 2 );
hold on
grid on
plot( freq/1e6, imag(Zin), 'r-', 'Linewidth', 2 );
title( 'feed point impedance' );
xlabel( 'frequency f / MHz' );
ylabel( 'impedance Z_{in} / Ohm' );
legend( 'real', 'imag' );

% 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}|' );

drawnow

%% export sparameter to touchstone file
##write_touchstone('s',freq,shiftdim(s11,-1),['export.s1p']);

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

%%
%disp( 'calculating 3D far field pattern and dumping to vtk (use Paraview to visualize)...' );
thetaRange = (0:2:360); %extend theta to 360 to get full polar plots
phiRange = (0:2:360) - 180;

nf2ff = CalcNF2FF(nf2ff, Sim_Path, f_res, thetaRange*pi/180, phiRange*pi/180,'Verbose',1,'Outfile','3D_Pattern.h5');

%% Resonance Mismatch and Efficiency Calculations
% get accepted power into port (this can be reduced by mismatch)
P_accepted = port.P_acc(f_res_ind);
% get incident power
P_incident = port.P_inc(f_res_ind);
% loss from mismatch
mismatch_loss = P_accepted / P_incident;
mismatch_loss_dB = 10*log10(mismatch_loss);

% display key results
disp( ['%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%']);
disp( ['Resonant Frequency: f_res = ' num2str(f_res/1e6) ' MHz']);

#disp( ['Radiated Power @ Resonance: ' num2str(nf2ff.Prad(1)) ' Watt']);
disp( ['Max Directivity @ Resonance: Dmax = ' num2str(nf2ff.Dmax(1)) ' (' num2str(10*log10(nf2ff.Dmax(1))) ' dBi)'] );
disp( ['Efficiency @ Resonance: ' num2str(100*nf2ff.Prad(1)./real(P_in(f_res_ind))) ' %']);
imp_res = 50*(1+s11(f_res_ind))/(1-s11(f_res_ind));
disp( ['Return Loss @ Resonance: ' num2str(20*log10(abs(s11(f_res_ind)))) ' dB (' num2str(imp_res) ')']) ;

% polar plot
figure
subplot(2,2,1)
polarFF(nf2ff,'xaxis','phi','param',[46],'logscale',[-20 10], 'xtics', 5);
subplot(2,2,2)
polarFF(nf2ff,'xaxis','theta','param',[91],'logscale',[-20 10], 'xtics', 5);
subplot(2,2,3)
polarFF(nf2ff,'xaxis','theta','param',[136],'logscale',[-20 10], 'xtics', 5);
drawnow

%figure
%plotFF3D(nf2ff)
%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);

[FolkSong]

Adding model image for good measure:

[VolkerMuehlhaus]

@FolkSong

I understand spherical coordinates but I extended theta to 360 degrees in order to get full circular plots when phi is the parameter. Is this valid? Is there any better way?

Yes, this is valid.

TI separates their plots into horizonal and vertical polarization, is there a way to produce this in the simulator? Otherwise it's hard to compare against real measurements.

The plot options using polarFF are somewhat limited. But you can always look at the source code and create your own evaluation for specific polarization instead of total field. We recently had a thread on that topic for the Python interface of openEMS.

For sweep in the xy plane (phi sweep at theta=90°), horizontal polarization is the Ephi component, and vertical polarisation is the Etheta component. For sweep in the xz and yz planes (theta sweep at phi=0° and phi=90°) it is more complicated to get horizontal and vertical component from the polar data.

My phi=0 plot is perfectly circular which seems unlikely, any ideas of what's happening there?

That result is indeed correct for that cutting plane, looking at absolute value over all polarizations.