# A Finite Element Approach for the Elastic-Plastic Behavior of a Steel Pipe Used to Transport Natural Gas

Authors: Miltiades C. Elliotis
Journal: Conference Papers in Energy
 (2013)

---

## Abstract

A finite element technique, with two-dimensional isoparametric elements, is developed, for the analysis of a steel pipe, with a circular cross-section. The pipe is installed above the ground and is used to transport natural gas at very high pressures. The material of the pipe is assumed to obey a bilinear elastic-plastic model. Double symmetry is considered when setting up the mathematical model problem and creating the finite element mesh. Step increments or decrements on the internal pressure are applied (cyclic loading). Yielding is detected by the Von Mises yield criterion. Also, aJ2 flow rule is employed to handle the plastic strain component. A four-point Gauss-Legendre quadrature is used to numerically perform all necessary integrations. Finally, the so-called “shakedown phenomenon” is studied when cyclic loading is applied after the commencement of plastic deformations on the pipe. Numerical results obtained by the method compare favorably with the analytic solution.


---

## Body

## 1. Introduction

The efficient and effective movement of natural gas from producing regions to consumption regions requires an extensive and elaborate transportation system. In many instances, natural gas produced from a particular well will have to travel a great distance to reach its point of use. The transportation system for natural gas consists of a complex network of pipelines, designed to quickly and efficiently transport natural gas from its origin to areas of high gas demand.In the USA and in Europe there are three major types of pipelines along the transportation route: the gathering system, the interstate pipeline system, and the distribution system. The gathering system consists of low pressure, small diameter pipelines that transport raw natural gas from the wellhead to the processing plant. Interstate highway system transports processed natural gas from processing plants to those areas with high natural gas requirements and can carry natural gas across state or country boundaries at pressures varying between 15 and 200 MPa. This reduces the volume of the natural gas being transported (by up to 600 times), as well as propelling of the gas through the pipeline. The distribution system consists of low pressure pipes which distribute the transported natural gas within any region with high natural gas requirements.Mainline transmission pipes, the principle pipelines in a given system, are usually between 200 and 1200 mm in diameter and with a wall thickness between 10 and 100 mm. The original part of a pipe, in which propelling of the gas is performed and pressures up to 200 MPa are often applied, is usually around 1000 mm in its diameter as shown in Figure1.Figure 1
An example of a steel pipe, used to carry natural gas, with a very large external diameter.The actual pipeline itself, commonly called “line pipe,” consists of a strong carbon steel material, engineered to meet high standards. Thus, there is usually a need to perform a high quality design of such a pipe. The engineering problem is to find the appropriate external diameter and wall thickness. Usually, the diameter of the pipe is chosen according to the needs in natural gas and the desired rate of supply. Therefore, it remains to estimate the wall thickness by considering the properties of the material used and the internal pressure of the gas.For the parts of the pipeline installed above ground we are interested in solving the problem of a steel pipe with internal pressure thousands of times greater than the external atmospheric pressure. Solutions for the elastic behavior (i.e., before the initiation of yielding) were given by many theoreticians like Timoshenko and Goodier [1], Beer and Johnston [2], and Fung [3], to mention some of them. But of great interest is the elastic-plastic case. The engineers who design such pipes consider that plastic deformations appear first on the internal surface and as the pressure increases yielding of the material progresses until it reaches the outer surface. They choose dimensions and material properties which allow either partial yielding of the cross-section or no yielding at all.The elastic-plastic problem is an old problem. It was first analytically approached by Prantdl (1925) and Reuss (1930). First they assumed that the changes of strain are to be divisible into static and plastic components. The elastic strain increment was expressed in terms of the deviatoric and hydrostatic stress components while the plastic strain component was assumed to be proportional to the stress gradient of the plastic potential. Also, by using the Von Mises yield criterion and simple equilibrium of stresses (Figure2), one obtains a hyperbolic system of three semilinear differential equations with first-order partial derivatives, in which the principal unknown is the radius of the elastic-plastic boundary. These equations were further simplified by Hencky who proposed the so-called stress-deformation law.Figure 2
Hoop stresses on the cross-section of a pipe in equilibrium with internal pressure (radial and longitudinal stresses also exist).But the resulting system of differential equations is still much complicated and cannot be solved theoretically to give a closed form of an exact solution. They can be solved numerically by replacing the partial derivatives with finite differences and obtaining a system of linear algebraic equations. Unfortunately, the numerical result is of poor accuracy. So, instead of numerically solving the original system of differential equations, researchers solved the system of equations for several values of the radius of the elastic plastic boundary and created a large number of graphs of the normalized radial stressσr/k in the wall of the tube (k being the yield stress) versus the normalized distance r/a from the center (a being the internal radius). But the necessity to apply a numerical scheme to estimate the ultimate loads, as well as the stresses, strains, and displacements at any position and with adequate accuracy, was still there. It was the finite element method (FEM) which satisfied this need. This powerful method was developed in parallel with the mathematical theory of plasticity mainly during the second half of the 20th century.The tremendous progress in the field of computer engineering alongside with the development of the FEM by Argyris [4], Zienkiewicz and Taylor [5], and most recently by Babuška et al. [6] and Bathe [7] and others (to mention some of them) gave to the engineers a powerful tool to solve many elastic-plastic problems by running codes on very fast computing machines. At the same time a more systematic approach and formulation of the mathematical theory of plasticity was developed by Hill [8], Prager [9] and others. Later on, in the 70s and 80s, the works of Bath and Wilson [10], of Hinton and Owen [11], and of others inaugurate a new era in the treatment of engineering problems concerning materials with elastic-plastic behavior.In this work the problem of a thick pipe carrying natural gas is solved. It is subjected to a gradually increasing internal pressure, followed by a decrease in pressure. It is also subjected to a repeated increase and decrease in the value of internal pressure (cyclic loading and unloading), with plane strain conditions being assumed in any cross section along the axis of the pipe. In Section2, a brief presentation of the mathematical theory of plasticity is made which includes a presentation of the Von Mises yield criterion, the mathematical models of strain hardening behavior, and the elastic-plastic stress/strain relation. In Section 3 a general presentation of the FEM is made, and in Section 4, the model problem is presented and numerical results are given together with the results of the analytical solution. This paper ends with conclusions and references which are given in Section 5 and References section, respectively.


