# A combined space-time extended Finite Element Method

Julien Réthoré, Anthony Gravouil, Alain Combescure

► To cite this version:

Julien Réthoré, Anthony Gravouil, Alain Combescure. A combined space-time extended Finite Element Method. International Journal for Numerical Methods in Engineering, 2005, pp.260-284. 10.1002/nme.1368. hal-00373829

HAL Id: hal-00373829

https://hal.science/hal-00373829v1

Submitted on 27 Apr 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# A combined space-time extended finite element method

J. Réthoré, A. Gravouil and A. Combescure*,†

LaMCoS, Laboratoire de Mécanique des Contacts et des Solides, UMR 5514, INSA Lyon,
Bat. Jean d'Alembert, 18, 20 rue des Sciences 69621 Villeurbanne, France

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.

KEY WORDS: extended finite element method; space-time formulation; accuracy; stability; energy conservation; dynamic fracture

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

Contract/grant sponsor: Office National d'Etudes et de Recherches Aérospatiales, centre de Lille; contract/grant number: F30075

<footnote>*Correspondence to: A. Combescure, LaMCoS, INSA Lyon, Bat. Jean d'Alembert, 18, 20 rue des Sciences 69621 Villeurbanne, France.</footnote>

<footnote>†E-mail: Alain.Combescure@insa-lyon.fr</footnote>
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–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).

## 2. 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 mass-spring system (parameters $k, m$) submitted to an external loading $f$. Here $x, \dot{x}, \ddot{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.

$$m\ddot{x} + kx - f = 0 \qquad (1)$$

### 2.1. 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, $\beta_N$ for $x$ and $\gamma_N$ for $\dot{x}$:

$$x_{n+1} = x_n + \Delta t \dot{x}_n + \frac{\Delta t^2}{2} \ddot{x}_n + \frac{\Delta t^3}{6} \left( 6\beta_N \frac{\ddot{x}_{n+1} - \ddot{x}_n}{\Delta t} \right) \quad (2)$$

$$\dot{x}_{n+1} = \dot{x}_n + \Delta t \ddot{x}_n + \frac{\Delta t^2}{2} \left( 2\gamma_N \frac{\ddot{x}_{n+1} - \ddot{x}_n}{\Delta t} \right) \quad (3)$$

The elastodynamics equation (Equation (1)) at time level $n$ is

$$m \ddot{x}_n + k x_n - f_n = 0 \quad (4)$$

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

$$\begin{aligned} & [m + \beta_N \Delta t^2 k] x_{n+1} + [-2m + (\frac{1}{2} + \gamma_N - 2\beta_N) \Delta t^2 k] x_n \\ & \quad + [m + (\frac{1}{2} - \gamma_N + \beta_N) \Delta t^2 k] x_{n-1} \\ & = \Delta t^2 [(\frac{1}{2} - \gamma_N + \beta_N) f_{n-1} + (\frac{1}{2} + \gamma_N - 2\beta_N) f_n + \beta_N f_{n+1}] \end{aligned} \quad (5)$$

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.

## 2.2. Literature survey in T-FEM

Zienkiewicz and then Wood presented other approaches using the weighted residual method and finite element concepts.

### 2.2.1. The Zienkiewicz approach [2].
In Reference [2], the displacement is interpolated in time as follows:

$$x(t) = \sum_{n=1}^{n+1} L_i(t) x_i \quad (6)$$

where $L_i$ are quadratic shape functions. Then, substituting Equation (6) into the weighted residual form of the elastodynamics equation (Equation (1)) integrated from $-\Delta t$ to $\Delta t$ ($W$ being any appropriate weight function), we have

$$\int_{-\Delta t}^{\Delta t} W(t) [m \sum \ddot{L}_i(t) x_i + k \sum L_i(t) x_i] dt = \int_{-\Delta t}^{\Delta t} W(t) f(t) dt \quad (7)$$

to obtain

$$\begin{aligned} & [m + \beta_Z \Delta t^2 k] x_{n+1} + [-2m + (\frac{1}{2} + \gamma_Z - 2\beta_Z) \Delta t^2 k] x_n \\ & \quad + [m + (\frac{1}{2} - \gamma_Z + \beta_Z) \Delta t^2 k] x_{n-1} \\ & = \Delta t^2 [(\frac{1}{2} - \gamma_Z + \beta_Z) f_{n-1} + (\frac{1}{2} + \gamma_Z - 2\beta_Z) f_n + \beta_Z f_{n+1}] \end{aligned} \quad (8)$$
where

$$2\beta_Z \frac{1}{\Delta t} \int_{-\Delta t}^{\Delta t} W(t) dt = \frac{1}{\Delta t^2} \int_{-\Delta t}^{\Delta t} W(t) t dt + \frac{1}{\Delta t^3} \int_{-\Delta t}^{\Delta t} W(t) t^2 dt \quad (9)$$

$$\gamma_Z \frac{1}{\Delta t} \int_{-\Delta t}^{\Delta t} W(t) dt = \frac{1}{2} + \frac{1}{\Delta t^2} \int_{-\Delta t}^{\Delta t} W(t) t dt$$

