Finite element approach of the buried pipeline on tensionless foundation under random ground excitation

Available online at www.sciencedirect.com

ScienceDirect Mathematics and Computers in Simulation 169 (2020) 149–165 www.elsevier.com/locate/matcom

Original articles

Finite element approach of the buried pipeline on tensionless foundation under random ground excitation Yang Shang Hsu Department of Mechanical Engineering, Polytechnic School at the Pontifícia Universidade Católica do Paraná, Brazil Received 30 September 2018; received in revised form 15 June 2019; accepted 9 September 2019 Available online 23 September 2019

Abstract This work presents an investigation in the numerical performance of finite element approach in the dynamic elastoplastic analysis of buried pipeline, which is subjected to random ground motion. The surrounding soil is modeled by Winkler and Pasternak type foundation, while the random ground motion is generated synthetically by generalized nonstationary Kanai–Tajimi model. The governing equation is formulated by Euler–Bernoulli beam theory and discretized by beam-pipe element. Moreover, the von Mises isotropic hardening model is employed for material behavior modeling. The Hilber–Hughes– Taylor (HHT) method and the Newton–Raphson method are adopted as the time incremental iterative algorithm to solve the global equilibrium equation. Several applications are carried out to investigate the numerical performance of the present numerical model in dealing with the dynamic elastoplastic analysis of buried pipe. The norm L 2 of stress and displacement error are determined for different time intervals and the factors that contribute to the error reduction are investigated in present work. c 2019 International Association for Mathematics and Computers in Simulation (IMACS). Published by Elsevier B.V. All rights ⃝ reserved. Keywords: Buried pipeline; Dynamic elastoplastic; Kanai–Tajimi model; Finite element; Norm L 2 of displacement and stress errors; Pasternak foundation; Winkler foundation

1. Introduction The seismic responses of the buried pipelines involve several factors and different areas of research [12,18]. Firstly, it should consider the characteristics of ground motion such as peak value, wave frequency, and type of wave. Secondly, the soil stiffness should be considered. Since the buried pipeline is surrounded by soil, therefore, the soil–structure interaction is also a crucial point to be included in the analysis. Finally, the buried pipeline analysis should be capable to predict the magnitude of the maximum strain and stress as they are generated by the seismic response. These factors should be considered in a numerical modeling of random seismic response of buried pipeline [13]. Moreover, there are 3 types of randomness that could be considered in a numerical modeling, such as the random variation of soil properties, the structural vibration driven by stochastic dynamic loads, and the random mechanical parameters of structural components. The last two types of randomness are extensively studied, and they are known as random vibration and stochastic structural analysis, respectively. In this work, only the randomness of structure response driven by stochastic dynamic loads is analyzed. E-mail address: [email protected]. https://doi.org/10.1016/j.matcom.2019.09.004 c 2019 International Association for Mathematics and Computers in Simulation (IMACS). Published by Elsevier B.V. All rights 0378-4754/⃝ reserved.

150

Y.S. Hsu / Mathematics and Computers in Simulation 169 (2020) 149–165

The present work aims to present an investigation concerning the numerical performance of finite element approach in the nonlinear time history analysis of buried pipeline, subjected to random ground motion. The seismic motion produced by an earthquake is modeled numerically as random ground motion in the transversal direction. To generate artificially the random ground motion, this work employs the generalized nonstationary Kanai–Tajimi model ([2,17,24,42,49]) and the ground acceleration profile is created by Shinozuka method. The soil–structure interaction is also included in the present work, and the present numerical approach includes two types of foundations, such as Winkler and Pasternak foundation [52]. The tensionless elastic–perfect-plastic discrete spring is implemented in the present work for Winkler foundation, and shear layer modulus is included for Pasternak foundation. Finally, the buried pipelines are discretized by using 3 nodes beam — pipe element, and Euler–Bernoulli beam theory [4] is employed in this work. The dynamic elastoplastic behavior of the buried pipeline is modeled by using von Mises isotropic hardening theory [32] by considering small displacement and infinitesimal strain, while the geometrical nonlinearity is not considered. A brief literature review is presented below and it can be divided into three parts. Firstly, the studies that address the seismic analysis of buried pipe are presented. A significant amount of studies on this area has already been developed, and this part of the review starts with the presentation of numerical approaches that were proposed to analyze different situations, which involve the seismic wave action. Such as pipeline crossing slip-fault area [22,51,53] or subjected to random ground motion. Secondly, a brief literature review is presented concerning the developed model in the soil–structure interaction. Several models were proposed to better represent the surrounding soil behavior in the ground motion. These models are formulated by considering different physical effects of soil on the buried pipeline. The solution of these model’s formulation can be carried out either by numerical approach or analytical approach. Finally, the literature review ends with the methodologies that are employed for the synthetic generation of ground random acceleration, such as improved Kanai–Tajimi method. Many numerical models were developed to carry out the time history seismic analysis of the buried pipeline. In the pioneer researches, the stress analyses were addressed, and the surrounding soil is modeled by a series of the spring element, which has the stiffness defined as a function of wavelength and frequency of the seismic wave [33,36]. The finite element approach also can be adopted by using shell element and beam element, and it is observed that the axial stress in the pipe is larger than hoop stress in the case of axial load. However, for the transversal loading, when the ratio between soil stiffness and pipe stiffness is small, the hoop stress is larger than the axial stress, and the 3D finite element is needed in this case. When the pipe stiffness is larger than soil stiffness, the axial stress has a more significant role than hoop stress, therefore, the employment of beam element in the numerical modeling is convenient [10,48], with the surrounding soil modeled by spring and dashpot element. In some studies, the random field theory is employed to simulate the non-homogeneous soil surrounding the buried pipe [35]. In this case, the Kelvin–Voigt model is adopted for soil–structure interaction. The interface element is constituted by spring and dashpots, representing stiffness and damping effect of soil. In addition, the simplified soil–structure interaction model was also proposed by using Winkler model, while the influences of pipe embedment and soil surrounding the buried pipeline response are analyzed [6,11]. Furthermore, there is a numerical model developed for the segmented buried pipe that is subjected to seismic wave, and it is modeled by the beam element [46]. In this study, the surrounding soil is modeled by elastoplastic spring. The proposed numerical model provides pipe strain and joint relative displacement with reasonable accuracy. Researches that adopt ground motion records are also developed by considering continuous pipeline. Further, the fragility function and the seismic intensity value of the failure were proposed [26]. Still, in this area of research, other studies propose inelastic cubic fiber element to simulate the pipeline, and the Winkler foundation is employed for surrounding soil [28]. Finite element approach was also developed in wave-seabed-pipeline interaction [19]. Moreover, the Navier–Bernoulli beam element formulated by Updated Lagrangian formulation was also proposed and it has considered geometrically nonlinear behavior of floating pipeline laying on a seabed [39]. Besides the finite element approach, the wave scattering mechanism was also developed, with the buried pipeline modeled as a waveguide, and the ground motion was converted to a waveform [54]. The stochastic seismic analysis based on the finite element approach is also widely developed. The power spectral density method can be adopted in the numerical approach for a buried pipeline that is rested on two parameters elastic foundation. While the stress developed in the pipeline is expressed in the domain of frequency [5]. Moreover, the joint element can be used to simulate the segmented pipes that are rested on a nonlinear spring element [31]. The probability density evolution method is also developed to evaluate the reliability in the pipeline network. The seismic stress on the pipe is used as a parameter to carry out the failure analysis [30]. This stochastic

