- Email: [email protected]
Multi-Scale Modelling in Time-Domain Electromagnetics
© Urban & Fischer Verlag http://www.urbanfischer.de/journals/aeue
Multi-Scale Modelling in Time-Domain Electromagnetics Christos Christopoulos Dedicated to Professor Peter Russer on the occasion of his 60th birthday Abstract: The paper reviews the techniques used to describe efficiently in time-domain electromagnetics the presence of fine physical features in an otherwise coarse numerical mesh. This is a major area of concern in electromagnetic modelling and simulation as it affects crucially modelling capabilities in dealing with complex practical configurations. Two broad classes of modelling methods are identified, those where the basic approach is to distort and/or refine the numerical mesh, or, interface the mesh to local solutions. A typical example of a technique in the former category is multi-gridding and in the latter mesh-to-wire interfaces. The emphasis of the paper is to describe the essence of the different approaches, their strengths and weaknesses and to identify areas of further work. Keywords: Computational electromagnetics, Multi-scale modelling, Time-domain simulation
1. Introduction The modelling of large-scale electromagnetic problems involves several stages. It assumes knowledge of the basic system configuration and components to proceed in establishing a numerical model which can be computed to offer the required outputs. It is relatively straightforward to devise efficient models in cases where all physical dimensions are comparable in electrical size (size relative to the wavelength). In such cases a numerical mesh with the appropriate spatial resolution can be devised to model the entire problem space. Several generic numerical techniques are available for general electromagnetic modelling including, Finite Difference TimeDomain (FDTD) [1], Finite Element Method (FEM) [2], Method of Moments (MoM) [3] and Transmission-Line Modelling (TLM) [4]. However, practical problems often do not conform to this requirement. This is particularly true in Electromagnetic Compatibility (EMC) problems where it is very common to have, in the same problem, very different physical scales (relative to the wavelength). A typical example is a wire inside an equipment cabinet. The wire diameter is typically much smaller than the wavelength while the cabinet size may be comparable to the wavelength. This and similar problems are described
Received September 1, 2002. Revised January 2, 2003. C. Christopoulos, Electromagnetics Research Group, University of Nottingham, Nottingham, NG7 2RD, U.K., Tel. +44 115 951 5557, Fax: +44 115 951 5616. E-mail: [email protected] ¨ 57 (2003) No. 2, 100−110 Int. J. Electron. Commun. (AEU)
as multi-scale problems and present severe difficulties to the modeller. It has to be stated that no amount of increase in the size and power of computers is likely to overcome this problem as for many systems full explicit knowledge of all details is impractical or impossible even if the computational power were available. An additional consideration, which is often overlooked, is the advantage of keeping models as simple as possible, thus offering intellectual economy and scope for creativity. Therefore for practical and methodological reasons it is essential that, in a large problem, features of a very small electrical scale or complex details are dealt with separately from the rest of the problem which is on a different scale. This requirement applies more or less irrespective of the actual numerical method used for field simulation be it FEM, FDTD, TLM etc. In tackling this problem two broad approaches have emerged. First, the mesh is refined or distorted locally to increase spatial resolution in areas where fine features are present and therefore fields may vary rapidly in space. In this way, increased resolution is applied only where it is needed thus maintaining computational efficiency. Difficulties however remain at the interface between meshes of different resolution and/or shape as will be discussed in the following sections. Second, the mesh is left unchanged at a more or less uniform resolution set with global considerations in mind. A local solution is obtained for the fields near the fine feature and this then is interfaced to the mesh. The difficulties here are both in obtaining local solutions and in devising appropriate interfaces. Both these approaches will be described in some detail in the following sections. As regards the types of fine feature that may be encountered, one should include, fine wires, wire looms, thin panels with/without perforations, fine conducting tracks, entire assemblies of such components etc. The objective is to present what has been done so far and to stimulate further research work in this area.
2. Mesh refinement and distortion A regular Cartesian mesh is shown in Fig. 1. Let us consider that there is a fine wire of radius a as shown in crosssection. In a mesh with spatial resolution ∆ and when the radius of the wire is smaller than the mesh size, the impact of the wire as regards its diameter and actual position cannot be accurately represented. One option is to reduce ∆ until it matches the radius of the wire. The however entails a heavy computational cost as an unnecessarily fine 1434-8411/03/57/2-100 $15.00/0
C. Christopoulos: Multi-Scale Modelling in Time-Domain Electromagnetics 101
Fig. 1. Schematic of various mesh arrangements: (a) regular mesh, (b) multigrid mesh and (c) hybrid mesh.
mesh is used in areas where it is not needed. Halving ∆ in a three-dimensional (3D) mesh will result to an increase in storage by a factor of eight. In practical problems a factor of two is rarely sufficient and very often an increase in resolution by an order of magnitude is necessary. Similar unacceptable increases are also inevitable in run-time making it impossible to tackle in this way any but the simplest problems. One possible way out of the dilemma is to refine the mesh as shown in Fig. 1b. This is known as a multigrid, multiple grid mesh or subgrid. The other option is shown in Fig. 1c and is described as a hybrid mesh. Both examples are for structured meshes of the type used in the main time-domain simulation techniques such as FDTD and TLM. When unstructured meshes are used (e.g. in FEM) there is greater flexibility in configuring the mesh to accommodate fine features but such techniques have not found favour in the time-domain. We therefore limit the discussion here to structured meshes. Each approach is examined separately. The obvious advantages of the scheme in Fig. 1b are that mesh refinement is implemented only where is needed, and that no other model developments are necessary to describe the fine features. Its perfectly conducting surface is registered in the mesh as zero tangential electric field (eg in FDTD) or as a boundary with a reflection coefficient of −1 (eg in TLM). There are however disadvantages and difficulties. Amongst the disadvantages one could mention the fact that the simplicity of the single regular mesh is lost and that each fine feature must be laboriously mapped onto the mesh. The difficulties are more profound. At the interface between the fine and coarse mesh regions information is available at different spatial and temporal rates. This therefore implies some form of averaging/filtering at the interface. What are the physical principles that must be enforced in order to construct such an interface? To-date several solutions to this problem have been proposed. To illustrate the complexity of the interface, let us consider the example shown in Fig. 1b. For this example, the spatial reduction ratio in 3D is 16:1 and over the time step of the coarse mesh there are 4 × 16 pulses from the fine mesh. Even higher reduction ratios may be needed in practice illustrating the complexity of the required interface. One of the earliest attempts to introduce local grids was described in [5] where, however, two separate FDTD simulations were required. The simultaneous coupling at the coarse-fine mesh interface, thus requiring only one FDTD solution, was described in [6]. Here fields are
passed onto the local grid after time and space interpolation. The authors reported two difficulties. First, the abrupt change in mesh size causes spurious reflections. Second, since a wave experiences a different phase error in each mesh, non-physical interactions may take place and therefore they advise that the coarse mesh size is still kept much smaller than the shortest wavelength of interest. A subgridding scheme was also described in [7], referred to as the Variable Step-Size Method (VSSM). A more efficient version known as the Mesh Refinement Algorithm (MRA) was described in [8]. It is shown there that the reflection coefficient at the coarse-fine mesh interface is less than 1%. However, in both these papers the interface must be placed in a homogeneous region. An alternative scheme is described in [9]. The authors have tried several algorithms and interpolation techniques and found that a pulsing overlapping scheme of computations yields the best results. Their interface is usually in homogeneous regions where a reflection coefficient below 60 dB can be expected. They state that if the interface is placed in inhomogeneous regions, somewhat higher reflections can be expected. Placing the interface in inhomogeneous regions was described in [10]. Here the refinement ratio must be an odd number and only the magnetic field is interpolated. Reflections at the interface of less than 1% are reported and no instabilities, although no formal demonstration of stability is offered. Results are shown for a refinement ratio of 3:1, although it is stated that higher ratios (5:1, 7:1) were tried successfully. It was reported in [11] that both the VSSM and MRA techniques suffer from long time instabilities. The same authors have developed an equivalent circuit representation to ensure a stable algorithm [12]. A similar algorithm is also described in [13]. A multigrid current method (MGCM) technique is described in [14], suitable for placing the interface in an inhomogeneous region using any value of the refinement factor. It is stated that the reflection coefficient at the interface is less than 2% and that there is improved stability. The wealth of papers on this topic over a number of years illustrates the inherent difficulties. Parallel developments to those described above in connection with FDTD have also taken place in TLM. Multigrid schemes were proposed in [4, 15, 16]. Since TLM is an equivalent circuit field representation, provided that the conversion across the coarse-fine interface is consistent with circuit theory, stability is guaranteed. The conditions that conversion must satisfy are: charge and energy conservation, no reflections and no delays. It was found that strict adherence to all four conditions was not possible and therefore a compromise is necessary. Stable schemes can be devised but at a slight energy loss at high frequencies. A TLM multigrid scheme that has been found acceptable is based on the number of transmission lines and time steps in each region of the mesh (fine and coarse). The notation for a region ‘i’ is as follows: • si , pi number of lines across the two polarization directions • ti number of time steps before conversion
102 C. Christopoulos: Multi-Scale Modelling in Time-Domain Electromagnetics • n i = si pi ti total number of pulses before conversion Then, in converting pulses from region 2 into region 1, the new pulses to be inserted into region 1 are given by
V1 j =
n2 1 V2k s1 p2 t2
, j = 1, 2 . . . n 1
(1)
k=1
where, V2k are the pulses in region 2 travelling towards the interface. All inserted pulses are taken to be the same to avoid introducing high frequencies. Similarly, for pulses inserted into region 2,
V2 j =
n1 1 V1k s2 p1 t1 k=1
, j = 1, 2 . . . n 2
(2)
where, V1k are pulses in region 1 travelling towards the interface. For the example shown in Fig. 2, s1 = 1, p1 = 1, t1 = 1, n 1 = 1, for the coarse mesh 1 s2 = 2, p2 = 2, t2 = 2, n 2 = 8, for the fine mesh 2 Hence, from equations (1) (2),
V1 =
1 V2 j = V1 2
1 V2k 4 k=1 8
, j = 1, 2 . . . 8
(3)
(4)
An example of a multigrid application is shown in [17]. An improved multigrid interface was described in [18]. Undoubtedly, multigrid schemes offer significant savings (several orders of magnitude) in computation, and in many cases can make the difference between being able to simulate a problem and abandoning the effort on the grounds of insufficient computational power. However, stability is always an active issue and care must be taken in the way these methods are applied. Irrespective of these benefits, it must be questioned whether it is the correct approach to refine the mesh indefinitely, to cope with finer details, rather than seek different approaches based
Fig. 2. Schematic showing conditions at the interface between coarse and fine mesh regions.
on local solutions and/or signal processing techniques, to couple such features to the surrounding coarse mesh. This issue is tackled in the next section. An alternative approach to the multigrid is the hybrid mesh depicted in Fig(1c). Here the nodes do not represent a cube of space, but can be shaped according to the required resolution in each direction. A full theory for this mesh is described in [4] where references to original papers can be found. Various schemes are available with different time step requirements as discussed in [19]. It is difficult to assess the loss in accuracy for high grading ratios. Comparative studies showing mesh behaviour for abrupt changes in grading for TLM and FDTD show that the former exhibits more robust behaviour [20]. The TLM hybrid mesh is unconditionally stable and is used extensively in commercial solvers. As an introduction to the next section the approach of making a local mesh modification is mentioned here. The problem considered is the accurate representation of edges and corners. Near such features fields exhibit singularities and variation which is not linear [21, 22]. However, the sampling of fields in FDTD and TLM is implicitly based on a linear field variation. This inconsistency results in extra energy storage leading to lower frequency estimates at resonance. This problem manifests itself in TLM by longer propagation delays near edges. Various approaches have been proposed to deal with this so called coarseness error. Naturally, the error may be reduced by a multigrid mesh with fine resolution near the edge, or by a hybrid mesh. However, modellers are bound to seek simpler, more efficient, solutions. One approach, described in [23], is to speed up propagation on the nodes immediately adjacent to the edge. A more efficient way of achieving the same objective is described in [24].
3. Mesh interface to local solutions We start in this section by revisiting the edge singularity problem mentioned at the end of the previous section as it illustrates the interfacing of local solutions to the mesh. A more systematic and rigorous approach is described in [25]. Here the field near the edge is obtained from quasi-static approximations and is used to derive an equivalent circuit representing the fields at the edge. This circuit, which in general 2D problems is a threeport, is then connected to the three TLM nodes which interact directly with the edge as shown in Fig. 3. The lumped components in this circuit are obtained from the Z-matrix representation which in turn is obtained from the quasi-static fields. Field singularities near the edge have also been dealt with using wavelets as described in [26]. Work along similar lines has also been done in connection with FDTD as discussed in [27–32]. Published work is limited in two-dimensions. For a discussion of some of the issues involved in 3D the reader is referred to [33]. Another area which has received considerable attention, is the modelling of thin wires by interfacing local solutions to the mesh without increasing mesh resolution. This
C. Christopoulos: Multi-Scale Modelling in Time-Domain Electromagnetics 103
After some manipulation one obtains E z (r0 ) = L d
Fig. 3. Three-port network connected to the three TLM nodes adjacent to an edge in a two-dimensional problem.
is a major area of concern in differential equation (DE) methods. The situation as regards models of thin wires is less critical in integral equation (IE) methods [34]. In the remainder of this section we thus focus on two major classes of problems, thin wires and thin panels with complex features, both in connection with TD-DE methods.
3.1 Thin wire models
∂I 1 ∂Q + ∂t Cd ∂t
(7)
where, L d , Cd are the wire inductance and capacitance per unit length with respect to some reference return conductor of radius r0 and E z (r0 ) is the value of the axial electric field at r0 . Equation (7) shows how the electric field and the conductor current and charge are related. This equation was obtained on the assumption that the fields near the wire are quasi-static. A second assumption is on the choice of the radius r0 of the reference conductor. This is typically chosen to be half the mesh size ∆/2. However, this is not always a straightforward choice. This basic technique was developed further to deal with arbitrarily oriented wires as described in [36, 37]. The extension to multi-conductor wire looms has been reported in [38]. Parallel developments have also taken place in TLM. The interface between the wire and the fields is configured as an equivalent circuit which is then connected to the TLM nodes. This is shown schematically in Fig. 4. In this figure the equivalent voltage source E z ∆ is obtained from the TLM nodes contributing to the E z field component i.e. V10 (x, y, z) and V6 (x + ∆, y, z). The wire equation is modelled by the circuit in Fig. 5, where the stub and line parameters are chosen to model the correct values of
As a broad generalization, thin wire models fall into three categories: i) ii)
special thin wire interfaces models based on contour−path integration formulations iii) hybrids between IE formulations for the wire and DE formulations for the rest of the problem. Within these categories there are various overlapping strands of work. We place particular emphasis on the thin wire formulations as they appear to offer the greatest flexibility in being able to incorporate complete circuits including lumped components into the mesh. 3.1.1 Thin wire formulations The first significant work in this area was reported in [35]. It was aimed at accounting for the presence of the wire without resolving in the mesh its radial dimensions. According to this formulation the problem is split into wire and field parts which are coupled. In order to achieve this coupling the quasi-static assumption is employed i.e. that if the current and the charge per unit length on the wire are I and Q respectively, then the corresponding azimuthal magnetic and radial electric field components are Hϑ (r) =
I 2πr
(5)
Er (r) =
Q 2πεr
(6)
Fig. 4. A z-directed wire in a three-dimensional mesh of TLM nodes.
Fig. 5. Model of the wire and its coupling to the field for the arrangement of Fig. 4.
104 C. Christopoulos: Multi-Scale Modelling in Time-Domain Electromagnetics the wire capacitance and inductance and to achieve synchronism with the rest of the mesh [39]. This interface is suitable for placing thin wires between nodes and it is therefore operated as part of the ‘connection’ process in TLM. An alternative interface for thin wires threaded through the centre of the node, which is operated as part of the ‘scattering’ was described in [40]. The treatment of multiconductor wire looms in TLM was fully described in [41]. The multiconductor interface follows closely the principles outlined for a single wire and is shown in Fig. 6. In this Figure U and V indicate the vectors of the incident and reflected voltages on the wire loom and z n , u z represent the coupling to the field. The construction of the interface proceeds as follows: • We choose the link line impedance matrix Zw to model all the required wire loom capacitance • We calculate the modelled inductance from the known time step ∆t and thus obtain the deficit in modelled inductance L stub . • Then we introduce a short-circuited stub with impedance Z stub = 2L stub /∆t. Thus the interface in schematic form is as shown in Fig. 7. In this figure the four parallel lines represent the coupling to the field i.e. i V i + V6i + V7i + V10 uz = 5 , 2
z n = Z 0 /4
(8)
where Z 0 is the impedance of free space. The Thevenin equivalent circuit of the interface is shown in Fig. 8, where 1v = [1, 1, . . . 1]T ,
1h = [1, 1, . . . 1]
Us , Ul , Uh represent the incident voltage vectors for the stub, low and high sides of the wire respectively. From Fig. 8 the loop current is obtained as Itot =
2(Uh − U − Us + 1v u z ) 2Z w + z s + 1v Z n 1h
(10)
The calculation then proceeds as in standard TLM with ‘scattering’, Vh = Uh − Z w Itot
(11)
V = U + Z w Itot
(12)
Vs = Us − Z s Itot
(13)
vz = u z − z n 1h Itot
(14)
(9)
Fig. 6. The structure of an n-wire bundle in a three-dimensional mesh.
Fig. 7. The structure of the interface for the bundle in Fig. 6.
Fig. 8. Thevenin equivalent of the circuit in Fig. 7 for the calculation of the loop current.
followed by ‘connection’. Multiwire junctions can also be treated as described in [41]. Wires inclined relative to the coordinate axis may be described by a piecewise approximation with some adjustment to obtain the correct propagation delays. All the approaches based on the model described in [35] suffer from several weaknesses. These are: • A single wire cannot be precisely located inside each numerical cell-it is assummed to be at the centre • In the case of multiwire looms the location of each wire within the loom cannot be taken fully into account and hence important timing information is missing from the coupling to each wire • Restrictions are imposed on the shape and size of each wire • In some cases, empirical factors are used to improve performance. It is believed that these weaknesses are a consequence of the fact that the local field expansions in (5) and (6) represent only the first (symmetrical) mode of the field near the
C. Christopoulos: Multi-Scale Modelling in Time-Domain Electromagnetics 105
wire. If more modes can be included a much better representation of the wire could be made. The field near an infinitely long z-directed wire of radius a is given by the expression [42], ∞ J(k0 a) E z (r, ϕ) = Bn e jnϕ Jn (k0r) − Nn (k0r) Nn (k0 a) n=−∞ (15) Hϕ (r, ϕ) =
1 ∂E z jωµ0 ∂r
(16)
where, J and N are Bessel and Neumann functions. In numerical schemes such as FDTD and TLM the fields are sampled at a finite number of mesh points in space. As an illustration, in a 2D mesh fields are sampled at four positions in each node. We can therefore include the four lowest order modes n = 0, 1, −1, 2. We now have the superposition of four rather than one modes and therefore, in principle, a far better accuracy and positional resolution should be expected [43]. The mapping between the modes and the sampled voltages or fields can then be implemented in TLM or FDTD. We sketch here the construction of the model in TLM [44]. In the case of a 2D calculation the model can be formulated in terms of a simple circuit equivalent which is interfaced to the standard TLM mesh in a manner similar to [39]. This however is neither necessary nor always the most efficient way in more complex 3D configurations. The computation proceeds by first calculating the impedance seen at the edges of the node. We make use of the small argument expansions for Bessel and Neumann functions as this is consistent with the normal modelling requirement of ten segments per wavelength. This approximation simplifies considerably the model without significant loss in accuracy. We find, ∆ E z0 jωµ0 ∆ ln , n=0 (17) Hϕ0 a E zn jωµ0 ∆2|n| − a2|n| ∆ 2|n| , |n| Hϕn ∆ + a2|n|
n = 0
(18)
where Z L is the impedance of the lines in the standard TLM formulation. For the case when n = 0, a cascade of the link line of impedance Z and a shortcircuit stub Z s is used. In this case, we obtain as the input impedance of this combination, √ Z in = jω µ0 ε0 ∆(Z + Z s ) (21) Equating this expression with (17) we obtain the required stub impedance, ∆ −Z (22) Z s = Z L ln a All these expressions can be brought into the 2D model of a wire node shown schematically in Fig. 9. The wire current is obtained directly from the stub current. It can be seen from this figure that modes n = 1 (cos ϕ and sin ϕ) and n = 2 with sample values of (0, 1, 0, −1), (1, 0, −1, 0), (1, −1, 1, −1) respectively at the node ports, do not see the stub (Vs = 0). Only mode n = 0 sees the stub. We have accounted fully in the model for modes n = 0, 1. Mode 2 should have a different impedance and phase constant, but for simplicity we do not account for this in the model as it does not appear to significantly degrade accuracy. It should be emphasized that the same calculation can be done without resort to an equivalent circuit. This simply requires the transform of mesh voltage quantities to modal quantities which after they see the correct modal impedances are transformed back into mesh quantities for propagation to the rest of the mesh. This indeed is the best way of proceeding when more complex configurations are studied. The outline of this calculation is as follows: Let us call [T ] the matrix transforming modal to mesh quantities and A the column vector of the four modal amplitudes. Then, the field quantities at the four ports of a 2D mesh are, E z = [T ]A,
Hϕ = [T ][Y ]A
(23)
where [Y] is a diagonal matrix with elements the modal admittances. From the normal TLM equivalences, Ez = V i + V r
,
Hϕ = Y L (V i − V r )
(24)
where in theses formulae ∆ is half the mesh size. The network is constructed by choosing link lines and stubs of suitable parameters so that the same impedance is seen as in [17, 18]. For the case n = 0, we choose a link line of impedance Z and phase constant k L which is terminated by a short-circuit. Employing (18) we obtain, jZ tan(k L ∆) = using
√ k L = ω µ0 ε0 ,
jωµ0 ∆2|n| − a2|n| |n| ∆2|n| + a2|n|
(19)
tan(k L ∆) k L ∆
the required impedance of the link line is obtained, Z=
Z L ∆2|n| − a2|n| |n| ∆2|n| + a2|n|
(20)
Fig. 9. Wire interface in 2D using the Modal Expansion Technique (MET) (after [43]).
106 C. Christopoulos: Multi-Scale Modelling in Time-Domain Electromagnetics Combining equation (23)(24), (Y L [I] + [Y ]) [T ]T V r = (Y L [I] − [Y ]) [T ]T V i
(25)
This equation relates the reflected to the incident voltages and represents the scattering process in the frequency domain. Matrix, (Y L [I] + [Y ])−1 (Y L [I] − [Y ]) is diagonal with elements the reflection coefficients (Z n − Z L ) /(Z L + Z n ) for each mode. Each of these coefficients represents a shift equal to, −e− j2k L ∆/n Fig. 11. Comparison of the MET with analytical results and other wire models (mesh resolution 50 mm, wire radius 5 mm). An incident field is scattered from the wire and the total field is observed 100mm in front of the wire.
for modes n = 1, 2, and e− j2k L ∆ for mode n = 0. This is shown schematically in Fig. 10 where for simplicity the shift is kept the same for mode n = 2. The accuracy of this model is shown in Fig. 11 [44], where in addition to comparisons with analytical results, the results from two other wire models are shown [39, 40]. The model uses no empirical factors and there are no restrictions on the size of the wire relative to the mesh size. The case of an offset wire, i.e. a wire not centred at the node can also be dealt with using this approach by exploiting the addition Theorems for Bessel Functions [45], e
jma
Z m (k R) =
∞
Jn (kr1 )Z n+m (kr01 )e
jnϕ
(26)
n=−∞
where all the relevant quantities are defined in Fig. 12. The field around the wire is expanded into cylindrical modes as before and then a mapping is implemented into the ports of the node using equation (26) [46]. The excellent behaviour of the model for the case of an offset wire is shown in Fig. 13. Many features of the model have been
Fig. 10. Schematic of the scattering and transformation of modal to total field quantities in a 2D mesh (after [44]).
Fig. 12. Definition of quantities used for the Bessel function addition Theorem.
implemented in 3D [47, 48]. It is worth stressing that results are obtained using closed form local solutions for the field around the wire, unlike the techniques described in [49, 50] where a full MoM numerical solution is necessary to calculate the wire current. The major advantage of the local solutions is simplicity and a minimal computational overhead.
Fig. 13. Modelling of an offset thin wire using the MET: wire radius = 0.1, offset is along the 45 degree line and the field in front of the wire is observed.
C. Christopoulos: Multi-Scale Modelling in Time-Domain Electromagnetics 107
3.1.2 Contour-path models
3.2 Models of thin panels with complex features
Contour path models have been used extensively in FDTD. The basic work is described in [49]. The following steps are implemented: • The calculation of an equivalent radius to replace a wire loom (a bundle) by a single equivalent wire • Faraday’s Law contour integral is applied to modify the FDTD algorithm and thus account for the presence of the wire • The equivalent electric and magnetic current sources are calculated on a virtual surface surrounding the wire • The EFIE and MoM are applied to calculate the wire current. A similar approach has been used to account for other features such as curved surfaces [51]. Further extensions of this technique may be found in [52, 53], where attention is given to end effects.
A problem similar in complexity and importance as that of fine wires is the description in numerical models of thin conducting panels containing a large number of holes, as are typically present in many equipment for cooling etc. Several possible problems can be envisaged e.g. thin conducting panels with/without holes, composite panels consisting of overlapping layers with different properties (carbon fibre composites), sandwiches of conducting and absorbing layers etc. In all these cases the challenging modelling feature is that the panel thickness and hole size are smaller than the nodal size. As explained in connection with wires, it is unprofitable to reduce the coarseness of the mesh to accommodate these features. We thus have another multiscale problem of the type already encountered. Early work in this area in connection with thin CFC panels was described in [57, 58]. In this approach onedimensional ladder networks connected between nodes were used to describe propagation through the panel. However the technique cannot deal with other features such as holes and lacks the generality and systematic parametrisation required in most problems. Taking the example of a perforated panel placed between nodes, it is clear that there will be both reflections and transmission through the panel. The relevant coefficients are in general frequency dependent and therefore scattering takes place described by a scattering matrix, R(ω) T(ω) [S(ω)] = T(ω) R(ω) (27)
3.1.3 Hybrid models As a general illustration of the hybridizing of methods to exploit their best features, we discuss briefly the hybridization of Integral Equation methods with TLM. An example is shown in [53, 54], where two separate regions with complex features are modelled by TLM and the intervening space, which is essentially free space, is modelled very efficiently by an IE method. However, we are interested here in the inverse process, where a wire inside a surface is modelled by an IE technique and the surrounding space by TLM. Effectively, this is another way of accounting for fine wires in TLM by hybridizing with a technique which is well suited to modelling wires. This approach is described in [55, 56]. The calculation proceeds as follows: • A surface SB is defined so that it completely surrounds the wire and is coincident with the boundaries of TLM cells • The total fields on SB are defined as a superposition of the incident fields from all sources located outside SB , and the fields which are radiated by the wire • The fields radiated by the wire, calculated at SB , are used to obtain the equivalent surface currents which are sources of radiation outside SB • The TLM technique, the equivalent current sources and other incident fields from sources outside SB are used to find the total fields outside SB • Inside SB , we enforce the condition that the tangential electric field (incident plus radiated) is zero at the wire surface. This conditions yields the current along the wire segments • From the current distribution along the wire, the radiated fields at the surface SB are obtained, which are then mapped onto the TLM pulses and superimposed on other pulses there to excite the total field outside SB . This process works satisfactorily but considerably more effort is involved compared to interfacing TLM to local solutions.
This scattering at the panel must now be included as part of the connection process in TLM, after it has been converted into a time-domain algorithm. This conversion results essentially to a digital filter algorithm which is embedded into the ‘connection’ process. The procedure is described in [60–62]. The basic steps are as follows: • Each coefficient F in (31) is put into standard Pade form NP
F(s) =
i=0 NP
bi s i (28) ai s i
i=0
where NP is the number of poles, decided from accuracy and computational cost considerations. • The frequency-domain Prony method [63] is employed to estimate ai and bi above from a set of NS measurements. After these coefficients have been estimated equation (28) is put into the form F(s) =
b NP (s − sz0 )(s − sz1 ) . . . (s − sz(NP−1) ) (s − s p0 )(s − s p1 ) . . . (s − s p(NP−1) )
where s pi , szi are the poles and zeroes.
(29)
108 C. Christopoulos: Multi-Scale Modelling in Time-Domain Electromagnetics • Equation (29) is converted into the discrete timedomain by the application of the bilinear transformation and after some manipulation we obtain NP
F(z) = B0 +
Bi z −i
i=1 NP
1+
(30) Ai
z −i
i=1
• As an illustration we assume that equation (30) represents a reflection coefficient, then F(z) = R(z) = k+1 V i /k V r
(31)
After some manipulation equation (31) can be put in a state space form, where the output equation is V i = B0 V r + B X
(32)
4. Discussion
and the state equation X = z −1 [A ]X + z −1 1 V
where B0 , B , [A ] are known and 1 = [1, 0 . . . 0]T . A signal flow diagram representing equations (32) and (33) is shown in Fig. 14. Typical results for the Shielding Effectiveness (SE) of a cabinet with a perforated wall are shown in Fig. 15. This simulation shows the accuracy to be expected and it is done at a fraction of the computational costs of a full-field model, where details of each perforation would have to be specifically included in the mesh. The process described is very general, as the characterization of the panel in the frequency domain (equation 27) can be done experimentally, analytically, or by a numerical model of a single perforation. Work is in progress to generalise this models to allow for arbitrary placement and to deal with many different generic objects that may be inserted into the mesh.
r
(33)
Fig. 14. Signal flow diagram of the ‘connection’ process in the presence of fine features using the Digital Filtering Interface (DFI).
Fig. 15. Shielding Effectiveness of an enclosure using the DFI and comparison with experiments: cabinet size (292 × 292 × 118mm3 ), area of perforations (212 × 78 mm3 ), hole radius 4 mm arranged in an equilateral triangle lattice of sides 10 mm. Simulations are shown for two different mesh resolutions (after [61]).
Several techniques developed to tackle multiscale problems were described. The number of publications, the variety of approaches and the period of time covered by the review testify to the importance and difficulty of this subject. It is not an understatement to say that without continuing progress in this area, it will be impossible to envisage electromagnetic computation in large practical systems of the type encountered in electromagnetic compatibility and signal integrity applications. Due to space limitations not all techniques could be covered in the same depth and detail. However, some conclusions are possible and should help point the way to future work. Multigridding schemes have a role to play but on their own cannot address fully the main modelling challenges. Some mesh refinement is desirable in certain areas and the modeller should have the facility of using one of the schemes described including the hybrid nodes. Even assuming that an entirely satisfactory multigrid scheme is developed, the effort of defining in detail the geometry of the problem in a very fine mesh appears to be out of proportion to any benefits, considering the limited resolution of the surrounding coarse mesh. It would be interesting to see developments which allow more flexibility in defining the structure of the mesh near particular features (as is possible in FEM) without sacrificing ease and simplicity of computation. Thin-wire nodes are highly desirable and effective in reducing complexity and are set to be an interesting and rewarding area of research. The Modal Expansion Technique (MET) appears to offer greater flexibility and accuracy compared to other approaches and could be generalised to deal with other types of problem, where modal expansion is the link between particular fine features and the surrounding coarse mesh. One can envisage MET acting as the ‘glue’ between the mesh and particular complex attachments to it. Without doubt, in addition to significant gains in run-time efficiency, major benefits are also possible in ease of problem definition. The treatment of essentially planar structures with complex details (e.g. perforated/composite panels) can be
C. Christopoulos: Multi-Scale Modelling in Time-Domain Electromagnetics 109
dealt with in a systematic and efficient way by the Digital Filtering Interface (DFI). The ease of characterization, the scope for automation and the possibility of developing a library of often used components and sub-assemblies, described in large-scale EM models by a small number of pre-calculated parameters as a DFI, is an exciting prospect. One can envisage work to expand this technique to encompass a wider range of components and structures. In conclusion, a long-term effort is required to develop further these techniques and to deal effectively with complexity reduction. In this context the pioneering work of Professor Russer should be mentioned [64, 65]. Acknowledgement. I wish to thank my colleagues Drs. P. Sewell, D.W.P. Thomas, J. Paul, V. Podlozny and Y.K. Choong for supplying me with material for this article and for their major contributions to my understanding of this topic. I also thank the EPSRC (UK) for support received over a number of years for multi-scale research. Thanks also to the Agency for Science Technology and Research (Singapore) and the Nanyang Technological University for their support and hospitality during the writing of this article while on attachment to NTU.
References [1] Taflove, A.: Computational electrodynamics: The Finite Difference Time-Domain Method. Norwood MA: Artech House, 1995. [2] Volakis, J.L.; Chaterjee, A.; Kempel, L.C.: Finite Element Method for Electromagnetics. New York: IEEE Press, 1998. [3] Harrington, R.F.: Field Computation by Moment Methods. New York: McMillan, 1968. [4] Christopoulos, C.: The Transmission-Line Modeling MethodTLM. New York: IEEE Press, 1995. [5] Kunz, K.S.; Simpson, L.: A technique for increasing resolution of finite-difference solution to maxwell’s equations. IEEE Trans. on EMC 23 (1981), 419–422. [6] Kim, I.S.; Hoefer, W.J.R.: A local mesh refinement algorithm for the FDTD method using maxwell’s curl equations. IEEE Trans. on MTT 38 (1990), 812–815. [7] Zivanovic, S.S.; Yee, K.S.; Mei, K.K.: A subgridding method for the time-domain finite-difference method to solve maxwell’s equations. IEEE Trans. on MTT 39 (1991), 471–479. [8] Prescott, D.T.; Shuley, N.V.: A method for incorporating different sized cells into the finite-difference time-domain analysis technique. IEEE Microwave and Guided Wave Letters 2 (1992), 434–436. [9] Okoniewski, M.; Okoniewska, E.; Stuchly, M.A.: Threedimensional subgridding algorithm for FDTD. IEEE Trans, on Antennas and Propagation 45 (1997), 422–429. [10] Chevalier, M.W.; Luebbers, R.J.; Cable, V.P.: FDTD local grid with material traverse. IEEE Trans. on Antennas and Propagation 45 (1997), 411–421. [11] Krishnaiah, K.M.; Railton, C.J.: Passive equivalent circuit of FDTD:AN application to subgridding. Electronics Letters 33 (1997), 1277–1278. [12] Krishnaiah, K.M.; Railton, C.J.: A stable subgridding algorithm and its application to eigenvalue problems. IEEE Trans. on Microwave Theory and Techniques 47 (1999), 620–628.
[13] Thoma, P.; Weiland, T.: A consistent subgridding scheme for the FDTD method. Int. J. Numer. Modeling 9 (1996), 359– 374. [14] White, M.J.; Yun, Z.; Iskander, M.F.: A new 3-D FDTD multigrid technique with dielectric traverse capabilities. IEEE Trans. on Microwave Theory and Techniques 49 (2001), 422– 430. [15] Herring, J.L.; Christopoulos, C.: Multigrid TLM method for solving EM field problems. Electronics Letters 27 (1991), 1794–1795. [16] Herring, J.L.; Christopoulos, C.: Solving EM field problems using a multiple grid TLM method. IEEE Trans. on Antennas and Propagation 42 (1994), 1654–1658. [17] Herring, J.L.; Christopoulos, C.: The application of different meshing techniques to EMC problems. – In: Proc. 9th Annual Rev. of Progress in Applied Computational Electromagnetics. Monterey, 1993, 755–762. [18] Wlodarczyk, J.: New multigrid interface for the TLM method. Electronics Letters 32 (1996), 1111–1112. [19] Trenkic, V.; Christopoulos, C.; Benson, T.M.: On the time step in hybrid symmetrical condensed TLM nodes. IEEE Trans. on Microwave Theory and Techniques 43 (1995), 2172–2174. [20] German, F.J.; Svigeli, J.A.; Mitra, R.: A numerical comparison of dispersion in irregular graded TLM and FDTD meshes. – In: 12th Annual Review of Progress in Applied Computational Electromagnetics. Monterey, 1996, 270–278. [21] van Bladel, J.: Singular electromagnetic fields and sources. New York: Oxford University Press, 1991. [22] Mexner, J.: The behaviour of electromagnetic fields at edges. IEEE Trans. on Antennas and Propagation 20 (1972), 302– 307. [23] Duffy, A.P.; Benson, T.M.; Christopoulos, C.: New method for accurate modelling of wires using TLM. Electronics Letters 29 (1993), 224–226. [24] Herring, J.L.; Hoefer, W.J.R.: Improved correction for 3-D TLM coarseness error. Electronics Letters 30 (1994), 1149– 1150. [25] Cascio, L.; Tardioli, G.; Rozzi, T.; Hoefer, W.J.R.: A quasistatic modification of TLM at the knife edge and 900 wedge singularities. IEEE Trans. on Microwave Theory and Techniques 44 (1996), 2519–2524. [26] Fujii, M.; Hoefer, W.J.R.: Field singularity correction in 2-D time-domain haar-wavelet modeling of waveguide components. IEEE Trans. on Microwave Theory and Techniques 49 (2001), 685–691. [27] Shorthouse, D.B.; Railton, C.: The incorporation of static field solutions into FDTD algorithms. IEEE Trans. on Microwave Theory and Techniques 40 (1992), 986–994. [28] Beilenhoff, K.; Heinrich, W.: Treatment of field singularities in the finite-difference approximation. Proc. of IEEE Symp. On Microwave Theory and Techniques, Atlanta, 1993. 979– 982. [29] Lindenmeier, S.; Russer, P.; Heinrich, W.: Hybrid dynamicstatic finite-difference approach for MMIC design. Proc. IEEE Symp. On Microwave Theory and Techniques, San Francisco, 1996. 197–200. [30] Lindenmeier, S.; Russer, P.; Heinrich, W.: Application of the hybrid dynamic-static finite-difference approach on 3DMMIC structures. Proc. of Progress on Applied Computational Electromagnetics, Monterey, 1997. 1182–1188. [31] Craddock, I.J.; Railton, C.J.: Stable inclusion of a priory knowledge of field behaviour in the FDTD algorithm:application to the analysis of microstrips lines. Proc. IEEE Symp. On Antennas and Propagation, Baltimore, 1996. 1300–1303.
110 C. Christopoulos: Multi-Scale Modelling in Time-Domain Electromagnetics [32] Esselle, K.P.; Okoniewski, M.; Stuchly, M.A.: Analysis of sharp metal edges at 45◦ to the FDTD grid. IEEE Microwave and Guided Wave Letters 9 (1999), 221–223. [33] Farina, M.; Rozzi, T.: Numerical investigation of the field and current behavior near lossy edges. IEEE Trans. on Microwave Theory and Techniques 49 (2001), 1355–1358. [34] Cui, T.J.; Chew, W.C.: Accurate analysis of wire structures from very-low frequency to microwave frequency. IEEE Trans. on Antennas and Propagation 50 (2002), 301–307. [35] Holland, R.; Simpson, J.W.: Finite-difference analysis of EMP coupling to thin struts and wires. IEEE Trans. on Electromagnetic Compatibility 23 (1981), 88–97. [36] Ledfelt, G.: A stable subcell model for arbitrarily oriented thin wires for the FDTD method. Int. J. Numer. Modelling 15 (to appear 2002). [37] Edelvik, F.: A new technique for accurate and stable modeling of arbitrarily oriented thin wires in the FDTD method. Uppsala University Technical Report 2002-016, www.it.uu.se/research/reports/2002-016/. [38] Berenger, J.-P.: A multiwire formulation for the FDTD method. IEEE Trans. on Electromagnetic Compatibility 42 (2000), 257–264. [39] Wlodarczyk, J.; Johns, D.P.: New wire interface for graded 3D TLM. Electronics Letters 28 (1992), 728–729. [40] Porti, J.A.; Morente, J.A.; Khalladi, M.; Gallego, A.: Comparison of thin wire models for the TLM method. Electronics Letters 28 (1992), 1910–1911. [41] Wlodarczyk, J.; Trenkic, V.; Scaramuzza, R.A.; Christopoulos, C.: A fully integrated multiconductor model for TLM. IEEE Trans. on Microwave Theory and Techniques 46 (1998), 2431–2437. [42] Dudley, D.G.: Mathematical Foundations of Electromagnetic Theory. New York: IEEE Press, 1994. [43] Choong, Y.K.; Sewell, P.; Christopoulos, C.: Accurate wire representation in numerical models for high-frequency simulation. Electronics Letters 37 (2001), 280–282. [44] Choong, Y.K.; Sewell, P.; Christopoulos, C.: New thin wire formulation for time-domain differential equation models. Int. J. Numer. Modelling 15 (2002 to appear). [45] Watson, G.N.: A Treatise on the Theory of Bessel Functions. Cambridge: CUP, 1958. [46] Choong, Y.K.; Sewell, P.; Christopoulos, C.: Accurate modelling of arbitrarily placed thin wire in a coarse mesh. Proc IEE Sci. Meas. Technnol. (to appear). [47] Choong, Y. K.; Sewell, P.; Christopoulos, C.: Hybrid node decsrption of thin wires based on analytical field representation in EMC problems. Proc. IEEE EMC Symp., Minneapolis, 2002. [48] Sewell, P.; Choong, Y.K.; Christopoulos, C.: An accurate thin wire model for 3D TLM simulations. IEEE Trans. on Electromagnetic Compatibility (accepted). [49] Umanshankar, K.R.; Taflove, A.; Bekker, B.: Calculation and experimental validation of induced currents on coupled wires in an arbitrarily shaped cavity. IEEE Trans. on Antennas and Propagation 35 (1987), 1248–1257. [50] Lindenmeier, S.; Christopoulos, C.; Russer, P.: Thin wire modelling with the TLM IE method. Proc. of Progress in Applied Computational Electromagnetics, Monterey, 2000. 587–593. [51] Jurgens, T.G.; Taflove, A.; Umanshankar, K.: FDTD modeling of curved surfaces validation of induced currents on coupled wires in an arbitrarily shaped cavity. IEEE Trans. on Antennas and Propagation 40 (1992), 357–366. [52] Douglas, M.; Okoniewski, M.; Stuchly, M.A.: Accurate modelling of thin wire antennas in the FDTD method. Microwave and Optical Technology Letters 4 (1999), 261–265.
[53] Makinen, R.M.; Juntunen, J.S.; Kivikovski, M.A.: An improved thin-wire model for the FDTD method. IEEE Trans. on Microwave Theory and Techniques 50 (2002), 1245–1255. [54] Lindenmeier, S.; Pierantoni, L.; Russer, P.: Hybrid space discretising-integral equation method for the numerical modelling of transient interference. IEEE Trans. on Electromagnetic Compatibility 41 (1999), 425–430. [55] Pierantoni, L.; Lindenmeier, S.; Russer, P.: Efficient analysis and modelling of the radiation of microstrip lines and patch antennas by the TLM-integral equation (TLM-IE) method. Int. J. Numer. Modelling 12 (1999), 329–340. [56] Lindenmeier, S.; Christopoulos, C.; Russer, P.: Thin wire modelling with the TLM-IE method. Proc. Progress in Applied Computational Electromagnetics, Monterey, 2000. 587– 593. [57] Lindenmeier, S.; Christopoulos, C.; Russer, P.: Methods for the modelling of thin wire structures with the TLM method. Proc. IEEE Symp. On Microwave Theory and Techniques, Boston, 2000. 387–390. [58] Trenkic, V.; Duffy, A.P.; Benson, T.M.; Christopoulos, C.: Numerical simulation of penetration and coupling using the TLM method. Proc. EMC Symp., Rome, 1994. 321–326. [59] Trenkic, V.; Christopoulos, C.; Benson, T.M.: Numerical simulation of polymers and other materials for electronic shielding applications. Proc. POLYMAT, London, 1994. 384– 387. [60] Paul, J.; Podlozny, V.; Thomas, D.W.P.; Christopoulos, C.: Time domain simulation of thin materials boundaries and thin panels using digital filters in TLM. Turkish J. of Electrical Eng. And Computer Sci. 10 (2002), 185–198. [61] Podlozny, V.; Paul, J.; Christopoulos, C.: Efficient calculation of shielding effectiveness of equipment cabinets in full-field numerical models. Proc. EMC Europe, Sorrento, 2000. [62] Paul, J.; Podlozny, V.; Christopoulos, C.: The use of digital filter techniques for the simulation of fine features in EMC problems solved in the time-domain. IEEE Trans. on Electromagnetic Compatibility (accepted). [63] Brittingham, J.N.; Miller, E.K.; Williams, J.L.: Pole extraction from real frequency information. Proc. IEEE 68 (1980), 263–273. [64] Russer, P.; Cangellaris, A.: Network oriented modeling, complexity reduction and system identification techniques for electromagnetic systems. Proc. 4th Int. Workshop on Computational electromagnetics in the Time-Domain: TLM/FDTD and Related Techniques, Nottingham, 2001. 105–122. [65] Russer, P.: Network oriented modeling of radiating electromagnetic structures. Turkish J. Elec. Engin. 10 (2002), 147– 162.
Christos Christopoulos was born in Patras, Greece in 1946. He received his first degree from the National Technical University of Athens and his MSc and DPhil from Sussex University. He is currently Professor of Electrical Engineering at the University of Nottingham. His interests include, computational electromagnetics, electromagnetic compatibility and the protection of large networks using fast transients. He is the author of more than 250 scientific articles and five books. He is current Chairman of the IEE Professional Network on EMC and an Associate Editor of the IEEE EMC Transactions.
Multi-Scale Modelling in Time-Domain Electromagnetics Christos Christopoulos Dedicated to Professor Peter Russer on the occasion of his 60th birthday Abstract: The paper reviews the techniques used to describe efficiently in time-domain electromagnetics the presence of fine physical features in an otherwise coarse numerical mesh. This is a major area of concern in electromagnetic modelling and simulation as it affects crucially modelling capabilities in dealing with complex practical configurations. Two broad classes of modelling methods are identified, those where the basic approach is to distort and/or refine the numerical mesh, or, interface the mesh to local solutions. A typical example of a technique in the former category is multi-gridding and in the latter mesh-to-wire interfaces. The emphasis of the paper is to describe the essence of the different approaches, their strengths and weaknesses and to identify areas of further work. Keywords: Computational electromagnetics, Multi-scale modelling, Time-domain simulation
1. Introduction The modelling of large-scale electromagnetic problems involves several stages. It assumes knowledge of the basic system configuration and components to proceed in establishing a numerical model which can be computed to offer the required outputs. It is relatively straightforward to devise efficient models in cases where all physical dimensions are comparable in electrical size (size relative to the wavelength). In such cases a numerical mesh with the appropriate spatial resolution can be devised to model the entire problem space. Several generic numerical techniques are available for general electromagnetic modelling including, Finite Difference TimeDomain (FDTD) [1], Finite Element Method (FEM) [2], Method of Moments (MoM) [3] and Transmission-Line Modelling (TLM) [4]. However, practical problems often do not conform to this requirement. This is particularly true in Electromagnetic Compatibility (EMC) problems where it is very common to have, in the same problem, very different physical scales (relative to the wavelength). A typical example is a wire inside an equipment cabinet. The wire diameter is typically much smaller than the wavelength while the cabinet size may be comparable to the wavelength. This and similar problems are described
Received September 1, 2002. Revised January 2, 2003. C. Christopoulos, Electromagnetics Research Group, University of Nottingham, Nottingham, NG7 2RD, U.K., Tel. +44 115 951 5557, Fax: +44 115 951 5616. E-mail: [email protected] ¨ 57 (2003) No. 2, 100−110 Int. J. Electron. Commun. (AEU)
as multi-scale problems and present severe difficulties to the modeller. It has to be stated that no amount of increase in the size and power of computers is likely to overcome this problem as for many systems full explicit knowledge of all details is impractical or impossible even if the computational power were available. An additional consideration, which is often overlooked, is the advantage of keeping models as simple as possible, thus offering intellectual economy and scope for creativity. Therefore for practical and methodological reasons it is essential that, in a large problem, features of a very small electrical scale or complex details are dealt with separately from the rest of the problem which is on a different scale. This requirement applies more or less irrespective of the actual numerical method used for field simulation be it FEM, FDTD, TLM etc. In tackling this problem two broad approaches have emerged. First, the mesh is refined or distorted locally to increase spatial resolution in areas where fine features are present and therefore fields may vary rapidly in space. In this way, increased resolution is applied only where it is needed thus maintaining computational efficiency. Difficulties however remain at the interface between meshes of different resolution and/or shape as will be discussed in the following sections. Second, the mesh is left unchanged at a more or less uniform resolution set with global considerations in mind. A local solution is obtained for the fields near the fine feature and this then is interfaced to the mesh. The difficulties here are both in obtaining local solutions and in devising appropriate interfaces. Both these approaches will be described in some detail in the following sections. As regards the types of fine feature that may be encountered, one should include, fine wires, wire looms, thin panels with/without perforations, fine conducting tracks, entire assemblies of such components etc. The objective is to present what has been done so far and to stimulate further research work in this area.
2. Mesh refinement and distortion A regular Cartesian mesh is shown in Fig. 1. Let us consider that there is a fine wire of radius a as shown in crosssection. In a mesh with spatial resolution ∆ and when the radius of the wire is smaller than the mesh size, the impact of the wire as regards its diameter and actual position cannot be accurately represented. One option is to reduce ∆ until it matches the radius of the wire. The however entails a heavy computational cost as an unnecessarily fine 1434-8411/03/57/2-100 $15.00/0
C. Christopoulos: Multi-Scale Modelling in Time-Domain Electromagnetics 101
Fig. 1. Schematic of various mesh arrangements: (a) regular mesh, (b) multigrid mesh and (c) hybrid mesh.
mesh is used in areas where it is not needed. Halving ∆ in a three-dimensional (3D) mesh will result to an increase in storage by a factor of eight. In practical problems a factor of two is rarely sufficient and very often an increase in resolution by an order of magnitude is necessary. Similar unacceptable increases are also inevitable in run-time making it impossible to tackle in this way any but the simplest problems. One possible way out of the dilemma is to refine the mesh as shown in Fig. 1b. This is known as a multigrid, multiple grid mesh or subgrid. The other option is shown in Fig. 1c and is described as a hybrid mesh. Both examples are for structured meshes of the type used in the main time-domain simulation techniques such as FDTD and TLM. When unstructured meshes are used (e.g. in FEM) there is greater flexibility in configuring the mesh to accommodate fine features but such techniques have not found favour in the time-domain. We therefore limit the discussion here to structured meshes. Each approach is examined separately. The obvious advantages of the scheme in Fig. 1b are that mesh refinement is implemented only where is needed, and that no other model developments are necessary to describe the fine features. Its perfectly conducting surface is registered in the mesh as zero tangential electric field (eg in FDTD) or as a boundary with a reflection coefficient of −1 (eg in TLM). There are however disadvantages and difficulties. Amongst the disadvantages one could mention the fact that the simplicity of the single regular mesh is lost and that each fine feature must be laboriously mapped onto the mesh. The difficulties are more profound. At the interface between the fine and coarse mesh regions information is available at different spatial and temporal rates. This therefore implies some form of averaging/filtering at the interface. What are the physical principles that must be enforced in order to construct such an interface? To-date several solutions to this problem have been proposed. To illustrate the complexity of the interface, let us consider the example shown in Fig. 1b. For this example, the spatial reduction ratio in 3D is 16:1 and over the time step of the coarse mesh there are 4 × 16 pulses from the fine mesh. Even higher reduction ratios may be needed in practice illustrating the complexity of the required interface. One of the earliest attempts to introduce local grids was described in [5] where, however, two separate FDTD simulations were required. The simultaneous coupling at the coarse-fine mesh interface, thus requiring only one FDTD solution, was described in [6]. Here fields are
passed onto the local grid after time and space interpolation. The authors reported two difficulties. First, the abrupt change in mesh size causes spurious reflections. Second, since a wave experiences a different phase error in each mesh, non-physical interactions may take place and therefore they advise that the coarse mesh size is still kept much smaller than the shortest wavelength of interest. A subgridding scheme was also described in [7], referred to as the Variable Step-Size Method (VSSM). A more efficient version known as the Mesh Refinement Algorithm (MRA) was described in [8]. It is shown there that the reflection coefficient at the coarse-fine mesh interface is less than 1%. However, in both these papers the interface must be placed in a homogeneous region. An alternative scheme is described in [9]. The authors have tried several algorithms and interpolation techniques and found that a pulsing overlapping scheme of computations yields the best results. Their interface is usually in homogeneous regions where a reflection coefficient below 60 dB can be expected. They state that if the interface is placed in inhomogeneous regions, somewhat higher reflections can be expected. Placing the interface in inhomogeneous regions was described in [10]. Here the refinement ratio must be an odd number and only the magnetic field is interpolated. Reflections at the interface of less than 1% are reported and no instabilities, although no formal demonstration of stability is offered. Results are shown for a refinement ratio of 3:1, although it is stated that higher ratios (5:1, 7:1) were tried successfully. It was reported in [11] that both the VSSM and MRA techniques suffer from long time instabilities. The same authors have developed an equivalent circuit representation to ensure a stable algorithm [12]. A similar algorithm is also described in [13]. A multigrid current method (MGCM) technique is described in [14], suitable for placing the interface in an inhomogeneous region using any value of the refinement factor. It is stated that the reflection coefficient at the interface is less than 2% and that there is improved stability. The wealth of papers on this topic over a number of years illustrates the inherent difficulties. Parallel developments to those described above in connection with FDTD have also taken place in TLM. Multigrid schemes were proposed in [4, 15, 16]. Since TLM is an equivalent circuit field representation, provided that the conversion across the coarse-fine interface is consistent with circuit theory, stability is guaranteed. The conditions that conversion must satisfy are: charge and energy conservation, no reflections and no delays. It was found that strict adherence to all four conditions was not possible and therefore a compromise is necessary. Stable schemes can be devised but at a slight energy loss at high frequencies. A TLM multigrid scheme that has been found acceptable is based on the number of transmission lines and time steps in each region of the mesh (fine and coarse). The notation for a region ‘i’ is as follows: • si , pi number of lines across the two polarization directions • ti number of time steps before conversion
102 C. Christopoulos: Multi-Scale Modelling in Time-Domain Electromagnetics • n i = si pi ti total number of pulses before conversion Then, in converting pulses from region 2 into region 1, the new pulses to be inserted into region 1 are given by
V1 j =
n2 1 V2k s1 p2 t2
, j = 1, 2 . . . n 1
(1)
k=1
where, V2k are the pulses in region 2 travelling towards the interface. All inserted pulses are taken to be the same to avoid introducing high frequencies. Similarly, for pulses inserted into region 2,
V2 j =
n1 1 V1k s2 p1 t1 k=1
, j = 1, 2 . . . n 2
(2)
where, V1k are pulses in region 1 travelling towards the interface. For the example shown in Fig. 2, s1 = 1, p1 = 1, t1 = 1, n 1 = 1, for the coarse mesh 1 s2 = 2, p2 = 2, t2 = 2, n 2 = 8, for the fine mesh 2 Hence, from equations (1) (2),
V1 =
1 V2 j = V1 2
1 V2k 4 k=1 8
, j = 1, 2 . . . 8
(3)
(4)
An example of a multigrid application is shown in [17]. An improved multigrid interface was described in [18]. Undoubtedly, multigrid schemes offer significant savings (several orders of magnitude) in computation, and in many cases can make the difference between being able to simulate a problem and abandoning the effort on the grounds of insufficient computational power. However, stability is always an active issue and care must be taken in the way these methods are applied. Irrespective of these benefits, it must be questioned whether it is the correct approach to refine the mesh indefinitely, to cope with finer details, rather than seek different approaches based
Fig. 2. Schematic showing conditions at the interface between coarse and fine mesh regions.
on local solutions and/or signal processing techniques, to couple such features to the surrounding coarse mesh. This issue is tackled in the next section. An alternative approach to the multigrid is the hybrid mesh depicted in Fig(1c). Here the nodes do not represent a cube of space, but can be shaped according to the required resolution in each direction. A full theory for this mesh is described in [4] where references to original papers can be found. Various schemes are available with different time step requirements as discussed in [19]. It is difficult to assess the loss in accuracy for high grading ratios. Comparative studies showing mesh behaviour for abrupt changes in grading for TLM and FDTD show that the former exhibits more robust behaviour [20]. The TLM hybrid mesh is unconditionally stable and is used extensively in commercial solvers. As an introduction to the next section the approach of making a local mesh modification is mentioned here. The problem considered is the accurate representation of edges and corners. Near such features fields exhibit singularities and variation which is not linear [21, 22]. However, the sampling of fields in FDTD and TLM is implicitly based on a linear field variation. This inconsistency results in extra energy storage leading to lower frequency estimates at resonance. This problem manifests itself in TLM by longer propagation delays near edges. Various approaches have been proposed to deal with this so called coarseness error. Naturally, the error may be reduced by a multigrid mesh with fine resolution near the edge, or by a hybrid mesh. However, modellers are bound to seek simpler, more efficient, solutions. One approach, described in [23], is to speed up propagation on the nodes immediately adjacent to the edge. A more efficient way of achieving the same objective is described in [24].
3. Mesh interface to local solutions We start in this section by revisiting the edge singularity problem mentioned at the end of the previous section as it illustrates the interfacing of local solutions to the mesh. A more systematic and rigorous approach is described in [25]. Here the field near the edge is obtained from quasi-static approximations and is used to derive an equivalent circuit representing the fields at the edge. This circuit, which in general 2D problems is a threeport, is then connected to the three TLM nodes which interact directly with the edge as shown in Fig. 3. The lumped components in this circuit are obtained from the Z-matrix representation which in turn is obtained from the quasi-static fields. Field singularities near the edge have also been dealt with using wavelets as described in [26]. Work along similar lines has also been done in connection with FDTD as discussed in [27–32]. Published work is limited in two-dimensions. For a discussion of some of the issues involved in 3D the reader is referred to [33]. Another area which has received considerable attention, is the modelling of thin wires by interfacing local solutions to the mesh without increasing mesh resolution. This
C. Christopoulos: Multi-Scale Modelling in Time-Domain Electromagnetics 103
After some manipulation one obtains E z (r0 ) = L d
Fig. 3. Three-port network connected to the three TLM nodes adjacent to an edge in a two-dimensional problem.
is a major area of concern in differential equation (DE) methods. The situation as regards models of thin wires is less critical in integral equation (IE) methods [34]. In the remainder of this section we thus focus on two major classes of problems, thin wires and thin panels with complex features, both in connection with TD-DE methods.
3.1 Thin wire models
∂I 1 ∂Q + ∂t Cd ∂t
(7)
where, L d , Cd are the wire inductance and capacitance per unit length with respect to some reference return conductor of radius r0 and E z (r0 ) is the value of the axial electric field at r0 . Equation (7) shows how the electric field and the conductor current and charge are related. This equation was obtained on the assumption that the fields near the wire are quasi-static. A second assumption is on the choice of the radius r0 of the reference conductor. This is typically chosen to be half the mesh size ∆/2. However, this is not always a straightforward choice. This basic technique was developed further to deal with arbitrarily oriented wires as described in [36, 37]. The extension to multi-conductor wire looms has been reported in [38]. Parallel developments have also taken place in TLM. The interface between the wire and the fields is configured as an equivalent circuit which is then connected to the TLM nodes. This is shown schematically in Fig. 4. In this figure the equivalent voltage source E z ∆ is obtained from the TLM nodes contributing to the E z field component i.e. V10 (x, y, z) and V6 (x + ∆, y, z). The wire equation is modelled by the circuit in Fig. 5, where the stub and line parameters are chosen to model the correct values of
As a broad generalization, thin wire models fall into three categories: i) ii)
special thin wire interfaces models based on contour−path integration formulations iii) hybrids between IE formulations for the wire and DE formulations for the rest of the problem. Within these categories there are various overlapping strands of work. We place particular emphasis on the thin wire formulations as they appear to offer the greatest flexibility in being able to incorporate complete circuits including lumped components into the mesh. 3.1.1 Thin wire formulations The first significant work in this area was reported in [35]. It was aimed at accounting for the presence of the wire without resolving in the mesh its radial dimensions. According to this formulation the problem is split into wire and field parts which are coupled. In order to achieve this coupling the quasi-static assumption is employed i.e. that if the current and the charge per unit length on the wire are I and Q respectively, then the corresponding azimuthal magnetic and radial electric field components are Hϑ (r) =
I 2πr
(5)
Er (r) =
Q 2πεr
(6)
Fig. 4. A z-directed wire in a three-dimensional mesh of TLM nodes.
Fig. 5. Model of the wire and its coupling to the field for the arrangement of Fig. 4.
104 C. Christopoulos: Multi-Scale Modelling in Time-Domain Electromagnetics the wire capacitance and inductance and to achieve synchronism with the rest of the mesh [39]. This interface is suitable for placing thin wires between nodes and it is therefore operated as part of the ‘connection’ process in TLM. An alternative interface for thin wires threaded through the centre of the node, which is operated as part of the ‘scattering’ was described in [40]. The treatment of multiconductor wire looms in TLM was fully described in [41]. The multiconductor interface follows closely the principles outlined for a single wire and is shown in Fig. 6. In this Figure U and V indicate the vectors of the incident and reflected voltages on the wire loom and z n , u z represent the coupling to the field. The construction of the interface proceeds as follows: • We choose the link line impedance matrix Zw to model all the required wire loom capacitance • We calculate the modelled inductance from the known time step ∆t and thus obtain the deficit in modelled inductance L stub . • Then we introduce a short-circuited stub with impedance Z stub = 2L stub /∆t. Thus the interface in schematic form is as shown in Fig. 7. In this figure the four parallel lines represent the coupling to the field i.e. i V i + V6i + V7i + V10 uz = 5 , 2
z n = Z 0 /4
(8)
where Z 0 is the impedance of free space. The Thevenin equivalent circuit of the interface is shown in Fig. 8, where 1v = [1, 1, . . . 1]T ,
1h = [1, 1, . . . 1]
Us , Ul , Uh represent the incident voltage vectors for the stub, low and high sides of the wire respectively. From Fig. 8 the loop current is obtained as Itot =
2(Uh − U − Us + 1v u z ) 2Z w + z s + 1v Z n 1h
(10)
The calculation then proceeds as in standard TLM with ‘scattering’, Vh = Uh − Z w Itot
(11)
V = U + Z w Itot
(12)
Vs = Us − Z s Itot
(13)
vz = u z − z n 1h Itot
(14)
(9)
Fig. 6. The structure of an n-wire bundle in a three-dimensional mesh.
Fig. 7. The structure of the interface for the bundle in Fig. 6.
Fig. 8. Thevenin equivalent of the circuit in Fig. 7 for the calculation of the loop current.
followed by ‘connection’. Multiwire junctions can also be treated as described in [41]. Wires inclined relative to the coordinate axis may be described by a piecewise approximation with some adjustment to obtain the correct propagation delays. All the approaches based on the model described in [35] suffer from several weaknesses. These are: • A single wire cannot be precisely located inside each numerical cell-it is assummed to be at the centre • In the case of multiwire looms the location of each wire within the loom cannot be taken fully into account and hence important timing information is missing from the coupling to each wire • Restrictions are imposed on the shape and size of each wire • In some cases, empirical factors are used to improve performance. It is believed that these weaknesses are a consequence of the fact that the local field expansions in (5) and (6) represent only the first (symmetrical) mode of the field near the
C. Christopoulos: Multi-Scale Modelling in Time-Domain Electromagnetics 105
wire. If more modes can be included a much better representation of the wire could be made. The field near an infinitely long z-directed wire of radius a is given by the expression [42], ∞ J(k0 a) E z (r, ϕ) = Bn e jnϕ Jn (k0r) − Nn (k0r) Nn (k0 a) n=−∞ (15) Hϕ (r, ϕ) =
1 ∂E z jωµ0 ∂r
(16)
where, J and N are Bessel and Neumann functions. In numerical schemes such as FDTD and TLM the fields are sampled at a finite number of mesh points in space. As an illustration, in a 2D mesh fields are sampled at four positions in each node. We can therefore include the four lowest order modes n = 0, 1, −1, 2. We now have the superposition of four rather than one modes and therefore, in principle, a far better accuracy and positional resolution should be expected [43]. The mapping between the modes and the sampled voltages or fields can then be implemented in TLM or FDTD. We sketch here the construction of the model in TLM [44]. In the case of a 2D calculation the model can be formulated in terms of a simple circuit equivalent which is interfaced to the standard TLM mesh in a manner similar to [39]. This however is neither necessary nor always the most efficient way in more complex 3D configurations. The computation proceeds by first calculating the impedance seen at the edges of the node. We make use of the small argument expansions for Bessel and Neumann functions as this is consistent with the normal modelling requirement of ten segments per wavelength. This approximation simplifies considerably the model without significant loss in accuracy. We find, ∆ E z0 jωµ0 ∆ ln , n=0 (17) Hϕ0 a E zn jωµ0 ∆2|n| − a2|n| ∆ 2|n| , |n| Hϕn ∆ + a2|n|
n = 0
(18)
where Z L is the impedance of the lines in the standard TLM formulation. For the case when n = 0, a cascade of the link line of impedance Z and a shortcircuit stub Z s is used. In this case, we obtain as the input impedance of this combination, √ Z in = jω µ0 ε0 ∆(Z + Z s ) (21) Equating this expression with (17) we obtain the required stub impedance, ∆ −Z (22) Z s = Z L ln a All these expressions can be brought into the 2D model of a wire node shown schematically in Fig. 9. The wire current is obtained directly from the stub current. It can be seen from this figure that modes n = 1 (cos ϕ and sin ϕ) and n = 2 with sample values of (0, 1, 0, −1), (1, 0, −1, 0), (1, −1, 1, −1) respectively at the node ports, do not see the stub (Vs = 0). Only mode n = 0 sees the stub. We have accounted fully in the model for modes n = 0, 1. Mode 2 should have a different impedance and phase constant, but for simplicity we do not account for this in the model as it does not appear to significantly degrade accuracy. It should be emphasized that the same calculation can be done without resort to an equivalent circuit. This simply requires the transform of mesh voltage quantities to modal quantities which after they see the correct modal impedances are transformed back into mesh quantities for propagation to the rest of the mesh. This indeed is the best way of proceeding when more complex configurations are studied. The outline of this calculation is as follows: Let us call [T ] the matrix transforming modal to mesh quantities and A the column vector of the four modal amplitudes. Then, the field quantities at the four ports of a 2D mesh are, E z = [T ]A,
Hϕ = [T ][Y ]A
(23)
where [Y] is a diagonal matrix with elements the modal admittances. From the normal TLM equivalences, Ez = V i + V r
,
Hϕ = Y L (V i − V r )
(24)
where in theses formulae ∆ is half the mesh size. The network is constructed by choosing link lines and stubs of suitable parameters so that the same impedance is seen as in [17, 18]. For the case n = 0, we choose a link line of impedance Z and phase constant k L which is terminated by a short-circuit. Employing (18) we obtain, jZ tan(k L ∆) = using
√ k L = ω µ0 ε0 ,
jωµ0 ∆2|n| − a2|n| |n| ∆2|n| + a2|n|
(19)
tan(k L ∆) k L ∆
the required impedance of the link line is obtained, Z=
Z L ∆2|n| − a2|n| |n| ∆2|n| + a2|n|
(20)
Fig. 9. Wire interface in 2D using the Modal Expansion Technique (MET) (after [43]).
106 C. Christopoulos: Multi-Scale Modelling in Time-Domain Electromagnetics Combining equation (23)(24), (Y L [I] + [Y ]) [T ]T V r = (Y L [I] − [Y ]) [T ]T V i
(25)
This equation relates the reflected to the incident voltages and represents the scattering process in the frequency domain. Matrix, (Y L [I] + [Y ])−1 (Y L [I] − [Y ]) is diagonal with elements the reflection coefficients (Z n − Z L ) /(Z L + Z n ) for each mode. Each of these coefficients represents a shift equal to, −e− j2k L ∆/n Fig. 11. Comparison of the MET with analytical results and other wire models (mesh resolution 50 mm, wire radius 5 mm). An incident field is scattered from the wire and the total field is observed 100mm in front of the wire.
for modes n = 1, 2, and e− j2k L ∆ for mode n = 0. This is shown schematically in Fig. 10 where for simplicity the shift is kept the same for mode n = 2. The accuracy of this model is shown in Fig. 11 [44], where in addition to comparisons with analytical results, the results from two other wire models are shown [39, 40]. The model uses no empirical factors and there are no restrictions on the size of the wire relative to the mesh size. The case of an offset wire, i.e. a wire not centred at the node can also be dealt with using this approach by exploiting the addition Theorems for Bessel Functions [45], e
jma
Z m (k R) =
∞
Jn (kr1 )Z n+m (kr01 )e
jnϕ
(26)
n=−∞
where all the relevant quantities are defined in Fig. 12. The field around the wire is expanded into cylindrical modes as before and then a mapping is implemented into the ports of the node using equation (26) [46]. The excellent behaviour of the model for the case of an offset wire is shown in Fig. 13. Many features of the model have been
Fig. 10. Schematic of the scattering and transformation of modal to total field quantities in a 2D mesh (after [44]).
Fig. 12. Definition of quantities used for the Bessel function addition Theorem.
implemented in 3D [47, 48]. It is worth stressing that results are obtained using closed form local solutions for the field around the wire, unlike the techniques described in [49, 50] where a full MoM numerical solution is necessary to calculate the wire current. The major advantage of the local solutions is simplicity and a minimal computational overhead.
Fig. 13. Modelling of an offset thin wire using the MET: wire radius = 0.1, offset is along the 45 degree line and the field in front of the wire is observed.
C. Christopoulos: Multi-Scale Modelling in Time-Domain Electromagnetics 107
3.1.2 Contour-path models
3.2 Models of thin panels with complex features
Contour path models have been used extensively in FDTD. The basic work is described in [49]. The following steps are implemented: • The calculation of an equivalent radius to replace a wire loom (a bundle) by a single equivalent wire • Faraday’s Law contour integral is applied to modify the FDTD algorithm and thus account for the presence of the wire • The equivalent electric and magnetic current sources are calculated on a virtual surface surrounding the wire • The EFIE and MoM are applied to calculate the wire current. A similar approach has been used to account for other features such as curved surfaces [51]. Further extensions of this technique may be found in [52, 53], where attention is given to end effects.
A problem similar in complexity and importance as that of fine wires is the description in numerical models of thin conducting panels containing a large number of holes, as are typically present in many equipment for cooling etc. Several possible problems can be envisaged e.g. thin conducting panels with/without holes, composite panels consisting of overlapping layers with different properties (carbon fibre composites), sandwiches of conducting and absorbing layers etc. In all these cases the challenging modelling feature is that the panel thickness and hole size are smaller than the nodal size. As explained in connection with wires, it is unprofitable to reduce the coarseness of the mesh to accommodate these features. We thus have another multiscale problem of the type already encountered. Early work in this area in connection with thin CFC panels was described in [57, 58]. In this approach onedimensional ladder networks connected between nodes were used to describe propagation through the panel. However the technique cannot deal with other features such as holes and lacks the generality and systematic parametrisation required in most problems. Taking the example of a perforated panel placed between nodes, it is clear that there will be both reflections and transmission through the panel. The relevant coefficients are in general frequency dependent and therefore scattering takes place described by a scattering matrix, R(ω) T(ω) [S(ω)] = T(ω) R(ω) (27)
3.1.3 Hybrid models As a general illustration of the hybridizing of methods to exploit their best features, we discuss briefly the hybridization of Integral Equation methods with TLM. An example is shown in [53, 54], where two separate regions with complex features are modelled by TLM and the intervening space, which is essentially free space, is modelled very efficiently by an IE method. However, we are interested here in the inverse process, where a wire inside a surface is modelled by an IE technique and the surrounding space by TLM. Effectively, this is another way of accounting for fine wires in TLM by hybridizing with a technique which is well suited to modelling wires. This approach is described in [55, 56]. The calculation proceeds as follows: • A surface SB is defined so that it completely surrounds the wire and is coincident with the boundaries of TLM cells • The total fields on SB are defined as a superposition of the incident fields from all sources located outside SB , and the fields which are radiated by the wire • The fields radiated by the wire, calculated at SB , are used to obtain the equivalent surface currents which are sources of radiation outside SB • The TLM technique, the equivalent current sources and other incident fields from sources outside SB are used to find the total fields outside SB • Inside SB , we enforce the condition that the tangential electric field (incident plus radiated) is zero at the wire surface. This conditions yields the current along the wire segments • From the current distribution along the wire, the radiated fields at the surface SB are obtained, which are then mapped onto the TLM pulses and superimposed on other pulses there to excite the total field outside SB . This process works satisfactorily but considerably more effort is involved compared to interfacing TLM to local solutions.
This scattering at the panel must now be included as part of the connection process in TLM, after it has been converted into a time-domain algorithm. This conversion results essentially to a digital filter algorithm which is embedded into the ‘connection’ process. The procedure is described in [60–62]. The basic steps are as follows: • Each coefficient F in (31) is put into standard Pade form NP
F(s) =
i=0 NP
bi s i (28) ai s i
i=0
where NP is the number of poles, decided from accuracy and computational cost considerations. • The frequency-domain Prony method [63] is employed to estimate ai and bi above from a set of NS measurements. After these coefficients have been estimated equation (28) is put into the form F(s) =
b NP (s − sz0 )(s − sz1 ) . . . (s − sz(NP−1) ) (s − s p0 )(s − s p1 ) . . . (s − s p(NP−1) )
where s pi , szi are the poles and zeroes.
(29)
108 C. Christopoulos: Multi-Scale Modelling in Time-Domain Electromagnetics • Equation (29) is converted into the discrete timedomain by the application of the bilinear transformation and after some manipulation we obtain NP
F(z) = B0 +
Bi z −i
i=1 NP
1+
(30) Ai
z −i
i=1
• As an illustration we assume that equation (30) represents a reflection coefficient, then F(z) = R(z) = k+1 V i /k V r
(31)
After some manipulation equation (31) can be put in a state space form, where the output equation is V i = B0 V r + B X
(32)
4. Discussion
and the state equation X = z −1 [A ]X + z −1 1 V
where B0 , B , [A ] are known and 1 = [1, 0 . . . 0]T . A signal flow diagram representing equations (32) and (33) is shown in Fig. 14. Typical results for the Shielding Effectiveness (SE) of a cabinet with a perforated wall are shown in Fig. 15. This simulation shows the accuracy to be expected and it is done at a fraction of the computational costs of a full-field model, where details of each perforation would have to be specifically included in the mesh. The process described is very general, as the characterization of the panel in the frequency domain (equation 27) can be done experimentally, analytically, or by a numerical model of a single perforation. Work is in progress to generalise this models to allow for arbitrary placement and to deal with many different generic objects that may be inserted into the mesh.
r
(33)
Fig. 14. Signal flow diagram of the ‘connection’ process in the presence of fine features using the Digital Filtering Interface (DFI).
Fig. 15. Shielding Effectiveness of an enclosure using the DFI and comparison with experiments: cabinet size (292 × 292 × 118mm3 ), area of perforations (212 × 78 mm3 ), hole radius 4 mm arranged in an equilateral triangle lattice of sides 10 mm. Simulations are shown for two different mesh resolutions (after [61]).
Several techniques developed to tackle multiscale problems were described. The number of publications, the variety of approaches and the period of time covered by the review testify to the importance and difficulty of this subject. It is not an understatement to say that without continuing progress in this area, it will be impossible to envisage electromagnetic computation in large practical systems of the type encountered in electromagnetic compatibility and signal integrity applications. Due to space limitations not all techniques could be covered in the same depth and detail. However, some conclusions are possible and should help point the way to future work. Multigridding schemes have a role to play but on their own cannot address fully the main modelling challenges. Some mesh refinement is desirable in certain areas and the modeller should have the facility of using one of the schemes described including the hybrid nodes. Even assuming that an entirely satisfactory multigrid scheme is developed, the effort of defining in detail the geometry of the problem in a very fine mesh appears to be out of proportion to any benefits, considering the limited resolution of the surrounding coarse mesh. It would be interesting to see developments which allow more flexibility in defining the structure of the mesh near particular features (as is possible in FEM) without sacrificing ease and simplicity of computation. Thin-wire nodes are highly desirable and effective in reducing complexity and are set to be an interesting and rewarding area of research. The Modal Expansion Technique (MET) appears to offer greater flexibility and accuracy compared to other approaches and could be generalised to deal with other types of problem, where modal expansion is the link between particular fine features and the surrounding coarse mesh. One can envisage MET acting as the ‘glue’ between the mesh and particular complex attachments to it. Without doubt, in addition to significant gains in run-time efficiency, major benefits are also possible in ease of problem definition. The treatment of essentially planar structures with complex details (e.g. perforated/composite panels) can be
C. Christopoulos: Multi-Scale Modelling in Time-Domain Electromagnetics 109
dealt with in a systematic and efficient way by the Digital Filtering Interface (DFI). The ease of characterization, the scope for automation and the possibility of developing a library of often used components and sub-assemblies, described in large-scale EM models by a small number of pre-calculated parameters as a DFI, is an exciting prospect. One can envisage work to expand this technique to encompass a wider range of components and structures. In conclusion, a long-term effort is required to develop further these techniques and to deal effectively with complexity reduction. In this context the pioneering work of Professor Russer should be mentioned [64, 65]. Acknowledgement. I wish to thank my colleagues Drs. P. Sewell, D.W.P. Thomas, J. Paul, V. Podlozny and Y.K. Choong for supplying me with material for this article and for their major contributions to my understanding of this topic. I also thank the EPSRC (UK) for support received over a number of years for multi-scale research. Thanks also to the Agency for Science Technology and Research (Singapore) and the Nanyang Technological University for their support and hospitality during the writing of this article while on attachment to NTU.
References [1] Taflove, A.: Computational electrodynamics: The Finite Difference Time-Domain Method. Norwood MA: Artech House, 1995. [2] Volakis, J.L.; Chaterjee, A.; Kempel, L.C.: Finite Element Method for Electromagnetics. New York: IEEE Press, 1998. [3] Harrington, R.F.: Field Computation by Moment Methods. New York: McMillan, 1968. [4] Christopoulos, C.: The Transmission-Line Modeling MethodTLM. New York: IEEE Press, 1995. [5] Kunz, K.S.; Simpson, L.: A technique for increasing resolution of finite-difference solution to maxwell’s equations. IEEE Trans. on EMC 23 (1981), 419–422. [6] Kim, I.S.; Hoefer, W.J.R.: A local mesh refinement algorithm for the FDTD method using maxwell’s curl equations. IEEE Trans. on MTT 38 (1990), 812–815. [7] Zivanovic, S.S.; Yee, K.S.; Mei, K.K.: A subgridding method for the time-domain finite-difference method to solve maxwell’s equations. IEEE Trans. on MTT 39 (1991), 471–479. [8] Prescott, D.T.; Shuley, N.V.: A method for incorporating different sized cells into the finite-difference time-domain analysis technique. IEEE Microwave and Guided Wave Letters 2 (1992), 434–436. [9] Okoniewski, M.; Okoniewska, E.; Stuchly, M.A.: Threedimensional subgridding algorithm for FDTD. IEEE Trans, on Antennas and Propagation 45 (1997), 422–429. [10] Chevalier, M.W.; Luebbers, R.J.; Cable, V.P.: FDTD local grid with material traverse. IEEE Trans. on Antennas and Propagation 45 (1997), 411–421. [11] Krishnaiah, K.M.; Railton, C.J.: Passive equivalent circuit of FDTD:AN application to subgridding. Electronics Letters 33 (1997), 1277–1278. [12] Krishnaiah, K.M.; Railton, C.J.: A stable subgridding algorithm and its application to eigenvalue problems. IEEE Trans. on Microwave Theory and Techniques 47 (1999), 620–628.
[13] Thoma, P.; Weiland, T.: A consistent subgridding scheme for the FDTD method. Int. J. Numer. Modeling 9 (1996), 359– 374. [14] White, M.J.; Yun, Z.; Iskander, M.F.: A new 3-D FDTD multigrid technique with dielectric traverse capabilities. IEEE Trans. on Microwave Theory and Techniques 49 (2001), 422– 430. [15] Herring, J.L.; Christopoulos, C.: Multigrid TLM method for solving EM field problems. Electronics Letters 27 (1991), 1794–1795. [16] Herring, J.L.; Christopoulos, C.: Solving EM field problems using a multiple grid TLM method. IEEE Trans. on Antennas and Propagation 42 (1994), 1654–1658. [17] Herring, J.L.; Christopoulos, C.: The application of different meshing techniques to EMC problems. – In: Proc. 9th Annual Rev. of Progress in Applied Computational Electromagnetics. Monterey, 1993, 755–762. [18] Wlodarczyk, J.: New multigrid interface for the TLM method. Electronics Letters 32 (1996), 1111–1112. [19] Trenkic, V.; Christopoulos, C.; Benson, T.M.: On the time step in hybrid symmetrical condensed TLM nodes. IEEE Trans. on Microwave Theory and Techniques 43 (1995), 2172–2174. [20] German, F.J.; Svigeli, J.A.; Mitra, R.: A numerical comparison of dispersion in irregular graded TLM and FDTD meshes. – In: 12th Annual Review of Progress in Applied Computational Electromagnetics. Monterey, 1996, 270–278. [21] van Bladel, J.: Singular electromagnetic fields and sources. New York: Oxford University Press, 1991. [22] Mexner, J.: The behaviour of electromagnetic fields at edges. IEEE Trans. on Antennas and Propagation 20 (1972), 302– 307. [23] Duffy, A.P.; Benson, T.M.; Christopoulos, C.: New method for accurate modelling of wires using TLM. Electronics Letters 29 (1993), 224–226. [24] Herring, J.L.; Hoefer, W.J.R.: Improved correction for 3-D TLM coarseness error. Electronics Letters 30 (1994), 1149– 1150. [25] Cascio, L.; Tardioli, G.; Rozzi, T.; Hoefer, W.J.R.: A quasistatic modification of TLM at the knife edge and 900 wedge singularities. IEEE Trans. on Microwave Theory and Techniques 44 (1996), 2519–2524. [26] Fujii, M.; Hoefer, W.J.R.: Field singularity correction in 2-D time-domain haar-wavelet modeling of waveguide components. IEEE Trans. on Microwave Theory and Techniques 49 (2001), 685–691. [27] Shorthouse, D.B.; Railton, C.: The incorporation of static field solutions into FDTD algorithms. IEEE Trans. on Microwave Theory and Techniques 40 (1992), 986–994. [28] Beilenhoff, K.; Heinrich, W.: Treatment of field singularities in the finite-difference approximation. Proc. of IEEE Symp. On Microwave Theory and Techniques, Atlanta, 1993. 979– 982. [29] Lindenmeier, S.; Russer, P.; Heinrich, W.: Hybrid dynamicstatic finite-difference approach for MMIC design. Proc. IEEE Symp. On Microwave Theory and Techniques, San Francisco, 1996. 197–200. [30] Lindenmeier, S.; Russer, P.; Heinrich, W.: Application of the hybrid dynamic-static finite-difference approach on 3DMMIC structures. Proc. of Progress on Applied Computational Electromagnetics, Monterey, 1997. 1182–1188. [31] Craddock, I.J.; Railton, C.J.: Stable inclusion of a priory knowledge of field behaviour in the FDTD algorithm:application to the analysis of microstrips lines. Proc. IEEE Symp. On Antennas and Propagation, Baltimore, 1996. 1300–1303.
110 C. Christopoulos: Multi-Scale Modelling in Time-Domain Electromagnetics [32] Esselle, K.P.; Okoniewski, M.; Stuchly, M.A.: Analysis of sharp metal edges at 45◦ to the FDTD grid. IEEE Microwave and Guided Wave Letters 9 (1999), 221–223. [33] Farina, M.; Rozzi, T.: Numerical investigation of the field and current behavior near lossy edges. IEEE Trans. on Microwave Theory and Techniques 49 (2001), 1355–1358. [34] Cui, T.J.; Chew, W.C.: Accurate analysis of wire structures from very-low frequency to microwave frequency. IEEE Trans. on Antennas and Propagation 50 (2002), 301–307. [35] Holland, R.; Simpson, J.W.: Finite-difference analysis of EMP coupling to thin struts and wires. IEEE Trans. on Electromagnetic Compatibility 23 (1981), 88–97. [36] Ledfelt, G.: A stable subcell model for arbitrarily oriented thin wires for the FDTD method. Int. J. Numer. Modelling 15 (to appear 2002). [37] Edelvik, F.: A new technique for accurate and stable modeling of arbitrarily oriented thin wires in the FDTD method. Uppsala University Technical Report 2002-016, www.it.uu.se/research/reports/2002-016/. [38] Berenger, J.-P.: A multiwire formulation for the FDTD method. IEEE Trans. on Electromagnetic Compatibility 42 (2000), 257–264. [39] Wlodarczyk, J.; Johns, D.P.: New wire interface for graded 3D TLM. Electronics Letters 28 (1992), 728–729. [40] Porti, J.A.; Morente, J.A.; Khalladi, M.; Gallego, A.: Comparison of thin wire models for the TLM method. Electronics Letters 28 (1992), 1910–1911. [41] Wlodarczyk, J.; Trenkic, V.; Scaramuzza, R.A.; Christopoulos, C.: A fully integrated multiconductor model for TLM. IEEE Trans. on Microwave Theory and Techniques 46 (1998), 2431–2437. [42] Dudley, D.G.: Mathematical Foundations of Electromagnetic Theory. New York: IEEE Press, 1994. [43] Choong, Y.K.; Sewell, P.; Christopoulos, C.: Accurate wire representation in numerical models for high-frequency simulation. Electronics Letters 37 (2001), 280–282. [44] Choong, Y.K.; Sewell, P.; Christopoulos, C.: New thin wire formulation for time-domain differential equation models. Int. J. Numer. Modelling 15 (2002 to appear). [45] Watson, G.N.: A Treatise on the Theory of Bessel Functions. Cambridge: CUP, 1958. [46] Choong, Y.K.; Sewell, P.; Christopoulos, C.: Accurate modelling of arbitrarily placed thin wire in a coarse mesh. Proc IEE Sci. Meas. Technnol. (to appear). [47] Choong, Y. K.; Sewell, P.; Christopoulos, C.: Hybrid node decsrption of thin wires based on analytical field representation in EMC problems. Proc. IEEE EMC Symp., Minneapolis, 2002. [48] Sewell, P.; Choong, Y.K.; Christopoulos, C.: An accurate thin wire model for 3D TLM simulations. IEEE Trans. on Electromagnetic Compatibility (accepted). [49] Umanshankar, K.R.; Taflove, A.; Bekker, B.: Calculation and experimental validation of induced currents on coupled wires in an arbitrarily shaped cavity. IEEE Trans. on Antennas and Propagation 35 (1987), 1248–1257. [50] Lindenmeier, S.; Christopoulos, C.; Russer, P.: Thin wire modelling with the TLM IE method. Proc. of Progress in Applied Computational Electromagnetics, Monterey, 2000. 587–593. [51] Jurgens, T.G.; Taflove, A.; Umanshankar, K.: FDTD modeling of curved surfaces validation of induced currents on coupled wires in an arbitrarily shaped cavity. IEEE Trans. on Antennas and Propagation 40 (1992), 357–366. [52] Douglas, M.; Okoniewski, M.; Stuchly, M.A.: Accurate modelling of thin wire antennas in the FDTD method. Microwave and Optical Technology Letters 4 (1999), 261–265.
[53] Makinen, R.M.; Juntunen, J.S.; Kivikovski, M.A.: An improved thin-wire model for the FDTD method. IEEE Trans. on Microwave Theory and Techniques 50 (2002), 1245–1255. [54] Lindenmeier, S.; Pierantoni, L.; Russer, P.: Hybrid space discretising-integral equation method for the numerical modelling of transient interference. IEEE Trans. on Electromagnetic Compatibility 41 (1999), 425–430. [55] Pierantoni, L.; Lindenmeier, S.; Russer, P.: Efficient analysis and modelling of the radiation of microstrip lines and patch antennas by the TLM-integral equation (TLM-IE) method. Int. J. Numer. Modelling 12 (1999), 329–340. [56] Lindenmeier, S.; Christopoulos, C.; Russer, P.: Thin wire modelling with the TLM-IE method. Proc. Progress in Applied Computational Electromagnetics, Monterey, 2000. 587– 593. [57] Lindenmeier, S.; Christopoulos, C.; Russer, P.: Methods for the modelling of thin wire structures with the TLM method. Proc. IEEE Symp. On Microwave Theory and Techniques, Boston, 2000. 387–390. [58] Trenkic, V.; Duffy, A.P.; Benson, T.M.; Christopoulos, C.: Numerical simulation of penetration and coupling using the TLM method. Proc. EMC Symp., Rome, 1994. 321–326. [59] Trenkic, V.; Christopoulos, C.; Benson, T.M.: Numerical simulation of polymers and other materials for electronic shielding applications. Proc. POLYMAT, London, 1994. 384– 387. [60] Paul, J.; Podlozny, V.; Thomas, D.W.P.; Christopoulos, C.: Time domain simulation of thin materials boundaries and thin panels using digital filters in TLM. Turkish J. of Electrical Eng. And Computer Sci. 10 (2002), 185–198. [61] Podlozny, V.; Paul, J.; Christopoulos, C.: Efficient calculation of shielding effectiveness of equipment cabinets in full-field numerical models. Proc. EMC Europe, Sorrento, 2000. [62] Paul, J.; Podlozny, V.; Christopoulos, C.: The use of digital filter techniques for the simulation of fine features in EMC problems solved in the time-domain. IEEE Trans. on Electromagnetic Compatibility (accepted). [63] Brittingham, J.N.; Miller, E.K.; Williams, J.L.: Pole extraction from real frequency information. Proc. IEEE 68 (1980), 263–273. [64] Russer, P.; Cangellaris, A.: Network oriented modeling, complexity reduction and system identification techniques for electromagnetic systems. Proc. 4th Int. Workshop on Computational electromagnetics in the Time-Domain: TLM/FDTD and Related Techniques, Nottingham, 2001. 105–122. [65] Russer, P.: Network oriented modeling of radiating electromagnetic structures. Turkish J. Elec. Engin. 10 (2002), 147– 162.
Christos Christopoulos was born in Patras, Greece in 1946. He received his first degree from the National Technical University of Athens and his MSc and DPhil from Sussex University. He is currently Professor of Electrical Engineering at the University of Nottingham. His interests include, computational electromagnetics, electromagnetic compatibility and the protection of large networks using fast transients. He is the author of more than 250 scientific articles and five books. He is current Chairman of the IEE Professional Network on EMC and an Associate Editor of the IEEE EMC Transactions.