A Fully Coupled Wave-to-Wire Model of an Array of Wave Energy Converters

This paper describes a fully coupled, wave-to-wire time-domain model that can simulate the hydrodynamic, mechanical, and electrical response of an array of wave energy converters. Arrays of any configuration can be simulated to explore both the effects of the array on the electricity network and of network events on the devices within the array. State-space modeling of the hydrodynamic radiation forces enables fast and accurate prediction of the interacting response of multiple devices, including the effects of wave climate, control strategies, and network power flow. Case studies include the demonstration of the bidirectional interaction of the array and the network.


I. INTRODUCTION
C ONNECTING large arrays of wave energy converters (WECs) to the electricity network will increase the variability of the power flow in the network, causing time variation in the voltage profile and potentially an adverse impact on local supply quality. It is therefore necessary to be able to predict this variation in power injection from arrays of WECs and the consequent effects on the electricity network. Additionally, there has been little work reported on the effects of electricity network events on the individual and collective response of the WECs (compared with, say, fault impacts on wind turbines). Insight into how the array control and buoy motions may react to such network perturbations is required for reliable operation of the array. This defines the need for a wave-to-wire model that can simulate the bidirectional hydrodynamic and electromechanical interaction of the devices in an array with one another and with the network. This paper describes the development of such a model. Section II presents a literature review of existing models of WECs and WEC arrays and identifies gaps in the literature that this paper aims to fill. Section III describes the development of the new wave-to-wire model of arrays of WECs, including the hydrodynamic model of the collectors and the modeling of the hydraulic PTOs and the generators. Section IV specifies conditions and parameters for the test-case studies. Section V discusses the results of the time-series modeling and the wholesystem response from the sea to the network and back again. Conclusion from the results is drawn in Section VI.

A. Single WEC Studies
Many wave-to-wire models of single WEC devices have been published and these examine the nature of the primary energy conversion and the effects of different electrical and hydraulic control strategies on the power delivered to the network [1]- [10]. References [1]- [4] describe models of oscillating water column WECs with air turbine generators, overtoppingtype WECs with water turbine generators, WECs with direct-drive rotating generators, and WECs with direct-drive linear generators, respectively. Wave-to-wire models of waveactivated body type WECs are presented in [5]- [10]. Significant understanding has been obtained through these publications about the control of a single WEC and about the electricity network effects of the wave energy conversion process. Unfortunately, these single WEC wave-to-wire models cannot be duplicated to model WEC arrays because the hydrodynamic interactions between the devices (the radiated and diffracted waves) are not taken into account.

B. Arrays of WECs for Hydrodynamic Control and Array Layout Studies
Models of arrays of WECs, with differing levels of complexity, have been developed to examine different areas of interest [11]- [25]. The effects of hydrodynamic interactions, control in arrays, and the layout of arrays have been analyzed using frequency-domain hydrodynamic models in [11]- [19]. Frequency-domain simulations are extremely fast and they are very useful for initial studies; however, they are restricted to linear problems and can only incorporate linear external forces such as linear PTO forces.
Time-domain models have been used to simulate the response of WEC arrays in [20]- [25]. In [20], a method was presented which can be used to investigate the optimal control of WEC arrays. In that approach, the WEC motions and PTO forces are expanded in a set of suitable basis functions and then the time-domain equations of motion are solved using a Galerkin method. In [21], the time-domain equations of motion are solved for an array of WECs but with just linear PTO forces. This work is licensed under a Creative Commons Attribution 3.0 License. For more information, see http://creativecommons.org/licenses/by/3.0/ The solution procedure does not use a state-space approach but instead is based on the frequency-dependent transfer functions which link the motion of each device to the free-surface elevation at the origin. An interesting aspect of the presented method is the way in which it solves the problem caused by the highfrequency oscillations which occur in these transfer functions when a device is located far from the origin.
The commercial hydrodynamics software WaveDyn incorporates a time-domain array code, which can also model all the linear hydrodynamic interactions between all the WECs [22]. WaveDyn directly computes the radiation convolution terms necessary to calculate the radiation forces but for large arrays, this approach is computationally inefficient.
The hydrodynamic time-domain array model used in this paper was also used by Nambiar et al. [23] to investigate the coordinated/global control of an array of closely spaced point absorbers. However, unlike in this paper, only a brief description of the time-domain model was presented in [23].
All the above publications either used simplified representations of the electrical generator or used the power extracted by the array as the generated power. None modeled the power takeoff (PTO) system and the generator in detail. The hypothetical electricity network to which the arrays were connected was assumed to have infinite capacity and consequently would not exhibit any supply quality variations or represent the effects of constraints in weak networks.
Two time-domain wave-to-wire models of arrays are presented in [24] and [25]. Reference [24] suggests how direct drive WECs may be arranged in arrays for enhanced power quality and the development of an advanced hydraulic PTO system to optimally control an array of WECs is described in [25]. Both these papers incorporated simplified models of the hydrodynamic interactions between the WECs. In [24], a regular wave power correction factor was used to obtain the power output of each WEC in the array based on the power output of that WEC in isolation. Superposition was then used to obtain the total power from the array in irregular seas. While this approach is sufficient for linear PTO forces, the method reported herein can also take into account nonlinear PTO forces. In [25], the self-induced radiation forces on the WECs are modeled, but the radiation forces between WECs are not.

