Influence of seabed proximity on the vibration responses of a pipeline accounting for fluid-structure interaction

Mechanical Systems and Signal Processing 114 (2019) 224–238

Contents lists available at ScienceDirect

Mechanical Systems and Signal Processing journal homepage: www.elsevier.com/locate/ymssp

Influence of seabed proximity on the vibration responses of a pipeline accounting for fluid-structure interaction H. Ribeiro Neto a, A. Cavalini Jr. a, J.M. Vedovoto a,⇑, A. Silveira Neto a, D.A. Rade b a b

UFU – Federal University of Uberlândia, Av. João Naves de Ávila, 2121, Campus Santa Mônica, Uberlândia, MG CEP 38400-902, Brazil ITA – Technological Institute of Aeronautics, Praça Marechal Eduardo Gomes, 50, Vila das Acacias, São José dos Campos, SP 12228-900, Brazil

a r t i c l e

i n f o

Article history: Received 23 September 2017 Received in revised form 20 April 2018 Accepted 3 May 2018

Keywords: Fluid-structure-interaction Vortex-induced-vibration Timoshenko beam theory

a b s t r a c t Cylindrical bodies subjected to external flow can vibrate due to the fluctuations of the forces induced by vortex shedding. The way these coherent fluid flow structures are formed and how they excite the structure depends on parameters, such as the Reynolds number, the reduced velocity, and the geometry e.g., the proximity of the structure to other bodies. These vibrations change the drag and lift forces by means of a nonlinear interaction. In addition, vibrations can cause crack nucleation and propagation in the structure. This is especially important when oil or natural gas is being transported in pipe-like structures, subjected to waves and sea currents. The present paper aims to characterize the influence of the proximity of the seabed on the fluid–structure interaction, considering horizontal pipes anchored by dunes. The simulations were undertaken for a nominally horizontal, elastic pipeline, 42 m in length and 0.273 m in diameter, with a mid–span static sag of 1.06 m due to self-weight. Seven different distances between the pipeline and the seabed were tested. The structural and fluid-dynamic models were coupled numerically, which allows the simulation and analysis of the flow using a single computational tool. The equations modeling the flow were solved in an Eulerian domain, while the surface of the immersed body was represented by a set of Lagrangian points. The immersed boundary method was used to impose a Dirichlet boundary condition on the Eulerian domain at the boundary between the structure and the fluid. It was also used to determine the fluid dynamic forces acting on the structure. An in-house three-dimensional computational framework was developed to simulate the turbulent incompressible flow subjected to fluid–structure interaction in conjunction with a beam modeled according to Timoshenko’s theory. The obtained results are consistent, as expected for this problem. Ó 2018 Elsevier Ltd. All rights reserved.

1. Introduction The interaction between flows and structures is a complex and recurring problem in engineering applications. This phenomenon can be found in aircraft [1,2], pipelines [3–5], wind turbines [6–8], bridges [9,10], offshore platforms [11–13], compressor valves [14,15], and others [16]. Computational fluid dynamics, together with techniques for the numerical solution of the equations that model the movement of structures, have been frequently used to solve fluid–structure interaction prob-

⇑ Corresponding author. E-mail address: [email protected] (J.M. Vedovoto). https://doi.org/10.1016/j.ymssp.2018.05.017 0888-3270/Ó 2018 Elsevier Ltd. All rights reserved.

H. Ribeiro Neto et al. / Mechanical Systems and Signal Processing 114 (2019) 224–238

225

lems. It is a typical multidisciplinary problem, since it involves disciplines such as fluid mechanics, structure mechanics, pure and applied mathematics, software engineering, and computer science. Flows over cylindrical structures can be a source of vibrations induced by instabilities. These vibrations can induce increase in fluid dynamic forces (i.e., drag and lift), leading to increased stress on the structures [17]. In addition, vibrations can cause nucleation and crack propagation in the structure leading to fatigue failure. In some cases, the RMS value of the lift coefficient can vary from 0.3 to 1.75, depending on the operating regime [17]. This justifies the concern with fluid–structure interaction in cylinders. This is especially important when the cylinders are pipelines, excited by waves and/or currents, used to transport oil or natural gas. The maintenance of this type of structure is generally expensive, since they can be immersed at a depth of thousands of meters, and failure has the potential to cause environmental disasters and great financial losses. Therefore, it is important to understand how the fluid–structure interaction process acts on the pipe dynamics to prevent failures. The main objective of this study is to acquire and expand the understanding of the influence of the seabed on the fluid–structure interaction process in horizontal pipelines anchored by dunes. This study was done through the numerical solution, in a distributed processing environment, of the equations that model the phenomenon. These simulations were done in a pipe of length L ¼ 42 m and diameter / ¼ 0:273 m in flows dynamically characterized by Re/ ¼ 8:64  104 . The seabed was modeled as a plane wall. Different mid–span gap distances have been tested, namely: 0.1/, 0.2/, 0.3/, 0.5/, 1/, 2/, and 5/. The results presented herein were obtained from the in-house numerical framework Fluids3D, which has been being developed in the MFLab for more than nine years. MFLab is the laboratory of fluid mechanics specialized in computational fluid dynamics located at the Federal University of Uberlândia in Brazil. In this program it is possible to perform simulations of incompressible flows taking into account the movement of structures by using a single computational tool [18]. The Immersed Boundary Method, used in the present research, is particularly suitable for problems involving fluid–structure interaction, as it allows to treat the fluid and structure domains independently [19]. The equations that model the flows are solved in a fixed Cartesian Eulerian domain, while the surface of the immersed body is represented by a set of Lagrangian points [20]. Through this method, the forces at the interface between the structure and the fluid are evaluated and used both to impose the non-slip boundary condition and to calculate the displacements and velocities of the structure. 2. Methodology The seabed geometry in deep water can be irregular. Pipes used to transport oil and natural gas laying on this uneven ground may have unsupported sections, known as free spans. These may be subjected to sea currents which, when interacting with the pipes, induce fluctuations of the fluid dynamical forces that excite the structure. The purpose of this section is to present the physical, mathematical, and numerical models used to simulate the problem of interest. 2.1. Physical model In the physical model, the problem of interest is evaluated and physical assumptions are adopted to make the solution feasible. The assumptions for each subsystem will be presented separately. 2.1.1. Fluid subsystem In this subsection, the physical assumptions for the fluid subsystem are presented. They are separated into three subgroups: domain, flow, and physical properties. The domain of the physical problem is the ocean, in which the structure is immersed. It is important to note that the domain is delimited by the seabed, which although static, influences the flow and, as a consequence, the vibration of the structure. Since the simulation of the ocean in all its extension is impractical, it is necessary to choose a reduced domain for the analysis of the problem. A domain with length 42 m (154/), 6 m wide (22/), and 4 m high (14.6/) was chosen, where, / is the diameter of the pipe. An illustration of the fluid domain with the structure deformed by its self weight is presented in Fig. 1 for the case in which the mid–span gap is 5/. For the values of the material properties used in the computations, the mid–span pipe static sag is 1.076 m or 3.941/ for all mid–span gap values. A supplementary view of Fig. 1 is presented in Fig. 2. A view looking downstream of the pipeline for a mid–span gap of 0.1/ showing mid–span deflection to scale. It is possible to see how little of the pipeline is actually in close proximity to the seabed. The flow inside the structure is not considered here. The internal fluid is modeled as a rigid body and accounted on the inertia and self-weight calculation of the structure. The external flow and its influence on the structure is considered. A preliminary analysis of the external flow is required to make the appropriate assumptions for the case study. The flow may be characterized dynamically by the Reynolds number, Re/ ¼ u/=m, where u is the mean fluid velocity at the inlet of the domain and

