^{1}School of Aerospace Engineering and Applied Mechanics,Tongji University, Shanghai 200092, China

The elastoplastic field near crack tips is investigated through finite element simulation.A refined mesh model near the crack tip is proposed. In the mesh refining area, element size continuously varies from the nanometer scale to themicrometer scale and the millimeter scale. Graphics of the plastic zone, the crack tip blunting, and the deformed crack tip elements are given in the paper.Based on the curves of stress and plastic strain, closely near the crack tip, the stresssingularity index and the stress intensity factor,as well as the plastic strain singularity index and the plastic strain intensity factor are determined.Thestress and plastic strainsingular index vary with the load, while the dimensions of the stress and the plastic strain intensity factorsdependon the stress and the plastic strain singularity index, respectively. The singular field near the elastoplastic crack tip is characterized by the stress singularity index and the stress intensity factor, or alternativelythe plastic strain singularity index and the plastic strain intensityfactor.At the end of the paper, following Irwin’s concept of fracture mechanics,${{\sigma }_{\delta K}}$criterion and${{\varepsilon }_{\delta Q}}$criterion are proposed.Besides, crack tip angle criterion is also presented.

Since the end of World War II, the fracture mechanics has developed steadily and fruitfully, including linear elastic fracture mechanics (LEFM), elastoplastic fracture mechanics (EPFM), and interfacial fracture mechanics (IFFM), etc. Many researchers have worked in this field. Among them, Irwin[1], Rice[2, 3, 4, 5], Huchinson [6, 7], and Wells[8]contributed many important advances. Thereupon, the application of fracture mechanics in structural design and analysis has increased significantly [9].

On account of that thefracture parameters of criterionsforLEFM and IFFM are extracted from the rigorous solutions in linear elasticity[2, 3, 6, 10, 11, 12], the theoretical background issolid. But, thefracture parameters of criterionsforEPFM, such as J integral, crack-tip opening displacement(CTOD) and crack-tipopening angle(CTOA), are based on the simplified solutions in nonlinear plasticity [4, 5, 7, 8], without the supports from the reliable theoretical verifications. Continuingeffortshave been devoted to improving the theory for more than half a century [13].Due to the mathematical difficultyof nonlinearity, up to date, advances in EPFM are limited.

Therefore, peopleswitch their attention to numerical analysis [14].Considering the inexactness of elastic-plastic finite element analysis (EPFEA) of the crack tip field, it is most important to improve the accuracy of EPFEA. Although, someprogressin EPFEAhas beenmade [15], the satisfactory results have not been achieved yet.

Inthe presentpaper, a finite element numerical study for the stress field near elastoplastic crack tips is carried out. To increase the accuracy, a refined mesh model near the crack tip is proposed. In the mesh refining area, element size continuously varies from the nanometer scale (nm) to the micrometer scale (μ m) andthe millimeter scale (mm). Two crack tip elements (1st and 2nd CTE) are specifically designed. The 1st CTE at the left side of the crack tip is an isosceles triangle having a center angle of 120^{o}, and the 2nd CTE at the right side of the crack tip is a right-angled triangle having a center angle of 60^{o}. Graphics of the plastic zone, the crack tip blunting, and the deformed crack tip elements are given in the paper. Based on the curves of stress and plastic strain, closely near the crack tip, the stress singularity index${{\lambda }_{p}}$and the stress intensity factor ${{K}_{p}}$ , as well as the plastic strain singularity index$\gamma _{p}^{{}}$and the plastic strain intensity factor ${{Q}_{p}}$ are determined. It is seen that, the stress and the plastic strain singular index vary with the load, while the dimensionsof the stress and the plastic strain intensity factors depend on the stress and the plastic strain singularity index, respectively. Therefore, the singular field near the elastoplastic crack tip is characterized by the stress singularity index and the stress intensity factor, or alternativelythe plastic strain singularity index and the plastic strain intensityfactor.The stress intensity factor ${{K}_{p}}$ or plastic strain intensity factor ${{Q}_{p}}$ cannot be used alone as a fracture parameter.At the end of the paper, following Irwin’ s concept of fracture mechanics, the ${{\sigma }_{\delta K}}$ criterion andthe ${{\varepsilon }_{\delta Q}}$criterion are proposed.In addition, thecrack tip angle criterion is also presented.

