Fatigue crack growth simulations of 3-D linear elastic cracks under thermal load by XFEM

Abstract

This paper deals with the fatigue crack growth simulations of three-dimensional linear elastic cracks by XFEM under cyclic thermal load. Both temperature and displacement approximations are extrinsically enriched by Heaviside and crack front enrichment functions. Crack growth is modelled by successive linear extensions, and the end points of these linear extensions are joined by cubic spline segments to obtain a modified crack front. Different crack geometries such as planer, non-planer and arbitrary spline shape cracks are simulated under thermal shock, adiabatic and isothermal loads to reveal the sturdiness and versatility of the XFEM approach.

Keywords

3-D cracks / fatigue life / Paris law / thermal load / XFEM

Introduction

Engineering components are often subjected to thermo-elastic loads. Few well-known examples of such components are automobiles, aircraft, and nuclear/thermal power plants. Hence, the failure of the engineering components is not only due to mechanical loads but also due to thermo-mechanical loads. Cracks are unavoidable in all the structural components. Hence, the studies of thermo-elastic cracks [ 1− 5] under thermo-elastic loads become quite important. The design against fatigue in 3-D domain under thermo-elastic load becomes even more important. These problems are difficult to solve because of challenges involved in the modeling and precise evaluation of fracture parameters such as stress intensity factors (SIFs). The present work deals with the fatigue crack growth in 3-D subjected to thermo-elastic loads. The analysis of crack growth by analytical methods is a very tedious job. Hence, numerical method would be better option to simulate fatigue cracks as it brings flexibility, and requires less simulation time as compared to the duration of experiment. So far, the finite element method (FEM) [ 6− 10] has been largely used to solve stationary and propagating crack problems. Boundary element method (BEM) [ 11− 13] is also used by researchers to solve fracture problems. FEM is found most successful and powerful numerical method for solving variety of engineering and science problems. However, the accurate modelling of crack and crack growth remains a challenging task for FEM. In FEM, the geometry is generally modelled by an adequate mesh, and a crack must coincide with the edges of the finite elements i.e. a conformal mesh is required apart from the requirement of special elements to handle crack tip asymptotic stress fields. FEM experiences difficulties in remeshing and adaptive analysis. Hence, the modelling and simulation of discontinuities by FEM becomes a cumbersome task, and often leads to inaccuracy in the numerical simulation of crack growth problems. To deal with these issues, several methods have been developed based on Partition of Unity (PU) concept to handle discontinuous domain problems namely: element free Galerkin method [ 14, 15], extended finite element method (XFEM) [ 16, 17], particle methods [ 18, 19], cracking particle method [ 20− 22], extended bridging domain [ 23], XIGA [ 24, 25], meshless adaptive multiscale method [ 26, 27], virtual node XFEM [ 28]. Among these methods, XFEM has been grown as an effective and efficient tool to solve crack growth problems [ 17, 29, 30]. XFEM removes the burden associated with mesh generation for the problems involving voids, cracks and interfaces, and thus provides the precise modelling of cracks. XFEM does not require a conformal mesh for crack growth modelling. In XFEM, the standard displacement based approximation is enriched by additional functions using the partition of unity [ 16, 31]. The enrichment functions are obtained from the theoretical background of the problem under consideration. Primarily two enrichment functions [ 32, 33] are required to model a crack in XFEM: first one is discontinuous on the crack surface while the second one is asymptotic at the crack front. Lots of efforts have been made in past by the computational community to solve 3-D crack problems using XFEM. The level set method [ 32, 34] is combined with XFEM to model the crack growth problems. Möes et al. [ 35, 36] presented 3-D XFEM formulation using vector level set. Sukumar et al. [ 30] coupled the XFEM with fast marching method for fatigue crack growth simulations. Shi et al. [ 37] presented 3-D XFEM for cracks in ABAQUS using fast marching and level set methods. Till date, only a limited number of papers have been published on the fatigue crack growth by XFEM such as cohesive zone fatigue crack growth modeling [ 38− 40], fatigue crack growth simulation [ 41, 42], three-dimensional crack growth [ 29] and fatigue life estimation of 2-D homogenous [ 43], FGM [ 44− 47] and heterogeneous plate [ 48]. Pathak et al. [ 49] extended XFEM for three-dimensional cracks simulation in linear elastic materials. Recently, it has been implemented to simulate crack problems in piezoelectric domain under mechanical and thermal loading environment [ 50− 56].

In above mentioned studies, so far, the fatigue crack growth is mainly limited to two-dimensional simulations under mechanical load. Therefore, the main aim of this work is to evaluate the fatigue life of 3-D cracked structures under thermal load. Both regular and arbitrary shape cracks have been taken for the simulation. A non-iterative scheme is adopted to define the geometry of the crack surface. The crack growth is modelled by successive line segments and the crack front is modelled by segments of cubic splines [ 57]. The vector level set approach [ 35, 36] is used for analyzing the non-planer 3-D cracks. Mixed mode stress intensity factors are computed by domain based interaction integral approach [ 49]. The maximum principal stress criterion has been used to find the crack growth direction, and Paris law of fatigue crack growth has been used to evaluate fatigue life under thermal load.

Problem formulations

Consider a three dimensional cracked domain ( \mathrm{\Omega } ) consisting of {\mathrm{\Gamma }}_c , {\mathrm{\Gamma }}_T , {\mathrm{\Gamma }}_t and {\mathrm{\Gamma }}_u (Fig. 1). All boundary conditions like displacement, temperature and traction are imposed on the boundaries {\mathrm{\Gamma }}_u , {\mathrm{\Gamma }}_T and {\mathrm{\Gamma }}_t . Crack surfaces are assumed to be traction free. The equilibrium and boundary conditions for linear thermo-elastic problem may be described as

(1(a)) \begin{array}{l}-\nabla \mathbf{q}+\mathrm{Q}=\mathrm{0}\hfill \\ \mathbf{q}=-k\nabla T\hfill \\ \nabla :\mathbf{\sigma }+\mathbf{b}=\mathrm{0} \mathrm{in} \mathrm{\Omega }\hfill \\ \mathbf{\sigma }=\mathbf{C}:(\mathbf{\epsilon }-{\epsilon }_T)\hfill \end{array}

(1(b)) \begin{array}{l}\mathbf{\sigma }\cdot \mathbf{n}=\overline{\mathbf{t}} \mathrm{on} {\mathrm{\Gamma }}_t\hfill \\ \mathbf{\sigma }\cdot \mathbf{n}=\mathrm{0} \mathrm{on} {\mathrm{\Gamma }}_{c^{+}}\hfill \\ \mathbf{\sigma }\cdot \mathbf{n}=\mathrm{0} \mathrm{on} {\mathrm{\Gamma }}_{c^{-}}\hfill \end{array},

and thermo-elastic strain tensor is defined as \epsilon ={\nabla }_s\mathbf{u} and {\epsilon }_{\mathbf{T}}=\alpha (T-T_{ref})\mathbf{I} .

In above governing equation, q is heat flux vector, Q represents the heat source, k is the thermal conductivity of material, n is the unit outward normal on the boundaries, u is displacement field, σ is Cauchy stress tensor, b is the body forces per unit volume, C is isotropic fourth order tensor, {\nabla }_s is the symmetric gradient operator, T is the temperature field within the domain, ε_T is the strain vector, ε_T is the thermal strain vector with respect to reference temperature, T_{ref} , \alpha is the thermal expansion coefficient and I represents the second order identity tensor.

Using constitutive relation, the weak form of elastic part of governing Eq. (1) can be written as [ 58],

(2) {\displaystyle \underset{\mathrm{\Omega }}{\int }\mathbf{\epsilon }(\mathbf{u}):\mathbf{C}:\mathbf{\epsilon }(\mathbf{V})d}\mathrm{\Omega }-{\displaystyle \underset{\mathrm{\Omega }}{\int }\mathbf{b}.\mathbf{v}}d\mathrm{\Omega }-{\displaystyle \underset{{\mathrm{\Gamma }}_t}{\int }\overline{\mathbf{t}}\cdot \mathbf{v}d\mathrm{\Gamma }}=\mathrm{0}\hspace{1em}\forall \mathbf{v}\in u_o,

where

(3) u_o=\{\mathbf{v}\in \Re \mathbf{:}\mathbf{v}=\mathrm{0}\mathrm{on}{\mathrm{\Gamma }}_{\mathrm{c}}, \mathbf{v}\mathrm{discontinuous}\mathrm{on}{\mathrm{\Gamma }}_{\mathrm{c}}\}.

The discrete equations are obtained by displacement approximation, weight and shape functions in Eq. (2).

In thermo-elastic adiabatic crack problems, both temperature and displacement fields are discontinuous across the crack surface whereas heat flux is singular at the crack tip [ 59]. Therefore, the displacement and temperature fields can be defined as,

