A 3D Unstructured Mesh FDTD Scheme for EM Modelling

The Yee finite difference time domain (FDTD) algorithm is widely used in computational electromagnetics because of its simplicity, low computational costs and divergence free nature. The standard method uses a pair of staggered orthogonal cartesian meshes. However, accuracy losses result when it is used for modelling electromagnetic interactions with objects of arbitrary shape, because of the staircased representation of curved interfaces. For the solution of such problems, we generalise the approach and adopt an unstructured mesh FDTD method. This co-volume method is based upon the use of a Delaunay primal mesh and its high quality Voronoi dual. Computational efficiency is improved by employing a hybrid primal mesh, consisting of tetrahedral elements in the vicinity of curved interfaces and hexahedral elements elsewhere. Difficulties associated with ensuring the necessary quality of the generated meshes will be discussed. The power of the proposed solution approach is demonstrated by considering a range of scattering and/or transmission problems involving perfect electric conductors and isotropic lossy, anisotropic lossy and isotropic frequency dependent chiral materials.


Introduction
In this paper, we consider the simulation of problems in electromagnetics using unstructured co-volume staggered meshes and a generalisation of the Yee finite difference time domain (FDTD) method [1]. In its standard form, the FDTD algorithm is implemented on structured staggered spatial meshes. The use of staggered meshes enables the interleaving of the unknown nodal values of the electric and magnetic field components, resulting in an explicit algorithm that is second order accurate in both space and time. The algorithm has the additional advantages of simplicity, low computational cost and, as no linear algebra is required, there is no inherent limit to the size of the simulation. In addition, the scheme preserves the energy and amplitude of waves, is divergence free and can readily be parallelised. To retain these benefits, for simulations involving geometries of complex shape, non-uniform and unstructured mesh FDTD implementations have been suggested, such as the generalised Yee algorithm [2] and the Yee-like algorithm [3]. Unfortunately, these methods do not retain the efficiency of the original scheme [4]. An alternative approach to generalising the FDTD algorithm to unstructured meshes is to employ a primal unstructured Delaunay mesh and its orthogonal Voronoi dual graph [5]. With this method, hybrid meshes can be naturally incorporated, without requiring the inter-mesh interpolation, and consequent non-physical diffraction effects mesh interfaces, that plagues many hybrid approaches. The proposed implementation is edge and face based, which means that it is not restricted to specific mesh element types and can handle any form of polyhedron.
This paper is structured as follows: In Sect. 2, we outline the problem of interest and employ a scattered field formulation for Maxwells equation in integral form for problems involving an isotropic lossy dielectric material and/ or a perfect electric conductor (PEC). Section 3 introduces the FDTD method and illustrates the discretisation of the Maxwell equations on both a structured and an unstructured mesh. Section 4 is devoted to the mesh generation process, which is one of the most crucial parts for the successful implementation of the proposed scheme. The treatment of the different kinds of boundary conditions is discussed in Sect. 5. Section 6 describes the method adopted for the calculation of the distribution of the radar cross section (RCS). In Sect. 7, we validate our algorithm by comparing the numerical and analytical solutions for the RCS distribution for problems involving, dielectric lossy and coated spheres of different electrical lengths. The calculation of the RCS distribution of a more complex PEC aircraft geometry and the simulation of transmission of an electromagnetic pulse through a radome are presented in Sect. 8. In Sects. 9-11, we move on to consider anisotropic lossy materials, where the material parameters, electric permittivity , the magnetic permeability and the electric and magnetic conductivites , m , become second order tensors. This allows us to consider the simulation of problems involving crystals, or composite materials, where the orientation of fibres plays an important role. In Sects. 12-13, we consider a further extension that allows for the modelling of bi-isotropic materials. Materials of this type are frequency dependent and have an additional parameter, the chirality, which couples the electric to the magnetic field, making the simulation even more challenging. Brief conclusions are presented in Sect. 14.

Problem Formulation
We consider problems involving the simulation of the interaction of an incident electromagnetic wave, generated by a source located in the far field, and a general medium surrounded by free space. We assume, initially, that the medium might represent a dielectric or a closed region of very high electrical conductivity, as illustrated in Fig. 1. A closed region of very high electrical conductivity is generally approximated as a PEC. With this approximation, electromagnetic fields do not penetrate the conducting surface and the outer boundary of the PEC forms the inner boundary of the problem domain. For such simulations, it is convenient to split the total electric and magnetic fields into incident and scattered components, according to with the subscripts tot, inc and scat referring to the total, incident and scattered fields respectively. To enable a numerical solution of this problem, the governing Maxwell equations are expressed in terms of the scattered field components [6] and written in the integral form provided by the (1) tot = inc + scat tot = inc + scat Laws of Ampère and Faraday [5]. For a three dimensional lossy dielectric medium, the resulting equations become and Here, A denotes the closed curve bounding a surface A, d is an element of area of the surface A, directed normal to the surface, and d is an element of contour length in the direction of the tangent to the curve A . In addition, denotes the electric permittivity of the medium, its magnetic permeability, its electric conductivity and m its magnetic conductivity. We will often consider plane incident waves, of free space wavelength 0 . In this case, if o and o denote the electric permittivity and magnetic permeability of free space respectively, we can write where is the unit wave vector, is the position vector, = 2 ∕ 0 is the angular frequency and = √ o ∕ o is the impedance of free space.

The Standard Algorithm
The standard FDTD scheme [1,7] is a popular solution algorithm in computational electromagnetics. In its original form, the algorithm is implemented on two mutually orthogonal staggered structured cartesian meshes in a straightforward manner. For example, on a general (x, y, z) cartesian mesh, with mesh spacing x , y , z in the three coordinate directions, the location of the unknown electric and magnetic field components is illustrated in Fig. 2. These components are advanced in time using a staggered leapfrog scheme. In the resulting algorithm, with the superscript n denoting an evaluation at time n t , where t is the time step, and subscripts (I, J, K) denoting an evaluation at the point (I x, J y, K z) , typical discretized forms of Eqs. (2) and (3) in the scattered field formulation [6] are and (5) Similar equations can be written directly for the other field components. For a stable implementation of this explicit scheme, the magnitude of the time step employed has to be restricted and, for a uniform cartesian mesh, this stability criterion may be written as c √ 3 t ≤ l , where c = ( ) −1∕2 is the speed of wave propagation through the medium and l is the edge length in the mesh [7]. It should be noted that, since the incident fields are defined analytically in Eq. (4), the time derivatives of E inc,x and H inc,x can be evaluated directly for use in Eqs. (5) and (6). In this form, this explicit algorithm is second order accurate, in both space and time, on a uniform mesh. The computational implementation of this FDTD method, characterized by its low storage requirements, of around 75 Mbytes per million nodes, and its low operation count, provides a highly efficient simulation process.

Unstructured Mesh Algorithm
Although the cartesian mesh algorithm is very efficient, unstructured meshes are more suitable for use in the solution of problems involving regions of complex shape. For the implementation of a corresponding unstructured mesh FDTD procedure, a suitable pair of mutually orthogonal meshes has to be identified. An obvious choice, given the widespread availability of automatic unstructured mesh generators, is to select a primal Delaunay tetrahedral mesh together with its Voronoi dual. These meshes are mutually orthogonal, in that each Voronoi face is a perpendicular bisector of the corresponding Delaunay edge, while each Delaunay face is perpendicular to the corresponding Voronoi edge. It is assumed that the Delaunay mesh employs N D e and the Voronoi mesh N V e edges. Each Delaunay edge has an associated closed loop of Voronoi edges and each Voronoi edge is surrounded by a closed loop of Delaunay edges. The unknowns in the unstructured mesh implementation of the FDTD scheme are located at the midpoints of these edges. The unknown at the centre of the ith Delauany edge corresponds to the projection, E scat,i , of the scattered electric field onto the direction of the edge. The unknown at the centre of the jth Voronoi edge corresponds to the projection, H scat,j , of   Fig. 3a. Similarly, the numbers i j,k , k = 1, … , M D j refer to the M D j edges of the Delaunay face corresponding to the jth Voronoi edge, as illustrated in Fig. 3b.
With these staggered equations, the magnetic field is advanced over the dual mesh at the half time step, using Eq. (14), and the electric field is advanced over the primal mesh at the full time step, using Eq. (13). It should be noted that, with the appropriate interpretation, the discrete approximations of Eqs. (13) and (14) reduce to the basic FDTD algorithm, when they are employed on a pair of staggered regular orthogonal cartesian meshes. For post-processing purposes, a local least squares problem has to be solved to determine all components of and at the vertices of the primal and dual meshes respectively. For a stable implementation of this explicit unstructured mesh scheme, practical experience has shown that a relationship of the form generally results in stability, where S f is a safety factor that depends on the mesh and is normally given a value in the range 0.8-2.0.

