Numerical Modelling of a Mussel Line System by Means of Lumped-Mass Approach

Abstract

This paper describes a numerical model to simulate the behavior of a mussel longline system, subjected to environmental loads such as waves and current. The mussel line system consists of an anchor, a mooring chain, a long backbone line, mussel collector lines and buoys. The lumped-mass open-source code MoorDyn is modified for the current application. Waves are modelled as a directional spectrum, and the current as a homogeneous velocity field with an exponential vertical distribution. A Coulomb model is implemented to model the horizontal friction between nodes and the seabed. Cylindrical buoys with three translational degrees-of-freedom are modelled by extending the simplified hydrodynamic model in use for line’s internal nodes with additional properties like cylinder height, diameter and mass. Clump weights are modelled in a similar way. For validation purposes, the results of the present software are compared with the commercially available lumped-mass based mooring dynamic software, OrcaFlex.

1. Introduction

Two different systems for mussel aquaculture are commonly described in the literature: Bottom culture, and suspension culture. In bottom culture the mussels are placed and grown directly on the seafloor without the need of any infrastructures [1]. The choice of location is a very sensitive aspect for bottom culture, due to the influence of predators, temperature, soil deposition and other organisms living on the sea bed. On the contrary, in suspension culture mussels are attached and grown on suspended ropes, away from the seafloor. The suspension of crop lines can be achieved with different strategies and different infrastructures. There are two main types of suspension culture: Raft culture and longline system [2]. In modern raft system, mussel crop lines are put within a structure made from a combination of wooden and steel frames to make a stable platform along with large floaters [2]. However, it can also be made from an old wooden boat or a catamaran boat to accommodate up to 1000 crop lines [1]. On the other hand, the longline suspension culture has a typical configuration consisting of a line that supports suspended ropes with mussels connected to buoys to keep them afloat. Moreover, the supporting line is moored to the seabed with anchors [3,4]. The longline system offers the possibility to have a substantial amount of mussel lines suspended from one long line. This makes it an advantageous technique, since it requires a minimum amount of auxiliary infrastructure [5]. Furthermore, the submerged mussel crops would be less influenced by surface waves, lessening the risk of system break or mussel detachment. A study to determine the influence of crop lines’ submergence in scallop longline culture has been done by performing a dynamic simulation with a FEM (Finite Element Method) based software [6,7]. The conclusion of the study was that the decrease in the stress on the mooring line could achieve 55% as the crop lines are submerged deeper [7]. Lien and Fredheim [8] performed dynamic analysis of mussel longline and mussel longtube system using FEM based software [9] originally used for offshore riser analysis. The purpose was to analyze the two floating systems with the conclusion that longtube yields in the reduction of dynamic loads compared to traditional longline system. The dynamics of a model mussel longline system was studied by Raman-Nair [10] by means of a lumped mass approach. The purpose of the paper was to come up with a method to perform a design optimization of shellfish aquaculture by means of numerical simulations. Hydrodynamic forces are modelled using Morison Equation and waves are modelled according to second-order Stokes theory. No validation work was done; however, the experiment is suggested in their conclusion to obtain the hydrodynamic coefficients used in the model.

A longline aquaculture system is essentially an underwater system of interconnected lines of various materials. The dynamics of such a system can be studied with different methods: Finite-difference [11,12], lumped-mass [13,14,15,16] and the inclusion of bending and torsion on cable segments by finite elements analysis (FEA) [17,18,19]. In lumped-mass approach, the governing equations are discretized in space, whereas it is discretized in space and time in finite-difference models [15]. FEA approach is superior to lumped-mass if the goal is to get an accurate prediction across different conditions [15]. In a dynamic mooring analysis study done by Hall et al. [20] with FEA-based approach, it was found that bending and torsion can be neglected in predicting tension of barge, spar and Tension Leg Platforms (TLP). Based on that finding, Hall [15,16] developed a simple lumped-mass based mooring dynamic solver to predict line tensions of the typical offshore mooring system. The code has been validated against the model test of a semisubmersible floating wind turbine with good fairlead tensions agreement [15,16].

In Section 3, research methodology is presented, which consists of 5 sub-sections. In Section 3.1, the open-source code MoorDyn and is briefly introduced. The equation of motion and assumptions made to model the environmental induced load are explained in Section 3.2. The different assumptions implemented in adapted MoorDyn and OrcaFlex are briefly described in Section 3.3. Test cases setup used for five simulation cases performed in adapted MoorDyn and OrcaFlex are described in Section 3.4. In Section 4, simulation results are presented and discussed. Section 5 will conclude the study done in this paper, along with suggestions for future work.

2. Field Site

In 2017 a research project to examine the feasibility of mussel aquaculture within the Belgian North Sea Wind Farms was established under the name of Edulis. Belgian wind farms are located in a restricted area where the space can be utilized for wind power generation and aquaculture. Two test setups are already deployed in the Belgian wind farms Belwind and C-Power (see Figure 1) to investigate both the forces on the system and the biological properties of the mussels.

