Applied The application of the spectral domain modeling to the power take-off sizing of heaving wave energy converters

The power take-off (PTO) system is a main component in wave energy converters (WECs), and it accounts for a notable proportion in the total cost. Sizing the PTO capacity has been proven to be significant to the cost-effectiveness of WECs. In the numerical modeling, the PTO size is normally represented by a force constraint. Therefore, to accurately evaluate the power performance of WECs with various PTO sizes, it is necessary to take the PTO force limitation, a nonlinear effect, into consideration. In this paper, a computationally-efficient spectral domain model of the PTO force saturation is developed for a heaving point absorber, and the nonlinear term is included by statistical linearization. For comparison, a frequency domain and nonlinear time domain model are implemented, and the developed spectral model is verified with the results of the nonlinear time domain model. Compared with the frequency domain model, the spectral domain model remarkably reduces the relative errors in predicting the power performance of WECs with force constraints, while the computational demand is much lower than the nonlinear time domain model. Furthermore, a case study is conducted to size the PTO capacity for reducing the levelized cost of energy (LCOE) in a chosen wave site. Three different numerical models are applied respectively. The frequency domain model could lead to a misestimate of the optimal PTO capacity, with a maximum relative error on the prediction of the annual energy production (AEP) of 24%. In contrast, the spectral domain model indicates the same optimal PTO size with the time domain modeling, and its relative errors on the prediction of the AEP are within 4.3%.


Introduction
Ocean waves contain a huge amount of clean energy, which is attractive to renewable energy communities. Wave energy has been researched over decades, and a number of wave energy converters (WECs) have been proposed and even tested (Aderinto and Li, 2018;Lehmann et al., 2017). However, its levelized cost of energy (LCOE) is much higher than other competitive renewable technologies (De Andres et al., 2017). This is a big hurdle to the large-scale commercialization of WECs.
Recently, scientific literature has investigated the impact of the power take-off (PTO) system sizing on the economic performance of WECs. The PTO system is actually a core component of WECs. On the one hand, the PTO converts the absorbed mechanical energy to usable electricity, and it is therefore highly related to the power performance. On the other hand, it is an expensive system and its cost normally accounts for over 20% of the total capital expenditures (CAPEX) (Tokat, 2018). In Tan et al. (2021), the influence of PTO sizing on the LCOE was studied, and three realistic wave sites were considered. The results showed that the PTO sizing is of significance to the LCOE, and suitably * Corresponding author.
downsizing the PTO capacity would penalize the annual energy production (AEP) but improve the overall economic performance. In addition, a series of studies have also indicated that the PTO sizing plays an important role in the techno-economic performance of WECs (Tokat and Thiringer, 2018;Tan et al., 2020;Tai et al., 2012). Hence, it is necessary to include the PTO sizing in the early design stage of WEC concepts.
The PTO size implicitly indicates the maximum force it could sustain during normal operation. Therefore, it has to be considered as a physical constraint in the numerical modeling when calculating the power performance of WECs with differently sized PTO systems. In mild wave states with sufficiently large PTO sizes, the violation of the force constraint is limited, and the dynamic response of the WEC can be thought linear. For linear systems, the traditional frequency domain (FD) modeling is considered to be a suitable numerical tool. When the force constraint comes to play, the force saturation becomes relevant and the system then is nonlinear. In this case, the time domain (TD) modeling is required to include nonlinear behavior. However, the TD modeling is much more computationally time-demanding than the FD https://doi.org/10.1016/j.apor.2022.103110 Received 3 October 2021; Received in revised form 11 January 2022; Accepted 15 February 2022 models. In this sense, the TD modeling is not an appealing solution since the PTO size optimization requires a large number of iterations through different wave states. To deal with the force constraint and save the computational efforts simultaneously, some assumptions were made in previous studies. For instance, a PTO sizing method based on the FD modeling was presented in Tan et al. (2021), in which the irregular wave states were transferred to harmonic waves by equaling the wave power transport. Then, the PTO parameters were tuned to make the PTO force amplitude in the corresponding harmonic waves lower than the defined constraint. In Backer (2009), a hydrodynamic optimization of a point absorber was conducted subjected to different physical constraints based on the FD modeling. The constraints were handled based on probabilistic analysis. The probability of exceeding the constraints was calculated by assuming that the responses obey a Rayleigh distribution. In the formulation of the Rayleigh distribution, the probability can be related to the standard deviation of the considered response. Therefore, given a defined tolerance on the exceedance probability, the PTO parameters were optimized to maximize the extracted power within the allowed range of the standard deviation of the response. Although these two methods ease the FD modeling's problem of lacking the consideration of the force constraint, they are still insufficient to model the contribution of the constraints to the power absorption. As an alternative, the spectral domain (SD) modeling can be used to take nonlinear components into account, while its computational demands are significantly less than the TD modeling. The SD modeling could be regarded as an extension of the FD modeling, and the nonlinear terms are incorporated by statistical linearization. Therefore, the SD modeling is not able to predict the transient response, such as the peak displacement and extreme loads. However, it could effectively provide the estimates of statistical parameters, such as the mean absorbed power, which makes it highly suitable for the early-stage conceptual design and performance evaluation of WECs.
Although statistical linearization has been widely used in dynamic engineering for long time, the SD modeling was only introduced into the field of wave energy in 2010 (Folley and Whittaker, 2010). In the work, a SD model considering the quadratic damping and wave force decoupling was derived for a flap-type WEC. Afterwards, follow-up studies were conducted to develop the SD modeling for taking more types of nonlinear components into account. In Folley and Whittaker (2013a), a SD model including nonlinear damping terms for oscillating water column WECs was developed and validated by physical wave tank tests. In Gunawardane et al. (2017), a SD model of nonlinear hydrostatic stiffness was established for a semi-sphere point absorber. Recently, in Silva (2019), da Silva et al. (2020), Silva et al. (2020) and Spanos et al. (2018), the SD modeling was further developed to incorporate the linearized representations of the end stop force, mooring force, viscous drag force, Coulomb damping and partial overlap effect between the stator and translator of the direct-drive PTO.
To the authors' knowledge, a SD model considering the PTO force saturation has not been reported yet, but it is expected to serve as an powerful tool in the PTO sizing. In addition, most of previous papers regarding the SD modeling mainly focused on model development or validation, and there is limited practical application of the SD modeling reported in the field of wave energy. Hence, the objective of this paper is to develop a SD model of the PTO force saturation and demonstrate its reliability in the application of the PTO sizing.
In this paper, focusing on a spherical heaving point absorber, the representation of the PTO force saturation is derived for the SD modeling by statistical linearization. Next, for comparison, a FD model and nonlinear TD model are established. The proposed SD model is verified by the nonlinear TD modeling along various PTO force limits and wave states. Besides, the estimates of the power capture from three models are compared, in which the relative errors and computational efficiencies are presented. Furthermore, the proposed SD model is applied to a case study of reducing the LCOE by optimizing the PTO size. A realistic wave site and preliminary economic model are considered. The impact of the numerical modeling on the tuning of PTO damping coefficients is analyzed. The reliability of the SD model on the prediction of the AEP, LCOE and optimal PTO size is demonstrated.

