# Strength Analysis Modelling of Flexible Umbilical Members for Marine Structures

Authors: S. Sævik; J. K. Ø. Gjøsteen
Journal: Journal of Applied Mathematics
 (2012)

---

## Abstract

This paper presents a 3-dimensional finite element formulation for predicting the behaviour of complex umbilical cross-sections exposed to loading from tension, torque, internal and external pressure including bending. Helically wound armours and tubes are treated as thin and slender beams formulated within the framework of small strains but large displacements, applying the principle of virtual displacements to obtain finite element equations. Interaction between structural elements is handled by 2- and 3-noded contact elements based on a penalty parameter formulation. The model takes into account a number of features, such as material nonlinearity, gap and friction between individual bodies, and contact with external structures and with a full 3-dimensional description. Numerical studies are presented to validate the model against another model as well as test data.


---

## Body

## 1. Introduction

Flexible risers such as flexible pipes and umbilical cables represent a crucial element of a floating production system. The concept is an attractive alternative to a rigid riser since it does not require heave compensation and tensioning devices at the top or riser manifold at the seabed. At the same time, it offers ease of installation, retrieval and usage elsewhere.Flexible pipes and umbilical cables both work as composite pipes that are compliant and highly deformable in bending but strong and stiff in response to internal pressure, external pressure, tension, and torque.For the flexible pipe, concentric polymeric layers are used to provide sealing. These layers are supported by interlocked metallic layers to resist pressure loading and tensile armour layers to resist tension, torque, and the pressure end cap effect. Umbilicals have homogenous layers, cables, hydraulic hoses, and metal tubes at their central core. Both helical armours and helical multifunction tubes are used to resist tension and torque; see Figure1.Figure 1
Umbilical cross-section.Predicting the structural response in such composite structures is a problem that has been dealt with by many authors. Knapp [1] addressed the axisymmetric response of   aluminium core with steel reinforcement cables. Raoof and Hobbs [2] proposed to use an orthotropic material formulation to treat the helix armours as shells. Costello [3] developed a model for wire ropes that is suitable for cables with few large tendons in lateral contact. Witz and Tan [4] presented expressions for cables and flexible pipes under moderate loads. Computer models for the same range of application on flexible pipes were also presented by McNamara and Harte [5]. Feld [6] investigated the response of power cables, material nonlinearity included, but assuming a constant factor between tension and core radial contraction. Custódio and Vaz [7] studied the axisymmetric response of umbilical cables, introduced improvements to previously published models, and compared estimations with experiments. Sævik and Bruaseth [8] presented a fairly general 2-dimensional model, however, also assuming axisymmetric behaviour.The response of flexible pipes subjected to bending has been less focused on. Lutchansky [9] and Tan et al. [10] allow axial movement of the helical reinforcing only. Ferét and Bournazel [11] suggest that the unbonded helically wound tendon will follow the geodesic of the curved cylinder. Leclair and Costello [12] used Love's equations and assumed wire geometry to calculate the local and global response of wires in a bent rope. Sævik [13] developed a curved beam element to study the stress and deformation in one single armour tendon exposed to a given curvature distribution. Most models presented in the literature, however, assume that the curvature is constant and given from global analyses assuming elastic bending stiffness, excluding the effect of the friction work done due to relative sliding between layers. Sævik et al. [14] presented a model considering the full cross-section in bending for flexible pipes; however, the assumptions applied in this model utilised the benefit from the layered structure of such pipes, thus not very well suited for the umbilical cross-section, often consisting of a few individual structural elements with a significant bending stiffness.The main purpose of the present work has therefore been to formulate a model that can describe stresses and slip of umbilical structural elements for both axisymmetric and bending loads and to validate the model against available test data.


## 2. Model

### 2.1. General

Based on a comprehensive literature review of models, Custódio and Vaz [7] listed the most common assumptions. For clarity and comparison purposes, these assumptions are repeated here using their definition.(1)
Regularity of initial geometry: (a) the homogenous layers are long and uniform cylinders; (b) the wires are wound on a perfectly cylindrical helix; (c) the wires are equally spaced; (d) the wires of an armour are numerous, hence the forces they exert on the adjacent layers may be replaced by uniform pressure; (e) the structure is straight.(2)
Reduction to simple plane analysis: (f) there are no field loads such as self-weight; (g) end effects may be neglected; (h) the material points from any layer have the same longitudinal displacement and twist; (i) all wires of an armour present the same stress state; (j) the wires maintain a helicoidal configuration when strained; (k) the angle between the wire cross-section principal inertia axis and a radial vector linking the centre of structure's cross-section and the centre of wire's cross-section is constant; (l) there is no overpenetration or gap spanning.(3)
The effects of shear and internal friction are neglected: (m) the wires are so slender that the movements of the material points are governed only by their tangent strain and not by the change in curvature.(4)
Linearity of the response: (n) the materials have linear elastic behaviour; (o) the changes in armour radii and pitch angles are linearly small; (p) the wires in one armour never touch laterally or are always in contact; (q) there are no voids between layers nor among cables in the functional core; (r) the homogeneous layers are thin and made of soft material so they simply transfer pressure; (s) the umbilical's core responds linearly to axisymmetric loading; (t) both loading and response are not timedependent.The model presented by Custódio and Vaz avoids the assumptions (n)–(s). The present model, however,  avoids the assumptions (c), (d), (f), (g), (h), (i), (j), (k), (m), (n), (o) and (p)  enabling contact to be handled on component level as well as allowing end effects and friction to be included. However, in order to limit the number of unknown parameters in a 3D model, the governing equations are formulated to only include 1st order helices. It is also assumed that local deformations in the cross-section can be handled by a surface stiffness penalty parameter, so that the structural response can be described by beam theory. The finite element approach has further been selected as it allows structural elements and contact interaction effects to be handled on an individual basis.


### 2.2. Finite Element Formulation