C. Arrays of WECs for Electricity Network Interaction Studies
In the past few years, several publications have examined the impacts of single WECs and arrays of WECs on the electricity network [26]- [30]. In these publications, the focus was on power quality aspects of the electricity network interaction. An evaluation tool to assess voltage flicker due to a single grid-connected WEC was presented in [26]. Time series of the predicted power generated by the WEC were used to determine the flicker severity indices for all the sea states seen at the European Marine Energy Centre test site.
This work with a single WEC was extended to include arrays of WECs in [27]- [29]. Flicker and voltage fluctuations at the point of common coupling (PCC) caused by a 20-MW array of WECs connected to electricity networks of different strengths were studied in [27]. Effects of wave power injection at the Pacific Marine Energy Centre test site were discussed in [28]. In that paper, flicker and voltage fluctuations were examined as well as the harmonics introduced by arrays of WECs and the low-voltage ride-through capability of those arrays. In [29], a techno-economic analysis of different wave farm collection networks was presented. In addition to flicker and voltage fluctuations, the losses in the collection network and the power factor at the PCC were considered in the technical assessment. Additionally, the effect of having energy storage in each WEC on the flicker level was also explored.
In the last three publications, array hydrodynamics were not modeled. The power time series obtained from a single WEC, either from measurements or from simulations, were used, along with appropriate time delays, to estimate the power time series of all the WECs in the array. Uncertainty in the assumed wave propagation characteristics was included by adding a random time-delay component.
A wave-to-wire model of a point absorber WEC was cloned to build a wave-to-wire array model in [30]. Three control levels were investigated in this paper: wave-side control (hydrodynamic control of the WEC), grid-side control, and centralized farm control. The last two of these control strategies dealt with real and reactive power control for reducing voltage fluctuations. Short-term storage was included in the model as a low-pass filter. Mutual interactions between the three main subsystems in the wave-to-wire model of the WEC (the hydrodynamic module, the PTO (including generator) module, and the electricity network/grid module) were included in the full model.
A methodology to generate high-resolution time series of the power generated by large WEC farms, over long durations, for power system studies was demonstrated in [31]. In that work, the WEC displacement was assumed to be the modeled incident wave elevation at the location of each WEC.
In all the above publications, the focus was on the power system effects of connecting wave power farms and simplified models of the array hydrodynamics were used. This was done with the aim of reducing the computational intensity of simulating large farms and was thought to be fit for purpose for those studies.

