*Article* **Numerical Simulation and Experimental Validation of an Oil Free Scroll Compressor**

#### **Massimo Cardone \* and Bonaventura Gargiulo**

Department of Chemical, Materials and Production Engineering, University of Naples Federico II, 80125 Naples, Italy; bonaventura.gargiulo@unina.it

**\*** Correspondence: massimo.cardone@unina.it; Tel.: +39-081-7683675; Fax: +39-081-2394165

Received: 15 October 2020; Accepted: 6 November 2020; Published: 10 November 2020

**Abstract:** This paper presents a virtual model of a scroll compressor developed on the one-dimensional analysis software Simcenter Amesim®. The model is semi-empirical: it needs some physical details of the modelled machine (e.g., the cubic capacity), but, on the other hand, it does not require the geometrical features of the spirals, so it needs experimental data to calibrate it. The model also requires rotational speed and the outlet temperature as boundary conditions. The model predicts the power consumption and the mass flow rate and considers leakages and mechanical losses. After the model presentation, this paper describes the test bench and the obtained data used to calibrate and validate the model. At last, the calibration process is described, and the results are discussed. The calculated values fit the experimental data also in extrapolation, despite the model is simple and performs calculations within 7 s. Due to these characteristics, the model is suitable for being used in a larger model as a sub-component.

**Keywords:** scroll-compressor; experimental validation; numerical model

#### **1. Introduction**

Scroll compressors are widely used in applications where noise pollution and low vibrations are relevant factors, such as domestic refrigeration and domestic climatic control. Due to its unique properties, these machines are given much attention by scientific and industrial researchers. Some works, especially the oldest ones [1,2], are mainly focused on the theoretical functioning of those machines, even though the final aim has always been the improvement of the scroll efficiency. The scroll compressors efficiency-enhancing is pursued using different approaches, such as researching on the design of rotor profiles to reduce volumetric losses [3,4] or investigating the cooling effects on the compression work [5,6].

The latest most attractive research branches on scroll machines concern studies on the machine behaviour into thermodynamics cycles [7], works on the performance of injecting water (or vapour) compressors [8,9] and studies on the scroll expanders [10].

Many works on scroll compressors or expanders use a mathematical model that is generally virtualised with a low-level programming language. Among them, some use a geometrical approach [11–15] while others use semi-empirical methods [16–18]. Some details of one work for each group are briefly illustrated below. Blunier et al. [14] presented a model written in very-high-speed-integrated-circuits hardware description language (VHDL) code including the scroll's geometrical features in it; the model does not need any calibration. Winandy et al. [18] used the Engineering Equation Solver (EES) software to virtualise the model equations. The model needs seven parameters to be obtained through a calibration process on mass flow rate, power output and outlet temperature. The required data points are obtained by themselves through several experimental runs.

The use of commercial 0-1D fluid-dynamic code was not frequent for machines models, while they are often used to simulate a whole thermodynamic cycle [18–20]. Ziliani et al. [19] used the commercial software Amesim® ( Siemens PLM Software, Plano, TX, USA)to model an entire Organic Rankine Cycle (ORC) plant where it is supposed to be used a screw or a scroll expander. The machines' behaviour was deduced by a Computational Fluid Dynamics (CFD) simulation (for the scroll compressor) and a geometry-based simulation (for the screw compressor). Bracco et al. [20] used a similar approach for their ORC plant simulation: their machine model needs four parameters, that are based on a combination of some manufacturer information, experimental data and a self-made scroll simulation tool.

However, Bell et al. [21] recently (2020) developed an open-source platform (named PDSim) specific for modelling positive-displacement compressors. Tanveer et al. [22] compared different software on a reciprocating compressor finding that the PDSim suit has good potential. On the other hand, Rak et al. [23] relied on a CFD analysis to model in detail a cooled scroll compressor, thus considering heat transfer issues.

The main purpose of this work was the development of a fluid-dynamic model of a scroll compressor and its experimental validation. The model does not consider the internal geometry of the scroll, nor its kinematic behaviour. The model aim is to perform calculations in a few seconds so that it can be used as a sub-component of a whole plant model. In this paper, an oil-less commercial scroll compressor is modelled by a zero-dimensional semi-empirical model developed using the commercial software Simcenter Amesim® (version 19.2). The compressor is viewed as a pneumatic system made up by a series of peculiar reciprocating compressors with driven valves. Leakages and a keyed fan power consumption are considered. The experimental activities are performed on a commercial oil-less scroll compressor, at four rotational speed levels and six pressure levels.

#### **2. The Scroll Compressor**

The same timing scheme characterises all rotary machines. The rotors uncover ports and intercept cells carved in the stator, carrying the working fluid from the inlet side to the outlet. As the cells have a decreasing volume, a design pressure ratio β*<sup>i</sup>* is generated by the internal volume reduction (for roots compressors β*<sup>i</sup>* = 1). The compression ratio β, required by the application, may overlap or not with the compression ratio β*<sup>i</sup>* , that the machine can produce due to volume change of cells. Usually, regarding the scroll compressors, there is a moving spiral and a fixed spiral. The moving spiral describes an orbit around the base circle centre of the fixed. The spirals are thick circle involutes, and, in most cases, they are equals. Many authors described the scroll geometry and kinematic characteristics in details (e.g., Chen et al. [13]). The proposed model is not based on the geometry of the spirals, so this work will not examine these issues.

Figure 1a shows a scroll compressor working scheme, while Figure 1b shows a rotary volumetric compressor's ideal cycle diagram. Referring to Figure 1, these machines make a first internal compression (from state 0 to state 1) through volume reduction. Then, at the opening of the last internal contacts of the machine, two volumes are put into communication: the last cell filled with gas in state 1 (*p*1, *T*1) and the discharge volume filled with gas in state 2 (*p*2, *T*2). Therefore, an instantaneous mixing phase (at constant volume) starts: the state of the whole gas becomes an intermediate state (*px*, *Tx*). In the next step of the machine rotors, while the volume of the last cell decreases, state 2 in the volume *V<sup>m</sup>* is restored.

An authors' previous work [24] theoretically analysed the machine β ≥ β*<sup>i</sup>* field. It illustrated an alternative representation of a generic ideal rotary compressor's working scheme (see Figure 1). The theoretical analysis confirmed that the two theoretical models are equivalent. It was demonstrated that the proposed representation leads to a simpler but rigorous equation to calculate the ideal specific work consumption (see Equation (1)).

**Figure 1.** (**a**) Scroll functioning scheme; (**b**) rotary compressors' diagram.

In this work, Equation (1) is used to give a theoretical validation to the proposed model, before implementing real losses and calibrating it on experimental data.

$$w\_{0-2} = \frac{k}{k-1} \text{RT}\_0 \left(\boldsymbol{\beta}^{\frac{k-1}{k}}\_i - 1\right) + \frac{k}{k-1} \text{RT}\_X \left(\boldsymbol{\beta}^{\frac{k-1}{k}}\_{X2} - 1\right) \left(\frac{m\_m + m\_\mathbb{C}}{m\_\mathbb{C}}\right) \tag{1}$$
 
$$\text{contronic index } \mathcal{V} \text{ is the grid constant (for air) } T\_\Gamma \text{ is the inlet temperature } T\_\Gamma \text{ is the}$$

*β*

*β β*

*β*