In order to establish the equilibrium equations, the principle of virtual displacement is applied. For an arbitrary equilibrium state(2.1)∫V(σ-σ0):δεdV-∫St⋅δudS=0,
where σ is the stress tensor, σ0 is the initial stress tensor, ε is the strain tensor, t is the surface traction, and u is the displacement vector. σ0 may be obtained from the initial strain tensor by applying the material law. In the case of nonconservative loading such as pressure, the resulting load will change as a function of the area change. Hence, the change in the surface area S of the volume V has to be formulated as a function of the strain. The above equation is used for equilibrium control. However, the equation also needs to be formulated on incremental form to allow equilibrium iterations to be carried out. The incremental form is obtained as
(2.2)∫VC∶Δε∶δεdV+∫Vσ:δΔEdV-∫sΔtdS=0,
where E is the Green-St. Venant strain tensor; see Belytschko et al. [15]. Equation (2.2) gives the incremental equilibrium equation to be used as basis for the stiffness matrix. The first term gives the material stiffness matrix, whereas the second term gives the geometric stiffness matrix.With respect to choice of reference configuration from which deformations are measured, the so calledcorotational ghost  reference formulation has been chosen. The basic idea is to separate the rigid body motion from the local or relative deformation of the element. This is done by attaching a local coordinate system to the element and letting it continuously translate and rotate with the element during deformation; see Horrigmoe and Bergan [16]. The basis for this procedure is by attaching a local coordinate system to each node of the structure. Along the helix, the following transformation describes the position of the local coordinate system, i, relative to the right-handed Cartesian coordinate system, I, located at the cross-section centre; see Figure 2:
(2.3)[i1i2i3]=[cos⁡⁡αsin⁡⁡αcos⁡⁡φsin⁡⁡αsin⁡⁡φsin⁡⁡α-cos⁡⁡αcos⁡⁡φ-cos⁡⁡αsin⁡⁡φ0sin⁡⁡φ-cos⁡⁡φ][I1I2I3],
where α is the lay angle and φ is the angular position.Figure 2
Coordinate system.During deformation the helix nodes will translate and rotate and in order to update the position of the local node base vector system from one configuration,n, to the next, n+1 the assumption of small incremental rotations is applied. This gives the following relation between the updated local node base vector system relative to the global base vectors:
(2.4)iin+1=TiknΔTkjnIj,
where the incremental rotation matrix ΔT is defined by the rotation increments Δθi about the global axes Xi for each node and where the rotation is carried out in subsequent order 1-2-3, which by multiplication of the three associated rotation matrices gives
(2.5)ΔT=[cos⁡⁡Δθ2cos⁡⁡Δθ3cos⁡⁡Δθ1sin⁡⁡Δθ3+sin⁡⁡Δθ1sin⁡⁡Δθ2cos⁡⁡Δθ3A-cos⁡⁡Δθ2sin⁡⁡Δθ3cos⁡⁡Δθ1cos⁡⁡Δθ3-sin⁡⁡Δθ1sin⁡⁡Δθ2sin⁡⁡Δθ3Bsin⁡⁡Δθ2-cos⁡⁡Δθ2sin⁡⁡Δθ1cos⁡⁡Δθ2cos⁡⁡Δθ3],
where 𝔄 denotes (sinΔθ1sinΔθ3-sinΔθ2cosΔθ1cosΔθ3) and 𝔅 denotes (cosΔθ3sinΔθ1+sinΔθ2sinΔθ3cosΔθ1).The element types needed can be divided into two categories:(1)
beam elements to model copper conductors, tubes, fill materials or user defined structural elements,(2)
contact elements to model contact between different element combinations such as contact between core and 1st helix layer, between two helix layers and between outer helix, and other interacting structural elements are examples.With reference to (2.1) and (2.2), the following is needed in order to develop equations for each element type that can be implemented into a computer code:(1)
Kinematics description, that is, a relation between the displacements and rotations and the strains at a material point.(2)
A material law connecting the strain with resulting stresses.(3)
Displacement interpolation, describing the displacement and rotation fields by a number of unknowns on matrix format that can be programmed.In the following the above will be separately addressed.


### 2.3. Beam Kinematics

The beam contribution to the equilibrium equation is established assuming that the Bernoulli-Euler and Navier hypotheses apply. The Green-St. Venant strain tensor is used as strain measure when formulating the incremental equilibrium equations. The second order longitudinal strain term in the Green-St. Venant strain tensor is neglected to avoid shear locking. However, all terms related to coupling between longitudinal strain and torsion are included.The displacements of an arbitrary pointP, defined by local coordinates xi in the cross-section as seen in Figure 2, may be expressed as
(2.6)u1(x1,x2,x3)=u10-x2u20,1-x3u30,1,u2(x1,x2,x3)=u20-x3ω1,u3(x1,x2,x3)=u30+x2ω1.
Introducing the previously mentioned assumptions, the longitudinal Green-St. Venant strain is found to be
(2.7)E11=u10,1-x2u20,11-x3u30,11+12(u20,12+u30,12)+ω1,1(x2u30,1-x3u20,1)+12ω1,12(x22+x32),
where ui is the displacements along the respective axis i and  ωi is the rotation about the same axis i. The subscript “0” quantities represent the displacement field in the respective directions. The shear deformation, γ, is assumed according to St. Venant’s principle assuming the beam to be long and slender with no end section warping as
(2.8)γ=Rω1,1.
The beam deformations resulting from the motion of the nodal ii system are referred to the corotational ghost reference j-system. For an arbitrary equilibrium state the ghost reference system is defined by the following procedure, where the superscript index refers to element end node:
(2.9)j1=dxk|dx⋅dx|Ik,j2=12(i21+i22),j3=j1×j2|j1×j2|.
The local beam deformations at end 1 and referred to as the local ji system, is found by transformation of the ii node system into the ji element system:
(2.10)ii1=Tij1Ij=Tij1Tjkjk=T~ij1jj.
The elemental rotational deformations in node i  are determined by
(2.11)θ1i=12θsin⁡⁡θ(T~23i-T~32i),θ2i=12θsin⁡⁡θ(-T~13i-T~31i),θ3i=12θsin⁡⁡θ(T~12i-T~21i),
where θ represents the resultant rotation at node i.


### 2.4. Contact Kinematics