It is noteworthythat in seventies of the last century, researcherspaid significantattention to the EPFEA of finite deformation [16, 17] and the EPFEAof crack problems [18, 19, 20], to investigate EPFM by numerical approaches.Their work laid the foundation for research in computational elastoplastic fracture mechanics [21, 22].

In the following decades, the finite element analysis of large elastoplastic deformation near the crack tip has become the focus of research [23, 24, 25, 26, 27].In order to correctly analyze the elastoplastic crack tip field, finite element modeling with suitable mesh refining near the crack tip is a key problem.Lamain[28] summarized the main points that the modeling needs to follow. He pointed out that six-node triangularisoparameter elements are more suitable for modeling the area containing the crack tip, and the crack tip elements should describe the plastic incompressibility.He also proposed the iteration process for equilibrium correction, the iteration process for yield surface correction, as well as the combination of iteration and substepping to minimize the nonlinear error.

Recently, Ji et al [29] took the full advantages of the modelling capability provided by the finite element software Abaqus/Standard, and achieved a promising numerical result. An element mesh model of a central cracked plate is presented in Fig. 1 [31], and then the stress intensity factor is determined from the finite element analysis. In the mesh refining area, the size of the smallest element is 3.8 nm, and the largest element is 0.82 mm, renderingthe refining ratio of around 216000. Element size continuously varies from the nanometer scale(nm) to the micrometer scale (μ m) and the millimeter scale (mm).

Compared with the rigorous solution, the error of the obtained stress intensity factor is less than 1% in general[31].It is demonstrated that the stress singularity near the crack tip can be simulated accurately through classic FEM.

Therefore, it is feasible to utilize the FEM with the element size continuously varying from the nanometer scale to the micrometer scale and the millimeter scale to solve the elastoplastic crack problems.

Consider an elastoplasticplate containing a central crack under tension, subjectedto finite plastic deformation.Assume the plate size is 200mm, the crack length is $2a$=20 mm, and themaximum applied tensile load is 150MPa.The origin is set at the crack tip, the *x*-axis and *y*-axis areoriented along and perpendicular to the crack, respectively. The material obeys von Mises flow criterion and the associated flow rule.Theelastoplastic crack problem isanalyzed with EPFEM (Abaqus/Standard).In the calculation, material nonlinearity and geometric nonlinearity are considered.

As showed in Fig.2a, b, the mesh refining area is a circulardomain, with the center at the crack tip. The radii of the inner circle, the middle circle and the outer circle are ${{r}_{1}}$=300 nm, ${{r}_{2}}$=1000 nm and ${{r}_{3}}$ =2mm, respectively.

Inside the inner circle, the modified six-node triangularisoparametric elements are adopted.10 divisionsare uniformly allocated along the crack surface, 1division of 60nm and 8divisions of 30nm are allocated along the positive *x*-axis. The area within the inner circle of ${{r}_{1}}$=300 nm, are divided into10 layers radially. The most outer layer is uniformly subdivided into 20 divisions circumferentially.Each inner layer decreases 2 divisions compared to the adjacentouter layer.Two crack tip elements (1st and 2nd CTE) are specially created for capturingthe seriously nonuniform distribution of plastic deformation.Asshowed in Fig.2c, the 1st CTEat the left side of the crack tip is an isosceles triangle having a center angle of 120^{o}, and the 2nd CTEat the right side of the crack tip is aright-angled triangle having a center angle of 60^{o}.

In the annulus between the inner and the middle circles, modified six-node triangularisoparametric elements are adopted. 14 elements are distributed along the radial direction with the element size varying according to the bias ratio of 2.28 (Abaqus index). Along the circumferential direction, elements are uniformly distributed. From the inner to the middle circles, the number of elements along the circumference increases from 20 to 48.

In theannulus between the middle and the outer circles, eight-node quadrilateral isoparametric elements are adopted.112elements are distributed along the radial direction with the element size being proportional to the radius according to the bias ratio of 1869 (Abaqus index), and 48 elements are uniformly distributed in the circumferential direction.

