Correctness of solution for Helmholtz equation for open RF cavity?

General Discussion
1

I am studying a square RF cavity of side lengths L = 1 with a small aperture of radius a = 0.1 using the 2D Helmholtz equation:

(\nabla^2 + k^2)E = 0

I model the RF cavity walls as perfectly reflective, so I have Dirichlet boundaries E(x, y)\big|_{\partial \Omega_W} = E_0 = 1, where \Omega_W represents the semi-enclosed region inside the cavity. However, since I am interested in EM radiation propogating out of the aperture, I also use a first-order approximation to the Sommerfield radiation condition \left(\frac{\partial E}{\partial r} - ikE\right)_{\partial \Omega_C} = 0, where \Omega_C represents the circular enclosing region.

I use the following code for my simulation:

/* Written for FreeFEM++ v4.15 */

// This solves the Helmholtz equation
// for a rectangular RF cavity with a small
// aperture

int top = 1;
int right = 2;
int bottom = 3;
int left = 4;
int edge = 5;
int all = 6;

real a = 0.1; // aperture diameter
real b = (1 - a) / 2;
real R = 10.0; // This should be a "large" number to simulate infinite domain

 // use (n, m)-th mode
int n = 1;
int m = 0;
real k = sqrt(n^2 + m^2) * pi;

border boxtop(t=1, 0){ x=t; y=1; label=top; } // top
border boxrightup(t=0, b){ x=1; y=t; label=right; } // upper right
border boxrightdown(t=0, b){ x=1; y=(1 - t); label=right; } // lower right
border boxbottom(t=0, 1){ x=t; y=0; label=bottom; } // bottom
border boxleft(t=1, 0){ x=0; y=t; label=left; } // left
border circle(t=0,2*pi){x=(R*cos(t) + 0.5); y=(R*sin(t) + 0.5); label=edge; } // circular boundary to fake a boundary at infinity

int nSides = 50;
int nOpening = 25;
int nCircle = 200;

mesh Th = buildmesh(
    boxtop(nSides) + boxrightup(nOpening)
    + boxrightdown(nOpening) + boxbottom(nSides) 
    + boxleft(nSides) + circle(nCircle)
    );
plot(Th, value=1, wait=1);

fespace Vh(Th, P1);
Vh<complex> TestPhi, Emag;

macro grad(u) [dx(u), dy(u)] //

solve Helmholtz(Emag, TestPhi)
    = int2d(Th)(
        grad(Emag)' * grad(TestPhi)
    )
    - int2d(Th) (
        k^2 * Emag * TestPhi
    )
    // radiative BC at "infinite" edge
    - int1d(Th, edge)(
        1i * k * Emag * TestPhi
    )
    // reflective BC
    + on(top, right, bottom, left, Emag=1);

Vh Ereal = real(Emag);
Vh Eimag = imag(Emag);
Vh Eabs = abs(Emag);

// plot real & imaginary parts of the field
plot(Ereal, value=true, fill=true, dim=2, cmm="Real Part of E-field (k=" + k + ")", wait=1);
plot(Eimag, value=true, fill=true, dim=2, cmm="Imaginary Part of E-field (k=" + k + ")", wait=1);
plot(Eabs, value=true, fill=true, dim=2, cmm="Magnitude of E-field (k=" + k + ")", wait=1);

The result for the n = 1 mode is shown below:

image
image1192Ă—719 96.5 KB

I, however, am somewhat confused about the fact the solution shows symmetric propagating waves when the aperture is towards the right side. I would expect to see an anistrophic pattern. Is this solution correct? If not, why is it different from the physical solution?

2

Hello,
I think that a problem is that your square cavity does not separate correctly the inside and outside parts. In order to do that you have to give some width to these walls, and create an “inside” boundary and and “outside” boundary, with no mesh inside the walls.

1 Like
3

Thank you for the tip! After trying your suggestion, I have obtained a much more reasonable-looking result:

image
image1197Ă—717 51.9 KB

This is my updated code for reference:

/* Written for FreeFEM++ v4.15 */

// This solves the Helmholtz equation
// for a rectangular RF cavity with a small
// aperture

int cavity = 1;
int aperture = 2;
int edge = 3;

real a = 0.3; // aperture diameter
real b = (1 - a) / 2;
real R = 3; // This should be a "large" number to simulate infinite domain
// real R = 10.0;

 // use (n, m)-th mode
