Finite element simulation of plasticity induced crack closure with different material constitutive models

The study presented in this paper analyses the mechanical effects of material constitutive modelling on the numerical prediction of plasticity induced crack closure. With this aim, an elastoplastic stress analysis of a MT specimen was conducted using an implicit three dimensional finite element program. Two materials were studied: an Aluminium Alloy and a High Strength Steel. Several constitutive models were used to describe their cyclic behaviour, ranging from pure isotropic hardening or pure kinematic hardening models to combined isotropic plus kinematic hardening models. Numerical results showed clear differences in plastic behaviour and crack closure predictions for the different types of mechanical models used to describe the mechanical behaviour of the materials. The mechanisms of opening stress stabilization, usually observed in numerical simulations, are explained in this work by analysing the evolution of plastic deformation along the crack flanks. The same type of plastic deformation stabilization behaviour was observed independently of the hardening model in use.


Introduction
Crack closure concept is often used to explain some aspects of fatigue damage in cracked components being currently accepted that, when the crack is closed, the external load produces negligible fatigue damage in the component. Many analytical or semi-analytical models have been developed to explain crack closure and several agents were identified as promoting this phenomenon, such as: roughness (Roughness Induced Crack Closure), oxides (Oxide Induced Crack Closure) and plastic deformation (Plasticity Induced Crack Closure -PICC). Among these, plasticity was considered the primary crack closure mechanism under many conditions [1][2][3][4].
Since a detailed examination and quantification of PICC by using experimental testing procedures is a relatively complex task, finite element numerical simulation became one of the most important techniques of analysis considering this subject. In 1989, McClung and Sehitoglu [5] published one of the first papers giving some general guidelines for the numerical simulation of PICC problems. From the time of this publication, to our days, it becomes consensual that several aspects must be considered to correctly simulate PICC, namely: an accurate selection of finite element parameters (type and size of crack front elements, and numerical parameters characteristic of each type of numerical simulation program), the selection of a proper crack increment strategy (propagation load, size of crack increments, number of load cycles between crack increments), the definition of the closure level (last contact of the crack faces, inversion of stresses, etc.), and finally, a correct choice of the material models used to describe the cyclic behaviour of the material in study.
In fact, mechanical modelling of metals elastoplastic behaviour requires the definition of a flow rule, a yield surface and hardening laws that determine the evolution of the yield surface with plastic deformation [6][7][8]. For the specific case of cracks growing under cyclic loading, modelling of cyclic plasticity inside the crack tip area implies a correct description of both isotropic and kinematic hardening mechanisms, which makes mixed hardening constitutive models necessary.
A brief review concerning material modelling in the numerical simulation of crack closure reveals that most of the published works considers the metallic materials stress-strain response modelled as elastic-perfectly plastic [9][10][11][12][13][14][15][16][17][18][19][20][21][22] or bi-linear [5,9,12,20,[23][24][25][26][27][28][29][30][31]. Some of the published works employed material constitutive relationships that ignore kinematic hardening [10,[32][33][34]. Considering the relevant importance of work-hardening modelling in the accurate numerical simulation of PICC, numerical studies comparing results obtained with different mechanical models were also published. McClung and Sehitoglu [23] and Ashbaugh et al. [35] compared numerical results obtained using linear and power-law material models for the numerical simulation of PICC. These authors showed that power-law models conduct to higher opening values than the linear models. In 1999, Pommier and Bompard [33], published a plane strain 2D numerical work comparing crack closure behaviour registered in numerical simulations performed by using purely isotropic and purely kinematic hardening models. These authors found a strong interaction between the cyclic plastic behaviour of the material and the parameters associated with PICC, such as, stress ratio and overload ratio. Some years later, Roychowdhury and Dodds [4] also compared purely isotropic and purely kinematic hardening models numerical results in a 3D numerical investigation concerning small scale yielding fatigue crack growth. In this work they had shown that the pure isotropic hardening model strongly overestimates the computed crack opening load.
Studies comparing the ability of different mixed hardening models in describing PICC were also published. Jiang et al. [1] published a comparative study considering three different mixed hardening models: a pure kinematic hardening rule of Prager-Ziegler with an elastic-perfectly plastic (EPP) stress-strain relationship, a pure Kinematic hardening rule of Prager-Ziegler with a bilinear (BL) stress-strain relationship and a Kinematic hardening model developed by Jiang and co-workers [36][37][38]. These authors found a strong dependence of crack closure stabilized values on the mechanical model in use. They also found a severe dependence of crack closure results on crack tip finite element size for the first two mechanical models. Jiang et al. [1] also discussed the inability of the elastic-perfectly plastic model, with the pure kinematic hardening rule of Prager-Ziegler, and of the bilinear model to describe stress ratcheting and stress relaxation during cyclic loading. Elliyn and Ozah [8] also published a study concerning the effect of material model in describing PICC under variable amplitude cyclic loading. The authors compared numerical results from a pure kinematic hardening model with that from a multi-surface yield model, the Ellyin-Xia model [39], and found significant differences in the results.
However, to the author's knowledge, numerical simulation works comparing both non-linear isotropic and kinematic hardening models are still rare. In this work, numerical results from two different non-linear mixed hardening models will be compared with non-linear pure isotropic and pure kinematic hardening models. These mechanical models are currently strain in the y direction GP Gauss point K max maximum stress intensity factor used in the numerical simulation of large plastic deformation problems, such as metal forming processes, for which material constitutive modelling is also of crucial importance in order to achieve accurate predictions of the technological parameters in study. Several benchmarks are proposed every year in order to evaluate the performance of the metal forming numerical simulation programs and to discuss the latest developments in the constitutive modelling of the plastic behaviour of commercial metallic alloys [6]. These models are currently used to describe the elastoplastic behaviour of metallic materials under intense forming sequences that may involve complex strain paths and/or large accumulated strains. Considering the enormous amount of experimental work developed in metals plastic deformation modelling and the huge amount of material parameter data already available, the possibility of applying these models in the study of PICC would become very interesting. This paper focuses precisely on the numerical simulation of the crack closure phenomenon by using some of these advanced constitutive models to describe the cyclic behaviour of the materials. The PICC numerical results obtained considering two different non-linear pure isotropic hardening models, one non-linear pure kinematic hardening model and two mixed hardening constitutive models, for two different materials, an Aluminium Alloy and a High Strength Steel, will be compared in the next.

Material constitutive models
During the plastic deformation of metals, hardening mechanisms enforce a continuous evolution of the yield stress of the material during the loading process. From a macroscopic point of view, the initial yield behaviour of any material and its evolution during loading is modelled by establishing a yield criterion and isotropic and kinematic work-hardening laws that determine the evolution of the initial yield surface with loading. Both components of hardening are associated to the plasticity criterion through a plastic potential f that can be expressed as f ðr 0 À X; a k Þ ¼ rðr 0 À X; a k1 Þ À Yða k2 Þ; where r and Y are equivalent stresses associated with the yield criteria and isotropic work-hardening law, respectively, r 0 is the stress deviator, X is the back stress associated with kinematic hardening, a k1 is a set of scalar parameters related to the plasticity criteria and a k2 are internal parameters associated with the isotropic component of hardening.
In present work, an anisotropic yield criterion [40] was considered, which is expressed by the quadratic function Fðr yy À r zz Þ 2 þ Gðr zz À r xx Þ 2 þ Hðr xx À r yy Þ 2 þ 2Ls 2 where r xx , r yy , r zz , s xy , s xz and s yz are the components of the effective stress tensor ðr 0 À XÞ defined in the orthotropic frame and F, G, H, L, M and N are coefficients that characterize the anisotropy of the material (Table 1). However, it is important to enhance that despite anisotropy was considered in the numerical simulations, previous tests from present authors, performed using Von Mises criteria instead of Hill criteria, showed that there is no influence of materials anisotropy in PICC numerical results. To describe the isotropic component of hardening, two classical phenomenological models were considered in this work: the Swift law [43,44] Y ¼ Cðe 0 þ e p Þ n ð3Þ and the Voce law [42] Y ¼ Y 0 þ R sat ð1 À e Ànv e p Þ: ð4Þ In previous equations, Y is the equivalent flow stress, e p is the equivalent plastic deformation and C, n, e 0 and Y 0 , R sat and n v are material constants that can be determined experimentally for the Swift and Voce models, respectively. The kinematic component of hardening, that determines the back stress evolution, was modelled by the non-linear model [43,44] where C x and X sat are material parameters and _ e p is the equivalent plastic strain rate. Based in previous equations, several material models were tested in this work: Pure isotropic hardening behaviour described by the Swift equation (SIH). Pure isotropic hardening behaviour described by the Voce equation (VIH). Mixed hardening behaviour: isotropic hardening described by the Voce law combined with kinematic hardening (VIKH). Mixed hardening behaviour: isotropic hardening described by the Swift law combined with kinematic hardening (SIKH). Pure kinematic hardening behaviour (KH).
These constitutive models were used to describe the cyclic elastic-plastic response of two metallic materials: an Aluminium Alloy (AA 6016 -T4) and a High Strength Steel (DP600). These materials were previously tested (monotonic and Bauschinger tests were performed) and material parameters, listed in Table 1, were determined according to the five types of material models listed above by using a genetic based numerical algorithm [45,46]. From the experimental results it was concluded that the AA 6016 -T4 aluminium exhibit a pronounced saturation of the flow stress in monotonic loading and a rather small Bauschinger effect (<2% difference between the flow stresses in forward and reversed directions). On the other hand, the results obtained for the steel DP 600 showed continuous hardening in monotonic loading and strong Bauschinger effect. During the reversed deformation in Bauschinger tests, work-hardening stagnation followed by a resumption of workhardening was also observed. This behaviour will be favourable to reversed plasticity.
The material parameters determined for these two alloys, with very different mechanical behaviour, were used in present numerical study in order to determine the effect of work-hardening modelling on the numerical prediction of PICC.

Numerical simulation procedure
The numerical simulations were performed with a 3D elastoplastic finite element code that makes use of an updated Lagrangian formulation and follows a fully implicit time integration scheme [47][48][49]. In this work, a 3D center cracked plate subjected to mode I loading was investigated. Due to symmetry conditions only 1/8 of the sample was simulated and the opposite crack face was modelled by placing a symmetry plane behind the growing crack front. The boundary conditions used in this work enforce frictionless contact conditions over this symmetry plane. As illustrated in Fig. 1, both plane stress and plane strain conditions were simulated in 0.1 mm thick specimens with straight cracks. 3D models were used because the numerical simulation program (DD3IMP) was exclusively developed to work with solid finite elements, but this does not represent a disadvantage in the numerical simulation of PICC, on the contrary. The finite element mesh used in the discretisation of the specimens is shown in Fig. 2. This mesh is highly refined near the crack front, to enable the numerical simulation of the local stress and strain gradients, and large at remote positions, to reduce the numerical effort. A mesh sensitivity study was previously performed in order to determine the optimum near crack tip finite element size. Elements with radial sizes (L 1 ) of 16 lm and 8 lm were considered for plane stress and plane strain conditions, respectively. The plane stress model had 2485 linear isoparametric elements and 5194 nodes, while the plane strain model had 3230 linear isoparametric elements and 6650 nodes. Only one layer of elements was considered along the thickness of the specimen. In order to avoid deterioration of results due to the artificial hydrostatic stresses resulting from the locking effect a selective reduced integration scheme was adopted in the numerical simulations. More details concerning the locking effect and the type of integration used can be found in [47].
The specimens were submitted to constant amplitude loading with stress ratio R = 0.1 and r max /r ys % 0.3, except for the aluminium alloy in plane strain conditions, which was submitted to r max /r ys % 0.5. Crack propagation was simulated by releasing, every two load cycles, all crack front nodes over the thickness. Previous numerical results indicate that only one load cycle between crack increments clearly is not enough, and that for more then two load cycles the small influence observed does not compensate the increase in the numerical effort [50]. Each crack increment (Da) corresponded to one finite element, i.e., Da = L 1 = 16 lm or 8 lm. The increment at minimum load was adopted to overcome convergence difficulties [2]. Numerical simulations were stopped after 20 or 40 crack propagation increments for the models with L 1 = 16 and 8 lm, respectively. Crack opening was determined by evaluating, in each layer, the contact status of the first node behind current crack tip. The crack was considered completely open when, at each through-thickness crack front location, it was detected that the first node behind current crack tip had lose contact with the symmetry plane. In the next, the opening level results will

Crack tip plastic deformation
The major aim of this work was to analyse the influence of material modelling on the numerical predictions of PICC. Since the numerical simulation of both the forward and reversed plastic zones is currently accepted as one of the main concerns in the numerical prediction of PICC, plastic deformation behaviour ahead of the crack tip was carefully analysed. Fig. 3 plots the normalized stress-strain curves (r yy -e yy ) registered at a Gauss point (GP) situated ahead of the initial crack tip position, as the crack propagates towards him Fig. 1c schematises the location of the initial crack tip (Position 1) and of the GP (Position 14). The r yy -e yy curves plotted in Fig. 3a-e were obtained for the aluminium alloy samples using the KH, SIH, VIH, SIKH and VIKH models, respectively. From the graphs it is possible to observe that, independently of the material model in use, at the end of the first cycle some plastic deformation occurs at the GP position, which indicates that this point is within the initial forward plastic zone. Compressive stresses arise at the GP position after a number of crack increments that varies according to the material model in use: 4 crack increments for the VIH and KH models, 5 crack increments for the SIKH and VIKH models and 6 crack increments for the SIH model. In all situations, compressive stresses increase as the crack tip approaches the GP position and start producing reversed plasticity. For example, for the KH model (Fig. 3a) reversed plasticity started after 8 crack increments. Since crack tip reached the GP position after 14 crack increments, it is possible to conclude that 6 elements of 16 lm were comprised inside the reversed plastic zone. For the VIKH model (Fig. 3d) 10 crack increments were applied before reversed plasticity initiates at the GP position, and therefore, only 4 elements were comprised inside the reversed plastic zone. Finally, for the SIH, VIH and SIKH models (Fig. 3b-d) 12 crack increments were necessary before reversed plasticity starts at the GP position, which indicates that the reversed plastic zone only encompassed 2 elements. After 15 crack increments, crack propagation extends ahead of the GP position, the level of stress applied to the GP becomes relatively low and the plastic deformation stops increasing. Similar trends were obtained for the steel (Fig. 4). For this material, the GP analysed was positioned at the 20th element, which means that propagation extended ahead of the GP position after 20 crack increments. Like in previous results, independently of the material model in use, it was observed that at the end of the first cycle the GP suffered some plastic deformation, which indicates that it was within the forward plastic zone. For this material, compressive residual stresses started arising at the GP position after 12 crack increments for all the material models in use and reversed plasticity took place after  16 crack increments for the SIH and VIH models (4 elements within the reversed plastic zone), 15 crack increments for the VIKH and SIKH models (5 elements within the reversed plastic zone) and 14 crack increments for the KH model (6 elements within the reversed plastic zone). These results show how sensitive is the numerical simulation of the forward and reversed plasticity at the crack tip to the material models used in the numerical simulations. According to Solanki et al. [11], there must be 3-4 linear elements within the reversed plastic zone, while Roychowdhury and Dodds [4] suggested 2-3 linear elements, in order to accurately simulate reversed plasticity. As previous results show, when the kinematic component of hardening is considered, reversed plasticity at the crack tip becomes more important and larger meshes can be used, reducing significantly the computational effort.
Under plane strain conditions, plastic deformation behaviour at the crack tip is different from that registered for the plane stress samples, since the through-thickness stress imposes an increase of the hydrostatic stress in the sample and reduces the deviatoric stresses responsible for plastic deformation. Therefore, modelling of reversed plasticity may be difficult, particularly when using pure isotropic hardening models. Considering this issue, in the numerical simulation of plane strain cases, the crack tip mesh area was refined to 8 lm elements. Fig. 5 presents r yy -e yy curves, for the aluminium alloy, registered as the crack propagates in plane strain conditions. Since the numerical results were very similar for all the material models in use, in this figure are only plotted the r yy -e yy curves relative to the numerical simulations performed with the KH and VIH models, which represent extreme situations among all the results registered for the aluminium alloy. The GP in study was located at the 28th element ahead of initial crack tip position and 29 crack increments were necessary until the crack tip propagates ahead of it. Due to the large amount of load cycles and very small plastic deformation values, it is very difficult to identify in the figure when forward plastic deformation was initiated at the GP position. However, it is possible to see that, in both numerical simulations, plastic deformation at the GP position remained very small until the crack tip approaches it. Amplifying the figure, it was possible to identify the beginning of reversed plasticity at the GP position. It was determined that three and two finite elements were within the reversed plastic zone for the KH and VIH models, respectively. This shows that an individual analysis of reversed plasticity is recommended for each physical situation in order to guarantee an adequate modelling, particularly in plane strain conditions. Similar r yy -e yy curves and sizes of plastic zones were obtained for the steel samples. A minimum of 2 elements were observed inside the reversed plastic zone.

Crack opening load stabilization
Another crucial aspect in the numerical predictions of PICC is the stabilization of the numerical crack opening load. The crack opening versus crack extension numerical results, obtained in plane stress conditions and considering the different material models, are plotted in Fig. 6a and b for the aluminium and steel, respectively. It is possible to observe in the figure that, for each material, independently of the material model in use, the opening values increase and stabilize with crack propagation. The number of crack increments necessary for crack opening load stabilization was lower for the aluminium than for the steel. In the graphs it is also possible to observe that the value of crack opening stabilized load varies according to the material model in use.
For the aluminium alloy (Fig. 6a), the higher opening stabilized loads were registered with the pure isotropic hardening models (VIH and SIH) and the SIKH mixed model. The small differences between the SIH and the SIKH models can be explained by the fact that this material is not very sensitive to Baushinger effect. In fact, very similar behaviour concerning reversed plasticity was already depicted when comparing the numerical results relative to these models in Fig. 3. The strong difference in crack opening numerical results, between the Voce based models (VIH and VIKH) and the Swift based models (SIH and SIKH), can be attributed to the different nature of these two equations. In fact, the Voce equation (Eq. (5)) is usually recommended to describe the mechanical behaviour of materials that present strong saturation in flow stress, like the aluminium alloys. For the steel (Fig. 6b), similar results were found for both pure isotropic (VIH and KIH) and for both mixed hardening models (SIKH and VIKH). The differences in crack opening values obtained with these two types of hardening models can be attributed to the strong sensibility of this material to the Baushinger effect. For both materials ( Fig. 6a and b), the numerical simulations performed with the pure kinematic hardening model (KH), conducted to the lowest opening load previsions.
The results obtained by fitting the Voce type equation (Eq. (4)) to the crack opening load numerical results are also plotted in Fig. 6a and b (lines connecting the symbols). The excellent fitting that can be observed for all results plotted in the graphs, shows the adequacy of this model to extrapolate crack opening results under plane stress conditions. Two approaches have been used to check the stabilization of the crack opening load. One of this approaches consists in the analysis of the asymptote of the Voce fitting model, while the other, is based in the study of the equivalent plastic strain ð e p Þ evolution along crack front, given by the derivative ð@ e p =@xÞ. In order to calculate the derivatives, the equivalent plastic strain along crack flanks was evaluated after 20 crack increments took place. The plastic strain derivatives for the aluminium alloy, in plane stress conditions, are plotted in Fig. 7a. Since it was used the one-node-per-two-cycles debonding scheme, crack increment and node position along crack front can be directly related. Analysing the graphs it is possible to observe a strong gradient in plastic deformation in a small region that coincides with the first nodes released. After this initial period, independently of the material model in use, all the derivatives approach zero as the distance from the initial crack tip increases, indicating a stabilization of plastic deformation along the crack front. Therefore, the gradient of equivalent plastic strain can be used to check the stabilization of crack opening loads. According to Fig. 6a, the crack opening load corresponding do the numerical simulation with the KH model stabilized after a larger number of crack increments than for the other models. In Fig. 7a it is also possible to observe that, for the numerical simulation with the KH model, the plastic deformation along the crack front also stabilized after a larger number of crack increments. Similar results were obtained for the steel. Fig. 7b presents the variation of crack opening load obtained from the extrapolation of the numerical results by using the Voce equation. It is possible to observe in the figure that, independently of the material model in use, the variation of the consecutive extrapolated values becomes always close to zero, which confirms crack opening load stabilization. The comparison of both graphs of Fig. 7 shows that converged opening stress values are attained after a number of crack increments that coincides with plastic deformation stabilization. Fig. 8a and b present crack opening results obtained for the aluminium alloy and steel, respectively, in the numerical simulations performed with plane strain conditions. It is important to remember that these results were obtained by using a more refined mesh than that used for the numerical simulation of the plane stress cases and, for the aluminium alloy, a different external load was used: r max /r ys = 0.5. For each material, the shape of the crack opening curves in Fig. 8 is similar for almost all the material models in use. With initial crack propagation, a maximum value of the crack opening load is reached in all the curves, and after that, the opening load values start decreasing gradually. The only exception is the KH curve, registered for the aluminium alloy, which continuously increase with crack propagation.
The decrease in opening values after a maximum value is reached, is much more important for the steel (Fig. 8b), than for the aluminium alloy (Fig. 8a), which is explained by the different load level. As Fig. 8b illustrates, it was possible to apply with success the Voce model to the decreasing region of the curves obtained for the steel. The extrapolated values indicate a possible stabilization of the opening load with continuous crack propagation. Previous stabilization analysis can be explained by the formation of a significant plastic wedge during the first load cycle. Since the material ahead of the crack tip, at the beginning of the numerical simulation, was not plastically deformed, the plastic deformation in the forward plastic zone resulting from the first load cycle is much higher than that resulting from the subsequent load cycles. This results from the fact that the material in the forward plastic zone of subsequent load cycles will be already plastically deformed and work-hardened during the first load cycle. As soon as crack propagation moves this initial plastic wedge away, crack tip plastic deformation starts stabilizing. The application of an overload produces a similar effect [8].
For each material, the variation in crack opening load values according to the material models in use follows the same trend as in the plane stress numerical simulations, with the lower values being registered with the KH model. For all the material models in study, the crack opening load values are significantly smaller for the steel, which could be expected considering the lower load level. Like in the plane stress numerical simulations, the higher opening values for the aluminium were registered with the VIH model and similar results were obtained with SIH, SIKH and VIKH models. For the steel, similar results were obtained for both pure isotropic hardening models (SIH and VIH) and for both mixed hardening models (SIKH and VIKH).
As in the previous analysis, the derivatives ð@ e p =@xÞ were plotted in order to study the relation between the plastic deformation along crack front and the opening load stabilization mechanisms. These results are shown in Fig. 9a and b for the aluminium and steel samples, respectively. In both graphs it is possible to observe that none of the derivatives tends to zero, which means that plastic deformation evolution along the crack front did not stabilize during the crack propagation interval considered in the numerical simulation. These results are in accordance with the crack opening load evolution of Fig. 8. However, for the aluminium alloy, that presents almost stabilized values, the plastic strain derivatives approaches zero during some period of crack propagation.

Discussion of results
In Fig. 10 the crack opening stabilized values for both materials and loading conditions are compared. As it is possible to observe in the figure, the stabilized crack opening loads are always higher in plane stress than in plane strain. As already mentioned, the lower opening values were always registered when using the KH model. Differences in crack tip deformation mechanisms according to the external load and material model in use explain the results obtained. As it is well known, main phenomena in cyclic plasticity are forward and reversed plasticity, strain ratcheting, mean stress relaxation and cyclic hardening or softening. Forward plastic deformation depends on the maximum load and is usually related with K max and stress state. The increase of the size of the forward plastic zone and of the magnitude of the plastic deformation inside of it, are expected to increase the opening load. The results obtained in this work, namely those presented in Figs. 3-5, indicate that the strain levels inside the forward plastic zone are higher for the KH model than for all the other models. According to these results, the higher opening results would be also obtained with this mechanical model. However, the occurrence of reversed plasticity also has a significant influence in crack opening behaviour. The reversed plastic zone is formed by the material near the crack tip undergoing compressive yielding at the minimum load. The increase of reversed plasticity is expected to decrease the opening load. In plane stress conditions, the KH model was found to produce the largest reversed plastic zone, while the pure isotropic models (VIH and SIH) were found to produce the smallest reversed zones. This explains the relatively low values of PICC presented in Fig. 6 for the KH model. Simultaneously, for plane strain conditions, results in Fig. 5 suggest that the plastic deformation inside the forward plastic zone is very small and that the size of this zone is almost the same of the reversed plastic zone. These results explain that the crack opening loads in the plane strain tests are much smaller than in the plane stress tests. According to all previous results, the effect of reversed plasticity is dominant over the effect of the forward plastic zone.
Concerning the other phenomena associated with cyclic plasticity, it is expected that strain ratcheting increases plastic deformation at the crack tip, and therefore, it is expected to increase the level of closure. In the same way, mean stress   relaxation increases the residual stresses near crack tip, and therefore, the closure level. In order to analyse the influence of these two phenomena in the observed cyclic behaviour, new plane stress tests were performed in which the crack was propagated 18 increments of 16 lm, with 2 load cycles between increments, and after that, 20 load cycles were applied without crack propagation. Fig. 11 presents, for all the material models in study, the r yy -e yy curves obtained in these tests for a Gauss point placed immediately ahead of last crack tip position. As can be seen in the graph, strain ratchening was registered for all material models, but with particular significance for the KH model. On the other hand, for the pure isotropic hardening models (SIH and VIH), this phenomenon was less important. Additionally, a slight decrease of mean stress and a slight cyclic hardening can be observed for all material models. An additional phenomenon was also observed when, after 18 crack increments, 20 load cycles were applied without propagating the crack. Fig. 12 presents the crack profile in steel samples, before and after the 20 load cycles without crack propagation were applied to it. As it is possible to see in Fig. 12a, where is shown the crack profile obtained in the numerical simulation performed with the SIH model and plane stress conditions, after 20 cycles without crack propagation, a slight crack opening occurs near the crack tip. Actually, it is the first node behind crack tip that moves promoting this phenomenon. This behaviour can be responsible for a slight decrease of crack opening load in plane stress conditions. Similar crack tip opening behaviour was observed with all the material models in use.
For plane strain conditions, as it is possible to see in Fig. 12b, the movement of the first node behind crack tip also occurs and is much more significant than in plane stress conditions. Comparing all plane strain results, it was also observed that the amplitude of this movement strongly depends on the material model used in numerical simulation. To illustrate this aspect, in Fig. 12b are compared the results relative to the SIH and SIKH models. As it is possible to see in the figure, node displacement is much higher when using the SIKH model than with the SIH model. Higher nodal displacements were always obtained with the material models comprising the kinematic hardening component. It is also interesting to enhance that, whatever is the mesh size (8 lm in plane strain and 16 lm in plane stress), only the node immediately behind the crack tip node is moved.

Conclusions
Five mechanical models were used in this work in order to analyse the influence of work-hardening modelling in the numerical prediction of PICC. For both materials considered in the study, noticeable differences in crack tip forward and reversed plastic deformation behaviour and crack opening predictions were depicted when comparing the results obtained with the different mechanical models. The lowest values of crack opening load were found for the pure kinematic models, while the pure isotropic models presented the highest values. Mixed models presented intermediate values, but closer to the isotropic values than to the kinematic values. These variations are intimately related to the near crack tip deformation mechanisms. Reversed plasticity and strain ratcheting were found to be particularly affected by the material model in use. An additional effect, the movement of the node behind crack tip, was also found to occur. The summation of the different effects explains the trends of the numerical results.
Analysing crack opening load stabilization mechanisms, it was found that the variation of equivalent plastic deformation along crack flanks and the variation of extrapolated values can be used to identify the stabilization of crack closure level.