Structured Mesh
Structured meshes of hexahedra are readily generated and the meshes employed for wave propagation problems are generally constructed to consist of cubes of uniform side length, . For problems involving waves of characteristic wavelength, 0 , a value of in the range 0 ∕30-0 ∕10 is typically employed for practical applications of the FDTD scheme on regular cartesian grids [7].

Unstructured Mesh and the Ideal Unstructured Mesh
The unstructured mesh discretisation employed in Eqs. (13) and (14) will be second order accurate provided that each Voronoi face is a perpendicular bisector of the corresponding Delaunay edge. This criterion is satisfied automatically for meshes generated by the Delaunay method, but the need to recover a given boundary triangulation may result in the loss of this property. In addition, each Delaunay face must be a perpendicular bisector of the corresponding Voronoi edge and the intersection point, between the edge and the face, must be at the centroid of the face. This criterion can only be achieved for a perfect tetrahedron, which is a tetrahedron for which all faces are equilateral triangles. In this case, the centre of the circumsphere lies inside the tetrahedron and coincides with its centroid. Failure to ensure that this criterion is satisfied means that the corresponding contour integral in Eq. (13) is approximated over a surface that does not intersect the corresponding Voronoï edge, as illustrated in Fig. 4. In this case, the accuracy of the approximation of the integral cannot be guaranteed.
It follows that certain requirements should be imposed on a Delaunay primal mesh and its Voronoi dual mesh, to ensure the successful implementation of the unstructured mesh FDTD solution method. These requirements are: 1. all Voronoi and Delaunay edge lengths should be bounded from below; 2. the centre of the circumsphere for each Delaunay tetrahedron should lie within the tetrahedron; 3. all Delaunay and Voronoi edge lengths should be bounded from above by a value that is not significantly greater than the specified value of ; 4. any deviation in the location of the midpoint of a Voronoi edge from the actual point of intersection with the corresponding Delaunay face should be minimised; 5. any deviation in the location of the circumcentre of a tetrahedron from its centroid should be minimised.
In this list of requirements, the first two entries are the most important in securing a practical, stable implementation of the unstructured FDTD algorithm. For complex three dimensional geometries, it is not always possible to satisfy exactly the last two requirements, so that these criteria have to be relaxed in practice. The third requirement guarantees a sufficiant resolution for the propagation of the electromagnetic wave.
In two dimensions, the ideal unstructured mesh is the space filling mesh of equilateral triangles [8]. The number of nodes connected to each node is 6, the Voronoi edge length l V = l D ∕ √ 3 ≈ 0.56 l D and all the vertex angles are 60 • . In three dimensions, a tetrahedralmesh made of elements with faces in the form of equilateral triangles is not space filling. It has been shown [9][10][11] that an ideal space filling mesh consists of equal non-perfect tetrahedra, with each face in the form of an isosceles triangle, with one side of length l D long and two shorter sides of length l D short = √ 3∕2 l D long . Six such tetrahedra form a parallelepiped tiling the space, as illustrated in Fig. 5. This configuration can be shown to maximise the minimum Voronoi edge for a fixed element size. All Voronoi edges have the same length The number of nodes connected to each node is 14. All vertices have an acute angle of 71.5 • and, hence, for each element, the centre of the circumsphere is located inside the element. However, certain dihedral angles are equal to 90 • , which is larger than the value 70.5 • for the dihedral angle of the perfect tetrahedron. The elements in this three dimensional ideal mesh satisfy the requirements set out above, but do not, of course, fit general boundaries. It can be readily  demonstrated that the solenoidal preserving property of the original FDTD algorithm is also valid for the discretisation of Eqs. (13) and (14) on these ideal meshes.

Unstructured Mesh Generation
Traditional automatic unstructured mesh generation methods, such as the advancing front technique [12] and the Delaunay triangulation [13], or their combination [14], are not designed to guarantee the creation of an unstructured mesh meeting the FDTD requirements set out above, even in the two dimensional case. These methods generate meshes in which the element edge length is acceptable, but the corresponding Voronoı diagram is often highly irregular and can include some very short Voronoı edges. This means that these methods cannot guarantee the regularity of the edge lengths of the dual mesh and the absence of bad elements. Although methods based on swapping, reconnection and smoothing [15] can be used to improve mesh quality, experience has shown that boundary constraints usually mean that their use does not result in a suitable unstructured mesh for the FDTD method, even in two dimensions [8].
An alternative approach is to generate meshes using an iterative constrained centroidal Voronoi tessellation (CVT) [16], followed by mesh quality optimisation. The iteration process adopted is Lloyd's algorithm [17] which, given an initial mesh, relocates the nodes, to the mass centroids of the corresponding Voronoi cells, and then creates a new Voronoi tessellation of the relocated nodes. The process is repeated until all nodes are close enough to the corresponding centroids. The initial tetrahedral mesh is constructed using the vertices of the ideal mesh as the point distribution. In addition to relocating the nodes, the CVT scheme changes the mesh topology and, although the quality of final mesh is much higher than the quality of the initial mesh, typically around 10% of the elements in the final mesh may still be found to violate the mesh quality measures that are required in the current context.

Unstructured Mesh Optimisation
The unstructured mesh is optimised in an attempt to ensure that both the primal and the dual mesh are of the highest possible quality. The requirement that a dual edge must be a bisector of the corresponding Delaunay edge is relaxed. At the same time, the corresponding dual mesh vertex is moved to a point which still ensures orthogonality between the two grids and which lies inside the corresponding primal element. The procedure adopted is similar to that found in the concept of a power diagram [18], where the perpendicular bisector of an edge is moved parallel to itself by an amount that is proportional to a weight assigned to the edge nodes.
To illustrate how this optimisation process can be achieved, consider, in two dimensions, the case when the Voronoi vertex, located at the centre, O, of the circumcircle enclosing a triangular element ABC, lies outside the element. The point O is located at the intersection of the perpendicular bisectors of the three edges of the triangle. For edge AB, the perpendicular bisector is identical to the common cord of the two circles, each having the same radius as the circumcircle, centred at the points A and B, as illustrated in Fig. 6a. When the radii of these two circles are changed, the new common cord remains perpendicular to AB but loses its property of being a bisector. This is illustrated in Fig. 6b. In this case, with careful choice of the radii, the common cords of the circles centred at the nodes of each edge still intersect at a common point O. The vertices and edges of the dual mesh are modified accordingly, with the dual mesh vertex for the element being relocated to this common point. This same process can be applied to a tetrahedral element in three dimensions.
These ideas are employed as the basis of a mesh optimisation process that is used to locate the position O of the dual mesh node, O, for each element, as close as possible to the element centroid. This will ensure that each element contains its dual node and guarantees the essential criterion of mesh orthogonality. For a tetrahedron element, ABCD, the dual node is located at the intersection point of the common cords of the spheres centred at the vertices of each element, with the sphere at A having a radius of R A . If the radius employed at each vertex of a tetrahedron is the same, then the dual node will coincide with the circumcentre. Suppose that O = (x O , y O , z O ) and that the vertex coordinates are E = (x E , y E , z E ) for E = A, B, C, D . The optimisation is achieved, during the CVT process, by changing the radius associated with vertex A by an amount R A , where the new dual node location for each tetrahedron is computed by solving the system where the summation convention is employed for i = x, y, z.

