Time-Splitting Coupling of WaveDyn with OpenFOAM by Fidelity Limit Identified from a WEC in Extreme Waves

: Survivability assessment is the complexity compromising Wave energy development. The present study develops a hybrid model aiming to reduce computational power while maintaining accuracy for survivability assessment of a Point-Absorber (PA) Wave Energy Converter (WEC) in extreme Wave Structure Interaction (WSI). This method couples the fast inviscid linear potential ﬂow time-domain model WaveDyn [1] with the fully nonlinear viscous Navier–Stokes Computational Fluid Dynamics (CFD) code OpenFOAM [2]. The coupling technique enables the simulation to change between codes, depending on an indicator relating to wave steepness identiﬁed as a function of the conﬁdence in the linear model solution. During the CFD part of the simulation, the OpenFOAM solution is returned to WaveDyn via an additional load term, thus including viscous effects. Developments ensure a satisfactory initialisation of CFD simulation to be achieved from a ‘hot-start’ time, where the wave-ﬁeld is developed and the device is in motion. The coupled model successfully overcomes identiﬁed inaccuracies in the WaveDyn code due to the inviscid assumption and the high computational cost of the OpenFOAM code. Experimental data of a PA response under extreme deterministic events (NewWave) are used to assess WaveDyn’s validity limit as a function of wave steepness, in order to validate CFD code and develop the coupling. The hybrid code demonstrates the applicability of WaveDyn validity limit and shows promising results for long irregular sea-state applications.


Introduction
Wave energy is a promising sector for sustainable low-carbon energy production.However, despite the efforts made in the last decade, uncertainties in wave energy performance and cost remain.In 2016, survivability assessment has been identified as the main Wave Energy Converter (WEC) concern in the five years to come by the Collaborative Computational Project in Wave Structure Interaction (CCP-WSI) Working Group [3].Numerical modelling can play a key role in the assessment of extreme responses that are necessary to the understanding of WEC capacity to survive extreme conditions.
WECs are subjected to complex hydrodynamics phenomena, such as slamming waves or green water effect.Numerical models that are based on the weakly simplified Navier-Stokes equations, assumed of 'high-fidelity' (Figure 1), are necessary to model such hydrodynamics and the interaction with the device, whereas numerical models of 'lower fidelity' (Figure 1) are restricted by the assumptions made on the fluid properties (e.g., inviscid, or small displacement) or the equations solving procedure (e.g., linear potential flow).Therefore, their use is debatable for survivability assessment of WEC, which, in the meantime, justifies the use of Computational Fluid Dynamics (CFD) models ( [4][5][6], among others).
However, assumed high-fidelity codes like CFD are too time-consuming for routine design use in industry [3].Codes of lower fidelity are used instead, although being based on major physical assumptions inducing concerning uncertainties, especially in surge and pitch [7].Inaccurate estimations of extreme response often lead to conservative design decisions that lower the commercial viability of a device [8].In addition, lower fidelity models were developed by the Oil and Gas and shipping industries.WECs differ to ships and platforms by their size (i.e., WECs are typically smaller) and their motion, which should be accentuated to generate power.Viscous effects tend to become more important for WECs, whereas diffraction and radiation effects tend to decrease.Hence, when radiation diffraction effects dominate viscous effects, Oil&Gas low-fidelity numerical models that are based on simplified physics remain reliable; whereas, when viscous effects dominate, high-fidelity numerical models, like CFD, are justified.During a given simulation-especially for a long irregular sea-state-the regime is likely to change from one to the other.
Uncertainties remain regarding which numerical model to use for a given application.At best, numerical models are classified, depending on the physical processes represented (i.e., wave-breaking, wave non-linearity) [9].The absence of a parametric study function of a wave-parameter makes the choice of the selected model for survivability vague.For long extreme irregular sea-states-as recommended by WEC standards (IEC TS 62600-2 [10])-CFD only simulations remain unaffordable, and so the use of CFD is restricted to single short extreme events [3].However, in large waves up to the breaking point, Coe et al. [11] argue that lower fidelity models are the most efficient, since they give good agreement for a fraction of the computational cost.On the other hand, the accuracy of Degree of Freedoms (DoFs) other than heave is lower due to viscous effects [7] that have a greater impact on pitch and surge than on heave, hence requiring CFD for their modelisation.Due to this lack of consensus on code selection between low and high fidelity models and where each should be applied, a wave parametric definition for the range of use of each code category [12] is needed-this is an intrinsic objective of the present study.
Besides, the requirement for long extreme irregular sea-state simulations demonstrates that no numerical model is capable of accurately modelling extreme Wave Structure Interaction (WSI) of WEC at an affordable computational cost.Hybrid models (or coupled models) have emerged to deal with the issue.Inspired from [13] and [14].
Coupling between models may be classified as weak compared to strong, and loose compared to tight.When one model runs before the other in a unidirectional exchange of information, the coupling is weak [15].This strategy initialises one model by outputs from the other.When the exchange of information is two-way, the coupling is strong, and each model is run simultaneously, thus exchanging information and benefiting from each model's feature during the entire simulation.Tight-coupling is a strongly coupled model, where the two models converge to the same time-step before advancing to next [16].In loose coupling, the time-step of a model is not imposed on the other.
Three mains coupling strategies exist: toolbox, zonal-splitting, and time-splitting.The toolbox technique (or function-splitting) uses an external model to simulate a specific component, for example, the mooring line.Palm et al. [17] use the mooring software MooDy (in-house) to describe catenary lines, where the modelling of the Rigid Body Motion (RBM) and flow are performed using CFD in OpenFOAM.A similar development is done by de Lataillade et al. [18] between the CFD solver PROTEUS [19] and the multi-body dynamics ProjectCHRONO [20].However, although it improves solution accuracy, this strategy increases computational power due to communication between solvers.
A second strategy is a zonal or space-splitting approach.It has been given much attention recently, since it specifically aims to reduce computational cost.Li et al. [21], Hildebrandt et al. [22], Lachaume et al. [23], Paulsen et al. [24], and Biausser et al. [25] weakly coupled a two-phase Navier-Stokes (NS)-Volume of Fluid (VoF) model with a Fully Nonlinear Potential Theory (FNPT) based model solved using the Finite Element Method (FEM) (QALE-FEM) and the Boundary Element Method (BEM), respectively.In these studies, the NS model uses the inputs from the FNPT based model.Hildebrandt et al. [22] justify the use of the weak coupling since no wave reflection comes from the structure located in the NS sub-domain before the extreme event of interest reaches it.Hence, the 3-dimensions (3D) simulation made using the NS solver is limited to immediately before the impact and the impact itself.Strong coupling is often required for long simulations, since sub-domains will affect one other [26].A relaxation-zone (or sponge layer) or an interface provides the connection between sub-domains.Lachaume et al. [23] further develop the weakly coupled into a strongly coupled model by creating a region where both FNPT and NS-VoF sub-domains overlap.Kim et al. [27] use the relaxation-zone technique between a BEM-FNPT sub-domain and a NS-VoF originally developed by Yan et al. [28] for a similar two-way coupling.A drawback of the strongly coupled zonal strategy is that the time-consuming model (i.e., CFD) is used during the whole simulation regardless of the complexity of the event, even in cases where such fidelity is not required.Therefore, there is a period in time and a region on space, where WSI requires such a level of accuracy and, outside this period, lower fidelity models are accurate enough and should be used.
Therefore, the present study develops a coupling based on the time-splitting strategy: a hybrid method that swaps between models for specific periods, depending on an indicator function of the model validity.It is similar to the strategy developed by Wang et al. [29], where three wave models (QSBI, ESBI, ENLSE-5F) for large temporal and spatial scale simulations are ranked by order of validity assumed from underlying equations in order to use the most efficient model (i.e., accurate while minimising computational cost) as a function of the instantaneous wave information.The present hybrid model uses the computationally efficient WaveDyn [1] model (time-domain linear potential flow, i.e., Cummins Equation) to rule the simulation until the WSI requires CFD to maintain solution fidelity.The most appropriate model is always used, which results in a reduction of the overall computational power that is required for the simulation.
The aim of this work is, therefore, to develop a computational tool that is both capable of simulating accurate extreme WSI and is at an affordable computational cost.The strategy employed is the coupling of WaveDyn with OpenFOAM [2].The original contribution of the present coupling is the assessment and definition of WaveDyn validity limit as a function of a wave-parameter that rules the coupling (i.e., choose the solver).
Present work is presented first with the overall method, including the coupling strategy, the experimental set-up, and the numerical models.The result section quantifies the fidelity of the linear potential flow model WaveDyn when compared with OpenFOAM, demonstrating the existence of a validity limit and motivating the present coupling strategy.Once the validity limit is identified and parametrically characterised, the coupling strategy is developed.Fully detailed development is available in [30].

