# Simulation of Airborne Antenna Array Layout Problems Using Parallel Higher-Order MoM

Authors: Zhongchao Lin; Yu Zhang; Shugang Jiang; Xunwang Zhao; Jingyan Mo
Journal: International Journal of Antennas and Propagation
 (2014)

---

## Abstract

The parallel higher-order Method of Moments based on message passing interface (MPI) has been successfully used to analyze the changes in radiation patterns of a microstrip patch array antenna mounted on different positions of an airplane. The block-partitioned scheme for the large dense MoM matrix and a block-cyclic matrix distribution scheme are designed to achieve excellent load balance and high parallel efficiency. Numerical results demonstrate that the rigorous parallel Method of Moments can efficiently and accurately solve large complex electromagnetic problems with composite structures.


---

## Body

## 1. Introduction

As the most traditional and widely adopted method, the method of moments (MoM) is a numerically accurate method for solving electromagnetic problems [1]. However, its memory requirement and computational complexity grow rapidly with O (N
2) or even O (N
3), where N is the number of unknowns. This makes the method very difficult to deal with large complex objects, such as airborne platform. To overcome this difficulty, on one hand, higher-order polynomials over wires and quadrilateral plates are used as basis functions over larger subdomain patches to reduce the number of unknowns [2]. Polynomial expansions for the basis functions over larger subdomains lead to a good approximation of the current distributions over large surfaces using approximately ten unknowns per wavelength squared, which is much less than the use of piecewise Rao-Wilton-Glisson basis functions (RWGs) [2, 3]. On the other hand, with the rapid development in computer hardware capabilities, massively parallel computing on computer clusters and multicore processors has been the method of choice for solving modern engineering and science problems arising from extremely complicated real-life applications. In addition, to improve the computation efficiency of the higher-order MoM, the large dense MoM matrix is divided into a number of smaller block matrices that are nearly equal in size and distributed among all participating processes. The strategy of distributing the blocks is designed appropriately according to the parallel lower/upper (LU) decomposition solver to minimize the communication between processes.The parallel efficiency and scalability of the parallel higher-order MoM were analyzed in our previous studies [2]. The accuracy of algorithm was validated through comparison with measurement and MoM code with RWG basis function which can be found in [2, 4]. In this paper, the parallel higher-order MoM is used to analyze the changes in radiation patterns of an antenna array installed on different positions of a lager airplane, which have been simulated with a MoM hybrid with high frequency asymptotic methods [5, 6] rather than a rigorous MoM.The paper is organized as follows. In Section2, the basics of integral equations, higher-order basis functions, and the parallel strategies are listed, respectively. Section 3 lists some numerical results: Section 3.1 validates the accuracy of this paper’s method through the comparison with the finite element method (FEM) [7]; Section 3.2 describes the computation platforms; Section 3.3 presents radiation patterns of a microstrip patch array antenna with 37 × 9 Elements; Section 3.4 analyzes the changes in radiation patterns of the microstrip array installed on different positions of a lager airplane; Section 3.5 analyzes another case of airborne array that a microstrip array sticks on a lager airplane. Finally, the conclusions are presented in Section 4 followed by the acknowledgements.


## 2. Parallel Higher-Order Method of Moments

### 2.1. Formulation of the Integral Equations