(4a) u_o=\left\{\mathbf{u}\in \Re \mathbf{:}\mathbf{u}=\overline{\mathbf{u}}\mathrm{on}{\mathrm{\Gamma }}_{\mathrm{u}}, \mathbf{u}\mathrm{discontinuous}\mathrm{on}{\mathrm{\Gamma }}_{\mathrm{c}}\right\},

(4b) T_o=\left\{T\in \Re \mathbf{:}T=\overline{T}\mathrm{on}{\mathrm{\Gamma }}_{\mathrm{T}}, T\mathrm{discontinuous}\mathrm{on}{\mathrm{\Gamma }}_{\mathrm{c}}\right\}.

In case of thermo-elastic isothermal crack problem, a specific temperature is prescribed on the crack surface therefore essential boundary conditions are imposed at the crack surface. In this case, discontinuity in heat flux is found across the crack surface whereas continuous temperature profile is found at the crack surface [ 60]. Angular variation of temperature and heat flux field obtained in isothermal crack is different then the adiabatic crack. In case of an adiabatic crack, temperature varies perpendicular to the crack surface, whereas in case of an isothermal crack, it varies in radial direction.

The temperature field may be defined as,

(5) T_o=\left\{T\in \Re \mathbf{:}T=\overline{T} \mathrm{on} {\mathrm{\Gamma }}_{\mathrm{T}}/{\mathrm{\Gamma }}_{\mathrm{c}}, \mathbf{q}\mathrm{discontinuous}\mathrm{on}{\mathrm{\Gamma }}_{\mathrm{c}}\right\}.

The weak form of heat governing equation (1) using constitutive equation becomes [ 60],

(6) {\displaystyle \underset{\mathrm{\Omega }}{\int }\mathbf{q}(S)k \mathbf{q}(T)d}\mathrm{\Omega }+{\displaystyle \underset{\mathrm{\Omega }}{\int }SQ}d\mathrm{\Omega }-{\displaystyle \underset{{\mathrm{\Gamma }}_q}{\int }S \overline{q}d\mathrm{\Gamma }}=\mathrm{0},

(7) {\displaystyle \underset{\mathrm{\Omega }}{\int }\mathbf{\epsilon }(\mathbf{u})\mathbf{:}\mathbf{C}\mathbf{:}\mathbf{\epsilon }(\mathbf{v})\mathrm{d}}\mathrm{\Omega }-{\displaystyle \underset{\mathrm{\Omega }}{\int }\mathbf{b}\mathbf{.}\mathbf{v} \mathrm{d}\mathrm{\Omega }}-{\displaystyle \underset{{\mathrm{\Gamma }}_t}{\int }\overline{\mathbf{t}}\mathbf{.}\mathbf{v} \mathrm{d}\mathrm{\Gamma }}-{\displaystyle \underset{\mathrm{\Omega }}{\int }{\mathbf{\epsilon }}_T(\mathbf{u})\mathbf{:}\mathbf{C}\mathbf{:}\mathbf{\epsilon }(\mathbf{v})\mathrm{d}}\mathrm{\Omega }=\mathrm{0.}

Using the displacement and temperature approximations, trial and test functions and use of arbitrariness of nodal variation, a set of discrete equations can be written as,

(8(a)) [\mathbf{K}]\{\mathbf{u}\}=\{\mathbf{f}\},

(8(b)) [\mathbf{K}]\{T\}=\{\mathbf{f}\},

where, u and T are the nodal unknowns, K and f are the global stiffness matrix and external force vector. The stiffness matrix and force vector are computed at element level and then assembled into their global counterparts through usual finite element assembly procedure.

Thermo-elastic crack problems are decoupled into thermal and structural problems. First, discrete equation for heat conduction Equation (8(b)) has been solved to get temperature distribution in the cracked geometry, further these temperature data are used as body force to get solution for elastic discrete Equation (8(a)). Elastic equation gives the solution in the terms of deformation u. In post-processing phase, thermal interaction integral approach has been used to extract individual stress intensity factor under thermal mixed-mode loading environment. In this way, decoupled algorithm has been used to get solution for thermos-elastic crack problems. As steady state heat conduction equation is solved to get temperature distribution in the domain, and then this temperature is used as load input for elastic equation to get displacement in the solution.

The additional degrees of freedom arise due to the enrichment and are handled by considering the fictitious nodes. The elemental contributions for K are given as,

(9a) {\mathbf{K}}_{ij}^e={\displaystyle \underset{{\mathrm{\Omega }}^e}{\int }{({\mathbf{{\rm B}}}_i^{\alpha })}_{case}^T}\overline{\mathbf{C}}{({\mathbf{{\rm B}}}_j^{\alpha })}_{case}d\mathrm{\Omega }\{(\alpha =classical,\mathrm{enriched})\hspace{1em}(case=elastic,\mathrm{thermal})\},

(9b) {\mathbf{f}}_i^e={\left\{\begin{array}{ll}{\mathbf{f}}_i^{classical}& {\mathbf{f}}_i^{enrich}\end{array}\right\}}^{\mathrm{T}},

(10) \overline{\mathbf{C}}=\{\begin{array}{ll}Cond\hfill & \mathrm{f}\mathrm{o}\mathrm{r} \mathrm{h}\mathrm{e}\mathrm{a}\mathrm{t} \mathrm{e}\mathrm{q}\mathrm{u}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n},\hfill \\ \mathbf{C}\hfill & \mathrm{f}\mathrm{o}\mathrm{r} elastic \mathrm{e}\mathrm{q}\mathrm{u}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{.}\hfill \end{array}

The constitutive matrix \mathbf{C} for a three dimensional linear isotropic elastic material can be written as,

(11a) \mathbf{C}={\displaystyle \frac{E}{(\mathrm{1}+\nu )(\mathrm{1}-\mathrm{2}\nu )}}\left[\begin{array}{llllll}\mathrm{1}-\nu & \nu & \nu & \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \nu & \mathrm{1}-\nu & \nu & \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \nu & \nu & \mathrm{1}-\nu & \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0.5}-\nu & \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0.5}-\nu & \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0.5}-\nu \end{array}\right].

E is the Young’s modulus, and \nu is the Poisson’s ratio.

The thermal conductivity matrix Cond for three-dimensional linear isotropic materials is defined as,

(11b) Cond=\left[\begin{array}{lll}k& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& k& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& k\end{array}\right].

The sub-matrices and vectors that appear in the foregoing equations are given as,

(12a) {({\mathbf{f}}_i^{classical})}_{thermal}={\displaystyle \underset{{\mathrm{\Gamma }}_t}{\int }{\mathrm{{\rm N}}}_i\mathbf{q}\mathrm{d}\mathrm{\Gamma }+{\displaystyle \underset{{\mathrm{\Omega }}^e}{\int }{\mathrm{{\rm N}}}_{i}Q\mathrm{d}\mathrm{\Omega }}},

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

(12c) {({\mathbf{f}}_i^{enriched})}_{thermal}={\displaystyle \underset{{\mathrm{\Gamma }}_t}{\int }{\mathrm{{\rm N}}}_i(\psi (\mathbf{x})-\psi ({\mathbf{x}}_i))\mathbf{q}\mathrm{d}\mathrm{\Gamma }+{\displaystyle \underset{{\mathrm{\Omega }}^e}{\int }{\mathrm{{\rm N}}}_i(\psi (\mathbf{x})-\psi ({\mathbf{x}}_i))Q\mathrm{d}\mathrm{\Omega }}},

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

where {\mathrm{{\rm N}}}_i is finite element shape function and \psi is additional enrichment function (Heaviside function or crack tip branch functions).

(12e) {({\mathrm{{\rm B}}}_i^{classical})}_{thermal}=\left[\begin{array}{lll}{N}_{i,x}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& N_{i,y}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& N_{i,z}\end{array}\right],

(12f) {({\mathrm{{\rm B}}}_i^{classical})}_{elastic}=\left[\begin{array}{lll}{N}_{i,x}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& N_{i,y}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& N_{i,z}\\ \mathrm{0}& N_{i,z}& N_{i,y}\\ N_{i,z}& \mathrm{0}& N_{i,x}\\ N_{i,y}& N_{i,x}& \mathrm{0}\end{array}\right].

To recover Kronecker delta property, shifted enrichment [ 61] is used with partition of unity.

(12g) {({\mathrm{{\rm B}}}_i^{enriched})}_{thermal}=\left[\begin{array}{lll}{(N_i(\psi (\mathrm{x})-\psi ({\mathbf{x}}_i)))}_{,x}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& {(N_i(\psi (\mathbf{x})-\psi ({\mathbf{x}}_i)))}_{,y}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& {(N_i(\psi (\mathbf{x})-\psi ({\mathbf{x}}_i)))}_{,z}\end{array}\right],

