An efficient high order plane wave time domain algorithm for transient electromagnetic scattering analysis

An efficient high order plane wave time domain algorithm for transient electromagnetic scattering analysis

Engineering Analysis with Boundary Elements 85 (2017) 13–19 Contents lists available at ScienceDirect Engineering Analysis with Boundary Elements jo...

1MB Sizes 0 Downloads 34 Views

Engineering Analysis with Boundary Elements 85 (2017) 13–19

Contents lists available at ScienceDirect

Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound

An efficient high order plane wave time domain algorithm for transient electromagnetic scattering analysis G.S. Cheng, Z.H. Fan, D.Z. Ding∗, R.S. Chen Department of Communication Engineering, Nanjing University of Science and Technology, Nanjing 210094, China

a r t i c l e

i n f o

Keywords: Hierarchical vector basis functions Time domain combined field integral equation Plane wave time domain Parallel algorithm

a b s t r a c t An efficient high order plane wave time domain algorithm is presented for analyzing the transient scattering from three dimensional electrically large conducting objects. This method uses a set of hierarchical divergenceconforming vector basis functions to accurately represent the current distribution on the perfect electrically conducting (PEC) surface. The higher order functions can significantly reduce the number of unknowns without compromise on the accuracy. The time domain combined field integral equation (TD-CFIE) is then discretized using the hierarchical divergence-conforming vector basis functions and shifted Lagrange polynomial functions in spatial and time domain, respectively. The final matrix equation can be accelerated using the plane wave time domain (PWTD) algorithm. Finally, a parallel algorithm that can execute on a distributed-memory parallel cluster is developed, which provides an appealing avenue for analyzing the transient scattering from three-dimensional electrically large complex PEC objects. Numerical examples are given to demonstrate the accuracy and efficiency of the method. © 2017 Elsevier Ltd. All rights reserved.

1. Introduction With the continuous improvement of late time instability occurring in the marching-on-in-time (MOT) based time domain integral equation (TDIE) method [1–7], the TDIE method has been of considerable interest in the computational electromagnetic (CEM) community for their appealing avenue in analyzing transient and broadband electromagnetics problems [8–10]. The computational complexity and memory requirement for the analysis of transient scattering from an object using the classical MOT algorithm are 𝑂(𝑁𝑡 𝑁𝑠2 ) and 𝑂(𝑁𝑠2 ), respectively, where Nt and Ns are the numbers of temporal and spatial basis functions. The high computational burden greatly hinders the performance of TDIE algorithm to analyze the large-scale electromagnetic scattering problems. A great deal of attention has been focused on the development of timedomain fast solvers, among which, the plane wave time domain (PWTD) algorithm is the most representative one. Relying on a Whittaker-type expansion of transient fields in terms of propagating plane waves, the computational cost and memory requirements of the PWTD algorithm scale as O(Nt Ns log 2 Ns ) and 𝑂(𝑁𝑠1.5 ), respectively [11–13]. However, the PWTD algorithm can only reduce the computational complexity of the classical MOT algorithm. There is another approach that directly reduces the number of unknowns to be solved by using the higher order spatial basis functions [14–18]. If we combine the fast algorithm with the high order basis function in the framework of the MOT based TDIE,



Corresponding author. E-mail address: [email protected] (D.Z. Ding).

http://dx.doi.org/10.1016/j.enganabound.2017.09.004 Received 19 May 2017; Received in revised form 3 August 2017; Accepted 10 September 2017 0955-7997/© 2017 Elsevier Ltd. All rights reserved.

this will greatly reduce the computational time and memory consumption of the TDIE method, from the perspective of reducing the number of unknowns and the computational complexity. The similar approach has been extensively studied in the frequency domain method of moments (MOM) [21–24]. However, few works has been proposed in the time domain. In this paper, an efficient high order PWTD algorithm is presented for analyzing the transient scattering from three dimensional electrically large conducting objects. We first use a set of hierarchical divergenceconforming vector basis functions [19–20] and shifted Lagrange polynomial functions to discrete the time domain combined field integral equation (TD-CFIE). The hierarchical divergence-conforming vector basis functions can accurately represent the current distribution on the PEC surface and significantly reduce the number of unknowns without compromise on the accuracy. Then the final matrix equation can be accelerated using the PWTD algorithm. Finally, a parallel high order PWTD algorithm is developed, which provides an appealing avenue for analyzing the transient scattering from three-dimensional electrically large complex PEC objects. The described solver executes on a distributedmemory parallel cluster and uses the message passing interface (MPI) paradigm to communicate data between processors. Numerical examples are given to demonstrate the accuracy and efficiency of the proposed method.

G.S. Cheng et al.

Engineering Analysis with Boundary Elements 85 (2017) 13–19

