Finite element modeling of pipe-laying dynamics using corotational elements

Abstract A three-dimensional finite element model is built to compute the motions of a pipe that is being laid on the seabed. Corotational beam elements account for geometric nonlinearity. The pipe is subject to contact, hydrodynamic forces, gravity, and buoyancy. New in this article is the addition of nodal moments due to buoyancy and nodal correctional forces to compensate for a cross-sectional area mismatch. The results show a modest increase in accuracy due to these moments and a significant increase due to the correctional forces.


Introduction
Numerical simulations of offshore pipe-laying are typically performed using a finite element approach. Commercial software is available for this type of simulations, such as Offpipe [1], Pipe-Lay [2], and OrcaFlex [3]. Each of these software packages has its advantages as well as limitations.
In pipe-lay analyses, large rotations of the elements occur, as an initially horizontal pipe rotates to an almost vertical position. Several methods for large rotations can be found in literature. Reissner [4] and Simo and Vu-Quoc [5,6] developed the geometrically exact beam method. This method is called exact because no assumptions have been made on the size of the deformations.
An alternative method for large rotations is the absolute nodal coordinate formulation. This method was presented by Shabana et al. [7]. An advantage of this method is that it results in a constant mass matrix. For explicit time integration, this is a big advantage. But when a time integration method from the Newmark family [8] is used, the advantage that the mass matrix is not rotated is of little interest and is offset by the disadvantage that the derivation of the stiffness matrix is more complex. Furthermore, Romero [9] and Bauchau et al. [10] compared the absolute nodal coordinate and geometrically exact formulations. They concluded neither method to be superior and that the selection of method was application dependent.
In another paper, Shabana and Wehage [11] suggested a substructuring technique with a floating frame of reference. The deformations of the substructures are described by elastic vibration modes. This method, which is often used in multibody dynamics, is effective when a substructure rotates as a whole. However, in the pipe-laying simulation, such a substructure cannot be defined.
The chosen method to account for large rotations in the current numerical model is the corotational method, which has been well established in literature [12]- [15]. The formulation used in the current model is the method of Crisfield [16,17], which is explained in more detail in the next section.
Compared to the absolute nodal coordinate formulation [7], an advantage of the corotational method is that it has a smaller number of degrees of freedom [9]. In the corotational method, a lumped mass matrix can be chosen to prevent the necessity of rotating the mass matrix. However, the reduction in calculation time by not having to rotate the mass matrix is negligible. Furthermore, a comparison between the corotational method, a recursive method, and a substructuring technique with a floating frame of reference by Disveld [18] showed the corotational method to be the most efficient.
The pipe is subject to gravity. When submerged, the pipe is also subject to buoyancy forces and can be subject to water current. Most commercial software programs model buoyancy using only Archimedes law, while disregarding effects of the pipe's curvature. Then a postprocessing step is required where a correction is made to the axial force, known as effective tension [19]. In detailed analyses, such as the analysis of torsional rotation due to residual curvature [20,21], the correct axial force is required. Therefore, it is necessary to calculate the axial forces during the simulation. In the current article, buoyancy is modelled using the method of Yazdchi [22], which avoids this postprocessing step. New in this article is the addition of equivalent moments due to buoyancy, which are derived from a distributed pressure. Furthermore, a correctional nodal force is introduced, which is related to an error due to mismatch in cross-sectional area.
In finite element analyses, mesh refinement is often applied locally, resulting in a mesh with elements of unequal lengths. The novel contributions of this article to the buoyancy calculation are in particular important when unequal element lengths are used.
The following section presents the three-dimensional corotational beam element. The three sections thereafter elaborate on external loads. First, the external distributed loads are considered, then buoyancy loads and finally hydrodynamic forces. The explanation of the model finishes with a section on waves and vessel response and a section on contact and friction. Both static and dynamic examples will be presented in the results section. Preliminary results of this model have been presented in [23], which included static results of buoyancy, validation of the dynamic model, results of numerical damping, and results of hydrodynamic forces.

