Transparent boundary condition for the calculation of eigenmodes in transversely infinite waveguides

Waveguides play one of the key figures in today's electronics and optics for signal transmission. Corresponding simulations of electromagnetic wave transportation along these waveguides are accomplished by discretization methods such as the Finite Integration Technique (FIT) or the Finite Element Method (FEM). For longitudinally homogeneous and transversely unbounded waveguides these simulations can be approximated by closed boundaries. However, this distorts the original physical model and unnecessarily increases the size of the computational domain size. In this article we present a boundary condition for transversely open waveguides based on the Kirchhoff integral which has been implemented within the framework of FIT. The presented solution is compared with selected conventional methods in terms of computational effort and memory consumption.

1 Introduction

In order to quantify electromagnetic wave propagation in longitudinally homogeneous waveguides, simulations are performed on a two-dimensional domain. Transversely unbounded waveguides such as the microstrip-line displayed in Fig. 1 or an on-chip open waveguide are often physically not enclosed by shields. For numerical simulations a specifically tailored boundary condition has to be considered.

Figure 1 Two-dimensional cross section in a reference plane used to calculate the modal fields on a transversely open microstrip line.

The two-dimensional problem of transversely unbounded waveguide can be calculated in terms of an eigenvalue problem. The required calculation is performed on a computational domain that consists in the considered case of a face and corresponding boundary conditions, by which the solution is defined. In contrary to closed waveguides, numerically efficient and physically correct representations of transversely open boundary conditions are still subject to research.

2 Conventional solutions

Simulations of the wave propagation along a three-dimensional, longitudinally homogeneous waveguides can be reduced to a two-dimensional formulation using the following ansatz:

$$\begin{array}{}\text{(1)}& \underset{\mathrm{¯}}{\mathit{E}}(x,y,z)=\underset{\mathrm{¯}}{\mathit{E}}(x,y)e^{-j\underset{\mathrm{¯}}{k_z}z},\end{array}$$

with $\underset{\mathrm{¯}}{\mathit{E}}$ being the phasor of the electric field strength and $\underset{\mathrm{¯}}{k_z}$ being the propagation constant of a corresponding mode along the propagation direction z. On account of the underlying frequency-domain formulation, a time harmonic dependency of electromagnetic field is assumed. Inserting Eq. (1) in Maxwell's equations, an eigenvalue problem can be formulated, where the eigenvectors represent the phasors of the electric field strength in the transverse plane (also known as waveguide modes) and the eigenvalues describe the propagation constants $\underset{\mathrm{¯}}{k_z}$ for each individual mode. Each solution for these modes in transversely unbounded and longitudinally homogeneous waveguides can be approximated using closed boundary conditions.

2.1 Closed boundary conditions

Hollow waveguides with good electric conductive materials can be simplistically simulated using Perfect Electric Conductor (PEC) instead of the conductive material. In this case, an efficient implementation can use this boundary condition to represent the conducting material.

A Perfect Magnetic Conductor (PMC) can also be used to formulate a closed boundary. Dependent on the specific application, PMC or PEC also handles a suitable symmetry condition. With PMC or PEC boundary conditions, the solution contains propagating, evanescent or even complex modes, which appear as complex conjugate pairs.

If a transversely unbounded waveguide is artificially truncated with a PEC or PMC boundary condition, the resulting solution will be distorted, since either the tangential or the normal electric or magnetic field components will be missing. Furthermore, in transversely inhomogeneous waveguides, complex modes have to be considered (Strube and Arndt, 1985; Monsoriu et al., 2003; Clarricoats and Slinn, 1965). However, in case of a longitudinally homogeneous open waveguide complex eigenmodes are only present at very specific conditions (Jabłoński, 1994). Furthermore, it is also important that choosing between PEC or PMC requires an a priori knowledge of the field distribution of a specific mode. Thus, for TEM or Quasi-TEM modes the PMC might generally be suited better than the PEC boundary condition, since the PEC boundary condition itself introduces an additional conductor. In order to obtain an accurate solution, closed boundary conditions must be set further away from the original waveguide which naturally is accompanied by an increase of the computational cost due to the enlarged problem size.

2.2 Absorbing boundary conditions

