Mooring systems with submerged buoys: influence of buoy geometry and modelling fidelity

Abstract Mooring systems often make use of submerged buoys (SBs) in order to make the moorings compliant. In this paper we present the dynamic effects of changing the buoy geometry or the buoy model fidelity on the mooring system response. Three cylindrical SBs with increasing slenderness (height/diameter) are studied for a mooring leg with two polyester ropes and a SB. The results show a large impact of SB geometry on the mooring dynamics. A larger height/diameter ratio (with preserved mass and buoyancy) is shown to be beneficial as it gives both smaller tension force magnitudes and, more importantly, avoids slack-snap occurrence in the upper cable. We further present a comparison between four numerical methods for SB dynamics: (i) a high-fidelity model using computational fluid dynamics (CFD); (ii) the Morison equation with slender body drag force approximation using numerical quadrature; (iii) the Morison equation with an independent evaluation of the fluid drag due to translation and rotation; and (iv) a translating Morison model which simulates a vertical cylinder in three degrees of freedom with no rotation. All methods are used together with a high-order finite element mooring dynamics solver. The results show that the translating method is inadequate to model this mooring configuration. The remaining three methods agree moderately well, but the Morison formulations give larger motions and higher tensions compared to the CFD results. We show that the quadrature drag model is better suited to model the drag moment on SBs than the independent model, and that the improvement increases with increasing slenderness of the buoy. The uncertainty, sensitivity and importance of the hydrodynamic coefficients of the buoy are discussed and examined by a regression analysis from the CFD data.


Introduction
Chains are the traditional choice for slack-moored offshore structures in shallow waters, such as semi-submersibles and floating production storage and offloading (FPSO) units. The weight of the chain (kg/m) and its length governs the restoring stiffness of the mooring. But chains have some well known disadvantages, such as the large weight during installation and a large footprint [1]. For deep and very deep waters, the chain is installed together with sections of rope and/or submerged buoys (SBs) and clump weights (CWs) [2]. A SB is defined as a mooring element designed to increase the buoyancy of the system. It is usually made of polyethylene filled with polyurethane foam. Analogously, a CW is here defined as a mooring element providing negative buoyancy, commonly made of commonly made of steel or concrete. The emerging industry of marine energy has a requirement of reduced mooring footprint to facilitate farms of tightly placed devices. Thus for ocean renewable energy (ORE) applications, chain moorings are often replaced by taut moorings using ropes [3,4]. The rope stiffness is generally too high to provide the required mooring compliance, so SBs and CWs are then introduced to create the restoring stiffness of the system. Moorings of this type is often referred to as hybrid moorings [5], in which we also include multi-catenary mooring systems with chains instead of ropes between each SB/CW [6]. An example of an ORE application using a mooring system with submerged buoys is the Wa-ves4Power [7] full-scale device deployed at Runde, Norway (see Fig. 1), which will be the starting point of the analysis in this paper.
shown to have significant dynamic effects on the system and there was a strong frequency-dependence tied to the hydrodynamics of the buoys and their coefficients [10]. It was experimentally shown that properly placed SBs can give beneficial reduction of mooring tension for very deep waters [2]. A numerical investigation of different buoy sizes and masses [11] further strengthened the importance of the dynamic effects. The resulting dynamic tension decreased when buoys were included, compared with the pure chain mooring. A similar comparison was later done by Surendan and Goutam [12], where the addition of a SB on a catenary mooring leg gave an overall reduction in dynamic tension by 8-25 % for different cases of load amplitudes, sea current conditions, pretensions and buoy placements. Other parameter-studies include [13] in which SBs and CWs were placed on a taut mooring, and more recently [14] in which the effect of inserting a spherical buoy on the baseline catenary chain mooring legs of a semi-submersible platform was investigated.
In shallow water depths, hybrid mooring systems were studied for aquaculture in [15], for different surface-piercing buoy types. Hybrid mooring systems were also suggested by Fitzgerald and Bergdahl [5] for mooring wave energy converters (WECs). They showed numerically that catenary moorings gave very high loads during slack-snap behaviour and that the mooring loads from a taut-mooring with two lines and a submerged buoy gave significant reduction in the dynamic range of the mooring force. These results were very promising but the horizontal stiffness of the proposed SB-mooring was lower than the pure catenary, which made the system more sensitive to drift offset. The results of the different mooring configurations are therefore hard to compare [16]. Paredes et al. [1] studied three mooring leg types for a generic cylinder in an experimental campaign: (A) a SB with two lines; (B) the same as A but with a CW attached as well; and (C) a catenary chain. The moorings were in this case designed to have the same horizontal stiffness and the performance of the moorings were similar. The addition of the SB was not found to be as beneficial as reported in [5]. Vicente et al. [17] investigated the influence of buoy location for a slack moored WEC. It was shown that a buoy close to the surface and horizontally close to the anchor point yielded higher power absorption and smaller horizontal displacement of the WEC. Ortiz et al. [18] developed a surrogate model in order to optimize the mooring system for a heaving point-absorber with regard to power production. The size of the SBs was one of the optimized parameters.

