{'fraction_non_alphanumeric': 0.022456974954526376, 'fraction_numerical': 0.00430250454736253, 'mean_word_length': 4.9697222802255165, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Recursive thoughts on the simulation of the flexible multibody dynamics of slender offshore structures

In this work, the floating frame of reference formulation is used to create a flexible multibody model of slender offshore structures such as pipelines and risers. It is shown that due to the chain-like topology of the considered structures, the equation of motion can be expressed in terms of absolute interface coordinates. In the presented form, kinematic constraint equations are satisfied explicitly and the Lagrange multipliers are eliminated from the equations. Hence, the structures can be conveniently coupled to finite element or multibody models of for example seabed and vessel. The chain-like topology enables the efficient use of recursive solution procedures for both transient dynamic analysis and equilibrium analysis. For this, the transfer matrix method is used. In order to improve the convergence of the equilibrium analysis, the analytical solution of an ideal catenary is used as an initial configuration, reducing the number of required iterations.


Introduction
In offshore engineering, many long and slender structures are encountered, such as pipelines and risers. These structures are subject of study in a wide variety of applications such as design optimization, operability studies, real time control and so on. The fact that such studies require many design evaluations, long transient analyses or iteration times that are faster than real time, emphasizes the importance of efficient and accurate structural models.
The most important contribution of this paper is in providing an overview of how various methods, that individually are all well-developed and well-documented, can be combined to form a powerful solution strategy for the dynamic simulation of slender offshore structures specifically. These various methods will be briefly explained here, together with their suitability for offshore applications.
For the modelling of slender offshore structures, different approaches can be pursued. In particular for the purpose of pipe-laying, different methods are explained in textbooks such as [1]. For the equilibrium analysis of pipelines, analytical methods based on catenary theory are well-documented. The ideal catenary solution can be found in standard textbooks on statics such as [2] and can be expanded to an elastic catenary in which the effects of bending stiffness are taken into account [3]. It is still common to compare more elaborate models that do not possess an analytical solution, with catenarylike solutions, as for example shown in [4].
In order to study the dynamic behavior of more complicated structural geometries and topologies, and the effects of additional phenomena such as contact mechanics, plasticity and hydro-elastic coupling, a broad range of numerical models is developed. In rigid multibody dynamics (sometimes referred to as the rigid finite element method) the structure is modeled as a series of rigid bodies connected by ideal elastic elements in the form of (torsion) springs. The details of this formulation are explained excellently in the standard work by Haug [5] and an application of this method in offshore engineering is for instance found in [6].
Departing from linear finite element formulations, co-rotational finite element formulations have been developed to cope with the geometric nonlinearities that result from the large displacements [7]. A benefit of this method is the strong relation with linear finite element formulation. However, in particular for three-dimensional problems and systems that move at high velocities, many elements may be required in order to produce accurate results. This is due to the fact that the effects of quadratic velocity terms and elastic deformation are not taken into account in the inertia forces.
Flexible multibody dynamics is a natural extension of rigid multibody dynamics. Within the framework of flexible multibody dynamics, the floating frame of reference formulation is used in this work, which is described clearly in the standard work by Shabana [8]. The key concept is that for each body the rigid and flexible motion are separated and described by distinct generalized coordinates: a floating frame attached to each flexible body describes the body's rigid motion with respect to an inertial coordinate frame, whereas a linear combination of elastic modes describes the local deformations with respect to the floating frame. These elastic modes can be obtained directly from a linear finite element model of the body, which allows the possibility of well-developed model order reduction techniques.
The decision to use the floating frame of reference formulation in this work is a fundamental one, as it emphasizes the natural possibility of including a pipeline or riser model in a multibody dynamics model of the vessel and its on-board equipment and a possible seabed model. How to efficiently include hydrodynamic forces into this flexible multibody model using potential theory is described in [9].
The efficiency of recursive solution procedures in systems that have a chain-like or tree-like topology lies in the fact that the number of operations to solve the equations of motion at a certain iteration is only of ( ), with the size of the system matrix. The majority of the recursive methods is based on the successive elimination of bodies from the chain by condensation. An example of this method for the purpose of dynamic simulations of space structures is found in [10]. Alternatively, the transfer matrix method can be used as a recursive method. The fundamentals of this method for linear elastic problems is well explained in [11]. How to extent this method for geometrically nonlinear problems is described in [12], in which also a comparison has been made with the condensation method. It is found that the transfer matrix method seems to be advantageous for closed-loop systems, as it does not require the use of Lagrange multipliers for loop-closure constrains.
The outline of this paper is as follows: In Sections 2 and 3, the details of the floating frame of reference formulations that are necessary for our present purposes are given. This concerns the kinematics of an elastic body and the development of the constrained equations of motion of a flexible multibody system. In Section 4, it is described how for the particular case of a chain-like structure that contains many bodies with two interface points only, the equations of motion can be rewritten in terms of the absolute interface coordinates. In this way, the Lagrange multipliers can be eliminated from the equations of motion and a formulation in terms of minimal coordinates is obtained. In Sections 5 and 6, it is shown how these equations of motion can be solved recursively using the transfer matrix method for transient dynamic analysis and equilibrium analysis, respectively. To this end, the equation of motion of a single flexible body is rewritten to transfer matrix form after which it is explained how to couple it with neighboring bodies. Section 7, summarizes the basics of catenary theory and explains how the analytical solution of the ideal catenary can be used in the initialization of equilibrium analysis. A numerical example is developed for the demonstration of the equilibrium analysis. The paper is closed with the conclusions in section 8.

Kinematics of an elastic body
Consider a material point located on a flexible body. A coordinate frame is rigidly attached to this material point. The combination of the location of and the orientation of the coordinate frame located in defines a coordinate system. In this work, this will be referred to as the position of . of can be expressed relative to a coordinate frame in another point by the (3 × 1) position vector , and the (3 × 3) rotation matrix . The position vector , defines the location of (lower index) with respect to the location of (second upper index) and its components are expressed in the coordinate system of (first upper index). The rotation matrix defines the orientation of (lower index) with respect to the orientation of . The rotation matrix is an orthogonal matrix of the proper type, i.e. its determinant is equal to 1 and its transpose equals its inverse: From the orthogonality property, it follows that the time derivative of the rotation matrix can be expressed as: where the overdot operator �̇� denotes partial differentiation with respect to time. � , is the skew symmetric matrix constructed from the angular velocity vector , . Here, the tilde operator � � � is introduced such that when applied on a (3 × 1) vector , it yields a skew symmetric matrix �: Similarly as above, the virtual change or variation in the rotation matrix can be expressed as: where delta operator ( ) denotes the variation of its argument. , is referred to as the virtual rotation of with respect to with its components expressed in frame .
In the floating frame of reference formulation, the absolute position of any material point on a flexible body with respect to a fixed inertial reference frame is decomposed in the absolute position of the body's floating frame, denoted here by , with respect to the inertial frame and the relative position of with respect to the floating frame . This decomposition is shown graphically in Figure 1. With the use of the floating frame, the location of with respect to is decomposed as: The orientation of with respect to can be decomposed as: The body's flexible behavior is described locally, with respect to the floating frame. Hence, the position vector , can be expressed as the sum of the position vector of relative to on the undeformed body , and an elastic displacement due to flexibility , : In case of small elastic deformations, the rotation matrix can be linearized: in which only the terms up to the first order in the Taylor series of the rotation matrix are retained. For this, the elastic rotation , that depends on the gradient of the displacement field must be small. Following the theory of linear elasticity, the displacement field , can be expressed as a linear combination of a finite number of modes: where is the matrix of mode shapes and the vector of generalized elastic coordinates. The floating frame of reference formulation allows in principle for any combination of modes. For bodies with a complex geometry, these modes could be for instance the body's natural modes as obtained using finite element software. Modal order reduction techniques, such as the Guyan-Iron method, Craig-Bampton method, etc. can be applied to reduce the number of modes. For bodies with a simple beamlike geometry, standard cubic Hermitian interpolation functions can be used, similarly to the finite element method. In fact, it can be shown that in this case the floating frame of reference formulation for a body is equivalent to the co-rotational formulation for a finite element [13].
In order to obtain an expression for the linear velocity of with respect to , Equation (5) is differentiated with respect to time. Using the appropriate mathematical properties of the rotation matrix and skew symmetric matrix, this can be rewritten to: For the angular velocity of with respect to holds: Equation (10) and (11) are combined to a compact notation: in which the (6 × 1) velocity vectors contain both the linear and angular velocity vectors, [ ] denotes a (6 × 6) compound rotation matrix and [�] is defined as: where sub-and superscripts are omitted above for brevity. The expression for the virtual displacement of relative to can be written in a similar form as Equation (12).

Equations of motion in the floating frame of reference formulation
The equation of motion of a flexible body is derived using Hamilton's principle: where is the virtual kinetic energy, is the virtual internal elastic energy and is the virtual work by externally applied body forces and surface tractions. Substitution of the appropriate kinematic relations in the expressions for the virtual energies in Equation (14) and using the argument that Hamilton's principle must hold for arbitrary intervals [ 1 , 2 ] and for all virtual displacements, the floating frame of reference equation of motion in standard form is obtained: where , is the absolute acceleration of the floating frame. is the rigid body mass matrix, the mass matrix due to elastic modes and are the modal participation factors that couple rigid and elastic is the generalized stiffness matrix due to elastic modes. The vector contains the inertia forces that are quadratic in the velocity and contains the externally applied generalized forces and modal forces. In a multibody system, the equations of motion of all individual bodies can be combined and written in the following compact form: in which the quadratic velocity forces and elastic forces are contained in , as well as the unknown constraint forces that result from kinematic constraints between distinct bodies or between a body and the fixed world. In order to obtain a solvable set of equations, the corresponding kinematic constraint equations need to be taken into account. Holonomic constraint equations are considered in the following form: Differentiation with respect to time twice yield the acceleration equation of the constraints: where is the Jacobian matrix of the constraint equations. The combination of Equation (16) and (18) forms the constraint equations of motion in Lagrange multiplier form, which is of the differential-algebraic type: in which the constraint forces are removed from and are the Lagrange multipliers.

Transformation to interface coordinates
The constrained equations of motion in Lagrange multiplier form Equation (19) are valid for any flexible multibody system subjected to holonomic constraints. In co-rotational finite element formulations, no Lagrange multipliers are present. This is due to the fact that the equations of motion are written in terms of the absolute nodal coordinates. Since constrains are applied at the nodes, they result in a direct relation between nodal coordinates or they prescribe nodal coordinates to have a certain value. Based on co-rotational finite element formulations, it is observed that it is possible to eliminate the Lagrange multipliers from the floating frame of reference equation of motion, if it is possible to establish a coordinate transformation from the floating frame of reference coordinates , to the interface coordinates corresponding to the interface points at which the kinematic constraint equations are applied.
In general, such a coordinate transformation is not straightforward. Previously, it was established by locating the floating frame in an interface point [14] or by assuming the position of floating frame to be the weighted average of the interface coordinates [15]. More generally, this coordinate transformation can be established by demanding zero elastic deformation at the location of the floating frame [13].
However, the bodies that form slender offshore structures in a typical chain-like topology, often have two interface points only: connecting the body to its two immediate neighbors. In this specific case, the desired coordinate transformation can be established easily for the case that the floating frame coincides with one of the interface points, provided that the elastic behavior is described by an appropriate set of modes. To this end, consider a flexible body with interface points and . Let the body's floating frame coincide with interface point . Using Equation (12), the absolute velocity of interface point can be expressed as: In the specific case that the elastic displacement field is described using the 6 static interface modes of a Craig-Bampton reduction, the local velocities , equal the time derivative of the generalized elastic coordinates ̇. This is an immediate consequence of the fact that a static interface mode evaluated at Since the absolute coordinates of the floating frame of reference coincide with the absolute nodal coordinates of node , the following coordinate transformation is obtained on the velocity level: Differentiating Equation (22) with respect to time yields the coordinate transformation on the acceleration level: ̈=̈+̇̇. ( Applying this coordinate transformation on the equations of motion Equation (16) results in: in which � is the transformed mass matrix. � is the corresponding generalized force vector in which also the squared velocity term due to the transformation Equation (23) is included.
With the equations of motion in the form of Equation (24), constrains between bodies can be applied direcly, eliminating dependent generalized coordinates from the equations. The resulting equations of motion are expressed in the minimal coordinates, i.e. the number of equations of motion equals the number of independent generalized coordinates.
For the sake of completeness, it is mentioned that � is not constant, but depends on both the rotation matrix as well as the body's elastic deformation. For small deformations, it is often reasonable to neglect the dependency on the elastic deformation. In [12] it is shown that under this assumption � can be related to the body's finite element mass matrix that is condensed to the interface points: Note that in co-rotational finite element formulations, the global mass matrix is also assembled by rotating the local finite element mass matrices, which are constant, to the global frame. In the floating frame of reference formulation, as presented in this work, the dependency of the mass matrix on the elastic deformation can, but does not need to be neglected: the terms follow naturally from the formulation. Also the quadratic velocity terms are still included in � , these terms are typically not taken into account in a co-rotational finite element formulation.

Transfer matrix method for dynamic analysis
As a consequence, of the chain-like topology of many offshore structures, recursive solution procedures are well-suited for dynamic simulations. In this work, the transfer matrix method is used for solving the geometrically nonlinear equations of motion. To this end, consider the equation of motion Equation (16) for a single body, expressed in its interface coordinates, in the following partitioned form: where and denote the constraint forces at interface points and respectively and are separated from the remaining right hand side forces. In the transfer matrix method, Equation (26) is rewritten such that the acceleration and interface force at are expressed in terms of the acceleration and interface force at : This is of the form: where and are the so-called state vectors of interface points and respectively, is the transfer matrix from with respect to and � is an internal state vector due to elastic and quadratic velocity forces.
When the equation of motion of each flexible body in the chain is written in the form Equation (28), the bodies can be coupled by considering the appropriate coupling conditions. For example, in case of a rigid coupling, the accelerations of the connected interface points are equal and the interface forces are equal but of opposite sign, following Newton's third law. Mathematically, such a coupling can be realized by a permutation matrix : where ( ) and ( +1) are the state vectors of the interface point that bodies and +1 have in common.
In this way, the state vector of last interface point of the chain, say can ultimately be related to the state vector of the first interface point 1 : At this point, the boundary conditions at both ends of the chain can be applied. In the case of pipelaying one could for instance prescribe zero deformation at the location of the sea bed, prescribe the water depth and the force applied by the tensioner system. Then, from Equation (30) the unknown displacements and constraint forces at both ends of the chain are determined. All intermediate states can be solved recursively by back substitution. The solving of all state vectors is then followed by numerical time integration to the next time step.
Because the transfer matrix method is a recursive solution procedure that requires ( ) computations, it is potentially a very efficient algorithm. A benefit of the transfer matrix method as presented in this work over many other recursive solution procedures, is that no knowledge of boundary conditions is required to set up the recursive equations. In particular, no loop-closure constraints are required in order to solve structures with a closed-loop topology. Such a loop-closure constraint would be required in a condensation-based recursive method if one prescribes the water depth.

Transfer matrix method for equilibrium analysis
Analyses that are performed with the purpose of e.g. operability studies, require transient dynamic simulations. In order to prevent unwanted start-up behavior, the transient simulation is started from a static equilibrium configuration. The equilibrium equations of a single body are obtained directly by removing all inertia related terms from the left hand side of Equation (15). Also for this static case, the assembly of the equilibrium equations of the entire system can be done using Lagrange multipliers. Due to the geometric nonlinearities on the position level, it is not possible to establish an explicit coordinate transformation to interface coordinates. As a consequence, it is not possible develop a state transfer matrix on the position level.
It is common however to solve nonlinear equilibrium equations using a total Lagrangian or updated Lagrangian formulation, in which the equations are solved iteratively using a Newton-Raphson scheme. Within each iteration, the incremental change in the degrees of freedom Δ is solved from a set of linearized equations in the tangent space of the equilibrium equations. For increment + Δ these equations can be written as [16]: where is the tangential stiffness matrix, which consists of the material stiffness matrix, as determined in Equation (15), as well as the geometric stiffness matrix that follows from taking the variation of the nonlinear equilibrium equations. The right hand side is the difference between the externally applied forces of the next iteration ( +Δ ) and the internal elastic forces of the current iteration ( ) . If at the current iteration, the nonlinear equations are solved exactly, the right hand side equals the incremental change in the externally applied forces Δ . However, in the notation used in Equation (31), a possible residual at the current iteration is taken into account during the next iteration.
The linearized equilibrium equation of the entire structure can be assembled using Equation (31) for each body in the chain. At this point it is possible to set up a transfer matrix method for solving the equilibrium equations. In this procedure, a state vector would consists of the position of an interface point and the interface forces. The remainder of the solution procedure is identical as the procedure for a transient dynamic simulation.
In order to determine the equilibrium configuration of an offshore pipeline that is suspended from a vessel and supported by the sea bed, an initially straight pipe can be considered. This pipe can be incrementally loaded, solving Equation (31) at each load increment, until equilibrium is reached.
To demonstrate this procedure, a numerical example is performed on pipe that is subjected to its own weight in downward direction and a vertical tip force equal to its own weight in upward direction. The pipe is divided into 10 bodies. The total length of the pipe is 1 (m). The cross section has radius 0.01 (m) and the wall thickness is 0.001 (m). The Young's modulus is 70 (GPa). The applied weight per unit length is 100 (kN/m). This load is such that the deformation of the pipe will be large, such that geometrical nonlinear effects need to be included. The pipe is initially straight and the equilibrium configuration is obtained using 5 equal load increments. Within each load increment, iterations are performed until convergence is reached. Figure 2 shows the converged equilibrium configurations of each load increment.

Equilibrium analysis using catenary theory
To solve the equilibrium configuration of a suspended pipeline from a pipe that is initially straight, large displacements are required. As a consequence, the use of load increments and the Newton-Raphson procedure are required, which may need quite some iterations before convergence is reached. It could be said that an initially straight pipe is a very poor initial estimation of the equilibrium configuration.
It is found that the equilibrium configuration of a pipe line is reasonable well approximated by the analytical solution from catenary theory. For instance, in [4] results are presented that in case for a J-lay process show the striking similarity between the ideal catenary solution and the pipeline. In order to benefit from the fact that the catenary theory yields an analytical solution, many references can be found in which this theory is expanded to include effects such as longitudinal straining, finite bending stiffness and inertia effects. In this work, this semi-analytical approach is not followed on purpose, in order to emphasize the possibility to include the pipeline into a multibody simulation. However, using the ideal catenary solution for the initialization of the equilibrium analysis could greatly reduce the effort required in this stage. The nonlinear differential equation of a standard catenary is [ in which ( ) is the shape of the catenary, the weight of the cable per unit length and 0 the constant horizontal component of the pretension in the cable. The solution of Equation (32) is: At every position, the slope of the catenary can be determined from: The horizontal component of the tension force equals 0 at every position. The vertical component of the tension force equals the weight of the cable up to that point, which can be expressed as: The internal bending moment in the catenary is zero by definition. For a given water depth and pretension 0 , the solution of the catenary is determined. Then, Equation (32) to (34) are evaluated at the interface points to establish the state vectors at these interface points. These state vectors are used as an initial guess in the iterative Newton-Raphson procedure that solves the equilibrium configuration of the actual multibody system. Also in the case of the numerical example used in the previous section, it is found that the catenary solution is close to the equilibrium configuration of the pipe. To demonstrate this, Figure 3 shows the final equilibrium configuration of the pipe again (solid line), together with the analytical catenary solution that passes through the end point (dashed line). The close resemblance of the catenary to the pipe line suggests that it would indeed be beneficial to use the catenary solution as an initial guess in the equilibrium analysis of the pipeline. In this way, the use of load increments could be prevented and only few Newton-Raphson iterations are required.

Conclusion
In this work, thoughts have been presented on how to model slender offshore structures using the floating frame of reference formulation, suitable for the fast simulation of the structure's flexible multibody dynamics. Starting with the standard floating frame of reference formulation, a coordinate transformation to absolute interface coordinates is applied in order to remove the Lagrange multipliers from the equations of motion. This allows for the use of potentially efficient recursive solution procedures, such as the transfer matrix method. For the equilibrium analysis, it is explained how the use of the ideal catenary solution could be used to prevent the use of load increments. In this way, the number of iterations can be reduced and convergence improved. Future work is aimed on the validation of the proposed approach with simulations published in literature. In this, the effect of using a recursive scheme in combination with the use of the ideal catenary for the equilibrium analysis on computation times will be investigated. Once reliable results have been obtained, a pipeline model can be coupled to a multibody dynamics model of the pipe-laying vessel. At this point, full hydrodynamic coupling can be included. In this way the transient dynamic behavior of the combined system subjected to various sea-states can be studied with great numerical efficiency.