Concept description
This section describes the considered WEC concept. A generic heaving point absorber is used as the WEC reference in this study. The schematic of the WEC is shown in Fig. 1. The geometry of the floating buoy is a sphere with a radius of 2.50 m and the buoy is semisubmerged in still water. Given the water density of 1025 kg∕m 3 , the mass of the buoy is calculated as 33543 kg according to Archimedes's principle. Besides, the WEC is assumed to operate in deep water.
The spherical buoy is connected to the moving part of the PTO, and the moving part could be a piston or generator translator depending on the PTO type. Furthermore, a passive control strategy is implemented for the studied point absorber, which implies that only a force proportional to the buoy velocity is applied by the PTO system.

Frequency domain modeling
In this subsection, the FD model of the studied WEC is established based on linear wave theory. As the device in this paper is assumed to oscillate only in heave motion, then the model is only discussed for this degree of freedom. According to Newton's second law, the motion of the WEC as a rigid body can be described as where m is the mass of the oscillating body, is the hydrostatic force, is the wave excitation force, is the wave radiation force, is the PTO force and is the velocity of the rigid body. If the body is assumed to perform harmonic motion under regular waves and a linear PTO model is used to simulate the behavior of the PTO system, (1) could be rewritten in the form of complex amplitudes (Johannes Falnes, 2003), aŝ where ( ) is the hydrodynamic damping coefficient, is the PTO damping coefficient, is the angular wave frequency, ( ) is the added mass of the buoy,̂is complex amplitude of the vertical velocity, and is the hydrostatic stiffness. Then, by solving (2), the complex amplitude of velocitŷcould be obtained aŝ If the ocean wave inputs are assumed to be Gaussian processes, then the stochastic waves can be represented by the linear superposition of a set of frequency components with a random phase. Hence, the wave elevation is expressed as where is time, ( ) and ( ) are the wave amplitude and phase of the regular wave component corresponding to . In a predefined wave spectrum, the amplitude of the wave component is related to the wave energy spectrum , as The variance of the wave elevation 2 is calculated as where is the standard deviation of the wave elevation. Similarly, as the velocity amplitude of WEC corresponding to each wave component can be obtained by (3), the standard deviation and spectral density of J. Tan et al. the WEC response can be calculated. Then, the mean absorbed power can be derived as where and denote the spectral density and standard deviation of the velocity of the WEC.

Time domain modeling
The above FD modeling is built based on the assumption that the system is fully linear. In the TD modeling, nonlinear terms can be taken into account. In the TD modeling, the Cummins equation (Cummins et al., 1962) is applied to describe the dynamic response of the floating body by the expression: where (∞) is the added mass evaluated at the infinite frequency and is the radiation impulse function, and they can be calculated based on the results of ( ) and ( ) (Cummins et al., 1962;Ogilvie, 1964). , and ℎ represent the excitation force, PTO force and hydrostatic force respectively. The representation of wave excitation force is expressed according to the predefined incident wave spectrum and random phase assumption, as where represents the excitation force coefficient which is the excitation force normalized to the wave amplitude of the harmonic wave component, is wave amplitude of each frequency component, and ( ) is the phase angle of the excitation force to the corresponding wave component.
The PTO force saturation can be included, and a typical force saturation effect is depicted in Fig. 2. Once the PTO force limit is set to , the is expressed as The time averaged absorbed power is calculated as The only nonlinear behavior addressed in this TD model is the PTO force saturation, and all the other force components are considered to be linear. It is realized that the nonlinear hydrodynamic and hydrostatic forces would also be relevant to the reliability of the modeling when wave heights are large relative to the diameter of the WEC. However, this paper is intended to demonstrate the relevance of PTO force saturation and therefore other nonlinearities are omitted.

