A Numerical Solution for Modelling Mooring Dynamics, Including Bending and Shearing Effects, Using a Geometrically Exact Beam Model

: During the operation of moored, ﬂoating devices in the renewable energy sector, the tight coupling between the mooring system and ﬂoater motion results in snap load conditions. Before snap events occur, the mooring line is typically slack. Here, the mechanism of energy propagation changes from axial to bending dominant, and the correct modelling of the rotational deformation of the lines becomes important. In this paper, a new numerical solution for modelling the mooring dynamics that includes bending and shearing effects is proposed for this purpose. The approach is based on a geometrically exact beam model and quaternion representations for the rotational deformations. Further, the model is coupled to a two-phase numerical wave tank to simulate the motion of a moored, ﬂoating offshore wind platform in waves. A good agreement between the proposed numerical model and reference solutions was found. The inﬂuence of the bending stiffness on the motion of the structure was studied subsequently. We found that increased stiffness increased the amplitudes of the heave and surge motion, whereas the motion frequencies were less altered. comparison to the analytical solution for the vibration of


Introduction
Floating devices are attracting increasing attention to meet the global demand for green and renewable energy sources. An important factor for the safety of devices such as floating offshore wind turbines (FOWTs) is the correct construction of the attached mooring systems. In comparison to traditional offshore systems, where mainly large but low-frequency motions occur [1], FOWTs require stronger fixing of the platforms in order to ensure safe operations with rotating turbines. This stronger coupling between the motion of the floating body and the mooring system leads to the occurrence of slack lines and subsequently snap loads.
The most efficient approach for modelling the influence of mooring systems on floating structures is based on the analytical catenary solution [1] for slack lines and the elastic spring solution for taught lines. Both approaches are generally restricted to their primary intended forms, and their validity is based on the assumption of small motions of the moored system. A more elaborate version of this class of simplified approaches is a quasi-static solution [2,3] which combines the flexibility of a numerical model with the efficiency of an analytical formulation. However, the quasi-static assumption is still a serious restriction that does not hold for floating offshore wind turbines. Here, the dynamics of the mooring lines, which are described by non-linear partial differential equations of second order in both space and time, have to be included in the analysis.
If bending stiffness is neglected, multiple numerical solutions are applicable. Amongst others, finite differences [4][5][6], finite elements [7][8][9] and lumped mass [10][11][12] methods have been proposed. Recently, Palm et al. [13] utilised a discontinuous Galerkin method for solving the system. Thus, snap loads could be modelled with a high order of accuracy. These loads necessarily follow a phase of local slack which is characterised by vanishing axial stresses, and thus an ill-posed system if no bending stiffness is considered [14]. As pointed out by Zhu and Meguid [15], the main mechanism of energy propagation in these low tension situations also changes from axial to bending dominant. Thus, the ability to capture the correct bending dynamics of slack cables in three dimensions is a required extension for accurate mooring modelling [8].
One of the first attempts to include bending stiffness into a mooring model was presented by Garrett [9]. He developed a three-dimensional finite element model assuming a linearly elastic and torque-free rod. The effects of rotary inertia and shear deformations were neglected. Similar assumptions were made by Burgess [16] for proposing a cable model based on finite difference methods. Here, Euler angles were applied to account for the rotation of the line. This implies numerical difficulties, known as the Gimbal lock effect, in practice. Palm and Eskilsson [14] extended their Galerkin method for bending stiffness by following the idea of Garrett [9] and Tjavaras [17]. Thus, the inertia effects of the rotation of the cross-section were neglected and no torsion or shear force was considered. In [8], a cubic cable element with the translational positions of the end nodes as unknown variables was proposed. This element choice and the assumption of no external moments allowed the reformulation of the conservation laws in terms of forces acting on the node points. The resulting formulation is more compact than traditional finite element formulations which are based on 12 state variables.
More elaborate approaches were developed based on beam models suitable for deformable one-dimensional elastic structures. Mooring lines, such as other marine cables, are considered as beams undergoing large displacements and rotations [18]. Euler-Bernoulli beam models can account for the axial stiffness, bending and torsion under the assumption of small displacements, and the fact that the cross-sections remain undeformed during bending. The linear Timoshenko-Reissner beam theory includes the shearing effects of Euler-Bernoulli models but holds validity only for geometrically linear cases [19]. The extension of the latter approach to the geometrically non-linear domain was introduced by Antman [20] and Reissner [21], which is also known as the Cosserat rod theory due to Cosserat and Cosserat [22]. These geometrically exact beam theories provide consistent strain measures, which arise from the equilibrium equations for the deformed configuration, irrespective of the magnitudes of the displacements and rotations. Simo [23] and Simo and Vu-Quoc [24,25] generalised the theory of Antman [20] to propose a geometrically exact beam model using a rotational vector for the parametrisation of the rotations and the displacement and rotational vectors as the primary unknowns of the system. Quan et al. [26] recently presented the application of this beam model to underwater cable dynamics, whereas Cottanceau et al. [18] introduced a quasi-static version for the simulation of flexible cables in air.
As pointed out by Zupan et al. [27], any choice of the rotational vector in the threedimensional space can ultimately lead to a singularity at some state of rotation. They overcame this numerical limitation using quaternions which represented a set of four singularity-free parameters. The proposed model solved the dynamics of geometrically non-linear beams in quaternion descriptions based on the Newmark integration scheme and the collocation method. Later, Weeger et al. [28] utilised NURBS curves and quaternions in combination with an isogeometric collocation method to solve the same set of equations. As an alternative, Lang and Arnold [29] discretised the quaternion-based set of equations using the finite difference method. Thus, a matrix-free formulation of the discrete system could be derived which can be easily solved with common ODE solvers. In [30], this model was successfully applied to the simulation of a large number of flexible rods in current flow using LES.
In this paper, the efficient Cosserat rod model of Lang and Arnold [29] is utilised to develop a mooring model with geometrically exact kinematics for the first time. Thus, the effects of shear force, torsion and bending on the behaviour of mooring systems can be analysed in contrast to previous attempts. The advantages of using their formulation and discretisation of the geometrically exact rod theory are the avoidance of the Gimbal lock effect and the efficient matrix-free formulation. However, solving the dynamics of the rotational degrees of freedom is necessarily connected with the disadvantage of increased computational costs compared to simpler mooring models. In particular, the small bending stiffnesses of mooring lines result in a stiff problem which requires small time steps. The model was implemented in the open-source CFD code REEF3D [31]. The code is based on a numerical wave tank which enables the application to fluid-structure interaction problems of complex free surfaces and moored, floating structures. The remainder of the paper starts by presenting the numerical model for the geometrically exact mooring model in Section 2. Then, multiple verification (Section 3) and validation cases (Section 4) are presented to analyse the accuracy of the approach and discuss its applicability. Next, we present the validated numerical model being applied to the simulation of a moored floating offshore wind platform using CFD in Section 5. Here, the influences of the bending stiffness on the tension forces in the mooring system and the motion of the structure are studied. The paper concludes with final remarks in Section 6.