The method is based on the solution of surface integral equations (SIEs) [8–10] in the frequency domain for equivalent currents over dielectric boundary surfaces and electric currents over perfect electric conductors (PECs). The set of integral equations obtained are solved by using MoM, specifically using the Galerkin method. The method is able to handle inhomogeneous dielectrics categorized by a combination of various homogeneous dielectrics. Therefore, any composite metallic and dielectric structure can be represented as an electromagnetic system consisting of a finite number of finite-size linear, homogeneous, and isotropic regions situated in an unbounded linear, homogeneous, isotropic environment.Let us consider an arbitrary regioni with a nonzero electromagnetic field. According to the surface equivalent theorem, the field outside region i becomes zero. Hence, the region outside region i can be homogenized with respect the region i, and then, a multiple-region problem may be decomposed into n single-region problems. The total electric field in region i may be expressed as a sum of scattered and incident fields as follows:
(1)
E
⃑
(
i
)
=
∑
j
=
0
,
j
≠
i
n
E
⃑
(
i
)
(
J
⃑
s
i
j
,
M
⃑
s
i
j
)
+
E
⃑
inc
(
i
)
,
where E
⃑
(
i
)
(
J
⃑
s
i
j
,
M
⃑
s
i
j
) represents the scattered field inside region i, which is calculated by the currents placed on the boundary surface between regions i and j, and E
⃑
inc
(
i
) is the corresponding incident field. The electric and magnetic scattered fields inside region i are given by
(2)
E
⃑
(
i
)
(
J
⃑
s
i
j
,
M
⃑
s
i
j
)
=
-
L
(
i
)
(
J
⃑
s
i
j
)
+
K
(
i
)
(
M
⃑
s
i
j
)
,
(3)
H
⃑
(
i
)
(
J
⃑
s
i
j
,
M
⃑
s
i
j
)
=
-
K
(
i
)
(
J
⃑
s
i
j
)
-
ε
(
i
)
μ
(
i
)
L
(
i
)
(
M
⃑
s
i
j
)
,
where L
(
i
) and K
(
i
) are linear operators given by
(4)
L
(
i
)
(
X
⃑
s
i
j
)
=
j
ω
μ
(
i
)
∫
S
i
k
(
X
⃑
s
i
j
(
r
⃑
i
j
)
g
(
i
)
(
r
⃑
,
r
⃑
′
)
+
1
ω
2
ε
(
i
)
μ
(
i
)
∇
s
i
j
m
m
m
m
m
m
m
m
m
·
X
⃑
s
i
j
(
r
⃑
i
j
)
∇
g
(
i
)
(
r
⃑
,
r
⃑
′
)
1
ω
2
ε
(
i
)
μ
(
i
)
)
d
S
i
j
,
K
(
i
)
(
X
⃑
s
i
j
)
=
∫
S
i
j
X
⃑
s
i
j
(
r
⃑
i
j
)
×
∇
g
(
i
)
(
r
⃑
,
r
⃑
′
)
d
S
i
j
,
where X
⃑
s
i
j  denotes the electric or magnetic current, r
⃑
′ is the vector position of the source point, r
⃑ is the vector position of the field point, and g
i
(
r
⃑
,
r
⃑
′
) is Green’s function for the homogeneous medium i. Making use of the boundary condition equations given by
(5)
n
^
i
j
×
(
E
⃑
(
i
)
-
E
⃑
(
j
)
)
=
0
,
(6)
n
^
i
j
×
(
H
⃑
(
i
)
-
H
⃑
(
j
)
)
=
0
and replacing (1), (2), and (3) into (4) and (6), the integral equations for regions i and k are obtained in the form
(7)
n
^
i
k
×
{
∑
j
=
0
,
j
≠
i
n
[
L
(
i
)
(
J
⃑
s
i
j
)
-
K
(
i
)
(
M
⃑
s
i
j
)
]
mmmm
-
∑
j
=
0
,
j
≠
k
n
[
L
(
j
)
(
J
⃑
s
k
j
)
-
K
(
j
)
(
M
⃑
s
k
j
)
]
}
=
n
^
i
k
×
(
E
⃑
inc
(
i
)
-
E
⃑
inc
(
k
)
)
,
n
^
i
k
×
{
∑
j
=
0
,
j
≠
i
n
(
ε
(
i
)
μ
(
i
)
L
(
i
)
(
M
⃑
s
i
j
)
+
K
(
i
)
(
J
⃑
s
i
j
)
)
m
m
m
m
-
∑
j
=
0
,
j
≠
k
n
(
ε
(
j
)
μ
(
j
)
L
(
j
)
(
M
⃑
s
k
j
)
+
K
(
j
)
(
J
⃑
s
k
j
)
)
}
=
n
^
i
j
×
(
H
⃑
inc
(
i
)
-
H
⃑
inc
(
k
)
)
.These two sets of equations represent a general form of the Poggio-Miller-Chang-Harrington-Wu (PMCHW) [2, 11] formulation. When the boundary surface of two different regions is PEC, the magnetic currents are equal to zero at the boundary surface and the equations degenerates into the electric field integral equation (EFIE) [12, 13].


### 2.2. Higher-Order Basis Functions