2. Formulation This section describes the implementation of proposed method in detail. Section 2.1 derives and establishes the time-domain electric field, magnetic field and combined field integral equation for solving the transient electromagnetic characteristics from PEC object. Section 2.2 describes the marching on-in-time scheme for solving these equations and the specific spatial basis functions adopted in this paper. Then, the high order PWTD algorithm is introduced in Section 2.3. Finally, a parallel implementation of the scheme is described in Section 2.4. 2.1. Time domain integral equation

Fig. 1. The mapping relationship between the quadratic curvilinear triangle and the parameter triangle.

When an arbitrary shaped PEC object residing in free space with permeability 𝜇 0 and permittivity 𝜀0 is illuminated by a transient incident field Einc (r,t) and Hinc (r,t). The induced current J(r, t) on surface S generates a scattered field Esca (r,t) and Hsca (r,t) which can be characterized as 𝐄𝑠𝑐𝑎 (𝒓, 𝑡) = − +

2.2. Discretization and marching on-in-time scheme To solve the above equations numerically, the PEC surface S is discretized into curved triangular elements firstly. Then, the surface current density J(r, t) can be expanded in terms of a set of vector spatial basis functions and scalar temporal basis functions as

𝜇0 1 𝜕𝐉(𝐫 ′ , 𝜏) 𝑑𝑆 4𝜋 ∬𝑠 𝑅 𝜕𝑡

𝜏 ∇′ ⋅ 𝐉(𝐫 ′ , 𝑡′ ) ′ ∇ 𝑑𝑆 𝑑𝑡 ∫0 4𝜋𝜀0 ∬𝑠 𝑅

(1)

𝐉(𝐫, 𝑡) ≅

𝑁𝑠 𝑁𝑡 ∑ ∑ 𝑛=1 𝑙=1

𝐇𝑠𝑐𝑎 (𝒓, 𝑡) = ∇ ×

𝐉(𝐫 ′ , 𝜏) 1 𝑑𝑆 4𝜋 ∬𝑠 𝑅

𝐙 0 𝐈𝑖 = 𝐕𝑖 −

(4)

[ 𝑖] ⟨ { }⟩| 𝑉 𝑚 = 𝒇 𝑚 (𝐫 ), 𝐕𝑐 𝐄𝑖𝑛𝑐 (𝐫, 𝑡), 𝐇𝑖𝑛𝑐 (𝐫, 𝑡) | |𝑡 = 𝑡 𝑖

𝜏 ∇′ ⋅ 𝐉(𝐫 ′ , 𝑡′ ) ′ ∇ 𝑑𝑆 𝑑𝑡 ∫0 4𝜋𝜀0 ∬𝑠 𝑅

𝐧̂ (𝐫) × 𝐇𝑖𝑛𝑐 (𝒓, 𝑡) 𝐉(𝐫 ′ , 𝜏) 𝐉(𝐫, 𝑡) 1 = − 𝐧̂ (𝐫 ) × 𝑑𝑆 ∇ × = 𝑳ℎ {𝐉(𝐫 , 𝑡)} 2 4𝜋 ∬𝑠 𝑅

(5)

(10) (11)

