Accurate modeling of the interaction of constrained floating structures and complex free surfaces using a new quasistatic mooring model

A new approach for the coupled simulation of moored floating structures in a numerical wave tank using CFD is subject of this article. A quasistatic mooring model is derived by dividing the cable into finite elements and solve force equilibria in each time step. A successive approximation is applied to solve the problem most efficiently. Here, the unknown internal and external forces are separated, and the system is corrected iteratively using the intermediate results for the unit vectors until convergence is reached. An improved algorithm based on the directional ghost‐cell immersed boundary method is utilized for modeling floating objects in a three‐dimensional numerical wave tank. It is based on an implicit representation of the body on a stationary grid using a level set function. The motion of the rigid body is described using Euler parameters and Hamiltonian mechanics. Dynamic boundary conditions are enforced by modifying the ghost points in the vicinity of the structure. A special procedure for staggered grids is included in the article. The mooring model is coupled to the fluid‐structure interaction solver to simulate physically constrained floating bodies in waves. Several validation cases provide an overview of the accuracy of the mooring model and the floating algorithm. In addition, the comparison of the numerical results for the motion of a moored floating barge in regular waves with experiments is presented.

Taking the Navier-Stokes equations as the basis of the calculations, several attempts have been proposed to solve fluid-structure interactions. In Arbitrary Lagrangian-Eulerian methods, 1 the fluid properties are described on a fixed Eulerian grid which is fitted to the moving body. As the body changes its position, the mesh has to be recalculated dynamically. This approach raises the theoretical possibility to resolve the boundary layer around the structure by applying fine layers of cells but is prone to loss of numerical accuracy and stability due to suboptimal remeshing, especially for large body motions. To prevent irregular grids, remeshing can be avoided using either completely Lagrangian methods such as meshless smoothed particle hydrodynamics models 2 or dynamic overset grids. 3,4 In overset methods, several moving meshes which are arranged on a fixed base grid are considered. Hence, body-fitted cell alignments remain possible. An interpolation mechanism handles the data transfer between the overlapping grid points, and an overall solution for the fluid properties can be found by including these interpolations into the system matrix. 5 The size of the overlapping region is related to the stencil size of the discretization. This complicates the implementation and affects the efficiency for high-order numerics. Alternative approaches are based on immersed boundary methods. 6 They require just one Eulerian grid, and the body is described implicitly. Direct forcing methods 7 incorporate the fluid-structure interaction as an additional source term in the Navier-Stokes equations. The term represents the amount of momentum between the actual fluid-solid interface and the nearest cell centers under consideration of the boundary conditions. Evaluating this term usually includes several iterations with multiple pressure solution processes. A noniterative and more efficient approach was recently presented by Yang and Stern. 8 Following the work of Calderer et al, 9 an easier way of respecting the boundary conditions can be found. Instead of an additional term, modified stencils are considered. In the vicinity of the body, interpolation stencils are modified by including their closest distance to the body. Berthelsen and Faltinsen 10 further simplified this procedure for fixed bodies by using the closest distances in the local directions of the stencils. In doing so, the necessary interpolations become straightforward evaluations of Lagrangian polynomials in the principal directions of the grid. Bihs and Kamath 11 utilized this advantage to develop a local ghost-cell immersed boundary method for moving bodies. Their research was mainly focused on falling objects, and the algorithm was based on first-order interpolations and force and geometry property calculations using the level set representation of the body. A possibly underresolved level set function has a significant effect on the accuracy of the results, which is particularly severe for floating objects depending on the equilibrium of forces. Therefore, this article presents several improvements of the original solver by reducing the dependency on the level set function to the most useful parts. Thus, the simulation of floating objects in complex waves is enabled for the first time using a local ghost-cell immersed boundary method.
The dynamics of mooring lines is described by a nonlinear partial differential equation of second order in both space and time even if bending stiffness is neglected. Therefore, the general solution can just be found numerically. Davidson and Ringwood 12 present a review of a wide range of possible discretization methodologies. They are based on finite differences, 13,14 finite elements, 15 or spring elements. 16 Palm et al 17 pointed out the hyperbolic nature of the equations and applied an hp-adaptive discontinuous Galerkin method. Their method is stable, of arbitrary order of accuracy and capable of simulating snap loads. The main disadvantage of fully dynamic mooring models is the necessity of stable initial conditions and the restriction of the time step imposed by the material stiffness. Thus, the time step due to the mooring dynamics can be several magnitudes smaller than the time step due to the CFL condition of the fluid solver. This either restricts the efficiency of the whole simulation or raises the need for interpolation in time. 18 If the purpose of the simulation is more concentrated on the FSI problem and less on the details of the mooring dynamics or the expected motion of the lines are small in any case, simplifications of the original equations are appropriate. By neglecting the dynamic effects, dependencies of mass, damping and fluid acceleration on the system are omitted. In this case, the mooring line shape and tension is often found analytically using the catenary solution. 19 This solution is restricted in its form and not suitable for tense or arbitrary shaped systems. Therefore, a quasistatic solution presents a suitable compromise because it combines the flexibility of a generically formulated numerical approach with similar efficiency and simplicity as an analytical solution.
This article aims at developing such a quasistatic mooring model, which is disregarded in research so far. The finite element method for tensile structures, 20 which is based on a discrete representation of the structure using knots connected with trusses, is taken as the basis of the presented developments. Here, a linear system of equations is generated using force equilibria at each knot of the discrete line and a geometrical constraint preventing unphysical correlation of tension forces and bar deformation. Important hydrodynamic effects on the lines are respected using Morison forces from the surrounding fluid solution. The presented mooring model is coupled to the new floating algorithm within the framework of the open-source CFD solver REEF3D. The coupling process is straightforward and adds significantly less computational time to the overall simulation than dynamic models due to the missing time step restriction. Furthermore, the transition from a hyperbolic problem based on initial and boundary conditions to a boundary value problem makes the new mooring model a more flexible approach.
In Section 2.1, the CFD model is briefly described. Afterward, details about the improved floating algorithm and new mooring model are given in Sections 2.2 and 2.3. Section 3 presents several verification and validation cases. First, Section 3.1 provides the validation of the new mooring model against measurements of the prescribed motion of a single line and wind channel experiments with a line mounted between two fixed points. Afterward, the rigid body dynamics solver is investigated and the complete FSI solver is validated for fluid-structure interactions of a floating and moored floating barge in regular waves in Sections 3.3 and 3.4. The article concludes with an application to a more complex moored semisubmersible platform in Section 4 and final remarks in Section 5.

