1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
|
%
% EXAMPLE / microstrip / MSL
%
% Microstrip line on air "substrate" in z-direction.
%
% This example demonstrates:
% - simple microstrip geometry
% - characteristic impedance
% - material grading function
% - geometric priority concept
%
%
% Tested with
% - Matlab 2009b
% - Octave 3.3.52
% - openEMS v0.0.14
%
% (C) 2010 Thorsten Liebig <thorsten.liebig@uni-due.de>
close all
clear
clc
%% setup the simulation
physical_constants;
unit = 1e-3; % all length in mm
% geometry
abs_length = 100; % absorber length
length = 600;
width = 400;
height = 200;
MSL_width = 50;
MSL_height = 10;
%% prepare simulation folder
Sim_Path = 'tmp';
Sim_CSX = 'msl.xml';
[status, message, messageid] = rmdir( Sim_Path, 's' ); % clear previous directory
[status, message, messageid] = mkdir( Sim_Path ); % create empty simulation folder
%% setup FDTD parameter & excitation function %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
max_timesteps = 2000;
min_decrement = 1e-5; % equivalent to -50 dB
f0 = 2e9; % center frequency
fc = 1e9; % 10 dB corner frequency (in this case 1e9 Hz - 3e9 Hz)
FDTD = InitFDTD( max_timesteps, min_decrement );
FDTD = SetGaussExcite( FDTD, f0, fc );
BC = {'PMC' 'PMC' 'PEC' 'PMC' 'PEC' 'PEC'};
FDTD = SetBoundaryCond( FDTD, BC );
%% setup CSXCAD geometry & mesh
% very simple mesh
CSX = InitCSX();
resolution = c0/(f0+fc)/unit /15; % resolution of lambda/15
mesh.x = SmoothMeshLines( [-width/2, width/2, -MSL_width/2, MSL_width/2], resolution ); % create smooth lines from fixed lines
mesh.y = SmoothMeshLines( [linspace(0,MSL_height,5) MSL_height+1 height], resolution );
mesh.z = SmoothMeshLines( [0 length], resolution );
CSX = DefineRectGrid( CSX, unit, mesh );
%% create MSL
% attention! the skin effect is not simulated, because the MSL is
% discretized with only one cell!
CSX = AddMaterial( CSX, 'copper' );
CSX = SetMaterialProperty( CSX, 'copper', 'Kappa', 56e6 );
start = [-MSL_width/2, MSL_height, 0];
stop = [ MSL_width/2, MSL_height+1, length];
priority = 100; % the geometric priority is set to 100
CSX = AddBox( CSX, 'copper', priority, start, stop );
%% add excitation below the strip
start = [-MSL_width/2, 0 , mesh.z(1)];
stop = [ MSL_width/2, MSL_height, mesh.z(1)];
CSX = AddExcitation( CSX, 'excite', 0, [0 -1 0] );
CSX = AddBox( CSX, 'excite', 0, start, stop );
%% fake pml
% this "pml" is a normal material with graded losses
% electric and magnetic losses are related to give low reflection
% for normally incident TEM waves
finalKappa = 1/abs_length^2;
finalSigma = finalKappa*MUE0/EPS0;
CSX = AddMaterial( CSX, 'fakepml' );
CSX = SetMaterialProperty( CSX, 'fakepml', 'Kappa', finalKappa );
CSX = SetMaterialProperty( CSX, 'fakepml', 'Sigma', finalSigma );
CSX = SetMaterialWeight( CSX, 'fakepml', 'Kappa', ['pow(z-' num2str(length-abs_length) ',2)'] );
CSX = SetMaterialWeight( CSX, 'fakepml', 'Sigma', ['pow(z-' num2str(length-abs_length) ',2)'] );
start = [mesh.x(1) mesh.y(1) length-abs_length];
stop = [mesh.x(end) mesh.y(end) length];
% the geometric priority is set to 0, which is lower than the priority
% of the MSL, thus the MSL (copper) has precendence
priority = 0;
CSX = AddBox( CSX, 'fakepml', priority, start, stop );
%% define dump boxes
start = [mesh.x(1), MSL_height/2, mesh.z(1)];
stop = [mesh.x(end), MSL_height/2, mesh.z(end)];
CSX = AddDump( CSX, 'Et_', 'DumpMode', 2 ); % cell interpolated
CSX = AddBox( CSX, 'Et_', 0, start, stop );
CSX = AddDump( CSX, 'Ht_', 'DumpType', 1, 'DumpMode', 2 ); % cell interpolated
CSX = AddBox( CSX, 'Ht_', 0, start, stop );
%% define voltage calc box
% voltage calc boxes will automatically snap to the next mesh-line
CSX = AddProbe( CSX, 'ut1', 0 );
zidx = interp1( mesh.z, 1:numel(mesh.z), length/2, 'nearest' );
start = [0 MSL_height mesh.z(zidx)];
stop = [0 0 mesh.z(zidx)];
CSX = AddBox( CSX, 'ut1', 0, start, stop );
% add a second voltage probe to compensate space offset between voltage and
% current
CSX = AddProbe( CSX, 'ut2', 0 );
start = [0 MSL_height mesh.z(zidx+1)];
stop = [0 0 mesh.z(zidx+1)];
CSX = AddBox( CSX, 'ut2', 0, start, stop );
%% define current calc box
% current calc boxes will automatically snap to the next dual mesh-line
CSX = AddProbe( CSX, 'it1', 1 );
xidx1 = interp1( mesh.x, 1:numel(mesh.x), -MSL_width/2, 'nearest' );
xidx2 = interp1( mesh.x, 1:numel(mesh.x), MSL_width/2, 'nearest' );
xdelta = diff(mesh.x);
yidx1 = interp1( mesh.y, 1:numel(mesh.y), MSL_height, 'nearest' );
yidx2 = interp1( mesh.y, 1:numel(mesh.y), MSL_height+1, 'nearest' );
ydelta = diff(mesh.y);
zdelta = diff(mesh.z);
start = [mesh.x(xidx1)-xdelta(xidx1-1)/2, mesh.y(yidx1)-ydelta(yidx1-1)/2, mesh.z(zidx)+zdelta(zidx)/2];
stop = [mesh.x(xidx2)+xdelta(xidx2)/2, mesh.y(yidx2)+ydelta(yidx2)/2, mesh.z(zidx)+zdelta(zidx)/2];
CSX = AddBox( CSX, 'it1', 0, start, stop );
%% write openEMS compatible xml-file
WriteOpenEMS( [Sim_Path '/' Sim_CSX], FDTD, CSX );
%% show the structure
CSXGeomPlot( [Sim_Path '/' Sim_CSX] );
%% run openEMS
openEMS_opts = '';
openEMS_opts = [openEMS_opts ' --engine=fastest'];
% openEMS_opts = [openEMS_opts ' --debug-material'];
% openEMS_opts = [openEMS_opts ' --debug-boxes'];
RunOpenEMS( Sim_Path, Sim_CSX, openEMS_opts );
%% postprocess
freq = linspace( f0-fc, f0+fc, 501 );
U = ReadUI( {'ut1','ut2','et'}, 'tmp/', freq ); % time domain/freq domain voltage
I = ReadUI( 'it1', 'tmp/', freq ); % time domain/freq domain current (half time step offset is corrected)
% plot time domain voltage
figure
[ax,h1,h2] = plotyy( U.TD{1}.t/1e-9, U.TD{1}.val, U.TD{3}.t/1e-9, U.TD{3}.val );
set( h1, 'Linewidth', 2 );
set( h1, 'Color', [1 0 0] );
set( h2, 'Linewidth', 2 );
set( h2, 'Color', [0 0 0] );
grid on
title( 'time domain voltage' );
xlabel( 'time t / ns' );
ylabel( ax(1), 'voltage ut1 / V' );
ylabel( ax(2), 'voltage et / V' );
% now make the y-axis symmetric to y=0 (align zeros of y1 and y2)
y1 = ylim(ax(1));
y2 = ylim(ax(2));
ylim( ax(1), [-max(abs(y1)) max(abs(y1))] );
ylim( ax(2), [-max(abs(y2)) max(abs(y2))] );
% calculate characteristic impedance
% arithmetic mean of ut1 and ut2 -> voltage in the middle of ut1 and ut2
U = (U.FD{1}.val + U.FD{2}.val) / 2;
Z = U ./ I.FD{1}.val;
% plot characteristic impedance
figure
plot( freq/1e6, real(Z), 'k-', 'Linewidth', 2 );
hold on
grid on
plot( freq/1e6, imag(Z), 'r--', 'Linewidth', 2 );
title( 'characteristic impedance of MSL' );
xlabel( 'frequency f / MHz' );
ylabel( 'characteristic impedance Z / Ohm' );
legend( 'real', 'imag' );
%% visualize electric and magnetic fields
% you will find vtk dump files in the simulation folder (tmp/)
% use paraview to visualize them
|