(12h) {({\mathrm{{\rm B}}}_i^{enriched})}_{elastic}=\left[\begin{array}{lll}{(N_i(\psi (\mathbf{x})-\psi ({\mathbf{x}}_i)))}_{,x}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& {(N_i(\psi (\mathbf{x})-\psi ({\mathbf{x}}_i)))}_{,y}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& {(N_i(\psi (\mathbf{x})-\psi ({\mathbf{x}}_i)))}_{,z}\\ \mathrm{0}& {(N_i(\psi (\mathbf{x})-\psi ({\mathbf{x}}_i)))}_{,z}& {(N_i(\psi (\mathbf{x})-\psi ({\mathbf{x}}_i)))}_{,y}\\ {(N_i(\psi (\mathbf{x})-\psi ({\mathbf{x}}_i)))}_{,z}& \mathrm{0}& {(N_i(\psi (\mathbf{x})-\psi ({\mathbf{x}}_i)))}_{,x}\\ {(N_i(\psi (\mathbf{x})-\psi ({\mathbf{x}}_i)))}_{,y}& {(N_i(\psi (\mathbf{x})-\psi ({\mathbf{x}}_i)))}_{,x}& \mathrm{0}\end{array}\right].

Crack growth tracking

In XFEM, geometrical discontinuities are tracked by using vector level set function [ 35, 36]. To reduce computational time and memory, vector level set is evaluated for those elements which have their mean centre lying in a narrow band of crack surface. Initially, vector level values are stored for particular crack surface, and as crack propagates, vector level set is updated by performing geometric operations [ 57]. Thus, a geometric representation of the crack is not needed, and a crack is entirely described by numerical data. Apart from crack tracking, vector level set is used to reduce the computational time [ 49],

• Stored level set functions are used for the calculation of the enrichment function in pre-processing.

• In post-processing, stored level set functions are used for the evaluation of auxiliary fields to calculate the stress intensity factors.

Since, the level sets are approximated by quadratic shape functions, the representation of the crack surface and crack front is more accurate. In this work, 8-node brick elements are used for the approximation of field variables and 20 node quadratic elements are used for the level set approximation at crack front [ 57]. The vector level set is assumed to be a compound function given by a vector and a Boolean. The signed distance function is used as a particular level set function.

(13) [\phi (\mathbf{x})\stackrel{\sim }{\psi }(\mathbf{x})]=\{f(\mathbf{x}),H(f(\mathbf{x}))\},

where, f(\mathbf{x}) is the vector between the point x and its projection on the crack surface {\mathrm{\Gamma }}_{cr} . It is oriented in such a way that it points from x to the crack surface, i.e. f=\overline{\mathbf{x}}-\mathbf{x} where \overline{\mathbf{x}} is the closest point projection of x on {\mathrm{\Gamma }}_{cr} . The signed distance function \phi (\mathbf{x}) can be computed from f(\mathbf{x}) using \phi (\mathbf{x})=\left|\left|f(\mathbf{x})\right|\right|H(f(\mathbf{x})) , where H(\mathbf{x}) is Heaviside function which takes+1 value above the crack surface and-1 below the crack surface. Crack surface between elements is represented by [ 49, 57],

(14a) {\phi }_{approx}(\mathbf{x})={\displaystyle \sum _{i=\mathrm{1}}^n{N}_i(\mathbf{x})\phi ({\mathbf{x}}_i)},

(14b) {\stackrel{\sim }{\psi }}_{approx}(\mathbf{x})={\displaystyle \sum _{i=\mathrm{1}}^n{N}_i(\mathbf{x})\stackrel{\sim }{\psi }({\mathbf{x}}_i)}.

The criterion for the finding whether an element belongs to crack front or crack surface, is given as,

(15a) \mathrm{F}\mathrm{o}\mathrm{r} \mathrm{c}\mathrm{r}\mathrm{a}\mathrm{c}\mathrm{k} \mathrm{f}\mathrm{r}\mathrm{o}\mathrm{n}\mathrm{t} \mathrm{e}\mathrm{n}\mathrm{r}\mathrm{i}\mathrm{c}\mathrm{h}\mathrm{e}\mathrm{d} \mathrm{e}\mathrm{l}\mathrm{e}\mathrm{m}\mathrm{e}\mathrm{n}\mathrm{t}\mathrm{s}:{\stackrel{\sim }{\psi }}_{approx}(\mathbf{x})=\mathrm{0} \mathrm{a}\mathrm{n}\mathrm{d} {\phi }_{approx}(\mathbf{x})=\mathrm{0.}

(15b)

Enrichment functions

Heaviside function is utilized for the modelling of strong discontinuity in displacement and temperature fields. The Heaviside function [ 62, 63] is defined as,