Three-dimensional corotational beam element
The corotational beam element has two nodes. Each node has six degrees of freedom: three translations and three rotations. The main idea of the method is that it uses a local coordinate system to calculate forces and moments based on a small strains assumption. For each element, the small strain results calculated in this local coordinate system are rotated to the global coordinate system. Because the corotational method is geometrically nonlinear, an iterative solution procedure is needed. It is chosen to use the Newton-Raphson method.

Local coordinate system
In a corotational formulation, a local "corotational" coordinate frame is introduced. The origin of this frame is placed on the first node of the element and the local coordinate frame is placed such that the local x l -axis runs through the other node. As a consequence, the displacements in local y l -and z l -direction are always zero at both nodes. This placement of the coordinate system results in only seven local degrees of freedom.
Here w i ; h i ; and / i are the rotations about the local x l -; y l -; and z l -axis, respectively. Subscript i refers to the first or the second node of the element and u l is the axial elongation of the element.
An important factor when rotating in three dimensions is that the rotations are not commutative, which means that a rotation cannot be added to another rotation. Therefore, the three-dimensional rotations are implemented using nodal triads. Detailed derivations can be found in [17]. For each element, two nodal triads, T and U, and one element triad, E, are defined. This can be seen in Figure 1.
Each triad contains three orthogonal unit vectors. The unit vectors of the three triads are named t i ; e i ; and u i : These unit vectors define the local rotations of the element nodes. The middle triad defines the element coordinate system. As shown by Crisfield [17], these unit vectors can be used to find the local rotations.
After each iteration, matrices T and U are updated as follows: Here Da and Db are the iterative spin variables. These variables are the increments of the rotations at the end of an iteration. They are nonadditive and represent the vector about which the nodal system has rotated. The length of the vector represents the angle of rotation. The rotation matrices DT and DU are calculated from the spin variables using Rodrigues rotation matrix [24]. The element triad E is updated using the mean rotation from the first to the second node, R ¼ r 1 ; r 2 ; r 3 ½ ; as explained in [17].

Internal forces
The three-dimensional stiffness matrix can be found in many textbooks, like [25], where it is presented as a 12 Â 12 matrix. The size of this matrix can be reduced when using the corotational method. Since the number of local degrees of freedom is reduced to seven, the local stiffness matrix is reduced to 7 Â 7: Using a linear elastic material model, the stiffness matrix multiplied by the local displacements vector yields the local element forces, F l ¼ K l d l . The forces in the global coordinate system are found by rotating the local forces using rotation matrix B: The rotation matrix is derived by taking the variation of the local variables with respect to the global degrees of freedom d: This derivation can be found in [17]. The resulting rotation matrix is: This is a 7 Â 12 matrix. Heret i is the skew symmetric cross-product matrix of vector t i : Furthermore, Q is the following 3 Â 12 matrix. and Here I is a 3 Â 3 identity matrix.

Consistent tangent stiffness
The tangential stiffness matrix is found by taking the variation of the forces in the global coordinate system.
The result is the rotated local stiffness matrix plus the geometric stiffness matrix. method used is a general method, which can also be used for other distributed loads such as water current.

Equivalent forces and moments
The forces and moments at the nodes are derived such that they are kinematically equivalent to a uniformly distributed load. The global equivalent forces at the first node (see Figure 1) are a function of the distributed force vector q; which contains the distributed loads in the global coordinate system.
The global forces at the second node are equal to the global forces at the first node. The moments at the first node are first calculated in the local coordinate system. The moment around the local x-axis is zero since the external forces do not contribute to torsional moments. This is exact for an undeformed element and an approximation for a deformed element. The moment around the local y-axis is due to all forces in e 3 -direction. The minus sign is due to the right-hand rule. The moment around the local z-axis is a result of the forces in e 2 -direction.
These moments can be rotated to the global coordinate system using nodal triad E from (2).
The moments at the second node are minus the moments at the first node. Now the element force vector due to distributed forces becomes as follows:

Consistent tangent stiffness of external distributed loads
Since these forces are a function of the displacements, a tangent stiffness matrix needs to be derived in order to maintain quadratic convergence when the Newton-Raphson method is used. This matrix is derived by taking the variation of F q with respect to the global degrees of freedom. The variation of the forces at each node is zero since dq ¼ 0: Therefore, the first three rows and rows 7-9 of the global stiffness matrix are zero. The variation of M q is calculated as follows: This yields the following contribution to the stiffness matrix:  Note that 0 nxm represents a n by m zero matrix. The matrix K q is of size 12 Â 12. Furthermore A is from (11), Q can be found in (9) and Y is as follows:

Buoyancy
The treatment of buoyancy forces is based on Yazdchi [22]. All elements are considered straight but connected to the adjacent elements at an angle. By integrating the pressure analytically over the outer surface of the pipe, Yazdchi finds three contributions to the buoyancy forces. The first contribution is due to the distributed pressure on the straight pipe, which is the contribution where this article adds the novel buoyancy moment. The second contribution is determined by the pipe's curvature and the third contribution, as described by Yazdchi, is due to capped ends. The novel correction on the cross-sectional area is presented as a fourth contribution.

Distributed pressure
The forces due to distributed pressure on a straight element can be calculated as follows. On a horizontal pipe section, the distributed force in N/m is: Here A o is the outer section area of the pipe, A i is the inner section area of the pipe, q w is the density of the sea water, q i is the density of the fluid in the pipe, and g is the standard acceleration due to gravity.
This equation can be found by integrating the pressure over the surface area of the pipe and the result is determined by the pressure difference in the water. However, if the pipe is not horizontal, the pressure difference between top and bottom of the pipe is smaller. Therefore, the factor e T n k is added: Here k is a unit vector in the global z-direction (see Figure 1). e n is a unit vector, normal to the local x l -axis and in the plane spanned by the local x l -axis and global z-direction. The derivation of e n can be found in [22].
If an element is only partially submerged, equivalent forces are calculated for both nodes. This is done by introducing factors b 1 and b 2 for the first and second node, respectively. If the pipe is completely submerged, b 1 and b 2 are 1=2: If only the first node is submerged, the parameters b i are determined from the position of the first and second node with respect to the water line, indicated by h 1 and h 2 ; which are positive when the node is submerged and negative when the node is above water.
Here h is calculated by Up to here, the buoyancy is identical to Yazdchi [22]. New in this article is the addition of buoyancy moments, which are determined using factors b m1 and b m2 : If the element is completely submerged, these terms are equal to 1=12 and À 1=12; respectively. This results in a nodal force vector that is statically and kinematically equivalent to the uniformly distributed load, comparable to the nodal loads of the general distributed loads.
If only the first node is submerged, the moment at the first node can be calculated using Here a is the distance between the first node and the water surface. The local displacement w l within the element is determined by a cubic Hermite interpolation function. By calculating this integral, factor b m1 is found. The same can be done for the moment at the second node to find b m2 : The equivalent nodal moments can now be found. With these moments, the buoyancy force vector due to distributed pressure becomes: Here e n is the direction of the equivalent forces and e 1 Â e n ð Þ is the direction of the equivalent moments.

Curvature effects
The forces due to curvature effects are implemented as given by Yazdchi [22]. They are derived from the change in surface area as a result form the connection to the adjacent element.
The forces at the second node are the negative of the forces at the first node, because of the direction of rotation at the second node. Note that the force vector F b 2 is dependent on the water depth.

Capped ends
Third, Yazdchi also gives a buoyancy term due to capped ends. This term is only applied to the end of the pipe. The following equations show the capped end term on the first and second node, respectively.

Buoyancy area mismatch
In the buoyancy calculation, the elements are assumed straight. At the nodes, the elements are connected at a certain angle, as illustrated in Figure 2. a and b indicate the local angles of rotation at the intersecting node, calculated using (3). Yazdchi [22] assumes that the cross-sectional area of both elements is equal at their common node; thus a ¼ b: However, this is not always true, for example, when two adjacent elements have different lengths. Therefore, a correctional force is derived, based on the difference between the crosssectional areas.
The total cross-sectional area when local rotation a ¼ 0 is the surface area of a circle. The total area when a 6 ¼ 0 is equal to: If the local rotation of the adjacent element b is unequal to the local rotation of the current element a; there is a difference in area equal to: In the calculation of the forces of the current element, the local rotation of the adjacent element is unknown. Therefore, the difference with respect to b ¼ 0 is calculated.
The force acting on this area is equal to the pressure multiplied by the area. To expand to three dimensions, cosa is substituted with t T 1 e 1 : This results in the following force vector for the area mismatch of one element: The minus one terms in this force vector could be removed, as it is cancelled out to the term of the adjacent element when constructing the system force vector. If it is left out, the capped end term needs to be removed to avoid calculating the capped end area twice. Here it is chosen to keep the minus one term, in order to illustrate that the correctional force is zero for a straight pipe.