Numerical modelling of submerged buoys
The motion response of the submerged buoy can be computed in several ways depending on the desired fidelity of the simulation. The most complete manner is, of course, to simulate the flow around the SB using CFD models solving the viscous Navier-Stokes equations. While moored submerged buoys have been modelled using CFD in the context of wave energy [19,20] this approach has not before been applied to SB/CW that are mooring elements. CFD comes with a substantial computational effort that generally is not justified as these components are considered to be secondary structures. On the other extreme is a quasi-static analysis of mooring lines with submerged bodies, as investigated in [21]. In this approach only the net buoyancy force of the body is included, which makes it suitable for rough screening design calculations only. The significant contribution of hydrodynamic forces [11] is neglected.
The general approach to modelling SBs/CWs is to treat them as Morison bodies [22]. This simplification is typically acceptable as the submerged bodies are small in relation to changes in the fluid velocity due to wave motion or ocean currents. For small bodies it is common to neglect the rotational modes and model a translating body only, see e.g. [13]. Mavrakos [11] first introduced the rotational modes of the buoys in mooring simulations by computing the quadratic drag force independently for the rotational and translational modes of motion of the SB, which was originally proposed for spherical buoys. When the buoy shape is cylindrical and its slenderness (height/diameter h/D) increases, the independent approach becomes increasingly inaccurate. Further, its accuracy is questionable for bodies experiencing significant wave loads or motions, which will be the case e.g. in intermediate water depth or in many moorings of wave energy devices. The short-comings of the independent drag modelling approach to SB dynamics motivates computing the Morison drag force along the cylinder with numerical quadrature, to consider the velocity variation along the height. With the computational resources available today, the increase in computational effort is considered to be minor in comparison with the increased accuracy and generality of the model. There are clear differences in how commercial software packages for offshore engineering treat SB/CWdynamics. Although there is always the option of meshing the buoy as a radiation-diffraction body, the standard mooring buoy elements use different model fidelity in different solvers. A buoy element in ANSYS AQWA [23] or DNV DeepC [24] is treated as a translating body without rotational effects, while a buoy in OrcaFlex can be modelled in six degrees of freedom but the drag is modelled using the independent assumption between translational and rotational forces. Finally a cylindrical buoy in ProteusDS [25] can be modelled with drag effects using numerical quadrature. The different approaches are of course suitable within a certain range of case specific conditions, however there are few studies investigating the effects of buoy model fidelity on the mooring reponse.

Paper contribution
This paper aims to investigate the importance of submerged buoy dynamics on mooring system response. We also analyse the numerical modelling fidelity required to predict the dynamic effects for mooring systems with submerged buoys. We present a comparison of four different methods: (i) CFD simulations, used to extend the modelling fidelity of SBs beyond the Morison approximation. We present unsteady Reynolds averaged Navier-Stokes (RANS) simulations of the submerged buoy. To the authors best knowledge, RANS simulations have not before been used for evaluating SB in the field of mooring analysis. (ii) quadrature Morison drag, where the drag forces and moments are computed from integration of the Morison drag over the SB using numerical quadrature. (iii) independent Morison drag, where the drag force is proportional to the square of the translating velocity and the drag moment is proportional to the square of the angular velocity. The total force and moment from hydrodynamic drag are thus computed independently. (iv) translating Morison body. This method neglects the rotational effects of the body altogether. It has the lowest fidelity of the methods compared.
The three Morison approaches (ii)-(iv) have been implemented in a rigid body library included in the discontinuous Galerkin (DG) mooring dynamics solver named Moody [26]. The first method (i) is run as a coupled CFD-mooring simulation using Moody's native coupling [27] with the rigid body solver of OpenFOAM [28]. Thus, the cable dynamics are solved in exactly the same manner in all four methods in order to isolate the effect of buoy model fidelity. The models are compared in terms of tension forces and body motions. The investigated test case is a simple mooring leg based on two lines and one submerged buoy under a range of cyclic load conditions at the fairlead. The mooring solution from Waves4Power's deployment at Runde, Norway [7,29], see Fig. 1, is the starting point of the test case analysed. We point out that the study is restricted to fully submerged non-spherical buoys that may be modelled as cylinders, and that some important changes have been made to the Runde trial case.
The paper is outlined as follows. Section 2 details the governing equations of the mooring cable and the complete equations of motion for a submerged cylinder. In Section 3 we explain the differences in buoy modelling fidelity applied in the paper. The detailed layout and properties of the mooring lines and cylinders used in the numerical investigations are presented in Section 4. The numerical settings used in Moody and in the CFD simulations made in OpenFOAM [28] are also outlined. The results are presented and discussed in two sections. Section 5 compares different buoy geometries, and Section 6 presents a comparison between the four modelling approaches. This is followed by conclusions in Section 7.

Governing equations
Throughout the paper we will with the superscripts c, f and b refer to cable, fluid and body, respectively. (1)

Mooring dynamics
where γ 0 is the cable mass per unit length, ϵ is the elongation of the cable, and ϵ˙the strain-rate. The cable tension force magnitude, T (ϵ, ϵ˙), contains the constitutive relation of the cable material. In the case of a linear visco-elastic cable the relation is = + T EA ξ ϵ ε 0 in which EA 0 is the axial stiffness and ξ is the internal damping coefficient. The variable → f c represents all external forces made up of: (i) added mass and Froude-Krylov forces; (ii) net force of gravity and buoyancy; (iii) contact forces, typically from sea-floor interaction; and (iv) drag forces acting on the cable. Eqs. (1) and (2) are solved with the Moody mooring solver. Moody uses an hp-adaptive discontinuous Galerkin method. For a full description of the method see [26].

Submerged buoy dynamics
SBs are most often used far from the surface to avoid the forces of the splash zone. Hydrodynamic forces may have a significant impact on the mooring response, but the buoy size is in general small in relation to changes in the fluid velocity. The Morison equation [22] is therefore a valid assumption and we will use it to derive the governing equations of motion for a submerged cylinder. Please note that only fully submerged buoys are discussed in this paper. The SBs are further assumed to be located deep enough so that no free surface waves are generated from the buoy motion, i.e. the radiation damping of the SB is assumed to be zero.

