Modeling and simulation of kinked cracks by virtual node XFEM

Accepted Manuscript Modeling and simulation of kinked cracks by virtual node XFEM Sachin Kumar, I.V. Singh, B.K. Mishra, Timon Rabczuk PII: DOI: Reference:

S0045-7825(14)00389-2 http://dx.doi.org/10.1016/j.cma.2014.10.019 CMA 10413

To appear in:

Comput. Methods Appl. Mech. Engrg.

Received date: 16 May 2014 Revised date: 16 October 2014 Accepted date: 17 October 2014 Please cite this article as: S. Kumar, I.V. Singh, B.K. Mishra, T. Rabczuk, Modeling and simulation of kinked cracks by virtual node XFEM, Comput. Methods Appl. Mech. Engrg. (2014), http://dx.doi.org/10.1016/j.cma.2014.10.019 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Modeling and Simulation of Kinked Cracks by Virtual Node XFEM Sachin Kumar#, I. V. Singh#, B. K. Mishra#, Timon Rabczuk$ #Department of Mechanical and Industrial Engineering, Indian Institute of Technology Roorkee, Uttarakhand, India, Phone: +91-1332-285888, Fax: +91-1332-285665, Email: [email protected] $Institute of Structural Mechanics, Bauhaus-Universitét Weimar, Germany, Email: [email protected]

ABSTRACT In the present work, virtual node extended finite element method (VNXFEM) is proposed for modeling and simulation of kinked cracks in a single element. The kinked crack growth in an element is performed using the concept of virtual nodes to model and improve the accuracy of the solution. The virtual nodes are enriched with additional degrees of freedom. In the proposed approach, an actual tip element is divided into one virtual tip element and one virtual split element with the help of virtual nodes. Special five-node transition elements are used to ensure the continuity in the displacement fields. A polygon law is proposed to determine the positions of the Gauss points with respect to a crack face in the virtual split elements for the purpose of integration. Several crack growth problems in homogeneous and bi-materials are solved to verify the proposed method. The results obtained by VNXFEM are compared with those available in literature and standard XFEM solutions. Keywords: Virtual node XFEM; virtual split element (VSE); virtual tip element (VTE); transition elements; discontinuities; stress intensity factors 1. INTRODUCTION Over the past few decades, the finite element method (FEM) has been widely used as a numerical tool for the simulation of fracture mechanics problems. Various software packages have been developed based on FEM to study the various problems of fracture mechanics. However, the modelling of complex problems (propagating cracks and moving discontinuities) becomes quite difficult in standard FEM as it requires conformal mesh to model the crack growth or moving discontinuities. To overcome these problems, several methods have been proposed by the researchers such as boundary elements method (Portela et al., 1991), boundary collocation method (Newman, 1971), meshfree method (Belytschko et al., 1994) and extended finite element method (Belytschko and Black, 1999; Daux et al.,

1

2000). Among these, extended finite element method (XFEM) is found to be more accurate and powerful method to analyze fracture mechanics problems. In the last decade, the XFEM has been used extensively by many researchers for the simulation of fracture mechanics problems. In XFEM, crack propagation problems have been solved without re-meshing by adding few enrichment functions into the standard finite element approximations. In 1999, Moes and his coworkers (1999) first used the partition of unity (Melenk and Babuska, 1996) to enrich the standard finite element approximation by incorporating discontinuous fields near the crack face and asymptotic stress fields near the crack tip. In 2000, Daux et al., (2000) applied the XFEM to model the cracks with multiple branches, holes and cracks emanating from the holes. They enriched the standard displacement approximation by discontinuous fields through the partition of unity. Zi and Belytschko (2003) modeled the elements intersected by the crack (split and tip elements) with the signed distance function so that no blending of the local partition of unity is required. In 2003, Sukumar and Prevost (2003) carried out the quasi-static crack growth modeling without remeshing in isotropic and bi-material media using XFEM. They also extended the XFEM to simulate the elastic-plastic fracture mechanics problems. Elguedj et al., (2006) introduced the new crack tip enrichment functions, based on Hutchinson-Rice-Rosengren (HRR) fields, to predict elastic-plastic fatigue crack growth in homogeneous two-dimensional solids. Anahid and Khoei (2008) coupled the XFEM with the Lagrangian formulation in order to simulate the arbitrary interfaces with large deformations. Furthermore, Khoei et al., (2009) extended this Lagrangian-XFEM to model the large deformation with contact nonlinearity for arbitrary interfaces in two and three dimensional problems. In 2014, Kumar et al., (2014a) used the crack tip enrichment functions (Elguedj et al., 2006) to simulate the stable crack growth in ductile materials. In recent years, a new set of orthotropic crack tip enrichment functions have been developed for XFEM and then effectively used for the static (Asadpoure et al., 2006 and Asadpoure and Mohammadi, 2007) and dynamic analysis (Motamedi and Mohammadi, 2010a; Motamedi and Mohammadi, 2010b) of fracture problems. In 2008, Fries (2008) proposed a new corrected XFEM to improve the accuracy of the standard XFEM by introducing the ramp function in the blending elements such that the partition of unity is satisfied everywhere in the domain. Since then, a lot of work has been carried out by many investigators (Gracie et al., 2008; Tarancon et al., 2009; Shibanuma and Utsunomiya, 2009) on blending elements. In 2008, Rabczuk et al., (2008) proposed a new crack tip element for 2

the phantom node method to solve arbitrary cohesive cracks using three-node triangular and four-node quadrilateral elements. They concluded that the proposed method can be effectively used for dynamic fracture problems with one-point quadrature scheme. Natarajan et al., (2009) suggested a new approach for the numerical integration over arbitrary polygon domains using Schwarz-Christoffel mapping and midpoint quadrature rule, which eliminate the need of two-level isoparametric mapping. In 2011, Ashari and Mohammadi (2011) proposed a set of crack tip enrichment functions for XFEM to analyze the interlaminar cracks in orthotropic bi-material domains. In 2011, Fries et al., (2011) proposed the concept of hanging nodes for the local mesh refinement in XFEM. They performed the simulations using hanging nodes with additional degrees of freedom and hanging nodes without additional degrees of freedom. They achieved higher convergence rate for hanging nodes with additional degrees of freedom. Further, Cheng and Fries (2012) explored the concept of hanging nodes for the simulation of twophase incompressible flow problems. The concept of phantom-node is also applied for threenode shell elements to model the cracks by Chau-Dinh et al., (2012). Chen et al., (2012) proposed an edge based strain smoothing technique to eliminate the requirement of elements subdivision by transforming the interior integration into boundary integration. In 2013, Bac et al., (2013) proposed a combined approach of phantom node with edge-based strain smoothing for the solution of linear elastic fracture mechanics problems. In the proposed method, the crack is modeled by adding phantom nodes. In 2012, extended finite element method has been extended to simulate moving crack problems in composites (Motamedi and Mohammadi, 2012), orthotropic functionally graded materials (Bayesteh and Mohammadi, 2013) and carbon nanotube reinforced concrete (Eftekhari et al., 2014). In the same year, Singh et al., (2012) implemented the XFEM to evaluate the fatigue life of homogeneous plate in the presence of multiple discontinuities under cyclic load, and found that the presence of multiple discontinuities significantly affect the fatigue life of materials. Prange et al., (2012) proposed a technique to control the discretization error in XFEM based on the approach of Zienkiewicz and Zhu (1987). Loehnert et al., (2012) proposed a discretization error controlled adaptive multiscale approach for the simulation of micro-defects in a macroscopic component. Bouhala et al., (2013) determined the crack tip enrichment functions for a crack terminating at the bimaterial interface. In 2013, Pathak et al., (2013) used the XFEM for the fatigue crack growth simulations of bi-material interfacial cracks.

3

Thus, it can be stated that XFEM provides a mesh independent approximation for non smooth solutions. In XFEM, crack propagation problems were modeled with an assumption that the crack tip must fall into the new element during crack growth on a curved path (Singh et al., 2012; Pathak et al., 2013). This assumption imposes limitations on the crack growth modeling by XFEM. In XFEM, very small crack growth increments are mandatory to model actual crack growth. Therefore, a very fine mesh is required for the crack growth simulation, which is computationally quite expensive. On the other hand, for a coarse mesh, the crack growth increment size becomes large then the piecewise linear segments cannot accurately represent a real crack path. Thus, there is a need to develop a new scheme in XFEM for the modeling of crack growth within an element. Therefore, in the present work, a new approach based on virtual node concept i.e. VNXFEM has been proposed to model the crack growth within an element. In past, Fries et al., (2011) obtained the higher convergence rate using a concept of hanging nodes with additional degrees of freedom. The additional degrees of freedom associated with virtual/hanging nodes contribute to the system matrix, hence improves the accuracy by local mesh refinement. In VNXFEM, to model the kink within an element, actual crack tip element has been divided into one virtual split element and one virtual tip element with the help of virtual nodes. These virtual nodes are enriched with additional degrees of freedom, depending upon their positions. The shape functions for all regular and virtual nodes are found in such a way that they follow partition of unity. The special five-node transition elements (Gupta 1978; Huang and Xie, 2010; Kumar et al., 2014b) are used to connect the virtual elements with the standard finite elements to ensure the continuity in the displacement field. The virtual nodes, virtual split element and virtual tip element keep changing with the advancement of crack tip position. For the integration of virtual split elements, level set method cannot decide the position of the Gauss points accurately with respect to the kinked crack. Therefore, a new scheme based on Polygon law, is proposed to obtain the location of Gauss points with respect to the crack in the split elements. In this scheme, a polygon is formed in anticlockwise direction for all the split elements with the help of kinks and nodes coordinates. Gauss point location is obtained with respect to the polygon. If the Gauss point is lying inside the polygon, means, it is above the crack; otherwise it is below the crack. Thus, the proposed approach can handle complex curved crack growth problems with extremely coarse mesh.

4

This paper is structured as follows: Section 2 describes the governing equations and XFEM formulation of linear elastic materials. Section 3 describes the modelling of crack growth within an element and a procedure to divide the tip element into virtual split elements (VSE) and virtual tip element (VTE). The shape function generation for the transition elements and numerical integration is also explained in Section 3. The formulation to compute the bi-material stress intensity factors is explained in Section 4. Section 5 describes the procedure to compute the crack growth direction and fatigue life. Several numerical examples of crack growth in homogeneous and bi-materials are presented in Section 6 to verify the VNXFEM. Section 7 provides conclusions and an outline of the future work. 2. NUMERICAL FORMULATION 2.1 Governing Equations Consider a domain ( Ω ) as shown in Figure 1, which is partitioned into displacement ( Γ u ) , traction ( Γ t ) and crack surface ( Γ c ) boundaries. The strong form of the equilibrium equations and boundary conditions are given as (Sukumar et al., 2001; Singh et al., 2012),

~ ∇. σ + b = 0 in Ω

(1a)

σ ⋅ nˆ = t

on Γ t

(1b)

σ ⋅ nˆ = 0

