Abstract

The extended isogeometric analysis (X-IGA) is the combination of the extended finite element method (X-FEM) and the isogeometric analysis (IGA), so the X-IGA possesses the advantages of both methods. In this paper, the X-IGA is extended to investigate the dynamic stress intensity factors of cracked isotropic/orthotropic media under impact loading. For this purpose, a corresponding dynamic X-IGA model is developed, the Newmark time integration scheme is used to achieve a dynamic response, and the dynamic stress intensity factors are evaluated through the contour interaction integral technique. Numerical simulations show that the X-IGA results agree with other available reference solutions, and accurate results can be obtained by using the X-IGA with a relatively coarse mesh.

1. Introduction

Among existing numerical methods, the extended finite element method (X-FEM), which was developed in 1999 by Belytschko and his coworkers [1, 2], is the most effective method for solving discontinuity problems. In the X-FEM, the geometry of the discontinuity is independent of the computational mesh; therefore, the computational mesh does not need to be updated in the simulation of the propagation of discontinuities. In the past decades, the extensive literature on improving and applying the original X-FEM in modeling discontinuities has been published [315]. The X-FEM aims to enrich standard finite element approximation by using discontinuous basis functions in the framework of partition of unity.

The standard finite element approximation is element-based polynomial approximation, and it yields discretization errors in complex geometry. In order to overcome the drawback, the isogeometric analysis (IGA) was proposed by Hughes et al. [16] in 2005. The principle of the IGA is that nonuniform rational B-splines (NURBS) basis functions are employed as shape functions for geometric description and field approximation. The IGA has some unique advantages [16], so it has been developed and applied in many fields including structural mechanics [1724], solid mechanics [2527], fluid mechanics [28], and contact mechanics [29].

Recently, the IGA has been improved with enrichment functions to solve linear elastic fracture mechanics problems [3032] and curved interface problems [33]. The extended isogeometric analysis (X-IGA) was proposed in [31] and contains the inherent advantages of both IGA and X-FEM. Currently, the X-IGA is further developed to analyze cracked orthotropic media [34]. In this study, the X-IGA is extended to investigate the dynamic fracture behavior of stationary cracks in isotropic/orthotropic media. A corresponding dynamic X-IGA model is developed, and the Newmark time integration scheme is used to achieve a dynamic response. The dynamic stress intensity factors (DSIFs) are evaluated by using the contour interaction integral technique.

The paper is organized as follows. Section 2 briefly reviews the fundamentals of NURBS-based IGA. The X-IGA for dynamic cracked isotropic/orthotropic media is described in Section 3. The DSIFs are derived by using the contour interaction integral technique in Section 4. Section 5 presents the numerical results obtained via X-IGA and compares the results with other solutions. Finally, conclusions and prospects are drawn in Section 6.

2. NURBS-Based Isogeometric Analysis

In the CAD, a two-dimensional NURBS surface can be constructed as [35] where is the NURBS basis function; and are, respectively, the orders of basis functions in the direction and the direction; represents the coordinate of control point; and and are the numbers of basis functions in the direction and the direction, respectively.

In NURBS-based isogeometric analysis [16], the displacement field is approximated similarly to the isoparametric finite element method as follows: where is the number of the control points; is the parametric coordinate; denotes the NURBS basis function at control point ; and is the displacements of control point .

3. X-IGA for Dynamic Cracked Media

3.1. Displacement Approximation

Similar to the X-FEM, the X-IGA aims to enrich the standard IGA approximation by using additional functions on the basis of the partition of unity to model discontinuities. For crack problems, the enriched displacement approximation can be expressed as where , , and are the shape functions of standard IGA (NURBS basis functions); , , and are, respectively, the displacement and enrichment variable vectors at control point; , , and are, respectively, the set of all control points in the computational domain, the set of control points enriched with a modified Heaviside step function , and the set of control points enriched with the crack-tip branch enrichment functions ; ; and the basis function support of the control point in is completely split by the crack, whereas that of the control point in is partly split by the crack.

The modified Heaviside step function is given by

The crack-tip branch enrichment functions for orthotropic materials are defined as [36]where and are the crack-tip local polar coordinates, and where and are crack-tip material parameters and    are the roots of where is the compliance coefficients.

Equation (5) cannot be directly used for isotropic materials because of the presence of 0/0 in the equation. The crack-tip branch enrichment functions for isotropic materials are defined as [1]

3.2. Discrete Equilibrium Equations

Considering a body with an initial traction-free crack in the state of dynamic equilibrium, the weak form of the momentum equation is where is the displacement vector; and are the stress and strain tensors, respectively; and are the body force and external traction vectors, respectively; and is the mass density.

The discretized form of (9) using the X-IGA approximation (3) can be written as where is the vector of the unknown variable at control points and , , and are, respectively, the global mass matrix, stiffness matrix, and external force vector at control points.

The element contribution to is where and the element contribution to is where with

Additionally, the element contribution to is with where and are the external traction and body force vectors, respectively.

3.3. Numerical Integration Scheme