int n = 1;
int m = 0;
real k = sqrt(n^2 + m^2) * pi;
real h = 0.05; // thickness of the RF cavity walls

// top wall
border boxtopInner(t=1, 0){ x=t; y=1; label=cavity; }
border boxtopOuter(t=1, 0){ x=((1 + 2*h)*t - h); y=1+h; label=cavity; }

// right wall
border boxrightupInner(t=0, 1){ x=1; y=b*t; label=cavity; } // upper right
border boxrightupOuter(t=0, 1){ x=1+h; y=((b + h)*t - h); label=cavity; }
border boxUpperConnector(t=1, 1+h) { x = t; y=b; label=cavity; }

border boxrightdownInner(t=0, 1){ x=1; y=(1 - b*t); label=cavity; } // lower right
border boxrightdownOuter(t=1, 0){ x=1+h; y=(1-b + (h + b)*t); label=cavity; }
border boxLowerConnector(t=1, 1+h) { x = t; y=(1-b); label=cavity; }

border boxaperture(t=0, 1){ x=1; y=((t + b/a)*a); label=aperture; } // aperture

// bottom wall
border boxbottomInner(t=0, 1){ x=t; y=0; label=cavity; } // bottom
border boxbottomOuter(t=0, 1){ x=((1 + 2*h)*t - h); y=-h; label=cavity; } // bottom

// left wall
border boxleftInner(t=1, 0){ x=0; y=t; label=cavity; } // left
border boxleftOuter(t=1, 0){ x=-h; y=((1 + 2*h)*t - h); label=cavity; } // left

// circular boundary to fake a boundary at infinity
border circle(t=0,2*pi){x=(R*cos(t) + 0.5); y=(R*sin(t) + 0.5); label=edge; }

int nSides = 50;
int nOpening = 25;
int nCircle = 200;

mesh Th = buildmesh(
    boxtopInner(nSides) 
    + boxtopOuter(nSides)
    + boxrightupInner(nOpening)
    + boxrightupOuter(nOpening)
    + boxUpperConnector(3)
    + boxrightdownInner(nOpening)
    + boxrightdownOuter(nOpening)
    + boxLowerConnector(3)
    + boxbottomInner(nSides) 
    + boxbottomOuter(nSides) 
    + boxleftInner(nSides)
    + boxleftOuter(nSides)
    + circle(nCircle)
    + boxaperture(nOpening)
    );
plot(Th, cmm="Discretized domain", value=1, wait=1);

fespace Vh(Th, P1);
Vh<complex> TestPhi, Emag;

macro grad(u) [dx(u), dy(u)] //

solve Helmholtz(Emag, TestPhi)
    = int2d(Th)(
        grad(Emag)' * grad(TestPhi)
    )
    - int2d(Th) (
        k^2 * Emag * TestPhi
    )
    // radiative BC at "infinite" edge
    - int1d(Th, edge)(
        1i * k * Emag * TestPhi
    )
    // perfectly-reflective BC
    + on(cavity, Emag=0)
    // absorption (neumann) BC
    + on(aperture, Emag=1);

Vh Ereal = real(Emag);
Vh Eimag = imag(Emag);
Vh Eabs = abs(Emag);

// plot real & imaginary parts of the field
plot(Ereal, value=true, fill=true, dim=2, cmm="Real Part of E-field (k=" + k + ")", ps="rf-cavity.ps", wait=1);
plot(Eimag, value=true, fill=true, dim=2, cmm="Imaginary Part of E-field (k=" + k + ")", wait=1);

I am, however, still somewhat concerned, due to the fact that I have to manually set the Dirichlet boundary condition E = E_0 = 1 on the aperture with + on(aperture, Emag=1). As far as I understand, across this boundary, I should have \frac{\partial E}{\partial n} = 0 which leads to the boundary term \displaystyle \int_{\partial \Omega}\phi \nabla E in the variational form to become zero, so there should be no need to hardcode this (\phi is the trial function). Is there some reason why this is necessary?

4

I don’t know much on the Helmholtz equation, thus I cannot say what is the physically relevant boundary condition.
What you do is to put a source at the gate, by specifying the value of E. What I would have done is, instead, is to put a source at the center of the square, via a right-hand side in the equation, localized on a small area. And I would have removed the boundary “boxaperture”. Note that in your code, this boundary is not a real boundary since the mesh goes through it.
But as I have said, I don’t know what is the physically relevant way to describe the source.