Spectral domain modeling
As an extension of the FD modeling, the SD modeling is formulated based on the assumption that the system is linear and that the inputs and outputs are uncorrelated among different frequency components. Therefore the superposition principle is valid. Besides, all the dynamic responses of the system are assumed to have a Gaussian distribution.
The SD modeling introduces quasi-linear coefficients to include the considered nonlinear effects. The derivation of the quasi-linear coefficients is based on statistical linearization, in which the contributions resulting from the nonlinear effects of the considered frequency components are taken into account. The derivation of statistical linearization is demonstrated in Appendix A.
According to the derivation of statistical linearization, the quasilinear damping to represent the PTO force saturation can be calculated as where ⟨⋅⟩ denotes the expected value of a variable.
As the upper and lower limits of the PTO force in expressed in (10) are symmetric with regard to the horizontal axis, then substituting (10) into (12) gives where ( ) is the probability density function of Gaussian distribution, which is expressed as where is the standard deviation of the buoy velocity. Then, (13) can be solved as where ''erf'' represents the error function, and 1 is the ratio between the PTO force limit and the PTO damping coefficient, expressed as Therefore, (2) in the SD modeling is adjusted in order to include the quasi-linear damping coefficient , , aŝ However, (17) cannot be solved directly since the standard deviation in (15) is unknown, and an iteration process is required. The initial guess of is taken from the value calculated by the FD modeling, and the iteration continues until the difference between the previous and iterative value converges within a certain margin. The mean absorbed power in the SD modeling is then calculated as

Implementation of simulation
The simulation set-up is presented in this subsection. In this work, a JONSWAP spectrum together with peakedness factor of 3.3 is applied to represent the irregular waves. For each wave state, 500 individual harmonic wave components with a random phase between frequency components are considered. The angular frequencies of the wave components are uniformly spaced from 0.05 to 4 rad/s. Besides, the significant wave height considered in this paper ranges from 0.5 m to 5.5 m.
The hydrodynamic characteristics of the buoy, including ( ), ( ) and ( ), are calculated numerically using the Boundary Element Method through the open source software Nemoh (Penalba et al., 2017), and the results are presented in Appendix B. The identical hydrodynamic results are applied to the three different models.
In the TD modeling, directly computing the convolution integral of radiation force in (3) is computationally expensive (Folley, 2016). However, it can be represented using an state-space approximation. In this paper, the frequency domain identification method (Pérez and Fossen, 2008) is used to determine the state-space parameters and (∞) as well. The TD simulation was implemented in the commercial Matlab environment and solved using a numerical integration scheme based on the ODE45 solver. The initial displacement and velocity of the buoy are set to zero. The simulation time duration and time step are set to 200 times and 0.01 times the considered peak period respectively. A ramp function is used to avoid strong transient flow at earlier time steps, and the ramp time is chosen as 25 (Lawson et al., 2014). The duration of the ramp time is not included in the calculation of the average power absorption. To mitigate random errors, the TD model is re-run 10 times for each case, and the mean values are calculated. In the SD simulation, as the iterative process is required to solve the standard deviation of the buoy velocity, a convergence criteria of 0.01% is selected (Folley and Whittaker, 2013a;Silva et al., 2020).

