Extracting modal field profiles from 3D unstructured transmission line modelling meshes for use as sources and observers

: Large-scale three-dimensional (3D) electromagnetic simulations are commonly used as design and investigation tools in a wide variety of technological fields. It is not uncommon for both excitations and observation quantities to be expressed in terms of particular field profiles of feed waveguides. These may then be used to evaluate, e.g. scattering parameters. These field profiles must be obtained as a pre-processing task before the main simulation. Use of a theoretical field profile as an excitation to a discretised structure will typically cause a non-physical reflection. It is therefore more desirable in practice to use a field profile that is consistent with both the discretisation of the geometry and the 3D method of simulation. The authors present an approach to extracting these 2D field profiles from large-scale 3D unstructured meshes which are to be simulated with the unstructured transmission line modelling method. Discretised slices from the 3D mesh are extracted and incrementally extruded into a form suitable for consistent pre-processing. The impacts of all the parameters of the approach are investigated. Benchmarking is undertaken on both coaxial cable and microstrip waveguide feed structures showing that good quality results can be obtained straightforwardly.


Introduction
Large-scale numerical computational electromagnetic methods have rapidly evolved to become a mainstay of design and simulation activities in a broad range of disciplines [1][2][3][4]. In the authors' opinion, experienced users understand and accept that the accuracy of the results from simulations are dependent upon the refinement of their discretised representations. However, there is a difference between the perturbation of values such as the scattering parameters and observing a non-physical phenomenon such as a spurious reflection. Therefore, as simulation developers, it is preferable where possible, to carefully manage how the consequences of discretisation errors manifest themselves to the user.
Besides the core algorithms and their computationally efficient implementation, practical deployment of numerical methods also demands a supporting portfolio of techniques to provide a complete simulation capability. For example, it is common to develop radiation boundary conditions, thin panel, thin wire models and complex material models. Moreover, appropriate means for extracting observation quantities of interest and exciting the simulation in a physically meaningful manner must be available. In this work, we demonstrate how to effectively extract particular two-dimensional (2D) field profiles from a 3D meshed geometry which are necessary to specify a class of excitations and observations when using the unstructured transmission line modelling (UTLM) method [5]. In particular, we show how to remove a non-physical consequence of discretisation errors from the point of view of a user. Fig. 1 shows a pair of illustrative scenarios which we have used to describe and benchmark the method presented below. First, Fig. 1a shows a compact antenna excited by the fundamental TEM mode of a coaxial feed line and its performance is characterised in part by the return loss, S 11 [6]. Second, Fig. 1b shows a Lange coupler which is excited at one of its ports by a fundamental microstrip mode and its performance is characterised in terms of a four-port scattering matrix [7,8]. In both cases, the 3D simulation is being performed using an algorithm that discretises the problem space using an unstructured tetrahedral mesh. (Such 3D meshes may be hybridised for efficiency with a cubic grid in the relatively large free space regions.) A pre-processing task before the main time-stepping simulation is to obtain suitable 2D field profiles for the modes of the coaxial and microstrip lines. These are to be used to; (i) provide the source excitation to the simulation and, (ii) to evaluate the amplitude of the scattered modal fields for the purposes of identifying the circuit scattering parameters. The modes being sought are 2D waveguide-like eigensolutions of Maxwell's equations. For the coaxial cable, there exists a simple closed form expression for this field. For the microstrip line, a numerical mode solver would often be employed beforehand for this purpose. The use of an independent mode solver is at least inconvenient, but more importantly, may discretise the fields in a different manner to the representation used by the 3D algorithm requiring the use of interpolation. In both cases, it is important to recognise that solutions for the modes obtained using different spatial discretisations and/or different approximations to Maxwell's equations will be numerically different.
The key motivation for this work is remove a non-physical artefact from 3D simulations that compromises the usefulness of their results for design and understanding of the phenomena involved.
If, for example an excitation is imposed in the form of a field profile which is not a waveguide-like eigensolution of Maxwell's equations discretised in exactly the same manner as that used by the 3D algorithm, then a numerical impedance mismatch will occur. Specifically, both the characteristic impedance and the field distributions will be slightly different. This leads to non-physical reflections. The simulation user could then be confronted by a situation where exciting a straight section of supposedly impedance matched waveguide by its fundamental mode actually generates a spurious return loss. The authors believe that, even if an exact closed form expression for the waveguide mode is known, it is more valuable to use the closest approximation to it that does not cause such a reflection. The error associated with the approximation is then simply a further small contribution to the background discretisation error expected by a user. Fig. 2 presents a summary of the context. A 2D slice through a 3D mesh is identified as an excitation or observation port. This slice will have been discretised in the 3D mesh. A modal field profile is to be evaluated on the slice such that if the slice were extruded to form a An important point to clarify is the application context motivating this study. There is a long history of developing, refining and implementing specialised techniques for electromagnetic mode solving techniques in both the microwave and optical regimes that have been proved highly successful (illustrated by chapter 7 of [4]). However, in the context of large 3D simulations, for example on the scale typical of aerospace electromagnetic compatibility studies, it is unlikely that the geometry around feed waveguides will be meshed very finely. Simply demanding a finer discretisation of the feeds is not a viable option due to the overall demands on resources.
Hence, the key measures of success for this work are proposed as: (i) To obtain good representations of field solutions for feed structures which have relatively crude discretisations. (ii) That the field solutions are obtained with a method directly compatible with the technique used to perform the full 3D study in order to avoid spurious reflections.

