New high order FDTD method to solve EMC problem

In electromagnetic compatibility (EMC) context, we are interested in developing new accurate methods to solve efﬁciently and accurately Maxwell’s equations in the time domain. Indeed, usual methods such as FDTD or FVTD present important dissipative and/or dispersive errors which prevent to obtain a good numerical approximation of the physical solution for a given industrial scene unless we use a mesh with a very small cell size. To avoid this problem, schemes like the Discontinuous Galerkin (DG) method, based on higher order spatial approximations, have been introduced and studied on unstructured meshes. However the cost of this kind of method can become prohibitive according to the mesh used. In this paper, we ﬁrst present a higher order spatial approximation method on cartesian meshes. It is based on a ﬁnite element approach and recovers at the order 1 the well-known Yee’s scheme. Next, to deal with EMC problem, a non-oriented thin wire formalism is proposed for this method. Fi-nally, several examples are given to present the beneﬁts of this new method by comparison with both Yee’s scheme and DG approaches.


Introduction
In solving electromagnetism problems, Yee's scheme (FDTD) [1] is one of the most widely used numerical method [2]. However this method is quite a low order one (based on a second order spatial approximation) and presents numerical dispersive errors. The later may lead to erroneous solutions for some configurations. For instance, the attenuation of an electromagnetic field inside a cavity requests long time range computations and, in this case, the accumulation of dispersive errors degrades significantly the solution computed. In order to avoid this problem the main idea is to increase the spatial order of approximation of the numerical method. Several alternative schemes have been proposed in that way. One of the most popular is the Discontinuous Galerkin (DG) approach [7,8,9]. By using this kind of method, it has been proved that the dispersive numerical error is reduced. Unfortunately, DG schemes can present an important cost of computation. An alternative is thus to look for FDTD-like methods (preserving low computational costs) but with high order (to lower dispersion error). The two main ways were thus to develop finite differences on bigger stencils [3] [4], or introduce polynomial basis functions in cells [5] [6].
In this paper we propose a new high order spatial approximation method on cartesian grids which recovers FDTD at the first order but which presents, for a similar spatial order, lower computational costs than DG methods. Main specificity of our method lies into the particular choice of basis functions used. This method is shown to be able to obtain accurate solutions on inhomogenous meshes (ie with spatial step sizes varying in each direction). Indeed, the method authorizes the use of different spatial order of approximation for each cell in a given direction. Finally, to address EMC problems, a thin wire model is necessary to take into account the presence of cables. We thus propose in this paper an oblique thin wire formalism adapted to our high order spatial approximation method. This paper is divided as follows. In the first part we give a description of the principle of the method, the stability criterion, and we end with some comparisons between our scheme and both FDTD and a DG approach. Then, in the second part of the paper, we present the introduction of the thin wire model in the high order scheme and propose numerical validations to validate it on several examples.

Higher spatial order FDTD formulation
In this section we present a high order method based on previous works by G.Cohen [10]. We detail the variational formulation of the problem and the approximation space used to construct the scheme. Then, we give the CFL condition and finally we propose numerical examples showing the interest of the method.

Mathematical formulation
Let Ω be a bounded domain in IR 3 and T ∈ IR + a quantity similar to an observed time. The studied physical problem consists in finding E(t, x) and H(t, x), for x ∈ Ω and t ∈ [0, T ], such that: where n stands for the outgoing normal on the boundary of Ω noted ∂Ω, ε and µ are respectively the permittivity and the permeability, σ is the conductivity and J is a source term.