where *k* is the isentropic index, *R* is the gas constant (for air), *T*<sup>0</sup> is the inlet temperature, *T<sup>X</sup>* is the temperature in the state *X*, *m<sup>c</sup>* is the mass elaborated per cycle, *m<sup>m</sup>* is the mass in the discharge volume. The other abbreviations are collected in Equation (2).

$$\beta\_i = \frac{p\_1}{p\_0} \quad ; \quad \beta\_{X2} = \frac{p\_2}{p\_X} \; ; \; T\_X = \frac{p\_1 V\_1 + p\_2 V\_m}{\frac{p\_1 V\_1}{T\_1} + \frac{p\_2 V\_m}{T\_2}} \tag{2}$$

#### **3. Numerical Model**

The numerical model is developed on Simcenter Amesim® software.

Referring to Figure 2, two pairs (A, B) of variable volumetric chambers are used to simulate the scroll compressor. The volumetric variation does not follow the actual scroll chambers variation. The model uses a simple sinusoidal function of the shaft rotation. The compression process is modelled through the subsequent steps: the suction phase of chambers A1 and B1 stands for the scroll suction. They are in phase opposition so that the system has the suction phase 360◦ long and the sum of their cubic capacity is the scroll capacity. Three controlled valves for each pair of chambers are used to simulate the openings of the scroll contacts. These valves are opened every 180◦ . Focus on one pair of chambers (e.g., chambers A1 and A2). The first valve stands for the contacts that enclose half scroll suction chamber when it reaches its maximum capacity; the second valve separates the first chamber from the second and it stands for the contacts that separate the scroll suction chambers from the other scroll chambers. Then, the third valve stands for the contacts that separate the scroll compression chamber from the discharge chamber. Chambers A2 and B2 are smaller than the others and they are in phase opposition respective to chamber A1 and B1. When the second valve opens, the compression phase begins due to the total volumetric decreasing achieved after the second valve opening.

**Figure 2.** Model code in Amesim environment. **Figure 2.** Model code in Amesim environment.

An ideal electric engine imposes a constant speed to the whole model. In this work, electric losses are neglected, because the experimental work input is measured downstream of the electric motor.

After the compression phase, the third valve opens and the second close. The discharge phase and the suction phase start respectively in chambers 2 and 1. The other pair of chambers (e.g., B1 and B2) are just the same, but they are in opposition of phase. Therefore, there is an entire suction and discharge that lasts 2π radians long in a revolution. An ideal infinite volume simulates the compressor discharge tank. π

In this model, the heat losses are neglected, so the chambers are adiabatic. Therefore, the outlet temperature, as well as the inlet temperature, are set based on the experimental data. According to Yu Chen et al. [13] work there are two types of leakages: the flank leakages and the radial leakages. The flank leakages are simulated by an imperfect closing of the valves (the valves V1, V2, V3 that simulate the internal contacts). Therefore, when they should be closed, they are slightly opened. The radial leakages are simulated by internal by-pass (R). The mass flow rate through the valves is calculated via Equation (3) [25].

$$
\dot{m} = \underset{A}{A} \mathbf{C}\_{q} \mathbf{C}\_{m} \frac{P\_{up}}{\sqrt{T\_{up}}} \tag{3}
$$

where *A* is the actual orifice area, *C<sup>q</sup>* is Perry's coefficient from his correlation [26], *C<sup>m</sup>* is the flow coefficient as also described by Szente et al. [25], *Pup* and *Tup* are the upstream pressure and temperature. ඥܶ௨

The mechanical losses (ML) are modelled as a fixed frictional loss (one constant to be estimated). The keyed fan (KF) is simulated by a torque load function of the square of the rotational speed, so there are three parameters (*A, B, C*) to be estimated (see also Equation (4)).

$$Torque\ load = \overline{A}\omega^2 + B\omega + \mathbb{C}\tag{4}$$

At first, the model is tested on an ideal case, suppressing the elements KF, ML and R. The model results are compared with the theoretical results of Equation (1), at different internal and total pressure ratios. Figure 3 shows that there are no significant deviations between the scroll ideal performances and the simplified numerical model. The maximum deviation obtained is always under 0.8% of the theoretical values (the model calculation is always greater than the theoretical value). These differences are caused by some small pressure losses still present in the model, and, secondary, by the model outlet tank, that the code treats as an infinite volume (the theoretical calculation is under the assumption of an outlet tank of 10<sup>5</sup> times the scroll cubic capacity). Overall, the proposed model is congruent with the theory, so, enabling KF, ML and R, it should be able to simulate real scrolls. The values of these elements' parameters are estimated through a calibration process based on experimental data.

**Figure 3.** Comparison between theory and the proposed model. **Figure 3.** Comparison between theory and the proposed model.

#### **4. Results**

#### *4.1. Experimental Activities*

#### 4.1.1. Experimental Setup

Figure 4a shows the experimental test bench used, Figure 4b illustrates its scheme and Table 1 collects the details on its main components. The machine tested is a 2.1 kW oil-free scroll compressor (C) taken from the ATLAS COPCO SF2 clean air generator (ATLAS COPCO, Nacka, Sweden). The required power is given by a three-phase oscillating-casing electrical asynchronous motor (EM). The motor is connected to an inverter that controls the motor speed modifying the current frequency. A load cell measures the required torque (Tq): the motor oscillating-casing is constrained by the load cell through an arm of a known length, so the product of the sensed force and the arm length is equal to the torque given by motor. A trapezoidal belt connects the compressor to the motor with a unitary transmission ratio (TB). The encoder (RPM) is integral with the compressor axis, measuring the compressor speed. The compressor original cooling fan is a centrifugal fan (KF), and it is keyed on the scroll axis. A turbine flow meter (V) measures the air volumetric flow rate at the suction of the compressor. Two K-type thermocouples (T) (mounted through a T-joint) measure the temperatures, both at the inlet and outlet pipe of the machine. The outlet thermocouple is mounted at 15 cm from the outlet port, due to the compressor built-in external fins. A piezo-resistive sensor (p) is linked to the calm reservoir (TANK) at the outlet to measure the required average pressure imposed on the compressor. The circuit ends with a regulator valve (RV) to control the outlet pressure. The data are digitalized by data acquisition system made up by a NI-DAQ USB 6259 (National Instrument, Austin, TX, USA) and a NI-FieldPoint cFP 1808 coupled with a cFP-TC-120 module that provides built-in cold joint correction for the thermocouples. All data are processed by a self-made software realized in LabVIEW ™ (version 15.5, National Instrument, Austin, TX, USA) code.

**Figure 4.** (**a**) Test bench; (**b**) test bench scheme.


**Table 1.** Test bench elements.

#### 4.1.2. Experimental Plan

The test campaign considered in this work consists of four runs, each performed at total pressure ratios (β) from 2 to 7 and then from 7 to 2. Each run is performed at a different constant speed: 1000, 1500, 2000 and 2500 RPM. Once the compressor speed is set, the pressure in the outlet tank is controlled through the regulation valve. Starting from the lowest pressure ratio, when the desired value is reached, the acquisition system waits until the outlet temperature is stationary and then saves the data. Then, the regulation valve is tightened to achieve the next desired pressure ratio. The process is repeated until the pressure ratio reaches 7, and the system saves the data. Then, compressor state is modified reaching a pressure ratio of 7.5 (approaching the limits of the test bench). After 5 min, the process is repeated from pressure ratio 7 to 2. Despite the different thermal dynamics between the rising and the falling part of the run, the next paragraphs show that the measured data are close to each other. *β*