Transversely open waveguide problems can also be addressed using absorbing boundary conditions like perfectly matched layers (PML) (Bérenger, 1994). There are different possibilities to implement PML (Roden and Gedney, 2000; Johnson, 2010) but independent of the formulation the underlying idea is not changed. The PML implements a boundary condition that contains one or multiple layers of nonphysical materials, that absorb the energy of incidenting electromagnetic wave. If PML is used as a transverse boundary condition for the case of longitudinally homogeneous waveguides, the propagation direction of the wave is directed parallel to the boundary. Hence, the assumed conditions for a PML boundary to function properly are not met. Additionally, a PML introduces nonphysical modes (Talukder, 2009; Rogier and De Zutter, 2001) that unfortunately distorts the physical solution spectrum. However, these modes can be identified (Bandlow et al., 2010). Furthermore, if the PML is enclosed by PEC or PMC, similar modes as described in Sect. 2.1 can arise.

Figure 2 The energy of an incoming wave traveling towards the boundary at the angle α is partially absorbed within the PML region. The PML layer is terminated by a closed boundary condition such as PEC or PMC or by a simple open condition.

3 Enhanced solution approach

The core idea for the realization of a physically correct transversely open boundary condition within the framework of FIT resides in (i) the evaluation of the field in the exterior domain and in (ii) its implication on the location of the boundary of the computational domain. In the presented approach, the electromagnetic field is evaluated using the Kirchhoff Integral (Jackson, 1998). Contrary to specifically available closed or absorbing boundary conditions, this approach defines a generally available open boundary condition. After a short survey of the fundamentals of Kirchhoff Integral as well as the basic principles of FIT, the combination of both will be used to formulate an open boundary condition.

3.1 Fundamentals of the Kirchhoff integral

In Fig. 3, a region R₁, containing all the sources that are enclosed by a surface S₁ and another surface S₂ is considered. Hence, the region R₂ is a source-free and homogeneous domain. Electromagnetic fields in region R₂ generated by the sources in region R₁ can be efficiently described by corresponding field components on the surface S₁ (Jackson, 1998).

Figure 3 Sources in the region R₁ enclosed by the surface S₁. An additional surface S₂ is considered here, but will be moved to infinity afterwards for technical reasons.

Using the electric and magnetic field values on the surfaces S₁ and S₂, the electric or magnetic field strength in the region R₂ can be evaluated. By expressing the magnetic field strength $\underset{\mathrm{¯}}{\mathit{H}}$ through the electric field strength $\underset{\mathrm{¯}}{\mathit{E}}$, the calculation of the required fields inside the region R₂ can be conducted using the surface integral

$$\begin{array}{}\text{(2)}& \underset{\mathrm{¯}}{\mathit{E}}\left(\mathit{r}\right)=\underset{S_{\mathrm{1}}\cup S_{\mathrm{2}}}{∯}\left[\underset{\mathrm{¯}}{\mathit{E}}\right(\mathit{r}\left)\right({\mathit{n}}^{\prime }\cdot {\mathbf{\nabla }}^{\prime }G\left(\mathit{r}\right))-G(\mathit{r}\left)\right({\mathit{n}}^{\prime }\cdot {\mathbf{\nabla }}^{\prime }\left)\underset{\mathrm{¯}}{\mathit{E}}\right(\mathit{r}\left)\right]\mathrm{d}{a}^{\prime },\end{array}$$

with n^′ being the normal vector on the respective surfaces directed towards the region R₂ and G(r) being the three-dimensional Green's Function for the Helmholtz equation (Jackson, 1998). By placing the surface S₂ infinitely far away from S₁, the surface integral (Eq. 2) has to be evaluated only along S₁, as specified in Jackson (1998).

The required integration over the entire two-dimensional surface can be carried out partially on account of the a priori knowledge given in Eq. (1). As a result, a one-dimensional integral over a closed contour in the transverse plane of the investigated waveguide remains and can be summarized as

