# No solution for circular waveguide

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:

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

// 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

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

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

Regards.

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)

``````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);

``````

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:

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:

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.

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

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