Computational fluid dynamics solver
In this article, the incompressible continuity and Reynolds-averaged Navier-Stokes equations are solved in the whole domain. They are given as Here, u i with i = 1,2,3 are the velocity components in the coordinate directions, is the fluid density, p presents the pressure contribution, and t are the kinematic and turbulent viscosity, and g i is the gravity acceleration vector. The additional viscosity is evaluated from the k-turbulence model with a modified source term to account for free surfaces. 21 The material properties of the two phases are determined for the whole domain in accordance to the continuum surface force model of Brackbill et al. 22 The equations are discretized on a rectilinear and staggered grid using a finite difference method. Chorin's projection method for incompressible flows is applied to resolve the velocity-pressure coupling. 23 The pressure is obtained as the solution of a Poisson equation from the continuity equation. The fully parallelized BiCGStab algorithm 24 with PFMG preconditioner 25 from the iterative linear system solvers of HYPRE are exploited. Temporal terms are discretized applying the third-order accurate total variation diminishing Runge-Kutta scheme. 26 An adaptive time-stepping control according to the CFL condition is used. Diffusion terms are treated implicitly to overcome their restrictions on this condition. The Convection term in (2) is discretized in a nonconservative form to increase numerical stability. 27 The fifth-order accurate weighted essentially nonoscillatory (WENO) scheme of Jiang and Shu 28 is chosen to reconstruct the fluxes. The scheme is adapted to nonconservative terms under consideration of the work of Zhang and Jackson. 29 The free surface is implicitly captured by the level set method of Osher and Sethian. 30 A smooth signed distance function , defined as the closest distance to the interface, is convected in the fluid velocity field using the linear advection equation The convection term is discretized by the fifth-order accurate Hamilton-Jacobi WENO method of Jiang and Peng. 31 In order to conserve the signed distance property, the level set function is reinitialized after each time step using the PDE-based reinitialization equation of Sussman et al. 27 All simulations are executed in the numerical wave tank of the open-source CFD solver REEF3D. 21 Waves are generated by applying the relaxation method to the x-and z-component of the velocities and the level set function . For a general variable it is given as Reference 32 with Γ the relaxation function taken from Reference 33. The same method is applied to damp potential reflections near the outlet of the tank. Extensive validation of the described numerical wave tank can be found in Reference 21.

