Modelling of a wave energy converter array with non-linear power take-off using a mixed time-domain/frequency-domain method

A mixed time-domain/frequency-domain method is proposed for modelling dense wave energy converter (WEC) arrays with non-linear power take-off (PTO). The model is based on a harmonic balance method which describes the system response in the frequency domain, while evaluating the non-linear PTO force and solving the system equations of motion in the time domain. The non-linear PTO force is computed with Lagrange multipliers. In order to apply the proposed method for WEC array responses in real sea states, the time series is split into time windows and the simulation is carried out individually per window. The method is demonstrated by investigating the dynamics of the Ocean Grazer WEC array (OG-WEC) with an adaptable piston pumping system. The key parameters thought to possibly inﬂuence model accuracy, including the number of harmonic components, the length of the time window and overlay, are discussed. It is shown that the proposed model can signiﬁcantly reduce the computational cost with an acceptable accuracy penalty.


INTRODUCTION
Ocean wave energy is a sustainable and abundant energy source but extracting energy from ocean waves has not become a commercially viable technology yet. The costs, survivability, and power quality are still major obstacles. Wave energy converter (WEC) arrays may partly overcome the current obstacles, as it can potentially increase the overall energy production, reduce the operation and maintenance costs and smoothen the power output. WEC arrays, in particular dense ones, are prone to exhibiting complex and non-linear behaviours due to strong wave interactions between the buoys and coupling effects with the non-linear power take-off (PTO) system. Development of an accurate and computationally cost-effective numerical model for array configuration optimization and WEC array production prediction is therefore of great value. It is common to use non-linear PTOs (e.g. hydraulic PTO) in WEC design. Cargo et al. [1] presented a study of implementable active tuning methods for WECs with hydraulic PTOs, and they found that over-simplification of the PTO in simulations was inadequate for WEC design studies. Penalba and Ringwood [2] presented a high-fidelity wave-to-wire model for WECs This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. © 2021 The Authors. IET Renewable Power Generation published by John Wiley & Sons Ltd on behalf of The Institution of Engineering and Technology with various non-linear PTOs and including different conversion stages, which can be used for WEC design and optimization. Hansen et al. [3] developed a Simulink-based time domain model to investigate an array with 20 devices combined with hydraulic PTOs, but the computational times were prohibitive for real time operation. These studies have shown the significance of assessing the non-ideal efficiency of the PTO and taking into account the non-linear effects of the PTO for the WEC (array) performance evaluation. The development of accurate and highly efficient numerical models is of great value for such purposes.
Conventional frequency-domain models have been widely used for WEC array studies [4,5]. Accompanied by a boundary element method, they can account for arbitrary geometries of WECs and wave interactions with low computational effort. Despite their computational advantages, one main drawback with this approach is its limitation in dealing with non-linear effects. On the other hand, time-domain models incorporating, for example, the Cummins' equation with a convolution operator, can incorporate such non-linearities, albeit with much higher computational costs. Alternatively, the state-space representation can approximate the representations of the convolutions, which can dramatically increase computational speeds [6]. The state-space approach has been widely employed for single WECs and small size WEC array studies [7]. However, this approach is rarely applied to large size, dense WEC arrays because of the difficulty of dealing with cross-coupled terms [8]. The more advanced computational fluid dynamics (CFD)-based models can describe the comprehensive flow details of waves and WEC array interactions. However, their applications are still limited to cases with a small number of WECs subjected to regular waves [9] as the computational cost is extremely high or even unacceptable.
The Harmonic Balance Method (HBM) is a popular method to approximate the frequency response of non-linear systems, as it is a computationally efficient method to obtain steady-state responses for non-linear dynamics problems and is relevant for WEC systems; examples include non-linear electrical circuits [10,11], non-linear mass-spring-damper systems [12][13][14][15] and computational fluid dynamics problems [16][17][18]. Recently, an HMB-based frequency-domain representation of WECs subject to non-linearities was proposed by Bacelli and Ringwood [19] and Merigaud and Ringwood [20]. In these papers, the displacement and the non-linear force were approximated by a truncated Fourier series. Thus, instead of the time consuming convolution integral in time-domain modelling, the HBM solves a set of equations of motion in the frequency domain. The HBM was extended to WEC real-time control applications [21] due to its competitive computational performance. In their work, an explicit Jacobian was applied in order to speed up the convergence. The paper [22] demonstrated that the HBM can be used to model large dense WEC arrays with non-linear pumping forces. They applied a numerical Jacobian in the computation due to occurrence of discontinuity in the pumping force, hence the vector field is not differentiable everywhere. Although the classical HBM is simple in its principle, it becomes computationally inefficient when a large number of harmonics is required, for example for irregular wave simulations [23]. This is due to the large amount of numerical integration operations required when computing the harmonic coefficients of a nonlinear force.
The present work extends the classical HBM model and propose a mixed frequency-domain/time-domain method with a windowing technique for a large dense WEC array with nonlinear PTO. The appealing feature of the proposed method is that it can describe the non-linear dynamics of the complex WEC array in random sea states while significantly reducing the computational cost. In this model, the real sea state is split into time windows with uniform length; thus, the state of each time window can be described with a finite number of harmonic components. A mixed frequency-domain/timedomain method was applied for each time window. The nonlinear term is calculated with a time-marching procedure in the time domain, using the Lagrange multipliers method. The OG-WEC device was used to demonstrate the validity of the proposed approach. This device consists of a large array of single floaters, each connected to a PTO utilizing adaptable piston pumps that introduces a strong non-linearity into the system. The paper is organized as follows. The system equations of motion of the OG-WEC and the numerical solution are described in Section 2. The model validation and the discussion about the key parameters are presented in Section 3. Finally, conclusions are drawn in Section 4.

