# A space-time approach for the simulation of brittle fracture with phase-field models in elastodynamics

Fabrice Kevine Feutang, Stéphane Lejeunes, Dominique Eyheramendy

► To cite this version:

Fabrice Kevine Feutang, Stéphane Lejeunes, Dominique Eyheramendy. A space-time approach for the simulation of brittle fracture with phase-field models in elastodynamics. Computers & Structures, 2025, 307, pp.107668. 10.1016/j.compstruc.2024.107616. hal-04830169

HAL Id: hal-04830169

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

Submitted on 10 Dec 2024

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 Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# A space-time approach for the simulation of brittle fracture with phase-field models in elastodynamics

F.K. Feutang, S. Lejeunes *, D. Eyheramendy

LMA, UMR 7031, Aix-Marseille Univ, CNRS, Centrale Marseille, Marseille, France

## ARTICLE INFO

### Keywords:

IGA

Phase-field

Space-time

Time discontinuous Galerkin

Damage

Fragmentation

## ABSTRACT

A space-time approach is proposed to simulate the propagation of brittle cracks in an isotropic and elastic solid in dynamics. We adopt the so called phase-field description of crack that is based on a variational representation of fracture mechanics. Due to this variational structure, the crack initiation and propagation can be then described thanks to a well chosen potential. In this approach, we propose to consider a space-time potential to derive the appropriate Euler equations on the space-time domain. A time discontinuous Galerkin approach is used and adapted to damage and elastodynamics such as to be able to account of time singularities in the considered fields. This approach follows a previous work done on elastodynamics (see [47]) in which we have proposed a stabilized formulation with the help of least square terms. The proposed space-time potential is discretized with either standard finite-elements (ST-FE) or isogeometric analysis (ST-IGA). We apply this approach to different numerical examples including dynamic fragmentation.

## 1. Introduction

Space-time approaches for computational mechanics (non-relativistic) go back to the work of Oden, Argyris and Sharpf who were the first to formulate space-time finite-elements for solving linear time-dependent problems. The main idea was to consider time as an additional dimension of the problem (but still in a Lagrangian framework, time is independent of space and vis-versa, see [42,6]). Following these first works, one can find in the literature some applications of this method in particular in elastodynamics. In [46] the space-time finite-element is used to treat a 1D viscoelasticity problem. One can also find applications to heat transfer problems [14], free surface problems [28], or advection-diffusion problems [40]. Hughes and Hulbert [27], as well as Hulbert and Hughes [26], are notable for their groundbreaking contributions to time-discontinuous Galerkin methods. Their work introduced a way to build unconditionally stable formulations for elastodynamics on space-time unstructured meshes, all the while accommodating discontinuities in the primary variables within finite time intervals. In reference [29], it was demonstrated that time-discontinuous Galerkin methods have the capability to effectively represent discontinuities, especially when dealing with problems that exhibit a discontinuous solution (e.g. crack propagation, impact, etc). These methods offer the advantage of enabling the definition of higher-order time-stepping schemes.

Few years latter, with the introduction of the isogeometric analysis (IGA) by Hughes and his co-authors (e.g. [25]), space-time IGA approaches start to emerge. The first interest of IGA is to obtain an exact discrete representation of complex geometries, but the IGA presents also a great flexibility in the definition of the approximation spaces (because it allows to control both the order of the functions and their order of continuity at some points). ST-IGA allows to deal efficiently with different problematic such as parabolic problems [30], elliptic and parabolic nonlinear problems [22] and advection-diffusion equations [8]. Nowadays, very complex space-time simulations can be realized with such approaches (see for instance [50] and the previous works of the authors)

In line with the previously mentioned works, we have proposed to use IGA in a space-time framework for various linear and nonlinear problems (elastodynamics, thermomechanical couplings or viscoplasticity, see [48]). In particular, we have proposed in [47] a new space-time formulation in a continuous Galerkin framework that exploits a Hellinger-Reisner transform of the kinetic energy to define new stabilization terms to control spurious numerical oscillation due to the discretization space. We obtained a two fields formulation (velocity/displacement) but the stabilization terms are general and in particular independent of the mechanical behavior. Therefore the extension of this approach to dynamics gradient damage models seems natural. Contrarily or complementary to what we have proposed in [47], we focus

<footnote>* Corresponding author.
E-mail address: lejeunes@lma.cnrs-mrs.fr (S. Lejeunes).</footnote>
here on a space-time discontinuous Galerkin method for the aforementioned reasons (time discontinuities capturing). For the phase-field part of our work we use standard formulations from the literature. In the late 1990s, the phase field method was used to model brittle fracture in elastic solids, and since then it has attracted a great interest in cracking problems in general. One intriguing aspect of the phase-field method lies in its inherent capacity to initiate, propagate, and even bifurcate cracks seamlessly, without any extra ad hoc strategies. The work of Francfort and Marigo [18] and Frémond and Nedjar [19] laid the foundation for quasi-static brittle fracture using the phase field method. In [18], authors introduced the concept of phase field models in fracture, which approach fracture as a problem of minimizing potential energy and surface energy. This energetic formulation aligns with Griffith's theory and offers a powerful framework to increase complexity and develop new models. Subsequently, drawing inspiration from Ambrosio and Tortorelli's [4] elliptic regularization of the Mumford and Shah [39] functional in image segmentation, Bourdin et al. [11] introduced a regularized form of the mathematical description of the phase-field model for brittle fracture proposed by Francfort and Marigo. There exist an extensive literature on quasi-static phase-field formulations and gradient damage models for brittle fracture adapted to different materials and situations, interested readers can refer for instance to [35,44,5,1,15,51,13,34,7,9] and references therein. Under the assumption that the damage field does not influence explicitly the kinetic energy (due to the fact that mass is conserved), phase-field damage models have been effectively applied to the study of dynamic fracture in brittle materials, see for instance [10,23,54,36]. Nevertheless, phase-field is numerically very expensive especially for dynamics problems due both to the non-convexity of the potential to be minimized and to the requirement of capturing accurately the damage gradient (defined by its characteristic length) and the stress/strain waves due to dynamics. The algorithm chosen for the minimization, the space discretization and the time scheme are crucial ingredients to capture brittle cracks with a controlled numerical error and avoiding numerical instabilities (due to spurious oscillations in a dynamic regime for instance). Our proposal is therefore to investigate the interest of using ST-FE or ST-IGA methods both from a theoretical (i.e. the possibility to write new formulations) and a numerical point of view. In the literature there are some very recent work in this direction for the formulation of ST-FE adapted to phase-field, see for instance [45,16].

The paper is structured as follows. The phase-field modeling in a space-time context is presented in Section 2. Section 3 describes the numerical implementation of the space-time phase-field model and we also presented, the staggered algorithm which is chosen here to deal with the non-convexity in the framework of the space-time dynamic phase-field model. Section 4 presents several benchmark problems both with ST-FE and ST-IGA models. Finally, Section 5 concludes the work

## 2. Space-time formulation

### 2.1. A space and time potential adapted to phase-field damage models

We consider the space-time cylinder $Q = \Omega \times [0, T]$, where $\Omega \subset \mathbb{R}^n$ is a closed domain such that $n$ is the number of space dimensions and $T$ is the final time of interest. The boundary of $\Omega$ is denoted by $\Gamma$. Let $\Gamma_u$ (for imposed displacement) and $\Gamma_\sigma$ (for imposed forces) denote the non-overlapping subregion of $\Gamma$ such that: $\Gamma = \Gamma_u \cup \Gamma_\sigma$, $\Gamma_u \cap \Gamma_\sigma = \emptyset$. We denote by $P = \Gamma \times [0, T]$ a space-time boundary (which can also be decomposed in $P_u$ and $P_\sigma$). In the context of isotropic small strain damage models, the fracture response of a solid is described by the displacement field of material point $x \in \Omega$ and by a local damage field which are defined as follows:

$$ \mathbf{u} : \begin{cases} \Omega \times [0, T] \to \mathbf{R}^r \\ (\mathbf{x}, t) \mapsto \mathbf{u}(\mathbf{x}, t) \end{cases} \quad d : \begin{cases} \Omega \times [0, T] \to [0, 1] \\ (\mathbf{x}, t) \mapsto d(\mathbf{x}, t) \end{cases} $$

where $r$ is the number of dimension of the kinematics. The strains are assumed to be small and the strain tensor is defined as:

$$ \epsilon(\mathbf{u}) = \nabla_s \mathbf{u} = \frac{\nabla \mathbf{u} + \nabla^T \mathbf{u}}{2} \quad (1) $$

$d(\mathbf{x}, t)$ locally defines the damage state at a material point $\mathbf{x}$, which is fully broken for $d = 1$ and unbroken for $d = 0$. In this work, we follow the standard description of the crack regularization with the phase field approach. This approach could be motivated by considering the following crack surface density function (see for instance [35] for more details):

$$ \gamma(d, \nabla d) = \frac{1}{2l} d^2 + \frac{l}{2} |\nabla d|^2 \quad (2) $$

where $l$ is a regularization parameter or a characteristic length scale. The minimization of the integral over the space domain $\Omega$ of this function leads to the classical phase field equation for the damage:

$$ d - l^2 \Delta d = 0 \quad \text{in } \Omega \quad \text{and} \quad (\nabla d)\mathbf{n} = 0 \quad \text{on} \quad \Gamma \quad (3) $$

for a 1D problem with the conditions $d(0) = 1$, $d(\pm\infty) = 0$ the solution of the previous equation is $d(x) = \exp(-|x|/l)$ and $l$ can be clearly interpreted as a length scale for the regularization of the crack. This solution also emphasizes that the gradient of the damage could be singular at some points. Some authors have proposed to use higher order theories such as the one proposed by [9] to obtain more regular solutions. The time variation of the crack surface density function gives the dissipated energy. Obtaining the phase field equation from a minimization principle motivates the development of variational formulations for cracks propagation and is the starting point of the following developments.

We propose to consider the following space-time potential (for the sake of simplicity time derivatives are denoted by a dot):

$$ \Pi = \int_Q k(\dot{\mathbf{u}}) + \rho \psi(\epsilon, d) - \varphi(\dot{d}, \nabla d) dQ + \int_Q \mathbf{u} f dQ + \int_P \mathbf{u} t dP \quad (4) $$

