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, a 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.
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 Figure 1.
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 (Figure 2), 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.
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 in the wall of the tube ( being the yield stress) versus the normalized distance from the center ( 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 Section 2, 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 where and are the stress and strain components respectively, and is the tensor of elastic constants which, for an isotropic material, has the explicit form where and are the Lamé constants and is the so-called Kronecker delta, a useful symbol in mechanics defined by 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 where is a function of stresses and 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 is a material parameter to be determined experimentally. It is the so-called yield stress and usually is denoted as .
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 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 where and are the second and third invariants where 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 reaches a critical value as where, again, is the yield stress. The second deviatoric stress invariant, , can be explicitly written as Yield criterion (9) may be further written as Quantity is termed effective or equivalent stress.
Figure 3 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 . The two-dimensional plot of this yield surface has an elliptic form (Figure 4). Now, a physical meaning of the constant can be obtained by considering the yielding of materials under simple stress states. So, the case of uniaxial tension requires that value is the uniaxial yield stress.
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 Figures 5 and 6.
In the case where the yield surfaces are a uniform expansion of the original yield curve, without translation, as shown in Figure 5, 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 stress , to the plastic deformation by means of the hardening parameter . 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, , only. Then, where are the plastic components of strain occurring during a strain increment. Alternatively, can be related to a measure of the total plastic deformation termed the effective or equivalent plastic strain which is defined incrementally as For situations where the assumption that yielding is independent of any hydrostatic stress is valid, we have . Consequently, (13) can be rewritten as Then, the hardening parameter is defined as where is the result of integrating over the strain path. This behavior is termed as strain hardening.
By considering an isotropic hardening model, all the points on the yield surface represent plastic states, while elastic behavior is characterized by . At plastic state, , the incremental change in the yield function due to an incremental stress change is Then, if , elastic unloading occurs (stress point is inside the yield surface); , neutral loading takes place (stress point remains on the yield surface); , 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: upon loading, and when 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 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 is as follows: where and 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 potential , so that where 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 must be a function of and but as yet cannot be determined in its most general form. However, the relation 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 is a valid one since it has been postulated that both are functions of and and such an assumption gives rise to an associated theory of plasticity. In this case (21) becomes and is termed the normality rule. The derivative 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 -dimensional space to give a vector which is normal to the yield surface. For the particular case of , we have Then, (22) becomes 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
The hardening law could just be easily expressed in terms of the effective stress (since it is proportional to ) to give, for the strain hardening hypothesis (15), or differentiating, In (27) the derivative is the hardening function and it is the slope of the uniaxial stress/plastic strain curve for which and . Thus, it can be further expressed as follows: where is the elastoplastic tangent modulus (Figure 8).
The theoretical expressions developed so far can now be converted to matrix form. The yield function, first defined in (4), can be rewritten as where is the stress vector and is the hardening parameter which governs the expansion of the yield surface. In particular, from (12), we have for the work hardening hypothesis. Rearranging (29), we get By differentiating , we have which is expressed as with the definitions Zienkiewicz et al. [15] give an elegant proof of the following result: Thus, parameter 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 where for node , vector is the vector of nodal variables, is the vector of virtual nodal variables, is the matrix of global shape functions, and is the global strain-displacement matrix. For a 4-node isoparametric element shape function, is expressed in the local system of the element as . The total number of nodes on the whole mesh is .
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: where are the boundary tractions which are also denoted by (external applied forces) and are the distributed internal loads per unit volume. Since (36) must be true for any arbitrary value, we have Here, for the plane strain problem, the discretized elemental volume represents an element area which on the local coordinate system of the element is given as where 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 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 . For the purpose of evaluating the global tangential stiffness matrix at any stage, the incremental form of (4) must be employed. Thus, within an increment of load, we have where 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 where 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 4-point isoparametric elements having all together nodes. The problem was solved for an inner radius equal to 100 mm and an outer radius which was varied from 140 to 200 mm. The internal pressure 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 mm.
The material of the pipeline is carbon steel with Young’s modulus MPa, Poisson’s ratio , and a uniaxial yield stress MPa. The strain hardening parameter is 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: and 100 MPa. For the elastic-plastic case the values of the internal pressure were = 120, 140, 160, 180, and 190 MPa. Also, the following cycling loadings were examined:(i) (10 cycles)(ii) (10 cycles).
The pressure/radial displacement characteristics for the elastic and for the elastic-plastic loading are shown in Figure 10. 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.
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 Table 1 the numerical results obtained for 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).
In Table 2 we present the results obtained by running the program for 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 is found from . It is always less than .
Figure 11 shows the graph of versus for an internal pressure MPa, and Figure 12 presents the plot of versus for the same internal pressure. On the versus 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 was observed at the value . 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.
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. Table 3 shows the execution time and the number of iterations required to achieve convergence at each load. Convergence occurs when becomes less than .
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 a 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.