Computer Simulation of High‐Frequency Electromagnetic Fields

Abstract

High‐frequency and microwave electromagnetic fields are used in billions of various devices and systems. Design of these systems is impossible without detailed analysis of their electromagnetic field. Most of microwave systems are very complex, so analytical solution of the field equations for them is impossible. Therefore, it is necessary to use numerical methods of field simulation. Unfortunately, such complex devices as, for example, modern smartphones cannot be accurately analysed by existing commercial codes. The chapter contains a short review of modern numerical methods for Maxwell's equations solution. Among them, a vector finite element method is the most suitable for simulation of complex devices with hundreds of details of various forms and materials, but electrically not too large. The method is implemented in the computer code radio frequency simulator (RFS). The code has friendly user interface, an advanced mesh generator, efficient solver and post‐processor. It solves eigenmode problems, driven waveguide problems, antenna problems, electromagnetic‐compatibility problems and others in frequency domain.

1. Introduction

High‐frequency electromagnetic fields are used now in telecommunications and radar systems, astrophysics, plasma heating and diagnostics, biology, medicine, technology and many other applications. Special electromagnetic systems excite and guide these fields with given time and space distribution. A designer or a researcher of such systems ought to know in detail their electromagnetic field characteristics. This goal can be achieved or by experimental study, often too long in time and expensive, or by computer simulation. The last choice becomes more and more preferable with fast progress of computational electrodynamics and computer efficiency.

All macroscopic electromagnetic phenomena are governed by Maxwell's equations. Unfortunately, these remarkable equations have so many solutions, that choice of the one satisfying given initial and boundary conditions (BC) often becomes very difficult problem. A number of commercial computer codes based on numerical solution of Maxwell's equations are available at this time. These codes make possible high‐frequency electromagnetic field simulation.

Researches have not yet created a universal code, efficiently simulating electromagnetic field excited by an arbitrary technical or nature source in an arbitrary medium. Some codes are more suitable for solving one kind of problems and other codes—another kind. Hence, developing new, more universal and efficient computer codes is an actual task. On the other side, a designer of electromagnetic devices and systems has to choose most efficient computer code for solving his/her particular problem. He/she can do right choice only if he/she understands the basics of a numerical method used in the given code.

The goal of the presented chapter is to formulate electromagnetic problems and to describe in short prevailing numerical methods of its solving. As an example, the chapter also gives more detailed description of the radio frequency simulator (RFS) computer code, developed in collaboration of Saint‐Petersburg State Electrotechnical University and LG Russian R&D Centre. Some results obtained by means of this code demonstrate its accuracy and efficiency.

The author hopes that this chapter would be useful for researchers and designers of modern telecommunication devices and systems.

2. Basic equations

Maxwell's equations are the basic ones, describing macroscopic electromagnetic fields in an arbitrary medium. In modern notation, these equations have the form

$\nabla \times H=\frac{\partial D}{\partial t}+J$E1

$\nabla \times E=-\frac{\partial B}{\partial t}$E2

$\nabla \cdot D=\rho$E3

$\nabla \cdot B=0$E4

In these formulas, $J,\rho$ are the electric current and electric charge densities (field sources), $E,H$ are the electric and magnetic field intensities (or simply electric and magnetic fields), $D,B$ are the electric and magnetic flux densities, $\nabla$ is Hamilton's differential operator and $\times$ is the sign of vector and scalar products.

The constitutive relations couple flux densities and field intensities:

$D=\epsilon E={\epsilon }_0{\epsilon }_{r}E$E5

$B=\mu H={\mu }_0{\mu }_{r}H$E6

Here, $\epsilon ,\mu$ are the absolute permittivity and permeability, ${\epsilon }_0={10}^7/(4\pi c^2),{\mu }_0=4\pi \cdot {10}^{-7}$ are the dielectric and magnetic constants (we use the SI units in this chapter) and ${\epsilon }_r,{\mu }_r$ are the relative permittivity and permeability, which can be scalars or tensors, depending on medium properties. Equations (1), (2), (5) and (6) form a system of 12 scalar differential equations of the first order with 12 unknowns—components of $E,H,D,B$ vectors, which are functions of space coordinates and time. To solve this system, one needs to define initial and boundary conditions. These conditions together with the equations form an electrodynamics problem (EMP).

Initial conditions define electric and magnetic field intensities in the computational region at the initial moment of time.

The most frequently used boundary conditions (BC) are as follows:

• On the surface separating two dielectrics:

  $n\times (E_2-E_1)=0;n\times (H_2-H_1)=J_s$E7

  where $n$ is the unit normal to the surface directed from the first medium to the second and $J_s$ is the surface electric current density.

• On the surface of perfect electric conductor (PEC):

  $n\times E=0$E8

• Similarly, on the surface of perfect magnetic conductor (PMC)

  $n\times H=0$E9

