AN ACCURATE AND COMPUTATIONALLY EFFICIENT SMALL-SCALE NONLINEAR FEA OF FLEXIBLE RISERS

This paper presents a highly efficient small-scale, detailed finite-element modelling method for flexible risers which can be effectively implemented in a fully-nested (FE) multiscale analysis based on computational homogenization. By exploiting cyclic symmetry and applying periodic boundary conditions, only a small fraction of a flexible pipe is used for a detailed nonlinear finiteelement analysis at the small scale. In this model, using three-dimensional elements, all layer components are individually modelled and a surface-to-surface frictional contact model is used to simulate their interaction. The approach is applied on a 5-layered pipe made of inner, outer and intermediate polymer layers and two intermediate armour layers, each made of 40 steel tendons. The capability of the method in capturing the detailed nonlinear effects and the great advantage in terms of significant CPU time saving are demonstrated by comparing the results obtained on elements of pipe of different lengths, equal to one pitch length Lp as well as Lp/5, Lp/20 and Lp/40.


INTRODUCTION
Unbonded flexible risers have become the main means for transporting oil and gas between the seabed and surface in ultra-deep waters. They consist of several polymer and steel layers that can move internally relative to each other. This gives them low bending stiffness and makes them highly valuable tools for subsea oil and gas companies. Their ability of withstanding large displacements and rotations makes them ideal for floating platforms. Due to the complex internal dynamics of flexible risers related to the possible interlayer slip, however, the conventional stress prediction and fatigue analysis tools based on simplified analytical formulations and linear methods have a limited 2 degree of accuracy. In many problems of very significant industrial interest, including but not limited to fatigue-life prediction and life extension for existing risers or forensic analyses after major accidents, sufficient accuracy can only be obtained by the use of models that properly take into account contact and friction between layers and how these are related to internal and external pressure, bending and torsion of individual tendons, large displacements and rotations [1].
Major established industries in the field rely on computational mechanics software to analyze and design flexible pipes. In few cases the finite difference method (FDM) is used to discretize the model in space [2]. Here, the focus will be on the finite element method (FEM), which is the most widely used approach as it is capable of handling geometrically complicated domains, a variety of boundary conditions, nonlinearities, and coupled phenomena that are common in flexible risers.
A number of authors have developed various finite-element (FE) modelling approaches for flexible pipes. de Sousa et al. [3] and Merino et al. [4] describe a model using a commercial FE package with concentric solid layers discretised with thin-walled 4-noded shells. The carcass and pressure armour layers are modelled as equivalent cylindrical layers with orthotropic properties, for which analytical derivations are presented. Tendons are modelled as three-dimensional Euler-Bernoulli beams, with principal axes (for moment calculations) in the pipe radial direction, and a penalty method for contact constraint enforcement is used. A similar approach using a general purpose FE program is adopted by Le Corre and Probyn [5]. In this model of a single-core umbilical, three concentric sheaths are modelled as cylindrical shells. The annulus between each pair of sheaths contains helically wound tubes and cables modelled as circular-section beams.
Simulations are carried out using an explicit-dynamics solution procedure, with a general contact algorithm to include all parts with a friction coefficient derived from tests. Although FE models can account for the complex internal structure of flexible risers, their computational requirements limit their applicability to just a few meters in length at most. So, a more efficient methodology with lower computational cost is required to bridge the gap between nonlinear dynamic simulations at the large scale and detailed finite element models at the small scale.
One approach to reduce computational cost of the analysis of flexible risers is to develop constitutive laws for large-scale beam models, which link generalised stresses and strains to model the hysteresis loops occurring for flexible pipes subjected to cyclic loading, as shown by Tan et al. [6]. Using these constitutive laws a nonlinear bending stiffness of the pipe, which captures the transition from no-slipping to slipping between layers, can be included in the analysis. Due to the close analogy between the aforementioned hysteretic response and the elasto-plastic behaviour of metal beams, this hysteretic response can be modelled as a rate-independent elasto-plastic relationship between generalised strains and stresses. Following this approach, Saevik [7,8] presents a FE model for predicting stresses due to axisymmetric loads and bending loads in flexible pipes. 3 The interlayer stick-slip behaviour due to friction is taken into account by formulating a constitutive relation based on a plastic beam model, with nonlinear stiffness derived from an analytical formulation in terms of the moment resultants and the wire slips. Experimental studies on the model using strain measurements showed that the numerical distribution of the longitudinal stress predicted by the method is in reasonable agreement with the measured data. A similar approach is used by Alfano et al. [9], building on the analogy between frictional slipping between different layers of a flexible riser and frictional slipping between micro-planes of a continuum medium in non-associative elasto-plasticity. In this way, a linear elastic relationship was used for the initial response, in which no slip occurs, and a non-associative rule with linear kinematic hardening was then introduced to model the full-slip phase. Using this approach, the above authors conducted several analyses using a detailed FE model at the small scale to estimate the input parameters of the beam-constitutive model at the large scale.
Key to accurate solution of the small scale FE problem is the use of suitable boundary conditions. Various boundary conditions can be used, including zero fluctuations over the whole model, uniform displacement, uniform traction and periodic boundary conditions [10]. The periodic boundary conditions are the most effective and accurate for most cases involving a periodic microstructure or when the microstructure is not periodic but the small scale model is sufficiently statistically representative [11,12]. Leroy et al. [13], in their FE model of a flexible pipe, assume periodic solutions (given constant curvature in the pipe) and determine an analytical solution of equilibrium of wires on a torus (the bent pipe). The 3D periodic model consisted of a single layer of helical wires, with all internal and external layers represented by rigid kernels. They also developed a detailed FE model of a full pitch length of the pipe, solved using explicit dynamics. Crossvalidation between the models was carried out for cyclic bending and a comparison of the axial stress distributions showed good correlation for the inner armour layer.
In their analysis of one-pitch long segment of a flexible riser using a detailed FE model, which includes all layers in frictional contact between each other, Edmans et al. [14,15] studied the influence of boundary conditions by considering two cases. In one case fixed in-plane boundary conditions are used, whereby the displacements of the nodes of two end-cross sections are prescribed so that the two cross sections undergo a rigid relative displacement and rotation. In a second case, periodic boundary conditions are used, whereby, in addition to the relative displacement and rotation between the two end cross sections, additional displacement (fluctuations) are allowed for all nodes of the cross section, with the constraint that the fluctuations must be the same for two nodes of either end cross section with the same position in the cross section. The comparison between the two types of boundary conditions, in a case in which the onepitch long riser segment is subject to prescribed bending, shows that fixed in-plane boundary 4 conditions result in a much stiffer response [14] and significant spurious edge effects, with stress concentrations building close to the two end sections [15], whereas periodic boundary conditions lead to a much more realistic stress distribution that is uniform across the longitudinal direction.
One major challenge in using the constitutive law based on the non-associative elasto-plasticity analogy is the determination of the parameters of the constitutive law to bridge the small scale of the detailed FE simulations with the large scale of the model accurately. An alternative approach which does not have this limitation is a fully-nested multi-scale procedure [16], currently in widespread use for the modelling of composite materials. With this method, at each integration point (i.e. cross section) of the large-scale beam model, the stress resultants corresponding to assigned generalised strains are determined through the solution of the small-scale FE problem.
This requires recasting the computational homogenisation problem in a more general theory which can link different structural models at different scales [17,18].
When FE models are used at both scales, the fully nested procedure is also known as the FE 2 method. This name very clearly highlights the significant computational cost associated with fully nested computational homogenisation, because, in an implicit incremental solution strategy, the nonlinear small-scale FE model is to be solved at each integration point of the large-scale model, at each equilibrium iteration conducted within each increment of the analysis. For this reason, a fully nested approach is often used only for 'hot-spots', i.e. critical areas of the structure where significant accuracy is needed [16].
On the other hand, even a long riser can be modelled at the large scale with reasonably good accuracy using hundreds or few thousands nodes and integration points. Therefore, the solution of the linear system at each equilibrium equation is normally not an issue. The assembly of the residual vector and tangent stiffness matrix is where the small-scale detailed FE models have to be solved at each integration point, but this is the 'perfectly scalable' part of the process if parallel computing is used, because there is no need for any exchange of information between small-scale analyses at different integration points. Therefore, a fully nested procedure can be feasible if the computational cost of the small-scale model is made sufficiently small so that it can be effectively run within one single node of a cluster in a limited amount of time.
Hence, building on some preliminary work done for an extremely simplified riser made of only two polymer layers and one armour layer [19], this paper describes an efficient modelling approach for the small-scale analysis of a more realistic flexible riser made of 3 polymer layers and 2 armour layers. It also explains how cyclic symmetry of the riser can be exploited by writing periodic boundary conditions based on the multiscale theory derived by Edmans et al. [18] for the case when different structural models are used at different scales. This results in periodic boundary conditions, 5 written in terms of dummy nodes, which are different from those used in the aforementioned work by Leroy et al. [13].
We consider a riser for which the pitch length of the tendons is the same in both the inner and the outer armour layers, which is a condition very close to the real design of risers. This means that, apart from the carcass and the pressure armour layer, if the latter is present, denoting by the pitch length of each tendon and by the number of tendons, the geometries of two cross sections at a distance / , or multiple of it, is the same. The carcass and the pressure armour layer are typically made of a single tendon wound at a small pitch with an interlocking mechanism that prevents unwinding. In general, these layers do not present the same cyclic symmetry. However, while their radial stiffness is significant, their extensional and bending stiffness can be normally considered negligible. Furthermore, it is widely accepted that they can be effectively modelled as a continuum pipe made of orthotropic material, with the longitudinal axis as one of the material directions.
Therefore, with this assumption, by definition they do not violate the cyclic symmetry of the riser.
On the other hand, one problem that has not been considered in previous work is that uniqueness of the solution cannot guaranteed in the presence of friction [20]. The probability of bifurcations is expected to increase with the number of layers because of the increase in the surfaces where frictional slips occur. To address this issue, various riser models with different lengths starting from the smallest repeat of unit of length equal to / to a model pitch length are used.
The results have been compared and discussed.
The structure of the paper is as follows. Section 2 describes how the fully-nested (FE 2 ) computational homogenization method can be applied to flexible risers using the extended theory presented by Edmans et al. [18]. The FE models for the small-scale analysis of risers of different lengths are described in Section 3; in particular, the implementation of periodic boundary conditions is discussed in detail in Section 3.1. In Section 4, the differences in numerical results and CPU time obtained by using models of different lengths are reported and discussed to evaluate the accuracy and the computational saving entailed by the use of the smaller and smaller 'slices' of repeating units of riser. The effectiveness of the use of periodic boundary conditions is also assessed further, by comparing results obtained with periodic and fixed in-plane boundary conditions. Finally, conclusive remarks are made and future work discussed in Section 5.