$$\begin{array}{}\text{(3)}& {\displaystyle }\begin{array}{rl}\underset{\mathrm{¯}}{\mathit{E}}\left(\mathit{r}\right)=& \underset{C}{\oint }\left[\underset{\mathrm{¯}}{\mathit{E}}\left(\mathit{r}\right)\left({\mathit{n}}^{\prime }\cdot {\mathbf{\nabla }}^{\prime }G\left(\underset{\mathrm{¯}}{k_z},\mathit{r}\right)\right)-G\left(\underset{\mathrm{¯}}{k_z},\mathit{r}\right)\right.\\ & \left.({\mathit{n}}^{\prime }\cdot {\mathbf{\nabla }}^{\prime })\underset{\mathrm{¯}}{\mathit{E}}\left(\mathit{r}\right)\right]\mathrm{d}{s}^{\prime },\end{array}\text{(4)}& {\displaystyle }G\left(\underset{\mathrm{¯}}{k_z},\mathit{r}\right)=\int G\left(\mathit{r}\right)e^{-j\underset{\mathrm{¯}}{k_z}z}\mathrm{d}z,\end{array}$$

where $G(\underset{\mathrm{¯}}{k_z}$, r) is the two-dimensional Green's function that additionally depends on the propagation constant due to the integration along the longitudinal coordinate. In the presented method, the open boundary problem is transferred into a nonlinear eigenvalue problem with $\underset{\mathrm{¯}}{k_z}$ being the desired eigenvalue. Hence, the linear eigenvalue problem in case of PMC or PEC boundary condition becomes nonlinear for the specified open boundary condition due to the two-dimensional Green's function being dependent on the eigenvalue.

3.2 Finite integration technique

The volume discretization method FIT (Weiland, 1977), Weiland (2017) is applicable for calculations in the frequency as well as in the time domain. It allocates electric and magnetic voltages and fluxes on a primal and a dual mesh. These field quantities are defined in FIT as integrals as follows¹

$$\begin{array}{ll}{\displaystyle }& \textcolor{red}{}={\int }_{C_i}\underset{\mathrm{¯}}{\mathit{E}}\cdot \mathrm{d}\mathit{s},\textcolor{red}{}={\int }_{{\stackrel{\mathrm{̃}}{C}}_i}\underset{\mathrm{¯}}{\mathit{H}}\cdot \mathrm{d}\mathit{s},\\ {\displaystyle }& \textcolor{red}{}={\iint }_{{\stackrel{\mathrm{̃}}{S}}_i}\underset{\mathrm{¯}}{\mathit{D}}\cdot \mathrm{d}\mathit{A},\textcolor{red}{}={\iint }_{S_i}\underset{\mathrm{¯}}{\mathit{B}}\cdot \mathrm{d}\mathit{A},\textcolor{red}{}={\iint }_{{\stackrel{\mathrm{̃}}{S}}_i}\underset{\mathrm{¯}}{\mathit{J}}\cdot \mathrm{d}\mathit{A}.\end{array}$$

Voltages and fluxes are allocated on mesh edges and faces, respectively, and are stored as algebraic vectors of dimension 3N_p.

$$\textcolor{red}{}=\left(\textcolor{red}{},\textcolor{red}{},\mathrm{\dots },\textcolor{red}{}\right),\textcolor{red}{}=\left(\textcolor{red}{},\textcolor{red}{},\mathrm{\dots },\textcolor{red}{}\right),$$

with N_p being the number of mesh nodes. Material settings are expressed as material matrices, thus the corresponding material relations can be expressed as follows:

$$\textcolor{red}{}={\mathbf{M}}_{\mathit{\mu }}\textcolor{red}{},\textcolor{red}{}={\mathbf{M}}_{\mathit{\epsilon }}\textcolor{red}{},\textcolor{red}{}={\mathbf{M}}_{\mathit{\kappa }}\textcolor{red}{}.$$

The required curl, divergence and gradient operators are defined on the primal and dual mesh as band matrices C, $\stackrel{\mathrm{̃}}{\mathbf{C}}$, S, $\stackrel{\mathrm{̃}}{\mathbf{S}}$, G, $\stackrel{\mathrm{̃}}{\mathbf{G}}$, respectively, with following properties

$$\begin{array}{}\text{(5)}& \begin{array}{rl}& \mathbf{SC}=\mathrm{0},\mathbf{CG}=\mathrm{0},\\ & \stackrel{\mathrm{̃}}{\mathbf{S}}\stackrel{\mathrm{̃}}{\mathbf{C}}=\mathrm{0},\stackrel{\mathrm{̃}}{\mathbf{C}}\stackrel{\mathrm{̃}}{\mathbf{G}}=\mathrm{0},\end{array}\end{array}$$