on Γ c

(1c)

u=u

on Γ u

(1d)

~ where, σ is the Cauchy stress tensor, u is the displacement field vector, b is the body force vector and nˆ is the unit outward normal. For small strains and displacements, the straindisplacement relation can be written as, ε = ε(u) = ∇ s u

(2)

In the above equation ∇ s is the symmetric part of the gradient operator and ε is the linear strain tensor. The constitutive relation for linear elastic material is given by Hook's law,

σ = Dε

(3)

where, D is the Hook's tensor. The weak form of the equilibrium equation can be written as, ~

∫ σ(u) : ε( v) dΩ = ∫ b.v dΩ + ∫ t.v dΓ

Ωe

Ωe

(4)

Γt

5

After substituting the trial and test functions into Eq. (4) and using the arbitrariness of nodal variations, the following discrete system of equations are obtained

Kd =f

(5)

where, K is the global stiffness matrix, d is the vector of nodal unknowns and f is the external force vector. 2.2 XFEM Formulation In two-dimension, at a particular node of interest x i , the displacement approximation with holes, inclusions and cracks can be written as (Belytschko and Black, 1999; Singh et al., 2012), ⎡ ⎤ ⎧ ⎫ ⎢ ⎥ ⎪⎪ 4 ⎪⎪ uh ( x ) = N i ( x ) ⎢ ui + [ Η ( x ) − H ( x i )] ai + ⎨ [ βα ( x ) − βα ( x i )] bαi ⎬ R( x ) + [V ( x ) − V ( x i )] c i + ϕ ( x ) d i + χ ( x )e i ⎥ 144 42444 3 1442443 123 123 ⎥ ⎢ i =1 =1 ⎪ α1444 i∈nc i∈nh i∈ni i∈nb 424444 3⎪ ⎢ ⎥ i∈nt ⎩⎪ ⎭⎪ ⎣ ⎦ n

∑

∑

(6)

where, ui is a nodal displacement vector of the standard finite element part at node i; n is the set of all nodes in the mesh; nc is the set of nodes associated with those elements which are completely cut by the crack; nt is the set of nodes associated with those elements which are partially cut by the crack; nh is the set of nodes associated with those elements which are cut by the hole; ni is the set of nodes associated with those elements which are cut by the inclusions; nb is the set of nodes associated with those elements which are cut by the interface. R(x) is a ramp function (Fries, 2008) used for the blending elements, and can be defined as, R ( x ) =

∑

i∈nenr

N i ( x ) , nenr is the set of crack tip enriched nodes (a typical discretized

geometry showing all types of nodes is shown in Figure 2). H (x) is the discontinuous enrichment function or Heaviside function, used to model a discontinuity in the displacement due to the presence of a crack, which takes +1 on one side and -1 on other side of the crack;

a i is the nodal enriched degree of freedom associated with H (x) . c i is the nodal enriched degree of freedom associated with V (x) (+1 outside the hole and zero inside the hole). d i is the nodal enriched degree of freedom associated with ϕ (x) , where ϕ (x) is a level set function used to model the inclusions. e i is the additional nodal degrees of freedom vector associated with signed distance function χ (x) used for simulating the material discontinuity.

bαi is the nodal enriched degrees of freedom vector associated with crack tip enrichment 6

βα (x) . The virtual crack tip element is modeled by asymptotic enrichment functions. In bimaterials, the stress singularity in the vicinity of the interfacial crack tip is oscillatory in nature. Thus, twelve crack tip enrichment functions (Sukumar et al., 2004) are required to model the oscillatory behavior. The time period of oscillations is very high, which creates the complexity in the modeling without much improvement in the results. Therefore, in the present work, four crack tip enrichment functions are used to capture the crack tip singularity in interfacial crack problems (Pant et al., 2011; Pathak et al., 2013). This reduces the size of the overall matrix as additional degrees of freedom are required is less. Consider a local polar coordinates system r and θ at the crack tip, the asymptotic crack tip enrichment functions can be written as,

θ ⎡ βα ( x ) = ⎢ r sin , 2 ⎣

θ

r cos , 2

r sin

θ 2

cos θ ,

θ ⎤ r cos cos θ ⎥ 2 ⎦

(7)

The elemental matrices, K and f in Eq. (5) are obtained by using Eq. (6), and are provided in Appendix. 3. MODELING OF KINKS IN AN ELEMENT In real life, the crack growth in a component or structure occurs on a curved path under actual loading conditions. To represent a real crack path, the kinked crack modeling is needed apart from the requirement of small size crack growth increment. In standard XFEM, the accurate modeling of crack growth demands kinked crack modeling and small size crack increments. Thus, the crack growth modeling in XFEM is based on the assumption that the crack tip must fall in the new element under mixed mode loading, which limits the size of crack growth increment. If the size of crack growth increment is small, then a linear crack segment can also model a crack path accurately but it requires very fine mesh, resulting in the increase of overall computational cost. On the other hand, if the crack growth increment size is large then the piecewise linear segments cannot accurately represent a real crack path. For the modeling of a kink crack, a typical triangulation procedure to determine the Gauss points is depicted in Figure 3. Figure 3(a) shows the standard triangulation (seven Gauss point per triangle) for a split element with a straight crack segment. Figure 3(b) shows an actual kinked crack in an element along with required sub-triangles but this type of triangulation (division into sub-triangles) imposes difficulty in generation of Gauss points as it requires many level sets depending upon the number of kinks within an element. To avoid this situation, the crack path with a kink is replaced by an approximate line segment as shown 7

in Figure 3(b), but this assumption does not model a real crack path. The standard triangulation procedure and Gauss point generation for a kinked element is depicted in Figure 3(c). To avoid the limitations of kinked crack modeling in an element, the VNXFEM has been proposed in the present work based on the concept of virtual nodes. In the VNXFEM, an actual tip element is divided into virtual tip and virtual split elements with the help of virtual nodes. At each virtual node, few additional degrees of freedom are imposed. These virtual nodes divide the actual tip element into VSE and VTE. Therefore, these virtual nodes form the common edge between VSE and VTE as shown in Figure 4 by dotted line. In Figure 4, virtual nodes (common nodes between VSE and VTE) are enriched by asymptotic crack tip enrichment functions. The detail of VNXFEM is explained in the following subsections. 3.1 Procedure to Divide Actual Tip Element into Two Virtual Elements The procedure to divide the current tip element into VSE and VTE is illustrated in Figure 4, where circle and rectangle represent a standard node and virtual node respectively. For the initial crack when no kink is present in the tip element, the tip element is divided into two virtual elements (VSE and VTE) at a specified distance behind the crack tip as shown in Figure 4(a). The specified distance is opted iteratively with a constraint of aspect ratio of VSE and VTE. Figure 4(b) shows a kinked crack within an element, here the actual tip element is divided into VSE and VTE with the help of virtual nodes. Now, consider a case where crack is increased by some specified amount within the same element as shown in Figure 4(c). In this case, two kinks are present in the same element. To model this, the virtual nodes, which were at the kink 1, will move to the kink 2 (previously known crack tip) as shown in Figure 4(c), and the actual tip element will be divided into VSE and VTE. In the VNXFEM, many kinks (as per need) present in an element can be modeled by the above explained procedure. Figure 4(a-f) depicts the procedure to split the current tip element into VSE and VTE for various kink possibilities in an element. In all these cases, if a current tip element possesses a kink then it is divided into VSE and VTE at the latest known kink position. Now, consider the case, where no kink is present in the actual tip element but several kinks are present in the actual split element as shown in Figure 5. In such a case, there is no need to introduce the virtual nodes in the actual split element. Only the virtual nodes are required to divide the actual tip element into VSE and VTE at a specified distance behind the crack tip by following the procedure as shown in Figure 4(a). Figure 5(b) presents the procedure to divide the actual tip element into VSE and VTE along with 8

triangulation procedure. After dividing the actual tip element into VSE and VTE, the VTE behaves as a standard four node isoparametric element. The shape functions of VTE are the standard four node quadrilateral element shape functions, which fulfills the partition of unity. Furthermore, the displacement and stress contours over a virtual tip element are shown in Figure 6 and Figure 7 respectively. 3.2 Integration Schemes for Discontinuous Elements For the elements partitioned by discontinuities, the standard Gauss quadrature cannot perform the integration accurately. Therefore, these elements are divided into sub-triangles (Fries, 2008; Singh et al., 2012) for the purpose of integration as shown in Figure 3. In these subtriangles, the integrand remains continuous. The procedure used to divide the elements into sub-triangles is illustrated in Figure 4 for different cases. As the displacement approximation is different from element to element therefore, the following integration orders are adopted in the various elements: • Tip element: 7 Gauss points in each sub-triangle • Split element: 5 Gauss points in each sub-triangle • Tip blending element: 16 Gauss points • Standard element: 4 Gauss points 3.3 Shape Functions for Transition Elements The XFEM with virtual/hanging nodes provides quite accurate results for crack growth problems. However, the main challenge lies in the construction of conforming partition of unity interpolation/shape functions. In 2006, Chao and Im (2006) constructed the shape functions for non-matching interface meshes. Consider a case illustrated in Figure 8(a), where transition elements are needed to patch the domain approximations. To construct the functions for a five noded transition element (Huang and Xie, 2010), first develop the shape function for node 5. N 5 (ξ ,η ) =

1 (1 − ξ 2 )(1 − η ) 2

(8)

The addition of node 5 will affect the node 1 and node 2, therefore a correction in the standard finite element shape functions is required. The original bilinear shape functions at the corner nodes 1 and 2 are modified in the following way, N1 (ξ ,η ) =

1 (1 − ξ )(1 − η ) − α1 × N 5 (ξ ,η ) 4

(9)

9

N 2 (ξ ,η ) =

1 (1 + ξ )(1 − η ) − α 2 × N 5 (ξ ,η ) 4

(10)

Now for Figure 8(b), the shape functions of transition element can be written in the similar way, N 6 (ξ ,η ) =

1 (1 − ξ 2 )(1 + η ) 2

(11)

N 3 (ξ ,η ) =

1 (1 + ξ )(1 + η ) − α 3 × N 6 (ξ ,η ) 4

(12)

N 4 (ξ ,η ) =

1 (1 − ξ )(1 + η ) − α 4 × N 6 (ξ ,η ) 4

(13)

where α1 , α 2 , α3 and α4 are the correction factors which are used to ensure that the final modified shape function should be zero at the virtual node, and satisfy the partition of unity condition as illustrated in Figure 9. For the purpose of plotting, the virtual node is considered at the midpoint of the transition element edge, and shape function contour plot for virtual node is shown in Figure 10. The values of correction factors can be obtained as, 1 4

1 4

1 4

1 4

