Applied Ocean Research 88 (2019) 48–62

Contents lists available at ScienceDirect

Applied Ocean Research

journal homepage: www.elsevier.com/locate/apor

T
Feasible numerical method for analysis of offshore pipeline in installation
Pavel A. Trapper
Department of Structural Engineering, Faculty of Engineering Sciences, Ben-Gurion University of the Negev, P.O.Box 653, Beer-Sheva 84105, Israel

A R T I C LE I N FO A B S T R A C T

Keywords: Feasible numerical method for a structural analysis of a pipeline configuration during the installation process is
Finite difference method presented. The method considers the whole pipeline, which is partially suspended and partially laid-on a seabed,
Energy minimization as a single continuous segment, and is valid for a complete range of laying angles between 0°–90°, i.e., valid for
Nonlinear analysis both S-lay and J-lay configurations. The method accounts for a pipeline–seabed interaction and the pipeline is
Pipeline-seabed interaction
modeled by means of nonlinear large deformation beam theory. The numerical solution is carried out in an
Elastic seabed
Moving boundary conditions
incremental-iterative manner by following the actual pipeline installation process, and thus allowing efficient
treatment of pipeline-seabed interaction circumventing the further complexities with contact detection. At each
increment, the length of the pipeline is increased and new sequential equilibrium configuration is assessed by
direct minimization of a total potential energy approximated as a Riemann sum, which yields algebraic system of
nonlinear finite difference equations that is further solved by iterations with Newton-Raphson technique. The
simplicity, flexibility and robustness of the proposed method allow to enhance the efficiency of engineering
calculations and design. Accounting for a bending stiffness in a suspended part allows analyzing variations in
laying angle and lay tension independently. The method convergence is validated and compared with Abaqus.
The results are in an excellent agreement. Moreover, the comparison with Abaqus shows that for the selected
parameters the assumption that the pipeline is inextensible and unshearable is very reasonable. Representative
parametric study is conducted to demonstrate the feasibility of the method. Parametric study considers the
effects of laying angle (0°–90°), lay tension, laying water depth (up to 3000 m) and seabed stiffness.

1. Introduction should be performed. Such an analysis can be usually carried out by

using commercial software, such as OFFPIPE, OrcaFlex, Riflex, etc. In
A design of offshore pipelines is a laborious process. The purpose of more involved cases or for more refined analyses even heavier general
the design is to achieve an optimal pipeline sizes and material grades to purpose finite element packages such as Abaqus or Ansys, which are
satisfy all the possible loadings and configurations during the pipeline able to perform the full 3D Multiphysics analysis, may be required.
life-time. The most severe stage is usually the pipeline installation These codes can yield accurate results, but they are quite time-con-
(Fig. 1), where the pipeline is usually being subjected to the most severe suming, and therefore are usually appropriate to be applied, if required,
loading due to the lay effects [1,2]. Consequently, the pipeline con- only for refined analysis at the final stage of the design, where the
figuration during the installation should be properly addressed. detailed specification is defined [3].
The design is usually carried out in several stages starting from a Taking into consideration that the dominant first order effects,
general conceptual design, which outcome is a rough estimate of as- which define a general pipe-lay configuration are usually static and
sociated cost of the project; and down to a detailed design, which occur in 2D [3], offshore engineers tend to look for more simple and
outcome is the final detailed specification for the construction [1]. time-saving ways to perform the analyses, especially at the earlier
Since the design is an iterative process, during all its stages numerous stages of the design, where preliminary calculations are being per-
pipe-lay analyses, which aim at refining the configuration may be re- formed to define a general pipeline configuration and estimate a rough
quired. cost of the project. Moreover, simple methods can reveal the relevance
In general, pipeline installation is a highly nonlinear 3D process, of further more sophisticated analyses.
which involves coupled in-plane and out-of-plane phenomena, such as Over the past decades there was a series of such techniques pro-
full pipeline-seabed interaction, lateral buckling, lateral stability under posed, both analytical (e.g., [3–5]) and numerical (e.g. [6,7],). Those
hydrodynamic loads, torsion, sideway walking, etc. In order to capture techniques are usually specific to a particular laying method: S-lay or J-
all these phenomena accurately, a complete dynamic 3D analysis lay technique, which mainly differ by a laying angle (Fig. 1.). Moreover,

E-mail address: [email protected].

https://doi.org/10.1016/j.apor.2019.04.018

0141-1187/ © 2019 Elsevier Ltd. All rights reserved.

Received 25 January 2019; Received in revised form 4 March 2019; Accepted 17 April 2019
P.A. Trapper Applied Ocean Research 88 (2019) 48–62

Fig. 1. Pipeline installation techniques: (a) S-lay; (b) J-lay.