where the operators marked with a tilde are defined on the dual mesh. Assuming harmonic time dependence as before and substituting magnetic quantities by electrical quantities, a three-dimensional wave equation in the frequency domain can be formulated in form of the standard eigenvalue problem

$$\begin{array}{}\text{(6)}& {\mathbf{M}}_{\mathit{\epsilon }}^{-\mathrm{1}}\stackrel{\mathrm{̃}}{\mathbf{C}}{\mathbf{M}}_{{\mathit{\mu }}^{-\mathrm{1}}}\mathbf{C}\textcolor{red}{}={\mathit{\omega }}^{\mathrm{2}}\textcolor{red}{},\end{array}$$

where ω represents the angular frequency. The specified system matrix ${\mathbf{M}}_{\mathit{\epsilon }}^{-\mathrm{1}}\stackrel{\mathrm{̃}}{\mathbf{C}}{\mathbf{M}}_{{\mathit{\mu }}^{-\mathrm{1}}}\mathbf{C}$ is sparse and unsymmetrical, but can be symmetrized. The three-dimensional problem (Eq. 6) can be reduced to two dimensions and formulated in a way such that $\underset{\mathrm{¯}}{k_z}$ becomes the eigenvalue (Weiland, 1977; Schmitt, 1994). As a result, the two-dimensional eigenvalue problem

$$\begin{array}{}\text{(7)}& \mathbf{A}\left(\mathit{\omega }\right)\textcolor{red}{}={\underset{\mathrm{¯}}{k_z}}^{\mathrm{2}}\textcolor{red}{}\end{array}$$

can be derived. The system matrix A is sparse and unsymmetrical and can generally not be symmetrized.

3.3 Kirchhoff integral used as an open boundary condition for FIT

Consider the two-dimensional rectangular mesh specified in Fig. 4. The presented Kirchhoff Integral boundary (KIR) evaluates the electric field $\underset{\mathrm{¯}}{\mathit{E}}$ at the boundary B_c of the computational domain via Eq. (3). For this evaluation, the electric field values $\underset{\mathrm{¯}}{\mathit{E}}$ on the reference boundary B_r are used. In order to study the influence of the singularity in the Green's function on the modeling accuracy, a variable spatial distance between B_c and B_r is introduced. Since Eq. (3) depends on $\underset{\mathrm{¯}}{k_z}$, the two-dimensional eigenvalue problem with KIR as a boundary condition becomes nonlinear. FIT 2D eigenvalue matrix (Eq. 7) becomes denser and non symmetrizable if KIR is used as a boundary condition.

Figure 4 (a) Model of a longitudinally infinite and transversely open artificial optical waveguide with η₁, η₂ being different indices of refraction and R being the radius of the cylindrical waveguide. (b) Two-dimensional rectangular mesh with the reference surface B_r and computational domain boundary B_c. Parameters are as follows: Δx=0.4 mm, frequency = 16 GHz, radius R=6.5 mm, L=4 mm, η₁=2, η₂=1.5.

In order to integrate KIR into the framework of FIT multiple steps are necessary. First, the boundary field values of the computational domain must be replaced by the integral (Eq. 3) along the cell edges. After evaluating the integral, the relationship for a boundary field cell  at the boundary point r is obtained:

$$\begin{array}{ll}{\displaystyle }& \textcolor{red}{}\\ \text{(8)}& {\displaystyle }& \textcolor{red}{}\end{array}$$

with c_0…U being the KI weights and _0…U,q being the field values at the reference boundary.

In the next step, Eq. (8) will be rearranged so that the values of individual reference sources _i,q can be expressed as a function of the boundary value _x,r and other reference sources.

$$\begin{array}{}\text{(9)}& \textcolor{red}{}=f\left(\textcolor{red}{},\textcolor{red}{},c_{\mathrm{0}\mathrm{\dots }U}\right)\end{array}$$

Each row and column of the A matrix from Eq. (7) that corresponds to a specific reference source _i,q must be altered according to Eq. (9), thus a separate matrix A_KIR is built up, which consists of the KIR components. This matrix is not sparse, non symmetric, complex-valued and has the same dimensions as the two-dimensional FIT system matrix A from Eq. (7). The structure of the FIT and KIR matrix are shown in Fig. 5.

