Analysis of Dynamic Fracture Parameters in Functionally Graded Material Plates with Cracks by Graded Finite Element Method and Virtual Crack Closure Technique

Based on the finite element software ABAQUS and graded element method, we developed a dummy node fracture element, wrote the user subroutines UMAT and UEL, and solved the energy release rate component of functionally graded material (FGM) plates with cracks. An interface element tailored for the virtual crack closure technique (VCCT) was applied. Fixed cracks and moving cracks under dynamic loads were simulated. The results were compared to other VCCT-based analyses. With the implementation of a crack speed function within the element, it can be easily expanded to the cases of varying crack velocities, without convergence difficulty for all cases. Neither singular element nor collapsed element was required. Therefore, due to its simplicity, the VCCT interface element is a potential tool for engineers to conduct dynamic fracture analysis in conjunction with commercial finite element analysis codes.


Introduction
Functionally graded materials (FGMs) are composites formed of two or more constituent phases with a continuously variable composition.They are attractive in potential applications owing to their numerous advantages, including the reduction of in-plane and transverse through-thethickness stresses and stress intensity, and the improvement in residual stress distribution, thermal properties, and fracture toughness.There are a number of reviews dealing with various aspects of FGMs in recent years [1][2][3].For instance, an embedded crack in an orthotropic FGM layer is considered in the case of mechanical loading [4].Diverse areas relevant to the theory and applications of FGMs are reflected in [5], including homogenization of particulate FGM, heat transfer, stress, stability and dynamic analyses, manufacturing and design, applications, and fracture.Many aspects of FGMs such as free vibration [6], shear deformation [7], thermal buckling [8], and stress intensity factor [9] have been investigated.Moreover, the fracture toughness of functionally graded (FG) sections is of interest especially for a material with elastic behavior [10,11].Using a plate bending finite element based on FOST, Singha et al. studied the nonlinear behaviors of FG plates under transverse load, considering the physical/exact neutral surface position and assuming the power law gradation of material properties in the thickness direction [12].Isogeometric analysis was also very promising to be applied to a wide range of practical mechanics problems such as laminated composite and sandwich plates based on inverse trigonometric shear deformation theory, functionally graded plates based on generalized shear deformation theory [13].Until now, FGM is one predominant mode of material and has been investigated extensively.
In recent years, there are growing concerns on how cracked functional material body responds to collision under impulse loading.To accurately evaluate the fracture mechanics under dynamic loading, researchers proposed dynamic fracture parameters, such as dynamic stress intensity factor (DSIF) and strain energy release rate (SERR).The dynamic 2 Advances in Materials Science and Engineering fracture parameter of simple geometric model, ideal material model, or special load model can be determined by the analytical method.However, this method is not applicable to complex structure or boundary conditions, and its experimental measurements are very expensive and timeconsuming.Nevertheless, this type of problems can be well resolved by numerical calculations.
At present, the finite element method (FEM) is widely used for fracture analysis in FGMs.For instance, a pair of FEM-based elastodynamic contour integrals was developed to calculate the elastodynamic asymptotic mixed-mode stress field for plane elastic materials containing a stationary notch tip [16].Graded finite elements can be used in fracture analysis in FGMs where the elastic moduli are smooth functions of spatial coordinates, which are integrated into the element stiffness matrix.The stress intensity factors for mode I and mixed-mode two-dimensional problems can be comparatively evaluated through three FGMs-tailored approaches: path-independent J-integral, modified crack closure integral, and displacement correlation [17].The feasibility of FEM in cracked or uncracked FGM plates was studied.The J contour integral of ABAQUS was used to calculate stress intensity factors for an edge cracked FGM plate [18].Matthews used the finite element analysis (FEA) for large displacement J-integral test to analyze mode I interlaminar fracture in composite materials [19].The dynamic crack tip fields were determined, and the crack propagation of anisotropic materials was also characterized [20].These previous works are important; however, they only focus on the dynamic cracks of isotropic and orthotropic materials, but not on the direction of crack propagation.
The methods used to resolve the fracture parameters include J-integral, M-integral, extrapolation, and virtual crack closure technique (VCCT).Among all fracture parameters, SERR is used increasingly in conjunction with linear elastic fracture mechanics (LEFM) and can be computed by VCCT together with FEA.VCCT requires a preexisting crack with a sharp neat tip within a material for crack initiation as well as conditions of small-scale yield to hold.With material nonlinearity at the crack tip (small process zone) ignored, LEFM-based approaches were proven effective in predicting crack initiation and subsequent growth [21,22].
VCCT was proposed for 2D crack configurations [23] and extended to 3D-VCCT later [24].Recently, the VCCT formulation for kinking cracks was proposed [25].Krueger summarized historical developments and discussed different applications [26].Combining 2D-VCCT and FEA, Sun and Qian compared the SERRs of interfacial cracks between two isotropic materials [27].As reported, the SERRs for the delamination between the face sheet and the core material of sandwich structures were calculated [28,29].Glaessgen et al. calculated SERRs to evaluate the suppressing effect of stitching on debonding [30].VCCT was also applied to electronic packaging [31][32][33].Leski used VCCT to study the interface crack propagation [34].Ramu used a differential transform method (DTM) to study free transverse vibration of isotropic rectangular plates resting on a Winkler foundation [35].Modified crack closure integral technique was extended to the element-free Galerkin method [36].
A cohesive theory assumes the presence of a process zone in front of the crack tip whose fracture properties consist of upper and lower surfaces controlled by the cohesive traction-displacement discontinuity relationship and allows non-self-similar crack propagation [37].An automated fracture procedure implemented in the large-scale, nonlinear, and explicit, finite element code DYNA3D can be used to simulate dynamic crack propagation in arbitrary directions [38].Manolis et al. used boundary element method (BEM) to analyze the dynamic fracture of a smoothly inhomogeneous and defective plane [39].Solving crack growth problems, the recent approach on smoothed finite element methods is really a good candidate [40,41].The DSIF around the antiplane crack in an infinite strip FGM under impact loading was investigated [42].FG cracked plates under different loads and boundary conditions were numerically simulated using NURBS-based XIGA [43].XIGA has been applied to stationary and propagating cracks in 2D [44], plastic collapse load analysis of cracked plane structures [45], and cracked plate/shell structures [46].
At present, the emerging computing method is strongly pertinent, nonversatile, and difficult to promote.Analysis of dynamic crack problems based on secondary development of ANSYS, ABAQUS, and so on is mainly focused on homogeneous materials, but it should be further expanded into FGMs.
In this study, based on the commercial FEA software ABAQUS and graded element method, we developed a dummy node fracture element, wrote the user subroutines UMAT and UEL, and solved the energy release rate component of cracked FGM plates.