In the refined mesh area, the size of the smallest element is 30nm, and that ofthe largest element is 0.131mm, which gives the refining ratio of around 4367. Element size continuously varies from the nanometer scale (nm) to the micrometer scale (μ m) and the millimeter scale (mm). The finite element mesh for a quarter of the plate contains 7974 elements in total.

Our experience shows that the refined mesh model of the elastoplastic crackedplate, showed in Fig.2, works well in the whole loading process, until the maximum tensile load 150MPa is reached.

In Fig.2, inside of the circle of radius r_{2} =1000 nm, modified six-nodetriangular isoparametric elements (CPS6M for plane stress analysis, CPE6M for plane strain analysis)are adopted. Outside of the circle of radius r_{2}=1000 nm, eight-node quadrilateral isoparametric elements(CPS8 for plane stress analysis, CPE8 for plane strain analysis) are adopted.

Regular second-order triangular elements(CPS6, CPE6)aretypically preferable in Abaqus/Standard.However, these elements may exhibit “ volumetric locking” when incompressibility is encountered, for example, in problemswith large plastic deformation.The modified triangularisoparametric elements exhibit minimal shear andvolumetric locking, andare robust during finite plastic deformation.

The modified six-node triangular isoparametric element and eight-node quadrilateral isoparametric elements are incompatible. But in our case, the results show that the disturbance atthe boundary of these two kinds of elements isnot significant and may be negligible.

For the purpose of simplicity, the material is assumed as linear elastic-linear hardening plastic (LE-LHP). Assume that Young’ s modulus is 198000MPa, Poisson’ s ratio is 0.3, and yield stress is 2000MPa. The tangential plastic modulus is 20000MPa. The input data of plastic properties for the material model are listed in Table1.

The continuity hypothesis in theory of elastoplasticity holds that the constitutive relationship of the material remains unchanged when the volume element approaches zero. Therefore, in the finite element analysis, the finite element of nanometer scale size can be modeled by the conventional plastic medium’ s constitutive relationship. Only in this way, the results from finite element modeling can be determined to be the approximate numerical solution which is consistent with the classical elastic-plastic theory.

Because of the large plastic deformation, the option “ NLGEOM” is set at “ on” in Abaqus/Standard to deal with the geometric nonlinearity. The updated Lagrangian formulation is used when NLGEOM is specified.

The external load is divided into three steps (0~50MPa, 50~100MPa, 100~150MPa). To avoid the divergence in iteration, for each load step, the minimum time increment is set at 0.001 (corresponding 0.05MPa). At the end of each time increment, the nodes are updated to the current (deformed) configurations, and a new stiffness matrix is formed.

Then the force residual for the iteration*R _{a}* is calculated by Abaqus. If

*R*is less than the force residual tolerance at all nodes, the solution as being in equilibrium is accepted. By default, thistolerance value is set to 0.5% of the average force in the structure, which is averaged over time.And the last displacement correction

_{a}*c*is also checked. If

_{a}*c*is less than a fraction(1% by default) of the incremental displacement, the solution is accepted. Bothconvergence checks must be satisfied before a solution is said to have converged for that time increment. Otherwise, another iteration (attempt) would be performed. If the time increment becomes smaller than thedefined minimum value or more than 5 attempts are needed, Abaqus/Standard interrupts the analysis.

_{a}The total computer time to analyze a typical crack problem with a high-performance personal computer is generally less than 60min.

As the calculation process is strictly monitored by the Abaqus/Standard, we did not find any sign of the divergence in the numerical results due to the vast span ofelement size from nanometer scale to the millimeter scale.

Figures3and 4showthe plastic zone of the central cracked plate at different loading levels in the undeformed configuration, in plane strain and plane stress, respectively. The initial mesh is also plotted in the figures.

In plane stress, the plastic zone is extended along the ligament of the cracked plane. It is noted that, in plane strain, the plastic zone is confined to a shorter range along the *x*-axis. But it does not mean that the plastic zone does not extend, it is just extending along theoblique direction of about ± 70^{o}.

The plastic zone radius${{r}_{p}}$of the plane strain and the plane stress, as well asthe plastic zone radius along the*x*-axis of the plane strain$r_{p}^{x}$, under different loads, are given in Fig.5.

Figure 5 showsthat the radius of the plastic zone near the crack tip increases faster than the load as the plastic deformation dominated. The plastic zone radius in plane strain is much smaller than that in plane stress under the same load.