Coupling Strategy
The present study develops a weak, loosely coupled time-splitting hybrid model that aims to benefit from the high-fidelity of the CFD code OpenFOAM and the efficiency of the WaveDyn model.The coupling model overcomes the limitation that is induced by the inviscid assumption made in WaveDyn while improving the computational cost in comparison with a CFD only simulation.The strategy consists of switching between models, depending on the estimated fidelity of the WaveDyn solution defined by a wave parameter indicator function (Section 3.2), as illustrated in Figure 2. At the initial time, t = 0, WaveDyn is initiated to simulate the WSI, while OpenFOAM is not used until : An event occurs that lies outside of WaveDyn's estimated fidelity range.This is known as the hot-start time t hot . 2 The active solver changes to OpenFOAM since the limitation induced by the inviscid assumption is attained.Meanwhile, the WaveDyn solver waits for OpenFOAM inputs.(4).

3
Using a build-up period, t minus , the wave-field and the RBM-state (i.e., position, velocity and acceleration) at t hot − t minus are set up in the CFD OpenFOAM simulation. 4 OpenFOAM CFD simulation starts at t hot − t minus , and OpenFOAM predicted loads and motion-state are transferred to WaveDyn via the External Load Controller (ELC).5 Once the higher fidelity CFD simulation is no longer required and the simulation returns within the fidelity range of WaveDyn, OpenFOAM stops, and the simulation continues in WaveDyn The present study is developed around the first four steps presented above.
Step 5 is left for future work since the complexity of swapping back into WaveDyn is expected to require a dedicated study.The coupling is weak and loose, since OpenFOAM solution is started with the WaveDyn one, and that there is always only one solver ruling simulation (i.e., either WaveDyn or OpenFOAM, neither both) without two-ways exchange of information (i.e., loosely coupled).A potential advantage of a strongly coupled model would be to have WaveDyn calculating specific components which are difficult to implement in OpenFOAM in a function-splitting coupling strategy.As this development is the early stage of this coupling strategy, the two-ways strong coupling is left for future work.We choose a loose coupling to avoid problems that are induced by forcing the time-step of one model into the other, and especially to respect the Courant-Friedrichs-Lewy (CFL) condition in CFD simulations.
Throughout the CFD part of the simulation, from t hot − t minus to step 5, the OpenFOAM RBM solution is transferred to WaveDyn.The conservation of momentum-i.e., Newton's second law of motion, ∑ − → F = m.− → a and ∑ − → M = I. dΩ dt -that is applied in both models needs to agree with the other, which is determined by the following pair of equations: where, − → F are forces, − → M are moments; and, the External Load Controller (ELC) inputs a user-specified additional load (and moment) into WaveDyn model.
OpenFOAM forces and moments are defined in the global coordinate system.Moments are transferred to the centre of gravity within the coupling script before being applied in WaveDyn at a given time-step determined by linear interpolation of the two closest time-steps in the OpenFOAM solution.
This time-splitting coupling strategy requires two main developments.First, the assessment of the WaveDyn validity threshold will define the indicator to switch between models when the validity threshold for WaveDyn is crossed.A parametric assessment that is based on wave characteristics is made to indicate code validity as a function of the simulation time.The second development is the hot-start procedure for initiating the CFD simulation, which has been already described in Musiedlak et al. [31].

