4.4 FEM solution for the two-dimensional Helmholtz equation

Finite element methods (FEM), initially developed for structural mechanics [70], are widely used for numerical simulations of electromagnetic fields [15], [23], [71], [72], [73]. Some advantages of the finite element method are listed below.

• A sparse system matrix enables the simulation of problems with a large
  number of unknowns.
• Suitable for inhomogeneous simulation domains.
• Adaptability to arbitrary geometries by triangular meshes for two-dimensional
  simulations and tetrahedral meshes for three-dimensional problems.
• Suitable for inhomogeneous simulation domains.
• Existence of an unambiguous solution for elliptical partial differential equations.
  This means convergence of the solution error towards zero with mesh refinement.

The finite element method provides an approximate solution $\tilde{U}$ for the elliptical partial differential equation (4.13). The residuum (approximation error) of this solution is

$$Res=\Delta \tilde{U}+\omega^{2}\mu\epsilon\left(1-\frac{j}{Q(\omega)}\right)\tilde{U}+j\omega\mu hJ_{z}.$$ (4.24)

For consistence to the voltage and current directions, defined by the impedance matrix in (4.17), a negative sign has been used for the current density $J_{z}$, as opposed to the sign in (4.13). The residuum $Res$, weighted with a function $W_{g}$, will vanish on average over the simulation domain. This is expressed in the variation integral

$$\\ \int_{\mathcal{A}}\left(\Delta \tilde{U}+\omega^{2}\mu\epsi... ...mega)}\right)\tilde{U}+j\omega\mu hJ_{z}\right)W_{g}\,\textrm{d}\mathcal{A}=0.$$ (4.25)

$\mathcal{A}$ denotes the whole simulation domain, that is, the parallel plane surface. Applying Green's theorem, (4.25) is transformed to the weak formulation