Construction of the scheme
Let τ h be a cartesian mesh of Ω such that Ω = ∪ K∈τ h K. Any element K of τ h is a parallelepiped cell. Let Q r be the set of polynomials from IR 3 into IR such that each term has a degree lower or equal to r. We introduce the approximation spaces of H 0 (curl, Ω) and L 2 (Ω), respectively by U and V with: Then, we define two sets of vector basis functions ϕ ν l1l2l3,K ∈ U and ψ ν l1l2l3,K ∈ V , for ν ∈ {x, y, z}, indexed on (l 1 , l 2 , l 3 ) and any K ∈ τ h , by where ( e x , e y , e z ) is the orthonormal basis in IR 3 , (x,ŷ,ẑ) belongs to the unit elementK = [0, 1] 3 , F K is a Q 1 bijective mapping betweenK and K, L g li and L gl li are Lagrange polynomials defined on [0, 1] respectively by using r Gauss ((ξ g k ) k=1...r ) and r+1 Gauss-Lobatto ((ξ gl k ) k=1...r+1 ) points. These points define respectively the quadrature schemes (4) and (5): Then, by using these basis functions, we approximate E(t, x, y, z) by a local expansion A similar approximation is obtained for H replacing basis functions ϕ l1l2l3,K with ψ l1l2l3,K ones and making correspond the index ranges. To clarify the spatial approximation chosen, figure 1 shows the location of the electric and magnetic unknowns in a cell K for some components. We note on this figure that these degrees of freedom are similar, for a spatial order equal to one at the locations of the unknowns in FDTD. Hence, we obtain a "rise in spatial order" with this method of the classical FDTD method. For this reason, we call this approach a "higher order FDTD method". According to figure 1, the first Gauss-Lobatto degree of freedom of a given cell should be coincident with the last one of the previous cell. To grant the continuity of the tangential components of the fields, we force the corresponding degrees of freedom to be the same. Hence, for instance, if K and K are adjacent cells on the y boundary (K on the left of this boundary and K on the right), we have Source term J is discritized spatially according to the same basis functions as E.
Reporting these functions into (3), we now look for E = (E x , E y , E z ) ∈ U and H = (H x , H y , H z ) ∈ V verifying: ∀ν ∈ {x, y, z}, All integrations in (6)-(7) are performed with quadrature methods. By considering the orthogonality properties of the basis functions and with a smart choice for the quadrature integration formula, we can obtain a numerical scheme with diagonal mass matrices and very sparse stiffness matrices. Actually, in (6) when the variable x of a given unknown function E ν or H ν is associated to a Gauss basis function L g of order r we use (4) in the x direction. Conversely, if E ν or H ν involves a Gauss-Lobatto function L gl in x then we approximate the integration in x with (5). Same philosophy is applied to variables y and z. To illustrate the process, the x component of (6) is evaluated with: .
where ε K and σ K are the values of ε and σ restricted to the cell K, suppϕ x l1l2l3 denotes the support of the function ϕ x l1l2l3 , and we divided by the volume |K| of the cell K whose lengths according to x, y and z are respectively dx K , dy K and dz K . For equation (7), we use for all the terms the same rule: if the x (resp. y or z) part of ψ ν is L g then we use (4), else we apply (5). For instance, x component of (7) becomes: where µ K is the value of µ on the cell K.
Finally, discretization in time is led with a leap-frog scheme, as done for the FDTD. Noting E n and H n+ 1 2 the vectors of all unknowns respectively in E at time n∆t and H at (n + 1 2 )∆t, we obtain the numerical method given by: where M E , M σ , M J and M H are diagonal matrices, and R E and R H are very sparse stiffness matrices which do not need storage in memory. Hence, we obtain for any spatial order approximation a fast, accurate and not memory consuming method.

Stability and convergence
To show the stability of the scheme (8), we follow the demonstration developped in [9]. So, we study the conservation in time of the quantity n h = Ω εE n · E n + Ω µH n+1/2 · H n−1/2 , and we can show the positivity n h , and hence the stability of the method, under the CFL condition: where λ max (A) is the highest eigenvalue of a given matrix A,M andR are respectively the mass and stiffness matrices evaluated on the reference elementK, and Consistency of the method have been determined numerically. To do so we considered the propagation of a mode inside a cavity and compared the solutions obtained with the analytical one. The study was led with different spatial orders varying from 2 to 4 and two meshes with cell sizes equal to h and h/4. Time step was set much smaller than the value prescribed by (9) to neglect the error due to the time approximation. Thus, we determined a convergence rate in O(h r+1 ) for our method. We note that this value is equal, in the case of a spatial order taken to 1, to FDTD one.