Experimental Set-up
Experiments [32] were carried out at the 35 m long, 15.5 m wide, and 2.8 m deep Ocean wave-basin at the University of Plymouth COAST Laboratory [33].Waves are generated from 24 paddles and absorbed by the parabolic beach.Wave Gauges (WGs) measure the surface-elevation along and across the wave-tank.
The extreme loading of marine energy devices due to waves, current, flotsam and mammal impacts (X-MED) buoy, first used in the X-MED project [34], as shown in Figure 3 and with characteristics summarised in Table 1, is used to model a generic Point-Absorber (PA)-WEC in survival conditions, hence with no Power Take Off (PTO).The X-MED buoy is a 0.5 m high cylinder with a hemispherical bottom moored by a single mooring line of 66.3 N/m stiffness as in Hann et al. [32].The resulting mooring pre-load when the sinking the buoy is −20.Four extreme events of increasing wave steepness are modelled while using NewWave theory [35]; a design wave widely used for offshore and WEC applications (e.g., [34,36,37]) that defines the average shape of the highest wave with a specified exceedance probability [38].When compared with irregular sea-states, design waves do not include motion history, although this has a known influence on the dynamics of the structure ( [34,39]).However, using irregular sea-states, Rafiee et al. [40] identified the most extreme event for the PTO of their PA device as a NewWave event with a focus phase-shift.This demonstrates the relevance of these four extreme events for survivability assessment.
A 100-year storm hindcast at Wave Hub (T p = 14.1 s, H s = 14.4 m, [41]) is input to the Pierson-Moskowitz spectrum to generate the first design wave (named ST1) from NewWave theory.The amplitude A is defined as: where N is the number of waves in the given sea-state and m 0 the zeroth moment of the spectrum.Assuming 1000 waves in 3 h sea-state [42], amplitude is equal to 0.273 m at 1/50 scale.This wave-group is used to generate three steeper cases (ST2, ST3, and ST4; Figure 4 and Table 2) by multiplying the wave group spectrum peak frequency ( f p) by specific factors to obtain the desired theoretical steepness (kA, where wave-number k is obtained assuming linear theory) values while maintaining crest amplitude (A) ( [32,34]).However, due to non-linear wave-wave interactions [34], the measured amplitude increases and a trial and error process (described in [43]) is required to achieve focusing at the front edge of the buoy at the initial position with a ±10 cm precision.All experiemental surface-elevations that are shown in this study are the measure of the WG at the front edge of the buoy at its initial position during wave-only cases (i.e., no buoy in the tank).Twelve additional events that are based on ST1 are used [44].These wave cases have the same steepness (kA = 0.152) and differ only by the theoretical phase of the underlying components at the theoretical focus location and time; i.e., a focus-phase of 0 degrees represents a theoretically 'crest-focused' wave and a focus-phase of 180 degrees represents a 'trough-focused' wave.This set of events will be called ST1 − f ocused, and a specific event will be referred as ST1 − X where X is the phase (e.g., ST1 − 180 for 'trough-focused').Figure 5 shows four of the twelve events.Experiments on X-MED motion responses under an extreme wave-group demonstrated having acceptable repeatability while using a steeper (kA = 0.35) breaking case [34].Minor differences build-up after the main event for surge and pitch but not for heave.Thus, small changes in the WSI can be seen to induce significant differences in motion responses after a few periods.This repeatability test is considered to be valid for the present study, since it involves a more complex WSI (i.e., breaking wave).

OpenFOAM: laminar RANS
A Numerical Wave Tank (NWT) mirror to the physical one is developed under the CFD open-source code OpenFOAM-4.1.waveDyMFoam solver solves the Reynolds Averaged Navier-Stokes (RANS) equations for an isotropic, Newtonian, incompressible, two-phase flow using the Finite Volume Method (FVM) and the VoF method for the interface capturing scheme.VoF method, in OpenFOAM, uses the MULES method (developed by H. G. Weller [45]) with 3 iterations.Momentum (on the i-axis) and Continuity equation are given as: where, − → u is the velocity vector with u i its components on the i-axis; p is the pressure; ρ the fluid density; ν is the kinematic viscosity (ν = µ/ρ where µ is the dynamic viscosity); and, S i is a source term.At each time-step, one pressure-velocity coupling is solved, within which the pressure and velocity fields are calculated and corrected twice.Solvers schemes are detailed in Table 3.
The flow is considered laminar (i.e., no turbulence model is used) across the study, because the wave-groups considered are below the breaking limit; and because Reynolds number remains below the critical value (Re = 6.65 × 10 5 < 1.5 × 10 6 [46]) estimated from a cylinder of the same diameter in a flow of velocity calculated from linear theory while using the steepest case.An investigation on turbulence models and its effect on drag is left for future work, since the coupling development is the main interest.
As presented in Figure 6 and Table 4, the NWT is made of a 1 m inlet ( [47,48]) acting as a wave-maker (left hand-side); a numerical beach absorbing the waves (right hand-side), and a working section where the rigid-body is (middle).waves2foam [49] performs wave-generation and absorption using relaxation-zones where are imposed the analytical solutions of the wave-group (and the still water for the numerical-beach) gradually along the relaxation length.This solution is chosen over other techniques, because wave-reflections can be avoided in short wave-groups using a sufficiently long numerical beach.The linear superposition that is imposed at the inlet (using combinedWaves) is the filtered surface elevation data from the most upstream WG.The two-dimensional Numerical Wave Tank (NWT) was first developed to ensure precise prediction of the wave-only cases.Square cells (cube in 3D) are used to avoid the influence of cell aspect when grading is used.The time-step is governed by the Courant-Friedrichs-Lewy (CFL) condition to maintain numerical stability [50].The CFL number (C0 = |u|∆ t ∆ x ) is restricted to 0.5; a value commonly used in WSI [14], which shows no difference in the surface elevation at the focus location of a NewWave event compared to 0.25 [51].Mesh resolution (i.e., 20 cells per wave-height (CPWH) in all directions) and numerical beach length (13 m) are selected-in this order-from a convergence analysis and a compromise between speed and accuracy.The least steep event is used, since it is the least nonlinear case that requires a finer vertical resolution to minimize the numerical diffusion.A long enough NWT (here 60 m long) is used to avoid completely any wave reflections, because, in the first 20 s, the wave-group has not reached the far wall.In 3D cases, the NWT is made symmetric with a symmetry plane along the length (i.e., wave direction).Mesh resolution is maintained as compared to wave-only cases in the MWL region containing the free-surface.Using the octree refinement strategy [52] (snappyHexMesh), the rest of the mesh is made three times coarser, apart on the surface of the rigid-body (here the X-MED buoy), which is one level finer than at the MWL.The X-MED buoy has a 66.3 N/m linear spring function of the distance between buoy bottom (i.e., attachment point) and seabed acting as mooring.As for the experiment, the −20.5 N pre-load in the mooring line is found to sink the buoy from its unmoored to moored resting position.
Mesh is deformable to account for the rigid-body motion.A RBM library solves the rigid-body motion and deforms the mesh accordingly.rigidBodyDynamics library is selected over sixDoFRigidBodyMotion RBM one, since sixDoFRigidBodyMotionRBM solutions are dependent on the user-specified deformable area; and because rigidBodyDynamics demonstrated good performance in the reproduction of extreme motion responses of the CETO PA [40].
The width of the symmetric NWT is set to 1.5 m.The motion responses of the X-MED buoy under the least steep event converge (R ≥ 0.9995 for all DoF) to those predicted by a 4 m wide NWT for half the computational effort.

