Fast and Accurate Simulation Technique for Large Irregular Arrays

TL;DR

This paper introduces a fast, accurate simulation method for large irregular antenna arrays using macro basis functions, interpolatory techniques, and FFT acceleration, validated on SKA radio telescope arrays.

Contribution

It presents a novel combination of MBF, HARP models, and FFT-based acceleration for efficient large array analysis, addressing a gap in existing simulation methods.

Findings

Significant reduction in computation time for large arrays.

Validated accuracy against commercial software and experiments.

Effective analysis of SKA radio telescope array configurations.

Abstract

A fast full-wave simulation technique is presented for the analysis of large irregular planar arrays of identical 3-D metallic antennas. The solution method relies on the Macro Basis Functions (MBF) approach and an interpolatory technique to compute the interactions between MBFs. The Harmonic-polynomial (HARP) model is established for the near-field interactions in a modified system of coordinates. For extremely large arrays made of complex antennas, two approaches assuming a limited radius of influence for mutual coupling are considered: one is based on a sparse-matrix LU decomposition and the other one on a tessellation of the array in the form of overlapping sub-arrays. The computation of all embedded element patterns is sped up with the help of the non-uniform FFT algorithm. Extensive validations are shown for arrays of log-periodic antennas envisaged for the low-frequency SKA (Square Kilometer Array) radio-telescope. The analysis of SKA stations with such a large number of elements has not been treated yet in the literature. Validations include comparison with results obtained with commercial software and with experiments. The proposed method is particularly well suited to array synthesis, in which several orders of magnitude can be saved in terms of computation time.

Figures

Tables

Table 1. TABLE I: Computational Time Obtaining Embedded Element Patterns of a SKA1-Low station at 110 MHz

	CST	WIPL-D	HARP
Simulation Time	96 hours	97 hours	0.5 min.

Table 2. TABLE II: Analysis Time for a SKA station using HARP at 110 MHz.

Operations	Time (in seconds)
Matrix Filling	10
Solving MoM equations	12.6
Calculation EEPs	6
Total	28.6

Table 3. TABLE III: computational Complexity by Sparse Matrix and Tessellation approaches.

	Sparse Matrix	Tessellation
Filling Matrix	$𝒪(N_a{N}_I{N}_m^2$)	$𝒪\left(N_a{N}_s{N}_m^2\right)$
Solving	$𝒪\left(N_a{N}_I^2{N}_m^3\right)$	$𝒪\left(N_a{N}_s^2{N}_m^3\right)$

Equations

$$Z_{\mathrm{TS}}^{\mathrm{app}}(r_{mn},\hat{\alpha })\approx -j\omega \mu {\vec{F}}_{T,m}^{\hat{\alpha },\ast }.{\vec{F}}_{S,n}^{\hat{\alpha }}\frac{e^{-jk{r}_{mn}}}{4\pi r_{mn}}$$

$$B_{\mathrm{TS}}(r_{mn},\hat{\alpha })=\frac{Z_{\mathrm{TS}}^{\mathrm{exact}}(r_{mn},\hat{\alpha })-Z_{\mathrm{TS}}^{\mathrm{app}}(r_{mn},\hat{\alpha })}{e^{-jk{r}_{mn}}}$$

$$B_{\mathrm{TS}}(d,\hat{\alpha })=\sum _{p=-P}^P{e}^{jp\alpha }\sum _{q=0}^Q{c}_{pq}{d}^q$$

$$Z_{\mathrm{TS}}(r_{mn},\hat{\alpha })=Z_{\mathrm{TS}}^{\mathrm{app}}(r_{mn},\hat{\alpha })+B_{\mathrm{TS}}(r_{mn},\hat{\alpha })e^{-jk{r}_{mn}}$$

$${\vec{F}}_{array}(\theta ,\varphi )$$

$$=\sum _{j=1}^{N_m}{\vec{F}}_j(\theta ,\varphi )\sum _{i=1}^{N_a}{c}_{ij}{e}^{jk(u_x{x}_i+u_y{y}_i)}$$

$$C=\alpha \frac{N_a}{{{\mathrm{R}}_{\mathrm{in}}}^2}({\mathrm{R}}_{\mathrm{in}}+{\mathrm{R}}_{\mathrm{ext}}{)}^4$$