Figure 5 (a) Structure of the KIR matrix A_KIR. (b) Matrix structure of the two-dimensional FIT eigenvalue problem (Eq. 7).

After A_KIR has been set up, it can be integrated into the eigenvalue problem. For this, the A_KIR matrix has to be added to the system matrix A:

$${\mathbf{A}}_K=\mathbf{A}+{\mathbf{A}}_{\mathrm{KIR}},$$

with A_K being the new system matrix of the nonlinear eigenvalue problem.

4 Computational examples

In order to evaluate the numerical errors of KIR, a representative computational example will be described in this section. In the following examinations, a transversely open cylindrical optical waveguide according to Fig. 4 will be considered. The presented waveguide supports 4 propagating modes which can be calculated analytically (Yeh and Shimabukuro, 2008). Using FIT with KIR, the numerically obtained eigenvalue solution for the monopole TE mode is compared with the analytical solution. This specific monopole mode is chosen, since the azimuthal angle of its field distribution in the transverse plane doesn't change, hence the comparison between the analytical and numerical solutions become easier. The analytical solution of the chosen mode's propagation constant is given by $\underset{\mathrm{¯}}{k_z}=\mathrm{520.728}$ m⁻¹. In order to quantify the field error, the following L₂ norm is used:

$$\begin{array}{}\text{(10)}& \left|\right|\textcolor{red}{}-\textcolor{red}{}|{|}_{\mathrm{2}}=\sqrt{{\displaystyle \frac{\mathrm{1}}{A}}{\sum }_i|\textcolor{red}{}-\textcolor{red}{}{|}^{\mathrm{2}}{a}_i},\end{array}$$

with $A={\sum }_i{a}_i$, , being the analytical and numerical solution, respectively, A being the computational area inside B_r and a_i being the area of individual cells i. In this particular case of the equidistant mesh in both coordinate directions all a_i areas in the computational domain are equal.

5 Investigation of KIR properties

The boundary formulation presented here has two parameters that can be used to adjust the accuracy and computation cost of the calculation. Both parameters will be introduced in the following subsections. Additionally, a particular KIR for spatially symmetrical structures is discussed.

5.1 Solution of the nonlinear eigenvalue problem

The solution of linearized subsystems of the nonlinear eigenvalue problem can be obtained using a fixed point iteration (FPI), which is used here to calculate individual eigenvalue pairs. The first eigenvalue $\underset{\mathrm{¯}}{k_z}$ in Eq. (3) is estimated by the averaged $\underset{\mathrm{¯}}{k_{z,\mathrm{avg}}}$ within the iterative solution process, which leads to a sequence of linear eigenvalue problems. Knowing the material properties of the computational domain, known limits for the smallest and highest propagation constants can be specified as

$$\begin{array}{}\text{(11)}& \underset{\mathrm{¯}}{k_{z,\min }}=\min \left(\sqrt{{\mathit{\epsilon }}_i{\mathit{\mu }}_i}\right)k_{\mathrm{0}},\underset{\mathrm{¯}}{k_{z,\max }}=\max \left(\sqrt{{\mathit{\epsilon }}_i{\mathit{\mu }}_i}\right)k_{\mathrm{0}},\end{array}$$

where ε_i and μ_i are the relative permittivity and permeability of an individual cell, respectively. Considering a real-valued $\underset{\mathrm{¯}}{k_z}$ parameter,
$\underset{\mathrm{¯}}{k_{z,\mathrm{avg}}}$ is set to the mean value between $\underset{\mathrm{¯}}{k_{z,\min }}$ and $\underset{\mathrm{¯}}{k_{z,\max }}$ in the initial iteration of the FPI. On the subsequent iterations, $\underset{\mathrm{¯}}{k_{z,\mathrm{avg}}}$ is set to one of the obtained solutions. The amount of required iterations depends on the required relative accuracy between individual iterations.