UTLM method
The transmission line modelling (TLM) method was initially developed in the early 1970s, from a 2D through to a 3D algorithm, for predicting electromagnetic field behaviour [2,3,9]. Originally operating with a Cartesian discretisation on a cuboidal grid, more recently a powerful variant that uses unstructured tetrahedral meshes has been developed [5]. This has been thoroughly characterised and industrially deployed, e.g. [10][11][12][13]. Unstructured meshes offer both smoother material and geometrical boundary representations than their Cartesian counterparts. Furthermore, they easily permit the use of a highly graded meshes to deal with multiscale problems which otherwise demand the use of multi-grid methods to maintain computational efficiency [10,[14][15][16][17]. TLM algorithms were originally developed by making an analogy between the field solutions of Maxwell's equations and the behaviour of voltages and currents on an interconnected electrical network of commensurate transmission lines. In the unstructured case, it is more simple to regard each tetrahedral cell as an incremental volume of compact homogeneous space in which local canonical solutions to Maxwell's equations are available as interpolators between discrete field sampling points defined on the surfaces of each cell. Requiring continuity of the fields across the boundaries of adjacent cells yields a time-domain algorithm for repeatedly updating the field samples as time evolves. Note that this is the reverse of the steps used to derive finite element (FE) methods which also start with a local field interpolation for each cell, but which are additionally required to a priori satisfy suitable continuity conditions. The FE update algorithm then follows from the requirement that the total field is then a solution of Maxwell's equations in a variational sense. Finally, the discretisation approaches of both TLM and FE methods contrast with the dual grid (Delaunay and Voronoi) arrangement used by finite-difference time-domain (FDTD) techniques. In the UTLM case, the feature that maintains its link with the original Cartesian TLM is its implementation, rather than its derivation, in terms of commensurate lengths of transmission lines. Advantageous features of the TLM algorithms are (i) a spatial and temporal colocation of the electric and magnetic field samples, unlike with FE time-domain and FDTD and most importantly, (ii) a priori assurance of stability. Stability is provable on a cell-by-cell basis without the need to rely on guidelines such as the Courant condition or else an impractical eigenanalysis of the whole mesh. For large-scale simulations, this is a critical advantage and late time instability has never been observed. One further enabler is the use of automatic cell clustering. In common with all explicit time stepping algorithms, the choice of the UTLM time step is often dictated by the smallest mesh feature, i.e. cell, present in the simulation. Mesh generation is not able to guarantee that small cells are only present where absolutely necessary to represent the geometry. Indeed, even in the presence of fine geometry, one usually wishes to set a notably larger time step based upon the physically frequencies of interest. Cell clustering is a technique whereby many thousands of TLM tetrahedral cells are automatically coalesced into larger entities [10]. A simple preprocessing of the characteristics of these clusters inherently separates their behaviour into that which is physically meaningful and that which is attributable to noise caused by the non-regular distribution of the mesh across the cluster. The latter can be explicitly excluded from the simulation without compromising the physical validity; the only effect being to change the particular nature, but not the magnitude, of the dispersion errors. Overall, cell clustering, which may be regarded as constructing a higher order scheme on the fly, allows significantly larger time steps to be used which reduce the time-stepping run time.

Extracting modal profiles
This section presents the core method used to achieve the objective of the work. The principles of the approach are drawn from established background material and the novel contribution is to apply them for the first time to the particular discretisation found in 3D UTLM simulations. Fig. 3 shows the field sampling on the surface of a single UTLM cell and how this is implemented using propagating voltage impulses on a transmission line network. (The complete details are given in [5,10].) Referring to Fig. 2, a surface triangulation is extracted from the 3D mesh upon which a 2D modal field profile is to be evaluated using the same sampled field representation as in the 3D mesh. Consistent with the formulation of a propagating mode eigenproblem, a stand-alone prismatic mesh is extruded from this surface triangulation and tetrahedralised. It is noted that the tetrahedra appear in degenerate (i.e. have the same circumcentre) groups of three, each of which forms a small prismatic cluster. Fig. 4 shows the sampled fields on such a structure, noting that there are three sets, those on the top and the bottom, which are tangential to the prism cross-section and those midway between. On the sides of the prism a suitable impedance or more complex boundary condition is to be used to relate the electric and magnetic  fields and this is typically drawn from the geometry of the 3D problem.
For each UTLM scattering cell indexed by n, a scattering matrix, S n is defined by the UTLM algorithm such that V nr = S n e − jωΔ t V ni (1) where the term e − jωΔ t embodies the time stepping delay of the algorithm in the Fourier domain and the superscripts nr and ni denote the reflected and scattered quantities of the nth cell, respectively. Combining the expressions (1) from each cell gives an overall expression for all cells A sparse connection matrix, C enforces field continuity between adjacent cells, the boundary conditions on the sides of problem and the reflection coefficients at the end of the transmission line stubs that form part of the cell scattering elements, S, as shown in Fig. 3 V The subscript m, t and b denote the samples on the midpoint, top and bottom of the prisms shown in Fig. 4b, respectively.
Combining (2) and (3) yields an expression which is conveniently expressed in the form The z-directed modal solutions being sought are defined by augmenting the relationships for the midpoint samples, (3), with for the top and bottom samples, where z is the direction of extrusion and Γ = e − jβΔ z (6) β is a propagation constant, which may be complex in order to account for attenuation, and Δ z the extrusion length of the prismatic mesh (see Fig. 2). From (4) and (5) which can be rearranged to present a generalised eigenvalue problem for Γ.
Although symmetry, losslessness and reciprocity will provide a number of relationships regarding the structure of the matrices A and B, in this work we do not manipulate (8) further in order to improve its numerical characteristics for such special cases. Furthermore, the dependence upon frequency does not appear in a manner amenable to its explicit extraction. Equation (8) can be solved using the shifted inverse power method based upon an estimate of the value of β for the solution being sought which is usually easy to provide in context. The matrices are sparse and hence the complex linear equation solutions required are obtained using general minimum residual (GMRES) employing a sparse incomplete LU threshold (ILUT) preconditioner [18]. Once solved, the sampled electric and magnetic fields, E and H, can be processed to yield a value for the modal wave impedance, where meaningful, using, e.g.
where the integrals are evaluated over the cross-section. This corresponds to use of a power-voltage definition of impedance [19] 4 Results

Coaxial cable
Fig . 5 shows the convergence with mesh size of the approach for a coaxial cable operating at 3 GHz. The cable has core and shield radii of 1 and 3.495 mm, respectively, and an insulator relative permittivity of 2.25, [6]. Both the core and the shield are assumed to be perfect conductors. As described above, a triangulated crosssection mesh was extracted from a 3D mesh (produced using an inhouse 3D Delaunay tetrahedral mesh generator) and extruded to provide a uniform prismatic slice. Meshes were generated with two measures of quality, specifying (i) a maximum area and (ii) a maximum Q-factor for the cross-section's triangles. The Q-factor for triangles is defined as the circumradius to minimum edge length ratio and provides a convenient quantification of the shape of the triangle [20]. Equilateral triangles have a Q-factor of 0.58 and generally speaking, a value of <2 is regarded as indicative of well-shaped triangles [20]. Further degrees of freedom are the extrusion length, Δ z , and the time step Δ t . Except where stated, the value of Δ t is set by the lowest prismatic cell cluster response time, as described in [10], such that in all cases at 3 GHz, ωΔ t ≪ 1. The use of small values of Δ t does not incur a run time penalty for mode identification as it does for time stepping. However, Fig. 5 shows one example of the effect of explicitly controlling the value of Δ t . Fig. 5a demonstrates that the errors in both the propagation constant, in terms of an effective relative permittivity ϵ eff = c o β/ω and the modal wave impedance, Z w , exhibit second-order convergence with triangle size. c o and ω are the free space speed of light and angular frequency. The triangular element size is quantified in terms of the square root of the maximum triangle area, denoted by A max , in order to provide a measure which is dimensionally a length. Fig. 5b considers the effect of the extrusion length, Δ z . It can be seen that the propagation constant remains accurate even as Δ z approaches the length scale of the triangles and in fact, remains constant until the sampling limit with respect to the propagation constant is reached. Typically, the value for Δ t is set automatically as described above, but it was also argued in the introduction that being able to use exactly the same parameters as for a subsequent 3D simulation is important to avoid numerical impedance mismatches. Fig. 5c shows that the convergence with the scaled measure, 2c o Δ t can conservatively be set using a criterion of 10% of A max .
The formation of higher order cell clusters has proved a powerful enabling technique for the UTLM method. Without clustering, the time step can be required to be significantly smaller than demanded by the physics, incurring a significant computational penalty when time stepping. The clusters are formed by specifying a threshold distance, D T , and cells whose scattering centres are closer than this are coalesced as described in [10]. In Fig. 5d, the use of larger thresholds is clearly seen to allow larger time steps to be used. The limit on the thresholding occurs when too many cells (typically >1000 in a 3D mesh) coalesce such that  the pre-processing of the cell clusters becomes prohibitively large. The rise in error as larger clusters are formed in Fig. 5d is due to the fact that the spatial resolution in the transverse plane eventually becomes diminished. Fig. 6 assesses the sensitivity of the approach to mesh quality: Good quality triangles are quasi-equilateral (Q-factor <2) and poor triangles are those possessing either a small (<10°) or a large (>170°) angle [20,21]. The Delaunay mesh refinement process recursively splits poor triangles at their circumcentre to achieve either better quality or smaller elements, so it is actually difficult to generate small element area, poor quality, meshes as splitting for area naturally improves the quality as well. Fig. 6 shows four cases with different worse case element qualities. It is remarkable that even with such crude meshes, reasonable results can be obtained. For comparison, Fig. 6b shows a field profile obtained using a more refined mesh and it is noted that one can see that the circular cable has been discretised by a polygon (of 16 sides) which supports the proposition that such geometrical resolution will not necessarily be high in very large-scale simulations.

Microstrip line
The second example presented is driven by the need to provide a source to and observe the outputs from the Lange coupler example shown in Fig. 1b. The relevant parameters for mode extraction are an operating frequency of 1 GHz, a substrate relative permittivity of 9.6 and thickness of 2 mils (0.0508 mm) and a strip with and a thickness of 10 mils (0.254 mm) and 2 mils (0.0508 mm), respectively [7]. For each of the microstrip ports, a small 'selecting cube' is used to extract a part of the surface triangulation from the 3D mesh from which an extrusion is generated as shown in Fig. 2. In general, this may encompass known boundary conditions, e.g. the ground plane here, but also leave (jagged) edges on which a boundary condition must be imposed for the purpose of mode identification. In this case, an open-circuit condition is used to close the top and sides of the problem creating a boxed microstrip structure. Fig. 7a examines the convergence of the propagation constant with mesh size. The reference value taken is the result obtained using a fine mesh ( A max = 0.003 mm). Again, second-order convergence is observed. Fig. 7b shows a selection of field profiles obtained. In the context of exciting very large-scale 3D simulations, it is useful that an error of ∼6% is available even with the crudest mesh shown.
Finally, it is commented that the computational effort required to obtain the results is wholly dependent upon the effective implementation of a sparse preconditioned matrix eigensolver. In the above examples, results were obtained in a matter of seconds when run serially on a Sandy Bridge CPU platform and never required more than a minute. This is to be compared to the subsequent 3D simulation times often extending to many hours or days.

Conclusion
This paper has presented for the first time a methodology for extracting waveguide mode parameters using the 3D UTLM approach. The value of the work to the simulation developer and user community is to be able to provide field profiles on 2D slices that can be used for the purpose of excitation or observation which   do not give rise to spurious non-physical reflections. If such artefacts are not avoided where possible, the interpretation of simulation results is less straightforward. Slices containing triangular elements are taken from a 3D tetrahedral mesh and extruded to produce a prismatic structure. This structure is then analysed to identify the eigensolutions corresponding to the propagating modes. The analysis is developed by applying well known waveguide and matrix theory to the particular discretisation used by the UTLM algorithm.
The method demonstrates second-order convergence with mesh size and an investigation of the other degrees of freedom has shown that the method is robust. The key contribution is considered to be that the solutions obtained are (i) fully consistent with both the field sampling method of the full 3D mesh upon which they are intended to interact and (ii) able to provide acceptable approximations to the true mode profiles even when the mesh is crude which is important in the context of resource limited very large-scale simulations.