The Stable Explicit Time Stepping Analysis with a New Enrichment Scheme by XFEM

Xue-cong Liu, Qing Zhang,2 and Xiao-zhou Xia

1 Introduction

The fracture analysis of structures and components has been widely applied and highly valued in recent years, and modeling discontinuities like crack is one of the important parts in the simulation of failure. In order to model the crack and crack growth behavior, the way of remeshing is used by classic finite element method (FEM) in order to align the mesh with discontinuities. In addition, other solutions such as meshfree method, boundary element method and extended finite element method (XFEM) are available.

As proposed by Belytschko and Black (1999); Noes, Dolbow and Belytschko (1999);Belytschko and Noes (2001), XFEM based on the concept of partition unity becomes a dominant numerical scheme. The crack can be modeled independent of finite element mesh.All the elements are divided into the normal parts and the enriched parts. Since the elements can be influenced directly by crack, the enrich functions are introduced. Heaviside function and the Westergaard stress function are used frequently for the discontinuities and the tip’s stress singularity, respectively. XFEM is used to simplify the discontinuous problems and perform well in stress analysis concerned with fracture mechanics.

Dealing with the dynamic fracture, Belytschko, Chen and Xu et al. (2003) proposed a tip element in which the crack opens linearly and developed a propagation criterion with loss of hyperbolic. Later on, a singular tip enrichment function is proposed for the elastodynamiccracks with explicit time integration scheme [Belytschko and Chen (2004)]. In order to deeply study the stability and energy conservation to get a more accurate result, Ré thoré,Gravouil and Combescure (2005a, 2005b) combined Space and Time XFEM (STX-FEM)to obtain a unified space-time discretization, and concluded that the STX-FEM is a suitable technique for dynamic fracture problems. On the other hand, a new lumping technique for mass matrix was proposed in order to be more suitable for dynamic problems by Menouillard, Ré thoré and Combescure et al. (2006); Menouillard, Ré thoré and Noes et al.(2008) and the robustness and stability of the approach has been proved.

As we noticed, the previous study is all based on the classical enrichment scheme, and a large number of additional degrees of freedom (DOFs) are required. In the meantime,various improved enrichment methods have been studied. Song, Areias and Belytschko(2006) has reinterpreted the conventional displacement field, described discontinuities by using phantom nodes and superimposed extra elements onto the intrinsic grid for dynamic fracture problems. The method doesn’t require subdomain integration for the discontinuous integrand and has a highly efficient but nevertheless quite accurate formulation. Further,Duan, Song and Menouillard et al. (2009) has shown its practicability on the shell problem as well as three-dimension problem [Song and Belytschko (2009)]. Besides, changing the basic enrichment function is another solution. Menouillard, Song and Duan et al. (2010)proposed a new enrichment method with only a singular enrichment, which shows great accuracy for stationary cracks. The similar research has been done by Rabczuk, Zi and Geretenberger et al. (2008). Without crack tip enrichment, the Heaviside function has also been improved. Nistor, Pantale and Caperaa (2008) used only Heaviside function to model the dynamic crack growth. Kumar, Singh and Mishra et al. (2015) presented a new approach based Heaviside function along with a ramp function which contains information like crack length and angle. A similar method was proposed by Wen and Tian (2016); Tian and Wen (2016), which is based on an extra-dof-free partition of unity enrichment technique, and no more extra DOFs are added in the dynamic crack growth simulation.

For all the study discussed above, the stability of the method is always concerned.Generally, using explicit scheme for dynamic problem, one goes through a very small time stepping that leads to high computation cost, while with a larger one the numerical result may be divergent. So, in the present paper, we will focus on the stability of the numerical scheme. A new enrichment scheme is used for the elements influenced by crack tip based on the analytical solution of the displacement fields near crack tip. The Newmark scheme is adopted for time integration, and different parameters are tested to investigate their influence on the stability. In addition, DSIF is calculated as an important parameter which represents the variation of the stress field around the crack tip, and also can determine the stability of the simulation.

This paper is organized as following: Section 2 illustrates the governing equations and the XFEM; The explicit time algorithm and the lumping technique are introduced in Section 3; The DSIF is shown in Section 4; In Section 5, in order to verify the feasibility and accuracy of the simulation, several numerical examples are provided.

2 Governing equations and XFEM formulation

2.1 Governing equations

In order to develop the equations governing the problem, a homogeneous two-dimensional domainwith cracks is considered. As described in Fig.1, the domainis bounded bywhich is composed ofThroughout this paper, prescribed displacements are imposedand assumed to be traction-free. The strong form of the linear momentum equation and the boundary conditions are