Reference is given to Figure3.Figure 3
Geometrical relations for contact element.A contact element with two nodes is developed by considering two bodiesA and B1 for this case body B2 is ignored. Each of the bodies occupies a region Vl and has a boundary Sl where l=A or B1. For a virtual time interval [t,t+Δt], the displacement fields are denoted by ul=ul(xl), where xlVl. When bodies A and B1 are brought into contact, let Sc be the unknown contact surface, which satisfies the relationship Sc=VA∩VB1 and Sl=Sσl∪Sul∪Scl, where Sσl denotes part of the surface with prescribed surface tractions Tl and   Sul denotes part of the surface with prescribed surface displacements. Also let n be the outward surface normal vector of body A at x∈Scl and t be the corresponding tangent vector. At the beginning of a time increment, an initial gap g0 at x∈Scl in the direction of n is defined as
(2.12)g0=(xB1-xA)⋅n=[(RB1+rB1)-(RA+rA)]⋅n,
where xl, l=A or B1,  represents the updated coordinates of a point at time t.In order to describe contact between two helical layers, contact is considered between 3 bodies. One body is from the inner layer, and is denotedA. The remaining two bodies are from the layer outside, and are denoted B1 and B2. With reference to Figure 3, the initial gap can be expressed by:
(2.13)g0=(1-ξ)⋅(xB1-xA)⋅n+ξ⋅(xB2-xA)⋅n,
where the non-dimensional parameter ξ is defined as:
(2.14)ξ=(RA+rA-RB1)⋅(RB2-RB1)|RB2-RB1|2.
Note that in this case, the direction of the surface normal vector n is from the center of the umbilical to the center of body A. It is directed outwards, to the next layer containing bodies B1 and B2.The current gap at timet+Δt in the direction of n can be described by:
(2.15)g=(ΔuB1-ΔuA)⋅n+g0≥0
for a two-body contact and by
(2.16)g=(1-ξ)(ΔuB1-ΔuA)⋅n+ξ(ΔuB2-ΔuA)⋅n+g0≥0
for a three-body contact.Two contact conditions may occur:(1)
gap opening, ifg≥0,(2)
contact, ifg<0.Further, if contact has been established, relative slippage including friction work will occur for a two-body contact when(2.17)Δγt=(ΔuB1-ΔuA)⋅t≠0,Δγs=(ΔuB1-ΔuA)⋅s≠0,
where s is pointing along the centre line of body A and t is obtained by taking the cross-product between s and n. The three-body contact is treated in the same manner:
(2.18)Δγt=(1-ξ)(ΔuB1-ΔuA)⋅t+ξ(ΔuB2-ΔuA)⋅t≠0,Δγs=(1-ξ)(ΔuB1-ΔuA)⋅s+ξ(ΔuB2-ΔuA)⋅s≠0.


#### 2.4.1. Patch Test

The meshing in the longitudinal direction is restrained to be identical for all the helical components. Hence, the 2-noded contact elements operate between coincident nodes, and the contact element will pass the patch test [17]. For the 3-noded contact elements, the force balance will be satisfied within each separate contact element. A contact force that is acting on body A, Figure 3, will have a reaction force that will be linearly distributed with ξ between body B1 and B2. The patch test is therefore also passed for this contact element.


### 2.5. Material Models

#### 2.5.1. Beams and Tubes

Denoting the longitudinal direction by index 1 the circumferential direction by index 2 and assuming that the normal direction is governed by a prescribed value,p, Hooke’s law for linear elastic materials reads:
(2.19)[σ11σ22τ12]=ν1-ν[pp0]+E1-ν2[1ν0ν10001-ν22(1+ν)][ε11ε22γ12].
For the beam and tube elements (2.19) is further simplified by prescribing σ22  by application of thin shell theory for the tube element.When the stresses exceed the elastic limit of the stress/strain relation, an elastoplastic formulation is required. In this case, it is of great importance to use an elastoplastic formulation that takes into account the two-dimensional stress state, that is, both the stresses in the axial and hoop directions for the tubes. The applied elastoplastic material model is based on expressing(i)
the yield criterion as a yield surface specified by the scalar functionf,(ii)
the 2nd Piola-Kirchhoff stress tensorSas a measure for stress together with the energy conjugate Green strain tensor E as a measure of strain.J
2-flow theory of plasticity is applied, assuming that the yielding is independent of the first and third deviatoric stress invariant, and any combination of kinematic and isotropic hardening is allowed for. The constitutive relation reads:
(2.20)S˙=[C(e)-αC(e):(∂f/∂S)(∂f/∂S):C(e)(∂f/∂S):C(e):(∂f/∂S)-(∂f/∂Seq)(dSeq/dW(p))S:(∂f/∂S)]:E˙,
where
(2.21)α=1if∂f∂S:S˙≥0,f=0,α=0if∂f∂S:S˙≥0orf<0,
where C(e) is the elastic material law, W(p) is the plastic work done, and Seq is the equivalent stress. For further details, references are made to Levold [18], McMeeking and Rice [19] and Belytschko et al. [15].


#### 2.5.2. Contact Element

In nonlinear finite element analysis, there are three commonly used principles when dealing with contact problems see Shyu et al. [20, 21]. These are(i)
the lagrange multiplier method (LM),(ii)
the penalty method (PM),(iii)
the mixed method (MM).In LM, the constraint conditions for a contact problem are satisfied by introducing Lagrange parameters in the variation statement, where both the displacements and Lagrange multipliers are treated as unknowns which again leads to an increased number of finite element equations.In PM, the contact pressure is assumed proportional to the amount of penetration by introducing a pointwise penalty parameter. The final stiffness matrix does not contain additional terms. However, since the resulting contact forces are of the same order as the assumed displacement field, this may lead to violation of the local Babuska-Brezzi stability condition, see Chang et al. [20], and the success of using MM is highly dependent on the selected order of the contact pressure. A guideline on this item is given in [20] where an MM contact element with excellent numerical properties was presented.In this case, however, the focus has been on formulating a contact element that allows easy and flexible modelling of contact effects in helical structures. Therefore a node-to-node contact formulation based on PM has been adopted.The constitutive relation used to model friction in the contact elements consists of two major ingredients:(i)
a friction surface,(ii)
a slip rule.The constitutive model used is established based on the work by Shyu et al. [21] but is adjusted to include material hardening to improve numerical stability. The friction surface fs is assumed to be a function of the normal traction λn, the tangential tractions λs and λt, and the friction coefficient μ, that is:
(2.22)fs=fs(λn,λs,λt,μ).
Assuming a linear surface, we have
(2.23)fs=λs2+λt2-μλn=0.
The slip increment is divided into two parts:
(2.24)Δγ=Δγe+Δγp,
where Δγe is the 2-dimensional increment of elastic slip and Δγp is the corresponding increment of plastic slip. The elastic slip increment is determined as
(2.25)Δγ-e=Δγs2+Δγt2e=Δλs2+Δλt2eE=Δλ-eE,
where E is a proportional constant. It is seen that as E grows towards infinity, a classical Coulomb law is obtained. Assuming symmetry, that is, an associative slip rule, the following slip rule is postulated:
(2.26)γp=η[λsλt],
where η is a still unknown proportional constant. By assuming a constant normal traction, the proportional constant is determined from the consistency condition:
(2.27)dfs=∇λλ-:Δλ-ηλn∇γpμ:∇λλ-=0,
where the last term is associated with material hardening. The expressions are therefore developed including the hardening parameter β, that is:
(2.28)∂μ∂γ-p=β.By combining the previous equations the following constitutive relation is obtained:(2.29)[ΔλsΔλt]=E[1-λs2λ-2(1+β/E)-λsλtλ-2(1+β/E)-λsλtλ-2(1+β/E)1-λt2λ-2(1+β/E)][ΔγsΔγt].