Motion equations of OG-WEC system
The Ocean Grazer is a hybrid renewable energy device which combines adaptable WEC technology with on-site energy storage to harvest renewable energy offshore. An OG-WEC system comprises a finite number of buoys with a wind turbine in the centre and its storage system is housed in a gravity-based concrete structure. The OG-WEC is designed to react to the high variability of waves by utilising an adaptable pumping system -MP 2 PTO (multi-piston, multi-pump power take-off) systemand store potential energy in the storage system. The detailed description can be found in Wei et al. [24]. Penalba et al. [24] applied a numerical model to investigate the hydrodynamic interactions in three WEC arrays with different sizes, and suggested that arrays with about 20 devices were hydrodynamically more efficient. We applied the numerical model to an OG-WEC system with 18 buoys, arranged as depicted in Figure 1(a). All buoys have a uniform geometry (see dimension in Figure 1(b)) and are connected to the central pillar via trusses (not sketched) which restrict their motion in heave only. Each buoy is linked to a controllable transmission system which drives the hydraulic piston-type PTO system to pump the working fluid from the inner reservoir to the outside flexible bladder. The stored potential energy can be converted into electricity with hydro-turbines. In the model, we further assume that the cable linking the buoy and the pumping system is stiff enough so that we can simply describe the piston displacement as the buoy displacement adjusted via a transmission ratio.
The motion of the buoy is governed by the following equation: where M b is the mass of the buoy, X b is the heave displacement of the buoy, and F e , F r , F hs and F c represent the excitation, radiation, hydrostatic restoring and cable forces, respectively. The motion of the piston in an adaptable piston pump can be described by: where m p is the mass of the piston, is the transmission ratio, F p is the non-linear pumping forces which is expressed as:

FIGURE 1
Sketch of the OG-WEC system: (a) honeycomb arrangement of an array; (b) geometry of a cone-cylinder buoy; and (c) schematic of the power take-off system. The pillar and rigid reservoir are maintained the atmospheric pressure (ATM). Reprinted from Wei et al. [22] where is the density of the working fluid (assumed in this work to be the same as sea water density), D is the water depth, h r is the depth of the reservoir, l p is the length of the piping between the reservoir and the bladder,̇z p andz p are the piston velocity and acceleration, and A c is the area of the check valve (used to adapt the piston pump to the wave excitation). The pumping force F p is very large during the upstroke but becomes zero during the downstroke. In the upstroke component of the pumping force, the first term represents the hydraulic head which saves potential energy and is the largest of the three terms; the second term represents the inertia effect which can lower the resonant frequency of the system; the third term is the kinetic energy which is considered as loss in the system. The present paper focused on the energy extracted from the waves, while the energy losses in the pumping system, electricity generation, and electricity to wire were not accounted for in this study. These losses should be investigated in further studies, for example in the context of a power matrix assessment.
Combining (1) and (2), we obtain the governing equation of one OG-WEC unit: ( Assuming a steady state periodic response for the WEC system, X b , F e , F r , F hs and F p are approximated as a truncated Fourier series: where N h being the number of harmonic components, is the frequency step for the signal with period T p = 2 , X b,n corresponds to the nth harmonic complex amplitude of the buoy displacement,F e,n A r,n and B r,n are the excitation force, added mass and radiation damping coefficients , andF p,n is the nth harmonic pumping force coefficient. For a dense WEC array, interactions between WECs are considered via cross-coupling hydrodynamic coefficients obtained by NEMOH. Substituting Equations (5)-(9) into Equation (4), the equation can be further transformed to the frequency domain with a set of linearised equations over each harmonic. Sorting out the terms with the same frequency, the n-th harmonic motion equation for the ith WEC, is expressed as: The matrix form of the motion equations of the WEC array can be written as: (11) where Z r is the condensed dynamic stiffness assembled from (10). The vectorX b contains the unknown Fourier coefficients of displacement, with N b being the number of buoys. Solving (11) with unknowñ X b is not straightforward asF p is undetermined.

A mixed frequency-domain/time-domain model
According to (3), the pumping force is a piston state-dependent variable, which cannot be explicitly obtained in the frequency domain. Alternatively, it is more convenient to calculate the pumping force using a time-marching procedure. The procedure of the A mixed frequency-domain/time-domain (MFT) approach is structured as follows: a. With the input of the time series of wave elevation per time window (either waves pre-generated with a defined wave spectrum or real time measurements), the pre-processing is carried out to transfer the wave elevation into the frequency domain; then, the wave excitation force is obtained in frequency domain. The buoy displacement is initialized as zero. b. The excitation force is transferred into the time domain, and the time series of the pumping force is predicted with the Lagrange multiplier method (see Section. 2.3), using time domain buoy displacement as input. c. The system equations of motion are solved in MFT form, cf.
(15). d. An iterative procedure is applied to each time window to obtain converged results. The iteration terminates once the system residual is smaller than a given tolerance. Otherwise, the buoy displacement is updated, the procedure returns to step 2 and the iteration continues.
For temporally periodic oscillations, the Fourier coefficients and the time variation can be transferred back and forth via a discrete Fourier transform (DFT) and its inverse (iDFT), for where T is the DFT matrix, and T inv is the iDFT matrix, and X b is the vector of the buoy displacements at all the selected time steps. However, the present paper deals with a WEC array in real sea states where, due to the random nature of ocean waves, the incident wave series is non-periodic; this makes the situation more complicated. Ekici and Hall [17] suggested to determine the Fourier coefficients with non-square matrices to replace T inv in (13), where T inv is a block diagonal matrix with each block written as where T , t 1 -t Nt are equally spaced time steps per time window, and N t > 1.5(2N h + 1) was recommended to ensure that the aperiodic equivalents of DFT and iDFT are well-conditioned. Then, T is replaced with the Moore-Penrose inverse of T inv in (13).
By using (13), (11) is equivalent to Note that (15) becomes an overdetermined system when N t > 2N h + 1, which can be solved by using the non-linear least squares method in the Matlab optimization toolbox. The numerical Jacobian matrix is calculated with the same approach as presented in Wei et al. [22] in order to improve the computational efficiency. It is worth remarking that (15) mixes the frequency-domain displacement and the time-domain pumping force. The pumping force is obtained directly by Lagrange multipliers without any smoothing. Thus, no additional treatment is required to deal with the sums and differences of frequencies generated by the non-linear force in the conventional harmonic balance method. This is the major advantage of this approach compared to solving the system in the frequency domain with (11).

Lagrange multipliers
Nacivet et al. [25] introduced a 'dynamic Lagrangian' algorithm for non-linear contact problems. The pumping force is computed based on the dynamic state of the piston, that is sticking, upstroke or downstroke, which is similar to the contact problem. Hence, we applied a similar algorithm in the present study. At each iteration, the pumping force is updated with Lagrange multiplier (corresponding to pumping forces), which is formulated as a penalization of the equations of motion in time domain:̄=F where is a penalty coefficient that can be chosen arbitrarily as positive but influences the convergence speed: a good choice is to set this to half of the hydrostatic force of the pump, = 0.5 g(D − h r );X b,r is a new vector of relative displacements which satisfies the updated dynamic state of the piston. In order to calculatēat each iteration, (16) is separated into two parts, and rewritten as wherēI is determined by the displacement obtained from the non-linear equations (15), and̄I I is the corrective pumping force vector which is computed with a predictioncorrection procedure. At each time step increment, it is firstly assumed that the piston is in a sticking state, that is the piston does not move and with̄0 II = 0. The corrected pumping force vector̄k cor ∈ ℝ N b is determined by enforcing the following rules per buoy: p,up , the piston moves upstroke, 3. WhenF p,up , the piston sticks and does not move,̄k Here the sticking state of the piston is determined by the force rather than piston velocity as the sticking phenomenon cannot be properly captured with the reconstruction of the time series of piston velocity from frequency-domain inversion.

Time window and overlay
Theoretically, the HBM can be straightforwardly used for WEC arrays oscillating in real sea states. Some researchers [15,17] successfully demonstrated the applicability of the HBM for oscillating systems with multiple excitation forces. However, the challenge is that a large number of harmonic components is required to obtain a satisfactory description of the time series of signal data in the frequency domain, which would geometrically increase the computational cost and render the HBM as under-performing compared to a time integration approach.
To overcome this difficulty, we split the time series of the incident wave signal into a series of time windows. For each time window, satisfactory accuracy can be achieved with a finite number of harmonic components and the simulation can be effectively carried out. Since each time window is independent within the HBM framework, the computation can be simply accelerated by using a parallel implementation. This is an appealing feature that can be applied to the development of real time control strategies. The main drawback of using this scheme is that discontinuities may occur when reconstructing the time series based on the results from each time window.
One distinctive feature of hydrodynamic systems is the fluid memory effect. Due to the independence between each time window, the fluid memory effect from the previous time window cannot be transferred to the following one. In order to compensate for the loss of accuracy due to this discontinuity, an overlay between the time windows is included, and only the results in the middle of the time window are used for further analyses such as computing the predicted power output of the OG-WEC system. The energy generated per window E w is calculated by trapezoidal numerical integration of: where t w is the length of time window, t b is the beginning time of the time window, and t o is the length of overlay. From the description above, there are several parameters that influence the simulation results and computational cost, including the length of the time window and overlay and the number of harmonic components. We investigate the influence of these parameters in the following section of this paper.

Comparison with the time domain model
The proposed MFT model is validated by comparing its results to the results obtained from our previously developed time domain (TD) model [26]. The TD model is based on the open source package WEC-Sim, augmented with the in-house developed adaptable piston pump model. The hydrodynamic coefficients were obtained from NEMOH [27], the incident waves were generated based on the JONSWAP spectrum with 1001 wave frequencies, and the configuration of the WEC array was chosen to correspond to previous research in Wei et al. [22]. In the MFT model, the input time series of incident waves was replicated from the TD model using the same pseudo-random phase, and the implemented OG-WEC parameters were chosen to be equivalent to those applied to the TD model, in order to enable the comparison of input and output values of the MFT and TD methods. Figure 2 depicts the wave elevation ( ) and the excitation force (F e ) over the total simulation time for one floater in the OG-WEC array (floater 14 as enumerated in Figure 1). The chosen sea state is a JONSWAP spectrum with a significant wave height of 2m and a peak period of 7s. Figure 2 represents the computation of the MFT model considering a time window length of 19s, an overlap of 7.6s (40% of the time window length) and 9 harmonic components. It can be seen from the first plot in Figure 2 that the input wave elevation is the same in both models. In the second plot in Figure 2, the excitation force of each time window is described with a different colour to show the functionality of the overlap; at the point in which the excitation force of one time window matches the excitation force of the previous one, they connect to create the new starting point for the latter time window's calculation. This ensures a smooth and more accurate representation of the system dynamics. The vertical lines show the time window delimitation within the calculation. The last plot in Figure 2 gives a more detailed representation of the time windows 4 to 7 that show high excitation force in the plot above, and highlights the overlapping technique; the remaining time windows of the simulation show similar accuracy. Large discrepancy in the excitation force can be observed between the MFT and TD models in the first time window which is due to the differences in the initialization. In the TD model, the computation starts in a resting position with a time ramp to avoid strong transient flows at the beginning of the simulation, while, in the MFT model, the simulation is performed in a steady-state manner. The excitation in the following time window shows good agreement with that in the TD model, indicating that the given number of harmonic components is sufficient to approximate the real sea state.
The floater's response is captured through the buoy displacement (X b ), pumping force (F p ) and pumping or predicted extracted power (P p ) (Figure 3). These results are obtained with the given inputs shown in Figure 2. The floater responses by the MFT model are in accordance with those by the TD model to an acceptable degree. The predicted power in each time window shows only slight differences between the two methods except for the first time window's results, which are justified by the excitation force difference explained above. Time windows 5 and 14 show the same absolute power difference; when considering the difference in percentages in the case of 14, the discrepancy appears to be considerably higher.
The pumping force calculated by the MFT model is characterised by less radical peaks and valleys within short time spans. Furthermore, negative values can be observed within the pumping force and pumping power; for comparison purposes, negative power values were set to zero, as this reflects the TD method's approach. Negative power is attributed to the truncation of the Fourier series of displacement results in a Gibbs artifact, manifested by undershoot and overshoot oscillations. Figure 3 depicts the response of a single floater; the predicted power and behaviour of the entire array need to be considered in order to quantify the overall performance of the two models. Figure 4 depicts the predicted power output of each floater in kW (a) and the difference between the predicted power extracted by each floater (b) in the case and settings mentioned above. The incident wave travels from left to right. It can be seen that the difference of predicted power extracted from each floater ranges from a maximum of 9.9% to minimum −6%. The slightly asymmetric results in (b) are attributed to accumulated time integration error in the TD model. The differences are further described and statistically depicted later in the paper ( Figure 6). The total predicted power from the MFT method is 178.61 kW and the array is predicted to produce 177.26 kW based on the TD method; the difference is only approximately 0.7%. Beyond its accuracy, the main advantage of the MFT method is the significant lower CPU time it requires, which will be further analysed in the following sections.

Sensitivity analysis and computational cost
Two sea states were analysed to measure the influence of changes in parameters on simulation time and the difference in predicted power output between the MFT and TD methods with different input conditions and the computational costs related to the simulations (we used a workstation with an Intel Xeon Processor at 3.5 GHz and 64GB of RAM). The first sea state-introduced in the earlier section-has a significant wave height of 2m and a peak period of 7s (Case 1), while the second sea state is described by a significant wave height of 2m and a peak period of 15s (Case 2); this allows us to investigate the influence of the incoming wave characteristics on the choice of parameter settings for the MFT method. In both cases, the time window length and number of harmonic components were altered as necessary; in addition, within Case 1, the overlap length between time windows was analysed. Because of the longer peak period describing Case 2, the total simulation time was chosen to be greater in order to decrease randomness  Table 1 the parameter settings applied in each case are summarised. For all comparisons, the results from the first time window were excluded as a time ramp was applied in the TD model.

Time window length
With increasing time window length, the total simulation time is divided into fewer separately calculated windows, thereby decreasing the computational effort. The advantage of the shorter computation time needs to be weighted against the agreement of the predicted power with the TD method's results.  As can be seen in Figure 5(a,b), choosing the time window to be too small or too large increases the difference in the results. This figure shows the normalised difference of power extraction between the MFT and TD simulations; the fine purple line represents the TD results, while the remaining lines show the MFT results. The time window length has been normalised by the peak period of each case study to highlight their relationship. A considerable disagreement in results in Figure 5(b) can be seen at time window lengths that are shorter than the peak period of the defined sea state. Furthermore, the choice of an adequate time window is influenced by the number of harmonic frequency components defined. With conditions as defined in Table 1, in Case 1, for three harmonic components, a time window of 9s gives a difference in predicted power lower than 5%, while, for six and nine harmonic components, the time windows with 14s and 17-22s show a difference in total predicted power extraction lower than 5%. In Case 2 ( Figure 5(b)), with three, six and nine harmonic frequency components, the time window length ranges to achieve less than 5% difference in total predicted power outcomes between the two methods are 9-27s, 15-35s and 24-45s, respectively. In Case 1, keeping the number of harmonic components constant (e.g. six components) gives a range of computational times between 160s and 10s with the computational time decreasing as the time window length is increased; with the same settings, the TD method's CPU time was 2193s. In Case 2, the CPU time was recorded to be between 634s and 14s with six harmonic frequency components, while the TD method took 6687s to compute its results.

Number of harmonic components
The effect of changing the number of harmonic components can be seen from Figure 5(a,b). Increasing the number of har-monic components makes the fluctuation between the results calculated with different time window lengths smaller; it diminishes the influence of time window length choice on the results. More specifically, the graphs in Figure 5(a) show that the predicted power with six and nine harmonic components stays closer to the TD method's results for an extended range of time window lengths. It can also be seen that the larger the number of harmonic components, the longer the time window length in which they approximate the TD method's results more closely. For more insight into the comparability of results, the predicted power extracted by each floater was considered. It is possible that the total predicted power shows good agreement, although the difference between single floater predicted power is inadmissible. For future research purposes, such as optimizing the energy extraction by adapting the individual PTO, an accurate power prediction per WEC is required. Figure 6 gives a statistical representation of the difference in power extracted by each floater between TD and MFT results; the edges of the box show the 25th and 75th percentile with the line in the middle describing the median. The whiskers of the boxes define the minimum and maximum values excluding outliers that are shown with red crosses. From Figure 6 it becomes clear that even if the solution with three harmonic components has a smaller absolute median and mean value of difference in predicted power, the maximum and minimum values of difference for a single floater's power extraction ranges over a larger span, from 28% (visible as the outlier in Figure 6) to -6%. Increasing the number of harmonic components decreases this range of differences between single floater results considerably (maximum and minimum differences with nine frequency components lie between 9.9% and −6%), while the overall discrepancy is decreased to less than 2%. The result with nine harmonic components showing the least difference is achieved using a time window length of 19s (Case 1); the computation with these parameter settings takes 165s with the same parameter settings; by changing the number of components to three and six the computational effort decreases to 7s and 22s, respectively. Given the considerable decrease in computational effort from 2193s with the TD model, the proposed MFT model can significantly reduce the computational time while obtaining satisfactory results and offering great flexibility to balance computational time and accuracy by choosing an appropriate number of harmonic components.

Time window overlap
As the time window overlap is a function of the time window length, the overlap between adjacent windows increases when longer time window lengths are chosen. In Case 1, two different values for the overlap of time windows were analysed: a time window overlap of 40% of the total time window length and a time window overlap of 20%. Increasing the overlap length can minimize the error introduced by fluid memory effects, but will lead to a larger total number of time windows necessary to cover the total simulation time. Since the results with 40% overlap showed good agreement (as seen above), the results were compared with shorter time window overlaps in order to achieve comparable results with lower computational effort. At long time window lengths, the computational effort decreases up to three times (with nine harmonic frequency components and a time window length of 24s); however, a drawback is that the total predicted power difference with the TD model's results changes within a range of −5.2% to 8.8%.

Optimal settings for the MFT model
The sensitivity analysis clarifies the choice of parameters necessary to tune the MFT model in a way to approximate the TD model's results with satisfactory accuracy and reasonable computational cost. The time window length needs to be chosen carefully depending on the peak period of the sea state defined beforehand and the selected number of harmonic frequency components. Especially in sea states with short peak periods, the correct choice of time window length is crucial (see Figure 5). Harmonic oscillations in the results are noticeable when considering the number of harmonic components: in Case 1, at three components, the time window length with most accuracy is close to the peak period (7s); with six components, it is double (approximately 14s); and, with nine components, the accuracy is greatest at around 21s. For longer peak periods, this behaviour is not as noticeable anymore and other parameters are more influential to the overall accuracy. Increasing the number of harmonic components leads to a longer CPU time, but increases considerably the accuracy of predicted power output per floater, leading to a higher accuracy for detailed analyses.
Considering the significant decrease in CPU time in comparison to the TD model, setting the number of harmonic components to nine still provides a significant advantage in computational effort. Given a random sea state, a rule of thumb is that the time window length can be chosen as twice or three times the peak period and the number of harmonic components should be set as six or nine since, with these settings, we can obtain a good approximation of the incident wave profiles. A constant overlap of 40% of the time window length can be chosen as it balances the advantages of a long overlap for an appropriate memory effect and a shorter overlap to decrease the total CPU time. Furthermore, a long simulation time will decrease the influence of the random incident wave profile on the mean predicted power (Case 2). Tables 2 and 3 summarise the results of the two models' comparison for two case studies, highlighting the significant decrease in CPU time when using the MFT model, while the total predicted power extracted remains accurate. For each case study, three different settings for the number of harmonic components and their corresponding time window length as multiples of the peak period of the case are represented. Although the MFT model can be implemented in parallel, we ran all simulations in series to obtain a fair CPU time comparison. To balance the advantages of shorter CPU time and predicted power agreement, six harmonic components and a time window length of twice the peak period are considered acceptable and generally applicable. These parameter settings will be applied in future work to create accurate power matrices of WEC arrays by significantly decreasing the necessary computational effort to create them.

CONCLUSIONS
The work successfully demonstrates that the MFT method can effectively assess the dynamic response of large WEC arrays with a cluster of highly non-linear pumping units. A series of simulations is carried out to investigate the influence of the key parameters on the accuracy and computational cost.
The following conclusions can be drawn from the present study: • Comparison of results from the TD and MFT models show good agreement in both analysed cases. With an appropriate choice of settings, less than 1% of difference in the overall predicted power output can be achieved. • The total CPU Time of the MFT method is at least tenfold shorter than that of the TD method. • Increasing the number of harmonic frequency components improves the accuracy of the individual WEC results within the array and decreases the dependence of the results' accuracy on the time window length. • The choice of time window length shows significant impact on the accuracy of the predicted power extraction, especially in sea states with shorter peak periods, while increasing the total simulation time decreases the influence of randomness of the incident wave profile. Depending on the number of harmonic frequency components selected, the time window length should be chosen to approximately be a multiple of the peak period length (once, twice or thrice).
Due to the innate nature of the HBM, aliasing and leakage errors of the MFT method cannot be eliminated. Hence, the MFT method shows difficulties in capturing the microscale dynamics, for example the transitions between dynamic states of the piston. The application of the MFT model is limited to describing the global system response. By using the MFT method, extensive simulations of WEC arrays in individual sea states can be carried out at an affordable computational cost. Hence, the MFT method is suitable to be used to assess the power production of WEC arrays by calculating their power matrices, while taking into account non-linear PTO forces. This is a work in progress. In addition, the MFT method can potentially serve as the predictive model in the model predictive control (MPC) based on real-time measurements of incoming waves per time window.