• Approximate Leontovich boundary condition holds on the impedance surface:

  $e_n\times \left(e_n\times \mathbf{Ė}\right)-j\frac{Z_s}{k_0{\eta }_0{\mu }_r}{e}_n\times \left(\nabla \times \mathbf{Ė}\right)=0$E10

  where $Z_s$ is the surface impedance. For metals, $Z_s=(1+j)\sqrt{\omega \mu /(2\sigma )}$.

• We define radiation (adsorption) condition on the surface throw where radiation propagates without reflections. Higdon [1] proposed general absorption boundary conditions (ABCs) theory of arbitrary order of approximation. For a wave, propagating under arbitrary angle $\varphi$ to the border $x=\mathrm{Const}$, the first‐order Higton's ABC are

  $\left(\frac{\partial }{\partial x}+\frac{\cos \varphi }{u}\frac{\partial }{\partial t}+\xi \right)E_{\tau }=0$E11

where $\xi$ is constant, providing solution stability and $u$ is the phase velocity of the wave. Higher order ABC can be constructed by multiplying several first‐order operators (11).

Another expression for ABC can be derived from the wave equation. For the wave propagating along the normal $e_n$ to the border, we get

${\left(\nabla \times \mathbf{Ė}\right)}_{\tau }+j{k}_0\sqrt{{\epsilon }_r{\mu }_r}{e}_{n\times }\times \mathbf{Ė}=0$E12

There are two types of EMPs: an inner problem, when the solution is defined in a closed space region with certain boundary conditions on its border, and an outer problem presuming existence of the solution in the unbounded space excluding some regions with prescribed conditions on their boundaries. We ought also to define field sources in the computational region and initial conditions—electrical and magnetic fields at some moment of time. The proper initial and boundary conditions guarantee existence and uniqueness of the solution [2, 3].

Equations (1)–(4) presume arbitrary time dependence of sources and fields. But very often field dependence on time expresses by the harmonic law: $a=a_m\cos (\omega t+\varphi )$, where $a$ is any component of the field, $\omega =2\pi f$ is the angular frequency, $f$ is frequency and $\varphi$ is the initial angle. In this case, we can simplify Maxwell's equations, using notation

$E(r,t)=\mathrm{Re}\left[\mathbf{Ė}(r){\mathrm{e}}^{j\omega t}\right],$E13

where $\mathbf{Ė}(r)=E_x{e}^{j{\varphi }_x}{e}_x+E_y{e}^{j{\varphi }_y}{e}_y+E_z{e}^{j{\varphi }_z}{e}_z$ is the complex electric field amplitude (phasor). Magnetic field $H$ is presented similarly. Using these notations, we can write Maxwell's equations for phasors:

$\nabla \times \dot{H}=j\omega \dot{\epsilon }\mathbf{Ė}+\dot{J}$E14

$\nabla \times \mathbf{Ė}=-j\omega \dot{\mu }\dot{H}$E15

$\nabla \cdot (\dot{\epsilon }\mathbf{Ė})=\dot{\rho }$E16

$\nabla \cdot (\dot{\mu }\dot{H})=0$E17

Here, $\dot{\epsilon }={\epsilon }_0{\dot{\epsilon }}_r,\dot{\mu }={\mu }_0{\dot{\mu }}_r$ are the complex absolute permittivity and permeability and ${\dot{\epsilon }}_r=({{\epsilon }^{\prime }}_r+j{{\epsilon }^{″}}_r)$, ${\dot{\mu }}_r=({{\mu }^{\prime }}_r+j{{\mu }^{″}}_r)$ are the relative complex permittivity and permeability. The equation systems (14)–(17) are simpler than the original one because it does not contain time as independent variable and has only six unknown functions.

3. Basis stages of electromagnetic problem solution

Solution of a given electromagnetic problem can be divided on several subsequent steps (stages):

1. Problem formulation—defining the goal of the computation, necessary input and output data, admissible inaccuracy of the results.

2. Analytical treatment—formulating of equations, initial and boundary conditions, geometrical description of the computational region and filling medium properties. The choice of the numerical solution method, transforming equations to the form, most suitable for the chosen method, a‐priory analysis of equations and their solutions properties.

3. Problem discretization—transfer from continuous functions to discrete ones and from functional equations to the system of linear algebraic equations (SLAE), in the certain sense approaching the initial problem.

4. Algebraic solution—choosing the most efficient numerical method and solving the SLAE with the prescribed accuracy.

5. Post‐processing—calculation of fields, characteristics and parameters of the electromagnetic system and visualization of the results.

Each stage of the solution adds its own contribution to the total solution error. The first step adds the so‐called inherent error arising due to inaccuracy of input data. This error cannot be eliminated on the next stages of solution. The second stage adds mathematical model error caused by imperfect adequateness of the model to the real physical process. Problem discretization adds the so‐called numerical method error, the value of which depends on the quality of the discretization process. At last, computational error arises on stages 4 and 5 due to finite accuracy of numbers presentation in a computer and finite number of operations. With progress in the computational mathematics and computer, the main error sources move from the uncontrolled first stages to the fourth and fifth stages, where we can often predict the error value a‐priory.