### 2.6. Solution Procedure

The developed finite element equations have been implemented into a computer code using standard procedures for matrix operations. For the beam elements, linear interpolation is used in the axial direction, whereas cubic interpolation has been applied in the transverse direction. For the contact elements, linear interpolation is used for the 3-noded element used between helical layers to allow the helix to distribute forces if the helix is positioned between two bodies in the next layer. Between umbilical helix and 1st helix layer and between last helix layer and external structures, two-noded contact elements are used, and hence no interpolation is applied for these elements. With respect to numerical integration, closed form expressions are developed for the beam elastic case; that is, no numerical integration is needed. In the elastoplastic case, for pipe elements numerical integration is performed by dividing the beam into a number of subvolumes each characterized by its strain/stress state.The numerical procedure used to solve the finite element equations was based on a Newton-Rapson iterative scheme controlled by applying displacement, force, and energy norms to ensure convergence, see Belytshcko et al. [15].


## 2.1. General

Based on a comprehensive literature review of models, Custódio and Vaz [7] listed the most common assumptions. For clarity and comparison purposes, these assumptions are repeated here using their definition.(1)
Regularity of initial geometry: (a) the homogenous layers are long and uniform cylinders; (b) the wires are wound on a perfectly cylindrical helix; (c) the wires are equally spaced; (d) the wires of an armour are numerous, hence the forces they exert on the adjacent layers may be replaced by uniform pressure; (e) the structure is straight.(2)
Reduction to simple plane analysis: (f) there are no field loads such as self-weight; (g) end effects may be neglected; (h) the material points from any layer have the same longitudinal displacement and twist; (i) all wires of an armour present the same stress state; (j) the wires maintain a helicoidal configuration when strained; (k) the angle between the wire cross-section principal inertia axis and a radial vector linking the centre of structure's cross-section and the centre of wire's cross-section is constant; (l) there is no overpenetration or gap spanning.(3)
The effects of shear and internal friction are neglected: (m) the wires are so slender that the movements of the material points are governed only by their tangent strain and not by the change in curvature.(4)
Linearity of the response: (n) the materials have linear elastic behaviour; (o) the changes in armour radii and pitch angles are linearly small; (p) the wires in one armour never touch laterally or are always in contact; (q) there are no voids between layers nor among cables in the functional core; (r) the homogeneous layers are thin and made of soft material so they simply transfer pressure; (s) the umbilical's core responds linearly to axisymmetric loading; (t) both loading and response are not timedependent.The model presented by Custódio and Vaz avoids the assumptions (n)–(s). The present model, however,  avoids the assumptions (c), (d), (f), (g), (h), (i), (j), (k), (m), (n), (o) and (p)  enabling contact to be handled on component level as well as allowing end effects and friction to be included. However, in order to limit the number of unknown parameters in a 3D model, the governing equations are formulated to only include 1st order helices. It is also assumed that local deformations in the cross-section can be handled by a surface stiffness penalty parameter, so that the structural response can be described by beam theory. The finite element approach has further been selected as it allows structural elements and contact interaction effects to be handled on an individual basis.


## 2.2. Finite Element Formulation

In order to establish the equilibrium equations, the principle of virtual displacement is applied. For an arbitrary equilibrium state(2.1)∫V(σ-σ0):δεdV-∫St⋅δudS=0,
where σ is the stress tensor, σ0 is the initial stress tensor, ε is the strain tensor, t is the surface traction, and u is the displacement vector. σ0 may be obtained from the initial strain tensor by applying the material law. In the case of nonconservative loading such as pressure, the resulting load will change as a function of the area change. Hence, the change in the surface area S of the volume V has to be formulated as a function of the strain. The above equation is used for equilibrium control. However, the equation also needs to be formulated on incremental form to allow equilibrium iterations to be carried out. The incremental form is obtained as
(2.2)∫VC∶Δε∶δεdV+∫Vσ:δΔEdV-∫sΔtdS=0,
where E is the Green-St. Venant strain tensor; see Belytschko et al. [15]. Equation (2.2) gives the incremental equilibrium equation to be used as basis for the stiffness matrix. The first term gives the material stiffness matrix, whereas the second term gives the geometric stiffness matrix.With respect to choice of reference configuration from which deformations are measured, the so calledcorotational ghost  reference formulation has been chosen. The basic idea is to separate the rigid body motion from the local or relative deformation of the element. This is done by attaching a local coordinate system to the element and letting it continuously translate and rotate with the element during deformation; see Horrigmoe and Bergan [16]. The basis for this procedure is by attaching a local coordinate system to each node of the structure. Along the helix, the following transformation describes the position of the local coordinate system, i, relative to the right-handed Cartesian coordinate system, I, located at the cross-section centre; see Figure 2:
(2.3)[i1i2i3]=[cos⁡⁡αsin⁡⁡αcos⁡⁡φsin⁡⁡αsin⁡⁡φsin⁡⁡α-cos⁡⁡αcos⁡⁡φ-cos⁡⁡αsin⁡⁡φ0sin⁡⁡φ-cos⁡⁡φ][I1I2I3],
where α is the lay angle and φ is the angular position.Figure 2
Coordinate system.During deformation the helix nodes will translate and rotate and in order to update the position of the local node base vector system from one configuration,n, to the next, n+1 the assumption of small incremental rotations is applied. This gives the following relation between the updated local node base vector system relative to the global base vectors:
(2.4)iin+1=TiknΔTkjnIj,
where the incremental rotation matrix ΔT is defined by the rotation increments Δθi about the global axes Xi for each node and where the rotation is carried out in subsequent order 1-2-3, which by multiplication of the three associated rotation matrices gives
(2.5)ΔT=[cos⁡⁡Δθ2cos⁡⁡Δθ3cos⁡⁡Δθ1sin⁡⁡Δθ3+sin⁡⁡Δθ1sin⁡⁡Δθ2cos⁡⁡Δθ3A-cos⁡⁡Δθ2sin⁡⁡Δθ3cos⁡⁡Δθ1cos⁡⁡Δθ3-sin⁡⁡Δθ1sin⁡⁡Δθ2sin⁡⁡Δθ3Bsin⁡⁡Δθ2-cos⁡⁡Δθ2sin⁡⁡Δθ1cos⁡⁡Δθ2cos⁡⁡Δθ3],
where 𝔄 denotes (sinΔθ1sinΔθ3-sinΔθ2cosΔθ1cosΔθ3) and 𝔅 denotes (cosΔθ3sinΔθ1+sinΔθ2sinΔθ3cosΔθ1).The element types needed can be divided into two categories:(1)
beam elements to model copper conductors, tubes, fill materials or user defined structural elements,(2)
contact elements to model contact between different element combinations such as contact between core and 1st helix layer, between two helix layers and between outer helix, and other interacting structural elements are examples.With reference to (2.1) and (2.2), the following is needed in order to develop equations for each element type that can be implemented into a computer code:(1)
Kinematics description, that is, a relation between the displacements and rotations and the strains at a material point.(2)
A material law connecting the strain with resulting stresses.(3)
Displacement interpolation, describing the displacement and rotation fields by a number of unknowns on matrix format that can be programmed.In the following the above will be separately addressed.