The numerical errors of the eigenvector depend on the difference $|\underset{\mathrm{¯}}{k_{z,\mathrm{avg}}}-\underset{\mathrm{¯}}{k_z}|$. The convergence of the FIT model with KIR is shown in Fig. 6, with the L₂ error on the ordinate and the mesh cell size on the abscissa. Calculations with FIT result in a convergence order between 1 and 2 if the error is calculated using field integrals. If instead one uses field strength in order to calculate the error, a convergence order between 0.5 and 1 is obtained. In the numerical example specified in Fig. 4 a cylindrical figure is being discretized by a quadratic mesh. Additionally electric field strength is being used to calculate the error. These two reasons lead to the convergence order of 0.51, shown in Fig. 6. Two additional calculations were performed at different fixed $\underset{\mathrm{¯}}{k_{z,\mathrm{avg}}}$ and one was conducted using the FPI. In this particular example it is evident, that (i) the error stagnates if no FPI is used and (ii) the error offset is decreasing when $\underset{\mathrm{¯}}{k_{z,\mathrm{avg}}}$ is getting closer to the analytical value.

Figure 6 Convergence study of FIT with KIR for the numerical example specified in Fig. 4, with and without FPI. The solution obtained with FIT and KIR converges wit the order of 0.51.

5.2 Distance to the boundary

One of the important parameters is defined as the distance L between the reference contour B_r and the boundary B_c as specified in Fig. 4. Changing this value affects the ratio between R_min, R_max, which are the closest and furthest distances between the evaluated field on B_c and the reference elements on B_s.

As shown in Fig. 7, the evaluated numerical error increases in two cases: (i) if the distance to the boundary becomes too small, thus the evaluated Green's function gets closer to the singularity and (ii) if the ratio R_min∕R_max increases. The choice of an optimal value for L depends on the application.

Figure 7 Numerical error over the increasing ratio R_min∕R_max for the computational example shown in Fig. 4. The distribution of the L₂-Error is shown for two calculations in (a) and (b).

Apart from the numerical error, the distance parameter also affects the computational cost. In Fig. 7 the determined amount of floating point operations (FLOP) dependent on the ratio R_min∕R_max is shown. Although the adjustment of the distance parameter has a relatively low influence on the overall error, the distribution of the corresponding L₂ error shows significant differences.

5.3 KIR for symmetrical structures

Symmetrical waveguides may support modes with a symmetric field distribution with regard to one of the axes of the chosen coordinates. In that case, a dedicated boundary condition can be imposed for the numerical simulation. This reduces the size of the computational domain in the given two-dimensional case by a factor of two of even four depending on the amount of symmetry planes. In case of the example given in Fig. 8 the computational domain can be divided in two parts and only one of those must be taken into account for the calculation. For the selected domain, a dedicated boundary condition must be specified on the dashed symmetry line.

Figure 8 Example of a horizontally symmetrical structure.

The proposed boundary condition KIR can naturally be used in combination with symmetry-enforcing boundary conditions. However, some implementation details must be particularly considered. Since in case of Fig. 8 one half of the reference sources on B_r is missing, these values must be considered using the a priori knowledge of the symmetry. A simplified example in Fig. 9 is used to explain the application of the method. Without the symmetry boundary, for example the value of e_A is calculated using:

$$\begin{array}{}\text{(12)}& e_A={\sum }_{n=\mathrm{1}}^{\mathrm{12}}{e}_{n}f\left(r_n\right),\end{array}$$

with r_n as the distance between the edge n and the edge A and f(r_n) as the function calculating the weights of the Kirchhoff-Integral. Field values e₇–e₁₂ can be calculated using the knowledge of the values e₁–e₆ and the waveguide symmetry.

Figure 9 Schematic of a model with a possible symmetry plane and its implication on the field evaluation formulation.

Using a symmetry boundary condition leads to a reduction of the system matrix by the factor of two.

Figure 10 Visualization of the obtained error in the L₂ norm as a function of the considered reference sources in percent represented as parameter P as well as the corresponding amount of FLOP for the computational example shown in Fig. 4.

5.4 Reduction of reference sources

The evaluation of the Kirchhoff Integral requires the integration over a closed surface. If the reference point is close to the respective surface, only a small fraction of the entire surface may contribute to the value on account of the weighting with the Green's functions kernel. If this situation holds, considering a fraction of the entire surface is sufficient to obtain an adequate solution.

Another KIR parameter, further referred to as parameter P defines the amount of reference sources in percentage used for evaluating an individual field value. On the one hand, the overall computational cost decreases with the number of reference sources. On the other hand, lower values of P leads to a higher numerical error (see Fig. 10). Similar investigations in form of a convergence study also shows the dependency on the fraction of the considered reference sources as can be seen in Fig. 11. For this particular example it is recognizable, that considered reference sources between 25 % and 100 % show a similar convergence behavior. The intentionally introduced reduction of reference sources leads to a smaller amount of DoF in the system matrix specified in Eq. (7), which in turn decreases the overall computational cost as expected.