(16) H(\mathbf{x})=\{\begin{array}{ll}\mathrm{1},\hfill & if(\mathbf{x}-\mathbf{x}\mathbf{*}).n\ge \mathrm{0};\hfill \\ -\mathrm{1},\hfill & \mathrm{otherwise}.\hfill \end{array}

The following crack tip enrichment functions are used to model the stress singularity at the crack tip.

(17) \beta (\mathbf{x})=\left[\sqrt{r}\sin {\displaystyle \frac{\theta }{\mathrm{2}}},\sqrt{r}\cos {\displaystyle \frac{\theta }{\mathrm{2}}},\sqrt{r}\cos {\displaystyle \frac{\theta }{\mathrm{2}}}\sin \theta ,\sqrt{r}\sin {\displaystyle \frac{\theta }{\mathrm{2}}}\sin \theta \right].

In order to implement an enrichment function at crack front, a nearest point from a Gauss point (evaluation point) is selected by using a local coordinate system. Initially crack front is divided into many piecewise curve spline parts, and then a nearest point on the crack segment from a Gauss point is evaluated.

Consider an arbitrary three-dimensional crack [ 49, 57] as shown in Fig. 2. The geometry of the crack front can be represented by the position vector r(\mathbf{x}) pointing from the origin of the (X,Y,Z) global Cartesian coordinate system to some point t on the crack front. The mathematical description of the crack front and crack surface can be written as,

(18a) r(\mathbf{x})=\mathrm{0},

(18b) f(\mathbf{x})=\mathrm{0.}

A curvilinear coordinate system [ 49] on the curve is denoted by {\xi }_{\mathrm{1}},{\xi }_{\mathrm{2}} and {\xi }_{\mathrm{3}} . The direction of {\xi }_{\mathrm{1}} is taken outward normal to the crack front and lie in the plane of the crack, {\xi }_{\mathrm{3}} is the perpendicular to crack surface, and is named as the gradient of the crack surface, and {\xi }_{\mathrm{2}} is tangent to the crack front and is calculated by the vector multiplication of {\xi }_{\mathrm{1}} and {\xi }_{\mathrm{3}} .

(19) {\xi }_{\mathrm{2}}={\xi }_{\mathrm{3}}\times {\xi }_{\mathrm{1}}.

Let {\xi }_{\mathrm{1}}=a_{\mathrm{1}}\mathbf{i}+b_{\mathrm{1}}\mathbf{j}+c_{\mathrm{1}}\mathbf{k} , {\xi }_{\mathrm{2}}=a_{\mathrm{2}}\mathbf{i}+b_{\mathrm{2}}\mathbf{j}+c_{\mathrm{2}}\mathbf{k} , {\xi }_{\mathrm{3}}=a_{\mathrm{3}}\mathbf{i}+b_{\mathrm{3}}\mathbf{j}+c_{\mathrm{3}}\mathbf{k} where i, j, k are unit vectors in X, Y and Z directions respectively. The perpendicular from a point/node on the surface is denoted by x_f. Now, a local vector from the foot to a point P can be written as,

(20) \mathbf{f}(\mathbf{x})={\mathbf{x}}_f-\mathbf{x}.

Local coordinate of a point with respect to the crack front can be represented by

(21) \tilde{x}=\mathbf{f}(\mathbf{x})\cdot {\mathbf{\xi }}_{\mathbf{1}},\tilde{y}=\mathbf{f}(\mathbf{x})\cdot {\mathbf{\xi }}_{\mathbf{2}},\tilde{z}=\mathbf{f}(\mathbf{x})\cdot {\mathbf{\xi }}_{\mathbf{3}}.

In orthogonal coordinate system, \tilde{z} automatically becomes zero, thus r and \theta can be calculated in terms of \tilde{x} and \tilde{y} as,

(22) r=\sqrt{\left({\tilde{x}}^{\mathrm{2}}+{\tilde{y}}^{\mathrm{2}}\right)},\theta ={\tan }^{-\mathrm{1}}\left({\displaystyle \frac{\tilde{y}}{\tilde{x}}}\right),

where, \tilde{x} and \tilde{y} are the function of (X,Y,Z) . After finding the values of \tilde{x} and \tilde{y} at each node, it can be easily approximated using partition of unity. This process is completed in small bend of elements,

(23a) {\tilde{x}}_{approx}(\xi ,\eta ,\zeta )={\displaystyle \sum _{i=\mathrm{1}}^{i=n}{\tilde{x}}_i{N}_i(\xi ,\eta ,\zeta )},

(23b) {\tilde{y}}_{approx}(\xi ,\eta ,\zeta )={\displaystyle \sum _{i=\mathrm{1}}^{i=n}{\tilde{y}}_i{N}_i(\xi ,\eta ,\zeta )}.

( \xi ,\eta ,\zeta ) are the coordinates of a point in the natural coordinate system. After creating local polar coordinate system, r and \theta can be easily differentiated using the derivatives of shape functions [ 57]. In the present work, the quadratic shape functions are used to approximate the level set functions while the linear shape functions are used to approximate the solution.

(24) r=\sqrt{\left({\tilde{x}}_{approx}^{\mathrm{2}}(\xi ,\eta ,\zeta )+{\tilde{y}}_{approx}^{\mathrm{2}}(\xi ,\eta ,\zeta )\right)},\theta ={\tan }^{-\mathrm{1}}\left({\displaystyle \frac{{\tilde{y}}_{approx}(\xi ,\eta ,\zeta )}{{\tilde{x}}_{approx}(\xi ,\eta ,\zeta )}}\right).

Numerical integration

The problem domain is numerically integrated using standard Gauss quadrature. In cracked elements, numerical integration is performed by dividing them into several tetrahedrons above and below the crack surface. To achieve this, initially a point of intersection of crack front with element surfaces and edges is calculated [ 30], then higher order Gauss quadrature is used for the integration in these enriched elements. In the present simulations, three point Gauss quadrature has been used in the elements intersected by crack surface whereas seven point Gauss quadrature has been used in the crack front elements. The rest of the elements in the problem domain are numerically integrated by two point Gauss quadrature.

Evaluation of stress intensity factors

Computation of accurate stress intensity factors (SIFs) under mixed mode thermal loading condition is quite essential to analyze the behavior of cracked problems. In this work, individual SIFs (K_I,K_{II} \mathrm{a}\mathrm{n}\mathrm{d} K_{III}) have been obtained by modified form of domain based interaction integral approach. A domain form of J-integral [ 57] is provided in Fig. 3. Virtual domain extension approach [ 64] has been implemented for three-dimensional evaluation of J-integral. In this approach, a crack front contour integral is expressed in terms of volume integral over the domain surrounding the crack front. Two states of the stress are superimposed with each other to extract the individual SIFs. For convenience, one state (auxiliary state) is assumed to be known while other one is the actual state. The SIF of the auxiliary state is taken as one for the mode which is being evaluated in actual state while the SIFs of other two modes in auxiliary states, are assigned zero values. The path independent J-integral for a homogeneous cracked body is given as

(25) J={\displaystyle \underset{\mathrm{\Gamma }}{\int }\left(W{\delta }_{\mathrm{1}j}-{\sigma }_{ij}{\displaystyle \frac{\partial u_i}{\partial x_{\mathrm{1}}}}\right)}{n}_j\mathrm{d}\mathit{\Gamma },

where, \mathrm{\Gamma } is an arbitrary contour which encloses the crack tip, W_s is the strain energy density, n_j is the outward unit normal to the contour and {\delta }_{\mathrm{1}j} is the Kronecker delta function. A local coordinate which is parallel to the crack surface has been created.

Interaction integral

Two independent equilibrium states have been taken to evaluate the interaction integral. State 1 corresponds to an actual state while state 2 is taken as auxiliary state which can be obtained from analytical solutions of asymptotic stress and displacement fields. After introducing actual and auxiliary fields, the J-integral can be defined as

(26) J^{(\mathrm{1}+\mathrm{2})}=J^{(\mathrm{1})}+J^{(\mathrm{2})}+M^{(\mathrm{1},\mathrm{2})},

where, J^{(\mathrm{1}+\mathrm{2})} is the J-integral of the superimposed state, J^{(\mathrm{1})} is the J-integral in the actual state, J^{(\mathrm{2})} is the J-integral due to the auxiliary state and M^{(\mathrm{1},\mathrm{2})} is the interaction term.

The domain form of the interaction integral under mechanical loading can be written as

(27) M^{(\mathrm{1},\mathrm{2})}={\displaystyle \underset{\mathrm{\Omega }}{\int }\left[{\sigma }_{ij}^{(\mathrm{1})}{\displaystyle \frac{\partial u_i^{(\mathrm{2})}}{\partial x_{\mathrm{1}}}}+{\sigma }_{ij}^{(\mathrm{2})}{\displaystyle \frac{\partial u_i^{(\mathrm{1})}}{\partial x_{\mathrm{1}}}}-W^{(\mathrm{1},\mathrm{2})}{\delta }_{ij}\right]{\displaystyle \frac{\partial q}{\partial x_j}}\mathrm{d}\mathrm{\Omega }},

(28) W^{(\mathrm{1},\mathrm{2})}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}({\sigma }_{ij}^{(\mathrm{1})}{\epsilon }_{ij}^{(\mathrm{2})}+{\sigma }_{ij}^{(\mathrm{2})}{\epsilon }_{ij}^{(\mathrm{1})}),

where, q is a scalar weight function whose value is one at the crack front, zero at the contour and arbitrary between crack front and contour. {\sigma }_{ij}^{(\mathrm{1})} and {\epsilon }_{ij}^{(\mathrm{1})} , are the actual Cauchy stress and engineering strain respectively, while {\sigma }_{ij}^{(\mathrm{2})} and {\epsilon }_{ij}^{(\mathrm{2})} are the auxiliary Cauchy stress and engineering strain respectively. W^{(\mathrm{1},\mathrm{2})} represents the strain energy density in actual and auxiliary states. Similarly, the thermal interaction integral [ 65] is described as

(29) M^{(\mathrm{1},\mathrm{2})}={\displaystyle \underset{\mathrm{\Omega }}{\int }\left[{\sigma }_{ij}^{(\mathrm{1})}{\displaystyle \frac{\partial u_i^{(\mathrm{2})}}{\partial x_{\mathrm{1}}}}+{\sigma }_{ij}^{(\mathrm{2})}{\displaystyle \frac{\partial u_i^{(\mathrm{1})}}{\partial x_{\mathrm{1}}}}-W^{(\mathrm{1},\mathrm{2})}{\delta }_{ij}\right]}{\displaystyle \frac{\partial q}{\partial x_j}}\mathrm{d}A+\alpha {\displaystyle \underset{A}{\int }{\displaystyle \frac{\partial T}{\partial x_{\mathrm{1}}}}}{\sigma }_{kk}^{(\mathrm{2})}q\mathrm{d}\mathrm{\Omega }.

In thermal interaction integral expression, \alpha represents thermal expansion coefficient.

SIF calculation

For linear elastic problems, the interaction integral is related to the mixed mode stress intensity factors by the following relation:

(30) M^{(\mathrm{1},\mathrm{2})}={\displaystyle \frac{\mathrm{2}(\mathrm{1}-{\nu }^{\mathrm{2}})}{E}}[K_I^{(\mathrm{1})}{K}_I^{(\mathrm{2})}+K_{II}^{(\mathrm{1})}{K}_{II}^{(\mathrm{2})}]+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mu }}{K}_{III}^{(\mathrm{1})}{K}_{III}^{(\mathrm{2})}.

In this work, a hollow cuboidal contour is used for the calculation of M^{(\mathrm{1},\mathrm{2})} integral. The calculation of auxiliary field is same as the calculation of enrichment function. In this work, auxiliary field are approximated using partition of unity approach

(31a) {\sigma }^{(\mathrm{2})}(j,k)\approx {\displaystyle \sum _{i=\mathrm{1}}^n{\sigma }_i^{(\mathrm{2})}(j,k)N_i}(X,Y,Z),

(31b) u^{(\mathrm{2})}\approx {\displaystyle \sum _{i=\mathrm{1}}^n{u}_i^{(\mathrm{2})}{N}_i}(X,Y,Z),

where, i is the node number.

Fatigue crack growth

In this work, incremental crack growth based on concept of linear elastic fracture mechanics has been used under constant amplitude cyclic thermal loading. Standard Paris law of fatigue crack growth has been used for the life prediction of crack geometries under cyclic thermal loading. After evaluating the SIFs using domain based interaction integral approach, the range of SIFs for constant amplitude cyclic loading is defined as,

(32) \mathrm{\Delta }K=K_{\max }-K_{\min },