FULLY NESTED MULTISCALE ANALYSIS OF RISERS
A fully-nested computational homogenization scheme is essentially based on the construction of a micro-scale (or more generally small-scale) boundary-value problem (BVP) at each integration 6 point of a macro (or large-scale) model. In particular, this 'micro' problem is defined for a suitably defined representative volume element (RVE) of the micro scale and is solved numerically to determine the constitutive response of the material at each integration point.
For the flexible risers considered in this paper, a 3D continuum model is used at the small scale, which is to be linked to a large-scale beam model, where generalised strains and stress resultants are employed. This means that different structural models are used at different scales, which makes classical computational homogenisation, e.g. [16], not directly applicable as shown in detail by Edmans et al. [18], where an extension of that theory has been formulated. Referring to the original paper for the details of the derivation, a geometrically non-linear formulation is assumed at the large scale: in particular, it is assumed that displacements and rotations are large, while macro strains are small enough so that a geometrically linear formulation can be adopted at the small scale. The small-scale problem is then written in its general form as follows: In this equation, ̅ is a linear operator that translates the large-scale strain at the integration point into a corresponding small-scale displacement field , is a suitably defined operator that extracts the appropriate boundary values of a small-scale displacement field, so that Equation (1) Given a large-scale strain , the small-scale problem consists of finding a small-scale displacement field solution of the variational problem (1). In practice, for each RVE a dummy reference node R is introduced, whose degrees of freedom are collected in a vector = , which therefore represents the components of the macro strain. Therefore Equation (1) becomes: Once the micro-stresses are computed in (2), by using a 'generalised Hill condition' (GHC) enforcing the equality of the virtual works at the two scales [18], the macro-stress is computed via the following variational equation: which, in practice, means that is simply obtained as the nodal reaction vector at the dummy reference node.
The practical implementation of Equation (2) in the small-scale FE model of a segment of riser is described in the next section.

