Elasto-plastic fatigue crack growth analysis of plane problems in the presence of flaws using XFEM

Abstract

In this paper, elasto-plastic XFEM simulations have been performed to evaluate the fatigue life of plane crack problems in the presence of various defects. The stress-strain response of the material is modeled by Ramberg-Osgood equation. The von-Mises failure criterion has been used with isotropic hardening. The J-integral for two fracture modes (mode-I and mode-II) is obtained by decomposing the displacement and stress fields into their symmetric and antisymmetric parts, then individual stress intensity factors are extracted from J-integral. The fatigue life obtained by EPFM is found quite close to that obtained by LEFM.

Keywords

XFEM / von-Mises yield criterion / isotropic hardening / fatigue crack growth / J-integral

Introduction

The analysis of static and growing cracks becomes quite important for the health monitoring of the structures/components under fatigue loading. In general, many components e.g. the engine of an aircraft, turbine blades and rotating shafts subjected to cyclic loads may sometimes lead to the fatigue failure. Thus, the design analysis of structures subjected to fatigue load needs to be performed accurately. The fatigue life of a component is highly influenced by the crack closure phenomenon, i.e., the crack surfaces contact during cyclic loading. Several crack closure phenomena like plasticity induced crack closure (PICC), oxide induced crack closure (OICC) and roughness induced crack closure (RICC) are observed in fatigue analysis. The plasticity induced crack closure phenomenon investigated by Elber [ 1] has a major impact on the fatigue crack growth. Since then, several crack closure mechanisms have been identified, but the primary mechanism for the crack closure is the plasticity ahead of the crack tip. During cyclic loading, a significant amount of plastic strains develops near the crack tip upon loading, which is not fully reversible during unloading. This residual plastic strains leads to the formation of plastic wake behind the crack tip as the crack extends. This phenomenon develops compressive stresses at the crack tip, which tends to close the crack tip as illustrated in Fig. 1. Therefore, for the next loading cycle, some applied tensile load consumed to overcome the compressive stresses at the crack tip. This reduces the crack driving force for the further fatigue crack growth.

Since, the introduction of the crack closure concept, a lot of efforts have been made to characterize and predict the crack closure effects on the crack growth rate. Experimental studies [ 2− 4] have played a major role to understand the phenomenon of crack closure in the cyclic loading. Several methods have been proposed to study the effect of crack closure. An analytical crack closure model for center cracked specimen was proposed by Newman [ 5]. Budiansky and Hutchinson [ 6] proposed an analytical model based on elastic-perfectly plastic Dugdale-Barenblatt crack model to study the effect of crack closure. Many investigators have performed fatigue crack growth with plasticity induced crack closure using finite element analysis for both plane stress and plane strain cases. Newman [ 7] was the first who investigated the effect of crack closure phenomenon in center cracked specimen using finite element method (FEM) under plane stress condition. A combined experimental and numerical study of crack closure behavior in Al 2024-T3 was performed by Blom and Holm [ 8]. McClung et al. [ 9] numerically performed the study of plastic zone size and crack closure behavior for the middle cracked tension specimen under plane strain and plane stress conditions. Solanki and his coworkers [ 10] numerically investigated the effect of plasticity induced crack closure for different types of specimens. They observed that the plasticity induced crack closure is negligible for compact tension specimen under plane strain condition. Alizadeh et al. [ 11] performed plane stress and plane strain finite element analysis for two and three dimensional center cracked plate. They observed that the results obtained with crack closure effect were found in good agreement with the plane stress strip yield results. A three-dimensional finite element analysis was carried out by the Ellyin and Ozah [ 12] to analyze the PICC phenomenon in a cracked plate under variable amplitude loading. Toribio and Kharin [ 13] numerically investigated the influence of load range, load ratio and overload over the near tip stress and strain fields under mode-I loading.

Although, a lot of work has been done on fatigue crack growth analysis but it still remains a challenging issue for the scientific community. So far, the finite element method [ 14] has been widely used to simulate the crack growth problems but it not only requires a conformal mesh to model the crack problems but also needs remeshing at each stage of crack propagation along with special elements to handle asymptotic crack tip stresses. To overcome these issues, various numerical methods have been proposed such as boundary element method [ 15, 16], meshfree methods [ 17− 25], extended finite element method (XFEM) [ 26− 28], cracking particle method [ 29− 31], phantom node method [ 32], strain smoothing finite element method [ 33, 34], coarse-graining technique [ 35], coupled extended finite element and meshfree method [ 36], extended isogeometric analysis [ 37], virtual node XFEM [ 38].

The XFEM allows the mesh independent crack modeling, and avoids remeshing during crack growth. To model a crack in XFEM, standard FE approximation is enriched with additional functions [ 26], which are obtained from the theoretical background of the problem. The level set method [ 27, 28] is advantageously used in XFEM to track a moving discontinuity. Several LEFM problems were analyzed by the XFEM such as elastic fatigue crack growth [ 28, 39], crack growth with friction [ 40], arbitrary three dimensional fatigue crack growth [ 41, 42] and dynamic crack growth [ 43, 44]. In 2011, Natarajan et al. [ 45] used XFEM for the simulation of cracked FGM plates, and computed the natural frequencies and mode shapes of simply supported and clamped plates. Baiz et al. [ 46] employed XFEM to study the linear buckling of isotropic plates using quadrilateral elements. They evaluated the buckling coefficient and mode shapes of rectangular and square plates as a functions of crack length, location and plate thickness. Bac et al. [ 47] combined node-based smoothed finite element method (NS-FEM) with XFEM to form the node-based extended finite element method (NS-XFEM) for simulating 2-D fracture problems. In 2012, an edge-based strain smoothing technique was proposed by Chen et al. [ 48] in conjunction with XFEM for the simulation of crack growth problems. The XFEM has also been applied in many areas including large strain analysis (Legrain et al. [ 49],), phase transformation (Ji et al. [ 50],) and plasticity [ 51, 52].

In general, all materials show some kind of heterogeneity at some scale. Therefore, it is essential to consider the effect of the heterogeneity on the life of structure. In 2012, Singh et al. [ 39] performed the linear elastic crack growth analysis of center cracked plate in the presence of flaws (holes, inclusions and minor cracks) using XFEM under monotonic loading. An adaptive multiscale approach was used for the simulation of crack growth problems [ 53]. They employed phantom node method to model the continuum region, while the discontinuous region was modeled with molecular statics. Further, this adaptive mutliscale approach was extended in conjunction with meshless methods to model the fracture problems [ 54]. In 2014, Nanthakumar et al. [ 55] proposed an iterative procedure to solve the inverse problem of multiple voids in piezoelectric structures. They employed XFEM for determining the response of the structures. The results show that this approach was effectively used to determine the number of voids, locations and their shapes. A strain energy based homogenization approach was employed by Kumar et al. [ 56, 57] to model the response of edge crack plate in the presence of multiple defects under mode-I cyclic loading. The results showed that the homogenization approach significantly increases the computational efficiency along with good accuracy. A stochastic multiscale approach was proposed by Bac et al. [ 58] to quantify the correlated key-input parameters influencing the mechanical properties of polymer nanocomposites. They employed a hierarchical multiscale approach to find out the variations of parameters at nano-, micro-, meso- and macro-scales. Further, Bac et al. [ 59] performed a sensitivity analysis to quantify the key-input parameters influencing the mechanical properties of polymeric nanocomposites.

The linear elastic crack growth analysis of center cracked plate in the presence of flaws (holes, inclusions and minor cracks) is already performed by authors using XFEM under monotonic loading [ 39]. In the present work, XFEM has been extended to perform the nonlinear fatigue crack growth analysis of the edge crack plate in the presence of multiple defects. The key points of PICC modeling used in this work are given as,

The concept of contact forces is used to model the PICC and its effect is observed on the fatigue life of the cracked plate.

The effect of multiple flaws such as holes, inclusions and minor cracks is evaluated on the fatigue life of the cracked plate.

The J-integral values for mode-I and mode-II are obtained by decomposing the displacement and stress fields into symmetric and antisymmetric parts.

The crack under cyclic load remains closed in some part of the loading cycle, which results in an opening of stress intensity factor.

Various plane crack problems are solved in the presence of randomly distributed discontinuities/flaws of arbitrary sizes. The contact between the crack faces during unloading is modeled by fixing the nodal displacements. In this work, small strain plasticity is assumed; hence the size of plastic zone ahead of the crack tip is small. Therefore, the values of equivalent stress intensity factor are obtained from the J-integral using linear elastic relations only. The direction of crack growth is determined by using maximum principal stress criterion [ 19, 39]. The fatigue life obtained by XFEM using linear elastic and elasto-plastic theories are compared with each other.

Numerical formulation

Plasticity modeling

The mathematical theory of plasticity provides a theoretical description of the relationship between stress and strain for the elasto-plastic material. Before the onset of plastic yielding, the relationship between stress and strain is given by the standard linear elastic expression as [ 60],

