Method of Moments Solver for Metal and Dielectric Structures

Method of Moments computation technique for metal and dielectrics in a PCB.

A PCB consists of a metal part and a dielectric part. The first step in the computational solution of electromagnetic problems is to discretize Maxwell's equations. The process results in this matrix-vector system:

$$V=ZI$$

• V — Applied voltage vector. This signal can be voltage or power applied to the antenna or an incident signal falling at the input port of the PCB component.

• I — Current vector that represents current on the PCB component surface.

• Z — Interaction matrix or impedance matrix that relates V to I. For calculating the interaction matrix, the effect of metal and dielectric parts of the PCB antenna are taken separately.

RF PCB Toolbox™ uses method of moments (MoM) to calculate the interaction matrix and solve system equations.

MoM Formulation

The MoM formulation is split into three parts.

Discretization of Metals and Dielectrics

Discretization enables the formulation from the continuous domain to the discrete domain. This step is called meshing. In the MoM formulation, the metal surface of the PCB is meshed into triangles and the dielectric volume is meshed into tetrahedrons.

Basis Functions

Metals.   To calculate the surface currents on the PCB, first define basis functions. RF PCB Toolbox uses Rao-Wilton-Glisson (RWG) [2] basis functions. The arrows show the direction of current flow.

The basis function includes a pair of adjacent (not necessarily coplanar) triangles and resembles a small spatial dipole with linear current distribution. Each triangle is associated with a positive or negative charge.

For any two triangle patches, $$t_n^{+}$$ and $$t_n^{-}$$, having areas $$A_n^{+}$$ and $$A_n^{-}$$, and sharing common edge $$l_n$$, the basis function is

