Finite element simulation of elastoplastic field near crack tips and results for a central cracked plate of LE?LHP material under tension

X. Ji · F. Zhu

Abstract 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 the micrometer 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 stress singularity index and the stress intensity factor, as well as the plastic strain singularity index and the plastic strain intensity factor are determined. The stress and plastic strain singular index vary with the load, while the dimensions of the stress and the plastic strain intensity factors depend on 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 alternatively the plastic strain singularity index and the plastic strain intensity factor. At the end of the paper, following Irwin’s concept of fracture mechanics,ò Kcriterion and σ Qcriterion are proposed. Besides, crack tip angle criterion is also presented.

Keywords Crack · Fracture mechanics · Elastoplasticity · Finite element method · Geometry nonlinearity* X. Ji jixing@#edu.cn

1 Introduction

Since the end of World War II, fracture mechanics has developed steadily and fruitfully, including linear elastic fracture mechanics (LEFM), elastoplastic fracture mechanics(EPFM), and interfacial fracture mechanics (IFFM). Many researchers have worked in this field. Among the m, Irwin[1], Rice et al. [2–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].

Because the fracture parameters of criteria for LEFM and IFFM are extracted from the rigorous solutions in linear elasticity [2, 3, 6, 10–12], the the oretical background is solid. But, the fracture parameters of criteria for EPFM, such as J integral, crack-tip opening displacement (CTOD), and crack-tip opening angle (CTOA) are based on the simplified solutions in nonlinear plasticity [4, 5, 7, 8], without the supports from the reliable the oretical verifications. Continuing efforts have been devoted to improving the theory for more than half a century [13]. Because of the mathematical difficulty of nonlinearity, up to now, advances in EPFM are limited.

Therefore, the attention has been switched 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 some progress in EPFEA has been made [15], satisfactory results have not been achieved yet.

In the present paper, 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), and the millimeter scale (mm).Two crack tip elements (CTE) (1st and 2nd) are specifically designed. The 1st CTE at the left side of the crack tip is an isosceles triangle having a center angle of 120°, and the 2nd CTE at the right side of the crack tip is a right-angled triangle having a center angle of 60°. 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 λpand the stress intensity factor Kp, as well as the plastic strain singularity indexpand the plastic strain intensity factor Qpare determined. It is seen that the stress and the plastic strain singular index vary with the load, while the dimensions of 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 alternatively the plastic strain singularity index and the plastic strain intensity factor. The stress intensity factor Kpor plastic strain intensity factor Qpcannot be used alone as a fracture parameter.At the end of the paper, following Irwin’s concept of fracture mechanics, the ò Kcriterion and the Qcriterion are proposed. In addition, the crack tip angle criterion is also presented.

2 Overview of elastoplastic finite element analysis on crack tip fields

It is noteworthy that in the 1970s, researchers paid significant attention to the EPFEA of finite deformation [16, 17] and the EPFEA of crack problems [18–20], to investigate EPFM by numerical approaches. The ir 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–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.Work by Lamain [28] summarized the main points that the modeling needs to follow. The y pointed out that six-node triangular isoparameter elements are more suitable for modeling the area containing the crack tip, and the crack tip elements should describe the plastic incompressibility. The y 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.

3 Finite element analysis o f elastic p late containing a central crack