A very important requirement to be complied with is that the position of the mussel line system, under all circumstances, will not interfere with wind farms activities. Careful design of the mussel growing system is required for deployment in a critical environment, such as the one found at offshore locations in the North Sea [21]. The semi-submerged system used in the test setups (Figure 2) consists of a long backbone line connecting several buoys. The backbone line supports a total of 145 m long v-shaped mussels collector lines. Properties of the mooring configuration for the test setups are summarized in

Table 1. Line properties.

Line Type [-]	Dry Mass per Length [kg/m]	Nominal Diameter [m]	Line Length [m]
Chain (Grade 3 steel)	10.910	0.022	108
Backbone (Movline Plus 8 strands)	2.1	0.068	57
Mussel sock (fully grown mussels)	21.8	0.15	145

,

Table 2. Buoy properties.

Buoy Type [-]	Outer Diameter [m]	Dry Mass [kg]	Length [m]	Volume [m3]	Quantity [-]
SPAR buoy	0.790	2500	8.865	4.345	2

and

Table 3. Anchor properties.

Anchor Type [-]	Dry Mass [kg]	Quantity [-]
Gravity	15000	2
Danforth	2500	2

. This system with submerged crop lines is expected to experience less load on the mooring line as it is submerged deeper [7]. Despite the protection from the surface waves, the submerged longline system is still in high risks of system failures (anchor loss, buoy immersion, line breakage), due to storms [22]. Hence, a mathematical model is needed to predict the dynamic behavior of a mussel line system under the influence of environmental loads. To this end, modifications are implemented to an open-source lumped-mass mooring dynamic solver. The aim of this paper is to show the results of newly adapted open-source code compared to commercial software. A comparison with the commercial mooring dynamic solver OrcaFlex [23] is performed to assess the performances of the newly adapted code with respect to wave-induced loads acting on a three degrees-of-freedom (all translations, zero rotations) buoy.

3. Methodology

The open-source code used as the starting point in the development of the mathematical model is MoorDyn (Fortran version) developed by Matthew Hall [15,16]. In this paper, comparison with a finite element based commercial software OrcaFlex will be presented.

3.1. General Description of Numerical Models: MoorDyn

The original version of MoorDyn is capable of modelling the chain, backbone and mussel collector lines. The segments of all these lines are modelled as homogenous cylinders, characterized by a certain diameter. Concerning chains, the nominal shackle diameter needs to be converted to an equivalent diameter when assuming a cylinder shape [24]. The original MoorDyn code is modified in order to be able to simulate a mussel longline in waves and current. This requires the introduction of the effects of environmental loads in the original code, as well as three degrees-of-freedom buoys, clump weights and the effects of seabed friction. The adapted code allows modelling the complete Edulis test setup as it is shown in Figure 3. This model assumed a flat sea bottom across the entire domain.

The code was originally developed with the goal of having a computationally efficient solver with the focus on predicting tensions on a typical mooring line system. This is achieved by omitting rotational degrees of freedom, neglecting bending and torsional stiffness. MoorDyn uses a lumped-mass approach to model the dynamics of underwater mooring lines and Morison Equation to calculate the hydrodynamic forces acting on such lines. For the purpose of modelling the mussel longline system, adaptations to MoorDyn are introduced to model environmental loads, buoys, clump weights and the effects of seabed friction. According to the lumped-mass approach, a line is represented as a series of contiguous segments with homogenous properties. Each line is evenly divided into N number of segments, corresponding to N+1 nodes (see Figure 3). The mass of each segment is equally transferred to both its extremity nodes. The same principle applies to the external forces acting on each segment. Forces are calculated at each node and at each time step [15]. The first and last node in a line is defined as connection nodes. There are three types of connection nodes:

• Fixed: Fixed type restricts the node to prevent any displacement in position.
• Vessel: This type only allows the node to move according to a prescribed motion provided as an input to the code. External forces do not influence the displacement of vessel type connection nodes.
• Connect: This type allows the node to freely move in any directions according to all the forces acting on it.

A comparison of MoorDyn with OrcaFlex against model test data of OC3-Hywind spar buoy was investigated, and good agreement between the two codes was found [25].

3.2. Mathematical Model

This sub-section will describe the assumptions used in the numerical model of. The focus will be put on the adaptations made to the original MoorDyn.

3.2.1. Equation of Motion

In MoorDyn, the equation of motion for each node can be expressed as follows [15]:

$$\left(m+a_i\right){\ddot{r}}_i\left(t\right)=F_i\left(r_i\left(t\right), {\dot{r}}_i\left(t\right), v_i\left(t\right), {\dot{v}}_i\left(t\right)\right),$$

(1)

where:

• $r_i\left(t\right)$ is the node’s position at instantaneous time t [m];
• ${\dot{r}}_i\left(t\right)$ is the node’s velocity at instantaneous time t [m/s];
• ${\ddot{r}}_i\left(t\right)$ is the node’s acceleration at instantaneous time t [m/s²];
• i represents the 3 degree of freedom in translation (i = 1, 2, 3);
• $v_i\left(t\right)$ is the fluid velocity at instantaneous time t [m/s];
• ${\dot{v}}_i\left(t\right)$ is the fluid acceleration at instantaneous time t [m/s²];
• $m$ is the mass of the node [kg];
• $a$ is the hydrodynamic added mass of the node [kg];
• $F$ is the total force acting on the node [N].

In this equation of motion, the force acting on the node is expressed as a function of time, outlining the explicit dependencies on the position and velocity of the node and on the instantaneous values of the fluid velocity and acceleration at the node’s location. The dependency of the force on the fluid instantaneous velocity ${\dot{v}}_i\left(t\right)$ and acceleration ${\ddot{v}}_i\left(t\right)$ is introduced in the current modification of MoorDyn. The force acting on a node can be further split into different components induced by different physical phenomena:

$$F_i\left(r_i\left(t\right), {\dot{r}}_i\left(t\right), v_i\left(t\right), {\dot{v}}_i\left(t\right)\right)=D_p+D_q+T+C+W+B.$$

(2)

The right-hand side of Equation (2) are the contributions of internal forces and external forces. Line numerical damping $C$ is a function of strain rate and damping coefficient. Tension $T$ is calculated as a function of segment strain and Young modulus of the line. Net weight $W$ is calculated based on the segment’s dry weight and buoyancy. These are the internal forces. As for the external forces, $D_p$ and $D_q$ are transverse and tangential hydrodynamic drag forces modelled using the Morison Equation. Vertical bottom contact $B$ is modelled by a spring-damper system. The forces are illustrated in Figure 4. For additional details on how the aforementioned forces are modelled, the reader is referred to Reference [15].

In three dimensions, Equation (1) can be expressed in a matrix form as:

$$\left[\begin{array}{ccc}\left(m+a\right)& 0& 0\\ 0& \left(m+a\right)& 0\\ 0& 0& \left(m+a\right)\end{array}\right]\left[\begin{array}{c}{\ddot{r}}_x\\ {\ddot{r}}_y\\ {\ddot{r}}_z\end{array}\right]=\left[\begin{array}{c}F\left(r_x, {\dot{r}}_x, t\right)\\ F\left(r_y, {\dot{r}}_y, t\right)\\ F\left(r_z, {\dot{r}}_z, t\right)\end{array}\right]$$

(3)

The entire equation of motion for a line of N segments is expressed in a matrix which consists of N+1 sub-matrixes like the one shown in Equation (3). The differential equations of motion are solved using a constant time step Runge-Kutta second order (RK2) integration scheme [26]. The implementation of this integration scheme is further explained in Appendix A. In OrcaFlex, bending and torsional elasticity are taken into account [24]. However, this can be omitted in the calculation, which was done in the simulations performed for this study.

3.2.2. Current

Current is implemented as a uniform velocity field over the whole simulation domain. Both in OrcaFlex and adapted MoorDyn, a vertical distribution of the current velocity in three dimensions can be implemented according to the power law:

$$v_{cx}= s_0{\left(\frac{d+z}{d+z_0}\right)}^{\alpha }·cos\theta ,$$

(4)

$$v_{cy}= s_0{\left(\frac{d+z}{d+z_0}\right)}^{\alpha }·sin\theta ,$$

(5)

$$v_{cz}= 0,$$

(6)

where:

• $s_0$ is the current magnitude at a reference depth (defined in the input) [m/s];
• d is the water depth (defined in the input) [m];
• $z_0$ is the reference depth (defined in the input) [m];
• $\alpha$ is the exponent to control the shape of vertical distribution [-];
• $v_{cx}$ is the current velocity component in x-direction [m/s];
• $v_{cy}$ is the current velocity component in y-direction [m/s];
• $v_{cz}$ is the current velocity component in z-direction [m/s].

3.2.3. Waves

The wave climate is defined by a directional wave spectrum. For each spectral component, Airy linear wave theory [27] is used to calculate the free surface elevation and the orbital velocities and accelerations induced by the spectral component. The total elevation and kinematic quantities, due to the wave spectrum can be calculated by a superposition of the contributions of each spectral component [28].

$$\eta ={\displaystyle \sum }_{j=1}^M{\displaystyle \sum }_{l=1}^P{a}_{j,l }cos\left(k_{j}xcos{\theta }_l+ k_{j}ysin{\theta }_l- {\omega }_{j}t+ {\alpha }_{j,l}\right),$$

(7)

$$v_{wx}={\displaystyle \sum }_{j=1}^M{\displaystyle \sum }_{l=1}^P{a}_{j,l}{\omega }_j\frac{cosh\left(k_j\left(d+z^{\prime }\right)\right)}{sinh\left(k_{j}d\right)}\mathit{cos}\left(k_{j}xcos{\theta }_l+ k_{j}ysin{\theta }_l- {\omega }_{j}t+ {\alpha }_{j,l}\right) cos{\theta }_l,$$