FGMs
FGMs are often formed by two or more materials whose volume fractions change continuously along certain dimensions of the structure (Figure 1) [22].The effective moduli of two constituents are homogenized by the rule of mixture or the Mori-Tanaka models which are used to evaluate the effective elastic properties of the grade composite.The effective property is expressed by a power law of volume fraction exponent as follows: where subscripts  and  refer to the metal and ceramic components, respectively;  is the thickness coordinate and varies from −ℎ/2 to ℎ/2;  is the power law index;   and   denote the material properties of ceramic and metal, respectively, including Young's modulus, Poisson's ratio, and density.Equation ( 1) denotes the volume fraction variation versus nondimensional thickness (/ℎ) with different .However, the rule of mixture does not reflect the interactions among the two materials [47,48].Meanwhile, the Mori-Tanaka model [49] is assumed to calculate their interactions through the effective bulk and shear modulus given by The effective Young modulus   and Poisson's ratio V  are, respectively, now written as Figure 2 illustrates comparison of the effective Young modulus of Al/ZrO 2 FGM plate calculated by the rule of mixture and the Mori-Tanaka scheme via the power index .Note that with homogeneous material the two models produce the same values.For inhomogeneous material, the effective property through the thickness of the former is Advances in Materials Science and Engineering higher than that of latter.Moreover, increasing in power index  leads to decrement of the material property due to the rise of metallic volume fraction.In this paper the power law distribution of constituent materials along the plate thickness is assumed, and the effective homogeneous properties are calculated by the rule of mixture [15].Graded elements were implemented by directly sampling the properties at the Gauss points of each element [50,51].The graded finite element stiffness matrix relations can be written as follows [52]: where U  is a nodal displacement vector, F  is a load vector, and where B  is the strain-displacement matrix which contains the gradients of the interpolating functions; D  (x) is a constitutive matrix variable; Ω  is the domain of element .
In the present work, the elasticity matrix D  (x) = D  (, ) was assumed to be a function of spatial coordinates.The integral in (6) was evaluated by Gauss quadrature, and D  (x) was specified at each Gaussian integration point.Thus the integral for two-dimensional problems becomes where the subscripts  and  refer to the Gaussian integration points, |  | is the determinant of the Jacobian matrix, and   is the Gaussian weight.D  (x) can be determined by interpolation, wherein the elastic modulus  and Poisson's ratio ] can be expressed as where   is the shape function of FEM.This part was written by UMAT in ABAQUS5.The file .inpshould include the following statements:

VCCT Interface Element
Figure 3 shows the definition and node numbering of a typical VCCT interface element for 2D fracture problems.The details for the VCCT interfacial element can be found in [36][37][38].Specifically, each element has five nodes.Such an element is placed in a way that nodes 1 and 2 are located at the crack tip, while nodes 3 and 4 are behind and node 5 is ahead of the crack tip.The element contains two sets of nodes: a top set (nodes 1, 3, and 5) and a bottom set (nodes 2 and 4).A very stiff spring is placed between nodes 1 and 2 to compute the crack tip nodal forces as follows: where ( 1 , V 1 ) and ( 2 , V 2 ) are the displacement components of nodes 1 and 2, respectively, under the global coordinate system (, );   and   are the and -direction spring stiffness, respectively.Initially, they are assigned with large numbers [37], but once the crack is predicted to grow, they are set to zero.Dummy nodes 3, 4, and 5 do not contribute to the stiffness matrix and are introduced only to extract the information of displacement opening behind the crack tip and the crack jump length ahead of the crack tip.For nodes 3 and 4 behind the crack tip, the displacement openings are where ( 3 , V 3 ) and ( 4 , V 4 ) are the displacement components of nodes 3 and 4, respectively, under the global coordinate system (, ).Therefore, the crack jump length, which is the distance between nodes 1 and 5, is calculated as follows: where ( 1 ,  1 ) and ( 5 ,  5 ) are the global coordinates of nodes 1 and 5, respectively.If they are updated at each step, the crack orientation is also updated.This is of particular interest when large deformation cannot be neglected [45].
In order to separate the fracture modes (modes I and II), we computed the SERRs ( I and  II ) with respect to the local coordinate system (, ) attached to the crack tip as shown in Figure 1.The included angle between  and  is determined as follows: Then the nodal forces and the displacement openings in ( 9) and ( 10 Based on 2D-VCCT, the SERRs can be approximated as follows: where  is the body thickness.
For the fixed crack problem under dynamic loading, the relationship between DSIF where plane stress  =  tip ; plane strain  =  tip /(1 − (] tip ) 2 );  tip and ] tip are the modulus of elasticity and Poisson's ratio in the crack tip, respectively.For dynamic running cracks with constant velocity , their DSIFs are related to the corresponding SERRs as follows: where  is shear modulus.The crack speed functions are where where where   is the dilatational speed;   is the shear wave speed;  tip = 3 − 4] tip for plane strain which is the case in examples 2 and 3;  tip is material density at the crack tip.
The above procedure for the VCCT interface element regarding dynamic crack propagation was implemented into ABAUQUS with its user subroutine UMAT and UEL.DSIFs are directly outputted through ABAQUS state variable "SDV" and thus can be determined simultaneously during FEA on ABAQUS without any postprocessing.This part was written by UEL in ABAQUS.The .inp file should include the following statements:

Example 1.
The rectangular panel with a central crack is the first example of a convergence study in this field.As shown in Figure 4, two FGMs Y-TZP and -TiAl were used.The components were subjected to ( 1)-( 2), where the shape factor   = 0, 0.2, 0.5, 1.0, and 5.0 and the elastic moduli of Y-TZP and -TiAl are   = 200 GPa and   = 186 GPa, respectively; Poisson's ratios ] of two materials are both 0.3 and the mass density  = 5 × 10 3 kg/m 3 .The load form is plotted in Figure 5, the amplitude  0 = 0.4 GPa, undamped.
In order to study the convergence of the present method, a discrete model of four grids (I: 100 × 200 elements, II: 200 × 400 elements, III: 400 × 800 elements, and IV: 500 × 1000 elements) is simulated to the normalized DSIFs  I .The body geometry was modeled with two-dimensional ABAQUS standard plane strain elements such as "CPE4."Neither special singular elements nor the collapsed element technique was used at the crack tip.
The SERR solved by VCCT is converted to DSIFs, which are normalized by The dimensionless DSIFs  I of FGM plate of four different grid models solved by our method at  = 7 s are plotted in Figure 6, when the shape factor  = 0.0, 0.2, 0.5, 1.0, and 5.0.At  = 0, material gradient does not change at all, and the FGM is only Y-TZP.Then the result was compared with FDM and [53].Clearly, at  = 0, the result of our method is basically the same as in [53], and well consistent with the calculation of FDM.The result proves the accuracy of our method.