$$ϵ=\sqrt{{\left|E_{\theta }-E_{\theta }^{ref}\right|}^2+{\left|E_{\varphi }-E_{\varphi }^{ref}\right|}^2}$$

Full text

I Introduction

Modern radars [1], massive MIMO [2] and radio astronomy arrays [3, 4] are in need of accurate modeling of their main electromagnetic properties. These are strongly impacted by mutual coupling and include the array impedance matrix and all embedded element patterns (EEP), i.e. the patterns of each individual elements in the presence of all other antennas. This can become a challenging task, knowing that arrays can consist of hundreds or thousands of elements, as considered in the SKA project [5, 6] and in modern radar systems [1].

The arrays can be regular as well as non-regular. Irregular arrays exhibit a series of advantages over regular arrays, namely the randomization of the side-lobes and of the mutual coupling effects [7]. Both of these factors allow the operation of the arrays over a much larger frequency band compared to traditional regular arrays [8], with a limited number of elements. The accurate knowledge of the antenna and array beams, as well as of the elements’ impedance in the presence of mutual coupling [9] is essential for the design of these systems, i.e. for the design of the antenna elements, the array configuration and the electronics connected to the antennas, especially the first stage amplification. The number of elements in the arrays, the geometrical complexity of the elements and the electrical size of these systems represent a challenge regarding their computer modeling. In particular, the irregularity of the elements’ positions in the array prevents us from using certain approximations, such as periodic boundary conditions possibly combined with corrections for truncation [10]. The present paper proposes a possible solution based on two novel contributions: (1) a fast full-wave technique able to accurately simulate large irregular arrays of complex structures in a matter of minutes; and (2) an analysis scheme that exploits this fast simulation technique for the analysis of extremely large arrays consisting of hundreds or thousands of elements.

Integral equation techniques solved with the method of moments (MoM) are very convenient for antenna problems, since they inherently account for radiation. However, when one uses Gaussian elimination, the MoM solution has a complexity of $O(({N_{a}}{N_{e}})^{3})$, where ${N_{a}}$ is the number of antennas in the array, and ${N_{e}}$ is the number of elementary basis functions per element. For large arrays of complex antennas, the multilevel fast multipole (MLFMA) technique [11], and its hybridization with other techniques (e.g. HOMoM-MLFMA [12]) are often exploited in an iterative scheme. Unfortunately, the calculation of the array impedance matrix and all EEPs in large finite arrays requires a new series of iterations for excitation at each antenna port. This often leads to an unaffordable computational effort for large arrays, besides the uncertainty related to the number of iterations. A non-iterative technique is then favorable, provided that it can benefit from the proper compression of unknowns and from further accelerations. That is the case for techniques such as Macro Basis Functions (MBF) [13], Characteristic Basis Functions (CBF) [14, 15], Synthetic Function Expansion (SFX) [16, 17], or entire/sub-entire-domain (SED) basis functions [18, 19].

MBF-type techniques achieve a compression of the MoM impedance matrix by replacing the original set of elementary basis functions by a new set of functions obtained through the solution of smaller problems. Once the MoM matrix size is reduced by means of the MBF technique, the $O((N_{a}N_{e})^{3})$ complexity for solving the system of equations is reduced to $O((N_{a}N_{m})^{3})$, where $N_{m}$ is the number of MBFs and $N_{m}<<N_{e}$. Nevertheless, the impedance matrix filling time remains $O((N_{a}N_{e})^{2})$, the same as for the original MoM, which is prohibitive for large arrays. To overcome this obstacle, different techniques have been proposed, such as the multipole expansion algorithm [20, 21], the interpolatory technique [22], the perturbation technique based on the domain Green’s function [23], or the combination of multipole and Contour-FFT approaches for printed structures [24]. The multipole expansion, however, suffers from the instability of near-field interaction, which appears when antennas are too close to each other, i.e. for distances smaller than about one wavelength. For dense wideband arrays, such as those envisaged for the SKA, antennas are very close to each other at low-frequencies, which requires additional efforts to explicitly calculate the near-field interactions. In this respect, the interpolatory technique presented in [22] is very effective in providing highly accurate solutions even in the near-field [25]. Moreover, MBF interactions are rapidly evaluated from a model, which is built beforehand and is independent of the array configuration. The method, renamed later as HARP [26], has shown great capability and has been extended to the analysis of antennas made of wires, in order to model the SKA1 - Low arrays [8, 27].