Verification of the SD modeling
The TD modeling is used as the verification reference in this paper. The SD modeling is established based on Gaussian assumption, but the addition of nonlinear force components makes the assumption less valid. Hence, it is necessary to have an insight into the validity of the assumption when the considered nonlinear components take effect. Figs. 3 and 4 show the normalized probability density of the buoy velocity calculated by the nonlinear TD modeling for different wave heights and PTO force limits. Besides, in these figures, the probability density of the Gaussian distribution with the standard deviation calculated by the TD modeling accordingly is also shown as the reference. In general, the responses of the WEC remain Gaussian to a suitable extent. However, the Gaussian assumption tends to be relatively less valid when the wave height turns higher or the PTO force limits become stricter, which implies a reduction of the reliability of the SD modeling. For instance, in Fig. 4, the difference of the probability density for the PTO force limit of 30 kN is larger than those for the PTO force limit of 50 kN and 100 kN. This happens since the highly constrained force limit makes the PTO force frequently saturated, in which the nonlinearity is highly intensified. This results in a decrease of the validity of the Gaussian assumption, even though the quantitative difference is hardly noticeable at the low PTO force limit of 30 kN. Nevertheless, the SD modeling is commonly applied for the performance evaluation or earlystage design of WECs, where the contribution of the highly powerful wave states and strong force saturation is limited. Therefore, the use of Gaussian assumption is thought to be reasonable.
To verify the SD modeling, its calculated power spectral density, standard deviation of the responses are compared with those obtained by the nonlinear TD modeling. For the results of the TD modeling, Fourier analysis is used to obtain the response components in different frequencies. Fig. 5 shows the power spectral density of the buoy velocity obtained by three different models. The shaded area in Fig. 5 represents the standard deviation of TD results by running the model multiple times. It can be seen that the SD modeling has a good agreement with the TD modeling, and its accuracy is clearly better than the FD modeling. Figs. 6 and 7 show the standard deviations of the buoy velocity calculated by the TD and SD modeling, and different wave states and PTO force constraints are considered. In Figs. 6 and 7, the standard deviations predicted by the SD and TD modeling tend to converge as the PTO force limit keeps increasing to the point when the system can be considered linear. Besides, it can be found that changes in the standard deviations predicted by the SD modeling are more notable from those calculated by the TD modeling with the increase of and the decrease of the PTO force limit. In Fig. 6, the maximum difference is 2.4% occurring at of 7.28 s and of 5 m, and it is 3.2% at the PTO force limit of 20 kN and of 3.5 m in Fig. 7. In mild wave states and PTO force constraints, the difference is even more limited. Therefore, the proposed SD modeling is thought to be reliable for the use in the PTO sizing.

Comparison of three models on the prediction of the power absorption
In this section, the accuracy and computational efficiency of three models are compared. Figs. 8 and 9 depict the capture width ratio (CWR) predicted by the three models and the average relative errors of the FD and SD modeling with respect to the TD modeling. Different wave heights and PTO force constraints are taken into account. In Figs. 8(b) and 9(b), the resulting errors are averaged over the three considered wave periods. The CWR is defined as the absorbed power divided by the wave power per unit length of wavefront and characteristic length of the WEC. It can be seen from Figs. 8 and 9 that the FD modeling has an acceptable reliability at low wave heights or large PTO force limits even though the force constraint effect is not included at all. This is because the PTO force constraint is not significant in these cases. When the force limit turns tighter or the wave height becomes higher, the force saturation occurs more often, and then the FD modeling can hardly give a reliable prediction. For instance, when is 5 m and is 7.28 s in Fig. 8(a), the CWR predicted by the TD modeling is approximately only half of that given by the FD modeling. In Figs. 8(b) and 9(b), the average relative error of the FD modeling even reaches 100% and 200% when the is 5 m and the PTO force limit is 20 kN respectively. As for the SD modeling, it presents comparable results with the TD modeling throughout all the considered conditions. At relatively low wave heights and large  force limits, the prediction difference on CWR between the SD and TD modeling is negligible. It can be noted that the accuracy of the SD modeling tends to decrease with the wave height and the strictness of the PTO force limit where the nonlinearity is stronger. However, the relative errors of the SD modeling are dramatically lower than those of the FD modeling, and the maxima is around 20% which only occurs at the extreme PTO force limit of 20 kN. This is an obvious improvement on the accuracy compared with the FD modeling. Furthermore, it can noticed that the CWR shown in these two figures are relatively lower than the reported values for point absorbers in literature (Babarit et al., 2012). This is because the PTO damping coefficient is not optimized here for the considered operating conditions. The relative computational time of the three models is compared in Table 1, and all the simulations are run in the same machine with an Intel i7/2.80 GHz processor. The TD modeling computational time is counted as the time of running one single simulation. It is seen that the computational efficiency of the SD modeling remains in a similar level with the FD modeling, and the TD modeling is dramatically more timedemanding. Therefore, compared with the FD and TD modeling, the SD modeling could clearly improve the accuracy of power prediction and still retain sufficient computational efficiency. It should be clarified that the computational time of the TD modeling depends on the run length and time step length of the simulation. However, the setup of these two parameters is important to the accuracy of the modeling. One comparable example is given in Lawson et al. (2014) where the run length was set to 125 with a time step of 0.01 in the TD modeling for calculating the motion response and power production while the run length in this paper is set to 200 . Nevertheless, the computational time of the TD modeling is several orders of magnitude higher than the SD modeling.

Case study: Determining the optimal PTO size for yeu island
In this section, three models are implemented to do the PTO sizing for a chosen wave site respectively. It is intended to demonstrate the performance of the SD modeling in a practical application. The impacts of the numerical modeling on tuning of the PTO damping coefficients, techno-economic performance and determination of the PTO size are identified.

Wave resource
As a realistic wave site, Yeu island is considered in this case study. It is located in the oceanic territory of France (Babarit et al., 2012), and its hours of occurrence of each wave states is given in Fig. 10. In Fig. 10, represent the mean zero-crossing wave period. The occurrence of wave states with higher than 5 m is not depicted, because the WEC is stopped from operating in order to protect the structure in severe wave states. Given the JONSWAP spectrum, the zero crossing period is transferred to by multiplying with a factor of 1.287 (Myrhaug and Kjeldsen, 1987;Journée et al., 2015).