those techniques usually divide the pipeline into number of segments We also leave out hydrodynamic forces due to underwater currents.
assuming different idealized structural behavior in each segment. Sus- And finally, in this paper we consider only horizontal flat seabed, still
pended part is usually treated as a cable with zero bending stiffness, and water and constant temperature.
laid-on part, is treated like a beam, which undergoes small deforma- The paper is organized as follows. In the following section we pre-
tions on a Winkler foundation (e.g., [6]). In some models the number of sent the proposed methodology, discretization and the solution tech-
such segments may become even higher (e.g. [5,7],). nique. Then we validate the model with Abaqus general purpose finite
Beside reduced accuracy due to various simplifications and re- element package and check numerical convergence. In what follows,
stricting assumptions, those models are very case-specific and will not representative parameter analyses are conducted to demonstrate the
simply allow adaptation of further modifications, such as additional feasibility of the method. In particular, the effects of laying angle, top
loadings or geometric constraints. In addition, ignoring the bending lay tension, laying water depth and seabed stiffness are considered.
stiffness in a suspended part disables the option to analyze laying angle
and lay tension variations independently, which significantly limits the
2. Model overview
application range of those methods since local perturbations of a laying
angle (e.g., due to waves) can lead to significant rise in internal forces
The structural analysis of the pipeline during the installation aims at
(Section 5.2.1).
determination of the overall configuration, including the horizontal
The aim of this paper is to find the balance between the simplicity
distance of the touchdown point (TDP) from the lift-off point (LOP) and
and the robustness of the calculations, thus enhancing the efficiency of
the embedment into the seabed (Fig. 1), and determination of the in-
engineering design. In this paper we present a simple numerical method
ternal forces (bending moments, shear forces and tension).
for the structural analysis of a pipeline configuration during the in-
While the internal forces and the embedment are more severe at the
stallation process. The method considers the whole pipeline, which is
neighborhood of the TDP, at some point away from the TDP, in opposite
partially suspended and partially laid-on a seabed, as a single con-
direction to the laying direction, the laying effects fade away and only
tinuous segment, and is valid for a complete range of laying angles
the embedment due to submerged self-weight of the pipeline remains
between 0°–90°, i.e., valid for both S-lay and J-lay configurations. The
(see Fig. 2). Let denote this point as SSP (steady state point). Beyond
method accounts for a pipeline–seabed interaction and the pipeline is
that point the forces and the embedment do not vary with the distance
modeled by means of nonlinear large deformation beam theory.
any longer and the remaining as-laid part of the pipeline, LAL , is not
The numerical solution is carried out in an incremental-iterative
relevant for the analysis, provided proper boundary conditions are
manner by following the actual pipeline installation process, and thus
specified at SSP.
allowing efficient treatment of pipeline-seabed interaction cir-
Two challenges do emerge. The first one is that the location of the
cumventing the further complexities with a contact detection. At each
above point SSP relative to LOP, which is consistent with pipeline
increment, the length of the pipeline is increased and new sequential
segment of length LAN (Fig. 2) that should be actually analyzed, is not
equilibrium configuration is assessed by direct minimization of a total
known in advance. The second challenge is that even for a given LAN ,
potential energy discretized as a Riemann sum, which yields algebraic
the location of TDP is still not known in advance either, i.e., it is still not
system of nonlinear finite difference equations that further solved with
known in advance which portion of LAN is already in a contact with a
Newton-Raphson iterative technique.
seabed, LTD , and which portion of it is still suspended, Lsus (Fig. 2).
The relative flexibility and simplicity of the model, allows easily
The first challenge is likely to be settled if long enough pipeline is
conducting sensitivity analyses and re-analyses in order to define the
assumed for the analysis. Otherwise, iterative techniques may be ap-
key roles of various parameters involved and assisting the offshore
propriate. On the other hand, the second challenge, in general, can be
engineer to decide which pipe-lay option is preferable in the certain
treated by techniques proposed by the theory of constrained optimi-
conditions. Accounting for a bending stiffness in a suspended part al-
zation, which apply special algorithms to identify the active set of
lows analyzing variations in laying angle and lay tension independently
material points to which the contact condition in the form of equality
of each other.
constraints should be applied [9].
The proposed model adopts the following assumptions and limita-
In reality, pipeline installation is a sequential process, where pipe-
tions. The model is limited to a static planar 2D problem. 3D out of
line material points are sequentially come into a contact with a seabed,
plane effects, such as torsion, lateral buckling, walking and etc. are not
i.e., it is known in advance which next material point will enter the
considered. The pipeline is modeled as inextensible unshearable beam,
active set of equality constraints. Setting numerical solution to follow
such that only the deformation due to bending is considered.
this reality allows us to circumvent both the abovementioned chal-
In addition, in order to demonstrate direct energy minimization
lenges in a simple way.
technique, which is valid in general only for conservative and path-
Consequently, in this paper the solution is being carried out as
independent systems, we make additional assumptions as follows. We
follows (Fig. 3). First, we set some initial configuration of a pipeline of
consider only an elastic pipe: stinger effects, which may involve plastic
some initial length L0 , such that it will have no more than one point in a
deformations at the top of suspended part [8], are ignored. We consider
contact with a seabed. The most convenient way to do this is probably
only linear elastic frictionless seabed in a form of a Winkler foundation.
to provide a catenary shape (e.g., [10]) taking into account given water

49
P.A. Trapper Applied Ocean Research 88 (2019) 48–62

Fig. 2. Pipeline configuration in installation: suspended and touchdown segments.

Fig. 3. Pipe-lay analysis technique: advance of numerical solution.

depth, lay tension and boundary condition at the fixed point FP (green
dashed line in Fig. 3). The above catenary shape ignores bending
stiffness of the pipeline and does not satisfy the fixed laying angle at Fig. 4. Scheme of the pipeline: geometry, loads, boundary conditions and co-
ordinate systems.
LOP. This initial configuration serves as initial estimate for numerical
solution.
From then and on, the solution is carried out in an incremental- a distributed submerged self-weight w and a horizontal tension point
iterative way, following the analogy to the real pipe-lay process, where load P0 controlled by the laying vessel at LOP (aka lay tension). The lay
the pipeline, with a fixed end FP on the bottom (only in x-direction), tension P0 is actually a horizontal component of axial tension force
being gradually released from the laying vessel into the water (Fig. 3.). within the pipe at LOP.
During this process, pipe material points are laid down and, one by one, The model is established in the global coordinate system xy where
they come into contact with the seabed soil. Technically, for each the origin is set at top left corner of the domain (Fig. 4) with x-axis
loading step j , we increase the length of the pipe by ∆L , such that directed to the right in the direction of the laying vessel propagation
Lj = Lj − 1 + ∆L , and solve static equilibrium to obtain a new config- and y-axis directed downwards to the seabed. The pipeline configura-
uration. For each consequent length increment as initial estimate serves tion is given by θ (s ) that defines the local orientation of the pipeline,
the configuration from the previous length increment (Fig. 3). Within where s is the curvilinear distance along the pipeline (0 ≤ s ≤ L )
each increment, the configuration is being iteratively updated by measured from the fixed point FP.
Newton-Raphson technique. Active set [11] of contact points beyond The geometric boundary conditions for the above configuration are
the TDP is sequentially updated with a new point. Namely, at each as follows. At point FP the pipeline is restrained in x-direction and
iteration, the contact force is being applied to material point if it falls against the rotation, while is free to move in y-direction:
beyond the water depth, i.e., y ≥ h w . The process should quit when the
x (0) = 0 (1)
steady state configuration is achieved where the effect of laid end (point
FP) boundary conditions is negligible, i.e., the internal forces and the θ (0) = 2π (2)
embedment at any point on the pipeline segment between LOP and SSP
do not change any longer with the increments. It should be noted that On the other hand, at LOP the pipeline is rotated with a prescribed
this problem involves moving boundary conditions and with each angle, φ0 , and restrained in y-direction, while is free to move in x-di-
length increment the size of the problem grows. rection:
In what follows we present mathematical model and numerical y (L) = 0 (3)