Preliminary results of HARP on SKA1-Low antennas have been presented in recent conference papers [25, 26, 28]. Nonetheless, limited details were presented about the technique itself and its performance. The present paper further extends this method including a detailed description of the proposed technique, new improvements and extended validation. HARP is now capable of simulating SKA stations over the whole frequency band (50 MHz to 350 MHz) in a matter of minutes for each frequency point.

Additionally, two approaches are presented to exploit the fast simulation method for the analysis of extremely large arrays, namely, the sparse matrix approach and the tessellation approach. These methods make use of the so-called radius of influence (RI) concept, which is often applied to the analysis of large arrays [29]. The first approach calculates only the interaction between antennas whose relative distance is smaller than the RI, which leads to a sparse interaction matrix that can be effectively LU-decomposed and solved exploiting sparse-matrix libraries. In the tessellation approach, the full array is first partitioned into overlapping subarrays. The EEPs are retained only for the internal (or non-overlapping) part of the subarray. The size of the internal part and of the extended surrounding area of the subarray are determined so as to minimize the computational cost. As each subarray can be analyzed independently, this approach enables the simulation of very large arrays of complex elements on a normal desktop computer; it also facilitates parallel implementation. Finally, as explained above, the MBF approach allows the fast calculation of induced currents for excitation at each port. The calculation of corresponding EEPs is then expressed as superposition of pattern multiplication problems, one per MBF [30], and the array factors calculation is accelerated using the non-uniform FFT (NFFT) [30] – [32].

In this paper, the SKA1-Low arrays are taken as a case study to extensively validate the performance of HARP. Upon its completion, the SKA1-Low instrument will consist of 512 stations, each with 256 log-periodic antennas (SKALA) [8, 27], (see Fig. 1), arranged in a pseudo-random configuration, as shown in Fig. 2. With an ultra-wideband characteristic (7:1 fractional bandwidth) and its complex structure, the modeling and simulation of the SKA1-Low frequency array presents a challenge as the array is quite dense at low frequencies, while at high-frequencies, the antenna needs to be meshed precisely. The accurate modeling of the SKALA antenna has a strong impact on the results obtained with a given simulation software. The results obtained using HARP will be validated with those calculated using commercial software, and with measurements of the isolated elements and a small 3$\times$3 array.

The remainder of the paper is organized as follows. In Section II, the full-wave technique and improvements made to HARP are explained. The sparse matrix and tessellation approaches for extremely large arrays are detailed as well. In Section III, numerical results on the SKA1-Low case are presented. Data obtained using HARP are compared with those from CST Studio [33], WIPL-D [34], and measurements for some occasions. In Section IV, the performance of the tessellation scheme on very large arrays is examined and compared with those from the traditional RI approach on different arrays of SKALA antennas. Finally, conclusions are drawn in Section V.

II Fast Full-wave Technique

The full-wave approach proposed in this paper is based on the interpolatory method presented in [22], renamed as HARP in [26], which uses an interpolatory approach to effectively deal with the analysis of metallic antenna arrays in free space. In this section, we will explain in detail the HARP method, and present improvements that make HARP able to provide with the properties of complex radiators, such as the SKALA antenna. The tessellation approach, which exploits HARP, is also introduced; it enables the effective simulation of extremely large arrays. For the sake of clarity, the interpolatory method is first recalled. Interested readers are referred to [22] for detailed explanations regarding the interpolatory technique. In the present work, the MBFs are obtained using the “primary-and-secondaries” approach [14], where the primary MBF corresponds to the current distribution on an isolated element and the secondary MBFs are the currents induced on a set of neighboring antennas by the element that supports the primary MBF.

II-A Background of HARP

The interpolatory approach was proposed in [22] to quickly compute the interaction matrix between MBFs, which are the entries of the reduced MoM matrix. The interaction between a source MBF $S$ on antenna $n$ placed at the origin and a testing MBF $T$ on antenna $m$, positioned at $(r_{mn},\hat{\alpha})$, was first approximated as

[TABLE]