Economic modeling
The PTO size is optimized with the objective to minimize the LCOE. The economic model proposed in Tan et al. (2021) is applied to calculate the LCOE. The CAPEX of WECs is divided into two portions, namely mass related cost and power related cost , expressed as where and can be calculated as where , & and are the cost related to the structure, foundation & mooring and the installation respectively, and are the cost of the PTO and the grid connection; and , & , , and are their corresponding percentages in the total CAPEX, and they are determined as 38.2%, 19.1%, 10.2%, 24.2% and 8.3% respectively (Tan et al., 2021).
Following Tan et al. (2021), a steel price of 1.95 Euros/kg is considered and the buoy structure cost is calculated by assuming that it totally comes from the steel cost. The density of the buoy is assumed to be half of the water. The PTO system is assumed to be a direct drive generator and all PTO costs come from the generator. The amount of required active material depends on the PTO force limit and the force density of the generator. As in Polinder (2013), the maximum force density per surface area in this work is assumed as 44 kN∕m 2 , which generally ranges from 30 to 60 kN∕m 2 depending on the design. The cost of active material of this generator in series production is estimated as 14,600 Euros∕m 2 (Tan et al., 2021). The PTO cost is considered as twice J. Tan et al.  the active material cost to account for manufacturing. The detailed parameters of the economic model are presented in Appendix C.
In this paper, the annual OPEX is assumed as 8% of the CAPEX, and the discount rate is assumed as 8% with a lifespan of 20 years, based on De Andres et al. (2016). Then, the LCOE of WECs is calculated as where represents the total years of the lifespan, and represents the evaluated year. The AEP (Annual Energy Production) of the WEC producing at the sea site is calculated as where is the overall conversion efficiency from the annual absorbed energy to the AEP and is assumed 0.7 (Chozas et al., 2014); is the availability of WECs to work and it is set as 0.9 due to the necessary operation and maintenance (Kramer et al., 2011); represents the total hours of the occurrence of a certain sea state, which is presented in Fig. 10; represents the sea state and is the number of concerned wave states.

Incorporation of viscous drag force
In the AEP calculation, the viscous drag force is taken into account in the SD and TD modeling for better accuracy of the technoeconomic assessment. In the TD simulation, the viscous drag force can be represented by a quadratic damping term J. Tan et al. Fig. 11. Comparison of capture width ratio predicted by different models, = 30 kN, = 70 kNs∕m and = 5.15 s. In the legend, 'saturation' represents the models including PTO force saturation as the only nonlinearity, and 'saturation + vis' represents the models with both nonlinear effects of the PTO namely force saturation and viscous drag force.
where is the water density, is the drag coefficient and is the characteristic area of the buoy perpendicular to the moving direction. The drag coefficient is selected as 0.6 to minimize the error of the power estimate based on the investigation reported in Giorgi and Ringwood (2017), in which the research reference is also a sphere.
For the SD modeling, the linearized viscous damping has been derived in Folley andWhittaker (2010, 2013a) and Silva et al. (2020), and it is expressed as Then (17) After including the viscous drag force into the SD and TD modeling, other set-ups and the solution process remain identical as demonstrated in Section 3. The CWR calculated by the models with and without considering the viscous drag forces are compared in Fig. 11. It is clearly reflected in the SD and TD modeling that the incorporation of the viscous drag force decreases the power absorption. In addition, the SD modeling still presents fairly close results to the TD modeling even though two nonlinear effects, namely the PTO force saturation and nonlinear viscous drag force, are taken into account at the same time. It also indicates that more nonlinear effects can be effectively considered in the SD modeling to further improve the accuracy. However, this is out of the scope of this work which is intended to focus on the effect of the PTO force saturation. The derivation of other commonly concerned nonlinear effects in the SD modeling can be found in Folley and Whittaker (2010), Gunawardane et al. (2017) and Silva et al. (2020).

The tuning of PTO parameters
The PTO damping coefficients are significant to the power absorption of WECs. To provide an insight into the dependence of the optimization of the PTO damping coefficients on the numerical modeling, a comparison is made among the models. Besides maximizing the power capture, the probability of occurrence of the PTO force saturation also needs to be considered in the determination of the PTO damping coefficients. Taking the direct-drive generator as an example, the PTO force saturation is associated with large currents, and the highly frequent PTO force saturation might lead to overheating conditions. Two methods are introduced below to incorporate the probability of occurrence when tuning the PTO parameters of WECs.

