Dynamic Crack Analysis in Isotropic/Orthotropic Media via Extended Isogeometric Analysis

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.


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 [3][4][5][6][7][8][9][10][11][12][13][14][15].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 elementbased 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 [17][18][19][20][21][22][23][24], solid mechanics [25][26][27], fluid mechanics [28], and contact mechanics [29].
Recently, the IGA has been improved with enrichment functions to solve linear elastic fracture mechanics problems [30][31][32] 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.

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; B , 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 u 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 u  is the displacements of control point .

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   (x), The crack-tip branch enrichment functions   (x) for orthotropic materials are defined as [36] {  (x)} where  and  are the crack-tip local polar coordinates, and where   and   ( = 1, 2) are crack-tip material parameters and   =   +   ( = 1, 2) 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 cracktip branch enrichment functions   (x) for isotropic materials are defined as [1] {  (x)} where u is the displacement vector;  and  are the stress and strain tensors, respectively; f and f 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  = [u a b]  is the vector of the unknown variable at control points and M, K, and F are, respectively, the global mass matrix, stiffness matrix, and external force vector at control points.
The element contribution to M is where and the element contribution to K is where with Additionally, the element contribution to F is with where t and t are the external traction and body force vectors, respectively.

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].

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  ≥ 0.5 and  ≥ 0.25( + 0.5) 2 .In this study,  = 0.5 and  = 0.25 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).
( (5) Repeat the loop for the next time step until the maximum time step is reached.

DSIF Calculations
The DSIFs are evaluated by using the domain form of the contour interaction integral technique.The following states are considered: state 1 ( (1)   ,  (1)   ,  (1)   ), which corresponds to the actual state; state 2 ( (2)   ,  (2)   ,  (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]  (1,2) = ∫  [− (1,2)  1 +  (1) where  (1,2) =  (1)    (2)   =  (2)    (1)   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.

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,  =  = 3. 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.

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/m 3 .The tensile stress  0 = 500 MPa.The theoretical solution of mode-I DSIF for the stationary crack is expressed as [38]  dyn where   = (/  ) (  : the dilatational wave speed).

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 1.5  ≤  ≤ 3  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.

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 meshfree 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 cracktip enrichments are very accurate during 1.5  ≤  ≤ 3  ; however, the results are relatively poor during  ≤ 1.5  .The X-IGA results are more accurate than other numerical results during 1.5  ≤  ≤ 3  .

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 2 = 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/m 3 .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.

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  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.

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 2 = 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 Δ = (ℎ/50  ) is selected.Two different inclination angles between the  1 materialaxis 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.

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/m 3 .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.

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 Tsplines 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.

Figure 1 :
Figure 1: Geometry and loading of a semi-infinite crack.

Figure 2 :
Figure 2: Control points and elements for degree 3 X-IGA with 27 × 9 elements (the blue cross represents the control points enriched by Heaviside function; the blue square represents the control points enriched by crack tip enrichment functions; the red line represents the crack).

Figure 3 :
Figure 3: Convergence of the normalized mode-I DSIFs versus mesh densities.

Figure 4 :
Figure 4: Convergence in percentage errors for the normalized mode-I DSIFs versus mesh densities during 1.5  ≤  ≤ 3  .

5 Figure 5 :Figure 6 :
Figure 5: Comparison of the normalized mode-I DSIFs among different methods.

5 Figure 7 :
Figure 7: Comparison of the normalized mode-I DSIFs obtained by the X-FEM, the sES-FEM, and the X-IGA for three crack inclination angles.

Figure 8 :
Figure 8: Comparison of the normalized mode-II DSIFs obtained by the X-FEM, the sES-FEM, and the X-IGA for three crack inclination angles.

Figure 9 :
Figure 9: Geometry and loads of an edge crack in a rectangular orthotropic plate.

Figure 10 :
Figure 10: The X-IGA discretization for an edge crack in a rectangular orthotropic plate.

5 Figure 11 :
Figure 11: Comparison of the normalized mode-I DSIFs among different methods.

Figure 12 :
Figure 12: Comparison of the normalized mode-I DSIFs for different meshes.

Figure 15 :Figure 16 :
Figure 15: Normalized mode-II DSIFs for a single central crack in a rectangular plate.

Figure 17 :Figure 18 :
Figure 17: Meshes and control points.(The cross represents the control points enriched by Heaviside function; the square represents the control points enriched by crack tip enrichment functions.)

Figure 19 :
Figure 19: Comparison of the normalized mode-I DSIFs between the X-FEM and the X-IGA.
(x), and   (x) are the shape functions of standard IGA (NURBS basis functions); u  , a  , and b  are, respectively, the displacement and enrichment variable vectors at control point;   ,  cut , and  tip are, respectively, the set of all control points in the computational domain, the set of control points enriched with a modified Heaviside step function (x), and the set of control points enriched with the crack-tip branch enrichment functions   (x); ( = 1, . .., 4) ; and the basis function support of the control point in  cut is completely split by the crack, whereas that of the control point in  tip is partly split by the crack.The modified Heaviside step function (x) is given by