The electric and magnetic currents are approximated by higher-order polynomials, which reduce the number of unknowns compared with the rational piecewise basis functions. The code makes use of truncated cones for wires and bilinear patches to characterize other surfaces.For bilinear surfaces, the surface current is decomposed into itsp- and s-components, p and s being the two parametric coordinates of the unit quadrangle, p
,
s
∈
[
-
1
,
1
]. The approximation for the s-component of the electric current is (analogous expressions stand for the p-component and for the magnetic current)
(8)
J
⃑
s
(
p
,
s
)
=
∑
i
=
0
N
p
[
c
i
1
E
⃑
i
(
p
,
s
)
+
c
i
2
E
⃑
i
(
p
,
-
s
)
+
∑
j
=
2
N
s
a
i
j
P
⃑
i
j
(
p
,
s
)
]
,
where N
p and N
s are the degrees of the approximations along the coordinates, and a
i
j, c
i
1, and c
i
2 are the unknown coefficients.Thus, expression (8) stands for the representation of the current in terms of edge basis functions E
⃑
i
(
p
,
s
) and interior or patch basis functions P
⃑
i
j
(
p
,
s
) which can be compactly expressed as
(9)
E
⃑
i
k
(
p
,
s
)
=
α
s
|
α
p
×
α
s
|
{
p
i
N
(
s
)
,
k
=
1
p
i
N
(
-
s
)
,
k
=
2
,
P
⃑
i
j
(
p
,
s
)
=
α
s
|
α
p
×
α
s
|
p
i
S
j
(
s
)
,
N
(
s
)
=
1
-
s
2
,
S
i
(
s
)
=
{
s
i
-
1
i
is
even
s
i
-
s
i
is
odd
,
where symbols a
s and a
p denote the unitary vectors along the transformed s and p coordinates.Edge basis functionsE
⃑
i
1 and patch basis functions P
⃑
i
j  are zero along the first edge (s
=
-
1), E
⃑
i
2 and P
⃑
i
j being zero along the second edge (s
=
1). Thus, the continuity equation can easily be imposed on a given mesh made of patches.


### 2.3. Parallel Schemes of Higher-Order MoM

The parallelization of the higher-order MoM solution involves parallel matrix filling followed by a parallel solution of the dense matrix equation. To ensure the load-balancing in terms of data and CPU operations, it is necessary to divide the matrix into blocks and distribute those blocks to all processors through a block-cyclic distribution procedure. The modes of distribution of blocks are related to matrix equation solution methods and we select the LU factorization as the solution method.As an example, consider a matrix [A] of Figure 1(a), which is distributed to different processes in the 3 × 3 process grid. Figure 1(b) shows to which process the blocks of [A] are distributed using the block-cyclic distribution methodology [2, 14]. In Figure 1(a), the outermost numbers denote the row and column indices of the process coordinates. The top and bottom number in any block of Figure 1(b) denotes the process ID (rank) and the process coordinates of a certain process, respectively, corresponding to the block of the matrix shown in Figure 1(a). By varying the dimensions of the blocks of [A] and those of the process grid, different mappings can be obtained.Block-cyclic distribution of a matrix: (a) a matrix consisting of 9 × 9 blocks (b) rank and coordinates of each process owning the corresponding blocks in (a).
(a)
(b)During the procedure of finding the solution to the matrix equation, balancing the computational load between all the processes is also important. However, it is less easy to track compared with the parallel matrix filling procedure, since an increase in the number of unknowns or the number of processes executing the solution can increase the amount of communication required between the processes. This increase in communication will decrease the gains of parallel speedup. However, as a rule of thumb, more processes typically means less wall clock time for solving large problems, even though the overall parallel efficiency may sometimes be increased by using fewer processes to solve the same problem.The solution of the matrix equation is essentially the same regardless of the type of basis functions used for the MoM. The parallel LU decomposition solver based on the ScaLAPACK library package is described in Chapter 2 in [2]. Also, one can refer to Chapter 6 in [2] for a detailed discussion about the iterative solvers based on the conjugate gradient (CG) method. The computational complexity of LU decomposition, scales with O (N
3), is much higher than that of CG type methods with O (N
2), where N is the number of unknowns. However, in this paper, the parallel LU decomposition is utilized as the parallel equation solver rather than the parallel iterative CG method due to the fact that the iterative CG method may encounter a divergence problem when dealing with complex antennas composed of thin structures and various materials.


## 2.1. Formulation of the Integral Equations

The method is based on the solution of surface integral equations (SIEs) [8–10] in the frequency domain for equivalent currents over dielectric boundary surfaces and electric currents over perfect electric conductors (PECs). The set of integral equations obtained are solved by using MoM, specifically using the Galerkin method. The method is able to handle inhomogeneous dielectrics categorized by a combination of various homogeneous dielectrics. Therefore, any composite metallic and dielectric structure can be represented as an electromagnetic system consisting of a finite number of finite-size linear, homogeneous, and isotropic regions situated in an unbounded linear, homogeneous, isotropic environment.Let us consider an arbitrary regioni with a nonzero electromagnetic field. According to the surface equivalent theorem, the field outside region i becomes zero. Hence, the region outside region i can be homogenized with respect the region i, and then, a multiple-region problem may be decomposed into n single-region problems. The total electric field in region i may be expressed as a sum of scattered and incident fields as follows:
(1)
E
⃑
(
i
)
=
∑
j
=
0
,
j
≠
i
n
E
⃑
(
i
)
(
J
⃑
s
i
j
,
M
⃑
s
i
j
)
+
E
⃑
inc
(
i
)
,
where E
⃑
(
i
)
(
J
⃑
s
i
j
,
M
⃑
s
i
j
) represents the scattered field inside region i, which is calculated by the currents placed on the boundary surface between regions i and j, and E
⃑
inc
(
i
) is the corresponding incident field. The electric and magnetic scattered fields inside region i are given by
(2)
E
⃑
(
i
)
(
J
⃑
s
i
j
,
M
⃑
s
i
j
)
=
-
L
(
i
)
(
J
⃑
s
i
j
)
+
K
(
i
)
(
M
⃑
s
i
j
)
,
(3)
H
⃑
(
i
)
(
J
⃑
s
i
j
,
M
⃑
s
i
j
)
=
-
K
(
i
)
(
J
⃑
s
i
j
)
-
ε
(
i
)
μ
(
i
)
L
(
i
)
(
M
⃑
s
i
j
)
,
where L
(
i
) and K
(
i
) are linear operators given by
(4)
L
(
i
)
(
X
⃑
s
i
j
)
=
j
ω
μ
(
i
)
∫
S
i
k
(
X
⃑
s
i
j
(
r
⃑
i
j
)
g
(
i
)
(
r
⃑
,
r
⃑
′
)
+
1
ω
2
ε
(
i
)
μ
(
i
)
∇
s
i
j
m
m
m
m
m
m
m
m
m
·
X
⃑
s
i
j
(
r
⃑
i
j
)
∇
g
(
i
)
(
r
⃑
,
r
⃑
′
)
1
ω
2
ε
(
i
)
μ
(
i
)
)
d
S
i
j
,
K
(
i
)
(
X
⃑
s
i
j
)
=
∫
S
i
j
X
⃑
s
i
j
(
r
⃑
i
j
)
×
∇
g
(
i
)
(
r
⃑
,
r
⃑
′
)
d
S
i
j
,
where X
⃑
s
i
j  denotes the electric or magnetic current, r
⃑
′ is the vector position of the source point, r
⃑ is the vector position of the field point, and g
i
(
r
⃑
,
r
⃑
′
) is Green’s function for the homogeneous medium i. Making use of the boundary condition equations given by
(5)
n
^
i
j
×
(
E
⃑
(
i
)
-
E
⃑
(
j
)
)
=
0
,
(6)
n
^
i
j
×
(
H
⃑
(
i
)
-
H
⃑
(
j
)
)
=
0
and replacing (1), (2), and (3) into (4) and (6), the integral equations for regions i and k are obtained in the form
(7)
n
^
i
k
×
{
∑
j
=
0
,
j
≠
i
n
[
L
(
i
)
(
J
⃑
s
i
j
)
-
K
(
i
)
(
M
⃑
s
i
j
)
]
mmmm
-
∑
j
=
0
,
j
≠
k
n
[
L
(
j
)
(
J
⃑
s
k
j
)
-
K
(
j
)
(
M
⃑
s
k
j
)
]
}
=
n
^
i
k
×
(
E
⃑
inc
(
i
)
-
E
⃑
inc
(
k
)
)
,
n
^
i
k
×
{
∑
j
=
0
,
j
≠
i
n
(
ε
(
i
)
μ
(
i
)
L
(
i
)
(
M
⃑
s
i
j
)
+
K
(
i
)
(
J
⃑
s
i
j
)
)
m
m
m
m
-
∑
j
=
0
,
j
≠
k
n
(
ε
(
j
)
μ
(
j
)
L
(
j
)
(
M
⃑
s
k
j
)
+
K
(
j
)
(
J
⃑
s
k
j
)
)
}
=
n
^
i
j
×
(
H
⃑
inc
(
i
)
-
H
⃑
inc
(
k
)
)
.These two sets of equations represent a general form of the Poggio-Miller-Chang-Harrington-Wu (PMCHW) [2, 11] formulation. When the boundary surface of two different regions is PEC, the magnetic currents are equal to zero at the boundary surface and the equations degenerates into the electric field integral equation (EFIE) [12, 13].


## 2.2. Higher-Order Basis Functions

The electric and magnetic currents are approximated by higher-order polynomials, which reduce the number of unknowns compared with the rational piecewise basis functions. The code makes use of truncated cones for wires and bilinear patches to characterize other surfaces.For bilinear surfaces, the surface current is decomposed into itsp- and s-components, p and s being the two parametric coordinates of the unit quadrangle, p
,
s
∈
[
-
1
,
1
]. The approximation for the s-component of the electric current is (analogous expressions stand for the p-component and for the magnetic current)
(8)
J
⃑
s
(
p
,
s
)
=
∑
i
=
0
N
p
[
c
i
1
E
⃑
i
(
p
,
s
)
+
c
i
2
E
⃑
i
(
p
,
-
s
)
+
∑
j
=
2
N
s
a
i
j
P
⃑
i
j
(
p
,
s
)
]
,
where N
p and N
s are the degrees of the approximations along the coordinates, and a
i
j, c
i
1, and c
i
2 are the unknown coefficients.Thus, expression (8) stands for the representation of the current in terms of edge basis functions E
⃑
i
(
p
,
s
) and interior or patch basis functions P
⃑
i
j
(
p
,
s
) which can be compactly expressed as
(9)
E
⃑
i
k
(
p
,
s
)
=
α
s
|
α
p
×
α
s
|
{
p
i
N
(
s
)
,
k
=
1
p
i
N
(
-
s
)
,
k
=
2
,
P
⃑
i
j
(
p
,
s
)
=
α
s
|
α
p
×
α
s
|
p
i
S
j
(
s
)
,
N
(
s
)
=
1
-
s
2
,
S
i
(
s
)
=
{
s
i
-
1
i
is
even
s
i
-
s
i
is
odd
,
where symbols a
s and a
p denote the unitary vectors along the transformed s and p coordinates.Edge basis functionsE
⃑
i
1 and patch basis functions P
⃑
i
j  are zero along the first edge (s
=
-
1), E
⃑
i
2 and P
⃑
i
j being zero along the second edge (s
=
1). Thus, the continuity equation can easily be imposed on a given mesh made of patches.


## 2.3. Parallel Schemes of Higher-Order MoM

The parallelization of the higher-order MoM solution involves parallel matrix filling followed by a parallel solution of the dense matrix equation. To ensure the load-balancing in terms of data and CPU operations, it is necessary to divide the matrix into blocks and distribute those blocks to all processors through a block-cyclic distribution procedure. The modes of distribution of blocks are related to matrix equation solution methods and we select the LU factorization as the solution method.As an example, consider a matrix [A] of Figure 1(a), which is distributed to different processes in the 3 × 3 process grid. Figure 1(b) shows to which process the blocks of [A] are distributed using the block-cyclic distribution methodology [2, 14]. In Figure 1(a), the outermost numbers denote the row and column indices of the process coordinates. The top and bottom number in any block of Figure 1(b) denotes the process ID (rank) and the process coordinates of a certain process, respectively, corresponding to the block of the matrix shown in Figure 1(a). By varying the dimensions of the blocks of [A] and those of the process grid, different mappings can be obtained.Block-cyclic distribution of a matrix: (a) a matrix consisting of 9 × 9 blocks (b) rank and coordinates of each process owning the corresponding blocks in (a).
(a)
(b)During the procedure of finding the solution to the matrix equation, balancing the computational load between all the processes is also important. However, it is less easy to track compared with the parallel matrix filling procedure, since an increase in the number of unknowns or the number of processes executing the solution can increase the amount of communication required between the processes. This increase in communication will decrease the gains of parallel speedup. However, as a rule of thumb, more processes typically means less wall clock time for solving large problems, even though the overall parallel efficiency may sometimes be increased by using fewer processes to solve the same problem.The solution of the matrix equation is essentially the same regardless of the type of basis functions used for the MoM. The parallel LU decomposition solver based on the ScaLAPACK library package is described in Chapter 2 in [2]. Also, one can refer to Chapter 6 in [2] for a detailed discussion about the iterative solvers based on the conjugate gradient (CG) method. The computational complexity of LU decomposition, scales with O (N
3), is much higher than that of CG type methods with O (N
2), where N is the number of unknowns. However, in this paper, the parallel LU decomposition is utilized as the parallel equation solver rather than the parallel iterative CG method due to the fact that the iterative CG method may encounter a divergence problem when dealing with complex antennas composed of thin structures and various materials.


## 3. Numerical Results

### 3.1. Comparison with the Finite Element Method

To validate the accuracy and efficiency of the proposed parallel higher-order MoM, a rectangular microstrip patch antenna element is simulated, which is shown in Figure2. The dimensions of the patch element are 205.6 mm × 154.8 mm and the material parameter of the substrate is ε
r
=
4.2. The operation frequency is 440 MHz and the simulation is performed on a personal computer, with a Quad-core Intel I5 processor (3.00 GHz) and 8 GB of RAM. The radiation pattern obtained using the parallel higher-order MoM and the finite element method are plotted together in Figure 3. The results agree with each other very well.Figure 2
Model of the microstrip patch antenna element.2D gain of microstrip patch antenna: (a)y
o
z plane and (b) x
o
y plane.
(a)
(b)


### 3.2. Description of the Parallel Computational Platform

The computing cluster used in this section is the Magic Cube with 64 nodes at Shanghai Supercomputer Center (SSC), and each node has four quad-core 1.9 GHz processors and 64 GB memory. The nodes are connected with Infiniband switches.


### 3.3. A Microstrip Patch Array Antenna with37
×
9 Elements

In this section, we present the radiation patterns of the algorithm for a microstrip array. Consider a rectangular microstrip patch array antenna with37
×
9 elements. The microstrip element parameters are the same as the element analyzed in Section 3.1 and the array model is shown in Figure 4. The dimensions of the array are 10 m × 2.5 m × 0.018 m and the distances between any two neighboring elements are 57.9 mm and 83.3 mm along the length and width directions. The amplitude at the feed of the array is designed by a −35 dB Taylor distribution, and the operation frequency of the array is 440 MHz. The three-dimensional (3D) gain pattern is given in Figure 5 and the two-dimensional (2D) gain patterns are given in Figure 6. The side-lobe levels in both the y
o
z and x
o
y planes are below −35 dB, which meet the design target. Note that θ coordinate is measured from x
o
y plane to z axis and φ coordinate is measured from +
x axis to y axis in this paper.Figure 4
Model of the microstrip patch array antenna with 37 × 9 elements.Figure 5
3D gain patterns of the microstrip patch array antenna.2D gain patterns of the microstrip patch array antenna: (a)y
o
z plane and (b) x
o
y plane.
(a)
(b)


### 3.4. A Microstrip Antenna Array with37
×
9 Elements over a Large Airplane

Install the 37 × 9 microstrip patch array that was analyzed in Section3.3 over a large airplane (http://www.scalecraft.com/737aewandc.aspx), as shown in Figure 7. The dimensions of the airplane are 55 m × 47.6 m × 15.8 m. The array elements are very thin compared with the airplane platform and the airborne antenna array is a composite and multiscale problem. Consider a realistic problem that the microstrip antenna array installed on nine positions of the airplane with different x coordinates or different z coordinates, as shown in Figure 8, the mainlobe direction is directed towards the wing. The operation frequency of the airborne array is 440 MHz and the number of unknowns is 308,371.Figure 7
Microstrip patch array over an airplane.Figure 8
Nine positions of the array in thexoz plane.In this Section, 1024 CPU cores are used and the direct solver is observed in solving this kind of problems to avoid the slow convergence of iterative CG type methods. The matrix filling time is about 4700 seconds, and the matrix equation solving time is about 18500 seconds. The computed three-dimensional (3D) gain patterns of the antenna array installed on different positions of the airplane are shown in Figure9. To make the comparison more clearly, the computed two-dimensional (2D) gain patterns of the antenna array installed on different positions of the airplane are compared in two ways, as shown in Figure 10. One is the same x coordinates and different z coordinates, as shown in Figure 11. One is the different x coordinates and same z coordinates, as shown in Figure 12. Furthermore, the patterns of the array alone are also given for comparison.3D gain patterns of the nine installed positions of the array. (a) Array at position 1. (b) Array at position 2. (c) Array at position 3. (d) Array at position 4. (e) Array at position 5. (f) Array at position 6. (g) Array at position 7. (h) Array at position 8. (i) Array at position 9.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)Two comparison ways for the 2D gain pattern.
(a)
(b)2D gain patterns of the array installed on samex coordinates and different z coordinates. (a) y
o
z plane of positions 1, 2, and 3. (b) x
o
y plane of positions 1, 2, and 3. (c) y
o
z plane of positions 4, 5, and 6. (d) x
o
y plane of positions 4, 5, and 6. (e) y
o
z plane of positions 7, 8, and 9. (f) x
o
y plane of positions 7, 8, and 9.
(a)
(b)
(c)
(d)
(e)
(f)2D gain patterns of the array installed on differentx coordinates and same z coordinates. (a) y
o
z plane of positions 1, 4, and 7. (b) x
o
y plane of positions 1, 4, and 7. (c) y
o
z plane of positions 2, 5, and 8. (d) x
o
y plane of positions 2, 5, and 8. (e) y
o
z plane of positions 3, 6, and 9. (f) x
o
y plane of position 3, 6, and 9.
(a)
(b)
(c)
(d)
(e)
(f)Form the comparison, it is clearly seen that the gain patterns of the airborne antenna array in they
o
z planes (elevation planes) seriously deteriorate compared with the antenna array alone, especially in the mainlobe region. And, in contrast, the patterns in the x
o
y plane (azimuth plane) show slight changes. From the comparisons results of the array installed on same x coordinates and different z coordinates shown in Figure 11, one can see that, with the z coordinates of the installation positions increasing (i.e., form position 1 to 2 and 3, form position 4 to 5 and 6, and form position 7 to 8 and 9), the mainlobe region of gain patterns (i.e., −30°~30° in the y
o
z plane) gets better mainly due to the reflection and diffraction effects of airplane wing becoming weak. Furthermore, the sidelobe region of gain patterns (i.e., −150°~−30° in the y
o
z plane, 115°~120° in the x
o
y plane) becomes low with the installation height increasing. Form the comparisons results of the array installed on different x coordinates and same z coordinates shown in Figure 12, we can see that gain patterns (i.e., −150°~30° in the y
o
z plane) get better with the x coordinates of the installation positions increasing (i.e., form position 1 to 4 and 7, form position 2 to 5 and 8, and form position 3 to 6 and 9). However, the overall design of airplane is a composite and multiscale problem, it is necessary to take many factors into account, such as the aerodynamics. Therefore, the array is installed not too high and not close to the head of airplane; taking these factors into consideration, the better installation position is position 5 in this paper.Moreover, the model of airborne array is a whole in practical engineering. The array is connected to airplane with a bracket, which is installed below the array. When the main beam of the array is direct to the wings of the airplane, the effect of the bracket can be neglected. So in the simulations of the airborne array in this section, the bracket is removed from the simulation model.


### 3.5. A Microstrip Antenna Array with9
×
1 Elements stick on an Airplane

To demonstrate the ability of the method to deal with the problem of an array connected to the airplane, a microstrip antenna array with 9 × 1 elements sticks on the surface of a large airplane is considered in this section [15]. The model of the array and the airborne array are shown in Figures 13(a) and 13(b) and the position of the array is given in Figure 13(c). The parameters of the array and the airplane are same to the previous sections. The computed 3D gain patterns of the array alone and airborne array are shown in Figure 14. The computed 2D gain patterns of the array alone and airborne array are compared in Figure 15.Simulation models: (a) microstrip array; (b) airborne array; (c) the position of the array.
(a)
(b)
(c)3D gain patterns of the array alone and airborne array: (a) array alone; (b) airborne array.
(a)
(b)Figure 15
Comparison result of 2D gain patterns inxoz plane.Form the comparisons, we can see that when the array sticks on the airplane, the back gain of the array is significantly reduced, while the front gain shows slight changes. It is mainly due to that when the array is stuck on the airplane, it can be considered as enlarging the ground of the array, so the result in Figure15 is reasonable.


## 3.1. Comparison with the Finite Element Method

To validate the accuracy and efficiency of the proposed parallel higher-order MoM, a rectangular microstrip patch antenna element is simulated, which is shown in Figure2. The dimensions of the patch element are 205.6 mm × 154.8 mm and the material parameter of the substrate is ε
r
=
4.2. The operation frequency is 440 MHz and the simulation is performed on a personal computer, with a Quad-core Intel I5 processor (3.00 GHz) and 8 GB of RAM. The radiation pattern obtained using the parallel higher-order MoM and the finite element method are plotted together in Figure 3. The results agree with each other very well.Figure 2
Model of the microstrip patch antenna element.2D gain of microstrip patch antenna: (a)y
o
z plane and (b) x
o
y plane.
(a)
(b)


## 3.2. Description of the Parallel Computational Platform

The computing cluster used in this section is the Magic Cube with 64 nodes at Shanghai Supercomputer Center (SSC), and each node has four quad-core 1.9 GHz processors and 64 GB memory. The nodes are connected with Infiniband switches.


## 3.3. A Microstrip Patch Array Antenna with37
×
9 Elements

In this section, we present the radiation patterns of the algorithm for a microstrip array. Consider a rectangular microstrip patch array antenna with37
×
9 elements. The microstrip element parameters are the same as the element analyzed in Section 3.1 and the array model is shown in Figure 4. The dimensions of the array are 10 m × 2.5 m × 0.018 m and the distances between any two neighboring elements are 57.9 mm and 83.3 mm along the length and width directions. The amplitude at the feed of the array is designed by a −35 dB Taylor distribution, and the operation frequency of the array is 440 MHz. The three-dimensional (3D) gain pattern is given in Figure 5 and the two-dimensional (2D) gain patterns are given in Figure 6. The side-lobe levels in both the y
o
z and x
o
y planes are below −35 dB, which meet the design target. Note that θ coordinate is measured from x
o
y plane to z axis and φ coordinate is measured from +
x axis to y axis in this paper.Figure 4
Model of the microstrip patch array antenna with 37 × 9 elements.Figure 5
3D gain patterns of the microstrip patch array antenna.2D gain patterns of the microstrip patch array antenna: (a)y
o
z plane and (b) x
o
y plane.
(a)
(b)


## 3.4. A Microstrip Antenna Array with37
×
9 Elements over a Large Airplane

Install the 37 × 9 microstrip patch array that was analyzed in Section3.3 over a large airplane (http://www.scalecraft.com/737aewandc.aspx), as shown in Figure 7. The dimensions of the airplane are 55 m × 47.6 m × 15.8 m. The array elements are very thin compared with the airplane platform and the airborne antenna array is a composite and multiscale problem. Consider a realistic problem that the microstrip antenna array installed on nine positions of the airplane with different x coordinates or different z coordinates, as shown in Figure 8, the mainlobe direction is directed towards the wing. The operation frequency of the airborne array is 440 MHz and the number of unknowns is 308,371.Figure 7
Microstrip patch array over an airplane.Figure 8
Nine positions of the array in thexoz plane.In this Section, 1024 CPU cores are used and the direct solver is observed in solving this kind of problems to avoid the slow convergence of iterative CG type methods. The matrix filling time is about 4700 seconds, and the matrix equation solving time is about 18500 seconds. The computed three-dimensional (3D) gain patterns of the antenna array installed on different positions of the airplane are shown in Figure9. To make the comparison more clearly, the computed two-dimensional (2D) gain patterns of the antenna array installed on different positions of the airplane are compared in two ways, as shown in Figure 10. One is the same x coordinates and different z coordinates, as shown in Figure 11. One is the different x coordinates and same z coordinates, as shown in Figure 12. Furthermore, the patterns of the array alone are also given for comparison.3D gain patterns of the nine installed positions of the array. (a) Array at position 1. (b) Array at position 2. (c) Array at position 3. (d) Array at position 4. (e) Array at position 5. (f) Array at position 6. (g) Array at position 7. (h) Array at position 8. (i) Array at position 9.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)Two comparison ways for the 2D gain pattern.
(a)
(b)2D gain patterns of the array installed on samex coordinates and different z coordinates. (a) y
o
z plane of positions 1, 2, and 3. (b) x
o
y plane of positions 1, 2, and 3. (c) y
o
z plane of positions 4, 5, and 6. (d) x
o
y plane of positions 4, 5, and 6. (e) y
o
z plane of positions 7, 8, and 9. (f) x
o
y plane of positions 7, 8, and 9.
(a)
(b)
(c)
(d)
(e)
(f)2D gain patterns of the array installed on differentx coordinates and same z coordinates. (a) y
o
z plane of positions 1, 4, and 7. (b) x
o
y plane of positions 1, 4, and 7. (c) y
o
z plane of positions 2, 5, and 8. (d) x
o
y plane of positions 2, 5, and 8. (e) y
o
z plane of positions 3, 6, and 9. (f) x
o
y plane of position 3, 6, and 9.
(a)
(b)
(c)
(d)
(e)
(f)Form the comparison, it is clearly seen that the gain patterns of the airborne antenna array in they
o
z planes (elevation planes) seriously deteriorate compared with the antenna array alone, especially in the mainlobe region. And, in contrast, the patterns in the x
o
y plane (azimuth plane) show slight changes. From the comparisons results of the array installed on same x coordinates and different z coordinates shown in Figure 11, one can see that, with the z coordinates of the installation positions increasing (i.e., form position 1 to 2 and 3, form position 4 to 5 and 6, and form position 7 to 8 and 9), the mainlobe region of gain patterns (i.e., −30°~30° in the y
o
z plane) gets better mainly due to the reflection and diffraction effects of airplane wing becoming weak. Furthermore, the sidelobe region of gain patterns (i.e., −150°~−30° in the y
o
z plane, 115°~120° in the x
o
y plane) becomes low with the installation height increasing. Form the comparisons results of the array installed on different x coordinates and same z coordinates shown in Figure 12, we can see that gain patterns (i.e., −150°~30° in the y
o
z plane) get better with the x coordinates of the installation positions increasing (i.e., form position 1 to 4 and 7, form position 2 to 5 and 8, and form position 3 to 6 and 9). However, the overall design of airplane is a composite and multiscale problem, it is necessary to take many factors into account, such as the aerodynamics. Therefore, the array is installed not too high and not close to the head of airplane; taking these factors into consideration, the better installation position is position 5 in this paper.Moreover, the model of airborne array is a whole in practical engineering. The array is connected to airplane with a bracket, which is installed below the array. When the main beam of the array is direct to the wings of the airplane, the effect of the bracket can be neglected. So in the simulations of the airborne array in this section, the bracket is removed from the simulation model.


## 3.5. A Microstrip Antenna Array with9
×
1 Elements stick on an Airplane

To demonstrate the ability of the method to deal with the problem of an array connected to the airplane, a microstrip antenna array with 9 × 1 elements sticks on the surface of a large airplane is considered in this section [15]. The model of the array and the airborne array are shown in Figures 13(a) and 13(b) and the position of the array is given in Figure 13(c). The parameters of the array and the airplane are same to the previous sections. The computed 3D gain patterns of the array alone and airborne array are shown in Figure 14. The computed 2D gain patterns of the array alone and airborne array are compared in Figure 15.Simulation models: (a) microstrip array; (b) airborne array; (c) the position of the array.
(a)
(b)
(c)3D gain patterns of the array alone and airborne array: (a) array alone; (b) airborne array.
(a)
(b)Figure 15
Comparison result of 2D gain patterns inxoz plane.Form the comparisons, we can see that when the array sticks on the airplane, the back gain of the array is significantly reduced, while the front gain shows slight changes. It is mainly due to that when the array is stuck on the airplane, it can be considered as enlarging the ground of the array, so the result in Figure15 is reasonable.


## 4. Conclusion

An efficient parallel higher-order MoM is developed to solve extremely large and complex electromagnetics problems. The parallel strategies based on block matrices and the higher-order basis functions make it possible to solve airborne antenna array problems using rigorous MoM, and the simulation time of one installed position is about seven hours using 1024 CPU cores. Numerical results confirm the accuracy and efficiency of the method and its applicability to the analysis of radiation from airborne platform systems with complex antennas.



---

*Source: 985367-2014-07-17.xml*