Improved rigid body dynamics and fluid-structure coupling
In the following, the rigid body dynamics solver and its coupling to the fluid model is described. In comparison to the previous algorithm, 11 the improved model calculates geometrical properties of a rigid body from a triangulated surface instead of a level set function approximation. This can be computed more efficiently by surface integrations if the assumption of a homogeneous material with density s holds. Given the approximation of the bodies' surface S as where S i is assumed to be triangular, the properties are defined as Here, ∇ ⋅ f (⃗ x) = 1 returns the body mass, ∇ ⋅ f (⃗ x) = x, y, z present the center of mass coordinates in an inertial system and ∇ ⋅ f (⃗ x) = x 2 , y 2 , z 2 , xy, yz, xz compute the moments of inertia. Furthermore, ⃗ n i is the face normal vector of triangle S i , and ⃗ k is a unit vector in the x, y, or z-direction dependent on the chosen integration of the function. Numerical integrations in (6) can be avoided by parameterizing the triangles and analytically integrating over a standard triangle as given in Reference 34.
The translational motion of the rigid body is described by Newtons second law as The position ⃗ x s and velocitẏ⃗ x s of the bodies center of gravity can be determined by numerically integrating (7). The orientation in a fixed coordinate system is described by the four Euler parameters ⃗ e = (e 0 e 1 e 2 e 3 ) T which are related to the physically more relevant Tait-Bryan angles Φ, Θ, Ψ via (c is cos and s is sin) 35 In practice, this choice overcomes eventual issues due to the Gimbal lock effect. The back-transformation can be calculated using The four Euler parameters are depend as ⃗ e T ⃗ e = 1. The transformation of a vector in the body-fixed coordinate system to a corresponding vector in the inertial system is described by the orthogonal rotation matrix Based on this, the kinematic equations for the rotational motion of the body in terms of the Euler parameters are given aṡ⃗ with ⃗ the components of the angular velocity vector in the body-fixed coordinate system and Introducing the momentum vector ⃗ h = I ⃗ with I the moment of inertia tensor, (16) can be rewritten aṡ Following the derivation of Shivarama and Fahrenthold, 36 a first-order ODE can be derived foṙ⃗ h using a system Hamiltonian. Here, the constraint for the Euler parameters is fulfilled automatically. By setting the potential energy function to zero and assuming imposed moments ⃗ M ′ ⃗ x in the body-fixed system, the equations reaḋ The moments are obtained from the inertia system using the transformation matrix (15). In total, (7), (16), and (19) forms a system of 13 first-order ODEs which can be solved using any type of numerical integration method. An advantage of the derivation based on Hamiltonian mechanics is the energy-based formulation. Thus, conservation of energy can be tracked easily in the absence of external forces. This is used later to choose a suitable numerical integration method.
An implicit description of the floating body is found efficiently by defining a level set function such that the zero level set describes the surface of the body using the triangulated body representation. A ray-tracing algorithm 37 is applied first to get the shortest distances in each coordinate direction between the zero level set and each grid point. The signed distance property for the level set function is ensured by the above-mentioned reinitialization process. The method avoids explicit calculations for getting the intersections between the contour and the grid but lacks accuracy for calculating the surface area. Thus, the algorithm determines the body forces by applying a linear interpolation between the grid point values and integrating over the discrete surface with ⃗ n the surface normal vectors, the dynamic viscosity, ⃗ the viscous stress tensor, and N the number of surface representing triangles. Furthermore, ΔS is the area of each triangle. This is in contrast to the previous force calculation based on the level set and the Dirac delta function. The proposed FSI model is a weakly coupling algorithm with the possibility for subiterations. As already indicated in Reference 11, the overall accuracy or stability is only marginally influenced by these subiterations and so can usually be omitted. Key steps in the new algorithm are the ghost-cell interpolations, which are adapted from the ghost-cell immersed boundary method of Berthelsen and Faltinsen, 10 to incorporate the boundary conditions of the solid. For this purpose, the staggered evaluation points for the velocity are defined as either fluid, interpolation, or ghost points as indicated in Figure 1A. Pressure points, which coincide with the cell centers of the original grid, are defined as fluid points if they are in the fluid and as ghost points elsewhere. The pressure values of ghost points are extrapolated in the different coordinate directions from calculated fluid values and the boundary condition for the pressure gradient at the interface which reads with u i,Γ the bodies' velocities at the interface which are given in Reference 11. Equation (21) simplifies to the common zero gradient boundary condition for nonmoving bodies. The pressure ghost cells are used during the solution of the Poisson equation and force calculations. A slightly different approach is implemented for the velocity components on the staggered grids. As can be seen in Figure 1A, interpolation points are defined as an additional class of points. The velocity values at these points are determined from a linear interpolation using a fluid point and the known velocity of the solid at the interface in the principal direction (see Figure 1B). Thus, the correct boundary conditions are implicitly respected by the momentum calculations in the predictor step. In comparison to other immersed boundary methods, this directional approach avoids modification of the fluid stencils or the setup of additional stencils in the vicinity of the body. The indicated problems of solid cells becoming fluid cells and the resulting missing velocity information 38,39 are circumvented by assigning a velocity to the ghost points in the body. This velocity equals the body velocity at the corresponding ghost point. The presented approach is easily extendable to higher order interpolations but preliminary tests did not show significant improvements in accuracy or stability.