Figures6 and 7 illuminate the crack tip blunting of the central cracked plateunder different loading levels, in plane strain and plane stress, respectively.Blunting greatly alters the shape of the crack surface near the crack tip.

Figures8 and 9describe the shapes of the crack tip elements andtheir neighboring elements subjected to plastic deformation under different loading levels, in plane strain and plane stress, respectively. A comparison of the deformed crack tip element with its initial shapedemonstrates the complication of the blunting phenomena.

The plastic deformation process of the 1st CTE has unique features. As loading increases from 0 to 150MPa, it is seen that thebottom edge (crack flank) rotates from horizontal to a certain angle ($\beta $) and elongated in*y*-direction.

It is especially interesting to draw attention to the stress and strain state of the 2nd CTE. Since the bottom edge of the 2nd CTE is located at thestarting portionof the ligament, the extension of the crack will likely start at this place.

Duringloading, the upper and lower crack flanks open to a concave angle $2\beta $ ($0< 2\beta < \text{ }\!\!\pi\!\!\text{ }$). It is seen that a rounded end does not appear at the blunted crack tip. The blunted crack tip is just a sharp vertex of a concave angle $2\beta $. Sharp vertex of a concave angle is an origin of stress singularity in elastoplastic field.

Angle $\beta $may be regarded as the “ crack tip angle (CTA)” .Crack tip angle$\beta $of the plane strain and plane stress under different loads areplotted in Fig. 10.

5.4.1 Plane strain

Figures11 and 12 describe the stress components(${{\sigma }_{y}}, {{\sigma }_{x}}, {{\tau }_{xy}}, {{\sigma }_{z}}$)and the plastic strain components ($\varepsilon _{y}^{p}, \varepsilon _{x}^{p}, \gamma _{xy}^{p}, \varepsilon _{z}^{p}$), respectively, along the ligament near the crack tip, under different loading levels in plane strain.

It is seen that, in plane strain case, the material element at the crack tip is in the state of tri-axial tensile stress, which could account for the plastic zone radius in plane strain is much smaller than that the plastic zone radius in plane stress under the same load, and why the plane strain more tends to brittle fracture.

5.4.2 Plane stress

Figures 13 and 14 describe the stress components (${{\sigma }_{y}}, {{\sigma }_{x}}, {{\tau }_{xy}}$) and the plastic strain components ($\varepsilon _{y}^{p}, \varepsilon _{x}^{p}, \gamma _{xy}^{p}$), respectively, along the ligament near the crack tip, under different loading levels in plane stress.

The curves${{\sigma }_{y}}\tilde{\ }r$in Figs.11a and 13a from point A to point Bare replotted onto the log-log coordinates in Fig.15a, b, respectively. Points A and B are labelled in Figs.11a and 13a, as ${{r}_{A}}$ and ${{r}_{B}}$ respectively.It is seen that the lines in Fig. 15a, b are nearly straight. That means an asymptotic expression existing near the crack tip,

${{\sigma }_{y}}={{K}_{p}}{{r}^{-{{\lambda }_{p}}}}, \text{ }r\to 0, \text{ }$ (1)

Where ${{\lambda }_{p}}$ is the stress singularity index, and ${{K}_{p}}$ is the stress intensity factor. The stress singularity index and stress intensity factor for plane strain and plane stress are listed in Table 2.