A Hybrid Mesh Approach
The unknowns in the FDTD algorithms are associated with the edges of the primal and dual meshes. It follows that, for a fixed number of points per wave length, an unstructured mesh made of ideal tetrahedral elements will over-resolve the wave, by about 30%, compared to a structured cartesian mesh. To obtain a resolution similar to that achieved with a structured cartesian mesh, the number of points per wavelength must be reduced when (16) generating the ideal tetrahedral mesh. Unfortunately, such a reduction can result in the under-resolution of boundary geometries, particularly in regions of high curvature. This suggests that a hybrid mesh approach could result in a computationally effective compromise. This has been demonstrated for two dimensional electromagnetic scattering simulations, which have been performed accurately and efficiently using such a hybrid approach, with a mesh of triangles used near complex boundaries and a cartesian mesh used elsewhere [19]. Initial attempts at extending these ideas to three dimensions, for problems involving simple geometries, have also been reported [11,20]. With such an approach, the resolution employed is selected to be that appropriate for the operation of the FDTD algorithm on a regular cartesian mesh. This will result in an unstructured mesh that over-resolves the solution in the free space region immediately adjacent to the scatterer. However, as it can be estimated that, typically, only around 10% of the total number of nodes will lie in the unstructured mesh region in three dimensions, following this approach will only result in a small increase in the total number of degrees of freedom, which is regarded as being acceptable given the improved boundary resolution.

Hybrid Mesh Generation
Consider the problem of discretising a general domain, using tetrahedral cells in the vicinity of the scatterer and regular cartesian hexahedral cells elsewhere. To provide a consistent hybrid mesh, these two groups of cells will be connected through a set of interface polygonal cells. The process of generating a hybrid mesh for the computational domain surrounding a general scattering obstacle is accomplished in four stages. In the first stage, an unstructured triangulation of the surface of the scatterer is produced [21] and the triangulation is then placed inside a hexahedral box. The surface of the box forms the outer surface of the computational domain. The region inside the box is discretised using a regular cartesian mesh of cubes, of a prescribed edge length , as illustrated in Fig. 7a. Cubes within a prescribed distance of the scatterer, or lying internal to the scatterer, are removed in a second stage, to create a staircase shaped surface that completely encloses the scatterer, as shown in Fig. 7b. In the third stage, a point distribution is specified to completely cover the unmeshed region, as in Fig. 7c, with those points from the distribution that lie either inside the scatterer or outside the staircase surface being removed, as shown in Fig. 7d. The fourth stage consists of using the points that remain to generate an unstructured tetrahedral mesh in the region between the surface triangulation of the scatterer and the staircase surface. However, this mesh generation process cannot start directly, as the exposed faces of the staircase surface are all squares. The traditional approach to handling the problem of fitting a hexahedral mesh to a tetrahedral mesh is either to place a pyramidal element on each exposed square face, leading naturally to a consistent mesh, or, at the expense of introducing a hanging edge, to divide each exposed square face into two triangles. The use of pyramids is not adopted here, because the geometrical constraints of the staircasing imply that the centre of the circumsphere of each pyramid would lie outside the pyramid, which is something that should be avoided. Instead, the hanging edge approach is employed, but is modified to produce a consistent mesh of arbitrary shaped polyhedral cells that are appropriate for use with the co-volume method. In this approach, each cube that forms part of the staircase surface is divided into six tetrahedra. For a cube with a single exposed face, a tetrahedron that contains one of the staircase surface triangles is removed from the mesh, as shown in Fig. 8a and b. A similar procedure is employed for cubes with either two or three exposed faces, producing the configurations shown in Fig. 8c and d. This process will replace the original staircase surface by a set of triangular faces that are appropriate as the starting point for the generation of the tetrahedral mesh. Each polyhedral interface cell is stored as a set of tetrahedra. The tetrahedra in each set have the same circumsphere and they are merged again during the solution process.

Element Merging
The primal Delaunay mesh will be degenerate if two, or more, adjacent tetrahedral elements share the same circumsphere. In this case, certain Voronoi edges will have zero length, certain Voronoi faces will have zero area and the discrete approximations of Eqs. (13) and (14) become invalid. However, the approximations will remain valid if, at the outset, each set of such adjacent degenerate tetrahedra is merged, to create a polyhedron of the appropriate shape, and the unknowns and the lists of Delaunay and Voronoï edges are modified accordingly. The primal mesh is then regarded as a hybrid of tetrahedra and polyhedra and the dual mesh is updated appropriately. To avoid creating a new Fig. 7 Illustration of the stages involved in the generation of an initial hybrid mesh: a the regular cartesian mesh covering the interior of the computational domain; b creation of the staircase surface at a prescribed distance from, and completely enclosing, the scatterer; c the introduction of the point distribution; d the final point distribution, retaining only those points lying in the region between the scatterer and the staircase surface Fig. 8 Treatment of the cubes forming the staircase surface to enable a consistent hybrid mesh to be generated: a typical cube; b cube with one face exposed; c cube with two exposed faces; d cube with three faces exposed data format to cater for the shapes that are created in this way, a polyhedron is stored as a set of tetrahedra, each with the same circumsphere. These tetrahedra are merged again during the solution process.
Voronoi edges having a 'small' length may still be present in the mesh and this will have a detrimental effect on the efficiency of the explicit time stepping process. To alleviate this problem, certain 'small' Voronoi edges are also collapsed and the corresponding primal elements are merged to form a polyhedron. For the examples that are presented here, this merging is performed for all Voronoi edges with a length that is less than 0.01 . Note that, if two merged polyhedra are adjacent, it is important to ensure that the triangles forming the common face are co-planar. Violation of this criterion implies that the original merging process should not be allowed. Here, two triangles are taken to be co-planar provided their dihedral angle in less than 1 • .

Material Interfaces
At interface boundaries, the material parameters change over the area of integration. In this case, the values of , , and m are replaced by appropriately averaged values av , av , av and mav . For the standard FDTD algorithm, the three most common averaging techniques, at an interface between materials with generic parameters a 1 and a 2 , are the arithmetic mean, the harmonic mean and the geometric mean. It has been reported [22] that the convergence rates obtained with the geometric and harmonic means deteriorate as the mesh size decreases, whereas second order convergence is maintained with the arithmetic mean. For this reason, the arithmetic mean average is adopted here and modified for use with unstructured meshes, where the cell size is not necessarily uniform.
To illustrate how an arithmetic mean is employed in the construction of the material parameters, consider the situation at an interface, between material 1 and material 2, illustrated in Fig. 9. At the interface, material property a on Delaunay edge i is assigned the averaged value a avDel,i , computed as In this equation, the term a cell(k) denotes the material parameter assigned to cell k, surrounding Delaunay edge i. Cell k, with volume w k , is defined as the region spanned by the end points of Delaunay edge i, the point of intersection of the Voronoi edge with the Delaunay face and the position of the circumcentre of cell k. For example, in Fig. 9a, w 1 would be the volume spanned in material 1 by the points P1, P5, P9, P8, and a cell(1) = 1 or a cell(1) = 1 . Similarly, w 2 would correspond to the volume spanned in material 2 by the points P1, P5, P9, P10, so that a cell(2) = 2 or a cell (2) Similarly, material property b, on a Voronoi edge, j, is assigned the average value b avVor,j and is computed using the weighted average formulae In this equation, the term b cell ( ) denotes the material parameter assigned to cell . The length of the Voronoi edge inside cell is g and this corresponds to the distance between the point of intersection of the Voronoi edge j with the Delaunay face and the circumcentre of the cell. For example, in A PEC boundary may also be simulated in this fashion, treating it as a resistive sheet. In this case, the Delaunay and Voronoi edges forming the interface are assigned a very high conductivity value in the order of 10 6 S/m. It can be readily demonstrated that this averaging process on an unstructured mesh corresponds to the arithmetic mean on a structured mesh [23].

Far Field Boundary
In the far field, the scattered fields should consist of outgoing waves only. This radiation condition is imposed by considering the outermost layers of the structured hexahedral mesh to be in the form of an artificial absorbing perfectly matched layer (PML) [24]. The PML region is defined to be of a constant thickness and a standard PML formulation is adopted [7,25].