New quasistatic numerical approach for solving mooring dynamics
The development of the new quasistatic mooring model based on finite elements which were originally developed for quasistatic modeling of more complex tensile structures. 20,40 The dynamics of a mooring line neglecting bending stiffness is described as Reference 17 with the specific weight of the material, F T the magnitude of the tension force, ⃗ f the unit vector pointing in the direction of this force, and ⃗ F e external forces including gravitation and hydrodynamic effects. Assuming small line motion in time and steady-state flow of the fluid, respectively, (22) simplifies to In order to solve this force equilibrium, each mooring line is split into N equally distanced bars of length l t with knots P in between. As shown in Figure 2, the first and last knot, P (0) and P (N) , are attached to the bottom and the floating object.
The mass of the line is equally distributed on the adjacent knots which gives a gravity force contribution of at any knot P (j) . Here, m is the density of the mooring material and ⃗ g is an unit vector pointing in negative z-direction. Hydrodynamic forces ⃗ F H , arising from the slowly varying relative motion between structure and surrounding fluid, are calculated as drag forces using Morison's formula at each bar with c n and c t the drag coefficients in normal and tangential direction. The hydrodynamic forces are assigned to knots by equally distributing the net amount. As indicated in Reference 20, hydrodynamic forces can also be written as the partial sum of hydrodynamic drag, lift, and shear forces. The coefficients then depend particularly on the local angle between the fluid and bar unit vector and have to be found experimentally for each mooring configuration. However, using (25) with coefficients for slender cylindrical shapes is more practical for general applications. Choo and Casarella 41 derived a formula for c n as a function of the Reynolds number For c t , typical values can be found in literature. Furthermore, it is given that the tension forces act at the knots in the direction of the adjacent bars such that (23) can be solved locally at each inner knot P (j) (see also Figure 3): A solution for the shape of the line and the distribution of tension forces can be found by gathering (27) into a suitable linear system of equations. Both sought unknowns are generally dependent on the direction of the bar unit vectors. However, ⃗ F H also depends on ⃗ f according to (25), whereas ⃗ F G is related to the length of the bars which is a function of tension forces. This elasticity of the material is respected by representing l t as a polynomial of F T such as F I G U R E 3 Force equilibrium at inner knot P : Inner knot (white filled circle), bars (thin vectors), forces (thick vectors) with K p the stiffness constants and k the iteration index. Hooke's law is considered in the applications shown below by only using K 1 . An iterative method has to be chosen for solving the system because all components of (27) are related to each other. As proven by Hackmann, 42,43 a successive approximation can be found to ensure a converged solution. For this purpose, (27) is rearranged such that internal and external forces are separated Figure 2 indicates that the discrete line consists of N bars but just N − 1 inner knots. Therefore, (29) presents an undetermined system if solved for ⃗ f . A geometrical constraint is introduced in order to close the system. This constraint can be chosen in such a way that it represents both the physical boundary condition of fixed end points P (0) and P (N) , and the physical coherence of the line during deformation: This means that the vector of the sum of all bar vectors has to be equal to the distance vector between the two end points ⃗ x(P (0) ) and ⃗ x(P (N) ), which is a conditional equation for a physical coherent solution of the problem. The resulting linear system of equations can now be formulated with the unit bar vectors as the (N × 3) solution matrix F because of the dependencies in (30). Under consideration of (29) and (30), the system reads The solution process of the mooring model is shown in Figure 4. As input parameters, the fluid velocity field and the current locations of the mounting points are communicated to the algorithm. The geometrical constraint vector is calculated from these and an initial system of equations is generated. For the first time step, less restrictive directions for the unit vectors and initial tension forces are set to fill the left-and right-hand sides. Thus, this model is independent of a predefined initial form of the cable and highly suitable for applications without this requirement. Based on the given values, (28) returns the initial lengths of the bars and (24) yields the gravity force matrix G. The hydrodynamic force matrix H is initialized considering the fluid velocity field and (25). The velocities at the bars are calculated from the surrounding velocity cells and a weighted interpolation. This process is fully parallelized using MPI. The solution of (31) at any iterative step k is determined by Gaussian elimination with pivoting as The lengths of bar unit vectors have to be equal to one by definition after being calculated from (32) which corresponds to the conservation of all bar unit vectors within a tolerance typically chosen as 10 −4 . In this article, the influence of the mooring system on the fluid is neglected because hydrodynamic transparency is assumed. This is justified by the focus on the motion of the body, which is not affected by the fluid disturbances due to the presence of the mooring lines. Thus, the communicated data from the mooring algorithm to the FSI solver only includes the tension  Figure 5). Alternatively, (29) can be solved using Newton-Raphson iterations. As shown in Reference 40, the solution vector then contains a single column with directions of the unit bars and tension forces. Under consideration of the multiple inversions to reach convergence, the successive approximation arises as the more efficient algorithm due to a significantly decreased matrix size. In addition, the boundary condition is not directly included in a Newton-Raphson method such that its fulfillment is not guaranteed without sufficient convergence.