The curves$\varepsilon _{y}^{p}\tilde{\ }r$in Figs.12a and 14a from point A to point B(for the curve of 50 MPa in Fig.12a, point B is instead by point B')are replotted onto the log-log coordinates in Fig. 16a, b, respectively. Points A and B are labelled in Figs. 12a and 14a as ${{r}_{A}}$and ${{r}_{B}}$, respectively. Point B' in Fig. 12a for the curve of 50 MPa is labelled as${{r}_{{{B}'}}}$.It is seen that the lines in Fig. 16a, b are nearly straight. That means an asymptotic expression existing near the crack tip.

$\varepsilon _{y}^{p}={{Q}_{p}}{{r}^{-{{\gamma }_{p}}}}, \text{ }r\to 0, $ (2)

aPlane strain.bPlane stress

Where${{\gamma }_{p}}$istheplastic strain singularity index, and ${{Q}_{p}}$ is the plastic strain intensity factor. The plastic strain singularity index and plastic strain intensity factor for plane strain and plane stress are listed in Table 3.

The resultsgivenin Tables 2 and 3show that the singularity near the elastoplastic crack tip is much more complicated than the singularity near the linear elastic crack tip. Thestresssingular index${{\lambda }_{p}}$andplastic strain singular index$\gamma _{p}^{{}}$vary with the load, while the dimensions of the stress intensity factor ${{K}_{p}}$ and plasticstrain intensity factor ${{Q}_{p}}$ are related to thesingularity index${{\lambda }_{p}}$and$\gamma _{p}^{{}}$, respectively.Moreover, the magnitude of the stress intensity factor and plastic strain intensity factor do not vary monotonically with the loading level. Therefore, the strength of the elastoplastic singular field near the crack tipcannotbe characterized by the intensity factor ${{K}_{p}}$ or ${{Q}_{p}}$ separately.Equations (1) and(2)show that the singular field near the elastoplastic crack tip is characterized by two coupling parameters:${{\lambda }_{p}}$and ${{K}_{p}}$ , oralternatively${{\gamma }_{p}}$and ${{Q}_{p}}$ .

Irwinestablished the SIF-based criterion in LEFM[1].He found that the singular stress field near the crack tip is described by the asymptotic expression [10, 11], i.e.

${{\sigma }_{y}}=K{{r}^{-\lambda }}, \text{ }\lambda \text{=-}\frac{1}{2}\text{, }r\to 0.\text{ }$ (3)

Thus, the singular stress field in the vicinity of the crack tip is wholly determined by the stress intensity factor *K* and the singularity index$\lambda $. As$\lambda $ is a constant, the singular stress field near the crack tip is characterized by the stress intensity factor*K*. Therefore, the fracture criterion is written as

$K={{K}_{cr}}.$ (4)

where*K _{cr}*is the critical value of

*K*.

The asymptotic expression may be further interpreted. At the crack tip$(r\to 0)$, ${{\sigma }_{y}}\to \infty $. Let $\delta $ be an arbitrary specified distance from the crack tip, ${{\sigma }_{\delta K}}$ be the stress at$\delta $,

${{\sigma }_{\delta K}}=K{{\delta }^{-\lambda }}.$ (5)

Then, the infinite stress at the crack tip and the finite stress atthe point$\delta $both are correspondingto the same *K*and$\lambda $.The infinite stress at the crack tip cannot be used as a fracture parameter, but the finite stress${{\sigma }_{\delta K}}$ can be used as a fracture parameter. Therefore, the fracture criterion also may be written as,

${{\sigma }_{\delta K}}={{({{\sigma }_{\delta K}})}_{cr}}.$ (6)

where (*σ _{δ K}*)

*is the critical value of*

_{cr}*σ*.

_{δ K}The criterion (4) and (6) are equivalent. In LEFM, only the stress intensity factor (SIF)-based criterion (4) has been widely accepted.

*β*criterion in EPFM

Following Irwin’ s criterion concept of fracture mechanics, the criterion in EPFM can be established.

Let $\delta $be an arbitrary specified distance from the crack tip, ${{r}_{A}}< \delta < {{r}_{B}}$(Figs.12a and14a). Define${{\sigma }_{\delta K}}$as the stress at point$\delta $, and${{\varepsilon }_{\delta Q}}$as the plastic strain at point$\delta $.${{\sigma }_{\delta K}}$and${{\varepsilon }_{\delta Q}}$are determined by the asymptotic expressions (1) and (2), respectively. Therefore, the applied load, (${{\lambda }_{p}}$, ${{K}_{p}}$ ) and (${{\gamma }_{p}}$, ${{Q}_{p}}$ ), ($\delta $, ${{\sigma }_{\delta K}}$)and ($\delta $, ${{\varepsilon }_{\delta Q}}$) reach their respective critical value at the same time. That means ($\delta $, ${{\sigma }_{\delta K}}$) and ($\delta $, ${{\varepsilon }_{\delta Q}}$) may be adopted as a fracture parameter.Taking into accountthat$\delta $is a prescribed fixed constant, the fracture criterion may be written as

$({{\sigma }_{\delta K}})={{({{\sigma }_{\delta K}})}_{cr}}$ , (7)

$({{\varepsilon }_{\delta Q}})={{({{\varepsilon }_{\delta Q}})}_{cr}}$ . (8)

Therefore, Irwin’ s criterion concepthasbeen generalized to the ${{\sigma }_{\delta K}}$ criterion and ${{\varepsilon }_{\delta Q}}$ criterion in EPFM.The${{\sigma }_{\delta K}}$ criterion is suitable for thematerials of high strain hardening.On the contrary, the${{\varepsilon }_{\delta Q}}$ criterion is suitable for thematerials of low strain hardening, including ideal plastic materials.

Besides, CTA can also be adopted as a fracture parameter. Such that, CTA criterion may be written as

$\beta ={{\beta }_{cr}}$ (9)

In criterions (7) -(9), ${{({{\sigma }_{\delta K}})}_{cr}}$, ${{({{\varepsilon }_{\delta Q}})}_{cr}}$and ${{\beta }_{cr}}$ are thecriticalvaluesof${{\sigma }_{\delta K}}$, ${{\varepsilon }_{\delta Q}}$and $\delta $, respectively, which controllelastoplasticcrack extension, and should be determined from fracture experiments ofplane strain and plane stress.

${{\sigma }_{\delta K}}$criterion and${{\varepsilon }_{\delta Q}}$ criterion in EPFM are valid under the condition of $\delta $being an arbitraryspecified constant.Once$\delta $is fixed, it cannot be changed again. Therefore, these criterions are likely to be transient ones.A more reliable version of the criterion could be developed after more efforts are paid.

On the other hand, because the definition of CTA is more precise than that of CTOA, CTA criterionis an improved version of CTOA criterion.

Furthermore, the above proposed criterions have anobvious advantage. They are applicable to small scale yieldingand medium scale yielding, even to large scale yielding.

The work was supportedbythe National Natural Science Foundation of China (Grant 11572226).

**Reference**

1. |
Analysis of stresses and strains near the end of a crack traversing aplate |

2. |
Plane problems of cracks in dissimilar media |

3. |
Elastic fracture mechanics concepts for interfacialcracks |

4. |
A path independent integral and the approximate analysis of strainconcentration by notches and cracks |

5. |
Plane Strain Deformation near a Crack Tip in aPower-Law Hardening Material |

6. |
Mixed mode cracking in layered materials |

7. |
Singular behavior at the end of a tensile crack tip in a hardening material |

8. |
Application of fracture mechanics at and beyond general yielding |

9. |
Fracture mechanics - fundamentals and applications, 3rd ed. CRC Press, Boca Raton |

10. |
Bearing pressures and cracks |

11. |
On the stress distribution at the base of a stationary crack |

12. |
SIF-based fracture criterion for interface cracks |

13. |
Review of fracture toughness (G, K, J, CTOD, CTOA) testing and stand ardization |

14. |
Comparison of finite element solutions with analytical and experimental data for elastic-plastic cracked problems |

15. |
A calculational round robin in elastic-plastic fracture mechanics. |

16. |
A finite element formulation for problems of large strain and large displacement |

17. |
Finite-element formulations for problems of large elastic-plastic deformation |

18. |
Small scale yielding near a crack in plane strain: a finite element analysis |

19. |
The role of large crack tip geometry changes in plane strain fracture |

20. |
Recent finite element studies in plasticity and fracture mechanics |

21. |
Elastic-Plastic Fracture Mechanics |

22. |
Computational fracture mechanics |

23. |
Elastic-plastic fracture analysis of finite bodies—I. Description of the stress field |

24. |
Constraint effects on the elastic-plastic fracture behaviour in strain gradient solids |

25. |
Elastic-plastic finite element analysis for double-edge cracked tension (DE(T)) plates. Eng |

26. |
Modelling of crack tip blunting using Finite Element Method |

27. |
Elastic-plastic finite element analysis of the effect of compressive loading on crack tip parameters and its impact on fatigue crack propagation rate |

28. |
Numerical analysis in EPFM |

29. |
Determination of stress intensity factor with direct stress approach using finite element analysis |