Numerical Model
In the following, the numerical implementation of the proposed mooring model is described in detail.

Governing Equations
The deformation of a Cosserat rod is described by the translational motion of its centreline r(X, t), with X ∈ [0; L] being the Lagrangian coordinate along the line of length L, in the inertia system and the conservation of the angular momentum expressed through the propagation of the rotation matrices of the cross sections R(X, t) in time. The latter presupposes that each cross section remains plane but freely rotatable during the deformation. This allows for shear force and torsion through the time and space-dependency of R. On the basis of this description, a geometrically exact representation of the three-dimensional conservation laws of momentum can be combined to two one-dimensional equations for the translational and rotational accelerations,r andω [29]: with ρ s being the material density, A the area of the rod and f ext and m ext the external forces and moments. Further, the internal forces f and moments m are defined under the assumption of linear viscoelastic material as with the strain and curvature vectors and the constitutive matrices in the Lagrangian frame of the rod E is the Young's modulus; G is the shear modulus; I = diag(I X , I Y , I Z ) are the second moments of area; and c = (c t , c s,Y , c s,Z ) are the torsion and shear correction factors.
It is elaborated in [19] that an efficient representation of the rotational motion without Gimbal locking is given in terms of the unit quaternions: with q ∈ H, Re(H) = R, Im(H) = R 3 and ||q|| 2 = 1. Within the Hamiltonian quaternion algebra H, the conjugate q of the quaternion q is defined as and * indicates the quaternion multiplication with the quaternion matrix For a vector v ∈ R 3 , the multiplication with a quaternion is defined as Finally, the rotation matrix is expressed in terms of quaternions using with Thus, the system of Equations (1) and (2) can be reformulated as [19] r = 1 with the inverse mass matrix and the quaternion matrix of inertia and its inverse The curvatures in the internal moments (4) are expressed using quaternions as