4. Numerical methods classification

The solution of real‐life electromagnetic problem is a very complicated task. There are no universal methods, capable to solve efficiently an arbitrary problem. Hence, a number of numerical methods were elaborated, each of them is the most efficient for its particular range of problems. The numerical methods divide on two large groups.

The first group solves problems involving Maxwell's equation (1)–(4), containing time as an independent variable. This group makes possible to find the solution in time domain (TD) with arbitrary time dependence of fields. This group of methods is the most suitable for solving non‐linear problems. If the problem is linear, Fourier transform can be used to find frequency spectrum of the solution.

The most popular method of this group is the Finite Difference Time Domain (FDTD) method, proposed by Yee [4, 5]. Another method, successfully used in time domain, is the Finite Integration Technique (FIT), firstly proposed by Weiland [6]. The Transmission Line Matrix (TLM) method, proposed by Johns and Beurle [7], also works efficiently in time domain. Its detailed description can be found in [8]. These methods were implemented in a number of computer codes, such as SPEAG SEMCAD™ (FDTD), CST Studio Suite™ (FIT, TLM) and many others.

The second group deals with Eqs. (14)–(17) for field phasors, supposing harmonic time dependence of fields. This group provides solution in Frequency (or Spectral) Domain (FD). The finite element method (FEM) is one of the most efficient representatives of this group. The first application of this method to the solution of mechanical problems refers to the year 1943 [9]. The book [10] contains detailed description of the method, which is rather universal and accurate. The Method of Moments (MoMs) and its varieties are also frequently used in FD. In contrast with FEM, MoM uses integral form of basic equations, where electric current density distribution on conducting surfaces excited by external sources is unknown. The book [11] reflects modern state of the MoM. Mentioned methods were implemented in codes ANSIS HFSS™ (FEM), Altair FEKO™ (MoM, FEM) and other commercial computer codes. Of course, there exist many other numerical methods and various implementations.

Most of modern methods allow implementation both in TD and in FD.

Because it is impossible to give detailed representation of all numerical methods in a limited space, we give here more detailed information about the FEM method, which was implemented in the computer code radio frequency simulator (RFS) [12].

5. Finite element method for electromagnetics

5.1. Main features of the method

The finite element method belongs to the variational methods of solving partial differential equation (PDE). It presumes formulating a functional, which is stationary (has minimum or maximum value) on the equation solution. In order to find functional extremum, the computational region is partitioned on a number of subregions (finite elements). After that, we approximate an unknown function in each finite element (FE) by superposition of basis functions. Basis functions have to be simple and form a linearly independent system. Applying Ritz method or Galerkin algorithm, we fulfil discretization of the problem, that is, transition from the partial difference equation to the system of linear algebraic equations (SLAE). Numerical solution of the SLAE gives unknown coefficients of basic functions, which are used to restore electromagnetic field. At last, needed components of electromagnetic field and system parameters are calculated.

5.2. Mesh generation

Partitioning of the computational region (mesh generation) is the first stage of problem solution by FEM. It supposes dividing the computation region on a set of subregions—finite elements. FEs must densely fill the region and be nearly conformal to its border. In contrast to FDTD or FIT methods, a FEM mesh can be irregular and contain FEs of different forms. Tetrahedron FEs are used most commonly, because they allow dense packing and quite correctly approximate curvilinear borders. However, mesh generation for regions with complex forms filled with different materials is a very complex task.

Delaunay tessellation is a common way of mesh building. It includes several stages. Firstly, a surface triangle mesh is generated. Then, a volume mesh based on the surface mesh and covering the whole computational region is being built. The method can build a rather good mesh without self‐intersections. Some commercial codes, such as SYMMETRIX MeshSym™, can be used to generate mesh by this procedure.

Unfortunately, this method cannot be applied to tessellation of complex geometrical models consisting of hundreds of parts made from various materials. Usually, an electromagnetic simulation code imports such models from various CAD systems (see, e.g. Figure 1). As a rule, imported models contain a number of errors caused by insufficient attention of a designer or arising in the process of graphic formats transforms. As a result, mesh generator fails to build the mesh.

Figure 1.

A part of the handset CAD model.

On the other side, CAD models are often excessively detailed. They contain peculiarities, not influenced on electromagnetic field. Figure 2 shows an example of such extra detailed CAD model. The application of standard mesh generator to such a model can result in excess mesh size. Correction and simplification of the model manually needs several working days of the qualified engineer.

Figure 2.

A part of the CAD model with excessive details.

An advanced method of mesh generation was developed and implemented in the RFS simulation code. The algorithm begins from assigning materials and attributes to the model parts. Most important objects, such as ports, printed circuit board (PCB) and antennas, are labelled as ‘electrically important'. This is the only stage of the algorithm, which is made manually. This stage can be omitted for simple models. Then the code builds surface meshes for each detail of the system. The third stage presumes elimination of individual meshes interceptions and building united surface mesh. Electrically important, metal parts and parts with higher permittivity have priority in this process. At last, a global volume mesh is generated [13]. The user can control mesh quality (maximum to minimum ratio of tetrahedron's dimensions and other mesh parameters).