WaveDyn: Cummins Equation
WaveDyn is a commercial code that was developed by the consultancy DNV-GL.WaveDyn solves the Cummins Equation [53], which is the conservation of momentum assuming the small displacement of a structure in an incompressible, inviscid, and irrotational flow (i.e., potential): where x is the vector body displacement, ẋ the velocity, and ẍ the acceleration; m is the rigid-body mass and m a the added mass; f hs (t) is the hydrostatic force, f e (t) is the excitation force, f ext (t) gathers additional forces that can be applied on the structure (e.g., mooring, PTO, mechanical joints); and, k(t) is the body impulse response function obtained from potential flow solvers (e.g., WAMIT [54], NEMOH [55]), along with the added masses and excitation force.As detailed in Table 5, the X-MED model is represented as in OpenFOAM with a 66.3 N/m linear spring to model the mooring line and joints to provide the 3 DoF, surge, heave, and pitch.Hydrodynamics coefficients were calculated while using WAMIT at the centre of gravity of the moored case using a wetted surface described by 1214 panels of 0.025 m size.The excitation force and body impulse response were imported into WaveDyn using the fine 0.01 Hz frequency sampling, since the WaveDyn computation times are negligible when compared with OpenFOAM.WaveDyn imposes the surface-elevation from a linear superposition of the measure of the WG aligned with the centre of the buoy at initial position.−20.5

Decay Tests Calibration
The heave and pitch decay tests are used for the calibration of each numerical model RBM.Decay responses are assumed to follow an exponentially damped oscillatory motion defined as: where f p is the resonance frequency; A is the amplitude of the first oscillation; T τ is the decrement; and, θ a phase constant.Only the first 10 s of the heave decay time-series are considered in order to avoid the effect of radiated waves reflecting on the physical tank sidewalls.Specifically for decay tests, the NWT width is set to the width of the physical one in order not to affect decay responses [30].
Only heave decay experiment time-series are available.A comparison of pitch decay is restricted to resonance frequencies.Table 6 shows that the two numerical models are correctly calibrated according to the decay tests.The lower accuracy obtained for the heave unmoored case might be partly explained by inaccuracies within the experimental set-up.As the heave motion is not driven by the mooring line, a small inaccuracy in the set-up might more significantly affect the response than for the moored decay test.Figure 8a shows that the reproduction of the four extreme events and the prediction of X-MED buoy resulting motion responses is accurately performed by OpenFOAM for all DoF (R > 0.9), while WaveDyn shows concerning inaccuracies in surge responses (R < 0.4).

Numerical Validation
Heave motion responses are more accurately reproduced by WaveDyn (--) when compared to OpenFOAM (----), as shown by Figure 8a.Below the breaking limit, the heave motion responses are mainly dictated by the surface-elevation.WaveDyn uses the linear superposition of a WG measurement aligned with the buoy centre at the initial position, while OpenFOAM tries to reproduce the propagation of a non-linear wave-group from a linear superposition at the position of the most upstream WG imposed at the inlet, which inevitably introduces inaccuracies, compared to the physical measurements, as the linear description does not properly account for nonlinearities that are present in the wave-group.It also means that OpenFOAM's surface elevation solution at the position of the buoy appears to be weaker when compared to WaveDyn, as the simulation of the wave propagation is imperfect.Figure 8a shows that the difference in the accuracy of the surface-elevation between OpenFOAM (----) and WaveDyn (--) explains WaveDyn's success in predicting heave motions well.
OpenFOAM predicts pitch motion with a higher level of accuracy (OpenFOAM R ≥ 0.95, WaveDyn R > 0.87-Figure 8a) due to the damping of the oscillations induced by viscosity.The inviscid assumption in WaveDyn also affects the capture of surge response, since the viscous drag force is omitted.In the experiment, Figure 7, the surge response shows a motion of low frequency and large amplitude (from 11 s to the end) when the buoy returns to equilibrium after reaching the extrema due to the passage of the main crest.This motion is not captured by WaveDyn, but it is captured by OpenFOAM.The formulation of drag force that is proportional to the square of velocity generally used in potential models is expected to generate higher frequency response.The assumption of small displacement-intrinsic to the linear potential flow solver-is deemed to be responsible for WaveDyn's inaccuracy in surge.
Similar levels of accuracy are obtained for the ST1 − f ocused events, as shown in Figure 8b (motion responses time-series are available in details in Musiedlak [30]).This second set of experiment confirms OpenFOAM accuracy (R > 0.9) in each DoF for all events; and, WaveDyn concerning predictions of surge responses (R < 0.3 (-×-)).