where $k(\dot{\mathbf{u}})$ is the kinetic energy, $\rho$ is the material density, $\psi$ is the Helmholtz free energy, $\varphi$ is a dissipation potential, $f$ and $t$ are the external loads applied on the space domain or its boundary. The kinetic energy is simply defined by

$$ k(\dot{\mathbf{u}}) = \frac{1}{2} \rho \dot{\mathbf{u}}^2 \quad (5) $$

For the free energy we use a decomposition into compression and tension (damage only affects the tension part):

$$ \psi(\epsilon, d) = g(d) \psi_0^+(\epsilon) + \psi_0^-(\epsilon) \quad (6) $$

where $g(d)$ is a degradation function and $\psi_0$ the non damaged free energy of the elastic media and the + and - exponent respectively indicate tension or compression.<fnref n="1" /> In this paper we have used the so-called spheric/deviatoric decoupling in the numerical examples. The splitting of the energy represents an important choice that must be made on the basis of experimental observations and material considerations. Interesting readers may refer to [31,49] and references therein for more information about this aspect of phase-field modeling. For the dissipation function we follow [35] and we adopt the following form (see also [41] for a discussion on dissipation potential within non-local approaches):

$$ \varphi(\dot{d}, \nabla d) = Y_d \dot{d} + \mathbf{Y}_{\nabla d} \nabla d + I_d \quad (7) $$

where $Y_d$ and $\mathbf{Y}_{\nabla d}$ are the thermodynamics forces associated to the damage and its gradient, $I_d$ is the indicator function:

$$ I_d = \begin{cases} 0 & \text{if } \dot{d} \ge 0 \\ \infty & \text{otherwise} \end{cases} \quad (8) $$

<footnote n="1">¹ In the literature there are several proposals for this decomposition, [35] propose to use the spectral decomposition of the strain while [5] propose to use the spheric/deviatoric decoupling. Each proposal is adapted for a given material/situation. Nevertheless this choice does not affect the space-time formulation.</footnote>
We also assume that the mass is locally conserved, which leads to the fact that $\rho = cte$ in a Lagrangian framework.

In a previous work on a space-time approach for (undamaged) elastodynamics (see [48]) we have shown that in the case of two fields (displacement and velocity) formulation, simple stabilization terms can be introduced that are independent of the material behavior. Therefore, in line with these ideas we propose to consider the following modified potential instead of the one at eq. (4):

$$ \bar{\Pi} = \int_Q \dot{k}_{HR}(\dot{\mathbf{u}}, \mathbf{v}) + \rho\psi(\epsilon, d) - \varphi(d, \nabla d) - p_{GAC}(\ddot{\mathbf{u}}, \dot{\mathbf{v}})dQ \\ + \int_Q \dot{\mathbf{u}}fdQ + \int_P \dot{\mathbf{u}}tdP \quad (9) $$

where $\mathbf{v}$ is the velocity field and $k_{HR}$ is the Hellinger-Reisner form of the kinetic energy:

$$ k_{HR}(\dot{\mathbf{u}}, \mathbf{v}) = \rho\dot{\mathbf{v}} - \frac{1}{2}\rho\dot{v}^2 \quad (10) $$

and $p_{GAC}$ is a least square stabilization term expressing the consistency condition on the acceleration ($\tau$ is a stabilization parameter with no dimension):

$$ p_{GAC}(\ddot{\mathbf{u}}, \dot{\mathbf{v}}) = \frac{1}{2}\tau\rho(\ddot{\mathbf{u}} - \dot{\mathbf{v}})^2 \quad (11) $$

Furthermore, as done by several authors, we consider a regularized version of the indicator function such that $I_d$ is approached by the following function:

$$ I_d^\gamma = \frac{\gamma}{p+1} \langle d \rangle^{p+1} \quad (12) $$

where $\gamma$ and $p$ are regularization parameters. The stationary conditions of the previously defined potential can now be studied. The variation of the potential can be written as (we assume a dead load case for which $f$ and $t$ are independent of the primary fields $\{\mathbf{u}, \mathbf{v}, d\}$):

$$ \begin{aligned} \delta \bar{\Pi} = & \int_Q \frac{\partial \dot{k}_{HR}}{\partial \dot{\mathbf{u}}} \delta \dot{\mathbf{u}} + \frac{\partial \dot{k}_{HR}}{\partial \mathbf{v}} \delta \mathbf{v} + \rho \frac{\partial \dot{\psi}}{\partial \epsilon} : \delta \epsilon + \rho \frac{\partial \dot{\psi}}{\partial d} \delta d - \frac{\partial p_{GAC}}{\partial \dot{\mathbf{u}}} \delta \dot{\mathbf{u}} \\ & - \frac{\partial p_{GAC}}{\partial \dot{\mathbf{v}}} \delta \dot{\mathbf{v}} - \frac{\partial \dot{\phi}}{\partial d} \delta d - \frac{\partial \dot{\phi}}{\partial \nabla d} \delta \nabla d + \frac{\partial I_d^\gamma}{\partial d} \delta d dQ \\ & + \int_Q \delta \dot{\mathbf{u}}fdQ + \int_P \delta \dot{\mathbf{u}}tdP \end{aligned} \quad (13) $$

where $\{\mathbf{u}, \mathbf{v}, d\} \in \{\mathcal{H}^u, \mathcal{H}^v, \mathcal{H}^d\}$ and their variations are such that $\delta \mathbf{u} \in \mathcal{H}^u_0$, $\delta \mathbf{v} \in \mathcal{H}^v_0$, $\delta d \in \mathcal{H}^d_0$ such as:

$$ \mathcal{H}^v = \{v \in \mathcal{H}^{0,1}(Q), v(x, t=0) = v_0(x)\} $$

$$ \mathcal{H}^u = \{u \in \mathcal{H}^{1,2}(Q), u = g(t) \text{ on } \Gamma, u(x, t=0) = u_0(x)\} $$

$$ H^d = \{w \in \mathcal{H}^{1,1}(Q), w(x, t=0) = d_0(x)\} $$

$$ H_0^v = \{v \in \mathcal{H}^{0,1}(Q), v = 0 \text{ on } \Gamma, v(x, t=0) = 0, v(x, t=T) = 0\} $$

$$ H_0^u = \{u \in \mathcal{H}^{1,2}(Q), u = 0 \text{ on } \Gamma, u(x, t=0) = 0, u(x, t=T) = 0\} $$

$$ H_0^d = \{w \in \mathcal{H}^{1,1}(Q), w(x, t=0) = 0, w(x, t=T) = 0\} $$

$\mathcal{H}^{l,k}(Q)$ is the Sobolev spaces of functions defined in the space-time cylinder such as:

$$ \{u \in \mathcal{L}^2(Q) : \partial_x^\beta u \in \mathcal{L}^2(Q) \forall \beta \text{ with } 0 \le |\beta| \le l, \partial_i^l u \in \mathcal{L}^2(Q), i = 0, \dots, k\} $$

where $\mathcal{L}^2(Q)$ denotes the space of square-integrable functions, $\beta = (\beta_1, \dots, \beta_d)$ is a multi-index with non-negative integers, $|\beta| = \beta_1 + \dots + \beta_d$, $\partial_x^\beta u := \partial_{x_1}^{\beta_1} \partial_{x_2}^{\beta_2} \dots \partial_{x_d}^{\beta_d}$ with $\partial_{x_i}^{\beta_i} := \frac{\partial}{\partial x_i^{\beta_i}}$ and $\partial_i^l u = \frac{\partial^l u}{\partial t^l}$.

Using integration by part on the time and the homogeneous conditions on the time boundary for the fields variation, one can obtain the following set of stationary conditions of the previously defined potential:

$$ \begin{aligned} \delta_u \bar{\Pi} &\equiv \int_Q \rho \dot{\mathbf{v}} \delta \dot{\mathbf{u}} + \sigma : \delta \dot{\epsilon} dQ - \int_Q \delta \dot{\mathbf{u}}fdQ - \int_P \delta \dot{\mathbf{u}}tdP \\ &\quad + \int_Q \rho \tau(\ddot{\mathbf{u}} - \dot{\mathbf{v}}) \delta \ddot{\mathbf{u}} dQ = 0 \end{aligned} \quad (15) \\ \delta_v \bar{\Pi} &\equiv \int_Q \rho(\mathbf{v} - \dot{\mathbf{u}}) \delta \dot{\mathbf{v}} dQ + \int_Q \rho \tau(\ddot{\mathbf{u}} - \dot{\mathbf{v}}) \delta \ddot{\mathbf{v}} dQ = 0 \\ \delta_d \bar{\Pi} &\equiv \int_Q (\rho g'(d) \psi_0^+ (\epsilon) + Y_d - \gamma \langle d \rangle^p) \delta d\dot{\mathbf{d}} + \mathbf{Y}_{\nabla d} \delta \nabla d dQ = 0 \end{aligned} $$

The local form of the previous equations, with $\tau = 0$ and $\gamma = 0$, corresponds to the conservation of the linear momentum, to the definition of the velocity and to the phase field equation, which are defined as (we assume that $\mathbf{Y}_{\nabla d} \cdot \mathbf{n} = 0$ and we use the condition $\sigma \cdot \mathbf{n} = \mathbf{0}$ on $\Gamma_u$):

$$ \left\{ \begin{aligned} \operatorname{div} \boldsymbol{\sigma} + \mathbf{f} &= \rho \dot{\mathbf{v}} \\ \rho(\dot{\mathbf{u}} - \mathbf{v}) &= 0 \\ \operatorname{div} \mathbf{Y}_{\nabla d} &= \rho g'(d) \psi_0^+ + Y_d \end{aligned} \right. \quad (16) $$

where $\boldsymbol{\sigma} = g(d) \partial \psi_0^+ / \partial \epsilon + \partial \psi_0^- / \partial \epsilon$ and $\boldsymbol{\sigma} \cdot \mathbf{n} = \mathbf{t}$ on $\Gamma_\sigma$.

Unfortunately the penalization function often leads to numerical difficulties, the stabilization parameter is tricky to defined and problem dependent (see [21] for a discussion about it). Furthermore this regularization does not guaranty that the irreversibility condition is fulfilled everywhere at anytime in $Q$. An alternative formulation can be defined by following the proposal of Miehe with the introduction of the history function, $\mathcal{H}$, that stores the maximum value of the argument over the time (see [35]). In this approach, the system of equations at eq. (15) is then replaced by:

$$ \begin{aligned} \delta_u \bar{\Pi} &\equiv \int_Q \rho \dot{\mathbf{v}} \delta \dot{\mathbf{u}} + \sigma : \delta \dot{\epsilon} dQ - \int_Q \delta \dot{\mathbf{u}}fdQ - \int_P \delta \dot{\mathbf{u}}tdP \\ &\quad + \int_Q \rho \tau(\ddot{\mathbf{u}} - \dot{\mathbf{v}}) \delta \ddot{\mathbf{u}} dQ = 0 \end{aligned} \quad (17) \\ \delta_v \bar{\Pi} &\equiv \int_Q \rho(\mathbf{v} - \dot{\mathbf{u}}) \delta \dot{\mathbf{v}} dQ + \int_Q \rho \tau(\ddot{\mathbf{u}} - \dot{\mathbf{v}}) \delta \ddot{\mathbf{v}} dQ = 0 \\ \delta_d \bar{\Pi} &\equiv \int_Q (\rho g'(d) \mathcal{H}(\psi_0^+(\epsilon)) + Y_d) \delta d\dot{\mathbf{d}} + \mathbf{Y}_{\nabla d} \delta \nabla d dQ = 0 \end{aligned} $$

The local system of equations defined at eq. (16) is obviously unchanged except for the last one which is replaced by:

$$ \operatorname{div} \mathbf{Y}_{\nabla d} = \rho g'(d) \mathcal{H}(\psi_0^+(\epsilon)) + Y_d \quad (18) $$

As remarked by many authors formulations at eqs. (15) and (17) could not be proven to be equivalent in general. However the history function approach does not require the introduction of a penalization parameter and allows to satisfy the fundamental thermodynamic requirement $\dot{d} \ge 0$. Unfortunately, it has been shown by some authors that the use of the history function with the traction/compression splitting of the energy could result in the loss of the $\Gamma$-convergence, e.g. [32]. Nevertheless we decided to use this formulation as it represents the simplest alternative for guaranteeing respect of the irreversibility condition.

To be consistent with the form of the previously defined crack surface density function, at eq. (2), we assume that the thermodynamics forces associated to $d$ and $\nabla d$ are chosen as:

$$ Y_d = \frac{G_c}{l} d, \quad \mathbf{Y}_{\nabla d} = G_c l \nabla d \quad (19) $$
<figure><img src="image_3.png" /><figcaption>Fig. 1. Time decomposition domain in 4 slabs.</figcaption></figure>

where G_c is the material toughness which is here assumed to be constant (not rate dependent for instance). Furthermore, we take the standard quadratic form for the degradation function g(d):

$$g(d) = (1 - d)^2 \qquad (20)$$

these choices correspond to the so-called AT2 model (e.g. [3,11]).

## 2.2. A time-discontinuous Galerkin formulation

The previously defined space-time potential, $\bar{\Pi}$, is expressed in a time continuous form. However, both phase field brittle fracture problems and elastodynamics problems could lead to highly singular solution in time. Furthermore, it is not trivial to account of the history function, $H$, in an explicit manner in a discrete version of the time continuous formulation. For all these reasons, we propose to adopt a time-discontinuous formulation.

We consider a finite discretization, denoted $Q_h$, of the space-time cylinder $Q$. This discretized domain is partitioned in the time direction into N subdomains $Q_h^i$, or time slabs, such that $t_1 = 0 < t_2 < \dots < t_i < T_{N+1} = T$. $Q_h^i$ is defined such that: $Q_h = \cup_i Q_h^i$, see Fig. 1.

To shorten the notations, we denote the $i^{th}$ time slab $Q_h^i$ by $Q_i$, such that $Q_i = \Omega \times ]t_i, t_{i+1}[$ and $P_i = \Gamma \times ]t_i, t_{i+1}[$. We assume that the fields can be discontinuous over the time slabs, such that the temporal jumps are defined by:

$$ \begin{aligned} [[ \mathbf{u}(\mathbf{x}, t_i) ]] &= \mathbf{u}(\mathbf{x}, t_i^+) - \mathbf{u}(\mathbf{x}, t_i^-) \\ [[ \mathbf{v}(\mathbf{x}, t_i) ]] &= \mathbf{v}(\mathbf{x}, t_i^+) - \mathbf{v}(\mathbf{x}, t_i^-) \\ [[ d(\mathbf{x}, t_i) ]] &= d(\mathbf{x}, t_i^+) - d(\mathbf{x}, t_i^-) \end{aligned} \qquad (21) $$

where $\mathbf{u}(\mathbf{x}, t_i^+), \mathbf{v}(\mathbf{x}, t_i^+), d(\mathbf{x}, t_i^+)$, denote the values, at time $t_i$, on the time slab $Q_i$ and $\mathbf{u}(\mathbf{x}, t_i^-), \mathbf{v}(\mathbf{x}, t_i^-), d(\mathbf{x}, t_i^-)$, denote the values, at time $t_i$, on the time slab $Q_{i-1}$.

We propose to consider here the following modified space-time potential:

$$ \begin{aligned} \Pi_{T DG} = & \sum_{i=1}^{N} \int_{Q_i} (\dot{k}_{HR}(\dot{\mathbf{u}}, \mathbf{v}) + \rho \dot{\psi}(\varepsilon, d) - \varphi(\dot{d}, \bar{\nabla} d) - p_{GAC}(\ddot{\mathbf{u}}, \dot{\mathbf{v}})) dQ \\ & + \sum_{i=2}^{N} \int_{B_i} [[k_{HR}(\dot{\mathbf{u}}, \mathbf{v})]] dB + \sum_{i=2}^{N} \int_{B_i} [[\rho \psi(\varepsilon, d)]] dB \\ & + \frac{1}{2} \int_{B_1} \alpha_u (u(\mathbf{x}, t_i^+) - u_0(\mathbf{x}))^2 dB + \frac{1}{2} \int_{B_1} \alpha_v (v(\mathbf{x}, t_i^+) - v_0(\mathbf{x}))^2 dB \\ & + \frac{1}{2} \int_{B_1} \alpha_d (d(\mathbf{x}, t_i^+) - d_0(\mathbf{x}))^2 dB + \sum_{i=1}^{N} \int_{Q_i} \dot{\mathbf{u}} f dQ + \sum_{i=1}^{N} \int_{P_i} \dot{\mathbf{u}} t dP \end{aligned} \qquad (22) $$

where $B_i$ stands for the space domain at time $t_i^+$, which is here perfectly equivalent to the space domain: $B_i = \Omega(t_i^+) = \Omega$. It can be seen that

compared to the previous continuous space-time potential we have supplementary terms that express the initial conditions in the weak sense and the time continuity of the kinetic and strain energies. The coefficients $\alpha_u, \alpha_v, \alpha_d$ are introduced to obtain the correct dimensionality of the initial condition terms. These coefficients must be chosen with care and can have an impact on the solution if not defined correctly. One can choose $\alpha_v = \rho$ and relate $\alpha_u$ and $\alpha_d$ to the elastic modulus and time increment (and possibly to a characteristic length, such as for instance $\alpha_d \approx E/\Delta t$, $\alpha_u \approx E/L^2\Delta t$ where L is some characteristic length and E the Young modulus). In the following, after some preliminary numerical tests, we have taken these coefficients to be equal to one (see Appendix A for more details).

As previously, computing the variation of $\Pi_{T DG}$ and regrouping terms it can be obtained the stationarity conditions of the space-time potential. It can be shown that these conditions can be expressed separately on each time slab, and the problem becomes incremental (the solution of the previous time slab or the initial conditions are required to compute the next one).

On the first time slab $Q_1$ we have:

$$ \begin{aligned} \delta_u \Pi_{T DG} &\equiv \int_{Q_1} \rho \dot{\mathbf{v}} \delta \dot{\mathbf{u}} + \sigma : \delta \dot{\varepsilon} dQ - \int_{Q_1} \delta \dot{\mathbf{u}} f dQ - \int_{P_1} \delta \dot{\mathbf{u}} t dP \\ &\quad + \underbrace{\int_{Q_1} \rho \tau (\ddot{\mathbf{u}} - \dot{\mathbf{v}}) \delta \ddot{\mathbf{u}} dQ}_{\text{Initial BC Weak form}} - \underbrace{\int_{B_1} (\mathbf{u}(\mathbf{x}, t_i^+) - \mathbf{u}_0(\mathbf{x})) \delta \dot{\mathbf{u}} dB}_{\text{= 0}} \\ \delta_v \Pi_{T DG} &\equiv \int_{Q_1} \rho (\mathbf{v} - \dot{\mathbf{u}}) \delta \dot{\mathbf{v}} dQ + \int_{Q_1} \rho \tau (\ddot{\mathbf{u}} - \dot{\mathbf{v}}) \delta \dot{\mathbf{v}} dQ \\ &\quad + \underbrace{\int_{Q_1} (\rho g'(d) H(\psi_0^+(\varepsilon)) + Y_d) \delta d + \mathbf{Y}_{\bar{\nabla} d} \delta \bar{\nabla} d dQ}_{\text{Initial BC Weak form}} - \underbrace{\int_{B_1} (\mathbf{v}(\mathbf{x}, t_i^+) - \mathbf{v}_0(\mathbf{x})) \delta \dot{\mathbf{v}} dB}_{\text{= 0}} \end{aligned} \qquad (23) $$

$$ \begin{aligned} \delta_d \Pi_{T DG} &\equiv \int_{Q_1} (\rho g'(d) H(\psi_0^+(\varepsilon)) + Y_d) \delta d + \mathbf{Y}_{\bar{\nabla} d} \delta \bar{\nabla} d dQ \\ &\quad + \underbrace{\int_{B_1} (d(\mathbf{x}, t_i^+) - d_0(\mathbf{x})) \delta d^+ d B}_{\text{Initial BC Weak form}} - \underbrace{\delta d^+ d B}_{\text{= 0}} \end{aligned} $$

and on $Q_i$ with $2 \le i \le N$, we have:

$$ \begin{aligned} \delta_u \Pi_{T DG} &\equiv \int_{Q_i} \rho \dot{\mathbf{v}} \delta \dot{\mathbf{u}} + \sigma : \delta \dot{\varepsilon} dQ - \int_{Q_i} \delta \dot{\mathbf{u}} f dQ - \int_{P_i} \delta \dot{\mathbf{u}} t dP \\ &\quad + \underbrace{\int_{Q_i} \rho \tau (\ddot{\mathbf{u}} - \dot{\mathbf{v}}) \delta \ddot{\mathbf{u}} dQ}_{\text{Galerkin discontinuous Terms}} - \underbrace{\int_{B_i} [[\sigma]] : \delta \dot{\varepsilon}^+ d B}_{\text{= 0}} \end{aligned} $$
$$
\begin{equation}
\begin{aligned}
\delta_{\mathbf{v}} \Pi_{TDG} \equiv & \int_{Q_i} \rho (\mathbf{v} - \dot{\mathbf{u}}) \delta\dot{\mathbf{v}} dQ + \int_{Q_i} \rho \tau (\dot{\mathbf{u}} - \dot{\mathbf{v}}) \delta\dot{\mathbf{v}} dQ \\
& \qquad \text{Galerkin discontinuous Terms} \\
& - \underbrace{\left[ \left[ \rho \mathbf{v} \right] \delta\dot{\mathbf{v}}^+ d\mathbf{B} \right]}_{\mathbf{B}_i} = 0
\end{aligned}
\tag{24}
\end{equation}
$$

$$
- \underbrace{\left[ \left[ \rho \mathbf{v} \right] \delta \dot{\mathbf{v}}^+ d \mathbf{B} \right]}_{\mathbf{B}_i} = 0
$$

$$
\delta_d \Pi_{TDG} \equiv \int_{Q_i} \left[ (\rho g'(d) H(\psi_0^+(\epsilon)) + Y_d) \delta d + \mathbf{Y}_{\nabla d} \delta \bar{\nabla} d \right] dQ
$$

$$
- \underbrace{\left[ \left[ \left( \frac{\partial \psi}{\partial d} \right) \right] \delta d^+ d B \right]}_{B_i} = 0
$$

It can be remarked that we neglected the continuity terms linked to **u**
when computing the variation of the kinetic energy: variations of discon-
tinuity terms at time *t*ᵢ are computed by frozen the time derivatives at *t*ᵢ.
In the previous expressions we have used the notation δ**u**⁺, δ**v**⁺, δ**d**⁺ to
denote the test functions associated to each field at time *t*ᵢ on the time
slab *Q*ᵢ. It can be remarked that this time discontinuous formulation
does not require the introduction of supplementary numerical parame-
ters and the time continuity terms are correctly dimensioned for each
equation.

3. Space-time finite-element (ST-FE) or space-time IGA (ST-IGA)

3.1. Space-time discretization

The finite discretization of the time slab $Q_i$ is realized with either
Lagrange Finite-Elements or Isogeometric Analysis (NURBS based ele-
ments). The restriction of $Q_i$ to a single element is denoted by $Q_e^h$ and
the primary fields are approximated with the same basis functions that
are used to describe the geometry (i.e. Lagrange polynomials or NURBS):

$$
\begin{align}
\mathbf{u}_e^h(\mathbf{x}, t) &= \mathbf{N}_u(\mathbf{x}, t) \mathbf{d}_e^u \notag \\
\mathbf{v}_e^h(\mathbf{x}, t) &= \mathbf{N}_v(\mathbf{x}, t) \mathbf{d}_e^v \notag \\
d_e^h(\mathbf{x}, t) &= \mathbf{N}_d(\mathbf{x}, t) \mathbf{d}_e^d \tag{25}
\end{align}
$$

where N_u, N_v, N_d denote the matrices<fnref n="2" /> made of approximation functions
build from an element and d_e^u, d_e^v, d_e^d the restriction of the global vector
of freedom, d, to an element for each field. The matrices are build such
that the vector/matrices products give the correct approximation fields
depending on the type of kinematic considered. It has to be remarked
that we do not necessarily assume the approximations to be the same
for each field. Therefore, we can choose to use different approximation
orders (or continuity orders for IGA) for each field. Furthermore, we do
not assume a particular restriction regarding the meshing of the time
slab, therefore we can consider cases with more than one element in
the time direction which can be viewed as a mixture of time contin-
uous and time discontinuous approaches. This may have an interest if
one considers space-time parallelization (which we use in our numerical
implementation), but we won't discuss this issue in this paper.

From the previous approximations, one can define time derivatives
of the primary fields. These derivatives on an element are given by:

$$
\begin{align}
\dot{\mathbf{u}}_e^h(\mathbf{x}, t) &= \mathbf{B}_u^t(\mathbf{x}, t) \mathbf{d}_e^u, & \ddot{\mathbf{u}}_e^h(\mathbf{x}, t) &= \mathbf{B}_u^{tt}(\mathbf{x}, t) \mathbf{d}_e^u \nonumber \\
\dot{\mathbf{v}}_e^h(\mathbf{x}, t) &= \mathbf{B}_v^t(\mathbf{x}, t) \mathbf{d}_e^v, & & \nonumber \\
d_e^h(\mathbf{x}, t) &= \mathbf{B}_d^t(\mathbf{x}, t) \mathbf{d}_e^d \tag{26}
\end{align}
$$

where **B'** and **B''** are the matrices formed from the time derivatives
of the previously defined shape functions. It has to be noted that due
to the presence of the second time derivatives we have to use at least
functions of order 2 for the kinematic field. In a similar manner we can
define space gradients or space-time gradients

$$
\begin{align}
\epsilon_e^h(\mathbf{x}, t) &= \mathbf{B}_u^x \mathbf{d}_e^u, \qquad \dot{\epsilon}_e^h(\mathbf{x}, t) = \mathbf{B}_u^{xt} \mathbf{d}_e^u \nonumber \\
\nabla_x d_e^h(\mathbf{x}, t) &= \mathbf{B}_d^x \mathbf{d}_e^d, \qquad \dot{\nabla}_x d_e^h(\mathbf{x}, t) = \mathbf{B}_d^{xt} \mathbf{d}_e^d \tag{27}
\end{align}
$$

As for the time derivatives, the matrices B^x and B^t denote the space or
space and time derivatives of the shape functions N. For the same rea-
sons as previously mentioned for the kinematic field, the approximation
of the damage field should be at least of order 1 in space and time due
to the occurrence of space and time derivatives.

By incorporating the elemental approximations defined at eqs
(25)-(27) into eqs. (23) and (24), we obtain the following discrete sys-
tem of equations, for the first time slab:

$$
\left\{
\begin{aligned}
\mathbf{r}_u &= A_e \left\{ \int_{Q_{e1}^h} (\rho \mathbf{B}_u^T \dot{\mathbf{v}}_e^h + \mathbf{B}_u^{xT} \sigma) dQ - \int_{Q_{e1}^h} \mathbf{B}_u^T \mathbf{f} dQ \right\} - A_e \int_{P_{e1}^h} \mathbf{B}_u^{T} \mathbf{t} dP \\
&\quad + A_e \int_{Q_{e1}^h} \rho \tau \mathbf{B}_u^{tt} (\dot{\mathbf{u}}_e^h - \dot{\mathbf{v}}_e^h) dQ - A_e \int_{B_{e1}^h} \mathbf{N}_u^{+T} (\mathbf{u}^+ - \mathbf{u}_0) dB = 0 \\[1em]
\mathbf{r}_v &= A_e \left\{ \int_{Q_{e1}^h} \mathbf{B}_v^T (\dot{\mathbf{u}}_e^h - \mathbf{v}_e^h) + \rho \tau \mathbf{B}_v^T (\ddot{\mathbf{u}}_e^h - \dot{\mathbf{v}}_e^h) \right\} dQ \\
&\quad - A_e \int_{B_{e1}^h} \mathbf{N}_v^{+T} (\mathbf{v}^+ - \mathbf{v}_0) dB = 0 \\[1em]
\mathbf{r}_d &= A_e \left\{ \int_{Q_{e1}^h} \mathbf{B}_d^T (\rho g'(d_e^h) H(\psi_0^+(\epsilon_e^h)) + Y_d^h) dQ \right. \\
&\quad \left. + A_e \int_{Q_{e1}^h} \mathbf{B}_d^{xT} \mathbf{Y}_{\nabla d} dQ - A_e \int_{B_{e1}^h} \mathbf{N}_d^{+T} (d^+ - d_0) dB = 0 \right\}
\end{aligned}
\right.
\tag{28}
$$

and the following ones for the other time slabs:

$$
\left\{
\begin{aligned}
\mathbf{r}_u &= A_e \left\{ \int_{Q_{e1}^h} (\rho \mathbf{B}_u^T \dot{\mathbf{v}}_e^h + \mathbf{B}_u^{xT} \sigma) dQ - \int_{Q_{e1}^h} \mathbf{B}_u^T \mathbf{f} dQ \right\} - A_e \int_{P_{e1}^h} \mathbf{B}_u^T \mathbf{t} dP \\
&\quad + A_e \int_{Q_{e1}^h} \rho \tau \mathbf{B}_u^{tt} (\dot{\mathbf{u}}_e^h - \dot{\mathbf{v}}_e^h) dQ - A_e \int_{B_{e1}^h} \mathbf{B}_u^{x+T} [[\sigma]] dB = 0 \\[1em]
\mathbf{r}_v &= A_e \left\{ \int_{Q_{e1}^h} \mathbf{B}_v^T (\dot{\mathbf{u}}_e^h - \mathbf{v}_e^h) + \rho \tau \mathbf{B}_v^T (\ddot{\mathbf{u}}_e^h - \dot{\mathbf{v}}_e^h) \right\} dQ \\
&\quad - A_e \int_{B_{e1}^h} \mathbf{N}_v^{+T} [[\rho \mathbf{v}_e^h]] dB = 0 \\[1em]
\mathbf{r}_d &= A_e \left\{ \int_{Q_{e1}^h} \mathbf{B}_d^T (\rho g'(d_e^h) H(\psi_0^+(\epsilon_e^h)) + Y_d^h) dQ + A_e \int_{Q_{e1}^h} \mathbf{B}_d^{xT} \mathbf{Y}_{\nabla d} dQ \right. \\
&\quad \left. - A_e \int_{B_{e1}^h} \mathbf{N}_d^{+T} \left[ \rho \frac{\partial \psi}{\partial d} \right] dB = 0 \right\}
\end{aligned}
\right.
\tag{29}
$$

<footnote n="2">² For the sake of simplicity we do not use the *e* exponent for the approximation functions. For FE, these functions are only defined elements by elements, for IGA these functions can be defined on the all parametric space. For the sake of brevity, the notion of reference element is not presented here and assumed as familiar to the readers.</footnote>
The previous system of equations is non-linear for each time slab. More precisely, $r_u$ is non-linear due to the decoupling of the energy into compression and traction part, $r_v$ is linear and $r_d$ is also non-linear in general (but linear regarding the damage field with the choice made in previous section for $Y_d$, $Y_{\nabla d}$, $g(d)$ if the displacement is kept fixed). Furthermore, it can be shown that the space-time potential $\Pi_{TDG}$ is non convex in general (but convex separately for $\{\mathbf{u}, \mathbf{v}\}$ and for $d$). For the numerical resolution we have adopted an alternate minimization strategy that is described here after. This leads to the fact that we used two schemes: a standard Newton iterative scheme for $\{\mathbf{u}, \mathbf{v}\}$ and a linear resolution for $d$. One iteration of the Newton scheme can be summarized as:

$$ \left\{ \begin{matrix} \mathbf{d}^u \\ \mathbf{d}^v \end{matrix} \right\} \left\{ \begin{matrix} \mathbf{d}^u \\ \mathbf{d}^v \end{matrix} \right\} - \left[ \begin{matrix} K_{uu} & K_{uv} \\ K_{vu} & K_{vv} \end{matrix} \right]^{-1} \left\{ \begin{matrix} \tilde{\mathbf{r}}_u \\ \tilde{\mathbf{r}}_v \end{matrix} \right\} \quad (30) $$

the tilde over the residuals designates original quantities defined at eqs. (28) or (29) but modified by boundary conditions of the problem (the tangent operator is only assembled for non constrained degrees of freedom, see [24]). The solution of the linear system for $d$ is simply:

$$ \{ \mathbf{d}^d \} \left\{ \begin{matrix} \mathbf{d}^d \end{matrix} \right\} - \left[ K_{dd} \right]^{-1} \{ \tilde{\mathbf{r}}^d \} \quad (31) $$

Due to this staggered resolution we have observed that the time continuity term for damage equation is inefficient or more precisely incomplete (because we are missing a continuity term for the gradient of $d$). Therefore, by analogy with what can be done with standard materials (for which the thermodynamics forces derive from the strain energy) we propose to replace the term $[[\rho \partial \psi / \partial d]]$ by a continuity condition on the thermodynamics forces $Y_d$ and $Y_{\nabla d}$. The last equation of the system defined at eq. (29) is then replaced by:

$$ \mathbf{r}_d = A_e \left\{ \begin{matrix} \int_{Q_{e_i}^h} \mathbf{B}_d^{xT} (\rho g'(d_e^h) H(\psi_0^+(\epsilon_e^h)) + Y_d^h) \\ -A_e \int_{B_{e_i}^h} \mathbf{N}_d^{+T} [[Y_d]] + \mathbf{B}_d^{+T} [[\mathbf{Y}_{\nabla d}]] \end{matrix} \right\} dQ + A_e \int_{Q_{e_i}^h} \mathbf{B}_d^{xT} \mathbf{Y}_{\nabla d} dQ \\ -A_e \int_{B_{e_i}^h} \mathbf{N}_d^{+T} [[Y_d]] + \mathbf{B}_d^{+T} [[\mathbf{Y}_{\nabla d}]] dB = 0 \quad (32) $$

The proposed formulation can be understood as a force time-continuity formulation (in the weak sense) and the continuity conditions for $\mathbf{u}$, $d$ and $\nabla d$ are similar. The terms of the tangent matrix<fnref n="3" /> are defined from (for all time slabs except the first one for which the continuity terms are slightly different):

$$ \begin{align} K_{uu} &= A_e \int_{Q_e^h} \mathbf{B}_u^{xT} \mathbf{C}_{uu}^d \mathbf{B}_u^x dQ + A_e \int_{Q_e^h} \rho \tau \mathbf{B}_u^{tt} \mathbf{B}_u^{tt} dQ \nonumber \\ &\quad - A_e \int_{B_e^h} \mathbf{B}_u^{x+T} \mathbf{C}_{TDG} \mathbf{B}_u^{x+} dB \nonumber \\ K_{uv} &= A_e \int_{Q_e^h} \rho \mathbf{B}_u^T \mathbf{B}_v^t dQ - A_e \int_{Q_e^h} \rho \tau \mathbf{B}_u^{tt} \mathbf{B}_v^t dQ \nonumber \\ K_{vu} &= A_e \int_{Q_e^h} \rho \mathbf{B}_v^T \mathbf{B}_u^t dQ + A_e \int_{Q_e^h} \rho \tau \mathbf{B}_v^{tt} \mathbf{B}_u^{tt} dQ \tag{33} \\ K_{vv} &= A_e \int_{Q_e^h} -\rho \mathbf{B}_v^T \mathbf{N}_v dQ - A_e \int_{Q_e^h} \rho \tau \mathbf{B}_v^{tt} \mathbf{B}_v^t dQ - A_e \int_{B_e^h} \rho \mathbf{N}_v^{+T} \mathbf{N}_v^+ dB \nonumber \end{align} $$

$$ K_{dd} = A_e \int_{Q_e^h} \mathbf{B}_d^{xT} C_{dd} \mathbf{N}_d dQ + A_e \int_{Q_e^h} \mathbf{B}_d^{xT} C_{\nabla d \nabla d} \mathbf{B}_d^x dQ \\ - A_e \int_{B_e^h} \mathbf{N}_d^{+T} \frac{\partial Y_d}{\partial d} + \mathbf{N}_d^+ dB - A_e \int_{B_e^h} \mathbf{B}_d^{xT} \frac{\partial \mathbf{Y}_{\nabla d}}{\partial \nabla d} + \mathbf{B}_d^x dQ $$

Using the definition at eqs. (6), (20) and (19), and the proposal of [5] for the decoupling into compression and tension part of the strain energy:

$$ \rho \psi_0^+(\epsilon) = \mu(\epsilon^D : \epsilon^D) + \frac{k}{2} (\langle tr(\epsilon) \rangle_+)^2, \quad \rho \psi_0^-(\epsilon) = \frac{k}{2} (\langle tr(\epsilon) \rangle_-)^2 \quad (34) $$

where $\epsilon^D$ denotes the deviatoric part of the strain, $\mu, k$ are Lamé parameters (shear and bulk modulus), the terms $C_{uu}$, $C_{dd}$, $C_{\nabla d \nabla d}$ and $\partial Y_d / \partial d$, $\partial \mathbf{Y}_{\nabla d} / \partial \nabla d$ are given by:

$$ \begin{align} C_{uu} &= (1-d)^2 (2\mathbf{K} + k\mathbf{I} \otimes \mathbf{I}h_+(tr(\epsilon))) + k\mathbf{I} \otimes \mathbf{I}h_-(tr(\epsilon)) \nonumber \\ C_{dd} &= \mathcal{H}(\psi_0^+(\epsilon)) + \frac{G_c}{l} \nonumber \\ C_{\nabla d \nabla d} &= G_c \mathbf{I} \nonumber \\ \frac{\partial Y_d}{\partial d} &= \frac{G_c}{l} \nonumber \\ \frac{\partial \mathbf{Y}_{\nabla d}}{\partial \nabla d} &= G_c l \tag{35} \end{align} $$

with $\mathbf{K}$ the deviatoric projector, $\mathbf{I}$ the second order identity tensor, $h_+(tr(\epsilon))$ and $h_-(tr(\epsilon))$ the Heaviside functions<fnref n="4" /> related to the sign of $tr(\epsilon)$.

It can be remarked that from the point of view of the elemental formulation there are no differences between FE and IGA except in the fact that we can control the order of continuity (between elements) with IGA. In this paper we only use quadratic space-time approximations for the 3 primary fields, $\{\mathbf{u}, \mathbf{v}, d\}$, both for ST-FE ($C_0$ continuity between elements) and ST-IGA ($C_1$ continuity).

## 3.2. Alternate minimization algorithm

As already mentioned, the space-time potential defined in eq. (22) is non-convex in general with respect to its three unknowns. Non-convex problems can pose challenges in optimization and numerical computations, as finding the global minimum of such functional can be difficult. This non-convexity problem can be tackled numerically, by capitalizing on the convexity of each equation when one variable is held constant. This fundamental principle forms the core of the alternate algorithm designed to address the issue (obviously with no guaranty of global convergence). By utilizing this property, the algorithm is constructed to effectively navigate the non-convex nature of the problem, providing a potential solution approach. This algorithm has been suggested and used by numerous authors (e.g. [34,36,38,5,12,21,2]).

We have implemented this algorithm and the formulation presented above in in-house code (which can be used both with IGA or FE). We first solve the dynamic mechanical part iteratively on a time slab with the standard Newton-Raphson algorithm. During this resolution the values of the degrees of freedom for the damage are hold constant to their previous converged values. Once the convergence of the Newton-Raphson is achieved for a given precision related to a given norm (we use the infinite norm) we freeze displacement and velocity and we solve the linear problem for the damage. This process is repeated in the sense of a fixed point algorithm, alternatively both problems are resolved until we reach the convergence of the fixed point algorithm or a maximum number of iterations. If the convergence is reached we just pass to the next time

<footnote n="3">³ For the sake of simplicity, we denote by a + exponent the terms evaluated at time $t_i$ from the time slab $Q_i$.</footnote>

<footnote n="4">⁴ $h_+(tr(\epsilon)) = \{1 \text{ if } tr(\epsilon) > 0 \text{ and } 0 \text{ elsewhere}\}$, $h_-(tr(\epsilon)) = \{1 \text{ if } tr(\epsilon) < 0 \text{ and } 0 \text{ elsewhere}\}$.</footnote>
<figure><img src="image_4.png" /><figcaption>Fig. 2. ST-Alternate algorithm.</figcaption></figure>

slab. In the opposite case (non-convergence of the fixed point algorithm) we retry the current time slab but we previously adapt the mesh in time by using a simple contractile homotopy operation that affect only the size of the mesh in the time direction (so the mesh is reduced in time but the number of elements or nodes is unaffected). This process is repeat until the mesh size in time is not below a given value. We have also (heuristic) criteria to (re)extend the mesh size in time if needed. This algorithm is schematized at Fig. 2. Within this algorithm, the history function consists in finding and storing the maximum value of the strain energy (traction part) over all the previous converged time slab at a given position in space:

$$ \mathcal{H}(\psi_0^+(\epsilon_e^h(\mathbf{x},t))) = \max_{t<t_j} \psi_0^+(\epsilon_e^h(\mathbf{x},t)) \quad (36) $$

The previous values of the maximum of the strain energy are updated once we have converged for the current time slab. This is therefore very similar to the proposal of [35] but not strictly identical as the history evolves here time slab by time slab. Therefore the history at time $t_{i+1}$ depends both on a past (fixed) history that comes from the previous time slab (for $t \le t_i^-$) and an evolving history that comes from the current time slab (i.e. for $t_i^+ < t < t_{i+1}$).

There are many convergence criteria for the fixed point algorithm in the literature as this strategy can be inefficient from the numerical point of view and obviously there are no guarantees to reach convergence in a fixed point process. Miehe et al. [34,35] introduced the so-called “one pass criterion”, which is taken up by [36,33]. In this criterion, the coupling between damage and mechanics is treated sequentially, and for each loading increment, no global convergence is verified. In some cases, this scheme gives inaccurate and often erroneous results. In [20,52], the residual-based criterion is used. The main concept is to check both $\mathbf{r}^u$, $\mathbf{r}^v$ and $\mathbf{r}^d$ when the displacement, velocity and damage fields for the current state of loading are obtained. A further criterion used is the energy [2,20]. This criterion verifies that for a given load, between each iteration of the alternate resolution process, the variation in energy must be less than a threshold value. The damage criterion [5,11,12], ensures that the change in the damage field between two

<figure><img src="image_5.png" /><figcaption>Fig. 3. Bar impact.</figcaption></figure>

consecutive iterations, under a given loading state, remains below a pre-defined threshold value. Nevertheless, the alternate algorithm is robust and can handle snaps back instabilities. To improve convergence rate of this algorithm, Farrell et al. [17] propose another perspective which is to treat it as a nonlinear Gauss-Seidel iterative scheme. In addition, they recommend incorporating an over-relaxed parameter to improve iteration speed. In fact, the authors propose a mixed algorithm made up of alternate minimization, which seeks to get as close as possible to the pool of possible solutions, and gave way to the monolithic resolution algorithm. In such instances, the challenge lies in managing the transition between the two algorithms, which often necessitates an ad hoc adjustment. Nevertheless, we will not discuss these questions in the following section, we simply use the criteria on the damage field between two fixed point iterations as explained previously. In the following examples, we have used a relative precision of $1.e-4$ for the residuum norm of the Newton iterations (there is also an absolute criterion for the residuum norm of $1.e-6$). For the staggered iterations, the criterion on the damage increment between two iterations is defined in the sense of the infinite norm and a precision of $1.e-4$ was used.

# 4. Numerical applications

## 4.1. Impact of a bar with a mechanical fuse

In this test, we consider the case of a 1D elastic bar, of unit length, with an homogeneous initial velocity, $v_0$, that hits a wall on one side. The elastic bar is made of 2 materials which only differ by their material toughness. The material in the middle of the bar is therefore more fragile and acts as a fuse. The geometry and the material properties are given on Fig. 3. The initial velocity is chosen such that the elastic strain energy in the first material is not large enough to fracture when the stress wave due to the impact will travel along the bar (initially the stress wave is a compressive wave and do not damage the bar but it transforms into tension after hitting the top of the bar). When the tension wave travels the fragile material, it breaks immediately and the initial bar is separated into two bars. The idea of this test is to illustrate some fundamental properties of the proposed space-time formulation: 1- the impact of the stabilization, with least square terms, on the dynamics is not important regarding the evolution of the damage (but it does for the velocity), 2- the formulation allows a correct evolution of the damage (i.e. irreversibility is always satisfied not only at interpolation/approximation points), 3- the formulation allows to capture the appropriate balance of energy with a minimal impact of the supplementary numerical terms (stabilization) on the numerical solution. Other authors have also used similar test with a bar in dynamics condition to investigate the numerical properties of phase-field models (see for instance [31,53]). The balance of energy is an important point mentioned by some authors to correctly capture the expected solution.

The exact solution of this test is easy to compute until the damage occurs. Knowing the speed of sound in the undamaged material, $c = \sqrt{E/\rho}$, we can determine that the crack will instantaneously break
<figure><img src="image_6.png" /><figcaption>Fig. 4. Value of d upon time at different location in the bar with and without stabilization, ST-IGA p2, mesh size: $\Delta x = 8.e - 3$, $\Delta t_0 = 2.e - 2$.</figcaption></figure>

the bar at $t = 2.5$ if the damage in material 1 stay null or very small. If the damage in the first material can not be neglected, the time for which the tension wave hits the material 2 can be modified. The bar is then supposed to separate immediately, and the remaining strain/kinetic energies in the two bars lead to travailing stress waves that are independent from each other. If the strain energy is sufficiently small (such as no supplementary damage occurs) the problem becomes conservative in each bar. However, in this example we have chosen to be not exactly in these conditions (i.e. null damage in material 1) to illustrate the fact that we can also have a damage wave prior and after the instantaneous crack in material 2. Therefore the analytical solution after the initiation of damage is not known exactly.

Fig. 4 shows the results obtained with quadratic (in time and space) ST-IGA elements (the mesh size for the space is $h = \Delta x = l/3.7$). As expected the damage initiates in the middle of the bar close to $t = 2.5$, obviously due to the time approximation the instantaneous solution (Heaviside function in time) can not be obtained exactly but it is only approached while refining the mesh. It can also be remarked that the stabilization terms do not seem to play a great role for the damage (some oscillations can be found if one looks closer to the solution without stabilization but their amplitudes stay small for this test). Furthermore, the tensile stress wave leads to a damage wave that propagates along the bar (even after that the material 2 has been fully broken). Obviously this damage wave modifies locally the Young modulus and consequently the speed of sound. The results in Fig. 4 have been computed using the time approximation functions and do not correspond only to the values at the degree of freedom (due to the non interpolating behavior of NURBS at the control points except on the boundaries of the patch). As previously mentioned it can be remarked that the condition $\dot{d} \ge 0$ is always satisfied. More generally, the space-time formulation makes it possible to verify that the irreversibility condition is satisfied everywhere in space and time in a “continuous” manner.

Fig. 5 illustrates that, as for the time, the space discretization with regular NURBS quadratic functions of maximum continuity leads to a regular solution which can only approach the expected solution (which should be closed to $\exp(-|x - 0.5|/l)$ in material 2, if material 1 was not damaged). Again, it can be concluded that the stabilization does not affect that much the solution for the damage (at least for the chosen value of the stabilization parameter).

Similar observations can be made with ST-FE (quadratic but $C_0$) regarding the stabilization parameter and the evolution of the damage. We have already compared in [47] ST-FE and ST-IGA for elastodynamics (for similar impact tests) therefore we focus here on damage. Fig. 6 shows that ST-FE and ST-IGA converge to nearly the same damage solution when refining the mesh. However, it can be remarked that the solutions obtained with ST-FE models predict a discontinuity in the

<figure><img src="image_7.png" /><figcaption>Fig. 5. Profile of the damage along the bar at the end of the test ($t=3.5$) with and without stabilization, ST-IGA p2, mesh size: $\Delta x = 8.e - 3$, $\Delta t_0 = 2.e - 2$.</figcaption></figure>

gradient of the damage (and therefore a discontinuity in the thermodynamics force $Y_{\nabla d}$) at the interface of the two materials. Obviously, ST-IGA solutions are more regular due to the $C_1$ continuity property of the damage between elements (while only $C_0$ continuity is imposed to ST-FE). For ST-IGA, we have also tested to use a lower continuity for the field $d$ (i.e. $C_0$ continuity between elements) could give similar results than ST-FE which is not the case. This is due to the fact that the kinematic field and the velocity field should also have a lower continuity, but in this case ST-IGA and ST-FE are perfectly equivalent (and so are the numerical results obtained).

The impact of the stabilization can be clearly seen from Fig. 7. The solution without stabilization (at $\tau = 0$) oscillates. The separation of the bar in two can be clearly seen after $t = 2.5$, the velocity in the two parts becomes of opposite sense (the bottom part get again compressed and the top part moves vertically with a constant velocity). One of the main advantages of the proposed space-time approach can be seen from Fig. 8 and relates to the energy conservation property. It can be seen from this figure that the total energy is nearly constant until damage starts to occur in material 1 (after $t = 2$) and propagates in material 2. After the full crack of the bar the total energy is again nearly constant. Obviously the stabilization with Galerkin least square terms induces some artificial dissipation but it remains small as long as the stabilization parameter, $\tau$, is well chosen (see [48] for a discussion on $\tau$).

## 4.2. Fragmentation of a bar

For this example we consider a similar test as proposed in [37], a bar is initially submitted to a velocity field (linear in $x$) such as the initial strain rate, $\dot{\epsilon}$, is constant in the bar. One part of the bar is here fixed and the opposite part is subjected to a constant velocity upon time (see Fig. 9). The expected solution of this test is a fragmentation in pieces of the bar as can be found in the literature. The number of pieces is expected to depend on the initial strain rate. As proposed by the authors in [37], we adopt a Weibull space<fnref n="5" /> distribution of the Young's modulus:

$$E(\alpha) = E_0 \sqrt{-\ln(\alpha)} + E_{min} \quad (37)$$

where $\alpha \in [0, 1]$ is a uniformly distributed random variable, $E_0$ and $E_{min}$ are defined from:

$$E_{min} = E(1 - 1.9130584e - 2), \quad E_0 = E(2.1586552e - 2) \quad (38)$$

<footnote n="5">⁵ The initial space distribution is constant regarding the time such as a material point in a bar do not see a change in the Young modulus at any time.</footnote>
<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Profile of the damage along the bar at the end of the test ($t=3.5$), comparison between quadratic ST-FE and ST-IGA for different mesh size for the space ($\tau = 0.01, \Delta t_0 = 2.e-2$).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 7.** Velocity history upon time at different location in the bar with and without stabilization, ST-IGA p2, mesh size: $\Delta x = 8.e-3, \Delta t_0 = 2.e-2$.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 8.** Energies integrated over the space domain upon the time with and without stabilization, ST-IGA p2, mesh size: $\Delta x = 8.e-3, \Delta t_0 = 2.e-2$.</figcaption></figure>

As in the previous example this test is unsized for space and time, the bar is of unit length, Young modulus and density are also equals to 1. For

the damage we use the following parameters: $G_c = 5.e-2, l = 2.e-3$. The results found in [37] tend to show that the damage evolves more or less homogeneously in the bar until a critical value is reached at which isolated cracks develop, leading to fragmentation of the bar. Fig. 10 illustrates some results obtained with quadratic ST-IGA and ST-FE, each corresponding to a different realization of the Young modulus. The other numerical parameters are perfectly identical between the simulations. It can be remarked that the results obtained with ST-IGA are more close to the one obtained in [37]. For ST-FE we have more difficulties in obtaining realizations for which we can achieve fragmentation and, in many cases these solutions lead to divergence as damage is close to 1. The ST-IGA model also permits to obtain expected results for very coarse meshes for the space (until $\Delta x = l/2$ as illustrate in Fig. 10(c), as a comparison in [37] very fine meshes are used).

From Fig. 10(d) it can be remarked that the energies obtained with ST-IGA for $\Delta x = l/2$ and $\Delta x = l/10$ are similar even if the strain energy is lowered and the crack energy<fnref n="6" /> is overestimated with the coarser mesh. The same figure also shows that for this initial strain rate, the appearance of cracks releases kinetic energy (locally) after fragmentation, which has an impact on the total kinetic energy of the bar.

Fig. 11, shows the results obtained at a higher strain rate for ST-IGA. As shown in the literature, the number of fragment increases with the

<footnote n="6">⁶ The term crack energy corresponds here to the crack surface density defined at eq. (2) multiplied by $G_c$.</footnote>
<figure><img src="image_11.png" /><figcaption>(a) Damage along the bar at different time, ST-IGA p2 with (b) Damage along the bar at different time, ST-FE p2 with
$\Delta x = l/10$</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(c) Damage along the bar at different time, ST-IGA p2 with
$\Delta x = l/2$</figcaption></figure>

<caption>(d) Energies upon time, ST-IGA p2</caption>

<figure><img src="image_13.png" /><figcaption>(a) Damage along the bar at different time, ST-IGA p2 with
$\Delta x = l/10$</figcaption></figure>

<caption>(b) Energies upon time, ST-IGA p2</caption>

<caption>Fig. 11. Fragmentation results of the bar for a initial strain rate $\dot{\epsilon} = 100$, with $\tau = 1.e-3$ and $\Delta t_0 = 1.e-3$, ST-IGA p2 with $\Delta x = l/10$.</caption>
<figure><img src="image_14.png" /><figcaption>**Fig. 12.** Geometry and boundary condition of the L-shape example, plane strain case.</figcaption></figure>

strain rate. It is also interesting to remark that the balance of energy is different. The kinetic energy is no more impacted by the occurrence of cracks. The crack energy also increases until total fragmentation even if the strain energy as started to decreases. It can therefore be supposed that further increases the strain rate will not increase the number of fragments indefinitely. It will only modify the total kinetic energy and the time it takes for cracks to appear (this effect is also shown in [37]). From this simple test it can be concluded that the higher degree of continuity of ST-IGA may have an interest for fragmentation problems, in particular due to the fact that relatively coarse meshes can be used.

### 4.3. L-shape impact

In the following example, we propose a mixed mode dynamic fracture of a L-specimen. This test is a simplified version of the one presented in [43]. The geometry and the boundary condition are given in Fig. 12 and we assume a plane strain state.

The material parameters are taken such as the Young modulus is $E = 1 \text{ N/m}^2$, the mass density is $\rho = 1 \text{ kg/m}^3$, the Poisson ratio is $\nu = 0.3$ and the material toughness is $g_c = 3 \times 10^{-3} \text{ J/m}^2$. The length parameter is $l = 5 \times 10^{-3}$ m. The space-time domain is discretized such that the element size in space $h$ respects the condition $h < l/2$. The stabilization parameter is chosen such that $\tau = 1 \times 10^{-3}$ and the time slab depth is chosen such as $\Delta t_{slab} = 10^{-2}$ s. The aim of this test is to highlight the dynamic influence on crack propagation. To accomplish this, the specimen is exposed to varying loading rates such as $v(t) = v_0$ and $v(t) = 8v_0$, with $v_0 = 0.1 \text{ m/s}$. For this example we consider only ST-FE (which are again quadratic in space and time).

Fig. 13 illustrates the influence of dynamics on the crack pattern. As expected it can be seen that if the loading rate increases, the angle of the crack path relative to the horizontal axis also increases.

In Fig. 14, we show the evolution of the elastic and dissipate (crack) energy during the crack evolution for different loading rates. By examining these curves, we can observe that they are consistent with expectations. In the case where $v = 8v_0$, we observe that the dissipated energy increases rapidly, in contrast to the case where $v = v_0$, where a completely elastic phase (no dissipation) occurs before it starts to evolve. The evolution of the dissipated energy indirectly reflects the evolution of damage, which is accompanied by a decrease in elastic energy. For the kinetic energy, see Fig. 15, the case $v = v_0$ is very close to the quasi-static case. It is interesting to remark that when the geometry is fully broken for the case $v = v_0$, the strain energy and the kinetic energy stay nearly constant, illustrating the fact that the model is broken in two part with one part that is free to move without generating any stress. As in

the previous examples, we see that this space-time approach is interesting to capture dynamic solutions without adding significant undesired numerical dissipation. These observations are a consistent demonstration of the impact of dynamics on crack propagation.

## 5. Conclusion

This paper proposes an original space-time approach for dynamics brittle fracture. There are two main contributions, the first one is theoretical and comes from the derivation of a three fields (displacement, velocity, damage) space-time formulation from a space-time potential that incorporates stabilization terms (with Galerkin least squares) to limit the spurious oscillations that come from the discretization of the continuous problem. The second one is the numerical form with the help of a time discontinuous Galerkin method. We have derived consistent continuity terms that allow to use a staggered resolution of the problem (displacement/velocity and damage are resolved separately) without adding supplementary numerical parameters. The use of this discontinuous Galerkin method allowed us to use the history function such as to guaranty the irreversibility of the damage. Furthermore for the discretization we can use either ST-IGA or ST-FE. The numerical examples illustrated some of the possibilities offered by such an approach. In particular we have shown that we can approach with a good qualitative agreement the expected solutions while getting rid of oscillations often observed in dynamics and having a correct balance of energy.

As a perspective we would like to emphasize the possibilities offers by this approach. From earlier studies (without damage) we know that ST-IGA and ST-FE can reach optimum rate of convergence in space and time and parallelization strategies can be used to solve a space-time slab. These aspects pave the way for the development of more efficient solvers for cracks propagation problems which are computationally intensive. Moreover, we obtain much more data from a ST approach as we have an approximate solution for all times and this could have an interest to feed machine learning algorithm for instance. The space-time formulation is also a powerful tool to account of different physics in a single numerical scheme. We can therefore extend the proposed approach to the case of dynamics cracks with a thermomechanical coupling (or other couplings, electrical, magnetical, etc) or to ductile fracture models. We can also account of more tight couplings between damage and other physics (for instance if the material toughness is dependent on a state variable). This study is only a step towards more complex space-time models.

### CRediT authorship contribution statement

F.K. Feutang: Writing – review & editing, Software, Methodology, Investigation. S. Lejeunes: Conceptualization, Methodology, Software, Supervision, Writing – original draft. D. Eyheramendy: Conceptualization, Supervision, Writing – review & editing.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Appendix A. Influence of the $\alpha$ coefficient for the imposition of initial conditions

As mentioned previously, the coefficients $\alpha_u$, $\alpha_v$ and $\alpha_d$ should not be chosen arbitrarily. We report here in Table A.1 the impact of these coefficients for the first application of the paper: the bar fragmentation test. It can be seen that due to the weak imposition of the continuity terms, the imposed initial conditions are only satisfied up to a certain accuracy
<figure><img src="image_15.png" /><figcaption>(a) Damage evolution for $v = v_0$</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>(b) Damage evolution for $v = 8v_0$</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>(a) Strain energy upon space domain</figcaption></figure>

<caption>(b) Dissipate energy upon space domain</caption>

<caption>Fig. 14. Evolution of strain energy and dissipate energy for different loading rate.</caption>

(the value for the velocity has not been printed with a sufficient precision and is also not perfectly accurate with respect to the imposed one). This behavior is expected, as is the fact that the initial field will not be perfectly homogeneous in space even if specified from the initial condi-

tion. As a consequence, the discrepancy from the exact values can play a role of a perturbation of the initial condition that could be problematic for the convergence if not sufficiently small (which was not observed for the test we have done). It can also be seen that the coefficients have
<figure><img src="image_18.png" /><figcaption>Fig. 15. Evolution of Kinetic energy for different loading rate.</figcaption></figure>

<table><caption>**Table A.1**

Influence of the coefficients for the initial condition: test case fragmentation of a bar with $\tau = 0.01$, $\Delta x = 8.e-3$, $\Delta t_0 = 2.e-2$.</caption><thead><tr><th></th><th>u</th><th>v</th><th>d</th></tr></thead><tbody><tr><td>imposed initial conditions (\&x)</td><td>0.</td><td>-0.5</td><td>0.</td></tr><tr><td>α<sub>u</sub>, α<sub>d</sub>, α<sub>d</sub> = 1 (x = 0.3)</td><td>-8.47e-19</td><td>-0.5</td><td>2.78e-24</td></tr><tr><td>α<sub>u</sub>, α<sub>d</sub>, α<sub>d</sub> = 1.e-10 (x = 0.3)</td><td>3.74e-19</td><td>-0.5</td><td>5.90e-24</td></tr><tr><td>α<sub>u</sub>, α<sub>d</sub>, α<sub>d</sub> = 1.e-4 (x = 0.3)</td><td>4.18e-19</td><td>-0.5</td><td>1.04e-24</td></tr><tr><td>α<sub>u</sub>, α<sub>d</sub>, α<sub>d</sub> = 1.e+4 (x = 0.3)</td><td>1.03e-31</td><td>-0.5</td><td>5.62e-28</td></tr><tr><td>α<sub>u</sub>, α<sub>d</sub>, α<sub>d</sub> = 1.e+10 (x = 0.3)</td><td>1.00e-37</td><td>-0.5</td><td>5.24e-34</td></tr></tbody></table>

an impact but it is rather limited at least for the values we have tested and for the test considered.

## Data availability

The data that has been used is confidential.

## References

[1] Ambati M, Gerasimov T, De Lorenzis L. Phase-field modeling of ductile fracture. Comput Mech 2015;55:1017–40.

[2] Ambati M, Gerasimov T, De Lorenzis L. A review on phase-field models of brittle fracture and a new fast hybrid formulation. Comput Mech 2015;55:383–405.

[3] Ambrosio L, Tortorelli V. Approximation of functional depending on jumps by elliptic functional via gamma-convergence. Commun Pure Appl Math 1990;43:999–1036.

[4] Ambrosio L, Tortorelli VM. Approximation of functional depending on jumps by elliptic functional via t-convergence. Commun Pure Appl Math 1990;43:999–1036.

[5] Amor H, Marigo JJ, Maurini C. Regularized formulation of the variational brittle fracture with unilateral contact: numerical experiments. J Mech Phys Solids 2009;57:1209–29.

[6] Argyris J, Scharpf D. Finite elements in time and space. Aeronaut J 1969;73:1041–4.

[7] Bleyer J, Roux-Langlois C, Molinari JF. Dynamic crack propagation with a variational phase-field model: limiting speed, crack branching and velocity-toughening mechanisms. Int J Fract 2017;204:79–100.

[8] Bonilla J, Badia S. Maximum-principle preserving space-time isogeometric analysis. Comput Methods Appl Mech Eng 2019;354:422–40.

[9] Borden MJ, Hughes TJ, Landis CM, Verhoosel CV. A higher-order phase-field model for brittle fracture: formulation and analysis within the isogeometric analysis framework. Comput Methods Appl Mech Eng 2014;273:100–18.

[10] Borden MJ, Verhoosel CV, Scott MA, Hughes TJ, Landis CM. A phase-field description of dynamic brittle fracture. Comput Methods Appl Mech Eng 2012;217:77–95.

[11] Bourdin B, Francfort GA, Marigo JJ. Numerical experiments in revisited brittle fracture. J Mech Phys Solids 2000;48:797–826.

[12] Bourdin B, Francfort GA, Marigo JJ. The variational approach to fracture. J Elast 2008;91:5–148.

[13] Bourdin B, Larsen CJ, Richardson CL. A time-discrete model for dynamic fracture based on crack regularization. Int J Fract 2011;168:133–43.

[14] Bruch Jr JC, Zyvoloski G. Transient two-dimensional heat conduction problems solved by the finite element method. Int J Numer Methods Eng 1974;8:481–94.

[15] De Lorenzis L, Maurini C. Nucleation under multi-axial loading in variational phase-field models of brittle fracture. Int J Fract 2022;237:61–81.

[16] Denis K, Christian SM, Thomas W. Space-time mixed system formulation of phase-field fracture optimal control problems. J Optim Theory Appl 2023;199:1222–48. https://doi.org/10.1007/s10957-023-02272-7.

[17] Farrell P, Maurini C. Linear and nonlinear solvers for variational phase-field models of brittle fracture. Int J Numer Methods Eng 2017;109:648–67.

[18] Francfort GA, Marigo JJ. Revisiting brittle fracture as an energy minimization problem. J Mech Phys Solids 1998;46:1319–42.

[19] Frémond M, Nedjar B. Damage, gradient of damage and principle of virtual power. Int J Solids Struct 1996;33:1083–103.

[20] Gerasimov T, De Lorenzis L. A line search assisted monolithic approach for phase-field computing of brittle fracture. Comput Methods Appl Mech Eng 2016;312:276–303.

[21] Gerasimov T, De Lorenzis L. On penalization in variational phase-field models of brittle fracture. Comput Methods Appl Mech Eng 2019;354:990–1026.

[22] Hesch C, Schuß S, Dittmann M, Eugster S, Favino M, Krause R. Variational space-time elements for large-scale systems. Comput Methods Appl Mech Eng 2017;326:541–72.

[23] Hofacker M, Miehe C. A phase field model of dynamic fracture: robust field updates for the analysis of complex crack patterns. Int J Numer Methods Eng 2013;93:276–301.

[24] Hughes T. The finite element method. Prentice-Hall INC.; 1987.

[25] Hughes TJ, Cottrell JA, Bazilevys Y. Isogeometric analysis: cad, finite elements, nurbs, exact geometry and mesh refinement. Comput Methods Appl Mech Eng 2005;194:4135–95.

[26] Hughes TJ, Franca LP, Hulbert GM. A new finite element formulation for computational fluid dynamics: VIII. The Galerkin/least-squares method for advective-diffusive equations. Comput Methods Appl Mech Eng 1989;73:173–89.

[27] Hulbert GM, Hughes TJ. Space-time finite element methods for second-order hyperbolic equations. Comput Methods Appl Mech Eng 1990;84:327–48.

[28] Jamet P, Bonnerot R. Numerical solution of the Eulerian equations of compressible flow by a finite element method which follows the free boundary and the interfaces. J Comput Phys 1975;18:21–45.

[29] Johnson C, Pitkäranta J. An analysis of the discontinuous Galerkin method for a scalar hyperbolic equation. Math Comput 1986;46:1–26.

[30] Langer U, Moore SE, Neumüller M. Space-time isogeometric analysis of parabolic evolution problems. Comput Methods Appl Mech Eng 2016;306:342–63.

[31] Li T. Gradient-damage modeling of dynamic brittle fracture: variational principles and numerical simulations. Ph.D. thesis. Université Paris-Saclay; 2016.

[32] Linse T, Hennig P, Kästner M, de Borst R. A convergence study of phase-field models for brittle fracture. Eng Fract Mech 2017;184:307–18. https://doi.org/10.1016/j.engfracmech.2017.09.013. https://www.sciencedirect.com/science/article/pii/S0013794417307488.

[33] Liu G, Li Q, Msekh MA, Zuo Z. Abaqus implementation of monolithic and staggered schemes for quasi-static and dynamic fracture phase-field model. Comput Mater Sci 2016;121:35–47.

[34] Miehe C, Hofacker M, Welschinger F. A phase field model for rate-independent crack propagation: robust algorithmic implementation based on operator splits. Comput Methods Appl Mech Eng 2010;199:2765–78.

[35] Miehe C, Welschinger F, Hofacker M. Thermodynamically consistent phase-field models of fracture: variational principles and multi-field fe implementations. Int J Numer Methods Eng 2010;83:1273–311.

[36] Molnár G, Gravouil A. 2d and 3d abaqus implementation of a robust staggered phase-field solution for modeling brittle fracture. Finite Elem Anal Des 2017;130:27–38.

[37] Moës N, Lé B, Sterhicic A. Fragmentation analysis of a bar with the lip-field approach. Mech Mater 2022;172:104365. https://doi.org/10.1016/j.mechmat.2022.104365. https://www.sciencedirect.com/science/article/pii/S0167663622001399.

[38] Msekh MA, Sargadoh JM, Jamshidian M, Areias PM, Rabczuk T. Abaqus implementation of phase-field model for brittle fracture. Comput Mater Sci 2015;96:472–84.

[39] Mumford DB, Shah J. Optimal approximations by piecewise smooth functions and associated variational problems. Commun Pure Appl Math 1989.

[40] Nguyen H, Reynen J. A space-time least-square finite element scheme for advection-diffusion equations. Comput Methods Appl Mech Eng 1984;42:331–42.

[41] Nguyen QS, Andreux S. The non-local generalized standard approach: a consistent gradient theory. C R, Méc 2005;333:139–45. https://doi.org/10.1016/j.crme.2004.09.010.

[42] Oden JT. A general theory of finite elements. II. Applications. Int J Numer Methods Eng 1969;1:247–59.

[43] Özbolt J, Bede N, Sharma A, Mayer U. Dynamic fracture of concrete l-specimen: experimental and numerical study. Eng Fract Mech 2015;148:27–41.

[44] Pham K, Amor H, Marigo JJ, Maurini C. Gradient damage models and their use to approximate brittle fracture. Int J Damage Mech 2011;20:618–52.

[45] Philipp J, Thomas W. Space-time variational material modeling: a new paradigm demonstrated for thermo-mechanically coupled wave propagation, visco-elasticity, elasto-plasticity with hardening, and gradient-enhanced damage. Comput Mech 2023;73:365–402. https://doi.org/10.1007/s00466-023-02371-2.

[46] Podhorecki A. The viscoelastic space-time element. Comput Struct 1986;23:535–44.

[47] Saadé C, Lejeunes S, Eyheramendy D, Saad R. Space-time isogeometric analysis for linear and non-linear elastodynamics. Comput Struct 2021;254:106594.

[48] Saadé C. Space-time iso-geometric methods for multi-field equations in mechanics. Ph.D. thesis. Ecole Centrale de Marseille; 2020. https://www.theses.fr/2020ECDM0011.
[49] Storm J, Kaliske M. The origin of the energy split in phase-field fracture and eigen-fracture. PAMM 2023;23. https://api.semanticscholar.org/CorpusID:262210038.

[50] Tezduyar TE, Takizawa K. Space-time computational flow analysis: unconditional methods and first-ever solutions. Comput Methods Appl Mech Eng 2023;116137. https://doi.org/10.1016/j.cma.2023.116137. https://www.sciencedirect.com/science/article/pii/S004578252300261X.

[51] Vajari SA, Neuner M, Arunachala PK, Ziccarelli A, Deierlein G, Linder C. A thermodynamically consistent finite strain phase field approach to ductile fracture considering multi-axial stress states. Comput Methods Appl Mech Eng 2022;400:115467.

[52] Wick T. Modified Newton methods for solving fully monolithic phase-field quasi-static brittle fracture propagation. Comput Methods Appl Mech Eng 2017;325:577-611.

[53] Yin B, Storm J, Kaliske M. Viscoelastic phase-field fracture using the framework of representative crack elements. Int J Fract 2022;237:139-63. https://doi.org/10.1007/s10704-021-00522-1.

[54] Zhou S, Zhuang X, Rabczuk T. Phase-field modeling of fluid-driven dynamic cracking in porous media. Comput Methods Appl Mech Eng 2019;350:169-98.