where $\hat{\alpha}$ points from the reference point in the source antenna to the reference point in the testing antenna and $r_{mn}$ represents the inter-element distance, $\omega$ is the angular frequency, $\mu$ is the free-space permeability, $\vec{F}_{T,m}^{\hat{\alpha}},\vec{F}_{S,n}^{\hat{\alpha}}$ are the radiation patterns of the testing and source MBFs, respectively, and ($*$) stands for complex conjugate. The actual interaction is then obtained as a sum of the above far-field approximation and a complementary part, which is represented using a harmonic-polynomial model. The model is built by firstly calculating explicitly the exact interaction between MBFs over a grid with a limited number of points. A mathematical representation of the difference between the exact and far-field approximated interaction is then constructed after phase correction and an appropriate change of variables. The first step can be expressed in a compact form as

[TABLE]

where $Z_{\rm{TS}}^{\rm{exact}}(r_{mn},\hat{\alpha})$ is the exact interaction explicitly calculated on a pre-defined polar-radial grid. The result of Eq. ( 2) is then modeled in a harmonic-polynomial (HARP) form [22] with a proper change of variables, i.e. $d=1/r^{2}$:

[TABLE]

where $P$ and $Q$ are the orders of Fourier Series and polynomial, respectively, and $c_{pq}$ are coefficients calculated in the least-squares sense. Once the model has been built, the interaction between MBFs is rapidly calculated by adding up the far-field contribution and the complementary part evaluated from the model.

[TABLE]

II-B Improved HARP model

Improvements have been introduced to the original version of HARP to capture the interaction of complex structures such as the SKALA element. The first enhancement appears at the change of variables step, where $d=1/r$ is now implemented instead of $d=1/r^{2}$ as in Eq. (3). This modification is physically interpreted as taking into account not only the contribution of $1/r^{2}$, but also of $1/r^{3}$ field components [35], which are non-negligible in the case of SKALA. Such a change of variables is also consistent with the multipole-based model in [25], where a Laurent series versus the $r$ variable is obtained. The latter work, however, suffers from the stability issues in the very near field. Therefore, in this case, we followed [22] to build the HARP model using the new change of variables.

Additionally, improvements have been added to the implementation of HARP. The DFT approach has been used instead of a least-squares estimator to obtain the harmonic coefficients ($c_{pq}$). Also, the sampling points on the grid are determined by an appropriate sampling routine, i.e. a regular grid on angular variable, and inversely regular along the radial distance (regular in $1/r$ scale). As it will be shown in Section III, HARP is now able to accurately model the interaction between the MBFs of the SKALA over the whole frequency band, while satisfying the minimum distance between elements in the SKA-station.

Regarding the computational aspect using HARP, the complexities for filling the reduced matrix and solving the system of equations are $O(N_{a}^{2}N_{m}^{2}S)$ and $O((N_{a}N_{m})^{3})$, respectively, where $S$ is a small factor related to the calculation of Eq. (4) [22]. The number of MBFs is generally smaller by about two orders of magnitude than the number of basis functions used to describe the current on the antenna (related to meshing); the drastic reduction in the number of unknowns enables the acceleration of the analysis using HARP. In [22], it is shown that interpolatory (or HARP) technique is more efficient than the multipole expansion [20]. Moreover, HARP is able to accurately evaluate the near-field interaction for closely packed arrays.

II-C Pattern Calculations

The use of MBFs allows one to write the pattern of each antenna as a sum of MBF patterns multiplied by the corresponding coefficients for the MBF. Hence, the array pattern can be expressed as:

[TABLE]

where $(x_{i},y_{i})$ is the position of the $i^{\rm{th}}$ antenna, $\hat{u}=(u_{x},u_{y},u_{z})$ is a unit vector in the direction $(\theta,\phi)$ of interest, $\vec{F}_{j}(\theta,\phi)$ is the pattern of MBF $j$, and $c_{ij}$ is the coefficient for MBF $j$ on antenna $i$. The swapping of the summations is possible since only one set of MBFs is used for all-identical-antennas of the array. As the MBF patterns in Eq. (6) are factored out, the second summation becomes the array factor, which can be rapidly evaluated using the N-FFT [31, 30, 32] for irregular arrays. It is noted that since the MBF patterns are smooth functions, they are evaluated over a much sparser grid in advance. The higher resolution patterns then can be rapidly obtained using 2-D interpolation on the grid used for the N-FFT. This approach is also applied for the calculation of each EEP of the array, where the induced currents are quickly calculated for excitation at a given antenna port.

II-D Analysis of Very Large Arrays