2.2 The XFEM formulation

Consider a typical finite element mesh with four-node elements as shown in Fig.2, in which the geometry of crack is independent of the mesh. As in the classical XFEM, the nodes by Heavisde enrichment are enriched with two addition DOFs, and the shape functions are constructed from the Heaviside function）.）is defined as a unit magnitude for the elements cut completely by crack, and takes ±1 on the two sides of the crack. For the nodes with crack tip enrichment, they are enriched frequently by eight addition DOFs. The basic enrichment functions are inspired from the near tip displacement fields of mode I andmode II cracks in FEM, and can be written as the functions [Belytschko and Black (1999)].

Figure 1: Domain with cracks and prescribed boundary conditions

Figure 2: Typical discretization of a domain with crack and enriched nodes by XFEM

As mentioned above, the enrichment functions are developed based the asymptotic displacement fields of the crack tip, and can take different forms. In the present paper, a new form of enrichment functions is used, which derives from the asymptotic displacement fields directly. By shifting the enrichment functions, we are able to correspond the enriched nodes’ displacement to the true displacement with XFEM. The displacement can be written approximately as

Without concerning the damping effect, we substitute the displacement field Eq. (4) into the weak form Eq. (2). It then yields a system of linear algebraic matrix equations, which can be expressed as

where M is the mass matrix, K is the stiffness matrix and f is the force matrix:

sub-matrices and vectors that come in the foregoing equations are defined as below for four-node element ():

3 The explicit time integration

3.1 Time integration

As the most commonly used for dynamic problems, the Newmark scheme is chosen as the time integration algorithm. As we know, the time integration algorithm can be divided into two types: the explicit and the implicit. With the implicit method, there is no intrinsic limit to the time step. But we need to solve the global equations by iterating in each step. For dynamic problems, lots of iterations are needed, which has many disadvantages such as vast computation and low efficiency by implicit scheme. Compared with the implicit scheme, the explicit scheme solves the equations independently with no iterative, which is chosen in this paper.

The derived equation is given as

For a numeric scheme, the stability, consistency and convergence are the main reference standard. As instability is a sufficient condition for non-convergence, the Newmark scheme and their stability are discussed in this paper. The stability conditions of Newmark scheme are deduced in detail by Réthoré, Gravouil and Combescure (2004) with their custom notations. They can concluded as

Furthermore, due to such a restriction of stability condition, there must be a critical time step. With the time stepbeyond the critical value, the numerical instability and convergence problem will happen at some point. In contrast, the numerical results are very stable within the critical time.

In this paper, we will focus on figuring out the critical time step, and finding out the factors that can affect it. Thus, tests with different grid densities and different parameters in the Newmark scheme will be conducted.

3.2 The lumped mass

The matrix above in Eq. (8) known as the consistent mass matrix, includes standard terms,block-diagonal enriched terms, and coupling terms [Menouillard, Ré thoré and Combescure et al. (2006)]. However, for the problem of dynamic, the mass matrix lumped is used more frequently in order to simplify the numerical calculations. Due to the existence of additional DOFs which have no clear physical significance, the distribution of mass is not just as a simple average as in traditional FEM. Menouillard, Ré thoré and Noes et al. (2008)had in-depth study of the lumping technique for the mass matrix based on the conservation of mass and momentum, and proved its effectiveness with explicit scheme for dynamics by XFEM. Besides, the lumping technique was also researched by Zi, Chen and Xu et al.(2005); Elguedj, Gravouil and Maigre (2009); Song and Belytschko (2009); Jim, Zhang and Fang et al. (2016). In this paper, the lumped mass proposed by Menouillard, Ré thoré and Combescure et al. (2006) is used

4 DSIF

As the relevant quantities of crack tip may be questionable on accuracy such as stress fields,the SIF based on energetic consideration is used as a parameter of the strength of singularity.There are a few schemes to calculate the SIF, such as the displacement extrapolation method, the virtual crack extension method, the virtual crack closure method and the interaction integral method. The interaction integral method is used here which has the highest accuracy according to the research of Nagashima, Omoto and Tani (2003). In the interaction integral method, the auxiliary fields are introduced and superimposed onto the actual fields.

For dynamic loading case, an item related to inertia is added, and the interaction integral with force-free on crack surface can be given as

The basic algorithm used here for the DSIF is concluded as following:

(1) Give an integral rangR, then search for all the integral elements;

(2) Loop through all the integral elements;

(3) Loop through all the Gauss points in each integral element;

(4) Calculate the actual state and the auxiliary state of each Gauss point;