• Probabilistic method
As the WEC is subjected to random wave inputs, the probability of occurrence can be estimated by the probabilistic analysis. If the dynamic process is assumed to be Gaussian with zero mean, the amplitude of the variables can be characterized by Rayleigh distribution (Journée et al., 2015). Hence, the probability of PTO force saturation can be calculated by where represents the standard deviation of PTO force. Then, a tolerance on the probability should be defined, and it is tackled as a constraint in the optimization of PTO parameters. This method has been used in Backer (2009) for the hydrodynamic optimization of WECs. However, the tolerance on the probability is dependent on plenty of factors, such as the PTO type and design, scale of the WEC, severity of the load, material strength, maintenance frequency and so on. It is expected that the tolerance on the probability of the PTO force saturation would differ from one project to another in practice. Therefore, it is challenging to properly determine the tolerance due to the lack of information from practical WEC projects. The relation of the absorbed power and probability of the force saturation to the PTO damping coefficient is calculated by the three models respectively, as depicted in Fig. 12. It is noticed that the selection of the numerical model makes a notable difference to the PTO damping optimization. For instance, the PTO damping coefficient with the highest power absorption is 75 kNs∕m in the FD modeling. Comparatively, in the SD and TD modeling, the power absorption tends to be saturated with the increase of the PTO damping without reaching a certain peak point. This is because the PTO force is highly limited by the saturation effect when the PTO damping coefficient is large. The further increase of the PTO damping has a limited effect on the standard deviation of the PTO force or the system dynamics. Besides, it is observed that the maximum power predicted by the FD modeling is higher than other models, and this can be explained by the existence of the viscous drag force in the SD and TD modeling. In addition, compared with the SD and TD modeling, the FD modeling overestimates the probability of the PTO force saturation. The reason is that the PTO force saturation does not take effect in the FD modeling, and hence the resulting PTO forces are accordingly higher. Moreover, the difference between the probabilities calculated by the SD and TD modeling is negligible, which implies the good reliability of the SD modeling on predicting the standard deviation of the PTO force. Fig. 13 shows the PTO damping and absorbed power optimized by the three models, and different allowable probabilities of the PTO force saturation are considered. It can be seen that both the optimal PTO damping and absorbed power increase with the allowable probability. In addition, the PTO damping optimized by the FD modeling is lower than those corresponding to the SD and TD modeling, and the difference becomes larger with the increase of the probability. For instance, the optimal PTO damping of the SD and TD modeling is around twice the value of the FD modeling when the probability is 40%. The reason is that the FD modeling overestimates the standard deviation of the PTO force as well as the probability of the exceedance, thus its optimized PTO damping should be correspondingly lower under certain tolerance. Comparatively, the SD and TD modeling result in similar optimal PTO damping coefficients throughout all the considered allowable probabilities. It is noticed in Fig. 13 that the FD modeling only presents a slightly higher power estimate than the SD and TD modeling, although its optimal PTO damping is very different. This observation can be explained by analyzing how the optimal PTO damping is found during the optimization. In most of cases, the optimal absorbed power is obtained when the allowable probability of the PTO force is reached. Then, the standard deviations of the PTO force at the optimal conditions are identical for the SD and FD modeling, and the standard deviations of the PTO force are expressed as As the PTO force is significantly more dominating than the viscous drag force in the SD modeling, the equivalent PTO damping , in (29) is close to in (28) after the PTO damping optimization. In this way, the dynamics estimated by the two models tend to be similar. Then, the absorbed power of the SD and FD modeling are also comparable, and the slight difference results from the viscous drag force. This explanation can be verified by Fig. 14, which depicts the , for the SD modeling and for the FD modeling after the PTO damping optimization. It is seen that the equivalent viscous damping is highly limited in this case. Besides, the equivalent PTO damping , is very similar with the PTO damping optimized by the FD modeling, even though the value of the optimized PTO damping by the SD modeling is much larger. This explains why the values of the optimal absorbed power of the three models are close after the PTO damping optimization. However, it should be noted that this observation is limited to this particular case. When other nonlinear effects are considered or the viscous damping is more influential, the difference in dynamics estimated by the SD and FD modeling could become higher with the optimized PTO damping coefficients. As a result, the values of the optimal absorbed power from the FD and SD modeling are also expected to have higher variations. However, further investigation is out of the scope of the present paper.
• Transferred method A simplified method can also be used to determine the PTO coefficients without defining the tolerance and calculating the exceedance, and this method has been used in Tan et al. (2021) and Tedeschi et al. (2011). Specifically, in this method, the irregular wave states are transferred to corresponding regular wave states by equaling their time-averaged power transport per unit length of wave front. Then, the PTO parameters are selected to suit the transferred regular wave states, namely energy period and ∕ √ 2. As the PTO force amplitude and time average power can be explicitly derived for the regular wave states, the optimal PTO damping coefficient can be obtained to maximize the absorbed power and guarantee the PTO force amplitude to be lower than the PTO force constraint. For convenience, it is called the transferred method in the following text. It is realized that this tuning method is only an approximation. Therefore, in irregular wave states, the maximization of the power absorption cannot be guaranteed and the violation probability of the force constraint is not necessarily within a certain tolerance either. In Fig. 15, the optimal values of PTO damping obtained by the probabilistic method are normalized to those obtained by the transferred method. For comparison, the SD and FD modeling are applied in the probabilistic method respectively. Different wave periods and tolerances on the probability of the force saturation are considered. It is seen from Fig. 15 that the tolerance has a remarkable impact on the PTO damping determination, and easing the tolerance increases the selected optimal damping coefficients. In this case, the transferred method implies a tolerance on the probability of approximately 20%, if the results of the SD modeling are taken as the reference. In addition, the PTO damping optimization is highly dependent on the numerical modeling,   and the difference between the optimal damping coefficients selected by the FD and SD modeling increases proportionally to the tolerance. In particular cases, the difference reaches up to 100%. Fig. 15. The optimal PTO damping coefficients with different probabilities of force saturation obtained by the FD and SD modeling normalized to those obtained by the transferred method, = 3.5 m, = 50 kN.

