Application of dynamic optimisation to stabilise bending moments and top tension forces in risers

The study discusses the problem of determining vertical displacements of a riser’s ends, which, despite its horizontal displacements induced by waves, mitigate stresses. A spatial model of riser dynamics is presented that considers the geometric nonlinearity due to large deﬂections. The Rigid Finite Element Method was used for riser discretisation. Analyses are reported that enabled riser’s vibration frequency determination according to the positions of its upper and lower ends. Then a dynamic optimisation task was formulated and solved. It consists of the selection of riser’s vertical displacements that provide bending moments’ stabilisation at its selected points, tension forces, despite the deﬁned horizontal movements of the riser’s upper end induced by waves. The calculations were performed for variable amplitudes of riser end’s horizontal movements.

cal motion of the base and its horizontal motion caused by sea weaving. In all cases using the HCS, it means using moving vertically the upper end of the raiser. An important issue in the riser design is to stabilise the load despite the movements of the unit (platform or vessel) to which the riser's upper end is connected. While compensation of the vertical movements caused by waves is relatively easy to implement, the horizontal motion compensation is much more difficult. To compensate vertical movement are often used solution that are combination of passive and active hydraulic systems [1][2][3][4]. The active heave compensation system based on linear control techniques, which have been mounted on drillship was object of paper [1]. A nonlinear controller was proposed for active compensator based on an electric hydraulic drive in [2]. In both papers the vertical movements caused by waving sea are compensated by AHC (Active Heave Compensation). In order to avoid a collision between two raisers, the [3] have been proposed control based on adjusting the top tension. Effective length of raisers is selected by the controller which is taking into account the forces of water environment impact. One the proposals of application of active compensation system is to minimalise of raiser angle deviation in connection with the well head and at the top joint. The paper [4] presents a mathematical model of the hydraulic system that allows to reduce of the upper and the lower raiser angle. The purpose of this study is to demonstrate that riser end's vertical movements can mitigate the impact of the unit's horizontal move-ments on bending moments at selected riser points. This may be important if some areas along the riser length need to be protected against overloads. A similar problem with regard to planar was considered in the study [5]. To formulate a riser dynamics model in this study, a modification was applied of the Rigid Finite Element Method (RFEM) presented in [6,7] for planar systems, and in [8] for spatial systems. The subject of the analysis was the riser presented in [9], which features large deflections due to substantial horizontal distance between the positions of its upper and lower ends. The statics calculation results were compared with the result of the aforementioned study [9]. Also analysed was the impact of riser's length and of the horizontal distance between its upper and lower ends directly affecting its curvature and natural vibration frequency. Calculations were also completed to demonstrate how the number of the elements into which a riser is divided affects the accuracy of vibration frequency results.
The nonlinear model of riser dynamics was used to formulate and solve the optimisation task. Harmonic motion was assumed in the horizontal plane with known amplitude and frequency. A spline function was sought that determines the riser end's additional vertical displacements, which ensure that the bending moment at a selected point shall be the least different from its initial (static) value. At each optimisation step, the riser dynamics equations needed to be integrated. The calculations were performed for various amplitudes of the base's horizontal movement.
An extensive discussion on the state of the art in modelling and on problems in total dynamic analysis of offshore systems is presented by Chakrabarti [10]. Although Starossek [11], Jain [12], Patel and Seyed [13] report many different modelling methods, the most popular methods used are the finite element methods [14,15] and the lumped mass methods [16,17]. Those methods are used in commercial packages such as Riflex and Orcaflex. The Flexible Segment Model (FSM) presented in [18,19] is another similar approach. The RFEM used in this paper has been also successfully used for dynamic analysis of flexible systems containing beam-like links [20,21], plates and shells [22,23]. Figure 1 shows the analysed riser.