$$\begin{array}{l}{\overrightarrow{f}}_n(\overrightarrow{r})=\{\begin{array}{cc}\frac{l_n}{2A_n^{+}}{\overrightarrow{\rho }}_n^{+S},& \overrightarrow{r}\in t_n^{+}\\ \frac{l_n}{2A_n^{-}}{\overrightarrow{\rho }}_n^{-S},& \overrightarrow{r}\in t_n^{-}\end{array}\\ \end{array}$$

• $${\overrightarrow{\rho }}_n^{+}=\overrightarrow{r}-{\overrightarrow{r}}_n^{+}$$ — Vector drawn from the free vertex of triangle $$t_n^{+}$$ to observation point $$\overrightarrow{r}$$

• $${\overrightarrow{\rho }}_n^{-}={\overrightarrow{r}}_n^{+}-\overrightarrow{r}$$ — Vector drawn from the observation point to the free vertex of the triangle $$t_n^{-}$$

and

$$\nabla \cdot {\overrightarrow{f}}_n\left(\overrightarrow{r}\right)=\{\begin{array}{cc}\frac{l_n}{A_n^{+}},& \overrightarrow{r}\in t_n^{+}\\ -\frac{l_n}{A_n^{-}},& \overrightarrow{r}\in t_n^{-}\end{array}$$

The basis function is zero outside the two adjacent triangles $$t_n^{+}$$ and $$t_n^{-}$$. The RWG vector basis function is linear and has no flux (no normal component) through its boundary.

Dielectrics.  Basis functions are used to represent unknown quantities. In the case of a PCB component, the unknown quantities are the surface current on the metal structure and flux density due to dielectric volume. RF PCB Toolbox uses Rao-Wilton-Glisson (RWG) [2] basis functions.

For the dielectric volume, RF PCB Toolbox uses a zeroth order edge basis function to model the flux density.

The figure shows an edge-based basis function. The vector variation is perpendicular to the base edge AB (or $$\overline{l}$$). The vector of the edge CD (or $$\overline{p}$$) defines the basis function. Within a tetrahedron, the basis function is a constant field given by

$$\overrightarrow{f}=c\overrightarrow{p}$$

• c – normalization coefficient.

• p – vector of the edge defining the basis function.

Interaction Matrix

Metals.  The interaction matrix is a complex dense symmetric matrix. It is a square N-by-N matrix, where N is the number of basis functions, that is, the number of interior edges in the structure. A typical interaction matrix for a structure with 256 basis functions is shown:

To fill out the interaction matrix, calculate the free-space Green's function between all basis functions on the PCB surface. The final interaction matrix equations are:

$$Z_{mn}=\left(\frac{j\omega \mu }{4\pi }\right){\displaystyle \underset{S}{\int }{\displaystyle \underset{S}{\int }{\overrightarrow{f}}_m}}\left(\overrightarrow{r}\right).{\overrightarrow{f}}_m\left(\overrightarrow{r^{\prime }}\right)gd\overrightarrow{r^{\prime }}d\overrightarrow{r}-\left(\frac{j}{4\pi \omega \epsilon }\right){\displaystyle \underset{S}{\int }{\displaystyle \underset{S}{\int }\left(\nabla .{\overrightarrow{f}}_m\right)}}\left(\nabla .{\overrightarrow{f}}_m\right)gd\overrightarrow{r^{\prime }}d\overrightarrow{r}$$

where

• $$g(\overrightarrow{r},\overrightarrow{r^{\prime }})=\frac{\exp (-jk\left|\overrightarrow{r}-\overrightarrow{r^{\prime }}\right|)}{\left|\overrightarrow{r}-\overrightarrow{r^{\prime }}\right|}$$ — Free-space Green's function

To calculate the interaction matrix, excite the PCB by a voltage of 1 V at the input port. So the voltage vector has zero values everywhere except at the feeding edge. Solve the system of equations to calculate the unknown currents. Once you determine the unknown currents, you can calculate the field and surface properties of the PCB.

Dielectrics.  The interaction matrix is a complex dense symmetric matrix. To fill out the interaction matrix, calculate the free-space Green's function between all the basis functions on the PCB surface. The final interaction matrix equations are:

• ZMM – metal to metal interaction. For a pure metal structure, you only calculate this symmetric square matrix.

  $$Z_{mn}^{MM}=\left(\frac{j\omega \mu }{4\pi }\right){\displaystyle \underset{S}{\int }{\displaystyle \underset{S}{\int }{\overrightarrow{f}}^M{}_m}}\left(\overrightarrow{r}\right).{\overrightarrow{f}}^M{}_n\left(\overrightarrow{r^{\prime }}\right)gd\overrightarrow{r^{\prime }}d\overrightarrow{r}-\left(\frac{j}{4\pi \omega \epsilon }\right){\displaystyle \underset{S}{\int }{\displaystyle \underset{S}{\int }\left(\nabla .{\overrightarrow{f}}^M{}_m\right)}}\left(\nabla .{\overrightarrow{f}}^M{}_n\right)gd\overrightarrow{r^{\prime }}d\overrightarrow{r}$$

• ZDD – dielectric to dielectric interaction. For pure dielectric structures, you only calculate this symmetric square matrix.

  $$\begin{array}{l}{\widehat{Z}}_{mn}^{DD}={\displaystyle \sum _{p-1}^P{\displaystyle \sum _{p^{\prime }-1}^{P^{\prime }}\frac{K_p}{{\widehat{\epsilon }}_p}}}{\displaystyle \underset{V_D}{\int }{\overrightarrow{f}}_{mp}(\overrightarrow{r})\cdot }{\overrightarrow{f}}_{n{p}^{\prime }}(\overrightarrow{r})d\overrightarrow{r}\\ -\frac{{\omega }^2{\mu }_0}{4\pi }{\displaystyle \sum _{p-1}^P{\displaystyle \sum _{p^{\prime }-1}^{P^{\prime }}{K}_p{K}_{p^{\prime }}}{\displaystyle \underset{V_D}{\int }{\displaystyle \underset{V_{D^{\prime }}}{\int }g(}}}\overrightarrow{r},\overrightarrow{r^{\prime }}){\overrightarrow{f}}_{mp}(\overrightarrow{r})\cdot {\overrightarrow{f}}_{n{p}^{\prime }}(\overrightarrow{r^{\prime }})d\overrightarrow{r}d\overrightarrow{r^{\prime }}\\ -\frac{1}{4\pi {\epsilon }_0}{\displaystyle \sum _{q-1}^Q{\displaystyle \sum _{q^{\prime }-1}^{Q^{\prime }}{\widehat{K}}_q}}{\widehat{K}}_{q^{\prime }}{\displaystyle \underset{{\Omega }_q}{\int }{\displaystyle \underset{{\Omega }_{q^{\prime }}}{\int }g(\overrightarrow{r},\overrightarrow{r^{\prime }})f_{\perp mq}(\overrightarrow{r})}}{f}_{\perp n{q}^{\prime }}(\overrightarrow{r^{\prime }})d_s{d}_{s^{\prime }}m,n=1,\mathrm{....},N\end{array}$$

• ZMD and ZDM – These matrices calculate the interaction between metal and dielectric. This matrix is not a symmetrical square matrix.

  $$\begin{array}{l}{Z}_{mn}^{MD}=-\frac{{\omega }^2{\mu }_0}{4\pi }{\displaystyle \sum _{p-1}^2{\displaystyle \sum _{p^{\prime }-1}^{P^{\prime }}{K}_p}}{\displaystyle \underset{t}{\int }{\displaystyle \underset{V_{D^{\prime }}}{\int }{\overrightarrow{f}}_n^M(\overrightarrow{r})\cdot }{\overrightarrow{f}}_{m{p}^{\prime }}(\overrightarrow{r^{\prime }})g(\overrightarrow{r},\overrightarrow{r^{\prime }})d\overrightarrow{r^{\prime }}}{d}_s\\ -\frac{1}{4\pi {\epsilon }_0}{\displaystyle \sum _{p-1}^2{\displaystyle \sum _{q-1}^Q{\widehat{K}}_q}}{\displaystyle \underset{t_D}{\int }{\displaystyle \underset{{\Omega }_{q^{\prime }}}{\int }({\nabla }_s\cdot {\overrightarrow{f}}_n^M(\overrightarrow{r}))f_{\perp mq}(\overrightarrow{r^{\prime }})g(\overrightarrow{r},\overrightarrow{r^{\prime }})}}{d}_{{\Omega }^{\prime }}{d}_s\\ m=1,\mathrm{....},N_D;n=1,\mathrm{...},N_M\end{array}$$

  $$\begin{array}{l}{Z}_{mn}^{DM}=-\frac{j\omega {\mu }_0}{4\pi }{\displaystyle \sum _{p-1}^2{\displaystyle \sum _{p^{\prime }-1}^{P^{\prime }}{K}_{p^{\prime }}}}{\displaystyle \underset{V_D}{\int }{\displaystyle \underset{S_{D^{\prime }}}{\int }{\overrightarrow{f}}_{np}(\overrightarrow{r})\cdot {\overrightarrow{f}}_m^M(\overrightarrow{r^{\prime }})g(\overrightarrow{r},\overrightarrow{r^{\prime }})d_{s^{\prime }}d\overrightarrow{r}}}\\ +\frac{1}{4\pi {\epsilon }_0\omega }{\displaystyle \sum _{p-1}^2{\displaystyle \sum _{q-1}^Q{\widehat{K}}_q}}{\displaystyle \underset{{\Omega }_q}{\int }{\displaystyle \underset{S_{D^{\prime }}}{\int }{f}_{\perp nq}(\overrightarrow{r})\cdot ({\nabla }_s\cdot {\overrightarrow{f}}_m^M(\overrightarrow{r^{\prime }}))g(\overrightarrow{r},\overrightarrow{r^{\prime }})}}{d}_{s^{\prime }}{d}_{\Omega }\\ m=1,\mathrm{....},N_D;n=1,\mathrm{...},N_M\end{array}$$

where

• $$g(\overrightarrow{r},\overrightarrow{r^{\prime }})=\frac{\exp (-jkR)}{R},R=\left|\overrightarrow{r}-\overrightarrow{r^{\prime }}\right|$$ is the free space Green's function.

• $$K=\frac{{\widehat{\epsilon }}^{\pm }-{\epsilon }_0}{{\widehat{\epsilon }}^{\pm }}$$ is the complex dielectric constant within every tetrahedron.

• $${\widehat{K}}_q=K_{+}-K_{-}$$ is the differential contrast on every face of the tetrahedron.

For a composite metal structure, you must calculate all four matrices.

Neighbor Region

From the interaction matrix plot, you observe that the matrix is diagonally dominant. As you move further away from the diagonal, the magnitude of the terms decreases. This behavior is same as the Green's function behavior. The Green's function decreases as the distance between r and r' increases. Therefore, it is important to calculate the region on the diagonal and close to the diagonal accurately.

This region on and around the diagonal is called neighbor region. The neighbor region is defined within a sphere of radius R, where R is in terms of triangle size. The size of a triangle is the maximum distance from the center of the triangle to any of its vertices. By default, R is twice the size of the triangle. For better accuracy, a higher-order integration scheme is used to calculate the integrals. The figure shows a typical interaction matrix for a metal structure ZMM with 256 basis functions.

The dielectric interaction matrix is also diagonally dominant. As you move further away from the diagonal, the magnitude of the terms decreases. This behavior is same as the Green's function behavior. The Green's function decreases as the distance between r and r' increases. Therefore, it is important to calculate the region on the diagonal and close to the diagonal accurately. For a dielectric the neighbor region is based on the average size of the tetrahedron.

Singularity Extraction

Along the diagonal, r and r' are equal and defines Green's function becomes singular. To remove the singularity, extraction is performed on these terms.

$$\begin{array}{l}{\displaystyle \underset{t_p}{\int }{\displaystyle \underset{t_q}{\int }\left({\overrightarrow{\rho }}_i.{{\overrightarrow{\rho }}^{\prime }}_j\right)}}g(\overrightarrow{r},\overrightarrow{r^{\prime }})ds\text{'}ds={\displaystyle \underset{t_p}{\int }{\displaystyle \underset{t_q}{\int }\frac{\left({\overrightarrow{\rho }}_i.{{\overrightarrow{\rho }}^{\prime }}_j\right)}{\left|\overrightarrow{r}-\overrightarrow{r^{\prime }}\right|}}}ds\text{'}ds+{\displaystyle \underset{t_p}{\int }{\displaystyle \underset{t_q}{\int }\frac{\left(\exp \left(-jk\left|\overrightarrow{r}-\overrightarrow{r^{\prime }}\right|\right)-1\right)\left({\overrightarrow{\rho }}_i.{{\overrightarrow{\rho }}^{\prime }}_j\right)}{\left|\overrightarrow{r}-\overrightarrow{r^{\prime }}\right|}ds\text{'}ds}}\\ {\displaystyle \underset{t_p}{\int }{\displaystyle \underset{t_q}{\int }g(\overrightarrow{r},\overrightarrow{r^{\prime }})ds\text{'}ds}}={\displaystyle \underset{t_p}{\int }{\displaystyle \underset{t_q}{\int }\frac{1}{\left|\overrightarrow{r}-\overrightarrow{r^{\prime }}\right|}}}ds\text{'}ds+{\displaystyle \underset{t_p}{\int }{\displaystyle \underset{t_q}{\int }\frac{\left(\exp \left(-jk\left|\overrightarrow{r}-\overrightarrow{r^{\prime }}\right|\right)-1\right)}{\left|\overrightarrow{r}-\overrightarrow{r^{\prime }}\right|}ds\text{'}ds}}\end{array}$$