(1) {\sigma }_{ij}=C_{ijkl}{\epsilon }_{kl},

where {\sigma }_{ij} and {\epsilon }_{kl} are stress and strain components respectively and C_{ijkl} is elastic constants tensor.

von-Mises yield criterion has been used to determine the stress level at which plastic deformation begins. The yield condition can be written as,

(2) F({\sigma }_{ij},k)=\mathrm{0},

where F is yield function and k is hardening parameter. The relation between material and hardening parameter defines the progressive development of yield surface.

Elasto-plastic stress-strain relations

The stress-strain relations are derived assuming the isotropic hardening. The change in strain is assumed to be divisible into elastic and plastic components for any increment of stress such that,

(3) \mathrm{d}\mathbf{\epsilon }=\mathrm{d}{\mathbf{\epsilon }}^e+\mathrm{d}{\mathbf{\epsilon }}^p.

From the constitutive law, the elastic component of incremental strain can be written as,

(4) \mathrm{d}{\mathbf{\epsilon }}^e={[\mathbf{D}]}^{-\mathrm{1}}\mathrm{d}\mathbf{\sigma },

where \mathbf{D} is the elastic constitutive matrix which is symmetric, hence {\mathbf{D}}^{\mathrm{T}}=\mathbf{D} .

To obtain the relationship between the plastic incremental strain component and stress increment, a further assumption on the material behavior must be made. Thus, the plastic strain increment is proportional to the stress gradient of a quantity termed as plastic potential function Q, and is written as,

(5) \mathrm{d}{\epsilon }^p=\mathrm{d}\lambda {\displaystyle \frac{\partial Q}{\partial \sigma }},

where \mathrm{d}\lambda is a proportionality constant and termed as plastic multiplier. For associated theory of plasticity F\equiv Q , thus Eq. (5) can be written as,

(6) \mathrm{d}{\mathbf{\epsilon }}^p=\mathrm{d}\lambda {\displaystyle \frac{\partial F}{\partial \mathbf{\sigma }}}.

Governing equations

The governing equilibrium equation in elasto-statics [ 28, 61] is defined as,

(7) \nabla .\mathbf{\sigma }+\mathbf{b}=\mathbf{0} \mathrm{i}\mathrm{n} \mathrm{\Omega },

where \mathbf{\sigma } is Cauchy stress tensor and \mathbf{b} is body force per unit volume.

The yield function in Eq. (2) must lie on the yield surface and hence the following consistency conditions must be hold,

(8) \mathrm{d}F={\displaystyle \frac{\partial F}{\partial \mathbf{\sigma }}}\mathrm{d}\mathbf{\sigma }+{\displaystyle \frac{\partial F}{\partial k}}\mathrm{d}k=\mathrm{0.}

The above Eq. (8) can be written as,

(9) {\mathbf{a}}^{\mathrm{T}}\mathrm{d}\mathbf{\sigma }-A\mathrm{d}\lambda =\mathrm{0},

where,

(10) {\mathbf{a}}^{\mathrm{T}}=\left({\displaystyle \frac{\partial F}{\partial \sigma }}\right)=\left[{\displaystyle \frac{\partial F}{\partial {\sigma }_{xx}}},{\displaystyle \frac{\partial F}{\partial {\sigma }_{yy}}},{\displaystyle \frac{\partial F}{\partial {\sigma }_{zz}}},{\displaystyle \frac{\partial F}{\partial {\sigma }_{yz}}},{\displaystyle \frac{\partial F}{\partial {\sigma }_{zx}}},{\displaystyle \frac{\partial F}{\partial {\sigma }_{xy}}}\right] \mathrm{a}\mathrm{n}\mathrm{d} A=-{\displaystyle \frac{\mathrm{1}}{\mathrm{d}\lambda }}{\displaystyle \frac{\partial F}{\partial k}}\mathrm{d}k.

The vector \mathbf{a} is termed as the flow vector. Using Eq. (4) and Eq. (6), the Eq. (3) can be written as,

(11) \mathrm{d}\mathbf{\epsilon }={[\mathbf{D}]}^{-\mathrm{1}}\mathrm{d}\mathbf{\sigma }+\mathrm{d}\lambda {\displaystyle \frac{\partial F}{\partial \mathbf{\sigma }}}.

After pre-multiplying both sides of Eq. (11) by {\mathbf{a}}^{\mathrm{T}}\mathbf{D} and eliminating {\mathbf{a}}^{\mathrm{T}}\mathrm{d}\mathbf{\sigma } using Eq. (9), we obtain the plastic multiplier \mathrm{d}\mathit{\lambda } as,

(12) \mathrm{d}\mathit{\lambda }={\displaystyle \frac{\mathrm{1}}{[A+{\mathbf{a}}^{\mathrm{T}}\mathbf{D}\mathbf{a}]}}{\mathbf{a}}^{\mathrm{T}}\mathbf{D}\mathrm{d}\mathbf{\epsilon }.

The stress-strain relations are obtained by substituting Eq. (10) and Eq. (12) into Eq. (11) as,

(13) \mathrm{d}\mathbf{\sigma }=\left(\mathbf{D}-{\displaystyle \frac{\mathbf{D}\mathbf{a} {\mathbf{a}}^{\mathrm{T}}\mathbf{D}}{[A+{\mathbf{a}}^{\mathrm{T}}\mathbf{D}\mathbf{a}]}}\right)\mathrm{d}\mathbf{\epsilon },

(14) \mathrm{d}\mathbf{\sigma }={\mathbf{D}}_{ep}\mathrm{d}\mathbf{\epsilon },

where {\mathbf{D}}_{ep}=\mathbf{D}-{\displaystyle \frac{{\mathbf{d}}_D{\mathbf{d}}_D^{\mathrm{T}}}{A+{\mathbf{d}}_D^{\mathrm{T}}\mathbf{a}}} and {\mathbf{d}}_D^{\mathrm{T}}={\mathbf{a}}^{\mathrm{T}}\mathbf{D} .

As stated above, for the associated theory [ 62] of plasticity, F\equiv Q i.e., yield function and potential function are identical, thus the elasto-plastic constitutive matrix {\mathbf{D}}_{ep} becomes symmetric. However, for non-associated theory, the elasto-plastic constitutive matrix {\mathbf{D}}_{ep} becomes un-symmetric.

Weak formulation

The weak form of governing equations can be written as [ 26, 28],

(15) {\displaystyle \underset{\mathrm{\Omega }}{\int }\mathbf{\sigma }(\mathbf{u}):\mathbf{\epsilon }(\mathbf{v})\mathrm{d}\mathrm{\Omega }={\displaystyle \underset{\mathrm{\Omega }}{\int }\mathbf{b}.\mathbf{v}\mathrm{d}\mathrm{\Omega }+{\displaystyle \underset{{\mathit{\Gamma }}_t}{\int }\overline{\mathbf{t}}.\mathbf{v}\mathrm{d}\mathit{\Gamma }}}}.

On substituting the trial and test functions and using the arbitrariness of nodal variations, the following discrete system of equations are obtained

(16) [\mathbf{K}]\{\mathbf{d}\}=\{\mathbf{f}\},

where \mathbf{K} is the global stiffness matrix, \mathbf{d} is the vector of nodal unknowns and \mathbf{f} is the external force vector.

Displacement approximation for crack, inclusions and holes

At a particular node of interest {\mathbf{x}}_i , the displacement approximation for 2-D body having multiple discontinuities such as cracks, inclusions and holes can be written as Refs. [ 19, 28, 39],

(17) {\mathbf{u}}^h\left(\mathbf{x}\right)={\displaystyle \sum _{i=1}^n{N}_i\left(\mathbf{x}\right)\left[{\overline{u}}_i+\left\{\begin{array}{l}\underset{i\in n_r}{\underbrace{\left[\mathit{{\rm H}}\left(\mathbf{x}\right)-\mathit{{\rm H}}\left({\mathbf{x}}_i\right)\right]{\mathbf{a}}_i}}+\underset{i\in n_A}{\underbrace{{\displaystyle \sum _{\alpha =1}^4\left[{\beta }_{\alpha }\left(\mathbf{x}\right)-{\beta }_{\alpha }\left({\mathbf{x}}_i\right)\right]{\mathbf{b}}_i^{\alpha }}}}+\mathrm{.......}\hfill \\ \underset{i\in n_i}{\underbrace{\left[\phi \left(\mathbf{x}\right)-\phi \left({\mathbf{x}}_i\right)\right]{\mathbf{c}}_i}}+\underset{i\in n_h}{\underbrace{\left[\psi \left(\mathbf{x}\right)-\psi \left({\mathbf{x}}_i\right)\right]{\mathbf{d}}_i}}\hfill \end{array}\right\}R\left(\mathbf{x}\right)\right]},