We proved this algorithm on more than 190 models of mobile handsets and showed its 100% reliability. Generated meshes contained more than 1.5 billion tetrahedrons. Figure 3 shows a CAD model of the handset (a), surface mesh (b) and volume mesh (c), built by the described algorithm.

Figure 3.

CAD model of the handset (a), united surface mesh (b) and volume mesh (c).

5.3. Basic equation

Consider a region $V$, delimited by a surface $S$, where electromagnetic field has to be calculated. We suppose that the region $V$ is filled by linear medium. Applying curl operator to Eq. (15) and substituting into result Eq. (14), we get PDE of the second order

$\nabla \times ({\dot{\mu }}_r^{-1}\nabla \times \mathbf{Ė})+j\sigma {\eta }_{0}k\mathbf{Ė}-{\dot{\epsilon }}_r{k}^2\mathbf{Ė}=-jk{\eta }_0{\dot{J}}^{imp},$E18

where ${\dot{J}}^{imp}$ is the imposed current density, $k=\omega /c$ is wave number, $c={({\epsilon }_0{\mu }_0)}^{-1/2}$ is the light velocity in free space, $\sigma$ is medium conductivity and ${\eta }_0=\sqrt{{\mu }_0/{\epsilon }_0}=120\pi$ ohm is the intrinsic impedance of free space. After solving this equation, we can find magnetic field by means of Eq. (15):

$\dot{H}=\frac{j}{k{\eta }_0{\dot{\mu }}_r}\nabla \times \mathbf{Ė}.$E19

To solve Eq. (16), we divide the computation region on finite elements and approximate the field in each FE by superposition of basis functions. It is natural to approximate vector unknown function in Eq. (16) by a set of vector basis functions. We name such variant of FEM as vector finite element method (VFEM).

Nedelec [14] proposed a set of vector basis functions for triangle and tetrahedral FE. These functions are associated with tetrahedron edges and have the form

$w_m={\zeta }_{m1}\nabla {\zeta }_{m2}-{\zeta }_{m2}\nabla {\zeta }_{m1},$E20

where ${\zeta }_{mn}$ is the barycentric function of the vertex (node) node $n$ of the edge $m$ ($m=1,\dots ,6,n=1,2$). The barycentric functions are polynomials of the first order. For a given Descartes coordinate system

${\zeta }_p=a_p+b_{p}x+c_{p}y+d_{p}z,$E21

where $p$ is vertex number. As a tetrahedron has six edges, Nedelec's basis consists of six functions. Figure 4 shows the numbering scheme for nodes and edges, needed for the correct use of Nedelec's functions.

Figure 4.

Numbering order of a tetrahedron's vertices and edges.

Coefficients in Eq. (17) can be found from the equation

$\left[\begin{array}{cccc}{b}_1& c_1& d_1& a_1\\ b_2& c_2& d_2& a_2\\ b_3& c_3& d_3& a_3\\ b_4& c_4& d_4& a_4\end{array}\right]={\left[\begin{array}{cccc}{x}_1& x_2& x_3& x_4\\ y_1& y_2& y_3& y_4\\ z_1& z_2& z_3& z_4\\ 1& 1& 1& 1\end{array}\right]}^{-1}.$E22

where $x_p,y_p,z_p,p=1\dots 4$ are the coordinates of the p‐th tetrahedron's node.

Nedelec's basis functions have the following properties:

1. Dimension of these functions is $L^{-1}$.

2. They have zero divergence and their curl is constant inside the FE. This property provides the absence of spurious solutions of the problem.

3. Projection of the function $w_m$ on the edge $m$ is constant, and its projection on the other edges is equal to zero. Hence, these functions provide continuity of tangential electric field on the border between neighbouring tetrahedrons.

4. Circulation of function $w_m$ along its edge is equal to unity.

Full first‐order polynomial set for electric field approximation contains 12 functions (3 E‐vector projections, each needs four functions). However, Nedelec's set contains only six basis functions. Hence, it does not form full first‐order polynomial system. The order of approximation by these functions is lower than one. Finite elements with these functions are called FE of order 1/2 or CT/LN (constant tangential/linear normal) according to the behaviour of these functions along edges and in normal direction to them.

Striving for increasing approximation accuracy led to creating another set of basis functions with higher order of approximation. A set of 20 functions (12 associated with edges and eight with faces) forms FE of order 1.5 or LN/QN (linear tangential/quadratic normal) [15, 16]. There exist basis functions of order 2.5 and higher, but they are rarely used.

With the aim of simplicity, we consider further only Nedelec's basis functions. So, electric field intensity in q‐th tetrahedron is approximated as

${\mathbf{Ė}}^q(r)={\displaystyle \sum _{n=1}^6{\dot{x}}_n^q{w}_n^q}(r),q=1,\dots ,N_q,$E23

where $N_q$ is the total number of tetrahedrons.