For extremely large arrays, comprising thousands of antennas, the direct implementation of HARP on the whole array would require too much memory and computation time. An effective approach exploiting HARP for such arrays is then needed. Generally, the so-called radius of influence (RI) concept is often applied [29, 23] for the analysis of very large arrays; here the interactions are only considered when the separation between elements is smaller than the RI. This results in a sparse matrix for the reduced MoM system of equations. The filling time of this sparse matrix and the solution time of the system of equations are significantly reduced owing to the use of sparse operations. More precisely, sparse matrix libraries [36] are applied to obtain a sparse LU-decomposition of the matrix. Fast solutions can then be computed for any excitation law. The approach, however, requires a large RI [29] in order to obtain a highly accurate pattern, which implies an increase of the computational cost.

An alternative technique is proposed below; it relies on the tessellation of the array into subarrays. The RI concept is then applied to each subarray individually. The idea is illustrated in Fig. 3, where the elements in the array are clustered into subarrays, i.e. hexagonal shapes with the inner non-overlapping part of radius $\rm{R_{in}}$ (blue dashed line), and an extended area between the inner radius $\rm{R_{in}}$ and the outer one $\rm{R_{in}}+\rm{R_{ext}}$ (green area). The EEPs in the array are obtained by analyzing separately each subarray including the inner and the outer parts, and retaining only the EEPs of elements within the inner part. The tessellation approach deals with subarrays having much smaller size than the whole array, which will ease the computational effort of solving the MoM system of equations. The use of an outer part with a width of $R_{\rm{ext}}$ at least equal to $R_{\rm{in}}$ helps to increase the accuracy of the results.

An optimal choice of $\rm{R_{ext}}$ can be made by analyzing the computational complexity of the approach. Assuming that a random array has a constant antenna density, the number of antennas on a subarray (inner and outer parts) will be proportional to $\rm{(R_{in}+R_{ext})^{2}}$. The number of subarrays will be of order of $N_{a}/\rm{{R_{in}^{2}}}$. Following the calculation in [22], the computational effort for filling the matrix can be expressed as a function of $\rm{R_{in}}$:

[TABLE]

where $\alpha$ is a constant taking into account all other constant factors such as the number of MBFs. $C$ presents a minimum when $\rm{R_{in}=R_{ext}}$. The same minimum can be found regarding the complexity of solving the system of equations. This gives a rule of thumb for the choice of $\rm{R_{in}}$ and $\rm{R_{ext}}$ so as to optimize the quality of the solution while minimizing the computational complexity.

One of the important advantages of HARP is the re-usability of the model, which is established only once and can be exploited to analyze any array made of the same element. This feature is very helpful for the tessellation approach, as in this scenario, several subarrays can be simulated separately, which straightforwardly allows parallelization. The implementation and evaluation of HARP using the sparse matrix and tessellation approaches will be studied in Section IV.

III Numerical and Experimental Validation

Preliminary data at a single frequency has been presented in recent conference papers [26], [28] with limited results on relatively small arrays and on a SKA1-Low station. In the present paper, the performance of HARP is studied in depth on several examples, which include an isolated antenna, a small compact array, and the SKA station, shown in Fig. 2, at 50 MHz, 110 MHz, 200 MHz and 350 MHz, covering the whole SKA1-Low band. Simulated results from two commercial software, CST Microwave Studio [33] and WIPL-D [34], are taken as a reference to examine the performance of HARP, while measured data is presented whenever they are available. Different parameters are analyzed, including input impedance, coupling coefficients, EEPs and array patterns. In all cases, the data for one polarization is shown, while the results for the other polarization have the same quality and are not shown. The SKALA model used in HARP consists in a wire mesh comprising 1218 basis functions. This HARP implementation is based on a MoM code that considers wires, with the Poklington approximation [37]. The 4 cm diameter tube that supports each arm does not satisfy that approximation in the upper part of the frequency range. Therefore, those tubes have been replaced by wire cages with a square cross-section and a perimeter equal to that of the tube. Fig. 1 shows a global view of the structure, while Fig. 4 shows a zoom near the feed points and near the cage. The numbers of MBFs are 15, 20, 45, 70 at 50 MHz, 110 MHz, 200 MHz and 350 MHz, respectively (linearly scaled versus frequency).

III-A Self and Mutual Impedance