Assessment of Numerical Models
Trends in numerical models accuracy are drawn using both data-sets, and each model is assessed for a given DoF.
Below the breaking limit, heave responses are dictated by surface-elevation.Therefore, for all present cases, the accuracy of heave prediction is ensured by a good representation of the surface-elevation.WaveDyn's accuracy in heave is not expected to decrease significantly until the event cannot be described by a linear superposition, that is, above the breaking limit.In Figure 8a,b, OpenFOAM's heave prediction (----) directly follows surface-elevation accuracy (----).The three least steep events show that small improvements in the surface-elevation lead to larger improvements in heave responses.However, this is not the case for ST1 − f ocused events (Figure 8b), where surface-elevation accuracy varies with heave accuracy.Accordingly, provided that the surface-elevation is correctly captured (as for ST1 to ST3 but not ST4), the accuracy of heave response increases with wave-steepness.Both of the numerical models are validated for extreme heave assessment for kA ≤ 0.263.Breaking events will be required to differentiate the models in heave prediction.
For all events, the pitch responses look similar and they consist of a first oscillation, followed by the damping induced by viscous effects.Increasing event steepness increases the amplitude of the first peak and the damping rate.WaveDyn's accuracy decreases as steepness increases from kA = 0.179 (i.e., ST2) to kA = 0.263 (i.e., ST4).Hence, for kA = 0.179 (i.e., ST2), viscous damping effects probably start to become too significant for the inviscid assumption within WaveDyn.On the other hand, OpenFOAM accuracy increases with steepness, as viscous effects become more important provided that surface-elevation remains correctly reproduced.That is why OpenFOAM pitch accuracy drops for ST4 event and seems unrelated to surface-elevation for ST1 − f ocused events, probably because viscous effects are of similar magnitude for these twelve events.
For ST1 − f ocused events, WaveDyn accuracy oscillates from R = 0.95 to R = 0.85 depending on the capture of the natural frequency response of the oscillations, rather than the amplitude or the phase-shift between signals.A bias induced by the correlation R. Since all ST1 − f ocused events have the same event steepness, it seems that WaveDyn accuracy in pitch is more influenced by the shape of the crest rather than event steepness.In conclusion, WaveDyn remains validated for extreme pitch assessment as OpenFOAM.
Surge response is crucial in determining the overall extreme motion response assessment, because it can introduce an incorrect estimation of the rigid-body position in space.The experimental repeatability tests demonstrate that relatively small changes in the wave body interaction can result in significant differences in body motion and load after only a few periods [32].For a survivability assessment using an irregular sea-state, the error in the rigid-body position is likely to induce a different response (i.e., an extreme event not leading to extreme motion or, a moderate event generating extreme motion).Accordingly, surge motion is a crucial indicator of the usability of numerical models when dealing with long irregular sea-state.Therefore, despite the correct assessment of heave and pitch responses, the WaveDyn validity limit is defined below the least steep event, hence kA < 0.152.
OpenFOAM demonstrates the correct prediction of surge response for all events, hence validating the model for extreme motion response assessment.However, OpenFOAM surge accuracy still has room for potential improvement.Interestingly, surge accuracy inversely follows surface-elevation accuracy.In ST events, the over-estimation of surge motion increases with non-linearity probably as the laminar assumption tends to over-estimate drag.The increase in accuracy for ST4 is due to the overall similitude in shape and it identifies a limit in the use of the correlation coefficient to assess motion responses.For ST1 − f ocused events, surge motion accuracy decreases with surge maximum amplitude that decreases when the wave-group crest gets smaller.This confirms the complexity of representing cases of smaller wave-groups, where a given error will have larger proportion.Additionally, this inverse trend between surface-elevation and surge suggests that focusing purely on the quality of the surface-elevation prediction might not improve surge response, since the flow field below the surface might be of greater importance.

Deviation Points
The fact that WaveDyn exhibits concerning inaccuracies in surge motion responses is taken to identify the point where WaveDyn solution deviates from experiment, called the deviation point.An indicator function of a time-dependent wave-parameter is used to track the solution validity during the simulation.
In the following, the Root Mean Square (RMS) of the error between WaveDyn (x w dy) and the experiment (x e xp) surge responses up to time t would be used.It is normed by the experiment maximum surge, hence called NRMS(t), and it is defined as: where, N is the number of time-steps up to time t.At the deviation point, the instantaneous error between WaveDyn solution and experiment becomes suddenly large and remains large for a period.At this point, the accumulation of error increases drastically, generating a step in NRMS(t) time-series (left hand side Figure 9).The lower angle of the step identifies the deviation time; i.e., WaveDyn last point of validity.The angle generates a peak in the NRMS second derivative, d 2 RMS dt 2 , of amplitude related to the severity of the angle, as shown in Figure 9.At the time of the peak, wave-parameters deemed to be responsible for WaveDyn's deviation are identified.A low-pass butter-worth filter with a 4 Hz cut-off frequency is used to avoid high frequency noise due to small changes in error.ST cases are used to assess numerical models prediction compared to event steepness kA.A time-dependent parameter replaces the event steepness kA to track steepness at each time-step; the variation of the surface elevation relative to time: η is a time-dependent non-linear wave-parameter related to wave steepness.η is opposite in sign to the derivative with respect to space, ∂η ∂x ; but they have similar behaviour in absolute value.In potential flow theory, η is the vertical velocity of the free-surface at a given location in space (i.e., w z=0 ), which is proportional to the horizontal acceleration (i.e., ∂u ∂t | z=0 ) in the linear limit: Because the WaveDyn model ignores viscous effects, the WEC motion is dominated by inertia forces (radiation, diffraction calculated within the model).Hence, evaluating the surface-elevation vertical velocity-η -is equivalent to evaluating inertia forces acting on the WEC in the WaveDyn model.
Because WaveDyn deviates from experiment data for all studied cases, 16 points of deviations of different amplitude are identified with their corresponding steepness at that time (i.e., η ).The causality on deviation is assessed by the coefficient of determination-R 2 -between NRMS second derivative peaks and identified steepness at that time.Figure 10   Minor peaks in deviation are added to the analysis to consider points where the error increases, while the WaveDyn solution remains valid (i.e., non-deviating points); vertical grey line (|) in Figure 10.The evolution of η may be characterised as a function of deviation.The analysis is limited to the three largest peaks before the maximum (included), as smaller peaks will add points very close to zero.This methodology is not applicable after the maximum because surge motion response is no longer predicted well.ST1 − f ocused events improve the assessment, as they identify deviation points of lower severity.A total of 48 points are identified while using the present methodology.