The two integrals on the right-hand side of the equations, called potential or static integrals are found using analytical results [3].

The equations for the singularity extraction of the ZDD matrix are:

$$\begin{array}{l}{\displaystyle \underset{V_D}{\int }{\displaystyle \underset{V_{D^{\prime }}}{\int }g(\overrightarrow{r},\overrightarrow{r^{\prime }})d\overrightarrow{r}d\overrightarrow{r^{\prime }}}}={\displaystyle \underset{V_D}{\int }{\displaystyle \underset{V_q}{\int }\frac{1}{\left|\overrightarrow{r}-\overrightarrow{r^{\prime }}\right|}}}d\overrightarrow{r}d\overrightarrow{r^{\prime }}+{\displaystyle \underset{V_D}{\int }{\displaystyle \underset{V_q}{\int }\frac{\left(\exp \left(-jk\left|\overrightarrow{r}-\overrightarrow{r^{\prime }}\right|\right)-1\right)}{\left|\overrightarrow{r}-\overrightarrow{r^{\prime }}\right|}d\overrightarrow{r}d\overrightarrow{r^{\prime }}}}\\ {\displaystyle \underset{{\Omega }_q}{\int }{\displaystyle \underset{{\Omega }_{q^{\prime }}}{\int }g(\overrightarrow{r},\overrightarrow{r^{\prime }})d\Omega }}d\Omega \text{'}={\displaystyle \underset{{\Omega }_D}{\int }{\displaystyle \underset{{\Omega }_q}{\int }\frac{1}{\left|\overrightarrow{r}-\overrightarrow{r^{\prime }}\right|}}}d\Omega d\Omega \text{'}+{\displaystyle \underset{S_D}{\int }{\displaystyle \underset{S_q}{\int }\frac{\left(\exp \left(-jk\left|\overrightarrow{r}-\overrightarrow{r^{\prime }}\right|\right)-1\right)}{\left|\overrightarrow{r}-\overrightarrow{r^{\prime }}\right|}d\Omega d\Omega }}\end{array}$$

References

[1] Harringhton, R. F. Field Computation by Moment Methods. New York: Macmillan, 1968.

[2] Rao, S. M., D. R. Wilton, and A. W. Glisson. “Electromagnetic scattering by surfaces of arbitrary shape.” IEEE. Trans. Antennas and Propagation, Vol. AP-30, No. 3, May 1982, pp. 409–418.

[3] Wilton, D. R., S. M. Rao, A. W. Glisson, D. H. Schaubert, O. M. Al-Bundak. and C. M. Butler. “Potential Integrals for uniform and linear source distribution on polygonal and polyhedral domains.” IEEE. Trans. Antennas and Propagation. Vol. AP-30, No. 3, May 1984, pp. 276–281.