(8)

$$v_{wy}={\displaystyle \sum }_{j=1}^M{\displaystyle \sum }_{l=1}^P{a}_{j,l}{\omega }_j\frac{cosh\left(k_j\left(d+z^{\prime }\right)\right)}{sinh\left(k_{j}d\right)}cos\left(k_{j}xcos{\theta }_l+ k_{j}ysin{\theta }_l- {\omega }_{j}t+ {\alpha }_{j,l}\right) sin{\theta }_l,$$

(9)

$$v_{wz}={\displaystyle \sum }_{j=1}^M{\displaystyle \sum }_{l=1}^P{a}_{j,l}{\omega }_j\frac{sinh\left(k_j\left(d+z^{\prime }\right)\right)}{sinh\left(k_{j}d\right)}sin\left(k_{j}xcos{\theta }_l+ k_{j}ysin{\theta }_l- {\omega }_{j}t+ {\alpha }_{j,l}\right),$$

(10)

$${\dot{v}}_{wx}={\displaystyle \sum }_{j=1}^M{\displaystyle \sum }_{l=1}^P{a}_{j,l}{\omega }_j{}^2\frac{cosh\left(k_j\left(d+z^{\prime }\right)\right)}{sinh\left(k_{j}d\right)}sin\left(k_{j}xcos{\theta }_l+ k_{j}ysin{\theta }_l- {\omega }_{j}t+ {\alpha }_{j,l}\right) cos{\theta }_l,$$

(11)

$${\dot{v}}_{wy}={\displaystyle \sum }_{j=1}^M{\displaystyle \sum }_{l=1}^P{a}_{j,l}{\omega }_j{}^2\frac{cosh\left(k_j\left(d+z^{\prime }\right)\right)}{sinh\left(k_{j}d\right)}sin\left(k_{j}xcos{\theta }_l+ k_{j}ysin{\theta }_l- {\omega }_{j}t+ {\alpha }_{j,l}\right) sin{\theta }_l,$$

(12)

$${\dot{v}}_{wz}={\displaystyle \sum }_{j=1}^M{\displaystyle \sum }_{l=1}^P-a_{j,l}{\omega }_j{}^2\frac{cosh\left(k_j\left(d+z^{\prime }\right)\right)}{sinh\left(k_{j}d\right)}cos\left(k_{j}xcos{\theta }_l+ k_{j}ysin{\theta }_l- {\omega }_{j}t+ {\alpha }_{j,l}\right),$$

(13)

where:

• $\eta$ is the free surface elevation [m];
• $M$ is the number of wave frequency [-];
• $P$ is the number of wave direction [-];
• $a$ is the wave amplitude [m];
• $\omega$ is the wave angular frequency [rad/s];
• $t$ is the instantaneous time [s];
• $k$ is the wavenumber [rad/m];
• $d$ is the water depth [m];
• $\theta$ is the wave direction [rad];
• $\alpha$ is the phase angle [rad];
• $v_{wx}$ is the wave-induced velocity component in x-direction [m/s];
• $v_{wy}$ is the wave-induced velocity component in y-direction [m/s];
• $v_{wz}$ is the wave-induced velocity component in z-direction [m/s];
• ${\dot{v}}_{wx}$ is the wave-induced acceleration component x-direction [m/s²];
• ${\dot{v}}_{wy}$ is the wave-induced acceleration component y-direction [m/];
• ${\dot{v}}_{wz}$ is the wave-induced acceleration component z-direction [m/s²].

Since the Airy linear theory formulation does not extend above free surface, the kinematic properties of each spectral component are extended above the calm water surface by using Wheeler’s kinematic stretching approach [29]. As for the calculation of dispersion relation, Fenton approximation is used [30]. This method yields exact results in shallow water and deep water and produces errors of no more than 0.05% throughout all wavelengths in intermediate water depth [30].

3.2.4. Wave-Current Interaction

The interaction between waves and current is a very complex hydrodynamic phenomenon, and the literature which covers it is extensive. In general, the current has an effect on both the kinematic and dynamic properties of waves. The effect of current on wave kinematics is essentially a change of the dispersion relation, due to the Doppler effect [31]. In the case of a steady uniform current, the Doppler effect can be trivially calculated based on the uniform current speed. For vertically varying currents, more complex relations need to be used. An approximate solution valid for current profiles of arbitrary shapes up to the second order was given in Reference [32]. On the other hand, the effects of current on wave dynamics account for changes in the wave amplitude, due to wave action conservation [33]. An expression to calculate the modifications to a wave spectrum, due to the presence of current, was first provided in Reference [34]. Through Doppler effect and considerations deriving from wave action conservation, the effect of current on waves can be included in a phase averaged representation of the wave field.

