ZIP with MATLAB scripts and note:
exercise 6.10
exercise 6.10 notes:
RECTANGULAR CAVITY
STATIONARY FIELDS
STORED
ENERGY
for any resonant structure
DISSIPATED POWER ON CONDUCTOR
Q FACTOR FROM CONDUCTOR
DISSIPATED POWER IN DIELECTRIC
Q FACTOR FROM DIELECTRIC
TOTAL (UNLOADED) Q0
RELATING STORED inside SURFACE S, RADIATED (or RECEIVED) THROUGH SURFACE S, and DISSIPATED (LOSS) in VOLUME V and on S
TEmn TMmn RECTANGULAR WAVEGUIDE PARAMETERS:
When resonating or transmitting device z size may not be a multiple of lambda/2, wave travelling left-to-right. If waveguide load or/and generator not impedance matched, then there is an additional wave term exp(+1j*z)
mu0=1.256637061435917e-06
eps0=8.854e-12
etha0=120*pi
c0=1/(mu0*eps0)^.5 % light speed
mur=1;er=2.1
etha1=etha0/er^.5
a=.02;b=.01;d=.02 % rectangular cavity dimensions
dx=.01*a % spatial resolution
a_range=[0:dx:a];b_range=[0:dx:b];d_range=[0:dx:d];
[X,Y,Z]=meshgrid(a_range,b_range,d_range);
n=1;m=1;p=1; % TM n=1 m=1 p=1
mu0 = 1.256637061435917e-06
eps0 = 8.853999999999999e-12
etha0 = 3.769911184307751e+02
c0 = 2.997956376932163e+08
er = 2.100000000000000
etha1 = 2.601485870070325e+02
d = 0.020000000000000
dx = 2.000000000000000e-04
If I chose for instance 5GHz
f0=5e9 % [Hz]
lambda0=c0/f0
k0=2*pi/lambda0 % wave number
lambda1=lambda0/er^.5
k1=2*pi/lambda1
beta=(k1^2-(m*pi/a)^2-(n*pi/b)^2)^.5
kcxy=((m*pi/a)^2+(n*pi/b)^2)^.5
f0 = 5.000000000000000e+09
lambda0 = 0.059959127538643
k0 = 1.047911396497575e+02
lambda1 = 0.041375728882633
k1 = 1.518567884327186e+02
beta = 0.000000000000000e+00 + 3.167168622295275e+02i
kcxy = 3.512407365520363e+02
Complex beta, 'pushing' the wave down a waveguide in transmission mode reaches a lot further and is far easier on battery consumption compared to when in cut-off.
If I pull up to for instance 12GHz.
f0=12e9 % [Hz]
lambda0=c0/f0
k0=2*pi/lambda0 % wave number
lambda1=lambda0/er^.5
k1=2*pi/lambda1
beta=(k1^2-(m*pi/a)^2-(n*pi/b)^2)^.5
kcxy=((m*pi/a)^2+(n*pi/b)^2)^.5
E0=1; % [V/m] test value
Ztm=beta*etha1/k1
mu0 = 1.256637061435917e-06
eps0 = 8.853999999999999e-12
etha0 = 3.769911184307751e+02
c0 = 2.997956376932163e+08
er = 2.100000000000000
etha1 = 2.601485870070325e+02
d = 0.020000000000000
dx = 2.000000000000000e-04
TM11x travelling wave in rectangular guide has field expressions (3.101 a/b/c/d) on cutting wave, metal plates on both sides of waveguide (assume correct feed ..) one stores energy in stationary waves generating We and Wm,
Ex=-1j*beta*m*pi/(a*kc^2)*E0/Ztm*cos(m*pi*X/a).*sin(n*pi*Y/b).*sin(p*pi/d);
Ey=-1j*beta*m*pi/(b*kc^2)*E0*sin(m*pi*X/a).*cos(n*pi*Y/b).*sin(p*pi/d);
Ez=E0*sin(m*pi*X/a).*sin(n*pi*Y/b).*sin(p*pi/d);
Hx=1j*2*pi*f0*eps0*er*n*pi/(b*kc^2)*E0*sin(m*pi*X/a).*cos(n*pi*Y/b).*cos(p*pi/d);
Hy=-1j*2*pi*f0*eps0*er*m*pi/(a*kc^2)*E0*cos(m*pi*X/a).*sin(n*pi*Y/b).*cos(p*pi/d);
Hz=zeros(size(Hx)); % otherwise it wouldn't be TM
E=(Ex.*conj(Ex)+Ey.*conj(Ey)+Ez.*conj(Ez)).^.5;
H=(Hx.*conj(Hx)+Hy.*conj(Hy)+Hz.*conj(Hz)).^.5;
We=sum(sum(sum(E)))
Wm=sum(sum(sum(H)))
We = 2.169007368554939e-10
Wm = 9.610363550946333e+02
But Wm >> We a device to electromagnetically resonate has We=Wm.
So we have to start with the correct f0, given the fact that nothing else but [a b d] sizes and materials determine resonant frequencies, for the given [a b d] TM111 can be exactly calculated:
f0_TM111=c0/(2*pi*(mur)^.5)*((m*pi/a)^2+(n*pi/b)^2+(p*pi/d)^2)^.5;
f0=f0_TM111
lambda0=c0/f0
k0=2*pi/lambda0 % wave number
lambda1=lambda0/er^.5
k1=2*pi/lambda1
f0 = 1.835865848651688e+10
lambda0 = 0.016329931618555
k0 = 3.847649490485592e+02
lambda1 = 0.011268723396380
k1 = 5.575773835391055e+02
where do We Wm cross
df=1e7
storeWe=[]
storeWm=[]
f_range=[10e9:df:21e9];
for s1=1:1:numel(f_range)
lambda0=c0/f_range(s1);
k0=2*pi/lambda0; % wave number
lambda1=lambda0/er^.5;
k1=2*pi/lambda1;
% a TM11x mode travelling wave in rectangular guide has field expressions (3.101 a/b/c/d)
kcxy=((m*pi/a)^2+(n*pi/b)^2)^.5;
beta=(k1^2-(m*pi/a)^2-(n*pi/b)^2)^.5;
Ztm=beta*etha1/k1;
f_range(s1);
Ex=(-2*1j)*(-1j)*beta*m*pi/(a*kcxy^2)*E0*cos(m*pi*X/a).*sin(n*pi*Y/b).*sin(beta*Z);
Ey=(-2*1j)*(-1j)*beta*n*pi/(b*kcxy^2)*E0*sin(m*pi*X/a).*cos(n*pi*Y/b).*sin(beta*Z);
Ez=(-1j*2)*E0*sin(m*pi*X/a).*sin(n*pi*Y/b).*sin(beta*Z);
Hx=(-2)*1j*2*pi*f_range(s1)*eps0*er*n*pi/(b*kcxy^2)*E0*sin(m*pi*X/a).*cos(n*pi*Y/b).*cos(beta*Z);
Hy=(-2)*(-1j)*2*pi*f_range(s1)*eps0*er*m*pi/(a*kcxy^2)*E0*cos(m*pi*X/a).*sin(n*pi*Y/b).*cos(beta*Z);
Hz=zeros(size(Hx)); % otherwise it wouldn't be TM
E2=(Ex.*conj(Ex)+Ey.*conj(Ey)+Ez.*conj(Ez));
H2=(Hx.*conj(Hx)+Hy.*conj(Hy)+Hz.*conj(Hz));
We=er*eps0/4*sum(sum(sum(E2)))*a*b*d;
Wm=mu0/4*sum(sum(sum(H2)))*a*b*d;
storeWe=[storeWe We];
storeWm=[storeWm Wm];
end
figure;
plot(f_range,storeWe)
hold on
plot(f_range,storeWm);grid on
ax21=gca
legend({'We','Wm'})
zooming on
We1=We([70:end]);Wm1=We([70:end]);
[f_int,Wint]= ...
intersections(f_range,storeWm,f_range,storeWe)
f_int =
1.0e+10 *
1.092784337414282
1.186029940547057
1.261769656709135
1.397966131817036
1.538079082900475
1.742384794882003
1.915024459913416
Wint =
1.0e-10 *
0.518815395505571
0.050295951414071
0.057039778144585
0.070215676707580
0.085177740516018
0.109552914062239
0.132516196060094
plot(f_int,Wint,'rs')
If attempting simplification as done for TE101: all that is not sin or cos, to different E0', and then applying Ztm to get H0 out of E0
df=1e8
f_range=[10e9:df:21e9];
storeWe=[]
storeWm=[]
for s1=1:1:numel(f_range)
lambda0=c0/f_range(s1);
k0=2*pi/lambda0; % wave number
lambda1=lambda0/er^.5;
k1=2*pi/lambda1;
kcxy=((m*pi/a)^2+(n*pi/b)^2)^.5;
beta=(k1^2-(m*pi/a)^2-(n*pi/b)^2)^.5;
Ztm=beta*etha1/k1
f_range(s1)
Ex=E0*cos(m*pi*X/a).*sin(n*pi*Y/b).*sin(beta*Z);
Ey=E0*sin(m*pi*X/a).*cos(n*pi*Y/b).*sin(beta*Z);
Ez=E0*sin(m*pi*X/a).*sin(n*pi*Y/b).*sin(beta*Z);
Hx=E0/Ztm*sin(m*pi*X/a).*cos(n*pi*Y/b).*sin(beta*Z);
Hy=E0/Ztm*cos(m*pi*X/a).*sin(n*pi*Y/b).*sin(beta*Z);
Hz=zeros(size(Hx)); % otherwise it wouldn't be TM
E2=(Ex.*conj(Ex)+Ey.*conj(Ey)+Ez.*conj(Ez));
H2=(Hx.*conj(Hx)+Hy.*conj(Hy)+Hz.*conj(Hz));
We=eps0*er/4*sum(sum(sum(E2)))*a*b*d;
Wm=mu0/4*sum(sum(sum(H2)))*a*b*d;
storeWe=[storeWe We];
storeWm=[storeWm Wm];
end
We Wm do not cross anywhere near 12GHz, at least with the field expressions supplied by [POZAR] here used.
figure;
plot(f_range,storeWe)
hold on
plot(f_range,storeWm);
legend({'We','Wm'});grid on
Where is the electric charge? Maxwell 1st equation: Gauss
divE=divergence(X,Y,Z,Ex,Ey,Ez);
figure;
colormap('jet')
h=slice(X,Y,Z,abs(divE),[a_range(1) a_range(end)],b_range(1),[d_range(1) d_range(end)])
shading interp;
daspect([1 1 1]);
axis tight;
camlight;
set([h(1),h(2)],'ambientstrength',.6);
campos([0.0900 0.1200 0.0700]);
title('TM111 charge distribution');
xlabel('x [0 a]');ylabel('y [0 b]');zlabel('z [0 d]');
if lossy dielectric, then, as done in chapter 3, we could shift an XY surface along Z axis, and each plane would show certain charge distribution inside the dielectric, but because the question states lossless dielectric, then Q0 = Qc.
For any arbitrary metal shape storing electromagnetic energy one could resort to symbolically solving 3D integrals hoping to obtain more or less simple expressions, like the one supplied in the solutions manual for this particular TM111, or in [POZAR] pg286 (6.43a) for TE10x.
Following, pg286 (6.43a) for TE10x We Wm expressions are NOT SAME MODE ! FOR TM111: DON'T USE THESE EXPRESSIONS TO SOLVE 6.10
However let's see what is the result of using such clever formula directly.
f0_TE101=c0/(2*pi*(mur)^.5)*...
((1*pi/a)^2+(0*pi/b)^2+(1*pi/d)^2)^.5
lambda0=c0/f0_TE101;
k0=2*pi/lambda0; % wave number
lambda1=lambda0/er^.5;
k1_=2*pi/lambda1;
kcxy_=((m*pi/a)^2+(n*pi/b)^2)^.5;
beta_=(k1^2-(m*pi/a)^2-(n*pi/b)^2)^.5;
E1=(-2*1j)*(-1j)*beta_*m*pi/(a*kcxy_^2)*E0
We_TE101=eps0*er*a*b*d/16*E1^2;
% this E1 wouldn't not same E0,
% but E1=(-2*1j)*(-1j)*beta*m*pi/(a*kcxy^2)*E0
Z_TE=k1_*etha1/beta_
Wm_TE101=eps0*er*a*b*d/16*E1^2/etha1
% Z_TE=k1*etha/beta not the same as Z_TM=beta*etha/k1
f0_TE101 = 1.059937641915093e+10
E1 = -1.355666846843143
Z_TE = 1.573088414223621e+02
2 orders of magnitude mismatch between We_TE101 Wm_TE101 that are supposed to be equal at resonating frequency.
f0_TE101
% so for TM111 the solutions manual expression for Wm is
A0=E0
a*b*d*mu0*(A0)^2*1/32*(1/a^2+1/b^2)
% closest Wint to 12GHz is Win(2),
%although not exactly at 12GHz has the following value
Wint(2)
A0=E1
a*b*d*mu0*(A0)^2*1/32*(1/a^2+1/b^2)
f0_TE101 = 1.059937641915093e+10
A0 = 1
= 1.963495408493621e-09
= 5.029595141407139e-12
A0 = -1.355666846843143
= 3.608575870952672e-09
Using table 3.2 field expressions seem more accurate than the simplification in chapter 6 for TE101 calculations of We Wm .
Now let's calculate Pc dissipated power on metal sheets:
sigma_Cu=5.813e7 % [S/m] 90% Cu 10% Ni
Rs=(pi*mu0/sigma_Cu)^.5
sigma_Cu = 58130000
Rs = 2.606031776056232e-07
For now assume amount repeated points among adjacent surfaces of volume V very small in comparison to total amount of points in envelope S type cell containing metal surfaces. Next version [S nS]=peel(V) where nS is going to contain indices to surface points, and resulting S does not contain repeated surface points among adjacent surfaces.
V=reshape([1:numel(X)],numel(a_range),numel(b_range),numel(d_range));
S=peel(V)
S =
6×1 cell array
{ 51×101 double}
{ 51×101 double}
{101×101 double}
{101×101 double}
{101×51 double}
{101×51 double}
I wrote function peel so eventually have a general solution that can be applied to diverse problems, not just rectangular and cylindrical wave guides.
But from the visualisation of the charge, the top plate removed, but if one checks, there's only charge where the figure shows, on z=0 and z=range_d(end) coloured surfaces.
Z=0
Hx1=(-2)*1j*2*pi*f_range(s1)*eps0*er*n*pi/(b*kcxy^2)*E0*...
sin(m*pi*X/a).*cos(n*pi*Y/b).*...
cos(beta*Z);
Hy1=(-2)*(-1j)*2*pi*f_range(s1)*eps0*er*m*pi/(a*kcxy^2)*E0*...
cos(m*pi*X/a).*...
sin(n*pi*Y/b).*...
cos(beta*Z);
Hz1=zeros(size(Hx));
Pc1=Rs/2*sum(sum(Hx1*conj(Hy1)))
Z=d
Hx2=(-2)*1j*2*pi*f_range(s1)*eps0*er*n*pi/(b*kcxy^2)*E0*...
sin(m*pi*X/a).*...
cos(n*pi*Y/b).*...
cos(beta*Z);
Hy2=(-2)*(-1j)*2*pi*f_range(s1)*eps0*er*m*pi/(a*kcxy^2)*E0*...
cos(m*pi*X/a).*...
sin(n*pi*Y/b).*...
cos(beta*Z);
Hz2=zeros(size(Hx));
Pc2=Rs/2*sum(sum(Hx1*conj(Hy1)))
Pc1 = 1.657315870636257e-07
Pc2 = 1.929389401953894e-08
Also, instead of using f0=12e9 let's use the frequency closest to 12GHz that really allows this device to store electric field and magnetic field energies equally.
Pc=Pc1+Pc2
Pc = 1.850254810831646e-07
Here one could either define the complete volume as 3D points where E and/or H exist, or define different volumes and surfaces, for instance the dielectric volume and the metal faces. Often it's easier to define the entire volume, 1st option, and then extract the surface as the metal container of the enclosed dielectric or air/vacuum.
Again let's check against the 'smart' f(a,b,d, .. ) expression (6.45) pg287
L=d
Rs*E0^2*lambda1^2/(8*etha1^2)*(L^2*a*b /d^2 + b*d/a^2 + L^2*a/(2*d) + d/(2*a))
Rs*E1^2*lambda1^2/(8*etha1^2)*(L^2*a*b /d^2 + b*d/a^2 + L^2*a/(2*d) + d/(2*a))
L = 0.020000000000000
= 1.834387796175174e-16
= 3.371297692173485e-16
solutions manual expression for Pc
= 4.560555608098407e-07
E0^2*Rs/4*(a^3*d+b^3*d+a^3*b+b^3*a)/(a^2*b^2)
definitely,on a couple metal surfaces total area 2*a*b
2*a*b
= 4.000000000000000e-04
in 4 cm2, with electron charge order -1.e-19, definitely there are more than 100 electrons the Q would be in turn
2*2*pi*f_int(2)*We/5.682595842627020e-17
= 9.129040637719834e+15
technology is currently, for GHz THz hitting Q values of 1e6, but not 1e16
Therefore we carry on with Pc obtained from We=Wm
Qc=2*2*pi*f_int(2)*We/Pc
Qd=Inf
Q0=1/(1/Qc+1/Qd)
Qc = 2.803756978303032e+06
Qd = Inf
Q0 = 2.803756978303032e+06
Qc expression (6.46) pg287
(k0*a*d)^3*b*etha0/(2*pi^2*Rs)*1/(2*L^2*a^3*b + 2*b*d^3 + L^2*a^3*d + a*d^3 )
Q0(a,b,d, ..) expression in solutions manual
k0*etha0/(4*Rs)*(a*b*d)*(a^2+b^2)/(a^3*d+b^3*d+a^3*b+b^3*a)
COMMENT 1: limitation of ind2sub
V=reshape([1:48],4,4,3)
this
indV = 1:48;
szV = [4 4 3];
[row,col,page] = ind2sub(szV,indV)
works ok, but the output expression to correctly collect sub-indices produced by ind2sub has to have the exact amount of fields, otherwise ind2sub produces 'merged' indices, for instance if
nV=ind2sub(szV,indV)
ind2sub does not
nV=[row,col,page]
it just produces
nV=indV
what was then the point of invoking ind2sub anyway
indV=[1:prod(szV)]
[nV1,nV2,nV3]=ind2sub(szV,indV)
so for varying szV and indV the call for ind2sub has to be done with evalin or eval .
and then, generating a string [nV1,nV2 .. ,nVN] to collect output of ind2sub with eval, is not going to work because nV1 do not have dimension > 1, nV-ith is just a linear vector, while what is needed is a matrix with same structure as V
V
but containing indices
nV1 =
Columns 1 through 10
1 2 3 4 1 2 3 4 1 2
Columns 11 through 20
3 4 1 2 3 4 1 2 3 4
Columns 21 through 30
1 2 3 4 1 2 3 4 1 2
Columns 31 through 40
3 4 1 2 3 4 1 2 3 4
Columns 41 through 48
1 2 3 4 1 2 3 4
nV2 =
Columns 1 through 10
1 1 1 1 2 2 2 2 3 3
Columns 11 through 20
3 3 4 4 4 4 1 1 1 1
Columns 21 through 30
2 2 2 2 3 3 3 3 4 4
Columns 31 through 40
4 4 1 1 1 1 2 2 2 2
Columns 41 through 48
3 3 3 3 4 4 4 4
nV3 =
Columns 1 through 10
1 1 1 1 1 1 1 1 1 1
Columns 11 through 20
1 1 1 1 1 1 2 2 2 2
Columns 21 through 30
2 2 2 2 2 2 2 2 2 2
Columns 31 through 40
2 2 3 3 3 3 3 3 3 3
Columns 41 through 48
3 3 3 3 3 3 3 3
V(:,:,1) =
1 5 9 13
2 6 10 14
3 7 11 15
4 8 12 16
V(:,:,2) =
17 21 25 29
18 22 26 30
19 23 27 31
20 24 28 32
V(:,:,3) =
33 37 41 45
34 38 42 46
35 39 43 47
36 40 44 48
COMMENT 2: limitation of attempting symbolic
[POZAR] starts stationary waves expressions
Hx=B/b*sin(pi*m*x/a)*cos(n*pi*y/b)*exp(-1j*beta*z)
Hy=B/a*cos(m*pi*x/a)*sin(n*pi*y/b)*exp(1j*beta*z)
The rectangular waveguide is then cut at lambda/2 multiple to have current peak on XY closing plane.
m=1;n=1
syms a_ b_ d_ x y z
Hx = @(x,y,z) A/b_*sin(pi*x/a_).*cos(pi*y/b_).*cos(pi*z/d_)
Hy = @(x,y,z) A/a_*cos(pi*x/a_).*sin(pi*y/b_).*cos(pi*z/d_)
prim_Wm=@(x,y,z) int(int(int(Hx^2+Hy^2,x),y),z)
Wm=prim_Wm(a,b,d)-prim_Wm(0,0,0)
the obtained symbolic function of the stored magnetic energy prim_Wm as function of the cavity dimensions and field peak amplitude is supposed to be the same as
prim_Wm =
(d*((A^2*(sin((2*pi*z)/d)/4 + (pi*z)/(2*d))*(2*b^2*x*y*pi^2 - 2*b^3*x*pi*cos((pi*y)/b)*sin((pi*y)/b)))/8 - (A^2*a*(sin((2*pi*z)/d)/4 + (pi*z)/(2*d))*(b^3*cos((pi*y)/b)*sin((2*pi*x)/a)*sin((pi*y)/b) - b^2*y*pi*sin((2*pi*x)/a)))/8))/(a^2*b^2*pi^3) + (A^2*d*(sin((2*pi*z)/d)/4 + (pi*z)/(2*d))*(2*x*y*pi^2 + 2*b*x*pi*cos((pi*y)/b)*sin((pi*y)/b)))/(8*b^2*pi^3) - (A^2*a*d*(sin((2*pi*z)/d)/4 + (pi*z)/(2*d))*(y*pi*sin((2*pi*x)/a) + b*cos((pi*y)/b)*sin((2*pi*x)/a)*sin((pi*y)/b)))/(8*b^2*pi^3)
cul-de-sac. The symbolic toolbox has improved a lot but sometimes, when getting to symbolic expressions that cannot be simplified or even matched to expected simple relations it's better to attempt numerical solving.
Power loss along cavity wall: integrating |H(on walls)|^2 across all walls where there is tangential magnetic field.
PL=Rs/2*..
again, after integration the expected result is
And the Q factor: Q=2*pi*f0*(We+Wm)/PL=2*2*pi*f0*Wm/PL