α1 (ξ ,η ) = (1 − ξ )(1 − η ) , α 2 (ξ ,η ) = (1 + ξ )(1 − η ) α 3 (ξ ,η ) = (1 + ξ )(1 + η ) , α 4 (ξ ,η ) = (1 − ξ )(1 + η ) where, ξ and η are the local co-ordinates of the virtual node in the element. In transition elements, the displacement field variation is quadratic along the virtual node edge whereas, in virtual split and virtual tip elements, the displacement field variation is linear along the element edge. Therefore, displacement field is consistent at the nodes (virtual and standard) while it is inconsistent at the other points along the edge as shown in Figure 8(c). 3.4 Polygon Law to Find the Location of Gauss Points in Split Elements For the purpose of integration in discontinuous elements (split and tip elements), subtriangles are developed in discontinuous elements to generate the Gauss points. In split elements, several types of kinks may present as shown in Figure 3, hence the simplest way to model a kink is the use of a straight line segment as depicted in Figure 3(b). This assumption introduces some error into solution. The other option of kink modeling is the use of triangulation as given in Figure 3(c), but in this case, the level set function i.e. signed distance function cannot be easily applied to decide the position of the Gauss point with respect to crack face. Therefore, a new scheme based on the polygon law is proposed to 10

decide the position of the Gauss point in the split elements with respect to a crack face. A simple way to locate whether a particular Gauss point is inside or outside the polygon is to find the intersections of a ray (started from a Gauss point) with the polygon as shown in Figure 11. If the intersection count is odd number then a Gauss point lies inside the polygon, and if this count is even number then the Gauss point lies outside the polygon. For example in Figure 11, the ray originated from the point O cuts the polygon only once, hence point O lies inside the polygon, and the ray originated from point P cuts the polygon twice, hence the point P lies outside the polygon. Thus, a polygon is formed in anticlockwise direction for all the split elements with the help of kinks and nodes coordinates. The polygon formation for the kinked cracks in the split elements is illustrated in Figure 12 under various situations. For the modeling of a kink in these elements, the value of Heaviside function is chosen as,

⎧ +1, if point isinside the polygon H (x) = ⎨ ⎩ −1, if point is outside the polygon 4. COMPUTATION OF STRESS INTENSITY FACTORS For the computation of individual stress intensity factors (SIFs), the domain based interaction integral approach is an effective tool under mixed mode loading conditions. Johnson and Qu (2006) investigated that the J-integral remains path independent for bi-material interface crack if the material inhomogeneity exists in the direction parallel to crack surface. For two independent equilibrium states of the homogeneous cracked body, the interaction integral is defined as,

M

(1,2 )

=∫

A

(1) ⎡ (1) ∂ui( 2 ) ⎤ ∂q ( 2 ) ∂ui dA + σ ij − W (1,2 ) δ1 j ⎥ ⎢ σ ij ∂x1 ∂x1 ⎣ ⎦ ∂x j

(14)

where W (1, 2 ) interaction term associated with the actual and auxiliary states, q is the scalar quantity which is one at the crack tip and zero on the contour. σ ij and ε ij are the stress and strain fields respectively, 1 and 2 signify the actual state and auxiliary state. For the bimaterial interface cracked body, the interaction integral can be written as (Sukumar et al., 2004; Pathak et al., 2013),

M

(1,2 )

2

=∑∫ m =1

Am

(1) ⎡ (1) ∂ui( 2 ) ⎤ ∂q ( 2 ) ∂ui dA + σ ij − W (1,2 ) δ1 j ⎥ ⎢ σ ij ∂x1 ∂x1 ⎣ ⎦ ∂x j

(15)

where, m represents a particular material in the bi-material body. In LEFM, the interaction integral is related to the SIFs through the relation,

11

M (1,2 ) =

(

)

(16)

⎛1− β ⎞ 1 log ⎜ ⎟ 2π ⎝1+ β ⎠

(17)

2 K I(1) K I( 2 ) + K II(1) K II( 2 ) 2 E cosh πε *

where, 2 1 1 = + * E E1 E2

and

ε=

where, β is the second Dundurs parameter defined as,

β =

μ1 (κ 2 − 1) − μ2 (κ1 − 1) μ1 (κ 2 + 1) + μ2 (κ1 + 1)

⎧ 3 −ν i ⎪ κi = ⎨ 1 +ν i ⎪3 − 4ν i ⎩

(18)

( plane stress)

(19)

( plane strain )

where, μi , ν i and κ i are the shear modulus, Poisson's ratio and Kolosov constant respectively for corresponding material i (i =1,2). The mixed mode stress intensity factors can be obtained form Eq. (16) using K I(2) = 1 and K II(2) = 0 and vice versa. KI =

E * cosh 2 (πε ) (1) E * cosh 2 (πε ) (2) M , K II = M 2 2

(20)

5. COMPUTATION OF CRACK GROWTH DIRECTION AND FATIGUE LIFE In the present work, fatigue crack growth simulations are performed by VNXFEM under constant amplitude cyclic loading. The discrete equations are solved to obtain the displacements, and then values of stress intensity factors are extracted from Eq. (20). The range of SIF for constant amplitude cyclic loading is defined as,

ΔK = K max − K min

(21)

where, K max and K min are the stress intensity factors corresponding to maximum and minimum applied loads respectively. In real life problems, a crack path is curved in nature. Therefore, in the present simulations, it is modeled through many small line segments. The maximum principal stress criterion is used to determine the direction of the crack growth. The equivalent mode-I SIF and the direction of crack growth θ c at each crack increment is obtained by the following expression (Erdogan and Sih, 1963), ⎛θ ⎞ ⎛θ ⎞ ⎛θ ⎞ ΔK Ieq = ΔK I cos3 ⎜ c ⎟ − 3ΔK II cos2 ⎜ c ⎟ sin ⎜ c ⎟ ⎝2⎠ ⎝2⎠ ⎝2⎠

12

(22)

⎛ K − K2 + K2 I I II ⎜ 4 K II ⎝

θ c = 2 tan −1 ⎜

⎞ ⎟ ⎟ ⎠

(23)

In the present work, it is assumed that the interface toughness is relatively higher than the bulk materials, hence a crack may kink into any of the bulk material depending on their material properties. The computed ΔK Ieq is compared with the local fracture toughness of both the materials to decide in-which material the crack growth will occur. For this, two ratios Z1 and Z 2 are calculated as (Bhattacharya et al., 2013), Z1 =

( ΔK Ieq )m1 ( ΔK Ieq )m2 and Z 2 = ( K IC )m1 ( K IC )m2

(24)

where, m1 and m2 represents the material-1 and material-2 of the bi-material. If the Z1 > Z 2 , the crack will grow in the first material along angle θ = θc otherwise it will propagate into the second material parallel to the interface. For quasi-static crack propagation, the fatigue life is obtained using the generalized Paris law (Paris et al., 1961), which is defined as, m da = C ( ΔK Ieq ) dN

(25)

where, a is the crack length, N is the number of loading cycles, C and m are material constant for Paris law. ΔK Ieq can be obtained with ΔK I and ΔK II from Eq. (22). The failure takes place whenever ΔK Ieq > K IC for any material, and the number of cycles elapsed till failure gives the fatigue life of the material. 6. NUMERICAL RESULTS AND DISCUSSION

In this section, eight crack growth problems have been solved by the proposed virtual node XFEM and the results are compared with the literature and/or standard XFEM. First, a plate with an edge crack is considered under mode-I cyclic loading. The plate exhibits a strong discontinuity in the displacement field along the crack face. For this problem, the results obtained by VNXFEM are compared with the experimental and standard XFEM. In the second problem, a plate with an edge crack and circular hole is considered for the simulation. Third problem considers a plate with an edge crack and circular inclusion. The fourth problem is the plate with a bi-material interfacial edge crack. This problem possesses discontinuous displacement along the crack face and discontinuous derivative along the bi13

material interface. In fifth problem, a plate with an interfacial edge crack and circular hole is considered for the simulation. In the sixth problem, a plate with an interfacial edge crack and circular inclusion is taken for the simulation. In seventh problem, a dog-bone specimen along with an inclined edge crack is taken for the simulation. Finally, a miniature dumb-bell specimen with an inclined edge crack is taken for the simulation. These simulations are performed under plane stress condition. 6.1 Plate with an Edge Crack

A rectangular plate of size 90 mm × 108 mm with an initial edge crack length a = 45 mm is taken for the simulation as shown in Figure 13. The plate is subjected to a cyclic tensile load of Fmin = 8 kN and Fmax = 16 kN. The thickness of the plate is taken as 6 mm while the material properties (Ma et al., 2006) of the plate are taken as, Young modulus, E = 200 GPa Poisson's ratio, ν = 0.30 Paris exponent, m = 2.1 Paris constant, C = 7 × 10−8 In the proposed VNXFEM, transition elements are used to obtain the consistency in displacement field at the junction nodes. The displacement field is consistent at the nodes (virtual and standard) while it is inconsistent at the other points on the edge as shown in Figure 8(c). To check the effect of edge inconsistency, a convergence study is performed for

SIF. The variation of SIF with number of nodes is shown in Figure 14 for VNXFEM and standard XFEM. From the results presented in Figure 14, it is observed that the convergence rate is nearly same for both VNXFEM and standard XFEM. In the proposed VNXFEM, the problem has been solved for two cases; in first case, a uniform mesh of size 30 × 40 nodes is used with crack growth increments of 2.0 mm and 4.0 mm, and in second case, a uniform crack growth increment of 2.0 mm is used with two different mesh sizes of 30 × 40 nodes and 60 × 80 nodes. In standard XFEM, the same problem has been solved for two uniform mesh sizes i.e. 30 × 40 nodes (element size: 3.10 mm × 2.77 mm) with a crack growth increment of 4.0 mm and 60 × 80 nodes (element size: 1.52 mm × 1.37 mm) with a crack growth increment of 2.0 mm. Standard XFEM requires that the crack tip must fall in the next element after each crack growth increment. Therefore, two different crack growth increments i.e. 4.0 mm (>3.10 mm) for a mesh size of 30 × 40 nodes and 2.0 mm (>1.52 mm) for a mesh size of 60 × 80 nodes, are taken for the simulations. A comparison of the results obtained by the VNXFEM with standard XFEM and 14

experimental (Ma et al., 2006) is presented in Table 1(a) for a mesh of 30 × 40 nodes using different crack growth increments, and in Table 1(b) for a crack growth increment of 2.0 mm using different mesh sizes. A contour plot of σ yy obtained by VNXFEM is shown in Figure 15 for the final crack length. As expected in Figure 15, the maximum tensile stress occurs

near the crack tip and a compressive stress zone is developed near the far end of the unbroken ligament due to the nature of specimen loading and boundary conditions. Figure 16(a) presents a variation of SIF with the crack length for a uniform mesh size of 30 × 40 nodes. In Figure 16(a), standard XFEM results are obtained for a crack growth increment of 4.0 mm,