Validation of the mooring model
The validation of the proposed mooring model includes comparisons of line shape and maximum tension forces with experimental data. First, the maximum tension force is analyzed for a single mooring line attached to a buoy with prescribed heave motion. The mooring line consists of chains and wires of different material properties which can be found in the work of Chenga et al. 44 Here, linear material is assumed. A preliminary static test is conducted to validate the forces for different positive offsets. The resulting distribution is shown in Figure 6A. The deviation between the numerical model and experiments is generally small but increases slightly with increasing offset. The error is between −0.2% for offsets smaller than 5 m and increases to ≈1.1% for 30 m offset. Next, a simple harmonic heave motion with an amplitude of 15 m and a period of 10 s is prescribed to the top of the line. Good agreement with the measurements can be stated from comparing the distributions shown in Figure 6B. Generally, the numerical model slightly under predicts the troughs and over predicts the crests by ≈2.5% in peak height. A Fourier analysis of the time series indicates an error of less than 1% for the amplitude of the fundamental frequency. Next, the mooring model is validated by comparison to experimental results from wind channel measurements at the Institute of Ocean Engineering at the University of Rostock, Germany. The considered mooring system consists of a single rope with a length of L m = 1.82 m, a diameter of d t = 0.004 m, a specific weight of = 0.089 kg/m, and Young's modulus of 7.9 GPa. The line is fixed at the coordinates (x,y) = (0m, 0m) and mounted to a moveable load cell on the top to measure maximum tension forces. The wind machine at the inlet of the channel generates a laminar flow of predefined velocity which the mooring system is exposed to. A touch probe moves laterally along the deformed rope to record coordinates during static tests. Dynamic deformations are recorded using side view pictures.
In the first set of experiments, three different locations of the upper mounting point are investigated without inflow. Numerical results are produced using the same input data as given above. The results in Figure 7 show the distribution of the discrete line with 5, 10, and 50 bar elements and the comparison to the experimental data. Regardless of the considered distribution, the numerical model can accurately reproduce the experimental measurements. Small under and over predictions can be found for the coarsest line discretization and the hanging rope distribution in Figure 7A which is overcome by increasing the number of elements slightly. The second set of experiments is conducted by applying inflow velocities of v = 8 and v = 18 m/s on the stretched rope distribution of Figure 7C. The drag coefficients are calculated from (35), and c t = 0.35 is chosen from literature. 45 The results, shown in Figure 8, indicate that the numerical model is capable of representing the physical deformation of ropes in stationary flow conditions. All line discretizations coincide with the experiment accurately. Further information on the accuracy is provided by presenting the numerical and experimental results for the maximum tension forces over different inflow velocities in Table 1. It can be seen that the model converges toward the experimental data as the number of bar elements increases. This is demonstrated by calculating an extrapolated value without discretization error from the numerical results using Richardson extrapolation. 46 The remaining deviation from experiments is related to the model error. Here, a good agreement with the experiments can be stated because the maximum error is below 3% for v = 8 m/s.