Riser dynamics model
The riser was divided into rfe and sde as described in detail in [8,24], into n + 1 rfes (numbered 0 . . . n) and n sdes (numbered 1 . . . n). The riser's lower end was connected with a semi-rigid spherical joint to the sea bottom, and the upper end with an articulated ball joint to the base (platform or vessel). The movement of each rfe is described by the vector components: where: (Fig. 2), determining the element's orientation, with torsional deflection neglected.
The equations of individual elements' motion, derived from the Lagrange equations, employ the homogeneous transformation method [24,25]. The Rigid Finite Element Method's modification is described in [8], so the description of the mathematical model which follows is reduced to major determinations and relationships. In paper [8] each rfe has 6 degrees of freedom, while in the approach presented the torsional flexibility is omitted. The coordinates of the point defined in the local coordinate system {i} combined with rfe i, may be expressed in the inertial system according to the following formula: where:  components of the vector q i . The kinetic energy of rfe i can be defined with following dependence [24]: where: H i = mi r i r T i dm i . The Lagrange operators from kinetic energy of rfe i may be presented as follows:  where: The high number of elements present in (4) is at the same time equal to nil. Therefore, in this paper, a similar method was used to that referred to in [9]: products and traces of matrices, referred to in (4), were calculated analytically.
In order to connect rfes in kinematic chain reflecting considered raiser, the constrain equation is formulated and reaction forces in points A i are introduced in equations of motion. Following procedure described in [8], the generalised forces acting on ith rfe are as follows ( Fig. 3): The constraint equations are involved in acceleration form: (q n ,q n ) are presented in paper [8]. The potential energy of gravity forces of refs can be exposed in the form: where m is mass of ith rigid element, g is gravity acceleration,y c,i is y coordinate of mass centre of i-th element. From (7) it is obtained: The elastic deformation energy sdei is defined by expression: coefficients c i were determined from: where: E r -stiffness modulus of riser material, I rmoment of inertia of riser cross section, L-length of the original element, the bendability of which is reflected by sdei. Since α i and angles ψ i , θ i are interrelated according to: then: Provided that the angles ψ i and θ i are small, formulas (5) are reduced to well-known expressions: 2) The water environment impact forces can be introduced using Morison equations [26,27]. The forces caused by flow inside the element are disregarded. Then into account are taken the following: -buoyancy force, -viscous resistance forces, -inertia force.
After calculations presented in details in paper [8], the equations of motion of the system considered can be written in the form: where: Q we i is vector of all generalised forces arising from water environment. In general matrix form the above equation can be written as: the solution of which enables determination of forces in sdes connection, (b) generalised acceleration formulas: F(t, q,q). (16.2) With vector R known, the longitudinal and transverse forces may be determined in rfe connection with adjacent elements and the base in point A n+1 . The system motion equations were integrated by the Runge-Kutta IVth order method. To stabilise the constrain equations, the Baumgart method was applied.

Model validation, natural vibration frequency determination
Subject to numerical analysis was one of the risers presented in [9]. Listed in Table 1 are the riser's geometric and material parameters and the parameters needed to determine the water impact forces. Table 2 shows how increasing the number of elements affects the minimum and maximum forces at point A n+1 , and the maximum and minimum bending moments for s = 70 m. These values were obtained assuming a forced movement of the base in the direction of axis x according to: where: a = 2 m, ω = 2π T f , T f = 12 s. The results shown in [9] were obtained by the finite difference method.
The difference between the values obtained in [9] and from the present model does not exceed 5% already at n = 60.
The impact of number n on the results accuracy was also analysed for the riser's natural vibration frequency calculation. The frequencies were calculated using Normal drag coefficient C n 1 Tangential drag coefficient C t 0 Fluid inside riser Mass density p F 200 kg/m 3 the procedure described in [28]. With the assumption that:   and the approximation of the right side of (16.2) by: and with the assumption thaẗ then after simple transformations the problem of eigenvalues and eigenvectors takes the form: Gradients matrix C may be calculated by the finite difference method. In this study the central five-point differences were employed. Table 3 shows the impact of number n on ω for a riser with the parameters from Table 1, for H = 1800 m. It is seen that already the differences between n = 160 and n = 40 do not exceed 1%. Table 4 shows the vibration frequencies of the same riser for different parameters of H . This parameter was so chosen, in order to eliminate the riser's contact with the sea bottom    A 0 ), which enabled elimination of the "sea bed interaction" problem. The vibration frequencies of the vertical riser (without bending) from Fig. 1b are shown in Table 5.
Comparison of ω from Tables 3 and 5 shows significant differences between them. The frequencies of the vertical riser's vibration in planes XY and Y Z are identical due to the circular cross section.