Coordinate systems
Let = x y z {^,^,^} be the inertial frame, and = k {î, ĵ ,^} denote the body fixed coordinate system with origin in the centre of mass, see in which where M b is the 6 × 6 mass matrix, We stress that → O is a position vector expressed in the inertial frame , while all other state vectors have components in the Lagrangian frame of reference . Expressing Eq. (4) in and separating the rotational and translational degrees of freedom gives where m b is the body mass, I is the 3 × 3 inertia matrix of the body and V b is the volume of the body. The fluid density is denoted by ρ f and → g is the acceleration of gravity. The expressions for added mass and drag damping are a bit more complicated. The derivation below is valid for a cylindrical buoy, with Lagrangian frame located in the geometrical center and k along the symmetry axis. We assume the size of the body to be small in comparison with changes in the surrounding flow field. The fluid velocity and acceleration are thus treated as constant over the body domain and evaluated at the origin of , i.e.
the added mass force and moment vector → F b a (6 × 1) is written as 0^ .
The index ⊥ denotes vector projection onto the î, ĵ plane and h b is the cylinder height. C M1 and C M2 are the coefficients of added mass perpendicular and parallel to the symmetry axis respectively. The factor h 2 /12 comes from integrating the moment equation along the cylinder. A complete derivation of the added mass matrix for the case of ballasted buoys ( → ≠ → OB 0 ) is presented in Appendix A. For drag forces, the quadratic term makes the resulting integral expressions exceedingly expensive to evaluate, and we therefore leave the expression in integral form as (13) where C D1 , and C D2 are the in-plane and out of plane drag coefficients of a circle respectively, C D3 is the sectional shear coefficient of tangential drag and D b is the cylinder diameter.
→ Q (1) and → Q (2) denote the integrals over the cylinder height as a function of the local section velocity We rewrite Eqs. (4)-(6) in terms of a state vector with 13 degrees of freedom For submerged bodies, the effective mass matrix M b e describing inertial effects of the body and the surrounding fluid is constant in the bodyfixed coordinate system. Further, if the centre of mass coincides with the centre of buoyancy M b e reduces to a diagonal matrix written as which are used as boundary conditions for any cable connected to attachment point → OM i .

Methods for submerged buoy motion
This paper compares four numerical methods for computing the dynamics of a submerged buoy (SB). In all methods, the SB is connected to the same mooring dynamics solver, Moody. Thus, the difference between the methods is isolated to the treatment of the rigid body motion of the buoy itself.

CFD Simulations
The CFD results in this paper come from viscous Reynolds averaged Navier-Stokes (RANS) simulations, which here represent the model with the highest fidelity for SB dynamics. The rigid body library of OpenFOAM (v1806) [28] is used to simulate the buoy response in the RANS interDyMFoam solver. The force and moment from the fluid acting on each computational face of the body surface is computed from the fluid pressure p i and the shear stress τ i of the surrounding flow. Thus, Eq. (4) becomes where r i is the position vector from the centre of gravity to the cell facecentre or the mooring attachment point respectively, n i is the unit outward-pointing normal of face i and A i is the face area. The total number of cells on the body surface is N cell . The mooring force computed by Moody acts as a restraint to the rigid body solver in the CFD domain. For details and validation of the coupling see [27,30]. The connection with Moody was made using the interpolation scheme of Moody API [31]. The typical time step size for the mooring cables in these simulations is at least an order of magnitude smaller than the time step used in the fluid part. A quadratic interpolation of attachment point position was used to sub-step the position boundary condition of the mooring solver. Please note that the mooring cables exist in the mooring solver only, and that interaction between the CFD fluid state and the mooring cables is neglected. The coupling is between the mooring point position in the CFD domain, and the net mooring force at the attachment point in the mooring solver.

Quadrature Morison drag
This method simulates a six degree of freedom submerged body with the Morison equation, as it is written in Eq. (4). Please note that the integrals in the drag force (Eqs. (12) and (13)) are computed with numerical quadrature using a seven point Gauss-Lobatto-Legendre quadrature (see e.g. [32]) on a 5th order Legendre polynomial. For cases when the buoy is both translating and rotating, the quadratic drag force gives a net drag on the buoy due to its rotation, and a net moment due to the translation. Hence, it will resist rotation more strongly than in the independent Morison formulation. That said, we are still subject to the limitations of strip theory. The sectional drag coefficients are from standard tables [33] and are typically obtained from translational tests, while they are used for rotational motion as well. The true forces from the quadrature Morison drag are more complicated, especially for short and fat cylinders, as was pointed out in [25].

Independent Morison drag force
This method computes the drag force of the Morison equation [22] independently for the rotational and translational degrees of freedom. Originally proposed by Mavrakos et al. [11], it differs from the quadrature Morison drag approach only in the evaluation of → Q (1) and → Q , (2) see Eqs. (12), (13). Due to the independent assumption, the integrals where is the relative translational velocity of the buoy.

Translating Morison body
We also simulate the buoy response using a translating body approach. The buoy dynamics are thus reduced to three degrees of freedom, for motion in x, ŷ, and ẑ respectively. We include it here to isolate the effect of the buoy rotation. In this paper we simulate Eq. (4) with the additional constraints to remove the rotation of the cylinder. Please note that typical threedegree-of-freedom bodies are defined with single values of C M and C D [23,24], while we here keep the initial configuration of the cylinder and its differentiation between vertical and horizontal hydrodynamic properties.

