Waveguide Transmission Benchmark

Model Geometry

Figure 1: Two-dimensional vacuum-filled waveguide geometry, with an incoming plane wave at left and exiting at right.

Governing Equations and Boundary Conditions

where

• is the complex electric field vector in V/m,

• is the vacuum magnetic permeability of the medium in H/m,

• is the electric permittivity of the medium in F/m, and

• is the operating frequency in rad/s.

Harrington in (Harrington, 1961) showed that the transverse electric field distributions could be determined from the component in the direction of wave travel, so only the scalar component will need to be modeled for comparison. Further, because the real and imaginary components of the field are 90 out of phase with each other, verifying both the form of the real component and the phase shift should be adequate to check for proper performance.

Given Gauss's Law with an absence of charge density (as the waveguide is filled with vacuum), a standard diffusion-reaction style equation (without a source) is the final form of the equation being modeled, with the Diffusion and ADMatReaction objects being used for each real and imaginary component.

With respect to boundary conditions, is set to zero on the waveguide walls, to satisfy the perfect electrical conductor boundary condition for this case. For the entry port, the EMRobinBC object will be used in a "port" configuration, while the exit port will be set in an "absorbing" configuration. The incoming wave shape across the waveguide port needed by EMRobinBC will be set to

Analytic Solution for Comparison

where

• (or just for a 2-D case),

• is the propagation constant of the wave,

• is the free space wavenumber, with being the speed of light,

• is the operating frequency in rad/s,

• is the peak electric field amplitude,

• is the width of the waveguide in the direction, and

• is the width of the waveguide in the direction.

Because our two-dimensional case is effectively infinite into and out of the page, the factor can be removed to render the analytic function for this verification.

Mesh

gmsh -2 waveguide.geo -clscale 0.12 -order 1 -algo del2d

As of Gmsh 4.6.0, this command produces a first order 2D mesh with an element size factor of 0.12 using a Delaunay algorithm. The unstructured, triangular mesh contains 1202 nodes and 2806 elements. An image of the result is shown in Figure 2.

Figure 2: Mesh used in the waveguide transmission benchmark study.

// Use global scaling factor = 0.12 to duplicate current saved MSH file.

width = 10;
length = 80;

Point(1) = {0, 0, 0};
Point(2) = {length, 0, 0};
Point(3) = {length, width, 0};
Point(4) = {0, width, 0};

Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};
Line(4) = {4, 1};
Line Loop(5) = {1, 2, 3, 4};

Plane Surface(6) = {5};

Physical Line("bottom") = {1};
Physical Line("exit") = {2};
Physical Line("top") = {3};
Physical Line("port") = {4};

Physical Surface(7) = {6};

Input File

# Test for EMRobinBC in port and absorbing modes with simple electric plane wave
# 2D, vacuum-filled waveguide with conducting walls
# u^2 + k^2*u = 0, 0 < x < 80, 0 < y < 10, u: R -> C
# k = 2*pi*freq/c, freq = 20e6 Hz, c = 3e8 m/s

[Mesh]
  [fmg]
    type = FileMeshGenerator
    file = waveguide.msh
  []
[]

[Variables]
  [E_real]
    order = FIRST
    family = LAGRANGE
  []
  [E_imag]
    order = FIRST
    family = LAGRANGE
  []
[]

[Functions]
  [inc_y]
    type = ParsedFunction
    expression = 'sin(pi * y / 10)'
  []
[]

[Kernels]
  [diffusion_real]
    type = Diffusion
    variable = E_real
  []
  [coeffField_real]
    type = ADMatReaction
    reaction_rate = kSquared
    variable = E_real
  []
  [diffusion_imaginary]
    type = Diffusion
    variable = E_imag
  []
  [coeffField_imaginary]
    type = ADMatReaction
    reaction_rate = kSquared
    variable = E_imag
  []
[]