m is the kinematic viscosity of the fluid. For the determination of the Reynolds number, u ¼ 0:5 m s1, m ¼ 1:58  106

m2 s1 and / ¼ 2:73  101 m were adopted, resulting in Re/ ¼ 8:64  104 . It is well known that flows over cylinders become unstable over Re/ ¼ 47:5 [21]. Therefore, it is safe to assume that the flow is turbulent downstream from the pipe. This Reynolds number is found on usual industrial cases. If the flow velocity increases, the boundary layer on the surface of the struc-

226

H. Ribeiro Neto et al. / Mechanical Systems and Signal Processing 114 (2019) 224–238

Fig. 1. Fluid domain with the structure deformed by its self weight and sea bed modeled as flat wall. Fluid reference axes (X f ; Y f ; Z f ) and the structure (X s ; Y s ; Z s ) visible.

Fig. 2. Downstream view of the fluid domain with the structure deformed by its self weight and sea bed modeled as flat wall. Case using mid–span gap of 0.1 /.

ture may transition to turbulence, changing the dynamic of the fluid–structure-seabed interaction. On stationary cylinder the transition to turbulent boundary layer occurs on flows dynamically characterized by Reynolds numbers of Re ¼ 2:50  105 [22]. In turbulent flows, the imposed initial and boundary conditions are important, since errors in these conditions can be amplified exponentially by nonlinear interactions, generating instabilities [23]. Any variation in the initial conditions may determine different states. The thickness and the velocity gradient in the boundary layer formed near the seabed are factors that modify the pressure distribution and the hydrodynamic forces on stationary cylinders [24]. This effect is not taken into account in the present paper. The inlet condition defined on the plane X f ¼ 0:0 m was imposed as a uniform flow of (0.5 m s1, 0.0 m s1, 0.0 m s1) and the initial condition was also imposed as (0.5 m s1, 0.0 m s1, 0.0 m s1) throughout the computational domain. Further details of the boundary conditions will be provided in Section 2.2. The physical properties of the fluid are considered constant. In addition, it is assumed that saltwater behaves as a Newtonian fluid. The dynamic viscosity was held at

l ¼ 1:62  103 Pa s and the density at qf ¼ 1:025  103 kg m3.

2.1.2. Structure subsystem The physical and geometrical properties of the pipeline and the dunes will now be presented. Pipes used for oil and natural gas transportation can be kilometers long. However, the calculation domain must be reduced to enable numerical simulation. Therefore, the domain chosen for the structure subsystem has length 42 m, which is the length of the computational domain chosen for the fluid. The pipe has a circular section of external diameter / ¼ 2:73  101 m and internal diameter /i ¼ 2:37  101 m, which are constant. The pipeline is made of steel and is covered with a layer of anti-corrosive material with density adopted as qac ¼ 950 kg m3. The pipe is also filled with oil with density adopted as qoil ¼ 940 kg m3. The layer of anti-corrosive and the oil are considered to be rigid bodies and are accounted in the effects of the inertia of the pipe and the self-weight calculation. The total mass of anti-corrosive material and oil are 110.628 kg and 1734.684 kg respectively. The physical properties of the structure are considered constant over the domain. It is assumed that the steel behaves as a linear material with modulus of elasticity E ¼ 2:07  1011 Pa, Poisson coefficient ms ¼ 3  101 , and transversal modulus G ¼ 7:962  1010 Pa. The steel density adopted is qs ¼ 7:850  103 kg m3.

227

H. Ribeiro Neto et al. / Mechanical Systems and Signal Processing 114 (2019) 224–238

The dunes on which the structure is supported may consist of solid rocks or accumulated sediments. They are modeled as linear springs with an elastic constant of 1  1010 N m1. Residual lay tension and external and internal fluid pressure causes a resultant traction axial force of magnitude F a ¼ 58:584 kN. This force increases the stiffness of the structure. 2.2. Mathematical model In this mathematical model, the differential equations take into account the simplifications discussed in Section 2.1. The models for the fluid and structure subsystems will now be presented. 2.2.1. Fluid subsystem For the fluid subsystem, the simplified mass balance equations for incompressible flow, and the simplified linear momentum balance equations for isothermal, constant physical properties, and incompressible flows were used. Eq. (1) expresses the simplified mass balance equations for the case of incompressible flow in Cartesian coordinates, using Einstein’s summation convention:

@uj ¼ 0; @xj

j ¼ 1; 2; 3;

ð1Þ

where xj are the coordinate directions X f ; Y f and Z f , respectively, and uj is the component of the velocity of the fluid in the direction j. Eq. (2) expresses the balance of linear momentum for a Newtonian fluid written in divergent form, simplified for the case of constant specific mass and written in Cartesian coordinates:

@ui @ui uj 1 @P @ ¼ þ þ @t @xj qf @xi @xj


 @ui @uj ; m þ @xj @xi

ð2Þ

where i; j ¼ 1; 2; 3; P is the pressure, and t is the time. Eqs. (1) and (2) are considered sufficient for the simulation of the fluid dynamics of the flows of interest. However, to use these equations, it is necessary to resolve all scales or eddies of the flow. This is impracticable, since the number of degrees of freedom for a turbulent flow is estimated at approximately 1:28  1011 /3 , present in each unit volume /3 . The computational domains of the fluid subsystems have 49,515/3 . Hence, it would be necessary to solve 6:2279  1015 simultaneous algebraic equations, which is an unpractical task nowadays. Therefore, it is necessary to use a closure model designed for turbulent flows. In the present paper, the Large Eddy Simulation method is chosen. In this method, the large scales are resolved and the interaction between these and the small scales is modeled. In the present research, the subgrid model described in [25] was used. In order to use this method, it is necessary to filter Eqs. (1) and (2) twice, obtaining the transport equations of the filtered velocities, which correspond to the large scales that are resolved. The following equation is obtained by modeling the subgrid Reynolds tensor under the assumption of turbulent viscosity of Boussinesq

@ui @ui uj 1 @p @ þ ¼ þ @t @xj qf @xi @xj



mef


 @ui @uj : þ @xj @xi

ð3Þ

Here, mef is the effective kinematic viscosity, evaluated as the sum of the molecular kinematic viscosity and the turbulent kinematic viscosity mt . The turbulent kinematic viscosity is obtained from the Germano model [25]. Eqs. (1) and (3) have been solved numerically. The discretization procedures adopted will be discussed in Section 2.3. The initial and boundary conditions used are shown in Tables 1 and 2, respectively. In addition to the boundary conditions described in Table 2, it is necessary to set the fluid boundary condition at the interface between the fluid and the structure. An imposed velocity condition (Dirichlet) was employed. The fluid velocity at the surface of the pipe matches the structure velocity itself, given by the structure finite element model at the time analyzed. 2.3. Numerical model In the numerical model, the problem formulation under study is presented taking into account the simplifications discussed in the physical model (Section 2.1) and in the mathematical model represented by differential equations (Section 2.2). The numerical methods used to obtain the discretized equations of the fluid and the structure subsystems are presented next.

Table 1 Initial conditions for the fluid subsystem. Velocities (m s1)

Initial conditions

Pressure (Pa)

u

v

w

0.5

0.0

0.0

0

228

H. Ribeiro Neto et al. / Mechanical Systems and Signal Processing 114 (2019) 224–238

Table 2 Boundary conditions for the fluid subsystem. Velocities (m s1) and gradients (s1)

Boundary condition location

Pressure gradient (Pa m1)

u

v

w

X f = 0 m (Inlet) (Dirichlet)

0.5

0.0

0.0

X f = 6 m (Outlet) [26]

Wave equation

Wave equation

Wave equation

Y f = 0 m (bottom plane) (Dirichlet)

0.0

0.0

0.0

Y max (Top plane) (Symmetry) f

@u @y

¼0

Z f = 0 m (Symmetry)

@u @z @u @z

¼0

Z f = 42 m (Symmetry)

¼0

0.0 @v @z @v @z

@w @y

¼0

¼0

0.0

¼0

0.0

@p @x @p @x @p @y @p @y @p @z @p @z

¼0 ¼0 ¼0 ¼0 ¼0 ¼0

2.3.1. Fluid subsystem In general terms, the numerical method chosen for solving the momentum and Poisson equations for pressure correction is based on a three-dimensional, staggered, finite-volume discretization. The central difference scheme (CDS) is applied to express the diffusive contributions of the transport equations. For the advective terms, the discretization was done through the CUBISTA method [27]. This method is a Total Variation Diminished (TVD) scheme, which is a concept introduced in [28]. This is an important property of a discretization method for advective terms, as it helps to control and prevent numerical instabilities. A semi-implicit approach using the Modified Crank–Nicolson–Adams–Bashfort (MCNAB) method has been adopted [29,30], and the resulting linear systems are solved using the Modified Strongly Implicit Procedure (MSIP) [31] for velocity components. Since the numerical code developed adopts a pressure–velocity coupling based on pressure variations, an appropriate algorithm for pressure–velocity coupling is needed. Here, a projection method based on the fractional steps technique [32] is used, resulting in a variable-coefficient Poisson equation which is solved with the Bi-Conjugate Gradient Stabilized solver, BICGSTAB [33]. Eq. (4) shows Eq. (3) discretized in time with the MCNAB method:

   1  nþ1 1 @pnþ1 1  ui  ui n ¼  þ ð2 þ xnþ1 Þf ðun Þ  xnþ1 f un1 2 Dtnþ1 qf @xi þ

     1 ð8xnþ1 þ 1Þg unþ1 þ ð7xnþ1  1Þg ðun Þ þ xnþ1 g un1 ; 16xnþ1

ð4Þ