The Gauss quadrature scheme is employed in the X-IGA. To obtain an accurate integration for crack tip elements and elements cut by crack, the triangular subdomain technique is used in the same way as that of the X-FEM [1]. For additional details, refer to [31].

3.4. Time Integration Scheme

The Newmark method is adopted for the time integration of (10). At time step , the discrete simultaneous equations are described as where is the time step and the unconditionally stable conditions are and . In this study, and are chosen.

The implementation procedure of the Newmark time integration scheme applying into the X-IGA is outlined as follows.(1)The initial velocities and initial displacements are set to be zero; then the initial accelerations are obtained with (10).(2)Compute the constants , , , , , , , and .(3)Compute .(4)Loop over the time steps as follows:(a)compute ;(b)compute the displacements by solving ;(c)solve for the value , and the velocities are calculated by using .(5)Repeat the loop for the next time step until the maximum time step is reached.

4. DSIF Calculations

The DSIFs are evaluated by using the domain form of the contour interaction integral technique. The following states are considered: state 1 (, , ), which corresponds to the actual state; state 2 (, , ), which is an auxiliary state that will be selected as the asymptotic field for model I or II. The interaction integral may be written as [37] where is the interaction strain energy and is a weighing function. The second term of the right-hand side of the equation is the contribution of the inertia forces.

For isotropic materials, the relation between the interaction integral and the stress intensity factors is expressed as [37] where is for plane stress and is for plane strain.

The relation between the interaction integral and the stress intensity factors for orthotropic materials is expressed as [36] with

Choosing state 2 as mode I or mode II leads to mode-I or mode-II DSIF in terms of the interaction integral.

For isotropic materials,

For orthotropic materials,

5. Numerical Simulations

In this section, several examples of stationary cracks in isotropic/orthotropic media with available reference solutions are investigated to assess the accuracy of the proposed approach. The first two examples demonstrate the efficiency of the X-IGA for isotropic material problems, and the next two illustrate the application of the X-IGA to orthotropic material problems, whereas the last example demonstrates the capability and versatility of the X-IGA in modeling the complicated geometries. In all examples, the Heaviside step loading is considered and degree 3 X-IGA is adopted; that is, . Meshes are generated with linear parameterization method. In the case of investigating the effects of different meshes on the results, different uniform meshes are adopted; in other cases, the fine meshes are used around the crack, while the coarse meshes are used in other domains; thus excellent accuracy can be achieved at a low cost.

5.1. Semi-Infinite Crack in an Infinite Isotropic Plate

For the first example, we consider an infinite isotropic plate with a semi-infinite crack subjected to a tensile stress wave perpendicular to the crack face, as shown in Figure 1. The material properties of the media are as follows: Young’s modulus 210 GPa, Poisson’s ratio 0.3, and mass density 8000 kg/m3. The tensile stress 500 MPa. The theoretical solution of mode-I DSIF for the stationary crack is expressed as [38] where   (: the dilatational wave speed).

5.1.1. Convergence Study of the DSIF versus Meshes

In this study, five meshes with 27 × 9, 47 × 17, 67 × 25, 87 × 33, and 107 × 41 elements are considered. Figure 2 illustrates the typical regular mesh of 27 × 9 elements.

Figure 3 presents the normalized mode-I DSIFs as a function of time for the five considered meshes. In the mesh refinement process, the mode-I DSIFs obtained by the X-IGA converge well to the theoretical solution. The maximum error increases beside for the coarse mesh and decreases significantly for the fine mesh. Similar observation was noted in the simulation of the singular edge-based smoothed finite element method (sES-FEM) [37]. The same simulations were conducted by other researches by using different methods. However, the same conclusion is recorded; that is, the error is lower in later stages than in earlier stages. The relative errors during are presented in Figure 4, in which the error is unchanged as the mesh is refined to a certain extent. In this study, the errors are the same for meshes with 67 × 25, 87 × 33, and 107 × 41 elements.

5.1.2. Comparison with X-FEM Results

The X-IGA results using 67 × 25 mesh elements are compared with three available X-FEM solutions, including the X-FEM using mass lumping (120 × 60 quadrilateral elements) [39], the mesh-free enriched X-FEM with crack-tip enrichments (78 × 39 quadrilateral elements) [38], and X-FEM with a new enrichment function (120 × 59 quadrilateral elements) [40]. Figure 5 shows the comparative study of the normalized mode-I DSIFs. The results obtained with X-FEM using mass lumping [39] are unstable, and the accuracy of this method is the worst recorded among the four approaches. The results obtained with the X-IGA, the X-FEM using an enrichment function, and the mesh-free enriched X-FEM with crack-tip enrichments are very accurate during ; however, the results are relatively poor during . The X-IGA results are more accurate than other numerical results during .

5.2. Arbitrarily Oriented Central Crack in a Rectangular Isotropic Plate

We consider a rectangular isotropic plate with an arbitrarily oriented central crack subjected on the top and bottom of a uniform impact loading, as shown in Figure 6(a). The length of the crack is 4.8 mm. The material properties of the plate are as follows: Young’s modulus 200 GPa, Poisson’s ratio 0.3, and mass density 5000 kg/m3. Three crack inclination angles are investigated, and the total time of the simulation is 20 μs. In the X-IGA, the computational mesh is independent of the crack; therefore, the computation meshes for the three crack inclination angles are the same, as shown in Figure 6(b).

