No solution for circular waveguide

General Discussion
1

I am using FreeFEM++ to simulate the cylindrical waveguide, but no reasonable solution is returned. Here I’m using the 94 GHz light (𝜆 = 3.17mm) as input, and radius of the waveguide is 1.27mm.

If I uncomment the savevtk lines to seed the solutions, nothing appeared:

image
image1200×1244 23.8 KB

Need I provide more information to help debugging?

Here is the code:

// heavily copied from https://github.com/FreeFem/FreeFem-sources/blob/master/examples/plugin/waveguide.edp

load "MUMPS"
load "iovtk"
load "medit"
load "msh3"
load "Element_Mixte3d"

// The boundary value problem:
// (sigma = 0, here k is the wavenumber) 
//  -k^2*E + curl(curl E) = 0   in Omega
//  E x n = 0   on the side
//  Curl(E) x n + i*beta n x (E x n) = Ginc   on z = 0
//  Curl(E) x n + i*beta n x (E x n) = 0      on z = c

bool debug = true;

real r    = 0.00127;
real d    = r;
int  nx   = 20;
int  nz   = 20;
real zmin = 0.0;
real zmax = 0.01;

// build the mesh
border C(t=0.0, 2.0*pi) {
    x=r*cos(t);
    y=r*sin(t);
    label=1;
};
mesh Thx = buildmesh(C(nx));

int [int] rup = [0, 1], rdown = [0, 2], rmid = [1, 3];
int guide = 3, in = 1, out = 2;
mesh3 Th = buildlayers(Thx, nz, zbound=[zmin, zmax]
        , coef = 1.
        , labelup = rup
        , labelmid = rmid
        , labeldown = rdown
        );
//medit("Cyl", Th, wait=debug);
//plot(Th, wait=debug);

cout << "labels(Th) = " << labels(Th) << endl;

real f  = 94*10^9;
real er = 1;
real c0 = 299792458;
real k  = 2*pi*f*sqrt(er) / c0;
int m = 1, n = 1;
real beta = sqrt(k^2 - (m*pi/d)^2 - (n*pi/d)^2);
real Z = sqrt(er) * 120 * pi;

real ukb = 1 / (k^2 - beta^2);
func expbz = exp(-1i*beta*z);
func ExTE = (1i*k*Z) * ukb * (n*pi) / d * cos(m*pi*x/d) * sin(n*pi*y/d) * expbz;
func EyTE = -(1i*k*Z) * ukb * (n*pi) / d * sin(m*pi*x/d) * cos(n*pi*y/d) * expbz;

func Gix = +2*1i*beta*ExTE;
func Giy = +2*1i*beta*EyTE;
func Giz = 0;


fespace Nh(Th, Edge13d);
Nh<complex> [Ex, Ey, Ez], [vx, vy, vz];

// Macros
macro Curl(ux,uy,uz) [dy(uz)-dz(uy),dz(ux)-dx(uz),dx(uy)-dy(ux)] // EOM
macro Nvec(ux,uy,uz) [uy*N.z-uz*N.y,uz*N.x-ux*N.z,ux*N.y-uy*N.x] // EOM //uxN 
macro Curlabs(ux,uy,uz) [abs(dy(uz)-dz(uy)),abs(dz(ux)-dx(uz)),abs(dx(uy)-dy(ux))] //EOM