where xiþ1 ¼ DDtnþ1 ; Dtn is the time step evaluated in step n, the function f refers to the discretized diffusive term, and the tn function g represents the discretized advective terms. In Section 2.3.2, the Immersed Boundary Method will be described. 2.3.2. The immersed boundary method The immersed boundary method is necessary to impose the boundary condition on the fluid subsystem at the fluid–structure interface. In addition, it is used to determine the fluid dynamics forces acting on the structure. In this section we will discuss the method of imposing the boundary conditions and in Section 2.3.4 will show how the fluid dynamics forces are obtained. The immersed boundary method used here is an adaptation of that of [20]. The fluid domain is discretized in a uniform Cartesian mesh, while the surface of the structure is represented with Lagrangian points. These points are free to move within the fluid domain. The fluid velocity at the Lagrangian points (the interface between the fluid and the structure) should be equal to the velocity of the material points which belong to the structure. Both the velocity and position of these Lagrangian points are variable in time. Therefore, it is necessary to use a method that locates the Lagrangian points in the fluid domain and imposes the local velocity of the structure on the fluid. The surface of the pipe is discretized in 20,986 Lagrangian points represented by the same number of triangles given in an STL file [34]. This quantity remains constant in the course of the simulation. To take into account the presence of the immersed boundary, it is necessary to add a source term f i [N m3] to Eq. (3):

@ui @ui uj 1 @p @ ¼ þ þ @t @xj qf @xi @xj where



mef


 @ui @uj f þ i; þ @xj @xi qf

 Z

 ! ! ! ! ! fi x ¼ ðF k Þi xk dð x  xk Þdxk ; X

ð5Þ

ð6Þ

! ! ! ! with dð x  xk Þ being the Dirac delta function, xk location of the relevant Lagrangian point, and x the position of the relevant Eulerian point. Eq. (4) is rewritten taking into account the source term as follows:

   1  nþ1 1 @pnþ1 1  ui  ui n ¼  þ ð2 þ xnþ1 Þf ðun Þ  xnþ1 f un1 2 Dtnþ1 qf @xi      f 1 f þ ð8xnþ1 þ 1Þg unþ1 þ ð7xnþ1  1Þg ðun Þ þ xnþ1 g un1 þ i ¼ RHS þ i : 16xnþ1 qf qf

ð7Þ

H. Ribeiro Neto et al. / Mechanical Systems and Signal Processing 114 (2019) 224–238

229

An auxiliary term can be added and subtracted from Eq. (7) as follows:

 1  nþ1 f ui þ ui nþ1  ui nþ1  ui n ¼ RHS þ i : Dtnþ1 qf

ð8Þ

It can then be decomposed as

 1   nþ1 u  ui n ¼ RHS; Dtnþ1 i

ð9Þ

 f 1  nþ1 ui  ui nþ1 ¼ i : Dtnþ1 qf

ð10Þ

Now, Eq. (10) was written for the Eulerian domain, but it is also valid for any Lagrangian particle within the computational domain. Therefore, Eq. (10) can be rewritten for any Lagrangian particle:

 F 1 nþ1 i Ui  U i nþ1 ¼ ; Dtnþ1 qf

ð11Þ

where U i nþ1 is the velocity of a Lagrangian point, U i nþ1 is the Eulerian predicted velocity interpolated on the same point and F i is the source term to be used in Eq. (6). From the numerical point of view, a discrete version of the convolution integral in Eq. (6) is necessary. In the present paper, the hat function of [35] and implemented in the Fluids3D framework by [34] will be used. More details on its implementation and verification can be found in [34]. 2.3.3. Structure subsystem In this section, the discretized equations of a Timoshenko beam are presented in the physical domain. It will be also shown how to obtain the forces applied to the elements of the structure (Section 2.3.4). The pipe is modeled as a circular cross-section beam with constant internal and external diameters. The finite element discretized equations are based on [36]. In the case of the structure subsystem, the finite element discretized equations are obtained from the evaluation of the kinetic and potential energies of the beam elements. Here, the elementary stiffness matrix takes into account the increase of stiffness due to an imposed axial force. The resulting system of global equations can be written as:

€ þ C G q_ þ K G q ¼ W þ F; MG q

ð12Þ

where M G is the global mass matrix, C G is the global damping matrix, and K G is the global stiffness matrix. The vector of displacements is represented by q and the weight by W. F represents the fluid dynamics forces. The damping matrix C G is expressed as a linear combination of the mass and stiffness matrices:

C G ¼ cM þ bK;

ð13Þ

where c and b must be assumed. The adopted values are c ¼ 1:0 and b ¼ 1:0  104 . This assumption, although common in forced vibrations, is rather problematic for problems of stability. For these problems, the stability can be extremely sensitive to the structure damping matrix [37,38]. It was used the conversion of the second-order system of ordinary differential equations (ODEs), described in physical coordinates, to modal coordinates which are also time-dependent [39,40]. In order to numerically solve the resultant set of ODEs, the state-space form was used to reduce the order of the ODEs from second to first order. In the present paper, the Runge–Kutta–Fehlberg integrator is used to numerically integrate the resultant system of equations. The dynamic models represented in state-space and some related subjects are treated in detail in [41]. The initial condition is imposed using:

q ¼ K 1 G W; _q ¼ 0:0:

ð14Þ

The support conditions are modeled as linear springs with an elastic constant of 1  1010 N m1. Those springs are imposed on the X s and Y s directions. There is no degree of freedom on the axial direction (Z s ) therefore, no support condition is required. 2.3.4. Obtaining the fluid dynamics forces through the immersed boundary method In Section 2.3.2 there was presented the immersed boundary method for the imposition of the boundary condition on the fluid at the interface between the fluid and the structure. In the present subsection, the method for obtaining the force at the Lagrangian points will be presented. To do so, we use the source term F i [N m3] given by Eq. (11). In order to obtain the force itself, in Newtons, it is necessary to use a transformation volume F lag ¼ F i in which is the transformation volume

230

H. Ribeiro Neto et al. / Mechanical Systems and Signal Processing 114 (2019) 224–238