$$\int_{\mathcal{A}}\left(\vec{\nabla}\tilde{U}\cdot\vec{\nabla}W_{... ... =\int_{\mathcal{A}}\left(j\omega\mu hJ_{z}W_{g}\right)\,\textrm{d}\mathcal{A}.$$ (4.26)

$\mathcal{\partial A}$ is the boundary curve of the parallel plane surface $\mathcal{A}$, $\textrm{d}s$ is the line segment of this boundary curve and $\partial /\partial_{n}\tilde{U}$ is the normal derivation of $\tilde{U}$ at the boundary. $\tilde{U}$ is expressed by

$$\tilde{U}=\sum_{k=1}^{p}\alpha_{k}(x,y)\tilde{u}_{k}.$$ (4.27)

Where $\alpha_{k}$ denote the finite element base functions, $\tilde{u}_{k}$ are the values of the solution at the mesh points, and $p$ denotes the number of mesh nodes. According to the method of Galerkin, the weighting function is

$$W_{g}=\sum_{k=1}^{p}\alpha_{k}(x,y).$$ (4.28)

Linear base functions on a triangular mesh are selected for the finite element discretization. Barycentric coordinates are utilized according to Figure 4.5.

Figure 4.5: Barycentric (triangular) coordinates $\xi$ and $\eta$. Coordinates ($\xi$,$\eta$) become (1,0), (0,1), and (0,0) in the nodes 1, 2, and 3, respectively.

With this coordinate system definition the base functions are

$$\alpha_{1}=\xi,\,\alpha_{2}=\eta, \alpha_{3}=1-\xi-\eta,$$ (4.29)

with indexes according to the triangle node numbers in Figure 4.5. The relation of the barycentric to the cartesian coordinates of the triangular nodes is

$$\begin{pmatrix} \xi_{i}\\ \eta_{i}\\ \zeta_{i}\\ \end{pmat... ...2i}-x_{1i}\\ \end{pmatrix} \begin{pmatrix} 1\\ x\\ y\\ \end{pmatrix},$$ (4.30)

with

$$A_{i}=\frac{1}{2}\abs \begin{vmatrix} 1&x_{1i}&y_{1i}\\ 1&x_{2i}&y_{2i}\\ 1&x_{3i}&y_{3i}\\ \end{vmatrix} .$$ (4.31)

$A_{i}$ is the surface area of the triangle with index $i$. With equations (4.27) and (4.28), (4.26) becomes the linear equation system

$$\textbf{K}\tilde{\textbf{U}}-\textbf{R}\tilde{\textbf{U}}_{n}=\textbf{F},$$ (4.32)

for the voltages $\tilde{\textbf{U}}$ on the mesh nodes. $\tilde{\textbf{U}}_{n}$ denotes the vector of the normal derivatives of $\tilde{\textbf{U}}$ at the boundary. At open boundaries the PMC boundary condition requires vanishing $\tilde{\textbf{U}}_{n}$. On closed metallic edges a Dirichlet boundary has to be introduced, by setting the boundary node voltage to zero and reducing the order of the linear system accordingly. Since the weighting function is set to zero in the boundary nodes with Dirichlet condition, the term with $\tilde{\textbf{U}}_{n}$ vanishes also at Dirichlet boundaries and (4.32) is reduced to

$$\textbf{K}\tilde{\textbf{U}}=\textbf{F}.$$ (4.33)

The system matrix elements $K_{kl}$ are obtained from the solution of

$$K_{kl}=\int_{\mathcal{A}}\left[\vec{\nabla}\alpha_{k}\cdot\vec{\n... ...1-\frac{j}{Q(\omega)}\right)\alpha_{k}\alpha_{l} \right]\textrm{d}\mathcal{A}.$$ (4.34)

The term $\vec{\nabla}\alpha_{k}\cdot\vec{\nabla}\alpha_{l}$ in( 4.34) becomes

$$\vec{\nabla}\alpha_{k}\cdot\vec{\nabla}\alpha_{l}= \left(\frac{\... ...x}+\frac{\partial\alpha_{k}}{\partial y}\frac{\partial\alpha_{l}}{\partial y}.$$ (4.35)

With the base functions

$$\alpha_{k}=\xi \alpha_{l}=\eta$$ (4.36)

and (4.30), the partial derivatives of (4.35) become

$$\begin{split} \frac{\partial\alpha_{k}}{\partial x}&=\frac{\... ...ac{\partial\xi}{\partial y}=\frac{c_{li}}{2A_{i}}. \end{split}$$ (4.37)

Using these partial derivatives expressions (4.34) becomes

$$K_{kl}=\int_{\mathcal{A}}\left[\frac{b_{ki}b_{li}+c_{ki}c_{li}}{4... ...epsilon\left(1-\frac{j}{Q(\omega)}\right)\xi\eta \right]\textrm{d}\mathcal{A}.$$ (4.38)

Utilizing

$$\int_{\mathcal{A}}\xi^a\eta^b\zeta^c\textrm{d}\mathcal{A}_{i}=2\mathcal{A}_{i}\frac{a!b!c!}{(a+b+c+2)!},$$ (4.39)

from [23], the integral of (4.38) is solved analytically and the coefficients of the FEM matrix $K_{kl}$ finally become

$$K_{kl}=\sum_{i=1}^{p}\left[\frac{b_{ki}b_{li}+c_{ki}c_{li}}{4A_{i}} -\omega^{2}\mu\epsilon\left(1-\frac{j}{Q(\omega)}\right)K_{kli}^{(1)}\right]$$ (4.40)

with

$$K_{kli}^{(1)}=\begin{cases} A_{i}/6&\text{,if $k=l$}\\ A_{i}/12&\text{,if $k\neq l$}. \end{cases}$$ (4.41)

At the node position $(x_{l},y_{l})$ the port excitation current is

$$J_{z}(x,y)=I_{l}\delta(x_{l},y_{l}).$$ (4.42)

$\delta(x_{l},y_{l})$ is the Dirac pulse function and $l$ is the node index running from one to the number of nodes $p$. With (4.26), (4.28), (4.36), and (4.42), the coefficients of the excitation vector $\textbf{F}$ in (4.33) become

$$F_{l}=\int_{\mathcal{A}}\left(j\omega\mu hI_{l}\delta(x_{l},y_{l})\eta\right)\textrm{d}\mathcal{A} =j\omega\mu hI_{l},$$ (4.43)

because the weighting function $\eta$ is unity at node position $(x_{l},y_{l})$ according to the base function definition in (4.29). With this port current excitation definition the impedance matrix, which relates the voltages on the mesh nodes to the node excitation currents, is

$$\textbf{Z}=j\omega\mu h\textbf{K}^{-1}.$$ (4.44)

Figure 4.6 depicts a comparison of the analytic results from (4.18) to the FEM results for the impedance $Z_{in}$ of a port at position $(x=10mm,\,y=10mm)$ on a rectangular cavity with dimensions $L=160mm$, $W=120mm$ and $h=7mm$. Although this FEM simulation was carried out using a very coarse mesh with $\Delta x=\Delta y=10mm$, the results of this simulation match the analytical solution well. Thus, the proposed FEM method provides an efficient and accurate solution for arbitrarily shaped parallel plane cavities. In the case of a large number of excitation ports, the FEM method is even more efficient than the analytical solution. While the analytical method requires the calculation of the double sum term in (4.18) for each coefficient in the port impedance matrix (4.17), the FEM solution requires only one inversion of the sparse system matrix $\textbf{K}$. Linear base functions and a triangular mesh enable an efficient, analytical composition of the system matrix $\textbf{K}$ and the excitation vector $\textbf{F}$ with (4.40) and (4.43) respectively. Since a high accuracy is achieved with the linear base functions even on a coarse mesh, there is no need for higher order base functions. Table 4.2 contains a summary of the FEM equations.

(a) Magnitude comparison.

(b) Phase angle comparison.

Figure 4.6: Comparison of $Z_{in}$ between the analytical solution (4.18) and the FEM solution. Cavity dimensions are ($L=160mm$, $W=120mm$, $h=7mm$), the position of the port for the $Z_{in}$ is $(x=10mm,\,y=10mm)$ and the mesh spacing is $\Delta x=\Delta y=10mm$.

In the following the FEM equations are summarized.

Table 4.2: Summary of the FEM equations.

C. Poschalko: The Simulation of Emission from Printed Circuit Boards under a Metallic Cover