Equation (8) is formally the same as Equation (5), but $\beta_Z$ and $\gamma_Z$ are obtained from Equation (9). Consequently $\beta_Z$, $\gamma_Z$ are consistent with the polynomial interpolation and with the choice of the weight function $W$. It is not the case for $\beta_N$, $\gamma_N$, which can be chosen arbitrarily in Equations (2) and (3). As a consequence, we can say that an arbitrary choice of $\beta_N$, $\gamma_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.

### 2.2.2. The Wood approach [3]
Wood chooses to write the displacement with a quadratic interpolation in terms of three parameters $x_n$, $x_{n+1}$ and $\dot{x}_n$ as follows:

$$x(t) = x_n + t \dot{x}_n + \frac{t^2}{\Delta t^2} (x_{n+1} - x_n - \Delta t \dot{x}_n) \quad (10)$$

This expression is then introduced into the weighted residual form of the elastodynamics equation and integrated from 0 to $\Delta t$:

$$[2m + \beta_W \Delta t^2 k](x_{n+1} - x_n - \Delta t \dot{x}_n) - \Delta t^2 k x_n + \gamma_W \Delta t^2 \dot{x}_n = \Delta t^2 F \quad (11)$$

where $F$ is the right-hand side of Equation (8) and $\beta_W$, $\gamma_W$ are given by

$$\beta_W \frac{1}{\Delta t} \int_0^{\Delta t} W(t) dt = \frac{1}{\Delta t^3} \int_0^{\Delta t} W(t) t^2 dt \quad (12)$$

$$\gamma_W \frac{1}{\Delta t} \int_0^{\Delta t} W(t) dt = \frac{1}{\Delta t^2} \int_0^{\Delta t} W(t) t dt$$

and then $\dot{x}_{n+1}$ is given by differentiating Equation (10):

$$\Delta t \dot{x}_{n+1} = 2(x_{n+1} - x_n) - \Delta t \dot{x}_n \quad (13)$$

This method is equivalent to the Zienkiewicz approach and also to the Newmark method, but once again only for appropriate $\beta_N$, $\gamma_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].

## 2.3. 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}]$:

$$v(t) = v_n \lambda_n(t) + v_{n+1} \lambda_{n+1}(t) \quad (14)$$
where

$$ \lambda_n(t) = \frac{t_{n+1} - t}{t_{n+1} - t_n}; \quad \lambda_{n+1}(t) = \frac{t - t_n}{t_{n+1} - t_n} \qquad (15) $$

The displacement is then obtained by the integration of Equation (14) considering the initial condition $u(t_n) = u_n$:

$$ u(t) = u_n + \int_{t_n}^{t} v(\tau) \, d\tau \qquad (16) $$