Numerical Discretisation Using Finite Differences
Following the idea of Lang and Arnold [29], the line is split into N segments and a staggered grid is introduced with r and m defined at the segment edges and q and f defined in the segment centres (compare Figure 1). Two ghost points are included to incorporate the rotational boundary conditions as explained in Section 2.3. The translational Equation (14) is solved at the grid points i = 1, ..., N + 1 using central differences: The rotational Equation (15) is solved in the grid points i = 2, ..., N + 1 so thaẗ with the central differences The discretisation is completed if a finite difference formulation for the curvatures is given. The difficulty is the choice of an appropriate interpolation of two adjoint quaternions in a segment edge. In this paper, the formulation is chosen. It represents a midpoint finite difference approximation on the circle between the adjoint quaternions in H (see [29] for an illustration and proper derivation).
The propagation of the discrete system from the old time n to the new time n + 1 is calculated by converting the equations into a system of first-order equations and applying the third-order accurate TVD Runge-Kutta scheme [32]: Here, ∆t is the size of the time step, Φ is the state vector and L represents the vector of the right-hand side.

Boundary Conditions
The staggered grid introduced above requires a more elaborate definition of the boundary conditions. For the translational equations defined at the segment edges, prescribed locations or velocities can be simply imposed on the boundary point. This is for example used to prescribe the motion of the fairlead of the mooring lines connected to floating objects in the CFD solver. In order to prescribe the rotation q b at e.g., r 1 , the quaternion at the ghost point q 1 is determined from extrapolating q 2 and q b to q 1 . Following the principal idea of the spherical linear interpolation (SLERP) for quaternions [33], this extrapolation is calculated as (compare [34]) In contrast, free end boundary conditions require vanishing internal forces and moments at the last grid point. Thus, free rotational ends are simply imposed by setting m = 0 at the corresponding point. For the translational equations, this requires f b = 0. A linear extrapolation of the internal forces to the fictive ghost point behind the last segment centre, e.g., f N+2 , is employed by using Inserting this expression into the finite difference formulation at the last grid point results in the modified equation which implicitly respects the free boundary condition for the internal forces.

External Forces
The external forces acting on the rod are calculated using the assumptions of Morison et al. [35] for hydrodynamic transparent structures. Thus, they are defined as the sum of hydrodynamic and static forces. The hydrodynamic forces are approximated as the inertia and drag forces acting on the centreline of the structure. The velocity dependent drag forces are given as [36] with c d,n and c d,t the drag coefficients in normal and tangential direction, d the diameter of the rod, ρ the fluid density, v the relative velocity vector between fluid and structure, and n the unit normal vectors along r(X). The normal component is calculated as a function of the Reynolds number Re [37]: whereas the tangential component is taken as 0.5 [13]. The inertia forces are calculated using [38] with a f being the fluid acceleration, a the relative acceleration between fluid and structure and c m the added mass coefficients in normal and tangential directions, taken as 3.8 and 0 in this study [13]. The static gravity and buoyancy forces are calculated under the assumption of an equally distributed mass: g represents the gravitational acceleration vector typically pointing in negative z direction. In addition, the effect of the bottom on the mooring line is modelled using the contact force defined by Palm et al. [13].