## 2. The Mathematical Theory of Plasticity

The object of the mathematical theory of plasticity is to provide a theoretical description of the relationship between stress and strain for a material which exhibits an elastic-plastic response. Plastic behavior is characterized by an irreversible straining which is independent of time and which can only be sustained once a certain level of stress has been reached. In this section we outline the basic assumptions and associated theoretical expressions for a general continuum. An extensive treatment of this theory can be found in the works of Hill [8], Prager [9], and others.In order to formulate a theory which models elastoplastic material deformation, three requirements must be met:(i)
An explicit relationship between stress and strain must be formulated to describe material behavior under elastic conditions, that is, before the onset of plastic deformation.(ii)
A yield criterion, indicating the stress level at which plastic flow commences, must be postulated.(iii)
A relationship between stress and strain must be developed for postyield behavior, that is, at the stage where the deformation consists of both elastic and plastic components.Before the onset of plastic yielding the relationship between stress and strain is given by the standard linear elastic expression(1)σij=Cijklεkl,
where σij and εkl are the stress and strain components respectively, and Cijkl is the tensor of elastic constants which, for an isotropic material, has the explicit form
(2)Cijkl=λδijδkl+μ(δikδjl+δilδjk),
where λ and μ are the Lamé constants and δij is the so-called Kronecker delta, a useful symbol in mechanics defined by
(3)δij={1,ifi=j,0,ifi≠j.
Equation (1) will suffice in the ideal case of an elastic and homogeneous material. But, in general, this is not the case for the metals and especially for carbon steel and other steel alloys which exhibit both elastic and elastic-plastic behaviors under certain combinations of loadings. The steel pipes used for the transportation of natural gas are not an exception.


### 2.1. The Yield Criterion

The yield criterion determines the stress level at which plastic deformation begins and can be written in the general form(4)f(σij)=Y(k),
where f is a function of stresses σij and Y expresses the yield surface of a material (Figure 3). It has the meaning that a state of stress which gives a point inside the surface represents elastic behavior. In (4) parameter k is a material parameter to be determined experimentally. It is the so-called yield stress and usually is denoted as σy.Figure 3
The mathematical representation of the Von Mises and Tresca yield surfaces, suitable for three-dimensional stress states.On physical grounds, any yield criterion should be independent of the orientation of the coordinate system employed, and therefore it should be a function of the three stress invariants only(5)J1=σii,J2=12σijσij,J3=13σijσikσki.
Experimental observation, notably by Bridgman [12], indicate that plastic deformation of metals is essentially independent of hydrostatic pressure. Consequently, the yield function can only be of the form
(6)f(J2′,J3′)=Y(k),
where J2′ and J3′ are the second and third invariants
(7)J2′=12σij′σij′,J3′=13σij′σik′σki′,
where
(8)σij′=σij-13δijσkk.
Most of the various yield criteria that have been suggested for metals are now only of historic interest, since they conflict with experimental predictions. The two most simple criteria which do not have this fault are the Tresca criterion and the Von Mises criterion.In this work we adopted the Von Mises yield criterion because it has been proved, so far, to work very well for metals. It was first proposed by Von Mises [13] in 1913. According to this criterion yielding occurs when J2′ reaches a critical value as
(9)J2′=Y(k),
where, again, k is the yield stress. The second deviatoric stress invariant, J2′, can be explicitly written as
(10)J2′=12σij′σij′=16[(σ1-σ2)2+(σ2-σ3)2+(σ3-σ1)]=12[(σ11′)2+(σ22′)2+(σ33′)2]+σ122+σ232+σ132.
Yield criterion (9) may be further written as
(11)σeff=3J2′=3Y(k).
Quantity σeff is termed effective or equivalent stress.Figure3 shows the geometrical interpretation of the Von Mises yield surface to be a circular cylinder whose projection onto the π-plane is a circle of radius k√2. The two-dimensional plot of this yield surface has an elliptic form (Figure 4). Now, a physical meaning of the constant k can be obtained by considering the yielding of materials under simple stress states. So, the case of uniaxial tension (σ2=σ3=0) requires that value k√3 is the uniaxial yield stress.Figure 4
The mathematical representation of the Von Mises and Tresca yield surfaces, suitable for two-dimensional stress states.


### 2.2. Mathematical Models of Strain Hardening Behavior

After initial yielding, the stress level at which further plastic deformation occurs may be dependent on the current degree of plastic straining. Such a phenomenon is termed work hardening or strain hardening. Thus the yield surface will vary at each stage of the plastic deformation, with the subsequent yield surfaces being dependent on the plastic strains in some way. Some alternative models, which describe strain hardening are illustrated in Figures5 and 6.Figure 5
Mathematical model for the representation of isotropic strain hardening.Figure 6
Mathematical model for the representation of kinematic strain hardening.In the case where the yield surfaces are a uniform expansion of the original yield curve, without translation, as shown in Figure5, the strain hardening model is said to be isotropic. On the other hand if the subsequent yield surfaces preserve their shape and orientation but translate in the stress space as a rigid body, as shown in Figure 6, kinematic hardening is said to take place. Such a hardening model gives rise to the experimentally observed Bauschinger effect on cyclic loading.The progressive development of the yield surface can be defined by relating the yield stressY, to the plastic deformation by means of the hardening parameter k. This can be done in two ways. Firstly, the degree of work hardening can be postulated to be a function of the total plastic work, Wp, only. Then,
(12)k=Wp=∫Ωσijd(εij)p,
where (dεij)p are the plastic components of strain occurring during a strain increment. Alternatively, k can be related to a measure of the total plastic deformation termed the effective or equivalent plastic strain which is defined incrementally as
(13)dεp(eff)=23(dεij)p(dεij)p.
For situations where the assumption that yielding is independent of any hydrostatic stress is valid, we have (dεij′)p=(dεij)p. Consequently, (13) can be rewritten as
(14)dεp(eff)=23(dεij′)p(dεij′)p.
Then, the hardening parameter k is defined as
(15)k=εp(eff),
where εp(eff) is the result of integrating dεp(eff) over the strain path. This behavior is termed as strain hardening.By considering an isotropic hardening model, all the points on the yield surfacef=Y(k) represent plastic states, while elastic behavior is characterized by f<Y(k). At plastic state, f=Y(k), the incremental change in the yield function due to an incremental stress change is
(16)df=∂f∂σijdσij.
Then, if df<0, elastic unloading occurs (stress point is inside the yield surface);df=0, neutral loading takes place (stress point remains on the yield surface);df>0, plastic loading occurs (stress point remains on the expanding yield surface).It can also be shown [8] that, for a stable material, the initial and all subsequent yield surfaces must be convex.According to Drucker’s definition [14], strain hardening in a material, which already exhibits plastic deformations, occurs when the following important inequality is true:
(17)dσijdεij>0,
upon loading, and when
(18)dσij(dεij)p≥0,
a completed cycle of loading is indicated.


### 2.3. Elastic-Plastic Stress/Strain Relation

After initial yielding the material behavior will be partly elastic and partly plastic. During any increment of stress, the changes of strain are assumed to be divisible into elastic and plastic components, so that(19)dεij=(dεij)e+(dεij)p.
The elastic strain increment is related to the stress increment by (1), and by decomposing the stress terms into their deviatoric and hydrostatic components the expression for d(εij)e is as follows:
(20)(dεij)e=dσij′2μ+1-2vEδij(dσkk),
where E and v are, respectively, the elastic modulus and Poisson’s ratio of the material.In order to derive the relationship between the plastic strain component and the stress increment, a further assumption on the material behavior must be made. In particular, it will be assumed that the plastic strain increment is proportional to the stress gradient of a quantity termed the plastic potentialQ, so that
(21)(dεij)p=dλ∂Q∂σij,
where dλ is a proportionality constant termed the plastic multiplier. A theoretical basis for this assumption was developed by Hill [8]. Equation (21) is termed the flow rule. It governs the plastic flow after yielding. The potential Q must be a function of J2′ and J3′ but as yet cannot be determined in its most general form. However, the relation f=Q has a special significance in the mathematical theory of plasticity, since for this case certain variational principles and uniqueness theorems can be formulated. The identity f≡Q is a valid one since it has been postulated that both are functions of J2′ and J3′ and such an assumption gives rise to an associated theory of plasticity. In this case (21) becomes
(22)(dεij)p=dλ∂f∂σij
and is termed the normality rule. The derivative ∂f/∂σij is a vector directed normal to the yield surface at the stress point under consideration, as shown in Figure 7. It is seen that the components of the plastic strain increment are required to combine vectorially in an n-dimensional space to give a vector which is normal to the yield surface. For the particular case of f=J2′, we have
(23)∂f∂σij=∂J2′∂σij=σij′.
Then, (22) becomes
(24)(dεij)p=dλσij′
which is the so-called set of the Prandtl-Reuss equations [8]. Experimental observations indicate that the normality condition is an acceptable assumption for metals and especially for carbon steel with which pipelines transporting natural gas are made of. Now, on use of (19), (20), and (22), the complete incremental relationship between stresses and strains for elastic-plastic deformation is found to be
(25)dεij=dσij′2μ+1-2vEδij(dσkk)+dλ∂f∂σij.Figure 7
Geometrical representation of the normality rule of plasticity.The hardening lawY=Y(k) could just be easily expressed in terms of the effective stress σeff (since it is proportional to √J2′) to give, for the strain hardening hypothesis (15),
(26)σeff=H(εp(eff)),
or differentiating,
(27)dσeffdεp(eff)=H′(εp(eff)).
In (27) the derivative H′ is the hardening function and it is the slope of the uniaxial stress/plastic strain curve for which σeff=σ11 and dεp(eff)=dεp. Thus, it can be further expressed as follows:
(28)H′(εp(eff))=dσ11dεp=dσ11dε-dεe=⋯=ET1-(ET/E),
where ET is the elastoplastic tangent modulus (Figure 8).Figure 8
Elastic-plastic strain hardening behavior for the uniaxial case.The theoretical expressions developed so far can now be converted to matrix form. The yield function, first defined in (4), can be rewritten as
(29)f(σ)=Y(k),
where σ is the stress vector and k is the hardening parameter which governs the expansion of the yield surface. In particular, from (12), we have dk=σTdεp for the work hardening hypothesis. Rearranging (29), we get
(30)F(σ,k)=f(σ)-Y(k)=0.
By differentiating F(σ,k), we have
(31)dF=∂F∂σdσ+∂F∂kdk=0,
which is expressed as
(32)αTdσ-Adλ=0,
with the definitions
(33)αT=∂F∂σ=[∂F∂σ11,∂F∂σ22,∂F∂σ33,∂F∂σ23,∂F∂σ31,∂F∂σ12],A=-1dλ∂F∂kdk.
Zienkiewicz et al. [15] give an elegant proof of the following result:
(34)A=H′.
Thus, parameter A is obtained to be the local slope of the uniaxial stress/plastic strain curve and can be determined experimentally from (28).


## 2.1. The Yield Criterion

The yield criterion determines the stress level at which plastic deformation begins and can be written in the general form(4)f(σij)=Y(k),
where f is a function of stresses σij and Y expresses the yield surface of a material (Figure 3). It has the meaning that a state of stress which gives a point inside the surface represents elastic behavior. In (4) parameter k is a material parameter to be determined experimentally. It is the so-called yield stress and usually is denoted as σy.Figure 3
The mathematical representation of the Von Mises and Tresca yield surfaces, suitable for three-dimensional stress states.On physical grounds, any yield criterion should be independent of the orientation of the coordinate system employed, and therefore it should be a function of the three stress invariants only(5)J1=σii,J2=12σijσij,J3=13σijσikσki.
Experimental observation, notably by Bridgman [12], indicate that plastic deformation of metals is essentially independent of hydrostatic pressure. Consequently, the yield function can only be of the form
(6)f(J2′,J3′)=Y(k),
where J2′ and J3′ are the second and third invariants
(7)J2′=12σij′σij′,J3′=13σij′σik′σki′,
where
(8)σij′=σij-13δijσkk.
Most of the various yield criteria that have been suggested for metals are now only of historic interest, since they conflict with experimental predictions. The two most simple criteria which do not have this fault are the Tresca criterion and the Von Mises criterion.In this work we adopted the Von Mises yield criterion because it has been proved, so far, to work very well for metals. It was first proposed by Von Mises [13] in 1913. According to this criterion yielding occurs when J2′ reaches a critical value as
(9)J2′=Y(k),
where, again, k is the yield stress. The second deviatoric stress invariant, J2′, can be explicitly written as
(10)J2′=12σij′σij′=16[(σ1-σ2)2+(σ2-σ3)2+(σ3-σ1)]=12[(σ11′)2+(σ22′)2+(σ33′)2]+σ122+σ232+σ132.
Yield criterion (9) may be further written as
(11)σeff=3J2′=3Y(k).
Quantity σeff is termed effective or equivalent stress.Figure3 shows the geometrical interpretation of the Von Mises yield surface to be a circular cylinder whose projection onto the π-plane is a circle of radius k√2. The two-dimensional plot of this yield surface has an elliptic form (Figure 4). Now, a physical meaning of the constant k can be obtained by considering the yielding of materials under simple stress states. So, the case of uniaxial tension (σ2=σ3=0) requires that value k√3 is the uniaxial yield stress.Figure 4
The mathematical representation of the Von Mises and Tresca yield surfaces, suitable for two-dimensional stress states.


## 2.2. Mathematical Models of Strain Hardening Behavior

After initial yielding, the stress level at which further plastic deformation occurs may be dependent on the current degree of plastic straining. Such a phenomenon is termed work hardening or strain hardening. Thus the yield surface will vary at each stage of the plastic deformation, with the subsequent yield surfaces being dependent on the plastic strains in some way. Some alternative models, which describe strain hardening are illustrated in Figures5 and 6.Figure 5
Mathematical model for the representation of isotropic strain hardening.Figure 6
Mathematical model for the representation of kinematic strain hardening.In the case where the yield surfaces are a uniform expansion of the original yield curve, without translation, as shown in Figure5, the strain hardening model is said to be isotropic. On the other hand if the subsequent yield surfaces preserve their shape and orientation but translate in the stress space as a rigid body, as shown in Figure 6, kinematic hardening is said to take place. Such a hardening model gives rise to the experimentally observed Bauschinger effect on cyclic loading.The progressive development of the yield surface can be defined by relating the yield stressY, to the plastic deformation by means of the hardening parameter k. This can be done in two ways. Firstly, the degree of work hardening can be postulated to be a function of the total plastic work, Wp, only. Then,
(12)k=Wp=∫Ωσijd(εij)p,
where (dεij)p are the plastic components of strain occurring during a strain increment. Alternatively, k can be related to a measure of the total plastic deformation termed the effective or equivalent plastic strain which is defined incrementally as
(13)dεp(eff)=23(dεij)p(dεij)p.
For situations where the assumption that yielding is independent of any hydrostatic stress is valid, we have (dεij′)p=(dεij)p. Consequently, (13) can be rewritten as
(14)dεp(eff)=23(dεij′)p(dεij′)p.
Then, the hardening parameter k is defined as
(15)k=εp(eff),
where εp(eff) is the result of integrating dεp(eff) over the strain path. This behavior is termed as strain hardening.By considering an isotropic hardening model, all the points on the yield surfacef=Y(k) represent plastic states, while elastic behavior is characterized by f<Y(k). At plastic state, f=Y(k), the incremental change in the yield function due to an incremental stress change is
(16)df=∂f∂σijdσij.
Then, if df<0, elastic unloading occurs (stress point is inside the yield surface);df=0, neutral loading takes place (stress point remains on the yield surface);df>0, plastic loading occurs (stress point remains on the expanding yield surface).It can also be shown [8] that, for a stable material, the initial and all subsequent yield surfaces must be convex.According to Drucker’s definition [14], strain hardening in a material, which already exhibits plastic deformations, occurs when the following important inequality is true:
(17)dσijdεij>0,
upon loading, and when
(18)dσij(dεij)p≥0,
a completed cycle of loading is indicated.


## 2.3. Elastic-Plastic Stress/Strain Relation

After initial yielding the material behavior will be partly elastic and partly plastic. During any increment of stress, the changes of strain are assumed to be divisible into elastic and plastic components, so that(19)dεij=(dεij)e+(dεij)p.
The elastic strain increment is related to the stress increment by (1), and by decomposing the stress terms into their deviatoric and hydrostatic components the expression for d(εij)e is as follows:
(20)(dεij)e=dσij′2μ+1-2vEδij(dσkk),
where E and v are, respectively, the elastic modulus and Poisson’s ratio of the material.In order to derive the relationship between the plastic strain component and the stress increment, a further assumption on the material behavior must be made. In particular, it will be assumed that the plastic strain increment is proportional to the stress gradient of a quantity termed the plastic potentialQ, so that
(21)(dεij)p=dλ∂Q∂σij,
where dλ is a proportionality constant termed the plastic multiplier. A theoretical basis for this assumption was developed by Hill [8]. Equation (21) is termed the flow rule. It governs the plastic flow after yielding. The potential Q must be a function of J2′ and J3′ but as yet cannot be determined in its most general form. However, the relation f=Q has a special significance in the mathematical theory of plasticity, since for this case certain variational principles and uniqueness theorems can be formulated. The identity f≡Q is a valid one since it has been postulated that both are functions of J2′ and J3′ and such an assumption gives rise to an associated theory of plasticity. In this case (21) becomes
(22)(dεij)p=dλ∂f∂σij
and is termed the normality rule. The derivative ∂f/∂σij is a vector directed normal to the yield surface at the stress point under consideration, as shown in Figure 7. It is seen that the components of the plastic strain increment are required to combine vectorially in an n-dimensional space to give a vector which is normal to the yield surface. For the particular case of f=J2′, we have
(23)∂f∂σij=∂J2′∂σij=σij′.
Then, (22) becomes
(24)(dεij)p=dλσij′
which is the so-called set of the Prandtl-Reuss equations [8]. Experimental observations indicate that the normality condition is an acceptable assumption for metals and especially for carbon steel with which pipelines transporting natural gas are made of. Now, on use of (19), (20), and (22), the complete incremental relationship between stresses and strains for elastic-plastic deformation is found to be
(25)dεij=dσij′2μ+1-2vEδij(dσkk)+dλ∂f∂σij.Figure 7
Geometrical representation of the normality rule of plasticity.The hardening lawY=Y(k) could just be easily expressed in terms of the effective stress σeff (since it is proportional to √J2′) to give, for the strain hardening hypothesis (15),
(26)σeff=H(εp(eff)),
or differentiating,
(27)dσeffdεp(eff)=H′(εp(eff)).
In (27) the derivative H′ is the hardening function and it is the slope of the uniaxial stress/plastic strain curve for which σeff=σ11 and dεp(eff)=dεp. Thus, it can be further expressed as follows:
(28)H′(εp(eff))=dσ11dεp=dσ11dε-dεe=⋯=ET1-(ET/E),
where ET is the elastoplastic tangent modulus (Figure 8).Figure 8
Elastic-plastic strain hardening behavior for the uniaxial case.The theoretical expressions developed so far can now be converted to matrix form. The yield function, first defined in (4), can be rewritten as
(29)f(σ)=Y(k),
where σ is the stress vector and k is the hardening parameter which governs the expansion of the yield surface. In particular, from (12), we have dk=σTdεp for the work hardening hypothesis. Rearranging (29), we get
(30)F(σ,k)=f(σ)-Y(k)=0.
By differentiating F(σ,k), we have
(31)dF=∂F∂σdσ+∂F∂kdk=0,
which is expressed as
(32)αTdσ-Adλ=0,
with the definitions
(33)αT=∂F∂σ=[∂F∂σ11,∂F∂σ22,∂F∂σ33,∂F∂σ23,∂F∂σ31,∂F∂σ12],A=-1dλ∂F∂kdk.
Zienkiewicz et al. [15] give an elegant proof of the following result:
(34)A=H′.
Thus, parameter A is obtained to be the local slope of the uniaxial stress/plastic strain curve and can be determined experimentally from (28).


## 3. The Finite Element Method

The finite element method (FEM) is a numerical technique for finding approximate solutions to partial differential equations and their systems, as well as integral equations. The basic concept of this method is the discretization of the domain of the problem into finite elements. There are several tricks about how to create an efficient finite element mesh in order to have results as accurate as possible and spend a computational time as small as possible.Several types of elements have been developed of which the isoparametric elements were found to be the most suitable for the work presented in this paper. So, in a finite element representation with isoparametric elements, the displacements and strains and their virtual counterparts may be expressed by the relationships(35)u=∑i=1nNidi=N·d,δu*=∑i=1nNiδdi*=N·δd*,ε=∑i=1nBidi=B·d,δε*=∑i=1nBiδdi*=B·δd*,
where for node i, vector d is the vector of nodal variables, δd* is the vector of virtual nodal variables, N is the matrix of  global shape functions, and B is the global strain-displacement matrix. For a 4-node isoparametric element shape function, Ni is expressed in the local system Oξζ of the element as Ni(ξ,ζ)=(1/4)·(1+ξξi)·(1+ζζi). The total number of nodes on the whole mesh is n.The basic expressions required for solution can be obtained by use of the principle of virtual work. So, the virtual work expression for the whole domain, represented by the finite element mesh, is as follows:(36)∑i=1n[δdi*]T{∫Ω[Bi]TσdV-∫Ω[Ni]TbdV-∫∂Ω[Ni]TtdS}=0,⟶∫Ω(δd*)T(BTσ-NTb)dV-(δd*)f=0,
where t are the boundary tractions which are also denoted by f (external applied forces) and b are the distributed internal loads per unit volume. Since (36) must be true for any arbitrary δd* value, we have
(37)∫ΩBTσdV-∫ΩNTbdV-f=0.
Here, for the plane strain problem, the discretized elemental volume dV(e) represents an element area which on the local coordinate system of the element is given as
(38)dV(e)=|J(e)|dξdζ,
where J(e) is the Jacobian matrix of the transformation from the local system of the element to the global Cartesian system of the problem’s domain.For the solution of nonlinear problems, with elastic-plastic behavior and strain hardening, (38) will not generally be satisfied at any stage of the computation and
(39)ψ=∫ΩBTσdV-(∫ΩNTbdV+f)≠0,
where vector ψ is the residual force vector of the nodal forces. For an elastic-plastic situation the material stiffness is continually varying and instantaneously the incremental stress/strain relationship is given by dσ=Depdε. For the purpose of evaluating the global tangential stiffness matrix KT at any stage, the incremental form of (4) must be employed. Thus, within an increment of load, we have
(40)Δψ=∫ΩBT(Δσ)dV-(Δf+∫ΩNT(Δb)dV)=KTd-(Δf+∫ΩNT(Δb)dV),
where
(41)KT=∫ΩBTDepBdV.
The procedure for evaluating Δψ in (40), after the initiation of yielding, involves the implementation of suitable iteration methods (refined process) in order to reduce a stress point to the current yield surface, at positions where loading in the plastic region occurs. Therefore, vector ψ of the nodal forces must be equivalent to the stress field satisfying elastoplastic conditions. During the application of an increment of load, an element, or part of an element, may yield. All stress and strain quantities are monitored at each Gaussian integration point, and therefore it is possible to determine whether or not plastic deformation has occurred at such points. Consequently an element can behave partly elastically and partly elastoplastically if some, but not all, Gauss points indicate plastic yielding. For any load increment it is necessary to determine what proportion is elastic and which part produces plastic deformation and then adjust the stress and strain terms until the yield criterion and the constitutive laws are satisfied.


## 4. Model Problem and Numerical Results

The problem studied in this work is that of a thick cylinder subjected to a gradually increasing internal pressure, with plane strain conditions being assumed in the radial direction of the cross section. In fact it is a part of a pipeline carrying natural gas which must be designed for a given variation of the internal gas pressure. This means that having considered the self-weight of the pipe, which is installed horizontally above the ground and having estimated the bending moments and the normal stresses on the cross sections, at several positions along the pipe, one then examines the stress state at selected Gaussian points on the surface of the isoparametric elements, to check if yielding occurs. It must be mentioned that the normal stresses due to bending moments are negligible compared to the values of the radial and of the hoop stresses in the plane of the cross section.The Von Mises yield criterion is assumed and the numerical solutions, obtained by implementing the procedure presented in the previous sections, are compared with the theoretical results of the elastic solution and the elastic-plastic solution developed by Prandtl and Reuss.A 4-point Gauss-Legendre quadrature rule is used in order to estimate the values of the elements of the stiffness matrix of an element of the structure(42)Kij(e)=∬-1+1[Bi(e)]TD(e)Bj(e)dV(e),
where dV(e) is the elemental volume (here is an elemental surface) given by (38). Using double symmetry of the problem we work only on a quartet of the cross section, and the solid wall is discretized into NE 4-point isoparametric elements having all together n nodes. The problem was solved for an inner radius ri equal to 100 mm and an outer radius ro which was varied from 140 to 200 mm. The internal pressure Pi was also varied from 0 MPa to 192 MPa. The atmospheric pressure (about 0.1 MPa) is neglected being extremely small compared with the internal pressure on the pipe. In Figure 9 a finite element mesh is presented for a quartet of the cross section with an outer radius equal to 140 mm. But the solution presented here is for ro=200 mm.Figure 9
A finite element mesh for a quarter of the cross section.The material of the pipeline is carbon steel with Young’s modulusE=210 MPa, Poisson’s ratio ν=0.3, and a uniaxial yield stress σy=240 MPa. The strain hardening parameter is H′=0 MPa; that is, there is an elastic-perfectly plastic behavior for the material. A FORTRAN 77′ code was used to solve the problem on the computer. Only the nonzero elements of the global stiffness matrix are stored in the computer’s memory. The elastic problem was solved for two values of the internal pressure: Pi=80 and 100 MPa. For the elastic-plastic case the values of the internal pressure were Pi= 120, 140, 160, 180, and 190 MPa. Also, the following cycling loadings were examined:(i)
140→0→140→0→140→0⋯    (10 cycles)(ii)
180→135→90→45→0→180→⋯ (10 cycles).The pressure/radial displacement characteristics for the elastic and for the elastic-plastic loading are shown in Figure10. A good agreement between the numerical and analytical solution, given by Prandtl-Reuss, is evident. Indeed, in using the Von Mises criterion the analytic approach gives the value of 103.75 MPa as the value for the initiation of yielding. The corresponding value on the graph of Figure 10 is 103.33 MPa. Thus, a very good agreement is observed between the numerical value and the Prandtl-Reuss analytic solution. In fact, the numerical and analytic curves practically fall one upon the other in Figure 10.Figure 10
Plot of internal pressurePi versus δ (displacement of the internal surface of the pipe).In the numerical studies, collapse was deemed to have occurred if the iterative procedure diverged for an incremental load increase. For the specific problem the collapse load given by the analysis previously described was found to be equal to 192 MPa. Indeed, for this value of internal pressure or for greater values, there was no convergence in the numerical procedure and the program on the computer had to be stopped.In Table1 the numerical results obtained for Pi=80 MPa (elastic case) are compared with the analytical solution. One may observe that in almost all positions the theoretical results coincide with those obtained by the FEM (isoparametric element approach).Table 1
Numerical and analytical results for the elastic case.
r (mm)
σ
r (MPa)
σ
θ (MPa)
FEM
Analytic solution
FEM
Analytic solution
104.21
−71.55
−71.55
124.89
124.88
115.58
−52.89
−52.93
106.22
106.27
124.21
−42.40
−42.47
95.80
95.80
135.76
−31.18
−31.20
84.52
84.54
146.32
−23.16
−23.15
76.49
76.49
163.64
−13.14
−13.17
66.48
66.50
176.32
−7.64
−7.64
60.98
60.98
193.64
−1.77
−1.78
55.10
55.10In Table2 we present the results obtained by running the program for Pi=169.36 MPa. On the same table the analytical results given by the theoretical graphs of the Prandtl-Reuss method, for the same value of the internal pressure, are also presented. Comparison between the two sets of values indicates, in general, a good agreement between numerical and theoretical results. Some deviations from the analytical results are mainly due to errors in reading the values on the Prandtl-Reuss theoretical graphs and the assumption of linear interpolation for intermediate values. Since we have a plane strain condition the axial stress σz is found from σz=ν(σr+σθ). It is always less than σθ.Table 2
Numerical and analytical results for the elastic-plastic case for an internal pressurePi = 169.36 MPa.
r (mm)
σ
r (MPa)
σ
θ (MPa)
FEM
Prandtl-Reuss
FEM
Prandtl-Reuss
104.21
−157.89
−163.20
118.79
113.93
115.76
−128.83
−135.48
148.28
143.18
124.21
−109.27
−113.93
167.79
166.28
135.76
−84.69
−89.30
191.73
184.75
146.32
−63.59
−67.74
210.00
203.23
198.64
−4.86
−6.16
151.30
160.12Figure11 shows the graph of σr versus r for an internal pressure Pi=180 MPa, and Figure 12 presents the plot of σθ versus r for the same internal pressure. On the σθ versus r graph one may clearly see the position of the elastic-plastic boundary, that is, at a distance 165 mm from the center of the pipe. Figure 13 shows the distribution of the tangential residual stresses which appear on the cross section after the end of the last cycle of the second loading program which is discussed earlier. At the last two cycles of this loading program a solidification of the maximum effective plastic strain (max|εp(eff)|) was observed at the value 0.1885×10-2. This observation is in accord to the theory which has been developed, the last three decades, around the so-called “shakedown phenomenon.” It is a phenomenon which takes place when cyclic loading and unloading occur elastically (i.e., within the yield surface of the chosen yield criterion) after the initiation of the plastic deformation in the structure.Figure 11
Plot ofσr versus r for an internal pressure Pi=180 MPa. Notice the discontinuity in entering form the gas to the solid material.Figure 12
Plot ofσθ versus r for an internal pressure Pi=180 MPa.Figure 13
Distribution of residual hoop stressesσθ versus r after cyclic loading.Since this is a deep theory which is much discussed in the literature (e.g., [16]), there will be no further discussion here.Finally, it is of interest to have a look at the total CPU time used in running the FORTRAN code on the multiuser Cyber. Table3 shows the execution time and the number of iterations required to achieve convergence at each load. Convergence occurs when |Δψ| becomes less than |error|.Table 3
Execution time for different values of internal pressure.
Internal pressure (MPa)
Number of iterations until convergence occurs
CPU time (sec)
100
1
7.91
120
2
13.97
140
3
19.50
160
3
20.35
180
4
27.40
190
6
39.06


## 5. Conclusions

In this work the problem of a steel pipe carrying natural gas was solved by using the FEM with isoparametric elements. The elastic-plastic behavior of the material of the pipe was studied by implementing the Von Mises yield criterion and aJ2′ flow rule of the mathematical theory of plasticity. The numerical results given by the FEM for both the elastic and the elastic-plastic loadings were compared favorably with the analytical results of the classical theory of elasticity and with the theoretical results of the Prandtl-Reuss solution.For cyclic loading and unloading performed elastically, after the commencement of yielding in the structure, the FEM was able to predict the “shakedown phenomenon.” Thus, residual stresses in the wall of the pipe after a complete cycle of loading give an acceptable state of stress which is within the safety limits. The prediction of this phenomenon and of its residual stresses is probably the most important capability of the FEM for this class of problems together with its utility as a valuable tool for the design of such pipes to avoid collapse under high internal pressures.



---

*Source: 267095-2013-05-29.xml*