Assessment of WaveDyn Validity Limit
Figure 11 shows the evolution of the absolute value of η identified at deviation ( ) and non-deviation points (×, ) as a function of the severity of deviation (i.e., second derivative of RMS).A linear (--) and a logarithmic fit (-) forced at the origin try to quantify the evolution.As expected, Figure 11 confirms that WaveDyn deviates for large values of |η |.The overall evolution seems more logarithmic (R 2 = 0.79) than linear (R 2 = 0.68).Accordingly, for low |η | values, a small increase in |η | will result in a small increase of deviation; whereas, for large |η | value, the same increase in |η | will result in a large increase of deviation.This justifies the existence of a threshold in |η |.Two asymptotes (--, --) of the logarithmic fit-respectively, at the origin and at large deviation value (0.3)-in Figure 11 highlight this difference in the rate of increase.Accordingly, as long as |η | amplitude remains under the limit, an increase (in |η |) results in a low increase of deviation, which means that the solution should remain valid.Above the given limit, deviation increases significantly for a small increase in |η |.Therefore, this limit is the validity limit of WaveDyn.The intersection of the two asymptotes defines the time in |η | from which the deviation changes from a slow increase to a sensibly larger one: WaveDyn's validity limit is |η | < 0.549 m/s.

Test of WaveDyn Validity Limit
Table 7 presents the use of the WaveDyn validity limit (i.e., |η | < 0.549 m/s) for all cases as compared to the actual deviation, as previously identified by the second derivative maxima: it is correct if the limit identifies the deviation peak, restrictive when the limit identifies the second largest deviation peak, and late when the limit identifies a time after deviation.Table 7 shows that WaveDyn validity limit is mostly accurate but restrictive and sometimes late.The applicability of the current limit to precede WaveDyn deviation is deemed correct since the limit is restrictive.The current limit is suitable for identifying WaveDyn deviation preceding a large surge low frequency response that will be not captured by WaveDyn, because the limit assessment was conducted using only such cases.From Figure 11, we could have predicted that the limit would be sometimes not accurate at identifying the deviation peak.The three cases where limit is late are the cases that generate deviation of low amplitude where the experiment surge low frequency response has the smallest amplitude.These deviation cases are the hardest to identify, since the closest to the limit of WaveDyn deviate or not.A limit specified as the |η | value obtained for the deviation point of lowest deviation amplitude (i.e., |η | < 0.38 m/s at 0.07 deviation) would be too restrictive when applied to other cases.
The limit in |η | is now applied to an extreme irregular sea-state case not used for its assessment.Figure 12 compares surge motion responses for the experiment (-) with WaveDyn (-).The vertical red line (|) marks the first value of |η | above the limit, while (-) is the deviation.
Figure 12 shows that the time identified by WaveDyn validity limit just precedes a large surge response of low frequency.Therefore, this confirms the applicability of the present limit to predict WaveDyn deviation preceding a large surge low frequency response that will be not captured by WaveDyn.It also confirms that large |η | values do generate a large surge low frequency response.
Figure 12 also shows that the limit in |η | is late compared to the estimation of the maximum of deviation which occurs around t = 20 s.This deviation generates a large peak because NRMS is calculated from t = 0 s to the considered time; hence, the proportion of increase of cumulative error is significant at that time.It is the same with the peak just following, and peaks around 30 s.At the identified limit, the rate of increase is slower since it is more constant due to the low frequency response.Accordingly, this example indicates an issue in the use of the methodology for long time-series, and it suggests that a second validity limit can be investigated to characterise these deviation points.A preceding study [31] develops the hot-start process for a heave decay case and a wave-only case-taken separately.In each, CFD simulation is initiated at an advanced time with either the wave-field or the RBM being in motion.These solutions are compared to a conventional start (i.e., still water or zero rigid-body motion and acceleration).This procedure is called 'hot-start' and it aims to obtain the same accuracy as a conventional start while reducing the computational cost.The present study merges these two developments by extending the hot-start procedure to an extreme WSI simulation.Subsequently, the coupling is developed using the WaveDyn validity limit.
For wave-only hot-started cases in previous study, the wave-field is imposed along the NWT from the linear superposition of harmonics used at the inlet, and builds-up towards the conventional solution.The results show that building up 4 s prior to the crest of the NewWave events (ST1, ST2, ST3, and ST4) is sufficient to reproduce (R > 0.999) surface-elevation, water column velocity, and pressure at model initial position of the same wave-field starting conventionally.This accuracy is considered to be related to the representation of a short wave-group, since this builds-up in a few seconds preceding the event and so the imposed linear solution is already close to the actual one.Applicability is questionable as the overall event length increases, but it is expected to remain accurate for NewWave event of larger peak frequency, hence f p ≥ 0.356 Hz.
For the hot-started heave decay test in previous study, RBM-state (position, velocity, and acceleration) is imposed from a previously run CFD simulation, while the wave-field is set as still water.A library called deformDyMMesh [31] allows for deforming the mesh to position the rigid-body at the 6-DoF of hot-start, hence maintaining mesh quality when rigid-body returns to equilibrium.RBM was found to require 5 build-up time-steps (i.e., a time-step where RBM solution is corrected) to impose the correct motion to the buoy.
Therefore, in order to merge these two previous developments, the event is now hot-started at the focusing time, t hot = t f ocus , while using the 4 s build-up period (i.e., t minus = 4, s) from a linear solution; and, RBM-state is imposed from the previously run simulation using 5 build-up time-steps.The least steep event is used as a demonstrator.Steeper cases are believed to generate similar results, because the 4 s wave-field build-up period was found to accurately reproduce any of the steeper events (i.e., up to ST4).
Figure 13 shows surge (a), heave (b), and pitch (c) responses for the conventional simulation (--), and the hot-started one (--).The agreement obtained (R > 0.997 for each DoF) demonstrates the reliability of the hot-start technique and validates its parameters (i.e., t minus = 4 s, and 5 build-up time-steps for the RBM) for NewWave events of peak frequency f p ≥ 0.356 Hz.The good agreement is probably due to the hot-start happening when rigid-body motion is very small.Diffracted and radiated waves are negligible at that time, and wave-group profile is mainly linear.Similar agreement is expected for the steeper cases for the same reasons.However, the application to an irregular sea-state needs further investigation, with the definition of a build-up period function of the previous wave-train.Further investigations where the rigid-body motion is large are necessary in order to assess the effect of neglecting radiation and diffraction.Because the non-linearity of ST1 is small, the wave-field build-up period can be reduced while maintaining the accuracy of the generated wave-field.At these times nearer to the focus, the rigid-body motion is large and radiated/diffracted effects are of more importance.Hot-starting the simulation at that time might be used to assess the radiated/diffracted effects.
All coupling prerequisites are now fulfilled.The following section describes the use of a WaveDyn solution to hot-start OpenFOAM and achieve coupling.