Figures 7 and 8, respectively, present the normalized mode-I and mode-II DSIFs for the three considered crack angles computed by the X-IGA, sES-FEM [37], and X-FEM. Figure 6(b) is also the mesh for the X-FEM, whereas the mesh of 60 × 120 triangular elements is used in the analysis of the sES-FEM.

A good agreement can be observed for the three considered crack inclination angles among the three numerical approaches.

5.3. Edge Crack in a Rectangular Orthotropic Plate

As a third example, a rectangular orthotropic plate with a 12 mm edge horizontal crack is considered (Figure 9). The plane stress condition is considered, and the material properties of the plate are as follows: 118.3 GPa, 54.8 GPa, 8.79 GPa, 0.083, and 1900 kg/m3. The time step is selected, where is the wave velocity along the material-axis.

Figure 10 shows the mesh with 75 × 57 elements used in the simulation. The normalized mode-I DSIFs achieved with the X-IGA (Figure 11) were compared with the reference boundary element method (BEM) (21 elements for the external boundary and eight elements for the crack) [41], the X-FEM (78 × 60 quadrilateral elements) [42], and the conventional finite element method (FEM) results (via ANSYS) [41]. Good agreement is observed among the methods. The X-IGA results more closely match the FEM and the BEM solution compared with the X-FEM results. In [42], the DSIFs are evaluated by using the domain separation integral method, which could be the reason that the X-FEM results are different.

To analyze the effects of different meshes on DSIFs, five different meshes with 27 × 19, 47 × 39, 75 × 57, 97 × 81, and 121 × 101 elements are adopted to simulate the problem. The results are shown in Figure 12. Accurate results can be obtained by using the X-IGA with a relatively coarse mesh.

5.4. Central Crack with Different Orientations of the Axes of Orthotropy

A rectangular orthotropic plate with a single central horizontal crack was subjected to a Heaviside step tensile distributed load, as shown in Figure 13(a). The length of the crack is 4.8 mm and 20 mm. The material properties of the plate are the same as those in the example in Section 5.3. The plane stress state is assumed, and the time step is selected.

Two different inclination angles between the material-axis and the crack-face are considered. Figures 14 and 15 show the comparison of the numerical results of BEM (24 elements for the external boundary and 10 for the crack) [41] and X-FEM (50 × 100 quadrilateral elements) [42]. General trends are in good agreement across all methods; the X-IGA results more closely match the BEM solution compared with the X-FEM; the reason for this phenomenon is the same as the above example.

5.5. Edge Crack in an Annular Isotropic Plate

To illustrate the capability and versatility of the X-IGA in modeling the complicated geometries, the last example deals with an annular isotropic plate under a uniform pressure. The geometry and loads of the plate are depicted in Figure 16. The length of the crack is 2 m; the material properties of the plate are as follows: Young’s modulus 210 GPa, Poisson’s ratio 0.3, and mass density 8000 kg/m3. The total time of the simulation is 0.004s. The plane strain state is assumed, and the time step 8 × 10−5 s is selected.

The meshes are depicted in Figure 17 for 285 control points and 221 elements. To perform a comparison, this example is solved by the X-FEM, and the X-FEM meshes are shown in Figure 18 for 1292 nodes and 1221 quadrilateral elements.

The DSIFs are computed with the X-FEM and the X-IGA, and the results are shown in Figure 19. It is found that the results from the X-IGA with a relatively coarse mesh are in good agreement with those from the X-FEM with a very fine mesh. The amplitude of the normalized mode-I DSIF is negative in some time ranges, which may be due to the contact between two crack-faces which is not taken into consideration in the present study.

6. Conclusions and Prospects

The NURBS-based extended isogeometric analysis (X-IGA) was extended to investigate the dynamic fracture behavior of stationary cracks in isotropic/orthotropic media. A corresponding dynamic X-IGA model was developed, the Newmark time integration scheme was used to achieve the dynamic response, and the DSIFs were evaluated with the contour interaction integral technique. Numerical results indicate that accurate results can be obtained by using the X-IGA with a relatively coarse mesh, and the X-IGA is suitable to model the cracked complicated geometries.

The NURBS-based isogeometric analysis has several drawbacks including the handling trimmed geometries and the local refinement. The isogeometric analysis based on T-splines can effectively overcome these drawbacks [43]. A key advantage of the extended isogeometric analysis over isogeometric analysis is the ability to model crack propagation. Therefore, the application of extended isogeometric analysis based on T-splines to dynamic crack propagation analysis in isotropic/orthotropic media is an area of considerable promise. On the other hand, the smoothed finite element method (SFEM) developed by combining the smoothing technique with the finite element method is an efficient and accurate numerical simulation tool for the dynamic fracture problems [44]; therefore combining the smoothing technique with the X-IGA for fracture analysis is another area of considerable promise.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgment

This work is supported by the National Natural Science Foundation of China (51179063). The financial support is gratefully acknowledged.