where, K_{\max } and K_{\min } are the values of SIFs corresponding to maximum and minimum applied thermal boundary conditions respectively.

To get crack growth direction, the maximum principal stress criterion is used. Thus, at each crack increment point, crack growth direction {\theta }_c is obtained by equating the local shear stress equal zero [ 66],

(33) K_I\sin {\theta }_c+K_{II}(\mathrm{3}\cos {\theta }_c-\mathrm{1})=\mathrm{0.}

The solution of above equation gives,

(34) {\theta }_c=\mathrm{2}{\tan }^{-\mathrm{1}}\left({\displaystyle \frac{K_I-\sqrt{K_I^{\mathrm{2}}+\mathrm{8}{K}_{II}^{\mathrm{2}}}}{\mathrm{4}{K}_{II}}}\right),

where, angle {\theta }_c is calculated in the orthogonal plane.

For stable crack propagation model, the generalized Paris’ law is written as

(35) {\displaystyle \frac{da}{dN}}=C{(\mathrm{\Delta }{K}_{eq})}^m,

(36) \mathrm{\Delta }{K}_{eq}=\sqrt{\mathrm{\Delta }{K}_I^{\mathrm{2}}+\mathrm{\Delta }{K}_{II}^{\mathrm{2}}+(\mathrm{1}-{\upsilon }^{\mathrm{2}})\mathrm{\Delta }{K}_{III}},

where, C and m are material constant, expression for crack front velocity becomes [ 36].

(37) \mathbf{V}=C{(\mathrm{\Delta }{K}_{eq})}^m({\xi }_{\mathrm{1}}\times \cos {\theta }_c+{\xi }_{\mathrm{2}}\times \sin {\theta }_c),

where, {\xi }_{\mathrm{1}} and {\xi }_{\mathrm{2}} are defined in the Fig. 2. Crack growth is modelled by successive linear extensions of crack surface [ 57]. The ends of these linear extensions are joined by cubic spline segments to obtain a new crack front. Since, in three-dimensional mixed-mode crack simulation, the stress intensity factors may vary along the crack front. The fatigue crack growth is governed by Eq. (38), the increments along the crack front must be applied accordingly. The crack propagation length of the crack curve at different point on the crack front is given by this equation.

(38) {\displaystyle \frac{a_{\max }}{a_I}}={\left({\displaystyle \frac{{(\mathrm{\Delta }{K}_{eq})}_{\max }}{\mathrm{\Delta }{K}_{eq}}}\right)}^m.

a_{\max } , a_I are the user specified maximum crack increment and actual crack increment at crack front in the orthogonal plane. All points after crack propagation are joined by a cubic spline. This spline curve is further divided into several curved pieces for the calculation of the nearest points on the crack front. The tracking of moving surface in crack growth modeling is one of the major issues. Initially, numerous points are generated on the original crack surface and then these points are used for triangulation of crack surface. As crack grows, new points for the extended crack surface are calculated. A new approximated surface is created, and connectivity among the points is developed using triangulation [ 57] as shown in Fig. 4. These triangles are used for the computation of the level set functions. The triangulation is a judicious choice because it makes the computations easier. After the crack growth, crack front becomes arbitrary and the points on the crack front become random, these points are joined again by a cubic spline crack segments. This updated crack front again divided into hundreds of segments and then used to evaluate nearest point on crack front from an arbitrary point. Above mentioned steps for tracking of crack geometry has been repeated until the crack becomes unstable i.e. {(\mathrm{\Delta }{K}_{eq})}_{\max }>K_{IC} . A complete solution algorithm is given below:

(a) Meshing of complete domain by 8-noded brick element without considering any discontinuities in the geometry.

(b) Evaluation of vector Level set function for those elements which have their mean center lying in a narrow band of crack surface.

(c) Segmentation of crack front by finite numbers of spline curve.

(d) Creation of orthogonal coordinate system at each crack front segment.

(e) Delaunay Triangulation of crack surface.

(f) Find enriched nodes using vector level set function.

(g) Creation of hanging nodes within crack front elements.

(h) Create virtual nodal points at enriched nodes to handle additional DOFs.

(i) Creation of Gauss point and Jaccobian for each element.

(j) Loop over elements to calculate global conductive matrix for heat equation.

(i) If non-enriched element

• Standard conductivity matrix calculation for each element.

(ii) Else

• Find nearest Gauss point to each crack front segment by vector projection method.

• Local enrichment for crack discontinuities.

• Enriched conductivity matrix calculation for enriched element.

(iii) End loop over elements.

(k) Get global Heat matrix.

(l) Enforce temperature boundary condition.

(m) Obtain solution in terms of temperature.

(n) Loop over elements to calculate global stiffness matrix and body force vector for elastic equation.

(i) If non-enriched element

• Standard stiffness matrix calculation for each element.

• Using obtained temperature data calculate force matrix for each element.

(ii) Else

• Find nearest Gauss point to each crack front segment by vector projection method.

• Local enrichment for crack discontinuities.

• Enriched stiffness matrix calculation for enriched element.

• Using obtained temperature data calculate force matrix for each element.

(iii) End loop over elements.

(o) Get solution in terms of displacements.

(p) Create J-domain position along the crack front.

(q) Evaluate mixed mode SIFs (K_I, K_II and K_III) using thermal interaction integral approach.

(r) Check for crack growth

(i) If (K_eq<K_IC)

• Find crack growth direction at each J-domain position.

• Crack increment along the crack front at each J-domain position.

• Enclose the incremented J-domain position to form updated crack front.

• Level set update.

(ii) Elseif (K_eq>K_IC)

• Failure of component.

• Implement Paris law to obtain the fatigue life.

(s) End of the program.

Numerical results and discussion

The fatigue failure analysis of three-dimensional cracks in cuboid has been performed by XFEM under cyclic thermo-elastic load. Extrinsic partition of unity enrichment has been implemented to capture geometrical discontinuities present in the domain. Heaviside function and crack tip enrichment functions are used to represent a crack in the domain. Three separate cases involving cuboid with different geometries: penny shape crack, inclined elliptical crack and arbitrary shape are solved by XFEM. All problems have been solved under three different loading conditions namely thermal shock, adiabatic and isothermal. In thermal shock load, uniform temperature change is imposed throughout the problem domain. Top and bottom surfaces are constrained which results in the thermal stress. In adiabatic case, crack surfaces are assumed as thermally insulated, and a uniform heat flux is maintained in the direction perpendicular to the crack surface. In isothermal case, crack faces are maintained at specific temperature, and some different temperatures are imposed on the outer surface of problem domain, which creates radial thermal flux from outer surface to the crack face. Both adiabatic and isothermal crack problems are decoupled into thermal and elastic problems. First, the temperature field is obtained for the whole domain, then it used as input for the elastic problem. The displacement based discrete equations are solved for the field variables.

A benchmark problem, penny shaped crack in cylinder is solved to check the accuracy of the developed XFEM code. The complete geometry (Fig. 5) of cylinder having height 6 m and radius 1 m contains a penny shaped crack of radius 0.1 m, is modelled by XFEM. A uniform temperature 100°C is applied on the cylindrical surfaces while 0°C is imposed on the crack surfaces. The cylinder is constrained to move in vertical direction. The normalized SIFs along crack front for different J domain position are presented in the Fig. 6, which shows a comparison of results obtained by XFEM with those available in literature using BEM [ 4] and analytical [ 67]. The results presented in the Fig. 6 shows that XFEM results are in good agreement with analytical and BEM results.

Three different crack geometries (penny, inclined ellipse, arbitrary) in a cuboid are simulated by XFEM. Before the start of these simulations, a convergence study is performed on penny shape crack in cuboid (Fig. 7) under thermal shock loading to obtain the optimal mesh size. To check the convergence, problem is simulated for five different set of nodes i.e., (\mathrm{16}\times \mathrm{16}\times \mathrm{16}) , (\mathrm{20}\times \mathrm{20}\times \mathrm{20}) , (\mathrm{24}\times \mathrm{24}\times \mathrm{24}) , (\mathrm{28}\times \mathrm{28}\times \mathrm{28}) and (\mathrm{32}\times \mathrm{32}\times \mathrm{32}) . The values of mode-I SIF obtained from these simulations are presented in the Fig. 8. These results clearly show a convergence for mesh size of (\mathrm{24}\times \mathrm{24}\times \mathrm{24}) . Hence, the further simulations are performed for a mesh size of (\mathrm{24}\times \mathrm{24}\times \mathrm{24}) . A generalized MATLAB code has been developed to obtain the results. The material properties used in the present simulations are tabulated in Table 1:

The results are presented in the form of number of life cycles and crack growth contours. In all these simulations, the geometries are uniformly discretized by eight node brick element for a mesh of \mathrm{24}\times \mathrm{24}\times \mathrm{24} .