while the VNXFEM results are obtained for two different crack growth increments of 2.0 mm and 4.0 mm. Figure 16(b) presents a variation of SIF with the crack length for a crack growth increment of 2.0 mm. In Figure 16(b), standard XFEM results are obtained for a mesh size of 60 × 80 nodes, while the VNXFEM results are obtained for two different mesh sizes of 30 × 40 nodes and 60 × 80 nodes. The values of da/dN obtained by the VNXFEM are compared with the experimental and standard XFEM as shown in Figure 17. In Figure 17(a), VNXFEM and standard XFEM results are obtained for a mesh size of 30 × 40 nodes with different crack growth increments whereas in Figure 17(b), VNXFEM and standard XFEM results are obtained using a crack growth increment of 2.0 mm for different mesh sizes. The computation time required for a mesh size of 30 × 40 nodes is found to be 61.5 seconds and 59.6 seconds in case of VNXFEM and standard XFEM respectively. These simulations show

that the results obtained by VNXFEM are in good agreement with those obtained by standard XFEM and experimental results. Moreover, the results obtained by proposed VNXFEM with course mesh are quite close to the results of standard XFEM with fine mesh without much affecting the computational efficiency. 6.2 Plate with an Edge Crack and Circular Hole

A rectangular plate of size 50 mm x 100 mm with an initial edge crack of length a = 10 mm and a circular hole of radius r = 10 mm is taken for the simulation as shown in Figure 18. The center of hole is placed at x = 20 mm and y = 70 mm. The plate is subjected to a tensile load of σ = 80 MPa at the top edge while the bottom edge of the plate is constrained in ydirection. The following material properties of the plate are taken from the Ref. (Singh et al., 2012): Young modulus, E = 74 GPa Poisson's ratio, ν = 0.30 Fracture toughness, K IC = 1897.36 N/mm3/2 15

This problem has been solved by proposed VNXFEM and standard XFEM using coarse mesh size of 15 × 30 nodes (element size: 3.57 mm × 3.45 mm). In VNXFEM, the crack growth within an element can be modeled easily, while in standard XFEM, the modeling of crack growth within the element is not possible. Hence, to model the crack growth in standard XFEM, the crack increment should be taken in such a way that the new crack tip always falls in the next element. Therefore, in standard XFEM, the crack growth increment is taken as 4.0 mm (>element size: 3.57 mm × 3.45 mm) while in the proposed VNXFEM, the crack growth increments are taken as 0.50 mm and 4.0 mm. A contour plot of σ yy obtained using VNXFEM is shown in Figure 19. The crack growth trajectories obtained by VNXFEM and standard XFEM are compared in Figure 20. The values of stress intensity factor obtained by proposed VNXFEM are compared with those obtained by standard XFEM as illustrated in Figure 21. For this problem, Figure 22 shows the points taken to obtain the stresses along

and across the kinked crack in the tip element. The stress variations along and across the kinked crack at these points are shown in Figure 23(a) and Figure 23(b) respectively. In Figure 23(a), the normal stress decreases exponentially ahead of the crack tip.

In the previous example of Section 6.1, a comparison between proposed VNXFEM and standard XFEM could not be made properly due to the requirement of larger crack growth increment in standard XFEM. To further check the accuracy and effectiveness of the proposed approach, several crack growth problems are solved for two cases given in Table 2. In the first case, the mesh size is same while the crack growth increment is different. In the second case, the crack growth increment is same while the mesh size is different. For curved crack growth modeling, the standard XFEM requires that the crack tip must fall in the next element (Singh et al., 2012; Pathak et al., 2013) after each crack growth increment. Therefore, the crack growth increments are taken more than the size of the element i.e. 2.0 mm (>1.72 mm) and 1.0 mm (>0.85 mm) for the mesh sizes of 30 × 60 nodes and 60 × 120 nodes respectively. Further, to check the effectiveness of the proposed approach, an edge crack plate with a circular hole as explained above is also solved for two different cases as given in Table 2. For the first case (different crack growth increments), the crack growth trajectories obtained by proposed VNXFEM are compared with the crack growth trajectory obtained by standard XFEM as shown in Figure 24(a). In the second case (different mesh sizes), the crack growth trajectories obtained by VNXFEM are compared with the crack growth trajectory obtained by standard XFEM as shown in Figure 24(b). The values of stress intensity factor obtained by 16

VNXFEM are compared with those obtained by standard XFEM as illustrated in Figure 25(a) and Figure 25(b) for the first and second cases respectively. From these simulations, it

has been observed that the results obtained by VNXFEM are in good agreement with standard XFEM even for small crack growth increments. In the presence of hole, the crack deviates towards the hole even under mode-I loading. 6.3 Plate with an Edge Crack and Circular Inclusion

In this problem, a rectangular plate of size 50 mm × 100 mm with an initial edge crack of length a = 10 mm and a circular inclusion of radius r = 10 mm is taken for the simulation as shown in Figure 26. The inclusion centre is taken at x = 18 mm and y = 65 mm. The loading, boundary conditions and material properties are taken same as in Section 6.2. The material properties for the inclusion are taken as, Young's modulus, E = 10 GPa, and Poisson's ratio,

ν = 0.33. This problem is also solved for two cases as given in Table 2. The contour plot of σ yy obtained by VNXFEM for a mesh size of 30 × 60 nodes is shown in Figure 27. For the first case, the crack growth paths predicted by VNXFEM using crack growth increments of 1.0 mm and 2.0 mm are compared with the crack growth path predicted by standard XFEM using crack growth increments of 2.0 mm as shown in Figure 28(a). In the second case, the crack growth paths predicted by VNXFEM using mesh sizes of 30 × 60 nodes and 60 × 120 nodes are compared with the crack growth path predicted by standard XFEM using the mesh size of 60 × 120 nodes as shown in Figure 28(b). From the simulations, it has been observed that the crack deviates towards the soft inclusion. The values of SIF evaluated by VNXFEM are compared with standard XFEM as depicted in Figure 29(a) and Figure 29(b) for the first and second case respectively. From these simulations, it is also noticed that the results obtained by the VNXFEM using coarse mesh are in good agreement with those obtained by standard XFEM using fine mesh. 6.4 Plate with an Interfacial Edge Crack

A bi-material rectangular plate of size 50 mm x 100 mm with an interfacial edge crack of length a = 10 mm is taken for the simulation as shown in Figure 30. A tensile load of σ = 50 MPa has been applied at the top edge of the plate, while the bottom edge of the plate is constrained in y-direction. The thickness of the plate is taken as 1.0 mm. The material properties for both the materials are given in Table 3.

17

The contour plot of σ yy obtained by proposed VNXFEM for a mesh of 30 × 60 nodes is shown in Figure 31. For the first case data as per Table 2, the crack path trajectories of a bimaterial interfacial edge crack obtained by VNXFEM for crack growth increments of 0.5 mm and 1.0 mm are compared with the crack growth trajectory obtained by standard XFEM as illustrated in Figure 32(a). In the second case as per Table 2, the crack path trajectories obtained by VNXFEM for mesh sizes of 20 × 40 nodes and 30 × 60 nodes are compared with the standard XFEM as given in Figure 32(b). The crack path trajectories show that the crack kinks out from the interface towards the softer material. From the simulations, it has been observed that the crack paths obtained by proposed VNXFEM and standard XFEM are quite close to each other for different meshes and crack growth increments. The SIF values obtained by VNXFEM and standard XFEM for two cases of Table 2 are shown in Figure 33(a) and Figure 33(b) respectively. 6.5 Plate with an Interfacial Edge Crack and Circular Hole

In this problem, a bi-material rectangular plate of size 50 mm x 100 mm with an interfacial edge crack of length a = 10 mm and a circular hole of radius r = 10 mm is taken for the simulation as shown in Figure 34. The co-ordinates of the hole centre are taken as x = 20 mm and y = 70 mm. The loading, boundary conditions and material properties are kept same as explained in Section 6.5. The stress contour plot of σ yy obtained by VNXFEM for a mesh size of 30 × 60 nodes is shown in Figure 35. The crack growth paths obtained by VNXFEM and standard XFEM for the first and second cases of Table 2 are shown in Figure 36(a) and Figure 36(b) respectively. The SIF values obtained by VNXFEM and standard XFEM are shown in Figure 37(a) for the first case and in Figure 37(b) for second case. From the simulations, it has been

noticed that the results obtained by VNXFEM using coarse mesh are quite close to the results obtained by standard XFEM using fine mesh. 6.6 Plate with an Interfacial Edge Crack and Circular Inclusion

In this case, a bi-material rectangular plate of size 50 mm x 100 mm with an interfacial edge crack of length a = 10 mm and a circular inclusion of radius r = 10 mm is taken for the simulation as shown in Figure 38. The co-ordinates of the centre of inclusion are taken as x = 20 mm and y = 70 mm. The loading, boundary conditions and material properties are kept same as given in Section 6.5. The material properties for the inclusion are taken as, Young's modulus, E = 20 GPa and Poisson's ratio, ν = 0.33.

18

Figure 39 shows the contour plot of σ yy obtained by VNXFEM for a mesh size of

30 × 60 nodes. For first case (Table 2), the crack growth paths obtained by VNXFEM for crack growth increments of 0.5 mm and 1.0 mm are compared with the standard XFEM in Figure 40(a) whereas the crack growth paths obtained by VNXFEM for the mesh sizes of

20 × 40 nodes and 30 × 60 nodes are compared with the standard XFEM in Figure 40(b) for second case (Table 2). The SIF values obtained by VNXFEM and standard XFEM are plotted in Figure 41(a) and Figure 41(b) for the first and second cases respectively. From these simulations, it has been observed that the results obtained by VNXFEM and standard XFEM are found quite close to each other. 6.7 Dog-bone Specimen with an Inclined Edge Crack

In this case, a dog-bone specimen along with an inclined edge crack at an angle of 40o from horizontal is taken for the simulation as shown in Figure 42. The dimensions of the specimen are taken as; L = 50 mm, H = 100 mm, D = 55 mm and a = 7 mm. The loading, boundary conditions and material properties are taken same as Section 6.2. Figure 43 shows the contour plot of σ yy obtained by VNXFEM for a mesh size of

30 × 60 nodes. In this case, the elements are skewed in the central region of the specimen due to the geometry of the specimen as illustrated in Figure 43. The virtual nodes divide the actual tip element into VSE and VTE, and the shape of VSE and VTE is also skewed. For the first case of Table 2, the crack growth paths obtained by VNXFEM for the crack growth increments of 1.0 mm and 2.0 mm are compared with the standard XFEM as shown in Figure 44(a), whereas for the second case of Table 2, the crack growth paths obtained by

VNXFEM for mesh sizes of 30 × 60 nodes and 60 × 120 nodes are compared with the standard XFEM as shown in Figure 44(b). The SIF values obtained by VNXFEM and standard XFEM are plotted in Figure 45(a) and Figure 45(b) for the first and second case respectively. From these simulations, it has been observed that the results obtained by VNXFEM and standard XFEM are found quite close to each other. 6.8 Miniature Dumb-bell Specimen with an Inclined Edge Crack

