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

Xue-cong Liu, Qing Zhang,*, 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 order to model the behavior of discontinuities, the way of re-meshing is used which is to align the element edges with the discontinuities by classic finite element method (FEM). However, for the case of crack arbitrary growth, the mesh changes at every step and incurs high computation cost. Besides, other solutions such as meshfree method, boundary element method and extended finite element method (XFEM) are available, and XFEM is considered to be the most applicable.

As proposed by Belytschko and Black (1999); Moes, Dolbow and Belytschko (1999);Belytschko and Moes (2001), XFEM becomes a dominant numerical scheme nowadays.XFEM is based on the concept of partition unity, and the crack can be modeled independent of finite element mesh. All the elements are divided into the normal partsand the enriched parts. For the enriched parts, the elements can be influenced directly by crack, so the enrich functions are introduced. As the most used enrich functions,Heaviside function and the Westergaard stress function are used frequently for the cracks and cracks’ tips, respectively.

Dealing with the dynamic fracture, Belytschko, Chen, Xu et al (2003) proposed a kind of tip element in which the crack opens linearly and developed a propagation criterion with loss of hyperbolic. Later on, a singular enrichment function for tips is proposed for the elastodynamic cracks with explicit time integration scheme [Belytschko and Chen(2004)]. In order to deeply study the stability and energy conservation and 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é, Combescure et al (2006); Menouillard, Ré thoré,Noes et al (2008) and the robustness and the stability of this approach has been proved.

As we noticed, the previous studies are 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, 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,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, 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, 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. Based on the analytical solution of the displacement fields near crack tip, a new enrichment scheme is used for the elements influenced by 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 formulation; The explicit time algorithm and the lumping technique are introduced in Section 3; The DSIF is shown in Section 4; In Section 5, several numerical examples are provided, the feasibility and stability of the simulation are discussed.

2 Governing equations and XFEM formulation

2.1 Governing equations

2.2 The XFEM formulation

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 displacement fields near the crack tip, and can have different forms. In the present paper, a new form of enrichment functions is used, which derives from the 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

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

The 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 it needs to solve the global equations by iterating in each step.For dynamic problems, lots of iterative calculations are needed, using the implicit scheme can bring many disadvantages such as vast computation and low efficiency. Compared with the implicit scheme, the explicit scheme which is chosen in the present paper can solve the equations independently with no iterative.

For explicit scheme, two parametersandare considered. Combine the Newmark scheme and Eq.(6), the derived equation can be given as

For a numeric scheme, the stability, consistency and convergence are the main reference standards. From now on, we will focus on the stability of Newmark scheme because of the instability is a sufficient condition for non-convergence. The stability conditions of Newmark scheme are deduced in detail by Ré thoré, Gravouil and Combescure (2004)with their custom notations:

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 is very stable within the critical time. We will focus on figuring out the critical time step, and finding out the factors that can affect it in this paper. 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) is known as the consistent mass matrix, which includes standard terms, block-diagonal enriched terms, and coupling terms [Menouillard, Ré thoré,Combescure et al (2006)]. However, for the problem of dynamic, the lumped mass matrix 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é, 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,Xu et al (2005); Elguedj, Gravouil and Maigre (2009); Song and Belytschko (2009); Jim,Zhang, Fang et al (2016). In this paper, the lumped mass proposed by Menouillard,Ré thoré, Combescure et al (2006) is used

where Ω is the element being considered, m is the element’s mass， mes(Ω)is the area of element in 2D, nnodeis the number of nodes in element, and H is the Heaviside function.

4 DSIF

In this section, the DSIF is illustrated. Based on energetic consideration, the SIF is used as a parameter of the strength of singularity, and some questionable relevant quantities of crack tip such as stress fields are avoided. 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. Due to the research of Nagashima, Omoto and Tani (2003), the interaction integral method has the highest accuracy and is used here. 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 rang R, 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).

where R is the ratio between the actual integral radius r and the minimum size Lminof all elements as shown in Fig.3.

Figure 3: The integral elements for DSIF

5 Numerical examples

5.1 Stationary mode I crack

First, in order to illustrate the effectiveness of the new enrichment scheme presented above, we study the problem of an infinite plate contains a semi-infinite crack as shown in Fig.4. A prescribed vertical loading of=500MPa is applied to the upper surface.The evolution of the loading is a type of Heaviside step function. For the geometry of configuration, a rectangular plate of size L2H=10m4m with an initial edge crack of length a=5m is used. Some material parameters are: Young’s modulus=210Gpa,Poisson’s ratio=0.3, the densityused for tests. A theoretical solution of this problem with a stationary crack was obtained by [Freund (1990)] and it is given by