calculated by ¼ At d, where At is the area of the triangular element of the Lagrangian mesh and d is the average of the lengths of the edges of the triangular element. Each Lagrangian point will have a force applied to it, which must be transferred to the nodes of the finite element described in Section 2.3.3. In order to do this, only the Lagrangian points contained in the same structural element are considered. The force along the direction X s of each Lagrangian element is equally divided between the two nodes of the structural element. The same is true for the forces along the direction Y s . This procedure is performed for all Lagrangian points that define the immersed structure. In this way, the algorithm identifies the position of each Lagrangian point, identifies the element to which it is inserted, as well as the structural nodes that define this element, and distributes the forces to the structural nodes. 2.3.5. Fluid–structure coupling In Sections 2.3.1 and 2.3.3 there was presented the numerical methods used for the computation of the solutions for the fluid and structural subdomains, respectively. It is now necessary to couple the subsystems so that the complete solution of the fluid–structure interaction problem can be achieved. There are basically two ways of coupling the subsystems. The monolithic coupling (outside the scope of the present paper) and partitioned coupling [19]. In the monolithic coupling, the fluid and structure subsystems are treated as a single nonlinear system, while in the partitioned coupling, the equations are solved separately. We use the weak coupling, that is a type of partitioned coupling, which consists of treating the fluid subsystem with the frame position at the time n (Dn ), obtaining the force acting on the structure at time n þ 1 (F nþ1 ). With this force, the structure subsystem is then treated, obtaining Dnþ1 . This position is then used to solve the equations for the fluid subdomain and obtain F nþ2 , and so on. This type of coupling has an error of OðDtÞ, according to [42]. In addition, according to [43], when used in an incompressible flow formulation, it may exhibit numerical instability, called the artificial added mass effect. 3. Results This section compares some results obtained in the simulations. Seven cases with different distances to the seabed (henceforth referred to as the mid–span gap) have been simulated: 0:1/; 0:2/; 0:3/; 0:5/; 1/; 2/, and 5/. The centerline of the pipe was positioned in X f ¼ 1:72 m for all mid–span gaps and in Y f = 1.2238 m, 1.2510 m, 1.2784 m, 1.3330 m, 1.4696 m, 1.7426 m and 2.5618 m for mid–span gaps 0:1/; 0:2/; 0:3/; 0:5/; 1/; 2/, and 5/ respectively. The Large Eddy Simulation (LES) method provides the transient behavior of the flow, which can be visualized in three dimensions in Fig. 3. The Q criterion [44] colored with u component of velocity (component in the direction X f described in Fig. 1) was used as a numerical tracer to visualize the flow in the case with Re/ ¼ 8:64  104 and mid–span gap 5/. It is worth mentioning that the main advantage of resorting to a method such as Large Eddy Simulation is its ability to predict

Fig. 3. Three-dimensional visualization of the flow over a horizontal pipeline of aspect ratio L=/ ¼ 154. This flow is dynamically characterized by Re/ ¼ 8:64  104 . The iso-Q criterion is iso  Q ¼ 0:514971 colored by u-component of velocity at t ¼ 50 s, mid–span gap = 5/.

H. Ribeiro Neto et al. / Mechanical Systems and Signal Processing 114 (2019) 224–238

231

the large coherent fluid flow structures that are responsible for the structural excitation. Due its three-dimensionality, in Fig. 3 the Kelvin–Helmholtz vortices, interspersed with longitudinal vortices, are clearly visible. The structural vibration can be visualized by the orbit generated by following the displacement of one point in the space over time. The orbits in Z s =L ¼ 0:2; Z s =L ¼ 0:4 and Z s =L ¼ 0:5 for the cases where Re/ ¼ 8:64  104 and mid–span gaps 0.1/ and 5/ are presented in Fig. 4. The influence of the fluid–structure–seabed interaction is remarkable. Fig. 2 illustrated how little of the pipeline is actually in close proximity to the seabed and the results of the orbits showed how profound the effect of this limited region of close proximity is on the structural dynamics. There is little difference between the orbit shape on different points of the structure. Fig. 5 presents a comparison between the maximum values of the in-line mean displacement and Fig. 6 in the cross direction. Fig. 7 presents a comparison between the maximum standard deviation (STD) of the in-line displacement and Fig. 8 in

Fig. 4. Orbit of the pipeline, for t from 50 s to 300 s. Positions relative to the reference axis of the structure (Fig. 1).

232

H. Ribeiro Neto et al. / Mechanical Systems and Signal Processing 114 (2019) 224–238

Fig. 5. Comparison between the mean displacement in the in-line direction of the industrial case simulations. Values normalized by diameter and sampled at the central point of the structure (Z s ¼ 21 m).

Fig. 6. Comparison between the mean displacement in the cross direction of the industrial case simulations. Values normalized by diameter and sampled at the central point of the structure (Z s ¼ 21 m).

Fig. 7. Comparison between the standard deviation values of the displacement in the in-line direction of the industrial case simulations. Values normalized by diameter and sampled at the central point of the structure (Z s ¼ 21 m).

H. Ribeiro Neto et al. / Mechanical Systems and Signal Processing 114 (2019) 224–238

233

Fig. 8. Comparison between the standard deviation values of the displacement in the transversal (cross) direction of the industrial case simulations. Values normalized by diameter and sampled at the central point of the structure (Z s ¼ 21 m).

Fig. 9. Fourier transform of velocity components and displacements at the central point of the structure (Z s ¼ 21 m). Gap 0.1/.