Coefficients ${\dot{x}}_m^q$ have the sense of voltage on edge (with reverse sign). Really, taking integral along an edge

$V_m^q=-{\displaystyle {\int }_{L_q}{\mathbf{Ė}}^q\cdot d{l}_m}=-{\displaystyle \sum _{n=1}^6{\dot{x}}_n^q{w}_n^{q}d{l}_m}=-{\dot{x}}_m^q,$E24

due to the fourth property of basis functions.

Let us consider q‐th tetrahedron. According to Galerkin's method, we multiply Eq. (16) by basis function $w_m^q$, belonging to this tetrahedron and integrate the product over the tetrahedron volume $V_q$. As a result, we obtain the so‐called week form of Eq. (16):

$\begin{array}{l}\underset{V_q}{{\displaystyle \int }}\left[\nabla \times \left({\dot{\mu }}_r^{-1}\nabla \times E^q\right)\right]\cdot w_m^{q}dV+j{\eta }_{0}k\underset{V_q}{{\displaystyle \int }}\sigma E^q\cdot w_m^{q}dV-k^2\underset{V_q}{{\displaystyle \int }}{\dot{\epsilon }}_r{E}^q\cdot w_m^{q}dV\hfill \\ =-jk{\eta }_0\underset{V_q}{{\displaystyle \int }}{J}^{imp}\cdot w_m^{q}dV.\hfill \end{array}$

Transformation of the first term of this equation and subsequent substitution into result approximation (19) gives a system of linear algebraic equations:

$\begin{array}{l}{\displaystyle \sum }_{n=1}^6{\dot{x}}_n^q[\underset{V_q}{{\displaystyle \int }}\left({\dot{\mu }}_r^{-1}\nabla \times w_n^q\right)\cdot \left(\nabla \times w_m^q\right)dV+jk\underset{V_q}{{\displaystyle \int }}{\eta }_0\sigma w_n^q\cdot w_m^{q}dV-k^2\underset{V_q}{{\displaystyle \int }}{\dot{\epsilon }}_r{w}_n^q\cdot w_m^{q}dV]\hfill \\ =-\underset{S_q}{{\displaystyle \oint }}\left[\left({\dot{\mu }}_r^{-1}\nabla \times {Ė}^q\right)\times w_m^q\right]dS-jk\underset{V_q}{{\displaystyle \int }}{\eta }_0{\dot{J}}^{imp}\cdot w_m^{q}dV,m=1,2,\dots ,6;q=1,2,\dots ,N_q,\hfill \end{array}$E25

where $N_q$ is the total number of FE in the computation region.

The system of Eq. (25) can be written in matrix form:

${\overline{Q}}^q{X}^q\equiv \left({\overline{T}}^q+jk{\overline{U}}^q-k^2{\overline{R}}^q\right){\overline{X}}^q=-{\overline{S}}^q-jk{\overline{B}}^q,$E26

where $Q^q,{\overline{T}}^q,{\overline{U}}^q,{\overline{R}}^q$ are square matrices having dimension 6 × 6, ${\overline{S}}^q,{\overline{B}}^q$ are column vectors having six components, dependent on the given boundary conditions and implicit current density in the q‐th finite element. Depending on BC, $S^q$ can sometimes be a square 6 × 6 matrix.

Equation (26) is called the local matrix equation, and matrix $Q^q$ is the local matrix of q‐th FE. Such matrices should be built for every FE in the computational region.

Components of matrices ${\overline{T}}^q,{\overline{U}}^q,{\overline{R}}^q$ and vectors ${\overline{S}}^q,{\overline{B}}^q$ are defined by the expressions

$t_{mn}={\dot{\mu }}_r^{-1}{\displaystyle {\int }_V(\nabla \times w_m)\cdot (\nabla \times w_n)dV;}$E27

$u_{mn}={\eta }_0\sigma {\displaystyle {\int }_V{w}_m\cdot w_{n}dV;}$E28

$r_{mn}={\dot{\epsilon }}_r{\displaystyle {\int }_V{w}_n\cdot w_{m}dV;}$E29

$s_m={\dot{\mu }}_r^{-1}{\displaystyle {\oint }_S\left[\left(\nabla \times \mathbf{Ė}\right)\times w_m\right]dS;}$E30

$b_m={\eta }_0{\displaystyle {\int }_V{\dot{J}}^{imp}\cdot w_{m}dV.}$E31

Medium in these expressions is supposed to be homogeneous inside the FE. Index ‘q' denoted FE number is omitted. We can see that in order to calculate elements (27)–(31), it is sufficient to calculate two integrals:

$f_{mn}={\displaystyle {\int }_V(\nabla \times w_n)\cdot }(\nabla \times w_m)dV\text{ and }{h}_{mn}={\displaystyle {\int }_V(w_n}\cdot w_m)dV.$E32

We use the next algorithm for the calculation of matrices and vectors elements:

For each tetrahedron

• Coefficients of the barycentric functions are calculated by means of Eq. (22).