Due to the complexity of the design, an accurate model of the SKALA in Fig. 1 is critical for analyzing the performance of the array, e.g. self and mutual impedances are required to study parameters such as noise coupling [38, 39]. The input impedance is then the first parameter to study. Fig. 5 shows the self-impedance simulated using different software, along with the measured data over the low-frequency band. A small deviation of the order of $20\,\Omega$ is seen around 200–300 MHz between the simulations and the measurements. This can be attributed to the fact that different software packages model and mesh the SKALA differently, and that the SKALA in the measurement had a metallic ground limited to 2 m, whereas the simulations assume an infinite ground plane. These variations may also explain some of the differences among patterns in the following sections.

Another validation example carries on a 3$\times$3 regular array of SKALAs with element spacing of 1.5 m; as shown in Fig. 6. The coupling coefficient is calculated between the center antenna and another element at the middle of the edge, as shown in Fig. 7. This example is a repetition of an example in [26] with an improved version of HARP used here. A clear improvement is seen for HARP w.r.t. the result in [26] as the HARP and CST simulations follow well the measurement across the whole band.

The mutual impedance is studied for the SKA station as well. Fig. 8 plots the real and imaginary parts of the mutual impedance $Z_{\rm{ij}}$ for different frequencies in a SKA station between an arbitrary antenna ${i}$ (a red dot in Fig. 2) and another antenna ${j}$ running from 1 to 256; the obtained values are sorted versus distance to antenna ${i}$. A good agreement is found between HARP and WIPL-D for coupling between antennas. The main difference regards the self-impedance, i.e. zero distance, which is also illustrated in Fig. 5. The results confirm the capability of the current version of HARP to capture the interaction between SKALA antennas in their array environment.

III-B Embedded Element Patterns

In applications such as the SKA telescope, the EEPs are of primary interest as they allow one to form different beams to get information of signal from sky. An accurate and effective calculation of all EEPs is then required, while taking into account the mutual coupling. In this section, the quality of EEPs obtained using HARP is evaluated. A small testbed of 16 SKALA elements is first analyzed, where the EEPs of all elements are calculated at all frequencies with 1 MHz frequency steps. Then, the SKA station in Fig. 2 is simulated at four different frequencies, i.e. 50 MHz, 110 MHz, 200 MHz and 350 MHz, using HARP and WIPL-D, and the calculated EEPs are compared.

III-B1 EEPs of A Small Array

The performance of HARP has been tested against the brute force MoM solution of a small array of 16 SKALA elements for the frequency band (50 – 450 MHz). The array corresponds to a testbed for SKA1-Low technology and is situated at the Mullard Radio Astronomy Observatory (MRAO) near Cambridge, UK [8]. This array with only 16 elements is in a pseudo-random configuration, such as that of SKA stations, and has the same expected element density, with an average inter-element spacing of half wavelength at 77 MHz. The error between the patterns is calculated as:

[TABLE]

where $E$ and $E^{ref}$ being the EEP calculated by HARP and brute-force MoM with the electric fields being normalized by the maximum of the isolated antenna pattern. This error is calculated in all directions for all involved antennas. The condensed result is shown in Fig. 9 by calculating the percentiles of the error and plotting them over frequency. The HARP method shows very good accordance to the brute-force MoM results over the whole frequency range. The error increases at the highest end of the band, but it is still below -30 dB.

III-B2 EEPs of SKA station

The performance of HARP on the SKA station is now validated with patterns simulated using the WIPL-D software [34]. The EEP of an arbitrary element (red dot in Fig. 2) obtained using HARP is taken as an example. Fig. 10 shows the EEPs at four different frequencies and their differences w.r.t. WIPL-D, where the error patterns are computed using (8) with $E$ and $E^{ref}$ being patterns obtained using HARP and WIPL-D, respectively. An excellent agreement is shown at low-frequencies, where the difference level is found below -30 dB at 50 MHz. The variation becomes more visible as the frequency increases. A detailed study of the performance of different software at 350 MHz is carried out on an array of 2 antennas positioned along the $x-$axis with varying distance. The EEPs obtained using CST, WIPL-D, brute-force MoM and HARP are plotted together in Fig. 11 for a separation equal to 8 m. For this 2 element array, the EEPs are similar to the isolated pattern, and mutual coupling clearly appears as ripples in the main beam. The results show that while the mutual coupling is accounted for in all software packages, a small difference is observed in the EEPs. The reason is attributed to the difference in the geometrical modeling of the SKALA, i.e. the use of different meshes in each software, which leads to differences in the isolated pattern. At the highest frequency, a small difference in antenna modeling will result in a significant variation in the radiation patterns. This explains the discrepancies in the EEPs visible in Fig. 10(h), as well as the self-impedance behavior in Fig. 5. A stronger indicator that the discrepancies mainly come from the mesh is that the errors in the isolated pattern are almost as large as in the EEPs. Moreover, it is important to note that, in Fig. 11, patterns calculated using brute-force MoM and HARP overlap very well, which confirms the accuracy of HARP method.