Coupling to the Fluid Solver
The open-source CFD code REEF3D [31] was applied for the simulation of the fluid-structure interaction of the moored, floating offshore wind platform presented in Section 5. The code can solve the incompressible Navier-Stokes equations on a staggered rectilinear grid using finite differences. The free surface is modelled using the level set method of Osher and Sethian [39] and Sussman et al. [40]. The dynamics of rigid floating structures are coupled to the fluid solver using a continuous direct forcing method [41] and the incremental pressure-correction algorithm of Timmermans et al. [42]. The third-order accurate TVD Runge-Kutta scheme of Shu and Osher [32] is applied to propagate the momentum and free surface in time, and the fifth-order accurate WENO schemes by Jiang and Peng [43] and Jiang and Shu [44] are used for convection terms. The boundary conditions are enforced with a ghost point approach, and an n-halo decomposition strategy with three layers and the message passing interface (MPI) handles the inter-processor communication efficiently.
As previously described [3], the influence of the mooring system on the fluid was neglected due to the assumption of hydrodynamic transparency and the negligible disturbance of mooring lines on the fluid in its vicinity. As boundary conditions, free rotations but fixed translational motions were set. On the ground, the mooring line was fixed, whereas the fairlead followed the motion of the rigid floating body. The coupling between the mooring and fluid dynamics solver was implemented in a twoway explicit manner. In each time step, the fluid velocities of the old time step were linearly interpolated at the segment edges of the rod to calculate the external forces on the mooring line. External moments were not included in this study. The motion of the line was then advanced to the new time step, and the tension forces acting at the fairlead were added to the rigid body motion solver as external loads. Sub-iterations were required because the time step of the mooring solver is much smaller than the time step of the fluid solver. Following Palm [45], the fairlead boundary condition was interpolated in time to ensure a smooth transition between the old location and the new location. An illustration of the complete algorithm can be found in Figure 2.

Three-Dimensional Static Deformation of a Rod Under Constant Loads
The convergence of the implemented Cosserat rod model was first assessed with the static test case proposed in [46]. As the initial condition, a squared rod was defined along one quarter of a circle with a radius of 1 m without pretension. The geometrical and material properties are stated in Table 1. One end was fixed for both translational and rotational motion, whereas the other end was free. Two constant loads of 3 and 6 N acted on the free end perpendicularly to the initial rod configuration. Thus, a three-dimensional nonlinear deformation of the rod was initiated at t = 0 s, and a converged static solution was calculated using a small portion of internal damping to speed up the rate of convergence. Figure 3 shows the L 2 norms of error for N between 5 and 160 and the two load cases. An even finer solution was taken as the reference for this calculation. The log-log plot reveals a second-order accurate solution in all directions, which is in accordance with the theoretical accuracy of the chosen discretisation scheme.

Rigid Pendulum with Continuous Mass Distribution
The presented formulation for the dynamics of a rod is capable of predicting the motions of rigid bodies without further modifications. In order to verify this feature, the motion of a rigid pendulum with a continuous mass distribution in a gravity field is compared to the theoretical solution given in [34]. Very high values for E and G ensured that the pendulum moved as a rigid body. The pendulum was initially located on the horizontal plane perpendicular to the direction of the gravity vector. At t = 0 s, the pendulum was released, with one end being pinned and the other moving freely. A very small time step was chosen to remove its influence on the solution. Figure 4a shows the differences between the theoretical and numerically predicted angle of the pendulum over the first two periods. The error increases over time due to numerical dissipation but is reduced with an increasing number of grid points. In Figure 4b, the corresponding convergence of the L 2 norm of the time-integrated error is visualised. As before, the nominal order of accuracy is reached.

Forced Vibration of a Pinned Rod Under Static Force
The vibration of a rod under a static force is verified next. Assuming small deflections, the analytical solution for the position of the centre of the rod under a force on this point is found from the solution of the Euler-Bernoulli beam equations [47]: The geometrical and material properties were taken from [27] and are presented in Table 2. The initial rod was straight, and a static force of 1 N acted perpendicular on its mid-point. Both ends were pinned. The motion of the mid-point over 25 s for different numbers of grid points is shown in Figure 5. The maximum deflection of the beam was well captured by the numerical model for N > 32, whereas a small deviation of the frequency was present after about 6 periods even for the finest solution with N = 128. Following the argument of Tschisgale [34], this observation was due to the coupling of shear and bending modes in the Cosserat rod equations. This theoretical difference from the Euler-Bernoulli beam equations has to be taken into account when comparing a fully non-linear solution to (32). However, Zupan et al. [27] obtained better agreement with the analytical solution for their geometrically exact model using a finite element method. This indicates that the chosen staggered finite difference method might also cause the observed deviations.