SMALL-SCALE FINITE-ELEMENT MODELS OF SEGMENTS OF RISER
A simplified 5-layer flexible pipe, made of three polymer layers and two armour layers, was considered. Both the inner and the outer armour layers are made of 40 steel tendons, with rectangular cross section, which are wound with the same pitch length equal to 320mm. The model is created using the FE package ABAQUS, version 6.13.1. All components are modelled with fully-integrated 8-noded 3D solid elements with incompatible strains [21], with surface-to-surface frictional contact between all components.
As previously discussed, it is widely accepted that the use of periodic boundary conditions in multiscale computational homogenization provides the most accurate results, at least at sufficient distance from the real boundary of the structure. On the other hand, the solution for a segment of riser whose length is any multiple of For each model, Table 1 reports the length and the number of nodes. In Table 2, for each layer of the smallest model, the number of elements in the axial and circumferential directions and the total number of elements are given. For each layer, the inner and outer radii, 0 and 1 , the material, the Young's modulus and Poisson's ratio, are reported in Table 2. The cross section of the model is 8 shown in Figure 2(a) and tendon layers arrangements for the smallest model are shown in Figure   2(b). The  The study in this paper focuses on cases in which the segments of pipes are subject to bending, as well as internal and external pressure, the latter being balanced to produce a relatively small and outward radial displacement of the inner surface of the inner polymer layer (inner liner). While a carcass layer is present in most designs of flexible risers, it is typically not a leak proof layer, whereby the internal pressure is applied to the inner liner. As a result, the carcass would remain undeformed under the load cases considered in our analysis and for this reason it has been not 9 considered in the models. Furthermore, we did not consider the presence of a pressure armour layer because this is not present in all risers and we wanted to limit the complexity of the analysis, by focusing on proof-of-concept, yet realistic case studies.
A fully-implicit nonlinear static analysis based on the Newton-Raphson method is used to solve the models. The analyses were carried out in parallel on a computer cluster with two dual-core 1.8 GHz processors (32 processors in total) using 8MB of RAM. Whilst the run time for the smallest model is only few minutes, it takes over 26 hours to complete the analysis of a model with one pitch length.   The contact pressure, , between two surfaces at a point is the sum of two components, and [22]: where is a function of the 'over closure', ℎ, of the surfaces (the interpenetration of the surfaces) and is regularization term.
A 'hard' contact model in the normal direction was used whereby, if ℎ < 0, there is no contact between the surfaces = 0. If ℎ ≥ 0, the surfaces are in contact, see Figure 3 and is then calculated as a Lagrange multiplier, using a penalty method [22]. A Coulomb friction model is used, assuming for the friction coefficient a value of 0.16.
The regularization term, , is used to avoid convergence difficulties arising due to the sudden violation of contact constraints. It is obtained as a 'viscous' pressure that is a function of the normal relative velocity, ℎ̇ , at which the surfaces approach or separate from each other, and of a damping coefficient , which in turn is a function of ℎ: The damping coefficient, , was generated automatically by the program and then scaled down by a factor ranging between 1 and 10. It was checked that varying this scaling factor does not influence the results.