Case description
We analyse a single mooring leg consisting of two mooring lines interconnected with a relatively large SB. The geometry is shown in Fig. 3. The mooring leg is adapted from a full-scale installation of the Waves4Power buoy at Runde, Norway, with details presented in Yang et al. [34], which was further studied in Lang et al. [29].
The mooring leg consists of one buoy and two ropes of polyester type (see Table 1 for material data) in 90 m water depth. As is depicted in Fig. 3a, Line 1 (L 1 ) extends from the anchor point A to the buoy connection point B, and Line 2 (L 2 ) connects point B to the fair-lead point C. Point F represents the centre of gravity (cog) of the buoy,  Table 2 for details.
coinciding with its centre of buoyancy. The mooring force at point B is taken from the end-point ( = s L) of Line 1, and the force at point C is sampled at = s L of Line 2, being the fair-lead position.

Buoy characteristics
The original submerged buoy used in the prototype demonstration of the Waves4Power WEC consisted of a cylinder diameter of  Table 2. We will refer to the cylinder with = D 2.1 m as the baseline buoy (D 0 ). The different buoy geometries are each computed from the diameter, D, and the volume, V 0 , as where h tot is the height of the buoy with the cones, h eff is the height of the cylinder, z B is the z-coordinate of the mooring attachment point B, V c is the volume of the cone and = k (0.937/2.1) is the cone slope, which is also assumed constant over the different buoy diameters. Please note that the mooring attachment of the baseline test is [34]. This is a consequence of the conical lid assumption. For the purpose of this generic study, the difference is acceptable. Further, note that the buoy centre (point F) will be placed differently for each buoy, because of their different heights. Due to the preserved buoyancy however, the mooring lines and point B will remain at the same initial equilibrium location as shown in Fig. 3 for all buoy variations.
Making a suitable choice of hydrodynamic coefficients for added mass and drag requires some special attention. When the body is far below the free surface, the effect of radiation damping can be neglected [35,36] and the added mass coefficient is therefore constant across the frequency range. This assumption is valid in this case study. The hydrodynamic mass (added mass) per m (i.e. sectional coefficient), of an infinitely long cylinder is 1.0 [37]. However for shorter cylinders, end effects play an important role [33]. In this work we use the tabulated formulas from [33] as a basis of choosing the C M1 coefficient. In the vertical direction (axis of symmetry of the cylinder) the coefficient is not listed in [33]. Therefore, we apply the empirical correction formula derived in [38] to get an estimate of C M2 as a function of the slenderness h/D. The sectional added mass coefficients chosen for all three buoys are presented in Fig. 4a. The drag coefficients are also taken from the DNV recommendations in [33] (for Re ≈ 1 × 10 5 ) and are presented for different h/D in Fig. 4b. The drag coefficient of a cylinder is relatively constant in the range of Re ∈ [10 3 , 10 5 ] which is the expected working range of the buoy. The hydrodynamic coefficients presented in Fig. 4 are used in Eqs. (10) and (11) to compute the total added mass force and drag force on each cylinder respectively. We will use the coefficients chosen for each buoy as constant throughout the simulations.

Load cases
To highlight the importance of the dynamic effects of the submerged buoy on the resulting mooring tension, we will look at the mooring response using a prescribed circular fair-lead motion of small amplitude = a 0.5 m. This amplitude represents a very small excitation in relation to the water depth ( = a h / 0.55 %) and is in a wave-load context corresponding to a wave-rider moving with a = H L / 2.5 % steep wave at

Table 2
Submerged buoy properties of the three different diameters. All buoys have the same net buoyancy. The cylinders are labelled by their indexed diameter, D 0 , D 1 and D 2 respectively.  Fig. 4. Hydrodynamic (sectional) coefficients for the different buoy designs studied. The DNV curves were interpolated from data in [33] and the h/D-dependence formula for C M2 from [38]. = T 5 p s, which is close to the resonance frequencies of the mooring system. We chose a small excitation amplitude to isolate the nonlinear response of the hybrid mooring system and the SB motion. The investigated periods are T p ∈ [1, 10] s for the comparison between floaters of different geometry, and T p ∈ [2, 10] s for the model fidelity investigation. The fairlead oscillation amplitude is increased from zero over 2 loading periods using a cosine ramp. All simulations were made in still water conditions, to maintain comparability between the methods. Any hydrodynamic effect on the SB dynamics is thus exclusively triggered by its own motion response.

Numerical settings
The mooring lines L 1 and L 2 were discretised in Moody using a uniform polynomial order of = P 4 with 5 and 10 elements respectively. The high-order formulation provides sufficient resolution of the dynamics despite the seemingly low element count [26]. A CFL-condition of 0.5 was applied to control time step size for an explicit third order strong stability preserving Runge-Kutta scheme. A mesh sensitivity analysis showed that this setting deviated less than 1 % compared with refined solutions, see B.1.
The CFD simulations were made on a structured hexahedral cell mesh in a cylindrical domain with diameter 10D and height 5h.