numerical validation.
To validate and quantify the advantages of our method we consider two examples. Results are compared with those obtained by two other numerical schemes: a DG method [9] and usual FDTD [1]. In all the sequel, the notation DIFOE Or will define our method used with a polynomial order r.

Transverse electric mode inside a cavity
For the first example, we compute the evolution of a T E mnk mode, with (m = 3, n = 3, k = 0), inside a perfectly metallic cubic cavity with a side of 1 meter length. Figure 2 shows a comparison between the analytical solution and FDTD for various cell sizes and spatial order. Corresponding computational costs are summarized in We can notice on the figure 2 that FDTD scheme is not accurate unless we use a mesh with cell size lower or equal to λ/50. Compared to this last result, we obtained a better solution in our method when using a spatial order of 5 but for a mesh with cells enlarged up to λ/2. This traduces into an important gain in CPU time as can be observed in the table 1. Now, we consider the same example solved with DG method on the same meshes and with the same orders as used for our higher spatial order FDTD method. Results are presented on figure 3 and show a similar accuracy between these two high order methods. However, computational costs for the DG approach are 3 times greater than those of our new FDTD method. This example shows once again the interest of our method.

Simulation of a building structure
We are now interested in a more complex and realistic test case given by the building structure presented on figure 4 illuminated by a plane wave. To take into account for this kind of complex geometries we use non-homogeneous meshes with a strategy of spatial steps variable per direction and various  approximation orders per component and per cell. Thus, for this example the mesh used is made with space step varying from 5cm up to 30cm. In the same way, the spatial approxi-  Numerical results for a point inside the building are presented on figure 5 and compared to FDTD simulations on a uniform mesh with step size equal to 2.5cm. Note that, in this test case, FDTD on coarser meshes did not seemed to be satisfactory. In this condition, the CPU time required to perform FDTD simulation is twice the one of our method.

Oblique Thin wire model
To address EMC problems, an oblique thin wire model suited to our higher order FDTD method is required. This section is thus devoted to propose a such formalism and validate it on some examples.

Mathematical aspects of the wire model
Following works of EldelviK [11] and Guiffaut [12], we obtain our formalism from the model proposed by Holland [13] adapting it to our higher order FDTD approach. First, we consider that wires are described by a set of segments. According to [13], on these segments currents I and charges q are supported and given by: where − → u s is the unit vector along the considered segment of wire s, and c 0 the speed of light in vacuum. L is an averaged inductance per length value obtained, as for the Holland formalism, by: where a stands for the radius of the wire and R is a length equal to the half distance between two unknowns of the discretization of I along the wire. The choice for this last quantity has been obtained by numerical experiments. Figure 6: Example of a segment going from Node(j,1) to Node(j,2) crossing two cells and divided in two parts.
The key point to use (10) is the evaluation of the source term L −1 E · − → u s on a given segment s of the wire. To obtain it, we first split the segment into several parts such that a part belongs to only one cell at a time. This step is pictured on figure 6. The set of cells K ∈ τ h intersting the segment s is noted Λ(s), and the part of s belonging to any K ∈ Λ(s) is noted s K . Hence we have s = ∪ K∈Λ(s) s K . This leads to the following evaluation of L −1 E · − → u s : Then, we evaluate each integral over s K by using the Gauss quadrature formula (4) and we obtain: where X i s K is the coordinate in IR 3 of the point corresponding to the i-th quadrature point on s K . The number of quadrature points r 0 is chosen big enough to evaluate accurately the integral. Finally, by reporting the approximation of the electric field, we can write the coupling term between field and wire (11) under the form: where coefficients are given by: ∀ν ∈ {x, y, z} Last part of the wire model is to introduce the action of the wire on electromagnetic fields. This is done in the Ampere equation by addition of a source term J [13]. Thus, this is performed in (6) by evaluating the term Ω J · ϕ ν l1l2l3,K for any K belonging to Λ(s). We assume that the wire is a cylinder with radius R around the segment s and that the current I on s, noted I s , is homogeneous inside this cylinder. Introducing the same quadrature formula as done in (11), we thus obtain on each K ∈ Λ(s) ∀ν ∈ {x, y, z}, whith C R the disk of radius R, = s K I s u s · ϕ ν l1l2l3,K by using C R Jds = I s u s . Finally, by using quadrature formula , we obtain for the previous integral: Note that the coupling coefficients C ν l1l2l3,K involved in (11) and (13) are identical. This particularity is very important as it permits to prove the stability of the coupled system Mawell -wire when the discretizations of (3) and (10) are separately led with stable schemes.

