Application of the finite segment method to stabilisation of the force in a riser connection with a wellhead

This paper presents a spatial model of riser dynamics formulated using the segment method and its applications. The model has been validated by comparison of the authors’ own results with those obtained from experimental measurements and Abaqus on the basis of forced vibration with large amplitude for the riser submerged in water. The influence of the sea environment is considered. Correctness and numerical effectiveness of the model enable us to formulate and solve the force stabilisation problem. A dynamic optimisation problem is formulated and solved. As a result vertical courses of movement of the upper end of the riser are obtained which compensate the horizontal movement of the base and stabilise the force in the connection of the riser with a wellhead.

carried out before production. Risers are long elements with bending and longitudinal flexibilities; these risers are loaded with significant forces acting at their ends and large transverse forces caused by sea currents and hydrodynamic forces. Deflections of risers can be very large, and they may undergo elastoplastic deformations [1].
In this paper, we present the authors' own formulation of the finite segment method for discretisation of a link with bending flexibility. The finite segment method has been used in modelling risers and cables for more than 50 years [7,10,17]. Its main idea lies in dividing a slender link into rigid elements connected by spherical joints treated as massless spring-damping elements in order to reflect torsional and bending flexibilities, and sometimes also longitudinal flexibility of the discretised links [18]. A similar approach derived from the rigid finite element method [19,20] is used by Drąg [15,16].
The finite segment method belongs to a group of methods used for dynamic analysis of slender systems, in which a continuous system is replaced by a system of rigid bodies reflecting inertial features of the system modelled. These rigid bodies are not point masses as in the lumped-mass method, but elements (segments) for which masses and moments of inertia are accurately determined. Using this method, it is not necessary to formulate partial differential equations for dynamics, but the equations of motion are usually derived from the Lagrange equations by means of multibody-system techniques. The choice of generalised coordinates is an important issue. Drąg in his papers [15,16] uses absolute coordinates, which are three translational coordinates (describing the position of the chosen point of each rigid element) and up to three rotational coordinates describing the orientation of the rigid element. The translational coordinates are then eliminated by introducing reactions and constraint equations. The advantage of using the absolute coordinates is a pseudo-diagonal mass matrix, while the large number of coordinates and the necessity of stabilisation of constraint equations are disadvantages of the method. In the paper, we propose to use joint coordinates. In this case, the vector of generalised coordinates has a minimum number of components; these components are three coordinates of a chosen point of the link discretised and up to a three coordinates describing the orientation of each segment. This means that it is not necessary to introduce reactions and formulate constraint equations. Consequently, there is no need to stabilise the geometrical constraint equations as in the case of absolute coordinates. However, the disadvantage of the method is a full mass matrix. When large deformations are analysed, the mass matrix obtained using both approaches depends on generalised coordinates. In general, the finite segment method with the use of joint coordinates is more numerically effective for a moderate number of segments (less than 100) than when the absolute coordinates are used.
The finite segment method is admittedly not in all respect more efficient than the finite element method applied in commercial software packages, including those for simulation of statics and dynamics of offshore appliances. However, as far as analysis of slender links (cables, risers) is concerned, the finite segment method poses important advantages, such as simple physical interpretation, the ease of describing large deformations of links, and the possibility of special applications, including those presented in the paper.
The authors' formulation of the finite segment method is implemented and validated by comparison of their own results with those presented by other researchers, especially those using experimental measurements, as well as with results obtained using professional software. Two experiments are consideredone examines natural frequencies of a riser under different tensions [21,22], and the other is concerned with measurements of vibrations caused by oscillatory motion of the upper end of the riser [23][24][25]. Another method of validation is to compare the results of authors' own calculations for forced vibrations of a riser with those obtained from Abaqus. In all cases of validation, good compatibility of results has been achieved. The validated model is then used to stabilise the force in the connection between the riser and a wellhead. The problem is formulated as a dynamic optimisation problem. The aim is to choose vertical displacements of the upper end of the riser which, despite the horizontal movement of the base (platform or vessel) and sea current, ensure that the force in the connection of the riser and the wellhead is kept as close as possible to the static position. It is assumed that the base follows a path on a surface [26]. The considerations are limited to compensation of horizontal movements of the upper end of the riser, since compensation of vertical movements is not difficult [16]. Optimisation calculations are carried out for an almost straight riser. The results obtained are interesting due to the possibility of the method presented being used for preparation of learning sets for ANN (artificial neural network) and its potential application to control of HCS (heave compensation system).
Analysis of literature concerned with modelling of risers shows that most attention is devoted to analysis of static deflections and bending vibrations both in 2D and 3D analysis. Some research, however, shows that it is necessary to consider torsion. For example, [5,27] demonstrate how to take into account the bendingtorsion coupling effect, especially in the case of seabed interaction. The finite segment formulation presented in this paper, similar to [25,[28][29][30], omits the torsional flexibility of the riser model.
Most papers dealing with modelling of risers consider longitudinal flexibility. Consideration of elongation of slender systems is especially necessary when the analysis is concerned with top-tensioned risers or catenary lines. The elongation of the riser (500 m long) considered in this paper does not exceed 3 cm. This means that the influence of longitudinal deformations on displacements of the riser is negligible, which is why considerations in this paper are limited to analysis of spatial bending vibrations. Consequently, forces acting at the ends of the riser can be calculated in a different way (to be explained in Sect. 2).

Model of a slender link using the finite segment method
The basic idea of the finite segment method lies in the discretisation of a flexible link into rigid elements (segments) reflecting inertial features of the link. The segments are connected by means of joints consisting of massless and nondimensional spring-damping elements reflecting stiffness characteristics of the link (Fig. 1). When the longitudinal flexibility is omitted, the following relation is true: This means that the length of the deformed link along its curvature is constant and it enables quantity shown in Fig. 1a, caused by large deformations of the flexible link, to be taken into account. In order to define the position of any point lying on a neutral axis of the link in segment i, notation presented in Fig. 2 is introduced. The position of point P of segment i lying at distance ξ from point P i depends on orientation of axes of reference coordinate system {i} with respect to axes of inertial frame of reference {}.
In this paper, Euler angles ZY X presented in Fig. 3 are used. These angles define absolute rotations, in contrast to other papers [7,18], in which Euler angles ZY X define orientation of system {i} with respect to the preceding segment {i − 1} .
The order of rotations which leads the end of segment i from position P 0 i+1 (when segment P i P 0 i+1 is parallel to axis x of the inertial frame) to position P i+1 is as follows: -first rotation about axis z parallel to z of angle ψ i , -then rotation about axis y of angle θ i , -finally rotation of angle ϕ i about axis x i if torsion is considered.
It should be noted that the rotations are performed about axes assigned to the segment. Such angles of rotations are often used in airplane mechanics and are called heading, attitude and bank angles [31].
Coordinates of point P of segment i can be defined by the formula: where r P i are coordinates of point In further consideration, torsion is omitted (angles ϕ i ≡ 0) and thus rotation matrices R i take the follow-I. Adamiec-Wójcik, S. Wojciech If we assume that points P 0 , . . . , P n lie on the neutral axis of the link, then coordinates of point P i can be defined as follows: where φ j = ⎡ ⎣ cψ j cθ j sψ j cθ j −sθ j ⎤ ⎦ is the first column of matrix R j , and coordinates of point P are defined by the relation: The system of segments 0, . . . , n is defined by components of the vector with N = 3 + 2 (n + 1) elements: where is the vector of coordinates describing relative motion of the ith segment.
The equations of motion are derived from the Lagrange equations: where T is kinetic energy of the link, V is the potential energy of the link, D is the function describing dissipation of energy (damping), q j ,q j are generalised coordinates and velocities, respectively, Q j are generalised forces and N is the number of generalised coordinates. The kinetic energy of the i-th segment is defined by the following expression: where ρ is the material density, A i is the cross-sectional area of segment i. When the riser is filled with fluid, and the added mass of water for a segment submerged is considered, the following takes place: (9) where A i,r is the cross-sectional area of the riser, A i,inn is the internal cross-sectional area of the riser, A i,out is the external cross-sectional area of the riser, r is the material density of the riser, F is the density of the fluid filling the riser, w is the density of water, c M is the inertial coefficient.
Having taken into account (5) after some transformations, energy T i can be expressed in the form: The kinetic energy of the whole link can be written in the following form: where m = n i=0 m i is the mass of the link. The following identities can be proved by induction: In view of (12), the kinetic energy can be expressed as follows: It should be noted that coefficientsb i andd i j are constant.
After differentiating φ i defined by (4) with respect to time t, the following is obtained: The above relations enable us to define the following: The potential energy of gravity forces of the link is defined by the following expression: where g is the gravity acceleration, y Ci = y 0 + i j=0 b i j sψ i cθ i are coordinates of the centre of mass of segment i.
Having used (12a), we obtain: The energy of spring deformation and its derivatives can be calculated in a similar way to that presented by Wittbrodt et al. [20] and Drąg [16].
It follows from Fig. 4 that the energy of spring deformation of a link with circular cross section can be expressed as: where c i is the stiffness coefficient of bending, α i is the angle between axes x i−1 and x i dependent on q i−1 and q i . Angles ψ i and θ i can be calculated in the form presented in [16]. Let where and since due to orthonormality of matrix R i , the following takes place: then elements of this matrix can be calculated when q i−1 and q i are known. From (19), the following results are obtained: Angles ψ i and θ i can be calculated from the relations: Quantities α i , ψ i and θ i are related in the following way: and thus [20]: Having taken into account (15), (17) and (24), the equations of motion of the link without consideration of the influence of the sea can be presented in the form:

Influence of the water environment
The influence of the sea on the riser consists of the following forces: -buoyancy forces, -hydrodynamic forces, -added mass forces.
Buoyancy forces acting only in the vertical direction can be calculated as: Having assumed that this force acts at the centre of mass of the segment, its coordinates are expressed as follows: From: we can obtain generalised forces resulting from buoyancy forces acting on segment i in the form: In order to calculate hydrodynamic damping forces, it is assumed that they can be expressed as in [26]: T is the relative velocity of water at point P of the riser (Fig. 2), Velocity v can be calculated from the formula: where v =ṙ 0 + i−1 j=0 l jφ j + ξ iφi is the absolute velocity of point P, v w is water velocity (sea current).
Force dF d,i in the inertial frame of reference is expressed as: and its resulting generalised forces are: where r = r i defined by (2), j = 0, . . . , i − 1, and I 3x3 is the identity matrix, Forces Q d,ir and Q d,i j are the sum of forces dQ d,ir and dQ d,i j : If we introduce the following notation: then after some transformations the following is obtained: Integrals from (35) are calculated by means of Gauss formulae.
Having taken into account that: the following formula for hydrodynamic damping forces is obtained: Having taken into account (25), (29) and (38), the dynamic equilibrium equations of the riser can be written in the form: where vector of external forces acting at point P 0 (upper end of the riser). When the motion of the upper end of the riser is known as a function of time: which means that vectorsṙ 0 (t) ,r 0 (t) are also known, and Eq. (39b) can be written as follows: and forces causing this motion F ext 0 can be calculated from (39a): Equations (39) or (41) form a system of m = 3 + 2 (n + 1) or m = 2 (n + 1) nonlinear differential equations of the second order with respect to t. The equations of motion are solved using the Runge-Kutta method with a constant integration step.
In the applications presented further in the paper, constraints are imposed on the riser, which involve placing its bottom end in a strictly defined position (Fig. 5) by means of the formula: In such a case, the system considered includes additional reactions as components of vector F E and con- where D i = − l i B T i , and Eq. (43) can be written in the acceleration form as: The Baumgart method is used for stabilisation of constraint equations. Initial conditions for the dynamics problem are calculated from Eqs. (39) or (41) and (43) as a solution of the static problem (ṙ 0 =r 0 = 0;q i = q i = 0, i = 0, 1, . . . , n). The system of nonlinear equations has been solved by means of Newton's iterative method. Natural frequencies with and without tension at the upper end of the riser are calculated in Sect. 3.1. A procedure similar to that presented in [4,12] is applied in order to calculate frequencies and avoid linearisation of Eqs. (39) or (41).

Analysis of frequencies of free vibrations
Analysis of vortex-induced vibrations of risers is essential due to the possible fatigue damage caused by the vibrations. Thus, prediction of natural frequencies of risers is one of the crucial points in the design process of risers. The method presented enables us to calculate frequencies of free vibrations of the riser. In order to validate the model, the natural frequencies obtained using our model are compared with those from the experiment explained in [21], where Chaplin and coauthors describe experimental measurements taken on a model of a riser placed in a water tank and exposed to a stepped current. The data of the riser are presented in Table 1.
The installation of the riser consisted of universal joints at each end of the riser and a tensioning system at the upper end. Measurements of the first seven natural frequencies were taken for five different values of the top tensions.
The results of any discretisation method depend on the number of elements into which the system under consideration is divided. Thus, first the influence of the number of finite segments on the results obtained is examined for two different sets of tension applied at both ends of the riser. The results presented in Tables 2  and 3 show convergence of the results obtained for the segment method and its efficiency.
The comparative analysis is carried out not only with the experimental measurements given in [21] but also with results obtained by means of the variational iteration method (VIM) presented in [22], whose authors use the same experiment. It is assumed that the number of segments used for numerical simulations is n = 200. Table 4 presents comparison of the first seven natural frequencies calculated for five cases of different values of tension. The results calculated using the finite segment method are in very good agreement with both those measured in the experiment and calculated using the variational iteration method. Figure 6 shows the relative percentage error of the simulation results compared to the measurements calculated according to the formula: where ω i,exp , ω i,calc are values of the i-th natural frequency measured and calculated, respectively. Analysis of the diagram shows that the error does not exceed 2% for most calculated values. The largest differences are obtained for the first case with the lowest tension, yet still the error is less than 6.5%, and in one case the second natural frequency is calculated with an error of about 10.5%. The distribution of errors is similar to that presented in [22], but it can be seen that in most cases the results from the finite segment method Fig. 11 Courses of the end of the riser in y direction a in the air, b in water are closer to the experimental measurements than those obtained using the variational iteration method. The analysis presented generally confirms the applicability of the finite segment method for prediction of natural frequencies of rises.

Regular oscillatory vibrations
In this section, validation of the method is carried out on the basis of experiment used as a reference by many researchers [23][24][25]. A laboratory test stand consists of a long thin tube filled with water simulating a riser submerged in water with its upper end attached to a mechanism generating oscillatory motion along the x axis (Fig. 7). Additional weight W is attached at the bottom of the riser in order to keep the model straight in still water.
Cameras installed along the y axis measured x displacements at chosen points. Parameters of the riser are presented in Table 5.
Oscillatory motion exerted on the riser at point P 0 = P with amplitude 0.1 m is shown in Fig. 8 as x P . Numerical simulations were carried out for a riser divided into n = 20 elements. Comparison of displacements in x direction at points A, C and E (Fig. 7) is presented in Fig. 8. The results of calculations are carried out using the finite segment method are compared with experimental measurements and simulation results presented in [23].
The comparison presented in Fig. 8 demonstrates good compatibility of the results obtained using the method presented in the paper with those from experiment and numerical simulation based on Hamilton's principle [23].

Irregular forced vibrations
Validation of the model presented has also been carried out by comparing the results obtained from the authors' own software with those from Abaqus. The compari- Fig. 12 Courses of the end of the riser in z direction a in the air, b in water son is concerned with forced vibrations of the riser both with and without consideration of the sea. The parameters of the riser analysed are presented in Table 6. Displacements of the free end of the riser loaded ( Fig. 9a) with two forces changing in time and acting along x and z axes and with constant force F y = − 30 kN are analysed. Courses of forces changing in time are presented in Fig. 9b.
Calculations have been carried out using the model derived by means of the segment method (SM) dividing the riser into n = 50 segments. First order, threedimensional beam elements B31H have been used for discretisation in Abaqus; these elements enable longitudinal and transversal deformations to be taken into account. The same number of elements for both methods has been assumed. Displacements of the end of the riser in directions x, y, and z both without and with consideration of the sea are presented in Figs. 10, 11 and 12, respectively.
Analysis of displacement courses of the free end of the riser in interval <0, 100 s> in all directions shows that results obtained by the finite segment method are compatible with those obtained using the finite element method (Abaqus). Enlargements presented in each figure show maximal and minimal displacements, and their values are presented in Table 7.
Differences between maximal and minimal displacements obtained by both methods for all directions when the riser is divided into 50 segments are listed in Table 8.
The largest difference (20 cm) is related to displacements along the z axis when vibrations in air are considered; however, it has to be noted that displacements in this direction amount to 100 m, which means that the relative error is 0.2%. Projections of the trajectory of the end of the riser onto plane xz are also compared, and it can be seen from Fig. 13 that the shape is similar for vibrations both without and with consideration of the sea.
It can be seen that the amplitudes when hydrodynamic damping is considered are almost an eighth of those obtained for vibrations in the air.

Stabilisation of the force in the connection with the wellhead
Forces occurring in the connection of a riser with a wellhead can considerably change in heavy seas and strong sea currents. Offshore equipment must be durable and reliable; therefore, it is desirable to ensure that values of these forces are as stable as possible within an accepted range despite the vertical and horizontal motion of the base (platform or vessel) .  The model proposed is numerically very effective and allows dynamic analysis for large displacements to be carried out. For example n = 25, t ∈ 0, 50s , for with integration step h = 0.01s, calculation time does not exceed 4 s. In order to illustrate the possibilities of the model, a dynamic optimisation task is formulated and solved. The problem involves selecting vertical movements of the upper end of the riser (Fig. 14) so that for a given horizontal motion of the base the forces in the connection with the wellhead are stable.
The upper end of the riser (point P 0 ) moves in the xz plane, reaching ± 7 m and ± 3 m alternately and following the trajectory presented in Fig. 15, while displacements and velocities in x and z directions change as in Fig. 16.
Two cases are considered. When simulations start, the upper end of the riser is placed at point P 0 (7 m, 0 m, 0 m), while coordinates of its bottom end are E(− 10 m, − 499.39 m, 0 m) in case A, which means that at the beginning the riser is initially slightly bent (Fig. 17a) or E(− 100 m, − 477.16 m, 0 m) in case B when the riser is initially bent (Fig. 17b).
In case B, a rotary spring with stiffness coefficient equal to 1 × 10 6 Nm is placed at point P 0 in order to keep the deflection angle of the upper end of the riser within the limits recommended in the literature [32].
The dynamic optimisation task is formulated as minimisation of the functional: where p i = y 0 (t i ), i = 0, . . . , m are values of the function which defines the vertical motion of the upper end of the riser at equally spaced points within time period <0, T >, F E is the value of the resultant force occurring at the bottom end of the riser (E), F 0 E is the initial value of the force, and T is the calculation time. Moreover, it is assumed that the minimum of functional  ical effectiveness of the model is especially important because the solution of the dynamic optimisation problem formulated in such a way requires calculation of the functional for each combination of parameters p 1 , . . . , p m−1 , and thus integration of the equations of motion. Numerical calculations are carried out assuming T = 48 s, and interval is <0, T > divided into subintervals m = 12, the first three of 1s, the next seven of 3 s and last two 6 and 12 s. This means that the displacements of upper end of the riser are analysed not only during the motion of the base when t ∈ <0, 36s> , but also when the base is stationary for t ∈ <36s, 48 s>. It is assumed that m = 11 and minimal and maximal values, respectively, are equal to the following: It is assumed that at the initial approximation (before optimisation) there are no vertical displacements at P 0 i.e. p i = 0 for i = 0, 1, . . . , m. The comparison of the force courses at point (E) in the connection of the riser with the wellhead, obtained before and after solving the optimisation task, is presented in Fig. 18a. Due to the vertical motion of the upper end of the riser, the force in connection E decreases from 1150 kN to 500 kN; thus it is less than half of the initial value. Figure 18b shows force F ext 0 acting at the upper end of the riser and corresponding to initial approximation y 0 ≡ 0 also before and after optimisation.
This considerable decrease in the force in connection E has been obtained as a result of applying vertical motion at the upper end of the riser in the form presented in Fig. 19.
The results presented are illustrative. Shapes of the riser at different moments of time before and after optimisation are presented in Fig. 20.
Results of calculations for case B are presented in Figs. 21 and 22. In this case, reduction of the force acting at the bottom end of the riser (point E) is not as great as in case A. The reason is that when the riser is almost vertical (case A) forces acting at point E, when the upper end of the riser move about 14 m in x, z directions in a relatively short time, lead to a large increase in longitudinal forces. In case B, however, when the riser is significantly bent, the longitudinal forces are considerably smaller. It is important to note that as a result of the optimisation procedure for case A (slight bending of the riser) force F E is about a half of its initial value,   yet still its optimal value is larger than the initial value of force F E for case B (large initial bending of the riser) .

Conclusions
The paper presents a different formulation of the finite segment method which uses absolute angles in contrast to the classical finite segment method, in which relative angles are used in order to describe the geometry of a slender system. Our method enables us to take into consideration the shortening of the distance between the joints, which means that the length of the slender link measured after deformation along the curvature is constant. The influence of the sea is taken into account. The method is validated by comparing the results obtained with those obtained from experimental measurements and Abaqus. The validation process shows the good correspondence of the results obtained by SM, both in natural frequencies and oscillatory vibrations, with results of measurements and calculations when different methods are applied. In addition, the comparison of results obtained with the method presented and those obtained using professional finite element method software confirms the correctness of the model and its numerical effectiveness. Even when the equations of motion are nonlinear, the integration time for less than 50 segments does not exceed 20 seconds over a time interval of 50 s with an integration step of 0.01s. Good numerical effectiveness of the method allows it to be used in solving the dynamic optimisation problem. In order to solve the optimisation problem, the equations of motion are integrated at each optimisation step (for each combination of parameters describing the objective function). In the application presented, the disturbances caused by horizontal motion of the base are neutralised by vertical motion of the upper end of the riser; this solves the problem of stabilising the force in the connection between the riser and the wellhead. The results of calculations demonstrate that even for unfavourable initial approximation of the optimisation problem (such as assuming that the upper end moves horizontally and slight initial bending) the force in the connection of the riser and wellhead is reduced considerably as a result of optimisation. Such a solution is likely to be feasible by using HCS (heave compensation system). Due to the very short calculation time, the method can be applied when teaching a neural network, which then can be used in control of the system in real time.