Large scale validation of a new isotropic hardening law for the barcelona plastic damage model under ultra low cycle fatigue conditions

IS - Modeling of Plasticity and Damage Under Cyclic LoadingLarge scale validation of the Barcelona plastic damage model with a new hardening law under Ultra Low Cycle Fatigue conditions XIII International Conference on Computational Plasticity. Fundamentals and Applications COMPLAS XIII E. Oñate, D.R.J. Owen, D. Peric and M. Chiumenti (Eds) LARGE SCALE VALIDATION OF A NEW ISOTROPIC HARDENING LAW FOR THE BARCELONA PLASTIC DAMAGE MODEL UNDER ULTRA LOW CYCLE FATIGUE CONDITIONS LUCIA G.BARBU*, XAVIER MARTINEZ†, SERGIO OLLER†AND ALEX H. BARBAT† *, † Centre Internacional de Metodes Numerics en Enginyeria (CIMNE) Universidad Politécnica de Cataluña Campus Norte UPC, 08034 Barcelona, Spain e-mail: [email protected] Key words: Plasticity, Ductile Fracture, Ultra Low Cycle Fatigue, Pipelines, Life Prediction, Constitutive Modelling. Abstract.This paper presents the results of finite element simulations made on a bent pipe subjected to an in-plane variable cyclic displacement combined with internal pressure. The results of the numerical analyses will be compared to experimental ones. The constitutive model used for the simulation of Ultra Low Cycle Fatigue (ULCF) loading and the hardening–softening law used are only briefly touched upon. Special emphasis is put on the correlation between the failure mode and the internal pressure applied. 1INTRODUCTION The present work is centred on the large scale validation of the Barcelona plastic damage model, proposed by Lubliner et al.[1]. An innovative application is given to this formulation by considering it for the cyclic loading case and incorporating a Friederick-Armstrong kinematic hardening law that allows the description of phenomena like cyclic ratcheting (under stress control conditions) or cyclic stress relaxation (under strain control or elastically constrained conditions). A new isotropic hardening law is also developed especially for steel materials, designed to reproduce their hardening and softening performance under monotonic and cyclic loading conditions. The exact expression of the constitutive law and the thermodynamic formulation of the model are presented in Martinez et al. [2][3]and Barbu et al.[4]. 2ULCF CONSTITUTIVE MODEL This work will not describe the complete plastic damage model, as it can be obtained from references [1] and [2]. The yield surface is defined by a function Fthat accounts for the residual strength of the material, which depends on the current stress state, the temperature and the plastic internal variables. This Ffunction has the following form, taking into account isotropic and kinematic 248 Lucia G. Barbu, Xavier Martinez, Sergio Oller and Alex H. Barbat plastic hardening (Bauschinger effect - Lemaitre and Chaboche [5]), 0) ,,()() ,,( ≤− −= θκαθγ p ijijij p ij SKSf SF (1) where )( ijij Sf α − is the uniaxial equivalent stress functions depending of the current value of the stresses ij S , ij α the kinematic plastic hardening internal variable, ) ,,( θκ p ij SK is the plastic strength threshold, p κ is the plastic isotropic hardening internal variable, and θ is the temperature at current time t ([1], [2],[6],[7] and [8]). 2.1 Kinematic Hardening Kinematic hardening accounts for a translation of the yield function and allows the representation of the Bauschinger effect in the case of cyclic loading. This translation is driven by the kinematic hardening internal variable ij α which, in a general case, varies proportionally to the plastic strain of the material point [5].There are several laws that define the evolution of this parameter. Current work uses a non-linear kinematic hardening law, which can be written as: pdE c ijk P ijkij    αα −= (2) Where k c and k d are material constants, p ij E is the plastic strain, and p  is the increment of accumulative plastic strain, which can be computed as: p kl p ij E Ep   :3 / 2 ⋅= . 2.2 Isotropic Hardening The evolution of isotropic hardening is controlled by the evolution of the plastic hardening function K , which is often defined by an internal variable p κ . The rate equation for these two functions may be defined, respectively: pp p kk Eh S G hH h HK       ⋅=       ∂ ∂ ⋅=⋅= ⋅=⋅ = kkk : λλ κ κλ (3) where k denotes scalar and k states for a tensor function. Depending on the functions defined to characterize these two parameters different solid performances can be obtained. 3NEW ISOTROPIC HARDENING LAW In the Barcelona model defined in [1], the laws defined are driven by the fracture energy of the material. This work presents a new law, especially developed for steel materials, that has been designed to reproduce their hardening and softening performance under monotonic and cyclic loading conditions. This law also depends on the fracture energy of the material and is derived from the hardening softening law presented in [1] and [2]. 249 Lucia G. Barbu, Xavier Martinez, Sergio Oller and Alex H. Barbat 3.1 Fracture Energy Classical fracture mechanics defines the fracture energy of a material as the energy that has to be dissipated to open a fracture in a unitary area of the material. Thus, the fracture energy per unit volume is obtained as the fracture energy of the material divided by the fracture length. This fracture length corresponds to the distance, perpendicular to the fracture area, in which this fracture propagates. In a real section, this length tends to be infinitesimal. However, in a finite element simulation, in which continuum mechanics is applied to a discrete medium, this length corresponds to the smallest value in which the structure is discretized: the length represented by a gauss point. Therefore, in order to have a finite element formulation consistent and mesh independent, it is necessary to define the hardening law in function of the fracture energy per unit volume ([1], [8]). This value is obtained from the fracture energy of the material, f G , and the size of the finite element in which the structure is discretized. 3.2 Hardening Function and Hardening Internal Variable The hardening function defines the stress of the material when it is in the non-linear range. There are many possible definitions that can be used to fulfil the rate equations for the plastic strength threshold (3). Here the use of a function that describes the evolution of an equivalent uniaxial stress state is proposed, like the one shown in Figure 1. Figure 1. Evolution of the equivalent plastic stress This equivalent stress state shown in Figure 1 has been defined to match the uniaxial stress evolution described by most metallic materials. This curve is divided into three different regions. The hardening internal variable, p κ , accounts for the evolution of the plastic hardening function, K . In current formulation p κ is defined as a normalized scalar parameter that takes into account the amount of volumetric fracture energy dissipated by the material in the actual strain-stress state. This is: d tES g t t p f p∫ = = 0 : 1 κ (4) Using the definition of the hardening internal variable defined in equation (7), it is possible to define the expression of the hardening function as: 250 Lucia G. Barbu, Xavier Martinez, Sergio Oller and Alex H. Barbat )( peq SK κ = (5) It can be easily proven that the hardening function and internal variable defined in equations (7) and (8) fulfil the rate equations (3). And the k h and k h functions defined in this expression become: t p eq k g S h S h = ∂ ∂ = k κ (6) The exact numerical expressions used to define the new hardening law can be found in [2]. The constitutive law described in Sections 2 and 3 has been implemented in the in-house code PLCd [9],[10],[11]. The code was programmed to allow OpenMP parallelization which greatly reduced the computational cost of the large scale FE simulations. 4MATERIAL CALIBRATION The material characteristics for the numerical simulations will be obtained by conducting a calibration analysis on small scale specimens. The hardening –softening law presented in section 3 requires of the following material parameters: •εp – σppoints obtained from uniaxial monotonic tensile tests necessary for curve fitting. They are important for a correct representation of the tendency of the monotonic curve. •Kinematic coefficients in accordance with the type of hardening chosen. They are important for the exact adjustment of the monotonic curve and for an accurate description of the hysteresis loop. •Equivalent plastic deformation, εp1eq , at which the linear region starts. This parameter is important both in the monotonic curve and in the overall cyclic behaviour as it ensures a stable behaviour throughout the fatigue life. •Equivalent plastic deformation, εp2eq , at which softening starts. •Fracture energy Gfrequired in the monotonic curve for adjusting the slope of the softening behaviour and in the cyclic behaviour for correctly calibrating the fatigue life of the specimen. When the entire energy is spent the specimen is considered completely fractured. In order to exemplify the calibration process, a smooth X60 specimen was chosen from the experimental program ran by Pereira et al. [12] (Figure 2). Figure 2. Geometry of the specimen used for calibration 251 Lucia G. Barbu, Xavier Martinez, Sergio Oller and Alex H. Barbat The specimen was meshed into 3456 quadratic hexahedral elements with 20 nodes each and 27 integration points, adding to a total of 17165 mesh nodes. At this point in the calibration procedure the series of chosen points have to follow the general tendency of the monotonic curve without reaching the same level of stress. The density of the points is recommended to be constant and quite high so that the polynomial interpolation can be effective. For this simulation the points were interpolated by a 5th order polynomial function. The ArmstrongFrederick kinematic hardening function was used for this simulation. The kinematic coefficients chosen were 10 1 1 06 ×= k and 400 2=k . In Figure 4the effect of the kinematic hardening on the monotonic curve can be observed. In order to obtain this behaviour a value of 0.1 was used for the εp1eq parameter. Figure 3. Monotonic and cyclic (Δε=5%) stressstrain curve. Numerical vs. experimental. It can be seen that taking into account the kinematic hardening causes the resulting stressstrain curve to elevate until it reaches the experimental monotonic ones. Also, the exact shape of the transition zone from linear to nonlinear is determined by the kinematic coefficients. In choosing the kinematic coefficients a compromise must be made between the accuracy of the monotonic behaviour and of the cyclical one. The experimental result chosen as the calibration test is the Δε=5% case, that had an experimental fatigue life of 100 cycles. With a value of 13 for the εp2eq and a fracture energy of 2.7 x 106 Nm/m2a total fatigue life of 100.35 cycles has been obtained from the numerical simulation. With these values, softening started in the 86th cycle, close to the end of the experimental life and a very low amount of energy was left for the softening branch so that it could be spent in a reduced number of cycles. After the adjustment of the Δε=5% case, simulations were ran with the exact same material and with strain amplitudes of 2% and 8%. The results can be seen in Figure 4. Table 1. Polynomial coefficients for the X60 material as obtained from curve fitting Coefficient no. 1 380000000,00 Coefficient no. 2 326947332,25 Coefficient no. 3 -861244568,93 Coefficient no. 4 1103673406,83 Coefficient no. 5 -657861660,66 Coefficient no. 6 147925577,48 252 Lucia G. Barbu, Xavier Martinez, Sergio Oller and Alex H. Barbat In Table 1 and Table 2a summary of the material properties as resulting from the calibration can be seen. A polynomial of the 5th degree was chosen for the curve fitting zone of the hardening function. Table 2. Material parameters for the numerical model for an X60 steel Young Modulus 1.95·10 5 MPa Poisson Modulus 0.30 Elastic Stress ( eq Y σ ) 380 MPa Plastic Strain Limit for region 1( p E 1 )10 % Plastic Strain Softening ( p E 2 ) 1300 % C1 kinematic hardening 6.0·10 4 MPa C2 kinematic hardening 400 Fracture Energy 2.7 MN·m/m 2 Figure 4. Comparison between numerical and experimental fatigue life for different strain amplitudes when calibrating with Δε=5% 5BENT PIPE UNDER VARIABLE CYCLIC LOADING 5.1 Geometry of the model Following the validation of the constitutive model made in [2] on small scale specimens, the model is to be applied to large scale numerical simulations of a bent pipe. The geometry of the model, boundary conditions and the sequence of loading are in accordance to the experiment made by Schaffrath et al. [13]. The specimen consists of a bended middle section (elbow pipe) and a straight pipe section at each end of the elbow. The fillet radius of the elbow pipe is three times the pipe diameter (R=3 x D). For the length of both straight pipe sections a value of five times the diameter (L=5 x D) was used, whereby the influence of the load introduction can be neglected. For the numerical model a specimen made of X60 steel has been chosen from the experimental program, with a diameter of 16 inches and a wall thickness of 9.5 mm. The pipe has an elbow angle of 90°. In Figure 5a view of the mesh is shown. For this simulation quadratic hexahedral elements were used, each with 20 nodes and 27 integration points. The mesh consisted of 42853 elements and 213415 nodes. Three elements were considered in the pipe thickness. 253 Lucia G. Barbu, Xavier Martinez, Sergio Oller and Alex H. Barbat Figure 5. Mesh of hexahedral quadratic elements 5.2 Loading history The loading history is based on the actual load history of the experimental test done by Schaffrath et al. [13]. The loading scheme was decided by the authors in accordance with the ECCS procedure ECCS-Nr. 45-1986 [14]. For practical reasons it was decided to neglect the mostly small difference between the compressive and tensile yield strain by choosing an average value 2 / )( −+ += y yy e ee as the reference amplitude. In Table 5the experimental loading sequence is described as a function of y e .The value adopted for this parameter was set by Schaffrath et al. [13] at ±82mm. Table 3. Loading sequence for SP2 specimen Step Amplitude Number of cycles 1 0.25 y e 1 2 0.50 y e 1 3 0.75 y e 1 4 1.00 y e 1 5 1.50 y e 3 6 2.00 y e 3 7 2.50 y e 3 8 3.00 y e 3 9 3.50 y e 3 10 4.00 y e 3 11 4.40 y e 27 The entire loading sequence is comprised of 49 cycles with increasing amplitude, 44 of which have amplitudes in the plastic range. The reversion factor of the applied displacement is -1. The pipe is also submitted to internal pressure. First, it is loaded until a level of internal 254 Lucia G. Barbu, Xavier Martinez, Sergio Oller and Alex H. Barbat pressure equal to 20 bars. Afterwards, it is submitted to the varying cyclic displacement presented in Table 5. The experimental test has shown that the internal pressure also oscillates when the cyclic displacement is applied. The boundary conditions of the model were chosen in accordance with the setting of the experiment. One end of the model as presented in Figure 5has its displacement blocked in the x, y and z direction while in the other end the cyclic displacement is applied in the z in-plane direction. In the numerical simulation the loads have been applied in two stages. First, the internal pressure was applied and, in this stage, one end of the pipe was clamped and on the other end the pipe was only allowed in-plane gliding. The variable displacement was applied in the second stage on the deformed geometry obtained from applying the internal pressure. The movement was restrained in the two directions perpendicular to the in-plane one. 5.3 Material characteristics The exact calibration procedure for the X60 material can be seen in Section 4, as well as the properties defined, that are presented in Table 1 and Table 2. 5.4 Results and discussion In Figure 6the comparison between the experimental forcedisplacement curve and the numerical one can be seen. The numerical curve is in very good agreement with the experimental one taking into consideration that the material calibration was done on small scale specimens with different experimental results. It can be seen how in compression the constitutive equation tends to underestimate the maximum force level, while in traction the opposite tendency is present. Figure 6. Force-displacement curve. Experimental vs. numerical Regarding the fatigue life, the simulation lasted a total of 41.75 plastic cycles as compared to the experimental life of 44 complete plastic cycles. This result also shows a good 255 Lucia G. Barbu, Xavier Martinez, Sergio Oller and Alex H. Barbat agreement between the experiment and the numerical simulation. Figure 7illustrates the deformed shape of the geometry in the last step of the analysis and presents the distribution of the plastic internal variable of the model. The deformed shape is represented with a scale factor of 2 in order to better reflect the general tendency. Only the central zone of the elbow is shown, as this is the zone where nonlinear effects appear. It can be seen that the failure mode resulting from the numerical simulation is by cross-sectional ovalization with a crack opening in the longitudinal direction of the elbow, at its flank. The comparison of the failure mode obtained with the numerical simulation with the failure mode obtained in the experimental test (Figure 8) shows that the model has been able to capture the number of cycles to failure but has not been able to capture the failure mode shown in the experiment. Figure 7. Distribution of the plastic internal variable of the model on the deformed shape (x2) Under extreme loading conditions, such as the high repeated incursions in the nonlinear zone that the imposed displacement in this case causes, elbows exhibit two different failure modes. These are either significant cross-sectional ovalization or local buckling, as reported by the experimental work described in references [15] to [23]. An important conclusion can be drawn from the work above mentioned. The first failure mode can be generally found when the elbow internal pressure is relatively low compared to the yield pressure, as is the case in the numerical simulation presented above. The second failure mode, occurring due to local buckling is habitual in the cases where internal pressure is significantly high, of the order of the yield pressure. This is the failure mode yielded by the experiment made by Schaffrath et al. [13] as it can be seen in Figure 8. For the case considered, the internal pressure applied to the elbow is 20 bars, which leads to a stress value, according to Barlow’s formula: Mpa tp D 7 7.4 22 == σ . This is less than 10% of the yield strength for an X60 steel and, consequently, the pressure applied is less than 10% of the internal yield pressure. This puts us in the first yield mode. Summarizing, the number of cycles the simulation lasted and the force displacement curve are in good accordance between the numerical model and the experiment and the numerical failure mode is different from the experimental one but justifiable given the low internal pressure. 256