% cutoff_fem_rect.f
% Cut off wave number of a rectangular waveguide for TM11 mode.
% a = 1m, b = 0.5m
% FEM is used with 8 rectangular elements and 15 global nodes
clear all;
clc;
clf;
er = 1; %dielectric filling of waveguide
a0=1; % width of the waveguide
b0=0.5; % height of the waveguide
N0=25; % number of global nodes
Ne=16; % number of elements
Nfix=16; % number of fixed nodes
Nfree=N0-Nfix; % number of free nodes, should be less than fixed nodes
% element dimensions
ae=a0/4;
be=b0/4;
LF=[7 8 9 12 13 14 17 18 19]; % list of free nodes
% rectangle connectivity table
% the sequence is according to the local node numbers 1,2,3,4
% the row number defines the element number
el = [1 2 7 6;2 3 8 7;3 4 9 8;4 5 10 9;6 7 12 11;7 8 13 12;8 9 14 13;9 10 15 14;11 12 17 16;12
13 18 17;13 14 19 18;14 15 20 19;16 17 22 21;17 18 23 22;18 19 24 23;19 20 25 24];
% initialize the global matrix
Ag=zeros(N0); Bg=zeros(N0);
% Elemental matrix elements
A11=(1/3)*((be/ae)+(ae/be)); B11=er*(ae*be)/9;
A12=(1/3)*((ae/(2*be))-(be/ae)); B12=er*(ae*be)/18;
A14=(1/3)*((be/(2*ae))-(ae/be)); B14=er*(ae*be)/18;
A13=-(1/6)*((be/ae)+(ae/be)); B13=er*(ae*be)/36;
Ae=[A11 A12 A13 A14; A12 A11 A14 A13; A13 A14 A11 A12; A14 A13 A12 A11];
Be=[B11 B12 B13 B14; B12 B11 B14 B13; B13 B14 B11 B12; B14 B13 B12 B11];
%Global matrix
for e=1:Ne
gn=el(e,:); % reads the global node number for each element
A=zeros(N0);
B=zeros(N0);
A(gn,gn)=Ae; % forming augmented(global) matrix out of elemental matrix
B(gn,gn)=Be;
Ag=Ag+A; %adding augmented matrices for assembly
Bg=Bg+B;
end
% System coefficient matrix defined by the free nodes
Agfree(1:Nfree,1:Nfree)=Ag(LF(1:Nfree),LF(1:Nfree))
Bgfree(1:Nfree,1:Nfree)=Bg(LF(1:Nfree),LF(1:Nfree));
Aij=inv(Bgfree)*Agfree;
Kc= eig(Aij).^0.5; % minimum value of cut-off wavenumber
K = [];
for i=1:Nfree % finds the nonzero eigenvalues
if round(Kc(i))>0
K=[K Kc(i)];
end
end
disp('FEM value of cut-off wave number = ');
disp(min(K));
K_ana = sqrt((pi/a0)^2+(pi/b0)^2); % analytical cutoff wavenumber for 11 mode
%K_ana=pi/a0; % analytical cutoff wavenumber for 10 mode
disp('Analytical value of cut-off wave number = ');
disp(K_ana);
err=100*(min(K)-K_ana)/K_ana;
disp('% Error');
disp(err);