Validation of the rigid body dynamics solver
The derived solver for the rigid body dynamics simulation is validated. A major advantage of the chosen Hamiltonian approach is the direct availability of the kinetic energy of the system T in runtime. In canonical form, it is calculated as (compare Reference 36): In absence of external forces, (36) is constant in time. Considering the time evolution of the phase space of fixed energy states, the volume of this space should remain constant according to Liouville's theorem. This condition is ensured by using time-reversible, symplectic methods. The second-order Verlet algorithm conserves energy in the mean 47 and has, hence, advantages over forward integration methods such as Runge-Kutta schemes. For system (16) and (19), the advance from time step n to n + 1 is calculated using Both, (37) and (39) are implicit equations which need time-consuming subiterations. Therefore, a fast forward method might be the better option in terms of efficiency. In the following, the motion of a rigid body 36  First, torque-free motion is assumed. The Verlet scheme benefits from its time-reversibility as can be seen in Table 2. It shows the error in energy conservation for different time step sizes after a 20 s simulation. The method conserves the original energy of the system up to machine precision irrespective of the chosen step size. In comparison, the Runge-Kutta method converges with the expected rate of p ≈ 4 and needs smaller time steps to conserve energy.
In a second example, a constant torque of 10 Nm around the x-axis is imposed. This implies a change of the kinetic energy as indicated in Table 3. Both schemes show the expected rate of convergence. Therefore, the Runge-Kutta scheme converges faster than the Verlet scheme which leads to a bigger discrepancy of the latter method to the converged energy of 528.520 J for small time steps. The time consumption of both methods is compared in Table 4. Here, the time consumed with the biggest time steps is taken as the reference. As expected, the Runge-Kutta scheme is not just faster than the Verlet scheme but also scales much better for smaller time-stepping due to its explicit nature. In combination with the results from Table 3, it can be concluded that the Runge-Kutta scheme is the more efficient choice for the calculation of the rigid body motions. In practice, the time steps are relatively small due to the velocity restriction in the fluid so that the lack of energy conservation is not seen as critical.

Validation of the floating algorithm
A validation case for the floating algorithm without mooring is presented. The numerical simulation of a free floating barge in waves is compared with the experimental data of Ren et al. 48 The two-dimensional (2D) setup is shown in Figure 9 and In a first step, the time series for the wave elevation at x = 5.5 m and the surge motion , the heave motion , and the pitch motion of the free floating barge are compared with the measurements for the time period between t/T = 6 and t/T = 12 ( Figure 10). It shows the results for the medium grid and CFL = 0.1. In Figure 10A, the wave elevation in front of the barge is presented. Good agreement with the experiments can be stated. Similarly, the heave motion in Figure 10B reproduces the measurements but shows minor under-and overshoots. The pitch motion is shown in Figure 10C. The numerically predicted frequency follows the experiment, but some deviations can be observed for the amplitudes of the motion. Both, experiments and the simulation show minor irregularities for the amplitudes at different time instances. Figure 10D presents the surge motion over time. A very good agreement of the numerical simulation with the experiments can be seen.
As a further step, the spatial and temporal convergence of the model is calculate based on the procedure described above. In case of an oscillatory converging behavior, the average of the two extreme values is taken as the extrapolated value. The presented error then corresponds to the discrepancy between this value and the experimental results. A constant CFL number of 0.1 is chosen for the spatial convergence test, and the medium grid is considered for the temporal convergence test. The mean frequency and amplitude are considered to be the variables of interest for this study. Special attention is given to the surge motion because a high-frequency motion superimposes the mean drift of the body due to Stokes drift forces. The latter presents the most important parameter in practice and is investigated in terms of its gradient in the space-time domain.
The spatial convergence study is presented in Tables 5 and 6. For the mean periods, all variables show oscillatory convergence. The resulting extrapolated values under predict the measurements by up to 4%. In comparison, the main amplitude tends to over predict the experiments by up to 10%. For the mean drift in the surge motion, the coarsest grid clearly underresolves the physics leading to a large over prediction. Thus, the extrapolated value and the resulting error is misleading. If the average value of the two finer grids for the surge motion is considered instead, the error would be 9% which is in a similar order as the predicted errors of other motions.
The temporal convergence study is shown in Tables 7 and 8. All variables converge except the heave period which shows an oscillatory diverging behavior. However, all the calculated periods are very close to the experiment and, therefore, the average value is still assumed to be valid. The general agreement of the numerical values with the measurements  is within 8%. Hence, mean periods tend to be predicted more accurately than mean amplitudes. It might also be noticed that the proposed model computes stable results even for higher CFL numbers.