Coupling of RBM: Moored Heave Decay
The heave decay test is first used to develop the coupling around the RBM.The hot-start time, t hot = 0.1 s, is chosen, because the radiated wave-field is negligible at this point limiting the error in initializing the NWT as still water at hot-start.There is no build-up period for the wave-field (t minus = 0 s) and, therefore, the coupled model imposes WaveDyn solution at t hot = 0.1 s.Besides, WaveDyn and OpenFOAM heave solutions are identical at this time.Therefore, the fully coupled solution should converge to OpenFOAM solution.
Figure 14 illustrates the change from the WaveDyn only simulation (-), to the fully coupled-model (-), which matches the OpenFOAM only solution (--).It demonstrates the correct definition of the ELC and the capacity to hot-start the RBM in OpenFOAM from a WaveDyn solution.The 0.1 s of WaveDyn simulation saves 16% of the computational effort of the OpenFOAM-only simulation.This saving is considered to be due to the time-steps size ruled by the Courant number, which generates small time-steps at times close to 0 due to the high acceleration of the rigid-body.

Deterministic Case: ST1
The coupled model is tested against the ST1 case hot-starting at focusing time, t hot = t f ocus , with a build-up period, t minus = 4 s.OpenFOAM hot-starts with the RBM-state evaluated by the WaveDyn-only solver at t hot − t minus .The coupled model imposes the OpenFOAM solution from that time (i.e., t hot − t minus ) into WaveDyn.
In Figure 15, the coupled-model (-) differs from the WaveDyn only model (--) after t hot − t minus , from where it closely follows the OpenFOAM only solution (--).WaveDyn and OpenFOAM-only RBM states at t hot − t minus are very similar for surge and pitch (around 0.1%), and slightly different (3%) for pitch.The small difference in pitch can explain the difference between OpenFOAM-only and the coupled model in Figure 15.The difference in surge, especially the one appearing after t = 16 s, might be due to differences appearing in the wave-field towards the end of the simulation; or that the small differences in surge builds-up, making surge response very sensitive to the RBM at t hot − t minus .Differences remain negligible, so the coupling between WaveDyn and OpenFOAM is considered to be demonstrated.
The coupling is also validated against experimental data for the assessment of extreme response.The accuracy of the coupled model in each DoF, measured by the correlation (i.e., R) between time-series, is similar to the CFD only: surge is 0.941 against 0.956; heave is 0.919 against 0.92; and, pitch is 0.974 against 0.95.The benefit in using the coupled model results in a 26% reduction in computational effort.Therefore, the coupled model demonstrates a level of accuracy comparable to CFD for extreme motion response, while minimising the computational effort, hence demonstrating the concept.
This result is deemed to be applicable to the three steeper events, ST2, ST3, and ST4.The wave-field was found to be accurately reproduced while using the 4 s build-up period.RBM-states are similar between OpenFOAM and WaveDyn four seconds prior to the focusing of the event (Figure 7).ST1 − f ocused events are also expected to be reproduced up to a similar accuracy, since their wave-fields are based upon the ST1 event (i.e., same peak frequency f p).
It should be noted that the WaveDyn validity limit identified by |η | limit is at the trough before the main crest (just before 10 s Figure 10).If the hot-start time would have been chosen by WaveDyn validity limit, the coupled model would have swapped to OpenFOAM at an earlier time that done for Figure 15 (just before 6 s due to t m inus = 4 s).This is equivalent to the coupled hot-starting at focus where a longer wave-field build-up period is used (i.e., t minus > 4).Because the wave-field converges using t minus = 4 s, it converges for larger values [31].Accordingly, it is expected that very similar motion responses would have been obtained in using WaveDyn validity limit, hence illustrating the full use of the coupling to model extreme WSI.