Data analysis and presentation
The results of the simulations will be analysed based on five key signals: the three motion responses of the buoy centre of gravity in surge (η x ), heave (η z ) and pitch (θ), as well as the dynamic tension force magnitude ( = − τ T T 0 ) in cable 1 and 2 respectively. T 0 is here the pretension magnitude of the mooring lines. We limit our analysis to view forces in cable 1 at the buoy connection point B, and forces in cable 2 at the fair-lead point C. As the cables themselves are relatively light-weight, stiff and under strong pretension, the force in each cable is relatively homogeneous. For the buoy motion properties (point F) we present results in terms of the first order Fourier amplitude based on the last 5 cycles of the 10 cycle test period to avoid start-up transients, which coincides with the motion amplitude at the excitation frequency. These amplitudes are also labelled with^. The phase shift of the motion at point B is also introduced as where the scalar product and the norm operators operate on the position time histories at point B (η (B) ) and point C (η (C) ) respectively. Both are sampled at a constant frequency. Further, the dynamic tension range τ is used to quantify the mooring forces. The definition used in this paper is

Influence of buoy geometry
In this section we compare the responses of the three SB geometries. The results were computed using the quadrature Morison drag formulation only.

Decay tests
Decay tests were made as a basis for understanding the system and as an initial sensitivity check on system parameters. We assign an initial buoy velocity of s. The high-frequency peak amplitude decreases with increasing slenderness, and the D 2 buoy spectrum (Fig. 6d) is basically single-peaked. Fig. 6 also shows a sensitivity study on the added mass coefficients C M1 and C M2 , where the coefficients were modified by ± 25 %. We notice a sharp decrease in the high-frequency response peak (around = T s 2.2 p ) due to changes in C M2 coefficient, which is consistent with the changes due to buoy geometry. The longer and thinner D 2 buoy has a smaller hydrodynamic mass in the vertical direction. We can thereby ascertain that the high-frequency peak is connected to the vertical offset of point B, while the low-frequency peak is the natural frequency of buoy pitch.

Dynamic load response
The three SB designs D 0 , D 1 and D 2 were analysed using a cyclic excitation of the fairlead. The motion amplitude is set to = a 0.5 m. Fig. 7 shows the envelope amplitude responses for all three buoy geometries. We make the following observations: • The surge motion is strongly amplified ( > η â 2.5 x ) for all geometries. Displacement amplitude decreases with increasing buoy slenderness. The peak motion response also occurs at increasing periods as the slenderness increases. See Fig. 7a and b.
• The D 2 cylinder outperforms D 0 and D 1 in that the dynamic loads are smaller for T p ≥ 2 s, see Fig. 7f and e. Also, by increasing the minimum tension the D 2 buoy avoids slack-snap response in the cables around the horizontal resonance of the buoy (for this loading amplitude). Several load cases were slacking in cable 1 for the D 0 and D 1 buoys.
• The high-frequency increase in tension for T p < 2 s is to be expected and is an artefact of keeping the load amplitude a constant. For high frequencies, the accelerations of the system grow very large, which is an unrealistic loading scenario for offshore moorings. Therefore, we leave these cases outside our analysis.
• As can be expected from the geometric setup of the mooring system, the fairlead force (at point C) envelope behaves similar to the pitch response and surge response of the buoy, while the response at point B is more related to the heave motion of the buoy. We also notice the similar behaviour of τ C and φ x .
• The dynamic tensions at point B and point C are overall in the same range, however we highlight that the tension range at point B can exceed that of the fairlead at higher frequencies.

Coefficient sensitivity
The same coefficient sensitivity analysis as in the decay tests was made for the full envelope of dynamic tests. Fig. 9 show results of τ , B τ , C η x and θ for variations of ± 25 % in the radial hydrodynamic coefficients for added mass (C M1 ) and drag (C D1 ) respectively. The plots of Fig. 9 have been normalised by the maximum value in the envelope for each parameter. Please note that the vertical coefficients (C M2 and C D2 ) were also analysed but gave much smaller dependency of coefficients. This was to be expected as the surge and pitch modes of motion are significantly more amplified in this load scenario (see Fig. 7). Some key aspects of the coefficient sensitivity are pointed out: • Fig. 9b and d reveal that drag coefficient variations affect the amplitude of the buoy response. The sensitivity to drag coefficient is relatively low (25% coefficient variation gives approximately 10% motion difference) and the differences are concentrated around the peak responses.
• Fig. 9 shows that the buoy responses (surge and pitch) are overall more affected by variations in hydrodynamic coefficients than the tension forces (τ B , τ C ).

Discussion
The results highlight several important aspects of mooring dynamics with submerged buoys (SBs). Previous studies have mainly focused on choosing a suitable buoy size [18], the number of buoys [11], or position [17], choices that affect both the pre-tension of the mooring system and the quasi-static force response curve. The three SBs analysed in this paper are all of the same net buoyancy, which makes them equivalent in a quasi-static mooring design. However, a comparison of the envelope response (Fig. 7) and time traces of tension (Fig. 8) shows that the shape of the buoy has a significant impact on the resulting mooring dynamics. For this load case the D 2 buoy effectively avoids the slack-snap response in the fairlead cable, see Fig. 8d. Slack-snap response is by contrast evident in the D 0 case, and is very close to happen also in the D 1 geometry at = T 4 p s. Avoiding regularly occurring slacksnap behaviour is a critical criteria of the design [3,39]. Another important aspect is the number of load cycles that the moorings are subjected to. In Fig. 8 we stress that the tension force τ B is bi-harmonic while the fairlead force is approximately of single harmonic shape. Rainflow counting is still the primary tool to determine partial damage in fatigue analysis [34], and the expected lifetime of the lower line will clearly be affected by the change in load history due to buoy geometry. Consequently, we recommend that field-tests or laboratory experiments of moorings using this type of system should measure also the force in the lower line, as well as the motion of the buoy.
The results consistently point to that the D 2 buoy outperforms the D 0 and D 1 buoys, based on that it generates the smallest loads per excitation. However, in real applications we must also consider the installation, maintenance and operation of the moored structure, and the restraining function of the mooring system must be evaluated in that context. For a moored body such as a wave energy converter subjected to environmental loads from wind, waves and currents, there is a strong coupling between mooring loads and mooring stiffness. Therefore, a complete judgement of performance of a particular buoy geometry cannot be made from a mooring analysis in isolation but requires coupled simulations with the moored object. Understanding the mechanisms of buoy dynamics enables a design parameter that the authors believe can provide interesting features to mooring systems of wave energy converters. Submerged buoys have been predominantly considered as spherical objects used for deep waters [11], or have been modelled as point-bodies without the rotational effects [3,34]. By analysing the full dynamics of the SB-mooring system, one can tune not only the mooring line angle and the net buoyancy but also the ballasting and shape of a buoy to best suit the purpose of the moored structure across the full frequency range.