Validation of the complete moored floating algorithm
The  The predicted and measured mean period and amplitude, as well as the discrepancy for the wave elevation and the three degrees of motion, are given in Tables 9 and 10. For the wave with T = 1.2 s, it is noticed that the numerical model generally under predicts both period and amplitude, for all motions. For the mean period, very good agreement with the experiments can be observed as the error band is below 4%. At the same time, the numerical amplitudes deviate from the measurements by up to 29% for the heave motion. As the period of the waves increases, improved prediction of the periods can be observed (below 2%), and the maximum amplitude error is 35% for the surge motion.
Both simulations indicate a reduction of accuracy in comparison to the free floating case. A source of error is possibly small measurement errors in one of the subsystems, which then adds up for this complicated fluid-structure interaction problem. In order to analyse this suspected sensitivity, reported parameters are varied slightly and compared with the original results. First, the location of the center of gravity is shifted in the positive and negative z-direction by 3 mm which might be linked to imperfect manufacturing of the barge. Tables 11 and 12 present the resulting errors. For T = 1.2 s, the chosen modification has no significant effect on the mean periods but considerably affects the simulated mean amplitudes. The error for the heave motion is reduced from 28% to 21%, and the mean wave amplitude prediction is improved by 50% if the center of gravity is slightly higher. By contrast, a lower center of gravity improves the pitch amplitude prediction from 26% to below 5% deviation from the experiments. For T = 1.6 s, considering a higher center of gravity has no significant effect on the computations. However, the decrease of the center improves the pitch amplitude which is also observed for the shorter wave. These findings indicate a strong sensitivity of the motion to the correct location of the center of gravity. At the same time, the strong coupling between the different motions is highlighted. An additional TA B L E 11 Error (%) between numerical and experimental mean period and amplitude for the simulation with a shifted center of gravity in z-direction with T = 1. investigation of the influence of the stiffness constant is shown in Table 13 for the shorter wave. The results indicate no significant difference if the value is increased by 5% but a higher difference is noticed if the stiffness constant is reduced by 5%.
The time series of the tension forces in the top of the front and back mooring line are shown in Figure 12 for both waves. As the experiment did not consider force measurements only the numerical results are provided. The alternating occurrence of force peaks can be observed and is explained through a strong back and forth motion of the barge as indicated in Figure 13. The relatively tense initial configuration of the mooring lines results in snap-like reaction forces as the barge is moving in and against wave direction in wave crest and trough situations. Generally, the forces increase with the increase of the wave period due to higher wave forces acting on the moored floating body.

