## Abstract

We review non-conformal domain decomposition methods (DDMs) and their applications in solving electrically large and multi-scale electromagnetic (EM) radiation and scattering problems. In particular, a finite-element DDM, together with a finite-element tearing and interconnecting (FETI)-like algorithm, incorporating *Robin transmission conditions* and *an edge corner penalty term*, are discussed in detail. We address in full the formulations, and subsequently, their applications to problems with significant amounts of repetitions. The non-conformal DDM approach has also been extended into surface integral equation methods. We elucidate a non-conformal integral equation domain decomposition method and a generalized combined field integral equation method for modelling EM wave scattering from non-penetrable and penetrable targets, respectively. Moreover, a plane wave scattering from a composite mockup fighter jet has been simulated using the newly developed multi-solver domain decomposition method.

## 1. Introduction

The accurate simulation of electrically large and geometrically complicated electromagnetic (EM) problems is of vital importance in many engineering applications. The scopes of traditional approaches, such as the finite-element method (FEM), the boundary-element method (BEM) and the finite-difference method, are limited to problems of moderate electrical size and simplified complexity. These limitations stem from the vast computational resources required by conventional numerical methods. Recently, a non-conformal domain decomposition method (DDM) has been proposed (Lee *et al.* 2005; Vouvakis *et al.* 2006) and has emerged as a promising technique for the numerically rigorous solution of Maxwell’s equations. This paper documents the extensions of the non-conformal DDM to solutions of EM radiation and scattering problems involving composite objects. The central idea of the DDM can be summarized as follows. Once the solution of a subregion is obtained, its Huygens’ equivalent sources (i.e. electric and magnetic current sources) are computed. These sources generate incident fields that subsequently impinge on all other subregions and make each subregion a well-defined scattering problem. This process continues until no significant information is exchanged between subregions.

Lee *et al.* (2005) and Vouvakis *et al.* (2006) have exploited a finite-element-based non-overlapping DDM to solve large- and multi-scale EM radiation and scattering problems with significant repetitions. The proposed non-conformal DDM enables each of the subdomains to be discretized independently, and thus greatly simplify the mesh generations for problems with multi-scaled features. Additionally, for many engineering applications, usually only small portions of the device geometry will be modified/changed during the design and optimization processes. The possibility of only remeshing the parts that have been modified without consideration of the rest of the problem geometry will then be highly desirable. Furthermore, an integral equation domain decomposition method (IE-DDM; Peng *et al.* (2011)) has been proposed for the computations of EM scatterings from non-penetrable objects. Instead of solving the original scattering problem as one large computational domain, the IE-DDM algorithm decomposes the perfect electric conductor (PEC) target into many smaller and non-overlapping subdomains. Each local subdomain is enclosed by a closed surface and solved individually through the combined field integral equation with excitations that include radiations from all the other subdomains. Subsequently, these subdomains are coupled to one another via the multi-level fast multipole algorithm (MLFMA) using the Stratton–Chu representation formulas for non-touching subdomains. The continuities of tangential fields, both electric and magnetic fields, on the interfaces between touching subdomains are implicitly enforced via Robin-type transmission conditions (TCs).

## 2. Non-conformal finite-element domain decomposition methods

Consider a finite three-dimensional region with its exterior boundary ∂*Ω*. The region *Ω* is further partitioned into *M* non-overlapping subdomains such that *Ω*=*Ω*_{1}∪*Ω*_{2}⋯∪*Ω*_{M}. They are mainly of cubic or parallelepiped shapes (figure 1), and the interfaces between subdomains are assumed to be planar. Moreover, we denote by and the interfaces between *Ω*_{m} and *Ω*_{n}, with the surface seen from *Ω*_{m} and the one from *Ω*_{n}. Furthermore, we associate with each subdomain interface a unit normal , which points from subdomain *Ω*_{m} towards subdomain *Ω*_{n}. Additionally, the exterior boundary of subdomain *Ω*_{m} is denoted by . We write the outward-directed unit normal to ∂*Ω*_{m} as , and the subscript *m*∈{1,2,…,*M*} connotes the restriction of a quantity to *Ω*_{m}.

We will also use the following tangential trace and twisted tangential trace operators on subdomain interface : 2.1