III. NEW MODEL DESCRIPTION
From the literature review presented in the preceding section, it becomes clear that there has been little work published on the controlled responses of devices within arrays connected to constrained electricity networks. Further, little has been reported on bidirectional (wave-to-wire and concurrently wire-to-wave) models of hydrodynamically interacting arrays of individually controlled WECs. This section describes the development of a generic, fully coupled, bidirectional, time-domain wave-towire model for arrays of interacting single-bodied, and sea-bed referenced WECs. The model is generic in that the different subsystems in the model can be changed (e.g., type and hull form of the buoys, the generators, and PTO systems). In order to have realistic models of the PTO systems and the generators, it is necessary to allow nonlinear PTO forces. If nonlinear PTO forces are included, then the hydrodynamic model must be a time-domain model.
Including all the hydrodynamic interactions (all the radiated and diffracted waves) is critical to accurately model the individual responses and the power output from the array. This means that all the radiation convolution terms must be included. Since computing these terms is computationally inefficient, it is better to replace them with an equivalent state-space model. However, unlike for a single device problem, it is not straightforward to obtain a numerically stable state-space model for an array of devices. Section III-A describes the development of a new hydrodynamic time-domain model with radiation state-space models which are designed to have increased stability.
The PTO, generator and the electricity network subsystems within the model can simulate dynamic or transient electromechanical events. The whole model allows the study of the controlled responses of devices within arrays connected to constrained electricity networks and also the concurrent effects of the network power flow and transient events on array control and in turn buoy motions.
The primary input to the model, shown schematically for a single WEC in Fig. 1, is a time-series of simulated or measured wave elevation. The wave excitation force f (t) is obtained from the wave elevation and drives the hydrodynamic model: a time-domain equation solver for the motion of the WEC. The PTO could be a linear/rotating generator-based direct-drive system or a hydraulic system driving a rotating generator. It applies a resistive force f PTO to the collector body that is controlled, based on current conditions of the body [e.g.,ẋ(t)], to meet real power and/or speed control objectives. For systems with hydraulic PTOs, as modeled in this paper, the controlled PTO unit determines the torque T m (t) on the generator shaft, while the generator sets the rotational speed. The generator is connected to the electricity network either directly or through power converters. Depending on this choice, different voltage and/or reactive power control strategies can be applied. Sections III-A-III-C describe the three main subsystems of the wave-to-wire model presented here.

