A combined space–time extended finite element method

The Newmark method for the numerical integration of second order equations has been extensively used and studied along the past fifty years for structural dynamics and various fields of mechanical engineering. Easy implementation and nice properties of this method and its derivatives for linear problems are appreciated but the main drawback is the treatment of discontinuities. Zienkiewicz proposed an approach using finite element concept in time, which allows a new look at the Newmark method. The idea of this paper is to propose, thanks to this approach, the use of a time partition of the unity method denoted Time Extended Finite Element Method (TX‐FEM) for improved numerical simulations of time discontinuities. An enriched basis of shape functions in time is used to capture with a good accuracy the non‐polynomial part of the solution. This formulation allows a suitable form of the time‐stepping formulae to study stability and energy conservation. The case of an enrichment with the Heaviside function is developed and can be seen as an alternative approach to time discontinuous Galerkin method (T‐DGM), stability and accuracy properties of which can be derived from those of the TX‐FEM. Then Space and Time X‐FEM (STX‐FEM) are combined to obtain a unified space–time discretization. This combined STX‐FEM appears to be a suitable technique for space–time discontinuous problems like dynamic crack propagation or other applications involving moving discontinuities. Copyright © 2005 John Wiley & Sons, Ltd.


INTRODUCTION
For the past fifty years, the Newmark method (see Reference [1]) has been considered as a reference for the numerical integration of second order equations in structural dynamics or various fields of mechanical engineering. The easy implementation and the nice properties concerning numerical damping and energy conservation of this method and its derivative (Fox and Goodwin, HHT and others) are appreciated. Zienkiewicz [2] and Wood [3] propose a finite element (FEM) view of the Newmark method and extend this concept in a unified set of single step algorithms (see Reference [4]). During the same period, space FEM was developed in various ways. Such space or Time FEM is quite easy to use and has a lot of advantages, but the main drawback concerns the treatment of discontinuities (and even more if those discontinuities are moving, e.g. cracks, material interface, etc.). For space discontinuous problems, methods have been developed, and particularly the Partition of the Unity method (see Reference [5]), which have inspired several extensions of classical FEM for treating discontinuities (see References [6][7][8]).
This paper proposes to use an extended finite element basis for time interpolation: enriched shape functions are added to a classical Time FEM to take into account time discontinuities. The formulation we propose is velocity based and uses discontinuous enrichment at each time level. The resulting method presents nice stability and accuracy properties as well as capabilities in the treatment of time discontinuities. In a particular case, this approach allows a new look at the time-discontinuous Galerkin methods (T-DGM) proposed in References [9,10]. The proposed combined space-time extended finite element method (STX-FEM) gives a unified space time discretization and accurate results when space and/or time discontinuities have to be simulated. In Reference [11], Chessa and Belytschko proposed a different approach to treat discontinuities in space-time FE. They consider a space discontinuity that is moving in space-time FE using an enriched approximation in space. Here, the interest is on the modelling of space and time discontinuities and the approximation will be enriched in space as well as in time.
This paper is organized as follows: Section 2 summarizes the original time FEM (T-FEM) and describes the formulation we use. Section 3 presents a time extended FEM (TX-FEM) framework and a particular case: the Heaviside step function is used to model time discontinuities. Stability and accuracy are studied for this particular case. Then the combination of space and time X-FEM is developed and several examples illustrate the capabilities of the proposed method (Section 4).

TIME FINITE ELEMENT METHOD
In this section, the numerical integration of the elastodynamics equation (Equation (1)) is studied. For easier understanding, the problem is written for a single degree of freedom massspring system (parameters k, m) submitted to an external loading f . Here x,ẋ,ẍ denote the displacement, the velocity and the acceleration, respectively. The original Newmark method, and Zienkiewicz [2] and Wood [3] approaches are summarized. An alternative formulation is presented.

Reference time-stepping formula: the Newmark method [1]
In the Newmark method, no assumption is made on the type of function that can describe displacement or velocity. The assumption is on the continuity and the approximation of the rest in the Taylor formulae. This approximation uses two parameters, N for x and N forẋ: The elastodynamics equation (Equation (1)) at time level n is Considering Equation (4) at time levels n − 1, n and n + 1 and Equations (2) and (3) at time levels n and n + 1, the problem can be written in terms of x. The resulting two-steps method is as follows: given x n−1 and x n , find x n+1 such that This formula is interesting because it will allow us to establish relations between all presented methods in the following paragraphs. Various combinations of Equations (2)-(4) can give other time-stepping formulae.