[BCs]
  [top_real]
    type = DirichletBC
    value = 0
    variable = E_real
    boundary = top
  []
  [bottom_real]
    type = DirichletBC
    value = 0
    variable = E_real
    boundary = bottom
  []
  [port_real]
    type = EMRobinBC
    coeff_real = -0.27706242940220277  # -sqrt(k^2 - (pi/10)^2)
    sign = positive
    profile_func_real = inc_y
    profile_func_imag = 0
    field_real = E_real
    field_imaginary = E_imag
    variable = E_real
    component = real
    mode = port
    boundary = port
  []
  [exit_real]
    type = EMRobinBC
    coeff_real = 0.27706242940220277
    sign = negative
    field_real = E_real
    field_imaginary = E_imag
    variable = E_real
    component = real
    mode = absorbing
    boundary = exit
  []
  [top_imaginary]
    type = DirichletBC
    value = 0
    variable = E_imag
    boundary = top
  []
  [bottom_imaginary]
    type = DirichletBC
    value = 0
    variable = E_imag
    boundary = bottom
  []
  [port_imaginary]
    type = EMRobinBC
    coeff_real = -0.27706242940220277
    sign = positive
    profile_func_real = inc_y
    profile_func_imag = 0
    field_real = E_real
    field_imaginary = E_imag
    variable = E_imag
    component = imaginary
    mode = port
    boundary = port
  []
  [exit_imaginary]
    type = EMRobinBC
    coeff_real = 0.27706242940220277
    sign = negative
    field_real = E_real
    field_imaginary = E_imag
    variable = E_imag
    component = imaginary
    mode = absorbing
    boundary = exit
  []
[]

[Materials]
  [kSquared]
    type = ADParsedMaterial
    property_name = kSquared
    expression = '0.4188790204786391^2'
  []
[]

[Executioner]
  type = Steady
  solve_type = 'PJFNK'
  petsc_options_iname = '-pc_type'
  petsc_options_value = 'lu'
[]

[Outputs]
  exodus = true
  print_linear_residuals = true
[]

Results and Discussion

Figure 3: Electric field result, (real component), of the 2-D waveguide verification case.

Figure 4: Electric field result, (imaginary component), of the 2-D waveguide verification case.

A comparison of the real component results to compared to that in Eq. (1) is shown in Figure 5. This data is taken down the centerline of the waveguide geometry (). There is decent agreement between the two solutions, but note the increasing phase error along the length of the waveguide. This is known as "numerical dispersion" or "numerical phase error" in the literature, and is a byproduct of the numerical discretization, where the velocity of the calculated wave differs from that of the exact velocity due to lower resolution. Jin as well as Warren and Scott showed that the use of higher order elements as well as smaller elements can help mitigate this effect. (Jin, 2014; Jin, 2010; Warren and Jr., 1996) Further, Warren and Scott later showed that the orientation of the elements in the mesh for a triangular mesh, as used here, had a notable impact on the calculated phase error. (Warren and Jr., 1995) Since the error accumulates in the direction of wave travel, alternating the direction of the triangular elements in a hexagonal configuration had the lowest phase error in a TM1 mode simulation. Interestingly, randomizing the triangular element direction did not have the same impact (though did achieve lower phase error overall). A recalculation of this result using a finer mesh is shown in Figure 6, and shows that the refined (18,102 elements compared to 2,219 in the original), simulated result is almost exactly in alignment with the analytic solution.

Figure 5: Comparison of the real electric field result with the analytic solution in the 2-D waveguide verification case.

Figure 6: Comparison of the real electric field result with the analytic solution in the 2-D waveguide verification case, with a refined mesh.

Finally to confirm the proper resolution of the phase between real and imaginary components (and to show that the port and absorbing boundary conditions have little reflection), a plot of the magnitude of the electric field is shown below in Figure 7. If the components are 90 out-of-phase with each other, then the magnitude of the field should be a profile, extending the length of the waveguide. Indeed, Figure 7 confirms that this is the case.

Figure 7: Electric field result for the magnitude of in the 2-D waveguide verification case, confirming that the real and imaginary components are in the proper phase relative to each other.

References