## 2.3. Beam Kinematics

The beam contribution to the equilibrium equation is established assuming that the Bernoulli-Euler and Navier hypotheses apply. The Green-St. Venant strain tensor is used as strain measure when formulating the incremental equilibrium equations. The second order longitudinal strain term in the Green-St. Venant strain tensor is neglected to avoid shear locking. However, all terms related to coupling between longitudinal strain and torsion are included.The displacements of an arbitrary pointP, defined by local coordinates xi in the cross-section as seen in Figure 2, may be expressed as
(2.6)u1(x1,x2,x3)=u10-x2u20,1-x3u30,1,u2(x1,x2,x3)=u20-x3ω1,u3(x1,x2,x3)=u30+x2ω1.
Introducing the previously mentioned assumptions, the longitudinal Green-St. Venant strain is found to be
(2.7)E11=u10,1-x2u20,11-x3u30,11+12(u20,12+u30,12)+ω1,1(x2u30,1-x3u20,1)+12ω1,12(x22+x32),
where ui is the displacements along the respective axis i and  ωi is the rotation about the same axis i. The subscript “0” quantities represent the displacement field in the respective directions. The shear deformation, γ, is assumed according to St. Venant’s principle assuming the beam to be long and slender with no end section warping as
(2.8)γ=Rω1,1.
The beam deformations resulting from the motion of the nodal ii system are referred to the corotational ghost reference j-system. For an arbitrary equilibrium state the ghost reference system is defined by the following procedure, where the superscript index refers to element end node:
(2.9)j1=dxk|dx⋅dx|Ik,j2=12(i21+i22),j3=j1×j2|j1×j2|.
The local beam deformations at end 1 and referred to as the local ji system, is found by transformation of the ii node system into the ji element system:
(2.10)ii1=Tij1Ij=Tij1Tjkjk=T~ij1jj.
The elemental rotational deformations in node i  are determined by
(2.11)θ1i=12θsin⁡⁡θ(T~23i-T~32i),θ2i=12θsin⁡⁡θ(-T~13i-T~31i),θ3i=12θsin⁡⁡θ(T~12i-T~21i),
where θ represents the resultant rotation at node i.


## 2.4. Contact Kinematics

Reference is given to Figure3.Figure 3
Geometrical relations for contact element.A contact element with two nodes is developed by considering two bodiesA and B1 for this case body B2 is ignored. Each of the bodies occupies a region Vl and has a boundary Sl where l=A or B1. For a virtual time interval [t,t+Δt], the displacement fields are denoted by ul=ul(xl), where xlVl. When bodies A and B1 are brought into contact, let Sc be the unknown contact surface, which satisfies the relationship Sc=VA∩VB1 and Sl=Sσl∪Sul∪Scl, where Sσl denotes part of the surface with prescribed surface tractions Tl and   Sul denotes part of the surface with prescribed surface displacements. Also let n be the outward surface normal vector of body A at x∈Scl and t be the corresponding tangent vector. At the beginning of a time increment, an initial gap g0 at x∈Scl in the direction of n is defined as
(2.12)g0=(xB1-xA)⋅n=[(RB1+rB1)-(RA+rA)]⋅n,
where xl, l=A or B1,  represents the updated coordinates of a point at time t.In order to describe contact between two helical layers, contact is considered between 3 bodies. One body is from the inner layer, and is denotedA. The remaining two bodies are from the layer outside, and are denoted B1 and B2. With reference to Figure 3, the initial gap can be expressed by:
(2.13)g0=(1-ξ)⋅(xB1-xA)⋅n+ξ⋅(xB2-xA)⋅n,
where the non-dimensional parameter ξ is defined as:
(2.14)ξ=(RA+rA-RB1)⋅(RB2-RB1)|RB2-RB1|2.
Note that in this case, the direction of the surface normal vector n is from the center of the umbilical to the center of body A. It is directed outwards, to the next layer containing bodies B1 and B2.The current gap at timet+Δt in the direction of n can be described by:
(2.15)g=(ΔuB1-ΔuA)⋅n+g0≥0
for a two-body contact and by
(2.16)g=(1-ξ)(ΔuB1-ΔuA)⋅n+ξ(ΔuB2-ΔuA)⋅n+g0≥0
for a three-body contact.Two contact conditions may occur:(1)
gap opening, ifg≥0,(2)
contact, ifg<0.Further, if contact has been established, relative slippage including friction work will occur for a two-body contact when(2.17)Δγt=(ΔuB1-ΔuA)⋅t≠0,Δγs=(ΔuB1-ΔuA)⋅s≠0,
where s is pointing along the centre line of body A and t is obtained by taking the cross-product between s and n. The three-body contact is treated in the same manner:
(2.18)Δγt=(1-ξ)(ΔuB1-ΔuA)⋅t+ξ(ΔuB2-ΔuA)⋅t≠0,Δγs=(1-ξ)(ΔuB1-ΔuA)⋅s+ξ(ΔuB2-ΔuA)⋅s≠0.