### (a) Boundary-value problem and interior penalty formulation

We designate the time-harmonic electric and magnetic fields by **E** and **H**, respectively. stands for the free-space wavenumber, where *ω*=2*πf* is the radial frequency of operation with frequency *f* in Hertz, and *ε*_{0} and *μ*_{0} represent the permittivity and permeability in free space. In a material region, the wavenumber is given by , where *ε* and *μ* are the permittivity and permeability of the material. We also define *ε*_{r,m}=*ε*_{m}/*ε*_{0} and *μ*_{r,m}=*μ*_{m}/*μ*_{0} as the relative permittivity and permeability of subdomain *Ω*_{m}. The intrinsic impedance in free space is given by , the relative wave impedance of subdomain *Ω*_{m} is subsequently defined as , and *λ*_{0} is the wavelength in free space. We partition the entire domain *Ω* into *M* subdomains, namely *Ω*=*Ω*_{1}∪*Ω*_{2}∪⋯*Ω*_{M}. Consequently, within each subdomain, the electric field **E**_{m}=**E**|_{Ωm} must satisfy the following vector wave equation:
2.2where and is the given impressed current density in *Ω*_{m}.

The trial function space that we have employed in this work is
2.3Subsequently, we have **E**_{m}=**E**|_{Ωm}∈**H**(curl;*Ω*_{m}). For any trial field **u**_{m}∈**H**(curl;*Ω*_{m}), there is associated with it a volume residual from equation (2.2),
2.4Proper dual pairing requires that the volume residual (2.4) is tested with **v**_{m}∈**H**(curl;*Ω*_{m}). We shall write the resulting volume penalty term for any test field **v**_{m}∈**H**(curl;*Ω*_{m}) as (after applying the Green identity)
2.5where the volume and surface integral terms are
2.6In (2.5), is the set that includes all subdomain interfaces in *Ω*_{m}. Moreover, the boundary of *Ω*_{m}, ∂*Ω*_{m}, has been decomposed into the exterior boundary part () and the domain interface part (). Furthermore, we count as two members of , and , for the interface(s) between *Ω*_{m} and *Ω*_{n}.

We shall omit our discussions on the treatment of the exterior boundary ∂*Ω* (or ) since it can be addressed using well-documented approaches such as the absorbing boundary conditions (Stupfel 1994) and the BEMs (Stupfel 2001). Therefore, from here onwards, we shall drop the term from the formulation. Finally, extending (2.5) to all the subdomains results in
2.7Additional notations are introduced in (2.7), which require further explanation. The set includes all the interfaces between subdomains. In addition, we have the obvious notations in (2.7),
2.8

The decomposed problem will not have the same solution as the original entire domain problem unless proper boundary conditions are imposed for every subdomain interface . Directly, the boundary conditions are
2.9which simply state that the tangential components of the electric and magnetic fields must be continuous across every subdomain interface . However, as pointed out repeatedly by many authors, the direct enforcements of (2.9) may not be profitable (Despres *et al.* 1992; Collino *et al.* 1997; Lee *et al.* 2005; Li & Jin 2007). Among the difficulties that might be encountered, we mention: (i) the subdomain problem may suffer the notorious internal resonances and (ii) the overall DDM system matrix could be ill conditioned and may suffer from slow convergence or even fail to converge during the DDM iteration.

In an effort to accelerate the convergence in the DDM iteration, Robin TCs (Despres *et al.* 1992) have been proposed to replace (2.9),
2.10Straightforward algebraic manipulations reveal that (2.10) and (2.9) are equivalent. To facilitate the implementation of the Robin TCs (2.10), we introduce an additional set of variables, similar to the implementation of the Robin TCs detailed in Zhao *et al.* (2007) and Peng *et al.* (2010), for each of the subdomain interfaces ,
2.11Namely, **f**_{mn} is a tangentially continuous vector field on . Substituting (2.11) into (2.7) results in
2.12Moreover, the Robin TCs can now be expressed as
2.13

Clearly, from equation (2.13), we see that for every trial pair fields and , we have the following residuals:
2.14with
2.15We shall test with and scale the equation by *k*_{0}. Note that , the same discontinuous curl-conforming function space on for **f**_{mn}. Therefore, one surface penalty term is added for each subdomain interface ,
2.16Extending the surface penalty term to include all the subdomain interfaces results in the following expression:
2.17Finally, we combine both the volume and the surface penalty terms to formally state the following Galerkin weak formulation.

Find ; with , such that
2.182.19
∀**v**={**v**_{m},*m*=1,2,…,*M*}, with **v**_{m}∈**H**(curl,*Ω*_{m}); , with .

### (b) Divergence-free condition and corner-edge penalty terms

Through eigenvalue analysis, Peng & Lee (2011) have shown that the DDM system matrix equation is singular (nearly singular) for conformal (non-conformal) meshes. Moreover, they have illustrated that it is mainly due to the redundancy of the cement variables (Lee *et al.* 2005) on the corner edges. For example, shows figure 2 a corner edge that is shared by four subdomains *Ω*_{1}, *Ω*_{2}, *Ω*_{3} and *Ω*_{4}. In total, there would be eight cement variables associated with the corner edge. As a consequence of this redundancy, there exist eigenmodes that violate the divergence-free condition, ∇⋅**B**=0, on the corner edges such as the one in figure 2. The corresponding eigenvalues of these eigenmodes are zeros and very small numbers for conformal and non-conformal DDMs, respectively. To mitigate this technical difficulty, a corner-edge penalty method was proposed in Peng & Lee (2011). Take for example, the corner edge in figure 2, the divergence-free condition ∇⋅**B**=0 will result in the following equation:
2.20Since the cement variable **f**_{mn} is closely related to (differing only by a scalar constant) on , we can rewrite (2.20) in terms of the cement variables as
2.21We shall denote **F**_{l} for each of the corner edges as the sum of all the cement variables associated with it. Specifically, for the corner edge shown in figure 2, we have
2.22Consequently, we suggest the addition to (2.18) of a corner-edge penalty term to enforce the required divergence-free condition for each of the corner edges. For the corner edge in figure 2, we have
2.23With the corner-edge penalty terms, the Galerkin weak statement for non-conformal DDM can be restated as the following.

Find **u**={**u**_{m},*m*=1,2,…,*M*}, with **u**_{m}∈**H**(curl,*Ω*_{m}); , with , such that
2.24∀**v**={**v**_{m},*m*=1,2,…,*M*}, with **v**_{m}∈**H**(curl,*Ω*_{m}); , with , where *C*_{l} stands for the set of all corner edges.

### (c) Finite-dimensional discretization

Let denote the tetrahedral mesh (with characteristic edge length *h*) of subdomain *Ω*_{m}, and and are the surface triangulations induced on and , respectively. The set of corner edges is collected as *C*_{l}. On each of the subdomains, we define discrete trial and test functions, and , respectively, by . Here, is taken to be the space of mixed-order curl-conforming vector basis functions defined in Sun *et al.* (2001), with order *p*=2 (with 20 vector basis functions within each tetrahedron). On the domain interfaces, we also define discrete trial and test functions, and , respectively, with . is taken as the two-dimensional tangential component trace of , but the degrees of freedom (d.f.) are assigned independently for each subdomain interface. Thus, it is tangentially continuous on each of the domain interfaces, but not tangentially continuous across the corner edges.

In summary, we have the finite-dimensional discretization of the proposed non-conformal DDM with the Robin TCs and the corner-edge penalty terms.

Find with ; with , such that 2.25 with ; with .

The matrix equation resulting from the finite-dimensional discretization (2.25) can be written explicitly as the following (using *M*=2 as an example):
2.26In (2.26), all the sub-matrices and the RHSs can be easily identified from the discrete Galerkin formulation (2.25). Moreover, the unknowns within subdomain *Ω*_{m} are divided into three groups: for the coefficients of the vector field **u**_{m} inside the subdomain; for the coefficients of **u**_{m} on the domain interface ; and for the coefficient vector of all the cement variables **f**_{mn} on .

To facilitate our discussions on the matrix solution procedure, we rewrite equation (2.26) in a compact form as
2.27with
2.28To solve the system matrix equation (2.27), we first apply a block-diagonal preconditioner, which is similar to the additive Schwarz procedure (Cai & Sarkis 1999), using only the inverse of the subdomain matrices (either by forward/back substitution or by an iterative generalized conjugate residual (GCR) Krylov method),
2.29We shall introduce a simple restriction operator for a coefficient vector **u**_{m}. Explicitly, we have
2.30Furthermore, it is easy to show that
2.31Similarly, we have
2.32Accordingly, we propose to solve the system matrix equation of the DDM (2.27) in two steps.

—

*Solving the surface unknowns.*With the help of (2.31) and (2.32), it is not difficult to derive the following matrix equation in terms of only the surface unknowns: 2.33By extending (2.33) to*M*subdomains, we have 2.34Of course, in (2.34), if subdomains*Ω*_{m}and*Ω*_{n}do not touch each other, then .—

*Recovering the subdomain solutions.*Once the surface unknowns ,*m*=1,2,…,*M*, are computed, the interior solution for each subdomain can be recovered via 2.35

### (d) Numerical example: a three-dimensional complex package benchmark

Recently, as pointed out by Dolean *et al.* (2009), Rawat & Lee (2010) and Peng & Lee (2011), the use of Robin TCs may not be effective for modelling complex multi-scale EM problems. Consequently, a group of higher order TCs, notably the optimal TC (Dolean *et al.* 2009) and the second-order TC (Rawat & Lee 2010; Peng & Lee 2011) were proposed to further improve the convergence in DDM iterations. Although, the higher order TCs are understandably more elaborate than the Robin TCs discussed here, the weak Galerkin statement and the finite-dimensional discretization can be formulated in much the same way. For details, we refer the interested readers to Rawat & Lee (2010) and Peng & Lee (2011).

The numerical results reported were computed using the non-conformal finite-element DDM with second-order TCs on a workstation with two quad-core 64 bit Intel Xeon X5450 CPUs and 48 GB of RAM. Double-precision arithmetic is used throughout the numerical studies. To demonstrate the potentials of the proposed approach, we applied the non-conformal finite-element DDM to study a real-life three-dimensional package benchmark example provided by IBM (Krauter *et al.* 2006; Shao *et al.* 2011). This package example includes an eight-layer structure: Ground/Mounting Pads (SURFACE), Signal (FC3), Ground (FC2), Signal (FC1), Signal (BC1), Power (BC2), Signal (BC3) and Ground/Mounting Pads (BASE). With this kind of stackup technique, routing layers are buried between metallic planes and are therefore shielded. Our target is to study the signal integrity of the highlighted 20 signal traces with all geometrical details. Subsequently, we cut the interested portion as in figure 3 and analyse the reduced model with dimensions 10.5 × 16.5 mm. We further partition the entire computational domain into 149 subdomains, as in figure 4. In the performed study, the meshes on the interfaces are non-matching and the adjacent subdomains are also geometrically non-conformal. In figure 5, we compare the computed *S*_{21} parameters against the measurements from DC to 10 GHz. As can be seen from the figure, very good agreement exists between the two results. Additionally, in figure 5, we plot the comparisons between the transient waveforms of the computed and measured results. Again, we observe overall good agreement. In particular, the computed transient output gives a delay time of 120.2 ps, showing excellent agreement with the 118.6 ps of the measurement.

### (e) Finite-element tearing and the interconnecting algorithm for structures with repetitions

Many practical EM applications involve many repetitions, examples include large finite antenna arrays, frequency selective surfaces, photonic band gap structures, metamaterials and graphene materials, to name just a few. In a straightforward manner, the solution of equation (2.34) would require solution of every subdomain problem for each DDM iteration. With significant repetitions, the majority of the subdomains are identical, and hence, would also have the same subdomain system matrix. Therefore, an elegant way to benefit from the overwhelming repetitions existing in the structure is to rewrite equation (2.34) even further,
2.36To emphasize the effect of repetitions, we have assumed in equation (2.36) that all the subdomains are the same. Let , which carries the physical significance of a *numerical Green function*. By computing the dense matrix in the preprocessing stage, we have it explicitly available during the DDM iteration.

As pointed out previously (Vouvakis *et al.* 2006), the matrix can be viewed as the numerical Green function of the subdomain. Subsequently, it possesses a number of interesting properties similar to the impedance matrix resulting from the BEM, such as symmetry for subdomains with only reciprocal materials and rank deficiency for off-diagonal coupling blocks. For general domain partitions, the finite-element tearing and interconnecting (FETI)-like algorithm may not be beneficial; however, in the case of repeating subdomains, the computational resources can be reduced drastically and consequently outperform the classical DDM algorithms. Furthermore, since the numerical Green function, the dense matrix, exhibits low-rank properties for many of its matrix blocks, the data sparse representation techniques, such as the adaptive cross approximation (Bebendorf 2000; Zhao *et al.* 2005) can be applied to further improve the computational efficiency. Figure 6 shows a sample example of applying the FETI-like algorithm to solve a large Vivaldi antenna array, with 548 Vivaldi antenna elements (table 1). As indicated in the figure, it took more than 8 h to compute the numerical Green function ; however, this process does not need to be repeated for computing the scanning or for different configurations of the Vivaldi antenna arrays. For any subsequent computations, whether it is the scanning scenario or arrays with different number of Vivaldi elements, the matrix can be recycled. Additionally, as shown in the figure, for the Vivaldi antenna array with 548 elements, it results in more than 360 million FEM unknowns. However, with the application of the FETI-like algorithm to take advantage of the repetitions, it took only 1.8 GB of memory.

## 3. Non-conformal integral equation domain decomposition methods for non-penetrable targets

Recently, the non-conformal domain decomposition approach has also been extended to surface integral equation (SIE) methods, particularly for computing EM wave scatterings from non-penetrable and composite targets (Peng *et al.* 2011; submitted *b*). Consider an EM wave scattering from a mockup PEC jet aircraft at the X-band, shown in figure 7. The dimensions of the jet aircraft are approximately 14.5 m long, 9.4 m wide and 4.5 m high. The proposed non-conformal decomposition scheme follows a hierarchical domain partitioning strategy, as illustrated in the following.

—

*Level 1. Geometrical partitioning.*The entire aircraft is first divided into 13 closed-surface groups, as depicted in figure 7, based on the computer-aided design-generated geometry employed in industrial applications.—

*Level 2. Computational partitioning based on local features.*Each of the 13 groups in figure 7 is further decomposed into smaller subdomains according to local geometrical complexities and available computational resources. For instance, the group*left wing*consists of seven subdomains, as shown in figure 8. The original model is finally partitioned into 45 subdomains, and everyone of them is a closed-surfaced PEC scatter.

We consider EM scattering from a non-penetrable target occupying a finite domain . Let **j**=−(1/*ik*_{0})*π*^{×}(∇×**E**)∈**H**(div_{Γ},∂*Ω*), the application of the combined field integral equation (CFIE) results in
3.1with
3.2In (3.2), p.v. stands for principal value evaluation, and *g*(**x**,**y**)=e^{−ik0|x−y|}/4*π*|**x**−**y**| is the free-space Green function, with **x** and **y** denoting the observation and the source points, respectively.

In figure 9, we show a sample decomposition of a domain (target), for simplicity, into four subdomains such that *Ω*=*Ω*_{1}∪*Ω*_{2}⋯∪*Ω*_{4}. Furthermore, we denote the interface between subdomains by *Γ*_{mn}=∂*Ω*_{m}∩∂*Ω*_{n}, with its unit normal and exterior boundaries as . We shall also designate the outward directed unit normal to as . Through this decomposition, each subdomain is a metallic scatter with a closed conducting surface. For instance, the boundary surface ∂*Ω*_{1} of subdomain *Ω*_{1} includes exterior boundary and interfaces *Γ*_{12} and *Γ*_{14}, such that . The next step is to derive an integral representation for each of the subdomains. Differently from a single scatterer in free space, the incident fields of subdomain *Ω*_{m}, and , need to incorporate the scattered fields contributed from the other three subdomains. For example, we express the incident fields and for subdomain *Ω*_{1} as
3.3Similar expressions of the incident fields existed for subdomains *Ω*_{2}, *Ω*_{3} and *Ω*_{4}. Subsequently, the CFIE formulation for subdomain *Ω*_{m}, *m*∈{1,2,3,4}, can be written as
3.4Straightforwardly combining (3.4) for *m*=1,2,3,4 results in a general CFIE formulation to solve for the decomposed problem. However, the numerical integrations where the source and observation points are located on the same interface between two touching subdomains can be further simplified.

Let us consider the incident fields on ∂*Ω*_{1} as stated in (3.3). For clarity, we shall rewrite (3.3) as
3.5However, in (3.5), the integrals , **x**∈*Γ*_{12}, and , **x**∈*Γ*_{14}, are somewhat troublesome since they are both integrals where the observation points and the excitation currents are located on the same artificial interfaces between touching subdomains. We propose, instead, the following TCs on *Γ*_{12} and *Γ*_{14}:
3.6It can be shown that (3.6) enforces the required field continuities across the interface between touching subdomains.

By testing (3.4), via dual-pairing (Houston *et al.* 2005; Peng *et al.* 2011), with div-conforming functions *λ*_{m}∈**H**^{−1/2}(div_{τ},∂*Ω*_{m}), *m*∈{1,2,3,4}, we have
3.7The RHS of (3.7) can be computed through the TCs on ∂*Ω*_{m}. Specifically, we have, for *m*=1, the following substitutions:
3.8The RHSs for the other three subdomains can be computed in much the same way.

We proceed by discretizing each subdomain independently with a triangular mesh, denoted as . The surface meshes on the exterior boundary are denoted by . On each of the subdomains, we define discrete trial and test functions, and , respectively, with . Here, is taken to be the space of surface div-conforming vector Rao–Wilton–Glisson (RWG) basis functions (Rao *et al.* 1982) over . Subsequently, after expanding the electric current in terms of the basis functions, the discrete system can be cast into a matrix equation of the following form with the coupling parameter :
3.9with
In writing (3.9), we have further classified the surface unknowns within each subdomain ∂*Ω*_{m} into two groups: the unknowns on the exterior boundaries and the unknowns on the artificial domain interfaces *Γ*. As expected, the subdomain matrix *A*_{m} is the usual CFIE matrix for subdomain *Ω*_{m}, *m*=1,2,3,4, and
Since the triangular meshes on and are non-conformal, the integrations of are implemented on a union mesh of two non-conformal triangulations to perform the numerical integrations accurately. We refer the interested readers to Peng & Lee (2011) and Xu & Hoppe (2005) for details.

A Krylov subspace iterative method, the GCR (Eisenstat *et al.* 1983), is adopted for the solution of (3.9). The convergence criteria for the GCR solver is defined as
3.10The surface mesh of the decomposed problem corresponds to 13 859 548 triangles and results in 22 195 872 d.f. For the simulation at the X-band, the IE-DDM requires nine iterations to reach *ϵ*=10^{−2}. Moreover, we plot in figure 10 the electric current distributions on the exterior surfaces of the PEC aircraft. From these figures, we have witnessed smooth current distributions, without noticeable discontinuities, across domain interfaces.

## 4. Generalized combined field integral equation method for penetrable targets

In §3, we illustrated the application of the IE-DDM algorithm using a mockup PEC jet aircraft as an example. However, *real* jet aircrafts are often made of composite metallic and dielectric structures, as shown in figure 11. Note that in figure 11, we have represented different material regions with different colours: grey for the PEC airframe, the dielectric radome at the nose is coloured cyan, and the magenta colour highlights the lossy thin coatings on the left and right wings, as well as vertical and horizontal stabilizers. To compute the full-wave solution of an EM wave scattering from such a composite aircraft effectively, we shall introduce an additional subdomain solver, a generalized combined field integral equation (GCFIE) method, and integrate it into the IE-DDM formulation. We remark that by adding the GCFIE solver into the IE-DDM formulation, we naturally form the multi-solver DDM (MS-DDM; Wang *et al.* 2011).

The radome region located at the nose of the aircraft is a dielectric shell with uniform thickness of 10 mm and a relative permittivity of 3.5. To compute EM wave scattering from a homogeneous penetrable target, the SIE method (Harrington 1989; Nédélec 2001) is an appealing full-wave approach. Indeed, quite a few SIE formulations exist in the literature (Harrington 1989; Nédélec 2001; Buffa *et al.* 2003; Yla-Oijala & Taskinen 2007). Herein, we summarize a new GCFIE formulation proposed by Peng *et al.* (submitted *a*) for modelling the EM scattering from the dielectric radome region. We consider a homogeneous dielectric sub-region *Ω*_{m} with closed surface ∂*Ω*_{m}, as illustrated in figure 12. We further introduce the superscripts − and +, which tag traces onto sub-region surface ∂*Ω*_{m} from the interior and the exterior regions, respectively. The surface unit normal on surface points from *Ω*_{m} to the exterior air region, and on surface points inward into sub-region *Ω*_{m}. The GCFIE formulation starts by defining four auxiliary variables, and , that represent the tangential electric field and (scaled) electric current on the exterior side of sub-region surface and on the interior side of sub-region surface, ,
4.1and
4.2

By incorporating multiple equivalent electric fields and electric currents on surfaces and , the entire can be decomposed into two almost decoupled regions: the interior region *Ω*_{m} and the exterior region *Ω*_{ext}, as shown in figure 12. For the exterior air region, we have the electric field integral equation and the magnetic field integral equation on as follows:
4.3and
4.4We shall extend our definition of the combined field integral operator to account for the material properties, namely
4.5where 0≤*α*≤1,
4.6with
4.7Subsequently, the GCFIE formulation for the exterior and the interior regions can be stated formally as
4.8and
4.9Naturally (4.9) applies the Green function of the homogeneous interior region *Ω*_{m}. Finally, the interior traces and the exterior traces are coupled together through the boundary conditions or TCs. In particular, in our current approach, Robin-type TCs are used to weakly enforce the required field continuities across the interfaces between the interior and exterior regions. Written explicitly, we have
4.10

To derive a weak formulation for (4.8)–(4.10), we again employ the interior penalty approach. Consider the following residuals:
4.11These residuals can be interpreted as the physical error currents that support the differences between the exact field solution and that obtained through the numerical method. Following the Galerkin procedure, we pair each of these residuals with proper test functions. We argue that the rotated magnetic current should be paired with a rotated magnetic field to give rise to a surface energy density of the form **H**⋅**M**. Similarly, we pair the electric current with a surface electric field to give the familiar surface energy density of the form **E**⋅**J**.

Additionally, we should also remark on one fairly unusual aspect of our proposed formulation. That is, for each of the residuals in (4.11), we test it twice (with different scalings as indicated later) with both surface curl- and div-conforming basis functions. Therefore, we have a linear combination of the weighted residuals as
4.12Naturally, , and we insert it into (4.12) to emphasize the physical significance of the proposed dual pairings. Note that the parameters *c*_{1}, *c*_{2}, *c*_{3}, *c*_{4} *c*_{5}, *c*_{6}, *c*_{7} and *c*_{8} should not be zero, other than that, they can be chosen freely. We use *c*_{1}=*c*_{5}=−1, *c*_{2}=*c*_{6}=1, and *α*=0.5. Therefore, we obtain the desired Galerkin weak statement as the following.

Seek and , such that 4.13In the discrete implementations, we have separate triangulations for the interior and the exterior surfaces. Subsequently, on each of the triangular meshes, we construct finite-dimensional discrete trial and test function spaces and . Moreover, we have and . Here, are taken to be the space of surface div-conforming vector RWG basis functions, and are the space of first-order curl-conforming vector basis functions. After expanding the electric field and electric current in terms of the corresponding basis functions, the discrete system results in a matrix equation of the following form: 4.14with The other sub-matrices and have similar expressions and can be easily identified by comparing (4.14) with the weak Galerkin statement in (4.13).

We applied the GCFIE to model a plane-wave scattering from the dielectric radome, with thickness 10 mm and relative permittivity *ε*_{r}=3.5 of the mockup composite jet aircraft shown in figure 11. The length of the radome is roughly 1.2 m, which corresponds to approximately 40 free-space wavelengths at 10 GHz. The surface current distributions, magnitude of the real part of **j**, are plotted in figure 13 for a plane wave incident from the nose at three different frequencies: 75 MHz, 1 GHz and 10 GHz. Moreover, the computational statistics at these three frequencies are summarized in table 2. Judging from table 2, we remark that the proposed GCFIE performs satisfactorily for different frequencies, even with quite different discretization sizes (0.001, 0.05 and 0.08 *λ*_{0}). We should also remark here, as can be seen from table 2, that the resources, both memory and CPU times, required for these three solutions were modest. This is mainly attributed to the application of the MLFMA (Song & Chew 1995) to reduce the required memory and to speed up matrix vector multiplications. Furthermore, we note that the computations at 75 MHz were much more costlier than the solution at 1 GHz, even though the number of unknowns is less at 75 MHz than the 1 GHz case. As pointed out by many authors in the literature (Andriulli *et al.* 2008), the performance of the MLFMA deteriorates when it is applied to problems with very small mesh sizes (0.001 free-space wavelength at 75 MHz).

## 5. Multi-solver domain decomposition method

Previously, a MS-DDM was proposed by Zhao *et al.* (2007) for well-separated targets, which has been extended by Wang *et al.* (2011) for treating multi-regional EM problems with touching subdomains. Specifically, in Wang *et al.* (2011), the MS-DDM has been applied to compute the isolations among multiple antennae mounted on a large air platform at the X-band. The fundamental strategy of the proposed MS-DDM is to decompose the entire computational domain into many subregions based on the local material properties and geometrical features. Subsequently, we employ the most suitable computational electromagnetic (CEM) technique for each of the subregions. Moreover, the coupling between well-separated subregions is implemented through Stratton–Chu representation formulae. However, for the touching interfaces between neighbouring subregions, a Robin TC is introduced to weakly enforce the required field continuities across domain interfaces. The proposed MS-DDM is therefore well suited for modelling multi-scale and electrically large EM radiation and scattering problems. Additionally, by using the proposed MS-DDM framework, multiple existing CEM solvers can be integrated with few minor modifications.

Presently, with GCFIE formulation readily available for modelling dielectric radome and partially coated PEC targets, it is natural to apply the MS-DDM technique to integrate them and compute EM scatterings from the composite jet aircraft in figure 11. The plane wave incidents upon the aircraft from the nose (radome) with the electric field polarized in the direction. The simulations are conducted at the X-band, or *f*=10 GHz to be specific. The material properties of the lossy thin coatings are assumed to be *ε*_{r}=15.0−*i*0.03 and *μ*_{r}=2.2−*i*0.7. Moreover, at 10 GHz, the dimensions of the mockup aircraft are 500*λ*_{0}×314*λ*_{0}×151*λ*_{0} (figures 12 and 13).

Figure 14*a* shows the computed surface current distribution on the composite aircraft model. As can be seen from the figure, with lossy thin coatings, the surface currents are obviously reduced on the leading edges of the left and right wings of the composite aircraft. Similar effects can be observed on the left and right tail stabilizers, where thin lossy materials were employed on the composite model. Figure 14*b* provides more detailed views of the comparisons between the surface current distributions on the composite and the PEC aircrafts. Finally, we elucidate in table 3 the computational statistics of the two simulations, the mockup PEC jet aircraft and the mockup composite jet aircraft.

## 6. Conclusions

A large Vivaldi antenna array, with 548 Vivaldi antenna elements, was analysed using a non-conformal finite-element DDM, which successfully incorporates Robin TCs and an edge-corner penalty term to achieve robust and speedy convergence in DDM iterations. Moreover, we presented modelling of EM scattering from a multi-scale conducting fighter jet, at the X-band, using a newly developed integral equation non-conformal DDM. Once again, the non-conformal nature of the proposed IE-DDM makes it possible to address systematically the inherent multi-scale feature of the complex target. Finally, by introducing a GCFIE method as a subdomain CEM solver, and integrating it into the IE-DDM to form naturally the multi-solver DDM framework, we are able to compute the full-wave solution of a plane wave scattering from a composite fighter jet at the X-band. Through electrically large and multi-scale EM wave radiation and scattering problems, the astounding potential of non-conformal DDMs has been fully illustrated.

- Received January 17, 2012.
- Accepted March 6, 2012.

- This journal is © 2012 The Royal Society