Second Order Finite Volume Scheme on Tetrahedral Meshes for Three-Dimensional Maxwell’s Equations with Discontinuous Dielectric Permittivity

An extension of a finite volume scheme for three-dimensional Maxwell’s equations with discontinuous dielectric permittivity on tetrahedral meshes is discussed. The scheme is second order accurate in time and space for regions with linear and curvilinear discontinuities. It was tested on problems with linear and curvilinear discontinuity configurations. The test results support the second order accuracy of the proposed scheme.


Introduction
The time-dependent Maxwell's equations describe the propagation and diffraction of the electromagnetic waves.For many problems of interest analytical solutions are not available and as a result a number of numerical algorithms were developed.
Arguably the most widely used approach for numerical solution is the finite difference time domain method (FDTD).It was first suggested in [1] and since then became very popular.FDTD uses structured cartesian grids and is very simple.FDTD has second order of approximation in space and time for uniform media.Its main shortcoming is its inability to accurately represent curvilinear boundaries.This leads to a reduced order of approximation for such cases.Another approach that does not have this shortcoming is the finite volume time domain method (FVTD).It relies on structured non-cartesian grids and unstructured meshes.As a result it can represent curvilinear boundaries more precisely and can preserve the second order of approximation for such cases.
Of special interest are the problems involving spatially variable electromagnetic properties.The development of second order numerical algorithms for such problems is a challenging task.For FDTD various ways of electromagnetic properties smoothing were considered [2].This resulted in better accuracy, but not the second order of convergence yet.In [3] an algorithm was suggested that allowed to preserve the second order of approximation.For FVTD the use of continuous variables for linear discontinuities in two dimensions was suggested in [4] and extended to three dimensions and curvilinear discontinuities in [5].This approach increased the accuracy but did not preserve the second order of approximation.Recently in [6,7] a simple approach to the gradient calculation and limitation was suggested that allowed to preserve second order of approximation for the Maxwell's equations with discontinuous electromagnetic properties.In this paper we extend the scheme of [6,7] to three dimensions.The resulting scheme is tested using problems with linear and curvilinear dielectric permittivity discontinuities.For both cases the computational results support the second order of approximation.

Maxwell's equations
The system of non-stationary three-dimensional Maxwell's equations can be written using nondimensional variables in vector form as: where U is the vector of conservative variables, F 1 , F 2 and F 3 are the flux vectors111 In these formulas E is the electric field, H the magnetic field, D = εE the electric flux density, B = μH the magnetic flux density, ε the dielectric permittivity, and μ the magnetic permeability (assumed 1.0 in this paper).The system of non-stationary Maxwell's equations can also be written using the flux variables V related to the conservative variables by U = QV as where and the matrices A 1 , A 2 , A 3 are written as

Numerical Scheme
By integrating the system of the Maxwell's equations over a tetrahedral cell Δ i with faces Γ k assuming constant dielectric permittivity in the cell, we can obtain the integral conservation law where (n 1 , n 2 , n 3 ) is the unit normal directed outside.For the approximation of this equation consider a finite volume Godunov scheme [8] where Ω Δ i is the volume of the ith cell, s k Δ i the area of its kth face and the flux F is calculated using the exact solution [9] to the Riemann problem are approximations of V from two adjacent cells at the face center X Γ .The matrices A + and A − can be written using the notation Y = √ ε as , where Y L,R are values in two adjacent cells.This scheme will have the second order of approximation in space and time if the values at the face center are approximated with the second order.To obtain such values the interpolation is used where X L,R are centers of two adjacent cells, and the derivatives are first order accurate [10,11].
To approximate and limit the derivatives we use an extended two-stage procedure [6,7].First we calculate gradients using values in the reference cell and in its neighbors with the same dielectric permittivity by means of the least squares method [12].Then we limit gradients.Two approaches can be used here: continuous limitation and independent limitation.To perform continuous limitation, continuous variables at vertices on the discontinuity are used.Values in cells are extrapolated on adjacent vertices using computed derivatives to calculate limit values.Extrapolated values are averaged.At the vertices located on the dielectric permittivity discontinuity, extrapolated valuers are transformed to Mathematical Modeling and Computational Physics 2015 02030-p.3continuous variables before averaging.Limited gradients are calculated from averaged values at vertices, values at vertices showing dielectric permittivity discontinuity need to be transformed back from continuous variables.To perform independent limitation, two limit values are computed at the nodes characterized by dielectric permittivity discontinuity at the two sides.Values in cells are extrapolated on the adjacent vertices using computed derivatives to calculate limit values.Extrapolated values are partitioned into groups according to the value of the dielectric permittivity in the cell from which extrapolation took place.Each group is averaged independently.Limited gradients are calculated from averaged values at vertices.If the vertex is located on the dielectric permittivity discontinuity and has two average values, the value that corresponds to the dielectric permittivity in the cell is used to calculate limited gradients.
The scheme is explicit and is conditionally stable subject to the Courant-Friederichs-Lewy (CFL) condition.The information propagation velocity is inversely proportional to √ ε.Therefore the stable timestep is proportional to √ ε.

Computational results
To test the order of approximation of the proposed scheme, a range of test computations for a number of problems were performed.We used tetrahedral meshes built using Gmesh software [13].The accuracy was evaluated by comparing the numerical solution to the analytic one.The error was computed in L 2 norm.Here we present results only for one problem with curvilinear discontinuity and continuous limitation.Results for other cases were similar.

Curvilinear discontinuity
Consider the propagation of a hybrid electromagnetic mode in a waveguide with a step-like profile of the dielectric permittivity.The discontinuity of the dielectric permittivity ε (or index of refraction In this case the system of Maxwell's equations has an analytic solution [14] that in cylindric coordinates can be written in the form where J 0 and J 1 are Bessel functions of the first kind, K 0 and K 1 Bessel functions of the second kind, and β can be obtained from the dispersion relation.We used the following test constants ε 0 = 1.0, ε 1 = 2.25, k = 6.0, a = 0.64, β = 8.402440923258, w = 3.764648073438, u = 2.063837416842.
Mathematical Modeling and Computational Physics 2015 02030-p.5As a computational region we used a cylinder of radius 1.0 and height 0.6.The computations were performed using a sequence of meshes composed of 18845, 51147, 144048, 405279, and 1250790 tetrahedra.The mesh was built in such a way so that the discontinuity of the dielectric permittivity was along the cell faces.The time steps were chosen proportional to the linear cell sizes.Fig. 1 shows a tetrahedral mesh composed of 51147 tetrahedra.Fig. 2 shows the distribution of the first component of the magnetic field H 1 at the time T = 1.49 for the cross-section x 3 = 0.3 obtained using a mesh of 144048 tetrahedra.The error evolution for a sequence of five meshes is presented in Fig. 3. Table 1 shows the maximum errors in the L 2 norm.The error behavior corresponds to the second order of approximation.

Conclusion
A finite volume scheme for the numerical solution of three-dimensional Maxwell's equations with dielectric permittivity on tetrahedral meshes was suggested.The scheme is second order accurate in space and time.A range of numerical tests for linear as well as curvilinear discontinuities confirm the second order of approximation.

Table 1 .
Maximum error