Penny shape crack in finite domain

A cuboid (\mathrm{2}\mathrm{m}\times \mathrm{2}\mathrm{m}\times \mathrm{2}\mathrm{m}) with penny shape crack as shown in Fig. 7 is solved for three different loading cases namely thermal shock, adiabatic and isothermal. The radius (r) of penny is taken as 0.3 m.

Thermal shock load

In this sub-section, a penny shape crack in cuboid has been analysed under cyclic thermal shock load of \mathrm{\Delta }T=-{\mathrm{20}}^{\mathrm{o}}\mathrm{C} . Top and bottom surfaces of cuboid are constrained in x, y and z directions. The uniform thermal shock load under the imposed boundary conditions induces mode-I loading on the crack. Applied thermal shock is cyclic in nature and imposed until the crack becomes unstable. Standard Paris law has been implemented to the predict fatigue life. The maximum value of the crack increment (a_max) is taken as 0.03 m as per Eq. (40). The fatigue life obtained by XFEM is presented in Fig. 9. The star marked in the figure indicates the final fatigue life of the cracked object. The fatigue life and critical crack length (crack surface diameter) are found as 21257×10² cycles and 0.8242 m respectively. A crack growth contour is presented in Fig. 10, which shows the initial and final crack surfaces. These contours also show that the crack grows in its plane, and maintain the circular shape until it fails. Front view and top view of crack growth are presented in Fig. 11, which show that the crack growth is planer, and remains circular after each crack increment. The magnified view of each propagation steps is presented in Fig. 11(c).

Adiabatic load

In adiabatic crack, crack surfaces are assumed as thermally insulated and specific temperatures are imposed on top and bottom faces of the cuboid. These temperatures are different in magnitude ( T_{\mathrm{1}} and T_{\mathrm{2}} ) so that a uniform heat flux is maintained across crack surfaces. In thermo-elastic adiabatic crack problems, both temperature and displacement fields are discontinuous across the crack surface whereas heat flux is singular at the crack tip. The approximation of temperature and displacement fields can be written as,

(39a) T^h(\mathrm{x})={\displaystyle \sum _{j=\mathrm{1}}^n{N}_j(\mathbf{x})\left[T_j+\underset{j\in n_r}{\underbrace{(H(\mathbf{x})-H({\mathbf{x}}_j))a_j}}+\underset{j\in n_A}{\underbrace{\sum _{\alpha =\mathrm{1}}^{\mathrm{4}}({\beta }_{\alpha }(\mathbf{x})-{\beta }_{\alpha }({\mathbf{x}}_j))b_j^{\alpha }}}\right]},

(39b) u^h(\mathrm{x})={\displaystyle \sum _{j=\mathrm{1}}^n{N}_j(\mathbf{x})\left[u_j+\underset{j\in n_r}{\underbrace{(H(\mathbf{x})-H({\mathbf{x}}_j))a_j}}+\underset{j\in n_A}{\underbrace{\sum _{\alpha =\mathrm{1}}^{\mathrm{4}}({\beta }_{\alpha }(\mathbf{x})-{\beta }_{\alpha }({\mathbf{x}}_j))b_j^{\alpha }}}\right]}.

Heaviside function is used to capture discontinuity in temperature and displacement fields. Four crack tip asymptotic functions are used to model the singularities in heat flux and stress fields.

A cuboid with penny shape crack is considered for the numerical simulation under cyclic adiabatic thermal load. The complete problem is shown in Fig. 7. A constant amplitude cyclic heat flux is applied perpendicular to the crack surface. Crack surfaces are assumed to be thermally insulated. For maximum cyclic thermal load, specified temperatures i.e., ( T_{\mathrm{1}}=\mathrm{200}°C , T_{\mathrm{2}}=-\mathrm{200}°C ) are imposed on the top ( T_{\mathrm{1}} ) and bottom ( T_{\mathrm{2}} ) surfaces of cuboid whereas for minimum thermal load, these specified temperatures become zero. In this way, constant cyclic heat flux is applied across the crack surface. Initially, heat equation is solved to obtain temperature profile over the entire domain; this temperature data is then used as load input to solve elastic equation. Top and bottom surfaces are constrained in x, y and z directions. Fatigue life obtained using Paris law is presented in Fig. 12. The fatigue life and critical crack length (maximum length between opposite points at crack front along crack surface) found as 48096×10² cycles and 0.9790 m respectively. Crack surface contours at first and final crack increments are presented in Fig. 13. Front and top views of crack contours are presented in Fig. 14. These contours shows that crack grows in non-planer direction and move towards the bottom surface of the cuboid.

Isothermal load

In isothermal crack modelling, temperature ( T=\bar{T} ) is imposed at the crack surface {\mathrm{\Gamma }}_c using Lagrange multiplier approach. First, elements intersected by the crack surface are identified using the level set approach, and then temperature is directly specified at these elemental nodes using Lagrange multiplier approach. These elements have been subdivided into several tetrahedrons above and below the crack surface to accurately perform the numerical integration. The approximation of temperature field and displacement field for isothermal crack model may be written as [ 60],

(40) T(\mathbf{x})={\displaystyle \sum _{j=\mathrm{1}}^n{N}_j}(\mathbf{x})\left[T_j+\underset{j\epsilon n_t}{\underbrace{\sqrt{r}\cos \frac{\theta }{\mathrm{2}}{c}_j}}\right],

(41) u(\mathbf{x})={\displaystyle \sum _{j=\mathrm{1}}^n{N}_j}(\mathbf{x})\left[u_j+\underset{j\epsilon n_s}{\underbrace{[H(\mathbf{x})-H({\mathbf{x}}_{\mathbf{I}})]a_j}}+\underset{j\epsilon n_t}{\underbrace{\sum _{\alpha =\mathrm{1}}^{\mathrm{4}}[{\beta }_{\alpha }(\mathbf{x})-{\beta }_{\alpha }({\mathbf{x}}_{\mathbf{I}})]b_j^{\alpha }}}\right].

A cuboid with penny shape crack is numerically solved under cyclic isothermal load. The complete problem domain is shown in Fig. 7. A constant amplitude cyclic radial heat flux is applied from domain outer surface to crack surface. For maximum cyclic thermal load, specified temperatures i.e., T_{\mathrm{1}}=\mathrm{30}°C are imposed on all outer surfaces of cuboid, and crack surface is imposed with T_{\mathrm{2}}=\mathrm{0}°C whereas for minimum thermal load, these specified temperatures become zero on cuboid outer surface and crack surface. In this way, constant amplitude cyclic radial heat flux is applied within the cracked domain. Initially, heat equation is solved to obtain temperature profile over the entire domain. This temperature data is then used as load input to solve elastic equation. Top and bottom surfaces are constrained in x, y and z directions. Domain based thermal interaction integral approach has been applied to obtain the mixed mode SIFs. These SIFs are further used to obtain fatigue life of cracked domain using Paris law. The life cycles verses crack length are presented in Fig. 15. The fatigue life and corresponding critical crack length (maximum length between opposite points at crack front along crack surface) are found as 33400×10² cycles and 0.8681 m respectively. The crack surface contours at initial and final crack increment are presented in Fig. 16. Front and top views of crack increment contours are presented in Fig. 17. These contours show that the crack grows in its plane and remains circular after each increment.

Inclined elliptical shape crack in finite domain

In this case, a cuboid (\mathrm{2}\mathrm{m}\times \mathrm{2}\mathrm{m}\times \mathrm{2}\mathrm{m}) with inclined elliptical surface crack has been solved under thermal shock, adiabatic and isothermal loads. The complete model of this problem is shown in Fig. 18. For this non-planer crack surface, semi minor axis and semi major axis are taken as 0.15 m and 0.3 m. The crack surface is tilted by {\mathrm{30}}^o with the horizontal. The complete cuboid domain uniformly discretized by eight node brick elements.

Thermal shock load

A uniform cyclic thermal shock of \mathrm{\Delta }T=-\mathrm{20}°C has been applied over the entire domain. The top and bottom surfaces are constrained in x, y and z directions. Standard Paris law is used to predict the fatigue life of the problem, and the results are presented in Fig. 19. The critical fatigue life is shown by star mark in the life cycle verses crack length plot. The fatigue life and critical crack length (maximum length between two opposite points at the crack front) obtained by XFEM are found as 36476×10² cycles and 0.9144 m respectively. The crack surface contours are presented in Fig. 20 for initial and critical crack sizes. These contours show that crack growth is changing to circular shape. The front view and top view of crack propagation contours for each increment steps are presented in Fig. 21. These views show that the major axis of crack surface grows in planer direction (x-direction) and minor axis grows to become circular in shape.

Adiabatic load