A. Hydrodynamic Model
The hydrodynamic time-domain model was designed specifically for arrays and it takes into account all of the hydrodynamic interactions, in terms of all the diffracted and radiated waves, between all of the WECs. The commercial frequencydomain code WAMIT is used once to generate the hydrodynamic input data, where each device can have any shape and the array can have any configuration, with each rigid body WEC moving in up to six degrees of freedom (DoF). The array size that can be efficiently modeled has been increased by an order of magnitude to several tens of WECs as a result of the fast and efficient state-space techniques developed.
The hydrodynamic time-domain equations of motion for an array of N freely floating rigid bodies are These linear equations were first proposed for a single body by Cummins [32]. In (1), x i (t) is the motion of mode (DoF) i. Modes 6(n − 1) + 1, 6(n − 1) + 2 to 6(n − 1) + 6 = 6n are, respectively, the surge, sway, heave, roll, pitch, and yaw modes of body n. Also in (1), , and C = [C ij ] are, respectively, the global mass matrix, the added mass matrix at infinite frequency, the matrix of radiation impulse response functions, and the matrix of hydrostatic and gravitational restoring coefficients. WECs that are constrained not to move in certain modes can be represented in a reduced system by removing the associated rows and columns in (1). Each radiation impulse response function k ij (t) is the inverse Fourier transform of the frequency-dependent radiation impedance function K ij (ω), which is defined as where A ij (ω) and B ij (ω) are the added mass and added damping coefficients. The global mass matrix M is obtained from the geometry and the mass distribution of the WECs. In preparation for use, the matrices A(ω), A(∞), B(ω), and C are calculated by WAMIT for a range of finite frequencies and the infinite frequency case. WAMIT also outputs wave excitation forces for each incident sinusoidal wave and superposition is used to recreate irregular seas.
The numerical solution of (1) is complicated by the integral radiation convolution terms. Although it is possible to compute them directly, this is computationally inefficient and so they are often replaced. Replacement usually involves system identification and this can be performed either in the time or frequency domain [33].
In this work, system identification is performed in the frequency domain. Each radiation impedance function K ij (ω) is approximated by a rational transfer function with numerator and denominator polynomials. For specified orders n and m of the numerator and denominator, respectively, a first approximation is obtained by a least-squares error fitting method [34]. This approximation is then used as the initial condition for a second algorithm which uses a damped Gauss-Newton iterative search method [35]. This additional algorithm guarantees the stability of the resulting transfer function, i.e., its poles (zeros of the denominator polynomial) are all in the left half-plane. The implementation of the above process is carried out using the invfreqs method available in MATLAB. The orders n and m are then systematically increased, with n < m, until the relative root-mean-square error between the transfer function and K ij (ω) is less than a prescribed value (<1%). Once the whole system of approximating transfer functions is obtained (i.e., for all i and j), the system is converted to a single equivalent state-space model. If none of the modes of the WECs are constrained, then this system will consist of (6N ) 2 transfer functions, although approximately half of them will be the same due to the symmetric nature of K ij (ω) [i.e., K ji (ω) = K ij (ω)]. A general state-space model has the forṁ x ss = Ax ss + Bu ss y ss = Cx ss + Du ss (3) where u ss , x ss , and y ss are the input, state, and output vectors and A, B, C, and D are the state, input, output, and feedthrough matrices, respectively. For the modeling in this paper, D will always be the zero matrix. The advantages, in terms of computational resource, of converting the radiation convolution terms to a state-space model are discussed in detail in [36]. The improvement in the computational efficiency is because of the Markovian property of state-space models, where the states of the model summarize all the past information [33].
One aim of this work was to derive a method which would obtain radiation state-space models for general arrays of wave energy devices. For each set of radiation convolution terms, there is an infinite number of equivalent state-space models and most of these will be numerically unstable, i.e., although a state-space formulation may be constructed to be mathematically stable, errors may arise with the numerical computation of the state trajectory using finite-length numbers. A multipleinput multiple-output (MIMO) system can be represented by a matrix of transfer functions and has a characteristic equation formed from the product of the denominators of the transfer functions. It is not uncommon to encounter number overflow and truncation to finite length when large-value coefficients are multiplied together. This truncation of the coefficients changes the roots of the characteristic equation (the eigenvalues of the system matrix) and can lead to computational instability. The objective of this work was to design a method which would produce numerically stable radiation state-space models.
State-space models can be constructed from transfer functions via any of the standard canonical realizations. Depending on the system being modeled, a particular canonical form results in numerically stable or unstable state-space systems. For example, here, using a modified version of the controllable canonical form tends to result in computationally unstable models through the number overflow and truncation problems mentioned above, i.e., when a large system of possibly high-order transfer functions, each of which is expressed in terms of expanded polynomials, is converted in this way to a MIMO state-space system, the resulting state matrix and its characteristic equation can contain a very wide range of elements/coefficients from the extremely small to the extremely large. Due to the finite precision of computer arithmetic and the extreme sensitivity of the roots of a polynomial to its coefficients, this can result in changes to the roots of the characteristic equation (eigenvalues of the state matrix) and a numerically unstable state-space system. A far more robust approach, and the one adopted here, is to first convert each transfer function to zero-pole-gain form by factorizing its numerator and denominator polynomials. The resulting system of zero-pole-gain models is then converted to a state-space system, which this time is in what is called modal canonical form and is found to be far more numerically stable. The MATLAB function ss is used to convert the system of transfer functions in the zero-pole-gain form to a single state-space model.
Another advantage of replacing the radiation convolution terms with a single state-space model is automation. Fig. 2 shows the simulation block diagram used for the solution process of the equations of motion. This block diagram describes the procedure for solving the following rearranged version of the matrix form of (1), where the radiation convolution terms have been replaced by the state-space system Since all the radiation forces are now calculated from the single state-space block, the single block diagram in Fig. 2 can be used to solve the hydrodynamic time-domain equations for any array, irrespective of how many DoFs there are. For a multiple device, multiple DoF problem, the (signal) lines connecting the blocks will be vector (signal) lines. The Simulink/MATLAB modeling environment was used for the full wave-to-wire model and the Simulink diagram/model for the hydrodynamics module is virtually identical to the block diagram in Fig. 2.