#### 4.1.3. Experimental Results

Figure 5 presents the working fluid (air) temperature versus the total pressure ratio and for various rotational speed. The experimental reproducibility is average: the higher error is below 5%. It is shown that the temperature is an increasing function of both compression ratio and rotational speed. The temperature is still increasing due to both the rising of the compressed mass flow rate and the rise of frictional losses, despite a higher rotational speed causes a higher cooling flow rate (generated by the keyed fan).

**Figure 5.** Experimental data: outlet temperature.

The mass flow rate and power consumption experimental data are shown together with the simulated data (the model results will be discussed in the next paragraphs): Figure 6 shows the compressor mass flow rate and Figure 7 shows the power consumption with the fan. According to Figure 6, the higher the pressure in the delivery tank, the lower the mass flow rate. This effect is heavier for low rotation speed. This tendency is caused by both the temperature effect (the whole machine temperature is higher at higher compression ratios) and the leakages (a higher-pressure gradient across the gaps leads to more air leaked). The experimental reproducibility is good: the higher error is below 3%, and the absolute errors are under 0.1 g/s.

**Figure 6.** Model validation: mass flow rate.

**Figure 7.** Complete model validation: power consumption.

Figure 7 shows the power consumption data of compressor for all experimental tests, calculated via Equation (5), where *Mexp* and ω*exp* are the measured torque and the measured rotational speed. The experimental reproducibility is average: the higher error is below 5%

$$P\_{\exp} = M\_{\exp} \omega\_{\exp} \tag{5}$$

*β β*

*β β*

*β*

*β*

Both the power and the temperature are increasing functions of both compression ratio and rotational speed.

#### *4.2. Model Results*

#### 4.2.1. Model Calibration and Validation: Leakages

In the model, the adiabatic process simplified hypothesis is done, supposing that the fan cooling affects only the final stage of the compression process (i.e., the curve X-2 shown in Figure 1). The first step of the calibration process is to determine the leakages. In the model (see Figure 3), the flank gaps are modelled as an imperfect closing of the internal contacts (the valves V1–V3) while the radial leakages are modelled as a by-pass. The model requires the orifice area for each of them. These areas are not the physical areas of the scroll gaps, but their model representation. It is possible to determine them using the experimental data on the mass flow rate. It can be postulated that both the centrifugal force and the increasing temperature would reduce the internal clearances. Therefore, the influence of the rotor speed over the leakages must be considered.

The calibration is performed using experimental data at β = 2 and β = 6 at 1000, 2000 and 2500 RPM. All the other experimental points are used to validate the model. In particular, the β = 7 points are used to verify the model consistency outside the calibration data field, while the 1500 RPM experimental points are used to evaluate the model consistency in simulating the scroll running at another rotational speed. The calibration points temperatures are set based on the experimental data. In other cases, the temperatures are based on a linear regression of the values measured at β = 2 and β = 6. As for the 1500 RPM temperatures, they are based on a linear regression of the values determined at 1000, 2000 and 2500 RPM. Similarly, the orifice areas are estimated anew for each rotor speed. For the validation at 1500 RPM, the orifice areas are calculated based on a linear regression of the values determined at 1000, 2000 and 2500 RPM. Figure 6 shows the model compliance with the real scroll in terms of mass flow rate after the calibration process (coloured dots for the experimental data and coloured lines the calculated one).

#### 4.2.2. Model Calibration and Validation: Mechanical Losses

A new set of experimental runs were performed to calculate the mechanical losses (modelled by the ML block). This experimental set-up is characterised by the absence of the keyed fan. In this way, we get the difference between the experimental power consumption and the simulated one (once the KF block is suppressed, see Equation (6)).

$$\mathcal{M}\_{\text{ML\\_block}} = \frac{\mathcal{P}\_{\text{experimental}}}{\omega} - \mathcal{M}\_{\text{model\\_scroll}} \tag{6}$$

The machine temperatures were higher than the previous case, so the higher pressure ratio tested is 4. As the thermal conditions were altered, the spirals thermal deformation is different; therefore, a new leakages calibration was needed. Figure 8 shows the accordance between the model and the new experimental set (the one without the fan) after the new calibration. It is possible to determine the value of the ML torque confronting the power consumption experimental data with the model output. The torque can be assumed as a constant independent by the rotational speed as the mechanical losses are just simple friction losses (there are no other auxiliaries nor inertial forces).

#### 4.2.3. Model Calibration and Validation: Keyed Fan

The speed effect on the torque losses (ML) was neglected, so the mechanical loss is set to constant. As previously said, the keyed fan (KF) power consumption is modelled by a torque load function of the square of the rotational speed, so there are three parameters (A–C) to be estimated (see also Equation (2)). Once these parameters are determined, they are set as constant, independent from both the speed and the pressure ratio. As previously, the data at β = 2 and β = 6 at 1000, 2000 and 2500 RPM are used to calibrate, the others are used to validate. Figure 7 shows the comparison between the experimental power consumption (coloured dots) and the calculated one (coloured lines), after the calibration process. The figure shows the power consumption as a function of the pressure ratio and

*β*

the rotational speed. In addition, the 1500 RPM and all β = 7 experimental data are well predicted by the model. The model complies with experimental data, but some deviations remain. These deviations are higher for lower rotational speed (that are further from the real compressor nominal speed and further from adiabatic behaviour).

− ܯௗ ௦

ܲ௫௧ ߱

= ெܯ

*β*

**Figure 8.** Comparison of mass flow rates and power consumption between the model and the real compressor without the keyed fan.

#### **5. Conclusions**

*β β β* In the first part of the paper, the key features of a scroll compressor are presented. Then a numerical model developed in Simcenter Amesim® is presented: the model includes leakages and mechanical losses. Despite that the integrated keyed cooling fan power consumption is considered, its cooling effect is neglected, and the compressor is considered adiabatic. To partially overcome this limitation, the outlet temperature is set as a boundary condition using experimentally derived values. The model was preliminary tested on ideal cases (no losses or leakages). Then, the experimental activity is presented: a series of experimental runs were performed on a commercial scroll compressor. Four levels of rotational speed (1000, 1500, 2000, 2500 RPM) and six levels of the pressure ratio (2–7) were considered. The experimental data are reproducible as the errors are under 5%. Finally, the model is calibrated on the experimental data. Six points are used to calibrate the model (at β = 2 and β = 6 at 1000, 2000 and 2500 RPM), while the other 26 are used to validate. After the proper calibration of the leakages, mechanical losses and fan power consumption, the model can follow the scroll's real behaviour.

The model does not consider the actual scroll geometry and kinematic. This feature could be both an advantage or a limitation: it is a limitation because it is not possible to calculate any inner quantity nor instantaneous quantities (e.g., instantaneous torque); it is an advantage when the geometrical features are unknown. Moreover, the calculations last less than seven seconds. Overall, the model is not suitable as a machine design aid, but, on the other hand, it can be used as a component for a whole plant simulation.

**Author Contributions:** The conceptualization, methodology and experimental design have been developed by M.C.; the data acquisition, software model development and model calibration have been performed by B.G. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Acknowledgments:** The authors would like to thank Salvatore De Cristofaro for his technical assistance in laboratory activities.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Abbreviations**


#### **References**


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **A Parametric Study of Wave Energy Converter Layouts in Real Wave Models**