Large Motion and Deformation of a Free-Free Rod
A typical validation case for nonlinear beam models is the two [24] and three-dimensional [25] for the two-dimensional case and as for the three-dimensional configuration. The geometrical and material properties are summarised in Table 3. The deformed rod at different time instances for N = 10, 20 and 40 is shown in the Figures 6 and 7 in comparison to the numerical solution of Hesse and Palacios [48]. Figure 6 shows the 2D solution at different time instances between 0 and 10 s. For N = 10, minor differences to the reference solution were observed after 8.0 s, whereas qualitatively good agreement is given for the finer grids. As stated in [48], potential differences were due to different time stepping because this test case imposed up to the fourth bending mode. In comparison, the results for the three-dimensional load case in Figure 7 show better agreement with the reference solution. An explanation for this might be that this case was only presented for the first 6.5 s of the simulation. Later time instances might reveal larger discrepancies as well. Table 3. Geometrical and material properties of the free-free rod undergoing large motions and deformations.

L [m]
A    Quantification of the simulation results is provided in Table 4. The z-position of the central point of the rod was measured at the end of the simulations for different grid sizes. As the forces and moments were spatially attached to the end of the rod, this point was theoretically not moving in z direction during the deformation. Thus, the corresponding error of the z-coordinate at this point can be used as a good measure for the performance of the model for long-duration integrations. It is evident from the table that the solutions converge towards the analytical result with an increasing number of grid points. For the two-dimensional case, the percentage error decreases from 0.76% for N = 4 to 0.32% for N = 128, whereas the minimum error is 2.8% for the three-dimensional case. This increase in error might be attributed to the increased complexity of the solution when adding the third dimension.

Swinging Motion of Pinned Elastic Rods with Free Ends
For the next validation case, a similar rod as in Section 4.1 was pinned at one of its ends and subject to gravity forces. Thus, a swinging motion with slack and tense situations was expected. The positions of this rod at different time instances are shown in Figure 8 for the geometrical and material properties given in Table 5. The simulated result was obtained with 10 grid points and compared to the solution using the same number of elements in the FEM toolbox Abaqus (reported in [29]). Overall, very good agreement with the reference results was achieved. This was quantified by plotting the temporal evolution of the tension forces at 0.4 and 0.8 m of the rod (see Figure 9). At X = 0.4 m, the tension force increased continuously while the rod was moving downwards due to the influence of gravity. At around t = 0.45 s, the forces started to oscillate due to the shock-like waves propagating along the rod. The tension force decreased significantly when the rod approached its point of return. Compression effects were present. All phenomena were well captured by the implemented model. Table 5. Geometrical and material properties of the swinging elastic rod.

L [m] r [m] A [m 2 ] ρ [kg/m 3 ] E [N/m 2 ] I Y [m 4 ] G [N/m 2 ] g [m/s 2 ]
1.0 0.005 πr 2 1100 5 × 10 6 Further assessment of the model's capabilities is provided by comparing to the time series of tension forces in a swinging rod measured by Koh et al. [49]. In comparison to the setup from above, the initial configuration of the rod was a hanging shape which was calculated from a catenary solution [50]. The setup for the simulation is given in Table 6. Palm and Eskilsson [14] pointed out that the experiment lacks the measurement of the damping properties for the bending modes. Therefore, they conducted multiple simulations with different values for Cκ ,0 (2, 2) and found the best agreement with the experiment when choosing a value of 0.02 Nsm 2 . Similarly to their results, the proposed model is to a large extent not sensible with regard to this value as long as it is not zero. As can be seen in Figure 10, the numerical result obtained with N = 15 grid points agrees well with the experiments with regard to the period and the peaks of the tension forces. At the minima, which are characterised by compression forces, the simulation tended to over-predict the peak values by about 20%. A possible explanation is given by Palm and Eskilsson [14], who suggested that the frictional forces at the pivot pin cause a reduction of the tension forces. In contrast, this influence was not accounted for in the simulations. We also noticed that a further increase of the bending damping coefficient to 0.04 Nsm 2 provided slightly better agreement with the experiments.