• Auxiliary values are calculated:

  $v_{ij}=\nabla {\zeta }_i\times \nabla {\zeta }_j=(c_i{d}_j-c_j{d}_i)\stackrel{⌢}{x}+(d_i{b}_j-d_j{b}_i)\stackrel{⌢}{y}+(b_i{c}_j-b_j{c}_i)\stackrel{⌢}{z};$E33
  ${\varphi }_{ij}=\nabla {\zeta }_i\cdot \nabla {\zeta }_j=b_i{b}_j+c_i{c}_j+d_i{d}_j.$E34

  It can be easily shown that $\nabla \times w_m=2v_{m1,m2}$.

• We introduce matrix $\overline{G}$ with components

  $g_{ij}=V^{-1}{\displaystyle {\int }_V{\zeta }_i{\zeta }_{j}dV.}$E35

  Using integration formula for three‐dimensional (3D) simplex, we can prove that

  $\overline{G}=\frac{1}{20}\left|\begin{array}{cccc}2& 1& 1& 1\\ 1& 2& 1& 1\\ 1& 1& 2& 1\\ 1& 1& 1& 2\end{array}\right|.$E36

• Integrals (32) are calculated:

  $f_{mn}=4V(V_{m1,m2}\cdot V_{n1,n2})dV;$E37
  $h_{mn}=V({\varphi }_{m2,n2}{g}_{m1,n1}-{\varphi }_{m2,n1}{g}_{m1,n2}-{\varphi }_{m1,n2}{g}_{m2,n1}+{\varphi }_{m1,n1}{g}_{m2,n2}),m,n=1,\dots ,6.$E38

In these expressions, $V$ is the volume of the tetrahedron.