1) Nonlinear Hydrostatic Forces:
The hydrostatic and gravitational forces acting on a WEC are represented, in a linear manner, by the terms C ij x j (t) in (1). This means that in heave motion, the water plane area is assumed to be constant as the WEC oscillates and there is no accounting for the motion of the free surface around the WEC.
In this work, an array of heaving buoy WECs was used to demonstrate the functionality of the numerical model. The buoys have hemispherical bottoms, of radius r = 3 m, and cylindrical tops. In equilibrium conditions, the hemispherical part of each buoy will be below the free surface and the cylindrical part will be above [ Fig. 3(a)].
The linear hydrostatic and gravitational terms C ij x j (t) in (1) are replaced by more accurate nonlinear representations, as described below. Fig. 3(b) shows the first heaving buoy, whose vertical motion is described by x 3 (t). The free-surface elevation is η, which is approximated to be constant in space around the WEC and equal to the instantaneous free-surface elevation the incident wave would have at the center of the WEC. S B is the wetted surface of the buoy and S WP is the waterplane surface whose area is A WP . Hence, a closed surface S can be formed as S = S B + S WP . The hydrostatic force on this WEC becomes (5) wheren is the unit normal pointing out of S. By simplifying this equation, applying a form of Gauss's divergence theorem to the integral over S and adding the gravitational force on the WEC, the total nonlinear hydrostatic and gravitational force on the buoy becomes where V is the volume inside the closed surface S andk is the unit vector in the z-direction. Equation (7) has three different cases: 1) x 3 − η < 0 (hemispherical part completely submerged); 2) 0 < x 3 − η < r (hemispherical part partially submerged); 3) x 3 − η > r (buoy out of the water). The introduction of this nonlinear stiffness was found not to affect the WEC motion significantly for the WEC geometry used here. This need not always be the case for other WEC types. For example, when heaving WECs are highly flared at the waterline, the waterplane area will change significantly with x 3 − η. Also, when heaving WECs are excited by sea states with peak frequencies close to their natural frequency, displacements bigger than the wave elevation may be seen that increase x 3 − η. In such cases, the nonlinear hydrostatic effects will become more significant.

B. PTO Model and Control
The relative motion between the collector bodies of the WECs and their fixed reference is converted into mechanical, hydraulic, and then electrical power by the PTO unit. The PTO unit applies a resistive force f PTO to the WEC collector body, which modifies (1) to The PTO force acting on a WEC can be controlled to meet different objectives (e.g., maximize power extraction and smooth cumulative power output). The advantages of hydraulic PTO systems-their robustness, speed control characteristics, opportunity for energy storage, and relatively lower costs-led to their inclusion in this model.
The PTO system modeled here (Fig. 4) is based on the hydraulic systems described in [37]- [39]. The motion of the collector body relative to a fixed reference point (such as the sea bed) drives a hydraulic ram, which forces the flow of highpressure (HP) oil through a hydraulic circuit. The hydraulic circuit consists of a HP accumulator for smoothing, a lowpressure (LP) accumulator to avoid cavitation, and a variable displacement oil motor. Both accumulators are Nitrogen gas charged. A rectifying valve system ensures that the fluid flow through the accumulator-motor set is unidirectional. The pressure difference between the HP and LP accumulators causes fluid flow through the hydraulic motor. Fig. 4 also shows schematically the coupling between the shaft of the hydraulic motor and the rotor of the generator and, finally, the connection of its stator to the electricity network.
The volume flow rate of the fluid from the piston to the accumulator Q p is given by A p dx/dt, where A p is the crosssectional area of the piston and x is the WEC displacement. Inertia effects and oil pressure losses in the hydraulic circuit are small enough to be omitted in this model, which means that the pressure difference between the lower and the upper parts of the cylinder is equal to that between the HP and LP accumulators. Assuming that the gas compression and expansion process in the accumulator is isentropic and the gas is ideal, the accumulators are modeled by where ΔQ(t) = Q p (t) − Q m (t) and ΔQ(t) = Q m (t) − Q p (t) for the HP and LP accumulators, respectively. In (9), V f (t) is the volume of the fluid within the accumulator and Q m is the flow through the hydraulic motor. The volume of gas within the accumulator is given by where V 0 is the accumulator capacity. It lies between V pre (maximum limit) and V pre (p pre /p rv ) 1/kg , where V pre is the precharged gas volume, p pre is the precharged pressure, p rv is the relief valve pressure, and k g is the specific heat ratio of Nitrogen.
The PTO force f PTO , in (8), is given by where p high (t) and p low (t) represent the pressures of the HP and the LP accumulators, respectively. These accumulator system pressures p acc are obtained using where p acc (t) = p high (t) and p acc (t) = p low (t) for the HP and LP accumulators, respectively. The oil flow rate Q m drives the variable displacement hydraulic motor that is coupled to the generator. The flow rate through the motor is given by Q m = DV m ω m where D, V m , and ω m are, respectively, the instantaneous displacement (in per-unit relative to the maximum displacement), the volume capacity of the hydraulic motor, and the rotational speed of the motor.
The loadings on the bladder that separates the fluid chamber from the gas chamber in the accumulator and the loadings on the motor shaft due to friction are low enough to be omitted from the model. Fluid compressibility effects and leakage flow in the motor have also been neglected. Fig. 5 shows schematically the hydraulic motor control strategy. The time variation in the gas pressure difference (p high − p low ) averaged over an interval of time (here 30 s) is used to determine the reference torque T ref through a linear relationship. This approach helps to smooth the torque at the motor shaft. A PID controller determines the motor displacement D required to maintain the motor torque T m at the reference value. The torque at the motor output shaft T m is given by The WECs across the array, controlled with this torque regulation algorithm, were found to extract on average the same total amount of power as with the optimal damping force (with force constraints) applied to individual devices in the array.