Finally, a miniature dumb-bell specimen along with an edge crack of length a = 10 mm (making an angle of 30o from the horizontal) is taken for the simulation as shown in Figure 46. The dimensions, loading and boundary conditions of the specimen are also illustrated in Figure 46. The specimen is subjected to a tensile load of σ = 80 MPa. The material

properties for this problem are taken same as in Section 6.2. 19

In the proposed VNXFEM, the problem has been solved for two cases; in first case, a mesh size consist of 1188 elements is taken for the crack growth increments of 0.5 mm, 1.0 mm and 2.0 mm, while in second case, two different mesh sizes consisting of 1188 elements and 2794 elements are used for a uniform crack growth increment of 1.0 mm. In standard XFEM, the same problem has been solved for two mesh sizes i.e. 1188 elements with a crack growth increment of 2.0 mm and 2794 elements with a crack growth increment of 1.0 mm. Figure 47 shows the normal stress contour plot obtained by VNXFEM for a mesh size of

2794 elements. In this case, the elements are skewed in the central region of the specimen due to the geometry of the specimen. The virtual nodes divide the actual tip element into skewed VSE and skewed VTE. For the first case, the crack growth paths obtained by VNXFEM for the crack growth increments of 1.0 mm and 2.0 mm are compared with the standard XFEM as shown in Figure 48(a), whereas for the second case, the crack growth paths obtained by VNXFEM for mesh sizes of 1188 elements and 2794 elements are compared with the standard XFEM as shown in Figure 48(b). The SIF values obtained by VNXFEM and standard XFEM are plotted in Figure 49(a) and Figure 49(b) for the first and second case respectively. From these simulations, it is observed that the results obtained by VNXFEM and standard XFEM are quite close to each other. From the above simulations, it has been observed that the proposed VNXFEM can simulate crack growth in a single element using very small crack growth increments which is not possible in case of standard XFEM. The main advantages of the VNXFEM over standard XFEM are as follows, firstly it can model crack growth within the element using very small crack growth increments which is limited by element size in case of standard XFEM. Secondly, the VNXFEM provides quite accurate results for highly curved cracks even at coarse mesh whereas standard XFEM requires very fine mesh for the modeling of crack growth using small crack increments as standard XFEM requires that the crack tip must fall in the new element after each crack growth increment for curved crack growth modeling. 7. CONCLUSIONS

In the present work, a new virtual node XFEM (VNXFEM) has been proposed for the simulation of crack growth within an element. To model the kinks in the element, the concept of polygon law has been proposed to decide the position of the evaluation (Gauss) points with respect to the crack faces. The proposed VNXFEM provides a smooth crack growth within an element by dividing an actual crack tip element into virtual split element and virtual tip element with the help of virtual nodes. The VNXFEM leads to higher accuracy due to local 20

mesh refinement near the crack tip through virtual nodes. It models the crack growth accurately with coarse mesh using small crack growth increments whereas standard XFEM requires the fine mesh to model the same as crack tip must fall into next element during crack growth. Although, only eight kinks have been modeled in an element in the present work however, the number of kinks within an element is not limited to eight, and can be extended to any number as per the requirement. The results obtained by the VNXFEM using coarse mesh are found in good agreement with those obtained by standard XFEM using fine mesh. The VNXFEM not only simulates the crack path more accurately but also saves the computational time for the same accuracy. The numerical results also demonstrate that the combination of virtual nodes with the XFEM can handle a large class of problems with good accuracy. Hence, the concept of virtual node will be extended further for the simulation of crack growth in complex geometries. Acknowledgement: This work is not supported by any specific funding from the funding

agency. Appendix

The elemental stiffness matrices K and force vector f in Eq. (5) can be written as (Mohammadi, 2008), ⎡K uu ij ⎢ au ⎢K ij ⎢ bu ⎢K ij K ij = ⎢ c u ⎢K ij ⎢K d u ⎢ ij ⎢K e u ⎣ ij

{

f i = f iu

⎤ K ue ij ⎥ K ijae ⎥ ⎥ K be ij ⎥ ⎥ K ijce ⎥ K ijde ⎥ ⎥ K ijee ⎥⎦

(A1)

}

(A2)

where, r , s = u , a, b, c, d , e

(A3)

K ua ij

K ub ij

K uc ij

K ud ij

K ijaa

K ijab

K ijac

K ijad

K ba ij

K bb ij

K bc ij

K bd ij

K ijca

K ijcb

K ijcc

K ijcd

K ijda

K ijdb

K ijdc

K ijdd

K ijea

K ijeb

K ijec

K ijed

fia

fib

fic

K rsij = ∫ (B ir )T D B sj dΩ

fid

fie

T

Ωe

fiu =

∫ Ν b dΩ + ∫ Ν t d Γ i

Γt

Ω

fia =

(A4)

i

e

∫ Ν ( Η (x ) − H (x )) b d Ω + ∫ Ν ( Η (x ) − H (x )) td Γ

(A5)

{

(A6)

i

i

i

Ωe

fib = fib1

i

Γt

fib 2

fib 3

fib 4

}

T

21

fibα =

∫ Ν βα ((x) − (x )) b d Ω + ∫ Ν ( βα (x ) − (x )) t d Γ , i

i

i

Ωe

i

where, α = 1, 2,3, 4

Γt

(A7) fic =

∫ Ν (V(x ) − V (x )) b d Ω + ∫ Ν i

Ω

fid =

i

e

(V ( x ) − V ( x i )) td Γ

∫ Ν ϕ ( x ) b d Ω + ∫ Ν ϕ ( x ) td Γ i

∫Ν

Ω

e

(A8)

Γt

(A9)

i

Ωe

fie =

i

Γt

i

χ ( x ) b d Ω + ∫ Ν i χ ( x ) td Γ

(A10)

Γt

where, Ν i are finite element shape functions. Βiu , Βia , Βib , Βic , Βid and Β ie are the matrices of shape function derivatives, which are given below,

⎡ Νi,x 0 ⎤ ⎢ ⎥ Β = ⎢ 0 Νi, y ⎥ ⎢⎣ Ν i , y Ν i , x ⎥⎦

(A11a)

⎤ ⎡ ( Ν i ( Η ( x) − H (xi )) ), x 0 ⎥ ⎢ 0 Β ia = ⎢ ( Ν i ( Η (x) − H (xi )) ), y ⎥ ⎥ ⎢( Ν ( Η ( x) − H ( x )) ) i , y ( Ν i ( Η ( x ) − H ( x i )) ) ⎥ ⎣⎢ i ,x ⎦

(A11b)

Βib = ⎡⎣ Βib1 Βib 2

(A11c)

u i

Β ibα

Βib 4 ⎤⎦

Βib 3

⎡ ( Ν i ( βα ( x ) − βα ( x i )) ), x ⎢ 0 =⎢ ⎢ ( Ν ( β ( x ) − β ( x )) ) i α ,y ⎣ i α

⎤ ⎥ ( Ν i ( βα (x ) − βα (x i )) ), y ⎥ where, α = 1, 2,3, 4 ⎥ ( Ν i ( βα (x ) − βα (xi )) ), x ⎦⎥ 0