### 2.4.1. Patch Test

The meshing in the longitudinal direction is restrained to be identical for all the helical components. Hence, the 2-noded contact elements operate between coincident nodes, and the contact element will pass the patch test [17]. For the 3-noded contact elements, the force balance will be satisfied within each separate contact element. A contact force that is acting on body A, Figure 3, will have a reaction force that will be linearly distributed with ξ between body B1 and B2. The patch test is therefore also passed for this contact element.


## 2.4.1. Patch Test

The meshing in the longitudinal direction is restrained to be identical for all the helical components. Hence, the 2-noded contact elements operate between coincident nodes, and the contact element will pass the patch test [17]. For the 3-noded contact elements, the force balance will be satisfied within each separate contact element. A contact force that is acting on body A, Figure 3, will have a reaction force that will be linearly distributed with ξ between body B1 and B2. The patch test is therefore also passed for this contact element.


## 2.5. Material Models

### 2.5.1. Beams and Tubes

Denoting the longitudinal direction by index 1 the circumferential direction by index 2 and assuming that the normal direction is governed by a prescribed value,p, Hooke’s law for linear elastic materials reads:
(2.19)[σ11σ22τ12]=ν1-ν[pp0]+E1-ν2[1ν0ν10001-ν22(1+ν)][ε11ε22γ12].
For the beam and tube elements (2.19) is further simplified by prescribing σ22  by application of thin shell theory for the tube element.When the stresses exceed the elastic limit of the stress/strain relation, an elastoplastic formulation is required. In this case, it is of great importance to use an elastoplastic formulation that takes into account the two-dimensional stress state, that is, both the stresses in the axial and hoop directions for the tubes. The applied elastoplastic material model is based on expressing(i)
the yield criterion as a yield surface specified by the scalar functionf,(ii)
the 2nd Piola-Kirchhoff stress tensorSas a measure for stress together with the energy conjugate Green strain tensor E as a measure of strain.J
2-flow theory of plasticity is applied, assuming that the yielding is independent of the first and third deviatoric stress invariant, and any combination of kinematic and isotropic hardening is allowed for. The constitutive relation reads:
(2.20)S˙=[C(e)-αC(e):(∂f/∂S)(∂f/∂S):C(e)(∂f/∂S):C(e):(∂f/∂S)-(∂f/∂Seq)(dSeq/dW(p))S:(∂f/∂S)]:E˙,
where
(2.21)α=1if∂f∂S:S˙≥0,f=0,α=0if∂f∂S:S˙≥0orf<0,
where C(e) is the elastic material law, W(p) is the plastic work done, and Seq is the equivalent stress. For further details, references are made to Levold [18], McMeeking and Rice [19] and Belytschko et al. [15].


### 2.5.2. Contact Element

In nonlinear finite element analysis, there are three commonly used principles when dealing with contact problems see Shyu et al. [20, 21]. These are(i)
the lagrange multiplier method (LM),(ii)
the penalty method (PM),(iii)
the mixed method (MM).In LM, the constraint conditions for a contact problem are satisfied by introducing Lagrange parameters in the variation statement, where both the displacements and Lagrange multipliers are treated as unknowns which again leads to an increased number of finite element equations.In PM, the contact pressure is assumed proportional to the amount of penetration by introducing a pointwise penalty parameter. The final stiffness matrix does not contain additional terms. However, since the resulting contact forces are of the same order as the assumed displacement field, this may lead to violation of the local Babuska-Brezzi stability condition, see Chang et al. [20], and the success of using MM is highly dependent on the selected order of the contact pressure. A guideline on this item is given in [20] where an MM contact element with excellent numerical properties was presented.In this case, however, the focus has been on formulating a contact element that allows easy and flexible modelling of contact effects in helical structures. Therefore a node-to-node contact formulation based on PM has been adopted.The constitutive relation used to model friction in the contact elements consists of two major ingredients:(i)
a friction surface,(ii)
a slip rule.The constitutive model used is established based on the work by Shyu et al. [21] but is adjusted to include material hardening to improve numerical stability. The friction surface fs is assumed to be a function of the normal traction λn, the tangential tractions λs and λt, and the friction coefficient μ, that is:
(2.22)fs=fs(λn,λs,λt,μ).
Assuming a linear surface, we have
(2.23)fs=λs2+λt2-μλn=0.
The slip increment is divided into two parts:
(2.24)Δγ=Δγe+Δγp,
where Δγe is the 2-dimensional increment of elastic slip and Δγp is the corresponding increment of plastic slip. The elastic slip increment is determined as
(2.25)Δγ-e=Δγs2+Δγt2e=Δλs2+Δλt2eE=Δλ-eE,
where E is a proportional constant. It is seen that as E grows towards infinity, a classical Coulomb law is obtained. Assuming symmetry, that is, an associative slip rule, the following slip rule is postulated:
(2.26)γp=η[λsλt],
where η is a still unknown proportional constant. By assuming a constant normal traction, the proportional constant is determined from the consistency condition:
(2.27)dfs=∇λλ-:Δλ-ηλn∇γpμ:∇λλ-=0,
where the last term is associated with material hardening. The expressions are therefore developed including the hardening parameter β, that is:
(2.28)∂μ∂γ-p=β.By combining the previous equations the following constitutive relation is obtained:(2.29)[ΔλsΔλt]=E[1-λs2λ-2(1+β/E)-λsλtλ-2(1+β/E)-λsλtλ-2(1+β/E)1-λt2λ-2(1+β/E)][ΔγsΔγt].


## 2.5.1. Beams and Tubes