Calculation of the RCS
The RCS provides a measure of the power scattered from a target obstacle towards a receiver. In numerical simulation, we evaluate the bi-static RCS, which is a function of aspect angle and bi-static angle. The scattered wave is decomposed into two components. The first component has the same polarization as the transmitted signal and is referred to as co-polarized, while the second component has an orthogonal polarization state compared to the transmitted signal and is referred to as cross-polarized. For numerical simulations, the RCS distribution is of practical interest because an analytical solution is available for certain simple cases. To compute this distribution, a near to far field transformation procedure [26] is employed to obtain the required far field information from the computed near field data. To achieve this, a closed collection surface, S, completely enclosing the scatterer, is constructed and, when steady state conditions have been achieved, a further cycle is computed. During this cycle, the harmonic solution produced by the time domain solver is used to calculate the phasors, which enable fictitious electric and magnetic currents to be defined on S. Standard electromagnetic theory can then be used to determine, using the surface equivalence theorem, the components of the scattered electric field in the far field [27]. With these values, the distribution of the quantity is computed, where (r, , ) denotes a standard spherical polar coordinate system. For the results presented here, the collection surface, S, is taken to be the surface formed in the mesh by removing all but the first layer of tetrahedra that are attached to the scattering surface and it is the distribution of the quantity that is displayed.

Scattering by a 2 0 PEC Sphere
The first example involves modelling the interaction between a monochromatic plane wave, of wavelength 0 , and a PEC sphere of diameter 2 0 . Simulations are performed initially using the standard FDTD algorithm and regular structured meshes with mesh spacing = 0 ∕15 and = 0 ∕60 . Fig. 10 shows the staircased representation of the surface of the sphere on these structured meshes. The simulation is repeated using the unstructured FDTD algorithm, with a conforming unstructured mesh with global spacing = 0 ∕15 . For this example, cuts through a typical hybrid mesh are shown in Fig. 11. The analytical [26] and computed RCS distributions are compared in Fig. 12. It is apparent that the distributions computed on the fine structured mesh do not achieve the accuracy of the solution computed on the conformal unstructured mesh. This indicates that, for this problem, a significantly finer structured mesh must be used with the standard FDTD method to achieve accuracy comparable to that produced with the geometry conforming unstructured mesh FDTD method.
The comparison is repeated for the case of 2 0 dielectric sphere, with the relative values r = 2 and r = 1 assigned to the material parameters in the dielectric. In addition, it is assumed that = m = 0 . The corresponding RCS distributions are displayed in Fig. 13. In this case, the fine structured mesh is 64 times larger than the unstructured mesh and the corresponding allowable time step is four times smaller. The result is that, for this example, the fine structured mesh FDTD algorithm CPU requirement is 256 times larger than that of the unstructured FDTD method.
The same example is used to demonstrate the order of convergence of the method on a general unstructured mesh. Meshes with global spacing of 0 ∕4 , 0 ∕8 , 0 ∕15 and 0 ∕20 are employed and the L 2 norm of the errors in the computed distributions of the RCS are evaluated after 20 cycles of the incident wave. Fig. 14a shows a convergence rate of 1.59 for the co-polarized RCS and Fig. 14b shows a corresponding rate of 1.48 for the cross-polarized RCS, giving an averaged value of 1.54 for the order of convergence. For a lossy dielectric sphere, with the same values for r , r = 1 and m , but with the electric conductivity = 0.7 S/m, the bi-static RCS distributions, computed after 20 cycles of the incident wave, are again seen to be in visually in a good agreement by comparing it with the analytical solution as can be seen in Fig. 15.

Scattering by a Coated 2 0 Sphere
We now consider the scattering of a plane wave by a sphere, of electrical length 0 , that is coated with a uniform dielectric, of thickness 0 . A cut through the mesh used to represent the sphere and the coating is shown in Fig. 16a. For the simulations reported here, the PML consists of 10 layers of hexahedral elements and the complete mesh consists of 876, 116 cells, with 1, 673, 527 Delaunay edges and 2, 076, 019 Voronoi edges. Initially, the inner sphere Fig. 16a is assigned the material parameter values r = 1 , r = 1 , = 10 6 S/m and m = 0 S/m, so that it effectively is a PEC. The coating is assumed to consist of a conducting dielectric material, with parameters r = 2 , r = 1 , = 0.7 S/m and m = 0 S/m. The computed bistatic RCS distribution is compared with the analytical distribution in Fig. 17. An illustration of how the scattered electric field propagates through the sphere is given in Fig. 16b. Finally, the sphere is considered to be a dielectric, characterized by the parameters r = 2 , r = 1 , = m = 0 S/m, and the dielectric coating is characterized by the parameters r = 1.5 , r = 1 , = m = 0 . The computed bi-static RCS distribution is compared with the analytical distribution in Fig. 18. For these two examples, the computed and analytical RCS distributions are seen to be in excellent agreement.

Scattering by a 15 0 PEC Sphere
To illustrate the performance of the procedure when it is applied to problems involving bodies of larger electrical length, we consider the simulation of the interaction between a plane single frequency incident wave and a PEC sphere of electrical length 15 0 . The incident wave propagates in the x direction, the global element size is = 0 ∕10 and the outer surface of the PML region is, at its closest point, located at a distance of 2 0 from the surface of the sphere. The staircase boundary is placed at a distance of 0 ∕2 from the surface of the sphere and the thickness of the PML region is set equal to 0 . The generated mesh has 7, 006, 329 vertices, of which 2, 460, 486 vertices are located within the PML. The unstructured mesh region near the sphere has 848, 205 vertices and, before the application of the mesh enhancement procedure, 4.3% of the cells in this region are such that the centre of the cell circumsphere lies outside the cell. Following the application of mesh enhancement, the percentage of such elements is reduced to 1.1% . By allowing the merging of adjacent cells to form polyhedron cells, this percentage is further reduced to 0.2% . The majority of these undesirable cells are located adjacent to the staircase boundary. A view of the surface mesh and of a cut through the final mesh is shown in Fig. 19.
The solution is advanced through 30 cycles of the incident wave, using 180 time steps per cycle. A view of the distribution of the computed contours of the E scat2 component of the scattered electric field, on the scattering surface and on a cut through the mesh, is displayed in Fig. 20. For comparison, the simulation is repeated using an explicit nodal low order finite element time domain scheme (FETD) [28,29]. The FETD scheme is significantly less efficient requiring 8 times longer to advance the solution for one complete cycle of the incident wave. The RCS distributions computed with both schemes are seen to be in excellent agreement with the exact distribution in Fig. 21.
To demonstrate the effectiveness of the advocated mesh generation approach, an alternative mesh is generated, using a basic Delaunay algorithm with automatic point insertion [13], in the region between the staircase and the sphere. In this case, 47% of the tetrahedra are such that the centre of the circumsphere of the tetrahedron is located outside the tetrahedron. This percentage is reduced to 4.3% after applying the traditional mesh enhancement procedures of  edge swapping, edge collapse and Laplacian smoothing [15].
Using the FETD scheme on this mesh, it is found that 162 time steps are required to advance the solution for one cycle.
If the merging criterion is applied, stability of the FDTD scheme on the resulting mesh requires the use of 5491 time steps to advance the solution through one complete cycle of the incident wave. This shows quite clearly that meshing techniques that are only designed to produce a good quality Delaunay mesh will not, generally, be suitable for use with the unstructured mesh FDTD algorithm.

Basic Algorithm: Geometries of More Complex Shape
To demonstrate its predictive capabilities, the described approach is now applied to problems involving more complex geometries.

