- Email: [email protected]
On the multi-scale computation of un-bonded flexible risers
Engineering Structures 32 (2010) 2287–2299
Contents lists available at ScienceDirect
Engineering Structures journal homepage: www.elsevier.com/locate/engstruct
On the multi-scale computation of un-bonded flexible risers A. Bahtui a,∗ , G. Alfano a , H. Bahai a , S.A. Hosseini-Kordkheili b a
School of Engineering and Design, Brunel University, London, United Kingdom
b
Sharif University of Technology, Aerospace Engineering Department, Tehran, Iran
article
info
Article history: Received 5 August 2009 Received in revised form 10 February 2010 Accepted 1 April 2010 Available online 11 May 2010 Keywords: Non-associative plasticity Frictional slipping Parameter identification Generalized finite-element formulation Large-displacements
abstract The purpose of this paper is to model the detailed effects of interactions that take place between components of un-bonded flexible risers, and to study the three-dimensional motion responses of risers when subjected to axial loads, bending moments, and internal and external pressures. A constitutive law for un-bonded flexible risers is proposed and a procedure for the identification of the related input parameters is developed using a multi-scale approach. A generalized finite element structural model based on the Euler–Bernoulli beam theory is developed in which the constitutive law is embedded. The beam theory is enhanced by the addition of suitable pressure terms to the generalized stresses to account for the internal and external pressures, and it can therefore be efficiently used for large-scale analyses. The developed nonlinear relationship between generalized stresses and strains in the beam is based 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. The merit of the present constitutive law lies in its ability to capture many important aspects of the structural response of risers, including the energy dissipation due to frictional slip between layers and the hysteresis response. Crown Copyright © 2010 Published by Elsevier Ltd. All rights reserved.
1. Introduction Multilayered un-bonded flexible risers are slender marine structures that are widely used in offshore production to convey oil between the well-head and the offshore rig. The main advantage of flexible risers when compared to rigid steel risers is their much lower bending stiffness, leading to smaller radii of curvature with the same pressure capacity, due to the complex makeup of flexible risers [1]. This, in turn results in an increased ability to undergo large deformations under loads induced by the ocean current and wave, vortex-induced vibrations, the motion of the floating vessel and the installation procedure [1]. The detailed layout and description of layers and the behavior of un-bonded flexible risers are addressed in the literature [2–8]. A typical assumption made in the large-scale analyses of flexible risers is to replace all layers of a riser with an equivalent linearly elastic beam having a constant tangent stiffness matrix. The major advantage given by this assumption is that it makes the computational solution time significantly shorter, and renders it suitable for day-to-day industrial projects. A static analysis procedure for the numerical determination of nonlinear static equilibrium configurations of deep-ocean risers
∗
Corresponding author. Tel.: +44 7828487635. E-mail addresses: [email protected] (A. Bahtui), [email protected] (G. Alfano), [email protected] (H. Bahai), [email protected] (S.A. Hosseini-Kordkheili).
was developed by Felippa and Chung [9]. The riser was modeled by three-dimensional beam finite elements which include axial, bending, and torsional deformations. They extended their model by taking the deformations coupled through geometrically nonlinear effects [10]. The resulting tangent stiffness matrix included three contributions identified as linear, geometric (initial stress) and initial displacement stiffness matrices. For the solution, a combination of load–parameter incrementation, state updating of fluid properties and corrective Newton–Raphson iteration was used. The detailed implementation of their procedure is presented in Felippa and Chung [11] in which they discussed the importance of drag force along the pipe, varying significantly along the depth. The authors are of the opinion that the sensitivity of the results to environment characterization necessitates a review of many modeling and analysis practices. McNamara and Lane [12] studied the two-dimensional response of the linear and nonlinear static and dynamic motions of offshore systems such as risers and single-leg mooring towers. Their proposed method was based on the finite element approach using convected coordinates for arbitrary large rotations and includes terms due to loads such as buoyancy, gravity, random waves, currents, ship motions and Morison’s equation. The same authors extended their work by developing a three-dimensional computational dynamic analysis of deep water multi-line flexible risers in the frequency domain [13]. O’Brien et al. [14] presented the three-dimensional finite element modeling of marine risers, pipelines and offshore loading towers based on the separation of the rigid body motions and deformations of elements under
0141-0296/$ – see front matter Crown Copyright © 2010 Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.engstruct.2010.04.003
2288
A. Bahtui et al. / Engineering Structures 32 (2010) 2287–2299
conditions of finite rotations. This paper treats risers as a homogeneous material and includes all nonlinear effects including geometry changes, bending-axial and bending-torsional coupling and follower forces and pressures. The authors preceded their work by dividing the riser into a series of catenaries where the buoyancy modules were described as inverted catenaries [15]. Hoffman et al. [16] reviewed the design technique of deep and shallow water marine riser systems as well as their dynamic analytical and numerical analysis and the nonlinearities arising from hydrodynamic loading and dynamic boundary conditions. Their paper contains design methodology criteria, parameters and procedures of flexible riser systems while treating the riser as one homogeneous material layer. Atadan et al. studied the dynamic three-dimensional response of risers in the presence of ocean waves and ocean currents undergoing large deflections and rotations. They included shear effects based on the nonlinear elastic theory in their formulation. It was concluded that the length of the riser is the most important parameter which affects the deflections of the marine-riser [17]. Chai et al. proposed a three-dimensional lumped-mass formulation for the riser structure analysis which is capable of handling irregular seabed interaction. They adopted a simplified model by replacing the seabed surface with an elastic foundation with independent elastic springs having an arbitrary thickness [18]. Ong and Pellegrino [19] studied the nonlinear dynamic behavior of mooring cables in the frequency domain. They ignored the effects of friction and impact between the cable and the seabed. Their proposed method modeled the time-varying boundary condition at the touchdown by replacing the section of cable interacting with the seabed with a system of coupled linear springs. This work again ignored the hysteresis behavior of cables by assuming an equivalent linearly elastic pipe representing all layers. Zhang et al. [20] discussed analytical tools for improving the performance of flexible pipes. Their work used an equivalent linear bending stiffness which was derived from experimental data to calculate the maximum bending angle range. It contained reports on irregular wave fatigue analysis, collapse, axial compression and bird-caging for riser systems. They observed that the combined bending, axial compression and torsion could lead to the tendon being separated from the cylinder in a helix layer, and may lead to out of plane buckling. Yazdchi and Crisfield used a simple two-dimensional lowerorder beam element formulation for static nonlinear analysis of riser structures including the effects of buoyancy, steady-state current loading and riser top-tension. They assumed a linearly elastic material property for the riser by employing a constant modulus of elasticity, and studied the static behavior of flexible pipelines and marine risers, using finite elements developed for conventional nonlinear analysis [21]. The same authors continued their research by using a beam finite element formulation, based on the Reissner–Simo beam theory, for the static and dynamic nonlinear analysis of three-dimensional flexible pipes and riser systems in the presence of hydrostatic and hydrodynamic forces. Employing linearly elastic material properties, their work concentrates on the nonlinearities due to the fluid loading and the associated large deformations and considers hydrodynamic forces due to effects such as wave, drag and current action [22]. In all the above models, the assumption of perfect bonding between layers fails to take into consideration energy dissipation due to frictional slip between layers and the hysteresis response of the risers, generating inaccuracy in the results and making the analysis unsuitable for accounting for the hysteresis damping in vibration or fatigue. Some nonlinearities have been considered by Lacarbonara and Pacitti, who presented a geometrically exact formulation for the
dynamic analysis of cables with linearly elastic material behaviors when undergoing axis stretching and flexural curvature. Their formulation was based on nonlinearly viscoelastic constitutive laws for the tension and bending moment with an additional constitutive nonlinearity accounting for the no-compression condition [23]. Some works in the literature consider a bi-linear bending stiffness, taking into account the energy dissipation originating from frictional contact between layers in the bending moment analysis [24,25]. However the bi-linear bending stiffness proposed in these papers does not depend on internal and external pressures. This generates inaccuracies in the results when studying the response of a long riser floating in ultra-deep water. Moreover the study has been performed only on the hysteresis effect in bending analysis whereby hysteresis effects for axial strain or any coupling with radial displacement is ignored. On the other hand, the bi-linear model in these works is based on two different stiffnesses, one when perfect bonding is assumed between layers and one when slippage occurs between layers. The limitations outlined in the above contributions motivated the authors to study the effect of the interaction between layers of un-bonded flexible risers on their hysteresis behavior in more detail. The present work addresses the above issues by combining the accuracy of a detailed finite element model of a short part of riser, with the low computational cost of a three-dimensional Euler–Bernoulli beam element, in the framework of a novel, multiscale approach to conduct nonlinear three-dimensional largescale structural analyses of an entire un-bonded flexible riser. The ‘equivalent elasto-plastic’ constitutive law for the beam cross section proposed by the authors [2], which relates axial strain, torsion and bending curvatures to the conjugate stress resultants is enhanced by including additional degrees of freedom for the beam model which account for axial strain and radial displacement. This law models the small-scale frictional slip occurring between the layers, when a combination of values of the stress resultants and of the internal and external pressures overcomes a given threshold, at the macroscopic level of the cross section. The ‘equivalent yield function’ proposed by the authors [2], which depends on the current values of the stress resultants, on a set of internal history variables and on the values of the external and internal pressures acting on or within the pipe, is enhanced by taking into account the axial stress effects, making the equivalent yield function threedimensional. This last functional relationship simulates the effect that the pressure values have on the inter-layer normal stress, which in turn affects the inter-layer frictional sliding and the overall dissipation. The zero level set of the yield function represents the boundary of a ‘no-slip domain’, that is, the onset of slipping. In formulating this model a linear kinematic hardening law has been employed. The assumptions of linear hardening and of linear elastic behavior within the ‘no-slip domain’ result in a bilinear response, which can provide sufficient approximation of the structural response in many cases of interest. This simplified model entails that the single cross section of the riser can be either in a state in which ‘no-slip’ is found, or in a state in which ‘full-slip’ occurs. To identify the parameters of the beam model to be used in large-scale analyses of the entire riser, ideally one would need a great number of results of experimental tests conducted on a reasonably large model of a flexible pipe using an experimental rig which reproduces the real conditions as closely as possible. This is very difficult and extremely expensive, which is the main motivation for the development of the multi-scale approach proposed in this paper, in which the small-scale detailed FE model is employed as a virtual experimental rig, and the parameters of the large-scale model are then determined by combining engineering assumptions and careful curve fitting of the results of the smallscale and large-scale models.
A. Bahtui et al. / Engineering Structures 32 (2010) 2287–2299
2289
In addition, the internal and external pressures need to be considered in the model. If εr and ur are the average radial strain and the radial displacement of the average radius, respectively, Pε and Pu are the conjugate actions working for εr and ur . For a detailed derivation of these variables please refer to [2]. Hence, in the proposed model N , Mt , Mx , My , Pu and Pε are the generalized stresses. The corresponding work-conjugated generalized strains are εz , φ, χx , χy , ur , εr , where εz , χx , χy and φ are the axial strain, the curvatures about the x and y directions and the torsional curvature. The detailed definition of the foregoing is in Section 3. A compact notation is conveniently introduced denoting by {σ} and {ε} the generalized stresses and strains:
N
Fig. 1. Stress resultants.
A generalized finite element beam model supplemented by the constitutive law is used for large-scale analysis, unlike the detailed finite element model used to formulate it and to determine its parameters. This model is able to capture all the important aspects of the structural response of the riser and is suitable for the study of many practical dynamic problems, such as hysteresis damping in vibration, fatigue or service life analysis. Forthcoming papers will present results related to these situations. This paper is divided into nine sections. In Section 2 the constitutive law for the three-dimensional beam model is formulated. In Section 3 the comprehensive global analysis of the riser is formulated. In Section 4 the tangent operator matrix is derived. In Section 5 the detailed finite element model, mainly based on earlier work by the authors [4,7], is described. In Section 6 the key issue of the identification of the parameters of the beam constitutive model is addressed. In Section 7 one validation case is presented and discussed in order to determine the extent to which the results of the various detailed finite-element simulations for a short riser, subject to internal and external uniform pressure and uniform cyclic loading, agreed with those given by the proposed constitutive law. A practical case of a riser is then presented in Section 8. Results from a horizontal un-bonded flexible riser with a length equal to 500 m subject to top end bending moment due to vessel movements, are compared to a steel rigid riser with the same length and conditions. Finally, in Section 9, conclusions are drawn based on the results presented. 2. Formulation of the constitutive model The constitutive law of an un-bonded flexible riser for an Euler–Bernoulli three-dimensional beam model, proposed in [2], is significantly enhanced in this section to make it more suitable for practical usage. In particular, axial strain and radial displacement are taken into account and the ‘slip-onset’ function is made dependent on both the axial force and the bending moment resultants. It is assumed that locally the pipe develops along a straight line, so that an infinitesimal element of pipe can be represented as in Fig. 1. A local right-handed Cartesian system is introduced with the origin and the x and y axes located at the cross section and the z axis coincident with the centroid axis. 2.1. Generalized stresses and strains The stress resultants, shown in Fig. 1, are the axial force N, the torque Mt , and the bending moments around axes x and y, denoted by Mx and My , respectively.
Mt Mx {σ} = My Pu Pε
ε z
φ χx {ε} = χy .
(1)
ur
εr
2.2. Constitutive law As discussed in Section 1, the constitutive law for the beam model is based on the observation that the local (small-scale) frictional slip between the different layers results in a macroscopic (large-scale) relationship between generalized stresses and strains which has many analogies with the laws of elasto-plasticity. The formulation of a constitutive law for the bending stiffness was presented by the same authors [2] and hence will not be presented here. The derivation of a constitutive law for the axial stiffness behavior is based on a similar procedure, leading to a different ‘slip-onset’ function. Fig. 2 shows the response of the flexible riser under cyclic loading, evaluated using a detailed threedimensional finite element analysis which accurately models the interaction between all layers. The case study of axial compression loading has not been considered in this work. As a consequence, Fig. 2 or all other axial force figures have a minimum value of zero. Besides, the no-slip portion at unloading has approximately same length as the no-slip portion in the initial cycle. For the qualitative analysis in this section, presented in order to justify the assumptions of the proposed constitutive model, it is sufficient to state that the values of the applied loading (internal and external pressures included) are, within a meaningful range, typical of practical cases. For the first monotonic increase of the axial strain, the axial force increases with an almost linear stiffness between points (a) and (b). Then the stiffness gradually decreases between points (b) and (c) to a value which remains approximately constant up to point (d). Upon unloading, the axial stiffness is initially very close to the initial stiffness, between points (d) and (e), and then it rapidly decreases between points (e) and (f) to a value which is constant between points (f) and (a) and approximately equal to that between points (c) and (d). Replacing the rounded parts between points (b) and (c) and between points (e) and (f) with two sharp elbows, a bilinear curve response is obtained with good approximation. Such a bilinear curve represents the same type of response as that obtained using an elasto-plastic model with linear kinematic hardening, which is the model proposed here. In analogy with the elasto-plastic case, the generalized strains are additively decomposed into an elastic-like, i.e. ‘no-slip’, part {εe }, and a plastic-like, i.e. ‘full-slip’, part {εs }.
{ε} = {εe } + {εs } .
(2)
For full details of the formulation of the constitutive law, particularly the ‘slip-onset’ function f and the real-valued ‘slip potential’ g, please refer to [2].
2290
A. Bahtui et al. / Engineering Structures 32 (2010) 2287–2299
denoted by βz , βx , βy , βr and βε , so that the argument of the onsetslip function is given by: N − βz Mx − βx {σ} − {β} = My − βy . P − β u r Pε − βε
(4)
To determine functions f and g and the coefficients of the kinematic hardening matrix, many results from several numerical simulations made with the detailed FE model have been reviewed and the following expression has been found to be sufficiently valid: Fig. 2. Plot of the axial force against the corresponding strain for an un-bonded flexible riser during cyclic loading; the response is evaluated using a detailed finiteelement simulation.
2.2.1. Specialization to the case of torsion-less loading In this section, the constitutive law proposed previously by the same authors [2] is significantly enhanced by considering the very important issue of cyclic axial loading and its coupling effect on the slip-onset function. The number of degrees of freedom of the constitutive model is increased to five by incorporating the axial strain, curvatures, radial displacement and the average radial strain εr at each node. The same reasonable assumption made by the authors [2] that the torque Mt is constant throughout the deformation process is retained here to better focus on the factors which have a major influence on the structural response. Hence, its influence on the determination of functions f and g can be incorporated within the coefficients and then ignored. Furthermore, a linear elastic relationship is assumed between Mt and the conjugated variable φ , so that the hypothesis made that Mt is constant allows us to completely ignore this term. A linear elastic relationship is also assumed between Pu and the conjugated variable ur . On the other hand, the main influence on the response of the flexible pipe at the onset of and during frictional slipping is given by the other term Pε . Considering its corresponding work-conjugated generalized strains εr as a degree of freedom of the constitutive model allows for better modeling of the coupling effects between Pε and N. This enhancement to the previous work of the authors [2] allows us to consider the coupling terms between the axial and radial strains. As for the coupling effects between the axial strain and the curvature, they are due to the changes in the radial contact stresses between layers caused by changes in the axial force due to Poisson effects, which in turn affects friction. Hence, these effects are embedded into the slip-onset function, and thus are not considered in the stiffness matrix. Instead, no elastic coupling exists between these variables and between radial displacement and the curvature, which results in zero values of the corresponding terms in the stiffness matrix. Finally, symmetry of the cross section results in the same flexural response about the x and y axes, with no elastic coupling between the two curvatures, or between the curvatures and the axial strain. Hence, the initial stiffness is represented by:
N D11 Mx My = P u Pε
0 D22 sym.
0 0 D22
D41 0 0 D44
D51 εz 0 χx 0 χy . D54 ur D55 εr
(3)
Note that coupling exists between axial generalized strain and radial strain and displacement terms. The components of the back-stress vector {β} which correspond to the generalized-stress components N, Mx , My , Pu and Pε are
h b (N − βz )2 + c (Mx − βx )2 2 i f ({σ} − {β}) = − (Pε − βε ) − a, + My − βy N − βz > 0 − (N − βz ) − a, N − βz ≤ 0.
(5)
where a, b and c are material parameters to be identified, and symmetry of the cross section has been exploited. The slip-onset surface in this case has two surfaces. This is due to the fact that the riser is not assumed to be capable of tolerating axial compression. This indicates that axial force is non-negative. The slip potential g has to account for the finding that frictional slipping is not accompanied by any significant relative opening between layers. At the small-scale level this means that no significant dilatancy is observed. A suitable expression for the potential g which meets this requirement is as follows:
h b (N − βz )2 + c (Mx − βx )2 2 i g ({σ} − {β}) = + My − βy , N − βz > 0 − (N − βz ) , N − βz ≤ 0
(6)
which results in the following slip rule:
∂g ∂N ∂g ∂M ε˙ sz x ∂g χ˙ sx ∂g {˙εs } = λ˙ = λ˙ {m} ⇒ χ˙ sy = λ˙ ∂ My u˙ ∂ {σ} sr ∂g ε˙ sr ∂ Pu ∂g ∂ Pε 2bN 2cMx λ˙ 2cMy , (N − βz ) > 0 0 0 = −1 0 ˙ 0 , (N − βz ) ≤ 0 λ 0
(7)
0
where ε˙ sz , χ˙ sx , χ˙ sy , u˙ sr and ε˙ sr denote the full-slip components of the axial tension, the bending curvatures about the x and y axes, the radial deformation and of the radial strain. Overlapping the generalized stresses and strains spaces, Fig. 3 describes two crosssectional views (a, b) as well as a three-dimensional view (c) of the above three-dimensional equation and shows that the difference
A. Bahtui et al. / Engineering Structures 32 (2010) 2287–2299
2291
N – βz
Mx – βx m
m n
n
Pε – βε
Pε – βε N
Mx
b
a x 10 4 2 m
N - βz
1.5 n
1 0.5 0 2 x 104
1
6 0 Mx - βx
-2
x 10 6
Pε - βε
0
c Fig. 3. Non-associative behavior of the proposed model.
between functions f and g results in the full-slip strain rate {˙εs } being not normal to the slip-onset surface (a, b), i.e. not parallel to the normal {n}. This is in analogy with the non-associative plasticity models used for many materials exhibiting internal friction for which the normality rule is not satisfied. The following proportional hardening law is used:
βz βx {β} = [H] {εs } ⇒ βy β r βε H11
=
0 H22
sym.
0 0 H22
0 0 0 0
0 εsz 0 χsx 0 χsy . 0 usr 0 εsr
(8)
It is worth observing that because of the slip law in Eq. (7), the radial full-slip strains usr and εsr result to be constantly zero, so that βr and βε also do so. A final remark is that, having the same response about the x and y axes because of symmetry, a simple plane model can be used to study the structural response and to identify the model parameters, as will be done in Sections 3 and 4. Nevertheless, it was preferred to derive all the above equations for the general case in which components of moment and curvature about both axes are non-zero to clarify that the model can be used in a general threedimensional problem. 3. Generalized finite element formulation A common approach to the numerical analysis of riser structures is to use an equivalent modeling scheme to simulate the ex-
istence of the various layers which make up the structure of the flexible riser. This approach fails to account for the interaction between layers, the hysteresis effects due to the occurrence of interlayer slip-stick and also the nonlinear behaviors due to the friction between the un-bonded layers. In this section a generalized nonlinear finite element formulation is developed in order to analyze the un-bonded flexible risers, embedding the constitutive model presented in Section 2. The generalized nonlinear finite element formulation is based on the Euler–Bernoulli beam theory with two extra degrees of freedom in the radial and thickness directions. Referring to Fig. 4, displacement components ui , vi and wi for node i are defined as ui = u0i (z ) − y
∂w0i (z ) ∂v0i (z ) +x ∂z ∂z
(9)
vi = v0i (z ) wi = w0i (z )
(10) (11)
where i = 1, 2 indicates the node numbers, u0i (z ), v0i (z ) and w0i (z ) are the displacement of the centre-line in x, y and z directions, respectively. Also ∂w0i (z )/∂ z and ∂ u0i (z )/∂ z are rotation component about x and y axis for node i. Using Eqs. (9)–(11) the nonlinear strain–displacement relations, based on the von Karman assumption, are
εzz =
∂ u0 1 + ∂z 2
εxx = 0 εxz = 0.
∂w0 ∂v0 + ∂z ∂z
2 −y
∂ 2 w0 ∂ 2 v0 +x 2 2 ∂z ∂z
(12) (13) (14)
To develop a nonlinear finite element model for this Euler– Bernoulli beam element we interpolate the displacement values
2292
A. Bahtui et al. / Engineering Structures 32 (2010) 2287–2299
k+1 [B]tNL
2 0 N l 1 k {q}n+1 = 0
B12 N003 /τnk+1
0
2 02 N5 θ1 n + 1k l2 N005 /τnk+1
0
0
2 0 N2 l 0
0 0 l 2
0
N1
B16
B17
N004 /τnk+1
N006 /τnk+1
0
0
0
l 2
0 0 0
N2
(20)
Box I.
+1 {q}kn+1 is given in Box I where [J] is the Jacobean matrix and [B]ktNL (Eq. (20)), where
B12 = B16 = B17 =
2N30 l2 2N40 l2 2N60 l2
k 0 k 0 k 0 k N30 w1n +1 + 2N4 w2n+1 + 2N5 θ1n+1 + 2N6 θ2n+1 k 0 k 0 k N40 w2n +1 + 2N5 θ1n+1 + 2N6 θ2n+1
(21)
2N5 θ1n+1 + N6 θ2n+1 0 k
0 k
Fig. 4. Euler–Bernoulli beam element with extra degree of freedom to model the radial displacement.
as well as generalized strain components ur and εr in the element field using their nodal values as follows
v = N3 v1 + N4 v2 + N5 θ1 + N6 θ2 w = N3 w1 + N4 w2 + N5 θ1 + N6 θ2 ur = N1 ur1 + N2 ur2 εr = N1 εr1 + N2 εr2
(e)
∗ σg εg · εg dz = {F}ext · q˜ za ∗ εg = [B]NL q˜ =
(16)
where Wi is virtual work, superscript (e) shows the elemental value, subscript g uses to indicate the globalcoordinate system pa ∗ rameter, {F}ext is the external load vector, q˜ is the virtual displacement vector, and [B]NL is the nonlinear strain–displacement matrix. Using Eqs. (16), the incremental finite element equation of equilibrium can be written as
{R} ({q}) = {F}int ({q}) − {F}int = 0
derivation of the consistent tangent operator
d{σ} d{ε}
{ε}=[B]NL {q}kn+1
is
1 −1 k k {q}kn+ +1 = {q}n+1 + [K]tNL {q}n+1
+1 {F}next − {F}int {q}kn+1
(18)
where k is the iteration number, and n stands for increment number, [K]tNL is the nonlinear tangent stiffness matrix, and
Z
[B]TNL V
d {ε}
d {σ} {ε}=[B]NL {q}kn+1
× [B]NL det [J] dV
Alfano et al. [26] presented a general approach to the evaluation of the tangent stiffness matrix for associated rate-independent elasto-plasticity materials. This section presents the derivation of this coefficient which follows a similar approach extending it to the non-associated case of interest here. From a computational point of view, the time domain needs to be subdivided into a number of finite steps. Adopting a fully implicit backward-Euler time-integration scheme, the equations to solve in each step are as follows:
{σ}n+1 = [D] {ε}n+1 − {εs }n+1 ∂ g {ε } {ε } − = 1 λ s n s n +1 ∂ {σ} n+1 {β}n+1 = [H] {εs }n+1 n +1 ≤ 0 1λ ≥ 0 f {σ}n+1 − {β} 1λf {σ}n+1 − {β}n+1 = 0.
(22)
Differentiating Eq. (221 ) with respect to {ε}n+1 , we get
∂ {σ}n+1 ∂ {εs }n+1 = [D] [I] − . ∂ {ε}n+1 ∂ {ε}n+1
(23)
(17)
where {q} is the unknown matrix in the form {q} = (u1 w1 θ1 ur1 u2 w2 θ2 ur2 )T . Implementing Newton–Raphson method to solve Eq. (17) leads to
[K]tNL {q}kn+1 =
where l is the element length, N 0 is first derivative of the shape functions, and N 00 is second derivative of the shape functions. The
4. Consistent tangent operator
zb
2 3/2
(15)
where u1 and u2 are the axial nodal displacements, w1 and w2 are transitional displacements, θ1 and θ2 are nodal rotational components. Also in Eq. (15), N1 and N2 are linear shape functions, while N3 –N6 are cubic shape functions. For this nonlinear problem, the governing equation in rectangular Cartesian coordinates is derived using the principle of minimum total potential energy, i.e. Wi
k 0 k 0 k 0 k τnk+1 = 1 + N30 w1n +1 + N4 w2n+1 + N5 θ1n+1 + N6 θ2n+1
presented in the next section.
u = N 1 u1 + N 2 u2
Z
(19)
From Eq. (222 ) we get the expression
∂ 2 g d {σ}n+1 ∂ {σ}2 n+1 ∂ 2g d {β}n+1 . + 1λ ∂ {σ} ∂ {β} n+1
d {εs }n+1 = d1λ {m} + 1λ
(24)
The assumed definition (Eq. (5)) implies that
∂ 2g ∂ 2g ∂ 2g = = . ∂ {σ} ∂ {β} ∂ {β} ∂ {σ} ∂ {β}2
(25)
Substituting Eq. (25) into Eq. (24), and then substituting the result into Eq. (23) and the derivative of (223 ) provides the following
A. Bahtui et al. / Engineering Structures 32 (2010) 2287–2299
d {ε}
d {σ}
2293
= ([D] − λ [D] [Z] [D])
n +1
−
(31)
[[D] − λ [D] [Z] ([D] + [H])] {m} ⊗ [[D] − λ [D] [Z] ([D] + [H])] {n} [([D] + [H]) − λ ([D] + [H]) [Z] ([D] + [H])] {m} · {n} Box II.
system of equations:
−1 ∂ 2 g ∂ 2 g [D] + 1λ − 1 λ 2 ∂ {σ} ∂ {σ}2 n+1 n +1 2 ∂ 2 g − 1λ ∂ g [H]−1 + 1λ 2 ∂ {σ} n+1 ∂ {σ}2 n+1 d {σ}n+1 d {ε}n+1 − d1λ {m} . × = d {β}n+1 d1λ {m}
(26)
From which we obtain d {σ}n+1 = [D] d {ε}n+1 − 1λ [D] [Q] [F]−1 [D] d {ε}n+1
+ 1λd1λ [D] [Q] [F]−1 ([D] + [H]) {m} − d1λ [D] {m} d {β}n+1 = 1λ [H] [Q] [F]−1 [D] d {ε}n+1
(27)
− 1λd1λ [H] [Q] [F]−1 ([D] + [H]) {m} + d1λ [H] {m} where Fig. 5. Detailed geometry of riser (cross-sectional view).
∂ 2 g [Q] = . ∂ {σ}2 n+1
(28)
The fulfillment of Prager’s consistency condition f implies [27]
= f˙ = 0
d {σ}n+1 − d {β}n+1 · {n} = 0.
(29)
Substituting Eq. (27) into Eq. (29) results in 1λ as below [D] − 1λ ([D] + [H]) [Q] [F]−1 [D] d {ε}n+1 · {n} 1λ = . ([D] + [H]) − 1λ ([D] + [H]) [Q] [F]−1 ([D] + [H]) {m} · {n}
(30)
By substituting Eq. (30) into (271 ) we obtain the equation in Box II, where [Z] = [Q] [F]−1 .
(32)
Substituting matrices and performing cross and dot products yields to the final equation which is presented in Appendix. Fig. 6. Finite-element mesh.
5. Detailed finite element model 5.1. Modeling Several FE simulations for a typical 1.7 m long un-bonded flexible riser have been conducted using the FE code ABAQUS to estimate the parameters of the constitutive law for the beam model. The FE model includes a complex makeup of seven-layers of internal and external plastic sheaths, helical armors, carcass and anti-wear layers and is described in the cross-sectional view of Fig. 5. It has been created by adding one additional layer to the model considered by the authors in previous work [4,5]. Fig. 6 shows the finite element mesh of the riser and indicates the number of nodes and elements used. The full details of the finite element model are available in [2–5] and are not presented here for the sake of brevity.
5.2. Analysis For the identification of the parameters, constructing the constitutive model, six sets of analyses are performed and are summarized in Table 1. The second column indicates whether the analysis is performed for several different values of input pressures Pε , referring to each one as an individual ‘load case’. In these analyses (5 and 6) the simulation starts with an initial ‘pressure load step’ which accounts for the average pressures Pε . The pressure step is then followed by the load(s) from column three. The different considered values of Pε for each axial load case are reported in Table 2. The third column, ‘Loadings’, shows the type of load which is applied at each analysis. The fourth column determines the additional boundary conditions required for the derivation of a target
2294
A. Bahtui et al. / Engineering Structures 32 (2010) 2287–2299
Table 1 Detailed information about individual analyses. No.
Dependency on Pε
Loadingsa
BC at the TRP
Equation to be solved
Target parameter(s)
1 2 3 4 5 6
– – – No Yesb Yesb
Only Pε Only Pε Only Pε Only Pu N Mx
εz = 0, ur = 0
D55 εr = Pε , D51 εr = N D45 ur = Pε − D55 εr − D15 εz D41 ur = N − D51 εr D44 ur = Pu − D14 εz − D54 εr D11 εz = N − D51 εr − D41 ur D22 χx = Mx
D55 , D51 D45 D41 D44 D11 , H11 , a, b, c D22 , H22 , a, b, c
a b
Free
εz = 0
Free Free Free
Maximum and minimum values are given in Table 3. Values of Pε considered are given in Table 2.
Table 2 Load cases—steps and magnitudes. Loads applied to riser Pressure loading (MPa) Internal = 0 External = 0 Internal = 10 External = 7.80 Internal = 20 External = 15.61 Internal = 30 External = 23.412
1 2 3 4
Pε (kN)
Max. axial tension (kN)
Min. axial tension (kN)
0
100
0
145.6
100
0
291.3
100
0
436.9
100
0
Table 3 Maximum values or semi-amplitudes of the applied loads. Analysis Loading Units Minimum Maximum Amplitude
1 Pε kN 0 873.9 873.9
2 Pε kN 0 873.9 873.9
3 Pε kN 0 873.9 873.9
4 Pu M N/m 0 10 10
5 N kN 0 100 100
6 Mx k Nm −16 16 32
Analysis 1 1000
Pε (kN)
Analysis 2
1000
800
Pε (kN)
Load case
600
400
800
200
600
0 0
0.0005
0.001
0.0015
0.002
U r (m)
400
Fig. 8. Result of analysis 2.
200 0 0
0.02
0.04
0.06
0.08
0.1
εr Fig. 7. Result of analysis 1.
parameter, applied at top reference point (TRP). The fifth column shows Eq. (3) specialized to the case once the boundary condition has been applied. The sixth column shows the unknown target parameters which are to be found through each analysis. Analyses 1–3, 5 and 6, study the cyclic loading response, while analysis 4 studies the monotonic loading response of the riser. The details of the maximum and minimum values of the loads in column 3 of Table 1 are reported in Table 3. In analysis 1, movement of the average radius of the riser is restricted to the radial direction by fixing the inner surface of the outer anti-wear layer, in the radial direction. This is because the inner radius of the outer anti-wear layer coincides with the average radius of the riser excluding the carcass layer. Meanwhile fixing the
outer anti-wear layer allows the two adjacent helical armor layers to move freely, so that all layers can be compressed radially. Fig. 7 shows the applied Pε against the computed radial strain. In the second analysis the load is the same as in the first, except for the boundary condition at the TRP, which is not free. Fig. 8 shows the applied Pε against the radial displacement. Analysis 3 is the same as the first two, except for the boundary condition at the TRP, where displacement is restrained in the axial direction. Fig. 9 shows the applied Pε against the computed radial displacement. In the fourth analysis the riser is subject to a pure monotonic increase in Pu linearly with time. The TRP is free to move in all directions. Fig. 10 shows the applied Pu against the computed radial displacement. The fifth analysis corresponds to a cyclic axial tensile load simulation preceded by an initial ‘pressure load step’. The analysis starts with an initial ‘pressure load step’ which accounts for the internal fluid pressure and the external hydrostatic pressure. Four different input pressures are analyzed and referred to as Cases 1, 2, 3, and 4, respectively. The values of the applied internal and external pressures are given in Table 2. The pressure step is then
A. Bahtui et al. / Engineering Structures 32 (2010) 2287–2299
Analysis 3 1000
Pε ( k N )
800
600
400
200
0 0
0.0002
0.0004
0.0006
0.0008
Ur (m) Fig. 9. Result of analysis 3.
Analysis 4 12
Pu (MN/m)
9
6
2295
ternal and external pressure loading, the behavior becomes more stabilized, and a linear curve can be fitted to represent the overall behavior. Fig. 8 shows the radial displacement which resulted from analysis 2. The cyclic response of radial displacement due to cyclic internal and external loading is unlikely to be due to frictional dissipation between adjacent layers, but is potentially due to the gaps closure and opening at the initiation of each cycle. Therefore a linear curve fitting assumption for this figure makes sense. Radial displacements of analysis 3 are shown in Fig. 9. When radial displacements in this figure are compared with that of analysis 2, it can be seen that the latter are much lower and could be ignored, but are instead curve fitted linearly and considered in the constitutive model. Fig. 10 shows the radial displacement of analysis 4. This figure well illustrates the characteristics of the variable Pu . As mentioned in Section 2, an increase in this variable makes an inward or outward overall radial movement of all layers which neither creates any compression and thus frictional dissipation between layers, nor opens and closes the gaps between layers. Therefore the response is clearly linear and a linear curve fitting is well suited to this curve. All results of Fig. 11, except for the case Pε = 0, illustrate a hysteresis behavior of the riser subjected to cyclic axial tension (analysis 5), the response tending towards a stabilized cyclic response increasing gradually with the number of cycles. The first cycle represents the installation phase, whereas the stabilized results are more suitable for use in the calibration procedure, which will be discussed in the next section. This figure also shows that the higher the pressure, the greater the energy dissipated during the hysteresis loop.
6. Calibration of the constitutive model
3
0 0
0.003
0.006
0.009
0.012
Ur (m) Fig. 10. Result of analysis 4.
followed by a ‘tensile step’, in which a 100 kN axial tension is cyclically applied to the free end of the pipe. In each loading or unloading phase of the cycles the load is applied linearly with time. In Fig. 11(a)–(d) the applied axial force is plotted against the computed axial strain for the four different load cases defined in Table 2. Notice that the initial pressure load step results in initial axial strain, radial strain and radial displacement offsets, which are not shown in the figures. As expected, results are similar to those in Fig. 2. The sixth case analysis corresponds to a cyclic bending moment simulation preceded by an initial ‘pressure load step’. Analysis 6 was presented comprehensively by the authors in the previous work [2] and hence is not presented here. 5.3. Discussion Fig. 7 shows the radial strain for analysis 1. The initial jump corresponds to the first time step where all gaps between layers close. These gaps develop between the adjacent layers of the riser during mesh due to the limited number of elements in the circumferential direction of each layer. As the simulation goes through cyclic in-
Several FE simulations have been used to calibrate the constitutive model. Each result from the numerical simulations was curve-fitted either using a linear or bilinear curve, based on the stabilized cycle. The linear curve yields a constant slope (analysis 1–4). The bilinear curve can give three specific parameters for the constitutive model (analysis 5 and 6): the initial no-slip slope, a slip-initiation point, and a full-slip slope. For analyses 1–4 a linear response can be identified. These slopes directly provides the flexural stiffnesses D55 , D15 = D51 and D45 = D54 by using appropriate equations from the fifth column of Table 1. For all cases of Fig. 11, except those with zero or very small applied pressure, an initial approximately linear response can be identified. Therefore, an initial ‘no-slip slope’ can be determined and is found to be the same for each of the analyzed cases. This validates the hypothesis that an initial linear part of the force–displacement diagram, in which frictional slip between layers is absent or negligible, does exist. This slope directly provides the flexural stiffness D11 . After this initial linear response, a nonlinear part of the curve follows, which relatively rapidly tends to join a final full-slip slope. Upon repeated loading and unloading cycles, this final slope changes and stabilizes on a straight line, which is that used in the identification procedure. This final full-slip slope is not equal to the hardening parameter H11 , but it is related to it and therefore its determination allows the estimation of H11 . The intersection point between the initial slope and the final slope provides a ‘slip-initiation’ point, which is also a point on the boundary of the no-slip domain, i.e. a point at the zero level set of the slip-onset function in the N–Pε plane.
2296
A. Bahtui et al. / Engineering Structures 32 (2010) 2287–2299
a
b
d
c
Fig. 11. Axial force–axial strain curves for the four different load cases of analysis 6 and reported in Table 2. Cases are 1(a), 2(b), 3(c) and 4(d).
x 10 4 8
in many cases, it is important to underline that the method can be further refined, as discussed later. The slip-initiation points identified have been fitted with the quadratic expression of Eq. (5) and the resulting surface is shown in Fig. 12. Table 4 reports the parameters of the model determined with the proposed procedure.
N - βz
7 6 5 4 3
7. Verification of the model
2 1 0 5000 0 Mx - βx
-5000 0
1
2
4 3 Pε - βε
5
6 x 10 5
Fig. 12. Curve-fitting of the slip-initiation points.
6.1. Remarks In this proposed model, the nonlinear transitional regime between the initial no-slip response and the final full-slip response is replaced with a bilinear response. During the transition phase, some of the layers in contact (particularly the armor tendons) slip with respect to the others, while the others do not. Although the results of the proposed simplified model can be accurate enough
To evaluate the proposed method, the same tensile cyclic analyses performed using the detailed finite-element model in ABAQUS have been reproduced using the proposed constitutive model adopting the parameters of Table 4 and using the solution scheme described in Section 2. The results of this verification for four cases corresponding to four values of the applied pressure load Pε are reported in Fig. 13. For each case, a graph is reported in which the load–displacement curve obtained using the proposed constitutive model is plotted with a solid line and is compared with the dashed-line curve obtained using the detailed finite-element model. In Case 1 (Pε = 0), no friction occurs in both cases, which results in no hysteresis, i.e. no energy dissipation. The two curves overlap and cannot easily be distinguished. The difference between the stabilized response and that in the first one or two cycles increases with increasing values of Pε , and could be an issue of concern only in the analysis of the installation phase, while in the other cases it can be considered not important for the accuracy of the analysis.
A. Bahtui et al. / Engineering Structures 32 (2010) 2287–2299
2297
Table 4 Identified parameters. Stiffness D11 (N) D41 (N m−1 ) D51 (N) D22 (N m2 ) D44 (N m−2 ) D54 (N m−1 ) D55 (N)
2.5 · 108 1.28 · 108 −5.88 · 106 6.08 · 105 4.38 · 109 −1.36 · 108 1.52 · 107
Slip-onset function a b (N−1 ) c (N−1 m−2 )
0 1.2 · 10−4 2.3 · 10−2
Hardening H11 (N) H22 (N m2 )
7.3 · 107 2.5 · 105
b
a
d
c
Fig. 13. Comparison of the results obtained from the constitutive model (solid line) and from the detailed FE model (dotted line).
For verification related to bending-moment analysis, refer to [2]. 7.1. Discussion of the constitutive model The comparison between the curves provided by the constitutive model and the stabilized cycles obtained using the FE numerical simulations indicate a reasonable agreement in terms of the overall response. The model precisely captures the minimum and maximum magnitudes of all cases of different input pressures. The initial stiffnesses of all cases are well coincident with the stabilized initial slopes of the FE numerical simulations. The hardening stiffnesses of all cases also perfectly capture all the final full-slip slopes of the FE simulations.
It is worth noting that for all cases the same stiffness and hardening parameters have been used, that these parameters have been determined by interpolating the FE results with a bilinear curve only assuming a monotonic case, and that very good agreement is reached for cases with significantly different values of the pressure and for the entire stabilized cycles. This demonstrates that the assumptions on which the constitutive model is based are well justified and that the proposed approach could represent an excellent method of analysis of long flexible risers. The constitutive model is sufficiently robust to capture many practical loadings and has the capability of being easily improved, enhanced and extended. One possible development would be to introduce nonlinear kinematic hardening, which would then be able to match FE simulation or experimental results
2298
A. Bahtui et al. / Engineering Structures 32 (2010) 2287–2299
dσ
dε
= ε=Bˆ NL qˆ kn+1
dσ
dε
n+ 1
2 δθ111 − 4b2 τ12 Γ111 −4bc τ1 τ2 Γ111 ∆222 1 = −4bc τ1 τ3 Γ111 ∆333 δ δθ141 − 4b2 τ 2 Γ111 Γ141 1 δθ151 − 4b2 τ12 Γ111 Γ151
−4bc τ1 τ2 Γ111 ∆222 δ Σ2 − 4c 2 τ22 ∆2222 −4c 2 τ2 τ3 ∆222 ∆333 −4bc τ1 τ2 Γ141 ∆222 −4bc τ1 τ2 Γ151 ∆222
−4bc τ1 τ3 Γ111 ∆333 −4c 2 τ2 τ3 ∆222 ∆333 δ Σ3 − 4c 2 τ32 ∆2333 −4bc τ1 τ3 ∆141 ∆333 −4bc τ1 τ3 ∆151 ∆333
δθ411 − 2bτ1 Γ111 (2bτ1 Γ141 − Γ541 ) −2c τ2 ∆222 (2bτ1 Γ141 − Γ541 ) −2c τ3 ∆333 (2bτ1 Γ141 − Γ541 ) δθ441 − 2bτ1 Γ141 (2bτ1 Γ141 − Γ541 ) δθ451 − 2bτ1 Γ151 (2bτ1 Γ141 − Γ541 )
δθ511 − 2bτ1 Γ111 (2bτ1 Γ151 − Γ551 ) −2c τ2 ∆222 (2bτ1 Γ151 − Γ551 ) −2c τ3 ∆333 (2bτ1 Γ151 − Γ551 ) δθ541 − 2bτ1 Γ141 (2bτ1 Γ151 − Γ551 ) δθ551 − 2bτ1 Γ151 (2bτ1 Γ151 − Γ551 )
Box III.
the behavior properly, this paper presents a more accurate constitutive model for these structures. The constitutive model characterizes the hysteresis behavior of various layers of risers when subjected to different modes of loadings. The constitutive model was incorporated into a generalized finite element formulation which has been developed based on the Euler–Bernoulli beam theory together with two additional degrees of freedom. From the results it is noted that using the constitutive model together with the developed generalized finite element formulation, it is possible to study the nonlinear behavior of long flexible risers with significantly low computational costs. In future work, the constitutive law will be improved by including nonlinear hardening as well as more refined hardening models capable of simulating the stabilization of the response during cyclic loading. Further improvements can be made by accounting for the presence of the shear force (switching to a Timoshenko beam model) and by incorporating the influence of variable axial force and torque on the nonlinear response. Clearly, more sophisticated models will also imply more parameters to be identified and therefore call on the development of more refined identification procedures. Fig. 14. Cantilever un-bonded flexible riser subjected to an end bending moment.
Acknowledgements
much more accurately in terms of the gradual change of stiffness between the no-slip and full-slip phases.
The work has been funded by the UK Engineering and Physical Science research Council (EPSRC) under grant number EP/C510232/1. The authors would like to gratefully acknowledge their support. They would also like to thank Lloyds Register EMEA in London for the precious and fruitful discussions.
8. Large-scale analysis: one mechanical example The situation where a cantilever flexible riser with length 500 m, which is subjected to a 1.5 kN m end bending moment, is studied. Fig. 14 shows the result for this flexible riser for three different conditions: a linear solution with constant flexural stiffness EI = 608 108 Pa, a nonlinear solution with constant EI = 608 108 Pa and a nonlinear solution with the inclusion of the constitutive model with the initial value EI = 608 108 Pa (the first EI value which is taken from the constitutive model’s output). From Fig. 14 it is noted that the presented constitutive model is capable of capturing the energy dissipation in the flexible riser structure due to frictional slip between layers and the hysteresis response. The detailed finite element analysis of a flexible riser using commercial FEA packages on a cluster of parallel processors takes a significant time to arrive at the solution. Using the developed constitutive model and a generalized finite element formulation, as presented in Section 3, the computational time for studying the nonlinear behavior of several kilometers length of un-bonded flexible risers takes only a few seconds. This example reveals the capability and merit of the proposed constitutive model. 9. Conclusions In view of the fact that the equivalent modeling scheme (treating the riser as one homogeneous material layer) for the simulation of the various layers of an un-bonded flexible riser fails to model
Appendix According to the equation given in Box III,
δ = 4b2 τ12 Π11 + 4c 2 τ22 Ω22 + 4c 2 τ32 Ω33 − 2bτ1 Π15 2bλ Dij + Hij Dji + Hji Πij = Dij + Hij − , 1 + 2bλ (Dii + Hii ) 2c λ Dij + Hij Dji + Hji Ωij = Dij + Hij − 1 + 2c λ (Dii + Hii ) 2bλDkj Djk + Hjk , Γijk = Dij − 1 + 2bλ (Dkk + Hkk ) 2c λDkj Djk + Hjk ∆ijk = Dij − 1 + 2c λ (Dkk + Hkk ) 2bλDkj Dik θijk = Dij − , 1 + 2bλ (Dkk + Hkk ) 2c λD2ii Σi = Dii − . 1 + 2c λ (Dii + Dii ) References [1] API. American Petroleum Institute. Recommended Practice for Flexible Pipe. 2002. [2] Alfano G, Bahtui A, Bahai H. Numerical derivation of constitutive models for un-bonded flexible risers. Int J Mech Sci 2009;51(4):295–304.
A. Bahtui et al. / Engineering Structures 32 (2010) 2287–2299 [3] Bahtui A, Bahai H, Alfano G. Numerical and analytical simulation of un-bonded flexible risers subjected to combined modes of loading. In: 9th US national congress on computational mechanics. 2007. [4] Bahtui A, Bahai H, Alfano G. A finite element analysis for un-bonded flexible risers under axial tension. In: Proceedings of OMAE2008 27th international conference on offshore mechanics and arctic engineering. 2008. [5] Bahtui A, Bahai H, Alfano G. A finite element analysis for un-bonded flexible risers under torsion. J Offshore Mech Arctic Eng, OMAE 2008;130(4):041301. 4 pages. [6] Bahtui A, Bahai H, Alfano G. Constitutive modeling of un-bonded flexible risers under tension. In: 8th world congress on computational mechanics. 2008. [7] Bahtui A, Bahai H, Alfano G. Numerical and analytical modeling of un-bonded flexible risers. J Offshore Mech Arctic Eng, OMAE 2009;131(2):1010. [8] Kagoura T, Ishii K, Abe S, Inoue T, Hayashi T, Sakamoto T, et al. Development of a flexible pipe for pipe-in-pipe technology. Furukawa Rev 2003;24:69–75. [9] Felippa CA, Chung JS. Nonlinear static analysis of deep-ocean mining pipe EM DASH 1. Modeling and formulation. Amer Soc Mech Eng 1980;80(48): 6 pages. [10] Felippa CA, Chung JS. Nonlinear static analysis of deep-ocean mining pipe EM DASH 1. Modeling and formulation. Trans ASME, J Energy Resour Technol 1981;103(1):11–5. [11] Felippa CA, Chung JS. Nonlinear static analysis of deep-ocean mining pipe EM DASH 2. Numerical studies. Trans ASME, J Energy Resour Technol 1981;103(1): 16–25. [12] McNamara JF, Lane M. Practical modelling for articulated risers and loading columns. Trans ASME, J Energy Resour Technol 1984;106(4):444–50. [13] McNamara JF, Lane M. 3D frequency domain analysis of offshore structures. In: Proceedings of the 9th conference on engineering mechanics. 1992. p. 192–5. [14] O’Brien PJ, McNamara JF, Dunne FPE. Three-dimensional nonlinear motions of risers and offshore loading towers. J Offshore Mech Arctic Eng 1988;110(3): 232–7. [15] O’Brien PJ, McNamara JF. Analysis of flexible riser systems subject to threedimensional seastate loading. In: Proceedings of the international conference on behaviour of offshore structures, 6. 1988. p. 1373–88.
2299
[16] Hoffman D, Ismail N, Nielsen R, Chandwani R. Design of flexible marine riser in deep and shallow water. In: Proceedings 23rd annual offshore technology conference. 1991. p. 253–65. [17] Atadan AS, Calisal SM, Modi VJ, Guo Y. Analytical and numerical analysis of the dynamics of a marine riser connected to a floating platform. Ocean Eng 1997; 24(2):111–31. [18] Chai YT, Varyani KS, Barlt NDP. Three-dimensional Lump–Mass formulation of a catenary riser with bending, torsion and irregular seabed interaction effect. Ocean Eng 2002;29:1503–25. [19] Ong PPA, Pellegrino S. Modelling of seabed interaction in frequency domain analysis of mooring cables. In: Proceedings of 22nd international conference on offshore mechanics and arctic engineering ASME. 2003. p. 1–9. [20] Zhang Y, Chen B, Qiu L, Hill T, Case M. State of the art analytical tools improve optimization of un-bonded flexible pipes for deepwater environments. In: The 2003 offshore technology conference. 2003. [21] Yazdchi M, Crisfield MA. Buoyancy forces and the 2D finite element analysis of flexible offshore pipes and risers. Internat J Numer Methods Engrg 2002;54: 61–88. [22] Yazdchi M, Crisfield MA. Nonlinear dynamic behavior of flexible marine pipes and risers. Internat J Numer Methods Engrg 2002;54:1265–308. [23] Lacarbonara W, Pacitti A. Nonlinear modeling of cables with flexural stiffness. Math Probl Eng 2008;(370767):1–21. [24] Bech A, Skallerud B, Sodahl N. Structural damping in design analysis of flexible risers. In: Proceedings of the first European conference on flexible pipes, umbilicals and marine cables. 1992. 14 pages. [25] Smith R, O’Brien P, O’Sullivan T, Weibe C. Fatigue analysis of un-bonded flexible risers with irregular seas and hysteresis. In: Offshore technology conference. 2007. [26] Alfano G, Rosati L, Valoroso N. A tangent–secant approach to rate-independent elastoplasticity: formulations and computational issues. Comput Methods Appl Mech Engrg 1999;179:379–405. [27] Han W, Reddy BD. Plasticity—mathematical theory and numerical analysis. New York: Springer-Verlag; 1999 [Chapter 3].
Contents lists available at ScienceDirect
Engineering Structures journal homepage: www.elsevier.com/locate/engstruct
On the multi-scale computation of un-bonded flexible risers A. Bahtui a,∗ , G. Alfano a , H. Bahai a , S.A. Hosseini-Kordkheili b a
School of Engineering and Design, Brunel University, London, United Kingdom
b
Sharif University of Technology, Aerospace Engineering Department, Tehran, Iran
article
info
Article history: Received 5 August 2009 Received in revised form 10 February 2010 Accepted 1 April 2010 Available online 11 May 2010 Keywords: Non-associative plasticity Frictional slipping Parameter identification Generalized finite-element formulation Large-displacements
abstract The purpose of this paper is to model the detailed effects of interactions that take place between components of un-bonded flexible risers, and to study the three-dimensional motion responses of risers when subjected to axial loads, bending moments, and internal and external pressures. A constitutive law for un-bonded flexible risers is proposed and a procedure for the identification of the related input parameters is developed using a multi-scale approach. A generalized finite element structural model based on the Euler–Bernoulli beam theory is developed in which the constitutive law is embedded. The beam theory is enhanced by the addition of suitable pressure terms to the generalized stresses to account for the internal and external pressures, and it can therefore be efficiently used for large-scale analyses. The developed nonlinear relationship between generalized stresses and strains in the beam is based 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. The merit of the present constitutive law lies in its ability to capture many important aspects of the structural response of risers, including the energy dissipation due to frictional slip between layers and the hysteresis response. Crown Copyright © 2010 Published by Elsevier Ltd. All rights reserved.
1. Introduction Multilayered un-bonded flexible risers are slender marine structures that are widely used in offshore production to convey oil between the well-head and the offshore rig. The main advantage of flexible risers when compared to rigid steel risers is their much lower bending stiffness, leading to smaller radii of curvature with the same pressure capacity, due to the complex makeup of flexible risers [1]. This, in turn results in an increased ability to undergo large deformations under loads induced by the ocean current and wave, vortex-induced vibrations, the motion of the floating vessel and the installation procedure [1]. The detailed layout and description of layers and the behavior of un-bonded flexible risers are addressed in the literature [2–8]. A typical assumption made in the large-scale analyses of flexible risers is to replace all layers of a riser with an equivalent linearly elastic beam having a constant tangent stiffness matrix. The major advantage given by this assumption is that it makes the computational solution time significantly shorter, and renders it suitable for day-to-day industrial projects. A static analysis procedure for the numerical determination of nonlinear static equilibrium configurations of deep-ocean risers
∗
Corresponding author. Tel.: +44 7828487635. E-mail addresses: [email protected] (A. Bahtui), [email protected] (G. Alfano), [email protected] (H. Bahai), [email protected] (S.A. Hosseini-Kordkheili).
was developed by Felippa and Chung [9]. The riser was modeled by three-dimensional beam finite elements which include axial, bending, and torsional deformations. They extended their model by taking the deformations coupled through geometrically nonlinear effects [10]. The resulting tangent stiffness matrix included three contributions identified as linear, geometric (initial stress) and initial displacement stiffness matrices. For the solution, a combination of load–parameter incrementation, state updating of fluid properties and corrective Newton–Raphson iteration was used. The detailed implementation of their procedure is presented in Felippa and Chung [11] in which they discussed the importance of drag force along the pipe, varying significantly along the depth. The authors are of the opinion that the sensitivity of the results to environment characterization necessitates a review of many modeling and analysis practices. McNamara and Lane [12] studied the two-dimensional response of the linear and nonlinear static and dynamic motions of offshore systems such as risers and single-leg mooring towers. Their proposed method was based on the finite element approach using convected coordinates for arbitrary large rotations and includes terms due to loads such as buoyancy, gravity, random waves, currents, ship motions and Morison’s equation. The same authors extended their work by developing a three-dimensional computational dynamic analysis of deep water multi-line flexible risers in the frequency domain [13]. O’Brien et al. [14] presented the three-dimensional finite element modeling of marine risers, pipelines and offshore loading towers based on the separation of the rigid body motions and deformations of elements under
0141-0296/$ – see front matter Crown Copyright © 2010 Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.engstruct.2010.04.003
2288
A. Bahtui et al. / Engineering Structures 32 (2010) 2287–2299
conditions of finite rotations. This paper treats risers as a homogeneous material and includes all nonlinear effects including geometry changes, bending-axial and bending-torsional coupling and follower forces and pressures. The authors preceded their work by dividing the riser into a series of catenaries where the buoyancy modules were described as inverted catenaries [15]. Hoffman et al. [16] reviewed the design technique of deep and shallow water marine riser systems as well as their dynamic analytical and numerical analysis and the nonlinearities arising from hydrodynamic loading and dynamic boundary conditions. Their paper contains design methodology criteria, parameters and procedures of flexible riser systems while treating the riser as one homogeneous material layer. Atadan et al. studied the dynamic three-dimensional response of risers in the presence of ocean waves and ocean currents undergoing large deflections and rotations. They included shear effects based on the nonlinear elastic theory in their formulation. It was concluded that the length of the riser is the most important parameter which affects the deflections of the marine-riser [17]. Chai et al. proposed a three-dimensional lumped-mass formulation for the riser structure analysis which is capable of handling irregular seabed interaction. They adopted a simplified model by replacing the seabed surface with an elastic foundation with independent elastic springs having an arbitrary thickness [18]. Ong and Pellegrino [19] studied the nonlinear dynamic behavior of mooring cables in the frequency domain. They ignored the effects of friction and impact between the cable and the seabed. Their proposed method modeled the time-varying boundary condition at the touchdown by replacing the section of cable interacting with the seabed with a system of coupled linear springs. This work again ignored the hysteresis behavior of cables by assuming an equivalent linearly elastic pipe representing all layers. Zhang et al. [20] discussed analytical tools for improving the performance of flexible pipes. Their work used an equivalent linear bending stiffness which was derived from experimental data to calculate the maximum bending angle range. It contained reports on irregular wave fatigue analysis, collapse, axial compression and bird-caging for riser systems. They observed that the combined bending, axial compression and torsion could lead to the tendon being separated from the cylinder in a helix layer, and may lead to out of plane buckling. Yazdchi and Crisfield used a simple two-dimensional lowerorder beam element formulation for static nonlinear analysis of riser structures including the effects of buoyancy, steady-state current loading and riser top-tension. They assumed a linearly elastic material property for the riser by employing a constant modulus of elasticity, and studied the static behavior of flexible pipelines and marine risers, using finite elements developed for conventional nonlinear analysis [21]. The same authors continued their research by using a beam finite element formulation, based on the Reissner–Simo beam theory, for the static and dynamic nonlinear analysis of three-dimensional flexible pipes and riser systems in the presence of hydrostatic and hydrodynamic forces. Employing linearly elastic material properties, their work concentrates on the nonlinearities due to the fluid loading and the associated large deformations and considers hydrodynamic forces due to effects such as wave, drag and current action [22]. In all the above models, the assumption of perfect bonding between layers fails to take into consideration energy dissipation due to frictional slip between layers and the hysteresis response of the risers, generating inaccuracy in the results and making the analysis unsuitable for accounting for the hysteresis damping in vibration or fatigue. Some nonlinearities have been considered by Lacarbonara and Pacitti, who presented a geometrically exact formulation for the
dynamic analysis of cables with linearly elastic material behaviors when undergoing axis stretching and flexural curvature. Their formulation was based on nonlinearly viscoelastic constitutive laws for the tension and bending moment with an additional constitutive nonlinearity accounting for the no-compression condition [23]. Some works in the literature consider a bi-linear bending stiffness, taking into account the energy dissipation originating from frictional contact between layers in the bending moment analysis [24,25]. However the bi-linear bending stiffness proposed in these papers does not depend on internal and external pressures. This generates inaccuracies in the results when studying the response of a long riser floating in ultra-deep water. Moreover the study has been performed only on the hysteresis effect in bending analysis whereby hysteresis effects for axial strain or any coupling with radial displacement is ignored. On the other hand, the bi-linear model in these works is based on two different stiffnesses, one when perfect bonding is assumed between layers and one when slippage occurs between layers. The limitations outlined in the above contributions motivated the authors to study the effect of the interaction between layers of un-bonded flexible risers on their hysteresis behavior in more detail. The present work addresses the above issues by combining the accuracy of a detailed finite element model of a short part of riser, with the low computational cost of a three-dimensional Euler–Bernoulli beam element, in the framework of a novel, multiscale approach to conduct nonlinear three-dimensional largescale structural analyses of an entire un-bonded flexible riser. The ‘equivalent elasto-plastic’ constitutive law for the beam cross section proposed by the authors [2], which relates axial strain, torsion and bending curvatures to the conjugate stress resultants is enhanced by including additional degrees of freedom for the beam model which account for axial strain and radial displacement. This law models the small-scale frictional slip occurring between the layers, when a combination of values of the stress resultants and of the internal and external pressures overcomes a given threshold, at the macroscopic level of the cross section. The ‘equivalent yield function’ proposed by the authors [2], which depends on the current values of the stress resultants, on a set of internal history variables and on the values of the external and internal pressures acting on or within the pipe, is enhanced by taking into account the axial stress effects, making the equivalent yield function threedimensional. This last functional relationship simulates the effect that the pressure values have on the inter-layer normal stress, which in turn affects the inter-layer frictional sliding and the overall dissipation. The zero level set of the yield function represents the boundary of a ‘no-slip domain’, that is, the onset of slipping. In formulating this model a linear kinematic hardening law has been employed. The assumptions of linear hardening and of linear elastic behavior within the ‘no-slip domain’ result in a bilinear response, which can provide sufficient approximation of the structural response in many cases of interest. This simplified model entails that the single cross section of the riser can be either in a state in which ‘no-slip’ is found, or in a state in which ‘full-slip’ occurs. To identify the parameters of the beam model to be used in large-scale analyses of the entire riser, ideally one would need a great number of results of experimental tests conducted on a reasonably large model of a flexible pipe using an experimental rig which reproduces the real conditions as closely as possible. This is very difficult and extremely expensive, which is the main motivation for the development of the multi-scale approach proposed in this paper, in which the small-scale detailed FE model is employed as a virtual experimental rig, and the parameters of the large-scale model are then determined by combining engineering assumptions and careful curve fitting of the results of the smallscale and large-scale models.
A. Bahtui et al. / Engineering Structures 32 (2010) 2287–2299
2289
In addition, the internal and external pressures need to be considered in the model. If εr and ur are the average radial strain and the radial displacement of the average radius, respectively, Pε and Pu are the conjugate actions working for εr and ur . For a detailed derivation of these variables please refer to [2]. Hence, in the proposed model N , Mt , Mx , My , Pu and Pε are the generalized stresses. The corresponding work-conjugated generalized strains are εz , φ, χx , χy , ur , εr , where εz , χx , χy and φ are the axial strain, the curvatures about the x and y directions and the torsional curvature. The detailed definition of the foregoing is in Section 3. A compact notation is conveniently introduced denoting by {σ} and {ε} the generalized stresses and strains:
N
Fig. 1. Stress resultants.
A generalized finite element beam model supplemented by the constitutive law is used for large-scale analysis, unlike the detailed finite element model used to formulate it and to determine its parameters. This model is able to capture all the important aspects of the structural response of the riser and is suitable for the study of many practical dynamic problems, such as hysteresis damping in vibration, fatigue or service life analysis. Forthcoming papers will present results related to these situations. This paper is divided into nine sections. In Section 2 the constitutive law for the three-dimensional beam model is formulated. In Section 3 the comprehensive global analysis of the riser is formulated. In Section 4 the tangent operator matrix is derived. In Section 5 the detailed finite element model, mainly based on earlier work by the authors [4,7], is described. In Section 6 the key issue of the identification of the parameters of the beam constitutive model is addressed. In Section 7 one validation case is presented and discussed in order to determine the extent to which the results of the various detailed finite-element simulations for a short riser, subject to internal and external uniform pressure and uniform cyclic loading, agreed with those given by the proposed constitutive law. A practical case of a riser is then presented in Section 8. Results from a horizontal un-bonded flexible riser with a length equal to 500 m subject to top end bending moment due to vessel movements, are compared to a steel rigid riser with the same length and conditions. Finally, in Section 9, conclusions are drawn based on the results presented. 2. Formulation of the constitutive model The constitutive law of an un-bonded flexible riser for an Euler–Bernoulli three-dimensional beam model, proposed in [2], is significantly enhanced in this section to make it more suitable for practical usage. In particular, axial strain and radial displacement are taken into account and the ‘slip-onset’ function is made dependent on both the axial force and the bending moment resultants. It is assumed that locally the pipe develops along a straight line, so that an infinitesimal element of pipe can be represented as in Fig. 1. A local right-handed Cartesian system is introduced with the origin and the x and y axes located at the cross section and the z axis coincident with the centroid axis. 2.1. Generalized stresses and strains The stress resultants, shown in Fig. 1, are the axial force N, the torque Mt , and the bending moments around axes x and y, denoted by Mx and My , respectively.
Mt Mx {σ} = My Pu Pε
ε z
φ χx {ε} = χy .
(1)
ur
εr
2.2. Constitutive law As discussed in Section 1, the constitutive law for the beam model is based on the observation that the local (small-scale) frictional slip between the different layers results in a macroscopic (large-scale) relationship between generalized stresses and strains which has many analogies with the laws of elasto-plasticity. The formulation of a constitutive law for the bending stiffness was presented by the same authors [2] and hence will not be presented here. The derivation of a constitutive law for the axial stiffness behavior is based on a similar procedure, leading to a different ‘slip-onset’ function. Fig. 2 shows the response of the flexible riser under cyclic loading, evaluated using a detailed threedimensional finite element analysis which accurately models the interaction between all layers. The case study of axial compression loading has not been considered in this work. As a consequence, Fig. 2 or all other axial force figures have a minimum value of zero. Besides, the no-slip portion at unloading has approximately same length as the no-slip portion in the initial cycle. For the qualitative analysis in this section, presented in order to justify the assumptions of the proposed constitutive model, it is sufficient to state that the values of the applied loading (internal and external pressures included) are, within a meaningful range, typical of practical cases. For the first monotonic increase of the axial strain, the axial force increases with an almost linear stiffness between points (a) and (b). Then the stiffness gradually decreases between points (b) and (c) to a value which remains approximately constant up to point (d). Upon unloading, the axial stiffness is initially very close to the initial stiffness, between points (d) and (e), and then it rapidly decreases between points (e) and (f) to a value which is constant between points (f) and (a) and approximately equal to that between points (c) and (d). Replacing the rounded parts between points (b) and (c) and between points (e) and (f) with two sharp elbows, a bilinear curve response is obtained with good approximation. Such a bilinear curve represents the same type of response as that obtained using an elasto-plastic model with linear kinematic hardening, which is the model proposed here. In analogy with the elasto-plastic case, the generalized strains are additively decomposed into an elastic-like, i.e. ‘no-slip’, part {εe }, and a plastic-like, i.e. ‘full-slip’, part {εs }.
{ε} = {εe } + {εs } .
(2)
For full details of the formulation of the constitutive law, particularly the ‘slip-onset’ function f and the real-valued ‘slip potential’ g, please refer to [2].
2290
A. Bahtui et al. / Engineering Structures 32 (2010) 2287–2299
denoted by βz , βx , βy , βr and βε , so that the argument of the onsetslip function is given by: N − βz Mx − βx {σ} − {β} = My − βy . P − β u r Pε − βε
(4)
To determine functions f and g and the coefficients of the kinematic hardening matrix, many results from several numerical simulations made with the detailed FE model have been reviewed and the following expression has been found to be sufficiently valid: Fig. 2. Plot of the axial force against the corresponding strain for an un-bonded flexible riser during cyclic loading; the response is evaluated using a detailed finiteelement simulation.
2.2.1. Specialization to the case of torsion-less loading In this section, the constitutive law proposed previously by the same authors [2] is significantly enhanced by considering the very important issue of cyclic axial loading and its coupling effect on the slip-onset function. The number of degrees of freedom of the constitutive model is increased to five by incorporating the axial strain, curvatures, radial displacement and the average radial strain εr at each node. The same reasonable assumption made by the authors [2] that the torque Mt is constant throughout the deformation process is retained here to better focus on the factors which have a major influence on the structural response. Hence, its influence on the determination of functions f and g can be incorporated within the coefficients and then ignored. Furthermore, a linear elastic relationship is assumed between Mt and the conjugated variable φ , so that the hypothesis made that Mt is constant allows us to completely ignore this term. A linear elastic relationship is also assumed between Pu and the conjugated variable ur . On the other hand, the main influence on the response of the flexible pipe at the onset of and during frictional slipping is given by the other term Pε . Considering its corresponding work-conjugated generalized strains εr as a degree of freedom of the constitutive model allows for better modeling of the coupling effects between Pε and N. This enhancement to the previous work of the authors [2] allows us to consider the coupling terms between the axial and radial strains. As for the coupling effects between the axial strain and the curvature, they are due to the changes in the radial contact stresses between layers caused by changes in the axial force due to Poisson effects, which in turn affects friction. Hence, these effects are embedded into the slip-onset function, and thus are not considered in the stiffness matrix. Instead, no elastic coupling exists between these variables and between radial displacement and the curvature, which results in zero values of the corresponding terms in the stiffness matrix. Finally, symmetry of the cross section results in the same flexural response about the x and y axes, with no elastic coupling between the two curvatures, or between the curvatures and the axial strain. Hence, the initial stiffness is represented by:
N D11 Mx My = P u Pε
0 D22 sym.
0 0 D22
D41 0 0 D44
D51 εz 0 χx 0 χy . D54 ur D55 εr
(3)
Note that coupling exists between axial generalized strain and radial strain and displacement terms. The components of the back-stress vector {β} which correspond to the generalized-stress components N, Mx , My , Pu and Pε are
h b (N − βz )2 + c (Mx − βx )2 2 i f ({σ} − {β}) = − (Pε − βε ) − a, + My − βy N − βz > 0 − (N − βz ) − a, N − βz ≤ 0.
(5)
where a, b and c are material parameters to be identified, and symmetry of the cross section has been exploited. The slip-onset surface in this case has two surfaces. This is due to the fact that the riser is not assumed to be capable of tolerating axial compression. This indicates that axial force is non-negative. The slip potential g has to account for the finding that frictional slipping is not accompanied by any significant relative opening between layers. At the small-scale level this means that no significant dilatancy is observed. A suitable expression for the potential g which meets this requirement is as follows:
h b (N − βz )2 + c (Mx − βx )2 2 i g ({σ} − {β}) = + My − βy , N − βz > 0 − (N − βz ) , N − βz ≤ 0
(6)
which results in the following slip rule:
∂g ∂N ∂g ∂M ε˙ sz x ∂g χ˙ sx ∂g {˙εs } = λ˙ = λ˙ {m} ⇒ χ˙ sy = λ˙ ∂ My u˙ ∂ {σ} sr ∂g ε˙ sr ∂ Pu ∂g ∂ Pε 2bN 2cMx λ˙ 2cMy , (N − βz ) > 0 0 0 = −1 0 ˙ 0 , (N − βz ) ≤ 0 λ 0
(7)
0
where ε˙ sz , χ˙ sx , χ˙ sy , u˙ sr and ε˙ sr denote the full-slip components of the axial tension, the bending curvatures about the x and y axes, the radial deformation and of the radial strain. Overlapping the generalized stresses and strains spaces, Fig. 3 describes two crosssectional views (a, b) as well as a three-dimensional view (c) of the above three-dimensional equation and shows that the difference
A. Bahtui et al. / Engineering Structures 32 (2010) 2287–2299
2291
N – βz
Mx – βx m
m n
n
Pε – βε
Pε – βε N
Mx
b
a x 10 4 2 m
N - βz
1.5 n
1 0.5 0 2 x 104
1
6 0 Mx - βx
-2
x 10 6
Pε - βε
0
c Fig. 3. Non-associative behavior of the proposed model.
between functions f and g results in the full-slip strain rate {˙εs } being not normal to the slip-onset surface (a, b), i.e. not parallel to the normal {n}. This is in analogy with the non-associative plasticity models used for many materials exhibiting internal friction for which the normality rule is not satisfied. The following proportional hardening law is used:
βz βx {β} = [H] {εs } ⇒ βy β r βε H11
=
0 H22
sym.
0 0 H22
0 0 0 0
0 εsz 0 χsx 0 χsy . 0 usr 0 εsr
(8)
It is worth observing that because of the slip law in Eq. (7), the radial full-slip strains usr and εsr result to be constantly zero, so that βr and βε also do so. A final remark is that, having the same response about the x and y axes because of symmetry, a simple plane model can be used to study the structural response and to identify the model parameters, as will be done in Sections 3 and 4. Nevertheless, it was preferred to derive all the above equations for the general case in which components of moment and curvature about both axes are non-zero to clarify that the model can be used in a general threedimensional problem. 3. Generalized finite element formulation A common approach to the numerical analysis of riser structures is to use an equivalent modeling scheme to simulate the ex-
istence of the various layers which make up the structure of the flexible riser. This approach fails to account for the interaction between layers, the hysteresis effects due to the occurrence of interlayer slip-stick and also the nonlinear behaviors due to the friction between the un-bonded layers. In this section a generalized nonlinear finite element formulation is developed in order to analyze the un-bonded flexible risers, embedding the constitutive model presented in Section 2. The generalized nonlinear finite element formulation is based on the Euler–Bernoulli beam theory with two extra degrees of freedom in the radial and thickness directions. Referring to Fig. 4, displacement components ui , vi and wi for node i are defined as ui = u0i (z ) − y
∂w0i (z ) ∂v0i (z ) +x ∂z ∂z
(9)
vi = v0i (z ) wi = w0i (z )
(10) (11)
where i = 1, 2 indicates the node numbers, u0i (z ), v0i (z ) and w0i (z ) are the displacement of the centre-line in x, y and z directions, respectively. Also ∂w0i (z )/∂ z and ∂ u0i (z )/∂ z are rotation component about x and y axis for node i. Using Eqs. (9)–(11) the nonlinear strain–displacement relations, based on the von Karman assumption, are
εzz =
∂ u0 1 + ∂z 2
εxx = 0 εxz = 0.
∂w0 ∂v0 + ∂z ∂z
2 −y
∂ 2 w0 ∂ 2 v0 +x 2 2 ∂z ∂z
(12) (13) (14)
To develop a nonlinear finite element model for this Euler– Bernoulli beam element we interpolate the displacement values
2292
A. Bahtui et al. / Engineering Structures 32 (2010) 2287–2299
k+1 [B]tNL
2 0 N l 1 k {q}n+1 = 0
B12 N003 /τnk+1
0
2 02 N5 θ1 n + 1k l2 N005 /τnk+1
0
0
2 0 N2 l 0
0 0 l 2
0
N1
B16
B17
N004 /τnk+1
N006 /τnk+1
0
0
0
l 2
0 0 0
N2
(20)
Box I.
+1 {q}kn+1 is given in Box I where [J] is the Jacobean matrix and [B]ktNL (Eq. (20)), where
B12 = B16 = B17 =
2N30 l2 2N40 l2 2N60 l2
k 0 k 0 k 0 k N30 w1n +1 + 2N4 w2n+1 + 2N5 θ1n+1 + 2N6 θ2n+1 k 0 k 0 k N40 w2n +1 + 2N5 θ1n+1 + 2N6 θ2n+1
(21)
2N5 θ1n+1 + N6 θ2n+1 0 k
0 k
Fig. 4. Euler–Bernoulli beam element with extra degree of freedom to model the radial displacement.
as well as generalized strain components ur and εr in the element field using their nodal values as follows
v = N3 v1 + N4 v2 + N5 θ1 + N6 θ2 w = N3 w1 + N4 w2 + N5 θ1 + N6 θ2 ur = N1 ur1 + N2 ur2 εr = N1 εr1 + N2 εr2
(e)
∗ σg εg · εg dz = {F}ext · q˜ za ∗ εg = [B]NL q˜ =
(16)
where Wi is virtual work, superscript (e) shows the elemental value, subscript g uses to indicate the globalcoordinate system pa ∗ rameter, {F}ext is the external load vector, q˜ is the virtual displacement vector, and [B]NL is the nonlinear strain–displacement matrix. Using Eqs. (16), the incremental finite element equation of equilibrium can be written as
{R} ({q}) = {F}int ({q}) − {F}int = 0
derivation of the consistent tangent operator
d{σ} d{ε}
{ε}=[B]NL {q}kn+1
is
1 −1 k k {q}kn+ +1 = {q}n+1 + [K]tNL {q}n+1
+1 {F}next − {F}int {q}kn+1
(18)
where k is the iteration number, and n stands for increment number, [K]tNL is the nonlinear tangent stiffness matrix, and
Z
[B]TNL V
d {ε}
d {σ} {ε}=[B]NL {q}kn+1
× [B]NL det [J] dV
Alfano et al. [26] presented a general approach to the evaluation of the tangent stiffness matrix for associated rate-independent elasto-plasticity materials. This section presents the derivation of this coefficient which follows a similar approach extending it to the non-associated case of interest here. From a computational point of view, the time domain needs to be subdivided into a number of finite steps. Adopting a fully implicit backward-Euler time-integration scheme, the equations to solve in each step are as follows:
{σ}n+1 = [D] {ε}n+1 − {εs }n+1 ∂ g {ε } {ε } − = 1 λ s n s n +1 ∂ {σ} n+1 {β}n+1 = [H] {εs }n+1 n +1 ≤ 0 1λ ≥ 0 f {σ}n+1 − {β} 1λf {σ}n+1 − {β}n+1 = 0.
(22)
Differentiating Eq. (221 ) with respect to {ε}n+1 , we get
∂ {σ}n+1 ∂ {εs }n+1 = [D] [I] − . ∂ {ε}n+1 ∂ {ε}n+1
(23)
(17)
where {q} is the unknown matrix in the form {q} = (u1 w1 θ1 ur1 u2 w2 θ2 ur2 )T . Implementing Newton–Raphson method to solve Eq. (17) leads to
[K]tNL {q}kn+1 =
where l is the element length, N 0 is first derivative of the shape functions, and N 00 is second derivative of the shape functions. The
4. Consistent tangent operator
zb
2 3/2
(15)
where u1 and u2 are the axial nodal displacements, w1 and w2 are transitional displacements, θ1 and θ2 are nodal rotational components. Also in Eq. (15), N1 and N2 are linear shape functions, while N3 –N6 are cubic shape functions. For this nonlinear problem, the governing equation in rectangular Cartesian coordinates is derived using the principle of minimum total potential energy, i.e. Wi
k 0 k 0 k 0 k τnk+1 = 1 + N30 w1n +1 + N4 w2n+1 + N5 θ1n+1 + N6 θ2n+1
presented in the next section.
u = N 1 u1 + N 2 u2
Z
(19)
From Eq. (222 ) we get the expression
∂ 2 g d {σ}n+1 ∂ {σ}2 n+1 ∂ 2g d {β}n+1 . + 1λ ∂ {σ} ∂ {β} n+1
d {εs }n+1 = d1λ {m} + 1λ
(24)
The assumed definition (Eq. (5)) implies that
∂ 2g ∂ 2g ∂ 2g = = . ∂ {σ} ∂ {β} ∂ {β} ∂ {σ} ∂ {β}2
(25)
Substituting Eq. (25) into Eq. (24), and then substituting the result into Eq. (23) and the derivative of (223 ) provides the following
A. Bahtui et al. / Engineering Structures 32 (2010) 2287–2299
d {ε}
d {σ}
2293
= ([D] − λ [D] [Z] [D])
n +1
−
(31)
[[D] − λ [D] [Z] ([D] + [H])] {m} ⊗ [[D] − λ [D] [Z] ([D] + [H])] {n} [([D] + [H]) − λ ([D] + [H]) [Z] ([D] + [H])] {m} · {n} Box II.
system of equations:
−1 ∂ 2 g ∂ 2 g [D] + 1λ − 1 λ 2 ∂ {σ} ∂ {σ}2 n+1 n +1 2 ∂ 2 g − 1λ ∂ g [H]−1 + 1λ 2 ∂ {σ} n+1 ∂ {σ}2 n+1 d {σ}n+1 d {ε}n+1 − d1λ {m} . × = d {β}n+1 d1λ {m}
(26)
From which we obtain d {σ}n+1 = [D] d {ε}n+1 − 1λ [D] [Q] [F]−1 [D] d {ε}n+1
+ 1λd1λ [D] [Q] [F]−1 ([D] + [H]) {m} − d1λ [D] {m} d {β}n+1 = 1λ [H] [Q] [F]−1 [D] d {ε}n+1
(27)
− 1λd1λ [H] [Q] [F]−1 ([D] + [H]) {m} + d1λ [H] {m} where Fig. 5. Detailed geometry of riser (cross-sectional view).
∂ 2 g [Q] = . ∂ {σ}2 n+1
(28)
The fulfillment of Prager’s consistency condition f implies [27]
= f˙ = 0
d {σ}n+1 − d {β}n+1 · {n} = 0.
(29)
Substituting Eq. (27) into Eq. (29) results in 1λ as below [D] − 1λ ([D] + [H]) [Q] [F]−1 [D] d {ε}n+1 · {n} 1λ = . ([D] + [H]) − 1λ ([D] + [H]) [Q] [F]−1 ([D] + [H]) {m} · {n}
(30)
By substituting Eq. (30) into (271 ) we obtain the equation in Box II, where [Z] = [Q] [F]−1 .
(32)
Substituting matrices and performing cross and dot products yields to the final equation which is presented in Appendix. Fig. 6. Finite-element mesh.
5. Detailed finite element model 5.1. Modeling Several FE simulations for a typical 1.7 m long un-bonded flexible riser have been conducted using the FE code ABAQUS to estimate the parameters of the constitutive law for the beam model. The FE model includes a complex makeup of seven-layers of internal and external plastic sheaths, helical armors, carcass and anti-wear layers and is described in the cross-sectional view of Fig. 5. It has been created by adding one additional layer to the model considered by the authors in previous work [4,5]. Fig. 6 shows the finite element mesh of the riser and indicates the number of nodes and elements used. The full details of the finite element model are available in [2–5] and are not presented here for the sake of brevity.
5.2. Analysis For the identification of the parameters, constructing the constitutive model, six sets of analyses are performed and are summarized in Table 1. The second column indicates whether the analysis is performed for several different values of input pressures Pε , referring to each one as an individual ‘load case’. In these analyses (5 and 6) the simulation starts with an initial ‘pressure load step’ which accounts for the average pressures Pε . The pressure step is then followed by the load(s) from column three. The different considered values of Pε for each axial load case are reported in Table 2. The third column, ‘Loadings’, shows the type of load which is applied at each analysis. The fourth column determines the additional boundary conditions required for the derivation of a target
2294
A. Bahtui et al. / Engineering Structures 32 (2010) 2287–2299
Table 1 Detailed information about individual analyses. No.
Dependency on Pε
Loadingsa
BC at the TRP
Equation to be solved
Target parameter(s)
1 2 3 4 5 6
– – – No Yesb Yesb
Only Pε Only Pε Only Pε Only Pu N Mx
εz = 0, ur = 0
D55 εr = Pε , D51 εr = N D45 ur = Pε − D55 εr − D15 εz D41 ur = N − D51 εr D44 ur = Pu − D14 εz − D54 εr D11 εz = N − D51 εr − D41 ur D22 χx = Mx
D55 , D51 D45 D41 D44 D11 , H11 , a, b, c D22 , H22 , a, b, c
a b
Free
εz = 0
Free Free Free
Maximum and minimum values are given in Table 3. Values of Pε considered are given in Table 2.
Table 2 Load cases—steps and magnitudes. Loads applied to riser Pressure loading (MPa) Internal = 0 External = 0 Internal = 10 External = 7.80 Internal = 20 External = 15.61 Internal = 30 External = 23.412
1 2 3 4
Pε (kN)
Max. axial tension (kN)
Min. axial tension (kN)
0
100
0
145.6
100
0
291.3
100
0
436.9
100
0
Table 3 Maximum values or semi-amplitudes of the applied loads. Analysis Loading Units Minimum Maximum Amplitude
1 Pε kN 0 873.9 873.9
2 Pε kN 0 873.9 873.9
3 Pε kN 0 873.9 873.9
4 Pu M N/m 0 10 10
5 N kN 0 100 100
6 Mx k Nm −16 16 32
Analysis 1 1000
Pε (kN)
Analysis 2
1000
800
Pε (kN)
Load case
600
400
800
200
600
0 0
0.0005
0.001
0.0015
0.002
U r (m)
400
Fig. 8. Result of analysis 2.
200 0 0
0.02
0.04
0.06
0.08
0.1
εr Fig. 7. Result of analysis 1.
parameter, applied at top reference point (TRP). The fifth column shows Eq. (3) specialized to the case once the boundary condition has been applied. The sixth column shows the unknown target parameters which are to be found through each analysis. Analyses 1–3, 5 and 6, study the cyclic loading response, while analysis 4 studies the monotonic loading response of the riser. The details of the maximum and minimum values of the loads in column 3 of Table 1 are reported in Table 3. In analysis 1, movement of the average radius of the riser is restricted to the radial direction by fixing the inner surface of the outer anti-wear layer, in the radial direction. This is because the inner radius of the outer anti-wear layer coincides with the average radius of the riser excluding the carcass layer. Meanwhile fixing the
outer anti-wear layer allows the two adjacent helical armor layers to move freely, so that all layers can be compressed radially. Fig. 7 shows the applied Pε against the computed radial strain. In the second analysis the load is the same as in the first, except for the boundary condition at the TRP, which is not free. Fig. 8 shows the applied Pε against the radial displacement. Analysis 3 is the same as the first two, except for the boundary condition at the TRP, where displacement is restrained in the axial direction. Fig. 9 shows the applied Pε against the computed radial displacement. In the fourth analysis the riser is subject to a pure monotonic increase in Pu linearly with time. The TRP is free to move in all directions. Fig. 10 shows the applied Pu against the computed radial displacement. The fifth analysis corresponds to a cyclic axial tensile load simulation preceded by an initial ‘pressure load step’. The analysis starts with an initial ‘pressure load step’ which accounts for the internal fluid pressure and the external hydrostatic pressure. Four different input pressures are analyzed and referred to as Cases 1, 2, 3, and 4, respectively. The values of the applied internal and external pressures are given in Table 2. The pressure step is then
A. Bahtui et al. / Engineering Structures 32 (2010) 2287–2299
Analysis 3 1000
Pε ( k N )
800
600
400
200
0 0
0.0002
0.0004
0.0006
0.0008
Ur (m) Fig. 9. Result of analysis 3.
Analysis 4 12
Pu (MN/m)
9
6
2295
ternal and external pressure loading, the behavior becomes more stabilized, and a linear curve can be fitted to represent the overall behavior. Fig. 8 shows the radial displacement which resulted from analysis 2. The cyclic response of radial displacement due to cyclic internal and external loading is unlikely to be due to frictional dissipation between adjacent layers, but is potentially due to the gaps closure and opening at the initiation of each cycle. Therefore a linear curve fitting assumption for this figure makes sense. Radial displacements of analysis 3 are shown in Fig. 9. When radial displacements in this figure are compared with that of analysis 2, it can be seen that the latter are much lower and could be ignored, but are instead curve fitted linearly and considered in the constitutive model. Fig. 10 shows the radial displacement of analysis 4. This figure well illustrates the characteristics of the variable Pu . As mentioned in Section 2, an increase in this variable makes an inward or outward overall radial movement of all layers which neither creates any compression and thus frictional dissipation between layers, nor opens and closes the gaps between layers. Therefore the response is clearly linear and a linear curve fitting is well suited to this curve. All results of Fig. 11, except for the case Pε = 0, illustrate a hysteresis behavior of the riser subjected to cyclic axial tension (analysis 5), the response tending towards a stabilized cyclic response increasing gradually with the number of cycles. The first cycle represents the installation phase, whereas the stabilized results are more suitable for use in the calibration procedure, which will be discussed in the next section. This figure also shows that the higher the pressure, the greater the energy dissipated during the hysteresis loop.
6. Calibration of the constitutive model
3
0 0
0.003
0.006
0.009
0.012
Ur (m) Fig. 10. Result of analysis 4.
followed by a ‘tensile step’, in which a 100 kN axial tension is cyclically applied to the free end of the pipe. In each loading or unloading phase of the cycles the load is applied linearly with time. In Fig. 11(a)–(d) the applied axial force is plotted against the computed axial strain for the four different load cases defined in Table 2. Notice that the initial pressure load step results in initial axial strain, radial strain and radial displacement offsets, which are not shown in the figures. As expected, results are similar to those in Fig. 2. The sixth case analysis corresponds to a cyclic bending moment simulation preceded by an initial ‘pressure load step’. Analysis 6 was presented comprehensively by the authors in the previous work [2] and hence is not presented here. 5.3. Discussion Fig. 7 shows the radial strain for analysis 1. The initial jump corresponds to the first time step where all gaps between layers close. These gaps develop between the adjacent layers of the riser during mesh due to the limited number of elements in the circumferential direction of each layer. As the simulation goes through cyclic in-
Several FE simulations have been used to calibrate the constitutive model. Each result from the numerical simulations was curve-fitted either using a linear or bilinear curve, based on the stabilized cycle. The linear curve yields a constant slope (analysis 1–4). The bilinear curve can give three specific parameters for the constitutive model (analysis 5 and 6): the initial no-slip slope, a slip-initiation point, and a full-slip slope. For analyses 1–4 a linear response can be identified. These slopes directly provides the flexural stiffnesses D55 , D15 = D51 and D45 = D54 by using appropriate equations from the fifth column of Table 1. For all cases of Fig. 11, except those with zero or very small applied pressure, an initial approximately linear response can be identified. Therefore, an initial ‘no-slip slope’ can be determined and is found to be the same for each of the analyzed cases. This validates the hypothesis that an initial linear part of the force–displacement diagram, in which frictional slip between layers is absent or negligible, does exist. This slope directly provides the flexural stiffness D11 . After this initial linear response, a nonlinear part of the curve follows, which relatively rapidly tends to join a final full-slip slope. Upon repeated loading and unloading cycles, this final slope changes and stabilizes on a straight line, which is that used in the identification procedure. This final full-slip slope is not equal to the hardening parameter H11 , but it is related to it and therefore its determination allows the estimation of H11 . The intersection point between the initial slope and the final slope provides a ‘slip-initiation’ point, which is also a point on the boundary of the no-slip domain, i.e. a point at the zero level set of the slip-onset function in the N–Pε plane.
2296
A. Bahtui et al. / Engineering Structures 32 (2010) 2287–2299
a
b
d
c
Fig. 11. Axial force–axial strain curves for the four different load cases of analysis 6 and reported in Table 2. Cases are 1(a), 2(b), 3(c) and 4(d).
x 10 4 8
in many cases, it is important to underline that the method can be further refined, as discussed later. The slip-initiation points identified have been fitted with the quadratic expression of Eq. (5) and the resulting surface is shown in Fig. 12. Table 4 reports the parameters of the model determined with the proposed procedure.
N - βz
7 6 5 4 3
7. Verification of the model
2 1 0 5000 0 Mx - βx
-5000 0
1
2
4 3 Pε - βε
5
6 x 10 5
Fig. 12. Curve-fitting of the slip-initiation points.
6.1. Remarks In this proposed model, the nonlinear transitional regime between the initial no-slip response and the final full-slip response is replaced with a bilinear response. During the transition phase, some of the layers in contact (particularly the armor tendons) slip with respect to the others, while the others do not. Although the results of the proposed simplified model can be accurate enough
To evaluate the proposed method, the same tensile cyclic analyses performed using the detailed finite-element model in ABAQUS have been reproduced using the proposed constitutive model adopting the parameters of Table 4 and using the solution scheme described in Section 2. The results of this verification for four cases corresponding to four values of the applied pressure load Pε are reported in Fig. 13. For each case, a graph is reported in which the load–displacement curve obtained using the proposed constitutive model is plotted with a solid line and is compared with the dashed-line curve obtained using the detailed finite-element model. In Case 1 (Pε = 0), no friction occurs in both cases, which results in no hysteresis, i.e. no energy dissipation. The two curves overlap and cannot easily be distinguished. The difference between the stabilized response and that in the first one or two cycles increases with increasing values of Pε , and could be an issue of concern only in the analysis of the installation phase, while in the other cases it can be considered not important for the accuracy of the analysis.
A. Bahtui et al. / Engineering Structures 32 (2010) 2287–2299
2297
Table 4 Identified parameters. Stiffness D11 (N) D41 (N m−1 ) D51 (N) D22 (N m2 ) D44 (N m−2 ) D54 (N m−1 ) D55 (N)
2.5 · 108 1.28 · 108 −5.88 · 106 6.08 · 105 4.38 · 109 −1.36 · 108 1.52 · 107
Slip-onset function a b (N−1 ) c (N−1 m−2 )
0 1.2 · 10−4 2.3 · 10−2
Hardening H11 (N) H22 (N m2 )
7.3 · 107 2.5 · 105
b
a
d
c
Fig. 13. Comparison of the results obtained from the constitutive model (solid line) and from the detailed FE model (dotted line).
For verification related to bending-moment analysis, refer to [2]. 7.1. Discussion of the constitutive model The comparison between the curves provided by the constitutive model and the stabilized cycles obtained using the FE numerical simulations indicate a reasonable agreement in terms of the overall response. The model precisely captures the minimum and maximum magnitudes of all cases of different input pressures. The initial stiffnesses of all cases are well coincident with the stabilized initial slopes of the FE numerical simulations. The hardening stiffnesses of all cases also perfectly capture all the final full-slip slopes of the FE simulations.
It is worth noting that for all cases the same stiffness and hardening parameters have been used, that these parameters have been determined by interpolating the FE results with a bilinear curve only assuming a monotonic case, and that very good agreement is reached for cases with significantly different values of the pressure and for the entire stabilized cycles. This demonstrates that the assumptions on which the constitutive model is based are well justified and that the proposed approach could represent an excellent method of analysis of long flexible risers. The constitutive model is sufficiently robust to capture many practical loadings and has the capability of being easily improved, enhanced and extended. One possible development would be to introduce nonlinear kinematic hardening, which would then be able to match FE simulation or experimental results
2298
A. Bahtui et al. / Engineering Structures 32 (2010) 2287–2299
dσ
dε
= ε=Bˆ NL qˆ kn+1
dσ
dε
n+ 1
2 δθ111 − 4b2 τ12 Γ111 −4bc τ1 τ2 Γ111 ∆222 1 = −4bc τ1 τ3 Γ111 ∆333 δ δθ141 − 4b2 τ 2 Γ111 Γ141 1 δθ151 − 4b2 τ12 Γ111 Γ151
−4bc τ1 τ2 Γ111 ∆222 δ Σ2 − 4c 2 τ22 ∆2222 −4c 2 τ2 τ3 ∆222 ∆333 −4bc τ1 τ2 Γ141 ∆222 −4bc τ1 τ2 Γ151 ∆222
−4bc τ1 τ3 Γ111 ∆333 −4c 2 τ2 τ3 ∆222 ∆333 δ Σ3 − 4c 2 τ32 ∆2333 −4bc τ1 τ3 ∆141 ∆333 −4bc τ1 τ3 ∆151 ∆333
δθ411 − 2bτ1 Γ111 (2bτ1 Γ141 − Γ541 ) −2c τ2 ∆222 (2bτ1 Γ141 − Γ541 ) −2c τ3 ∆333 (2bτ1 Γ141 − Γ541 ) δθ441 − 2bτ1 Γ141 (2bτ1 Γ141 − Γ541 ) δθ451 − 2bτ1 Γ151 (2bτ1 Γ141 − Γ541 )
δθ511 − 2bτ1 Γ111 (2bτ1 Γ151 − Γ551 ) −2c τ2 ∆222 (2bτ1 Γ151 − Γ551 ) −2c τ3 ∆333 (2bτ1 Γ151 − Γ551 ) δθ541 − 2bτ1 Γ141 (2bτ1 Γ151 − Γ551 ) δθ551 − 2bτ1 Γ151 (2bτ1 Γ151 − Γ551 )
Box III.
the behavior properly, this paper presents a more accurate constitutive model for these structures. The constitutive model characterizes the hysteresis behavior of various layers of risers when subjected to different modes of loadings. The constitutive model was incorporated into a generalized finite element formulation which has been developed based on the Euler–Bernoulli beam theory together with two additional degrees of freedom. From the results it is noted that using the constitutive model together with the developed generalized finite element formulation, it is possible to study the nonlinear behavior of long flexible risers with significantly low computational costs. In future work, the constitutive law will be improved by including nonlinear hardening as well as more refined hardening models capable of simulating the stabilization of the response during cyclic loading. Further improvements can be made by accounting for the presence of the shear force (switching to a Timoshenko beam model) and by incorporating the influence of variable axial force and torque on the nonlinear response. Clearly, more sophisticated models will also imply more parameters to be identified and therefore call on the development of more refined identification procedures. Fig. 14. Cantilever un-bonded flexible riser subjected to an end bending moment.
Acknowledgements
much more accurately in terms of the gradual change of stiffness between the no-slip and full-slip phases.
The work has been funded by the UK Engineering and Physical Science research Council (EPSRC) under grant number EP/C510232/1. The authors would like to gratefully acknowledge their support. They would also like to thank Lloyds Register EMEA in London for the precious and fruitful discussions.
8. Large-scale analysis: one mechanical example The situation where a cantilever flexible riser with length 500 m, which is subjected to a 1.5 kN m end bending moment, is studied. Fig. 14 shows the result for this flexible riser for three different conditions: a linear solution with constant flexural stiffness EI = 608 108 Pa, a nonlinear solution with constant EI = 608 108 Pa and a nonlinear solution with the inclusion of the constitutive model with the initial value EI = 608 108 Pa (the first EI value which is taken from the constitutive model’s output). From Fig. 14 it is noted that the presented constitutive model is capable of capturing the energy dissipation in the flexible riser structure due to frictional slip between layers and the hysteresis response. The detailed finite element analysis of a flexible riser using commercial FEA packages on a cluster of parallel processors takes a significant time to arrive at the solution. Using the developed constitutive model and a generalized finite element formulation, as presented in Section 3, the computational time for studying the nonlinear behavior of several kilometers length of un-bonded flexible risers takes only a few seconds. This example reveals the capability and merit of the proposed constitutive model. 9. Conclusions In view of the fact that the equivalent modeling scheme (treating the riser as one homogeneous material layer) for the simulation of the various layers of an un-bonded flexible riser fails to model
Appendix According to the equation given in Box III,
δ = 4b2 τ12 Π11 + 4c 2 τ22 Ω22 + 4c 2 τ32 Ω33 − 2bτ1 Π15 2bλ Dij + Hij Dji + Hji Πij = Dij + Hij − , 1 + 2bλ (Dii + Hii ) 2c λ Dij + Hij Dji + Hji Ωij = Dij + Hij − 1 + 2c λ (Dii + Hii ) 2bλDkj Djk + Hjk , Γijk = Dij − 1 + 2bλ (Dkk + Hkk ) 2c λDkj Djk + Hjk ∆ijk = Dij − 1 + 2c λ (Dkk + Hkk ) 2bλDkj Dik θijk = Dij − , 1 + 2bλ (Dkk + Hkk ) 2c λD2ii Σi = Dii − . 1 + 2c λ (Dii + Dii ) References [1] API. American Petroleum Institute. Recommended Practice for Flexible Pipe. 2002. [2] Alfano G, Bahtui A, Bahai H. Numerical derivation of constitutive models for un-bonded flexible risers. Int J Mech Sci 2009;51(4):295–304.
A. Bahtui et al. / Engineering Structures 32 (2010) 2287–2299 [3] Bahtui A, Bahai H, Alfano G. Numerical and analytical simulation of un-bonded flexible risers subjected to combined modes of loading. In: 9th US national congress on computational mechanics. 2007. [4] Bahtui A, Bahai H, Alfano G. A finite element analysis for un-bonded flexible risers under axial tension. In: Proceedings of OMAE2008 27th international conference on offshore mechanics and arctic engineering. 2008. [5] Bahtui A, Bahai H, Alfano G. A finite element analysis for un-bonded flexible risers under torsion. J Offshore Mech Arctic Eng, OMAE 2008;130(4):041301. 4 pages. [6] Bahtui A, Bahai H, Alfano G. Constitutive modeling of un-bonded flexible risers under tension. In: 8th world congress on computational mechanics. 2008. [7] Bahtui A, Bahai H, Alfano G. Numerical and analytical modeling of un-bonded flexible risers. J Offshore Mech Arctic Eng, OMAE 2009;131(2):1010. [8] Kagoura T, Ishii K, Abe S, Inoue T, Hayashi T, Sakamoto T, et al. Development of a flexible pipe for pipe-in-pipe technology. Furukawa Rev 2003;24:69–75. [9] Felippa CA, Chung JS. Nonlinear static analysis of deep-ocean mining pipe EM DASH 1. Modeling and formulation. Amer Soc Mech Eng 1980;80(48): 6 pages. [10] Felippa CA, Chung JS. Nonlinear static analysis of deep-ocean mining pipe EM DASH 1. Modeling and formulation. Trans ASME, J Energy Resour Technol 1981;103(1):11–5. [11] Felippa CA, Chung JS. Nonlinear static analysis of deep-ocean mining pipe EM DASH 2. Numerical studies. Trans ASME, J Energy Resour Technol 1981;103(1): 16–25. [12] McNamara JF, Lane M. Practical modelling for articulated risers and loading columns. Trans ASME, J Energy Resour Technol 1984;106(4):444–50. [13] McNamara JF, Lane M. 3D frequency domain analysis of offshore structures. In: Proceedings of the 9th conference on engineering mechanics. 1992. p. 192–5. [14] O’Brien PJ, McNamara JF, Dunne FPE. Three-dimensional nonlinear motions of risers and offshore loading towers. J Offshore Mech Arctic Eng 1988;110(3): 232–7. [15] O’Brien PJ, McNamara JF. Analysis of flexible riser systems subject to threedimensional seastate loading. In: Proceedings of the international conference on behaviour of offshore structures, 6. 1988. p. 1373–88.
2299
[16] Hoffman D, Ismail N, Nielsen R, Chandwani R. Design of flexible marine riser in deep and shallow water. In: Proceedings 23rd annual offshore technology conference. 1991. p. 253–65. [17] Atadan AS, Calisal SM, Modi VJ, Guo Y. Analytical and numerical analysis of the dynamics of a marine riser connected to a floating platform. Ocean Eng 1997; 24(2):111–31. [18] Chai YT, Varyani KS, Barlt NDP. Three-dimensional Lump–Mass formulation of a catenary riser with bending, torsion and irregular seabed interaction effect. Ocean Eng 2002;29:1503–25. [19] Ong PPA, Pellegrino S. Modelling of seabed interaction in frequency domain analysis of mooring cables. In: Proceedings of 22nd international conference on offshore mechanics and arctic engineering ASME. 2003. p. 1–9. [20] Zhang Y, Chen B, Qiu L, Hill T, Case M. State of the art analytical tools improve optimization of un-bonded flexible pipes for deepwater environments. In: The 2003 offshore technology conference. 2003. [21] Yazdchi M, Crisfield MA. Buoyancy forces and the 2D finite element analysis of flexible offshore pipes and risers. Internat J Numer Methods Engrg 2002;54: 61–88. [22] Yazdchi M, Crisfield MA. Nonlinear dynamic behavior of flexible marine pipes and risers. Internat J Numer Methods Engrg 2002;54:1265–308. [23] Lacarbonara W, Pacitti A. Nonlinear modeling of cables with flexural stiffness. Math Probl Eng 2008;(370767):1–21. [24] Bech A, Skallerud B, Sodahl N. Structural damping in design analysis of flexible risers. In: Proceedings of the first European conference on flexible pipes, umbilicals and marine cables. 1992. 14 pages. [25] Smith R, O’Brien P, O’Sullivan T, Weibe C. Fatigue analysis of un-bonded flexible risers with irregular seas and hysteresis. In: Offshore technology conference. 2007. [26] Alfano G, Rosati L, Valoroso N. A tangent–secant approach to rate-independent elastoplasticity: formulations and computational issues. Comput Methods Appl Mech Engrg 1999;179:379–405. [27] Han W, Reddy BD. Plasticity—mathematical theory and numerical analysis. New York: Springer-Verlag; 1999 [Chapter 3].