The problem of wave current interaction is further complicated in shallow water, where both waves and currents interact with the bottom as well. This very complex problem is still the object of ongoing research, mainly focused on the use of phase resolving wave models. Recently, the problem was tackled in Reference [35] through a non-linear vertical 2D model which solves the wave-current-bottom interaction taking also into account the effects of vorticity induced by vertically non-uniform currents. In this case, the model only solves the interaction of waves on a vertical strip under the effects of collinear currents. A different approach was followed in Reference [36], where a turbulent jet current interacting with frontal waves was investigated by means of a 2DH non-linear model based on shallow water equations. A 2DH non-linear shallow water model was also used in Reference [37] to investigate the effects of turbulence and seabed friction on macrovortices generated at a river mouth.

While phase resolving models allow investigation into detail the hydrodynamics of wave-current interaction, such models are still rather heavy from a computational point of view. In the present study, a numerical model to investigate the dynamics of a moored system in waves and current is presented. For the moment, wave-current interaction is not considered in the code, and the orbital velocities induced by the wave components of the input wave spectrum are superimposed to the current velocity, which is defined as a vertical profile. The present code could be easily improved in the future by implementing a simplified wave-current interaction based on the Doppler effect and on a modification of the spectral density according to wave action conservation. On the other hand, major modifications to the code would need to be implemented in order to use the results of more computationally expensive codes like the ones outlined above.

3.2.5. Hydrodynamic Forces

The wave-induced modifications to the free surface act on the buoyancy of the line segments and of the connection nodes. On the other hand, wave and current induced fluid velocities, and wave-induced fluid accelerations all contribute to the total hydrodynamic forces acting on the system’s nodes. Wave-current interaction is not modelled, hence, total fluid velocity component $v_i$ is the summation of wave orbital velocity component $v_{wi}$ and current component $v_{ci}$. Furthermore, current acceleration component ${\dot{v}}_{ci}$ is not modelled $({\dot{v}}_{ci}=0)$, which makes the fluid acceleration component ${\dot{v}}_{\mathrm{i}}$ is only due to the wave-induced acceleration component ${\dot{v}}_{wi}$.

$$v_i =v_{wi}+ v_{ci} ,$$

(14)

$${\dot{v}}_i = {\dot{v}}_{wi}+ {\dot{v}}_{ci} ,$$

(15)

$${\dot{v}}_{ci} = 0.$$

(16)

In order to calculate the hydrodynamic forces, relative flow velocity is calculated for each node by subtracting the fluid from the node’s velocities $\dot{r}$. The relative fluid velocities are decomposed into a transverse and tangential direction. Drag forces on the aforementioned directions are then calculated through Morison Equation. Fluid acceleration ${\dot{v}}_i$ is used to calculate the inertial force, which is a contribution of Froude-Krylov force and hydrodynamic added mass. The force per unit length, according to the Morison Equation is expressed in the following equations [38]:

$$F_{hydrodinamic}\left(t\right)= F_{Inertia}+F_{Drag},$$

(17)

$$F_{Inertia}=\frac{\pi }{4}\rho D^2 {\dot{v}}_i\left(t\right)+\frac{\pi }{4}\rho C_a{D}^2 \left({\dot{v}}_i\left(t\right)-\ddot{r}\left(t\right)\right),$$

(18)

$$F_{Drag}= \frac{1}{2}\rho C_{D}D \left(v_i\left(t\right)-r\left(t\right)\right)\left|v_i\left(t\right)-\dot{r}\left(t\right)\right|,$$

(19)

where:

• $C_D$ is the drag coefficient [-];
• $C_a$ is the inertia coefficient [-];
• $D$ is the cylinder equivalent diameter [m];
• $\rho$ is the density of the water [kg/m].

3.2.6. Buoys and Clump Weights

The buoy is modelled as a cylinder and an extension to a connection node with additional properties of cylinder height, diameter and mass. The submerged volume of the buoy (orange colour in Figure 5) is calculated based on the free surface calculated at the exact position of the node. Wet volume fraction is calculated as the ratio of the submerged volume and the total volume of the cylinder. Buoyancy is then calculated using the total volume and the newly calculated wet volume fraction. The drag force is calculated through Morison Equation again by taking into account the total projected area, multiplied by the wet volume fraction calculated before.

3.2.7. Seabed Friction

When the absolute horizontal velocity $\left|V_a\right|$ of a node is zero, the static friction force is applied. The direction of the friction force is opposite to the direction of the node’s resultant force in the horizontal axis $\overrightarrow{F_H}$. The static friction coefficient ${\mu }_s$ is applied in this case. On the other hand, the kinetic friction force is applied when the node starts moving. In this case, the friction force acts in the opposite direction of the node’s velocity, with a magnitude equal to the normal force multiplied by the kinetic friction coefficient ${\mu }_k$.