Scattering by a PEC Aircraft
This example involves the interaction between a plane single frequency incident wave and a PEC full aircraft configuration. The wave frequency selected is such that the electrical length of the aircraft is 10 0 . The main axis of the aircraft is aligned with the x direction and the incident wave propagates in this direction, with the wave impinging directly onto the nose of the aircraft. The staircase surface is located at a distance 0 ∕2 from the scatterer and the thickness of the PML is taken to be 0 . Two meshes, with an element size = 0 ∕15 , are created. The first mesh is designed for use with the FETD method, with the region of tetrahedral elements generated using the standard Delaunay approach with automatic point insertion [13]. This mesh contains 4, 274, 438 vertices, of which 2, 129, 120 vertices are located within the PML region. The second mesh, for use with the FDTD co-volume algorithm, is generated using the approach outlined above and contains 4, 378, 740 vertices, with 2, 208, 512 vertices in the PML. For this second mesh, the percentage of tetrahedra having the property that the centre of a circumsphere of a tetrahedron is located outside the tetrahedron is reduced, from 16 to 2.2% , by using the proposed mesh enhancement procedure. This percentage is reduced further to 0.5% by using the advocated merging process. The majority of these undesirable elements  again lie in the vicinity of staircase interface. A view of the discretised aircraft surface is given in Fig. 22 and a detail of this view together with a view of the mesh on a cut through the domain is given in Fig. 22b. The solution is advanced for 40 cycles of the incident wave. The FETD computation on the first mesh takes 6.5 times longer per cycle compared to the FDTD scheme on the second mesh. It is interesting to note that the time required for the generation of the mesh for the co-volume scheme is 81 min. A detail of the computed distribution of contours of the E scat,x component of the scattered electric field on the aircraft surface is given in Fig. 23. The RCS distributions computed by the FETD method and the FDTD algorithm on these meshes are compared in Fig. 24. Note that, in the context of scattering by a PEC dodecahedron using a completely ideal mesh, we have demonstrated previously [11] the convergence properties of the co-volume algorithm and its improved performance on a given mesh, for problems involving singularities, over that of the FETD method. This means that the discrepancies that are apparent here between the two RCS distributions for this example are probably due to the inability of the nodal finite element scheme to adequately represent the solution in the vicinity of geometrical singularities, without the use of additional local mesh refinement [19].

Illumination of a Radome by a Narrow Band Pulse
This  Fig. 25a. To verify that one layer of elements is  where FT corresponds to the Fourier transform of the electric field. The value of the transmission obtained from the finer mesh is compared with that obtained from the coarser mesh in Fig. 27. The small differences between these results suggests that acceptable accuracy may be achieved here by using just one layer of cells to model the thickness of the radome.
To study the effect of different frequencies on the same mesh, two narrowband pulses are considered. Both have a pulse width of 1 × 10 −9 s , with the first centred at 1 GHz and the second centred at 1.2 GHz. In the frequency range where the pulse spectra overlap, the transmission is expected to be the same. The predicted variation of the computed transmission efficiency with frequency, for both pulses, is plotted in Fig. 28. The mesh generation for this problem is challenging due to the rather thin thickness of the layer and extra care needs to be taken when collecting the total field history inside the radome. To reduce the influence of the dispersion error the incident field values in the update Eqs. (13) and (14) are not taken from the known function of Eq. (21) but are interpolated from a 1D FDTD code running in parallel. This 1D code runs with the same time step and with a uniform mesh, with the mesh spacing equal to the spacing of the regular hexahedra in the actual problem.

Modelling Anisotropic Lossy Materials
For anisotropic materials, the electromagnetic parameters, such as permittivity, permeability and conductivity, vary in different directions, so that they must be treated as tensors.
Such materials offer many new and interesting possibilities in engineering, e.g. a thin anisotropic coating may significantly change the RCS distribution of an aircraft. With their additional mechanical strength and weight advantages, compared to metals, composite anisotropic materials, with applications initially limited to stealth aircrafts, satellites and space shuttles, are now part of everyday life [30]. In addition, anisotropy is at the heart of the new metamaterials, which have electromagnetic properties that do not occur naturally, such as a negative index of refraction [31]. It is clearly important to extend the capabilities of the FDTD method to enable the modelling of problems involving such anisotropic materials. The approach that will be followed was initially introduced within the context of a total field formulation [30], but our unstructured mesh extension employs a scattered field formulation.
For a three dimensional lossy anisotropic dielectric medium, the Laws of Ampère and Faraday are expressed, in a scattered field form, as and where denotes the electric flux density and the magnetic flux. In addition, Ī is the unit matrix, ̄ is the electric permittivity tensor, ̄ is the magnetic permeability tensor and ̄ and ̄m are the electric and magnetic conductivity tensors respectively. The constitutive equations define the relationships between the electric and magnetic fields and the electric flux density and the magnetic flux.

Unstructured Mesh Algorithm
For the solution of this equation set, the unknown at the centre of the ith Delaunay edge, illustrated in Fig. 3a, corresponds to the projection, (D scat,i , E scat,i ) , of the scattered electric field onto the direction of the edge. The unknown at the centre of the jth Voronoi edge, illustrated in Fig. 3b, corresponds to the projection, (B scat,j , H scat,j ) , of the scattered magnetic field onto the direction of that edge. Direct discretization of the Laws of Ampère and Faraday then leads to the equations and Here, the bracket notation ⟨ ,̂ i ⟩ is used to denote the projection of a field vector along the ith edge. At interface boundaries, the material parameters in the Eqs. (27) and (26) are not constant and are replaced by averaged values ̄a v , ̄a v , ̄a v and ̄m av . In this case, every component of the material parameter tensor is averaged following the approach employed earlier for the isotropic case [32]. To represent the vector-matrix multiplications that are required, we define and these definitions enable us to write Eqs. (26) and (27)

Obtaining Approximated Field Vectors from Edge Projections
To perform the necessary matrix-vector multiplications in Eqs. (32) and (33), we need to obtain the field vectors scat and scat and associate them with Delaunay edge i and Voronoi edge j respectively. Unfortunately, we cannot get exact full field vector components directly from edge projections, but we are able to approximate the full field components at any location in the mesh. To achieve this, we assume that, with a set of three orthogonal vectors 1 , 2 , 3 , a general vector can be reconstructed as in terms of the projections As the mesh is unstructured, we cannot use this method directly to obtain an exact field vector (Please note that Einstein notations is not employed in this paper). The trivial solution would be to consider all the edges connected to one node, as depicted in Fig. 39b, and solve a system of equations. This will give us a field vector at node n 1 . As one edge is always formed by two points, we can do the same for point n 2 and average these field vectors and project to the corresponding edge. For example, to obtain the displacement field vector scat on the edge i, denoted by the red line in Fig. 29a, we have to construct the set of equations where (31) and ̂ i is the normalized Delaunay edge vector corresponding to Delaunay edge i. The entries in the matrix ̄ , based upon the three components of the vector ̂ i and the projections scat ⋅̂ i , are all known. The displacement field vector, scat , at a node belonging to Delaunay edge i is approximated by considering the sum of the system in Eq. (36) for each Delaunay edge connected to the node, as depicted in Fig. 29b. In this case, we solve the system with Here, ̄ ′ and scat,q ⋅̂ q are known and N is the number of Delaunay edges connected to the node. This system of equations is solved locally, node by node, resulting in an approximated field vector at all nodes of the Delaunay mesh. The computed field vector is then linked to the corresponding edges. A Delaunay edge is assumed to connect the two nodes n 1 and n 2 and the displacement field vector, associated to the Delaunay edge i, is obtained by averaging the electric field vectors at nodes n 1 and n 2 , i.e. scat,i = scat,n 1 + scat,n 2 ∕2 . The same procedure is applied for approximating the magnetic flux vectors scat , but using now the Voronoi edges. However, unlike the standard FDTD scheme, where we have one update equation for each component of the field vectors, the approximations of the full vector fields obtained using

i x e i y e i x e i z e i y e i x e i y e i y e i y e i z e i z e i x e i z e i y e i z e i z
∑ m=1 e q l e q m this method are not good enough for time iterating the field components. It is found that error accumulation causes the algorithm to become unstable. We modified it in order to circumvent this difficulty.