#### **Erfan Amini <sup>1</sup> , Danial Golbaz <sup>1</sup> , Fereidoun Amini <sup>2</sup> , Meysam Majidi Nezhad <sup>3</sup> , Mehdi Neshat <sup>4</sup> and Davide Astiaso Garcia 5,\***


Received: 2 October 2020; Accepted: 18 November 2020; Published: 20 November 2020

**Abstract:** Ocean wave energy is a broadly accessible renewable energy source; however, it is not fully developed. Further studies on wave energy converter (WEC) technologies are required in order to achieve more commercial developments. In this study, four CETO6 spherical WEC arrangements have been investigated, in which a fully submerged spherical converter is modelled. The numerical model is applied using linear potential theory, frequency-domain analysis, and irregular wave scenario. We investigate a parametric study of the distance influence between WECs and the effect of rotation regarding significant wave direction in each arrangement compared to the pre-defined layout. Moreover, we perform a numerical landscape analysis using a grid search technique to validate the best-found power output of the layout in real wave models of four locations on the southern Australian coast. The results specify the prominent role of the distance between WECs, along with the relative angle of the layout to dominant wave direction, in harnessing more power from the waves. Furthermore, it is observed that a rise in the number of WECs contributed to an increase in the optimum distance between converters. Consequently, the maximum exploited power from each buoy array has been found, indicating the optimum values of the distance between buoys in different real wave scenarios and the relative angle of the designed layout with respect to the dominant in-site wave direction.

**Keywords:** layout assessment; wave energy conversion; renewable energy; real wave model

#### **1. Introduction**

Wave energy is expected to contribute towards the development of a carbon-free electricity generation. The theoretical computation of wave energy potential over the oceans is projected to be in the order of 1–10 TW [1], which can cover the current global energy demand [2]. This tremendous potential has attracted attention from research societies, which have proved that harnessing electric power from ocean waves is possible [3,4]. Wave energy converters (WEC) are planned to be stationed in an

array constituted of many converters -similar to offshore wind turbines. The initial developments of the analytical modeling of hydrodynamic forces on submerged buoys can be found in [5]. However, it has been further developed since then. The next stage of studies was focused on enhancing the design and power take off system of a single buoy [6]. The next studies bring the idea of WEC's array by conducting a comparative study on different configurations [7]. The proceeding research phase concentrated on finding the optimal value for WEC's array parameters (such as optimal position or layout) using either numerical, parametric or optimisation-based solutions [8–11], as this story falls in this category. The position of converters in the array which is scattered through the array has a direct relationship with the performance of the array because hydrodynamic interactions between them can be constructive or destructive. These interactions depend on the configuration of the array. Consequently, this is the main reason to investigate these interactions in order to apply them to reinforce the total power output. There are many relevant publications with this subject by several R&D units across Europe in the past by the pioneering works [12–16], and it is still an interesting research field, as several investigations have been published recently [17–20] . Furthermore, the identification of techno-economically feasible decarbonisation paths and sustainability transitions have been investigated by [21–23]. Some of the related research projects that considered the performance of arrays or converters' distance were undertaken by [9,24–26] and the effects of nonlinear mooring forces via a time-domain analysis and the influence of interactions between WECs are well described in [27,28], respectively. Table 1 demonstrates a brief survey of some of the recent literature on the various aspects of WECs including layouts, PTO and design optimisation. Some of the mentioned research has used hindcast wave models; however, different layout configurations were considered regarding real wave scenarios in this study.


**Table 1.** A briefly survey some of the recent literature on the layout, Power Take-Off (PTO) parameters and design optimisation of wave energy converters.


**Table 1.** *Cont.*

The CETO6 is a fully-submerged point absorber wave energy converter that is manufactured, installed, and updated by Carnegie Clean Energy Ltd. The prospective location of the WECs array is off the coast of Albany due to its exposure to open ocean wave conditions [58]. This study has been conducted based on the numerical simulation of this converter's array. Our concentration is particularly on the arrangement optimisation of WEC arrays and shows the effectiveness of the inter-distance among WECs to produce more power. In order to establish an array of WECs, an optimal layout is chosen to maximise the power conversion; however, the number of WECs is a significant factor. We evaluate various numbers of WECs as an array, arrangements and separations, and report the performance of the layouts using q-factor, power of each converter and total power output. The distances between the WECs, and the array size are constrained, which is a more realistic approach for studying WEC arrays. Finally, a landscape numerical analysis is performed with regard to evaluating the position effect of each WEC in the array's power output using a grid search approach.

It should be noticed that such research has not been investigated in the mentioned real wave scenarios (Perth, Adelaide, Sydney, and Tasmania) regarding this parametric study. Therefore, the main motivation of this study is to evaluate the output performance of the simulated CETO6 arrays to find a suitable layout with optimal distance, and the rotation angle to the dominant wave direction in these specific case studies. Moreover, the investigation coverage is more comprehensive than in other research studies by exploiting wave power using a ten-degree resolution covering the whole area of study, compared to, Bozzi et al. [24] presented the implementation and evaluation of a few numbers of WEC separation distances (5, 10, 20 and 30 buoy diameters) and incident wave directions (30◦ apart).

This paper is structured into five sections. Section 2 presents a brief description of the hydrodynamic WEC array interaction model, modeling the wave climate and the equations used to compute the produced power. Section 3 expresses the layout assessment routine and presents the strategy to explore the optimal position of the WECs in the array. Section 4 discusses the array layout investigation results in terms of performance and optimal array layout solutions. Subsequently, Section 5 summarises the principal finding of the paper.

#### **2. Numerical Modelling**

#### *2.1. Wave Energy Converter*

In this study, a CETO6 wave energy converter with a three-tethered mooring system is considered which has a fully submerged spherical buoy attached to the seabed by the tethers, as shown in Figure 1. This model is developed in MATLAB and was modified in 2020 [59]. The WEC details are:

buoy radius = 5 (m), submergence depth = 3 (m), water depth = 50 (m), buoy mass = 376 (t), buoy volume = 523.6 (m<sup>3</sup> ), tether angle = 55 (degree), PTO stiffness = 2.7 <sup>×</sup> <sup>10</sup><sup>5</sup> (N/m), PTO damping= 1.3 <sup>×</sup> <sup>10</sup><sup>5</sup> (Ns/m).

**Figure 1.** Schematic representation of the CETO6 modelled point absorber wave energy converter (adapted from [60]).

This buoy, which is floating at sea, moves in six degrees of motion. However, due to the converter's spherical shape, its displacement is in three degrees of freedom which are surge, heave, and sway. Based on these degrees, the motion equation can be written on the frequency domain.

$$\begin{aligned} \Sigma F &= m\Xi, \\ &= F\_m + F\_{hs} + \mathcal{W} + F\_R + F\_{PTO} + F\_{W\_k} + F\_{VD} \end{aligned} \tag{1}$$

where *F<sup>m</sup>* is the mooring force, *Fhs* is the hydro-static force resulting from buoyancy, *W* is the body weight, *F<sup>R</sup>* represents added mass and wave damping forces, force resulted by PTO system is *FPTO*, *FW<sup>k</sup>* represents the vertical components of the wave exciting force and *FVD* is the vertical viscous drag force [61]. This equation is used in order to describe a time-domain response of the WECs in waves, and can be rewritten as:

$$\dot{z}(m+A\_{\infty})\ddot{z} + \int\_{0}^{t} \mathbf{K}\_{rad}(t-\tau)\dot{z}(\tau)\mathbf{d}\tau + \mathbf{C}z = F\_{\text{exc}} + F\_{\text{pto}} + F\_{\text{hs}} \tag{2}$$