problem waveguide([Ex, Ey, Ez], [vx, vy, vz], solver=sparsesolver) =
      int3d(Th)(Curl(Ex,Ey,Ez)' * Curl(vx,vy,vz))
    - int3d(Th)(k^2*[Ex,Ey,Ez]' * [vx, vy, vz])
    + int2d(Th,in,out)(1i * beta * Nvec(Ex,Ey,Ez)' * Nvec(vx,vy,vz))
    - int2d(Th,in)([vx,vy,vz]' * [Gix,Giy,Giz])
    + on(guide,Ex=0,Ey=0,Ez=0);

waveguide;

cout << "r = " << r << " Integrated Ex = " << int3d(Th)(abs(Ex)^2) << endl;
medit("real", Th, [real(Ex),real(Ey),real(Ez)]);
//cout << (abs(Ex)^2) << end;
//cout << (abs(Ey)^2) << end;
//cout << (abs(Ez)^2) << end;

//int[int] Order = [1];
//savevtk("result_x.vtu", Th, abs(real(Ex)), dataname="RealEx", order=Order);
//savevtk("result_y.vtu", Th, abs(real(Ey)), dataname="RealEy", order=Order);
//savevtk("result_z.vtu", Th, abs(real(Ez)), dataname="RealEz", order=Order);
//medit("real", Th, [real(Ex),real(Ey),real(Ez)]);

and the log

2

Can you upload the code? Copying from the post is difficult.
Thanks.

4

Thanks for your reply!

Here is the source file:
cylinder_waveguide.edp (2.6 KB)

Regards.

5

I don’t think this mode propogates and beta is real. If you make beta et al
complex it may run. I’m also not sure if the boundary condition is right on
the guide surface assuming it is a metal. I was curious so I hacked up the code
and did get eveescent field apparently. fwiw.

d and r are equal, is that intended? ( although if your line wrap is wrong it looks like r is zero lol ).

cylinder_waveguide.edp (3.1 KB)

exanes
exanes600×600 43.2 KB

diff cylinder_waveguide.edp*
2d1
< // modified mjm 
33d31
< //medit("3d", Thx, wait=debug);
51c49
< complex k  = 2*pi*f*sqrt(er) / c0;
---
> real k  = 2*pi*f*sqrt(er) / c0;
53c51
< complex beta = sqrt(k^2 - (m*pi/d)^2 - (n*pi/d)^2);
---
> real beta = sqrt(k^2 - (m*pi/d)^2 - (n*pi/d)^2);
56c54
< complex ukb = 1 / (k^2 - beta^2);
---
> real ukb = 1 / (k^2 - beta^2);
58,66d55
< cout<<" ukb= "<<ukb; 
< cout<<" k= "<<k; 
< cout<<" beta= "<<beta; 
< cout<<" d= "<<d; 
< cout<<" r= "<<r; 
< cout<<endl; cout.flush;
< cout<<endl; cout.flush;
< cout<<endl; cout.flush;
< 
70,73c59,60
< func gate=1.0; //  (x*x+y*y<.2*r*r)?1e-20:0 ; 
< real sc=1.0; // 1e-20;
< func Gix = +2*1i*beta*ExTE*gate*sc;
< func Giy =  +2*1i*beta*EyTE*gate*sc;
---
> func Gix = +2*1i*beta*ExTE;
> func Giy = +2*1i*beta*EyTE;
78d64
< //fespace Nh(Th, P13D);
87d72
< //problem waveguide([Ex, Ey, Ez], [vx, vy, vz], solver=CG) =
90,91c75
<     + int2d(Th,out)(1i * beta * Nvec(Ex,Ey,Ez)' * Nvec(vx,vy,vz))
<     - int2d(Th,in)(1i * beta * Nvec(Ex,Ey,Ez)' * Nvec(vx,vy,vz))
---
>     + int2d(Th,in,out)(1i * beta * Nvec(Ex,Ey,Ez)' * Nvec(vx,vy,vz))
93,95c77
<     - int2d(Th,guide)(Ex*N.x*vx+Ey*N.y*vy)
<     //+ on(guide,Ex=0,Ey=0,Ez=0);
<     + on(guide,Ez=0);
---
>     + on(guide,Ex=0,Ey=0,Ez=0);
6

Thank you very much for the modification, it works like a charm!

You’re right, changing beta et al complex makes it run.

There are some issues with the parameters: the light frequency is too low thus the wave is cut, which is shown in your figure. If I set f to higher frequencies, the wave passes through as expected.

Here is the test result:

image
image1958×1042 106 KB

image
image1426×922 78.7 KB

Also, I referred other simulations on waveguide, and found this Mathematica Demonstration very helpful. Set the parameters according to it, I can reproduce the result, which confirms that your modification is right.

As for the boundary conditions, I’m still a newbie in FEM field, thus the BCs are copied from the example of FreeFEM++. I tried to understand it:

image
image1302×514 68 KB

I don’t see significant differences of results between your BCs and the original’s, both of them result in similar wave propagation.

All in all, I am gratitude to you for helping me.

7

I guess I was not sure what the boundary wall BC was supposed to be and
I put in e DOT n insteadof E x n but IIRC the original was something like E is zero.
The ends are more issues of excitation and I guess the idea is to excite one mode.
I guess you could try to excite it with a plane wave and see if you can decompose
the stuff on the other end in terms of modes :slight_smile:

If it works with beta complex, they should maybe change the FF example .