Based on the parametric triangular patch, higher order hierarchical vector basis functions can then be constructed. The lowest-order vector basis functions are the curvilinear RWG functions and the specific expression is [ ] ) 𝜕𝐫 1 ( 𝜕𝐫 𝒇 𝑒1,0 (𝐫) = + 𝜉2 𝜉1 − 1 𝐽 𝜕 𝜉1 𝜕 𝜉2 [ ] ( ) 1 𝜕𝐫 𝜕𝐫 𝒇 𝑒2,0 (𝐫) = 𝜉1 + 𝜉2 − 1 𝐽 𝜕 𝜉1 𝜕 𝜉2 [ ] 1 𝜕𝐫 𝜕𝐫 𝑒 𝒇 3,0 (𝐫) = 𝜉 + 𝜉2 (13) 𝐽 1 𝜕 𝜉1 𝜕 𝜉2

(6)

Then the time domain combined field integral equation (TD-CFIE) can be derived as a linear combination of the time domain electric and magnetic field integral equations: { } 𝐕𝑐 𝐄𝑖𝑛𝑐 (𝐫, 𝑡), 𝐇𝑖𝑛𝑐 (𝐫, 𝑡) = 𝛼 𝐧̂ (𝐫) × 𝐧̂ (𝐫) × 𝐄𝑖𝑛𝑐 (𝐫, 𝑡) + (1 − 𝛼)𝜂 𝐧̂ (𝐫) × 𝐇𝑖𝑛𝑐 (𝐫, 𝑡) = 𝛼𝑳𝑒 {𝐉(𝐫, 𝑡)} + 𝜂(1 − 𝛼)𝑳ℎ {𝐉(𝐫, 𝑡)} = 𝑳𝑐 {𝐉(𝐫, 𝑡)}

(9)

Here, the spatial basis functions are the higher order hierarchical basis functions defined on curved parametric triangles. We use software ANSYS to generate the quadratic curvilinear elements with six nodes [27]. Fig. 1 shows the mapping relationship between the quadratic curvilinear triangle and the parameter triangle. The quadratic curvilinear triangle can be transformed into the planar parameter triangle using the transformation ( ) ( ) 𝐫 (𝑥, 𝑦, 𝑧) = 𝐫200 𝜉1 2𝜉1 − 1 + 𝐫020 𝜉2 2𝜉2 − 1 ( ) + 𝐫002 𝜉3 2𝜉3 − 1 + 𝐫110 4𝜉1 𝜉2 + 𝐫011 4𝜉2 𝜉3 + 𝐫101 4𝜉1 𝜉3 (12)

(𝒓, 𝑡) 𝜇 1 𝜕𝐉(𝐫 ′ , 𝜏) = 𝐧̂ (𝐫) × 𝐧̂ (𝐫) × 0 𝑑𝑆 ∬ 4𝜋 𝑠 𝑅 𝜕𝑡

= 𝑳𝑒 {𝐉(𝐫, 𝑡)}

𝐙 𝑖 −𝑗 𝐈𝑗

The expression of matrix elements is [ 𝑖−𝑗 ] ⟨ { }⟩ = 𝒇 𝑚 (𝐫 ), 𝑳𝑐 𝒇 𝑛 (𝐫 )𝑇𝑖−𝑗 (𝑡) 𝑍 𝑚𝑛

𝑖𝑛𝑐

− 𝐧̂ (𝐫) × 𝐧̂ (𝐫) ×

𝑖−1 ∑ 𝑗=1

we can get the time domain electric and magnetic field integral equations: 𝐧̂ (𝐫) × 𝐧̂ (𝐫) × 𝐄

(8)

Here 𝐼𝑛𝑙 are the unknown coefficients. Nt and Ns denote the number of time steps and spatial basis functions, respectively. Substituting (8) into (7), Galerkin testing the Eq. (7) with the surface basis function fm (r) and point matching at time iΔt, leads to the system of equations that can be represented explicitly in a matrix form as

(2)

where J(r, t) denotes the induced current on the PEC surface S, r and r′ represent the observation point and source point, respectively, R = |r − r′|, c is the speed of light in free space, 𝜏 = t − R/c is the retarded time. According to the boundary condition on the PEC surface: ]| [ (3) 𝐧̂ (𝐫) × 𝐧̂ (𝐫) × 𝐄𝑖𝑛𝑐 (𝐫, 𝑡) + 𝐄𝑠𝑐𝑎 (𝐫, 𝑡) | = 0 |𝑆 ]| [ 𝐧̂ (𝐫) × 𝐇𝑖𝑛𝑐 (𝐫, 𝑡) + 𝐇𝑠𝑐𝑎 (𝐫, 𝑡) | = 𝐉(𝐫, 𝑡) |𝑆

𝐼𝑛𝑙 𝒇 𝑛 (𝐫)𝑇𝑙 (𝑡)

where J is the Jacobian. The first subscript denotes the edge and the second subscript denotes the order of the polynomials. The superscript denotes that they are edge-based basis functions. The higher order bases functions are obtained by forming the product of the curvilinear RWG basis functions with a set of polynomial functions. Here we give the edge-based and face-based basis functions up to

(7)

where 𝐧̂ (𝐫) is the outward directed unit normal, 𝜂 is the wave impedance of the free space. 𝛼is a positive real constant. The EFIE and the MFIE can be derived from the CFIE by setting 𝛼 = 1 and 𝛼 = 0, respectively. So the discussions hereafter will focus on the CFIE. 14

G.S. Cheng et al.

Engineering Analysis with Boundary Elements 85 (2017) 13–19

the order of 2.5 associated with edge 1. The remaining two subsets of edge-based associated with edge 2 and edge 3 and one subset of facebased functions associated with edge 2 can be obtained by rotating the simplex coordinates. The edge-based functions associated with edge 1 can be expressed as √ ( ) 3 2𝜉2 − 1 + 𝜉1 𝒇 𝑒1,0 √ ] )2 5[ ( 𝒇 𝑒1,2 (𝐫) = 3 𝜉2 − 𝜉3 − 1 𝒇 𝑒1,0 2 𝒇 𝑒1,1 (𝐫) =

(14)

The face-based functions associated with edge 1 can be expressed as √ 𝒇 𝑓1,01 (𝐫) = 2 3𝜉1 𝒇 𝑒1,0 √ ( ) 𝒇 𝑓1,02 (𝐫) = 2 3𝜉1 5𝜉1 − 3 𝒇 𝑒1,0 √ ( ) 𝒇 𝑓1,11 (𝐫) = 6 5 𝜉2 − 𝜉3 𝜉1 𝒇 𝑒1,0

(15)

Here, the first subscript denotes the edge and the sum of the second subscript denotes the order of the polynomials. The superscript denotes that they are face-based basis functions. In our implementation, the nonsingular part of the inner products in (10)–(11) is evaluated using Gaussian quadrature rules and the (near-) singular part is evaluated on the tangent plane triangle using the exact integration technique proposed in [5].

2.3. Higher order PWTD algorithm Although the higher order method can significantly reduce the number of unknowns, its computational complexity and memory requirement still scales as 𝑂(𝑁𝑡 𝑁𝑠2 ) and 𝑂(𝑁𝑠2 ), respectively. The classical HOTDIE method will quickly swamp unavailable due to its high computational resources for the analysis of large-scale scattering problem. Here, the PWTD algorithm is employed to reduce the memory requirement and accelerate the iterative solution of the HO-TDIE method. In order to embed the PWTD in a hierarchical framework, the scatterer needs to be enclosed in a virtual cube box and recursively subdivided into eight small boxes until the edge length of the lowest cube is about half wavelength. We use Nl to denote the number of levels, level 1 refers to the lowest level and level Nl is the coarsest level. For levels i = 1, 2, ..., Nl , Rs (i) denotes the size of the level ibox. Then two groups at level i are said to be far-field pairs if the distance between their centers is greater than a prescribed distance Rmin (i) = kRs (i) and if their parents are near-field pairs. As in the multilevel scheme, all nearfield groups are at the lowest level and evaluated using the classical MOT algorithms. The far field interaction can be reconstructed using the PWTD algorithms via the aggregation-translation-disaggregation procedure. The near- and far-field groups in a hierarchical framework are illustrated in Fig. 2. Specifically, consider the observation group Gi and source group Gj are the far-field pair. The observation and source groups contain the observation function fm (r) and source function fn (r′), respectively. ri and rj denote the centers of the two groups, respectively. rm and rn are the centers of the edges m and n for the edge-based functions and the centroids of the triangles Tm and Tn which contain the mth and nth edges for the faced-based functions, respectively. The source current density Jn (r′,t) can be represented as 𝐉𝑛 (𝐫 ′ , 𝑡) =

𝑁𝑡 ∑ 𝑗=1

𝐽𝑛,𝑗 𝐟𝑛 (𝐫 ′ )𝑇𝑗 (𝑡) = 𝐟𝑛 (𝐫 ′ )𝑔𝑛 (𝑡)

Fig. 2. The schematic diagram of the near- and far-field groups in a hierarchical framework.

𝐉𝑛 (𝐫 ′ , 𝑡) = =

𝐿 ∑ 𝑙=1 𝐿 ∑ 𝑙=1

𝐉𝑛,𝑙 (𝐫 ′ , 𝑡) = 𝐟𝑛 (𝐫 ′ )

𝐿 ∑ 𝑙=1

𝐟𝑛 (𝐫 ′ )𝑔𝑛,𝑙 (𝑡)

𝑙 𝑀𝑡 ∑ 𝑗=(𝑙−1)𝑀𝑡 +1

𝐽𝑛,𝑗 𝑇𝑗 (𝑡)

(17)

The tested scattering field generated by the subsignal Jn,l (r′, t) with the observation function fm (r) can be expressed as ⟨

{ }⟩ 𝒇 𝑚 (𝐫 ), 𝑳𝑐 𝐉𝑛,𝑙 (𝐫 ′ , 𝑡) =

𝜂0 8𝜋 2 𝑐 2

𝑖

𝑡

∫𝑇𝑙,𝑠𝑡𝑜𝑝

𝑑 𝑡′

𝑖

𝐾 𝐾 ∑ ∑ 𝑝=0 𝑞=−𝐾 𝑖

[ 𝜔𝑖𝑝𝑞 𝛼𝑆𝑚− (̂𝐤𝑖𝑝𝑞 , 𝑡′ , ̂𝐤𝑖𝑝𝑞 )

]𝑇 +(1 − 𝛼)𝑆𝑚− (̂𝐤𝑖𝑝𝑞 , 𝑡′ , 𝐧̂ ) ∗ 𝑇 (̂𝐤𝑖𝑝𝑞 , 𝑡′ ) 𝑖 ∗ [𝑆𝑛+ (̂𝐤𝑖𝑝𝑞 , 𝑡′ , ̂𝐤𝑖𝑝𝑞 )] ∗ 𝑔𝑛,𝑙 (𝑡′ )

(18)

where Ki represents the number of spherical harmonics, 𝜔𝑖𝑝𝑞 are quadrature weights on the plane wave directions ̂𝐤𝑖 , p = 0, ..., Ki and q = −Ki ,..., 𝑝𝑞

Ki represent directions of outgoing rays and incoming rays,∗ represents time convolution. The aggregation and disaggregation functions 𝑆𝑜± (̂𝐤𝑖𝑝𝑞 , 𝑡′ , 𝒗̂ ) are

(16)

𝑆𝑜± (̂𝐤𝑖𝑝𝑞 , 𝑡′ , 𝒗̂ ) =

The time signal gn (t) should be broken up into L consecutive subsignals gn,l (t). The duration of each subsignal is Ts = Mt Δt(LMt = Nt ) and the source current density can be represented as

∫𝑆𝑜

𝑑𝑠′ 𝛎̂

×𝒇 𝑜 (𝐫 ′ )𝛿(𝑡 ± ̂𝐤𝑖𝑝𝑞 (𝐫 ′ − 𝐫𝑜𝑐 )∕𝑐 ), 𝑜 ∈ {𝑚, 𝑛} 15

(19)

G.S. Cheng et al.

Engineering Analysis with Boundary Elements 85 (2017) 13–19

The translation function 𝑇 (̂𝐤𝑖𝑝𝑞 , 𝑡′ ) is

Table 1 The parameters of incident wave in different examples.

𝑖

𝑇̄ (̂𝐤𝑖𝑝𝑞 , 𝑡′ ) = −

𝐾 𝑐𝜕𝑡3 ∑

2𝐑𝑐

𝑘=0

(2𝑘 + 1)𝑃𝑘 (𝑐𝑡∕𝑅𝑐 )

× 𝑃𝑘 (cos 𝜃), |𝑡| ≤ 𝑅𝑐 ∕𝑐

(20)

where Pk (•) is the Legendre polynomial,𝜃 is the angle between vectors ̂𝐤𝑖 and Rc . 𝑝𝑞 The PWTD algorithm can efficiently accelerate the computation of the matrix vector product on the right-hand side of (9) as follows. First, for source group Gj , outgoing rays for all directions are constructed by convolving the aggregation function with the subsignal Jn,l (r′, t). Then the outgoing rays for source group Gj are translated into incoming rays for observation group Gi by convolving with the translation function 𝑇 (̂𝐤𝑖𝑝𝑞 , 𝑡′ ).Finally, the incoming rays are projected onto the observation function fm (r)by convolving with the disaggregation functions. Detailed description of the PWTD algorithm will be omitted, and the readers can refer to the literature [11–13]. There are both edge-based functions and face-based functions in the HO-TDIE method. It should be noted that in the aggregation and disaggregation procedure of the PWTD method, the corresponding aggregation factor and disaggregation factor for the face-based functions should use the center of the triangles where the face-based function is located.

Example

𝑖𝑛𝑐 𝒌̂

𝒑̂ 𝑖𝑛𝑐

f0

fbw

1 2 3 4

𝒛̂ 𝒙̂ 𝒛̂ 𝒙̂

𝒙̂ −𝒛̂ 𝒙̂ −𝒛̂

150 MHz 200 MHz 150 MHz 300 MHz

300 MHz 200 MHz 300 MHz 400 MHz

3. Numerical results This section presents several numerical examples to demonstrate the accuracy and efficiency of the parallel HO-PWTD solver. The scatterers considered here are all illuminated by the modulated Gaussian plane wave defined by [ ( )] 𝑖𝑛𝑐 𝐄𝑖𝑛𝑐 (𝒓, 𝑡) = 𝒑̂ 𝑖𝑛𝑐 cos 2𝜋𝑓0 𝑡 − 𝒓 ⋅ 𝑘̂ ∕𝑐 [ ] ( )2 ̂ 𝑖𝑛𝑐 ∕𝑐 ∕𝜎 2 (21) × exp −0.5 𝑡 − 𝑡𝑝 − 𝒓 ⋅ 𝒌 where f0 and fbw are the center frequency and bandwidth of the plane 𝑖𝑛𝑐 wave, respectively. 𝒌̂ and 𝒑̂ 𝑖𝑛𝑐 denote the incident and polarization direction, respectively, 𝜎 = 6∕𝜋𝑓𝑏𝑤 . It should be noted the power in the incident pulse is down by 160 dB at 𝑓 = 𝑓0 ± 𝑓𝑏𝑤 ∕2 relative to f0 . All computations were carried out on a 64-bit workstation Intel(R) Xeon(R) CPU E7-8850 @2.3 GHz, 2 TB RAM. Double precision computation is always used. First, the transient scattering from a sphere with radius 1.0 m is analyzed to validate the accuracy of the HO-TDIE solver. The parameters of incident field have been listed in Table 1. The time step size is Δt = 0.333ns and Nt = 300. In Table 2, various parameters affect the solution precision and computational efficiency for different order basis functions presented. Here the relative RMS error is used as a measuring tool for the comparison of the accuracy and the relative RMS error is defined as √ ∑𝑁 𝑎 ∑𝑁𝑎 𝑖 𝑖 𝑖 |2 ∕ 𝑖 |2 , 𝑅𝑀𝑆 = |𝑓𝐻 − 𝑓𝑀𝑖𝑒 |𝑓𝑀𝑖𝑒 𝑓𝐻 𝑂−𝑇 𝐷𝐼 𝐸 𝑂−𝑇 𝐷𝐼 𝐸

2.4. Parallel implementation In order to parallelize the HO-PWTD solver described in Section 2.3 efficiently, the “space partitioning” scheme is used to distribute the task load evenly and minimize the communications between the processors. All of the nonempty boxes at the level 1 are divided into P processors equally. After distributing the boxes between each processor, the near-field interactions are directly calculated using the classic MOT process and stored in memory without any communication; this operation can be parallelized in a straightforward manner. The far-field contributions are calculated through the PWTD algorithm via the aggregation-translation-disaggregation procedure.

𝑖=1

𝑖 𝑓𝑀𝑖𝑒

2.4.1. Aggregation The far-field interaction begins with aggregating basis functions at the lowest level to obtain the radiation pattern. Each processor aggregates the basis functions at the lowest level to its group center to obtain the outgoing rays. Then these outgoing rays information will be written into the random access memory (RAM). When it is necessary to construct the outgoing rays of the second level, the outgoing rays of the lowest level required is read out from the RAM, and then the interpolation and splicing procedure is carried out to obtain the rays for the second level. This procedure repeats until the coarsest levels.

𝑖=1

and are the radar cross section (RCS) values obtained by the HOTDIE method and Mie solution in the elevation angles between 𝜃 = 0o and 180o ,with the azimuthal angle set as 𝜙 = 0o , respectively. Na is the number of sampled points normally chosen as Na = 720. In Table 2, first the hierarchical basis functions of order 1.5 is employed and the mesh size is fixed as 0.3𝜆, gradually increasing the number of Gaussian integration points, we can find that the solution time increases while the RRMS error of the RCS is almost constant. This is because six-point Gaussian quadrature has reached the saturation accuracy for the mesh size. When the mesh size increases to 0.4𝜆, the computational time and memory consumption of the HO-TDIE method are further reduced, but the computational precision should be questioned. Then the hierarchical basis functions of order 2.5 is employed, the relationship between the mesh size, the number of Gaussian integration points, the memory/time consumption, and the accuracy of the proposed method is verified. It can be seen that when the 2.5-order basis function is used, the system matrix becomes ill-conditioned and the number of iterations and the iterative solution time increases obviously. Therefore, in the absence of special instructions, the hierarchical basis functions of order 1.5 with a mesh size of 0.3𝜆∼ 0.4𝜆 are employed in the following examples. Then, the efficiency of the HO-PWTD solver is demonstrated through analyzing the transient scattering from a simplified model of B2 plane as shown in Fig. 3(a). The plane fits within a cuboid of dimension 10.492m × 26.188m × 1.546m. The parameters of incident field have been listed in Table 1. We use the RWG-PWTD scheme and the HOPWTD scheme, respectively. For the RWG-PWTD scheme, the plane is discretized into 56884 flat triangular patches with a mesh size of

2.4.2. Translation Since all the outgoing rays’ information is stored in the RAM, the translation procedure can be implemented conveniently. For example, if cube i at processor a needs the outgoing rays of the cube j at processor b, the outgoing rays can be directly read out from the RAM and no communication is required. 2.4.3. Disaggregation The disaggregation stage is generally the inverse of the aggregation stage. The far-field contributions at the observers are reconstructed using the information contained in the incoming rays. Starting from the coarsest layer, the incoming rays are projected onto their children by truncation and anterpolation. This procedure repeats until the lowest levels. Consider all the steps of the multilevel algorithm, each processor is responsible for constructing the outgoing rays and incoming rays of the boxes belonging to it at each level. 16

G.S. Cheng et al.

Engineering Analysis with Boundary Elements 85 (2017) 13–19

Table 2 Computation time and memory requirements for the calculations of bi-RCS of the PEC sphere using the HO-TDIE solver. Bases order

Avg. patch size(𝜆)

Gauss integral point number

Total unknowns

0.5 1.5 1.5 1.5 1.5 1.5 1.5 2.5 2.5 2.5 2.5

0.1 0.3 0.3 0.3 0.4 0.4 0.5 0.3 0.4 0.5 0.6

6 6 12 25 6 25 25 25 25 25 25

3852 1400 1400 1400 730 730 570 2940 1533 1197 735

RRMS error (%) 30 MHz

150 MHz

270 MHz

0.013 0.105 0.175 0.111 0.145 0.161 0.214 0.054 0.087 0.285 0.464

0.923 0.617 0.697 0.711 1.047 0.762 1.216 0.837 1.015 1.320 2.078

2.131 1.490 1.873 1.621 2.361 2.333 4.526 1.490 1.598 2.908 3.740

Memory (MB)

Matrix filling time(s)

Iteration time(s)

Iteration step

757.1 145.9 151 156.9 52.7 55.1 39.9 634.3 211.7 137.8 65.8

227 10 24 72 4 20 13 213 61 39 15

716 292 318 323 46 63 37 2512 826 492 192

18 68 71 71 66 69 70 122 116 122 135

Total time(s) 968 303 344 399 50 85 52 2733 889 532 208

Table 3 Comparison of the time and memory requirements.

HO-PWTD RWG-PWTD

Unknowns

Total time (h)

Memory (GB)

33440 85326

29.8 51.9

15.4 39.3

Fig. 4. Parallelization efficiency for the solution of the transient scattering from a sphere of radius 4 m.

Next, in order to verity the parallel efficiency of the parallel HOPWTD solver, the parallel HO-PWTD solver is applied to the analysis of transient scattering from a sphere with radius 4.0 m. The parameters of incident field have been listed in Table 1. The sphere is discretized with 4614 curvilinear triangular patches for order 1.5 hierarchical basis functions with a mesh size of 0.3𝜆, giving rise to 23070 HO unknowns. The lowest group size is 0.4 m, the time step size is Δt = 0.333ns and Nt = 300. Fig. 4 depicts the efficiency with respect to the solution with a single processor when the solution is parallelized into 2, 4, 8, 16 and 32 processes. The parallel efficiency is defined as 𝜀p = T1 /(pTp ) where Tp is the solution time with p processes. Fig. 4 shows that the overall efficiency is about 80% when the number of processes is 32. Finally, in order to verify the capacity of the parallel HO-PWTD solver, the parallel HO-PWTD solver is applied to the analysis of transient scattering from the simplified model of B2 plane, which fits in a hypothetical box with dimensions 20.984m × 52.376m × 3.092m. The parameters of incident field have been listed in Table 1. For the HObased PWTD scheme, the plane is discretized with 62828 curvilinear triangular patches for order 1.5 hierarchical basis functions with a mesh size of 0.3𝜆, giving rise to 314140 HO unknowns. The time step size is Δt = 0.2ns, and Nt = 800. The solver requires about 387.5 GB of memory and 12.5 h of CPU time when the solution is parallelized into 80 processes. If we use the RWG-PWTD scheme, the plane is discretized

Fig. 3. (a) Simplified model of the B2 plane. (b) Comparison of the RCS data for the model of B2 plane obtained using HO-PWTD and RWG-PWTD method at different frequencies.

0.1𝜆and the total number of unknowns is 85326. For the HO-based PWTD scheme, the plane is discretized with 6688 curvilinear triangular patches for order 1.5 hierarchical basis functions with a mesh size of 0.3𝜆, giving rise to 33440 HO unknowns. The lowest group size is 0.4 m, the time step size is Δt = 0.333ns, and Nt = 500. RCS data at 114, 204 and 300 MHz obtained using the HO-PWTD method is compared against RWG-PWTD data in Fig. 3(b); again, all RCS results agree well with each other. In order to verify the efficiency of the HO-PWTD algorithm, the CPU and memory requirements for the HO-PWTD scheme and RWG-PWTD scheme are listed in Table 3. It can be seen that the significant reduction in computational time and memory requirements can be achieved through the HO-PWTD algorithm. 17

G.S. Cheng et al.

Engineering Analysis with Boundary Elements 85 (2017) 13–19

(MLFMA) [25–26] at 125 MHz, 300 MHz and 425 MHz, respectively in Fig. 5(a)–(c); again, all results are in good agreement. 4. Conclusion In this paper, an efficient high order PWTD algorithm is presented for analyzing the transient scattering from three dimensional electrically large conducting objects. The hierarchical divergence-conforming vector basis functions can accurately represent the current distribution on the PEC surface and significantly reduce the number of unknowns without compromise on the accuracy. Then the final matrix equation can be accelerated using the PWTD algorithm. Finally, a parallel high order PWTD algorithm using the MPI paradigm is developed. Numerical examples are given to demonstrate the accuracy and efficiency of the proposed method. Acknowledgment This work was supported in part by Natural Science Foundation of 61522108, 61431006, 61371037, and in part by the Fundamental Research Funds for the Central Universities of no. 30916011321. References [1] Rynne BP, Smith PD. Stability of time marching algorithms for the electric field integral equations. J Electromagn Waves Appl 1990;12:1181–205. [2] Sadigh A, Arvas E. Treating the instabilities in marching-on-in-time method from a different perspective. IEEE Trans Antennas Propag 1993;41(12):1695–702. [3] Sarkar TK, Lee W, Rao SM. Analysis of transient scattering from composite arbitrarily shaped complex structures. IEEE Trans Antennas Propag 2000;48:1625–34. [4] Weile DS, Pisharody G, Chen N-W, Shanker B, Michielssen E. A novel scheme for the solution of the time-domain integral equations of electromagnetics. IEEE Trans Antennas Propag 2004;52(1):283–95. [5] Shanker B, Lu M, Yuan J, Michielssen E. Time domain integral equation analysis of scattering from composite bodies via exact evaluation of radiation fields. IEEE Trans Antennas Propag 2009;57(5):1506–20. [6] Shi Y, Xia M-Y, Chen R-S, Michielssen E, Lu M. Stable electric field TDIE solvers via quasi-exact evaluation of mot matrix elements. IEEE Trans Antennas Propaga 2011;59:574–85. [7] Beghein Y, Cools K, Bagci H, De Zutter D. A space-time mixed Galerkin marching-on-in-time scheme for the time-domain combined field integral equation. IEEE Trans Antennas Propag 2013;61(3):1228–38. [8] Rao SM, Wilton DR. Transient scattering by conducting surfaces of arbitrary shape. IEEE Trans Antennas Propag 1991;39(1):56–61. [9] Chew W, Jin J-M, Michielssen E, Song J. Fast and efficient algorithms in computational electromagnetics. Boston, MA, USA: Artech House; 2001. [10] Shanker B, Ergin AA, Aygun K, Michielssen E. Analysis of transient electromagnetic scattering from closed surfaces using a combined field integral equation. IEEE Trans Antennas Propaga 2000;48:1064–74. [11] Ergin AA, Shanker B, Michielssen E. The plane wave time domain algorithm for the fast analysis of transient wave phenomena. IEEE Antennas Propaga Mag 1999;41:39–52. [12] Shanker B, Ergin AA, Aygun K, Michielssen E. Analysis of transient electromagnetic scattering phenomena using a two-level plane wave time domain algorithm. IEEE Trans Antennas Propaga 2000;48:510–23. [13] Shanker B, Ergin AA, Lu M, Michielssen E. Fast analysis of transient electromagnetic scattering phenomena using the multilevel plane wave time domain algorithm. IEEE Trans Antennas Propag 2003;51:628–41. [14] Graglia RD, Wilton DR, Peterson AF. Higher order interpolatory vector bases for computational electromagnetics. IEEE Trans Antennas Propaga 1991;45:329–42. [15] Weile DS, Ergin AA, Shanker B, Michielssen E. An accurate discretization scheme for the numerical solution of time domain integral equations. In: IEEE Antennas Propaga Int. Symp., Salt Lake City, UT; 2000. [16] Wildman RA, Pisharody G, Weile DS, Shanker B, Michielssen E. An accurate scheme for the solution of the time-domain Integral equations of electromagnetics using higher order vector bases and bandlimited extrapolation. IEEE Trans Antennas Propaga 2004;52(11):2973–84. [17] Andriulli FP, Bagci H, Michielssen E. Time-domain single-source integral equations for analyzing scattering from homogeneous penetrable objects. IEEE Trans Antennas Propaga 2013;61(3):1239–54. [18] Valdés F, Ghaffari-Miab M, Andriulli FP, Cools K, Michielssen E. High-order Calderón preconditioned time domain integral equation solvers. IEEE Trans Antennas Propag. 2013;61(5):2570–88. [19] Zha LP, Hu YQ, Su T. Efficient surface integral equation using hierarchical vector bases for complex EM scattering problems. IEEE Trans Antennas Propag 2012;60:952–7. [20] Graglia RD, Peterson AF, Andriulli FP. Curl-conforming hierarchical vector bases for triangles and tetrahedral. IEEE Trans Antennas Propag 2011;59(3):950–9.

Fig. 5. Comparison of the RCS data for the model of B2 plane obtained using HO-PWTD and MLFMA method at (a) 125 MHz, (b) 300 MHz, (c) 425 MHz.

into 640908 flat triangular patches with a mesh size of 0.1𝜆 and the total number of unknowns is 961362. The RWG-PWTD solver requires about 988.2 GB of memory and 21.8 h of CPU time when the solution is parallelized into 80 processes. The bistatic RCS of the plane is compared with the frequency-domain multilevel fast multipole algorithm 18

G.S. Cheng et al.

Engineering Analysis with Boundary Elements 85 (2017) 13–19

[21] Donepudi KC, Song J, Jin J-M, Kang G, Chew WC. A novel implementation of multilevel fast multipole algorithm for higher order Galerkin’s method. IEEE Trans Antennas Propag 2000;48(8):1192–7. [22] Ismatullah, Eibert TF. Surface integral equation solutions by hierarchical vector basis functions and spherical harmonics based multilevel fast multipole method. IEEE Trans Antennas Propag 2009;57(7):2084–93. [23] Nie ZP, Ma W, Ren Y, Zhao Y, Hu J, Zhao Z. A wideband electromagnetic scattering analysis using MLFMA with higher order hierarchical vector basis functions. IEEE Trans Antennas Propag 2009;57(10):3169–78.

[24] Borries O, Meincke P, Jorgensen E, Hansen C. Multilevel fast multipole method for higher order discretizations. IEEE Trans Antenna Propag 2014;62(9):4695–705. [25] Ergul O, Gurel L. The multilevel fast multipole algorithm (MLFMA) for solving large-scale computational electromagnetics problems. IEEE Press, Wiley; 2014. [26] Ergul O, Gurel L. Efficient parallelization of the multilevel fast multipole algorithm for the solution of large-scale scattering problems. IEEE Trans Antenna Propag 2008;56(8):2335–45. [27] ANSYS® Academic Research, Release 10.0.

19