An inclined elliptical shape crack (shown in Fig. 18) subjected to cyclic adiabatic thermo-elastic load is simulated by XFEM. Crack surfaces are assumed to be thermally insulated. To apply cyclic adiabatic load, specified temperatures i.e., T_{\mathrm{1}}=\mathrm{200}°C and T_{\mathrm{2}}=-\mathrm{200}°C are imposed on the top and bottom surfaces of the cuboid. For second set of thermal boundary condition, these specified temperatures become zero. Cyclic adiabatic thermal load is applied until crack becomes critical. A linear crack increment is used for crack growth simulation. Paris model is used to predict fatigue life of the model, and results are reported in Fig. 22. The fatigue life and critical crack length (maximum length between opposite points at crack front along crack surface) are found as 31624×10² cycles and 0.9303 m respectively. The crack surface contours at initial and final steps are presented in Fig. 23. Front and top views of crack contours are presented in Fig. 24. The crack growth contours shows that crack propagate in non-planer manner. The major axis end grows towards upper face of the cuboid whereas another end moves towards the bottom face of the cuboid.

Isothermal load

A cuboid with inclined elliptical shape crack as shown in Fig. 18 is subjected to cyclic isothermal crack load. A constant amplitude cyclic radial heat flux from outer surface of problem domain to the crack surface is applied. For maximum cyclic thermal load, specified temperature i.e., T_{\mathrm{1}}=\mathrm{30}°\mathrm{C} is imposed on all outer surfaces of the cuboid and T_{\mathrm{2}}=\mathrm{0}°\mathrm{C} is maintained at the crack faces whereas for minimum thermal load, these specified temperatures become zero on outer surfaces and crack faces. A constant amplitude cyclic radial heat flux has been applied until crack becomes critical. A linear crack increment is used for crack growth simulation. In post processing phase, Paris model is applied to predict fatigue life and the results are presented in Fig. 25. The fatigue life and critical crack length (maximum length between opposite points at crack front along crack surface) are found as 46904×10² cycles and 1.0780 m respectively. Crack surface contours at initial and final steps are presented in Fig. 26. Front and top views of crack contours are presented in Fig. 27. The crack growth contours shows that the crack grows in such a way that finally it becomes circular. In subsequent steps, the crack grows in x-y plane (horizontally).

Arbitrary spline shape crack in finite domain

In general, all real life crack problems are 3-D in nature, and do not have regular shapes, therefore, a cuboid with arbitrary spline shape crack lying at centre has been taken for simulation. The outer crack surface has been created manually. Total 26 arbitrary virtual points are created, and joined to each other to form arbitrary spline shape. After Delaunay triangulation of closed cracked boundary, an arbitrary crack surface has been formed as shown in Fig. 28.

Thermal shock load

The arbitrary shape crack geometry is subjected to uniform cyclic thermal shock load of \mathrm{\Delta }T=-\mathrm{20}°\mathrm{C} over the entire domain. The top and bottom surfaces of the domain are constrained in x, y and z directions as shown in Fig. 18, which induces mode-I effect loading on the crack. Standard Paris law has been implemented to predict the life cycle of the problem at each crack increment. The numerically obtained fatigue life is presented in Fig. 29. The fatigue life and critical crack length corresponding to maximum length between the opposite points at the crack front are found as 48021×10² cycles and 1.1181 m respectively. The crack surface contours for initial and final crack increments are presented in Fig. 30. The front and top views for all the increments are presented in Fig. 31. These contours show that crack grows in its own plane and try to become circular.

Adiabatic load

In this case, arbitrary crack surface is assumed as thermally insulated. To impose cyclic adiabatic load, specified temperatures i.e., T_{\mathrm{1}}=\mathrm{200}°\mathrm{C} is imposed on the top surface while T_{\mathrm{2}}=-\mathrm{200}°\mathrm{C} is applied on the bottom surface of the cuboid. For second set of thermal boundary condition, these specified temperatures become zero. The complete problem domain is shown in Fig. 28. A linear crack increment is given to all 26 virtual points using Eq. (40). Fatigue life predicted by Paris law is presented in Fig. 32. The fatigue life and critical crack length (maximum length between opposite points at crack front along crack surface) are found as 59123×10² cycles and 0.9654 m respectively. Crack surface contours for initial and final crack increments are presented in Fig. 33. The front and top view of crack contours are presented in Fig. 34. These crack growth contours show that the crack grows in its own plane towards the bottom surface of the cube.

Isothermal load

In the last sub-section, an arbitrary spline shape crack domain is subjected to constant amplitude cyclic isothermal crack load as shown in Fig. 28. For maximum cyclic thermal load, T_{\mathrm{1}}=\mathrm{30}°\mathrm{C} is applied on all the outer surfaces of cuboid and crack surface are kept at T_{\mathrm{2}}=\mathrm{0}°\mathrm{C} whereas to impose minimum thermal load, these specified temperature values become zero. Total twenty five J-domain points are created on crack front to calculate the SIFs. A linear crack increment is given to these points according to Eq. (40). After each crack increment, the new points are joined by spline which represents a new crack front. The numerically obtained SIFs are used to predict the fatigue life of cracked domain using Paris law, and the results are presented in Fig. 35. The fatigue life and corresponding critical crack length (maximum length between opposite points at crack front along crack surface) are found as 57172×10² cycles and 0.9891 m respectively. The crack surface contours for initial and final crack increments are presented in Fig. 36. The front and top views for all crack increments are presented in Fig. 37. The arbitrary shape crack simulation shows that the crack grows in own plane i.e., x-y plane and tries to become circular to attain the state of minimum potential energy.

Conclusions

The fatigue crack growth simulations of three-dimensional cracks have been performed by XFEM for three types of thermal boundary conditions namely thermal shock, adiabatic and isothermal. The results obtained by XFEM were found in good agreement with analytical/ reference solutions. The adiabatic load (except non-planer crack geometry) provides more fatigue life as compared to shock and isothermal loads. Thermal shock and isothermal loads induce mode-I effect whereas mode-II effect is more dominant in case of adiabatic load. Planer crack growth is observed in case of thermal shock and isothermal loads whereas adiabatic load produces for non-planer crack growth due to the presence of mode-II. It was also found that the division of crack front into several piecewise curve parts makes the computations easier as compared to iterative schemes. The capability and versatility of XFEM is well demonstrated by solving non-planer and arbitrary shape crack growth problems.

References

[1] Mukhopadhyay N K, Maiti S K, Kakodkar A. Effect of modelling of traction and thermal singularities on accuracy of SIFS computation through modified crack closure integral in BEM. Engineering Fracture Mechanics, 1999, 64(2): 141–159

[2] Dai H L, Wang X. Thermo-electro-elastic transient responses in piezoelectric hollow structures. International Journal of Solids and Structures, 2005, 42(3–4): 1151–1171

[3] Wang T H, Lai Y S. Submodeling analysis for path-dependent thermomechanical problems. Journal of Electronic Packaging, 2005, 127(2): 135–140

[4] Balderrama R, Cisilino A P, Martinez M. Boundary element method analysis of three-dimensional thermoelastic fracture problems using the energy domain integral. Journal of Applied Mechanics, 2006, 73(6): 959–969

[5] Zhong X C, Lee K Y. A thermal-medium crack model. Mechanics of Materials, 2012, 51: 110–117

[6] Barsoum R S. Triangular quarter-point elements as elastic and perfectly plastic crack tip elements. International Journal for Numerical Methods in Engineering, 1977, 11(1): 85–98

[7] Nash Gifford L Jr, Hilton P D. Stress intensity factors by enriched finite elements. Engineering Fracture Mechanics, 1978, 10(3): 485–496

[8] Nikishkov G P, Atluri S N. An equivalent domain integral method for computing crack tip integral parameters in non-elastic, thermo-mechanical fracture. Engineering Fracture Mechanics, 1987, 26(6): 851–867

[9] Rhee H C, Salama M M. Mixed-mode stress intensity factor solutions of a warped surface flaw by three-dimensional finite element analysis. Engineering Fracture Mechanics, 1987, 28(2): 203–209

[10] Carter B, Chen C S, Chew L P, Chrisochoides N, Gao G R, Heber G, Ingraffea A R, Krause R, Myers C, Nave D, Pingali K, Stodghill P, Vavasis S, Wawrzynek P A. Parallel FEM simulation of crack propagation—challenges, status and perspectives, Parallel and Distributed Processing Lecture Notes in Computer Science, 2000, 1800: 443–449

[11] Prasad N N V, Aliabadi M H, Rooke D P. Incremental crack growth in thermoelastic problems. International Journal of Fracture, 1994, 66(3): 45–50

[12] Prasad N N V, Aliabadi M H, Rooke D P. The dual boundary element method for transient thermoelastic crack problems. International Journal of Solids and Structures, 1996, 33(19): 2695–2718