Literature survey in T-FEM
Zienkiewicz and then Wood presented other approaches using the weighted residual method and finite element concepts. [2]. In Reference [2], the displacement is interpolated in time as follows:

The Zienkiewicz approach
where L i are quadratic shape functions. Then, substituting Equation (6) into the weighted residual form of the elastodynamics equation (Equation (1)) integrated from − t to t (W being any appropriate weight function), we have to obtain Equation (8) is formally the same as Equation (5), but Z and Z are obtained from Equation (9). Consequently Z , Z are consistent with the polynomial interpolation and with the choice of the weight function W . It is not the case for N , N , which can be chosen arbitrarily in Equations (2) and (3). As a consequence, we can say that an arbitrary choice of N , N may not be consistent with polynomial interpolation. This two-steps method shows the analogy between the Newmark method and T-FEM but does not allow us to change the time-step size easily. [3]. Wood chooses to write the displacement with a quadratic interpolation in terms of three parameters x n , x n+1 andẋ n as follows:

The Wood approach
This expression is then introduced into the weighted residual form of the elastodynamics equation and integrated from 0 to t: where F is the right-hand side of Equation (8) and W , W are given by and thenẋ n+1 is given by differentiating Equation (10): This method is equivalent to the Zienkiewicz approach and also to the Newmark method, but once again only for appropriate N , N . As soon as it is formulated as a one-step method it gives more flexibility for adaptive time-step size. Zienkienwicz et al. propose a generalization of this approach for p-degree polynomial approximation in Reference [4].

One-step velocity-based formulation
The formulation we are presenting is velocity based. For a better understanding we change the notation: u denotes the displacement, v the velocity and a the acceleration. Let us choose linear shape functions to discretize the velocity in [t n ; t n+1 ]: where The displacement is then obtained by the integration of Equation (14) considering the initial condition u(t n ) = u n : The problem is then expressed with the weighted residual form of Equation (1) integrated from t n to t n+1 ( t = t n+1 − t n ). As this formulation is velocity based the problem becomes: given v n and u n find v n+1 such that This equation can also be obtained from the Newmark method in the particular case 2 N = N .
Then we obtain u n+1 using Equation (16): This was obtained in the Wood approach (see Equation (13)) and can result from the Newmark approach combining Equations (2) and (3) assuming that 2 N = N . This one-step method is equivalent to the Wood approach, but, as it is velocity based using linear approximation, external loading is interpolated linearly (whereas quadratic shape functions are used in the Wood approach). This approach will be justified in the next section.