Figure 11 Convergence study of FIT with KIR as a function of the mesh cell size and the fraction of the reference sources considered. Calculations are performed for the representative example shown in Fig. 4, with $\mathrm{\Delta }x<R/\mathrm{10}$.

6 Sources of KIR approximations

The presented KIR boundary representation contains three types of approximations that will be addressed in this sections:

1. numerical derivative of the electric field strength;

2. numerical integration of the electric field strength;

3. numerical integration of the Green's function.

6.1 Numerical derivative of the electric field strength

The evaluation of the Kirchhoff Integral (Eq. 3) requires the calculation of a derivative of the electric field strength. In order to evaluate this quantity numerically within the framework of FIT, a central differential quotient is used:

$${\displaystyle \frac{\mathrm{\Delta }f}{\mathrm{\Delta }x}}:={\displaystyle \frac{f(x+\mathrm{\Delta }x)-f(x-\mathrm{\Delta }x)}{\mathrm{2}\mathrm{\Delta }x}}.$$

The order of this approximation is 𝒪(Δx²), which is the same for FIT in the used formulation with field integrals as eigenvector components. Since this approximation depends on the mesh size, it is not possible to analyze its sole influence on the overall convergence rate independently because changing the mesh size will also affect the approximations of FIT.

6.2 Numerical integration of the electric field strength

The Kirchhoff Integral (Eq. 3) allows the evaluation of the electric field strength at any position in a given region. However, for the used FIT formulation, integrals of electric fields are required, as described in Sect. 3.2. Thus, a numerical integration of the result is needed, which is simply achieved by the midpoint rule:

$${\int }_a^{b}f\left(x\right)\mathrm{d}x:=(b-a)\cdot f\left({\displaystyle \frac{a+b}{\mathrm{2}}}\right).$$

The approximation error here is also 𝒪(Δx²), which is again the same for FIT in the used formulation with field integrals as eigenvector components. Just like the approximation of the derivative, here it is also not possible to analyze the sole influence of this error on the overall convergence rate.

6.3 Numerical integration of the Green's function

In the one-dimensional Kirchhoff Integral (Eq. 3) both, the electric field strength and the Green's function are integrated over a contour C. The integration over an individual cell edge with an assumed constant field quantity allows to pull out the electric field from the actual integral. Hence, the remaining Green's function and its derivative can be integrated separately. For example by using the Gaussian quadrature (Stoer, 2005). The error for a Gaussian quadrature of the order 2n+2 is defined as

$$Q^{\left(n\right)}\left(f\right)={\displaystyle \frac{f^{(\mathrm{2}n+\mathrm{2})}\left(\mathit{\theta }\right)}{(\mathrm{2}n+\mathrm{2})}}{\int }_a^b{p}_{n+\mathrm{1}}^{\mathrm{2}}\left(x\right)\mathrm{d}x$$

with

$$f\in C^{\mathrm{2}n+\mathrm{2}},\mathit{\theta }\in [a,b]$$

and

$$p_{n+\mathrm{1}}={\prod }_{j=\mathrm{0}}^n\left(x-{\mathit{\lambda }}_j\right).$$

This approximation can be evaluated independent of the mesh cell size. Using the example described in Sect. 4, the L₂ error was calculated with Gaussian quadrature orders n=2, 4, 8. In Table 1, selected results are displayed as a relative difference to the midpoint rule integral. The accuracy gain between the midpoint rule and higher order Gaussian quadrature is relatively low, as this example shows.

Table 1 L₂ error for different orders of Gaussian quadrature for the numerical integration of the Green's function. The relative gain is defined as $|f_{\mathrm{g}}-f_{\mathrm{m}}|/\left|f_{\mathrm{g}}\right|$, with f_g and f_m being the solution using Gaussian quadrature and midpoint rule respectively.

7 Applications of KIR

In this section, the presented KIR boundary is compared with closed and absorbing boundaries PEC, PMC and PML. All required algorithms were implemented within the framework of FIT in a custom MATLAB^® code. The PML was realized as a uniaxial PML (UPML) as described in Tischler (2004). In order to compare different boundary conditions with each other, three different computational examples are examined.