Construction of the global matrix (assembling) for all tetrahedrons is a rather complex task, because many nodes, edges and faces are common for neighbour FEs (see Figure 5, which shows two neighbour tetrahedrons, separated for reader's convenience). They have three common vertices and three common edges. Hence, the united matrix for two FEs has only 9 dimensions instead of 12.

Figure 5.

Two neighbour tetrahedrons. Edges e1¹ and e2², e3¹ and e3², e5¹ and e6² and vertices v1¹ and v1², v2¹ and v3², v4¹ and v4² are common for both FEs.

The code RFS uses the next algorithm for building the global matrix:

1. Calculate a local matrix for the tetrahedron number one (q = 1). Give its edges global numbers from 1 to $N_c=6.$

2. Choose the tetrahedron with q = 2

   1. Choose the edge of this tetrahedron with m = 1. Check whether this edge coincides with any edge of the previous tetrahedron. If it does, calculate its matrix elements and add them to corresponding matrix elements of the previous tetrahedron. In the other case, give the edge global number $N_c+1$ and fill a corresponding element of the global matrix. Set $N_c=N_c+1$.

   2. Do this procedure for all edges of the tetrahedron.

3. Repeat point 2 for every tetrahedron in the region testing whether a given edge of the tetrahedron coincides with any edge, previously included in analysis.

Similar assembling procedure is valid for basis functions of higher order.

As a result, we get a global matrix of the problem, a global vector of the right side of Eq. (26) and a global vector of unknowns. Dimension of the global matrix is equal to the total number of coefficients (degrees of freedom) in formulas (23).

5.4. Numerical approximation of boundary conditions

• Border between two dielectrics. A finite element mesh is built so that every tetrahedron lies inside one of the media. Hence, the common face of two neighbour tetrahedrons lies on the separation surface or approximate it in the case of curvilinear border. Due to basis function properties, tangential electric field of both tetrahedrons is equal on the face edges, consequently electric fields are equal on the whole face area. As a result, BCs (7) for electric field are satisfied automatically. As for magnetic field intensity, its continuity on the border between media with different permeabilities is not guaranteed.

• Electrical wall. According to Eq. (8), tangential electric field on PEC surface is equal to zero. Hence, three coefficients of Eq. (25) belonging to edges, lying on the surface, are equal to zero. Corresponding rows and columns of local matrix of the given tetrahedron have to be eliminated.

• Magnetic wall. According to Eq. (9), surface tangential magnetic field on the PMC is equal to zero. We can see from Eqs. (27)–(31) that only term $s_m$ has to be changed:

  $s_m=j{k}_0{\eta }_0{\displaystyle {\int }_S\left(H\times w_m\right)e_n}=0.$E39

  This expression is equal to zero as the mixed product under the integral has two vectors ($H,e_n$) with the same direction. Therefore, PMC BC does not change local matrix.

• Absorption boundary conditions. Formula for matrix elements $s_{mn}$ can be derived from Eq. (11):

  $s_{mn}=j{k}_0\sqrt{{\mu }_r/{\epsilon }_r}{\displaystyle {\int }_S(e_n\times w_m)\cdot }(e_n\times w_n)dS.$E40

  Expressions for the higher order ABC can be found in [10].

• Impedance surface. Impedance surface has finite surface conductivity ${\sigma }_s$. Such surface is a good approximation of the metal surface in microwave frequency band. Using Leontovich BC (10), we can derive expression

  $s_{mn}=-j{k}_0{\eta }_0{\sigma }_s{\displaystyle {\int }_S{w}_m\cdot w_{n}dS}.$E41

  These values must be added to the corresponding elements of the local matrix. Of course, they are nonzero only for functions associated with edges lying on the impedance surface.

• Excitation sources. Electromagnetic field in the system excites by electric current or by electric or magnetic fields, defined on a part of the system border. In the first case, we use Eq. (26) to calculate right‐hand vector elements. In the second case, we define the so‐called ports—surfaces on which nonzero excitation field exists. Frequently, we define port as cross section of the regular transmission line (TL). Such ports are called wave ports. If TL cross section has simple form (rectangular, circular, etc.), we can define field on the port analytically. Otherwise, we have to solve two‐dimensional (2D) boundary problem.

  Together with wave ports, lumped ports (LPs) are often used. These ports are defined by current $I$ and intrinsic admittance $Z_i$. LPs are used when excitation source dimensions are very low compared to the whole‐system dimensions or wavelength. Geometrical model of LP is a segment of a straight line (linear LP) or a rectangle (planar LP). Linear and planar LPs are implemented in the RFS code. Arbitrary number of tetrahedron's edges or (and) faces can cover such LPs. We suppose that surface current density and electric field distributions on the PL are homogeneous. Values of right‐hand vector for a planar LP are calculated by Eq. (31), where ${\dot{J}}^{imp}=I/w$, $I$ is port current and $w$ is port width. Implemented model of the LP excludes limitations imposed on the mesh generator by other models.

• Lumped elements. Embedding of lumped elements (LEs) into 3D field model enlarges the capability of the code, as an LE excludes the necessity of generating very fine mesh inside such parts as resistors, capacitors and inductors. In the RFS code, linear and planar LEs are implemented [17].

  The presence of LEs modifies the global matrix by adding to diagonal matrix elements terms

  ${q^{\prime }}_{mm}=-j{k}_0{\eta }_{0}Yl,$E42

  where $Y$ is the admittance of the LE, $l$ is its length. The RFS code implements models of a resistor, capacitor and inductor, and their parallel connection. The number of LEs in the model is unlimited.

5.5. Fast‐frequency sweep

The finite element method in frequency domain solves the problem at one specific frequency. Calculation of the system frequency response (e.g. S‐matrix) needs solving the problem many times (often several hundreds or even thousands). This procedure takes much computation time.

Considerable economy of computational resources can be achieved by implementing the so‐called fast‐frequency sweep (FFS). The algorithm of FFS is based on model‐order reduction (MOR) technique. According to the MOR, we seek solution in one frequency point and then expand left and right sides of Eq. (21) in frequency power series. Terms with the same power are equated and full solution as a function of frequency is restored. The RFS code implements one of such algorithms—well conditioned asymptotic wave form evaluation (WCAWE) [18]. This method was adapted to the RFS procedures, particularly for LE models, and showed high efficiency, making possible to calculate system fields and parameters in three‐octave frequency band using full solution only in one frequency point. Hence, computation time for wideband systems decreases by tens and hundreds times without losing accuracy.

6. The RFS code description

The RFS code is designated to solve complex 3D electromagnetic problems, especially mobile handsets modelling. It contains graphic user interface for creating and editing geometric objects or for import CAD models. Geometric primitives, besides boxes, cylinders, spheres, pyramids and cones, include coaxial and strip lines and wizards for creating human head and hand models (phantoms), spiral antennas and other objects. All Boolean operations, such as extrusion, skinning and others, are also implemented. The code includes a rich material library.

Two types of meshing are available: exact automatic meshing based on Delaunay tessellation and hand meshing with automatic geometry corrections. We recommend the last one for meshing complex models, imported from CAD codes, for example, for full models of a handset. Two types of basis functions can be used: low order (CT/LN) and high order (LT/QN). The last is used as default, but for complex models with many small details, we recommend low‐order basis.

The code solves several types of problems: driven solution in one frequency point, driven solution in frequency band with given frequency step (solution in each point or FFS), eigenmode, EMS and EMI problems. The code can perform parametric solution when one or several geometric or (and) material parameters changes in a given order. The results can be represented as functions of these parameters.

Post‐processing includes the calculation of S‐, Y‐ and Z‐parameters, field distribution in a volume, on a surface or along a line, the calculation of cavities intrinsic impedance, specific absorption rate (SAR) and EMC/EMI parameters. Results are represented in one‐dimensional (1D), 2D and 3D graphs, tables and can be exported in a file.

7. Computation results

Firstly, we validated the code by solving problems having analytical solution. One of such problems is an eigenvalue problem for cylindrical cavity. A cylindrical cavity has curvilinear surface, so calculation can demonstrate the quality of surface approximation by a tetrahedron mesh. The analysed cavity had radius 10 mm and height 15 mm. Cavity walls were supposed to be perfectly conducting. Eight eigen frequencies and eigen fields were calculated using high‐order basis functions (LT/QN).

Computation results together with analytical data are shown in Table 1. The table also contains relative calculation error $\delta$. Azimuthal inhomogeneous modes are twice degenerated. Numerical errors remove this degeneration, so two eigen frequencies for each mode are given. Degenerated modes differ by their field, turned by 90° in azimuth direction and computation results confirm this phenomenon. Difference between degenerated eigen frequencies is no more than hundreds of per cent.

Mode	Analytical eigen frequency (EF), GHz	Mesh size
		6768		15,026		26,580
		EF, GHz	δ, %	EF, GHz	δ, %	EF, GHz	δ, %
TM010	11.474	11.502	0.235	11.4905	0.135	11.4848	0.085
TE111	13.302	13.317; 13.319	0.11; 0.135	13.3132; 13.3139	0.084; 0.089	13.3104; 13.311	0.063; 0.067
TM011	15.216	15.232	0.098	15.2287	0.077	15.2241	0.052
TE211	17.670	17.696; 17.705	0.147; 0.198	17.6853; 17.6877	0.086; 0.102	17.6802; 17.6811	0.058; 0.063
TM110	18.284	18.3246; 18.3287	0.222; 0.244	18.3081; 18.3102	0.131; 0.143	18.299; 18.291	0.082; 0.083

Table 1.

Comparison of analytical and computational eigen frequencies for a cylindrical cavity.

The table demonstrates that eigen frequency calculation error decreases with mesh refinement and not exceeds 0.1% on the last mesh.

Another example, an electric dipole, consisted of two perfectly conducting cylinders, connected to the linear‐lumped source (Figure 6). To solve this outer problem, we ought to constrain computational region by the so‐called air box with ABC on its surface. The code automatically determines air box size so that the minimal distance between the dipole and the air box is not less than $\lambda /4$, where $\lambda$ is wavelength in free space.

Figure 6.

A dipole model: a, dipole structure; b, mesh; c, 3D directivity chart; d, 2D directivity chart in the plane $\varphi =0$.

Cylinders have diameter 0.5 mm and length 9.5 mm. Together with a lumped port of 1‐mm length, the total dipole length is 20 mm. Excitation frequency was chosen to be 1.5 GHz, so the ratio $l/\lambda =0.1$. Such a dipole can be considered as elemental (Figure 6a). Its directivity is described by the formula $D(\theta ,\varphi )=1.5si{n}^2\theta .$

Figure 6b shows mesh used in simulation. The sphere inside cylinder contains refined local mesh, surrounding the dipole. The total number of tetrahedrons is 64553. The 3D directivity chart is shown in Figure 6c, and 2D directivity chart in plain $\varphi =0$ in Figure 6d. Calculated dipole directivity is 1.506 with relative error of 0.4 %. This error can be caused by the finite diameter of cylinders.

Figure 7a shows handset together with human head and hand phantoms. The handset has PIFA‐type antenna. Reflection coefficient of this antenna in the presence of the phantoms was calculated by RFS code and code SEMCAD [19], based on FDTD method.

Figure 7.

A handset with human head and hand phantoms: a, model; b, reflection coefficient chart.

Finite element mesh built by the RFS code in the phantom contains only 200,000 tetrahedrons. The total mesh size was about 450,000 tetrahedrons. Efficient algorithm of SAR calculation, implemented in the RFS code, provides using such comparatively small mesh size without the loss of accuracy. SEMCAD hexagonal mesh had about 1.5 billion of cells.

Figure 7b shows the module of reflection coefficient versus frequency for FRS and SEMCAD. As can be seen, both codes give similar results, especially for resonant frequencies.

We also calculated maximum SAR level in a brain tissue, averaged on 1 g of the tissue. The results are shown in Table 2. As can be seen, both codes give similar results, but SEMCAD solution time was nearly five times greater than RFS.

Frequency, MHz	SAR, W/kg
	RFS	SEMCAD
900	1.08	1.06
2000	1.18	0.99

Table 2.

SAR, averaged over 1 g of brain tissue

8. Conclusion

The problem of high‐frequency electromagnetic field simulation is formulated. A short survey of numerical methods for solving high‐frequency electromagnetic problems is presented. It was shown that one of the most efficient methods for solving inner EMP and outer EMP with moderate electrical size is the vector finite element method in frequency domain. The algorithm of this method, including mesh generation, building of the global matrix, boundary conditions approximation, SLAE solving and FFS technique was described. Some peculiarities of lumped ports and elements implementation into the RFS computer code were also depicted. Simulation results for a simple systems having analytical solution, as well as for complex problems, such as handset antenna analysis near human head show high accuracy and efficiency of the code. The code can be used for the solution of antenna problems, waveguide problems, PCB analysis and microwave resonator problems.

Acknowledgments

The author expresses his deep gratitude for valuable help provided by R.V. Salimov and R.I. Tikhonov from the Russian R&D Centre LG Electronics Inc.