III-C Array Patterns

In this section, the EEPs calculated using HARP and WIPL-D software are exploited to form the array pattern. A simple beamforming method is implemented with uniform excitation and scanning at broadside. Fig. 12 shows the radiation pattern of the SKA station at different frequencies in two principal planes calculated using HARP and WIPL-D. The two methods show a great agreement across the frequency band, as the difference in patterns, evaluated by (8), is only shown at a -30 dB level, except for the highest frequency, i.e. 350 MHz, where the error in the main beam is near -24 dB. As explained above, the visible difference in the patterns at higher frequency is due to the difference in modeling of the SKALA in the two software. To demonstrate the impact of the isolated element modeling, we show that it can be partially undone by properly adding the difference of isolated patterns between HARP and WIPL-D into the EEPs, as shown in Fig. 13. A better agreement is clearly seen w.r.t. Fig. 12(h), especially near the main beam. This exercise illustrates that the difference between WIPL-D and HARP at highest frequency is mainly due to the modeling of the SKA antennas in these two software packages. Expectedly, this fine modeling difference becomes more visible at high-frequency.

Finally, the performance of WIPL-D, CST and HARP for the analysis of a SKA station at 110 MHz is detailed in Table I. The CST simulations were carried out on a server with two processors and two cores each with 384 GB of RAM at the University of Cambridge, UK, whereas the WIPL-D ones were performed on a Windows Server 2008 R2 Enterprise with an Intel Xeon X7550 processor and 125 GB RAM at ASTRON, The Netherlands. All the simulations for HARP were carried out on a desktop computer with 16 GB RAM and Processor Intel Core $i7-3770$ CPU 3.4 GHz at UCL. For HARP, the accounted time stands for matrix filling, solving the system of equations and calculation of all EEPs once the HARP model is built, as reported in Table II. A great time saving factor is achieved using HARP, while providing highly accurate results, as demonstrated in Fig. 12. The establishment of the HARP model takes about 3 hours per frequency. It is important to recall that this preparation phase takes place only once and for all, regardless of the array configuration, and can be re-used to analyze any array made of the same element. This feature is very beneficial for simulating SKA arrays, where hundreds of stations are planned, as well as to analyze extremely large arrays, as will be discussed in the next section.

IV Simulation of Extremely Large Arrays

This section shows the performance of the sparse matrix and tessellation techniques on the analysis of very large arrays. The impact of the Radius of Influence (RI) on the array patterns is studied on arrays of different sizes containing 256, 500, 800, 1000 and 1500 elements. The latter has a diameter of 85 meters, as shown in Fig. 14, and the same array density has been considered for all the other arrays. The arrays are simulated at 110 MHz using the HARP approach; EEPs and the array pattern are calculated. It is noted that the array of 1500 elements cannot be handled using CST nor WIPL-D with our current hardware.

IV-A Sparse Matrix Approach

Fig. 15 shows the exact patterns and the error for different RIs in two main planes for the 1500 elements array. The errors are evaluated using (8), where $\overrightarrow{E}(\theta,\phi)$ and $\overrightarrow{E}^{\textrm{ref}}(\theta,\phi)$ stand for the exact pattern and those obtained using the sparse array approach, respectively. It is clear that for larger values of the RI, the errors are smaller, since more elements are included in the RI region. It is also seen that, using only the isolated pattern, the error level is higher than the sidelobe level at far-out angles. Even for RI = 30 m, the error is close to some sidelobes. The results again show the strong effect of mutual coupling in a random array made of SKALA elements. Fig. 16 displays the convergence of error curves versus the RI for different array sizes at 110 MHz. It shows the trend of reducing errors for larger RI and one can observe that to achieve a given error level, larger arrays require a smaller RI (normalized to array size). This is probably due to the slow averaging effect of randomization of the array on mutual coupling. As the RI is smaller (w.r.t. array size) for large arrays, the interaction matrix will be sparser, which implies the capability to analyze large arrays using the RI concept. The computational complexity of the approach is detailed in Table III, where $N_{I}$ is the average number of elements within the RI, which is proportional to $\rm{RI}^{2}$.