Local Orthogonal Unit Vectors
It is the projection, D scat,i = scat ⋅ i , of displacement field vector scat onto a Delaunay edge i , that is updated in the solution process. For an isotropic material, the electric permittivity and the magnetic permeability are scalars, so that updating the fields only involves scalar multiplication between the field projections and the scalar material properties. For an anisotropic material, matrix-vector multiplications, between the material tensors ( ̄ , ̄ , ̄,̄m ) and the fields ( scat , scat ), are required. To perform these multiplications, we create two linearly independent vectors for each Delaunay and Voronoi edge. Using the stabilized Gram-Schmidt orthonormalization procedure, we generate three orthonormal vectors, as illustrated in Fig. 30. The first vector , to which the two linearly independent vectors ( ′ 2 , ′ 3 ) are added, remains unchanged during the whole process. Each set of three orthogonal vectors represents one local coordinate system, leading to as many local systems as Delaunay and Voronoi edges. We have already seen how we can reconstruct approximated field vector components for each local coordinate system, using the field projections of the surrounding Delaunay or Voronoi edges. The field vectors, obtained in this way, are projected onto the two orthogonal vectors forming the local frame. The first vector of each subset can be immediately updated using the projection equation employed in the isotropic case. For this projection, no error is induced by the field averaging.

Coordinate Transformation
The material tensors are expressed in the global reference frame formed by the orthonormal vectors ̂ , ̂ and ̂ . Three orthonormalized vectors ̂ ′ ,̂ ′ ,̂ ′ form the basis of each local coordinate system. Here, the symbol ̂ denotes a unit vector and the symbol ′ refers to vectors or vector components in the local coordinate system. The Jacobian matrix J , defined by can be used to transform from the global (x, y, z) to a local frame. Each component of the Jacobian can be interpreted as an amplification factor, describing how one coordinate in a given reference frame stretches, shrinks or rotates with respect to another coordinate in another reference frame. In our case, the Jacobian is purely a rotation matrix J R , which can be directly calculated as This rotation matrix has unit determinant, i.e. det (J R ) = 1 . Using a simplified set of Maxwells equations the form invariance [33,34] is easiest explained in the following paragraph by expressing them in the local coordinate system, as The electric and magnetic fields, in local and global coordinates, are related as while a general material parameter tensor, M , is transformed to the local coordinate system as In this way, the coordinate transformation may be absorbed completely into the material properties.

Time Updating
The magnetic field is updated over the dual mesh at the half time step, using Eq. (33), and the electric field is updated over the primal mesh at the full time step, using Eq. (32). We will describe the updating of the electric field projection E n scat . The magnetic field H n scat is updated similarly. The coordinate transformation applied to the terms ā + ,ā − , defined in Eq. (28) (38). These vectors are projected in each of the three orthogonal directions of every local frame, e.g. Fig. 30, and the matrix vector multiplication is performed. In this way, Eq. (33) in the local frame takes the vector form The constitutive equation may then be used to obtain the electric field projection. In practice, we only have to consider the first line of ̄� −1 av , due to the fact that our data storage is based upon field projections along the edges and not field vectors. As a result, we can obtain

Anisotropic Material Modelling Algorithm: Initial Validation
A series of examples, involving scattering of an incident plane wave by an anisotropic sphere, is performed to validate the revised algorithm and to demonstrate its numerical performance. Validation of the implementation is achieved by comparing the computed RCS distributions with those obtained from the open source program Discrete Dipole Scattering (DDSCAT) [35], which is an implementation of the frequency domain discrete dipole approximation [36]. In each case, the incident wave, propagating in the x direction, has free space wavelength 0 and the electrical length of the sphere is 2 0 . The mesh employed is illustrated in Fig. 31 and has an average edge length of 0 ∕20 . The PML region is located at a minimum distance of 0 from the scatterer and is discretised using 10 layers of hexahedra. The complete mesh consists of 876, 116 cells, 1, 673, 527 Delaunay edges and 2, 076, 019 Voronoi edges [23]. It is worth noting that our unstructured mesh FDTD scheme only needs 8 points per wavelength to produce accurate results in free space [8]. However, inside a dielectric, the wavelength, Diel , is less than the free space wavelength, 0 . For the anisotropic case, it is estimated that a mesh spacing of 0 ∕20 mesh should be sufficient here. Despite the fact that coarser meshes could have been used in some cases, we decided to use the same mesh for all the examples.

Electrically Lossy Anisotropic Sphere
For the first example, the material parameters are employed and matrix multiplications are required in Eq. (27). Steady state conditions were attained after ten cycles of the incident wave. The computed RCS distributions are shown to be in excellent agreement with the DDA results in Fig. 32.

Magnetically Lossy Anisotropic Sphere
In the next example, the material parameters of the sphere are characterised, by an anisotropic permeability tensor, together with magnetic conductivity, as Now, matrix multiplication is required in Eq. (26) and steady state is again reached after ten cycles of the incident wave. The computed RCS distribution is compared with the DDA results in Fig. 33.

Analysis of the Full Anisotropic Model
To investigate the time penalty that results from the use of the anisotropic model, two simulations are performed on the same mesh. In the first example, the sphere is taken to be an isotropic lossy dielectric material. In the second example, the sphere is modelled as an anisotropic lossy dielectric material, characterised by the parameters  Comparison with DDA results is not possible for the second example, as the DDA code only allows for the modelling of non-magnetic materials. It is found that the computational cost for the simulation involving the anisotropic sphere is ten times the cost of the simulation for the isotropic sphere. This extra cost mainly arises from Eq. (38), which is a requirement to solve a system of three equations for each Voronoi and each Delaunay node. The y and z components of the electric and magnetic fields respectively are represented in Fig. 34 Given this additional cost, it is reasonable to ask if the current scheme is competitive when compared to the standard FDTD method for this class of problems. We have previously shown [23] that the accuracy of this scheme, with a 0 ∕15 unstructured mesh, compares favourable with that of the standard FDTD method, with a 0 ∕90 or even a 0 ∕120 structured mesh for objects of arbitrary shape. This means that we get comparable results when using a mesh that is 6-8 times coarser. A standard FDTD cell stores the electric and magnetic field components at 12 edges and 6 faces respectively, corresponding to 18 degrees of freedom per wavelength. One FDTD cell may be discretised using 6 tetrahedra, in the worst case scenario, leading to 19 edges and 14 faces and a total of 33 degrees of freedom per wavelength. Although the number of degrees of freedom has nearly doubled, this is more than compensated by the fact that we may use a mesh that is, at least, 6 times coarser. Furthermore, for the same volume in three dimensions, an unstructured mesh with 33 degrees of freedom produces the same accuracy as a structured mesh with 6 3 × 6 = 1296 degrees of freedom. The extra operations, needed for averaging and reconstructing vectors, imply a cost penalty. For each node, we require 18 averaging operations for each field, in addition to 3 × 3 operations of vector reconstruction for each edge. For the 6 tetrahedra lying inside one cube, the number of operations will then be 2 × 8 × 18 + 9 × 33 = 585 . For one cell of a cartesian mesh, 8 operations are required for averaging the two offset components of a field vector, leading to 48 extra operations on one cell node or 144 per FDTD cell. This implies that, when dealing with curved boundaries, to achieve the same level of accuracy, the standard FDTD scheme will actually require 6 3 × 48 = 10,368 operations per computational volume, whereas the current unstructured implementation has only 585 operations.