Dynamics of a Catenary Chain in Water
As a final validation case, the proposed mooring model is compared to the experiment of a moving catenary chain in water presented in [51]. The chain of L = 33 m was placed in a water tank with a depth of 3 m. One end of the chain was fixed to the bottom at (0, 0, 0), whereas the other end was subjected to a circular motion with a radius of 0.2 m around the point (32.554, 0, 3.3). The motion had periods of T = 1.25 and 3.5 s. The line is discretised by 20 points, and the geometrical and material properties are summarised in Table 7. It is to be noticed that the chain did not allow for compression effects. This was respected in the model by setting the internal forces to zero if the forces along the centreline of the rod were locally negative.  Figure 11a presents the time series of the tension force magnitudes for the motion with T = 1.25 s. The numerical model is capable of reproducing the physical frequency and peak values of the measured tension forces. In comparison to the experiment, the model predicted longer intervals of low tension due to the inclusion of bending stiffness and viscous damping in the tension force direction. The shock, induced by the sudden lifting of the completely slack chain, resulted in a strong increase of the tension forces. The numerical solution overshot at a saddle point close to the peak due to the additional momentum in this model, which altered the shock speed. Further insights into this behaviour were given by simulating this case with increased bending stiffness and without rotational motions. This ideal case was calculated by neglecting the rotational equations and replacing the calculation of the inner forces with the calculation provided in, e.g., [13]. A small amount of viscous damping was required to avoid unphysical oscillations of the numerical solution. As can be seen in Figure 12, the results without bending stiffness show less overshooting and undershooting of the experimental distribution, particularly at the saddle point. At the same time, oscillations at about half the maximum tension force arose due to the ground forces. In contrast, the increase of the bending stiffness introduced additional saddle points and also delayed the occurrence of the peak values. This substantiates the postulation that the bending stiffness is the main cause of the discrepancy between numerics and experimentation. Better agreement with the experiment is given for the fairlead motion with T = 3.5 s in Figure 11b. This was probably caused by the slower motion of the chain. Thus, the additional inertia of the rotational motions did not alter the principal propagation of the shocks alongside the chain.

Application to the Simulation of a Moored Floating Offshore Wind Turbine Support Structure
The validation case in Section 4.3 reveals that the proposed mooring model is capable of predicting the correct tension force propagation to a large extent. The main advantage of this model is, however, the capability of incorporating shear, bending and torsion effects into the analysis. The influences of these effects on the expected tension forces and motions of a moored floating offshore wind turbine support structure are analysed in the following.
The chosen design was the well-established DeepCwind OC4 semi-submersible FOWT [52] at a 1:50 Froude scale taken from [53]. The sub-structure consists of three vertical columns connected to a more slender central column using multiple thin braces. Heave plates are attached to the bottom of the columns to increase stability. The platform is held in place through a mooring system consisting of three catenary lines spread symmetrically about the vertical axis of the structure. The fairleads of the lines are attached to the top of the heave plates, whereas the bottom ends of the lines lie on the sea ground. The stiffness matrices required for the present model were calculated from the given information about the elasticity module and second moments of area calculated from the assumption of a cylindrical shape. The reference solution was taken as the mean of the various numerical models presented in [54].
The structure was placed in the middle of a numerical wave tank, as shown in Figure 13. The tank had a length of 24 m, a width of 8 m and a height of 6.5 m, and the water depth was 4 m. It can be noticed that the lower part of the mooring system was placed outside the computational domain to lower the computational costs of the simulations. This is also justified by the observation that the fluid velocities close to the bottom of the tank are small and a large portion of the mooring lines lies on the ground. Thus, the fluid velocities outside the computational domain but needed for the external force calculation of the mooring system were assumed to be zero. Near the inlet of the domain, a wave generation zone was defined using the relaxation method [31]. A numerical beach at the end of the domain damped the wave energy to avoid reflections. A rectilinear grid was used to discretise the numerical domain. As can be seen in Figure 14, a refinement box with uniform grid sizes was defined around the structure. Cells of linearly increasing size were placed between the box and the domain boundaries to ensure a smooth transition. The grid size in the inner domain was determined from decay tests in heave and surge, presented next.  At first, a free heave decay test was conducted for the moored structure using three different grid sizes in the refinement box. The grid size outside of this inner region was automatically calculated based on linear grid growth with a ratio of 1.1. Figure 15 presents the time series of the heave motion in comparison to the numerical solution in [54]. At the first peak after 1.2 s, the numerical solution converged towards the reference, whereas all grids predicted the correct amplitude at the subsequent peaks. Further, the present model tended to predict a slightly longer natural period than the reference. Next, the grid with a minimum distance of 0.03 m was chosen for the surge decay test (Figure 16). Both the amplitudes and the frequency of the motion were well predicted by the model. Therefore, this grid size was also used for the simulations below. The corresponding mesh had approximately 4.7 million grid points.   The accuracy of the numerical setup was assessed using a regular wave case with a wave height of H = 0.12 m and a wave period of T = 1.41 s. This corresponds to a wave length approximately twice the structural length, as can be seen in Figure 17. The wave was modelled as a second-order Stokes wave in the simulation. The qualitative comparison of the three body motions to the reference solution can be found in Figures 18-20, whereas the quantification of the mean amplitude and frequency is shown in Figures 21-23. The latter parameters were found from fast Fourier transformations of the steady-state time signals of the motions. Generally, good agreement between the presented and reference simulation can be stated. The frequencies of the three motions are close to the encountered wave frequency of 0.71 Hz. Further, the investigated model predicted lower motion amplitudes for the positive peaks. This is particularly prevalent for the surge motion, which indicates that the front mooring line might have caused this difference. The time series of the tension force magnitudes at the fairlead of the front mooring line is presented in Figure 24 to investigate this. In comparison to the reference solution, similar frequencies and crests of the forces but lower force troughs were predicted. This difference might be the reason for the slightly different motions of the structure.      Further insights into the influence of bending stiffness on the motion of moored, floating structures was provided by additional simulations with increased stiffness values. Figures 21-23 present the amplitudes and frequencies of the motions for different bending stiffnesses. The reference solution from above, which is only available for the original bending stiffness, is also shown for comparison's sake. It was first noticed that the frequencies of the heave and surge motions seem not to have been influenced by the increased stiffness. In contrast, the amplitudes of these motions increased by up to 40% and 70% respectively, and the pitch amplitude decreased by up to 30%. The reason for this was the additional rotational momentum, which restrained the rotation of the structure. Thus, the wave energy was rather converted into the translational modes. This also involved an increase of the mooring tension forces. It was finally observed that the bending stiffness led to an increased pitch frequency due to the larger system stiffness.