IV-B Tessellation Approach

Fig. 17 shows the radiation patterns and errors using the tessellation approach on the array of 1500 SKALA elements shown in Fig. 14 for $\rm{R_{in}}=10\,m$ and different $\rm{R_{ext}}$ values. It is interesting to see that tessellation is able to provide high-quality solutions w.r.t. the exact ones as the error levels are quite low even for small $\rm{R_{ext}}$. Despite sharing the same RI idea as the sparse array approach, the tessellation seems to perform better. The reason lies in the fact that even if the sparse array approach only accounts for the interaction between elements inside the RI, the solution of the resulting system of equations still provides non-zero currents for elements outside of the RI region, while the excited elements remain within the RI. This “artificial” current contributes to the radiation patterns, while probably not being highly relevant nor accurate. Conversely, the tessellation approach only solves for the subarray and calculates the contribution of those elements inside the subarray.

Fig. 18 shows the convergence of the tessellation approach for different array sizes. As $\rm{R_{ext}}$ increases, the error decreases almost linearly in logarithmic scale. This allows one to obtain a given accuracy by finding a proper $\rm{R_{ext}}$. Increasing $\rm{R_{ext}}$, however, results in increasing the computational complexity, both for matrix filling and for solving the system of equations. An optimum choice of $\rm{R_{ext}}$ can be found following Eq. (7), as discussed in Section II-D; it favors an extended region comprising half of all neighboring inner-hexagons. A detailed complexity of the tessellation approach is listed in Table III, where $N_{s}$ is the average number of elements in a subarray, which is proportional to $\rm{(R_{in}+R_{ext})^{2}}$ with $\rm{R_{in}=R_{ext}}$. By selecting $N_{I}$ and $N_{s}$, i.e. choosing $\rm{RI}$ and $\rm{R_{in}}$, respectively, one can obtain a same complexity for both methods. For extremely large arrays, $N_{I}$ and $N_{s}$ are at least one order of magnitude smaller than $N_{a}$, which results in a big saving factor compared to the direct solution over the whole array.

It is also interesting to see that the convergence curve in Fig. 18 is quite independent from the array size, which assures the performance of the tessellation approach for other (larger) arrays, and allows one to select $\rm{R_{in}}$ and $\rm{R_{ext}}$ in advance. Regarding the analysis, the tessellation approach offers the possibility of parallelizing the computation as the analysis of each subarray is totally independent. Moreover, the tessellation is not limited by the array size as each subarray can be solved effectively with a present-day desktop computer, providing enough temporary memory RAM to store the reduced MoM matrix and solutions, for example 16 GB has been proven sufficient in all tested cases. For the largest array, i.e. consisting of 1500 SKALA antennas, the analysis took about 15 minutes using the tessellation approach to simulate the array and calculate all the EEPs.

V Conclusion

A fast full-wave technique for the analysis of large arrays made of complex 3-D metallic antennas has been presented. The calculation of interactions between antennas has been accelerated using the HARP technique, where a non-expensive model is built for MBF interactions. The reduced matrix is then rapidly filled using such model. Several improvements have been added to HARP to cope with the complexity of solving very large arrays. An experimental example of compact array and simulations of SKALA arrays have demonstrated the performance of HARP. Furthermore, a tessellation approach has been proposed for analyzing very large arrays of SKALA elements, offering a solution to the analysis of extremely large arrays, such as the SKA telescope, for which hundreds of stations are foreseen. In particular, the proposed approach enables the analysis of very large arrays on present-day desktop computers, while providing highly accurate results.

The simulation technique opens the possibility to integrate the full-wave simulation in an optimization procedure [40], which is capable of taking into account mutual coupling to optimize the performance of the arrays (e.g. in terms of sensitivity, including noise coupling). Moreover, the accurate EEPs obtained using HARP can be incorporated in the calibration stage for modern radars and astronomical instruments such as the SKA. The use of the PMCHWT approach [41], combined with the surface equivalence principle [42], allows the extension of HARP to antennas comprising dielectric material.

Acknowledgment

The authors would like to thank the UK’s STFC funding that has supported this work. They are also grateful to their SKA colleagues for numerous discussions about SKA station modeling.