A hybridizable discontinuous Galerkin method for electromagnetics with a view on subsurface applications

TL;DR

This paper introduces two hybridizable discontinuous Galerkin (HDG) schemes for solving Maxwell's equations in the time domain, focusing on dielectric and conductive media with applications in subsurface electromagnetics.

Contribution

The paper develops two novel HDG methods for Maxwell's equations, incorporating conduction effects and electric field hybridization, with detailed convergence and validation studies.

Findings

Both methods achieve optimal convergence rates.

Numerical results validate the accuracy in dielectric and conductive media.

The schemes are suitable for subsurface electromagnetic applications.

Abstract

Two Hybridizable Discontinuous Galerkin (HDG) schemes for the solution of Maxwell's equations in the time domain are presented. The first method is based on an electromagnetic diffusion equation, while the second is based on Faraday's and Maxwell--Amp\`ere's laws. Both formulations include the diffusive term depending on the conductivity of the medium. The three-dimensional formulation of the electromagnetic diffusion equation in the framework of HDG methods, the introduction of the conduction current term and the choice of the electric field as hybrid variable in a mixed formulation are the key points of the current study. Numerical results are provided for validation purposes and convergence studies of spatial and temporal discretizations are carried out. The test cases include both simulation in dielectric and conductive media.

Tables

Table 1. Table 1: Error table of the spatial convergence study.

	P1		P2
$1/h$	Error	Order	Error	Order
$5$	$2.63e-2$	$-$	$2.05e-3$	$-$
$10$	$6.42e-3$	$2.04$	$2.67e-4$	$2.94$
$15$	$2.84e-3$	$2.01$	$7.96e-5$	$2.98$
$20$	$1.60e-3$	$2.01$	$3.37e-5$	$2.99$
$25$	$1.02e-3$	$2.00$	$1.73e-5$	$2.99$
$30$	$7.08e-4$	$2.00$	$1.01e-5$	$2.96$

Table 2. Table 2: Error table of the time convergence study.

$1/\Delta t$	Error	Order
$25$	$3.61e-2$	$-$
$50$	$1.89e-2$	$0.93$
$100$	$9.67e-3$	$0.97$
$200$	$4.89e-3$	$0.98$
$400$	$2.46e-3$	$0.99$
$800$	$1.24e-3$	$0.99$

Table 3. Table 3: Error table of the spatial convergence study of a cubic cavity with conductive media.

	P1
$1/h$	Error	Order
$5$	$3.34e-2$	$-$
$10$	$8.20e-3$	$2.03$
$15$	$4.40e-3$	$1.84$
$20$	$2.58e-3$	$1.85$

Equations

$$\nabla \times \mathbf{E}$$

$$\nabla \times \mathbf{H}$$

$$\nabla \cdot \mathbf{D}$$

$$\nabla \cdot \mathbf{B}$$

$$\nabla \cdot \mathbf{i}=-\frac{\partial \rho }{\partial t}$$

$$\mathbf{D}$$

$$\mathbf{B}$$

$$\mathbf{J}$$

$$\mathbf{J}=\mathbf{\sigma }\mathbf{E}+{\mathbf{i}}_{\mathrm{s}}.$$

$$\mathbf{n}\times \mathbf{E}=0,$$

$$\sqrt{\epsilon \mu }\left(\frac{\partial \mathbf{E}}{\partial t}\times \mathbf{n}\right)\times \mathbf{n}-\nabla \times \mathbf{E}\times \mathbf{n}=0,$$

$$\mathbf{H}\times \mathbf{n}+\sqrt{\frac{\epsilon }{\mu }}\left(\mathbf{E}\times \mathbf{n}\right)\times \mathbf{n}={\mathbf{H}}^{\mathrm{inc}}\times \mathbf{n}+\sqrt{\frac{\epsilon }{\mu }}\left({\mathbf{E}}^{\mathrm{inc}}\times \mathbf{n}\right)\times \mathbf{n}={\mathbf{g}}^{\mathrm{inc}},$$

$${\mathbf{E}}_0$$

$${\mathbf{H}}_0$$

$$\epsilon \frac{{\partial }^2\mathbf{E}}{\partial t^2}+\mathbf{\sigma }\frac{\partial \mathbf{E}}{\partial t}+\frac{1}{\mu }\nabla \times \nabla \times \mathbf{E}$$

$$\epsilon \frac{{\partial }^2\mathbf{B}}{\partial t^2}+\mathbf{\sigma }\frac{\partial \mathbf{B}}{\partial t}+\frac{1}{\mu }\nabla \times \nabla \times \mathbf{B}$$

$$\epsilon \frac{{\partial }^2\mathbf{H}}{\partial t^2}+\mathbf{\sigma }\frac{\partial \mathbf{H}}{\partial t}+\frac{1}{\mu }\nabla \times \nabla \times \mathbf{H}=\frac{1}{\mu }\nabla \times {\mathbf{i}}_{\mathrm{s}}.$$

$$\mathbf{\sigma }\frac{\partial \mathbf{E}}{\partial t}+\frac{1}{\mu }\nabla \times \nabla \times \mathbf{E}=-\frac{\partial {\mathbf{i}}_{\mathrm{s}}}{\partial t}.$$

$$\mathbf{\sigma }\frac{\partial \mathbf{H}}{\partial t}+\frac{1}{\mu }\nabla \times \nabla \times \mathbf{H}=\frac{1}{\mu }\nabla \times {\mathbf{i}}_{\mathrm{s}}.$$

$$\frac{1}{\mu }\nabla \times \nabla \times \mathbf{E}=-\frac{\partial {\mathbf{i}}_{\mathrm{s}}}{\partial t},$$

$$\mu \frac{\partial \mathbf{H}}{\partial t}+\nabla \times \mathbf{E}$$

$$\epsilon \frac{\partial \mathbf{E}}{\partial t}+\mathbf{\sigma }\mathbf{E}-\nabla \times \mathbf{H}$$

$$\mu \frac{\partial \mathbf{H}}{\partial t}+\nabla \times \mathbf{E}$$

$$\mathbf{\sigma }\mathbf{E}-\nabla \times \mathbf{H}$$

$${\mathbf{V}}^h=\{{\mathbf{v}}^h\in [L^2({\mathcal{T}}_h){]}^d:{\mathbf{v}}^h{|}_{{\Omega }_e}\in {\mathcal{P}}_m^d({\Omega }_e)\forall {\Omega }_e\in {\mathcal{T}}_h\},$$

$${\mathbf{M}}_t^h=\{{\mathbf{\eta }}^h\in [L^2({\mathcal{E}}_h){]}^{d-1}:{\mathbf{\eta }}^h{|}_F\in {\mathcal{P}}_m^{d-1}({\mathcal{E}}_h),({\mathbf{\eta }}^h\cdot \mathbf{n}){|}_F=0\forall {\Omega }_e\in {\mathcal{T}}_h\},$$

$${\mathbf{M}}_t^h(\mathbf{f})=\{{\mathbf{\eta }}^h\in {\mathbf{M}}_t^h:{\mathbf{\eta }}^h\times \mathbf{n}=\Pi (\mathbf{f}\times \mathbf{n})\text{ on }{\Gamma }_D\}.$$

$${\left({\mathbf{v}}^h,\mu \frac{\partial {\mathbf{H}}^h}{\partial t}\right)}_{{\Omega }_e}+{\left({\mathbf{v}}^h,\nabla \times {\mathbf{E}}^h\right)}_{{\Omega }_e}$$

$${\left({\mathbf{w}}^h,\epsilon \frac{\partial {\mathbf{E}}^h}{\partial t}\right)}_{{\Omega }_e}+{\left({\mathbf{w}}^h,\mathbf{\sigma }{\mathbf{E}}^h\right)}_{{\Omega }_e}-{\left({\mathbf{w}}^h,\nabla \times {\mathbf{H}}^h\right)}_{{\Omega }_e}$$

$${\left({\mathbf{v}}^h,\mu \frac{\partial {\mathbf{H}}^h}{\partial t}\right)}_{{\Omega }_e}+{\left(\nabla \times {\mathbf{v}}^h,{\mathbf{E}}^h\right)}_{{\Omega }_e}+{⟨{\mathbf{v}}^h\times \mathbf{n},{\hat{\mathbf{E}}}^h⟩}_{\partial {\Omega }_e}=0,$$

$${\left({\mathbf{w}}^h,\epsilon \frac{\partial {\mathbf{E}}^h}{\partial t}\right)}_{{\Omega }_e}+{\left({\mathbf{w}}^h,\mathbf{\sigma }{\mathbf{E}}^h\right)}_{{\Omega }_e}-{\left(\nabla \times {\mathbf{w}}^h,{\mathbf{H}}^h\right)}_{{\Omega }_e}-{⟨{\mathbf{w}}^h\times \mathbf{n},{\hat{\mathbf{H}}}^h⟩}_{\partial {\Omega }_e}$$

$$=-{\left({\mathbf{w}}^h,{\mathbf{i}}_{\mathrm{s}}^h\right)}_{{\Omega }_e},$$

$${\left({\mathbf{v}}^h,\mu \frac{\partial {\mathbf{H}}^h}{\partial t}\right)}_{{\mathcal{T}}^h}+{\left(\nabla \times {\mathbf{v}}^h,{\mathbf{E}}^h\right)}_{{\mathcal{T}}^h}+{⟨{\mathbf{v}}^h\times \mathbf{n},{\mathbf{\Lambda }}^h⟩}_{\partial {\mathcal{T}}^h}=0,$$

$${\left({\mathbf{w}}^h,\epsilon \frac{\partial {\mathbf{E}}^h}{\partial t}\right)}_{{\mathcal{T}}^h}+{\left({\mathbf{w}}^h,\mathbf{\sigma }{\mathbf{E}}^h\right)}_{{\mathcal{T}}^h}-{\left(\nabla \times {\mathbf{w}}^h,{\mathbf{H}}^h\right)}_{{\mathcal{T}}^h}-{⟨{\mathbf{w}}^h\times \mathbf{n},{\hat{\mathbf{H}}}^h⟩}_{\partial {\mathcal{T}}^h}$$

$$=-{\left({\mathbf{w}}^h,{\mathbf{i}}_{\mathrm{s}}^h\right)}_{{\mathcal{T}}^h}.$$

$${\hat{\mathbf{H}}}^h={\mathbf{H}}^h+\tau \left({\mathbf{E}}_{\mathrm{t}}^h-{\mathbf{\Lambda }}^h\right)\times \mathbf{n},$$

$$[[\mathbf{n}\times {\hat{\mathbf{H}}}^h]]=0\mathrm{on}{\mathcal{E}}_h^0,$$

$$[[\mathbf{n}\times {\hat{\mathbf{H}}}^h]]:={\mathbf{n}}^{+}\times {\hat{\mathbf{H}}}^{h+}+{\mathbf{n}}^{-}\times {\hat{\mathbf{H}}}^{h-}$$

Full text

1 Introduction

The field of computational electromagnetics has been gaining interest since the first finite difference method for Maxwell’s equations was developed by Yee in 1966 [1]. With the advent of more powerful computers, it was possible to implement alternative methods such as finite element methods enabling to overcome some of the shortcomings of the original finite difference method, such as problems related to using unstructured meshes and difficulties in achieving higher-order approximations.

The applications of computational electromagnetics are vast and can be found in many engineering fields. Both direct and inverse problems were investigated in the past, e.g., [2, 3, 4, 5, 6], to cite just a few. Examples of fields that benefit from such methods are telecommunications, aeronautics, geophysics, medical imaging and semiconductor industries.

Finite elements are generally more computationally expensive than finite differences, but their ability to represent complex geometries and the possibility to locally refine meshes, either via mesh refinement or by increasing the polynomial order of the approximation, makes them very attractive for the present application. The development of such methods greatly benefited from the work of Nédélec on edge-based elements in [7, 8]. Among others, these elements enable a correct representation of the null space of the curl operator; see, e.g., [9].

In continuous Galerkin (CG) methods, the continuity between elements is enforced strongly, such that it is therefore not possible to represent discontinuities in the normal component of the electric (or magnetic) field at interfaces between different materials or at domain boundaries. An alternative way to solve this problem, retaining the geometric flexibility and high-order capability, is to use discontinuous Galerkin (DG) methods. In such methods, the continuity between elements is relaxed, such that it is possible to represent discontinuities at element interfaces. However, DG discretizations lead to larger systems of linear equations as compared to CG discretizations. However, to retain the good properties of DG methods, while addressing their drawbacks, hybridized discontinuous Galerkin (HDG) methods were proposed in [10]. Such methods are obtained from DG methods using flux definitions that reduce the number of globally coupled unknowns and make the solution of the equations faster, as a result. These methods perform similarly to CG methods for matrix-based (see, e.g., [11]) and matrix-free (see, e.g., [12]) implementations, if appropriate preconditioners are available. Our research focus on the simulation of electromagnetic diffusion problems in underground strata formations, which might considerably benefit from the advantages of HDG methods. The topic has been investigated using finite differences, e.g., in [13, 14] and finite elements, e.g., [15, 4, 16, 17, 18]. While HDG methods were developed for other electromagnetic applications in [5, 6], to the best of our knowledge, they have never been applied to the direct (or inverse) problem of electromagnetic diffusion in geological formations. We focus on an efficient solution of the forward problem, which plays a fundamental role also in inverse problems. To achieve our goal, we first derive and implement a complete HDG formulation for a mixed formulation and validate it by computing test cases which solutions are available in analytical or experimental form. We also derive a formulation based on an electromagnetic diffusion equation that is typically the equation of choice for geophysical applications.

We would like to point out the differences of our formulation to previous ones proposed in [5, 6]. In [5], a formulation based on the displacement current term was proposed, while in this paper the conduction current is considered. Moreover, while a frequency-domain approach is used in [5], our approach is formulated in the time domain. A crucial difference between the mixed formulation presented in this paper and the one in [6] is related to the definition of the hybrid variable. In [6], the hybrid variable is defined as the tangential component of the magnetic field, and an equation is introduced to enforce the electric field to be continuous, while in this paper the opposite choice is made.

The paper is organised as follow: Section 2 introduces the constitutive equations in their original form and in the form that will be used as a starting point for our formulation. In Section 3, our formulation is derived in a form that is general enough to allow the validation of the method with classical test cases and to obtain the electromagnetic diffusion equation required for subsequent research steps. Section 4 gives an overview of the chosen implicit time-integration scheme and its implementation. In Section 5, numerical examples for validating the implementation, including convergence studies, are presented. Finally, Section 6 provides conclusions of this work.

2 Governing equations

Classical electrodynamics is mathematically described by electric and magnetic field equations. In combined form, these equations are known as Maxwell’s equations. A derivation of the equations as well as a thorough discussion of their mathematical properties can be found, e.g., in [19].

2.1 Maxwell’s equations

Maxwell’s equations are given as follows:

[TABLE]

The variables that appear are: electric field $\mathbf{E}$ ($\mathrm{V/m}$), magnetic field $\mathbf{H}$ ($\mathrm{A/m}$), electric displacement field $\mathbf{D}$ ($\mathrm{C/m^{2}}$), magnetic induction field $\mathbf{B}$ ($\mathrm{T}$), electric current density $\mathbf{J}$ ($\mathrm{A/m^{2}}$), and electric charge density $\rho$ ($\mathrm{C/m^{3}}$). Equations (1a)–(1d) are Faraday’s law, Maxwell–Ampère’s law, (electric) Gauss’ law, and (magnetic) Gauss’ law, respectively.

The equation of continuity, stating conservation of charge, can be used to relate the current $\mathbf{i}$ to the charge density:

[TABLE]

To close the system, it is necessary to add the following constitutive equations, valid for isotropic materials, to equations (1a)–(1d):

[TABLE]

with the following constitutive parameters: permittivity $\varepsilon$ ($\mathrm{F/m}$), permeability $\mu$ ($\mathrm{H/m}$), and conductivity tensor $\bm{\sigma}$ ($\mathrm{S/m}$), which is the inverse of the resistivity tensor $\bm{\rho}_{\mathrm{e}}$ ($\Omega\mathrm{m}$). Equation (5), which is a generalized form of Ohm’s law, may be extended by adding a current-density source $\mathbf{i}_{\mathrm{s}}$:

[TABLE]

In this paper, two different boundary conditions are considered. The simplest boundary condition is a Dirichlet-type boundary condition on the tangential component of the electric field, known as metallic or perfect electric conductor (PEC) boundary conditions:

[TABLE]

where $\mathbf{n}$ denotes the outward-pointing unit normal vector to the surface.

Another, more complex set of boundary conditions are absorbing boundary conditions (ABC). These boundary conditions are used to truncate physically unbounded domains to finite numerical domains, reducing unwanted reflections. ABC are different in every field of application in which equations enable wave propagation. We derive the ABC for the present problem in the form of the Silver–Müller condition according to [2] as

[TABLE]

obtaining

[TABLE]

where $\mathbf{H}$ and $\mathbf{E}$ represent the total fields in the problem domain, while $\mathbf{H}^{\mathrm{inc}}$ and $\mathbf{E}^{\mathrm{inc}}$ represent external incoming fields.

Furthermore, initial conditions as

[TABLE]

are prescribed, such that the initial fields are solutions of the Maxwell system; see, e.g., [20]. To study the response to step-on sources, the initial fields are set to zero everywhere. In geophysical applications, however, where usually the response to step-off sources is of interest, the fields need to be initialized to represent the solution of the static Maxwell equations by solving direct current and magnetometric resistivity problems, as outlined, e.g., in [4, 14, 21]. A discussion of such initializations is beyond the scope of this paper, and therefore, in this work, either zero fields or given analytical solutions are used as initial conditions.

2.2 Wave equations for electric and magnetic field

From equations (1a) and (1b), using constitutive equations (3), (4), and (6), damped (vector) wave equations for conductive media can be derived for both the electric field $\mathbf{E}$ and the magnetic flux density $\mathbf{B}$:

[TABLE]

Alternatively, the second wave equation may be formulated in terms of the magnetic field $\mathbf{H}$:

[TABLE]

For subsurface applications, where typically electromagnetic (EM) signals of very low frequency (up to $f=10\,\mathrm{Hz}$) are used, the conductivity is usually much larger than the permittivity (e.g., nine orders of magnitude in rocks); see, e.g., [22]. As a result, the displacement current becomes negligible with respect to the conduction current. In this case, the damped wave equation reduces to an EM diffusion equation:

[TABLE]

The analogous equation for the magnetic field reads

[TABLE]

When the frequency goes to zero, the EM diffusion equation (15) reduces to the Laplace/Poisson equation

[TABLE]

which describes potential-field methods. In this limit, intrinsic resolution becomes almost nonexistent; see [22].

In Section 3.2, the EM diffusion equation (15), along with boundary conditions (7) and (9), as well as initial condition (10), will be considered, and a hybridizable discontinuous Galerkin method will be derived.

2.3 Mixed equations for electric and magnetic field

From equations (1a) and (1b), using the constitutive equations (3), (4), and (6), the following mixed equations for the electric field $\mathbf{E}$ and the magnetic field $\mathbf{H}$ may also be derived:

[TABLE]

Assuming the conductivity to be much larger than the permittivity, as above, the equations can be reduced to

[TABLE]

However, to keep the formulation more general, in Section 3.1, equations (18a) and (18b), along with boundary conditions (7) and (9) as well as initial conditions (10) and (11), will be considered, and a hybridizable discontinuous Galerkin method will be derived.

3 Hybridizable discontinuous Galerkin method for electromagnetic diffusion equation

3.1 Formulation based on mixed equations

The domain $\Omega$ is discretized with $n_{\mathrm{el}}$ elements $\Omega_{e}$ such that the collection of all elements, $\mathcal{T}_{h}$, partitions $\Omega$. The boundary of the domain is divided in two disjoint portions, such that $\Gamma_{A}\cup\Gamma_{D}=\partial\Omega$. In these two portions, different boundary conditions will be applied, i.e., ABC on $\Gamma_{A}$ and PEC boundary conditions on $\Gamma_{D}$. The boundary of the collection, $\partial\mathcal{T}_{h}$, consists of all element boundaries $\partial\Omega_{e}$, where $e=1,...,n_{\mathrm{el}}$. The union of all interior faces $F$, where an interior face denotes the intersection of the boundaries of two neighboring elements $i$ and $j$ as $F=\partial\Omega_{i}\cap\partial\Omega_{j}$, is denoted by $\mathcal{E}_{h}^{0}$. A boundary face is defined as $F=\partial\Omega_{i}\cap\partial\Omega$, where $i=1,...,n_{\mathrm{el}}$, and the union of such faces is denoted by $\mathcal{E}_{h}^{\partial}$. The union of interior and boundary faces is then given by $\mathcal{E}_{h}=\mathcal{E}_{h}^{\partial}\cup\mathcal{E}_{h}^{0}$.

The space of polynomials of degree of at most $m$ on $D$ is denoted by $\mathcal{P}_{m}(D)$, where $D$ is an open domain in $\mathbb{R}^{d}$, and $L^{2}(D)$ is the space of square integrable-functions on $D$. Furthermore, we set $\mathcal{P}_{m}^{d}(D)\equiv[\mathcal{P}_{m}(D)]^{d}$. The approximation spaces for solution and weighting functions are then defined as

[TABLE]

where $\Pi(\mathbf{f}\times\mathbf{n})$ is a projection of $\mathbf{f}\times\mathbf{n}$ onto $\bm{\mathbf{M}}_{t}^{h}$. Note that the space $\bm{\mathbf{M}}_{t}^{h}$ has a smaller dimension than $\bm{\mathbf{V}}^{h}$ because the functions in the former are defined on the faces of the elements, $F$, while the functions in the latter are defined in the elements, $\Omega_{e}$.

Equations (18a) and (18b) are now multiplied by (discrete) weighting functions $(\mathbf{v}^{h},\mathbf{w}^{h})\in\bm{\mathbf{V}}^{h}\times\bm{\mathbf{V}}^{h}$, respectively, and integrated over one element $\Omega_{e}$ with boundary $\partial\Omega_{e}$:

[TABLE]

where $(\mathbf{E}^{h},\mathbf{H}^{h})\in\bm{\mathbf{V}}^{h}\times\bm{\mathbf{V}}^{h}$ denote approximations of the fields $\mathbf{E}\text{ and }\mathbf{H}$ in the element $\Omega_{e}$ and $\mathbf{i}_{\mathrm{s}}^{h}$ is a projection of $\mathbf{i}_{\mathrm{s}}$ onto $\bm{\mathbf{V}}^{h}$.

Integration by parts yields

[TABLE]

where the numerical traces $\hat{\mathbf{E}}^{h}$ and $\hat{\mathbf{H}}^{h}$ appear as additional variables.

As the next step, the hybrid variable $\bm{\Lambda}^{h}:=\hat{\mathbf{E}}^{h}_{\mathrm{t}}$ is introduced, defined as the tangential component of the numerical trace $\hat{\mathbf{E}}^{h}$, such that $\bm{\Lambda}^{h}\in\bm{\mathbf{M}}_{t}^{h}(\mathbf{0})$. Note that the tangential component can be obtained as $\hat{\mathbf{E}}^{h}_{\mathrm{t}}=\mathbf{n}\times\left(\hat{\mathbf{E}}^{h}\times\mathbf{n}\right)$ or as $\hat{\mathbf{E}}^{h}_{\mathrm{t}}=\hat{\mathbf{E}}^{h}-\hat{\mathbf{E}}^{h}_{\mathrm{n}}$, where $\hat{\mathbf{E}}^{h}_{\mathrm{n}}=\mathbf{n}\left(\hat{\mathbf{E}}^{h}\cdot\mathbf{n}\right)$. This formualation is similar to the one in [6], with the difference that here the hybrid variable is defined as the tangential component of the electric instead of the magnetic field. The main motivation for this choice lies in the fact that the tangential component of the electric field is always continuous. While charge density, surface currents and material heterogeneity can induce discontinuities in any component of the magnetic field or in the normal component of the electric field, $\mathbf{E}_{t}$ is always continuous, as emphasized in [19] (Chapter 4–5). There is not any notable difference in dielectric media, where surface currents are not allowed, but there might be one in conductive media, such as the ones encountered in geophysical applications.

Summing all contributions by individual elements, the following equation is obtained:

[TABLE]

The numerical trace of the magnetic field $\hat{\mathbf{H}}^{h}$ is defined here as a sum of $\mathbf{H}^{h}$ and a stabilization term, weighting the difference between the tangential component of the discrete electrical field $\mathbf{E}^{h}$ and the hybrid variable $\bm{\Lambda}^{h}$:

[TABLE]

where $\tau$ is a stabilization parameter.

Via an additional equation, the continuity of the tangential component of $\hat{\mathbf{H}}^{h}$ across inter-element boundaries is enforced:

[TABLE]

with the definition

[TABLE]

for interior faces, or simply

[TABLE]

for boundary faces.

Definition (23) and condition (24), are required to close the system of equations and ensuring the numerical trace $\hat{\mathbf{H}}^{h}$ to be single-valued. As a result, the method is locally conservative. To prove that (24) ensures a single-valued numerical trace $\hat{\mathbf{H}}^{h}$, equation (23) is inserted in (24), yielding

[TABLE]

from which it is possible to obtain an equation for the hybrid variable:

[TABLE]

Substituting (26) in (23), the following equation is obtained

[TABLE]

eventually yielding

[TABLE]

Note that to obtain (28), the relation $\mathbf{n}^{+}=-\mathbf{n}^{-}$ is used.

For an individual element boundary $\partial\Omega_{e}$, the variational form of (24), with weighting function $\bm{\eta}^{h}$, is given as

[TABLE]

To include the ABCs defined by (9), equation (29) is integrated over the whole domain and the boundary condition terms are added as

[TABLE]

where $\mathbf{g}^{h}_{\mathrm{inc}}$ is a projection of $\mathbf{g}_{\mathrm{inc}}$ onto $\mathbf{M}_{t}^{h}$.

Adding (30) to equations (22) yields

[TABLE]

Applying integration by parts to (31) and rearranging the equations leads to the final problem: find $\left(\mathbf{H}^{h},\mathbf{E}^{h},\bm{\Lambda}^{h}\right)\in\bm{\mathbf{V}}^{h}\times\bm{\mathbf{V}}^{h}\times\bm{\mathbf{M}}^{h}(0)$ such that

[TABLE]

for all $\left(\mathbf{v}^{h},\mathbf{w}^{h},\bm{\eta}^{h}\right)\in\bm{\mathbf{V}}^{h}\times\bm{\mathbf{V}}^{h}\times\bm{\mathbf{M}}^{h}(0)$.

From these equations, the following matrix-vector formulation is obtained:

[TABLE]

3.2 Formulation based on equations for electric field

The derivation of the second formulation is sketched here, and the final form is given. The complete derivation closely follows the steps outlined above for the mixed equations.

A new variable $\mathbf{u}=\mu^{-1}\nabla\times\mathbf{E}$ is introduced into the EM diffusion equation (15), such that it can be split up into two equations as follows: of equations is obtained:

[TABLE]

with PEC boundary conditions (7). The new variable $\mathbf{u}$ is related to the magnetic field $\mathbf{H}$ via equations (1a) and (4) as follows:

[TABLE]

and vice versa,

[TABLE]

Note the differences of the present system of equations (34) to (18), wich was the starting point for the derivation in the preceding section.

Based on the discretization already introduced above, the equations are multiplied by (discrete) weighting functions $(\mathbf{v}^{h},\mathbf{w}^{h})\in\bm{\mathbf{V}}^{h}\times\bm{\mathbf{V}}^{h}$, respectively, and integrated over one element $\Omega_{e}$ with boundary $\partial\Omega_{e}$:

[TABLE]

where $(\mathbf{u}^{h},\mathbf{E}^{h})\in\bm{\mathbf{V}}^{h}\times\bm{\mathbf{V}}^{h}$ are the approximations of the fields $\mathbf{u}$ and $\mathbf{E}$, respectively, and $\mathbf{i}_{\mathrm{s}}^{h}$ is a projection of $\mathbf{i}_{\mathrm{s}}$ onto $\bm{\mathbf{V}}^{h}$.

Summing all contributions by individual elements, integrating by parts and substituting the definitions

[TABLE]

yields

[TABLE]

To close the system of equations the continuity equation

[TABLE]

is added.

The final problem reads: find $\left(\mathbf{u}^{h},\mathbf{E}^{h},\bm{\Lambda}^{h}\right)\in\mathbf{V}^{h}\times\mathbf{V}^{h}\times\mathbf{M}^{h}_{t}(0)$ such that

[TABLE]

for all $\left(\mathbf{v}^{h},\mathbf{w}^{h},\bm{\eta}^{h}\right)\in\mathbf{V}^{h}\times\mathbf{V}^{h}\times\mathbf{M}^{h}_{t}(0)$.

From this system of equations, the following matrix-vector formulation is obtained:

[TABLE]

Note that, even if the same notation is used, apart from $\mathbb{I}\text{ and }\mathbb{J}$, the matrices in (41) are different from the ones in (33). Note further the different structure of the matrices in that a zero matrix apperas in the upper left block for the non-time-dependent $(\mathbf{H},\mathbf{E})$ $2\times 2$ system in (33), while it is observed for the time-dependent $(\mathbf{\dot{U}},\mathbf{\dot{E}})$ $2\times 2$ system in (41).

4 Implicit time-integration

Due to the frequency content and high variability of the constitutive parameters of the target application, the equation is very stiff in time [3]. Therefore, in our implementation, an implicit time-integration is chosen, to avoid time-step restrictions of explicit time-integration schemes. In this section, the implicit time-integration for the formulation in Section 3.1 is shown, by using a backward-Euler scheme for demonstration purposes. Previous matrix and vector indices are neglected to simplify the notation. The implementation of a higher-order time-integration scheme would follow a similar rationale.

Applying the backward-Euler scheme to (33), the following matrix-vector system is obtained

[TABLE]

where $\Delta t$ is the chosen time-step and the superscript $n$ indicates that we are referring to the solution at the time $t_{0}+n\Delta t$.

Condensation of the upper left $2\times 2$ block yields a system of equations for the degrees of freedom of the hybrid variable exclusively:

[TABLE]

with

[TABLE]

and

[TABLE]

After having solved (43), the interior degrees of freedom of the magnetic and electric fields are computed elementwise as

[TABLE]

The time-discretized formulation according to Section 3.2 can be obtained in a similar way, but the derivation is here omitted for brevity.

5 Numerical examples

In this section, four numerical examples demonstrate that our method is able to reproduce physical phenomena such as wave propagation and scattering correctly as well as the possibility to utilize unstructured meshes and absorbing boundary conditions to bound numerical domains. The chosen implicit time-integration scheme makes the solution of a system of linear equations necessary for every time step of the simulation. There are several approaches that can be used to solve a system of linear equations, which efficiency depends upon the properties of the system itself. In this study, we use SuperLU_DIST [23, 24] within Amesos interface [25], which is sufficient for the size of the systems of linear equations considered in the following four examples. For larger problem sizes, iterative solvers will become necessary. However, the development of specific solvers and preconditioners, respectively, for the problems such as the present one will be the subject of future work.

5.1 Cubic cavity with perfect electric conductor boundary conditions

In this subsection, the solution of Maxwell’s equations in a cubic cavity with perfectly conductive boundary conditions is presented. The domain is a unit cube, discretized with hexahedral elements. The permittivity and permeability are set to $\varepsilon=\sqrt{3}\,\text{F/m}$ and $\mu=\sqrt{3}\,\text{H/m}$, respectively.

The frequency of the resonant wave is given by

[TABLE]

where $m,n$ and $p$ are integers that define the wave pattern in the domain, and $a,b$ and $d$ are the dimensions of the cavity. The period of the wave can then be obtained using the relation $T_{mnp}=1/f_{mnp}$, which is $T_{111}=2\,$s in this test case.

The analytical solution for this mode reads

[TABLE]

The intial conditions are obtained by (48), when setting $t=0$.

The error considered here is a relative error computed on the electric field in the $L_{2}$-norm, which is defined as

[TABLE]

5.1.1 Spatial convergence

To make sure that the two sources of errors resulting from spatial and temporal discretization do not interfere with each other, we study the spatial convergence rates for a particularly very small time step. The chosen time-integration scheme is only of first-order accuracy and therefore small time steps are anyway recomended. The time-step length chosen here is $\Delta t=10^{-4}\,$s.

Fig. 1 and Table 1 show that the expected spatial convergence rate of the order $p+1$, where $p$ is the order of the polynomial interpolation, is indeed achieved.

5.1.2 Temporal convergence

In this convergence study, the spatial discretization is assumed constant and fine enough such that the spatial discretization error is negligible. Second-order hexahedral elements, with a characteristic element size of $h=1/15$, are used. Fig. 2 shows the expected first-order convergence rate, for a backward-Euler time-integration scheme. In Table 2, the computed errors and convergence orders are given.

5.2 Wave propagation

In this section, the propagation of a plane wave in free space is considered. The simulation domain is a spherical portion of free space with radius $R=1.5\,$m. The domain is discretised with about $40000$ tetrahedral elements and ABCs as defined in (9) are used to bound the numerical domain. The total time for this simulation is $T_{\mathrm{max}}=3.33\times 10^{-8}\,$s, and the time-step length is $\Delta t=10^{-11}\,$s. The vacuum electromagnetic properties are $\varepsilon_{0}=8.854\times 10^{-12}\,\text{F/m}\text{ and }\mu_{0}=1.257\times 10^{-6}\,\text{H/m}$.

The incoming wave propagates in the $y$-direction and has a frequency of $f=300\,\text{MHz}$. The respective electric and magnetic fields are given as

[TABLE]

where $\omega=2\pi f$ is the angular frequency of the wave. The velocity of the wave is given by $c_{0}=1/\sqrt{\varepsilon_{0}\mu_{0}}$ and coincides with the physical speed of light in vacuum.

As in Section 5.1, it is possible to compute the error with respect to an analytical solution of the plane wave as given by (50). The results of this simulation can be used as a reference state for the scattering of a plane wave. The error in the $L_{2}$-norm of the electric field is $\|\mathbf{E}^{\mathrm{inc}}-\mathbf{E}^{h}\|_{\tau_{h}}=0.16814$ at $t=10^{-8}$s, corresponding to about three periods of oscillation of the wave. Note that $c_{0}=1/\sqrt{\mu_{0}\varepsilon_{0}}\approx 3\times 10^{8}$ m/s, such that in $10^{-8}$ s the wave travels $3$ m, that is, the length of the domain.

To quantify the errors due to the spatial temporal discretizations, a second simulation is run, using a time-step length of $\Delta t_{2}=0.5\Delta t$. The error is computed at the same time $\|\mathbf{E}^{\mathrm{inc}}-\mathbf{E}^{h}\|_{\tau_{h}}=0.15471$, indicating that the error is mainly due to the spatial approximation.

In Fig. 3, the real part of the Discrete Fourier Transform (DFT) of the solution is shown for the frequency $f=300\,\text{MHz}$, allowing for a comparison with results in the frequency domain. The Fourier transform is computed by a pointwise temporal transformation using Numpy’s FFT implementation (see, e.g., [26]) and normalized by $1/n$, where $n$ is the number of time-steps used in the DFT. The normalization is required for enabling a comparison of the numerical results to the analytical Fourier transform of the plane wave, which is given by

[TABLE]

where $f_{0}$ is the frequency of the wave and $\delta(f_{0}-f)$ the Dirac delta function. The Fourier transform of $E_{z}$ for the frequency $f_{0}=300$ MHz reads

[TABLE]

and it can be observed that the Fourier transform of the electric field is related to the electric field itself by the relation

[TABLE]

The numerical error in the $L_{2}$-norm of $\hat{E}_{z}^{h}$ is $\|\hat{E}^{\mathrm{inc}}_{z}-\hat{E}_{z}^{h}\|_{\tau_{h}}=0.2763$.

Note that, in the simulations, the incoming wave is partially reflected by the boundaries, reducing the accuracy of the solution especially in the region of the positive $y$-axis.

5.3 Wave scattering

In this test case, the scattering of a plane wave by means of a dielectric sphere is computed. The domain $\Omega_{1}$ is spherical, with radius $R=1.5\,\text{m}$, and contains a sphere $\Omega_{2}$ of radius $r=0.5\,\text{m}$ and dielectric material, as shown in Fig. 4.

On the boundary of $\Omega_{1}$, ABCs of Silver–Müller type with the components of the wave defined in (50) are prescribed. The spatial discretization is the same as in Section 5.2. The material properties in $\Omega_{1}$ are also the same as in 5.2, while in $\Omega_{2}$, it is assumed that $\varepsilon_{2}=2\,\varepsilon_{0}$ and $\mu_{2}=\mu_{0}$.

For this test case, two different simulations frequencies are considered. Fig. 5 shows the scattering of a wave oscillating with a frequency $f=300\text{MHz}$ on the left-hand side and $f=600\text{MHz}$ on the right-hand side. Again, For $f=300\text{MHz}$, the time-integration parameters are the same as in Section 5.2, that is $T_{\mathrm{max}}=3.33\times 10^{-8}\,$s and $\Delta t=10^{-11}\,$s. For $f=600\text{MHz}$, due to the higher frequency of the wave, both simulation time $T_{\mathrm{max}}$, and time-step length $\Delta t$ are halved, such that $T_{\mathrm{max}}=1.67\times 10^{-8}\,$s and $\Delta t=5\times 10^{-12}\,$s.

5.4 Cubic cavity with perfect electric conductive boundary conditions in conductive media

To validate the formulation for problems in conductive environments, such as those that can be found in geophysical applications, a numerical example based on a manufactured solution is used. In this test case, permittivity and permeability values are set to $\varepsilon=\sqrt{3}\,\text{F/m}$ and $\mu=\sqrt{3}\,\text{H/m}$, respectively, and a non-zero conductivity value is introduced: $\sigma=1.429$ S/m. This conductivity is inspired by values used in [4] as exemplary conductivities of seabed models.

The manufactured solution aims at enforcing the analytical solution (48) in the mixed equations (18), yielding

[TABLE]

as a constraint for the forcing term. As a result, if $\mathbf{i_{\mathrm{s}}}$ satisfies (54), the analytical solution (48) will satisfy equations (18).

The time discretization is the same as in Section 5.1.1, that is, $\Delta t=10^{-4}$ and $T_{max}=5\times 10^{-4}$. Different spatial discretizations are used to simulate the proposed problem. First, hexahedral elements within structured meshes are used. A spatial convergence study is carried out by discretising the domain with the meshes used in Section 5.1.1. The obtained error values and spatial convergence rates are identical to the ones shown in Table 1, indicating that the conductivity term is not critical in terms of the accuracy of the discretization scheme.

Additionally, the domain is dicretised using tetrahedral elements within an unstructured mesh, as exemplary depicted in Fig. 6.

A convergence study is carried out using increasingly finer meshes, and the results are reported in Table 3. As expected, the convergence rate only approximates the convergence rate expected for linear elements, due to the non uniform element size across the domain.

6 Conclusions

Two implicit Hybridizable Discontinuous Galerkin formulations for the solution of Maxwell’s equations in three-dimensional domains have been presented. The formulations are designed to be used in domains with diffusive mediums, and therefore, the conductivity term is present in both formulations. While the formulation for the electric field only takes into account the conduction current, the formulation for the mixed equations is more general in that it retains both displacement and conduction currents.

The formulation for the mixed equations has been validated for three well-known test cases, such as the resonance in a cubic cavity, wave propagation and scattering, and for a manufactured solution to particularly take into account the conductivity term. Spatial convergence studies have been carried out for the wave-resonance test case in dielectric media, and the expected convergence rates have been obtained. Solutions of wave propagation and scattering have been proposed in time and frequency domain, allowing for comparisons with different solution methods. The case of a resonant cavity with conductive media has also been simulated, and the accuracy of the results turned out to be identical to the test case with dielectric media.

The diffusive term in the equations will be particularly investigated in future work on applications in the field of geophysics, that is, for transient electromagnetics in three-dimensional diffusive earth media, as were already considered before, e.g., in [4, 16]. Future research will also focus on solvers and preconditioners for HDG methods for electromagnetics. Moreover, a comparison between implicit and explicit time-integration will be investigated. Numerical evidence shows that, even though explicit schemes can have severe time-step restrictions, explicit implementations may still be more efficient than implicit ones in the end; see, e.g., [27], where such a comparison was carried out for HDG methods in the context of problems of acoustic wave propogation, which exhibit similarities to the present problems.

Acknowledgements

The support by the German Federal Ministry for Economic Affairs and Energy via the Central SME Innovation Program “Zentrales Innovationsprogramm Mittelstand (ZIM)” (Grant No. ZF4415001PO7) is gratefully acknowledged.