Hydrodynamic forces
Hydrodynamic forces are calculated using Morison's equation [26]. The force on a moving pipe in moving water is given by: Here q w is the density of the water, C d is the drag coefficient, C a is the added mass coefficient, D is the pipe diameter, and _ d rn and € d rn are the relative velocity and relative acceleration in normal direction of the pipe. The last term of (34) is called the Froude-Krylov force, which is dependent on the total acceleration of the water normal to the pipe, due to the pressure gradient in the water.
The relative velocity is calculated from the difference between the water velocity and the pipe velocity. The same can be done for the relative acceleration. The pipe velocity _ d is the velocity in global x-, y-and z-directions.
The relative velocity in the normal direction of the pipe is calculated by subtracting the tangential relative velocity from the relative velocity. Again, the same can be done for the relative normal acceleration.

Waves and vessel response
In this section, the vessel's motions are described, which are dependent on the encountering waves. First regular and irregular waves are described followed by the response of the vessel to the wave motion.

Regular wave
A standard method of parameterizing a wave g is with its height h in meters and frequency x in rad/s.
Here k is the wavenumber and t is time. The wave is propagating in positive x direction, that is determined by the encountering wave angle u; which is the counterclockwise angle between the global x-axis and the encountering wave direction. This is illustrated in Figure 3. At the center of motion x ¼ 0; at all other positions on the water surface, x is: Airy wave theory is used to describe the motion of the seawater. In this theory, a velocity potential is defined. The seawater is assumed incompressible, inviscid, and irrotational. Furthermore, the seawater is assumed to be deep. A water is said to be deep if the water depth is larger than half the wavelength (see, e.g., Dean and Dalrymple [27] or Faltinsen [28]). With these assumptions, the velocity potential associated with (37) is: Water velocity and acceleration can easily be determined from this velocity potential [27,28].

Irregular wave
In reality, a wave is not a perfect sine function. For irregular waves, represented by a collection of regular waves, a wave spectrum is introduced Different wave spectra can be found in literature, such as the JONSWAP (JOint North Sea WAve Project) spectrum [29], the Pierson-Moskowitz spectrum [30], and the Bretschneider spectrum [31]. Implemented in the current model is a Pierson-Moskowitz spectrum, representing a fully developed sea. It is based on a spectrum density S; which is a function of a significant wave height H s and a mean zero-up-crossing period T z [32]. The wave surface profile is determined by a summation of regular waves.
The term i is the random phase angle, which is between 0 and 2p: Figure 3. The global coordinate system related to pipe-lay vessel. The wave angle u is zero when encountering the ship at the stinger.

Vessel response
The vessel's motion is defined as a response to the wave surface profile at the center of motion. This is done by using the response amplitude operators, or RAO's. At a certain frequency, each degree of freedom has two RAO parameters: K for the amplitude and f for the phase. For a regular wave the response n is calculated as follows: For an irregular wave the contribution of each frequency is summed:

Contact and friction
The part of the pipe in contact with the seabed experiences reaction forces from it. Furthermore, friction forces are introduced when the pipe moves while in contact with the seabed. Contact and friction as implemented in the model are explained in this section. The Penalty method is used to describe contact. To start with, the equation of motion is extended with a contact term, N d ð Þ : Here M is the mass system matrix, F int and F ext are the internal and external force system vectors, and d; _ d; and € d are the pipe displacement, velocity, and acceleration system vectors. For each node, the penalty force is equal to the following: Here p > 0 is the penalty factor. The gap C z is the distance between the pipe and the seabed in the direction normal to the seabed. If there is no contact, the gap is larger than zero and the contact force remains zero. If there is contact, the force is equal to the penalty factor multiplied by the gap. This force is in the direction normal to the seabed and it can be seen as a spring embedded in the seabed at the point where contact is enforced.
Here n is the direction normal to the seabed. Within a single iteration, the gap is constant, meaning that a node cannot switch from being in contact to not being in contact or vice versa. The contribution to the stiffness matrix is derived to be: Note that all stiffness terms associated with rotations are zero.
If the pipe is in contact with the seabed, a friction force is present, which is represented by the Coulomb friction model. Since the first end of the pipe is fixed, the pipe cannot move in axial direction. Therefore, only friction forces perpendicular to the pipe axis are modeled. First, the slip force is introduced: Here l is the friction coefficient and F n ¼ minðpC z ; 0Þ is the magnitude of the normal force. The magnitude of the friction force can now be calculated as follows: Here C t is the "gap" tangential to the seabed, The friction force is in direction t f : The friction force therefore becomes F f ¼ t f F f : The stiffness matrix is found to be the following:

Results
The numerical model described in this article is implemented in Cþþ using the FeaTure framework for finite element codes [33]. This section shows static and dynamic results, including a comparison with and without the suggested buoyancy moment, an example that illustrates the importance of the novel buoyancy term, a validation with the industry standard numerical simulation software Offpipe [1] and a J-lay pipelaying simulation. All dynamic results are obtained using the Hilber-Hughes-Taylor time integration method, or HHT-a method [34].

Buoyancy with and without additional moment
It is expected that the additional buoyancy moment increases the accuracy of the results, which is investigated by comparing the results with and without the additional moment to a benchmark. The benchmark uses a refined mesh of 100 elements and includes the additional buoyancy moment. If the benchmark is created excluding the buoyancy moment, the results of the comparison are the same, because the difference between the benchmark with and without the buoyancy moment is very small. This small difference can be explained by the fact that the buoyancy moments are very small for small elements, due to the L 2 term in (27). Moreover, when the angle between two adjacent elements of the same length is zero, the summation of moments on their common node is zero. Thus, when using 100 equally sized elements where the angle between two adjacent elements is small, the summation of moments at their common node is almost zero.
In the example used for the comparison, a massless pipe with a length of 100 meters is placed horizontally 100 meters below the water surface. The pipe is inclined at the first pipe end, whereas the second pipe end is free. Only buoyancy forces are applied on the pipe. Pipe and water properties are presented in Table  1. The solution is found in two steps. In the first step, 20% of the total buoyancy force is applied to the pipe. In the second step, the full buoyancy force is applied. The static equilibrium of both steps is found using the Newton-Raphson method.
In Figure 4, the result of the example is depicted. The blue area is the water surface and the red line is the centerline of the pipe. Figure 5 shows the convergence of both steps, with iterations numbers on the horizontal axis and the displacement-based error on the vertical axis. The displacement-based error is defined as the ratio between the l 2 -norm of the solution update vector of the Newton-Raphson iteration and the l 2 -norm of solution vector: The result in Figure 4 is as expected, the buoyancy forces push the pipe towards the water surface, and there is no displacement in y-direction. The maximum displacement in z-direction is w ¼ 84:54 m: The convergence of both steps, as depicted in Figure 5, is quick due to the Newton-Raphson method.
Next, the example is recalculated with 10 equally sized elements, once with the additional moment and once excluding the additional moment. The same is done using 10 unequally sized elements. These results are compared to the benchmark. At each node that the result and the benchmark have in common, the error of the result is calculated using the vector u; which consists of the global displacements in x-, y-, and z-directions.    In Figure 6, this error is depicted for the calculations described above. The logarithm of the error is on the vertical axis and the s-coordinate on the horizontal axis, which is the coordinate along the pipe's axis: the first pipe end is at s ¼ 0 m and the second pipe end is at s ¼ 100 m: It can be observed that in the neighborhood of s ¼ 0 m; where the degrees of freedom are fixed, the error with and without additional buoyancy moments is equal. At the free pipe end, where s ¼ 100 m; the difference in error with and without the additional buoyancy moment is significant. Furthermore, these results show that when using unequally sized elements, the results can become more accurate without increasing the number of elements. At the free pipe end, this increase in accuracy is not obtained, when the additional buoyancy moment is neglected.
The improvement in accuracy seems only present at the free end of the pipe. At the left end of the pipe, the rotations are fixed, and thus the moment is not applied. Furthermore, if the elements are of equal length and the rotation between two adjacent elements is small, the nodal moments are very small. This is the case close to the fixed pipe end, which explains the equal accuracy in the neighborhood of this end. The node at the free end of the pipe has the largest moment applied to it, because it is only connected to one element. This causes the largest difference in accuracy at the free pipe end.