(5) Get the value of DSIF through Eq. (15) and Eq. (16).

whereRis the ratio between the actual integral radiusrand the minimum sizeLminof all elements as shown in Fig.3.

Figure 3: The integral elements for DSIF

5 Numerical examples

5.1 Stationary mode I crack

First, let us consider the problem of an infinite plate contains a semi-infinite crack whose geometry is shown in Fig.4. A theoretical solution of the problem was obtained by Freund(1990). To model this configuration, a rectangular plate of sizeL2H=10m4m with an initial edge crack of lengtha=5m under uniaxial tensile stress was used. The tensile stress was a typeof Heaviside step loading, and=500MPa. The material properties Young’s modulus =210Gpa, Poisson’s ratio=0.3 and the density=8000 kg/m3. A mesh of 3999 uniform square was used for tests. The theoretical DSIF of the problem with a stationary crack was given by Freund (1990):

Fig.6 presents the results of DSIF with different time steppingwhile=5. It shows good consistency and the results are not sensitive to the time step. So, this inspires us to improve the computational efficient with a larger step time which is less than the critical time. With a much larger time stepping, =20, the numerical result is rapidly divergent. As a consequence, there is a critical time stepping, which we will discuss it shortly.

Figure 4: The geometry and loading of a homogeneous material plate with crack

Figure 5: The DSIF with different integral path

Figure 6: The DSIF with different time step

5.2 Finite size edge plate with an arbitrarily oriented central crack

Secondly, the central cracks ofdifferent inclined angle are considered. The length of crack is the same, and the angles, =15°, 30°, 45°, 60°, 75° are examined. The problem has been discussed by Phan, Gray and Salvadori (2010) with Symmetric-Galerkin Boundary Element Method and by Liu, Bui and Zhang et al. (2012) with Smoothed Finite Element Method. The results are shown and compared with Phan, Gray and Salvadori (2010) in Fig.10. As depicted in Fig.10a, for the case of mode I, the values in the peak of DSIF curvesfor a small period of time after the stress wave arrive in the tip. Fig.10b reveals that DSIF in mode II are practically the same for the pair of=15°and=75°, and the pair of=30° and=60°. At the meanwhile, the curve of=45°has the highest peak value.

Figure 7: The rectangular plate with crack of different angle

Figure 8: The DSIF with different integral path of the left crack tip

Figure 9: The DSIF of the left crack tip with different time increments

Figure 10: The DSIF of crack tip with different rotation angle: (a) Mode I; (b) Mode II

5.3 The stable explicit time stepping analysis

This part focuses on the main factors that influence the critical time step. The grid density and iteration form are the two main subjects. The experiment configuration model is presented in Fig.7 with =0°. The material’ properties and the other parameters are the same as that used in last example. In order to get the critical time, the method of numerical approximation is used.

Firstly, the results were obtained with different grid densities. Three uniform meshes are considered, which are of CCT: 4999, CCT1: 2449, CCT2: 1324 elements. WithNewmark scheme, the critical time step of different meshesin Fig.11(a), the critical time we got is about=4.825×10-8s with 4999 elements. When the time stepcalculation results are completely consistent and do not produce divergence. Conversely,divergence is presented in the calculation when. The divergence occur at about 4.750, 7.154, 11.495, 15.141whenis 5.000×10-8s, 4.900×10-8s, 4.850×10-8s, 4.838×10-8s, respectively. As a comparison, Fig.11(b) is presented with the mesh of 2449. It is seen that the critical time is 10.025×10-8s which is improved than the one in Fig.11(a). The divergence occurs at about 4.095, 9.494, 12.090, 17.085whenis 10.500×10-8s, 10.100×10-8s, 10.075×10-8s, 10.050×10-8s, respectively.

Figure 11: Numerical stability with different time stepping (R=5,=0=1/2): (a) CCT:4999, (b) CCT1: 2449.

Table 1: The critical time stepping for different densities of grid (R=5=0,=1/2)

Table 1: The critical time stepping for different densities of grid (R=5=0,=1/2)

?

To clarify this case further, we repeated the above steps with CCT2: 1324 elements, and the comparison results are shown in the Table 1. The critical time step is about 16.568×10-8s in the case of CCT2, which is larger than the case of CCT1. It is hence concluded that the critical time step decreases with the increase of grid density. Besides, the critical time step of thestandard FEM for the lumped mass is also listed. With more elements, the critical time step is decreased, and this is consistent with the case of. The values ofare similar, which range from 88.064% to 91.344%. As Menouillard, Réthoré and Combescure et al. (2006); Elguedj, Gravouil and Maigre (2009) suggested,for the stationary crack, the value 2/3 is within the numerical range listed in this paper. So, the numerical stability can be guaranteed.