APPLICATION TO A MOORED FLOATING SEMISUBMERSIBLE PLATFORM
The range of applications of the proposed moored floating FSI algorithm is broad and can include both, mooring system analysis and restricted motion analysis. Exemplarily, the analysis of a moored floating semisubmersible platform is presented. The geometry and properties of the platform are taken from Reference 49 and shown in Figure 14. The model consists of two rectangular floaters and six cylindrical columns placed symmetrically on the floaters. The height of the cylinders is such that they pierce the free surface at any time. The platform has a mass of m s = 34.58 kg, and the principal moments of inertia are (I xx , I yy , I zz ) = (3.227 , 5.23, 6.0246 kg m 2 ). The center of gravity is 0.086 m above keel, and the draft is 0.2 m. The platform is placed in the middle of a three-dimensional wave tank with a length of 10 m, a width of 2 m, a height of 3 m, and a water height of d = 1.9 m (see Figure 15). Active wave generation, 50 which prescribes the velocities and the free surface as Dirichlet boundary conditions at the inlet, is utilized. A numerical beach is placed at the outlet to prevent reflections. Six mooring lines are attached centered under the cylinders: two front and back lines, respectively, with the same length and two shorter middle lines. Their properties can be found in Table 14. Six different cases are considered: Three head wave conditions with H = 0.1 m and = 0.9735, 1.947, and 3.894 m, and three beam wave conditions with the same wave input. Both, free floating and moored floating configurations are investigated. In Figure 16A, the mean periods and amplitudes of the wave elevation 1.5 m in front of the platform are shown. It can be seen that the influence of the body on the wave periods is small, but the amplitude is influenced by the wave direction and length. As the wavelength decreases, the amplitude of the wave decreases. For head waves, the amplitude decreases if a mooring system is attached, whereas the amplitude increases in long beam waves due to stronger reflections. The same investigation is presented in Figure 16B for the heave motion. Without any restriction, the heave amplitudes and periods decrease with decreasing wavelength. When the platform is moored, strongly nonlinear interactions lead to the smallest periods for the intermediate wavelength, whereas the amplitude decreases with the length of the incident wave. As a tendency, the heave amplitudes are higher in head than in beam waves. For the pitch motion in Figure 16C, a general increase of the amplitude is observed if the mooring system is attached. In head waves, the pitch period increases with mooring. By contrast, the period decreases for beam waves. Overall, the most significant influence of the mooring system is on the pitch motion of the platform in head waves and on the heave motion in beam waves. An explanation for these observations is shown in Figure 17 which depicts the moored floating platform in head and beam wave crest situations for a wavelength of 1.947 m. In head waves ( Figure 17A), the two front lines get tensed, whereas the middle line is tensed in beam waves ( Figure 17A). The angle of attack at the mounting point is larger at the front lines than at the middle lines. This leads to larger horizontal forces in the front lines and larger vertical forces in the middle lines, and, accordingly, to higher influences of rolling and heave motions.

CONCLUSION
A new approach for simulating the effect of mooring systems on floating structures is presented. It is based on quasistatic assumptions and a discretization into elastic bars and mass points connecting them. The introduction of successive approximation for this problem raises the possibility to converge to a physically relevant solution within a small number of iterations. Any restriction on time stepping is removed and no initial shape of the line has to be given. It is also noticed that the method was originally developed for nets 20 and direct coupling of both system matrices is now possible for applications in aquaculture technology. An improved noniterative weakly coupled model has been proposed for the accurate and efficient simulation of fluid-structure interaction problems including the physical constraint of mooring systems. The basis of these developments is a previously presented FSI algorithm 11 which adapted a directional ghost-cell immersed boundary method for moving bodies. The simplicity of this algorithm, due to the absence of an iterative algorithm and the usage of a ghost-cell method, is now combined with several improvements to increase the overall accuracy. The rigid body dynamics solver is based on Euler parameters and Hamiltonian mechanics. Thus, the gimbal lock leading to a loss of a degree of freedom for Euler angles is avoided and arbitrary bodies and motions can be modeled. The resulting system is of the first order and free of constraints or Lagrangian multipliers. This enables the application of simple but high-order numerical integration techniques. A further advantage of a Hamiltonian is the energy-based formulation. In the validation process, this could be utilized to investigate the applicability of explicit methods like the Runge-Kutta method. It was shown that even relatively large time steps lead to sufficient conservation of energy such that conservative but less efficient methods do not have to be used. Furthermore, the computation of areas, moments of inertia and forces is improved by considering the original triangulated surface of the body. The integration process of the geometrical quantities is simplified by mapping and solving it analytically. Similarly, the introduction of linear interpolations and the explicit distinction between pressure and velocity grids leads to a more reliable formulation. The new FSI solver is validated for a challenging free floating case in regular waves. Convergence in time and space is shown, and all physically important aspects of the motion are captured by the proposed model.
The comparison of the physically constrained motion of a floating barge in waves with experiments underlines the accuracy of the coupled solver. Very good agreement exists for the mean periods. The complexity of this problem is highlighted by investigating the influence of small inaccuracies in the body's center of gravity position on the results. It is shown that the prediction of the mean amplitudes for heave and pitch is highly dependent on this position. In addition, the interaction of the different body motions and the mooring lines changes eventually. The application to a semisubmersible platform indicates the possibility of time-domain analyses including the coupled dynamics of waves, rigid bodies, and mooring systems.