The problem is then expressed with the weighted residual form of Equation (1) integrated from $t_n$ to $t_{n+1}$ ($\Delta 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

$$ [m + \beta \Delta t^2 k] v_{n+1} = -\Delta t k u_n + [m + (\beta - \gamma) \Delta t^2 k] v_n + \Delta t [(1 - \gamma) f_n + \gamma f_{n+1}] \quad (17) $$

where

$$ \beta \frac{1}{\Delta t} \int_{t_n}^{t_{n+1}} W(t) \, dt = \frac{1}{\Delta t^3} \int_{t_n}^{t_{n+1}} W(t)(t - t_n)^2 \, dt \qquad (18) $$

$$ \gamma \frac{1}{\Delta t} \int_{t_n}^{t_{n+1}} W(t) \, dt = \frac{1}{\Delta t^2} \int_{t_n}^{t_{n+1}} W(t)(t - t_n) \, dt $$

This equation can also be obtained from the Newmark method in the particular case $2\beta_N = \gamma_N$. Then we obtain $u_{n+1}$ using Equation (16):

$$ u_{n+1} = u_n + \frac{\Delta t}{2}(v_{n+1} + v_n) \qquad (19) $$

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\beta_N = \gamma_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

$$ a(t) = a_{n+1} = \frac{1}{\Delta t}(v_{n+1} - v_n) \qquad (20) $$

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 $\ddot{u}(t_n) = \ddot{u}_n$ and $\dot{v}(t_n) = \dot{v}_n$:

$$ u_{n+1} = u_n + \Delta t v_n + \frac{\Delta t^2}{2} \ddot{u}_n \qquad (21) $$
$$v_{n+1} = v_n + \Delta t \dot{v}_n \qquad (22)$$

To be consistent with the Newmark method, one could define $\ddot{u}_n$ and $\dot{v}_n$ by

$$\ddot{u}_n = (1 - 2\beta_N)a_n + 2\beta_N a_{n+1} \qquad (23)$$

$$\dot{v}_n = (1 - \gamma_N)a_n + \gamma_N a_{n+1}$$

Those definitions are kinematically consistent if $\ddot{u}_n = \dot{v}_n$, i.e. if $2\beta_N = \gamma_N$. Indeed, in the Newmark method, kinematic relations between displacement velocity and acceleration are not exactly satisfied for arbitrary $\beta_N, \gamma_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.

### 3. TIME EXTENDED FINITE ELEMENT METHOD

#### 3.1. 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 {$\lambda_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:

$$v(t) = \sum_{i=0}^{N} \lambda_i(t) v_i^c + \sum_{j=0}^{M} \sum_{i \in \mathcal{N}_j} \lambda_i(t) \phi_j(t) v_{i,j}^c \qquad (24)$$

where $v_i^c$ are classical degrees of freedom, $\phi_j$ are enriched shape functions, $v_{i,j}^c$ additional degrees of freedom and $\mathcal{N}_j$ the set of time levels that are enriched with $\phi_j$. The $\phi_j$ can be chosen to capture the non-polynomial part of the solution with good accuracy.

#### 3.2. Heaviside enrichment

##### 3.2.1. 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 $\phi_j = H(t - t_j)$ centred at each time level, and $\mathcal{N}_j$ is reduced to the time level $t_j$. Because $\phi_j = H(t - t_j)$ is zero for all time $t$ earlier than $t_j$ and $\lambda_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:

$$v(t) = v^c(t) + v^e(t) \qquad (25)$$

$$u(t) = u^c(t) + u^e(t)$$
<figure><img src="image_2.png" /><figcaption>Figure 1. Linear shape function and activated Heaviside function in the time-slab $I_{n+1}$.</figcaption></figure>

where

$$ 
\begin{aligned}
v^c(t) &= v_n^c \lambda_n(t) + v_{n+1}^c \lambda_{n+1}(t) \\
v^e(t) &= v_{n+1}^e \lambda_n(t) H(t - t_n)
\end{aligned} 
\qquad (26) 
$$

and

$$ \lambda_n(t) = \frac{t_{n+1} - t}{t_{n+1} - t_n}; \quad \lambda_{n+1}(t) = \frac{t - t_n}{t_{n+1} - t_n} \qquad (27) $$

We abusively denote the additional degree of freedom $v_{n,n}^e$ by $v_{n+1}^e$ because it has to be determined when the time slab $I_{n+1}$ is solved.

The kinematic constraints between displacement and velocity ($\dot{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

$$ 
\begin{aligned}
u^c(t) &= u_n^c + \int_{t_n}^t v^c(\tau) \, d\tau \\
u^e(t) &= u_n^e + \int_{t_n^+}^t v^e(\tau) \, d\tau
\end{aligned} 
\qquad (28) 
$$

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$:

$$ \int_{t_n^+}^{t_{n+1}^-} W(t)[m\dot{v}(t) + ku(t)]dt + W(t_n^+)mv_{n+1}^e = \int_{t_n^+}^{t_{n+1}^-} W(t)f(t)dt \qquad (29) $$

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_n^c, u_n^e, v_n^c$ and $v_n^e$, find $u_{n+1}^c, u_{n+1}^e, v_{n+1}^c$ and $v_{n+1}^e$ such that

$$ \mathbf{H}_1 \mathbf{q}_{n+1} = \mathbf{H}_0 \mathbf{q}_n + \mathbf{F}_{n+1} \qquad (30a) $$
where

$$ H_1 = \begin{bmatrix} 1 & 0 & -\frac{\Delta t}{2} & 0 \\ 0 & 1 & 0 & -\frac{\Delta t}{2} \\ 0 & 0 & m + \beta_1 \Delta t^2 k & (\delta_1 - 1)m + (\gamma_1 - \beta_1) \Delta t^2 k \\ 0 & 0 & m + \beta_2 \Delta t^2 k & (\delta_2 - 1)m + (\gamma_2 - \beta_2) \Delta t^2 k \end{bmatrix} \quad (30b) $$

$$ q_n = \begin{bmatrix} u_n^c \\ u_n^e \\ v_n^c \\ v_n^e \end{bmatrix} \quad (30c) $$

$$ H_0 = \begin{bmatrix} 1 & 0 & \frac{\Delta t}{2} & 0 \\ 0 & 1 & 0 & 0 \\ -\Delta tk & -\Delta tk & m + (\beta_1 - \gamma_1)\Delta t^2 k & 0 \\ -\Delta tk & -\Delta tk & m + (\beta_2 - \gamma_2)\Delta t^2 k & 0 \end{bmatrix} \quad (30d) $$

$$ F_{n+1} = \begin{bmatrix} 0 \\ 0 \\ \Delta t((1 - \gamma_1)f_n + \gamma_1 f_{n+1}) \\ \Delta t((1 - \gamma_2)f_n + \gamma_2 f_{n+1}) \end{bmatrix} \quad (30e) $$

and

$$ \beta_i \frac{1}{\Delta t} \int_{t_n^+}^{t_{n+1}^-} W_i(t) \, dt = \frac{1}{\Delta t^3} \int_{t_n^+}^{t_{n+1}^-} W_i(t) (t - t_n^+)^2 \, dt \\ \gamma_i \frac{1}{\Delta t} \int_{t_n^+}^{t_{n+1}^-} W_i(t) \, dt = \frac{1}{\Delta t^2} \int_{t_n^+}^{t_{n+1}^-} W_i(t) (t - t_n^+) \, dt \quad (31) $$

$$ \delta_i \frac{1}{\Delta t} \int_{t_n^+}^{t_{n+1}^-} W_i(t) \, dt = W_i(t_n^+) $$

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 constraint $\dot{u} = v$ is always satisfied. This time integrator is defined by six parameters $\{\beta_i, \gamma_i, \delta_i\}_{i=1,2}$ which are computed from the weight functions $\{W_i\}_{i=1,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(t) = a_{n+1} = \frac{1}{\Delta t} (v_{n+1}^c - v_{n+1}^e - v_n^c) \quad (32)$$

### 3.2.2. A new look to the T-DGM.
As it is possible to have a further look at the Newmark method using T-FEM, XT-FEM proposes a new look at T-DGM. Let us write

$$\begin{aligned} v_n^+ &= v_n^c + v_{n+1}^e \\ v_{n+1}^- &= v_{n+1}^c \end{aligned} \quad (33)$$

and choose the weight functions to be

$$\begin{aligned} W_1 &= \lambda_n \\ W_2 &= \lambda_{n+1} \end{aligned} \quad (34)$$

The parameters are

$$\begin{aligned} \beta_1 &= \frac{1}{12}, & \gamma_1 &= \frac{1}{3}, & \delta_1 &= 2 \\ \beta_2 &= \frac{1}{4}, & \gamma_2 &= \frac{2}{3}, & \delta_2 &= 0 \end{aligned} \quad (35)$$

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

$$\begin{gathered} u_{n+1} = u_n + \frac{\Delta t}{2} (v_{n+1}^{-} + v_n^{+}) \\ \begin{bmatrix} m + \frac{\Delta t^2}{12} k & -\frac{\Delta t^2}{12} k \\ \frac{\Delta t^2}{3} k & m + \frac{\Delta t^2}{6} k \end{bmatrix} \begin{bmatrix} v_n^+ \\ v_{n+1}^- \end{bmatrix} = \begin{bmatrix} m v_n^- + \frac{\Delta t^2}{6} (f_n - f_{n+1}) \\ m v_n^- - \Delta t k u_n + \frac{\Delta t^2}{2} (f_n + f_{n+1}) \end{bmatrix} \end{gathered} \quad (36)$$

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 ($\dot{u} = v$) is strongly enforced, P3-P1 approximation allows high order approximation but the kinematic constraint is only weakly enforced.
### 3.2.3. 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 non-zero initial velocity. Weight functions for the numerical results are $W_1 = \lambda_n$, $W_2 = \lambda_{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

$$ q_{n+1} = Aq_n \quad (37) $$

$A = H_1^{-1}H_0$ is the amplification matrix and the stability depends only upon the eigenvalues of $A$. They are obtained solving

$$ \det(H_0 - rH_1) = 0 \quad (38) $$

It can be shown that this equation is equivalent to

$$ r(r-1)(\alpha_2 r^2 - \alpha_1 r + \alpha_0) = 0 \quad (39) $$

where

$$ \alpha_2 = \left[ \frac{\omega^2 \Delta t^2}{\zeta^2} \right] \left\langle \delta - \frac{\omega^2 \Delta t^2}{\zeta^2} (1 - \gamma \zeta^2) \right\rangle - \left\langle \frac{\omega^2 \Delta t^2}{\zeta^2} \right\rangle \left[ \delta - \frac{\omega^2 \Delta t^2}{\zeta^2} (1 - \gamma \zeta^2) \right] \quad (40) $$

$$ \alpha_1 = \left[ \frac{\omega^2 \Delta t^2}{\zeta^2} \left( 2 - \left( \gamma + \frac{1}{2} \right) \zeta^2 \right) \right] \left\langle \delta - \frac{\omega^2 \Delta t^2}{\zeta^2} \left( 1 - \left( \gamma - \frac{1}{2} \right) \zeta^2 \right) \right\rangle - \left\langle \frac{\omega^2 \Delta t^2}{\zeta^2} \left( 2 - \left( \gamma + \frac{1}{2} \right) \zeta^2 \right) \right\rangle \left[ \delta - \frac{\omega^2 \Delta t^2}{\zeta^2} \left( 1 - \left( \gamma - \frac{1}{2} \right) \zeta^2 \right) \right] \quad (41) $$

$$ \alpha_0 = \left[ \frac{\omega^2 \Delta t^2}{\zeta^2} \left( 1 - \left( \gamma - \frac{1}{2} \right) \zeta^2 \right) \right] \langle \delta \rangle - \left\langle \frac{\omega^2 \Delta t^2}{\zeta^2} \left( 1 - \left( \gamma - \frac{1}{2} \right) \zeta^2 \right) \right\rangle [\delta] \quad (42) $$

with $\omega = \sqrt{k/m}$, $\zeta^2 = 1 + \beta\omega^2\Delta t^2$ and $\langle b \rangle = \frac{1}{2}(b_2 + b_1)$, $[b] = b_2 - b_1$ for constant parameters. The eigenvalues 0 and 1 correspond to trivial modes of the updating equations. The interest is now on the second order equation $\alpha_2 r^2 - \alpha_1 r + \alpha_0$. In Equations (40)–(42), one can recognize, inside the operators $\langle \cdot \rangle$ and $[\cdot, \cdot]$, terms which govern the stability of the Newmark method. $r_1, r_2$ are the roots of this equation. The spectral radius of the amplification matrix is defined by $\rho = \max(|r_1|, |r_2|)$. Then the stability conditions are (see Reference [12])

(i) $\rho$ 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 ($\gamma_N = \frac{1}{2}, \beta_N = \frac{1}{4}$). Those spectral radii are plotted (see Figure 2) versus parameter $\omega\Delta t$. The first observation is that the method is unconditionally stable for this set of weight functions ($W_1 = \lambda_n, W_2 = \lambda_{n+1}$). We can observe that the spectral radius is not equal to one and the method will show numerical damping. Two regions of $\omega\Delta t$ ($\omega\Delta t \approx \pi$ and $\omega\Delta t \ge 10$) are particular. For these regions, we can observe bifurcations because
<figure><img src="image_3.png" /><figcaption>Figure 2. Spectral radii for the Newmark method ($γ_N = \frac{1}{2}$, $β_N = \frac{1}{4}$), HHT ($α = -0.3$) and TX-FEM.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 3. Numerical damping ratios for the Newmark method ($γ_N = \frac{1}{2}$, $β_N = \frac{1}{4}$), HHT ($α = -0.3$) and TX-FEM.</figcaption></figure>

the roots of Equation (38) separate into real distinct values (and not complex conjugate). This implies that for such values of $ωΔt$ the numerical response of the method is aperiodic:

$$q_n = \sum e^{-(n\Delta t/\tau_i)} \eta_i$$

with $η_i$ the eigenvector corresponding to the real eigenvalue $r_i = e^{-(\Delta t/\tau_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
<figure><img src="image_5.png" /><figcaption>Figure 4. Relative period errors for the Newmark method ($\gamma_N = \frac{1}{2}$, $\beta_N = \frac{1}{4}$), HHT ($\alpha = -0.3$) and TX-FEM.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 5. Errors in $L_2$ norm for the Newmark method ($\gamma_N = \frac{1}{2}$, $\beta_N = \frac{1}{4}$) and TX-FEM.</figcaption></figure>

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.

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

$$ \langle x \rangle = \frac{1}{2}(x_{n+1}^{-} + x_{n}^{+}) $$
$$ 
\begin{aligned}
[x] &= x_{n+1}^{-} - x_{n}^{+} \\
[|x|] &= x_{n}^{+} - x_{n}^{-}
\end{aligned}
\qquad (43) $$

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

$$ 
\begin{aligned}
& \left[\frac{1}{2}mv^2\right] + \left[\frac{1}{2}ku^2\right] + \left\langle\beta - \frac{\gamma}{2}\right\rangle \Delta t^2 \left[\frac{1}{2}kv^2\right] + 2\left\langle\gamma - \frac{1}{2}\right\rangle \frac{1}{2}k[u]^2 + \langle\delta\rangle m\langle v\rangle [|v|] \\
& = \left\langle\gamma - \frac{1}{2}\right\rangle [u][f] + [u]\langle f\rangle
\end{aligned}
\qquad (44a) $$

$$ \left[ \beta - \frac{\gamma}{2} \right] \Delta t^2 \left[ \frac{1}{2} k v^2 \right] + 2 \left[ \gamma - \frac{1}{2} \right] \frac{1}{2} k[u]^2 + [\delta] m \langle v \rangle [|v|] = \left[ \gamma - \frac{1}{2} \right] [u][f] \quad (44b) $$

The set of parameters for $W_1 = \lambda_n$, $W_2 = \lambda_{n+1}$ (Equation (35)) gives

$$ 
\begin{aligned}
\left\langle\beta - \frac{\gamma}{2}\right\rangle &= -\frac{1}{12}, & \left[\beta - \frac{\gamma}{2}\right] &= 0 \\
\left\langle\gamma - \frac{1}{2}\right\rangle &= 0, & \left[\gamma - \frac{1}{2}\right] &= \frac{1}{3} \\
\langle\delta\rangle &= -2, & [\delta] &= 1
\end{aligned}
\qquad (45) $$

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)). As soon as the energy created $m\langle v\rangle [|v|] = (2[\gamma - \frac{1}{2}]/[\delta])\frac{1}{2}k[u]^2$ comes from a dissipative term the formulation is not conservative. If we denote $W = \frac{1}{2}ku^2$, $T = \frac{1}{2}mv^2$, $N = \langle\beta - \frac{\gamma}{2}\rangle\Delta t^2 \frac{1}{2}kv^2$, the formulation would be energy preserving if

$$ [|W + T + N|] + [W + T + N] = 0 $$

Using Equation (44a), numerical energy conservation is ensured if

$$ [|W + T + N|] = 2\langle\gamma - \frac{1}{2}\rangle\frac{1}{2}k[u]^2 + \langle\delta\rangle m\langle v\rangle [|v|] $$

Using the continuity of $u$ and Equation (44b)

$$ 
\begin{aligned}
& \left[\left|\frac{1}{2}mv^2\right|\right] + \left\langle\beta - \frac{\gamma}{2}\right\rangle \Delta t^2 \left[\left|\frac{1}{2}kv^2\right|\right] \\
& = 2\left\langle\gamma - \frac{1}{2}\right\rangle \frac{1}{2}k[u]^2 - \frac{\langle\delta\rangle}{[\delta]} \left(\left[\beta - \frac{\gamma}{2}\right] \Delta t^2 \left[\frac{1}{2}kv^2\right] + 2\left[\gamma - \frac{1}{2}\right] \frac{1}{2}k[u]^2\right)
\end{aligned}
\qquad (46) $$
<figure><img src="image_7.png" /><figcaption>Figure 6. Energy balance of TX-FEM discretized formulation.</figcaption></figure>

Hence,

$$ \frac{1}{2} \left( m + \left( \beta - \frac{\gamma}{2} \right) \Delta t^2 k \right) [|v|^2|] \\ = \frac{1}{2} \Delta t^2 k \left( 2 \left( \left( \gamma - \frac{1}{2} \right) - \frac{\langle \delta \rangle}{[\delta]} \left[ \gamma - \frac{1}{2} \right] \right) \langle v \rangle - \frac{\langle \delta \rangle}{[\delta]} \left[ \beta - \frac{\gamma}{2} \right] [v] \right) \langle v \rangle \quad (47) $$

Consequently, for this set of weight functions ($W_1 = \lambda_n$, $W_2 = \lambda_{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 $\omega\Delta t \approx 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 = \lambda_n$, $W_2 = \lambda_{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.
<figure><img src="image_8.png" /><figcaption>Figure 7. Illustration of the numerical dissipation.</figcaption></figure>

## 4. 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:

$$U(x, t) = U^h(t) \qquad (48)$$

This stage is detailed in the following paragraph. Then a time approximation is made using the method proposed in Section 3 to write

$$U(x, t) = U^h, \Delta t \qquad (49)$$

### 4.1. 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]:

$$U^h = \sum_{i \in \mathcal{N}} N_i(x)u_i + \sum_{i \in \mathcal{N}_{\text{cut}}} N_i(x)\mathcal{H}(x)a_i + \sum_{i \in \mathcal{N}_{\text{branch}}} \sum_{\alpha} N_i(x)B_{\alpha}(x)b_{i, \alpha} \qquad (50)$$

where $\mathcal{N}$ is the set of all nodes in the mesh, $\mathcal{N}_{\text{cut}}$ the set of nodes whose support is completely cut by the crack and $\mathcal{N}_{\text{branch}}$ the set of nodes that belong to elements partially cut by the crack. $N_i$ are the classical shape functions and $\mathcal{H}$ is a discontinuous function, whose value is 1 if $x$ is above the crack surface and $-1$ if $x$ is below. $[B_{\alpha}]$ are branch functions $((r, \theta))$ are local cylindrical crack tip co-ordinates):

$$[B_{\alpha}] = \left[ \sqrt{r} \sin\left(\frac{\theta}{2}\right), \sqrt{r} \cos\left(\frac{\theta}{2}\right), \sqrt{r} \sin\left(\frac{\theta}{2}\right) \sin(\theta), \sqrt{r} \cos\left(\frac{\theta}{2}\right) \sin(\theta) \right] \qquad (51)$$
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.

## 4.2. Time discretization

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

$$V^{h, \Delta t} = \sum_{i=0}^{N} \lambda_i(t)(V_i^{h,c} + H(t - t_i)V_{i+1}^{h,e}) \quad (52)$$

## 4.3. 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_{\text{branch}}^{n+1}$. For discontinuous enrichment, new shape functions are only added on the set $N_{\text{cut}}^{\text{new}} = N_{\text{cut}}^{n+1} \setminus N_{\text{cut}}^n$. 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_x]$) in time-dependent problem. It has been validated 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 numerical 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.

## 4.4. Examples

### 4.4.1. 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 $\rho = 8000$ kg m$^{-3}$). It is submitted to a tensile stress $\sigma = 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
<figure><img src="image_9.png" /><figcaption>Figure 8. Strategy for enrichment when the crack is growing.</figcaption></figure>

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\Delta t$ is solicited. For such solutions, the corresponding value of the parameter $\omega\Delta t = 2\pi\Delta t/T$ is $\pi$. Looking at the spectral
<figure><img src="image_10.png" /><figcaption>Figure 9. Normalized displacement, velocity and acceleration for the Newmark average acceleration method and TX-FEM at the free end of the beam (left) and at the created interface of the fixed half-broken beam (right).</figcaption></figure>

radius of the TX-FEM, $\pi$ 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.

4.4.2. 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><img src="image_11.png" /><figcaption>Figure 10. Geometry and loading for the example of the infinite plate.</figcaption></figure>

geometry described in Figure 10, the analytical solution is only valid for time $t \le 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

$$K_1^{\text{dyn}}(0, t) = \frac{2\sigma_0}{1-\mu} \sqrt{\frac{c_1 t (1-2\mu)}{\pi}} \quad (53)$$

For a moving crack we can write

$$K_1^{\text{dyn}}(\dot{a}, t) = k(\dot{a}) K_1^{\text{dyn}}(0, t) \quad (54)$$

where $k$ is a universal function of the crack tip velocity $\dot{a}$, which can be approximated by

$$k(\dot{a}) = \frac{1 - (\dot{a}/c_r)}{1 - (\dot{a}/2c_r)} \quad (55)$$

And finally, we have

$$K_1^{\text{dyn}}(\dot{a}, t) = \frac{2\sigma_0}{1-\mu} \sqrt{\frac{c_d t (1-2\mu)}{\pi} \frac{1-(\dot{a}/c_r)}{1-(\dot{a}/2c_r)}} \quad (56)$$

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:

$$I^{\text{int}} = - \int_S q_{k,j} [(\sigma_{ml}^{\text{aux}} u_{m,l} - \rho \dot{u}_l \dot{u}_l^{\text{aux}}) \delta_{kj} - (\sigma_{ij}^{\text{aux}} u_{i,k} + \sigma_{ij} u_{i,k}^{\text{aux}})] dS \\ + \int_S q_k [(\sigma_{ij,j}^{\text{aux}} u_{i,k} + \ddot{u}_i u_{i,k}^{\text{aux}}) + (\rho \dot{u}_i^{\text{aux}} \dot{u}_{i,k} + \rho \dot{u}_i \dot{u}_{i,k}^{\text{aux}})] dS \quad (57)$$

The interaction integral (see Equation (57)) is computed using a J-domain which consists of an additional mesh (only used for numerical integration) following the crack tip. We are interested in the case when the crack propagates at a prescribed constant speed $v_0$ after time $t = 1.5t_c$. The following numerical results are obtained with $H = 2$ m, $L = 10$ m and $l = 5$ m for plate dimensions and $E = 210$ GPa, $v = 0.3$ and $\rho = 8000$ kg m$^{-3}$ for material properties. The tensile stress $\sigma_0$ is 500 MPa. $v_0$ is 1500 m s$^{-1}$ and stress intensity factors are normalized by $\sigma_0 \sqrt{H}$.
<figure><img src="image_12.png" /><figcaption>Figure 11. Numerical and analytical solutions $\bar{K}_1$ for a stationary then moving crack.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 12. Geometry and loading in the mixed mode configuration.</figcaption></figure>

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, $\gamma = \frac{1}{2}$, $\beta = \frac{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.

### 4.4.3. 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, $v = 0.25$ and $\rho = 8000$ kg m$^{-3}$). It is submitted to a
<figure><img src="image_14.png" /><figcaption>Figure 13. Mesh and initial geometry.</figcaption></figure>

compressive stress $\sigma = 500$ MPa at one end and fixed at the other. The dimensions are $L = 1$ m, $H = 0.5$ m, $l = L/2 = 0.5$ m, $r = 0.125$ and $e = 0.05$ m. The numerical simulation uses space enrichment for modelling the hole and the crack. The geometry is then not explicitly defined as shown in Figure 13.

To predict the crack path, we use the maximum hoop stress criterion for the crack propagation direction:

$$ \theta_c = 2 \arctan \left[ \frac{1}{4} \left( \frac{K_1}{K_2} - \operatorname{sign}(K_2) \sqrt{8 + \left( \frac{K_1}{K_2} \right)^2} \right) \right] \quad (58) $$

An equivalent mode 1 stress intensity factor $K_{1eq}$ is defined from this max hoop stress criterion by

$$ K_{1eq} = \cos^3 \left( \frac{\theta_c}{2} \right) K_1 - \frac{3}{2} \cos \left( \frac{\theta_c}{2} \right) \sin (\theta_c) K_2 \quad (59) $$

The crack tip speed is then obtained using the following criterion:

$$ \dot{a} = \begin{cases} 0 & \text{if } K_{1eq} \le K_{1c} \\ c_r \left( 1 - \frac{K_{1c}}{K_{1eq}} \right) & \text{otherwise} \end{cases} \quad (60) $$

where $K_{1c}$ is the dynamic fracture toughness that is assumed to be a constant whose value is $100 \text{ MPa}\sqrt{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–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><img src="image_15.png" /><figcaption>Figure 14. Crack front co-ordinate $X_f$ history for a circular hole.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure 15. Crack path when the crack arrests (circular hole).</figcaption></figure>

the compressive stress wave propagates through the specimen ($0 < t < 2t_c = L/c_d \approx 200$ µs). During this wave propagation, a negative mode I stress intensity factor develops once the wave reaches the crack front ($t \approx 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 \approx 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 \text{ ms}^{-1}$). At time $t \approx 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 \approx 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 \approx 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
<figure><img src="image_17.png" /><figcaption>Figure 16. Final crack path (circular hole).</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Figure 17. Crack front co-ordinate $X_f$ history for a square hole.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Figure 18. Final crack path (square hole).</figcaption></figure>
propagation more complicated with reflexions on the free surface and its interpretation much
more difficult. The crack propagation is also completely different.

5. CONCLUSION

After a brief literature survey of the first time finite element approaches, enriched shape func-
tions 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 equiv-
alent 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 dis-
cretization 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 nu-
merical 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.

REFERENCES

1. Newmark NM. A method of computation for structural dynamics. *Proceedings of the American Society of Civil Engineers* 1959; **85**:67–94.

2. Zienkiewicz OC. A new look at Newmark, houbolt and other time stepping formulae, a weighted residual approach. *Earthquake Engineering and Structural Dynamics* 1977; **5**:413–418.

3. Wood WL. A further look at Newmark, Houbolt, etc., time-stepping formulae. *International Journal for Numerical Methods in Engineering* 1984; **20**:1009–1017.

4. Zienkiewicz OC, Wood WL, Hine NW. A unified set of single step algorithms. Part i: general formulation and applications. *International Journal for Numerical Methods in Engineering* 1984; **20**:1529–1552.

5. Babuska I, Melenk JM. The partition of unity method. *International Journal for Numerical Methods in Engineering* 1997; **40**:727–758.

6. Wells GN, de Borst R, Sluys LJ. A consistent geometrically non-linear approach for delamination. *International Journal for Numerical Methods in Engineering* 2002; **54**:1333–1355.
7. Duarte CA, Hamzeh ON, Liska TJ, Tworzydlo WW. A generalized finite element method for the simulation of three-dimensional dynamic crack propagation. *Computer Methods in Applied Mechanics and Engineering* 2001; **190**:2227-2262.

8. Belytschko T, Moës N, Usui S, Parimi C. Arbitrary discontinuities in finite elements. *International Journal for Numerical Methods in Engineering* 2001; **50**(4):993-1013.

9. Li X, Yao D, Lewis RW. A discontinuous Galerkin finite element method for dynamic and wave propagation problems in non-linear solids and saturated porous media. *International Journal for Numerical Methods in Engineering* 2003; **57**:1775-1800.

10. Michler C, van Brummelen EH, Hulshoff SJ, de Borst R. The relevance of conservation for stability and accuracy of numerical methods for fluid-structure interaction. *Computer Methods in Applied Mechanics and Engineering* 2003; **192**:4195-4215.

11. Chessa J, Belytschko T. Arbitrary discontinuities in space-time finite elements by levelset and X-FEM. *International Journal for Numerical Methods in Engineering* 2004; **61**(15):2595-2614.

12. Hughes TJR, Belytschko T. *Nonlinear Finite Element Analysis*. ICE Division, Zace Services Ltd., 2000.

13. Black T, Belytschko T. Elastic crack growth in finite elements with minimal remeshing. *International Journal for Numerical Methods in Engineering* 1999; **45**:601-620.

14. Moës N, Dolbow J, Belytschko T. A finite element method for crack growth without remeshing. *International Journal for Numerical Methods in Engineering* 1999; **46**(1):133-150.

15. Chessa J, Wang H, Belytschko T. On the construction of blending elements for local partition of unity enriched finite elements. *International Journal for Numerical Methods in Engineering* 2003; **57**:1015-1038.

16. Réthoré J, Gravouil A, Combescure A. An energy conserving scheme for dynamic crack growth with the extended finite element method. *International Journal for Numerical Methods in Engineering* 2005; **63**(5):631-659.

17. Réthoré J, Gravouil A, Combescure A. A stable numerical scheme for the finite element simulation of dynamic crack propagation with remeshing. *Computer Methods in Applied Mechanics and Engineering* 2004; **193**:4493-4510.

18. Belytschko T, Chen H. Singular enrichment finite element method for elastodynamic crack propagation. *International Journal of Computational Methods* 2004; **1**(1):1-15.

19. Belytschko T, Chen H, Jingxiao X, Goangseup Z. Dynamic crack propagation based on loss of hyperbolicity and a new discontinuous enrichment. *International Journal for Numerical Methods in Engineering* 2003; **58**:1873-1905.

20. Freund LB. *Dynamic Fracture Mechanics*. Cambridge Monographs on Mechanics and Applied Mathematics, 1990.

21. Organ D. Numerical solutions to dynamic fracture problems using the element-free Galerkin methods. Ph.D. Thesis, Northwestern University, 1996.

22. Krysl P, Belytschko T. Dynamic propagation of arbitrary 3-d cracks. *International Journal for Numerical Methods in Engineering* 1999; **44**(6):767-800.

23. Carin T. Modélisation de la propagation dynamique de fissure. Ph.D. Thesis, École Nationale des Ponts et Chaussées, 2000.