Conclusions
A new approach for simulating the dynamics of mooring systems and their interactions with floating structures was presented in this paper. It was based on a geometrically exact beam model originally developed for slender flexible rods. An efficient solution for the arising system of equations was found from a quaternion description of the rotational deformations and a finite difference discretisation. Explicit time integration was utilised, in contrast to previous approaches. The fluid forces acting on the system were calculated from Morison's formula and the assumption of hydrodynamic transparency. The resulting mooring model has the advantage over existing models that it can account not just for bending but also shearing and torsion effects. Several verification and validation cases showed that the model is of second-order accuracy and can accurately represent the structural deformations and tension force distributions of rods and mooring lines.
The mooring model was then coupled to a two-phase CFD solver to investigate moored, floating structures in waves. The applicability of this approach was shown for the semi-submersible floating offshore wind support structure of the DeepCwind OC4 project. Initial decay tests in heave and surge were conducted to study the convergence of the numerical model and determine the required number of grid points. Then, the motion of the structure in a regular wave was compared to previous numerical solutions. Good agreement for the amplitude and frequency of the heave, surge and pitch motion was shown. A study of the influences of increasing bending stiffness on these motions revealed only a minor change to the motion frequencies but significant effects on the amplitudes. The additional rotational momentum decreased the rotation of the structure so that the wave energy acted more strongly on the translational modes.
The presented development of a dynamic mooring model is only one possible application for geometrically exact beam theory. Within further research, the same model will be applied to other slender marine structures, such as monopiles, offshore wind turbine towers and floaters. Additionally, electrical submarine cables might be of interest, as they do not only experience axial elongation but also torsion in the helical armour [55]. Further, the comparison to the analytical solution for the forced vibration of a pinned rod revealed a principal deviation of the chosen numerical approach. Previous research showed better results using finite element methods. Within further developments, it has to be investigated what causes the problems of the finite difference discretisation and whether the deviations can be avoided by switching to a finite element solution.