In addition, we took into account the effect of Newmark scheme for the critical time step.Four cases are concerned. Before studying the impact of iterative format on critical time step, all the cases are listed under the same conditions: CCT1, a mesh of 24×49 elements,R=5,=5×10-8s. We listed the first 30 microseconds with different parameter valuesin Fig.12. An approximately identical result can be obtained. The stability conditions of the Newmark scheme are also verified directly.

In Fig.11(b), we presented the test result with=1/2. As a comparison, the result with=2/3 is shown in Fig.13. The divergence occur at about 6.030, 12.848, 75.168whenis 9.000×10-8s, 8.800×10-8s, 8.700×10-8s, respectively. The critical time we obtained is about=8.685×10-8s, which is smaller than the case of=1/2. For further investigation, the cases of=3/4,=1 are tested. The results are listed in Table.2. The critical time step of the standard FEM for the lumped mass are also listed. In Table.2, it is seen that the critical time step decreases with the increase ofFurthermore, we observed that the values ofand the different parameter

Figure: 12 Numerical results with different parameters(CCT1: 2449, R=5, =0,=510-8s)

Figure: 13 Numerical stability with different time stepping (CCT1: 2449, R=5, =0,=2/3)

Table 2: The critical time stepping for different parameters (CCT1: 2449, R=5=0)

Table 2: The critical time stepping for different parameters (CCT1: 2449, R=5=0)

?

6 Conclusions

In the present paper, we carried out some numerical experiments of the stable explicit time stepping within the XFEM framework. A new enrichment scheme for crack tip is proposed and its applicability and availability has been sufficiently verified. The DSIF is used as an important parameter of the dynamic response and is also a parameter of judging the stability of numerical method. Objective to studying the factors that can affect the stability,different densities of grid and different parameters of Newmark scheme have been tested.The conclusions are shown as:

· The grid density and the form of iterative method have obvious effects on stability;

· The critical time steppingdecreases with the increase of grid density;

· The critical time steppingdecreases with the increase of the parameterbetween 0.5 and 1of Newmark scheme;

Furthermore, the simulation results are found in good agreement with each other when they are stable. Therefore, increasing time stepping appropriately in the range of critical value can improve the computational efficiency.

Acknowledgment:The authors are grateful to the National Natural Science Foundation of China(No.11672101, No.11372099), the 12th Five-Year Supporting Plan Issue (No.2015BAB07B10),Jiangsu Province Natural Science Fund Project (No. BK20151493) and the Postgraduate Research and Innovation Projects in Jiangsu Province (No.2014B31614) for the financial support.

Belytschko, T.; Black, T.(1999): Elastic crack growth in finite elements with minimal remeshing.International Journal for Numerical Methods in Engineering,vol. 45, no. 5,pp. 602-620.

Belytschko, T.; Chen, H.(2004): Singular Enrichment Finite Element Method for Elastodynamic Crack Propagation.International Journal of Computational Methods, vol.1, no. 1, pp. 1-15.

Belytschko, T.; Chen, H.; Xu, J.; Zi, G.(2003): Dynamic crack propagation based on loss of hyperbolicity and a new discontinuous enrichment.International Journal for Numerical Methods in Engineering,vol. 58, no. 12, pp. 1873-1905.

Belytschko, T.; Moes, N.; Usui, S.; Parimi, C.(2001): Arbitrary discontinuities infinite elements,International Journal for Numerical Methods in Engineering, vol. 50, no. 4, pp.993-1013.

Dolbow, J.; Moes, N.; Belytschko, T.(2000): Discontinuous enrichment in finite elements with a partition of unity method.Finite Elements in Analysis and Design,vol. 365, no. 3-4, pp. 235-260.

Duan, Q.; Song, J. H.; Menouillard, T.(2009): Belytschko. Element-local level set method for three-dimensional dynamic crack growth.International Journal for Numerical Methods in Engineering, vol. 80, no. 12, pp. 1520-1543.

Elguedj, T.; Gravouil, A.; Maigre, H.(2009): An explicit dynamic extended finite element method. Part 1: Mass lumping for arbitrary enrichment functions.Computer Methods in Applied Mechanics and Engineering, vol. 198 no. 30-32, pp. 2297-2317.

Freund, L. B.(1990): Dymamic fracture mechanics.Cambridge Monographs on Mechanics and Applied Mathematics.