where *m* is a buoy mass, *A*∞ is the infinite-frequency added mass coefficient, *C* is the hydro-static stiffness, *Krad*(*t*) is the radiation impulse response function, *Fexc* is the wave excitation force, *Fpto* is the load force exerted on the buoy from the power take-off system [62]. Free surface elevation height results from a linear superposition consisting of some wave characteristics in irregular waves. This is usually determined by a wave spectrum which describes the distribution of energy in a vast number of wave frequencies. Significant wave height and peak period are utilized as the basic identification of the wave in the spectrum. The irregular excitation force can be calculated as the real part of an integral term across all wave frequencies as follows.

$$F\_{\rm exc} = \mathbf{R} \left[ \int\_0^\infty \sqrt{2S(\omega\_r)} F\_\mathbf{x} \, e^{i(\omega\_r t + \phi)} \, \mathbf{d}\omega\_I \right] = \int\_{-\infty}^{+\infty} \eta\_{(\tau)} f\_\mathbf{e}(t - \tau) \mathbf{d}\_{(\tau)} \tag{3}$$

where **R** denotes the real part of the equation, *F<sup>x</sup>* is the excitation vector consists of amplitude and phase of the wave, *S* is the wave spectrum, *φ* is the stochastic phase angle, *η<sup>τ</sup>* represents water elevation and *f<sup>e</sup>* is the element of force vector [48]. The load force of PTO is modeled as a linear spring-damper system.

$$F\_{pto} = -B\_{pto}\dot{z} - K\_{pto}z \tag{4}$$

$$\begin{aligned} F\_{\rm hs} &= -K\_{\rm hs,min}(z - z\_{\rm min})u(z\_{\rm min} - z) \\ &- K\_{\rm hs,max}(z - z\_{\rm max})u(z - z\_{\rm max}) \end{aligned} \tag{5}$$

where in Equation (4) *Kpto* and *Bpto* are control parameters which represent stiffness and damping of PTO and in Equation (5) *u* is the Heaviside step function, *Khs*,min and *Khs*,max are the hard stop spring coefficients, and *z*min and *z*max are the stroke limits which are related to the nominal position of the converter. It is important to note that, for computing useful absorbed energy, the effect of this force is not considered [63].

In order to calculate the energy produced by each buoy, the sum of three forces is necessary: wave excitation(*Fexc*,*p*(*t*)), force of radiation (*Frad*,*<sup>p</sup>* (*t*)), and power take off force (*Fpto*,*p*(*t*)). The scattered irregular waves are included in the wave field when computing the excitation force. Furthermore, the stiffness and damping parameters of the PTO system at the end of each tether along with hydrodynamical parameters are taken in order to compute the total power output of an array. To calculate the average power absorbed by the array, several variables have to be taken in to account, as follows.

$$P\_{\rm n}(H,T) = \int\_0^{2\pi} \int\_0^{\infty} \mathcal{Z} S\_{\rm n}(\omega) D(\beta) p(\beta, \omega) \mathrm{d}\omega \mathrm{d}\beta \tag{6}$$

where *Pn*(*H*, *T*) is the average power absorbed by the array in a regular wave of unit amplitude, *Sn*(*ω*) is the irregular wave spectrum which is calculated with the Bretschneider spectrum and *D*(*β*) represents the directional spreading spectrum, particularly for this site which is come from the wave rose [64]. *ω* is the wave frequency and *p*(*β*, *ω*) is the power function of each submerged buoy defined by Equation (7).

$$p(\boldsymbol{\beta},\omega) = \frac{1}{2}D\_{pto}\omega^2\Gamma(\boldsymbol{\beta},\omega)^2\tag{7}$$

where Γ(*β*, *ω*) is the response amplitude operator (RAO) of the productive degree of freedom of the buoy obtained by solving the equation of motion from Equation (2), and *Dpto* is the Power Take-Off (PTO) damping. The wave angle is based on the *z*(*β*, *ω*) [40], which can be calculated by equation (2) at the beginning of this section. The array at a certain test site is generated by total mean annual power *Parray*, and to calculate that, the contribution of energy absorption from a wave climate in each state can be summarized as:

$$P\_{array} = \sum\_{n=1}^{N\_{\rm s}} O\_n(H\_{\rm s}, T\_p) P\_n(H, T) \tag{8}$$

where *N<sup>s</sup>* is a number of chosen sea state, *H<sup>s</sup>* is the significant wave height and *T<sup>p</sup>* is the peak wave period for each sea state, *On*(*H<sup>s</sup>* , *TP*) represents the probability of occurrence of sea state which stems from the wave scatter diagram and *Pn*(*H<sup>s</sup>* , *Tp*) is a power which the array produces in the *n*th sea state [40]. Significant wave height and peak wave period are statistics of a sea state which can refer to the condition of the ocean/sea surface. To calculate *Pn*(*H<sup>s</sup>* , *Tp*) in irregular waves, it is necessary to sum all power contributions in each frequency and significant wave direction.

#### *2.2. Wave Resource*

According to previous works [40,45], four different sea sites were chosen for this study. The wave height directional distribution (wave rose) can be seen in Figure 2, at the chosen locations (as an example, the wave condition at the Sydney site is shown). The wave rose shows that the significant wave directions are from 15 degrees to 190 degrees, where 90 percent of the incident waves travelled. Consequently, the dominant wave direction is from the south.

**Figure 2.** The wave rose plot at Sydney.

Each array is constrained by the maximum area and the minimum distance between WECs. Firstly the minimum separation between buoys (*R* ′ ) has to be 50 m to provide a safe pass for vessels. Secondly, although the area grows by increasing the number of buoys, it has to be constrained within the area Ω, where <sup>Ω</sup> = *<sup>l</sup>* × *<sup>ω</sup>*, *<sup>l</sup>* = *<sup>ω</sup>* = √ *N* × 120,000 m [45].

The Bretschneider spectrum is used for modeling irregular waves in this study. This spectrum is a modified Pierson-Moskowitz spectrum which is based on significant wave height and peak period. These two parameters are highly dependent on wind speed and its direction [65].

$$S(f) = \frac{H\_{\text{n0}}^2}{4} (1.057f\_p)^4 f^{-5} \exp\left[-\frac{5}{4} (\frac{f\_p}{f})^4\right] \tag{9}$$

where *Hm*<sup>0</sup> and *f<sup>p</sup>* are the significant wave height and the frequency of the peak wave period, respectively.

#### *2.3. Array Interaction Criteria*

The optimal designs of the array for four different locations in Australia use power matrices of various configurations (i.e., different layout geometry, WEC distance and relative angles regarding dominant wave direction wave directions). To be more precise, the goal is to select the best site for each layout configuration with the optimal separation among WECs and rotation angle, namely, the one that provides the highest annual energy output per each converter. For this aim, based on the number of WECs, different layouts can be deployed with various orientations and separation among converters. Note that there are certain constraints for distances, and this depends on the number of WECs. Similarly, the number of

converters in an array allows configurations to be chosen. Hence, it is crucial to measure the effectiveness of interactions between converters by the q-factor coefficient. The q-factor is shown to be an important evaluation criterium such that if *q* > 1, then it has a positive effect on the total energy of an array; otherwise, the interactions are destructive.

$$q = \frac{P\_{array}}{NP\_{isulated}}\tag{10}$$