$$F_{friction}=\{\begin{array}{c}-min\left({\mu }_s F_N, \left|\overrightarrow{F_H}\right|\right)\frac{\overrightarrow{F_H}}{\left|\stackrel{\rightharpoonup }{F_H}\right|}, | \left|V_a\right|=0\\ {\mu }_k{F}_N\frac{V_a}{\left|V_a\right|}, \left|V_a\right|\ne 0\end{array}.$$

(20)

3.3. Adapted MoorDyn and OrcaFlex Comparison

This sub-section will briefly describe the similarities and differences of the assumptions implemented in OrcaFlex in comparison to the ones in the adapted MoorDyn.

3.3.1. Seabed Friction

Sea bed friction in OrcaFlex is implemented differently. A modified version of the Coulomb friction model is implemented. The modification to the standard Couloumb model is the implementation of a ramping zone to create a transition from the node being static into moving to a certain deflection criterion. The friction force is increasing linearly throughout this ramping zone [24].

3.3.2. Line Theory

OrcaFlex modelled a line in a similar way to MoorDyn, which is made up of numerous segments with a node connecting contiguous segments. The mass, weight, buoyancy, drag and other properties of half-segment are lumped into the neighboring nodes. As for the axial and torsional properties, it is modelled in the segment as a massless spring following the FEA based approach [24]. The tension at the segment’s center is calculated as follows [24]:

$$T_e= T_w+\left(p_o{a}_o- p_i{a}_i\right),$$

(21)

$$T_w= EA\epsilon -2\upsilon \left(p_o{a}_o- p_i{a}_i\right)+ k_{tt}\frac{\tau }{l_0}+EAc\frac{dl}{dt}\frac{1}{l_0},$$

(22)

where:

• $T_e$ is the effective tension;
• $T_w$ is the wall tension;
• $p_o$, $p_i$ are the internal and external pressure;
• $a_o$, $a_i$ are the internal and external cross-section stress area;
• $EA$ is the axial stiffness of the line, which is Young’s modulus multiple by cross-section area;
• $l$ is the instantaneous length of the segment;
• $\epsilon$ is the strain;
• $l_0$ is the segment’s unstretched length;
• $\upsilon$ is Poisson ration;
• $k_{tt}$ is tension and torque coupling;
• $\tau$ is the segment twist angle;
• $c$ is damping coefficient;
• $dl/dt$ is the rate of increase of length.

3.3.3. Waves

Various wave theories are available to choose in OrcaFlex: Airy, Dean, Stokes’ 5th, Cnoidal [24]. Waves are modelled the same way as in the adapted MoorDyn, which was explained in Section 3.2.3. Waves in MoorDyn are generated through an external program, based on JONSWAP, which randomizes the phase for each wave component. OrcaFlex offers the possibility to have wave input based on JONSWAP spectrum as well. However, in order to have the exact same component, all the wave components need to be defined manually for the test cases used in this paper.

3.3.4. Integration Scheme

In OrcaFlex, the user might choose two different integration scheme: An explicit scheme, and an implicit scheme. The explicit scheme is essentially a semi-implicit Euler scheme, which does not have the same order of accuracy as the one implemented in MoorDyn (see Appendix A). As for the implicit scheme, it is implementing the generalized-$\propto$ integration scheme [24]. For the purpose of adapted MoorDyn and OrcaFlex comparison, the explicit scheme is chosen in OrcaFlex.

3.4. Test Cases Setup

A numerical model which consists of an anchor connected to a three degrees-of-freedom buoy through a chain is built in the adapted version of MoorDyn and in OrcaFlex. In both codes, the anchor is modelled as a fixed point, which restricts any movements of the node in all degrees of freedom. Simulations are performed in both codes to compare the forces at the anchor and buoy’s position. The properties of the model, summarized in

Table 4. Line properties.

Equivalent Diameter [m]	Dry Mass per Length [kg/m]	Axial Stiffness [N]	Chain Length [m]	Can [-]	Cat [-]	Cdn [-]	Cdt [-]
0.042	10.910	48884000	50	1.0	0.50	1.4	0.2

and

Table 5. Buoy properties.

Equivalent Diameter [m]	Dry Mass [kg]	Height [m]	Volume [m3]	Can [-]	Cat [-]	Cdn [-]	Cdt [-]
0.790	1200	8.865	4.345	0.94	0.50	0.81	0.40

, are kept the same in both codes.

In both codes, the 50 metres chain is discretized as an array of 12 contiguous segments of homogeneous cylinder. The inertia coefficients of the chain are taken from recommendations provided by Bureau Veritas [39] in BV-NR493. As for the drag coefficients of the chain, suggested values were taken from Det Norske Veritas Germanische Lyold (DNVGL) [40]. The suggested values of drag coefficients are adjusted according to the chain equivalent diameter used in the numerical calculations. However, since these coefficients are kept the same in both OrcaFlex and MoorDyn, the chosen values are not strictly relevant to the results of this comparative study. The same principle applies to the Morison coefficients used for the buoy.

The seabed is modelled as flat bottom/constant bathymetry with a depth of 30 m. The initial position of the system is summarized in

Table 6. The initial position of the system.

Object Type	X [m]	Y [m]	Z [m]
Anchor	0.00	0.00	−30.00
Buoy	30.00	0.00	−4.43

.

OrcaFlex gives the option to select an explicit or implicit time integration scheme. As Moordyn uses an explicit scheme for time integration of the equations of motion, this is chosen in OrcaFlex as well. Moreover, compression forces are not modelled in both codes.

3.5. Environmental Loads

There are five simulation cases performed in the adapted MoorDyn and OrcaFlex. The implementation of the regular wave, seabed friction, current and unidirectional spectrum are investigated in these test cases. The environmental input for the simulation cases is described further in Section 4.

4. Results and Discussions

4.1. Simulation Case 1: Regular Wave

The system is subjected to a regular wave amplitude of 5 m and period of 8.33 s (wavelength = 103 m) propagating in the direction of the positive x axis.

Figure 6 shows the time series of buoy’s position calculated in adapted MoorDyn and OrcaFlex. Good agreement can be found for all direction (x, y, z). Due to the limitation of the chain length, no wave drift apparent as the buoy cannot move further. Figure 7 shows the comparison of calculated tension at the anchor for the period of t = 100 s–114 s. High frequency oscillations are found in OrcaFlex results with overshoot values that are 10 times higher than the peak values predicted by MoorDyn. These values occur when the last of the chain segments which lay on the seabed (furthest from anchor) suddenly changes its state from slack to taut. This change of state eventually leads to the segment being lifted up from the seafloor. Further investigation shall focus on seabed contact model and line numerical damping. The seabed contact model influences the total forces acting on that specific segment as the seabed provides reaction force to the chain segment. When sea bed stiffness provides high reaction force over a very short period of time, this induces a very fast motion of the node, which is difficult to be damped out. Additionally, the line discretization could cause a high frequency oscillation in the segment, especially when the adjacent nodes are out of phase with each other [16]. It is worth mentioning that for time integration scheme of the equations of motion, OrcaFlex implements a semi-implicit Euler scheme, one order lower than the RK2 scheme used in MoorDyn. The fundamental difference between the two schemes is that RK2 calculates an additional state (position and velocity in this case) in between two time steps (see Appendix A). For this reason, with the same constant time step, OrcaFlex is expected to give a numerically less stable result.

4.2. Simulation Case 2: Regular Wave and Seabed Friction

In simulation case 2, environmental input is kept the same as in case 1: Regular wave with a wave amplitude of 5m and period of 8.33 s. The same friction coefficients are set in both adapted MoorDyn and OrcaFlex to study the influence of sea bed friction, which is 0.5. In his implementation of sea bed friction modelling, Hall [41] used three different friction coefficients—0, 0.05, and 0.1. It was concluded that even with the coefficient of 0.1, the influenced is well pronounced in changing the motion of the modelled mooring line. DNV [40] recommends a friction coefficient of 1.0 for a chain in contact with sea bottom. In order to assess the influence of such parameter on the results of the simulations, a value in between the highest one tested by Hall [41] and the one recommended by DNV [40] was chosen. In a real case scenario, this parameter needs to be calibrated according to the sea bed characteristics.

Good agreement of buoy’s position in all direction can be found between OrcaFlex and adapted MoorDyn as it is shown in Figure 8. Overall the results between simulation case 1 and simulation case 2 are similar. When comparing overshoot values in anchor tension, however, the implementation of sea bed friction gives a reduction to these values. This is consistent in both adapted MoorDyn and OrcaFlex as it is shown in Figure 9 and Figure 10. As it was discussed in Section 4.1, the addition of sea bed friction would provide damping to the whole system, which in turns reduce the tension on the line. However, high overshoot values are still found in OrcaFlex results even when sea bed friction is modelled.

4.3. Simulation Case 3: Irregular wave

In simulation case 3, the unidirectional spectrum is used as an input. A list of wave components with different frequencies are generated based on JONSWAP spectrum and used as input for simulation case 3, which can be found in Appendix B. Time series of wave sea elevation at the origin (0, 0) is shown in Figure 11. All the wave components are propagating to the direction of the negative y axis.

In Figure 12, the positions of the buoy calculated by OrcaFlex and adapted MoorDyn are compared. Some discrepancies are found in the calculated horizontal positions of the buoy. The discrepancies found are likely because of the implementation of wave-induced loads in OrcaFlex and adapted MoorDyn, which becomes more apparent in the results of simulation case 5. In terms of anchor tension prediction, OrcaFlex shows overshoot values (Figure 13), which is not the case in the adapted MoorDyn.

4.4. Simulation Case 4: Current

The system is subjected to an environmental load of current propagating to the direction of the positive x axis. Current is kept constant throughout simulation time with the magnitude along the depth modelled using power law. The current magnitude at the surface ($s_0$) is set to be 1 m/s, while the exponent alpha ($\alpha$) is 1/7. Figure 14 shown the current vertical profile used in both codes.

Perfect agreement between OrcaFlex and adapted MoorDyn are found when comparing positions of the buoy (Figure 15). There is a minor discrepancy when comparing anchor tension. As it is shown in Figure 16, OrcaFlex gives slightly higher tension.

4.5. Simulation Case 5: Irregular Wave and Current

In simulation case 5, wave input is kept the same as the one in the simulation case 3. An environmental load of current is now added propagating in the negative y-direction. Current magnitude along the depth is shown in Figure 14 of simulation case 4. The mooring state at t = 143.9s is visualized in the following Figure 17.

In simulation case 4, it was shown that both codes are in perfect agreement when predicting the buoy position, due to current-induced load. In general, good agreement of buoy positions is found between OrcaFlex and adapted MoorDyn as it is shown in Figure 18. The addition of current makes the chain reached its maximum length sooner in comparison with irregular wave only in simulation case 3. This makes the wave drift non-apparent in simulation case 5. In Figure 19, anchor tension calculated by OrcaFlex and adapted MoorDyn are compared. Overshoot values in OrcaFlex results are still apparent. However, compared to the simulation case 3, these values are much less. This is due to the current-induced load gives raised to tension in the chain, reducing the occurrence of snap loads.

4.6. Summary of Anchor Forces

The anchor forces for all simulation cases are summarized in the following

Table 7. Summary of anchor tension for all simulation cases.

Simulation Case	OrcaFlex min	OrcaFlex mean	OrcaFlex max	Adapted MoorDyn min	Adapted MoorDyn mean	Adapted MoorDyn max
1	0 N	2986 N	79,442 N	0 N	2355 N	16,597 N
2	0 N	2696 N	68,756 N	0 N	2609 N	12,672 N
3	0 N	1605 N	58,753 N	0 N	293 N	14,024 N
4	0 N	1406 N	1558 N	0 N	1380 N	1535 N
5	0 N	5058 N	71,766 N	0 N	3319 N	24,173 N

:

5. Conclusions and Future Work

The flexibility of open-source code is a big advantage to have an efficient code that can adapt to the needs of a user, in this case, aquaculture purposes. The original MoorDyn was partially capable of modelling a mussel line system, and was adapted to cope with 3DOF buoy objects and environmental induced loads (waves and current). The adaptation is focused on modelling the mussel line test setup described in the introduction, which uses a longline system. The modified code can calculate the reactions of an arbitrary mussel line system with slender cylindrical floaters to the environmental loads of wave and current. Based on the calculated forces and displacements, a mussel line system can be optimized for a real case scenario. From the viewpoint of placing mussel line systems in wind farms, the system has to be restricted in its movement based on very strict safety distances from the wind farm equipment. The model can be used to simulate the behavior of mussel line systems in an arbitrary location. In order to do that, an analysis of the wave and current characteristics of the location of interest has to be performed. The environmental loads will define the forces on the system and its maximum displacement, which can be used to design an optimal configuration. The modelling of a mussel line system with floaters different from slender cylinders would require further adaptations to the code.

OrcaFlex has been used to compare the numerical calculations with some discrepancies found. Simulation case 1 and 2 shows that the snap loads phenomenon is captured; however, the magnitudes calculated by OrcaFlex and adapted MoorDyn do not show the same results. Looking at the high overshoot values, OrcaFlex with bending and torsion not modelled and explicit integration scheme chosen appear to give numerically unstable results in the occurrence of snap loads. A comparison with experiment is suggested to determine the real magnitude of these snap loads. In simulation case 3, wave drift on the buoy, predicted by the two codes, shows some discrepancies. Concluding the simulation case 4, OrcaFlex and adapted MoorDyn predicted drag force, due to currents that are in perfect agreement. When the system is subjected to both irregular wave and current, the two codes show good agreement.

Some aspects of the code will be further improved in the near future: A more complex modelling of floating objects and seabed contact effects, a more efficient algorithm for wave calculations and a more computational efficient integration scheme. Finally, additional research effort will be put in the validation of the numerical model with full scale measurements. The adapted MoorDyn code is used in the framework of Edulis project to perform a validation study of the test setup line. Positions of the two SPAR buoys in the test setup line were recorded. A study is being done by comparing the positions of the buoys from measured data and simulation results of the adapted MoorDyn during a full tidal cycle of 13 h. The preliminary results showed that positioning the system in line with the direction of ebb and flood current gives the lowest maximum displacements of the system. Furthermore, forces are measured on one of the SPAR buoy’s shackle. A numerical calculation during 30 min of sea state with wave and current as the input will be compared with the measured forces. In the near future, the results of this on-going study will be published, and the validated code will be used to design a layout of mussel lines array for the location of the Belgian North Sea.