Numerical validations
To validate our model, we investigate in this section three typical examples of coupling between electromagnetic fields and wires: a straight line, a circular loop and a double loop. Obtained results are compared to FDTD with Holland model and the method of moment (MoM).

Straight line
The first example consists in the evaluation of the current on a straight wire illuminated by a plane wave. The direction of propagation for the plane wave is chosen perpendicular to the axis of wire and the electric component field is set tangential to the wire. The length of the wire is 0.525m and its radius is equal to 1mm. Simulations with our higher FDTD method are performed with a spatial approximation order equal to 3.  Comparisons have been done with the FDTD scheme using Holland formalism for wires [13] on two configurations. First, we consider the wire in line with the x axis of the grid. Then, the grid is turned with an angle of 26.56 degrees around the y axis. These configurations and corresponding mesh cells intersecting the wire are pictured on figure 7.
As parameters of the plane wave are only depending on the geometry of the wire, resulting solutions in these two configurations should be identical. Hence, they are proposed to both validate the coupling Maxwell -wire and ability of the scheme to deal with oblique wires. Note that FDTD computations are led on the first case and taken as a reference for the comparisons. Results obtained for a point in the vicinity of the wire are plotted on the figure 8. One can notice a good accordance of the results between our method and FDTD which validates our wire formalism. Moreover, our coupled method seems to deal efficiently with oblique wires, this last point usually requiring specific development to obtain acurate results in FDTD [12].

Circular loop in reception mode
The second example consists in the computation of the current at two points P 1 and P 2 on a circular loop. The loop is constituted by a thin wire with a radius equal to 1mm and is illuminated by a plane wave. We consider two configurations for this test case. The first one assumes that the axis of the spire is perpendicular to the y axis. In the second one, the axis of the spire is rotated with an angle of 30 degrees. These configurations are pictured on figure 9.
Once again, we chose to treat a configuration lined up with the mesh directions and another one turned to test the efficiency of our method to address complex wire configurations. On figure 10 we propose the comparison in the frequency domain of the results obtained between our method and Method of Moments(MoM). To treat this example, the spatial order of approximation have been fixed to 2 as the mesh used is not coarse. On the figure, we can observe the good agreement between the solutions obtained with the two methods. This confirms the validity of our wire model.

Double loop in reception mode
The last example we treat is a double loop which contains three wires connections. It is illuminated by a plane wave. The device is represented on figure 11). Figure 12 shows the solutions obtained with the three methods. It shows once again a good agreement between them. In particular, this example shows the ability of our model to take into account correctly the connections of wires.

Conclusion
In this paper, we have presented a time domain high order cartesian finite elements method which can be considered as a higher spatial order FDTD method. We have shown that this method is very efficient for problems having a long time of observation like for instance in cavity problems. More-23 (a) Loop with axis perpendicular to y axis.
(b) Loop rotated with 30 degrees. Figure 9: The two configuration studied for the circular loop over, the method offers the ability to adapt the spatial order per cells in each direction. This permits to reduce the CPU time and the memory storage for configuration involving localized spatial elements. An other interest of this adaptable spatial order is to reuse a given mesh in various frequency spectrum while keeping a good accuracy simply by increasing or decreasing the order in the cells. For these specificities, the method we proposed is a natural evolution of the usual FDTD scheme to be considered in future applications.
To address EMC problems, we have proposed a wire model permitting to take into account thin oblique wires. This model, coupled to our scheme, has been numerical validated by comparisons with FDTD and MoM.
In future works, we will look for an hybridization between this method and another higher order method but defined on unstructured meshes. Actually, the basis functions used to define our scheme seems to be well-suited to the com-(a) Currents obtained at P 1.    of the geometry are curved but keeping accuracy and good CPU performances.