Fig.5 presents the values of DSIF with different integral domains. The time step is chosen as=0.1. The DSIF is normalized byand the numerical time is normalized by. At the beginning, time, due to the stress wave has not reach the integral domain, the value of DSIF is 0. Then, after stress wave reaches the crack tip, the results with different integral domains are in good agreement with each other as shown in Fig.5. The theoretical solution for comparison has also been plotted, and shows good accuracy to the present result.

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. The results inspire us to improve the computational efficient with a larger step time which is less than the critical time. However, with a much larger time stepping for test,=20, the numerical result is rapidly divergent. As a consequence, it must have a critical time stepping, which we will discuss it later.

In addition, the result of DSIF with four crack tip enrichment functions of standard XFEM is also plotted for comparison in Fig.6. It can be seen that the new enrichment scheme offers the almost same accuracy as XFEM with standard enrichment functions.

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.2Finite size edge plate with an arbitrarily oriented central crack

In this example, A plate with a central crack under uniaxial tensions is considered. As shown in Fig.7, the dimensions of the plate is 2h=0.04 m and 2b=0.02 m, and the length of the central crack is 2a=0.0048m. The material’s properties are:E=199.99Gpa,=0.3,=5000kg/m3. A prescribed vertical loading of=100MPa is applied to the upper surface and the lower surface. A mesh of 4999 uniform elements is used.

First of all, a horizon central crack=0° is considered. For the cases of different integral domains, Fig.8 shows that the rangeRhas little effect on the crack tip’s DSIF. The results agree very well with the conclusion that the SIF under different integral path are the same.The results are also compared with the standard XFEM with four crack tip enrichment functions, and shows good accuracy to the present result.

Fig.9 presents the numerical results with different time step, the same conclusion can be drawn as in Fig.6. In additional, the numerical result given by Song, Areias and Belytschko (2006) is compared, and shows the similar accuracy.

Secondly, the central cracks of different inclined angles are considered. The length of crack is the same 2a=0.0048 m, and the angles,=15°, 30°, 45°, 60°, 75° are tested. The problem has been discussed by Phan, Gray and Salvadori (2010) with Symmetric-Galerkin Boundary Element Method and by Liu, Bui, Zhang, et al (2012) with Smoothed Finite Element Method. For the case of mode I, as depicted in Fig.10a, the values in the peak of DSIF curves decrease by the increase offor 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. The results are compared with Phan, Gray and Salvadori (2010) and shown the good accuracy.

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 considered here. The experiment configuration model is presented in Fig.7 with=0°. The material’ properties and the other parameters are the same as that usedin last example. The method of numerical approximation is used to get the critical time .

Firstly, the results with different grid densities were obtained. Three uniform meshes are considered, which are of CCT: 4999, CCT1: 2449, CCT2: 1324 elements. Withscheme, the critical time step of different meshes can be turned out. As shown in Fig.11(a), the critical time we got is about=4.825×10-8s with 4999 elements. When the time stepis less than, the numerical calculation 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)

CCT:49×99 CCT1:24×49 CCT2:13×24 ( ) 4.825×10-2 10.025×10-2 16.568×10-2 ( ) 5.479×10-2 11.012×10-2 18.138×10-2 88.064% 91.037% 91.344%

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 the standard FEM for the lumped mass is also listed. With moreis decreased, and this is consistent with the case ofare similar, which range from 88.064% to 91.344%. As Menouillard, Ré thoré, Combescure et al (2006); Elguedj, Gravouil and Maigre (2009)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 of. So does.Furthermore, we observed that the values ofare nearly the same (about 91%),and the parameterhave nearly no influence on the values of.

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 parameter

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. 2015 BAB07B10), Jiangsu Province Natural Science Fund Project (No. BK 20151493) and the Postgraduate Research and Innovation Projects in Jiangsu Province (No.2014B 31614) 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-5.

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 dynamics 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 edge-based 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.

Menouillard, T.; Ré thoré, J.; Combescure, A.; Bung, H.(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.; Bung, H.(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.;Belytschko, T.(2010): Time dependent crack tip enrichment for dynamic crack propagation.International Journal of Fracture, vol. 62,no. 1-2, pp. 33-69.

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.

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.; Wall, W.A.(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.