Modelling fidelity
In this section we compare the four methods for buoy dynamics used in this paper: (i) CFD simulations (cfd); (ii) the quadrature Morison drag model (quad); (iii) the independent Morison drag model (ind); and (iv) the translating Morison body model (xyz). The CFD method is considered the reference solution in the absence of experimental data. As the trend due to SB slenderness is consistent in Section 5, we consider only the D 0 and D 2 buoys. We also limit the CFD analysis to T p ∈ [2, 10].

Force response
We begin with analysing the total force acting on the buoy in the different methods. Fig. 10 shows the total surge and pitch force on the buoy, i.e. the right hand side of (4).
The overall impression is that the translation only (xyz) Morison model fails to represent the physics of the CFD simulations, while both the independent (ind) and the quadrature (quad) models provide similar and adequate approximations at most frequencies. However both Morison models (quad and ind) predict an increase in pitch moment around period = T 4 p s, of which there is no sign in the CFD simulations. The buoy simulated by CFD is overall experiencing smaller forces than the Morison methods, and around the peak surge force frequency, there is a tendency that the quad method is better describing the CFD results than the independent method. The predicted moments on the buoy are in good agreement for all methods except the xyz model where it is left outside the modelling. In this load scenario the ind and quad methods are relatively close to each other, however, in a more general simulation the two models may differ significantly. To highlight the potential effect of the modelling fidelity of the drag force, we use the baseline buoy D 0 as a numerical example and provide the resulting drag force for a range of linear (v î ) and rotational velocities (ω ĵ ). Fig. 11 a shows the drag force in the î direction for different pitch rotation speeds. Clearly the independent approach is a valid approximation for the drag force as the total force is only marginally affected by the rotation. However, in Fig. 11b the increasing mean velocity v î has a strong influence, resulting in a drag moment difference of up to an order of magnitude between the two methods. For SBs in strong currents or those excited by large amplitude motions, the impact of modelling fidelity may therefore be significant. Fig. 12 shows the envelope amplitude responses of the four methods for both SBs. First, we see that the translating method results deviate significantly from the other three methods, which was to be expected from the simplistic modelling assumptions. The other three methods agree moderately well on the response prediction for the baseline buoy (D 0 ). Mooring forces shown in Fig. 12e and f are in surprisingly good agreement given that the CFD predicts a lower surge and pitch motion amplitude, see Fig. 12c and a. There is however a difference between the CFD and Morison approaches in that the CFD tension at point B has only a single peak around = T 4 p s while the Morison results show a double peak. It is likely that this second peak in the Morison models originates from an interaction with the pitch motion, as it occurs at the same frequency as the pitch moment deviates from the cfd results in Fig. 10b. The results from the two Morison drag formulations are naturally in much closer agreement with each other than with the CFD results. However, around the peak of the surge motion of the buoy, the    quadrature model is somewhat closer to the CFD results than the independent method results, which is supported by the differences observed in Fig. 10.

Dynamic load response
The same trends between the methods discussed above for D 0 are seen also in the D 2 responses of Fig. 12, but with larger differences due to method fidelity, which indeed is to be expected. A significant difference in surge-and pitch-response can now be seen in Fig. 12c and a. Fig. 12 shows differences between the dynamic models for SB motion, but at the same time the Morison approaches have a sensitivity to the hydrodynamic coefficient values (see Fig. 9). Therefore it is of interest to separate the difference due to model fidelity from the difference due to a sub-optimal choice of hydrodynamic coefficients. The quadrature and the independent Morison drag models are therefore fitted to the time series of the CFD results to produce a best-fit set of hydrodynamic coefficients for each case. This is achieved by constructing an inverse model for the buoy dynamics, where the buoy state variables and the driving force from the moorings are known quantities. First we isolate drag and added mass effects from Eq. (4), and group remaining (known) variables as F*:

Coefficient analysis
We focus on the surge, heave and pitch equations and seek the best-fit set of coefficients → C from the linear system where A and Ã are the coefficient matrices of the quadrature and independent models respectively. We have here applied a small and constant tangential shear force coefficient , as well as still water conditions for brevity. We apply Eq. (32) on each of the n output times of the CFD results, forming the extended system of 3n equations The best fit coefficients are then computed by a least square error analysis, as with corresponding expression for the independent model using Ã . From the coefficient sensitivity analysis (see Fig. 9), we know that the drag coefficient has negligible influence in T p < 4 s. We therefore limit the analysis to the T p ∈ [4, 10] s window of the CFD results, where both coefficients have sufficient stiffness.
The resulting best-fit coefficients are presented in Table 3 with mean value and standard deviation. Both DNV added mass coefficients for D 0 match the CFD predictions very well, as does the C M2 value for the D 2 buoy. However, the CFD-predicted C M1 value for D 2 was lower than expected, but still matches within 15 %. All the added mass coefficients have a low standard deviation, verifying the low frequency dependence expected for deeply submerged buoys [11]. There are however considerable differences in drag coefficient prediction. Both the variation across the periods (σ) value, and the mean value differ significantly. Especially for the vertical drag (C D2 ), we see a very large coefficient predicted by the CFD, with a significant spread over the periods. Table 4 illustrates the effects of the best-fit coefficient analysis at the surge peak of each buoy respectively. We compare the phase-compensated relative difference ε, where Ψ is the phase difference from the cross correlation. The phase difference between the Morison models and the CFD data is small overall for the DNV coefficients (maximum 21 deg. for τ B ), but we also see a consistent reduction of Ψ (to maximum 6 deg. for θ) when applying the best-fit coefficients from Table 3. In particular, we highlight that also in cases where the relative difference ε is hardly affected by Table 3 Resulting mean coefficients including standard deviation. Based on fitting the quadrature and independent Morison models to the CFD results of T p ∈ [4, 10] s using Eq. (36). x denotes mean value and σ(x) is the standard deviation.  Table 4 Estimated difference between CFD data and Morison model simulations at the surge peak period of each buoy: = T 5.5 p for D 0 ; and = T 7.0 p for D 2 . Differences between CFD data and the quadrature (quad) and independent (ind) models using DNV coefficients are presented together with differences for best-fit (*) mean value coefficients presented in Table 3. Values show relative difference ε(x), and cross-correlation phase lag Ψ (deg). the coefficient variation (see η x , D 2 for the quadrature model), the temporal improvement is still substantial. Further, in the case of the fair-lead force τ C , it has a direct impact on the mooring induced damping [40]. Table 4 also quantifies the impression from Figs. 12 and 10 that the quadrature drag method performs better than the independent method for the D 2 case. The improvement of the average ε in Table 4 is 0.07 (25 %) when the two models are compared using DNV coefficients. However, the differences between the models are small for most cases, and as expected, both models arrive at almost the same accuracy when the best-fit coefficients are used. The inherit difference between the methods is instead absorbed in the different coefficient values of Table 3. The best-fit coefficient results still show a residual difference from the CFD results ( = ε 0.09 for D 0 and = ε 0.15 for D 2 ) which is then attributed to the difference in modelling fidelity.

Flow characteristics
Snap-shots of the vortical structures generated by the D 0 and D 2 cylinder movements are presented in Fig. 13 for = T 6 p s. As expected the leading corner of the buoy is the location of interest which dominates the flow. A vortex is built up over the surge amplitude and is released from the body during the large pitch accelerations when the buoy turns back into its own wake. The results show how the pulling of the upper cable induces a stronger response in the right-going surge than the buoyancy driven return leg (leftward in Fig. 13). It is likely that the damping effect of these vortices is not fully included in the drag coefficients of DNV [33] as they are an effect of combined pitch and surge velocities.
For the D 2 buoy the difference in pitch response and surge phase is clearly seen. It remains more vertical than the D 0 buoy, which induces The D 2 buoy is also subjected to vortex induced motion in sway, which breaks the xz symmetry of the problem. This effect was not modelled in the Morison methods and is another clear example of physical effects that lead to differences between the modelling methods compared.

Discussion
We begin by stressing the importance of the rotational motion of the buoy in this configuration. The translation only (xyz) model is by the results of Figs. 10 and 12 proven inadequate for this type of buoys.
The fact that the peak response frequencies agree rather well in the comparison between the different methods (see Fig. 12), points to that the added mass coefficients used in the Morison methods was an adequate approximation. This is further strengthened by the coefficient analysis of the CFD results in Table 3, where in particular the DNV D 0 added mass was closely reproduced in the model fitting.
The CFD method is of significantly higher fidelity and thereby captures many more case specific fluid effects than the Morison methods, which is evident from the complex flow shown in Fig. 13. One such flow characteristic is that the buoy is moving in its own wake, as we model the case without current. These wake effects can be a contributing factor to the method differences, and may in part explain the large standard deviation (frequency dependence) of the estimated drag coefficients in Table 3.
From the study example in Fig. 11 we expect to see large differences in pitch and surge between the independent and quadrature methods. This is demonstrated in the D 2 buoy response of Fig. 12c where the independent drag method predicts significantly larger pitch amplitudes at the peak frequency compared to the quadrature drag method, which is closer to the CFD results. The trend is consistently seen in Fig. 12a for the surge motion and in the surge forces in Fig. 10. At the surge peak period ( = T 7 p s) of D 2 , the quadrature method was 25 % closer overall to the CFD results, see Table 4. The fact that there is a 35 % reduction in overall difference to CFD results when the best-fit mean coefficients from Table 3 are used points to that the coefficients are important, but that there is a significant portion of the difference inherit in the Morison approximations compared with CFD. A possible extension of the standard Morison equation, with e.g. an additional linear damping term, could further increase its performance compared with CFD. However, this has not been analysed in the present contribution. We stress that although the best-fit coefficients were chosen from mean values, the phase difference match to CFD was excellent, suggesting that the overall damping properties of the system were very well captured, and the coefficients were close to the optimal choice. The phase difference is of particular importance for the fair-lead force (τ C ) as it directly affects the mooring-induced damping of a moored structure [40].
Finally, a note on computational effort. The CFD simulations are obviously orders of magnitude more demanding than the Morison approaches (at least in the order of 1000 times) and is prohibitively expensive for many engineering applications. As typical mooring simulations are made with several buoys and lines on a sea-state timescale, the Morison approaches are the only feasible choice. However, we do recommend that some simulations are validated with a CFD response to ensure that no significant physical effects are badly modelled by the Morison approximation. The three Morison methods studied in this paper require a similar computational effort. Therefore there is seemingly no reason to use the translating or independent Morison models, as the quadrature model has higher fidelity and generality.