where {\overline{\mathbf{u}}}_i = nodal displacement vector associated with the continuous part of the FE solution; n = set of all nodes in the mesh; n_r = set of nodes belonging to those elements which are completely cut by the crack; n_A = set of nodes belonging to those elements which are partially cut by the crack; H(\mathbf{x}) = Heaviside function, defined for elements completely intersected by the crack; {\mathbf{a}}_i = nodal enriched degrees of freedom associated with H(\mathbf{x}) ; {\beta }_{\alpha }(\mathbf{x}) = the asymptotic crack tip enrichment functions, which are defined as,

(18) {\beta }_{\alpha }(\mathbf{x})=\{{\beta }_{\mathrm{1}},{\beta }_{\mathrm{2}},{\beta }_{\mathrm{3}},{\beta }_{\mathrm{4}}\}=\left[\begin{array}{llll}{r}^{\ell }\cos {\displaystyle \frac{\theta }{\mathrm{2}}},& r^{\ell }\sin {\displaystyle \frac{\theta }{\mathrm{2}}},& r^{\ell }\sin {\displaystyle \frac{\theta }{\mathrm{2}}}\cos \theta ,& r^{\ell }\cos {\displaystyle \frac{\theta }{\mathrm{2}}}\cos \theta \end{array}\right],

where (r,\theta ) is a polar coordinate system with its origin defined at the crack tip. For linear analysis, the value of exponent \ell =\mathrm{0.5} , and for nonlinear analysis, the value of exponent \ell =\left({\displaystyle \raisebox{1ex}{$\mathrm{1}$}\!\left/ \!\raisebox{-1ex}{$\mathrm{1}+n$}\right.}\right) , where n is the strain hardening exponent which depends on the material. {\mathbf{b}}_i^{\alpha } = nodal enriched degrees of freedom vector associated with crack tip enrichment {\beta }_{\alpha }(\mathbf{x}) ; n_i = set of nodes belonging to those elements which are cut by inclusions; n_h = set of nodes belonging to those elements which are cut by holes; {\mathbf{c}}_i = nodal enriched degrees of freedom associated with \phi (\mathbf{x}) ; where \phi (\mathbf{x}) is level set function which can be defined as,

(19) \phi (\mathbf{x})=\pm \min \Vert \mathbf{x}\mathbf{-}{\mathbf{x}}_{\mathrm{\Gamma }}\Vert ,

where {\mathbf{x}}_{\mathrm{\Gamma }} is the nearest point on the interface from the point \mathbf{x} . {\mathbf{d}}_i = nodal enriched degrees of freedom associated with \psi (\mathbf{x}) ; where \psi (\mathbf{x}) is the Heaviside function, which takes the value+ 1 for the nodes lying outside the hole and 0 for the nodes lying inside the hole.

R(\mathbf{x}) = ramp function [ 63], used to model the blending elements. This can be defined as, R(\mathbf{x})={\displaystyle \sum _{i\in n_{\mathrm{e}\mathrm{n}\mathrm{r}}}{N}_i(\mathbf{x})} , where n_{\mathrm{e}\mathrm{n}\mathrm{r}} is the set of enriched nodes.

The elemental matrices \mathbf{k} and \mathbf{f} in Eq. (16), are obtained using the approximation function defined in Eq. (17),

(20) {\mathbf{k}}_{ij}^e=\left[\begin{array}{lllll}{\mathbf{k}}_{ij}^{uu}& {\mathbf{k}}_{ij}^{ua}& {\mathbf{k}}_{ij}^{ub}& {\mathbf{k}}_{ij}^{uc}& {\mathbf{k}}_{ij}^{ud}\\ {\mathbf{k}}_{ij}^{au}& {\mathbf{k}}_{ij}^{aa}& {\mathbf{k}}_{ij}^{ab}& {\mathbf{k}}_{ij}^{ac}& {\mathbf{k}}_{ij}^{ad}\\ {\mathbf{k}}_{ij}^{bu}& {\mathbf{k}}_{ij}^{ba}& {\mathbf{k}}_{ij}^{bb}& {\mathbf{k}}_{ij}^{bc}& {\mathbf{k}}_{ij}^{bd}\\ {\mathbf{k}}_{ij}^{cu}& {\mathbf{k}}_{ij}^{ca}& {\mathbf{k}}_{ij}^{cb}& {\mathbf{k}}_{ij}^{cc}& {\mathbf{k}}_{ij}^{cd}\\ {\mathbf{k}}_{ij}^{du}& {\mathbf{k}}_{ij}^{da}& {\mathbf{k}}_{ij}^{db}& {\mathbf{k}}_{ij}^{dc}& {\mathbf{k}}_{ij}^{dd}\end{array}\right],

(21) {\mathbf{f}}^h={\left\{\begin{array}{llllllll}{\mathbf{f}}_i^u& {\mathbf{f}}_i^a& {\mathbf{f}}_i^{b\mathrm{1}}& {\mathbf{f}}_i^{b\mathrm{2}}& {\mathbf{f}}_i^{b\mathrm{3}}& {\mathbf{f}}_i^{b\mathrm{4}}& {\mathbf{f}}_i^c& {\mathbf{f}}_i^d\end{array}\right\}}^{\mathrm{T}}.

Also, the sub-matrices and vectors are given as,

{\mathbf{k}}_{ij}^{rs}={\underset{{\mathrm{\Omega }}^e}{\int }({\mathbf{{\rm B}}}_i^r)}^{\mathrm{T}}\mathbf{C}{\mathbf{{\rm B}}}_j^{s}h\mathrm{d}\mathrm{\Omega },

where,

(22) r,s=u,a,b,c,d,

(23) {\mathbf{f}}_i^u={\displaystyle \underset{{\mathrm{\Omega }}^e}{\int }{\mathrm{{\rm N}}}_i\mathbf{b} \mathrm{d}\mathrm{\Omega }}+{\displaystyle \underset{{\mathrm{\Gamma }}_t}{\int }{\mathrm{{\rm N}}}_i\overline{\mathbf{t}} d\mathrm{\Gamma }},

(24) {\mathbf{f}}_i^a={\displaystyle \underset{{\mathrm{\Omega }}^e}{\int }{\mathrm{{\rm N}}}_i(\mathit{{\rm H}}(\mathbf{x})-\mathit{{\rm H}}({\mathbf{x}}_i))\mathbf{b} \mathrm{d}\mathrm{\Omega }}+{\displaystyle \underset{{\mathrm{\Gamma }}_t}{\int }{\mathrm{{\rm N}}}_i(\mathit{{\rm H}}(\mathbf{x})-\mathit{{\rm H}}({\mathbf{x}}_i))\overline{\mathbf{t}} \mathrm{d}\mathrm{\Gamma }},

(25) {\mathbf{f}}_i^{b\alpha }={\displaystyle \underset{{\mathrm{\Omega }}^e}{\int }{\mathrm{{\rm N}}}_i({\beta }_{\alpha }(\mathbf{x})-{\beta }_{\alpha }({\mathbf{x}}_i))\mathbf{b} \mathrm{d}\mathrm{\Omega }}+{\displaystyle \underset{{\mathrm{\Gamma }}_t}{\int }{\mathrm{{\rm N}}}_i({\beta }_{\alpha }(\mathbf{x})-{\beta }_{\alpha }({\mathbf{x}}_i))\overline{\mathbf{t}} \mathrm{d}\mathrm{\Gamma }},

where \alpha =\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{4.}

(26) {\mathbf{f}}_i^c={\displaystyle \underset{{\mathrm{\Omega }}^e}{\int }{\mathrm{{\rm N}}}_i\phi (\mathbf{x})\mathbf{b} \mathrm{d}\mathrm{\Omega }}+{\displaystyle \underset{{\mathrm{\Gamma }}_t}{\int }{\mathrm{{\rm N}}}_i\phi (\mathbf{x})\overline{\mathbf{t}} \mathrm{d}\mathrm{\Gamma }},

(27) {\mathbf{f}}_i^d={\displaystyle \underset{{\mathrm{\Omega }}^e}{\int }{\mathrm{{\rm N}}}_i(\psi (\mathbf{x})-\psi ({\mathbf{x}}_i))\mathbf{b} \mathrm{d}\mathrm{\Omega }}+{\displaystyle \underset{{\mathrm{\Gamma }}_t}{\int }{\mathrm{{\rm N}}}_i(\psi (\mathbf{x})-\psi ({\mathbf{x}}_i))\overline{\mathbf{t}} \mathrm{d}\mathrm{\Gamma }},

where {\mathit{{\rm N}}}_i are finite element shape function, {\mathbf{{\rm B}}}_i^u , {\mathbf{{\rm B}}}_i^a , {\mathbf{{\rm B}}}_i^b , {\mathbf{{\rm B}}}_i^c and {\mathbf{{\rm B}}}_i^d are the matrices of shape function derivatives given by

(28) {\mathbf{{\rm B}}}_i^u=[\begin{array}{l}{\mathit{{\rm N}}}_{i,x}\\ \mathrm{0}\\ {\mathit{{\rm N}}}_{i,y}\end{array}\begin{array}{l}\mathrm{0}\\ {\mathit{{\rm N}}}_{i,y}\\ {\mathit{{\rm N}}}_{i,x}\end{array}],

