Abstract

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.

1. 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-the-thickness 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 [13]. 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 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 time-consuming. 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 [3133]. 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.

2. 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 to ; 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 are, respectively, now written as

Figure 2 illustrates comparison of the effective Young modulus of Al/ZrO2 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 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 is a nodal displacement vector, is a load vector, andwhere is the strain-displacement matrix which contains the gradients of the interpolating functions; is a constitutive matrix variable; is the domain of element . In the present work, the elasticity matrix was assumed to be a function of spatial coordinates.

The integral in (6) was evaluated by Gauss quadrature, and was specified at each Gaussian integration point. Thus the integral for two-dimensional problems becomeswhere the subscripts and refer to the Gaussian integration points, is the determinant of the Jacobian matrix, and is the Gaussian weight.

can be determined by interpolation, wherein the elastic modulus and Poisson’s ratio can be expressed aswhere is the shape function of FEM.

This part was written by UMAT in ABAQUS®. The file  .inp should include the following statements:MATERIAL, NAME=FGMDepvar1,User Material, constants=2200000., 0.3where 200000 and 0.3 are the initial values of and , respectively.

3. 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 [3638]. 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 () and () 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 arewhere () and () 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 () and () 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 ( and ) 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) are projected into the local coordinate system as follows: 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 () and SERR () is expressed as follows:where plane stress ; plane strain ; and are the modulus of elasticity and Poisson’s ratio in the crack tip, respectively.

This part was written by UEL in ABAQUS. The  .inp file should include the following statements:USER ELEMENT, NODES=5, TYPE=U600, PROPERTIES=3, COORDINATES=2, VARIABLES=91,2ELEMENT, TYPE=U600, ELSET=CRACKTIP60001, 13, 9013, 12, 9014, 14UEL PROPERTY, ELSET=CRACKTIP200.0E6, 200.0E6, <B>

In the above statements, 60001 is the element number; 13, 9013, 12, 9014, and 14 are the node numbers, 200.0E6, 200.0E6, and <B> are the -direction spring stiffness, -direction spring stiffness, and thickness of the plate, 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 arewherewherewhere is the dilatational speed; is the shear wave speed; = for plane strain which is the case in examples and ; 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:USER ELEMENT, NODES=5, TYPE=U600, PROPERTIES=11, COORDINATES=2, VARIABLES=91,2ELEMENT, TYPE=U600, ELSET=CRACKTIP60001, 13, 9013, 12, 9014, 14UEL PROPERTY, ELSET=CRACKTIP74.588E12, 74.588E12, <B>, 0.0, <C>, 0.0, 2, 0.0, 32.0, 12.8, 44.8

In the above statements, 0.0, <C>, and 0.0, 2 are the amount of crack propagation, speed of crack propagation, initial crack length, and number of propagation steps; 0.0 and 32.0 are the - and -direction coordinates of the initial crack, respectively; 12.8 and 44.8 are the - and -direction coordinates of the final crack, respectively.

4. Numerical Results

4.1. 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.2, 0.5, 1.0, and 5.0 and the elastic moduli of Y-TZP and -TiAl are  GPa and  GPa, respectively; Poisson’s ratios of two materials are both 0.3 and the mass density = 5 × 103 kg/m3. The load form is plotted in Figure 5, the amplitude  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 . 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 of FGM plate of four different grid models solved by our method at μs are plotted in Figure 6, when the shape factor , 0.2, 0.5, 1.0, and 5.0. At , 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 , 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.

4.2. 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 (Figure 7). Two types of FGMs Y-TZP and -TiAl were used. The components were subjected to (1)-(2), where the shape factor , 0.2, 0.5, 1.0, and 5.0, the elastic moduli of Y-TZP and γ-TiAl are  GPa and  GPa, respectively; Poisson’s ratios ν of two materials are both 0.3, and the mass density ρ = 5 × 103 kg/m3. The load form is plotted in Figure 5, the amplitude  GPa, undamped. The normalized DSIFs and 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 and of FGM plate solved by our method are plotted in Figures 8 and 9, respectively, when the shape factor , 0.2, 0.5, 1.0, and 5.0. At , 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 , 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 deformation patterns of inclined-crack functionally graded plate in different methods and time steps are plotted in Figure 10. Visibly, the new method is correct and effective.

4.3. Example 3

The third example is a rectangular panel with a central crack as shown in Figure 11. Two FGMs Y-TZP and -TiAl were used. The components were subjected to (1)-(2), where the shape factor , 0.2, 0.5, 1.0, and 5.0 and the elastic moduli of Y-TZP and -TiAl are  GPa and  GPa, respectively. The initial crack starts to propagate in a self-similar manner along the dashed line shown in Figure 11. The plate was subjected to a time-independent tensile stress at the edges. The constant crack velocity was  m/s. The normalized DSIF solution computed here is shown in Figure 12. The results developed here at (Figure 12) are found to be nearly identical to the numerical solutions by Xie et al. [54]. The load form is plotted in Figure 5, at the amplitude  GPa, undamped. The normalized is plotted in Figure 12 and compared with Xie et al. [54]. The solutions obtained by the two methods are similar. Figure 13 shows the mesh before deformation and at several time sets.

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 VCCT interface 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  m/s are shown in Figure 12. Clearly, the normalized remains almost constant throughout the crack propagation, while the normalized gradually decreases. Moreover, the DSIFs drop immediately after the crack propagation begins. Figure 13 shows the deformation during the dynamic crack propagation.

4.4. 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 (,  mm) 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  Cs. The components were subjected to (1)-(2), where the shape factor , 0.2, 0.5, 1.0, and 5.0 and the elastic moduli of Y-TZP and -TiAl are  GPa and  GPa, respectively; Poisson’s ratios of two materials are both 0.3, and the mass density is = 5 × 103 kg/m3. 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 was employed, without encountering convergence difficulty. No damping was included in the model.

Variations of normalized DSIFs during the crack propagation at  Cs are shown in Figures 15 and 16. Clearly, the normalized remains almost constant throughout the crack propagation, while the normalized gradually decreases. Moreover, the DSIFs drop immediately upon the crack propagation. Figure 17 shows the deformation during the dynamic crack propagation. The results obtained by the new method agree well with those by Qian and Xie [38] using moving singular elements.

5. 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 propagation problems by resorting to the commercial FEA codes with user subroutines.

Competing Interests

The authors declare that there are no competing interests regarding the publication of this paper.

Acknowledgments

This work was financially supported by the National Natural Science Foundation of China (Grant no. 51305157), the National Key Scientific Instrument and Equipment Development Projects, China (Grant no. 2012YQ030075), and Jilin Provincial Department of Science and Technology Fund Project (Grant nos. 20160520064JH, 20130305006GX).