where *Pisolated* is the power that an isolated WEC generates, *N* is the number of converters [41]. As Equation (10) indicates, there is a direct relationship between q-factor and power output of each array; however, both of them have to be studied separately, due to different objectives of finding optimal values for each parameter. The maximum feasible amount of q-factor is investigated to achieve the best constructive effects of interactions of buoys in an array. On the flip side, since the power output of a WEC array plays a significant role in the assessment of the system's response to energy consumption needs in coastal areas, this parameter is considered along with the q-factor. The equation below calculates the mean q-factor by considering the number of converters, variety of wave directions and allowable distance within a 5 m interval.

$$\textbf{mean q-factor}\_{\text{(each wave scenario)}} = \frac{\sum\_{i=0}^{\text{maxa}} \sum\_{j=50}^{l} \textbf{q}\_{factor\_{(i,j)}}}{\text{total number of cases}} \tag{11}$$

where *α* is the direction of wave n 10-degrees resolution, except when *N* is 5, the interval changes to 9 degrees ranging from 0–63 degrees and *l* is the maximum allowable distance between buoys within the area. This mean q-factor coefficient will be considered to find the best location for the max q-factor over different angles and distances.

#### **3. Layout Assessment Routine**

According to the mentioned equations in Section 2, the following outcomes are obtained. Four different layouts are considered regarding the number of buoys, and they are thoroughly described in detail. By looking at the mentioned literature, it is evident that using more converters results in more potential destructive interaction. Specifically, the q-factor may decrease when the number of buoys rises to greater than five [54]. Therefore, we decided to choose the five buoy layout as the maximum complexity for the model; however, evaluating the larger wave farm characteristics is a part of our future research plan. Furthermore, the symmetric design is proposed for the layouts based on the following reasons: (i) To find a single variable to handle, the buoy-buoy distance must remain constant for each configuration. Thus, the symmetric configuration fulfills this requirement in our assessment, like in previous studies [66]. (ii) To cover the whole area of study in power absorption assessment process, the array rotates 10 degrees, regarding the dominant wave direction, in each evaluation. The asymmetric effect of the array configuration rules out the duplicate assessment, resulting in less computational cost.

In the first step, There are two buoys in this array, and one line connecting them. The dominant wave direction indicates the direction in which most waves travel. However, in the calculation of power output from each buoy or an array, the significant wave directions are used. These are directions from where 90 percent of the waves are traveling. Furthermore, the resolution size of the evaluation is 15 degrees, by which we detect the dominant direction of the waves.

The dominant wave direction is obtained by considering one-third of the maximum waves in the wave rose. The angle between dominant wave direction and a hypothetical line is considered to be alpha (*α*), which is clearly illustrated in Figure 3. The interval of alpha is chosen to be tested every 10 degrees; therefore, there would be 18 different alpha ranging from 0 to 170 degrees. This range has been considered to prevent the extra calculation of results that have already been calculated. When there are three buoys to

consider in an array, one of the most common geometries is the equilateral triangle. If some lines are used for connecting these buoys, angles between vertexes of the triangle will be 60. There is a line from this converter perpendicular to the line, which connects the two other buoys. The angle between the dominant wave direction and this perpendicular line is alpha, shown in Figure 3, and this parameter has twelve degrees from 0 to 110, which changes every 10 degrees. A regular quadrilateral is taken into account to configure four buoys. To describe alpha in this layout, firstly, the dominant wave direction needs to be determined. Secondly, a hypothetical line should be drawn from one converter to the furthest one. For example, if the converters are numbered clockwise and the closest buoy to the front wave is buoy number one, the line should be drawn from 1 to 3, exactly like in Figure 3. Finally, the angle between this line and dominant wave direction is alpha, and the range of this is from 0 to 80 degrees, which has nine different amounts with equal intervals.

In this layout, five similar converters form an array are shown in the shape of a regular pentagon. The dominant wave direction is illustrated in Figure 3 with a blue arrow. Converters are numbered clockwise, and the first number starts from the closest buoy to the front wave. As shown in Figure 3, each converter has the longest distance with two buoys. In this case, the furthest converter to buoy number 1 is number 3 and 4. If a perpendicular line is drawn from the first converter to the connecting line between furthest converters, the angle between the dominant wave direction and the perpendicular line represents alpha. The range for alpha is from 0 to 63 in 9-degree intervals, so there are eight alphas to test in this layout. Distances are also assumed to change every five meters between the allowable period. In the end, three measurements, which are the power output of each buoy, array power, and q-factor, are taken in each step for all layouts, separately. The details of all results are discussed comprehensively in Section 5.

**Figure 3.** Layout Setup of (**a**) 2 buoys array (linear), (**b**) 3 buoys array (triangle-shape), (**c**) 4 buoys array (square-shape) and (**d**) 5 buoys array (pentagon-shape), with regard to dominant wave direction.

#### **4. Results and Discussions**

This section represents the results of different array layouts when the number of buoys rises from 2 to 5 in considered locations on the Australian coast. The results demonstrate the sensitivity of the array power and the q-factor due to the changes in buoy-buoy distance and rotation angle. It is worth mentioning that there are 16 conditions in this study, which will be discussed in detail as follows. To choose the optimal degree in this section, one of the most important variables is alpha, whose optimum value leads to the average maximum amount of the power output.

#### *4.1. Sensitivity of Two-Buoy Array Performance To Distance*

Figure 4 shows the sensitivity analysis of the array power output to the different buoy-buoy distances. It can be seen that Tasmania has the most wave array power, which is almost 0.534 Mw where the *α* is 80 degrees and the buoy-buoy distance is 160 m. The 60-degree angle line, which has the most average power, rises with a sinusoidal trend from the beginning to 200 m. Then, it increases gradually. The second location is Sydney, which has a considerable array power. Although its array power is 0.218 Mw, which is far less than Tasmania, the maximum average array power occurs in 130 degrees with a similar trend to the mentioned location. The maximum power that Sydney's layout determines is achievable when the distance is around 400 m. Adelaide and Perth are similar in terms of the array power range, which is roughly from 0.18 to 0.196 Mw. As the green line in these figures shows, when the rotation angle and buoy-buoy distance are 40 degrees and 165 m, respectively, both reach the highest array power. The maximum average of array power can be witnessed in 20 degrees in Adelaide and 30 in Perth. In the two mentioned sites, it is apparent that figure lines follow different trends. To compare, when *α* is 20 degrees, and buoy-buoy distance is 100 m, the first peak of the array power is observed in Adelaide. Next, it falls until the distance goes over 150 m. The first peak in Perth happens when the distance is near 60 m. Then, it remains unchanged for the next 60 m. After significant growth, it reaches around 0.194 Mw with 160 m buoy-buoy distance. Overall, it is remarkable that the maximum power can be harnessed in three locations when the distance is 160 m.

**Figure 4.** Array power of the two-buoy layout over different distances in four wave models.

#### *4.2. Sensitivity of Three-Buoy Array Performance to Distance*