Example 2.
The second example is a 60 mm long and 30 mm wide rectangular plate containing a 14.14 mm long central crack inclined at an angle  = 45 ∘ (Figure 7).Two types of FGMs Y-TZP and -TiAl were used.The components were subjected to (1)-( 2), where the shape factor  = 0, 0.2, 0.5, 1.0, and 5.0, the elastic moduli of Y-TZP and -TiAl are   = 200 GPa and   = 186 GPa, respectively; Poisson's ratios ] of two materials are both 0.3, and the mass density  = 5 × 10 3 kg/m 3 .The load form is plotted in Figure 5, the amplitude  0 = 0.4 GPa, undamped.The normalized DSIFs  I () and  II () of the right crack tip are plotted in Figures 8 and 9, respectively.The solutions obtained by the two methods are similar to those of Xie and Biggers Jr. [53] who used the VCCT and finite differencing method (FDM).Figure 10 shows the mesh before deformation and at several time points.
The dimensionless DSIFs  I and  II of FGM plate solved by our method are plotted in Figures 8 and 9, respectively,   when the shape factor  = 0.0, 0.2, 0.5, 1.0, and 5.0.At  = 0, material gradient does not change at all, and the FGM is only Y-TZP.Then the result was compared with FDM and [53].Clearly, at  = 0, the result of our method is basically the same as in [53], and well consistent with the calculation of FDM.The result further proves the accuracy of our method.Besides high precision and simplicity, the new method puts forward relevant information to calculate DSIF with the help of ABAQUS, and it endows the program with strong commonality and large extension space.The elements were placed along the crack path (dashed line in Figure 11) and used as easy as the ABAQUS standard elements.An implicit dynamic analysis procedure was employed, without encountering convergence difficulty.No damping was included in the model.Variations of normalized DSIFs during the crack propagations at  = 1000 m/s are shown in Figure 12.Clearly, the normalized  I remains almost constant throughout the crack propagation, while the normalized  I gradually decreases.Moreover, the DSIFs drop immediately after the crack propagation begins.Figure 13 shows the deformation during the dynamic crack propagation.

Example 4.
An FGM plate with a slanted edge crack shown in Figure 14 was used as an example of mixed-mode dynamic crack propagation to illustrate the new method of VCCT interface element.The initial crack ( 0 cos  0 = 0.4,  = 32mm) starts to propagate in a self-similar manner along the dashed ling shown in Figure 14.The plate was subjected to a time-independent tensile stress at the edges.The crack velocity was constant at  = 0.2 Cs.The components were subjected to (1)-( 2), where the shape factor  = 0, 0.2, 0.5, 1.0, and 5.0 and the elastic moduli of Y-TZP and -TiAl are   = 200 GPa and   = 186 GPa, respectively; Poisson's ratios ] of two materials are both 0.3, and the mass density is  = 5 × 10 3 kg/m 3 .The body geometry was modeled with 2D ABAQUS standard plane strain elements such as "CPE4" and "CPE3."Neither special singular elements nor the collapsed element technique was used at the crack tip.The VCCT interface elements were placed along the crack path (dashed line in Figure 14) and used as easily as the ABAQUS standard elements.An implicit dynamic analysis procedure  the dynamic crack propagation.The results obtained by the new method agree well with those by Qian and Xie [38] using moving singular elements.

Conclusions
A 2D-VCCT interface element based on graded finite element method was presented for dynamic fracture analysis.The element was implemented into commercial FEA software ABAQUS via the user element subroutine UMAT and UEL.This new method was evaluated with three examples.The results agree well with the available numerical and experimental results in previous studies.With this interface element, fracture mechanics can be directly applied to crack growth problem on commercial FEA software.The element has several significant advantages, such as no need of extra postprocessing to extract fracture parameters, and no special burden on definition of body mesh.It can be applied in conjunction with conventional finite elements.Finally, the interface element proven reliable in a variety of cases can be employed potentially in other cases of dynamic fracture.In summary, the new 2D-VCCT interface element is simple, efficient, universal, and robust.This element method is a potential tool for structure-level analysis of dynamic crack Advances in Materials Science and Engineering 13 propagation problems by resorting to the commercial FEA codes with user subroutines.

Figure 6 :
Figure 6: The convergence of the normalized DSIFs  I at  = 7 s.

Figure 8 :
Figure 8: Dynamic stress intensity factor  I of inclined-crack functionally graded plate with different shape factors.

Figure 9 :
Figure 9: Dynamic stress intensity factor  II of inclined-crack functionally graded plate with different shape factors.
This part was written by UEL in ABAQUS.The .inp file should include the following statements: *