matlab examples: cleaned up MSL
Sebastian Held
parent
896c7f21f3
commit
3205d31654
|
|
@ -1,91 +1,187 @@
|
|||
%
|
||||
% 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 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
abs_length = 250;
|
||||
length = 1000;
|
||||
width = 500;
|
||||
height = 250;
|
||||
%% 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;
|
||||
mesh_res = [5 5 10];
|
||||
|
||||
EPS0 = 8.85418781762e-12;
|
||||
MUE0 = 1.256637062e-6;
|
||||
|
||||
%% define openEMS options %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
openEMS_opts = '';
|
||||
% openEMS_opts = [openEMS_opts ' --disable-dumps'];
|
||||
% openEMS_opts = [openEMS_opts ' --debug-material'];
|
||||
|
||||
%% prepare simulation folder
|
||||
Sim_Path = 'tmp';
|
||||
Sim_CSX = 'msl.xml';
|
||||
|
||||
mkdir(Sim_Path);
|
||||
[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 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
FDTD = InitFDTD(5e5,1e-6);
|
||||
FDTD = SetGaussExcite(FDTD,0.5e9,0.5e9);
|
||||
BC = [1 1 0 1 0 0];
|
||||
FDTD = SetBoundaryCond(FDTD,BC);
|
||||
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 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
%% setup CSXCAD geometry & mesh
|
||||
% very simple mesh
|
||||
CSX = InitCSX();
|
||||
mesh.x = -width/2 : mesh_res(1) : width/2;
|
||||
mesh.y = [linspace(0,MSL_height,11) MSL_height+1 MSL_height+3 MSL_height+mesh_res(2):mesh_res(2):height];
|
||||
mesh.z = 0 : mesh_res(3) : length;
|
||||
CSX = DefineRectGrid(CSX, 1e-3,mesh);
|
||||
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 );
|
||||
|
||||
%%% MSL
|
||||
CSX = AddMaterial(CSX,'copper');
|
||||
CSX = SetMaterialProperty(CSX,'copper','Kappa',56e6);
|
||||
start = [-0.5*MSL_width, MSL_height , 0];stop = [0.5*MSL_width, MSL_height+1 , length];
|
||||
CSX = AddBox(CSX,'copper',0 ,start,stop);
|
||||
%% 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 );
|
||||
|
||||
start = [-0.5*MSL_width, 0 , 0];stop = [0.5*MSL_width, MSL_height , mesh_res(1)/2];
|
||||
CSX = AddExcitation(CSX,'excite',0,[0 -1 0]);
|
||||
CSX = AddBox(CSX,'excite',0 ,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 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
finalKappa = 0.3/abs_length^4;
|
||||
%% 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,'pml');
|
||||
CSX = SetMaterialProperty(CSX,'pml','Kappa',finalKappa);
|
||||
CSX = SetMaterialProperty(CSX,'pml','Sigma',finalSigma);
|
||||
CSX = SetMaterialWeight(CSX,'pml','Kappa',['pow(abs(z)-' num2str(length-abs_length) ',4)']);
|
||||
CSX = SetMaterialWeight(CSX,'pml','Sigma',['pow(abs(z)-' num2str(length-abs_length) ',4)']);
|
||||
start = [mesh.x(1) mesh.y(1) length-abs_length];
|
||||
stop = [mesh.x(end) mesh.y(end) length];
|
||||
CSX = AddBox(CSX,'pml',0 ,start,stop);
|
||||
|
||||
%% define dump boxes... %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
CSX = AddDump(CSX,'Et_','DumpMode',2);
|
||||
start = [mesh.x(1) , MSL_height/2 , mesh.z(1)];
|
||||
stop = [mesh.x(end) , MSL_height/2 , mesh.z(end)];
|
||||
CSX = AddBox(CSX,'Et_',0 , start,stop);
|
||||
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 );
|
||||
|
||||
CSX = AddDump(CSX,'Ht_','DumpType',1,'DumpMode',2);
|
||||
CSX = AddBox(CSX,'Ht_',0,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 boxes %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
%voltage calc
|
||||
CSX = AddProbe(CSX,'ut1',0);
|
||||
start = [ 0 MSL_height length/2 ];stop = [ 0 0 length/2 ];
|
||||
CSX = AddBox(CSX,'ut1', 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 );
|
||||
|
||||
%current calc
|
||||
CSX = AddProbe(CSX,'it1',1);
|
||||
start = [ -MSL_width MSL_height/2 length/2 ];stop = [ MSL_width MSL_height*1.5 length/2 ];
|
||||
CSX = AddBox(CSX,'it1', 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 compatoble xml-file %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
WriteOpenEMS([Sim_Path '/' Sim_CSX],FDTD,CSX);
|
||||
%% write openEMS compatible xml-file
|
||||
WriteOpenEMS( [Sim_Path '/' Sim_CSX], FDTD, CSX );
|
||||
|
||||
%% cd to working dir and run openEMS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
savePath = pwd();
|
||||
cd(Sim_Path); %cd to working dir
|
||||
args = [Sim_CSX ' ' openEMS_opts];
|
||||
invoke_openEMS(args);
|
||||
cd(savePath);
|
||||
%% 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
|
||||
|
|
|
|||
Loading…
Reference in New Issue