IMPLEMENTATION OF PERIODIC BOUNDARY CONDITIONS
Due to cyclic symmetry and since the length of each model is a multiple of / , for each one of the considered models there is a one-to-one correspondence between the nodes of the two end cross sections that have the same position in the plane of the un-deformed cross-section. Therefore, denoting by the number of nodes on either end cross section, pairs of nodes are defined by this correspondence. Following [17,18], for each one of these pairs of nodes, a new 'dummy' projected node is introduced on a plane which is parallel to the end cross section in their undeformed configuration. This is shown in Figure 4, where, for the generic pair of nodes, and indicate the corresponding nodes on the left-hand and right-hand end cross section, respectively, and denotes the corresponding projected node. 13 A set of linear constraint equations was then generated, relating the degrees of freedom of each pair of nodes on the two end cross section to those of the corresponding projected node. The link is enforced for all displacement components and rotation components of a node , as follows: It is worth noting that, in ABAQUS, a linear constraint equations between three degrees of freedom, 1 , 2 and 3 , are written in the form: where 1 , 2 and 3 are real coefficients. Therefore, for example, Equations (6) are implemented by choosing 1 = , 2 = and 3 = , and by setting 1 = 1, 2 = −1 and 3 = −1, for = 1,2,3. An analogous procedure is used for Equations (7).
The displacement and rotation vectors, and of each dummy projected node are rigidly constrained to the displacement and rotation vectors, and , of a projected reference point at their center, using the following rigid-body constraint equations: where indicates the position vector of the projected node with respect to the reference point and × denotes the standard cross product of two vectors.
In this way, displacements and rotations of the reference point correspond to generalised strain components in the model [17,18]. For example, prescribing a rotation of the reference point about one axis in the cross section is equivalent to prescribing periodic boundary conditions corresponding to a bending curvature equal to / , being the length of the model.
In other words, the vector appearing in Equation (2) and, in its variation, in Equation (3), is here given by:

NUMERICAL RESULTS
The analyses were conducted by applying internal and external pressure in a first step, after which one symmetric cyclic history of bending curvature was prescribed, the maximum and    The key role played by the periodic boundary conditions is shown in the results reported in Figures 9-11, which show the contour plots of the stress for the armour layer and for the outer and 17 inner layers, respectively. With fixed in-plane boundary conditions, tendons cannot slip or rotate with respect to the polymer layers at the two ends of the model. This produces reactive bending moments on the tendons resulting in spurious stress concentrations close to the end sections, which make the stress profile vary between the end sections and the mid-section of the riser segment, as can be appreciate in Figures 9(a), 10(a) and 11(a). Instead, with periodic boundary conditions, the stress profiles in Figures 9(b), 10(b) and 11(b) is practically the same at each cross section.
Figures 12 show the difference in the bending moment vs bending curvatures for the longest model. In Figure 13, it can be seen that fixed in-plane boundary conditions result in a much stiffer response and significant spurious edge effects, with stress concentrations building close to the two end sections. The reason is that, since the tendons cannot slip or rotate with respect to the polymer layers at the two ends of the model, they do not have enough length between the two fixed ends to deform in the way they do in reality. Since edge effects are felt in a volume whose size is independent on the model length, the shorter the model the higher is the stiffness when fixed inplane boundary conditions are used. Instead, the models with periodic boundary condition capture the hysteretic response even with the smallest unit. Figure 14 highlights the huge saving in CPU time, from 1600 minutes for the longest model of one pitch length (considered in the past in [14,15]) to only few minutes for the smallest model. This confirms that the latter should be used in a nested multi-scale strategy.

CONCLUSIONS
This paper has presented an approach to minimize the computational cost of a detailed nonlinear Further analyses reported confirm the importance of using periodic boundary conditions, particularly for the shorter length. 22 It is worth noting that precise cyclic symmetry implies lay angles which are almost but not precisely equal and opposite in the two armour layers, because, for the two armour layers, the pitch lengths are the same but their radii are different. More generally, this issue is related to the necessity to balance the torsional response of the two armour layers to avoid, or at least minimise, coupling between torsion and axial tension. The axial-torsional balance is normally obtained with a suitable combined choice of the lay angles, that often (but not always) are assumed to be equal and opposite, as well as the average radii of the armour layers, the cross sectional area of the tendons and their number. The cyclic symmetry considered in our work requires the same pitch length and the same number of tendons. Assuming the radii are dictated by the thickness of the anti-wear polymer layer, torsional balance can be still achieved by proper sizing of the tendon cross sections. For these cases, exact cyclic symmetry can be achieved in the design and the methodology presented in this paper can be applied.
On the other hand, some designs used in the industry do not perfectly satisfy cyclic symmetry to achieve axial-torsional decoupling. Therefore, future work will include the study of how to best approximate the small-scale problem with the assumption of cyclic symmetry in cases where cyclic symmetry is not exact in the real geometry.
Furthermore, while some successful experimental validation of a modelling approach that is similar to the one used in this paper has been presented by Leroy et al. [13], further experimental validation of our proposed numerical model will be an important part of the future developments of this research.