C. Induction Generator Model
The opportunity for torque control provided by variable displacement hydraulic PTO systems, described in the preceding section, allows the use of three-phase cage induction generators that are directly connected to the electricity network. Such a system has been modeled here. This system provides both flexibility and controllability through the use of the variable displacement hydraulic PTO. This means that the use of expensive, fully rated power converters between the generator and the electricity network can be avoided. With no power converters between the generator and the network, the generator can contribute to power system damping and also does not add any harmonics into the network. The Pelamis device used fixed-speed induction generators with its digital hydraulic PTO system [40]. The use of fixed-speed induction generators, coupled with hydraulic PTO systems, has been discussed as an option for WECs in [41] and [42]. Note that models of any generator system can be incorporated into this generic wave-towire model and the induction generator was chosen here only for demonstrating the functioning of the bidirectional model.
The dynamics of the cage induction generator were modeled in Simulink using the established state equations [43] V qs = r s i qs + ωλ ds + d(λ qs )/dt where i, V , and λ represent current, voltage, and flux, subscripts d, q, and 0 refer, respectively, to the direct, quadrature, and zero axes components, subscripts s and r are used, respectively, for stator and rotor quantities, r is the resistance, ω is the angular speed of the reference frame, and ω r is the angular speed of the rotor. For cage induction generators, V dr and V qr are zero. The flux linkage expressions in terms of current are where L m is the mutual inductance, and L s and L r are the stator and rotor leakage inductances, respectively. The electrical torque of the machine is given by The speed ω r is related to the developed electrical and mechanical torques as where T e is the electrical torque, T m is the load torque or the torque at the hydraulic motor shaft, J is the rotor moment of inertia, and P is the number of poles in the stator winding.

D. Electricity Network Model
A power flow solver written in MATLAB was used to model three-phase balanced electricity networks. The currents i ds and i qs from the generators in the array were converted to stator real power P s and reactive power Q s using and The stator real and reactive power were inputs to the power flow solver. The network layout and demand data and the power generation profiles of all the other generators connected to the network were made available to the solver prior to running the simulation. The power flow solver feeds back the computed direct and quadrature axis voltages to the induction machines. The power flow solver captures the effects of the variable input power on the network and also translates network events back to the array. This approach, which was used to model the electricity network, allows the simulation of power flow in complex networks and also allows the investigation of demand and network transients.

A. WEC and the Array Configuration
The single-body WECs modeled as an array here are semisubmerged, vertical, hemispherically ended buoys. They are constrained to move only in heave under the action of the wave excitation forces. Fig. 6 shows the layout of the array of four WECs, which, in this instance, is similar to the Wavestar device [23]. The array was modeled at a site with a mean water depth of 20 m.

B. Details of the Derived Hydrodynamic Time-Domain Array Model
The details of the derived hydrodynamic time-domain model are described in the following. WAMIT was run for the 80 equally spaced frequencies ω = 0.05, 0.10, . . . , 4.00 rad/s. Table I shows the orders of the numerator and denominator polynomials of the transfer functions which approximate the radiation impedance functions K ij (ω).
Here, i and j refer to WECs 1 to 4, where WEC 1 is nearest to the incoming waves (see Fig. 6). Recall that the algorithm automatically obtains each approximating transfer function by systematically increasing the order of the numerator and denominator polynomials until the relative root-meansquare error between the transfer function and its associated K ij (ω) is less than a prescribed value (<1%). Only the unique radiation impedance functions K ij (ω) are listed in Table I. Due to the symmetry of the present array and the symmetric nature of K ij (ω) [i.e., K ij (ω) = K ji (ω)], all the other radiation impedances (all the other combinations of i and j) are each equal to one of those listed. This system of approximating transfer functions is then converted to a single equivalent state-space model, as described previously. The resulting state matrix A, input matrix B, and output matrix C are 168 × 168, 168 × 4, and 4 × 168 matrices, respectively.