Anisotropic Material Modelling Algorithms: Geometry of More Complex Shape
To illustrate the predictive modelling capability of the enhanced solution algorithm, we consider again a problem involving the evaluation of the transmission efficiency of a pulse through a radome.The dimensions of the radome remain unaltered, but an anisotropic dielectric material, containing conducting fibres, is used to construct the radome. The mesh employed for the simulation, with 40 degrees of freedom per wavelength, is the mesh that was used earlier for the isotropic material case. The incident plane wave is again the narrow band pulse defined in Eq. (21), polarised in the y direction. Figure 35a is a schematic of a cut through a slab of composite material, in which the fibres are oriented in the ̂ ′′ direction. Reinforced carbon-carbon is such a type of composite where carbon fibres are reinforcing a graphite matrix. This material is for example used as heatshields for space shuttles. For the illustrated slab orientation, the material frame (̂ �� ,̂ �� ,̂ �� ) and the global frame (̂ ,̂ ,̂ ) are identical. The situation is more complicated for the curved shell of the radome, as the orientation of the fibres changes in space. In this case, for each location on the shell, we have a material frame which, in general, differs from the global frame as illustrated schematically in Fig. 35b. Previously, we have only used the coordinate transformation to pass from a global frame ( (̂ ,̂ ,̂ ) coordinates) to a local frame ( (̂ � ,̂ � ,̂ � ) coordinates), linked to each Voronoi or Delaunay edge. For a composite it becomes more complicated as we first have to make a coordinate transformation from the material frame, linked to the orientation of the fibres (referred to as (̂ �� ,̂ �� ,̂ �� ) coordinates) to the global frame ( (̂ ,̂ ,̂ ) coordinates). We then pass as before from the global frame, to a local frame (̂ � ,̂ � ,̂ � ) which is linked to each Voronoi or Delaunay edge. Mathematically these transformations are performed as follows: where J R1 corresponds to the first transformation matrix describing the transformation from the material to the global frame, while J R2 corresponds to the second transformation matrix from the global to the local frame. This is therefore identical to the matrix J R . Both J Ri , i = 1, 2 are computed using Eq. (41), but with respect to the corresponding coordinate system leading to Eq. (53). We create an orthonormal material frame, which simplifies the coordinate transformation as the transformation matrix becomes a simple rotation matrix. The procedure here is illustrated in Fig. 36a. Initially, we only consider Voronoi edges at the dielectric interface. Each Voronoi edge (Vor1,blue) crossing Scattering of a plane wave by a fully anisotropic dielectric sphere of electrical length 2 0 : a contours of H scat,z shown on a crosssection through the mesh; b contours of E scat,y shown on a cross section through the mesh Fig. 35 Composite material: a orientation of fibres in the material frame; b orientation of the fibres inside the radome the face of a tetrahedron is surrounded by three Delaunay edges. We select two of these edges, Del 1 and Del 2 . By construction, these edges are parallel to the surface and taking the cross product between them the creation of the vector �� = Del 1 × Del 2 that is perpendicular to the surface. In the next step we assign �� = Del 1 . By construction ′′ and ′′ are perpendicular to each other. To create the third vector ′′ perpendicular to ′′ and ′′ , the cross product between �� × Del 1 = �� × �� = �� is taken. Following normalization, we have one orthonormal material frame (̂ �� ,̂ �� ,̂ �� ) for each face of the dielectric interface. In the next step, the material frame is linked to the cells inside the composite, by comparing the distances between the intersection points of a Voronoi edge with the interface and the circumcentre of all cells C 1 inside the dielectric. After finding the smallest distance, we link the local material frame (̂ � ,̂ � ,̂ � ) from the surface to the corresponding cell. This process is illustrated in Fig. 36b. To illustrate this process, consider a typical material parameter tensor M , defined with respect to the material frame (̂ �� ,̂ �� ,̂ �� ) . The first coordinate transformation converts the material parameters, from the material to the global frame, according to The second coordinate transformation then converts the parameters from the global to the local frame, a process that has been described earlier. If we consider, for example, Eqs. (28) and (29), the material parameters ̄,̄,̄,̄m have to be replaced by ̄′ ′ ,̄′ ′ ,̄′ ′ ,̄′ ′ m in each equation. The radome is characterized by the material parameters Due to the orientation of the fibres, the properties of the composite are the same in the ̂ �� =̂ �� and the ̂ �� =̂ �� directions, but differ along the ̂ �� =̂ �� axis (see Fig. 36). This means that for the material parameters M (1, 1) �� =M(2, 2) �� ≠M(3, 3) �� , and the off-diagonal components are 0, for ̄′ ′ ,̄′ ′ . The magnetic permeabiility tensor is set to the unit matrix Ī , ̄� � =Ī corresponding to free space and ̄� � m =0 . Typically, a material that minimally attenuates the electromagnetic signal transmitted or received by the antenna is used.
The computed transmission for a 1 GHz narrowband pulse through the composite radome is represented in Fig. 37.

Modelling of Bi-isotropic Materials
Innovative fabrication techniques are currently being employed to create new chiral materials for use in communication technologies. For these materials, a metallic structure, These inclusions need to be much smaller than the wavelength, so that the incoming wave sees the matrix and the inclusions as one homogeneous material. Often, these inclusions take the form of a combination of a straight wire and a ring, to which the electric and magnetic field respectively will couple, leading to a material (Fig. 38a) with a frequency dependent response. For such materials, the electric permittivity and the magnetic permeability are frequency dependent parameters. Furthermore, for specific shapes or arrangement of the metallic inclusions, e.g. the use of left or right handed helices, as shown in Fig. 38b, the electric and magnetic fields will be related by an additional frequency dependent material parameter in the constitutive equation, referred to as the chirality. The chirality allows for a negative index of refraction with, simultaneously, a positive relative permittivity and permeabilty.Chiral materials can also create phenomena such as optical rotatory dispersion (ORD) and circular dichroism (CD). In ORD, the polarization is continuously rotated inside the material, while CD refers to the change of polarisation of the incident wave, from linear to elliptical, due to the different absorption coefficients of a right and a left circularly polarized wave. Such materials are being used to create smaller antennas, super lenses [38], polarizers and radomes or RCS reducing materials.
To extend the modelling approach to enable the simulation of problems involving chiral dispersive materials, a Lorentz model is employed for the electric permittivity and magnetic permeability and a Condon model is adopted for the chirality [39]. Several methods have been developed for modelling frequency dependent materials within a time domain method, including the auxiliary differential equation method [40], the piecewise recursive convolution method [41] and the Z-transform technique [42]. As these three methods have been found to produce generally comparable results, we decided to implement the Z-transform technique. This approach, which demonstrates good convergence close to the resonant frequency [43], is widely used in signal processing and can be viewed as the discrete version of the Laplace transform. It was first used to model chiral materials in 3D with FDTD [44], using the Z-transform to handle the frequency dependence of the material parameters. First order approximations were used in the Z-domain and the Z-transform coefficients were derived analytically. In a different approach [45], second order accuracy in the Z-domain is achieved by employing a bilinear transformation for calculating Padé approximants [46]. In our approach, chiral materials are modelled by employing a generalization of this method.

Problem Formulation
It is not possible to employ a pure scattered field formulation when modelling problems involving frequency where c is the speed of light and the chiral parameter coupling the electric and magnetic fields. The Lorentzian model is employed to determine the variation of the permittivity and the permeability with frequency. Here, the subscript ∞ denotes high frequency values, while the subscript refers to the low frequency limit. The values e , h refer to the resonance frequency and e , h , to the damping coefficients. The chiral parameter is obtained, from the Condon model, as where k denotes the resonance frequency, while k corresponds to the damping coefficient for the Chiral material. The coupling constant, k , defines the magnitude of the chirality. Using these relations, Eqs. (56) and (57) may then be expressed as and where

Discrete Equations and the Z-Transform
To convert Eqs. (63) and (64) from the frequency to the time domain, we use the substitution j → t and then apply the FDTD time stepping scheme to discretize the resulting equations. The bilinear transformation is used, to convert a frequency dependent function to the Z-Domain [42]. Equations (65)  (63) is required on the Voronoi edges, but it includes the electric field projections, that are only available on the Delaunay edges. The solution to this problem is illustrated in Fig. 39. For a tetrahedron, we have six Delaunay edges and the electric field projection E i , i = 1, … , 6 is stored on each edge. Using these projections, the electric field vector Cell 1 can be reconstructed and linked to the centre of the corresponding cell, in this case cell 1. To achieve this, the equation set �� is constructed, where we again use the procedure described earlier. The resulting field vector is associated to the corresponding cell center. As a Voronoi edge connects two cells, we need to repeat this procedure for the second cell leading to Cell 2 . The field vectors linked to these cells are then averaged and projected on to the corresponding Voronoi edge according to ( Cell1 + Cell2 )∕2 ⋅ , as illustrated in Fig. 39b.
With the electric field defined on the Voronoi edges in this manner, we may then compute J n e,i av . A similar process is adopted for the magnetic field, but in this case we have a point from the Delaunay mesh inside a Voronoi cell [32].