The necessity of the new buoyancy term due to area mismatch
This example shows the importance of the new buoyancy term by looking at the summation of forces in global x-and z-direction. For this, a pipe is placed 2000m below the water surface. Its left end is fixed, and the right end is rotated to an angle of p=2 by a prescribed rotation about the y-axis. Pipe and water properties are presented in Table 1. Three elements are used, two of which have length L ¼ 25 m and one has length L ¼ 50 m: The node that connects two unequally sized elements has different local rotations for both elements, as can be seen in Figure 7.
The summation of buoyancy forces in the global xand z-directions are shown in Table 2. The expected result from the summation of buoyancy forces in xdirection is zero. In z-direction, the expected result is Archimedes force, which is 1=4pD 2 qgL ¼ 82:8 kN: From the results in the table, it can be concluded that the forces calculated without F b 4 are incorrect. The results calculated including F b 4 are very accurate.
If the example is repeated with elements that have equal sizes, the results without F b 4 are still inaccurate, as can be seen in Table 3. This is because the capped end forces are calculated with respect to the cross-section of a straight pipe. Using F b 4 increases the capped end force, such that the results are accurate.

Contact
In this static example, the pipe is partially in contact with the seabed. The solution of this example will be the starting point for the dynamic simulation of pipelaying with the J-lay method. Some properties of the pipe and the environment are presented in Table 4.
Before the first static step, the pipe starts as a straight and horizontal pipe from position x ¼ À 200 to x ¼ 0: The z coordinate of the pipe is À 50; equal to the depth of the water. This can be seen in Figure  8. The blue area at z ¼ 0 m is the water surface.
In the first static step, gravity, buoyancy, and contact forces are added to the pipe. Because the contact is at the outer surface of the pipe while the pipe coordinates are at the center of the pipe, the result of this first step is a displacement of all nodes in positive zdirection of approximately 0:16 m: The displacement is slightly smaller than the radius of the pipe due to the gap in the penalty method for contact.   The second end of the pipe is the end which is connected to the pipe-laying vessel. In the second static step, this end of the pipe is moved by a prescribed displacement of 60 m: The result is that, at the end of the second step, the second end of the pipe is positioned 10 m above the water surface. Since there is no friction in x-direction, the pipe displaces in x-direction such that all static forces are in balance. The result can be seen in Figure 9.
In preparation of the dynamic calculation, the second pipe end has to be moved to the position of the vessel at time t ¼ 0: This is done in the third and last step, which can be seen in Figure 10. The position of the vessel at time t ¼ 0 can be different for each calculation, due to the random phase angle in (40).
Due to the highly nonlinear on/off behavior of this contact method, several substeps are required to converge. In this example 5, substeps are used. The convergence, with respect to the error in (50), is depicted in Figure 11. In all steps, quadratic convergence can be recognized in the neighborhood of the static balance. It can be seen that the first step, where all distributed forces were added, converges in just a few iterations. Also, the last step converges quickly. The 5 (sub)steps in-between, where the second end is moved in z-direction, need a higher number of iterations. This is due to the highly nonlinear contact behavior. The number of iterations per substep can be reduced by increasing the number of substeps. However, this does not result in a reduction of the total number of iterations.