Conclusion
We have presented an analysis of the effects of buoy geometry and modelling fidelity for a mooring system with two polyester cables and a submerged buoy (SB). The effect of buoy geometry was studied by comparing three submerged cylinders of the same volume and mass but with varying slenderness, showing large differences in dynamic response of the mooring system. Previous studies have focused on different sizes, placement, and types of buoys, but in this paper we highlight the effect of the shape alone as a powerful design tool to change the dynamic mooring response. Overall D 2 , a more slender cylinder ( = h D / 3.9) out-performed the two shorter and wider cylinders (D 0 with = h D / 1.2 and D 1 with = h D / 1.8) giving both significantly smaller mooring force amplitudes as well as smaller displacements of the buoy. It also avoided slack-snap occurrences in the upper cable connected to the fairlead, which was seen for the other geometries. The lower line in the mooring leg is strongly affected by the buoy motion, often resulting in bi-harmonic load response when the fair-lead tension shows single harmonic loading, as was also mentioned in [11]. We therefore recommend that field-tests and experiments of moorings with an intermediate buoy should measure the force response in all lines, as well as tracking the motion of the buoy in six degrees of freedom.
We further compared four numerical methods of computing the dynamics of the SB: (i) viscous CFD simulations; (ii) the Morison equation with drag force using numerical quadrature; (iii) the Morison equation with independent treatment of the drag force in translation and rotation; and (iv) the Morison equation on a vertical cylinder in translation only. Our contribution compares four methods that encompass the full hierarchy of SB model fidelity, from the lowest (translation only) to the seldom used quadrature approach and further to the first (to the best of our knowledge) coupled mooring-CFD analysis of a mooring buoy. The paper shows that CFD as a high-fidelity simulation tool can be used to assess the performance of parametrised methods (such as the Morison equation in this case), and how improved hydrodynamic coefficients can be estimated from the CFD results. The translation only model is found inadequate for modelling mooring buoys of this type. The other two Morison methods are in relatively good agreement with the CFD results, with the quadrature model performing better than the independent model for slender SBs. Using hydrodynamic coefficients estimated from the CFD simulations reduced much of the phase difference of the fair-lead between the methods, which has a direct influence on the mooring-induced damping. The potential difference in drag moment between the methods can in cases with large currents or larger fair-lead excitation be significant. The computational costs of the Morison models are comparable, but the quadrature drag approach is more general and of higher fidelity. We therefore recommend to always use the quadrature-based Morison drag method for simulation of submerged mooring buoys.

Declaration of Competing Interest
The authors declare no conflict of interest Appendix A. Added mass for ballasted submerged buoys We will here present the derivation of the added mass forces on a slender cylinder-type rigid body, submerged in water. We define → = r δk b as the vector from the centre of gravity to the center of buoyancy of the cylinder. In the local tangential direction (k), the added mass force and torque due to linear acceleration and yaw acceleration respectively are independent of δ. We will therefore concentrate our analysis on the local normal direction of the cylinder.

A1. Accelerating flow contribution
The added mass forces and moments due to accelerating surrounding flow → a are added on the right hand side as: where C M1 and C M2 are the in-plane and out-of-plane added mass coefficients of a circle, V is the cylinder volume and h is the cylinder height.

A2. Inertial contribution
On the left hand side, the sectional force on a cylinder slice is integrated to give The mesh dependency of Moody was investigated by doubling the number of elements. As shown in Table B.5 the results are of a high-resolution and the mesh sensitivity is overall very small. The maximum difference in the first order amplitude response is 1.4%, seen for = T 5 p s in τ B . Hence, for the purpose of this paper we use the small resolution of = N 5, 1 = N 10 2 elements per cable 1 and cable 2 respectively, with an order = P 4 in each element. The exponential error decay of the high-order formulation, see [30], clearly provides sufficient resolution of the dynamics despite the seemingly low element count. The maximum difference in cases of no slack events were negligible for all three resolutions. We therefore focus the remainder of our analysis on the numerical treatment of the submerged buoy dynamics and how it affects the motions and tension response of the mooring leg.

B2. CFD Mesh sensitivity
Three different meshes were used to quantify the discretisation errors, the mesh data is presented in Table B.6. The results of a mesh independence study of the CFD simulations are shown in Table B.7. Table B.7 shows that the mesh sensitivity of the CFD analysis is relatively low, especially for pitch response θ which even for M1 is within 4 % of the M4 result. The results converge nicely for periods = T 6 p s and = T 10 p s in all quantities with a maximum difference of 3.5 % in τ B at = T 10 p s. At higher frequencies, larger differences are noted. In Fig. B.14 we present the time histories of the dynamic tension in cable 1 and 2 respectively for the = T 4 p s case. We notice an increasing amplitude response in τ B as the mesh refines. The difference is still relatively small compared to the cell count increase (M4 is 47 times larger than M1, see Table B.6). All meshes show the same trend so to maintain feasibility of the simulations we take note of the discrepancies in the mesh convergence but judge that the level of resolution is sufficient for use in our present effort to quantify the merits of parametrised methods using the CFD simulations as a reference.