(29) {\mathbf{{\rm B}}}_i^a=\left[\begin{array}{ll}{({\mathit{{\rm N}}}_i(H(\mathbf{x})-H({\mathbf{x}}_i)))}_{,x}& \mathrm{0}\\ \mathrm{0}& {({\mathit{{\rm N}}}_i(H(\mathbf{x})-H({\mathbf{x}}_i)))}_{,y}\\ {({\mathit{{\rm N}}}_i(H(\mathbf{x})-H({\mathbf{x}}_i)))}_{,y}& {({\mathit{{\rm N}}}_i(H(\mathbf{x})-H({\mathbf{x}}_i)))}_{,x}\end{array}\right],

(30) {\mathbf{{\rm B}}}_i^b=\left[\begin{array}{llll}{\mathbf{{\rm B}}}_i^{b\mathrm{1}}& {\mathbf{{\rm B}}}_i^{b\mathrm{2}}& {\mathbf{{\rm B}}}_i^{b\mathrm{3}}& {\mathbf{{\rm B}}}_i^{b\mathrm{4}}\end{array}\right],

(31) {\mathbf{{\rm B}}}_i^{b\alpha }=\left[\begin{array}{ll}{({\mathit{{\rm N}}}_i({\beta }_{\alpha }(\mathbf{x}))-({\beta }_{\alpha }({\mathbf{x}}_i)))}_{,x}& \mathrm{0}\\ \mathrm{0}& {({\mathit{{\rm N}}}_i({\beta }_{\alpha }(\mathbf{x}))-({\beta }_{\alpha }({\mathbf{x}}_i)))}_{,y}\\ {({\mathit{{\rm N}}}_i({\beta }_{\alpha }(\mathbf{x}))-({\beta }_{\alpha }({\mathbf{x}}_i)))}_{,y}& {({\mathit{{\rm N}}}_i({\beta }_{\alpha }(\mathbf{x}))-({\beta }_{\alpha }({\mathbf{x}}_i)))}_{,x}\end{array}\right],

where \alpha =\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{4} .

(32) {\mathbf{{\rm B}}}_i^c=\left[\begin{array}{ll}{({\mathit{{\rm N}}}_i(\phi (\mathbf{x})-\phi ({\mathbf{x}}_i)))}_{,x}& \mathrm{0}\\ \mathrm{0}& {({\mathit{{\rm N}}}_i(\phi (\mathbf{x})-\phi ({\mathbf{x}}_i)))}_{,y}\\ {({\mathit{{\rm N}}}_i(\phi (\mathbf{x})-\phi ({\mathbf{x}}_i)))}_{,y}& {({\mathit{{\rm N}}}_i(\phi (\mathbf{x})-\phi ({\mathbf{x}}_i)))}_{,x}\end{array}\right],

(33) {\mathbf{{\rm B}}}_i^d=\left[\begin{array}{ll}{({\mathit{{\rm N}}}_i(\psi (\mathbf{x})-\psi ({\mathbf{x}}_i)))}_{,x}& \mathrm{0}\\ \mathrm{0}& {({\mathit{{\rm N}}}_i(\psi (\mathbf{x})-\psi ({\mathbf{x}}_i)))}_{,y}\\ {({\mathit{{\rm N}}}_i(\psi (\mathbf{x})-\psi ({\mathbf{x}}_i)))}_{,y}& {({\mathit{{\rm N}}}_i(\psi (\mathbf{x})-\psi ({\mathbf{x}}_i)))}_{,x}\end{array}\right].

Modeling of crack face contact

A change in boundary conditions characterizes a crack under cyclic loading. To prevent the penetration of crack faces into each other at the minimum load under cyclic loading, some mechanism must be implemented for the crack growth simulations. Therefore, to model this phenomenon, several methods have been proposed such as, changing the stiffness of the spring elements attached to the crack faces, imposing crack face nodal constraints, use of truss elements on the crack face and use of contact elements. Newman [ 5] was the first who implemented the spring elements in FEM to simulate the change in boundary conditions. In this approach, the element is connected to each boundary node on the crack face. For the open nodes, the spring stiffness is set equal to zero, while for the closed nodes, the spring stiffness is assigned a large value. McClung et al. [ 64, 65] followed Newman's approach in their earlier studies and found that the imposition of large stiffness at constrained crack face nodes leads to numerical difficulties. Thereafter, to overcome such difficulty, a direct approach has been used by Solanki et al. [ 10, 66] and Blom and Holm [ 8]. In this approach, the crack face nodal displacements are monitored during each unloading steps, and if the nodal displacement between any two unloading steps becomes negative, which means the node is closed and the nodal fixity must be applied to prevent the crack face penetration during further unloading. Now, during loading steps, the reaction forces of the fixed nodes are monitored, and when the reaction force at these nodes becomes positive the nodal fixity is removed.

In the present work, this approach is used to model the contact behavior between the crack faces. In XFEM, the conformal mesh is not required hence the nodes do not lie on the crack faces. Therefore, to measure the displacements and reaction forces at the crack faces, fictitious nodes (points), as shown in Fig. 2, are generated. In Fig. 2, circles represent the standard nodes while the black dotes represent the fictitious nodes. These fictitious nodes are not the part of mesh, hence no need to initialize them. Fictitious nodes are used to compute the displacements and reaction forces through interpolation. During the unloading, crack closure phenomenon develops near the crack tip. The effect of crack closure is small, hence a small region behind the crack tip will observe this phenomenon. Thus, only six fictitious nodes are generated behind the crack tip. As the crack propagates, new fictitious nodes (Fig. 2(b)) are generated for interpolating the new displacements and reaction forces.

Computation of J-integral

In the present work, Ramberg-Osgood material model [ 67] is used to simulate the elasto-plastic behavior of the material.

(34) {\displaystyle \frac{\epsilon }{{\epsilon }_{\mathrm{0}}}}={\displaystyle \frac{\sigma }{{\sigma }_{\mathrm{0}}}}+\alpha {\left({\displaystyle \frac{\sigma }{{\sigma }_{\mathrm{0}}}}\right)}^n,

where \alpha is the strain hardening coefficient and n is the strain hardening exponent. In the presence of plasticity, the singularity in stress field near the crack tip is modified to Hutchinson-Rice-Rosengren (HRR) singularity [ 68, 69]. The displacement field in the presence of HRR singularity can be written as,

(35) \mathbf{u}=\alpha {\epsilon }_{\mathrm{0}}r{\left({\displaystyle \frac{J}{\alpha {\epsilon }_{\mathrm{0}}{\sigma }_{\mathrm{0}}{I}_{n}r}}\right)}^{\frac{n}{n+\mathrm{1}}}\overline{\mathbf{u}}(\theta ,n),

where J is the J-integral; I_n is a dimensionless parameter, depends on n; and \overline{\mathbf{u}} is a dimensionless function of \theta and n.

Li et al. [ 70] expressed the J-integral in the equivalent area integral form as,

(36) J={\displaystyle {\int }_A\left[{\sigma }_{ij}{\displaystyle \frac{\partial u_i}{\partial x_{\mathrm{1}}}}-W{\delta }_{\mathrm{1}j}\right]}{\displaystyle \frac{\partial q_{\mathrm{1}}}{\partial x_j}}\mathrm{d}A,

where W={\displaystyle {\int }_{\mathrm{0}}^{{\epsilon }_{ij}}{\sigma }_{ij}\mathrm{d}{\epsilon }_{ij}} is the strain energy density; {\sigma }_{ij} is the stress tensor; u_i is displacement field vector; {\delta }_{ij} is the Kronecker’s delta; A is the area which contains both top and bottom crack faces; q is a function which is unity along the inner boundary of A and zero along the outer boundary of A.

Ishikawa et al. [ 71] and Bui [ 72] proposed an effective technique to separate the J-integral into symmetric and antisymmetric modes. The modes are separated by decomposing the near crack tip fields into their symmetric and antisymmetric parts with respect to crack. Consider a coordinate system x₁ and x₂ centered at the crack tip as shown in Fig. 3. Then the total displacement field can be expressed as,