Dynamic optimisation task
It was assumed in the study that the movement of the riser's upper end A n+1 takes the cyclic form with periodic time T f = 12 s and amplitudes a = 3, 4 and 5 m. Figure 4 shows the motion forced in cycles 1 and 2. The deformation in the initial stage (t < 3 s) is due to the assumption of uniform initial conditions for vectorq| t=0 .
The bending moment course for a = 4 m at points s = 70 m and s = 170 m is shown in Fig. 5. For s ≈  Fig. 6. Figure 7, alternatively, shows the riser's positions at several selected moments in time during the forced motion's first and second cycles. The analysis of the above graphs shows that bending moments M and top tension forces T L change by approx. 10% relative to the baseline. Therefore, the dynamic optimisation task was formulated as: Sought for is the function: such as that functional where: p min i , p max i − minimum and maximum permissible limits of parameters p i , p i = y A i+1 (t i ), t i − equidistant points in range 0; T S for i = 0, 1, . . . , m. The downhill simplex method [29] has been applied in order to solve the optimisation task. Optimisation procedure used for calculating the values of function y A n+1 (p) is shown in Fig. 8. This procedure has been used for both cases considering stabilisation of the bending moment and torsional force. The criterion for stopping the optimisation process is the maximum acceptable difference between the following values:  2) or the number of optimisation steps j max . A characteristic feature of the nonlinear optimisation task so formulated is that for each combination of parameters p 1 ÷ p m−1 the system motion equations (23) needed to be integrated to calculate the value of functional . Hence, the numerical effectiveness of the methods of riser discretisation and of system motion equations integration is so important. Assumed in the calculations, the results of which are shown hereinafter, was n = 40. Under this assumption, one step of the optimisation procedure with integration step h = 0.02 s can be completed on a mid-range PC computer in 2 s or less.

Optimisation results
Numerical calculations involving the minimisation of functionals (23.1) and (23.2) were performed for amplitudes a = 3, 4 and 5 m of the riser's upper end horizontal displacement in XZ plane. Time interval 0; T S , where T S = 24 s, was divided into m = 10 subintervals, and parameters p i = y A i+1 (t i ), at the start and end of each sub-interval, were chosen considering In all analysed cases the riser top end's vertical displacements needed to stabilise the moments at the selected point, or the tension force, did not exceed 1.1 m. The choice of a better stabilisation requires direct comparison of the resulting courses. Compared below are the two ways considered to stabilise movement amplitude a = 4 m. It may be found by analysing courses M (s, t) in Fig. 12a that the bending moment at point s = 70 m is better stabilised in case O1. The difference M max − M min for the results after optimisation in case O1 is 1.33 kN, and in case O2 is 5.4 kN, whereas the resulting T L (t) courses differ only slightly (Fig. 12c). Similar findings were observed when analysing the other courses for amplitudes a = 3 and 5 m.  The riser's top end positions before and after the optimisation is shown in Fig. 13.
In case O1, when the bending moment was stabilised at the indicated point in the two motion cycles, similar trajectories were obtained of the riser top end's movement projected on XY and YZ planes. It is different when the riser's tension force is stabilised, i.e. in case O2, when the resulting trajectories in the first and second cycles vary considerably.
Both the stabilisation methods herein reported have effect of the reduction of the bending moment and the tension force alike. Shown in Table 6 are the calculated differences between the maximum and minimum M and T L .
In the first optimisation case (O1) differences M max − M min did not exceed 1.6 kN, which at the initial difference of approx. 34 kN, should be considered a very good result. Better results in the top tension minimising were obtained in case O1, where they were slightly better than in case O2. Difference T max L − T min L after the optimisation is on average higher by 5.2 kN in case O1 than in case O2.

Conclusions
The riser dynamics model formulated in this study features high numerical effectiveness. This makes it suitable for solving certain optimisation tasks. They consist in such a selection of riser top end's vertical displacements that despite its periodic movement in the horizontal plane the stresses from bending moments at certain points are stabilised. It is also indicated how the distance riser's top and bottom end effects its curvature and changes its natural vibration frequency. The optimisation results shown in this paper are fragmentary. But it is possible to carry out more calculations, which could provide a learning set for ANN (artificial neural network), as shown in [7]. Because AHC (active heave