C. Wave Resource
The wave elevation input was obtained from linear spectral representations of the sea. The two-parameter (significant wave height H s and mean zero-crossing wave period T z ) Bretschneider spectrum [44] was used. An irregular, multidirectional sea state with H s = 1.55 m and T z = 5.98 s, which is a commonly occurring sea state for the 20 m depth considered, was used in all the simulations.

D. Power Takeoff
The PTO system for each WEC consists of a piston-driving HP oil through a hydraulic circuit (see Fig. 4). The piston has a cross-sectional area of 0.0173 m 2 . The HP and LP accumulators have capacities of 1 and 0.5 m 3 , respectively. The precharge and relief valve pressures of the HP accumulator are 30 and 400 bar, respectively, while those of the LP accumulator are 10 and 20 bar, respectively. The hydraulic motor has a volume of 312 cm 3 .

E. Generator and Electricity Network
The hydraulic motor is coupled to a three-phase, fourpole, cage induction generator rated at 400 V and 110 kVA. The array is connected to the electricity network shown in Fig. 7 (all impedances are on a 100-MVA base). Each generator in the array is connected to an offshore substation through 30-m-long subsea cables. The offshore substation has a 500-kVA 0.4/11 kV transformer. The power is transmitted onshore through a 10-km-long subsea cable. The onshore network comprises a 10-km overhead line connected to the infinite grid. A fixed 300-kVAr reactive power compensator is available at bus 5. Onshore local demand (500 kVA at 0.95 lagging power factor) is supplied through a 500-kVA 11/0.4 kV transformer at bus 2, having a fixed tap ratio of 0.93. The electricity network modeled here is a generic representation of the offshore and onshore assets, including local demand.

V. RESULTS AND DISCUSSION
This section presents and discusses the results of simulated bidirectional power flow in the array of WECs, in the normal state and during electricity network events.

