Trying to simulate a finite dipole. #54
-
|
Hi everyone! I want to create a dipole with the following characteristics:
And output its impedance diagram and its current current distribution. I tried putting it together but I have some issues with the current distribution. I have the following script: The dipole seems fine, the impedance diagram is ugly, and the my current distribution does not exist. What am I doing wrong? |
Beta Was this translation helpful? Give feedback.
Replies: 2 comments 14 replies
-
|
Hi, there are 3 errors: Here is corrected script which give you right results: clc
clear all
close all
%% set up the simulation %%
postprocessing_only = 0;
physical_constants;
f = 1e9;
lambda = c0/f;
unit = 1e-3; % units in mm
dipole_length = lambda / 2 / unit;
dipole_orientation = 3; % 1,2,3: x,y,z
f_start = 5e8;
f_stop = 2e9;
delta_f = 1e7;
dipole_width = 2e-3;
%% prepare simulation folder %%
Sim_Path = 'tmp_simN1';
Sim_CSX = 'tmp_simN1.xml';
%% FDTD object initialization %%
FDTD = InitFDTD('NrTS', 20000);
FDTD = SetGaussExcite(FDTD, f_start, f_stop);
% apply PML-8 boundary conditions in all directions %
BC = {'PML_8' 'PML_8' 'PML_8' 'PML_8' 'PML_8' 'PML_8'};
FDTD = SetBoundaryCond(FDTD, BC);
%% setup CSXCAD geometry & mesh %%
CSX = InitCSX();
% create an equidistant mesh %
mesh.x = -dipole_length*10:dipole_length/2:dipole_length*10;
mesh.y = -dipole_length*10:dipole_length/2:dipole_length*10;
mesh.z = -dipole_length*10:dipole_length/2:dipole_length*10;
% dipole %
CSX = AddMetal(CSX, 'Dipole'); % create a perfect electric conductor (PEC)
%
% ORIGINAL SCRIPT CORRECTION: wire splitted into 2 parts and added air gap between them with size 2mm
%
CSX = AddBox(CSX, 'Dipole', 1, [0 0 -dipole_length/2], [0 0 -1]);
CSX = AddBox(CSX, 'Dipole', 1, [0 0 1], [0 0 dipole_length/2]);
% add space for PML %
mesh = AddPML( mesh, [8 8 8 8 8 8] );
%
% ORIGINAL SCRIPT CORRECTION: added lines for input port in z +1mm -1mm
%
mesh.z = [mesh.z -1 1];
% define the mesh %
CSX = DefineRectGrid(CSX, unit, mesh);
%% define dump box %%
CSX = AddDump(CSX, 'distribution', 'DumpType', 2, 'FileType', 1);
start = [-dipole_width, -dipole_width, -dipole_length/2];
stop = [dipole_width, dipole_width, dipole_length/2];
CSX = AddBox(CSX, 'distribution', 0, start, stop);
%% setup feeding %%
feed.R = 50; %feed resistance
%% port with a 50 ohms impedance %%
start = [-1 -1 -1];
stop = [1 1 1];
%
% ORIGINAL SCRIPT CORRECTION: excitation port defined as box sitting in air gap and have some XYZ dimensions
%
[CSX port] = AddLumpedPort(CSX, 5 ,1 ,feed.R, start, stop, [0 0 1], true);
%
% Also something like this could be used
%
%[CSX port] = AddCurvePort( CSX, 5, 1, feed.R, [0 0 -1], [0 0 1], true)
%% launching the simulation %%
if (postprocessing_only==0)
[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
CSXGeomPlot([Sim_Path '/' Sim_CSX]);
% run openEMS
RunOpenEMS(Sim_Path, Sim_CSX);
endif
%% post processing %%
freq = f_start:delta_f:f_stop;
port = calcPort(port, Sim_Path, freq);
% feed-point impedance diagram %
Zin = port.uf.tot ./ port.if.tot;
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' );
% current distribution %
figure
dump_file = [Sim_Path '/distribution.h5'];
PlotArgs.slice = {dipole_length/2*unit dipole_width/2*unit 0};
PlotArgs.pauseTime=0.01;
PlotArgs.component=0;
PlotArgs.Limit = 'auto';
PlotHDF5FieldData(dump_file, PlotArgs); |
Beta Was this translation helpful? Give feedback.
-
|
4nec results for comparison, most confusing for me that impedance graphs are different after 1GHz, maybe it's that created model in openEMS has some mistakes, S11 is correct but impedances are different so I'm not really sure if the openEMS model of dipole is right :/ maybe somebody else will know and can answer here |
Beta Was this translation helpful? Give feedback.
-
|
and here is how dipole is meshed |
Beta Was this translation helpful? Give feedback.
-
|
That's so impressive. I also tried solving the issue but I got nowhere near any of this. My idea was to look for a way to plot the incident current from the port over the surface of the dipole. Needless to say I ran into a wall pretty quickly after looking at the function reference. I'm new to CEM, so I've been struggling to understand how it all work. Thanks a bunch, really. |
Beta Was this translation helpful? Give feedback.
-
|
Hi, I was yesterday thinking how to plot that current and found easy solution: So here is final dipole simulation with current distribution in dipole good new is that that current distribution graph is similar to current density what it should be. clc
clear all
close all
%% set up the simulation
postprocessing_only = 0;
physical_constants;
f = 1e9;
lambda = c0/f;
unit = 1e-3; % units in mm
dipole_length = lambda / 2 / unit;
dipole_orientation = 3; % 1,2,3: x,y,z
f_start = 500e6;
f_stop = 2000e6;
delta_f = 10e6;
dipole_width = 2e-3;
%%prepare simulation folder
Sim_Path = 'tmp_simN1';
Sim_CSX = 'tmp_simN1.xml';
%%FDTD object initialization
%% CHANGE:
%% STOP SIMULATION AT -30dB to got fast simulation 'EndCriteria'
%% 0.001 = -30dB
%% 0.0001 = -40dB
%% 0.00001 = -50dB
%% ...
%%
FDTD = InitFDTD('NrTS', 300000, 'EndCriteria', 0.001);
FDTD = SetGaussExcite(FDTD, f_start, f_stop);
%f_stop = f_start;
%FDTD = SetSinusExcite(FDTD, f_start);
% apply PML-8 boundary conditions in all directions %
BC = {'PML_8' 'PML_8' 'PML_8' 'PML_8' 'PML_8' 'PML_8'};
FDTD = SetBoundaryCond(FDTD, BC);
%%setup CSXCAD geometry & mesh
CSX = InitCSX();
% create an equidistant mesh %
mesh.x = -dipole_length*10 : dipole_length/2 : dipole_length*10;
mesh.y = -dipole_length*10 : dipole_length/2 : dipole_length*10;
mesh.z = -dipole_length*10 : dipole_length/2 : dipole_length*10;
% dipole %
CSX = AddMetal(CSX, 'Dipole'); % create a perfect electric conductor (PEC)
% ORIGINAL SCRIPT CORRECTION: wire splitted into 2 parts and added air gap between them with size 2mm
%
CSX = AddBox(CSX, 'Dipole', 1, [-0.5 -0.5 -dipole_length/2], [0.5 0.5 -1]);
CSX = AddBox(CSX, 'Dipole', 1, [-0.5 -0.5 1], [0.5 0.5 dipole_length/2]);
%% add space for PML %
mesh = AddPML( mesh, [8 8 8 8 8 8] );
% ORIGINAL SCRIPT CORRECTION: added lines for input port in z +1mm -1mm
% added more fine mesh around dipole to be able plot fields
%
mesh.x(mesh.x >= -5 & mesh.x <= 5) = [];
mesh.y(mesh.y >= -5 & mesh.y <= 5) = [];
mesh.z(mesh.z >= -dipole_length & mesh.z <= dipole_length) = [];
mesh.x = [mesh.x -5:1:5 0.1 0.3 -0.1 -0.3];
mesh.y = [mesh.y -5:1:5 0.1 0.3 -0.1 -0.3];
mesh.z = [mesh.z -dipole_length:dipole_length/20:dipole_length];
mesh.z = [mesh.z -1 1]; %add mesh lines for air gap
% define the mesh %
CSX = DefineRectGrid(CSX, unit, mesh);
%%define dump box
CSX = AddDump(CSX, 'distributionJ','FileType', 1, 'DumpType', 13, 'Frequency', f);
start = [-20, -20, -dipole_length/2];
stop = [20, 20, dipole_length/2];
CSX = AddBox(CSX, 'distributionJ', 0, start, stop);
%
% Putting 10 current probes along dipole.
% - their names will be distributioCurrent_1, distributioCurrent_2, ...
%
a = 1
for k = -dipole_length/2 : dipole_length/10 : dipole_length/2
CSX = AddProbe(CSX, ['distributionCurrent_' num2str(a)], 1, 'frequency', f);
start = [-1, -1, k];
stop = [1, 1, k];
CSX = AddBox(CSX, ['distributionCurrent_' num2str(a)], 0, start, stop);
a = a + 1;
end
%% ORIGINAL SCRIPT CORRECTION: excitation port defined as box sitting in air gap and have some XYZ dimensions
%%
feed.R = 50;
start = [-0.5 -0.5 -1];
stop = [0.5 0.5 1];
[CSX port] = AddLumpedPort(CSX, 5 ,1 ,feed.R, start, stop, [0 0 1], true); %amplitude of E field will be 1
%%weight{1} = 0;
%%weight{2} = 0;
%%weight{3} = 100 %change amplitude to 100
%%CSX = SetExcitationWeight(CSX,'port_excite_1',weight);
%% launching the simulation
if (postprocessing_only==0)
[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
CSXGeomPlot([Sim_Path '/' Sim_CSX]);
% run openEMS
RunOpenEMS(Sim_Path, Sim_CSX);
endif
%% post processing
freq = f_start:delta_f:f_stop;
port = calcPort(port, Sim_Path, freq);
% feed-point impedance diagram %
Zin = port.uf.tot ./ port.if.tot;
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' );
figure
s11 = port.uf.ref./ port.uf.inc;
s11_dB = 20*log10(abs(s11));
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}|' );
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% PLOT CURRENT FROM PROBES
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure()
files = {
'distributionCurrent_1',
'distributionCurrent_2',
'distributionCurrent_3',
'distributionCurrent_4',
'distributionCurrent_5',
'distributionCurrent_6',
'distributionCurrent_7',
'distributionCurrent_8',
'distributionCurrent_9',
'distributionCurrent_10',
'distributionCurrent_11'
};
measuredCurrent = ReadUI(files, Sim_Path, 1e9);
yValues = []
for k = 1:1:11
currentValue = abs(measuredCurrent.FD{1,k}.val)*1e6
yValues = [yValues currentValue]
end
plot(yValues)
ylabel('current [uA]');
xlabel('probe [-]');
|
Beta Was this translation helpful? Give feedback.
-
|
It does look like it. I'm having a problem though I get the following message when running the script: It seems I need to write the files in the simulation folder is that correct? |
Beta Was this translation helpful? Give feedback.
-
|
aaah, sorry I'm looking I left this at the beginning of script: |
Beta Was this translation helpful? Give feedback.
-
|
Oh. My bad I didn't pay attention to that. |
Beta Was this translation helpful? Give feedback.
-
|
So I finaly got right simulation for dipole antenna. Before I provided answer but it wasn't quite correct as one of the target plots is current distribution along dipole and it must be same/similar to the one from MoM simulation from 4nec. As there is possible to specify in 4nec voltage source amplitude in middle of dipole let's say 1V current distribution must be same for results from FTDT and MoM simulation.
FTDT simulation at its core calculates E-field and H-field values over grid points which has some material parameters (permitivity, permeability, electric conductivity, magnetic conductivity). Values like voltage or current marked as U, I are scalar units that means there are measured between two points (voltage) or amount through surface (current). They can be calculated from E-field and H-field when taking some line (E-field -> voltage) and surface (H-field -> current). Now when you put into simulation current (plane) and voltage probe (line) openEMS data result is: In simulation for lumped port there is specified E-field amplitude in Z direction for active port with amplitude 1V/m and impedance 50Ohm, but voltage which you will measure in real world when put voltemeter on dipole input depends on it's physical size, if input port is 1mm you should measure U_meas = 1V/m * 1mm = 1mV. This leads you that to have simulation with real values you should put there instead 1V/m value 1V/mm = 1000V/m, but it's not so easy (I tried and tried but really it's not all) due input port has more gridlines in XY plane which excite all points inside excitation box with Ez = 1000V/m and also this gridlines don't have equidistance spacing, so changing Ez amplitude is not way how to get right value for input voltage and therefore to get right current distribution when Uin = 1V. Right way what I think is that whatever excitation Ez field value is set, simulation is calculated and E-field, H-field are computed so you will get right values for all Sxx params but not right values for U, I params due their values depends on probes placement. To get Uin as 1V you need to scale in postprocessing values from simulation. scaleFactor = 1 ./ U_middle;
for k = 2:1:numel(port)
ifVector = abs(port{k}.if.tot)*normFactor .* scaleFactor;
...This will provide you current along dipole for different frequencies. To get current distribution for specific frequency you need to get this value from each current array from same index: freqIndex = 32
yValues = [];
scaleFactor = 1 ./ U_middle(freqIndex);
for k = 2:1:numel(port)
ifVector = abs(port{k}.if.tot(freqIndex))*normFactor .* scaleFactor;
yValues = [yValues ifVector'];
end
yValues' * 1e3 %show current values as [mA]
plot(yValues);That's it, here is complete code:clc
clear all
close all
%% set up the simulation
postprocessing_only = 0;
opts = '';
physical_constants;
f = 920e6;
lambda = c0/f;
unit = 1e-3; % units in mm
dipole_length = lambda / 2 / unit;
dipole_orientation = 3; % 1,2,3: x,y,z
f_start = 500e6;
f_stop = 1200e6;
delta_f = 10e6;
dipole_width = 2e-3;
%%prepare simulation folder
Sim_Path = 'tmp_simN1';
Sim_CSX = 'tmp_simN1.xml';
%%FDTD object initialization
%% CHANGE:
%% STOP SIMULATION AT -30dB to got fast simulation 'EndCriteria'
%% 0.001 = -30dB
%% 0.0001 = -40dB
%% 0.00001 = -50dB
%% ...
%%
FDTD = InitFDTD('NrTS', 500000, 'EndCriteria', 0.001);
FDTD = SetGaussExcite(FDTD, f_start, f_stop);
% apply PML-8 boundary conditions in all directions %
BC = {'PML_8' 'PML_8' 'PML_8' 'PML_8' 'PML_8' 'PML_8'};
FDTD = SetBoundaryCond(FDTD, BC);
%%setup CSXCAD geometry & mesh
CSX = InitCSX();
% create an equidistant mesh %
mesh.x = -dipole_length*10 : dipole_length/2 : dipole_length*10;
mesh.y = -dipole_length*10 : dipole_length/2 : dipole_length*10;
mesh.z = -dipole_length*10 : dipole_length/2 : dipole_length*10;
% dipole %
CSX = AddMetal(CSX, 'Dipole'); % create a perfect electric conductor (PEC)
% ORIGINAL SCRIPT CORRECTION: wire splitted into 2 parts and added air gap between them with size 2mm
%
CSX = AddBox(CSX, 'Dipole', 1, [-0.5 -0.5 -dipole_length/2], [0.5 0.5 -0.5]);
CSX = AddBox(CSX, 'Dipole', 1, [-0.5 -0.5 0.5], [0.5 0.5 dipole_length/2]);
%% add space for PML %
mesh = AddPML( mesh, [8 8 8 8 8 8] );
% ORIGINAL SCRIPT CORRECTION: added lines for input port in z +1mm -1mm
% added more fine mesh around dipole to be able plot fields
%
mesh.x(mesh.x >= -5 & mesh.x <= 5) = [];
mesh.y(mesh.y >= -5 & mesh.y <= 5) = [];
mesh.z(mesh.z >= -dipole_length & mesh.z <= dipole_length) = [];
mesh.x = [mesh.x -5:1:5 0.1 0.3 0.4 -0.1 -0.3 -0.4];
mesh.y = [mesh.y -5:1:5 0.1 0.3 0.4 -0.1 -0.3 -0.4];
mesh.z = [mesh.z -dipole_length:dipole_length/20:dipole_length];
mesh.z = [mesh.z -1 -0.5 -0.4 -0.3 -0.2 -0.1 0.1 0.2 0.3 0.4 0.5 1]; %add mesh lines for air gap
% define the mesh %
CSX = DefineRectGrid(CSX, unit, mesh);
%% ORIGINAL SCRIPT CORRECTION: excitation port defined as box sitting in air gap and have some XYZ dimensions
%%
feed.R = 50;
delta1 = 0.0;
start = [-0.5+delta1 -0.5+delta1 -0.5];
stop = [0.5+delta1 0.5+delta1 0.5];
[CSX port{1}] = AddLumpedPort(CSX, 5 ,1 ,feed.R, start, stop, [0 0 1], true); %amplitude of E field will be 1
%
% Putting 10 current probes along dipole.
% - their names will be distributioCurrent_1, distributioCurrent_2, ...
%
b = 2
for k = -dipole_length/2 : dipole_length/10 : dipole_length/2
start = [-1, -1, k];
stop = [1, 1, k+1];
[CSX port{b}] = AddLumpedPort(CSX, 5, b, inf, start, stop, [0 0 1/0.001]);
b = b + 1;
end
%% launching the simulation
if (postprocessing_only==0)
[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
CSXGeomPlot([Sim_Path '/' Sim_CSX]);
% run openEMS
RunOpenEMS(Sim_Path, Sim_CSX, opts);
endif
%% post processing
freq = f_start:delta_f:f_stop;
port = calcPort(port, Sim_Path, freq);
% feed-point impedance diagram %
Zin = port{1}.uf.tot ./ port{1}.if.tot;
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' );
figure
s11 = port{1}.uf.ref./ port{1}.uf.inc;
s11_dB = 20*log10(abs(s11));
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}|' );
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% PLOT CURRENT FROM PROBES
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
calcPort(port, Sim_Path, f);
%values in freq domain must be unormn by max(port{k}.f)
normFactor = max(port{1}.f)
k = 1;
figure;
subplot(2,4,1);
plot(port{k}.it.tot);
title("Current in middle of dipole, Time Domain");
ylabel(['port ' num2str(k) ' current [A]']);
xlabel('time sample [-]');
##figure;
subplot(2,4,2);
plot(port{k}.ut.tot);
title("Voltage in middle of dipole, Time Domain");
ylabel(['port ' num2str(k) ' voltage [V]']);
xlabel('time sample [-]');
##figure;
subplot(2,4,3);
title("Current through probes, Freq Domain");
hold on;
grid on;
yValues = [];
for k = 2:1:numel(port)
currentValue = abs(port{k}.if.tot)*normFactor;
yValues = [yValues currentValue'];
end
plot(freq/1e6, yValues(:,1:floor(end/2)))
plot(freq/1e6, yValues(:,ceil(end/2):end), '--')
legend('2','3','4','5','6','7','8','9','10','11','12');
ylabel('current [A]');
xlabel('freq [MHz]');
##figure;
subplot(2,4,6);
hold on;
grid on;
yValues = [];
k = 1;
U_middle = abs(port{k}.uf.tot)*normFactor;
I_middle = abs(port{k}.if.tot)*normFactor;
plot(port{k}.f/1e6, U_middle);
title("Voltage in middle of dipole, Frequency Domain");
ylabel('voltage [V]');
xlabel('freq [MHz]');
##figure;
subplot(2,4,5);
title("Time Domain");
hold on;
for k = 2:1:numel(port)
plot(port{k}.it.tot);
end
legend('2','3','4','5','6','7','8','9','10','11','12');
ylabel(['port' ' current [A]']);
xlabel('time sample [-]');
##figure;
subplot(2,4,7);
title("Current through probes scaled for U_middle=1V, Freq Domain");
hold on;
grid on;
yValues = [];
scaleFactor = 1 ./ U_middle;
for k = 2:1:numel(port)
ifVector = abs(port{k}.if.tot)*normFactor .* scaleFactor;
yValues = [yValues ifVector'];
end
plot(freq/1e6, yValues(:,1:floor(end/2)))
plot(freq/1e6, yValues(:,ceil(end/2):end), '--')
legend('2','3','4','5','6','7','8','9','10','11','12');
ylabel('current [A]');
xlabel('freq [MHz]');
##figure;
subplot(2,4,8);
hold on;
grid on;
%for freqIndex = 10:1:35
for freqIndex = 32
yValues = [];
scaleFactor = 1 ./ U_middle(freqIndex);
for k = 2:1:numel(port)
ifVector = abs(port{k}.if.tot(freqIndex))*normFactor .* scaleFactor;
yValues = [yValues ifVector'];
end
yValues' * 1e3
plot(yValues);
end
title(["Current through probes at " num2str(freq(freqIndex)/1e6) "MHz, Freq Domain"]);
ylabel('current [A]');
xlabel('probe index [-]');Results
|
Beta Was this translation helpful? Give feedback.
-
|
Hi, sorry I couldn't answer earlier, I messed up my .profile file and it took a while to figure out how to revert the changes lol. |
Beta Was this translation helpful? Give feedback.
Hi, I was yesterday thinking how to plot that current and found easy solution:
- I put 10 current probes along dipole in frequency domain for 1GHz
- after simulation read current from each probe file and plot it in graph
So here is final dipole simulation with current distribution in dipole good new is that that current distribution graph is similar to current density what it should be.