The techno-economic analysis
To search for the optimal PTO size, the techno-economic performance of the WEC is calculated for different PTO force limits. The step length of the PTO force limit is selected to be 10 kN. In the AEP calculation of the techno-economic analysis, both of the above methods are used to tune the PTO parameters for comparison. The PTO parameters are updated for each sea state.

LCOE using the probabilistic method
In this method, the PTO parameters are optimized by each model independently, and hence the effect of the numerical modeling on tuning the PTO is reflected in the AEP and LCOE calculation. The allowable probability of the PTO force saturation acts as a constraint during optimization, and it is defined as 20% here for a preliminary assessment. Regarding the PTO damping optimization by the SD and FD modeling, the ''interior point'' algorithm in Matlab environment is used, and the tolerance of the function is set as 1e-5. The bound of the PTO damping is set as [0,20 ( )] for each sea state, in which is the peak frequency of the irregular wave state. As for the optimization using the TD modeling, an exhaustive searching process is used for a range around the optimal PTO damping optimized by the SD modeling , , namely [0.9 , , 1.1 , ], and the searching step is selected to be 0.01 , . The selected searching range is thought to be fair, since the deviation between the PTO damping optimized by the SD modeling and by the TD modeling is observed to be limited, as shown in Fig. 13. This is intended to save the optimization time given the low computational efficiency of the TD modeling. Table 2 and Fig. 16 present the AEP and LCOE of the WEC with different PTO force limits respectively. It can be seen that the relative errors of the SD modeling to the TD modeling are less than those of the FD modeling. Regarding the AEP prediction, the maximum error of the SD modeling is 3.5% while it is 5.5% for the FD modeling. However, the three models generally indicate similar results, and they also lead to the same optimal PTO sizing in this case. This can be expected, because the difference of the absorbed power of the three models is limited when using the probabilistic method to tune the PTO parameters. The detailed explanation can be found in Section 6.4. However, this conclusion could not be simply extended to other scenarios, particularly when other nonlinear effects are more influential than the PTO force saturation.
To provide an insight into the difference among the three models on the determination of the PTO parameters, the weighted mean values of the optimal PTO damping coefficients are calculated. In the calculation, Fig. 16. LCOE of the WEC predicted by three different models, and PTO parameters are tuned based on the probabilistic method.

Table 2
The comparison of AEP (MWh) predicted by three models and the relative errors of the FD and SD modeling to the TD modeling, and PTO parameters are tuned based on the probabilistic method. the value of hours of the occurrence of each sea state is considered as the weight. Fig. 17 shows the weighted mean of the optimized PTO damping coefficients estimated by the FD and SD modeling normalized to that estimated by the TD modeling. It is noticed that the SD modeling has a good reliability on tuning PTO parameters, and the discrepancies are limited to 2%. In contrast, the optimal PTO damping estimated by the FD modeling is much lower, with an error of more than 12% relative to the TD modeling when the PTO force limit ranges from 20 to 30 kN. This implies that the FD modeling is insufficient to identify the optimal PTO parameters.

LCOE using the transferred method
In the transferred method, the PTO damping coefficients applied in the three models are kept identical, and they are updated to each sea state. In the sense, the relative errors on the AEP and LCOE estimates of the FD and SD modeling only result from the modeling itself.
The AEP of the WEC predicted by different models is shown in Table 3. It can be seen that both the FD modeling and SD modeling result in an overestimation of the AEP compared with the TD modeling. However, the predication deviation between the FD modeling and TD modeling is significant, and the relative error increases as the PTO force limit turns stricter. When the PTO force limit is 20 kN, the deviation is over 24%. In this sense, the reliability of the FD modeling in tackling cases where the PTO force saturation is relevant is clearly insufficient. In contrast, the SD modeling presents a much better accuracy on the prediction of the AEP, the relative errors are observed to be no more than 4.3%, ranging from 2.0% to 4.3%, along all the considered cases.    Fig. 18 shows the LCOE predicted by these three models. The FD modeling shows an over-optimistic result, and its corresponding lowest LCOE is around 0.362 Euros/kWh which is 11.3% lower than that predicted by the TD modeling, namely 0.408 Euros/kWh. The SD modeling presents a close result with that obtained by the TD modeling, and its resulted lowest LCOE is 0.394 Euros/kWh, with a deviation of 3.4% from the TD modeling. In addition, it is seen that the FD modeling indicates a slightly lower optimal PTO size than the SD modeling and TD modeling do. The FD modeling corresponds to the optimal PTO size of 80 kN while the SD modeling and TD modeling are associated with 90 kN. In this case, considering the step length of the PTO force limit, the difference between the FD modeling and TD modeling on the PTO size determination is not significant. However, the PTO sizing differs with the wave resource and types of WECs, and the insufficient accuracy of the FD modeling might lead to an unacceptable misestimate of the optimal PTO size in other scenarios.
Regardless of the PTO tuning method, the SD modeling always presents a good agreement with the TD modeling on estimating the AEP and LCOE. At the same time, the application of the SD modeling contributes to a huge reduction in computational time in the technoeconomic assessment. Specifically, with the explicit PTO parameters, the FD and SD modeling consume 1.02 s and 2.11 s to complete the AEP calculation for the considered thirteen PTO force limits in this case. However, the TD modeling takes more than 1.6 h if the simulation is run single time for each sea state, and this is approximately 3000 times the duration for the SD modeling. Regarding the optimization time in the probabilistic method, it takes over 30 h with counting the singlerun time for the TD modeling to find all the necessary optimal PTO parameters, but it is only 75 s for the SD modeling.