Recently, Ji et al. [29] took full advantage 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 [29], and the n 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, rendering the refining ratio of around 216,000. 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[29]. It is demonstrated that the stress singularity near the crack tip can be simulated accurately through classic finite element method (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.

Fig. 1 a Finite element mesh of a quarter plate. b Refined mesh model near the crack tip [29]

4 Finite element simulation of elastoplastic plate containing a crack

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

4.1 Mesh refinement near the crack tip

As shown in Fig. 2a, b, the mesh refining area is a circular domain, with the center at the crack tip. The radii of the inner circle, the middle circle and the outer circle are r1= 300 nm, r2= 1000 nm, and r3= 2 mm, respectively.?

Fig. 2 a Finite element mesh of a quarter of an elastoplastic plate. b Refined mesh model near the elastoplastic crack tip. c Crack tip elements

Inside the inner circle, the modified six-node triangular isoparametric elements are adopted. Ten divisions are uniformly allocated along the crack surface, one division of 60 nm and eight divisions of 30 nm are allocated along the positive x-axis. The area within the inner circle of r1= 300 nm, are divided into 10 layers radially. The most outer layer is uniformly subdivided into 20 divisions circumferentially. Each inner layer decreases two divisions compared to the adjacent outer layer. Two crack tip elements (1st and 2nd CTE) are specially created for capturing the seriously nonuniform distribution of plastic deformation. As showed in Fig. 2c, the 1st CTE at the left side of the crack tip is an isosceles triangle having a center angle of 120°, and the 2nd CTE at the right side of the crack tip is a right-angled triangle with a center angle of 60°.

In the annulus between the inner and the middle circles,modified six-node triangular isoparametric elements are adopted. Fourteen elements are distributed along the radial direction with the element size varying according to the bias ratio of 2.28 (Abaqus index). A long 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 the annulus between the middle and the outer circles,eight-node quadrilateral isoparametric elements are adopted.One hundred and twelve elements 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 30 nm, and that of the largest element is 0.131 mm,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 cracked plate, showed in Fig. 2, works well in the whole loading process, until the maximum tensile load 150 MPa is reached.

4.2 Element type

In Fig. 2, inside of the circle of radius r2= 1000 nm, modified six-node triangular isoparametric elements (CPS6 M for plane stress analysis, CPE6 M for plane strain analysis) are adopted. Outside of the circle of radius r2= 1000 nm, eightnode quadrilateral isoparametric elements (CPS8 for plane stress analysis, CPE8 for plane strain analysis) are adopted.

Regular second-order triangular elements (CPS6, CPE6)are typically preferable in Abaqus/Standard. However, these elements may exhibit “volumetric locking” when incompressibility is encountered, for example, in problems with large plastic deformation. The modified triangular isoparametric elements exhibit minimal shear and volumetric locking, and are 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 at the boundary of these two kinds of elements is not significant and may be negligible.

4.3 Material property

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

The continuity hypothe sis 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 can the results from finite element modeling be determined to be the approximate numerical solution, which is consistent with the classical elastic–plastic theory.

4.4 Nonlinear analysis

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–50 MPa,50–100 MPa, 100–150 MPa). To avoid the divergence in iteration, for each load step, the minimum time increment is set at 0.001 (corresponding 0.05 MPa). At the end of each time increment, the nodes are updated to the current(deformed) configurations, and a new stiffness matrix is formed.

Table 1 Plastic property of the material

The n the force residual for the iteration Rais calculated by Abaqus. If Rais less than the force residual tolerance at all nodes, the solution being in equilibrium is accepted.By default, this tolerance value is set to 0.5% of the average force in the structure, which is averaged over time.The last displacement correction cais also checked. If cais less than a fraction (1% by default) of the incremental displacement, the solution is accepted. Both convergence 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 the defined minimum value or more than five attempts are needed, Abaqus/Standard interrupts the analysis.

The total computer time to analyze a typical crack problem with a high-performance personal computer is generally less than 60 m in.

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 of element size from nanometer scale to the millimeter scale.

5 Results from finite element analysis

5.1 Plastic zone

Figures 3 and 4 show the 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 the oblique direction of about ± 70°.

Fig. 3 Plastic zone in plane strain. a Load: 10–40 MPa. b Load: 50–90 MPa. c Load: 100–150 MPa

The plastic zone radius rpof the plane strain and the plane stress, as well as the plastic zone radius along the x-axis of the plane strainunder different loads, are given in Fig. 5.

Fig. 4 Plastic zone in plane stress. a Load: 10–40 MPa. b Load: 50–90 MPa. c Load: 100–150 MPa

Fig. 5 Plastic radius under different load

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

5.2 Blunting of the crack tip

Figures 6 and 7 illuminate the crack tip blunting of the central cracked p late under different loading levels, in plane strain and plane stress, respectively. Blunting greatly alters the shape of the crack surface near the crack tip.

Fig. 6 Blunting of the crack tip in plane strain

Fig. 7 Blunting of the crack tip in plane stress

5.3 Deformed crack tip elements

Figures 8 and 9 describe the shapes of the crack tip elements and their neighboring elements subjected to plastic deformation under different loading levels, in plane strain and plane stress, respectively. In Figs. 8 and 9, equivalent plastic strain (PEEQ) is expressed with color. A comparison of the deformed crack tip element with its initial shape demonstrates the complication of the blunting phenomena.

Fig. 8 Deformed crack tip elements in plane strain

Fig. 9 Deformed crack tip elements in plane stress

Fig. 10 CTA under different load

The plastic deformation process of the 1st CTE has unique features. As loading increases from 0 to 150 MPa,it is seen that the bottom edge (crack flank) rotates from horizontal to a certain angle () 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 the starting portion of the ligament,the extension of the crack will likely start at this place.

During loading, the upper and lower crack flanks open to a concave angle2(0<2