It can be seen that in Figure 5, the most obvious inferred outcome is that, the longer distance between WECs leads to more extracted power output. However, widening the area might not be proficient because the line only rises 0.2 Mw by increasing the distance from 250 to 500 m. The maximum harnessed power output can be seen when *α* is 10 degrees, except in Sydney, which is 110. The range of array power is quite narrow in the mentioned locations, where only a 0.05 Mw gap can be witnessed among 12 tested angles. The main reason for obtaining the similar results among different angles' experiments can be explained as follows. Where the three converters are placed in equilateral triangle geometry, there would always be two buoys in the zone of radiation. Therefore, the changes in *α* cannot produce considerable effects. By comparing the power output of WECs over the changes of distances, we can see the same overall trend has been followed in all studied locations. A sharp rise in array power can be achieved by increasing the distance up to 100 m, followed by a gradual rise by increasing the distance up to the maximum allowed size. It can be mentioned that differences between maximum distance in each layout relate to the area constraints, which have already been discussed in Section 2.

**Figure 5.** Array power of the three-buoy layout over different distances in four wave models.

#### *4.3. Sensitivity of Four-Buoy Array Performance to Distance*

The geometry chosen for four converters is a square shape. *α* values in this layout range from 0 to 80 degrees with 10-degree intervals. The highest harnessed array power is 1.05 Mw in Tasmania, while the lowest one observed in the Perth layout is around 0.387 Mw. Turning to the rotation angle, for Perth and Sydney, the angle is 40 degrees to extract the maximum average array power; however, for Tasmania and Adelaide, *α* is 0 and 80 degrees in order of appearance. Considering the distances between WECs, it is interesting that where the buoy-buoy distance is between 150 and 200 m, the maximum wave array power can be exploited in Perth, Adelaide, and Tasmania. However, in Sydney, it seems that in this *α*, the array power evens off after reaching the highest amount. Hence, the minimum distance between WECs is more cost beneficial for considering layout design 160 m of distance takes into account. In contrast, in Adelaide and Tasmania, a 160 or 170-m distance seems to be the best distance between converters. This amount is a bit greater in Perth, where the highest array power is firstly seen in the 180-m distance (Figure 6).

**Figure 6.** Array power of the four-buoy layout over different distances in four wave models.

#### *4.4. Sensitivity of Five-Buoy Array Performance to Distance*

In this study, the configuration of five converters is chosen as a regular pentagon array. This is because there is no difference between each WEC, and the range of rotation angles is restricted to be between 0 and 63 with eight different angles. The maximum averaged array power in all four case studies is witnessed when the *α* is either 18 or 63 degrees. To be more precise, in Tasmania and Perth, the *α* is 18 degrees, and for the other two, it is 63 degrees. Also, it is evident that when the distance is between 200 and 250 m, the maximum power output is harnessed in all wave scenarios; and the optimal choice can be found in the mentioned range consequently. As Figure 7 shows, the trend of all case studies are similar, except in the 18-degree's line in Tasmania and Perth, where the array power reduces gradually after the peak, instead of leveling at the peaks power. By comparing this result with recent similar studies, we can see the same trend of absorbed power by raising the distance between WECs in each layout up to an optimal value, after which the results were roughly stable [54,67].

**Figure 7.** Array power of the five-buoy layout over different distances in four wave models.

#### *4.5. Sensitivity Analysis of q-Factor to the Relative Angle of Rotation*

The interaction of converters is measured with a well-known parameter called q-factor. A sensitivity analysis has done by monitoring the q-factor distribution over different rotation angles of the WECs in each layout. Since the interaction of the buoys in each layout could be constructive or destructive, the q-factor is calculated per 5 m of distance between converters. Figures 8–11 show the distribution of results of calculated q-factors over different relative angles (*α*) in a box chart, describing the values as they spread across the entire range. In each angle, there is a box that reveals the amount of fifty percent of q-factor results. Also, the middle line indicates the mean value of all results. The other amounts of q-factor, which are far from the mean values (i.e., the greatest 25 percent and the least 25 percent of the results), are shown by two lines located above and below the rectangular box.

There is a lot of similarity between Adelaide and Perth in terms of their q-factors, but Sydney and Tasmania have different trends. In both Adelaide and Perth, as Figure 8 and 9 show in two-buoy layout, a fluctuating pattern is witnessed which indicates the importance of the rotation angle of the array power with respect to the dominant wave direction. Thus, when *α* is between 30 and 40 degrees, the highest q-factor is achieved, and the layout design process should be followed by choosing the best buoy-buoy distance in the mentioned angle. The distinction between q-factors in the three-buoy layout is negligible due to the effects of dominant wave direction on the equilateral triangle layout, in which one converter affects two others by radiated waves. The maximum q-factor can be seen in the 30 and 40 degrees area. In the four-buoy layout, when *α* is 40 degrees, the q-factor is around 1 in both locations. Among the 9 discussed dominant wave directions, the highest q-factors are seen when *α* is 20, 40, or 60 degrees, and the average q-factor in each direction is a bit over 0.96. In the five-buoy layout, it is clear that the average q-factors are between 0.95 and 0.98, and the highest q-factors happen when *α* is either 20, 30, 60, or 70 degrees, and the range of q-factors is from 0.85 to 0.98. It is important to note that in the mentioned degrees, q-factors are mostly close to the highest amount because the average line is on top of each box.

**Figure 8.** q-factor results distribution and mean value per five meters of the Wave Energy Converters (WECs)' distance over rotation angle due to significant wave direction in the Perth wave model. (Fifty percent of results near the mean value are plotted in a box, the range of other results is shown by a dashed line).

**Figure 9.** The q-factor results distribution and mean value every five meters of WECs' distance over rotation angle due to the significant wave direction in the Adelaide wave model. (Fifty percent of results near the mean value are plotted in a box, the range of the other results is shown by a dashed line).

In Tasmania, due to the symmetry between WECs and small effects of changing *α* in q-factor for the three-buoy and five-buoy layout, changes are not considerable. The average q-factor in each rotation angle is approximately 0.98 and 0.96, respectively. When *α* is between 40 and 90 degrees in the two-buoy layout, q-factors in each distance are between 0.95 and 1.005. The maximum averaged array power occurs at 20 degrees. In the four-buoy layout, maximum q-factors are observed in four rotation angles, which are 0, 20, 40, and 60 degrees. Further details can be found in Figure 10.

**Figure 10.** The q-factor results distribution and mean value per five meters of WECs' distance over rotation angle due to significant wave direction in Tasmania wave model. (Fifty percent of results near the mean value are plotted in a box, the range of other results are shown by a dashed line).

Looking at Sydney in Figure 11, in the two-buoy layout, it is evident that maximum q-factors happen when *α* is either between 110 and 120 degrees or 130 and 140. One of the distinctions compared to the mentioned locations is that the lowest q-factor is seen at 40 degrees. In the three-buoy and four-buoy layout, the average q-factor in each *α* is around 0.98 and 0.97, respectively, and its changes are not recognizable in all 12 tested angles. The q-factors in the five-buoy layout ranged from 0.84 to 0.98. The closest q-factor values to 1 are found at 10 and 50 degrees. Moreover, it can be inferred that when *α* is 30 or 70 degrees, the related q-factors are near 0.98. These results would help further feasibility studies of the WECs' array analysis by presenting the possible range of achievable q-factors in Perth, Adelaide, Tasmania, and Sydney ports.

**Figure 11.** q-factor results distribution and mean value per five meters of WECs' distance over rotation angle due to significant wave direction in the Sydney wave model. (Fifty percent of results near the mean value are plotted in a box, the range of the other results is shown by a dashed line).

#### *4.6. Landscape Analysis*