θ (L) = 2π − φ0
technique to solve static equilibrium at each length increment.
(4)

3. Mathematical formulation The reason for the choice of the coordinate system as defined in
Fig. 4 is that this way boundary conditions (1) and (3) can be
Let us consider an elastic pipeline of length L , constrained to deform straightforward adapted in the formulation in the easiest manner. Then,
in plane, partially laid on a seabed with a linear stiffness (per unit the position of the pipeline is given by Cartesian coordinates x and y as
length) ks , and partially suspended (Fig. 4). The pipeline is subjected to follows

50
P.A. Trapper Applied Ocean Research 88 (2019) 48–62

x (s ) = ∫0
s
cosθ (s ') ds' (5)
The second term in Eq. (7) represents the elastic energy stored in the

∫s
deformation of elastic seabed. Since there is only a fraction of the pi-

y (s ) = −
L peline that is in contact with a seabed, while the other one is being
sinθ (s ') ds' (6) suspended, the seabed stiffness function kSB (s ) and the embedment, e
Pay attention, Eq. (5) and (6) satisfy boundary conditions (1) and (s), are given piecewise by

ks, y (s ) ≥ h w
kSB (s ) = kSB (y (s )) = ⎧
(3) automatically. Establishing model in other coordinate system rather
than the one as defined in Fig. 4 may require more sophisticated way of ⎨
⎩ 0, y (s ) < h w (9)
boundary conditions (1) and (3)specification (by means of the Lagrange
multipliers [12], for example), which might significantly complicate the and

y (s ) − h w , y (s ) ≥ h w
e (s ) = e (y (s )) = ⎧
solution.

⎨
⎩ 0, y (s ) < h w (10)
3.1. Minimum total potential energy principle
respectively. We assume that there is no suction force in a case of a
The more common way to derive the field equations in engineering separation [14]. The seabed resistance (per unit length) is given by
Rs (s ) = kSB (s ) e (s )
is by considering balance of forces and moments acting on the differ-
ential part of the body (e.g., [3–7]). (11)
In this paper, we derive equations, which define the pipeline con- The third term in Eq. (7) represents the potential energy of sub-
figuration at equilibrium directly from Minimum total potential energy merged self-weight of the pipe, w . And the last term in Eq. (7) represent
principle. This principle says that among all the possible geometric the potential energy of lay tension P0 .
configurations consistent with the geometric constraints and forces Once θ (s ) is defined, it is possible to evaluate bending moments,
applied on the system, only the one that minimizes the total potential shear forces and tension using the usual differential relations given here
energy corresponds to the stable equilibrium configuration (e.g., [13]). in the above defined coordinate system
M (s ) = EIκ (s )
For static cases, the total potential energy consists of elastic strain en-
ergy of deformed body and potential energy of applied forces. The (12)

S (s ) = −
above approach is usually simpler than more widely accepted equili- dM (s )
brium considerations, since it naturally accounts for applied loads and (13)
ds
force boundary conditions, and often leads towards more convenient
= − w sinθ (s )
mathematical expressions, such as lower order equations, sometimes dT (s )
even with a lower number of variables [12]. Nevertheless, both ds (14)
methods are equivalent, and eventually yield the same results. By ap-
plying calculus of variations [12] it is possible to derive those equili-
4. Numerical implementation
brium field equations directly from total potential energy. Nevertheless,
in the present paper we do not apply the calculus of variations, but
4.1. Finite difference discretization
directly minimize the discretized total energy expression.
The limitation of the direct application of the above principle is that
The simplest way to tackle the above problem is by finite difference
it is valid only for conservative systems, where after the loading re-
method (nevertheless, it can be equally solved by finite elements).
moval both the pipeline and the seabed return to their initial un-
The total potential energy of the pipeline, Eq. (7), can be expressed
disturbed state, i.e., plastic seabed deformation, plastic pipeline de-
in terms of discrete variables (Fig. 5) written as a left Riemann sum:
formation and hydrodynamic loads (unless assumed to work similar to
1 ⎛ θi + 1 − θi ⎞2
N +1 N +1
Utot = ∆s ∑ + ∑ k SB (yi − h w )2 − ∑ wi yi − P0 xN +1
the body forces) are out of the scope. Consequently, in order to de- N
1 i
2 ⎝ ∆s ⎠
monstrate pure direct energy minimization approach, in this paper, we EI
i=1 i=1 i=1
2
limit ourselves to consider elastic pipeline and elastic seabed, leaving
out hydrodynamic forces and other non-conservative, path-dependent (15)
processes. Nevertheless, this method, in its present form, still yields where i and j are the indices (in what follows both indices are used
proper results provided the assumptions are appropriate for the ana- interchangeably). There are N subintervals of size ∆s = L/ N , and N + 1
lyzed case. Moreover, in general, the proposed technique, though not in nodes; wi = w∆s is a node gravity load and kSB
i
is a node seabed stiffness,
its direct form, can be extended to non-conservative cases. obtained by rewriting Eq. (9) as follows

ks ∆s, yi ≥ h w
=⎧
⎩ 0, yi < h w
⎨
3.2. Total potential energy i
kSB
(16)
The total potential energy of the system shown in Fig. 4 is given by

Utot =
1
∫0 EIκ (s )2ds +
1
∫0 kSB (s ) e (s )2ds − ∫0 wy (s ) ds − P0 x (L)
L L L

2 2
(7)
where E is the Young’s modulus of the pipe, I is the moment of inertia
of the pipe’s cross-section, κ is its planar curvature

κ (s ) =
dθ (s )
ds (8)
kSB (s ) is the seabed stiffness function and e (s ) is the embedment of the
pipeline into the seabed.
The first term in Eq. (7) represents the bending elastic energy stored
in a deformed pipeline. Since in this paper we consider only planar
deformation of inextensible unshearable pipeline, the strain energy is
contributed exclusively by bending. Fig. 5. Numerical scheme of the pipeline: finite difference discretization.

51
P.A. Trapper Applied Ocean Research 88 (2019) 48–62

At each step, the pipeline configuration is given by rewriting Eq. (5)

and (6) in a discrete form as follows
i−1
xi = ∑ cosθj ∆s
j=1 (17)

yi = − ∑ sinθj ∆s
N

j=i (18)
In order to achieve a proper pipeline configuration at static equili-
brium we look for extremum (which yields minimum) of total energy as
a function of N − 1 discrete variables θ2, …, θN :
∂Utot (θ1, …, θN + 1)
= 0 for i = 2, …, N
∂θi (19)
The conditions of stationarity with respect to θ1 and θN + 1 are not
required since these values are known from boundary conditions (2) Fig. 6. Scheme of the pipeline: depiction of the presented parameters and co-
and (4) and can be directly substituted ordinate systems.
θ1 = 2π (20)

θN + 1 = 2π − φ0
Table 1
(21) Basic parameters.
Substituting Eq. (15),(17) and (18) into Eq. (19) yields Parameters Symbol Value

⎡ ⎤
(−θi − 1 + 2θi − θi + 1) + ∑ − kSBj ⎢ ⎛⎜− ∑ ∆ssinθk ⎞⎟ − h w ⎥ ∆scosθi
i N
EI Pipe’s outer diameter [m] D 1.22

⎢
⎣⎝ ⎠ ⎥
⎦
Pipe’s wall thickness [m] t 0.0318
∆s j=1 k=j Elasticity modulus [N/m2] E 2.1·1011

∑ − wj ∆scosθi − P0 (−∆ssinθi) = 0
Steel density [kg/m3]
−
i ρs 7850
Seawater density [kg/m3] ρw 1030
j=1 Acceleration of gravity [m/s2] g 9.81
Laying angle [°] φ0 0°
Lay tension [N] P0 3·106
for i = 2,…, N (22) Laying water depth [m] hw 200
Seabed stiffness per unit length of the pipe [N/m2] ks 4·103
which is a set of difference equations in terms of θi . Finite difference subinterval size [m] ∆s 4
In order to obtain a more convenient expression for (22) it is rea-

Kθ (m) + F (θ (m) ) = 0
sonable to apply the Fubini Principle for sequences for exchange of the
order of summation in a second term (27)

⎛ ⎞
∑ kSBj ⎜∑ ∆ssinθk ⎟ = ∑ ∆ssinθj ⎛⎜ ∑ k ⎞
i N N min(i, j )
We solve Eq. (27) iteratively with Newton-Raphson method [9].
kSB ⎟
⎝ k=j ⎠ ⎝ ⎠
Linearizing Eq. (27) with a first order term of Taylor series expansion
j=1 j=1 k=1 (23)
yields the following linear equation for iteration n
and introduce the aggregated values as follow
⎛⎜K + dF ⎞⎟ ∆θ (m, n) = −Kθ (m, n − 1) − F (θ (m, n − 1))
= ∑ kSBk ⎝ ⎠
i
dθ θ (m, n − 1)
i
KSB (28)
k=1 (24)

∑ wk
with
Wi =
i

k=1 (25) θ (m, n) = θ (m, n − 1) + ∆θ (m, n) (29)

Substituting Eq. (23),(24) and (25) into Eq. (22) yields which is solved repeatedly. The iterative process quits when ∆θ (m, n)

∑ KSBmin(i,j) ∆ssinθj + h w KSBi ∆scosθi

meets the convergence criterion. The solution of previous increment,
(−θi − 1 + 2θi − θi + 1) + ∆scosθi θ (m − 1) , serves as initial estimate for subsequent increment,
N
EI
∆s j=1 θ (m,0) = θ (m − 1) .

+ Wi ∆scosθi + P0 ∆ssinθi = 0
The expressions, which appear in Eq. (28) are as follows
(26)
The Eqs. (20–21,26) is a nonlinear algebraic system of N − 1 ⋱ ⋱ ⋱
EI ⎡ ⎤
⎢ −1 2 −1 ⎥
equations with N − 1 unknowns θi , i = 2, …, N . K=
⎢ ⋱ ⋱ ⋱⎥
⎣ ⎦
∆s
(30)
4.2. Solution technique
and the other nonlinear terms are presented here in a form of matrix
The solution of nonlinear system of finite difference Eqs. (20–21,26) entries,
is carried out in an incremental - iterative manner as follows. At each
increment m , we increase the pipeline length by ∆L (=∆s) , i.e.,
∑ KSB
(m, n − 1)
{F (θ (m, n − 1) )}i = ∆scosθi(m, n − 1) ∆ssinθj(m, n − 1)
N

L(m) = L(m − 1) + ∆L . Consequently, at each increment, the dimensions

min(i, j )

of the whole system grow by one: N (m) = N (m − 1) + 1, since we add j=1

additional node.
Now, for the sake of clarity, let us rewrite Eq. (20–21, 26) for in- + h w KSB
i (m, n − 1) ∆scosθ (m, n − 1) + W ∆scosθ (m, n − 1) + P ∆ssinθ (m, n − 1)
i i i 0 i (31)
crement m in a matrix form as follows

52
P.A. Trapper Applied Ocean Research 88 (2019) 48–62

Fig. 7. Validation and comparison: (a) Bending moment vs. distance from the laying vessel for different mesh densities and dense Abaqus solution (dashed line); (b)
Maximum bending moment vs. varying number of subintervals/elements; (c) Shear force vs. distance from the laying vessel for different mesh densities and dense
Abaqus solution (dashed line); (d) Maximum shear force vs. varying number of subintervals/elements; (e) Tension force vs. distance from the laying vessel for
different mesh densities and dense Abaqus solution (dashed line); (f) Maximum tension force vs. varying number of subintervals/elements.

Table 2
The errors for the maximum values of bending moments, shear forces and tension.
Subinterval size, ∆s [m] Number of subintervals, N Maximum bending moment error [%] Maximum shear force error [%] Maximum tension force error [%]

8 104 11.0 18.1 8.4

4 209 5.2 8.8 4.3
2 418 2.2 4.0 2.1
1 837 0.6 1.6 1.0
0.5 1674 0.1 0.4 0.5

⎡ dF ⎤ = ∆scosθ (m, n − 1) K min(i, j) (m, n − 1) ∆scosθ (m, n − 1) and tension can be back-calculated by using Eq. (14–16) in their dis-
⎢ dθ ⎥
⎣ (m, n − 1)
⎦ij
i SB j crete forms as follows
θ

θi − θi − 1
+ δij ⎡ ⎤ Mi = EIκi = EI
*
dF
⎢ dθ ⎥
⎣ θ (m, n − 1) ⎦ii
∆s (35)
(32)
Mi − Mi − 1
Si = −
⎡ dF ⎤ = −∆ssinθ (m, n − 1) ∑ KSB
(m, n − 1)
*
∆ssinθj(m, n − 1)
N

⎢ dθ ⎥
min(i, j ) ∆s (36)
⎣ θ (m, n − 1) ⎦ii
i
j=1

⎡ ⎤
− h w KSB
i (m, n − 1) ∆ssinθ (m, n − 1) − W ∆ssinθ (m, n − 1) + P ∆scosθ (m, n − 1) Ti = P0cosθi − ⎢RLOP + ∑
N
(33) wj⎥ sinθi
⎣ ⎦
i i i 0 i
j=i+1 (37)
where δij is a Kronecker delta

0 if i ≠ j
δij = ⎧
where the vertical reaction at LOP (Fig. 6) is given by

⎩ 1 if i = j
⎨ (34) RLOP = − TN + 1sinθN + 1 + SN + 1cosθN + 1 (38)

Once the variables θi are calculated, bending moments, shear forces and tension at LOP by

53
P.A. Trapper Applied Ocean Research 88 (2019) 48–62

Fig. 8. (a) Configurations for different laying angles; (b) Distance from the laying vessel/LOP to TDP vs. varying laying angle.

Fig. 9. (a) Pipeline embedment into the seabed for different laying angles; (b) Maximum pipeline embedment vs. varying laying angle; (c) Seabed resistance for
different laying angles; (d) Maximum seabed resistance vs. varying laying angle; (e) Vertical reaction at LOP vs. varying laying angle; (f) Bending moment at LOP vs.
varying laying angle.

54
P.A. Trapper Applied Ocean Research 88 (2019) 48–62

Fig. 10. (a) Bending moment diagrams for different laying angles; (b) Maximum absolute bending moment vs. varying laying angle; (c) Shear force diagrams for
different laying angles; (d) Maximum absolute shear force vs. varying laying angle; (e) Axial tension force diagrams for different laying angles; (f) Maximum axial
tension force vs. varying laying angle.

P0 − SN + 1sinθN + 1 resistance at as-laid part is Rs, ∞ = w .

TN + 1 =
cosθN + 1 (39) Though we formulate the problem in the global coordinate system
xy, for the sake of convenience we adopt two additional abscissas to
present the results. The first one, xLOP , measured from LOP; and the
5. Validation and parameter analysis second one, xTDP , measured from TDP, both in negative global x di-
rections (Fig. 6).
In what follows we present the validation of the proposed method
with Abaqus [15], finite difference grid convergence and parameter
5.1. Validation and comparison
analysis. For all the forthcoming calculations we adopt the basic
parameters as summarized in Table 1. All the calculations, except for
The present method is first validated through the basic case
Abaqus reference solution, were performed with MATLAB.
adopting the characteristics as presented in Table 1 with finite element
We adopt the pipeline dimensions and materials following Ref. [3],
solution by Abaqus/Standard [15]. For the sake of comparison, a pi-
while other parameters were chosen in a range as proposed by Ref.
peline with a total length of 836 m was selected. The effect of finite
[1,2,14].
difference grid density on bending moments, shear forces and tension is
In the following simulations we do not apply the above-mentioned
analyzed. In order to demonstrate the reliability of the validation we
criterion to stop the analysis when the steady state is reached, since we
demonstrate the convergence of Abaqus solution either. Since Abaqus
would like to visualize the effect of boundary conditions on the overall
formulation adapts extensibility and shear-flexibility we can also test
configuration.
the authenticity of the assumptions.
Other parameters, which appear in Eq. (26) are calculated as fol-
lows. The moment of inertia of the pipe’s cross-section is given by
I = π /64[D 4 − (D − t ) 4], while its area by A = π /4[D 2 − (D − t )2]. A 5.1.1. Abaqus model
submerged (buoyant) density of the pipe is given by ρb = ρs − ρw . For a more reasonable comparison, in Abaqus model, two-node
Considering the pipe being filled with water during the installation, the hybrid linear beam elements (B21H) were selected to model the pipe-
submerged weight per unit length of the pipeline is given by w = ρb gA. line [15]. The hybrid elements improve the performance by including
In what follows, the results are presented for the quantities as illu- the axial and transverse shear forces as primary variables along with the
strated in Fig. 6, where emax is a maximum embedment, e∞ = w / ks is an nodal displacements and rotations. The seabed is modelled as 2D planar
as-laid embedment, and Rs, max is a maximum seabed resistance. Seabed analytical rigid shell body. The contact between the pipeline and the

55
P.A. Trapper Applied Ocean Research 88 (2019) 48–62

Fig. 11. (a) Configurations for different lay tension; (b) Distance
from the laying vessel/LOP to TDP vs. varying lay tension; (c)
Pipeline embedment into the seabed for different lay tension; (f)
Maximum pipeline embedment vs. varying lay tension; (e) Seabed
resistance for different lay tension; (f) Maximum seabed resistance
vs. varying lay tension.

Fig. 12. (a) Bending moment diagrams for different lay tension; (b) Maximum absolute bending moment vs. varying lay tension; (c) Shear force diagrams for
different lay tension; (d) Maximum absolute shear force vs. varying lay tension.

seabed is defined using softened contact pressure-overclosure re- 5.1.2. Comparison

lationship [15]. Comparisons of bending moments, shear forces and tension are
shown in Fig. 7a,c,e respectively. From these figures, it can be seen that
the results obtained by the proposed method and Abaqus are in a good

56
P.A. Trapper Applied Ocean Research 88 (2019) 48–62

Fig. 13. (a) Axial tension force diagrams for different lay tension; (f) Maximum axial tension force vs. varying lay tension; (c) Normalized axial tension force
diagrams for different lay tension; (d) Normalized maximum axial tension force vs. varying lay tension.

Fig. 14. (a) Configurations for different laying water depths; (b) Distance from the laying vessel/LOP to TDP vs. varying laying water depth; (c) Normalized distance
from the laying vessel/LOP to TDP vs. varying laying water depth. (d) Pipeline embedment into the seabed for different laying water depths; (e) Seabed resistance for
different laying water depths.

57
P.A. Trapper Applied Ocean Research 88 (2019) 48–62

Fig. 15. (a) Bending moment diagrams for different laying water depths ( xLOP ≥ 100 ); (b) Bending moment diagrams for different laying water depths
(100 ≥ xLOP ≥ 0 ); (c) Maximum absolute bending moment vs. varying laying water depth.

agreement. Fig. 7b,d,f show the convergence of those quantities vs. φ0 . Horizontal distance from the laying vessel/LOP to the TDP,
number of subintervals for the proposed method, and vs. number of LLOP − TDP , decreases almost linearly with growing laying angle (as the
elements for Abaqus results with hybrid elements. It can be seen that angle becomes steeper) as shown in Fig. 8b. Both pipeline embedment
while the bending moment in both models converges nearly at the same and seabed resistance look identical for all the laying angles (Fig. 9a–d).
rate, shear force and tension within Abaqus framework, due to hybrid As shown in Fig. 9e, vertical reaction between the pipeline and the
formulation, have close values already with a relatively coarse mesh, laying vessel at LOP does not change much for lower laying angles
while within the proposed method both shear and tension are back- while grows drastically towards 90°. On the other hand, the bending
calculated from bending moments. The maximum error of the solution moment at LOP (bending moment the pipeline imposes on the vessel)
for various meshes relatively to Abaqus reference solution (h el = 0.5 m ) increases nearly linearly with the angle throughout the whole angles
is shown in Table 2. It can be seen that the highest error for all the range (Fig. 9f). The curve passes through zero bending moment at
discussed quantities is at their maximum values. At other locations the 47.07°. Alternatively, this configuration could be received directly,
error is lower. The maximum errors for bending moments and shear simply by not prescribing boundary condition on φ0 at all, which con-
forces occur at LOP, while for tension about 50 m away from LOP to- sistent with prescribing at LOP zero bending moment (pinned condi-
wards TDP (Fig. 6). In this example, the solution with dense grid tion).
( ∆s = 0.5 m ) shows error less than 0.5% for all those quantities. The Figs. 10a,c,e depict bending moments, shear forces and axial tension
errors for maximum values are summarized in Table 2. for different laying angles. Laying angle highly affects both maximum
absolute bending moments and maximum absolute shear forces at LOP,
5.2. Parameter analysis while for the chosen parameters (Table 1) has no visible effect on those
values away from LOP at sagbend (Fig. 10b and d). It should be noted,
In what follows we present the parameter analysis. The basic set of that choice of other parameters rather than those given in Table 1 might
parameters is as summarized in Table 1. Each time we change values of lead to other conclusions. Maximum tension grows only slightly for
only one specific parameter. lower laying angles while increase growing towards 90° analogically to
The parametric study presented in this paper is by no means com- the vertical reaction (Fig. 10f). In reality, tension at laid-on part away
prehensive. The aim of this study is to demonstrate the general scope of beyond the TDP decreases due to friction forces applied by the seabed,
the proposed method and to show the general tendencies, rather than however, maximum values are still at suspended part.
provide specific design values. Consequently, the conclusions raised
about the results in the following subsections are related only to the 5.2.2. Effect of lay tension
used parameter ranges. In what follows we study the effect of lay tension (independently of
a laying angle). In this series we perform a set of five simulations for
5.2.1. Effect of laying angle different lay tension forces: 0.1 MN, 0.5 MN, 1 MN, 2 MN and 4 MN.
Sharper laying angles consistent with J-lay technique while more Other parameters adopted for this series kept constant as in Table 1.
gentle with S-lay one. The difference, from the mechanics standpoint is Fig. 11a shows the plot of configurations for different lay tension
only the angle. Since one of the distinctive features of the proposed force, P0 . Horizontal distance from the laying vessel/LOP to the TDP,
model is its ability to solve for any top laying angle varying from 0° to LLOP − TDP , increases with lay tension as shown in Fig. 11b. Both max-
90° independently of lay tension, we shall present this series of simu- imum pipeline embedment and maximum seabed resistance drop with
lations first. Other parameters adopted for this series kept constant as in increasing lay tension, approaching pipeline as-laid embedment and
Table 1. pipeline self-weight values respectively (Fig. 11c–f).
Fig. 8a shows the plot of configurations for different laying angle, Fig. 12a,c depict bending moments and shear forces for different lay

58
P.A. Trapper Applied Ocean Research 88 (2019) 48–62

Fig. 16. (a) Shear force diagrams for different laying water depths ( xLOP ≥ 100 ); (b) Shear force diagrams for different laying water depths (100 ≥ xLOP ≥ 0 ); (c)
Maximum absolute shear force vs. varying laying water depth; (e) Axial tension force diagrams for different laying water depths; (f) Maximum axial tension force vs.
varying laying water depth.

tension. It can be seen that maximum absolute bending moment at LOP reasonable.
only very slightly drops with the increase of lay tension, while at the
sagbend the drop is more substantial (Fig. 12b). On the other hand,
5.2.3. Effect of laying water depth
maximum absolute shear force at LOP rises with the lay tension almost
In this series we perform a set of six simulations for different laying
proportionally, while at the sagbend the drop is less substantial
water depths: 100 m, 200 m, 500 m, 1000 m, 2000 m, and 3000 m.
(Fig. 12d).
Other parameters adopted for this series kept constant as appear in
Fig. 13a depicts axial tension for different lay tension. Maximum
Table 1.
tension not only rises linearly with lay tension but also with the same
Fig. 14a shows the plot of configurations for different laying depths,
rate, i.e., the difference between the lay tension value and maximum
h w . Horizontal distance from the laying vessel/LOP to the TDP,
tension value is constant for all the tested lay tension values (Fig. 13b).
LLOP − TDP , increases with laying depth as shown in Fig. 14b, however the
However, in order to study this effect more explicitly, it is useful to look
above distance relative to laying depth exponentially drops (Fig. 14c).
at the normalized tension values T / P0 and Tmax / P0 (Fig. 13c and d).
Both pipeline embedment and seabed resistance look identical for all
Fig. 13c,d show that though the maximum tension rises with rising lay
the laying depths (Fig. 14d–e).
tension, maximum normalized tension drops, i.e., tension in all material
In order to better visualize bending moments for different laying
points along the pipeline approach the lay tension. Again, as friction
depths at both the sagbend and at the LOP, we split the bending mo-
ment diagram in two intervals, xLOP ≥ 100 m (Fig. 15a) and
forces are not considered, the residual tension at the laid-on segment
100 m≥ xLOP ≥ 0 m (Fig. 15b) respectively. It can be seen that max-
equals the lay tension. For higher lay tension values the assumption of
constant tension distribution along the pipeline [4,5] becomes more
imum absolute bending moment at LOP substantially grows with

59
P.A. Trapper Applied Ocean Research 88 (2019) 48–62

Fig. 17. (a) Configurations for different seabed stiffness; (b) Distance from the laying vessel/LOP to TDP vs. varying seabed stiffness; (c) Pipeline embedment into the
seabed for different seabed stiffness; (d) Maximum pipeline embedment vs. varying seabed stiffness; (c) Normalized pipeline embedment into the seabed for different
seabed stiffness; (e) Normalized maximum pipeline embedment vs. varying seabed stiffness.

increasing depth, while at the sagbend the maximum absolute bending Fig. 17a shows the plot of configurations for different seabed stiff-
moment is constant (Fig. 15c). ness, ks . Horizontal distance from the laying vessel/LOP to the TDP,
As with bending moments, for improved visualization, we split the LLOP − TDP , slightly grows with growing stiffness as shown in Fig. 17b.
shear forces diagram in two intervals, xLOP ≥ 100 m (Fig. 16a) and Fig. 17c depicts the plot of embedment for different seabed stiffness, ks .
100 m≥ xLOP ≥ 0 m (Fig. 16b) respectively. It can be seen that max- Maximum embedment, emax , exponentially drops with growing seabed
imum absolute shear force at LOP nearly linearly grows with increasing stiffness, nearly in a similar manner to as-laid (self-weight) embedment
depth, while at the sagbend the maximum absolute shear force is e∞ (Fig. 17d). In order to study the effect of seabed stiffness more ex-
constant (Fig. 16c). Fig. 16d depicts axial tension for different laying plicitly, it is useful to look at the normalized embedment values e / e∞
depths. Fig. 16e shows, that maximum tension grows linearly similar to and emax / e∞ (Figs. 17c,d), where the later relation is called touchdown
maximum absolute shear force at LOP. factor (TDF), which multiplies submerged weight of the pipe, and ty-
High growth in both maximum absolute shear force and maximum pically assumed to be 2–3 for the design [1]. Fig. 17c and d show that
tension with increasing laying depth is attributed to the increasing total though the maximum embedment drops with growing seabed stiffness,
weight of growing suspended part. maximum normalized embedment, which is TDF, substantially grows.
Fig. 18a depicts the plot of seabed resistance for different seabed
stiffness, ks . It can be seen, that the growth in seabed stiffness leads to
5.2.4. Effect of seabed stiffness growth in maximum seabed resistance, while as softer the seabed, the
In this series we perform a set of five simulations for different effect of sagbend becomes negligible and the resistance along the whole
seabed stiffness: 2 kN/m2, 4 kN/m2, 8 kN/m2, 16 kN/m2 and 32 kN/m2. laid part approaches the self-weight (Fig. 18b). Fig. 18c,d,e depict
Other parameters adopted for this series kept constant as in Table 1.

60
P.A. Trapper Applied Ocean Research 88 (2019) 48–62

Fig. 18. (a) Seabed resistance for varying seabed stiffness; (b) Maximum seabed resistance vs. varying seabed stiffness; (c) Bending moment diagrams for different
seabed stiffness; (d) Shear force diagrams for different seabed stiffness; (e) Axial tension force diagrams for different seabed stiffness.

bending moments, shear forces and axial tension for different seabed Riemann sum.
stiffness. It can be seen that all the curves are almost identical, which Parametric study demonstrated the ability of the method to provide
indicates that, in general, the internal forces within the pipeline do not the results for any laying angle in a range between 0°–90° in-
depend on a seabed stiffness. However, seabed stiffness highly affects dependently of lay tension, and this is due to accounting for a bending
pipeline embedment. stiffness in a suspended part. It was shown that variation in the laying
It should be noted, that in reality, accurately assessing pipeline angle can lead to significant increase in internal forces in the pipeline
embedment, is very important, since pipeline embedment into the close to LOP. Consequently, ignoring the bending stiffness in a sus-
seabed affects pipeline stability against hydrodynamic forces, laying pended part while keeping the laying angle fixed can bring to undesired
vessel movements sideways and lateral buckling, lateral and axial soil consequences due to vessel oscillations. On the other hand, away from
resistance and thermal insulation [14,16–21]. the LOP, the change in a laying angle has less substantial effect on the
internal forces. Nevertheless, it should be noted, that some works,
6. Conclusions which ignore the bending stiffness in a suspended part, reported quite a
wide variation in internal forces along the whole analyzed segment for
Simple numerical method for the structural analysis of a pipeline the variation of only 4° in a laying angle (e.g., [22]). The reason for this
configuration during the installation process is presented. The method is that in their works the lay tension force was automatically adapted so
considers the whole pipeline, which is partially suspended and partially that the condition of zero bending moment in the suspended segment
laid-on a seabed, as a single continuous segment, and is valid for both S- would be satisfied.
lay and J-lay configurations. The solution technique is based on a direct In the presented method, bending moment at LOP can be prescribed
minimization of a total potential energy of the system discretized as a by adding its potential energy to the total potential energy expression,

61
P.A. Trapper Applied Ocean Research 88 (2019) 48–62

while ignoring the specification of laying angle, φ0 (Eq. (4)). Conse- [6] Y. Zan, C. Yang, D. Han, L. Yuan, Z. Li, A numerical model for pipelaying on
quently, zero bending moment at LOP can be automatically prescribed nonlinear soil stiffness seabed, J. Hydrodyn. Ser. B. 28 (2016) 10–22, https://doi.
org/10.1016/S1001-6058(16)60603-0.
by just ignoring the Eq. (4). [7] L.-Z. Wang, F. Yuan, Z. Guo, Numerical analysis for pipeline installation by S-Lay
In addition parameter analysis has revealed that seabed stiffness method, ASME. Int. Conf. Offshore Mech. Arct. Eng. 29th Int. Conf. Ocean. Offshore
does not have a substantial effect on the internal forces within the pi- Arct. Eng. vol. 5, (2010), pp. 591–599, https://doi.org/10.1115/OMAE2010-
20552.
peline, however, it highly affects the pipeline embedment. Also it was [8] P. Xie, Q. Yue, A.C. Palmer, Cyclic plastic deformation of overbend pipe during
observed, that for higher lay tension values, the assumption of constant deepwater S-lay operation, Mar. Struct. 34 (2013) 74–87, https://doi.org/10.1016/
tension force along the pipeline [4,5] becomes more reasonable. j.marstruc.2013.08.003.
[9] J. Nocedal, S.J. Wright, Numerical Optimization, 2nd eds., World Scientific, 2006.
Comparison with Abaqus, beside the good agreement, showed that for [10] Z. Kang, L. Zhang, X. Zhang, Analysis on J lay of SCR based on catenary and large
the parameters given in Table 1, the assumption that the pipeline is deflection beam theory, Ocean Eng. 104 (2015) 276–282, https://doi.org/10.1016/
inextensible and unshearable is very reasonable. j.oceaneng.2015.05.017.
[11] P. Wriggers, Computational Contact Mechanics, 2nd eds., Springer-Verlag, Berlin
Direct minimization of a total potential energy is valid only for
Heidelberg, 2006, https://doi.org/10.1007/978-3-540-32609-0.
conservative systems, where after the loading removal both the seabed [12] B. Audoly, Y. Pomeau, Elasticity and Geometry: From Hair Curls to the Non-linear
and the pipeline return to their undisturbed initial state. Consequently, Response of Shells, OUP Oxford, 2010, https://books.google.co.il/books?id=
in order to demonstrate direct energy minimization, in this paper, we FMQRDAAAQBAJ.
[13] J.N. Reddy, An Introduction to Nonlinear Finite Element Analysis: With
have considered only conservative system. Nevertheless, the presented Applications to Heat Transfer, Fluid Mechanics, and Solid Mechanics, Oxford
technique can be extended to non-conservative cases as well, by University Press, 2015, https://books.google.co.il/books?id=6k5VBQAAQBAJ.
adapting special mechanisms to account for time-history and energy [14] M. Randolph, S. Gourvenec, Offshore Geotechnical Engineering, Taylor & Francis,
2011, https://books.google.co.il/books?id=qI8xPQAACAAJ.
dissipation. Moreover, for non-conservative cases (e.g., plastic seabed) [15] Dassault Systems, Abaqus Users’ Manual; Version 6.16. Providence, RI: Dassault
the presented technique of contact detection may become even more Systems Simulia Corp. (2016).
useful. [16] M.F. Randolph, D. White, Pipeline embedment in deep water: processes and
quantitative assessment, Off Shore Technol. Conf. (2008), https://doi.org/10.4043/
OTC-19128-MS.
References [17] H. Shiri, Response of steel catenary risers on hysteretic non-linear seabed, Appl.
Ocean Res. 44 (2014) 20–28, https://doi.org/10.1016/j.apor.2013.10.006.
[18] C.P. Aubeny, G. Biscontin, Seafloor-riser interaction model, Int. J. Rock Mech. Min.
[1] Q. Bai, Y. Bai, Subsea Pipeline Design, Analysis, and Installation, Elsevier Inc, 2014,
Sci. 9 (2009) 133–141, https://doi.org/10.1061/(ASCE)1532-3641(2009)
https://www.sciencedirect.com/science/book/9780123868886.
9:3(133)}.
[2] B. Guo, S. Song, A. Ghalambor, T.R. Lin, Offshore Pipelines: Design, Installation,
[19] A. Palmer, Touchdown indentation of the seabed, Appl. Ocean Res. 30 (2008)
and Maintenance, 2nd eds., Elsevier Science, Boston, 2014https://books.google.co.
235–238, https://doi.org/10.1016/j.apor.2008.09.004.
il/books?id=8r6tdixGNNwC.
[20] E.T.R. Dean, Offshore geotechnical engineering, Thomas telford ltd, London, 2010,
[3] S. Lenci, M. Callegari, Simple analytical models for the J-lay problem, Acta Mech.
https://doi.org/10.1680/oge.36413.
178 (2005) 23–39, https://doi.org/10.1007/s00707-005-0239-x.
[21] W.O. McCarron, Deepwater Foundations and Pipeline Geomechanics, J. Ross Pub.,
[4] F. Yuan, L. Wang, Z. Guo, Y. Xie, Analytical analysis of pipeline-soil interaction
2011, https://books.google.co.il/books?id=kdaEC_yuItIC.
during J-lay on a plastic seabed with bearing resistance proportional to depth, Appl.
[22] L. Wang, F. Yuan, Z. Guo, L. Li, Numerical analysis of pipeline in J-lay problem, J.
Ocean Res. 36 (2012) 60–68, https://doi.org/10.1016/j.apor.2012.02.004.
Zhejiang Univ. Sci. A. 11 (2010) 908–920, https://doi.org/10.1631/jzus.
[5] F. Yuan, M.F. Randolph, L. Wang, L. Zhao, Y. Tian, Refined analytical models for
A0900773.
pipe-lay on elasto-plastic seabed, Appl. Ocean Res. 48 (2014) 292–300, https://doi.
org/10.1016/j.apor.2014.10.003.

62