Jim, L.; Zhang, T.; Fang, E.; Song, J. H.(2016): Explicit phantom paired shell element approach for crack branching and impact damage prediction of aluminum structures.International Journal of Impact Engineering,vol. 87, pp. 28-43.

Kumar, S.; Singh, I. V.; Mishra, B. K.; Singh, A.(2015): New enrichments in XFEM to model dynamic crack response of 2-D elastic solids,International Journal of Impact Engineering, vol. 87, pp. 198-211.

Liu, P.; Bui, T.; Zhang, C.; Yu, T.T.; Liu, G. R.; Golub, M. V.(2012): The singular edgebased smoothed finite element method for stationary dynamic crack problems in 2D elastic solids.Computer Methods in Applied Mechanics & Engineering,vol. 233-236, no. 4, pp.68-80.

Moes, N.; Dolbow, J.; Belytschko, T.(1999): A finite element method for crack growth without remeshing.International Journal for Numerical Methods in Engineering, vol. 46,no. 1, pp. 131-150.

Menouillard, T.; Réthoré, J.; Combescure, A. et al.(2006): Efficient explicit time stepping for the eXtended Finite Element Method(X-FEM).International Journal for Numerical Methods in Engineering, vol. 68, no. 9, pp. 911-939.

Menouillard, T.; Réthoré, J.; Moes, N. et al.(2008): Mass lumping strategies for X-FEM explicit dynamics: Application to crack propagation.International Journal for Numerical Methods in Engineering, vol. 74, no. 3, pp. 447-474.

Menouillard, T.; Song, J. H.; Duan, Q. et al.(2010): Time dependent crack tip enrichment for dynamic crack propagation.International Journal of Fracture, vol. 62, no.1-2, pp.33-69.

Nagashima, T.; Omoto, Y.; Tani, S.(2003): Stress intensity factor analysis of interface cracks using X-FEM,International Journal for Numerical Methods in Engineering,vol.56, no. 8, pp. 1151-1173.

Nistor, I.; Pantale, O.; Caperaa, S.(2008): Numerical implementation of the extended finite element method for dynamic crack analysis.Advances in Engineer Software,vol. 39,no. 7, pp. 573-587

Phan, A. V.; Gray, L. J.; Salvadori, A.(2010): Transient analysis of the dynamic stress intensity factors using SGBEM for frequency-domain elastodynamics.Computer Methods in Applied Mechanics & Engineering,vol. 199, no. 45-28, pp. 3039-3050.

Rabczuk, T.; Zi, G.; Gerstenberger, A. et al.(2008): A new crack tip element for the phantom-node method with arbitrary cohesive cracks,International Journal for Numerical Methods in Engineering, vol. 75, no. 5, pp. 577-599.

Réthoré, J.; Gravouil, A.; Combescure, A.(2004): A stable numerical scheme for the finite element simulation of dynamic crack propagation with remeshing.Computer Methods in Applied Mechanics & Engineering, vol. 193, no. 42-44, pp. 4493-4510.

Réthoré, J.; Gravouil, A.; Combescure, A.(2005a): An energy-conserving scheme for dynamic crack growth using the eXtended finite element method,International Journal for Numerical Methods in Engineering, vol. 63, no. 5, pp. 631-659.

Réthoré, J.; Gravouil, A.; Combescure, A.(2005b): A combined space-time extended Finite Element Method.International Journal for Numerical Methods in Engineering, vol.64, no. 2, pp. 260-284.

Song, J. H.; Belytschko, T.(2009): Dynamic fracture of shells subjected to impulsive loads.Journal of Applied Mechanics, vol. 76, no. 5, pp. 051301-1-9.

Song, J. H.; Belytschko, T.(2009): Cracking ong node method for dynamic fracture with finite elements,International Journal for Numerical Methods in Engineering,vol. 77, no.3, pp. 360-385.

Song, J. H.; Areias, P. M. A.; Belytschko, T.(2006): A method for dynamic crack and shear band propagation with nodes phantom nodes.International Journal for Numerical Methods in Engineering, vol. 67, no. 6, pp. 868-893.

Tian, R.; Wen, L.(2016): Improved XFEM: An extra-dof free,well-conditioning, and interpolating XFEM.Computer Methods in Applied Mechanics & Engineering,vol. 285,no. 3, pp. 639-658.

Wen, L.; Tian, R.(2016): Improved XFEM: Accurate and robust dynamic crack growth simulation.Computer Methods in Applied Mechanics & Engineering,vol. 308, pp. 256-285.

Zi, G.; Chen, H.; Xu, J.; Belytschko, T.(2005): The extended finite element method for dynamic fractures.Shock and Vibration,vol. 12, no. 1, pp. 9-23.