Long Irregular Sea-State
The WaveDyn validity limit, |η | < 0.549 m/s, is implemented as a trigger in the ELC and applied to the long irregular sea-state previously used.The WaveDyn validity limit is identified for the long irregular sea-state at t hot = 32.86 s.Using a 4 s build-up period for the wave-field; hence, at t hot − t minus = 28.86 s, OpenFOAM starts from WaveDyn RBM solution using 5 build-up time-steps.The coupled model simulation stops at 40 s, since it is not realistic to run the entire remainder of the simulation with OpenFOAM.
Figure 16 compares the motion of the X-MED buoy under an extreme irregular sea-state between the WaveDyn-only (--) model, the coupled model (-), and the experiment (--).The time where the WaveDyn validity limit is attained, t hot , and the time where CFD is used, t hot − t minus , are highlighted as vertical green lines (|).
The coupled model performs a simulation not industrially realistic in CFD only due to the significant resulting computational cost.In this simulation, the mid-fidelity model WaveDyn is identified as outside its validity range by the previously developed limit.The coupled model solution improves the prediction of extreme motion responses.WaveDyn-only surge response does not reproduce the low frequency large surge motion happening from 32 s to 40 s, while the coupled model does.In addition, pitch response is correctly captured after the validity limit by the coupled model.Heave response accuracy decreases for the coupled model as for OpenFOAM only simulation for ST and ST1 − f ocused events due to the difference in wave modelling.In overall, the coupled model demonstrates its interest to model accurately long extreme irregular WSI.

Conclusions
In this study, the mid-fidelity time domain potential flow solver WaveDyn is weakly coupled while using a time-splitting strategy to the CFD-NWT validated for extreme WSI developed in OpenFOAM-4.1.A new time-splitting strategy is introduced, which switches between models at times, depending on WaveDyn's validity limit, caused by inviscid and small rigid-body displacement assumptions.During the coupling, WaveDyn receives the OpenFOAM RBM solution via an additional load output.The coupled model uses the computationally expensive but more accurate OpenFOAM model when necessary and the less expensive WaveDyn model when this one is valid.This coupling strategy overcomes WaveDyn's validity limit while reduces computational effort in comparison with CFD only and maintains accuracy.
Both numerical models are first assessed against experiment data for a PA-like buoy motion responses under NewWave extreme events of increasing steepness, and of various shapes.OpenFOAM is validated for all events and all DoFs; hence, for kA ≤ 0.263.WaveDyn deviates from experiment in surge since it does not capture low frequency large responses due to linear potential flow assumptions.Even so, heave and pitch responses are well captured, WaveDyn is not validated for extreme WSI because of the importance of rigid-body surge position on following responses for dynamics structures.
Times where WaveDyn surge response deviates from experiment-due to a large low frequency response not captured-are used to identify WaveDyn's validity limit.The vertical velocity of the surface elevation (i.e., derivative with respect to time, η = ∂η/∂t) is introduced as the wave parameter relating to instantaneous steepness and inertia forces since WaveDyn ignores viscous effects.Large instantaneous steepness values demonstrate a relation with large low frequency response.Therfore, WaveDyn's validity limit is assessed at |η | < 0.549 m/s.Above the limit, WaveDyn surge response is expected to deviate, meaning that CFD is required to accurately capture all DoF.
Finally, the proof of concept of the coupling is demonstrated.OpenFOAM starts when WaveDyn validity limit is estimated as attained, from WaveDyn RBM solution (i.e., position, velocity, and acceleration) and the imposed linear superposition for the wave-field (developed in [31]).Coupling demonstrates interest using the least steep event by reducing the computational effort by 26% while maintaining accuracy to the level of CFD-only.A result considered to be applicable to the other events.A long irregular sea-state is used to indicate the coupling potential for long simulations that would be prohibitively expensive for CFD alone.

Figure 2 .
Figure 2. Time-splitting strategy as function of time.
5 N. Its dry-mass is 43.2 kg.Resonance frequencies are obtained from experimental decay tests, and are 0.926 Hz for unmoored heave, 0.917 Hz for moored heave and 0.75 Hz for moored pitch.Motion-responses are recorded by optical tracking targets of a ±0.5 mm residual in precision.

Figure 3 .
Figure 3. Extreme loading of marine energy devices due to waves, current, flotsam and mammal impacts (X-MED) model [32].

Figure 5 .
Figure 5. Surface elevations time-series for four of the twelve ST1 − f ocused events.

3. 1 . 1 .Figure 7
Figure 7 shows physical and numerical surface-elevation (d) and motion responses-surge (a), heave (b), and pitch (c)-time-series.Figure 8a illustrates the evolution of accuracy quantified by the correlation (R = cov(A, B)σ A σ B, where A and B are two random variables[56]) as a function of the measured event steepness, kA, for each time-series: surge (×), heave ( ), pitch responses (•), and the surface elevation ( ).OpenFOAM is represented using coloured dashed lines (--) and WaveDyn coloured plain lines (-).WaveDyn surge motion accuracy (-×-) is set separately to highlight the difference with other DoF.

Figure 7 .
Figure 7. Surface-elevation and motion responses time-series of the X-MED buoy under ST events.

Figure 8 .
Figure 8. Correlation (R) between each numerical model and experiment.

Figure 9 .
Figure 9. Evaluation of the lower angle in the step of the NRMS time-series.
illustrates the methodology applied to case ST1, where the vertical black (|) line indicates the deviation point where |η | amplitudes are identified.

Figure 10 .
Figure 10.Identification of wave-parameters at deviation point for ST1.| is the deviation point, | is the second largest peak (i.e., a non-deviating point).

Figure 11 .
Figure 11.Evolution of |η | identified at deviation and non-deviating points for all events.

Figure 12 .
Figure 12.Estimation of the deviation point using the WaveDyn's validity limit (| identifies |η | < 0.549 m/s) applied to an irregular extreme sea-state.

Figure 13 .
Figure 13.CFD motion responses time-series when conventionally started and when hot-started at focus (Hot-started with t minus = 4 s and 5 time-steps for RBM build-up).

Figure 15 .
Figure 15.Comparison of numerical models motion responses prediction time-series using ST1 case.| indicate hot-start time (t hot ) and when coupling starts (t hot − t minus ).

Figure 16 .
Figure 16.Comparison of motion response time-series between experiment and numerical models prediction where the hybrid model is triggered by WaveDyn's validity limit.| indicate hot-start time (t hot ) and when coupling starts (t hot − t minus ).

Table 1 .
Wave-tank and X-MED details.

Table 5 .
Parameters of the X-MED model in WaveDyn.

Table 6 .
Decay tests resonance frequencies and decrements.