Validation
The results of the static contact example in Figure 9 are compared to results from industry standard numerical simulation software Offpipe [1]. Figure 12 shows the solution of both simulations in the same graph, where the individual graphs seem to be on top of each other. The difference is investigated more closely in Figure 13, where it is normalized to the length of the pipe and plotted versus the length of the pipe. The term "difference" is used instead of the term "error," since it cannot be said that the solution from the Offpipe simulation is more accurate than the solution from the current model. The difference is calculated as in (51), but the benchmark solution is replaced with the Offpipe solution.  Figure 13 and Table   Table 4    5, it can be concluded that the results of the current method correspond well to the results from Offpipe.

Dynamic J-lay example
This dynamic J-lay example starts from the result of the static example in Figure 10. Properties for this example are given in Table 4. The first end of the pipe, which is in contact with the seabed, is to remain at a constant position. This is enforced by applying the penalty method in all three directions at this node. The second end of the pipe is considered as were it attached to the center of motion of the pipelaying vessel. The center of motion responds to an irregular wave with H s ¼ 3 m and T z ¼ 7 s; with a wave encountering angle of 110 deg. In finite element terms, the pipe is modeled with 20 elements and 21 nodes. The position of node number 21 is determined by a prescribed displacement following the center of motion of the vessel, which can be seen in Figure 14.
In the dynamic simulation, using the HHT-a method, the time step is set to Dt ¼ 0:1 s: The simulation is done over a period of 180 seconds and thus requires 1800 steps. Each step consists of two or three iterations. The total computation time is approximately 18 seconds. In Figures 15 and 16, two snapshots of the results are depicted. These snapshots are chosen around the peak in the displacement of the center of motion (see Figure 14).
To show the dynamic result in more detail, the response of four interesting nodes is chosen to be depicted in Figures 17-20. It is noted that the displacement u has been set to zero at t ¼ 0 to make the plots better readable. Node A ( Figure 17) is in contact with the seabed. It can be seen that the pipe moves in y-direction in steps due to the stick-slip in the friction model. Node B ( Figure 18) is just above the seabed. Furthermore, node C ( Figure 19) is approximately 10 Figure 11. Convergence of static contact example. The different colored lines illustrate the different substeps.    Figure 14.
Response of the center of motion of a pipe-lay vessel to a wave with H s ¼ 3 m and T z ¼ 7 s: meters above the seabed. Of the four depicted nodes, node D ( Figure 20) is closest to the water surface. The displacements of node C and D have a peak around t ¼ 100 seconds. This is as expected, since the displacements of the vessel's center of motion also have a peak around this time. The responses of nodes A and B show that the horizontal displacement of the nodes is mostly in positive y-direction, making the solution unsymmetric. This is due to the encountering wave angle, which is set at 110 degrees. Numerical experiments have been done with a wave that encounters the ship from the other side, thus with a wave angle of 250 degrees. The responses of nodes A and B showed, as expected, horizontal displacements mostly in negative y-direction.

Conclusion
A nonlinear dynamic corotational finite element model for pipe-laying has been developed. The corotational method accounts for geometric nonlinearity.    Numerical time integration is done by the implicit HHT-a method. A buoyancy model with equivalent moments on all nodes has been implemented. Morison's equation is used for all hydrodynamic forces on the pipe. Contact with the seabed is modeled with the penalty method. Response amplitude operators are used to calculate the response of the pipe-laying vessel to the waves. The nonlinear corotational finite element model presented in this article can simulate J-lay pipe-laying. Validation showed that the results of the current method are in excellent agreement with industry standard software Offpipe. In contrast to Offpipe, the current method does not require a postprocessing step for correcting the axial force.
The implicit method used for numerical time integration works well with buoyancy load, hydrodynamic forces, gravity, and contact. The highly nonlinear contact phenomenon can cause bad convergence due to its on/off nature. In this article, this is solved by decreasing the step size.
A nodal buoyancy moment, resulting from a variational analysis, was introduced. It is an equivalent moment at both nodes of an element, based on the distributed pressure. This additional buoyancy moment showed an increase in accuracy, mainly for models with unequal element sizes.
It was shown that the buoyancy method of Yazdchi was incorrect in deep water, which was caused by a mismatch in cross-sectional area of two adjacent elements. The same error was found on the capped end force. A force was introduced to correct this mismatch and the results with this correctional force were shown to be accurate.

Disclosure statement
No potential conflict of interest was reported by the authors.