[13] Keppas L K, Anifantis N K. BEM prediction of TBC fracture resistance. Engineering against Fracture Proceeding of the 1st Conference, 2009, 551–560

[14] Pant M, Singh I V, Mishra B K. Evaluation of mixed mode stress intensity factors for interface cracks using EFGM. Applied Mathematical Modelling, 2011, 35(7): 3443–3459

[15] Pathak H, Singh A, Singh I V. Fatigue crack growth simulations of homogeneous and bi-material interfacial cracks using element free Galerkin method. Applied Mathematical Modelling, 2014, 38(13): 3093–3123

[16] Melenk J, Babuska I. The partition of unity finite element method: Basic theory and applications. Computer Methods in Applied Mechanics and Engineering, 1996, 139(1–4): 289–314

[17] 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

[18] Li S, Liu W K. Meshfree Particle Methods. Springer, ISBN 978-3-540-71471-2, 2004

[19] Moral PD, Doucet A. Particle methods: An introduction with applications. RR-6991, 2009, 46

[20] 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

[21] Rabczuk T, Song J H, Belytschko T. Simulations of instability in dynamic fracture by the cracking particles method. Engineering Fracture Mechanics, 2009, 76(6): 730–741

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

[23] Xu M, Gracie R, Belytschko T. Multiscale Modeling with Extended Bridging Domain Method. Bridging the Scales in Science and Engineering, Oxford Press, 2009, 1–32

[24] Bhardwaj G, Singh I V, Mishra B K, Bui T Q. Numerical simulation of functionally graded cracked plates using NURBS based XIGA under different load and boundary conditions. Composite Structures, 2015, 126: 347–359

[25] 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

[26] Holl M, Loehnert S, Wriggers P. An adaptive multiscale method for crack propagation and crack coalescence. International Journal for Numerical Methods in Engineering, 2013, 93(1): 23–51

[27] 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

[28] 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

[29] Sukumar N, Möes N, Moran B, Belytschko T. Extended finite element method for three-dimensional crack modeling. International Journal for Numerical Methods in Engineering, 2000, 48: 1549–1570

[30] Sukumar N, Chopp D L, Moran B. Extended finite element method and fast marching method for three-dimensional fatigue crack propagation. Engineering Fracture Mechanics, 2003, 70(1): 29–48

[31] Daux C, Möes N, Dolbow J, Sukumar N, Belytschko T. Arbitrary branched and intersecting cracks with the extended finite element method. International Journal for Numerical Methods in Engineering, 2000, 48: 1741–1760

[32] Sukumar N, Chopp D L, Möes N, Belytschko T. Modelling 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

[33] Liu X Y, Xiao Q Z, Karihaloo B L. XFEM for direct evaluation of mixed mode SIFs in homogeneous and bi-materials. International Journal for Numerical Methods in Engineering, 2004, 59(8): 1103–1118

[34] Ventura G, Budyn E, Belytschko T. Vector level sets for description of propagating cracks in finite elements. International Journal for Numerical Methods in Engineering, 2003, 58(10): 1571–1592

[35] Möes N, Gravouil A, Belytschko T. Non-planar 3D crack growth by the extended finite element and level sets-Part-I: Mechanical model. International Journal for Numerical Methods in Engineering, 2002, 53: 2549–2568

[36] Möes N, Gravouil A, Belytschko T. Non-planar 3D crack growth by the extended finite element and level sets-Part II: Level set update. International Journal for Numerical Methods in Engineering, 2002, 53: 2569–2586

[37] Shi J, Chopp D, Lua J, Sukumar N, Belytschko T. Abaqus implementation of extended finite element method using a level set representation for three-dimensional fatigue crack growth and life predictions. Engineering Fracture Mechanics, 2010, 77(14): 2840–2863

[38] Unger J F, Eckardt S, Könke C. Modelling of cohesive crack growth in concrete structures with the extended finite element. Computer Methods in Applied Mechanics and Engineering, 2007, 196(41–44): 4087–4100

[39] Asferg J L, Poulsen P N, Nielsen L O. A consistent partly cracked XFEM element for cohesive crack growth. International Journal for Numerical Methods in Engineering, 2007, 72(4): 464–485

[40] Zhang X D, Bui Q T. A fictitious crack XFEM with two new solution algorithms for cohesive crack growth modelling inn concrete structures. Engineering Computations, 2015, 32(2): 473–497

[41] Chopp D L, Sukumar N. Fatigue crack propagation of multiple coplanar cracks with the coupled extended finite element/fast marching method. International Journal of Engineering Science, 2003, 41(8): 845–869

[42] Giner E, Sukumar N, Denia F D, Fuenmayor F J. Extended finite element method for fretting fatigue crack propagation. International Journal of Solids and Structures, 2008, 45(22–23): 5675–5687

[43] 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

[44] Bhattacharya S, Singh I V, Mishra B K. Fatigue life estimation of functionally graded materials using XFEM. Engineering with Computers, 2013, 29(4): 427–448

[45] Bhattacharya S, Singh I V, Mishra B K. Mixed-mode fatigue crack growth analysis of functionally graded materials by XFEM. International Journal of Fracture, 2013, 183(1): 81–97

[46] Bhattacharya S, Singh I V, Mishra B K, Bui T Q. Fatigue crack growth simulations of interfacial cracks in bi-layered FGMs using XFEM. Computational Mechanics, 2013, 52(4): 799–814

[47] Bhattacharya S, Singh I V, Mishra B K. Fatigue life simulation of functionally graded materials under cyclic thermal load using XFEM. International Journal of Mechanical Sciences, 2014, 82: 41–59

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

[49] 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

[50] Bui Q T, Zhang Ch. Extended finite element simulation of stationary dynamic cracks in piezoelectric solids under impact loading. Computational Materials Science, 2012, 62: 243–257

[51] Bui Q T, Zhang Ch. Analysis of generalized dynamic intensity factors of cracked magnetoelectrostatic solids by XFEM. Finite Elements in Analysis and Design, 2013, 69: 19–36

[52] Sharma K, Bui Q T, Zhang Ch, Bhargava R R. Analysis of a subinterface crack in piezoelectric bimaterials with the extended finite element method. Engineering Fracture Mechanics, 2013, 104: 114–139

[53] Liu P, Yu T T, Bui Q T, Zhang Ch. Transient dynamic crack analysis in non-homogeneous functionally graded piezoelectric materials by the X-FEM. Computational Materials Science, 2013, 69: 542– 558

[54] Liu P, Yu T T, Bui Q T, Zhang Ch, Xu Y P, Lim C W. Transient thermal shock fracture analysis of functionally graded piezoelectric materials by the extended finite element method. International Journal of Solids and Structures, 2014, 51(11–12): 2167–2182

[55] Liu P, Bui Q T, Zhu D, Yu T T, Wang J W, Yin S H, Hirose S. Buckling failure analysis of cracked functionally graded plates by a stabilised discrete shear gap extended 3-noded triangular plate element. Composites. Part B, Engineering, 2015, 77: 179–193

[56] Yu T T, Bui Q T, Liu P, Zhang Ch, Hirose S. Interface dynamic impermeable cracks analysis in dissimilar piezoelectric materials under coupled electromechanical loading with the extended finite element method. International Journal of Solids and Structures, doi: 10.1016/j.ijsolstr.2015.03.037

[57] Pathak H, Singh A, Singh I V. Fatigue crack growth simulations of 3-D problems using XFEM. International Journal of Mechanical Sciences, 2013b, 76: 112–131

[58] Moës N, Dolbow J, Belytchsko T. A finite element method for crack growth without remeshing. International Journal for Numerical Methods in Engineering, 1999, 46: 131–150

[59] Sih G C. On singular character of thermal stress near a crack tip. Journal of Applied Mechanics, 1962, 29(3): 587–589

[60] Duflot M. The extended finite element method in thermo-elastic fracture mechanics. International Journal for Numerical Methods in Engineering, 2008, 74(5): 827–847

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

[62] Mohammadi S. Extended finite element method for fracture analysis of structures. Singapore: Blackwell Publishing, ISBN-978-1-4051-7060-4, 2008

[63] Singh I V, Mishra B K, Bhattacharya S. XFEM simulation of cracks, holes and inclusions in functionally graded materials. International Journal of Mechanics and Materials in Design, 2011, 7(3): 199–218

[64] Moran B, Shih C. A general treatment of crack tip contour integrals. International Journal of Fracture, 1987, 35(4): 295–310

[65] Banks-Sills L, Dolev O. The conservative M-integral for thermo-elastic problems. International Journal of Fracture, 2004, 125(1): 149–170

[66] 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

[67] Das B R. Thermal stresses in a long cylinder containing a penny shaped crack. International Journal of Engineering Science, 1968, 6(9): 497–516

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag Berlin Heidelberg