Denoting the longitudinal direction by index 1 the circumferential direction by index 2 and assuming that the normal direction is governed by a prescribed value,p, Hooke’s law for linear elastic materials reads:
(2.19)[σ11σ22τ12]=ν1-ν[pp0]+E1-ν2[1ν0ν10001-ν22(1+ν)][ε11ε22γ12].
For the beam and tube elements (2.19) is further simplified by prescribing σ22  by application of thin shell theory for the tube element.When the stresses exceed the elastic limit of the stress/strain relation, an elastoplastic formulation is required. In this case, it is of great importance to use an elastoplastic formulation that takes into account the two-dimensional stress state, that is, both the stresses in the axial and hoop directions for the tubes. The applied elastoplastic material model is based on expressing(i)
the yield criterion as a yield surface specified by the scalar functionf,(ii)
the 2nd Piola-Kirchhoff stress tensorSas a measure for stress together with the energy conjugate Green strain tensor E as a measure of strain.J
2-flow theory of plasticity is applied, assuming that the yielding is independent of the first and third deviatoric stress invariant, and any combination of kinematic and isotropic hardening is allowed for. The constitutive relation reads:
(2.20)S˙=[C(e)-αC(e):(∂f/∂S)(∂f/∂S):C(e)(∂f/∂S):C(e):(∂f/∂S)-(∂f/∂Seq)(dSeq/dW(p))S:(∂f/∂S)]:E˙,
where
(2.21)α=1if∂f∂S:S˙≥0,f=0,α=0if∂f∂S:S˙≥0orf<0,
where C(e) is the elastic material law, W(p) is the plastic work done, and Seq is the equivalent stress. For further details, references are made to Levold [18], McMeeking and Rice [19] and Belytschko et al. [15].


## 2.5.2. Contact Element

In nonlinear finite element analysis, there are three commonly used principles when dealing with contact problems see Shyu et al. [20, 21]. These are(i)
the lagrange multiplier method (LM),(ii)
the penalty method (PM),(iii)
the mixed method (MM).In LM, the constraint conditions for a contact problem are satisfied by introducing Lagrange parameters in the variation statement, where both the displacements and Lagrange multipliers are treated as unknowns which again leads to an increased number of finite element equations.In PM, the contact pressure is assumed proportional to the amount of penetration by introducing a pointwise penalty parameter. The final stiffness matrix does not contain additional terms. However, since the resulting contact forces are of the same order as the assumed displacement field, this may lead to violation of the local Babuska-Brezzi stability condition, see Chang et al. [20], and the success of using MM is highly dependent on the selected order of the contact pressure. A guideline on this item is given in [20] where an MM contact element with excellent numerical properties was presented.In this case, however, the focus has been on formulating a contact element that allows easy and flexible modelling of contact effects in helical structures. Therefore a node-to-node contact formulation based on PM has been adopted.The constitutive relation used to model friction in the contact elements consists of two major ingredients:(i)
a friction surface,(ii)
a slip rule.The constitutive model used is established based on the work by Shyu et al. [21] but is adjusted to include material hardening to improve numerical stability. The friction surface fs is assumed to be a function of the normal traction λn, the tangential tractions λs and λt, and the friction coefficient μ, that is:
(2.22)fs=fs(λn,λs,λt,μ).
Assuming a linear surface, we have
(2.23)fs=λs2+λt2-μλn=0.
The slip increment is divided into two parts:
(2.24)Δγ=Δγe+Δγp,
where Δγe is the 2-dimensional increment of elastic slip and Δγp is the corresponding increment of plastic slip. The elastic slip increment is determined as
(2.25)Δγ-e=Δγs2+Δγt2e=Δλs2+Δλt2eE=Δλ-eE,
where E is a proportional constant. It is seen that as E grows towards infinity, a classical Coulomb law is obtained. Assuming symmetry, that is, an associative slip rule, the following slip rule is postulated:
(2.26)γp=η[λsλt],
where η is a still unknown proportional constant. By assuming a constant normal traction, the proportional constant is determined from the consistency condition:
(2.27)dfs=∇λλ-:Δλ-ηλn∇γpμ:∇λλ-=0,
where the last term is associated with material hardening. The expressions are therefore developed including the hardening parameter β, that is:
(2.28)∂μ∂γ-p=β.By combining the previous equations the following constitutive relation is obtained:(2.29)[ΔλsΔλt]=E[1-λs2λ-2(1+β/E)-λsλtλ-2(1+β/E)-λsλtλ-2(1+β/E)1-λt2λ-2(1+β/E)][ΔγsΔγt].


## 2.6. Solution Procedure

The developed finite element equations have been implemented into a computer code using standard procedures for matrix operations. For the beam elements, linear interpolation is used in the axial direction, whereas cubic interpolation has been applied in the transverse direction. For the contact elements, linear interpolation is used for the 3-noded element used between helical layers to allow the helix to distribute forces if the helix is positioned between two bodies in the next layer. Between umbilical helix and 1st helix layer and between last helix layer and external structures, two-noded contact elements are used, and hence no interpolation is applied for these elements. With respect to numerical integration, closed form expressions are developed for the beam elastic case; that is, no numerical integration is needed. In the elastoplastic case, for pipe elements numerical integration is performed by dividing the beam into a number of subvolumes each characterized by its strain/stress state.The numerical procedure used to solve the finite element equations was based on a Newton-Rapson iterative scheme controlled by applying displacement, force, and energy norms to ensure convergence, see Belytshcko et al. [15].


## 3. Case Study

The purpose of the case study was to investigate the performance of the developed model compared to others. The case selected is the umbilical test case presented by Lutchansky [9], and the basic properties are summarized in Table 1.Table 1
Custódio and Vaz. umbilical cross-section.
Layer
Dext (mm)
Features
External sheath
94.0
Material HDPE,E=720 MPa, ν=0.42
Outer armour
83.8
N
=
56, round wires  d=4.1 mm, lay angle −20 deg.
Inner armour
75.6
N
=
50, round wires   d=4.1 mm, lay angle = 20 deg.
Intermediate polymer
67.4
Material HDPE,E=720 MPa, ν=0.42
8 hose bundle
52.4
N
=
8, telescopic settlement
Soft fller
27.0
Polyurethane,E=720 MPa, ν=0.420
Central hose
12.7
N
=
1, lay angel 0 deg.The contact between the beam and the armour was modelled by setting the surface stiffness penalty parameter to 7 MPa/mm, except for the armour to armour contact, where 100 MPa/mm was applied. Table2 shows the performance of the present model as compared to test data, the Custódio and Vaz model, and the 2-dimensional model presented by Tan et al. [10]. It is seen that the predicted values are similar to the other model values with respect to stiffness parameters. The present 3-dimensional model suggests stronger torsion coupling. However, few data points on coupling parameters were included in the measurements for this case study.Table 2
Comparison with Custódio and Vaz test data, simulation data from Custido and Vaz, and simulated results from Sævik and Bruaseth.
Layer
Test results
Custódio and Vaz
Sævik and Bruaseth
Present model
EA under tension, MN
71–101
82.9
100
82.9
Coupling  (m)(axial strain/torsion)
No data
0.19
0.21
0.11
GJ clockwise (kNm2) (tension = 0 kN)
14.5–56.3
44.7
44.3
44.4
Coupling clockwise (rad/m)(torsion/axial strain)
No data
182
188.0
202
GJ anticlockwise (kNm2) (tension = 0 kN)
15.9–17.22
19.1
19.5
16.5
Coupling anticlockwise (rad/m) (torsion/axial strain)
No data
−406
−330.0
−190For the test cases applied in this study a linear convergence was observed for the complete model.