(A11d) ⎤ ⎡ ( Ν i (V ( x ) − V ( x i )) ), x 0 ⎥ ⎢ 0 Βic = ⎢ ( Ν i (V (x ) − V (xi )) ), y ⎥ ⎥ ⎢( Ν (V ( x ) − V ( x )) ) i i Ν − x x ( ( ) ( )) V V ) ( y , ⎥ i i ⎣ ,x ⎦

(A11e)

⎡ ( Ν iϕ ( x ) ), x ⎢ d 0 Βi = ⎢ ⎢( Ν ϕ (x ) ) ,y ⎣ i

(A11f)

⎤ ⎥ ( Ν iϕ (x ) ), y ⎥ ⎥ ( Ν iϕ (x ) ), x ⎥⎦ 0

22

⎡ ( Ν i χ ( x ) ), x ⎢ 0 Βie = ⎢ ⎢( Ν χ (x) ) ,y ⎣ i

⎤ ⎥ ( Ν i χ ( x ) ), y ⎥ ⎥ ( Ν i χ (x ) ), x ⎥⎦ 0

(A11g)

REFRENCES Anahid, M. and Khoei, A.R. (2008), New development in extended finite element modeling of large elasto-plastic deformations, International Journal for Numerical Methods in Engineering, vol. 75, pp. 1133–1171. Asadpoure, A., Mohammadi, S. and Vafai, A. (2006), Crack analysis in orthotropic media using the extended finite element method, Thin-Walled Structures, vol. 44, pp. 1031– 1038. Asadpoure, A. and Mohammadi, S. (2007), Developing new enrichment functions for crack simulation in orthotropic media by the extended finite element method, International Journal for Numerical Methods in Engineering, vol. 69, pp. 2150–2172. Ashari, S.E. and Mohammadi S. (2011), Delamination analysis of composites by new orthotropic bimaterial extended finite element method, International Journal for Numerical Methods in Engineering, vol. 86, pp. 1507–1543. Bac, N.V., Xuan, H.N., Chen, L., Lee, C.K., Zi, G., Zhuang, X., Liu, G.R. and Rabczuk, T. (2013), A phantom-node method with edge-based strain smoothing for linear elastic fracture mechanics, Journal of Applied Mathematics, Vol. 2013, Article ID 978026, 12 pages. Bayesteh, H. and Mohammadi, S. (2013), XFEM fracture analysis of orthotropic functionally graded materials, Composites: Part B, vol. 44, pp. 8–25. Belytschko, T., Lu, Y.Y. and Gu, L. (1994), Element-free Galerkin methods, International Journal for Numerical Methods in Engineering, vol. 37, pp. 229–256. Belytschko, T. and Black, T. (1999), Elastic crack growth in finite elements with minimal remeshing, International Journal for Numerical Methods in Engineering, vol. 45, pp. 601–620. Bhattacharya, S., Singh, I.V., Mishra, B.K. and Bui, T.Q. (2013), Fatigue crack growth simulations of interfacial cracks in bi-layered FGMs using XFEM, Computational Mechanics, vol. 52, pp. 799–814. Bouhala, L., Shao, Q., Koutsawa, Y., Younes, A., Nunez, P., Makradi, A. and Belouettar, S. (2013), An XFEM crack-tip enrichment for a crack terminating at a bi-material interface, Engineering Fracture Mechanics, vol. 102, pp. 51–64. Chau-Dinh., T., Zi, G., Lee, P.S., Rabczuk, T. and Song, J.H. (2012), Phantom-node method for shell models with arbitrary cracks, Computers and Structures, vol. 92-93, pp. 242–256. Chen, L., Rabczuk, T., Bordas, S.P.A., Liu, G.R., Zeng, K.Y. and Kerfriden, P. (2012), Extended finite element method with edge-based strain smoothing (ESm-XFEM) for linear elastic crack growth, Computer Methods in Applied Mechanics and Engineering, vol. 209-212, pp. 250–265.

23

Cheng, K.W. and Fries, T.P. (2012), XFEM with hanging nodes for two-phase incompressible flow, Computer methods in applied mechanics and engineering, vol. 245246, pp. 290–312. Cho, Y.S. and Im, S. (2006), MLS-based variable-node elements compatible with quadratic interpolation. Part I: Formulation and application for non-matching meshes, International Journal for Numerical Methods in Engineering, vol. 65, pp. 494–516. Daux, C., Moes, N., Dolbow, J., Sukumar, N. and Belytschko, T. (2000), Arbitrary branched and intersected cracks with the extended finite element method, International Journal for Numerical Methods in Engineering, vol. 48, pp. 1741–1760. Eftekhari, M., Ardakani, S.H. and Mohammadi, S. (2014), An XFEM multiscale approach for fracture analysis of carbon nanotube reinforced concrete, Theoretical and Applied Fracture Mechanics, vol. 72, pp. 64–75. Elguedj, T., Gravouil, A. and Combescure, A. (2006), Appropriate extended functions for XFEM simulation of plastic fracture mechanics, Computer Methods in Applied Mechanics and Engineering, vol. 195, pp. 501–515. Erdogan, F. and Sih, G. (1963), On the crack extension in plates under plane loading and transverse shear, Journal of Basic Engineering, vol. 85, pp. 519–527. Fries, T.P. (2008), A corrected XFEM approximation without problems in blending elements, International Journal for Numerical Methods in Engineering, vol. 75, pp. 503–532. Fries, T.P., Byfut, A., Alizada, A., Cheng, K.W. and Schroder, A. (2011), Hanging nodes and XFEM, International Journal for Numerical Methods in Engineering, vol. 86, pp. 404– 430. Gracie, R., Wang, H. and Belytschko, T. (2008), Blending in the extended finite element method by discontinuous Galerkin and assumed strain methods, International Journal for Numerical Methods in Engineering, vol. 74, pp. 1645–1669. Gupta, A.K. (1978), A finite element for transition from a fine grid to a coarse grid, International Journal for Numerical Methods in Engineering, vol. 12, pp. 35–45. Huang, F. and Xie, X. (2010), A modified nonconforming 5-node quadrilateral transition finite element, Advances in Applied Mathematics and Mechanics, vol. 2, pp. 784–797. Johnson, J. and Qu, J. (2006), An interaction integral method for computing mixed mode stress intensity factors for curved bimaterial interface cracks in non-uniform temperature fields, Engineering Fracture Mechanics, vol. 74, pp. 2282–2291. Khoei, A.R., Biabanaki, S.O.R. and Anahid, M. (2009), A Lagrangian-extended finiteelement method in modeling large-plasticity deformations and contact problems, International Journal of Mechanical Sciences, vol. 51, pp. 384–401. Kumar, S., Singh, I.V. and Mishra, B.K. (2014a), XFEM simulation of stable crack growth using J-R curve under finite strain plasticity, International Journal of Mechanics and Material in Design, vol. 10, pp. 165–177. Kumar, S., Singh, I.V. and Mishra, B.K. (2014b), A multigrid coupled (FE-EFG) approach to simulate fatigue crack growth in heterogeneous materials, Theoretical and Applied Fracture Mechanics, Vol. 72, pp. 121–135.

24

Loehnert, S., Prange, C. and Wriggers, P. (2012), Error controlled adaptive multiscale XFEM simulation of cracks, International Journal of Fracture, vol. 178, pp. 147–156. Ma, S., Zhang, X.B., Recho, N. and Li, J. (2006), The mixed-mode investigation of the fatigue crack in CTS metallic specimen, International Journal of Fatigue, vol. 28, pp. 1780–1790. Melenk, J.M. and Babuska, I. (1996), The partition of unity finite element method: Basic theory and applications, Computer Methods in Applied Mechanics and Engineering, vol. 139, pp. 289–314. Moes, N., Dolbow, J. and Belytschko, T. (1999), A finite element method for crack growth without remeshing, International Journal for Numerical Methods in Engineering, vol. 46, pp. 131–150. Mohammadi, S. (2008), Extended Finite Element Method: for Fracture Analysis of Structures, Blackwell Publishing Ltd. Motamedi, D. and Mohammadi, S. (2010a), Dynamic analysis of fixed cracks in composites by the extended finite element method, Engineering Fracture Mechanics, vol. 77, pp. 3373–3393. Motamedi, D. and Mohammadi, S. (2010b), Dynamic crack propagation analysis of orthotropic media by the extended finite element method. International Journal of Fracture, vol. 161, pp. 21–39. Motamedi, D. and Mohammadi, S. (2012), Fracture analysis of composites by time independent moving-crack orthotropic XFEM, International Journal of Mechanical Sciences, vol. 54, pp. 20–37. Natarajan, S., Bordas, S. and Mahapatra, D.R. (2009), Numerical integration over arbitrary polygonal domains based on Schwarz-Christoffel conformal mapping, International Journal for Numerical Methods in Engineering, vol. 80, pp. 103–134. Newman, J. (1971), An improved method of collocation for the stress analysis of cracked plates with various shaped boundaries. Technical Report, TN D–6376, NASA. Pant, M., Singh, I.V. and Mishra, B.K. (2011), Evaluation of mixed mode stress intensity factors for interface cracks using EFGM, Applied Mathematical Modeling, vol. 35, pp. 3443–3459. Paris, P.C., Gomez, M.P. and Anderson, W.E. (1961), A rational analytic theory of fatigue, The Trend in Engineering, vol. 13, pp. 9–14. Pathak, H., Singh, A. and Singh, I.V. (2013), Fatigue crack growth simulations of bi-material interfacial cracks under thermo-elastic loading by extended finite element method, European Journal of Computational Mechanics, vol. 22, pp. 79–104. Portela, A.A., Aliabadi, M. and Rooke, D. (1991), The dual boundary element method: effective implementation for crack problems, International Journal for Numerical Methods in Engineering, vol. 33, pp. 269–1287. Prange, C., Loehnert, S. and Wriggers, P. (2012), Error estimation for crack simulations using the XFEM, International Journal for Numerical Methods in Engineering, vol. 91, 1459–1474.

25

Rabczuk, T., Zi, G., Gerstenberger, A. and Wall, W.A. (2008), A new crack tip element for the phantom-node method with arbitrary cohesive cracks, International Journal for Numerical Methods in Engineering, vol. 75, pp. 577–599. Shibanuma, K. and Utsunomiya, T. (2009), Reformulation of XFEM based on PUFEM for solving problem caused by blending elements, Finite Elements in Analysis and Design, vol. 45, pp. 806–816. Singh, I.V., Mishra, B.K., Bhattacharya, S. and Patil, R.U. (2012), The numerical simulation of fatigue crack growth using extended finite element method, International Journal of Fatigue, vol. 36, pp. 109–119. Sukumar, N., Chopp, D.L., Moes, N. and Belytschko, T. (2001), Modeling Holes and inclusions by level sets in the extended finite element method, Computer Methods in Applied Mechanics and Engineering, vol. 190, pp. 6183–6200. Sukumar, N., Huang, Z.Y., Prevost, J.H. and Suo, Z. (2004), Partition of unity enrichment for bimaterial interface cracks, International Journal for Numerical Methods in Engineering, vol. 59, pp. 1075–1102. Sukumar, N. and Prevost, J.H. (2003), Modeling quasi-static crack growth with the extended finite element method, Part I: computer implementation, International Journal of Solids and Structures, vol. 40, pp. 7513–7537. Tarancon, J.E., Vercher, A., Giner, E. and Fuenmayor, F.J. (2009), Enhanced blending elements for XFEM applied to linear elastic fracture mechanics, International Journal for Numerical Methods in Engineering, vol. 77, pp. 126–148. Zi, G. and Belytschko, T. (2003), New crack-tip elements for XFEM and applications to cohesive cracks, International Journal for Numerical Methods in Engineering, vol. 57, pp. 2221–2240. Zienkiewicz, O.C. and Zhu, J.Z. (1987), A simple error estimator and adaptive procedure for practical engineering analysis, International Journal for Numerical Methods in Engineering, vol. 24, pp. 337–357.

26

Figures

t

y 

c u

x

z Figure 1: Domain with a crack and boundary conditions

Enriched nodes

Hole

Inclusion

Blending elements

Enriched nodes Tip element

Tip nodes Split elements

Split nodes

Figure 2: A typical discretization of the 2-D domain with crack, hole and inclusion

27

(a) Decomposition of a split element with a straight crack into sub-triangles

(b) Simplest decomposition process of a split element with kink into sub-triangles (introduces error due to approximation)

Gauss point

(c) Standard decomposition of split element with kink into sub-triangles

Figure 3: Triangulation process for the evaluation of Gauss points

28

Sub-triangles in VSE

VSE

Sub-triangles in VTE

VTE

=

+

(a)

Sub-triangles in VTE

Sub-triangles in VSE

VSE

VTE

=

+

(b)

Sub-triangles in VTE

Sub-triangles in VSE

VSE VTE

=

+

(c)

29

Sub-triangles in VSE

VSE

Sub-triangles in VTE

VTE

=

+

(d)

Sub-triangles in VSE

VSE

Sub-triangles in VTE

VTE

=

+

(e)

Sub-triangles in VSE

Sub-triangles in VTE

VSE VTE

=

+

Standard Node (f)

Virtual Node

Figure 4: Generation of VSE and VTE in standard tip element and sub-triangles in VSE and VTE

30

VTE

VSE VSE

VTE

Element-1

(a)

Element-2

Element-1

(b)

Element-2

Figure 5: Schematic representation of kink modeling in an element: (a) if kink is present in the actual tip element; (b) if kink is not present in the actual tip element

x 10

Displacement, vyy

Displacement, uxx

-3

20

0.03

15

-0.0032

0.02

-0.0048

0.01

-0.0048

0

0.0036

0.0022

10

-0.0022

5

-0.0032 0

-0.0036

-0.02

(a) -5

-0.01

(b)

-0.03

Figure 6: Variation of displacement over a virtual tip element: (a) displacement component in xdirection; (b) displacement component in y-direction yy

xx

xy

100

30

150

80

20 60

10

100 40

0

50

-10

20

-20

0

(a)

-20

0

(b)

(c)

Figure 7: Variation of stress fields over a virtual tip element: (a) normal stress, σ xx ; (b) normal stress,

σ yy ; (c) shear stress  xy

31

-30

4

3 Transition Element

VE

5

1

VE

2

VE

4

Quadratic variation

3

6

Linear variation

Transition Element

VE

1

(a)

2

(b)

(c)

Figure 8: Schematic representation of node numbering in transition elements (a) and (b); (c) displacement field variation along the virtual node edge in the transition element

1

0.9

N1 N2 N3 N4 N5 Sum(N)

0.8

0.6

0.4

1

0.8

0.8

0.7

0.6

0.6

0.4

0.5 0.4

0.2

0.2

0.3

0 1

0

0.2 0.5

1 0.5

0

-0.2

-0.5 -1

-1

-0.5

0

0.5

1

0.1

0

-0.5 -1

Figure 10: Shape function contour plot of the virtual node in the transition element; virtual node is at the mid point of the edge

Figure 9: Shape functions for the transition element; virtual node is at the mid point of the edge

O

P Figure 11: Numbers of intersection of a ray with polygon are odd if a point O lies inside the polygon, and numbers of intersection are even if a point P lies outside the polygon

32

4

3

1

2

2

3

4

1

3

1

2

(c)

(b)

(a) 4

5

5

6

3

4

6

7

3 1

5

2

4

3 1 (d)

1

2

2 (f)

(e)

Figure 12: Generation of polygon in different cases for the integration in VSE

108 mm

45 amm

54 mm

90 mm Figure 13: Edge crack plate along with dimensions

33

yy 16.4 VNXFEM Standrad XFEM

300

16.2

200 SIF (MPa.m1/2)

16

100 15.8

0 15.6

-100 15.4

-200 15.2

0

1000

2000 3000 Number of nodes

4000

5000

Figure 15: Stress distribution of  yy for an

Figure 14: SIF variation with number of nodes for VNXFEM and standard XFEM

edge crack plate

50 45

45 40

35

SIF (MPa.m1/2)

SIF (MPa.m1/2)

40

50

Experimental (Ma et al., 2006) VNXFEM (da = 2 mm) VNXFEM (da = 4 mm) Standard XFEM (da = 4 mm)

30 25

35 30 25

20

20

15

15

10 40

45

50

55 60 Crack Length (mm)

65

10 40

70

Experimental (Ma et al., 2006) VNXFEM (Mesh size: 30x40) VNXFEM (Mesh size: 60x80) Standard XFEM (Mesh size: 60x80)

45

50

55 60 Crack Length (mm)

65

70

Figure 16b: SIF variation with crack length of an edge crack for a crack growth increment (da) of 2 mm

Figure 16a: SIF variation with crack length of an edge crack for a mesh size of 30x40

34

2.5

x 10

-4

2.5

Experimental (Ma et al., 2006) VNXFEM (da = 2 mm) VNXFEM (da = 4 mm) Standard XFEM (da = 4 mm)

-4

Experimental (Ma et al., 2006) VNXFEM (da = 2 mm) VNXFEM (da = 4 mm) Standard XFEM (da = 4 mm)

2

da/dN (mm/cycle)

da/dN (mm/cycle)

2

x 10

1.5

1

1.5

1

0.5

0.5

0 15

20

25

30

35

40

45

0 15

50

20

25

30

35

40

45

50

SIF (MPa.m1/2)

SIF (MPa.m1/2)

Figure 17b: da/dN vs SIF plot of an edge crack plate for a crack growth increment (da) of 2 mm

Figure 17a: da/dN vs SIF plot of an edge crack plate for a mesh size of 30x40

 yy 500 400

Hole (x,y)

300

r

200

a H

100 0

H 2

-100

L

-200

Figure 19: Stress distribution of  yy for an edge

Figure 18: Plate with an edge crack, hole and dimensions

crack plate with hole

35

Crack Trajectories

Hole

Standard XFEM (da = 4.0 mm) VNXFEM (da = 0.5 mm) Figure 20: Crack growth trajectory of an edge crack in the presence of hole for a mesh size of 15x30

2000 VNXFEM (da = 0.5 mm) VNXFEM (da = 4 mm) Standard XFEM (da = 4 mm)

1800

SIF (MPa.mm1/2)

1600

1400

1200

1000

800

8

10

12

14

16 18 20 Crack Length (mm)

22

24

26

28

Figure 21: SIF variation with crack length of an edge crack plate in the presence of hole for a mesh size of 15x30

36

yy 300 12

250

11 10

200

9 8

Crack 1

7

3

2

6

8

6

4

10

13 12

150 100

Crack tip

5

Kink

7

5

11

9

50

4 3

0 2 1

-50 -100

Figure 22: Schematic representation of the points selected for the stress plot on the sections along and across the kinked crack

300

300

250 250

yy (MPa)

yy (MPa)

200

150

200

150

100 100

50

0

0

2

4

6 8 Number of points

10

12

50

14

(a)

0

2

4

6 8 Number of points

10

12

(b) Figure 23: Normal stress (  yy ) variation: (a) along the kinked crack section; (b) across the kinked crack section

37

14

Crack Trajectories

Hole

Standard XFEM (da = 2 mm) VNXFEM (da = 1 mm) VNXFEM (da = 0.5 mm)

Figure 24a: Crack growth trajectory of an edge crack in the presence of hole for a mesh size of 30x60

Crack Trajectories

Hole

Standard XFEM (Mesh size: 60 x 120) VNXFEM (Mesh size: 30 x 60) VNXFEM (Mesh size: 20 x 40) Figure 24b: Crack growth trajectory of an edge crack in the presence of hole for a crack growth increment (da) of 1 mm

38

2000

2000 VNXFEM (da = 0.5 mm) VNXFEM (da = 1 mm) VNXFEM (da = 2 mm) Standard XFEM (da = 2 mm)

1800

1800

1600

SIF (MPa.mm1/2)

SIF (MPa.mm1/2)

SIF (MPa.mm1/2)

1600

1400

1400

1200 1800

1200

1750 1000

1000

1700 800

VNXFEM (Mesh size: 20 x 40) VNXFEM (Mesh size: 30 x 60) VNXFEM (Mesh size: 60 x 120) Standard XFEM (Mesh size: 60 x 120)

8

10

12

14

16 18 20 Crack Length (mm)

22

24

800

26

1650 Figure

25a: SIF variation with crack length of an edge crack in the presence of hole for a mesh size of 30x60 VNXFEM (da = 1 mm) VNXFEM (da = 0.5 mm) VNXFEM (da = 2 mm) Standard XFEM (da = 2 mm)

1600

1550

18.5

19

19.5

20 20.5 21 Crack Length (mm)

21.5

22

8

10

12

14

16 18 20 Crack Length (mm)

22

24

Figure 25b: SIF variation with crack length of an edge crack in the presence of hole for a crack growth increment (da) of 1 mm

22.5

 yy 600

Inclusion

500

(x,y)

400

r

300

a H

200 100

H 2

0 -100

L -200

Figure 27: Stress distribution of  yy for an edge Figure 26: Plate with an edge crack, inclusion and dimensions

crack plate with inclusion

39

26

Crack Trajectories

Inclusion

Standard XFEM (da = 2 mm) VNXFEM (da = 1 mm) VNXFEM (da = 2 mm) Figure 28a: Crack growth trajectory of an edge crack in the presence of inclusion for a mesh size of 30x60

Crack Trajectories

Inclusion

Standard XFEM (Mesh size: 60 x 120) VNXFEM (Mesh size: 30 x 60) VNXFEM (Mesh size: 60 x 120) Figure 28b: Crack growth trajectory of an edge crack in the presence of inclusion for a crack growth increment of da = 1 mm

40

1700

1700 VNXFEM (da = 0.5 mm) VNXFEM (da = 1 mm) VNXFEM (da = 2 mm) Standard XFEM (da = 2 mm)

1600 1500

1500 1400 SIF (MPa.mm1/2)

SIF (MPa.mm1/2)

1400 1300 1200 1100

1300 1200 1100

1000

1000

900

900

800

VNXFEM (Mesh size: 20 x 40) VNXFEM (Mesh size: 30 x 60) VNXFEM (Mesh size: 60 x 120) Standard XFEM (Mesh size: 60 x 120)

1600

8

10

12

14

16 18 20 Crack Length (mm)

22

24

800

26

Figure 29a: SIF variation with crack length of an edge crack in the presence of inclusion for a mesh size of 30x60

8

10

12

14

16 18 20 Crack Length (mm)

22

24

Figure 29b: SIF variation with crack length of an edge crack in the presence of inclusion for a crack growth increment (da) of 1 mm



400

Material-1

300

E1, 1 a

200

H 100

H 2

Material-2

E2,  2

0

L

-100

Figure 31: Stress distribution,  yy for a bi-

Figure 30: Bi-material plate with an interfacial edge crack and dimensions

material interfacial edge crack plate

41

26

Crack Trajectories Material-1

Material Interface

Standard XFEM (da = 2 mm) Material-2

VNXFEM (da = 1 mm) VNXFEM (da = 0.5 mm)

Figure 32a: Crack growth trajectory of a bi-material interface edge crack for a mesh size of 30x60

Material-1

Material Interface

Crack Trajectories

Material-2

Standard XFEM (Mesh size: 60 x 120) VNXFEM (Mesh size: 30 x 60) VNXFEM (Mesh size: 20 x 40)

Figure 32b: Crack growth trajectory of a bi-material interface edge crack for a crack growth increment (da) of 1 mm

42

1400

1400 VNXFEM (da = 0.5 mm) VNXFEM (da = 1 mm) VNXFEM (da = 2 mm) Standard XFEM (da = 2 mm)

1200

1200

1000 SIF (MPa.mm1/2)

1000

800 940 920

600

800 1300

600

400

900

1250 1/2

SIF (MPa.mm )

1/2

SIF (MPa.mm )

SIF (MPa.mm1/2)

VNXFEM (Mesh size: 20 x 40) VNXFEM (Mesh size: 30 x 60) VNXFEM (Mesh size: 60 x 120) Standard XFEM (Mesh size: 60 x 120)

880

400 860 840

200

8

10

12

14

16 18 20 Crack Length (mm) 820

18.5

19

22

24

200

26

VNXFEM 19.5 20(da = 0.5 mm)20.5 VNXFEM (da = 1 mm) Crack Length (mm) VNXFEM (da = 2 mm) Standard XFEM (da = 2 mm)

Figure 33a: SIF variation with crack length of an interface edge crack plate for a mesh size of 30x60

8

10

12

1200

1150

14

16 18 20 Crack Length (mm)

22

24

26

1100

Figure 33b: SIF variation with crackVNXFEM length(Mesh of an size: 20 x 40) VNXFEM (Mesh size: 30 x 60) interface edge crack for a crack growth increment (da) ofx 120) VNXFEM (Mesh size: 60 1050 21 21.5 23.5XFEM 24(Mesh24.5 size: 6025 x 120) 1 mm22 22.5Crack23Standard Length (mm)

 250

Material-1 E1, 1

Hole

200

(x,y) r

150 100

a H

50

H 2

0

Material-2

E2,  2

-50

L

-100

Figure 35: Stress distribution of  yy for a bi-

Figure 34: Bi-material plate with an interfacial edge crack, hole and dimensions

material interfacial edge crack plate with hole

43

Crack Trajectories Material-1

Hole

Material Interface

Standard XFEM (da = 2 mm) VNXFEM (da = 1 mm)

Material-2

VNXFEM (da = 0.5 mm) Figure 36a: Crack growth trajectory of an interface edge crack in the presence of hole for a mesh size of 30x60

Crack Trajectories

Material-1

Hole

Material Interface

Standard XFEM (Mesh size: 60 x 120) VNXFEM (Mesh size: 30 x 60)

Material-2

VNXFEM (Mesh size: 20 x 40) Figure 36b: Crack growth trajectory of an interface edge crack in the presence of hole for a crack growth increment (da) of 1 mm

44

1400

1400 VNXFEM (da = 0.5 mm) VNXFEM (da = 1 mm) VNXFEM (da = 2 mm) Standard XFEM (da = 2 mm)

1300 1200

1200 1100 SIF (MPa.mm1/2)

SIF (MPa.mm1/2)

1100 1000 900

1000 900

800

800

700

700

600

600

500

VNXFEM (Mesh size: 20 x 40) VNXFEM (Mesh size: 30 x 60) VNXFEM (Mesh size: 60 x 120) Standard XFEM (Mesh size: 60 x 120)

1300

500

8

10

12

14 16 Crack Length (mm)

18

20

Figure 37a: SIF variation with crack length of a bimaterial interface edge crack plate in the presence of hole for a mesh size of 30x60

8

10

12

14 16 Crack Length (mm)

18

20

Figure 37b: SIF variation with crack length of a bimaterial interface edge crack plate in the presence of hole for a crack growth increment (da) of 1 mm



200

Material-1 E1, 1 Inclusion

150

(x,y) r

100

a H

50

H 2

Material-2 0

E2,  2 L

-50

Figure 39: Stress distribution of  yy for a bi-

Figure 38: Bi-material plate with an interfacial edge crack, inclusion and dimensions

material interfacial edge crack plate with inclusion

45

Crack Trajectories Material-1

Inclusion

Material Interface

Standard XFEM (da = 2 mm) VNXFEM (da = 1 mm)

Material-2

VNXFEM (da = 0.5 mm)

Figure 40a: Crack growth trajectory of an interface edge crack in the presence of inclusion for a mesh size of 30x60

Crack Trajectories

Material-1

Inclusion

Material Interface

Standard XFEM (Mesh size: 60 x 120) Material-2

VNXFEM (Mesh size: 30 x 60) VNXFEM (Mesh size: 20 x 40)

Figure 40b: Crack growth trajectory of an interface edge crack in the presence of inclusion for a crack growth increment (da) of 1 mm

46

1400

1400 VNXFEM (da = 0.5 mm) VNXFEM (da = 1 mm) VNXFEM (da = 2 mm) Standard XFEM (da = 2 mm)

1200

VNXFEM (Mesh size: 20 x 40) VNXFEM (Mesh size: 30 x 60) VNXFEM (Mesh size: 60 x 120) Standard XFEM (Mesh size: 60 x 120)

1300 1200 1100 SIF (MPa.mm1/2)

SIF (MPa.mm1/2)

1000

800

600

1000 900 800 700 600

400

500 200

400 8

10

12

14 16 18 Crack Length (mm)

20

22

24

8

10

12

14 16 18 Crack Length (mm)

20

22

24

Figure 41b: SIF variation with crack length of an interface edge crack in the presence of inclusion for a crack growth increment (da) of 1 mm

Figure 41a: SIF variation with crack length of an interface edge crack in the presence of inclusion for a mesh size of 30x60



a H

D

L

Figure 42: Dog-bone specimen geometry along with an inclined edge crack

47

yy

350 

yy VTE

350

300

300

Virtual Nodes

250

250

200

200

150 100

150 VSE

50 0

100 50 0

Figure 43: Stress distribution,  yy in a dog-bone specimen in the presence of inclined edge crack

Crack Trajectories

Standard XFEM (da = 2 mm) VNXFEM (da = 1 mm) VNXFEM (da = 2 mm) Figure 44a: Crack growth trajectory of an inclined edge crack in the dog-bone specimen for a mesh size of 30x60

48

Crack Trajectories

Standard XFEM (Mesh size: 60 x 120) VNXFEM (Mesh size: 30 x 60) VNXFEM (Mesh size: 60 x 120)

Figure 44b: Crack growth trajectory of an inclined edge crack in a dog-bone specimen for a crack growth increment of da = 1 mm

2500

2200 VNXFEM (da = 0.5 mm) VNXFEM (da = 1 mm) VNXFEM (da = 2 mm) Standard XFEM (da = 2 mm)

2000

1800 SIF (MPa.mm1/2)

SIF (MPa.mm1/2)

VNXFEM (Mesh size: 20 x 40) VNXFEM (Mesh size: 30 x 60) VNXFEM (Mesh size: 60 x 120) Standard XFEM (Mesh size: 60 x 120)

2000

1500

1000

1600 1400 1200 1000 800

500

6

8

10

12 14 16 Crack Length (mm)

18

20

600

22

Figure 45a: SIF variation with crack length of an inclined edge crack in a dog-bone specimen for a mesh size of 30x60

6

8

10

12 14 16 Crack Length (mm)

18

20

Figure 45b: SIF variation with crack length of an inclined edge in a dog-bone specimen for a crack growth increment (da) of 1 mm

49

22

5 10 40 6 15

a

11.5

100

All dimensions are in mm

Figure 47: Stress distribution of  yy in a miniature

Figure 46: Miniature dumb-bell specimen geometry along with an inclined edge crack

dumb-bell specimen in the presence of inclined edge crack

Crack Trajectories

Standard XFEM (da = 2 mm) VNXFEM (da = 1 mm) VNXFEM (da = 2 mm)

Figure 48a: Crack growth trajectory of an inclined edge crack in a miniature dumb-bell specimen for a mesh size of 1188 elements

50

Crack Trajectories

Standard XFEM (No. of element=2794) VNXFEM (No. of element=1188) VNXFEM (No. of element=2794) Figure 48b: Crack growth trajectory of an inclined edge crack in a miniature dumb-bell specimen for a crack growth increment of da = 1 mm

3000

2500

VNXFEM (da = 0.5 mm) VNXFEM (da = 1 mm) VNXFEM (da = 2 mm) Standard XFEM (da = 2 mm)

2500 680

2000

SIF (MPa.mm1/2)

1/2

SIF (MPa.mm )

SIF (MPa.mm1/2)

660 VNXFEM (da = 0.5 mm) 640 2000 VNXFEM (da = 1 mm) VNXFEM (da = 2 mm) Standard XFEM (da620 = 2 mm) 600

1500

580 560

VNXFEM (No. of element: 1188) VNXFEM (No. of element: 2794) Standard XFEM (No. of element: 2794)

1000

1500

1000

540 520

500

500 17.6

17.7

0

8

17.8

10

17.9 18 18.1 18.2 Crack Length (mm)

12

14

18.3

18.4

18.5

16 18 20 Crack Length (mm)

22

24

26

0

28

Figure 49a: SIF variation with crack length of a miniature dumb-bell specimen in the presence of inclined edge crack for a mesh size of 1188 elements

51

8

10

12

14

16 18 20 Crack Length (mm)

22

24

26

Figure 49b: SIF variation with crack length of a miniature dumb-bell specimen in the presence of inclined edge crack for a crack growth increment (da) of 1 mm

Tables

Table 1(a): A comparison of VNXFEM results with standard XFEM and experimental results for mesh size of 30x40

da/dN (  105 ) (mm/cycle)

K eq (MPa m) Crack VNXFEM VNXFEM Standard XFEM length (da: 2.0 (da = 4.0 (da = (mm) mm) mm) 4.0 mm) 45 15.9 15.8 15.8

Ref.

VNXFEM VNXFEM (da = 2.0 (da = 4.0 mm) mm)

15.9

0

0

Standard XFEM Ref. (da = 4.0 mm) 0 0

47

17.4

–

–

17.5

2.8

–

–

2.9

49

19.1

19.0

19.0

19.2

3.4

3.4

3.4

3.5

51

21.1

–

–

21.2

4.2

–

–

4.3

53

23.3

23.3

23.2

23.5

5.2

5.2

5.1

5.3

55

25.8

–

–

26.2

6.4

–

–

6.7

57

28.7

28.6

28.5

29.2

8.1

8.0

7.9

8.4

59 61 63 65

32.2 36.2 41.0 46.8

– 35.9 – 46.8

– 35.8 – 46.9

32.8 37.1 42.1 48.3

10.2 13.1 17.1 22.6

– 12.9 – 22.5

– 12.8 – 22.6

11.0 14.0 18.0 24.0

Table 1(b): A comparison of VNXFEM results with standard XFEM and experimental results for crack growth increment size (da) of 2.0 mm

da/dN (  105 ) (mm/cycle)

K eq (MPa m)

45

15.9

16.0

Standard XFEM (Mesh size: 60x80) 16.0

47

17.4

17.5

17.5

17.5

2.8

2.9

2.8

2.9

49

19.1

19.1

19.2

19.2

3.4

3.4

3.4

3.5

51

21.1

21.2

21.2

21.2

4.2

4.3

4.3

4.3

53

23.3

23.4

23.4

23.5

5.2

5.2

5.2

5.3

55

25.8

26.0

26.0

26.2

6.4

6.5

6.5

6.7

57

28.7

28.9

28.9

29.2

8.1

8.2

8.2

8.4

59 61 63 65

32.2 36.2 41.0 46.8

32.3 36.5 41.1 47.1

32.3 37.6 41.3 47.2

32.8 37.1 42.1 48.3

10.2 13.1 17.1 22.6

10.4 13.4 17.1 22.8

10.3 14.2 17.3 22.9

11.0 14.0 18.0 24.0

Crack VNXFEM VNXFEM (Mesh (Mesh length size: size: (mm) 30x40) 60x80)

Ref.

VNXFEM VNXFEM (Mesh (Mesh size: size: 30x40) 60x80)

15.9

0

0

52

Standard XFEM (Mesh Ref. size: 60x80) 0 0

Table 2: Crack growth increment, mesh-size and element size for different cases CASE-1: Same mesh size with different crack growth increments Standard XFEM Crack growth increment (mm) Mesh size (nodes  nodes) Element size (mm  mm)

VNXFEM

2.0

0.5

1.0

2.0

30  60 1.72  1.69

CASE-2: Same crack growth increment with different mesh sizes Standard XFEM Mesh size (nodes  nodes) Element size (mm  mm) Crack growth increment (mm)

VNXFEM

60  120

20  40

30  60

60  120

0.85  0.84

2.63  2.56

1.72  1.69

0.85  0.84

1.0

Table 3: Properties of the constituents of bi-materials Young Modulus

E (GPa)

Poisson's ratio, 

Fracture Toughness K IC (M Pa m)

Material-1

74

0.30

40

Material-2

200

0.30

60

53

Highlights (for review)

Research Highlights 

A new virtual node extended finite element method (VNXFEM) has been proposed for modeling and simulation of kinked cracks in a single element.



The kinked crack growth in the element is performed using the concept of virtual nodes to model and improve the accuracy of the solution.



Special five-node transition elements are used to ensure the continuity in the displacement fields.



A polygon law is proposed to determine the positions of the evaluation points with respect to a crack face in the virtual split elements.



The VNXFEM models the crack growth accurately with coarse mesh and small crack growth increments.