DOI: 10.1002/fld.4894
R E S E A R C H A R T I C L E
Accurate modeling of the interaction of constrained floating structures and complex free surfaces using a new
quasistatic mooring model
Tobias Martin Arun Kamath Hans Bihs
Department of Civil and Environmental Engineering, Norwegian University of Science and Technology, Trondheim, Norway
Correspondence
Tobias Martin, Department of Civil and Environmental Engineering, Norwegian University of Science and Technology, Trondheim 7491, Norway.
Email: tobias.martin@ntnu.no
Summary
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 moor- ing 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 mod- eling 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 modi- fying 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.
K E Y W O R D S
CFD, fluid-structure interaction, immersed boundary method, mooring, numerical methods, rigid body dynamics
1 I N T RO D U CT I O N
Floating offshore and coastal constructions are exposed to a challenging environment with intense fluid-structure inter- action. Their safety is often related to fixing their position, which is usually achieved by mounting an appropriate mooring system. Therefore, the computer-aided design process requires not just the solution of a two-phase problem for the fluid and an accurate determination of the rigid body dynamics but also a mooring model.
This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.
© 2020 The Authors.International Journal for Numerical Methods in Fluidspublished by John Wiley & Sons, Ltd.
504 wileyonlinelibrary.com/journal/fld Int J Numer Meth Fluids. 2021;93:504–526.
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 dynam- ically. 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 models2 or dynamic overset grids.3,4 In overset methods, several mov- ing 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.5The 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.6They require just one Eulerian grid, and the body is described implicitly. Direct forcing methods7incorporate the fluid-structure interaction as an addi- tional 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.8Following the work of Calderer et al,9an easier way of respecting the bound- ary 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 Faltinsen10further sim- plified 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 Kamath11utilized this advantage to develop a local ghost-cell immersed boundary method for moving bod- ies. 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. David- son and Ringwood12present a review of a wide range of possible discretization methodologies. They are based on finite differences,13,14finite elements,15or spring elements.16Palm et al17 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.18If the purpose of the simu- lation 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.19This 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,20which 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.
2 N U M E R I C A L M O D E L S
2.1 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
𝜕ui
𝜕xi =0, (1)
𝜕ui
𝜕t +uj𝜕ui
𝜕xj
= −1 𝜌
𝜕p
𝜕xi
+ 𝜕
𝜕xj
(
(𝜈+𝜈t)⋅ (𝜕ui
𝜕xj
+ 𝜕uj
𝜕xi
))
+gi. (2)
Here,uiwithi=1,2,3 are the velocity components in the coordinate directions,𝜌is the fluid density,ppresents the pressure contribution,𝜈and𝜈tare the kinematic and turbulent viscosity, andgiis the gravity acceleration vector. The addi- tional viscosity is evaluated from thek-𝜔turbulence model with a modified source term to account for free surfaces.21The material properties of the two phases are determined for the whole domain in accordance to the continuum surface force model of Brackbill et al.22The 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.23The pressure is obtained as the solution of a Poisson equation from the continuity equation. The fully parallelized BiCGStab algorithm24 with PFMG preconditioner25 from the iterative linear system solvers of HYPRE are exploited. Temporal terms are dis- cretized applying the third-order accurate total variation diminishing Runge-Kutta scheme.26An 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.27The fifth-order accurate weighted essentially nonoscillatory (WENO) scheme of Jiang and Shu28is 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
𝜕𝜙
𝜕t +ui𝜕𝜙
𝜕xi =0. (3)
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.21Waves are gen- erated by applying the relaxation method to thex- andz-component of the velocities and the level set function𝜙. For a general variable𝜀it is given as Reference 32
𝜀(x) = Γ(x)𝜀analytical+ (1− Γ(x))𝜀computed, (4) with Γ the relaxation function taken from Reference 33. The same method is applied to damp potential reflec- tions near the outlet of the tank. Extensive validation of the described numerical wave tank can be found in Reference 21.
2.2 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 previ- ous algorithm,11the 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𝜌sholds. Given the approximation of the bodies’ surfaceSas
S=
⋃N i=1
Si, (5)
whereSiis assumed to be triangular, the properties are defined as 𝜌s∫V
∇⋅f(⃗x)dV=𝜌s∫S
(n⃗⋅k)f⃗ (⃗x)dS=𝜌s
∑N i=1
(n⃗i⋅k)⃗
∫Si
f(x)dS⃗ . (6)
Here,∇⋅f(⃗x) =1 returns the body mass,∇⋅f(x) =⃗ x,y,zpresent the center of mass coordinates in an inertial system and∇⋅f(⃗x) =x2,y2,z2,xy,yz,xzcompute the moments of inertia. Furthermore,n⃗iis the face normal vector of triangleSi, andk⃗is a unit vector in thex,y, orz-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
̈⃗
xs= F⃗⃗x
𝜌sV. (7)
The position⃗xsand velocityẋ⃗sof 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= (e0e1 e2e3)T which are related to the physically more relevant Tait-Bryan anglesΦ,Θ,Ψvia (cis cos andsis sin)35
e0=c(Φ 2
)⋅c(Θ 2
)⋅c(Ψ 2 )
+s(Φ 2
)⋅s(Θ 2
)⋅s(Ψ 2
), (8)
e1=s(Φ 2
)⋅c(Θ 2
)⋅c(Ψ 2 )
−c(Φ 2
)⋅s(Θ 2
)⋅s(Ψ 2
), (9)
e2=c(Φ 2
)⋅s(Θ 2
)⋅c(Ψ 2 )
+s(Φ 2
)⋅s(Θ 2
)⋅s(Ψ 2
), (10)
e3=c(Φ 2
)⋅c(Θ 2
)⋅s(Ψ 2 )
−s(Φ 2
)⋅s(Θ 2
)⋅c(Ψ 2
). (11)
In practice, this choice overcomes eventual issues due to the Gimbal lock effect. The back-transformation can be calculated using
Ψ =arctan 2(2⋅(e1⋅e2+e3⋅e0),1−2⋅(e2⋅e2+e3⋅e3)), (12)
Θ =arcsin(2⋅(e0⋅e2−e1⋅e3)), (13)
Φ =arctan 2(2⋅(e2⋅e3+e1⋅e0),1−2⋅(e1⋅e1+e2⋅e2)). (14) The four Euler parameters are depend as⃗eT⃗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
R=2
⎛⎜
⎜⎜
⎝
e20+e21−e22−e23
2 e1e2−e0e3 e0e2+e1e3
e0e3+e1e2 e20−e21+e22−e23
2 e2e3−e0e1 e1e3−e0e2 e0e1+e2e3
e20−e21−e22+e23 2
⎞⎟
⎟⎟
⎠
. (15)
Based on this, the kinematic equations for the rotational motion of the body in terms of the Euler parameters are given as
̇⃗
e= 1
2GT⃗𝜔, (16)
with ⃗𝜔the components of the angular velocity vector in the body-fixed coordinate system and
G=
⎛⎜
⎜⎜
⎝
−e1 e0 e3 −e2
−e2 −e3 e0 e1
−e3 e2 −e1 e0
⎞⎟
⎟⎟
⎠
. (17)
Introducing the momentum vectorh⃗=I⃗𝜔withIthe moment of inertia tensor, (16) can be rewritten as
̇⃗
e= 1
2GTI−1⃗h. (18)
Following the derivation of Shivarama and Fahrenthold,36a first-order ODE can be derived forḣ⃗using a system Hamil- tonian. Here, the constraint for the Euler parameters is fulfilled automatically. By setting the potential energy function to zero and assuming imposed momentsM⃗′⃗xin the body-fixed system, the equations read
ḣ⃗= −2GĠT⃗h+M⃗′⃗x. (19) 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 algorithm37is 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
⃗ Fx⃗=
∫Ω(−⃗np+𝜇⃗n⃗𝜏)dS(⃗x) =
∑N i=1
(−⃗np+𝜇⃗n⃗𝜏)i⋅ΔSi, (20)
withn⃗the surface normal vectors,𝜇the dynamic viscosity, ⃗𝜏the viscous stress tensor, andN the number of surface representing triangles. Furthermore,ΔSis 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,10to 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
𝜕p
𝜕xi||
||Γ= −1 𝜌⋅ Dui,Γ
Dt , i=1,2,3, (21)
(A) (B)
F I G U R E 1 Details of the directional ghost-cell interpolation method for one cell (blue): fluid points (white), interpolation points (black), ghost points (gray)
with ui,Γ 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 solu- tion 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 information38,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 extend- able to higher order interpolations but preliminary tests did not show significant improvements in accuracy or stability.
2.3 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 𝛾 𝜕2⃗r
𝜕t2 = 𝜕FT⃗f
𝜕s +F⃗e, (22)
with𝛾the specific weight of the material,FTthe magnitude of the tension force,⃗fthe unit vector pointing in the direction of this force, andF⃗eexternal forces including gravitation and hydrodynamic effects. Assuming small line motion in time and steady-state flow of the fluid, respectively, (22) simplifies to
𝜕FT⃗f
𝜕s = −F⃗e. (23)
F I G U R E 2 Discrete mooring lines: Inner knots (small black dots), Outer knots (big black dots), bars (red vectors)
In order to solve this force equilibrium, each mooring line is split intoNequally distanced bars of lengthltwith knots Pin between. As shown in Figure 2, the first and last knot,P(0)andP(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 F⃗G
(j)=𝛾⃗g⋅𝜌m−𝜌
𝜌m ⋅ l(j)t +l(j+1)t
2 , j=1,…,N−1, (24)
at any knotP(j). Here,𝜌mis the density of the mooring material and⃗gis an unit vector pointing in negativez-direction.
Hydrodynamic forcesF⃗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
F⃗H
(j)=l(j)t d(j)t 𝜌
2 ⋅[ct(⃗v⋅⃗f)|⃗v⋅f⃗|⋅⃗f +cn(⃗v− (⃗v⋅⃗f)⃗f)|⃗v− (⃗v⋅⃗f)⃗f|](j), j=1,…,N, (25) withcnandctthe 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 Casarella41derived a formula forcnas a function of the Reynolds number
cn(Re) =
⎧⎪
⎨⎪
⎩
8𝜋
sRe⋅(1.0−0.87s−2.0) if Re<1.0 1.45+8.55 Re−0.9 if 1.0≤Re<30.0 1.1+4.0 Re−0.5 else.
(26)
Forct, 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 knotP(j)(see also Figure 3):
⃗f(j+1)F(j+1)T −⃗f(j)FT(j)+F⃗(j)H +F⃗(j)G =0⃗, j=1,…,N−1. (27) A solution for the shape of the line and the distribution of tension forces can be found by gathering (27) into a suit- able 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), whereasF⃗Gis related to the length of the bars which is a function of tension forces. This elasticity of the material is respected by representingltas a polynomial ofFTsuch as
l(j)t,(k+1)= [
l(j)t ⋅ {
1.0+
∑P p=1
Kp(j)⋅
(FT(j)+FT(j−1) 2
)p}]
(k)
, j=1,…N, (28)
F I G U R E 3 Force equilibrium at inner knotP𝜈: Inner knot (white filled circle), bars (thin vectors), forces (thick vectors)
withKp the stiffness constants andkthe iteration index. Hooke’s law is considered in the applications shown below by only usingK1. 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,43a successive approximation can be found to ensure a converged solution. For this purpose, (27) is rearranged such that internal and external forces are separated
⃗f(j+1)F(j+1)T −⃗f(j)FT(j) = −(⃗F(j)H +F⃗(j)G), j=1,…,N−1. (29) Figure 2 indicates that the discrete line consists ofN bars but justN−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 pointsP(0)andP(N), and the physical coherence of the line during deformation:
∑N j=1
⃗f(j)l(j)t =⃗x(P(N)) −⃗x(P(0)). (30)
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))andx(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 Fbecause of the dependencies in (30). Under consideration of (29) and (30), the system reads
(T L )
⏟⏟⏟
A
⋅F=
( −(H+G)
⃗
x(P(N)) −⃗x(P(0)) )
⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟
B
, (31)
with the (N×N) matrixAcontaining the (N−1×N) submatrix of unknown tension forcesTand the (1×N) vector of the bar lengthsL. On the right-hand side,Bis filled with the (N−1×3) submatrices of hydrodynamic and static forces HandGand the vector connecting the two end points.
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 matrixG. The hydrodynamic force matrixHis 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 stepkis determined by Gaussian elimination with pivoting as
F(k)= (A(k))−1⋅B(k). (32)
F I G U R E 4 Illustration of the quasistatic mooring algorithm [Colour figure can be viewed at wileyonlinelibrary.com]
The lengths of bar unit vectors have to be equal to one by definition after being calculated from (32) maxj ||
||||
||||
||(⃗f(j))(k)||
||||
||2−1||
||=0, j=1,…,N. (33)
As this condition is generally not fulfilled, a correction step according to
(⃗f(j))(k∗)= (⃗f(j))(k)
|(⃗f(j))(k)|, j=1,…,N, (34) is performed. Consistency is given by multiplying the columns ofAby the Euclidean norm of the corresponding line of F(k). As the next iteration starts,F(k*)is inserted in (28), (24), and (25) to fill the system again. Convergence is formally proven in References 42,43 and typically reached within 100 iterations in the first time step and less than 50 afterward.
The algorithm stops after reaching the residual criterion maxj ||
||||
||||
||(⃗f(j))(k)||
||||
||2−1||
||<tol, j=1,…,N, (35)
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
F I G U R E 5 Illustration of the fluid-structure interaction algorithm for simulating moored floating structures in waves [Colour figure can be viewed at wileyonlinelibrary.com]
force and angle of attack at the upper mounting point. The force is added to the calculated fluid forces on the floating body (see also 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.
3 VA L I DAT I O N P RO C E S S
3.1 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 pre- scribed 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 numeri- cal 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.
(A) (B)
F I G U R E 6 Comparison of numerical and experimental results for the static analysis and the prescribed heave motion [Colour figure can be viewed at wileyonlinelibrary.com]
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 ofLm=1.82 m, a diameter ofdt=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 con- sidered 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 ofv=8 and v=18 m/s on the stretched rope distribution of Figure 7C. The drag coefficients are calculated from (35), andct=0.35 is chosen from literature.45The 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 exper- imental 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.46The remain- ing 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%forv=8 m/s.
3.2 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 systemTin runtime. In canonical form, it is calculated as (compare Reference 36):
T= 1
2(GTh)TGTI−1GGTh. (36)
F I G U R E 7 Numerical and experimental mooring line distributions without inflow [Colour figure can be viewed at wileyonlinelibrary.com]
(A) (B)
(C)
F I G U R E 8 Numerical and experimental mooring line distributions for different inflow velocities [Colour figure can be viewed at
wileyonlinelibrary.com] (A) (B)
v(m/s) Exp N=5 N=10 N=50 Extrapol. Error (%)
0.0 1.54 1.61 1.57 1.56 1.56 1.45
8.01 1.65 1.74 1.71 1.70 1.70 2.93
11.97 1.90 1.93 1.90 1.89 1.89 −0.43
14.08 2.01 2.06 2.03 2.02 2.02 1.22
16.08 2.21 2.21 2.19 2.18 2.18 −0.94
18.05 2.40 2.39 2.37 2.36 2.36 −1.64
Note: Values in (N) if not defined differently.
T A B L E 1 Numerical and experimental maximum tension forces for inflow velocities between v=0 andv=18m/s
Scheme 𝚫t=1 𝚫t=0.1 𝚫t=0.01 p
RK4 28.21 0.0043 5.57e-08 3.82
Verlet 1.99e-13 3.13e-13 3.13e-13 –
T A B L E 2 Convergence of the error in energy conservation in (J) using the fourth-order Runge-Kutta and the second-order Verlet scheme
Scheme 𝚫t=1 𝚫t=0.1 𝚫t=0.01 p
RK4 237.78 528.47 528.52 3.77
Verlet 574.25 529.37 528.56 1.74
T A B L E 3 Convergence of the total kinetic energy in (J) using the fourth-order Runge-Kutta and the second-order Verlet scheme
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 mean47and has, hence, advantages over forward integration methods such as Runge-Kutta schemes. For system (16) and (19), the advance from time stepnton+1 is calculated using
h∗=h(n)+Δt
2 L(e,h∗), (37)
e(n+1)=e(n)+Δt
2 L(e(n),h∗), (38)
h(n+1)=h(n)+Δt
2 L(e(n+1),h∗). (39)
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 body36with the principal moments of inertia of 400,307.808, and 200 kg m2and an initial momentum ofh⃗0= (346.410,0.0,200.0)TN ms shall be investigated with the second-order accurate Verlet and the fourth-order accurate Runge-Kutta scheme.
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 ofp≈4 and needs smaller time steps to conserve energy.
In a second example, a constant torque of 10 Nm around thex-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.
T A B L E 4 Time increase for different time step sizes using the
fourth-order Runge-Kutta and the second-order Verlet scheme Scheme 𝚫t=1 𝚫t=0.1 𝚫t=0.01
RK4 1 (0.021s) 1.372 4.26
Verlet 1 (0.070s) 5.67 51.72
F I G U R E 9 Numerical setup for the simulation of a two-dimensional barge in a wave tank [Colour figure can be viewed at wileyonlinelibrary.com]
3.3 Validation of the floating algorithm
A validation case for the floating algorithm without mooring is presented. The numerical simulation of a free float- ing 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 consists of a wave tank with a length of 20 m, a height of 0.8 m, and a water depth ofd=0.4 m. A 0.3 m×0.2 m barge with density𝜌s=500 kg/m3is placed in the tank at (x,z)=(7.0 m, 0.4 m). Regular waves are modeled as second-order Stokes waves in a wave generation zone at the inlet. The wave generation zone is one wavelength long.
The incoming waves have a height ofH=0.04 m, a period ofT=1.2 s, and a wavelength of𝜆=1.936 m. Wave reflec- tions are prevented by placing a numerical beach of two wavelengths at the outlet. The numerical methods described in Section 2.1 are applied. The convergence is investigated using three different mesh configurations with the smallest cell sizesΔx=0.012,0.01,0.007 m, leading to 80,120, and 200 thousand cells. Stretching functions are used to coarsen the cell size inx-direction toward the inlet and outlet andz-direction toward bottom and top.
In a first step, the time series for the wave elevation𝜂atx=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 betweent/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 con- stant 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, there- fore, the average value is still assumed to be valid. The general agreement of the numerical values with the measurements
(A)
(B)
(C)
(D)
F I G U R E 10 3DOF motion of the two-dimensional barge over time. Comparison of numerical and experimental results forΔx=0.01 m and CFL=0.1 [Colour figure can be viewed at wileyonlinelibrary.com]
T A B L E 5 Spatial convergence of the numerical mean period (dimensionless) in comparison to the experimental results
Motion Coarse Medium Fine Extrapol. Exp Error (%)
𝜂 0.951 1.046 0.997 0.974 0.998 2.4
𝜁 0.992 1.012 0.995 0.994 0.997 0.3
𝜃 1.018 1.051 1.027 1.023 1.059 3.5
Note: For the surge motion, the high-frequency component is analyzed here.
T A B L E 6 Spatial convergence of the numerical mean amplitude
(dimensionless) in comparison to the experimental results
Motion Coarse Medium Fine Extrapol. Exp Error (%)
𝜂 0.0567 0.0558 0.0546 0.0531 0.0513 −3.4
𝜁 0.0580 0.0593 0.0587 0.0584 0.0615 5.3
𝜃 0.0437 0.0394 0.0426 0.0432 0.0388 −10.2
𝜉 0.1008 0.0763 0.0812 0.0911 0.0715 −21.5
Note: For the surge motion, the mean drift is considered.
T A B L E 7 Temporal convergence of the numerical mean period (dimensionless) in comparison to the experimental results
Motion CFL 0.5 CFL 0.3 CFL 0.1 Extrapol. Exp Error (%)
𝜂 1.021 1.036 1.046 1.048 0.998 −4.7
𝜁 1.001 0.999 1.012 1.004 0.997 −0.7
𝜃 1.027 1.047 1.051 1.051 1.059 0.8
T A B L E 8 Temporal convergence of the numerical mean amplitude
(dimensionless) in comparison to the experimental results
Motion CFL 0.5 CFL 0.3 CFL 0.1 Extrapol. Exp Error (%)
𝜂 0.0587 0.0571 0.0558 0.0554 0.0513 −7.4
𝜁 0.0621 0.0604 0.0593 0.0591 0.0615 4.1
𝜃 0.0426 0.0409 0.0394 0.0387 0.0388 0.3
𝜉 0.0774 0.0683 0.0763 0.0769 0.0715 −7.0
Note: For the surge motion, the mean drift is considered.
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.
3.4 Validation of the complete moored floating algorithm
The three degrees of motion of a moored floating barge in waves is investigated. Experiments were conducted at the wave flume of the Ludwig-Franzius-Institute Hannover, Germany. The investigated barge Is made of a material with a uniform density of 680 kg/m3. It has a length of 0.3 m and a height of 0.15 m. The incoming waves have a height of H=0.03 m and periods ofT=1.2 s andT=1.6 s, respectively. Both waves are modeled using the second-order Stokes theory. Figure 11 illustrates the 2D numerical domain. The tank is 20 m long, and the water height is fixed tod=0.85 m.
The center of the barge is initially located at (x,z)=(10 m,0.823 m). Based on the results of the previous simulations, a cell size of 0.012 m and CFL=0.5 are considered. Two mooring lines are fixed to the barge at still water level. The distance between the center of the barge and the bottom mounting point of the lines is 4.15 m along thex-axis. Both lines have a length ofLm=4.07 m, a diameter ofdt=0.94 mm, and a stiffness ofK1=10e5 N−1. Their weight is neg- ligible in water. A wave gauge is placed 2 m in front of the center of the barge to calculate the incident and reflected waves.
F I G U R E 11 Numerical setup for the simulation of a two-dimensional moored floating barge in a wave tank [Colour figure can be viewed at wileyonlinelibrary.com]
T A B L E 9 Numerical and experimental results and error (%) between numerical and experimental mean period and amplitude for the original simulation of a moored floating barge in waves withT=1.2 s
Exp Num Error
Motion Period (s) Amplitude (m) Period (s) Amplitude (m) Period Amplitude
𝜂 1.22 0.015 1.20 0.013 1.505 14.012
𝜁 1.23 0.018 1.19 0.013 3.363 28.978
𝜃 1.21 8.857 1.18 5.605 2.154 26.717
𝜉 1.22 0.016 1.19 0.012 1.805 23.136
T A B L E 10 Numerical and experimental results and error (%) between numerical and experimental mean period and amplitude for the original simulation of a moored floating barge in waves withT=1.6 s
Exp Num Error
Motion Period (s) Amplitude (m) Period (s) Amplitude (m) Period Amplitude
𝜂 1.61 0.016 1.59 0.014 1.155 12.817
𝜁 1.62 0.016 1.59 0.014 1.798 11.524
𝜃 1.62 2.158 1.59 3.379 1.757 −26.595
𝜉 1.61 0.029 1.60 0.0185 0.769 35.903
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 withT=1.2 s, it is noticed that the numer- ical 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 negativez-direction by 3 mm which might be linked to imperfect manufacturing of the barge. Tables 11 and 12 present the resulting errors. ForT=1.2 s, the chosen modification has no significant effect on the mean periods but considerably affects the simulated mean ampli- tudes. 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 pre- diction from 26%to below 5%deviation from the experiments. ForT=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
T A B L E 11 Error (%) between numerical and experimental mean period and amplitude for the simulation with a shifted center of gravity inz-direction withT=1.2 s
CG+3 mm CG−3 mm
Motion Period Amplitude Period Amplitude
𝜂 1.591 7.151 1.941 32.452
𝜁 3.268 21.095 3.750 27.568
𝜃 1.164 41.296 2.894 4.665
𝜉 3.231 28.260 2.988 31.751
T A B L E 12 Error (%) between numerical and experimental mean period and amplitude for the simulation with a shifted center of gravity inz-direction withT=1.6 s
CG+3 mm CG−3 mm
Motion Period Amplitude Period Amplitude
𝜂 1.185 14.220 0.839 12.919
𝜁 1.692 12.670 1.849 10.060
𝜃 9.537 −30.673 0.789 −14.231
𝜉 1.571 35.449 1.097 41.024
T A B L E 13 Error (%) between numerical and experimental mean period and amplitude with varying stiffness constants withT=1.2 s
EA+5% EA−5%
Motion Period Amplitude Period Amplitude
𝜂 1.697 16.854 1.662 29.674
𝜁 3.288 27.910 3.571 30.415
𝜃 1.371 26.889 1.146 39.142
𝜉 3.393 25.794 2.963 27.619
F I G U R E 12 Tension forces in the top of the front and back mooring line over time. Forces forT=1.2 s are marked in red; forces for T=1.6 s are marked in blue. Solid lines indicate the front lines; dashed lines indicate the lines behind the body [Colour figure can be viewed at wileyonlinelibrary.com]
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 indi- cated 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.