the cross direction. In Fig. 5, it is evident that the mean in-line displacement grows as the structure moves away from the seabed. In Fig. 6, it can be seen that as the distance of the pipe from the seabed is increased, there is an increase in the mean in-line displacement, up to a mid–span gap of 0.5/. With a mid–span gap larger than that, there is a downward trend. The maximum standard deviation (STD) of the displacement in the in-line direction is given in Fig. 7. For those cases with mid–span gaps from 0.1/ to 0.5/, there is an increase of the standard deviation, whereas for those cases with mid–span gaps from 0.5/ to 5/, there is a decrease. The STD of the displacement in the cross direction is given in Fig. 8. There is an increase in this deviation for greater mid– span gaps. The u velocity component (X f direction) and v (Y f direction) were sampled using a numerical probe placed in the fluid domain. The probe’s position for the cases with mid–span gaps 0.1/, 0.5/ and 5/ is displayed in Figs. 13–15, respectively, represented by the symbol ⁄. These figures show the vorticity field at a plane positioned at Z f ¼ 21 m and will be discussed later in this paper. ^ and v ^ ) and from In Figs. 9–11 are shown the discrete Fourier transforms obtained from the data acquired by the probes (u

 b dc the displacement of the structure, non-dimensionalized by the diameter, both in-line /il and cross-directional d/c at the center of the structure (Z s ¼ 21 m). These results are for the cases of mid–span gaps 0.1/, 0.5/, and 5/, respectively. The Strouhal number is a dimensionless parameter that expresses the vortex shedding frequency in the wake of bluff bodies. It is defined as

St ¼

f w/ U1

ð15Þ

234

H. Ribeiro Neto et al. / Mechanical Systems and Signal Processing 114 (2019) 224–238

Fig. 10. Fourier transform of velocity components and displacements at the central point of the structure (Z s ¼ 21 m). Gap 0.5/.

Fig. 11. Fourier transform of velocity components and displacements at the central point of the structure (Z s ¼ 21 m). Gap 5/.

where f w is the vortex shedding frequency, U 1 is the fluid velocity and St is the Strouhal number. For a stationary cylinder submitted to a flow at Reynolds number Re/ ¼ 8:64  104 , the Strouhal number is approximately 0.2 [45]. This value was obtained through experiments with a rigid smooth cylinder positioned in the tank in such a way that the borders did not influence the wake. Therefore, the approximate vortex shedding frequency would be 0.366 Hz if there were no fluid–structure–seabed interaction. A numerical pluck test was performed in order to determine the natural frequencies of the structure under water, which are different from those found in vacuum, due to the added mass effect [46]. A force of (0.0, 500.0, 0.0) kN was applied at the center of the structure (Z s ¼ 21 m) on one time step. On the other time steps only the fluid dynamics forces and the self weight were present. The position was sampled at the same point(Z s ¼ 21 m). Initially the fluid is set still until the force is applied to the structure when the fluid starts to move due to the fluid–structure interaction. In Table 3, a comparison between the natural frequencies obtained in vacuum and the results of the pluck test is presented. The natural frequencies in vacuum were obtained through the solution of the eigenvalue problem [39] and the natural frequencies under water were obtained through the pluck test. The pluck test was performed for 67 s, therefore the frequency resolution is equal to 1:49  10 2 Hz. The mode shapes were obtained through the solution of the eigenvectors problem. Fig. 12 presents the first 5 mode shapes on plane (Y s ; Z s ). In this case, the mode shapes on planes (Y s ; Z s ) and (Z s ; X s ) are equal due to the considered structure boundary conditions. Although the natural frequencies changed in water due to the added mass effect, the mode shapes were considered to be similar. As the first natural frequency of the structure is close to the vortex shedding estimated by the Strouhal number, the phenomenon called lock-in [17,47] can occur. Such a phenomenon occurs when the vortex shedding frequency matches one of the natural frequencies of the structure. In Fig. 11 it is seen that the frequency of vortex shedding matches the frequency of

H. Ribeiro Neto et al. / Mechanical Systems and Signal Processing 114 (2019) 224–238

235

Table 3 Comparison between the frequencies obtained in the pluck test with the natural frequencies of the structure in vacuum. Vacuum [Hz]

Pluck test [Hz]

0:42 1:48 3:24 5:69 8:85 12:70 17:24

0:26 0:93 2:02 3:60 5:53 8:11 10:72

Fig. 12. Illustration of the first 5 mode shapes of the structure being analyzed.

Fig. 13. Vorticity field at Z f ¼ 21 m. Gap 0.1/ case at time = 170 s.

vibration of the structure at a frequency very close to the first natural frequency of the structure obtained in the pluck test and shown in Table 3. Therefore, lock-in is observed in this case. The frequencies of structure displacements are correlated with those of vortex shedding. In the case with a mid–span gap of 0.5/, there are two frequencies more energized in the in-line direction: 0.292 Hz and 0.584 Hz. In this case, the lower frequency has a greater amplitude than the second. In the case with a mid–span gap of 5/, these are also the most energized, however the higher frequency has a greater amplitude than the other. In the case with a mid–span gap of 0.1/, it is not possible to identify a frequency more energized than the other. The frequencies below 0.4 Hz have more energy than the others. The presence of the seabed in the fluid–structure interaction becomes even more evident. In Fig. 15 it can be seen that the presence of seabed has little influence on the vortex shedding. When a flow instability is being formed above the pipe, the force in the in-line direction increases in the positive direction of X s and the force in the transverse direction increases in the positive direction of Y s . When this structure detaches, there is a decrease in the forces. When a flow instability is being formed under the pipe, the force in the in-line direction increases, pulling the structure in

236

H. Ribeiro Neto et al. / Mechanical Systems and Signal Processing 114 (2019) 224–238

Fig. 14. Vorticity field at Z f ¼ 21 m. Gap 0.5/ case at time = 170 s.

Fig. 15. Vorticity field at Z f ¼ 21 m. Gap 5/ case at time = 170 s.

the positive direction of X s , and the force in the transverse direction increases in the negative direction of Y s . When this structure is detached, the forces decrease. Therefore, when there is vortex shedding above and below the pipe, there is a double frequency in the in-line direction 0.584 Hz. In the transverse direction, a complete cycle vortex shedding above and below causes only a simple frequency of movement of the structure 0.292 Hz. When there is an asymmetry in the formation of flow instabilities, causing a greater force above or below, the simple frequency appears in the in-line direction, as in the case with a mid–span gap of 0.5/ as can be seen in Fig. 10. As discussed in [48], in cases with small mid–span gaps the formation of traditional flow instabilities does not happen: the suppression of vortex shedding is observed. In the case with a mid–span gap of 0.1/, this can be verified. This makes the movement of the structure very different from the traditional movement without the influence of the seabed, as can be seen in Figs. 4a, c, e, 9 and 13.

4. Conclusions The objective of the present paper was to contribute to the study of an industrial problem of fluid–structure interaction considering pipelines in the proximity of the seabed. In order to achieve such a goal, the in-house computational framework Fluids3D was used. The study of an industrial problem was made. It was possible to show the strict correlation between the dynamics of vortex shedding, and the dynamics of the immersed structure, as expected. Interesting anomalies were observed for certain values of the mid–span gap (distance between the structure and the seabed) due to an inherent nonlinearity. Different dynamic behaviors were identified in the movement of the structure. The orbits vary noticeably with the mid–span gap. The study of problems involving fluid–structure interaction must be done in such a way as to take into account the nonlinear interaction between the two subsystems involved: the fluid and the structure. In addition, the determination of the velocity and pressure fields must be done in three dimensions and using closure turbulence models, especially in industrial

H. Ribeiro Neto et al. / Mechanical Systems and Signal Processing 114 (2019) 224–238

237

cases in which the Reynolds number is large. A two-dimensional simulation of the flow around the structure would not take into account the effects of the transition to turbulence, such as the three-dimensionalization of the flow. Acknowledgements The authors would like to thank PROPP-UFU, CNPq, CAPES, FAPEMIG, and Petrobras S.A. for financial and material support. References [1] C. Farhat, K.G. van der Zee, P. Geuzaine, Provably second-order time-accurate loosely-coupled solution algorithms for transient nonlinear computational aeroelasticity, Comput. Methods Appl. Mech. Eng. 195 (17) (2006) 1973–2001. [2] H.H. Mian, G. Wang, Z.Y. Ye, Numerical investigation of structural geometric nonlinearity effect in high-aspect-ratio wing using CFD/CSD coupled approach, J. Fluids Struct. 49 (2014) 186–201. [3] J.P. Pontaza, R.G. Menon, On the numerical simulation of fluid-structure interaction to estimate the fatigue life of subsea pipeline spans: effects of wall proximity, in: Proceedings of the ASME 2010 29th International Conference on Ocean, Offshore and Arctic Engineering, vol. 6, 2010, pp. 817–827. [4] A.S. Borges, F.P. Mariano, J.M. Vedovoto, A. da Silveira Neto, D.A. Rade, Fluid-structure interaction of cylinders by combining the Cosserat beams theory and immersed boundary methodology, in: Argentina Mechanical Computation Conference, vol. 29, 2010, pp. 15–18. [5] L. Lee, D. Allen, J.P. Pontaza, F. Kopp, V. Jhingran, In-line motion of Subsea pipeline span models experiencing vortex-shedding, in: ASME 2009 28th International Conference on Ocean, Offshore and Arctic Engineering, vol. 5, 2009, pp. 269–276. [6] L. Wang, X. Liu, A. Kolios, State of the art in the aeroelasticity of wind turbine blades: aeroelastic modelling, Renew. Sustain. Energy Rev. 64 (2016) 195– 210. [7] R. Rafiee, M. Tahani, M. Moradi, Simulation of aeroelastic behavior in a composite wind turbine blade, J. Wind Eng. Ind. Aerodyn. 151 (2016) 60–69. [8] L. Wang, R. Quant, A. Kolios, Fluid structure interaction modelling of horizontal-axis wind turbine blades based on CFD and FEA, J. Wind Eng. Ind. Aerodyn. 158 (2016) 11–25. [9] X. Ying, F. Xu, M. Zhang, Z. Zhang, Numerical explorations of the limit cycle flutter characteristics of a bridge deck, J. Wind Eng. Ind. Aerodyn. 169 (2017) 30–38. [10] N. Lee, H. Lee, C. Baek, S. Lee, Aeroelastic analysis of bridge deck flutter with modified implicit coupling method, J. Wind Eng. Ind. Aerodyn. 155 (2016) 11–22. [11] A.J. Dunbar, B.A. Craven, E.G. Paterson, Development and validation of a tightly coupled CFD/6-DOF solver for simulating floating offshore wind turbine platforms, Ocean Eng. 110 (Part A) (2015) 98–105. [12] T.T. Tran, D.H. Kim, The coupled dynamic response computation for a semi-submersible platform of floating offshore wind turbine, J. Wind Eng. Ind. Aerodyn. 147 (2015) 104–119. [13] Y. Liu, Q. Xiao, A. Incecik, C. Peyrard, D. Wan, Establishing a fully coupled CFD analysis tool for floating offshore wind turbines, Renew. Energy 112 (2017) 280–301. [14] X. Yu, Q. Tan, Y. Ren, X. Jia, L. Jin, Numerical study of the reed valve impact in the rotary compressor by FSI model, Energy Proc. 105 (2017) 4890–4897 (8th International Conference on Applied Energy, ICAE2016, 8–11 October 2016, Beijing, China). [15] F. Barbi, J.L. Gasche, A. da Silveira Neto, M.M. Villar, R.S. de Lima, Numerical simulation of the flow through a compressor-valve model using an immersed-boundary method, Eng. Appl. Comput. Fluid Mech. 10 (1) (2016) 255–271. [16] E. Naudascher, D. Rockwell, Flow-Induced Vibrations: An Engineering Guide, Dover Publications, 2012. [17] M.J. Chern, Y.H. Kuan, G. Nugroho, G.T. Lu, T.L. Horng, Direct-forcing immersed boundary modeling of vortex-induced vibration of a circular cylinder, J. Wind Eng. Ind. Aerodyn. 134 (2014) 109–121. [18] A.S. Borges, Desenvolvimento de procedimentos de modelagem de interação fluido-estrutura combinando a Teoria de Vigas de Cosserat e a metologia de fronteira imersa, Phd thesis, Universidade Federal de Uberlândia, Uberlândia, Minas Gerais, Brazil, 2010. [19] F. Sotiropoulos, X. Yang, Immersed boundary methods for simulating fluid–structure interaction, Prog. Aerosp. Sci. 65 (2014) 1–21. [20] Z. Wang, J. Fan, K. Luo, Combined multi-direct forcing and immersed boundary method for simulating flows with moving particles, Int. J. Multiph. Flow 34 (3) (2008) 283–302. [21] B. Kumar, S. Mittal, Effect of blockage on critical parameters for flow past a circular cylinder, Int. J. Numer. Meth. Fluids 50 (8) (2006) 987–1001. [22] I. Rodríguez, O. Lehmkuhl, R. Borrell, L. Paniagua, C.D. Pérez-Segarra, High performance computing of the flow past a circular cylinder at critical and supercritical Reynolds numbers, Proc. Eng. 61 (2013) 166–172. [23] M. Damasceno, J. Vedovoto, A. Da Silveira-Neto, Turbulent inlet conditions modeling using large-eddy, Simulations 104 (2015) 105–132. [24] C. Lei, L. Cheng, K. Kavanagh, Re-examination of the effect of a plane boundary on force and vortex shedding of a circular cylinder, J. Wind Eng. Ind. Aerodyn. 80 (3) (1999) 263–286. [25] M. Germano, U. Piomelli, P. Moin, W.H. Cabot, A dynamic subgrid-scale eddy viscosity model, Phys. Fluids A 3 (7) (1991) 1760–1765. [26] I. Orlanski, A simple boundary condition for unbounded hyperbolic flows, J. Comput. Phys. 21 (1976) 251–269. [27] M. Alves, P. Oliveira, F. Pinho, A convergent and universally bounded interpolation scheme for the treatment of advection, Int. J. Numer. Meth. Fluids 41 (2003) 47–75. [28] A. Harten, High resolution schemes for hyperbolic conservation laws, J. Comput. Phys. 135 (1997) 260–278. [29] D. Wang, S. Ruuth, Variable step-size implicit-explicit linear multistep methods for time-dependent partial differential equations, J. Comput. Math. 26 (2008) 838–855. [30] J.M. Vedovoto, Mathematical and Numerical Modeling of Turbulent Reactive Flows using a Hybrid LES/ PDF Methodology, Phd thesis, Universidade Federal de Uberlândia, Uberlândia, Minas Gerais, Brazil and the École Nationale Supérieure de Mécanique et d’Aérotechnique–ENSMA, Poitiers, France, 2011. [31] G.E. Schneider, M. Zedan, A modified strongly implicit procedure for the numerical solution of field problems, Numer. Heat Transf. 4 (1981) 1–19. [32] A.J. Chorin, The numerical solution of the Navier-Stokes equations for an incompressible fluid, Bull. Am. Math. Soc. 73 (6) (1967) 928–931. [33] H.A. van der Vorst, Iterative solution methods for certain sparse linear systems with a non-symmetric matrix arising from PDE-problems, J. Comput. Phys. 44 (1) (1981) 1–19. [34] J.M. Vedovoto, R. Serfaty, A.D. Silveira Neto, Mathematical and numerical modeling of turbulent flows, Anais Acad. Bras. Ciências 87 (2015) 1195–1232. [35] B.E. Griffith, C.S. Peskin, On the order of accuracy of the immersed boundary method: higher order convergence rates for sufficiently smooth problems, J. Comput. Phys. 208 (1) (2005) 75–105. [36] A.A. Cavalini, Detecção e Identificação de Trincas Transversais Inicipientes em Eixos Horizontais Flexíveis de Máquinas Rotativas, Doctoral thesis, Universidade Federal de Uberlândia, Minas Gerais, Brazil, 2013. [37] P. Hagedorn, M. Eckstein, E. Heffel, A. Wagner, Self-excited vibrations and damping in circulatory systems, J. Appl. Mech. 81 (10) (2014) 1–9. [38] P. Hagedorn, E. Heffel, P. Lancaster, P.C. Müller, S. Kapuria, Some recent results on MDGKN-systems, ZAMM Z. Angew. Math. Mech. 95 (7) (2015) 695– 702. [39] S. Rao, Vibrações mecânicas, Prentice-Hall, 2009.

238

H. Ribeiro Neto et al. / Mechanical Systems and Signal Processing 114 (2019) 224–238

[40] N. Maia, J. Silva, E. Almas, R. Sampaio, Damage detection in structures: from mode shape to frequency response function methods, Mech. Syst. Sign. Process. 17 (3) (2003) 489–498. [41] K. Ogata, Engenharia de controle moderno, Prentice-Hall, 1998. [42] H.G. Matthies, R. Niekamp, J. Steindorf, Algorithms for strong coupling procedures, Comput. Methods Appl. Mech. Eng. 195 (17–18) (2006) 2028–2049. [43] C. Förster, W.A. Wall, E. Ramm, The artificial added mass effect in sequential staggered fluid-structure Interaction algorithms, in: European Conference on Computational Fluid Dynamics ECCOMAS CFD, 2006. [44] J. Jeong, F. Hussain, On the identification of a vortex, J. Fluid Mech. 285 (1995) 69–94. [45] F.M. White, Fluid Mechanics, fourth ed., McGraw-Hill, New York, 1998. [46] J.N. Newman, Marine Hydrodynamics, MIT Press, 1977. [47] S. Kaneko, T. Nakamura, F. Inada, M. Kato, Flow Induced Vibrations, Elsevier, Amsterdam, 2008. [48] C. Lei, L. Cheng, S.W. Armfield, K. Kavanagh, Vortex shedding suppression for flow over a circular cylinder near a plane boundary, Ocean Eng. 27 (10) (2000) 1109–1127.