A. Normal Operation
The normal operation of the array is described here using a 100-s time series of the different variables associated with the first WEC (nearest to the incoming waves, Fig. 6). Fig. 8(c) shows the torque at the induction generator shaft. The torque is regulated by the hydraulic motor control strategy discussed earlier.
Between time t = 200 s and time t = 214 s, the mean absolute excitation force is relatively low [ Fig. 8(d)]. Therefore, to regulate the torque at the motor shaft, the hydraulic motor displacement (in per-unit relative to the maximum displacement) and flow rate through the hydraulic motor increase [ Fig. 8 Between time t = 214 s and time t = 230 s, the mean absolute excitation force is higher (when compared to the preceding period). To regulate the motor shaft torque, the hydraulic motor displacement reduces, which slightly reduces the fluid flow through the motor. The PTO force no longer balances the excitation plus the hydrostatic (and gravitational) force and the WEC moves. Thus, there is fluid flow from the ram, which at times is higher than the flow required through the motor. The excess fluid from the ram is stored in the HP accumulator, thereby reducing the gas volume in the HP accumulator and increasing the HP accumulator pressure. The PTO force appears to be capped in this interval but this is because its magnitude is limited to the differential pressure between the HP and LP accumulators multiplied by the area of the piston.

B. Hydrodynamic Coupling and Electromechanical Connection
Here, the array response of the four WECs to a three-phase balanced fault at the terminals of the first WEC's generator is analyzed, to demonstrate both the electromechanical and hydrodynamic coupling in the model. In this case, there is no electrical interactions between the four generators because they are connected through an equal impedance Z 1 to bus 5 in Fig. 7, which is taken to be the infinite grid. This isolates the hydrodynamic coupling between the four devices from the electrical effects of the fault. The closeup fault was chosen here because it represents the greatest electromechanical perturbation possible to the generator in the WEC, which is expected to radiate out waves to the neighboring buoy. The fault occurs at time t = 175 s and is cleared in 10 s. The fault duration was exaggerated to explore the model's functionality in terms of its effects on the PTO system and the WEC motions.
The generator terminal voltage during the fault is close to zero. The reactive power required by the generator for magnetization, which was obtained from the network prior to the fault, is no longer available due to the voltage collapse. The generator At time t = 176.5 s, the generator reaches the overspeed limit of 20% above synchronous speed. In-built over speed protection control in the hydraulic motor control system operates to prevent any further overspeeding of the generator and possible mechanical damage. This is accomplished by rapidly reducing and then controlling the flow through the hydraulic motor, by control of the motor displacement [ Fig. 9(d)]. This in turn keeps the power [ Fig. 9(b)] and torque at the motor shaft close to zero [ Fig. 9(c)], just sufficient to maintain a constant generator speed (188 rad/s here). Once the fault is cleared at time t = 185 s, the system returns to normal operation.
From time t = 175 s, when the fault occurs, the displacement of the first WEC with the fault at its generator terminal is different when compared to the normal operation case [see Fig. 10(a)]. The displacement difference is smaller than expected, even for a generator terminal fault that is the biggest perturbation possible in this system. This is because of the inherent inertia and damping terms within the WEC system. This simulation shows how electromechanical events affect  The difference in the motion of the neighboring WEC, between fault and normal conditions, is an order of magnitude lower when compared to the first WEC [see Fig. 9(b)]. This difference is purely because of the hydrodynamic interactions included in the model, since the WECs are not coupled electrically in this case. This simulation thus also illustrates the hydrodynamic interactions between the WECs.

C. Coupling Between Individual Generators and the Electricity Network
To demonstrate the coupling between the individual generators and with the wider electricity network, the performance of the array is analyzed when sudden disturbances are introduced.  Fig. 11(a) shows a step voltage increase of approximately 4% at the offshore bus 5 brought about by the sudden loss of load at bus 3. The generator speed oscillates at the moment the loss of load occurs and finally settles at a marginally lower speed (157.31 rad/s down to 157.28 rad/s).
To illustrate the coupling between the generators, WECs 1 and 2 in the array are forced, by by-passing the PTO control, to not generate any real power after a specific instant of time. Fig. 11(b) shows the 0.5% reduction in the voltage at bus 5 when these two WECs stop generating. Their rotational speeds drop and they start motoring. The rotational speeds of the two WECs (3 and 4) still generating oscillate due to the voltage step change and then stabilize. The generator speed change after the transients die out is not significant owing to the relatively small change in the voltage seen by the machines still generating.
The mechanical effects of energy exchange between the generators through the network interconnection are visibly demonstrated to be of the form expected, but less severe than envisaged, due to the relatively high energy storage in the inertia of the original and hydrodynamic added masses and, to a lesser extent the on-board storage. The simulations do, however, illustrate the sensitivity and completely coupled response of the system model and they also demonstrate the interaction between the generators in the array.

D. Effects of Integration of Wave Power to Electricity Networks
This section illustrates some of the modeled effects of injecting wave array output power into the edges of a distribution network. As an example, the effects of the varying power generated by the array on the voltage at the load bus (bus 3 in Fig. 7) is examined. The array was excited by a more energetic sea state (with H s = 4.75 m and T z = 14.5 s) to increase the power generated by the array, which in turn has a bigger influence on the network. The power generation profile of the array and the voltage at the load bus are shown in Fig. 12. The voltage profile resembles the power profile, as expected in weak networks. In addition to the average voltage rising within the statutory limits, the cyclic power production can introduce voltage flicker, as shown using the simplified hydrodynamic array models in [27]- [30]. These and other effects on supply quality can now be explored using the developed model.

VI. CONCLUSION
This paper has described a fully coupled, bidirectional, waveto-wire model able to predict the hydrodynamic, mechanical, and electrical response of an array of WECs of any configuration or collector form in complex mixed seas. A state-space model of the radiation forces has been used to accurately and efficiently simulate the hydrodynamic response of the devices in the array. Device-device hydrodynamic interactions are fully embodied across the array. The normal operation of an example array of four WECs and the influence of correctly modeling the hydrodynamic interaction between the devices were initially examined. Then, the fully coupled nature of the model was demonstrated using two sample cases: one involving the response of the wave-to-wire system during a fault close to the terminals of one generator and the second where the effects of sudden disturbances in the electricity network on the array were studied. The results allow the identification of the impact of individual and collective power production on the network power quality, as well as the investigation of the effects of network contingencies on individual buoy responses.