## 4. Experimental Studies and Analyses

The experimental studies were performed for validation of a 2D umbilical model and have been reported in this context by Tan et al. [10].Different umbilical cross-sections was tested with respect to the axial stiffness under tension, coupling between axial strain and torsion and torsion stiffness in both directions. The different umbilical cross-sections cover the range from torsion balanced dynamic umbilicals to torsion unbalanced static umbilicals and with variable complexity. The cross-section details required to define input to numerical models are attached in the appendix.Measurements, predictions, and comparison between experimental data and model results are presented in Table3. In all analyses a surface stiffness parameter of 100 MPa/mm has been applied to all metal-to-metal interfaces. For the other surfaces, a stiffness parameter of 7 MPa/mm has been applied. An exception was done for cross Section 1, which contains layers of fillers creating a soft core. For this cross section the 7 MPa/mm surface stiffness value was reduced to 3.5 MPa/mm.  When obtaining the torsional stiffness and corresponding coupling factors, 100 kN tension was applied in the analyses. The numerical model was 1.5 m long, and the element lengths were 1.5 cm. A typical element meshing is shown in Figure 4.Table 3
Comparing full-scale tests to analysis. Measurements denoted (M), simulation results (S) and measurement to simulation ratio (R). Average value and standard deviation are given for the ratio parameter.
Umbilical no.
Avrg.
Std.
1
2
3
4
5
6
7
8
9
EA under tension (MN)torsion fixed
M
609
—
—
206
194
209
466
136
366
—
—
S
630
—
—
202
191
232
529
168
406
—
—
R
0.97
—
—
1.02
1.02
0.90
0.88
0.81
0.90
0.93
0.08
EA under tension (MN)torsion free
M
565
466
261
190
—
201
474
114
366
—
—
S
606
538
281
172
—
215
529
141
406
—
—
R
0.93
0.87
0.93
1.10
—
0.93
0.90
0.81
0.90
0.92
0.08
Coupling  (rad/m)(torsion/axial strain)
M
—
−39
9.5
35.2
—
—
—
—
—
—
—
S
—
−36
7.1
38.2
—
—
—
—
—
—
—
R
—
1.08
1.34
0.92
—
—
—
—
—
1.11
0.21
GJ clockwise  (kNm2)
M
339
48
65
12.9
4.7
35.5
50.5
3.85
48.5
—
—
S
288
41
64
14.1
4.2
32.6
48.8
3.94
52.9
—
—
R
1.18
1.17
1.02
0.91
1.12
1.09
1.03
0.98
0.92
1.05
0.10
Coupling clockwise (10-3 m/rad) (axial strain/torsion)
M
—
−3.9
2.2
3.1
2.4
—
—
2.3
—
—
—
S
—
−3.0
1.8
3.5
2.8
—
—
2.2
—
—
—
R
—
1.30
1.22
0.89
0.97
—
—
1.05
—
1.09
0.17Figure 4
Typical finite element model.An example of measured and predicted data is shown in Figure5, where the axial force is plotted against axial strain for umbilical 1 and umbilical 4. Figure 5 shows three graphs, measured, computed, and analytical where the analytical is based on simply adding together the axial stiffness of all structural elements, taking the lay angle into account but neglecting the effect of radial interaction. It is seen that good correlation is obtained between the measured and the computed results in terms of stiffness. For umbilical 1, the analytical stiffness is considerable higher than the simulations and the tests. This is likely to be due to radial interaction effects, which are included in the model, but cannot be taken into account in the analytical solution. Umbilical 1 has a soft core, and radial contraction is therefore more significant than for a stiffer cross-section, as umbilical 4.Example of test results compared to analysis and theoretically obtained values.
(a)
(b)It is seen that the measured stiffness is low for the first load cycle of umbilical 1, which most probably is due to gap closure.In addition to test and analysis results, Table3 gives a summary of the relation tested value/predicted value for all umbilical cross-sections for all test parameters, including the statistics in terms of average correlation and standard deviation.It is seen that in terms of axial and torsion stiffness good correlation is found, with a standard deviation within a 10% range. This must be considered reasonable for such complex structures. It is noted that the axial stiffness in average is predicted to be 7-8% lower than the measurements. With regard to tension/torsion coupling the number of valid measurements was limited, and therefore the statistical basis is too weak to conclude. This was due to the fact that some of the specimens were torsion balanced; hence very small inherent test rig friction forces would influence the measured rotation. The reported values include only torsion unbalanced cross-sections where the amount of rotation is larger. For these cross-sections reasonable correlation was found with a standard deviation within a 20% range.It should be noted that a better match of analysis to the experiment could have been obtained if the surface stiffness parameter had been adjusted for each cross-section. This would also be a reasonable approach, as the surface stiffness represents the components ability to get a local depression at the contact point. This is clearly different for different materials and geometries. However, the results indicate that using the suggested stiffness is sufficient for most purposes.


## 5. Conclusions

In this paper, a 3-dimensional model for simulation of complex umbilical cross-sections were presented and tested with respect to published data for axisymmetric response. The model was based on the finite element approach, enabling contact interface effects to be handled on individual component level and in 3 dimensions. Material nonlinearities, gap formation, friction, lateral contact between wires, contact with external structures, and wires' curvature change are taken into account. The model can estimate the stresses and the displacements of individual structural element as well as the overall structural response. The influence of gap formation effects between structural elements was handled by a surface stiffness penalty parameter applied to the contact elements, and generally good correlation was found with test data. This demonstrates model robustness at least for the range of umbilical cross-sections used in the experiments. The model enables generality with respect to external loading, and work is ongoing with respect to investigating the model performance with respect to bending and friction stresses as well.



---

*Source: 985349-2012-08-02.xml*