First, the numerical solution and its error is obtained using KIR. Next, the distance to the boundary is varied such that the same eigenvector accuracy, defined by the L₂ norm (Eq. 10) is obtained for PEC, PMC and PML termination. Compared quantities are given by the distance between both boundaries, the amount of degrees of freedom (DoF), the amount of FLOP, the amount of nonzero elements (NZE) in the eigenvalue matrix and the amount of mesh cells.

7.1 Transmission line

A transversely open double cylindrical wire transmission line as show in Fig.  is used as the first application. For this configuration an analytical solution for the TEM mode with $\underset{\mathrm{¯}}{k_z}=k_{\mathrm{0}}$ exists and is used as a reference solution. In Table 2 the simulation results obtained are summarized and it is apparent that the KIR boundary uses less DoF and less cells then the other examined types of boundaries. For this calculation with KIR only one FPI step was needed, since the only existing propagation mode in this example is know in advance. The solution obtained with PEC shows a higher distance between boundaries B_r and B_c than in the case with PMC. This is due to the fact, that the PEC boundary introduces an additional conductor, which distorts the physical model.

Figure 12 Layout of an open two-wire transmission line. Parameters are as follows: mesh cell size Δx=0.3 mm, frequency = 16 GHz, diameter d=3 mm, distance S=12 mm.

Table 2 Comparison between different boundary conditions used to terminate the transmission line model from Fig.  as a computational example. This shows that the computational effort needed for the simulation is lower for the calculation performed using KIR.

^a Distance between boundaries B_r and B_c. ^b Number of Non-Zero-Elements. ^c Amount of cells.

7.2 Optical waveguide

An uncoated optical waveguide as shown in in Fig. 4 is considered in the second application example. The model itself is described in Sect. 4.

Results for the comparison between different boundary conditions with the KIR boundary is shown in the Table 3. For this calculation with KIR two FPI steps were needed. Two iterations were needed to reach the solution with KIR. Although the amount of FLOP for calculations with PEC and KIR boundaries are equal, the amount of DoF differs by the factor of 2.5. What is more important is the fact, that the obtained solution spectrum with the PEC boundary is polluted with modes that do not correspond to the initial physical model.

Table 3 Comparison between different boundary conditions used to terminate the optical waveguide from Fig. 4 as a computational example. This shows that the computational effort needed for the simulation is lower for the calculation performed using KIR.

^a Distance between boundaries B_r and B_c. ^b Number of Non-Zero-Elements. ^c Amount of cells.

7.3 Microstrip line

A microstrip line, displayed in Fig. 13 is used as a representative example. Since no analytical solution for this problem can be found, a numerical reference solution was calculated using the commercial CEM Software CST MICROWAVE STUDIO^®. For the comparison the quasi TEM mode is chosen here.

Figure 13 Microstrip line with a dielectric substrate and PEC elements. Following parameters for this model are used: h₁=0.12 mm, h₂=0.012 mm, l₁=1.2 mm, l₂=0.12 mm, ε_r=4, mesh width Δx=0.012 mm, frequency = 16 GHz.

In Table 4 the compiled results show that the model with the KIR boundary requires less DoF and less mesh cells than the other models calculated using classical boundaries on the microstrip line described above. For this calculation with KIR one FPI step was needed. The solution obtained with the PMC model shows computational costs in the same order of magnitude as the one, calculated using KIR. Nevertheless, with the remaining boundaries the eigenvalue spectrum contains modes that do not exist in the transversely open waveguide (e.g. complex modes).

Table 4 Comparison between different boundary conditions used to terminate the microstrip line from Fig. 13 as a computational example. This shows that the computational effort needed for the simulation is lower for the calculation performed using KIR.

^a Distance between boundaries B_r and B_c. ^b Number of Non-Zero-Elements. ^c Amount of cells.

8 Conclusions

An alternative boundary condition KIR for transversally open waveguides was defined and successfully implemented within the framework of FIT. Various parameters and approximations of KIR were introduced and discussed. Finally, a comparison between conventional boundary conditions and KIR was performed. In contrast to classical boundary conditions, KIR representation results in a better approximation of the physical problem. This significantly lowers the computation time.