Remarks
(i) As the velocity is piecewise linear, the acceleration is piecewise constant. Inside the time-slab I n+1 =]t n ; t n+1 [, using the differentiation of Equation (14), we have Consequently, the acceleration is not unique at the time levels in this approach. In the Newmark method, the assumptions are on the continuity of the displacement that enables to define the acceleration. Using Taylor formulae, we may here expand u and v in an excessive way, definingü(t n ) =ü n andv(t n ) =v n :  (22) To be consistent with the Newmark method, one could defineü n andv n bÿ Those definitions are kinematically consistent ifü n =v n , i.e. if 2 N = N . Indeed, in the Newmark method, kinematic relations between displacement velocity and acceleration are not exactly satisfied for arbitrary N , N if u is quadratic and v linear in time.
(ii) Such one-step methods present the following advantages: an easy implementation and adaptive time-step size can be used. The simplicity of the method is due to the polynomial interpolation of the displacement and the velocity. Like the Newmark method, these methods are consequently not well suited for the treatment of time discontinuities.

General formulation
In this section, we will use the partition of the unity method (PUM) (see Reference [5]). This is motivated by the fact that polynomial interpolation cannot well approximate time discontinuities.
As the set of functions { i } i=0..N is a partition of the unity in the interval [0; T ], we can enrich this basis of linear shape function using the PUM: where v c i are classical degrees of freedom, j are enriched shape functions, v e i, j additional degrees of freedom and N j the set of time levels that are enriched with j . The j can be chosen to capture the non-polynomial part of the solution with good accuracy.

TX-FEM formulation.
We now consider that the solution to be computed presents time discontinuities. Consequently, we choose to enrich the basis with the Heaviside step function H . We put one enriched function j = H (t − t j ) centred at each time level, and N j is reduced to the time level t j . Because j = H (t − t j ) is zero for all time t earlier than t j and j has a compact support [t j −1 ; t j +1 ], H (t − t n ) is the only enriched function that is active inside the time-slab I n+1 =]t n ; t n+1 [ (see Figure 1). Inside I n+1 , the velocity and the displacement are decomposed in a continuous and an enriched contribution: where and We abusively denote the additional degree of freedom v e n, n by v e n+1 because it has to be determined when the time slab I n+1 is solved.
The kinematic constraints between displacement and velocity (u = v) and the displacement continuity at the time levels (u(t + n ) = u(t − n )) are strongly enforced. Consequently, the continuous and enriched contributions of the displacement are We then write the weighted residual form of the elastodynamics equation in I n+1 and enforce the velocity continuity weakly at time level t n : As we added enriched degrees of freedom, two independent weight functions {W i } i=1, 2 must be chosen to solve the problem, which is formulated as follows: given u c n , u e n , v c n and v e n , find u c n+1 , u e n+1 , v c n+1 and v e n+1 such that This formulation has the following properties: the displacement is continuous with time, the velocity is linear inside the time slab and can have discontinuities at each time level (because the velocity continuity is only weakly enforced), kinematic constraintu = v is always satisfied. This time integrator is defined by six parameters 2 using Equations (31). The method is now reduced to find appropriate weight functions W 1 , W 2 that give the time integrator the best possible properties.
Concerning the acceleration, it can be defined as a piecewise-constant function and its value inside I n+1 is obtained by differentiation of Equation (28):

A new
and choose the weight functions to be The parameters are and the resulting system (Equations (36)) is the same as that obtained with a P3-P1 approximation in T-DGM (see for e.g. Reference [9]) or a velocity-based formulation from Reference [10]: The presented method has the same properties as the formulation proposed by Michler et al. [10]. Both mentioned formulations [9,10] are very close: the computed values of velocity and displacement at the time levels are the same and the displacement continuity is ensured. The difference concerns the displacement interpolation and the kinematic constraint: here the displacement approximation is P2 because the kinematic constraint (u = v) is strongly enforced, P3-P1 approximation allows high order approximation but the kinematic constraint is only weakly enforced.

Stability and accuracy.
All numerical results presented in this paragraph are obtained for a mass-spring single degree of freedom system without external load but with a nonzero initial velocity. Weight functions for the numerical results are W 1 = n , W 2 = n+1 . The results are consequently valid for this special TX-FEM and the above-mentioned T-DGM. Using Equations (30), the single degree of freedom problem without external loads can be recast in A = H −1 1 H 0 is the amplification matrix and the stability depends only upon the eigenvalues of A. They are obtained solving It can be shown that this equation is equivalent to (i) is less than 1.
(ii) The modulus of r i is strictly less than one if the order of multiplicity of r i is greater than one.
The spectral radius of the amplification matrix is compared to the spectral radius of the Newmark average acceleration method ( N = 1 2 , N = 1 4 ). Those spectral radii are plotted (see Figure 2) versus parameter t. The first observation is that the method is unconditionally stable for this set of weight functions (W 1 = n , W 2 = n+1 ). We can observe that the spectral radius is not equal to one and the method will show numerical damping. Two regions of t ( t ≈ and t 10) are particular. For these regions, we can observe bifurcations because with i the eigenvector corresponding to the real eigenvalue r i = e −( t/ i ) . A consequence of this property is that the numerical response for a multi-degree of freedom problem will not be disturbed by spurious numerical oscillations (in the corresponding region of t). This phenomenon will be illustrated by the examples of the next section.
As mentioned above, the method gives numerical damping. Figure 3 presents the evolution of the numerical damping ratio . As shown by Figures 4 and 5, the accuracy of the method  is not affected by this numerical damping: the relative period error is ten time less than that for the Newmark method and the method is third order accurate.

Discretized energy balance.
One of the interests of the TX-FEM is that the framework gives a way to study the properties of the resulting methods (e.g. T-DGM). For example, we can use the following notations for kinematic quantities: Then, using the definition of v + n and v − n+1 from Equation (33) and multiplying the last two equations obtained in Equation (30a) and (30d) by [u] we have The set of parameters for W 1 = n , W 2 = n+1 (Equation (35)) gives Considering a mass-spring system without external load we can have the following interpretation: in a simplified way, some energy is created at the time level (Equation (44b) and then dissipated during the time slab (Equation (44a Using the continuity of u and Equation (44b) Hence, Consequently, for this set of weight functions (W 1 = n , W 2 = n+1 ), looking at the discrete energy, the method is not energy preserving as illustrated in Figures 6 and 7 and as could be predicted by the numerical damping previously pointed out. The energy balance (normalized by T 0 the initial kinetic energy) of the discretized formulation is plotted in Figure 6. One can observe (in Figure 7 which is a zoom of Figure 6) the energy created at the time levels and the energy dissipated during the time slabs (results are presented for one period of the single degree of freedom mass-spring system with t ≈ 0.2). As the method dissipates more energy in average in the time slabs than it creates at the time levels, it is globally dissipative. Note that the dissipation of the discrete energy occurs during the time slabs whereas the conservation of the continuous energy is ensured by Equation (29).
Finally, whereas the method presents numerical damping and numerical energy dissipation, it is third order accurate and its advantages (particularly for time discontinuities integration) will be illustrated by the examples in the next section. As the proposed results are obtained with W 1 = n , W 2 = n+1 , they are also valid for the above-mentioned T-DGM. Such a choice of weight functions appears quite reasonable and gives interesting properties to the time integrator.

SPACE-TIME EXTENDED FINITE ELEMENT METHOD
Using a combined space-time X-FEM, the solution of the continuous problem U is first approximated in space: This stage is detailed in the following paragraph. Then a time approximation is made using the method proposed in Section 3 to write

Space discretization
In the method presented, we use the extended finite element method first introduced by Black [13]. In this method, an enrichment is added to the classical finite element approximation using the PUM developed by Babuska [5]. For the static problem, the displacement field can be written with enriched basis of shape functions as in Reference [14]: where N is the set of all nodes in the mesh, N cut the set of nodes whose support is completely cut by the crack and N branch the set of nodes that belong to elements partially cut by the crack. N i are the classical shape functions and H is a discontinuous function, whose value is 1 if x is above the crack surface and −1 if x is bellow.
[B ] are branch functions ((r, ) are local cylindrical crack tip co-ordinates): In this way the approximation in space is enriched by a local partition of unity with no particular treatment for the blending elements (see Reference [15]). In our approach, all fields (displacement, velocity and acceleration) are discretized with Equation (50). Consequently, the discrete enriched problem can be written as a classical dynamic problem in the scope of the finite element method: we define the mass matrix and the stiffness matrix in the usual manner.

Time discretization
If we use the approach presented in Section 3 we can write the space discretized velocity as follows:

Application to dynamic crack growth
New space shape functions are added to simulate the crack growth. As shown in Figure 8, new singular enrichments are added on the new set N n+1 branch . For discontinuous enrichment, new shape functions are only added on the set N new cut = N n+1 cut \N n cut . New degrees of freedom are initialized to zero. This strategy is stable and energy preserving (see Reference [16]) referring to theoretical studies of Newmark-type schemes for time-dependent discretization (see Reference [17]). This strategy allows the use of singular enriched functions ([B ]) in timedependent problem. It has been valitated by the authors in an earlier article [16]. Other references on dynamic crack propagation with X-FEM can be found in References [18,19]. The present strategy differs from Reference [18], where the stress intensity factors are considered to be additional degrees of freedom of the enriched approximation. The main advantage of the technique presented in Reference [18] is that the stress intensity factors are obtained directly without postprocessing of the mechanical fields. In this approach, the enriched zone must contain the crack path (to guarantee the numercial stability of the time integrator), which consequently has to be known in advance. Belytschko et al. [19] have proposed a new discontinuous enrichment. Combined with a cohesive zone model, they avoid the use of a singular enrichment and the difficulties encountered when it is used in time-dependent problems.

Brittle beam.
This first example has been chosen because of its simplicity and easy interpretation. An elastic beam is studied: it is 10 m long and its width is 0.5 m. The material is assumed to be homogeneous and isotropic with elastic behaviour (material parameters are E = 210 GPa and = 8000 kg m −3 ). It is submitted to a tensile stress = 500 MPa at one end and fixed at the other. To simulate the brittle fracture of the beam, we suddenly introduce an enriched discontinuous shape function. This discontinuity is located at the middle of the beam and occurs at time t/t c = 2.5 (t c is half the period of the vibration of the entire beam). Figure 9 compares the results obtained with the Newmark average acceleration method and the TX-FEM. The time-step size has been chosen to obtain the same accuracy on the displacement using both methods. For the Newmark method, it is two-third of the critical time-step size and twice with the TX-FEM.
The results are plotted for the free end of the beam and for the node just before the created interface of the fixed half of the broken beam. As a consequence of the choice of the time-step size, the displacement is the same in both methods. Looking at the velocity, the response of the Newmark method shows oscillations after time discontinuities. The results obtained with the TX-FEM are more accurate and not disturbed by numerical spurious oscillations. The difference between both methods increases at the interface because of the sudden introduction of the discontinuity (a time discontinuous space discontinuity). The capability of the presented method is also well illustrated by the computed acceleration that keeps a physical meaning after the 'rupture'. This example well illustrates the capability of the method in the treatment of time discontinuities. Indeed, the particular property of this method (illustrated by the spectral radius of the amplification matrix of the single degree of freedom problem) comes from the aperiodic response for typical spurious frequency that appears with the Newmark method. The frequency of the numerical oscillations with the Newmark method is about twice the chosen discretization frequency: i.e. spurious numerical vibration modes whose period T = 2 t is solicited. For such solutions, the corresponding value of the parameter t = 2 t/T is . Looking at the spectral radius of the TX-FEM, is one of the two particular regions mentioned in Section 3, where the response is aperiodic. Consequently, TX-FEM avoids spurious numerical high-frequency oscillations. This gives to the TX-FEM an attractive property for structural dynamics. The example has also been chosen because of its similarity with dynamic crack propagation, where space discontinuities appear suddenly in the numerical model.

Semi-infinite crack in an infinite plate submitted to a tensile stress wave.
The second example is the infinite plate with a semi-infinite crack, where the theoretical solution is known (see Reference [20]). Under these assumptions (infinite plate with a semi-infinite crack), for the Figure 10. Geometry and loading for the example of the infinite plate.
geometry described in Figure 10, the analytical solution is only valid for time t 3t c = 3H/c d (c d is the dilatational wave speed) when the reflected stress wave arrives on the crack tip. As the wave reaches the crack, mode I stress intensity factor can be written for a moving crack as For a moving crack we can write where k is a universal function of the crack tip velocityȧ, which can be approximated by And finally, we have To compute stress intensity factors, an interaction integral is used. It has been presented by authors in Reference [16]. This approach is based on the Lagrangian conservation law using auxiliary (denoted by a superscript aux) and virtual crack extension fields (q). It is valid for arbitrary 2D moving cracks: The interaction integral (see Equation (57)   Solutions are computed using a 40 × 80 quadrangle element with a linear approximation. The results are presented in Figure 11. As shown in Section 3, the convergence of the TX-FEM allows us to use an approximately four time larger time-step size for the same accuracy. Using a Newmark-type scheme (average acceleration method, = 1 2 , = 1 4 ), the slope is well captured, but oscillations appear when the crack is growing. Such oscillations are not observed with the TX-FEM and the solution remains satisfying. During the simulation of crack growth, new enriched shape functions are added to model the crack extension. Consequently, the numerical model is time discontinuous. As illustrated by the previous example, the TX-FEM is able to treat such phenomenon and the results are not disturbed by numerical oscillations.

Mixed mode dynamic crack growth and arrest.
With this last example, we want to show the capabilities of combined STX-FEM to simulate the dynamic propagation of a crack under mixed mode loading. The specimen consists of a plate with a circular hole and an off-centre crack (see Figure 12). As the crack is off centre, it will be submitted to a mixed mode loading. The material is assumed to be homogeneous and isotropic with elastic behaviour (material parameters are E = 210 GPa, = 0.25 and = 8000 kg m −3 ). It is submitted to a  Figure 13.
To predict the crack path, we use the maximum hoop stress criterion for the crack propagation direction: An equivalent mode 1 stress intensity factor K 1eq is defined from this max hoop stress criterion by The crack tip speed is then obtained using the following criterion: where K 1c is the dynamic fracture toughness that is assumed to be a constant whose value is 100 MPa √ m and c r is the Rayleigh wave speed that is also the theoretical maximum crack tip speed. More information about these criteria can be found in References [7,[20][21][22].
This particular configuration is very interesting because it has been experimentally observed by Carin and Maigre (see Reference [23]) that with an in-line initial crack, the crack runs at a constant speed after initiation, stops and reruns at the same constant speed. We suppose that in an experimental configuration the notch faces are separated by a gap so that the faces of the crack modelling the notch are allowed to overlap in the numerical model. As a consequence, a negative mode I stress intensity factor is developing when the compressive stress wave arrives on the crack tip.
The computational time is 900 s subdivided in 130 time steps. In Figure 14, the time history of the co-ordinate X f of the crack front in the loading axis is plotted. The first time, Figure 14. Crack front co-ordinate X f history for a circular hole. the compressive stress wave propagates through the specimen (0 < t < 2t c = L/c d ≈ 200 s). During this wave propagation, a negative mode I stress intensity factor develops once the wave reaches the crack front (t ≈ t c ). The wave is then reflected on the fixed end of the specimen, the crack faces open (2t c < t < 3t c ) and the crack starts growing (t ≈ 3t c ). Just after the initiation, the crack front reaches the loading axis of the specimen. Then the crack grows following this axis at an approximately constant speed (about 1300 ms −1 ). At time t ≈ 4t c , the main stress wave is reflected on the loaded face of the specimen and the crack stops when it reaches the crack front (t ≈ 5t c ). The crack starts again at the same constant speed when the main stress wave is reflected on the fixed face and arrives on the crack front (t ≈ 7t c ). Figures 15 and 16 show, respectively, the crack path when the crack arrests and the final crack path.
This geometrical configuration with a circular hole allows a simplified interpretation of the numerical results. As described in the previous paragraph, the propagation of the crack can be related to the wave propagation. Indeed, the circular hole does not disturb the main wave propagation a lot because of the shape of the free surface it produces. On the contrary, the square hole, which produces a plane free surface (see Figures 17 and 18), makes the wave   propagation more complicated with reflexions on the free surface and its interpretation much more difficult. The crack propagation is also completely different.

CONCLUSION
After a brief literature survey of the first time finite element approaches, enriched shape functions in time are added to the polynomial time approximation to give a TX-FEM. A particular case of enrichment function and choice of weight functions is studied and appears as an equivalent of T-DGM. A stability study is done for this case and the particular spectral radius of the amplification matrix of this method seems to be the key for its performances. The capabilities of the method to model time discontinuities without spurious numerical oscillations lie in an aperiodic response to numerical sollicitations of spurious high frequencies. The drawbacks of this over-damped behaviour are numerical damping and energy dissipation. Nevertheless, the method is third order accurate and gives less period error than the Newmark method. The cost of these time integration improvements in terms of numerical implementation and CPU time is balanced by an improved rate of convergence that allows a larger time-step size. For the particular obtained formulation in time (choice of enrichment function and weight functions), the method is an alternative approach of T-DGM. The framework of TX-FEM appears to be an attractive way for time integration improvements. Indeed, one could imagine to use other enriched functions for specific applications, to locate enrichment in the time slab and various ways of research to be explored. Combining space and time eXtended finite element methods gives a unified space-time discretization and accurate results when space and/or time discontinuities have to be simulated. Applications to dynamic crack propagation are presented in this paper. Avoiding spurious numerical oscillations, STX-FEM permits to simulate dynamic crack growth with a good accuracy. In the last example, the geometrical configuration gives interesting results concerning the crack behaviour. Results obtained using a square hole instead of a circular hole show that relating crack growth with a simplified wave propagation is not satisfying. Experimental results on such specimens would be interesting in order to improve crack physics knowledge. The contact between the crack faces would also be taken into account in the numerical simulations for more realistic models.