Solution Update Process
To update the solution over one timestep, the total field formulation is employed as an intermediate step. In this manner, the solution is updated over one time step as follows:  39 Determining the electric field on the Voronoi edges: a reconstruct the electric field vector from its projections on the Delaunay edges and link this vector to the corresponding cell center; b aver-age the electric field vectors from the cell centers linked to the same Voronoi edge and project them to the corresponding edge

Material Interfaces
Boundary conditions at material interfaces are applied by using a weighted average process that is similar to that employed earlier for simulations involving isotropic and anisotropic media. However, for the Padé approximants, which are required in expressions such as Eq. (68), we only average the coefficients in the numerator. To illustrate this process, consider two tetrahedral cells, separated by a boundary face, as illustrated in Fig. 40, and consider the calculation of the permeability for the corresponding Voronoi edge. Applying the bilinear transformation to the frequency dependent component ∞ − , from Eq. (61), leads to ( ∞ − P(z)) , where The magnetic flux is determined as where the subscript i(= 1, 2) refers to the medium 1 or 2. The average magnetic flux at the interface is calculated as where w i = l i ∕(l 1 + l 2 ) for i = 1, 2 and l i corresponds to the length of the segment of the Voronoi edge between the intersection with the Delaunay face and the center of the tetrahedron C i as shown in Fig. 40.

Rotation of the Plane of Polarization
To verify the implementation of the scheme, the transmission coefficients for an electromagnetic wave incident upon a chiral slab are computed, along with an evaluation of the variation of the angle of the plane of polarisation of the incident field. The dielectric slab has dimensions 0.3 m, 0.3 m and 0.1 m in the x, y, and z directions respectively and the mesh employed is such that 40 cells are used to discretise a length of 0.1 m. The electric permittivity and magnetic permeability of the slab are described by a Lorentz model and the chirality by a Condon model. The slab is located in free space, as depicted in Fig. 41. The slab has a real chirality, which induces an ORD of the incident field. In this case, the variation of the angle of polarisation, , can be computed as where is the chirality, is the angular frequency of the incident wave and L is the thickness of the slab [47]. To compute this angle, we use a Fourier transform, which allows us to evaluate the variation of the amplitude without time dependence. A point in the Delaunay mesh is selected, such that, with respect to the incident wave, it lies behind the chiral slab. Then, for a wave propagating in z direction, the numerical value for the angle of polarisation is computed as  Figure 42 shows the numerically computed rotation compared to the theoretical result, for different values of the chirality. For values of less than 0.00407, the numerically computed angles of rotation of the plane of polarisation remain constant over time and are in good agreement with the theoretical values. For higher chiralities, however, the numerical solution shows very strong fluctuations indicating an instability of the scheme. We performed a stability analysis with respect to variations in the time step, the spatial step and the real and imaginary part of the chirality [48]. The results showed that the stability is independent of the time step size. To investigate the effect of the spatial step size on the stability, a rotation per cell is defined by expressing the slab thickness, L, as L = n z , where n is the number of cells through the slab, in the z direction, and z is the edge length of a given cell. The cell normalised rotation of the plane of polarization is then defined as and the relative error in this quantity is computed. For a given real chirality, reducing the mesh size leads to a more stable scheme. A similar result is obtained if we use a given mesh and vary the chirality. Typically the lower the real part of the chirality the more stable the scheme. For a stable simulation and a real chirality, the results suggest that the scheme remains stable if cell ≤ 0.4 • . An imaginary component of the chirality also has a stabilizing effect on the scheme, as can be seen in Fig. 43. A coupling coefficient k = 2 ⋅ 10 −12 s leads to strong instabilities in the scheme without damping. We increase the damping, using k = 0.01, 0.025, 0.5, 0.1 in turn, and again plot the angles over 20 cycles. From the results, shown in Fig. 43, we see that the scheme is now stable well above the Cell = 0.4 • limit suggested for the non-damped case. It appears that, provided Im( ) ≥ Re( )∕15 , the scheme is stable, but, in this case, we are no longer dealing with a pure ORD.

Reflection/Transmission on a Chiral Slab
We consider the problem of simulating the interaction between an incident wave and a slab of bi-isotropic chiral material, with the object of determining the transmission and reflection coefficients. The chiral slab, which is assumed to be of infinite length and width, with a finite thickness L, fills the region 0 ≤ z ≤ L . The regions −∞ ≤ z ≤ 0 and L ≤ z ≤ ∞ are filled with free space. When the incident wave takes the form of a linearly polarised narrowband pulse, as in Eq. (21), propagating in the z direction, the problem has an analytic solution [47].
In the numerical simulation, the thickness of the slab L = 0.1 m is discretised with 40 cells and the chiral material is characterised by the values s = 1.8 0 , ∞ = 1.8 0 , s = 1.1 0 , ∞ = 1.0 0 , e = h = k = 2 × 3.5 G H z , e = 0.14 , h = 0.12 ,   cycles Φ/n=0.847˚,ξ k =0.1, κ=-0.0567-0.0157i Φ/n=0.895˚,ξ k =0.05, κ=-0.0599-0.0083i Φ/n=0.908˚,ξ k =0.025, κ=-0.0608-0.0042i Φ/n=0.911˚,ξ k =0.01, κ=-0.0610-0.0017i Fig. 43 for different damping factors k = 0.01, 0.025, 0.5, 0.1 , the relative error of the computed angle over 20 cycles is plotted k = 0.1 and k = 1 ps. A point p1, in the free space region just ahead of the slab, is chosen for the computation of the reflection coefficient, while the transmission coefficient is calculated at a point p2, in the free space region slightly behind the slab. The co-polarized reflection coefficient, R co , the co-polarized transmission coefficient, T co , and the cross polarized transmission coefficient, T cr , are then determined as For a simple isotropic material, T cr is equal to zero. Good agreement is obtained in the comparison between the analytical and computed distributions of the transmission coefficients shown in Figs. 44 and 45. The differences in the results arise because the analytical model assumes a slab which is infinite in both the x and y directions and, in addition, the numerical dispersion inherent in the scattered field formulation leads to a shift in the peaks for the reflection coefficient.

Chiral Sphere
This example involves the simulation of the interaction between a plane wave and a chiral sphere. We consider a GHz , e = 0.12 , h = 0.12 are employed for the material parameters. A time step size t = 6.18 ps is used and the numerical solution is advanced through 20 cycles of the incident wave. The results produced for the RCS are compared to those produced in an analytical solution of this problem, using the program developed by Demir et al. [49].
The analytical and numerical RCS distributions are compared in Figs. 46, 47, 48, and 49 where, due to the symmetry of the RCS distributions, only the results for angles lying between 180 • and 360 • are shown. All the results demonstrate a very good agreement between the numerical and analytical solution for different values of the coupling constants k and the wavelength 0 . In particular, the very weak signal at high angles, from 300 • to 360 • , is well captured.

Conclusion
We have demonstrated an unstructured mesh implementation of the classical Yee FDTD co-volume algorithm for computational electromagnetics. The complexities associated with generating an appropriate primal Delaunay mesh, together with its Voronoi dual, have been identified and addressed. The approach was used initially for the simulation of wave scattering by PEC and isotropic dielectric objects, where the efficiency of the solution process, and the accuracy of the results, has been demonstrated. An advantage of the approach is that,with the same set of equations, it allows for the modelling of isotropic lossy dielectric materials or PECs by simply increasing the conductivities to very high values. This does not lead to any instabilities in the scheme.A description has been given of the extension of the approach, to enable the modelling of problems involving anisotropic lossy materials and of problems involving frequency dependent bi-isotropic materials, with coupling between electric and magnetic fields. The incorporation of frequency dependent materials was accomplished by adopting the Z-Transform method. This solution procedure must now be parallelised, so that its performance may be fully compared with that of other numerical techniques when applied to problems that are computationally more challenging. When this has been accomplished, the performance of the method can then be compared with that of the powerful current generation of hybrid finite difference time domain/finite element time domain methods [50][51][52][53].