Discussion
This paper suggests that the PTO force saturation can be incorporated into the SD modeling of WECs. As the SD modeling combines the high efficiency and adequate reliability, it could make a contribution to the WEC optimization related to the PTO sizing. However, in this work, the PTO force saturation is only derived for the ideal PTO system behaving as a passive component or damper. In practice, the modeling of practical PTO systems could be more complicated depending on their designs and types, which might lead to different linearization solutions. For instance, in hydraulic PTO systems, the coupling of the hydraulic circuit with the rest of the components makes the PTO force profile highly nonlinear even without the force saturation (Plummer and Schlotter, 2009;António, 2007). Therefore, further development is needed to incorporate the effect of the force limit in the SD modeling for more realistic PTO systems, and then the accuracy and computational efficiency also remain to be investigated.
In this work, only a single WEC is considered in the case study for a preliminary assessment. Given the high computational efficiency when compared to other approaches, it is attractive to apply the SD modeling to the optimization or design of WEC farms, which are characterized by large computational loads. The reliability of the SD modeling in predicting the performance of a WEC farm has been validated in Folley and Whittaker (2013b). Thus, the next phase of this work is to extend the representation of the PTO force saturation to the SD modeling for WEC farms. Besides, the present SD model has only been developed for a point absorber with a single degree of freedom, but it is feasible to extend the modeling to multiple degrees of freedom, as demonstrated in Silva et al. (2020).
The current SD modeling is formulated based on linear wave theory with a Gaussian process, which limits the method to be used extensively. For instance, the relevance of the wave nonlinearity becomes more important in large wave heights, and then Gaussian assumption would be challenged. Although the assumption of the Gaussian distribution can be thought adequately valid in most operational regions of WECs, the source of the error resulting from the assumption is still unavoidable. To ease the limitation, as stated in Folley (2016), the consideration of non-Gaussian responses will be a promising future work, which could contribute to a more accurate estimation.

Conclusions
In this paper, a SD model is developed to include the effect of the PTO force saturation in the heaving point absorber WEC. The SD model is constructed based on the framework of the FD modeling, and the nonlinear effect is incorporated by a quasi-linear term. The quasilinear term is derived by statistical linearization, which is based on the assumption that all the responses of the system follow Gaussian distribution. The proposed SD model is verified against the nonlinear TD modeling, and a reliable accuracy is identified over a range of relevant wave conditions and force constraints.
Secondly, a comparison among the FD, SD and TD modeling is conducted regarding the prediction of power performance. In low significant wave heights or large PTO force limits, three models present similar power performance. With the increase of the significant wave heights and strictness of the PTO force limit, the FD modeling tends to be insufficiently reliable and the highest error can reach values up to 200% relative to the nonlinear TD modeling. However, the SD modeling shows a good reliability, and its error relative to the nonlinear TD modeling is only around 20% even at the strict PTO force constraint. At the same time, the computational efficiency of the SD modeling is comparable with the FD modeling and significantly higher than the TD modeling.
Thirdly, a case study of the PTO sizing is performed to demonstrate the potential of the SD modeling in practical applications. A realistic sea site and preliminary economic modeling are considered. The SD modeling proves to be a desired alternative for the PTO sizing and other similar applications of performance evaluation or early-stage design optimization. It combines the high computational efficiency and the high accuracy on the performance prediction. Its errors on the prediction of the AEP are no more than 4.3% relative to the TD modeling, and it shows an agreement on the optimal PTO size with respect to the TD modeling. In contrast, the FD modeling could lead to a misestimate on the PTO sizing, and its prediction of the AEP and LCOE is questionable, especially with strict PTO force limits. The relative error of the FD modeling on the predication of the AEP is over 24% in the particular case. In addition, with considering nonlinear effects, the SD modeling suggests an adequate reliability on tuning the PTO parameters, while the FD modeling is found to be insufficient.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix A. The derivation of statistic linearization
The derivation of the quasi-linear coefficients in the SD modeling is briefly introduced here, and more details can be found in Folley andWhittaker (2010, 2013a), Folley (2016). Assuming a general differentiable nonlinear function ( ) and its linear equivalent function ( ) are given as If the nonlinear function ( ) satisfies the zero-mean Gaussian distribution, its expected value is as zero and then only to be solved.

Appendix B. Hydrodynamic results
The hydrodynamic coefficients of the concept calculated by Nemoh are presented in Figs. B.19 and B.20.

Appendix C. Parameters in the economic modeling
The proportion of each cost component to the total CAPEX is depicted in Table C.4.
The parameters used in the techno-economic analysis are described in Table C.5.