Figures 12–15 reveal the power for each buoy in four different layouts. Overall, it is inferred that the asymmetry in arrays makes the power more predictable in each distance and rotation angle. The illustrated plots in this section indicate four kinds of buoy layout for each area in Adelaide, Tasmania, Sydney, and Perth. For each layout, the Colour-bar presents the amount of total power per buoy. This amount has been shown in the form of a contour. In general, there is a similar pattern for studied locations by considering the configuration of the arrays. The maximum amount of extracted energy is more likely to be found in the 5-buoy layout. Moreover, the power of each converter in the center area is not considerable, and the maximum amounts of the exploited energy found for higher buoy-buoy distances. The asymmetry of these energy distributions is high as well, but Perth is an exception. To observe abrupt changes, an increase in the resolution of distance and angles are needed. Take the 2-buoy layout in Adelaide as another example of non-asymmetric layout; the distribution of the incident waves over different angles implied such non-asymmetric contour of exploited energy. The reason behind the asymmetry in other layouts would be that by increasing the number of buoys, possible shadowing effects proportional to each rotation angle occur in every layout. In some placements of the 2-buoy layout, as the array experiences different rotation angles, the buoys may see the same dominant wave direction; however, in some angles, the shadowing effect of one buoy over the other can play a crucial role in reducing the

energy. This shadowing effect of buoys to each other becomes more drastic, as the number of buoys increases because, in each angle, there is more chance of interaction between at least two buoys. To describe blue spots it should be noted that in this research, the minimum separation between buoys is considered to be 50 m. Then, by assessing the different angles, the configuration of the array rotates over the center point, such that the lowest amounts of energy may be witnessed in the middle of the figures by interpolation.

Taking a look at Figure 12, it is interesting that the order of maximum power for each layout is around 0.1 to 0.105 Mw in most areas. However, this range increases in the two-buoy layout, so more areas with dark red colour can be seen. In this figure, asymmetry in each layout is more evident than in the other locations.

**Figure 12.** Exploited energy distribution of the WECs array over entire area in Sydeny wave scenario.

There are many similarities between Adelaide and Perth in Figures 13 and 14. For instance, in a two-buoy layout, their contour has resemblance, and the range of power is identical. Furthermore, more power is extracted in Perth based on these plots. To compare the four-buoy layout, Adelaide has more symmetrical power distribution, and the chance of reaching 1 Mw power is more in Adelaide in general. Finally, there is a small difference between choosing Perth or Adelaide as the installation site among the other surveyed sites.

**Figure 13.** Exploited energy distribution of the WECs array over the entire area in the Perth wave scenario.

**Figure 14.** Exploited energy distribution of the WECs array over the entire area in the Adelaide wave scenario.

The highest amount of power can be exploited from converters in Tasmania, which is obvious in Figure 15. Although seemingly green areas in the two-buoy layout cover the majority of the zone, the power that belongs to the green areas is very close to the other layouts of the orange ones. The chance of extracting over 0.27 Mw power is seen in the five-buoy and four-buoy layouts. It is worth considering the three-buoy layout power when the *X* axis is over 250 m, which reveals the potential of this layout under certain conditions. Likewise, when the *Y* axis is less than 150 or more than 350, this potential is met.

**Figure 15.** Exploited energy distribution of the WECs array over the entire area in the Tasmania wave scenario.

#### *4.7. Interaction Based Layout Selection*

Figure 16 shows the maximum and mean value of the q-factors in a given number of buoys in each wave model. The most significant observations inferred from Figure 16 are addressed as follows. The maximum q-factor in Tasmania and Sydney is less than in the other locations because of the lack of constructive interactions to compare to the other layouts. The mean q-factors are also higher in Perth and Adelaide in all locations for the same reason. Turning to the maximum q-factor, it is apparent that the highest constructive interactions in the two-buoy layout occur in Adelaide. However, in Sydney's wave scenario, installing buoys, whether separately or in an array, is almost the same because the q-factor equals 1 in the best-case scenario. Although this amount is over 1 in all locations of the four buoy layout, Sydney is an exception. In the three-buoy layout and five-buoy layout, only the q-factor of Perth is more than 1. It is worth considering that the latter has the least maximum q-factor. These results confirm results from previous studies on the decrease of q-factor after increasing the number of WECs arrays, specifically after incorporating more than five buoys [54]. Turning to the mean q-factor, it is evident that by increasing the number of buoys, this variable decreases, and a reduction of almost 0.07 is seen by adding a buoy. Also, this measure is observed to be a trend because constructive interactions are more likely to be seen in Perth and Adelaide. These interactions will occur if the *α*, buoy-buoy distance, and the geometry of layout are appropriately chosen.

Finally, it has to be noticed that further details and numbers are written in Table 2, which enables a comparison between each location. The bold numbers in Table 2 represent constructive interactions between converters. Therefore, in those cases, installing an array is more efficient.

**Figure 16.** Comparison of maximum and mean q-factor in different wave scenarios in each layout.


**Table 2.** best solutions related to maximum q-factor in each wave scenario.

#### **5. Conclusions**

Investigating for an appropriate arrangement an array layout constitutes a complicated problem in wave energy projects. Wave energy converters can reinforce each other to provide more power output in the form of an array if the distance among them is efficiently adjusted and the arrangement of the layout appropriately defined. In this paper, we analyzed the CETO6-project WECs separation in an array with different numbers of devices and arrangements. In order to assess the impact of various wave models, we perform and compare all numerical analyses in four real wave scenarios including the Sydney, Perth, Adelaide, and Tasmania sea sites. According to the numerical analysis, there is a direct relationship between the number of converters and optimal inter-distance among them and also relative angle to the significant wave direction. Greater separation between converters leads to more array harnessed power output. However, the most exploited energy can be achieved in 2 buoy layout with a 165 m buoy-buoy distance. A sensitivity analysis has revealed that the q-factor distribution differed due to different rotation angles of the WECs array. Moreover, the maximum q-factor output analysis showed that results in a two-buoy layout in all scenarios, tree-buoy layout in Perth and four-buoy layout in all scenarios (excluding Tasmania) are far higher than the other locations' q-factor, and this parameter is almost the same in the five-buoy layout sea sites. However, the landscape analysis-approved maximum amount in terms of the extracted net power output was found in the 5-buoy layout in the Tasmania wave scenario.

**Author Contributions:** Conceptualization, E.A.; Data curation, D.G., M.M.N., M.N.; Formal analysis, E.A., D.G., F.A.; Investigation, E.A., D.G., F.A., M.M.N. and M.N.; Methodology, E.A., D.G., F.A., M.M.N. and M.N.; Resources, M.N.; Supervision, D.A.G., F.A.; Validation, E.A., D.G., F.A., and M.M.N.; Visualization, E.A., D.G., and F.A.; Writing—original draft, E.A., D.G., M.M.N. and M.N.; Writing—review & editing, M.N. and D.A.G. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Acknowledgments:** The authors would like to express their gratitude to ODYSSEA project that received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement No. 727277. Furthermore, the authors would like to appreciate Nataliia Surgiienko from the University of Adelaide due to publishing the MATLAB source code of the wave energy simulator. This work has been supported by the High Performance Computing Research Center (HPCRC)—Amirkabir University of Technology under Contract No. ISI-DCE-DOD-Cloud-700101-4504.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Abbreviations**

The following abbreviations are used in this manuscript:


#### **References**

1. Barstow, S.; Mørk, G.; Mollison, D.; Cruz, J. The wave energy resource. In *Ocean Wave Energy*; Springer: Berlin/Heidelberg, Germany, 2008; pp. 93–132.





**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