Y.S. Hsu / Mathematics and Computers in Simulation 169 (2020) 149–165

151

methodology shows competitiveness in comparison with Monte Carlo method and the recursive decomposition algorithm. The soil–structure interaction plays a key role in the seismic analysis of the buried pipeline. Besides the analytical solution, the finite element approach of soil–structure interaction is also widely used to carry out the free vibration analysis of a beam resting totally or partially on a two parameters elastic foundation [55]. Moreover, this numerical approach is also applied for free vibration analysis of Euler–Bernoulli beam resting on a Winkler type stepped elastic foundation [50]. In addition, higher order plate element is also developed for analysis related to Winkler foundation [14,21,37,38]. Another numerical approach, such as boundary element based method is also developed to carry out the analysis of a shear deformable beam and inelastic Euler–Bernoulli beam resting on tensionless nonlinear viscoelastic foundation [23,43,44]. Many numerical approaches are also proposed using the differential quadrature element method [7] or the differential transform method [3]. Besides, the Runge–Kutta and Regula– Falsi methods are also developed to solve the problem of free vibration in prismatic Timoshenko beam resting on Pasternak foundation [27]. The time history of random ground motion can be obtained by using the Kanai–Tajimi model. This model has been investigated and developed during decades. Initially, the empirical formula of the spectrum is developed by considering the magnitude and the frequency of the spectrum in the statistical methodology [25]. After, the generalized nonstationary Kanai–Tajimi model is developed to represent the real accelerogram records in several earthquake sites [17,42]. And it is also applied to investigate the nonlinear response of the structure with multidegree-of-freedom [40]. Other spectral function, such as Normalized Power Spectral Density function (NPSD) is also proposed, which includes Kanai–Tajimi filter, Clough–Penzien filter and High-Pass filter. Recently, a modified Kanai–Tajimi model is proposed and it is named as fully nonstationary bimodal Kanai–Tajimi model. It is used to generate synthetically the time history of ground motion of Orion and Chi-Chi earthquake [8,9]. Besides the Kanai–Tajimi model, another methodology was also proposed [1,34], and it is based on wavelet packet procedure [29,45]. 2. Governing equations The structural model of the buried pipeline is presented schematically in Fig. 1(a) and (b), for seismic motion acting in the axial and transversal directions, respectively. Fig. 1(c) shows the force–displacement relationship of the spring element. The Euler–Bernoulli beam theory is employed to formulate the governing equation. The foundation model considers the spring element and the shear layer. The governing equation is shown below: ∂u ∂ 2u + ρ A 2 = f u (x, t) − qu (x, t) ∂x ∂t 2 ) 2 ( ∂ v ∂ 2v ∂ E (x, t) I 2 + ρ A 2 = f v (x, t) − qv (x, t) 2 ∂x ∂x ∂t

E (x, t) A

qu (x, t) = k1 u − G 1

∂ 2u ∂x2

(1a) (1b)

(2a)

∂ 2v (2b) ∂x2 where ρ is the pipeline mass density, u and v are the displacements in the axial and transversal directions, respectively. The A and I are area and moment of inertia, respectively. The E (x, t) is the material modulus, f u (x, t) and f v (x, t) are the loads derived from random ground motion in the axial and transversal directions, respectively. The qu (x, t) and qv (x, t) are reactions of foundations, k1 and k2 are the stiffnesses of spring element in the axial and transversal directions, respectively. Finally, G 1o and G 2o are the shear moduli of the shear layer, used in Pasternak type foundation. In this work, a conventional three nodes beam element is employed [4,15,56]. For dynamic elastoplastic behavior of the buried pipeline, the von Mises isotropic hardening is employed [32]. The weak form is adopted to develop the global equilibrium equation, which leads to: { } { } (3) [M] d¨ + [K ] {d} = {P(t)} = −M (1) d¨g qv (x, t) = k2 v − G 2

152

Y.S. Hsu / Mathematics and Computers in Simulation 169 (2020) 149–165

Fig. 1. Structural model of a buried pipeline (a) axial seismic motion and (b) transversal seismic motion, (c) force–displacement relationship of the spring element. Fy is foundation limit force.

where the element mass and stiffness matrices are defined as: ∫ ρ [N ]T [N ] d V [Me ] = ∫ ∫ Ve [ N L ]T [ ] [ ep ] T {σ } B1N L d V B2 [B] D [B] d V + [K e ] =

(4) (5)

Ve

Ve

In Eq. (4), [N ] is the shape function of Lagrange–Hermite polynomials. In Eq. (5), [B] is [ ] matrix,[ composed ] the strain–displacement matrix, B1N L and B2N L are nonlinear strain–displacement matrix and [D ep ] is the elastoplastic constitutive matrix [4]. The elastoplastic material behavior can be modeled by using incremental relations of stress. By using the indicial notation for Cartesian axes, the incremental stress is given by [4,32,47]: ep

dσi j = Di jkl dεkl

(6)

Eq. (6) is written incrementally for constitutive law, where the incremental strain components dεi j are defined as a ep sum of linear strain and nonlinear strain from the expansion of Taylor series. In Eq. (6), Di jkl is a tangential tensor defined by suitable state variables, such as constitutive tensor and effective stress, among others. Additionally, in ep this work, appropriate directions are also considered for Di jkl in unidimensional problems. Within the context of isotropic work hardening theory, the tangent constitutive tensor is defined as: ep

Di jkl = Di jkl − (1/γ ) Di jmn amn aop Dopkl

(7)

Where ( ) Di jkl = 2µν/ (1 − 2ν) δi j δkl + µ δik δ ji + δil δ jk

(8)

akl = ∂ σ¯ /∂σkl

(9a)

γ = ai j Di jkl akl + H

(9b)

H = ∂σo /∂ ε¯

(9c)

P

In Eqs. (9a)–(9c), σ¯ and ε¯ P are the effective stress and plastic strain, respectively; σ0 is the uniaxial yield stress; H is the plastic hardening modulus and µ and ν stand for the material shear modulus and the Poisson ratio, respectively. In the case of von Mises isotropic( strain-hardening material, the tensor of incremental elastoplastic material moduli ( )) ep takes the form Di jkl = Di jkl − 3µ/ σ02 (1 + H/3) si j skl , where si j = σi j − (1/3) δi j σkk is the stress deviator, and for the case of a perfectly plastic material H = 0. In the case of elastic analyses, the Cauchy stresses can be defined by σi j = Di jkl εkl , where Di jkl is the elastic constitutive tensor.

153

Y.S. Hsu / Mathematics and Computers in Simulation 169 (2020) 149–165

For the initial stress formulation, it is convenient to define a fictitious “elastic” stress increment as follows: dσiej = Di jkl dεkl

(10)

And to rewrite Eq. (6) as indicated below dσi j = dσiej − dσiPj

(11)

where the initial stress increments are computed by dσiPj = (1/γ ) Di jmn amn akl dσkle

(12)

The random ground acceleration is generated by generalized nonstationary Kanai–Tajimi model [17,42]. The global equilibrium equation is solved by using Hilber–Hughes–Taylor (HHT) time stepping algorithm [20]. In the present work, the elastoplastic criterion is evaluated in the cross-section of the pipeline. The numerical integration scheme is adopted in this work to evaluate stiffness and mass matrices, with seven points of numerical integration over each beam element. In each point of numerical integration, the elastoplastic criterion is evaluated in the crosssection, which is divided into layers and circumferential parts. The Newton–Raphson algorithm is adopted to solve the nonlinearity problem, which involves pipeline material isotropic hardening model and foundation perfect plastic effect. The flowchart, presented in Fig. 2, shows the numerical solution procedure. The solution procedure begins with the input data for buried pipeline and foundation mechanical properties, follows from input data to evaluate the synthetic accelerogram by using Kanai–Tajimi spectral density function. Then, the HHT algorithm is activated to calculate the nodal displacement vector. The nodal displacement vector is used to evaluate the strain, then the von Mises stress is calculated. After, the von Mises criterion is used to verify if the buried pipeline is undergoing material hardening. If affirmative, the Newton–Raphson algorithm is initiated to update the new stress state of the material, and the stiffness matrix is updated by changing the material modulus according to the stress–strain relationship. Otherwise, the algorithm proceeds the verification of loading and unloading conditions. The global convergence is verified in each numerical point of the beam element. This work uses a point mesh constituted by seven points of numerical integration over the beam element. In every numerical integration point, the hardening condition is verified in the cross-section, which is divided into one layer in the radial direction, and 36 parts in the circumferential direction. After this procedure, the residual force of the pipeline and foundation is calculated, and the convergence condition is verified. This work adopts the energy criterion to evaluate the convergence parameter with the tolerance value defined as 1 × 10−8 . This loop continues until the convergence condition is fulfilled. Finally, the time step is incremented, and Start In present work, an a-posteriori error estimator is adopted in a simple form. The norm in L 2 for displacement and stress errors is calculated, as stated in expressions (13) and (14) respectively. In addition, several refinements in time discretization are adopted in present work. Therefore, the error is measured between the different quantities of the degree of refinement in the finite element procedure discretization. This consideration is made due to the limitation to find literature results to be considered as reference results. [∫ ]1/2 ( )T ( ) ∥e∥ L 2 = u − uˆ u − uˆ dΩ (13) Ω

[∫ ∥eσ ∥ L 2 =

(

σ − σˆ

)T (

σ − σˆ dΩ )

]1/2 (14)

Ω

3. Generalized nonstationary Kanai–Tajimi model The synthetic accelerogram records in present work are generated by generalized nonstationary Kanai–Tajimi model. The original version of this model is proposed as a filtered white noise model in early 60 [24,25,49]. Since then, it is widely used in engineering application to evaluate synthetic accelerogram. The most attractive feature of this model is its capacity in the analysis of random response in seismically excited structures. It also provides a simple way to generate synthetic ground motion records with a single dominant frequency. In its original version, this model is a stationary random process, by this manner, the amplitude and frequency are both timeindependent. However, the available ground motion records show to be nonstationary in amplitude and frequency domain. Therefore, the generalized nonstationary model is proposed, that aims to transform the amplitude and

154

Y.S. Hsu / Mathematics and Computers in Simulation 169 (2020) 149–165

Fig. 2. Flowchart of the incremental-iterative solution procedure.

Y.S. Hsu / Mathematics and Computers in Simulation 169 (2020) 149–165

155

frequency used in the spectral function into the time-dependent variable [17,42]. The generalized nonstationary Kanai–Tajimi model for earthquake ground acceleration is given by: ( ) X¨ g = − 2ζg (t)ωg (t) X˙ f + ωg2 (t)X f e(t) (15) 2 ¨ ˙ X f + 2ζg (t)ωg (t) X f + ωg (t)X f = n(t) (16) where X¨ g is the ground acceleration, X f is the filter response, ωg (t) is the time-dependent predominant ground frequency, ζg (t) is the effective ground damping that is considered as a constant in present work, and e(t) is the amplitude of envelope function, time-dependent as well. Here, n(t) is a stationary Gaussian white noise process stated by Eq. (17): ⟨n(t)⟩ = 0, ⟨n (t1 ) n (t2 )⟩ = 2π S0 δ (t1 − t2 )

(17)

where S0 is the constant power spectral intensity of the noise, δ (t1 − t2 ) is the Dirac delta function, and the angular brackets stand for ensemble averages. Consider, for instance, that e(t) = 1, ωg and ζg being constants, then Eqs. (15)–(16) provide the original Kanai–Tajimi model. In this case, the power spectral intensity function is given by: ( )2 1 + 4ζg2 ω/ωg S0 (18) Sg = ( ( )2 )2 1 − ω/ωg + 4ζg2 ω2 /ωg2 To evaluate the ωg (t) and e(t) as a time-dependent amplitude envelope and ground frequency function, a recorded accelerogram is needed. Further, different earthquake records have different ground frequencies and amplitude envelope functions. The ζg is assumed as a constant, and the ground acceleration record is statically analyzed for evaluation of these functions, with some considerations to be followed as stated by several researches [17,42]. Firstly, the strong motion segment of synthetic accelerogram is selected in the way that it contains more than 95% energy of the original record. Secondly, the amplitude standard deviation of the ground acceleration is evaluated by using an appropriate time averaging procedure. Thirdly, the time averaged crossing rate is evaluated, and it is used for determination of time variation in the ground frequency. Finally, simple algebraic time functions are fitted to the data for the amplitude standard deviation and the mean zero-crossing rate. The amplitude envelope function, e(t), is assumed to be proportional to the amplitude standard deviation, and its coefficient is adjusted so that the mean energy of synthetic accelerograms is the same as that of the original record [17]. Based on these considerations, the amplitude envelop function and ground frequency function are proposed as [17,42]: e(t) = C0 σa (t) ωg (t) = π Fc (t)

(19) (20)

where C0 is a constant, σa (t) is a smooth standard deviation function, and Fc (t) is a smooth zero-crossing rate function. The σa (t) can be determined by using Eq. (21) [16,41]: σa (t) = c1 (σm − σw ) (t/ts )3 e−c2 t/ts + σw

(21)

where c1 and c2 are constants, σm is the maximum intensity of ground motion, and σw is the standard deviation of the weak motion, that corresponds to the tail of σa (t), and ts is the time duration of strong ground shaking, which is evaluated by the total energy of the smooth standard deviation. In addition, the zero-crossing rate function can be obtained by expression (22) [16], ( ) Fc (t) = F0 + F1 e−b1 t − e−b2 t (22) where F0 , F1 , b1 and b2 are constants to be obtained according to the conditions of earthquake records. The present work uses the smooth standard deviation function and zero-crossing rate function of Naghan [42] to generate synthetic accelerogram. For Naghan earthquake, consider c0 = 0, 074, ξg = 0, 35 and peak ground acceleration is 709,46 cm/s2 , Eqs. (20) and (21) become: ⎞ ⎛ 14, 71 )⎠ (23a) ωg = π ⎝−0, 998 + ( 1 1 + (t/0,558) 2,842

156

Y.S. Hsu / Mathematics and Computers in Simulation 169 (2020) 149–165

Fig. 3. Synthetic accelerogram profile Naghan.

2

σa = 49, 1 + 339, 84e−0,5(t−1,917/0,5134)

(23b)

The amplitude envelope function and ground frequency function are given by (19) and (20) respectively. Introducing these functions into (18), it is possible to evaluate the Kanai–Tajimi power density spectral for desired synthetic accelerogram. In the present work, S0 is assumed to be 1 cm2 /s3 . The Kanai–Tajimi power density spectral function is introduced into Shinozuka and Jan model, as stated by Eq. (24), and the synthetic accelerogram based on generalized nonstationary Kanai–Tajimi model is generated by the expression (24), N

as (t) =

f √ √ ∑ 2 s (ωk ) ∆ωcos (ωk t + φn )

(24)

k=1

where ωk is the frequency adopted in the present work with a gamma varies from 0 to 40 Hz, and φn is the random value of the phase angle, with the value lower or equal to unity. The synthetic acceleration is generated by the presented methodology and shown in Fig. 3. 4. Applications This section presents three examples to illustrate the effectiveness of the finite element approach adopted in present work and presents the investigation on the performance of the numerical model proposed. The first example deals with free vibration analysis of buried pipeline with finite length, resting on Winkler and Pasternak foundation. The reference results are obtained by the analytical solution proposed by Valsangkar and Pradhanang [52]. The natural frequencies are calculated for several values of spring stiffness and shear layer modulus. The second example presents the stress and strain analyses with the results obtained from dynamic elastic behavior of a buried pipeline. The error norm L 2 of displacement is determined by employing different values of the time interval. The last example presents a dynamic elastoplastic analysis of the buried pipeline, subjected to transversal ground motion. The error norm L 2 of stress is also determined by using different values of the time interval. The first example adopts mechanical properties employed by Valsangkar and Pradhanang [52], while the other two examples adopt the mechanical properties of API 5L grade X65. All examples adopt simply — supported Euler–Bernoulli beam as the boundary condition. 4.1. Free vibration analysis of pipe on Pasternak foundation This simple example aims to present the numerical results obtained from finite element approach and it has the purpose to validate the present numerical model by making a comparison between the numerical results and

157

Y.S. Hsu / Mathematics and Computers in Simulation 169 (2020) 149–165 Table 1 Value of frequency factor C for simply supported beam, literature results versus FEM3 results. Valsangkar and Pradhanang [52]

Present work FEM3

First mode α 1 1 × 102 1 × 104 1 × 106

First mode α

Sr 0

0,5

1

2,5

3,15 3,75 10,02 31,62

3,53 4,02 10,10 31,63

3,81 4,20 10,11 31,63

4,42 4,57 10,12 31,63

Second mode α 1 1 × 102 1 × 104 1 × 106

1 1 × 102 1 × 104 1 × 106

0

0,5

1

2,5

3,15 3,75 10,02 31,56

3,48 3,96 10,04 31,56

3,74 4,14 10,05 31,56

4,30 4,58 10,08 31,56

0

0,5

1

2,5

6,28 6,38 10,37 31,69

6,47 6,56 10,41 31,70

6,64 6,73 10,46 31,70

7,09 7,16 10,58 31,72

0

0,5

1

2,5

9,43 9,46 11,57 31,92

9,55 9,58 11,64 31,93

9,68 9,70 11,70 31,94

10,02 10,05 11,90 31,97

Second mode α

Sr 0

0,5

1

2,5

6,28 6,38 10,36 31,63

6,46 6,62 10,47 31,64

6,64 6,82 10,51 31,64

7,01 7,22 10,61 31,64

Third mode α

1 1 × 102 1 × 104 1 × 106

Sr

1 1 × 102 1 × 104 1 × 106

Sr

Third mode α

Sr 0

0,5

1

2,5

9,42 9,45 11,57 31,72

9,57 9,58 11,60 31,72

9,73 9,74 11,65 31,71

10,01 10,03 11,83 31,72

1 1 × 102 1 × 104 1 × 106

Sr

analytical results [52]. The structure model is represented schematically in Fig. 1(b). First three modes are calculated by using different values of spring stiffness and shear layer modulus. The mechanical and geometric parameters are: E = 210 GPa, ρ = 7830 kg/m3 , I = 7.3631 × 10−17 m4 , A = 2.3562 × 10−8 m2 and L = 1 m. ¯ 2 , with S¯ = G 2 L 2 /E I , and For the effect of a parametric study, consider for instance two parameters, Sr = S/π 4 α = k2 L /E finite element mesh is constituted by ten FEM3 elements. Table 1 shows the parametric ( I [52]. The)1/4 value C = m L 4 ω2 /E I , where ω is the natural frequency in rad/s. As it can be observed, the relative difference between analytical approach and finite element approach is less than 1% in almost all situations, solve the exception of the first mode, with α equal to 1 and Sr equal to 2,5. The numerical results show the competitiveness of the finite element approach adopted in the present work in the modeling of the buried pipeline. For α less than or equal to 1 × 102 , such as soft soil, the effects of the shear layer are significant. While for α higher than this value, the effects of the shear layer are insignificant. The improvement of the results can be done by making a h-refine in the finite element mesh. 4.2. Dynamic elastic analysis of buried pipeline surrounded by Winkler and Pasternak foundation This example investigates the effects of the time interval in the mechanical behavior of the buried pipeline, which adopts material API5L grade X65, with pipeline outer diameter 762 mm, wall thickness 17.5 mm, length 5 m and density 7850 kg/m3 . The stress–strain data is shown in Fig. 4a. The seismic wave propagates in the longitudinal direction, with components acting in the transversal direction of the pipeline. Two elastic foundation models are considered in this analysis, Winkler and Pasternak types. The Winkler spring stiffness is 36,24 MPa, and Pasternak shear layer modulus is 27,52 kN. Only Naghan synthetic accelerogram is considered in this analysis. The buried pipeline is modeled with finite element mesh, constituted by 20 FEM3 elements, with 51 nodes in the mesh. Moreover, a simply supported boundary condition is considered in this example, and HHT algorithm is employed with different values of ∆t, 1 × 10−3 s, 2 × 10−4 s, 1 × 10−4 s and 5 × 10−5 s. The error in norm L 2 is calculated for displacement and stress errors between a time discretization and another.

158

Y.S. Hsu / Mathematics and Computers in Simulation 169 (2020) 149–165

Fig. 4a. The stress–strain curve of material API5L — X65.

Fig. 4b. Schematic representation of finite element mesh.

Fig. 5 shows the displacement profile in node 26, Fig. 4b, located in the midpoint of the pipeline. The results are obtained by using different time intervals and different foundation types. The results related to Pasternak type foundation are depicted in Fig. 5(a), while the results related to Winkler type foundation are depicted in Fig. 5(b). It is possible to observe that the maximum displacement in both foundation types occurs at approximately 1,6 s, independent of time interval adopted. Moreover, the time interval influences the displacement magnitude, but not the frequency, as it can be observed in Fig. 5. In addition, the maximum displacement magnitude in Pasternak model is less than that produced by the Winkler foundation, because the Pasternak foundation has the shear layer in its model. The displacement profiles versus pipeline length are depicted in Fig. 6(a)–(b). The magnitude of maximum displacement and the instance that it occurs, are presented in each figure, for different time interval values and different foundation types. As it can be observed, the maximum value occurs approximately in 1,6 s, in almost all cases, which corresponds to the maximum variation in the PGA of Naghan synthetic accelerogram. Fig. 7(a)–(b) depicts the velocity versus displacement of node 26 when the time interval is ∆t = 1 × 10−5 s. As it can be observed, the displacement varies between −1.8 × 10−5 m and 2.1 × 10−5 m, while the velocity varies between −1.75 × 10−3 m/s and 1.7 × 10−3 m/s, for Pasternak type foundation. In the case of Winkler type foundation, the displacement varies between −2 × 10−5 m and 2.5 × 10−5 m, while the velocity varies between −2 × 10−3 m/s and 1.75 × 10−3 m/s. The norm L 2 is calculated for displacement error by making a comparison between the different values of time interval. The norm L 2 versus time is depicted in Fig. 8 for the different types of foundations. The denomination Error 1 in the legend symbolizes that the error is determined between ∆t = 1 × 10−3 s and 2 × 10−4 s. The Error 2 is the error determined between ∆t = 2 × 10−4 s and 1 × 10−4 s. Finally, the Error 3 is determined between ∆t = 1 × 10−4 and 5 × 10−5 s. The refinement in time discretization is made in this example, while the space discretization is maintained the same. It is possible to note in Fig. 8, that the maximum error in norm L 2 is remarkable in Error 1 during 1,5 to 2 s. However, with the time refinement, this error is reduced until it has almost the same value for whole time gamma, as it can be observed in the curve Error 3 in Fig. 8(a) and (b). Furthermore, the comparison between Fig. 8(a) and (b) leads to an observation that the foundation type has a significant influence on the error reduction. This fact can be observed by comparing the magnitude of the curve Error 2 and Error 3 of

Y.S. Hsu / Mathematics and Computers in Simulation 169 (2020) 149–165

159

Fig. 5. Displacement profile of node 26, Naghan synthetic accelerogram.

Fig. 6. Displacement profile along the pipeline length.

both figures. Finally, the mean error is also determined for each evaluated norm L 2 and the results are depicted in Fig. 9. Once again, it is possible to observe the remarkable reduction in the error between Error 1 and Error 2. The mean error reduces from 5 × 10−6 to 1 × 10−6 , approximately. However, the same situation is not repeated for Error 2 and Error 3. Nevertheless, based on the comparison of the results, it is possible to confirm that the time discretization refinement can reduce exponentially the error. 4.3. Nonlinear response of buried pipeline with different yield forces of the foundation In this example, an analysis of nonlinear response is carried out, with the same mechanical and geometrical parameters employed by the previous example. In addition, the finite element mesh is constituted by 20 elements. The limit reaction force is adopted as the foundation response force. The Naghan spectrum is employed in this example, however, with hypothetic PGA value equal to 4114.868 m/s2 to induce material hardening in the analysis.

160

Y.S. Hsu / Mathematics and Computers in Simulation 169 (2020) 149–165

Fig. 7. Velocity versus displacement in node 26, Naghan synthetic accelerogram, ∆t = 1 × 10−5 s.

Fig. 8. Error in norm L 2 of displacement calculated between the different values of time interval.

Fig. 9. Mean error of displacement for the different types of foundations.

Y.S. Hsu / Mathematics and Computers in Simulation 169 (2020) 149–165

161

Fig. 10. Stress distribution versus length of the pipeline, for different time intervals and different limit reaction forces.

By considering this value, the analyses are carried out by using different time intervals, and the error norm in L 2 is determined for stress. In the first place, the maximum von Mises stress distribution versus pipeline length for different ∆t and different values of foundation limit forces is presented in Fig. 10(a)–(b). It can be observed that the length of pipeline that undergoes the material hardening is significantly affected by the foundation limit force, as it is shown in Fig. 10(a). In the situation that involves Pasternak foundation and Fy = 80 kN/m, there is almost no presence of material hardening in the buried pipe. A similar situation is observed in Fig. 10(b). Moreover, with Fy = 50 kN/m, the material hardening is observed either in Pasternak foundation or in Winkler foundation, independent of the adopted time interval. In the case of ∆t = 1 × 10−3 s, the material hardening is not observed, even with F y = 50 kN/m or F y = 80 kN/m. To investigate the stress error, the norm L 2 is evaluated for different time discretization refinements. The results are depicted in Fig. 11(a)–(d). In the element domain, seven points of numerical integration are adopted, while in each integration point, the cross-section is divided into one layer and 36 circumferential parts. Therefore, the von Mises stress is evaluated in each point and in every circumferential part. The error norm L 2 of stress is then evaluated, as it can be observed in Fig. 11(a)–(d). The stress errors have higher value during material hardening when the time gamma varies from 1 s to 2 s, approximately. Which corresponds to the peak of Naghan accelerogram. In addition, Fig. 11(a)–(b) shows the results obtained by Pasternak and Winkler type foundations with Fy = 50 kN/m. It is possible to observe that the error 2 and error 3 increase significantly after the initiation of material hardening. However, a similar situation is not repeated in Figs. 11(c) and 11(d), which present the results obtained by Pasternak and Winkler type foundations,but with the value of Fy set as infinite. Finally, the mean error of stress is evaluated for the different foundation types and different limit forces considered in this example. The results are depicted in Fig. 12. As it can be observed, the errors are reduced significantly with the time discretization refinement. And the increase in the foundation limit force also contributes to this reduction. 5. Discussion and conclusions This work presents a finite element approach in the numerical modeling of dynamic elastoplastic response in a buried pipeline, subjected to seismic ground motion. In addition, an investigation is carried out concerning the effects of time discretization refinement in the numerical performance, which is measured by determining displacement and stress errors in norm L 2 . The numerical model is constituted by soil–structure interaction, which uses different foundations to simulate the behavior of surrounding soil. Further, the buried pipeline is discretized by Euler–Bernoulli 3 nodes beam — pipe element and it is subjected to the random ground motion generated by generalized nonstationary Kanai–Tajimi model. The governing equation is solved by using HHT time increment

162

Y.S. Hsu / Mathematics and Computers in Simulation 169 (2020) 149–165

Fig. 11. Norm L 2 of von Mises stress error for different foundation types and foundation limit forces.

Fig. 12. Mean error in norm L 2 of von Mises stress for different foundation types and different values of foundation limit forces.

Y.S. Hsu / Mathematics and Computers in Simulation 169 (2020) 149–165

163

method and Newton–Raphson algorithm is adopted to solve the nonlinear equation, that involves material isotropic hardening. The norm L 2 is evaluated for displacement and stress errors with a different time interval. In the free vibration analysis, the comparison is made between numerical results and reference results, and it demonstrates the competitiveness of the present model. In the dynamic elastic analysis, the buried pipeline is subjected to random ground motion, considering that the seismic wave travels in the longitudinal direction of the pipeline, but acts in the transversal direction. The time increment done by HHT algorithm uses a different value of time interval. For ∆t equal to 2 × 10−4 s, 1 × 10−4 s, and 5 × 10−5 s, the peak of displacement profiles versus pipeline length is similar. However, the same situation cannot be observed for velocity and acceleration. Moreover, the instance of peak value in the displacement, the velocity, and the acceleration is almost the same for different ∆t, and it occurs in the first moment of peak PGA. Therefore, it can be concluded that different values of ∆t only affect the magnitude of peak value, but not the instance that occurs the maximum value, neither oscillation frequency. The same conclusion is made for both Winkler type foundation and Pasternak type foundation. Moreover, the displacement error is reduced with time discretization refinement. A dynamic elastoplastic analysis is carried out with the different values of foundation limit forces in the Pasternak type and Winkler type foundations. In this case, the value of PGA is evaluated to ensure that it could lead the pipeline to material hardening. The maximum stress is determined along the pipeline for different ∆t. In the case of ∆t = 1 × 10−3 s, the pipeline does not suffer material hardening. However, with smaller ∆t, the material hardening is observed. Finally, with ∆t = 5 × 10−5 s, the foundation limit force affects the magnitude of the von Mises stress and contributes to reduce the error when the force is increased. Furthermore, the Pasternak type foundation contributes to reduce the numerical error more than Winkler type foundation for the same value of foundation limit force. Despite the HHT algorithm is characterized as an unconditionally stable time integration scheme, however, the error in norm L 2 of stress and displacement is reduced significantly only with ∆t ≤ 1 × 10−4 s. This observation leads to a suggestion to adopt ∆t ≤ 1 × 10−4 s in the dynamic elastoplastic analysis of buried pipeline subjected to random seismic ground motion. From the discussion presented above, some concluding remarks are listed below: • Different values of ∆t, adopted in the HHT, only affect the peak value in displacement and other dynamic variables, but not in the oscillation frequency of displacement versus time profile. • For the same value of limit force, Pasternak type foundation contributes more than Winkler type foundation in the reduction of stress and displacement errors in norm L 2 . • For dynamic elastoplastic analysis of buried pipeline, the HHT with ∆t ≤ 1 × 10−4 s is recommended, because it leads to very small stress and displacement errors in norm L 2 . References [1] A.M. Abbas, Critical seismic load inputs for simple inelastic structures, J. Sound Vib. 296 (2006) 949–967. [2] J.A. Abdalla, Y.M. Hag-Elhassan, Simulation of earthquake ground motion for generation of artificial accelerograms, in: Earthquake Resistant Engineering Structures, 2005. [3] M. Balkaya, M.O. Kaya, A. Saglamer, Analysis of the vibration of an elastic beam supported on elastic soil using the differential transform method, Arch. Appl. Mech. 79 (2009) 135–156. [4] K.J. Bathe, Finite Element Procedure, Prentice-Hall, 1996. [5] K. Bi, H. Hao, C. Li, H. Li, Stochastic seismic response analysis of buried onshore and offshore pipelines, Soil Dyn. Earthq. Eng. 94 (2017) 60–65. [6] Y. Chen, Simplified and refined earthquake analyses for buried pipes, Math. Comput. Modelling 21 (1995) 47–60. [7] C.N. Chen, Vibration of prismatic beam on an elastic foundation by the differential quadrature element method, Comput. Struct. 77 (2000) 1–9. [8] W.W. Chen, B.J. Shih, Y.C. Chen, J.H. Hung, H.H. Hwang, Seismic response of natural gas and water pipelines in the Ji-Ji earthquake, Soil Dyn. Earthq. Eng. 22 (2002) 1209–1214. [9] H. Chen, T. Zhong, G. Liu, J. Ren, Improvement of and parameter identification for the bimodal time varying modified Kanai – Tajimi Power Spectral Model, Shock Vib. (2017). [10] J.D. Colton, P.H.P. Chang, H.E. Lindberg, Measurement of dynamic soil-pipe axial interaction for full scale buried pipelines, Int. J. Soil Dyn. Earthq. Eng. 1 (1982) 183–188. [11] V. Corrado, B. D’Acunto, N. Fontana, M. Giugni, Estimation of dynamic strains in finite end-constrained pipes in seismic areas, Math. Comput. Modelling 49 (2009) 789–797. [12] T.K. Datta, Seismic response of buried pipelines: a state-of-the-art, Nucl. Eng. Des. 192 (1999) 271–284. [13] T.K. Datta, Seismic Analysis of Structures, John Wiley & Sons Pte Ltd, 2010.

164

Y.S. Hsu / Mathematics and Computers in Simulation 169 (2020) 149–165

[14] T. Dede, Y. Ayvaz, Y.I. Özdemir, Materially nonlinear free vibration analysis of reinforced concrete beam, J. Vib. Control 17 (2011) 1090–1098. [15] G. Dhatt, G. Touzot, G. Cantin, The finite element method displayed, John Wiley & Sons, 1984. [16] G.W. Ellis, A.S. Cakmak, Modelling Earthquake Ground Motions in Seismically Active Regions using Parametric Time Series Methods, Technical Report NCEER-87-0014, 1987. [17] F.C. Fan, G. Ahmadi, Nonstationary Kani-Tajimi models for El Centro 1940 and Mexico City 1985 earthquakes, Probab. Eng. Mech. 5 (1990) 171–181. [18] M. Fragiadakis, D. Vamvatsikos, M.G. Karlaftis, N.D. Lagaros, M. Papadrakakis, Seismic assessment of structures and lifelines, J. Sound Vib. 334 (2015) 29–56. [19] F.P. Gao, D.S. Jeng, H. Sekiguchi, Numerical study on the interaction between non-linear wave, buried pipeline and non-homogenous porous seabed, Comput. Geotech. 30 (2003) 535–547. [20] H.M. Hilber, T.J.R. Hughes, R.L. Taylor, Improved numerical dissipation for time integration algorithms in structural dynamics, Earthq. Eng. Struct. Dyn. 5 (1977) 283–292. [21] S. Imanzadeh, A. Denis, A. Marache, Simplified uncertainties analysis of continuous buried steel pipes on an elastic foundation in the presence of low stiffness zones, Comput. Geotech. 48 (2013) 62–71. [22] S. Joshi, A. Prashant, S.K. Jain, Analysis of buried pipelines subjected to reverse fault motion, Soil Dyn. Earthq. Eng. 31 (2011) 930–940. [23] A.E. Kampitsis, E.J. Sapountzakis, Dynamic analysis of beam-soil interaction systems with material and geometrical nonlinearities, Int. J. Non-Linear Mech. 90 (2017) 82–99. [24] K. Kanai, Semi-empirical formula for the seismic characteristics of the ground, Bull. Earthq. Res. Inst. Univ. Tokyo 35 (1957) 309–325. [25] K. Kanai, An empirical formula for the spectrum of strong earthquake motion, Bull. Eq. Res. Inst. 39 (1961) 85–95. [26] G. Lanzano, E. Salzano, F.S. De Magistris, G. Fabbrocino, Seismic vulnerability of natural gas pipelines, Reliab. Eng. Syst. Saf. 117 (2013) 73–80. [27] J.K. Lee, S. Jeong, J. Lee, Natural frequencies for flexural and torsional vibrations of beams on Pasternak foundation, Soils Found. 54 (2014) 1202–1211. [28] D.H. Lee, B.H. Kim, H. Lee, J.S. Kong, Seismic behavior of a buried gas pipeline under earthquake excitations, Eng. Struct. 31 (2009) 1011–1023. [29] Y. Li, G. Wang, Simulation and generation of spectrum compatible ground motions based on wavelet packet method, Soil Dyn. Earthq. Eng. 87 (2016) 44–51. [30] W. Liu, Z. Li, Z. Song, J. Li, Seismic reliability evaluation of gas supply networks based on the probability density evolution method, Struct. Saf. 70 (2018) 21–34. [31] W. Liu, Q. Sun, H. Miao, J. Li, Nonlinear stochastic seismic analysis of buried pipeline systems, Soil Dyn. Earthq. Eng. 74 (2015) 69–78. [32] J. Lubliner, Plasticity Theory, Pearson Education, 2006. [33] K. Matsubara, M. Hoshiya, Soil spring constants of buried pipelines for seismic design, J. Eng. Mech. 126 (2000) 76–83. [34] A. Moustafa, Critical earthquake load inputs for multi-degree-of-freedom inelastic structures, J. Sound Vib. 325 (2009) 532–544. [35] D. Nedjar, M. Hamane, M. Bensafi, S.M. Elachachi, D. Breysse, Seismic response analysis of pipes by a probabilistic approach, Soil Dyn. Earthq. Eng. 27 (2007) 111–115. [36] P.M. O’Leary, S.K. Datta, Dynamics of buried pipelines, Int. J. Soil Dyn. Earthq. Eng. 4 (1985) 151–159. [37] Y.I. Özdemir, Development of a higher order finite element on a Winkler foundation, Finite Elem. Anal. Des. 48 (2012) 1400–1408. [38] Y.I. Özdemir, Using fourth order element for free vibration parametric analysis of thick plates resting on elastic foundation, Struct. Eng. Mech. 65 (2018) 213–222. [39] J.G. Palacios, A. Samartin, V. Negro, A nonlinear analysis of laying a floating pipeline on the seabed, Eng. Struct. 31 (2009) 1120–1131. [40] Y. Peng, J. Chen, J. Li, Nonlinear response of structures subjected to stochastic excitations via probability density evolution method, in: Advances in Structural Engineering, vol. 17, 2014, pp. 801–816. [41] N.W. Polhemus, A.S. Cakmak, Simulation of earthquake ground motion using autoregressive moving average (ARMA) model, Earthq. Eng. Struct. Dyn. 9 (1981) 343–354. [42] F.R. Rofooei, A. Mobarake, G. Ahmadi, Generation of artificial earthquake records with a nonstationary Kanai-Tajimi model, Eng. Struct. 23 (2001) 827–837. [43] E.J. Sapountzakis, Kampitsis, Nonlinear response of shear deformable beams on tensionless nonlinear viscoelastic foundation under moving loads, J. Sound Vib. 330 (2011) 5410–5426. [44] E.J. Sapountzakis, A.E. Kampitsis, Inelastic analysis of beams on two-parameter tensionless elastoplastic foundation, Eng. Struct. 48 (2013) 389–401. [45] K. Sarkar, V.K. Gupta, R.C. George, Wavelet-based generation of spatially correlated accelerograms, Soil Dyn. Earthq. Eng. 87 (2016) 116–124. [46] P. Shi, Seismic wave propagation effects on buried segmented pipelines, Soil Dyn. Earthq. Eng. 72 (2015) 89–98. [47] D. Soares, Dynamic analysis of elastoplastic models considering combined formulations of the time-domain boundary element method, Eng. Anal. Bound. Elem. 55 (2015) 28–39. [48] H.O. Soliman, T.K. Datta, Response of overground pipelines to random ground motion, Eng. Struct. 18 (1996) 537–545. [49] H. Tajimi, A statistical method of determining the maximum response of a building structure during an earthquake, in: Proc. 2d World of Earthquake Eng., vol. 2, 1960, pp. 781–798. [50] D. Thambiratnam, Y. Zhuge, Free vibration analysis of beams on elastic foundation, Comput. Struct. (1996) 971–980.

Y.S. Hsu / Mathematics and Computers in Simulation 169 (2020) 149–165

165

[51] E. Uckan, B. Akbas, J. Shen, W. Rou, F. Paolacci, M. O’Rourke, A simplified analysis model for determining the seismic response of buried steel pipes at strike-slip fault crossings, Soil Dyn. Earthq. Eng. 75 (2015) 55–65. [52] A.J. Valsangkar, R. Pradhanang, Vibrations of beam–columns on two-parameter elastic foundations, Earthq. Eng. Struct. Dyn. 16 (1988) 217–225. [53] P. Vazouras, P. Dakoulas, S.A. Karamanos, Pipe-soil interaction and pipeline performance under strike-slip fault movements, Soil Dyn. Earthq. Eng. 72 (2015) 48–65. [54] Y. Yan, Response of pipeline structure subjected to ground motion excitation, Eng. Struct. 19 (1997) 679–684. [55] T. Yokoyama, Vibrations of timoshenko beam–columns on two-parameters elastic foundations, Earthq. Eng. Struct. Dyn. 20 (1991) 355–370. [56] O.C. Zienkiewicz, R.L. Taylor, The Finite Element Method. Volume 1: The Basis, Butterworth Heinemann, 2000.