(37) \mathbf{u}={\mathbf{u}}_{\mathrm{I}}+{\mathbf{u}}_{\mathrm{I}\mathrm{I}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left\{\begin{array}{l}{\mathbf{u}}_{\mathrm{1}}+{\overline{\mathbf{u}}}_{\mathrm{1}}\\ {\mathbf{u}}_{\mathrm{2}}-{\overline{\mathbf{u}}}_{\mathrm{2}}\end{array}\right\}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left\{\begin{array}{l}{\mathbf{u}}_{\mathrm{1}}-{\overline{\mathbf{u}}}_{\mathrm{1}}\\ {\mathbf{u}}_{\mathrm{2}}+{\overline{\mathbf{u}}}_{\mathrm{2}}\end{array}\right\},

where {\mathbf{u}}_{\mathrm{I}} and {\mathbf{u}}_{\mathrm{I}\mathrm{I}} are the symmetric and antisymmetric parts of displacement field about crack, and \overline{\mathbf{u}}(x_{\mathrm{1}},x_{\mathrm{2}})=\mathbf{u}(x_{\mathrm{1}},-x_{\mathrm{2}}) .

A similar decomposition can be used for stress field,

(38) \mathbf{\sigma }={\mathbf{\sigma }}_{\mathrm{I}}+{\mathbf{\sigma }}_{\mathrm{I}\mathrm{I}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left\{\begin{array}{ll}{\mathbf{\sigma }}_{\mathrm{11}}+{\overline{\mathbf{\sigma }}}_{\mathrm{11}}& {\mathbf{\sigma }}_{\mathrm{12}}-{\overline{\mathbf{\sigma }}}_{\mathrm{12}}\\ sym& {\mathbf{\sigma }}_{\mathrm{22}}+{\overline{\mathbf{\sigma }}}_{\mathrm{22}}\end{array}\right\}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left\{\begin{array}{ll}{\mathbf{\sigma }}_{\mathrm{11}}-{\overline{\mathbf{\sigma }}}_{\mathrm{11}}& {\mathbf{\sigma }}_{\mathrm{12}}+{\overline{\mathbf{\sigma }}}_{\mathrm{12}}\\ sym& {\mathbf{\sigma }}_{\mathrm{22}}-{\overline{\mathbf{\sigma }}}_{\mathrm{22}}\end{array}\right\}.

The mode separated J-integral values can be evaluated by Eq. (36) using the decomposed displacement and stress fields. In this work, small strain plasticity is assumed; hence the size of the plastic zone size developed ahead of the crack tip is small. The domain for the integration is taken in such a way that the plastic zone size remains inside it. Therefore, the following relations hold and can be used to evaluate stress intensity factors in mode-I and mode-II,

(39) J_{\mathrm{I}}=G_{\mathrm{I}}={\displaystyle \frac{K_{\mathrm{I}}^{\mathrm{2}}}{E}} \mathrm{o}\mathrm{r} K_{\mathrm{I}}=\sqrt{G_{\mathrm{I}}E}=\sqrt{J_{\mathrm{I}}E},

(40) J_{\mathrm{I}\mathrm{I}}=G_{\mathrm{I}\mathrm{I}}={\displaystyle \frac{K_{\mathrm{I}\mathrm{I}}^{\mathrm{2}}}{E}} \mathrm{o}\mathrm{r} K_{I\mathrm{I}}=\sqrt{G_{\mathrm{I}\mathrm{I}}E}=\sqrt{J_{\mathrm{I}\mathrm{I}}E}.

Calculation of fatigue life

In real situation, crack propagates along a curved path but in this analysis, it is assumed that the crack growth take place through successive linear increments. Thus, there is need to determine the length and direction of these increments. In this work, the crack growth increment size is kept fixed for each crack growth step, while the direction of crack growth is computed through maximum principal stress criterion. As per this criterion, the crack growth occurs in a direction perpendicular to the maximum principal stress [ 39, 73]. Thus, at each crack tip, the local direction of crack growth {\theta }_{\mathrm{c}} is determined by the condition that the local shear stress should be zero, i.e.,

(41) K_{\mathrm{I}}\sin {\theta }_{\mathrm{c}}+K_{\mathrm{I}\mathrm{I}}(\mathrm{3}\cos {\theta }_{\mathrm{c}}-\mathrm{1})=\mathrm{0.}

The solution of Eq. (41) gives

(42) {\theta }_{\mathrm{c}}=\mathrm{2}{\tan }^{-\mathrm{1}}\left({\displaystyle \frac{K_{\mathrm{I}}\pm \sqrt{K_{\mathrm{I}}^{\mathrm{2}}+\mathrm{8}{K}_{\mathrm{I}\mathrm{I}}^{\mathrm{2}}}}{\mathrm{4}{K}_{\mathrm{I}\mathrm{I}}}}\right),

where {\theta }_{\mathrm{c}} is the crack growth angle for each crack increment and measured in local crack tip coordinate system. If K_{\mathrm{I}\mathrm{I}} = 0 then {\theta }_{\mathrm{c}} = 0 (pure mode-I). By noting that if K_{\mathrm{I}\mathrm{I}} >0, then crack growth angle {\theta }_{\mathrm{c}} <0, and vice versa. According to this criterion, the equivalent mode-I SIF is given as,

(43) K_{\mathrm{I}\mathrm{e}\mathrm{q}}=K_{\mathrm{I}}{\cos }^{\mathrm{3}}\left({\displaystyle \frac{{\theta }_{\mathrm{c}}}{\mathrm{2}}}\right)-\mathrm{3}{K}_{\mathrm{I}\mathrm{I}}{\cos }^{\mathrm{2}}\left({\displaystyle \frac{{\theta }_{\mathrm{c}}}{\mathrm{2}}}\right)\sin \left({\displaystyle \frac{{\theta }_{\mathrm{c}}}{\mathrm{2}}}\right).

In cyclic loading, the residual compressive stresses occur along the crack faces at minimum load, which indicates the crack closure phenomenon near the crack tip as shown in Fig. 1. These residual stresses are used to calculate a negative residual stress intensity factor ( K_{\mathrm{r}\mathrm{e}\mathrm{s}} ). However, this residual stress intensity factor does not have physical sense, but by using a superposition argument, it can be used to compute the opening stress intensity factor ( K_{\mathrm{o}\mathrm{p}\mathrm{e}\mathrm{n}} ). Hence, the opening SIF is given as [ 74, 75],

(44) K_{\mathrm{o}\mathrm{p}\mathrm{e}\mathrm{n}}={(K_{\mathrm{I}\mathrm{e}\mathrm{q}})}_{\min }+K_{\mathrm{r}\mathrm{e}\mathrm{s}}.

Due to this phenomenon, the crack driving stress intensity factor known as the effective stress intensity factor ( \mathrm{\Delta }{K}_{\mathrm{e}\mathrm{f}\mathrm{f}} ) gets reduced as illustrated in Fig. 4. This effective SIF can be obtained from the difference of equivalent mode-I SIF at maximum load and the opening SIF at minimum load.

For linear analysis, the effective SIF can be obtained as,

(45) \mathrm{\Delta }{K}_{\mathrm{e}\mathrm{f}\mathrm{f}}={(K_{\mathrm{I}\mathrm{e}\mathrm{q}})}_{\max }-{(K_{\mathrm{I}\mathrm{e}\mathrm{q}})}_{\min },

whereas, for elasto-plastic analysis, the effective SIF can be obtained as,

(46) \mathrm{\Delta }{K}_{\mathrm{e}\mathrm{f}\mathrm{f}}={(K_{\mathrm{I}\mathrm{e}\mathrm{q}})}_{\max }-K_{\mathrm{o}\mathrm{p}\mathrm{e}\mathrm{n}},

where {(K_{\mathrm{I}\mathrm{e}\mathrm{q}})}_{\max } and {(K_{\mathrm{I}\mathrm{e}\mathrm{q}})}_{\min } are the equivalent stress intensity factors at maximum and minimum load.

The fatigue life obtained by generalized Paris law [ 19, 76] is expressed as,

(47) {\displaystyle \frac{\mathrm{d}a}{\mathrm{d}N}}=C{(\mathrm{\Delta }{K}_{\mathrm{e}\mathrm{f}\mathrm{f}})}^m,

where a is the crack length, N is the number of loading cycles, C and m are material properties. The failure takes place whenever {(K_{\mathrm{I}\mathrm{e}\mathrm{q}})}_{\max }>K_{\mathrm{I}\mathrm{C}} where K_{\mathrm{I}\mathrm{C}} is the fracture toughness of the material. A detailed procedure for the fatigue crack growth analysis is explained in Fig. 5.

Numerical results and discussion

In this paper, the various edge crack problems are simulated by XFEM in the presence of various flaws (holes, inclusions and minor-cracks). The volume fraction of these flaws i.e., holes or inclusions or holes and inclusions, is taken as 10%. Additionally, 60 minor cracks are also randomly distributed in the plate. The holes and inclusions are assumed circular in shape, and their radii vary from 2.0 to 4.0 mm. The minor cracks sizes vary from 2.0 to 4.0 mm with their orientation vary from 0° to 90°. In all the examples under consideration, the results are computed 10 times to obtain the average effect of the randomly distributed flaws. These average values are used for plotting the results. The elasto-plastic simulations for aluminum alloy 6061-T6 are performed using the generalized Ramberg-Osgood material model given in Eq. (34). The values of \alpha =\mathrm{1.169} and n=\mathrm{13.949} , obtained from the stress-strain diagram of the aluminum alloy 6061-T6, are used as input for all the examples. The material properties (ASM material data sheet, see Table 1) used for these simulations are,

In all the example problems, a plate of size 100 mm × 200 mm along with an edge crack of initial length a = 15 mm is taken for the simulations. The plate is subjected to a cyclic tensile load of {\sigma }_{\min } = 0 MPa and {\sigma }_{\max } = 50 MPa at the top edge of the plate. The bottom edge of the plate is constrained in the y-direction. A crack growth of 2 mm is given at each step until the final failure i.e., the simulation stops when {(K_{\mathrm{I}\mathrm{e}\mathrm{q}})}_{\max }>K_{\mathrm{I}\mathrm{C}} . The SIF values are calculated at each stage of crack growth to evaluate the fatigue life of the plate using Paris law. The effect of discontinuities/flaws on the fatigue life of the plate is evaluated by XFEM. A comparison of the results obtained through linear elastic and elasto-plastic analysis is discussed in detail. In the linear elastic approach, the loading and unloading follows the same path whereas in elasto-plastic analysis, the plastic region developed near the crack tip during loading does not regain its original shape during unloading, which results in the crack closure effect during unloading.

Plate with an edge crack

A plate of size 100 mm × 200 mm along with an edge crack of initial length a = 15 mm is taken for simulation as shown in Fig. 6. A convergence study for mesh size is performed in Table 2. From the table, it has been observed that the results are converged for a mesh size of 60 by 120 nodes. Therefore, for the further numerical simulations, a uniform mesh of 60 × 120 nodes discretized by four node isoparametric quadrilateral elements is used. In addition to this, the stress intensity factor also depends on the size of J-integral domain. Thus, a circular domain of radius 2 mm around the crack tip is taken for the SIF computation such that the sufficient number of elements lie inside the integral domain. In elasto-plastic analysis, the cyclic loading produces the crack closure phenomenon at the crack tip during unloading. The fatigue crack closure phenomenon is visualized here in terms of residual stresses at minimum load as shown in Fig. 7. As expected, the compressive normal stresses are developed at the crack tip and normal stresses are almost zero behind the crack tip. To satisfy the equilibrium condition and compensate for the plasticity, a large region of tensile stresses is generated around the compressive stress zone near the crack tip. For the next loading cycle, some applied tensile load will be consumed to overcome the residual compressive stresses at the crack tip, which leads to the reduction in crack driving force.

During constant amplitude cyclic loading, the maximum and opening values of SIF increase with the increase in number of loading cycles. Therefore, stabilized values of maximum and opening SIF are required for the fatigue crack growth analysis. To obtain stabilized values of SIF, eight load cycles are applied for a constant crack length. The variation of maximum and opening SIF values with number of cycles for initial crack length and final crack length is presented in Table 3.

The variation of stress intensity factor against crack length at maximum load is shown in Fig. 8(a) for both linear elastic and elasto-plastic analysis. Fig. 8(b) shows the variation of opening stress intensity factor with crack length at minimum load. In the elasto-plastic analysis, \mathrm{\Delta }{K}_{\mathrm{e}\mathrm{f}\mathrm{f}} is obtained by subtracting the opening SIF from the maximum SIF. The fatigue life obtained by XFEM is plotted with crack length as shown in Fig. 9. In case of linear elastic analysis, the fatigue life is found as 33836 cycles whereas in case of elasto-plastic analysis, the fatigue life is found as 34547 cycles. Thus, in case of elasto-plastic analysis, the fatigue life is increased nearly by 2.1% with respect to the linear elastic analysis. This increase in fatigue life is attributed to the presence of residual stresses during unloading.

Plate with an edge crack and holes

In this case, a plate of size 100 mm × 200 mm along with an edge crack of initial length a = 15 mm and multiple holes is taken for the simulation as shown in Fig. 10. Holes of arbitrary size, radii varying from 2.0 to 4.0 mm, are randomly distributed in the plate keeping the volume fraction constant at 10%. The loading, boundary conditions, mesh size and crack increments are kept same as explained above. The variation of stress intensity factor with crack length is shown in Fig. 11(a) and (b) at maximum and minimum load respectively. The variation of fatigue life with crack length is plotted in Fig. 12. Thus, in case of linear elastic analysis, the fatigue life is found as 18225 cycles whereas in case of elasto-plastic analysis, it is found as 19016 cycles. Thus, in case of elasto-plastic analysis, the fatigue life is increased by 4.3% with respect to linear elastic analysis, which is due to the presence of residual stresses during unloading.

Plate with an edge crack and inclusions

In this case, quasi-static fatigue crack propagation of an edge cracked plate with randomly distributed inclusions is performed by XFEM. A plate of size 100 mm × 200 mm with an edge crack of initial length a = 15 mm containing multiple inclusions (10% volume fraction) of radii vary from 2.0 mm to 4.0 mm is taken for simulation as shown in Fig. 13. The loading, boundary conditions, mesh size and crack increments are kept same as above. The variation of stress intensity factor with crack length at maximum load is shown in Fig. 14(a). Fig. 14(b) shows the variation of opening stress intensity factor with crack length at minimum load. The fatigue life is obtained using Paris law, and is plotted with crack extension in Fig. 15. The fatigue life is found to be 26348 and 27236 cycles for linear elastic and elasto-plastic analysis respectively. The fatigue life for elasto-plastic analysis is increased by 3.4% with respect to linear analysis.

Plate with an edge crack and minor cracks

In this case, the fatigue life of a rectangular plate having multiple minor cracks along with a major edge crack is evaluated by XFEM. A plate of size 100 mm × 200 mm having an edge crack of initial length a = 15 mm and 60 minor cracks of size varying from 2.0 to 4.0 mm is taken for the simulation as shown in Fig. 16. The inclination of the randomly distributed minor cracks varies from 0° to 90°. The modeling parameters i.e., loading, boundary conditions, mesh size and crack increments are kept same as defined above. The variation of SIF with crack length at maximum and minimum load is plotted in Figs. 17(a) and (b) respectively. A plot of fatigue life with crack length is shown in Fig. 18. In case of linear analysis, the fatigue life is found as 32054 cycles whereas in case of elasto-plastic analysis, the fatigue life is found as 32721 cycles. Thus, in case of elasto-plastic analysis, the fatigue life is increased by 2.1% with respect to linear elastic analysis.

Plate with an edge crack, holes and inclusions

In this case, fatigue crack growth is modeled in the presence of holes and inclusions. A plate of size 100 mm × 200 mm with an edge crack of initial length a = 15 mm containing randomly distributed holes and inclusions is taken for the simulation as shown in Fig. 19. The volume fraction of the flaws (holes and inclusions) is fixed to 10%. The radius of holes and inclusions vary from 2.0 mm to 4.0 mm. The loading, boundary conditions, mesh size and crack increments are kept same as explained earlier. The variation of stress intensity factor against crack length at maximum load is shown in Fig. 20(a) while, Fig. 20(b) shows the variation of opening stress intensity factor with crack length at minimum load. A plot of fatigue life against crack length is shown in Fig. 21. The fatigue life is found to be 23443 cycles and 24274 cycles for linear elastic and elasto-plastic analysis respectively. In elasto-plastic analysis, the presence of residual stresses at minimum load increase the fatigue life by 3.5%.

Plate with an edge crack, holes and minor cracks

In this case, an edge cracked plate along with randomly distributed holes and minor cracks is modeled by XFEM. The plate of size 100 mm × 200 mm having an edge crack of initial length a = 15 mm along with multiple holes (10% volume fraction) and 60 minor cracks is taken for simulation as shown in Fig. 22. The loading, boundary conditions, mesh size and crack increments are kept same as above. The variation of stress intensity factor against crack length is plotted in Fig. 23(a) and (b) at maximum and minimum load respectively. The variation of fatigue life with crack length is illustrated in Fig. 24. In case of linear elastic analysis, the fatigue life is found as 17388 cycles whereas in case of elasto-plastic analysis, the fatigue life is found as 18476 cycles. Thus, in case of elasto-plastic analysis, the fatigue life is increased by 6.3% with respect to linear elastic analysis, which is due to the presence of residual stresses at minimum load.

Plate with an edge crack, inclusions and minor cracks

In this case, a plate of size 100 mm × 200 mm with an edge crack of initial length a = 15 mm is taken for the simulation as shown in Fig. 25. The fatigue life of the plate is evaluated in the presence of inclusions (by 10% volume) and 60 minor cracks. The loading, boundary conditions, mesh size, flaws size and crack increments are kept same as above. The variation of SIF with crack length at maximum and minimum load is plotted in Figs. 26(a) and (b) respectively. A plot of the fatigue life with crack length is shown in Fig. 27. In case of linear analysis, the fatigue life is found as 25128 cycles whereas, in case of elasto-plastic analysis, the fatigue life is found as 26078 cycles. Thus, in case of elasto-plastic analysis, the fatigue life is increased by 3.8% with respect to linear elastic analysis. This increase in fatigue life is due to the crack closure phenomenon.

Plate with an edge crack, holes, minor cracks and inclusions

Finally, a plate of size 100 mm×200 mm with an edge crack of initial length a = 15 mm containing randomly distributed multiple defects is taken for simulation as shown in Fig. 28. The fatigue life of the plate is evaluated in the presence of holes and inclusions (by 10% volume) and 60 minor cracks. The modeling conditions i.e., loading, boundary conditions, mesh size, flaws size and crack increment are kept similar as earlier explained. Fig. 29(a) and (b) show the variation of SIF against crack length at maximum and minimum load respectively. A plot of the fatigue life with crack length is shown in Fig. 30. In case of linear analysis, the fatigue life is found as 21268 cycles whereas, in case of elasto-plastic analysis, the fatigue life is found as 22102 cycles. Thus, in case of elasto-plastic analysis, the fatigue life is increased by 3.8% with respect to linear elastic analysis. This increase in fatigue life is due to the compressive stresses developed during unloading.

The percentage reduction in fatigue life, due to the presence of flaws, is summarized in Table 4 for all cases discussed above. From these simulations, it is observed that the holes have got most severe effect on the fatigue life of the plate. The presence of minor cracks has got least effect on the fatigue life. The fatigue life of the plate obtained using elasto-plastic analysis is found more as compared that obtained using linear elastic analysis. This increase is due the presence of residual stresses ahead of the crack tip during unloading. Moreover, it is also noticed that the effect of plasticity is not significant in case of fatigue crack growth simulations.

Conclusions

In this work, elastic and elasto-plastic simulations have been performed by XFEM to evaluate the fatigue life of the plate in the presence of defects. A Ramberg-Osgood material model is combined with von-Mises failure criterion to simulate various problems under plane stress condition. The values of equivalent stress intensity factor are obtained from the J-integral values of two different modes. The plasticity induced crack closure phenomenon is modeled by fixing the nodal displacement of the crack faces and the effect of PICC on the fatigue life is obtained during constant amplitude cyclic loading. The results obtained using linear elastic analyses are compared with elasto-plastic analyses for various edge crack configurations. On the basis of the present simulations, the following conclusions have been drawn:

In cyclic loading, the plastic strains developed near the crack tip during loading lead to the formation of plastic wake behind the crack tip as the crack extends, which results in the development of compressive stresses at the crack tip.

The residual compressive stresses reduce the crack driving force for the next loading cycle as some portion of applied tensile load is consumed to overcome the compressive stresses developed at the crack tip. These compressive stresses results in the improvement of fatigue life by reducing \mathrm{\Delta }{K}_{eff} . Thus, plasticity developed at the crack tip improves the fatigue life.

The fatigue life obtained by linear elastic analysis is found quite close to small strain plastic analysis. Hence, it can be stated that the linear elastic analysis is sufficient to perform fatigue crack growth simulations for small scale yielding.

The presence of flaws significantly affects the fatigue life.

The effect of holes is found most severe in comparison to other defects.

References

[1] Wolf E. Fatigue crack closure under cyclic tension. Engineering Fracture Mechanics, 1970, 2(1): 37–45

[2] Guo W, Wang C H, Rose L R F. Influence of cross-sectional thickness on fatigue crack growth. Fatigue & Fracture of Engineering Materials & Structure, 1999, 22(5): 437–444

[3] de Matos P F P, Nowell D. Experimental and numerical investigation of thickness effects in plasticity-induced fatigue crack closure. International Journal of Fatigue, 2009, 31(11–12): 1795–1804

[4] Antunes F V, Branco R, Costa J D, Rodrigues M. Plasticity induced crack closure in middle-crack tension specimen: Numerical versus experimental. Fatigue & Fracture of Engineering Materials & Structure, 2010, 33(10): 673–686

[5] Newman J C Jr. A crack-closure model for predicting fatigue crack growth under aircraft spectrum loading. ASTM STP 748, Philadelphia, PA, 1981, 53–84

[6] Budiansky B, Hutchinson J W. Analysis of closure in fatigue crack growth. ASME Journal of Applied Mechanics, 1992, 42: 389–400

[7] Newman J C Jr. Finite element analysis of fatigue crack closure. ASTM STP 590, Philadelphia, PA, 1976, 281–301

[8] Blom A F, Holm D K. An experimental and numerical study of crack closure. Engineering Fracture Mechanics, 1985, 22(6): 997–1011

[9] McClung R C, Thacker B H, Roy S. Finite element visualization of fatigue crack closure in plane stress and plane strain. International Journal of Fracture, 1991, 50: 27–49

[10] Solanki K, Daniewicz S R, Newman J C Jr. Finite element modeling of plasticity-induced crack closure with emphasis on geometry and mesh refinement effects. Engineering Fracture Mechanics, 2003, 70(12): 1475–1489

[11] Alizadeh H, Hills D A, de Matos P F P, Nowell D, Pavier M J, Paynter R J, Smith D J, Simandjuntak S. A comparison of two and three-dimensional analysis of fatigue crack closure. International Journal of Fatigue, 2007, 29(2): 222–231

[12] Ellyin F, Ozah F. The effect of material model in describing mechanism of plasticity-induced crack closure under variable cyclic loading. International Journal of Fracture, 2007, 143(1): 15–33

[13] Toribio J, Kharin V. Finite-deformation analysis of the crack-tip fields under cyclic loading. International Journal of Solids and Structures, 2009, 46(9): 1937–1952

[14] Cheung S, Luxmoore A R. A finite element analysis of stable crack growth in an aluminium alloy. Engineering Fracture Mechanics, 2003, 70(9): 1153–1169

[15] Yan A M, Nguyen-Dang H. Multiple-cracked fatigue crack growth by BEM. Computational Mechanics, 1995, 16(5): 273–280

[16] Yan X. A boundary element modeling of fatigue crack growth in a plane elastic plate. Mechanics Research Communications, 2006, 33(4): 470–481

[17] Belytschko T, Gu L, Lu Y Y. Fracture and crack growth by element-free Galerkin methods. Modelling and Simulation in Materials Science and Engineering, 1994, 2(3A): 519–534

[18] Belytschko T, Lu Y Y, Gu L. Crack propagation by element-free Galerkin methods. Engineering Fracture Mechanics, 1995, 51(2): 295–315

[19] Duflot M, Nguyen-Dang H. Fatigue crack growth analysis by an enriched meshless method. Journal of Computational and Applied Mathematics, 2004, 168(1–2): 155–164

[20] Rabczuk T, Belytschko T. A three-dimensional large deformation meshfree method for arbitrary evolving cracks. Computer Methods in Applied Mechanics and Engineering, 2007, 196(29–30): 2777–2799

[21] Rabczuk T, Areias P M A, Belytschko T. A meshfree thin shell method for non-linear dynamic fracture. International Journal for Numerical Methods in Engineering, 2007, 72(5): 524–548

[22] Rabczuk T, Bordas S, Zi G. A three-dimensional meshfree method for continuous multiple-crack initiation, propagation and junction in statics and dynamics. Computational Mechanics, 2007, 40(3): 473–495

[23] Rabczuk T, Samaniego E. Discontinuous modelling of shear bands using adaptive meshfree methods. Computer Methods in Applied Mechanics and Engineering, 2008, 197(6–8): 641–658

[24] Bordas S, Rabczuk T, Zi G. Three-dimensional crack initiation, propagation, branching and junction in non-linear materials by an extended meshfree method without asymptotic enrichment. Engineering Fracture Mechanics, 2008, 75(5): 943–960

[25] Nguyen V P, Rabczuk T, Bordas S, Duflot M. Meshless methods: A review and computer implementation aspects. Mathematics and Computers in Simulation, 2008, 79(3): 763–813

[26] Belytschko T, Black T. Elastic crack growth in finite elements with minimal remeshing. International Journal for Numerical Methods in Engineering, 1999, 45(5): 601–620

[27] Stolarska M, Chopp D, Moes N, Belytschko T. Modeling crack growth by level sets in the extended finite element method. International Journal for Numerical Methods in Engineering, 2001, 51(8): 943–960

[28] Sukumar N, Chopp D L, Moes N, Belytschko T. Modeling of holes and inclusions by level sets in the extended finite element method. Computer Methods in Applied Mechanics and Engineering, 2001, 190(46–47): 6183–6200

[29] Rabczuk T, Belytschko T. Cracking particles: A simplified meshfree method for arbitrary evolving cracks. International Journal for Numerical Methods in Engineering, 2004, 61(13): 2316–2343

[30] Rabczuk T, Belytschko T. Application of particle methods to static fracture of reinforced concrete structures. International Journal of Fracture, 2006, 137(1–4): 19–49

[31] Rabczuk T, Zi G, Bordas S, Nguyen-Xuan H. A simple and robust three-dimensional cracking-particle method without enrichment. Computer Methods in Applied Mechanics and Engineering, 2010, 199(37–40): 2437–2455

[32] Rabczuk T, Zi G, Gerstenberger A, Wall W A. A new crack tip element for the phantom-node method with arbitrary cohesive cracks. International Journal for Numerical Methods in Engineering, 2008, 75(5): 577–599

[33] Bordas S P A, Rabczuk T, Hung N X, Nguyen V P, Natarajan S, Bog T, Quan D M, Hiep N V. Strain smoothing in FEM and XFEM. Computers & Structures, 2010, 88(23–24): 1419–1443

[34] Bordas S P A, Natarajan S, Kerfriden P, Augarde C E, Mahapatra D R, Rabczuk T, Pont S D. On the performance of strain smoothing for quadratic and enriched finite element approximations (XFEM/GFEM/PUFEM). International Journal for Numerical Methods in Engineering, 2011, 86(4–5): 637–666

[35] Budarapu P R, Gracie R, Yang S W, Zhuang X, Rabczuk T. Efficient coarse graining in multiscale modeling of fracture. Theoretical and Applied Fracture Mechanics, 2014, 69: 126–143

[36] Amiri F, Anitescu C, Arroyo M, Bordas S P A, Rabczuk T. XLME interpolants, a seamless bridge between XFEM and enriched meshless methods. Computational Mechanics, 2014, 53(1): 45–57

[37] Bhardwaj G, Singh I V, Mishra B K. Stochastic fatigue crack growth simulation of interfacial crack in bi-layered FGMs using XIGA. Computer Methods in Applied Mechanics and Engineering, 2015, 284: 186–229

[38] Kumar S, Singh I V, Mishra B K, Rabczuk T. Modeling and simulation of kinked cracks by virtual node XFEM. Computer Methods in Applied Mechanics and Engineering, 2015, 283: 1425–1466

[39] Singh I V, Mishra B K, Bhattacharya S, Patil R U. The numerical simulation of fatigue crack growth using extended finite element method. International Journal of Fatigue, 2012, 36(1): 109–119

[40] Dolbow J, Moes N, Belytschko T. An extended finite element method for modelling crack growth with frictional contact. Computer Methods in Applied Mechanics and Engineering, 2001, 190(51–52): 6825–6846

[41] Moes N, Gravouil A, Belytschko T. Non-planar 3D crack growth with the extended finite element and level set−Part I: Mechanical model. International Journal for Numerical Methods in Engineering, 2002, 53: 2549–2568

[42] Pathak H, Singh A, Singh I V, Yadav S K. A simple and efficient XFEM approach for 3-D cracks simulations. International Journal of Fracture, 2013, 181(2): 189–208

[43] Rethore J, Gravouil A, Combescure A. An energy-conserving scheme for dynamic crack growth using the extended finite element method. International Journal for Numerical Methods in Engineering, 2005, 63(5): 631–659

[44] Kumar S, Singh I V, Mishra B K, Singh A. New enrichments in XFEM to model dynamic crack response of 2-D elastic solids. International Journal of Impact Engineering, 2015

[45] Natarajan S, Baiz P M, Bordas S, Rabczuk T, Kerfriden P. Natural frequencies of cracked functionally graded material plates by the extended finite element method. Composite Structures, 2011, 93(11): 3082–3092

[46] Baiz P M, Natarajan S, Bordas S P A, Kerfriden P, Rabczuk T. Linear buckling analysis of cracked plates by SFEM and XFEM. Journal of Mechanics of Materials and Structures, 2011, 6(9–10): 1213–1238

[47] Vu-Bac N, Nguyen-Xuan H, Chen L, Bordas S, Kerfriden P, Simpson R N, Liu G R, Rabczuk T. A Node-based smoothed extended finite element method (NS-XFEM) for fracture analysis, CMES. Computer Modeling in Engineering & Sciences, 2011, 73: 331–356

[48] Chen L, Rabczuk T, Bordas S P A, Liu G R, Zeng K Y, Kerfriden P. Extended finite element method with edge-based strain smoothing (Esm-XFEM) for linear elastic crack growth. Computer Methods in Applied Mechanics and Engineering, 2012, 209–212: 250–265

[49] Legrain G, Moes N, Verron E. Stress analysis around crack tips in finite strain problems using the extended finite element method. International Journal for Numerical Methods in Engineering, 200<?Pub Caret?>5, 63: 209–314

[50] Ji H, Chopp D, Dolbow J E. A hybrid finite element/level set method for modelling phase transformation. International Journal for Numerical Methods in Engineering, 2002, 54(8): 1209–1233

[51] Elguedj T, Gravouil A, Combescure A. A mixed augmented Lagrangian-extended finite element method for modeling elastic-plastic fatigue crack growth with unilateral contact. International Journal for Numerical Methods in Engineering, 2007, 71(13): 1569–1597

[52] Kumar S, Singh I V, Mishra B K. XFEM simulation of stable crack growth using J-R curve under finite strain plasticity. International Journal of Mechanics and Materials in Design, 2014, 10(2): 165–177

[53] Budarapu P R, Gracie R, Bordas S P A, Rabczuk T. An adaptive multiscale method for quasi-static crack growth. Computational Mechanics, 2014, 53(6): 1129–1148

[54] Yang S W, Budarapu P R, Mahapatra D R, Bordas S P A, Zi G, Rabczuk T. A meshless adaptive multiscale method for fracture. Computational Materials Science, 2015, 96: 382–395

[55] Nanthakumar S S, Lahmer T, Rabczuk T. Detection of multiple flaws in piezoelectric structures using XFEM and level sets. Computer Methods in Applied Mechanics and Engineering, 2014, 275: 98–112

[56] Kumar S, Singh I V, Mishra B K. A multigrid coupled (FE-EFG) approach to simulate fatigue crack growth in heterogeneous materials. Theoretical and Applied Fracture Mechanics, 2014, 72: 121–135

[57] Kumar S, Singh I V, Mishra B K. A homogenized XFEM approach to simulate fatigue crack growth problems. Computers & Structures, 2015, 150: 1–22

[58] Vu-Bac N, Rafiee R, Zhuang X, Lahmer T, Rabczuk T. Uncertainty quantification for multiscale modeling of polymer nanocomposites with correlated parameters. Composites. Part B, Engineering, 2015, 68: 446–464

[59] Vu-Bac N, Silani M, Lahmer T, Zhuang X, Rabczuk T. A unified framework for stochastic predictions of mechanical properties of polymeric nanocomposites. Computational Materials Science, 2015, 96: 520–535

[60] Hinton E, Owen D. Finite Element in Plasticity. Pineridge Press Limited, 1980, 215–265

[61] Moes N, Dolbow J, Belytschko T. A finite element method for crack growth without remeshing. International Journal for Numerical Methods in Engineering, 1999, 46: 131–150

[62] Bland D R. The associated flow rule of plasticity. Journal of the Mechanics and Physics of Solids, 1957, 6(1): 71–78

[63] Fries T P. A corrected XFEM approximation without problems in blending elements. International Journal for Numerical Methods in Engineering, 2008, 75(5): 503–532

[64] McClung R C, Sehitoglu H. On the finite element analysis of fatigue crack closure-1. Basic modeling issues. Engineering Fracture Mechanics, 1983, 33(2): 237–252

[65] McClung R C, Sehitoglu H. On the finite element analysis of fatigue crack closure-2. Numerical results. Engineering Fracture Mechanics, 1983, 33(2): 253–257

[66] Solanki K, Daniewicz S R, Newman J C Jr. A new methodology for computing crack opening values from finite element analyses. Engineering Fracture Mechanics, 2004, 71(7–8): 1165–1175

[67] Elguedj T, Gravouil A, Combescure A. Appropriate extended functions for XFEM simulation of plastic fracture mechanics. Computer Methods in Applied Mechanics and Engineering, 2006, 195(7–8): 501–515

[68] Hutchinson J W. Singular behavior at the end of a tensile crack in a hardening material. Journal of the Mechanics and Physics of Solids, 1968, 16(1): 13–31

[69] Rice J R, Rosengren G F. Plane strain deformation near a crack tip in a power-law hardening material. Journal of the Mechanics and Physics of Solids, 1968, 16(1): 1–2

[70] Li F Z, Shih C F, Needleman A. A comparison of methods for calculating energy release rates. Engineering Fracture Mechanics, 1985, 21(2): 405–421

[71] Ishikawa H. A finite element analysis of stress intensity factors for combined tensile and shear loading by only a virtual crack extension. International Journal of Fracture, 1980, 16(5): 243– 246

[72] Bui H D. Associated path independent J-integrals for separating mixed modes. Journal of the Mechanics and Physics of Solids, 1983, 6(6): 439–448

[73] Erdogan F, Sih G. On the crack extension in plates under plane loading and transverse shear. Journal of Basic Engineering, 1963, 85(4): 519–527

[74] de Matos P F P, Nowell D. On the accurate assessment of crack opening and closing stresses in plasticity-induced fatigue crack closure problems. Engineering Fracture Mechanics, 2007, 74(10): 1579–1601

[75] Antunes F V, Chegini A G, Correia L, Branco R. Numerical study of contact forces for crack closure analysis. International Journal of Solids and Structures, 2014, 51(6): 1330–1339

[76] Paris P C, Gomez M P, Anderson W E. A rational analytic theory of fatigue. Trend in Engineering, 1961, 13: 9–14

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag Berlin Heidelberg