C Saadé

S Lejeunes
email: lejeunes@lma.cnrs-mrs.fr

D Eyheramendy

R Saad

Space-Time Isogeometric Analysis for linear and non-linear elastodynamics

Keywords: Space-Time, Isogeometric Analysis, Stabilization techniques, Impact

In this study we investigate a stabilized space-time formulation for linear and non-linear elastodynamics. We use Isogeometric Analysis (IGA) in order to benefit from its numerical qualities. We focus on two points: the formulation of stabilized weak-forms in a linear and non-linear context and the interest of using continuous Galerkin schemes in space and time with higher order and higher continuity basis functions. We illustrate the numerical performance of these methods through typical impact or vibration problems commonly encountered in the field of the elastodynamics of solids.

Introduction

The earliest ideas on space-time formulations for solving time-dependent mechanical problems came from the pioneering work of Oden [START_REF] Oden | A general theory of finite elements. II. Applications[END_REF] and Argyris and Scharpf [START_REF] Argyris | Finite elements in time and space[END_REF]. In Argyris and Scharpf [START_REF] Argyris | Finite elements in time and space[END_REF], the authors proposed a standard finite element formulation for space. The corresponding discrete time-dependent problem was handled using a second finite element interpolation over a time interval. The authors used a Hamiltonian formulation considering both velocity and displacement and showed that the elastodynamics problem leads to the resolution of sequential linear problems. This approach implies a separation of time and space and therefore different polynomial approximations can be used for space and time. In one section of his very general paper on finite elements, Oden [START_REF] Oden | A general theory of finite elements. II. Applications[END_REF] proposed another approach for which space and time are considered as an augmented dimensional finite element problem. Oden [START_REF] Oden | A general theory of finite elements. II. Applications[END_REF], proposed the concept of meshing the space-time domain with finite elements in a more general manner. In this approach, polynomial approximation should be the same for space and time, as the finite element can be arbitrarily positioned in the space-time mesh. These ideas of space-time methods did not gain much popularity during the next 20 years due to the fact that during the 70s more robust and efficient implicit algorithms, e.g. [START_REF] Hilber | Improved numerical dissipation for time integration algorithms in structural dynamics[END_REF], were developed and applied to engineering problems with an acceptable numerical cost. Some authors proposed some applications of space-time finite elements for elastodynamics such as linear viscoelasticity in the case of a 1D problem [START_REF] Podhorecki | The viscoelastic space-time element[END_REF], a nonlinear 1D problem [START_REF] Cella | Space-time elements for the shock wave propagation problem[END_REF] or in the case of a linear 1D beam and truss or 2D plane strain problem [START_REF] Bajer | Triangular and tetrahedral space-time finite elements in vibration analysis[END_REF][START_REF] Bajer | Notes on the stability of non-rectangular space-time finite elements[END_REF]. Other applications than elastodynamics were also proposed, such as free surface problems (see [START_REF] Bonnerot | Numerical computation of the free boundary for the two-dimensional stefan problem by space-time finite elements[END_REF]), compressible fluid flows [START_REF] Jamet | Numerical solution of the eulerian equations of compressible flow by a finite element method which follows the free boundary and the interfaces[END_REF], heat transfer problems [START_REF] Bruch | Transient two-dimensional heat conduction problems solved by the finite element method[END_REF][START_REF] Chung | The fourth-dimention concept in the finite element analysis of transient heat transfer problems[END_REF], advection-diffusion equations [START_REF] Nguyen | A space-time least-square finite element scheme for advection-diffusion equations[END_REF]. However, most of the previously mentioned studies assumed either a separation of space and time in the way the approximation and the mesh are built (always structured mesh), or a time integration done with simplified quadrature rules reducing in some cases to a Crank-Nicolson scheme for the time discretization.

New ideas started to emerge in the early 90s: a Time Discontinuous Galerkin (TDG) method was proposed for elastodynamics [START_REF] Hughes | Space-time finite element methods for elastodynamics: Formulations and error estimates[END_REF][START_REF] Hulbert | Space-time finite element methods for second-order hyperbolic equations[END_REF]. It combined a discontinuous Galerkin formulation, a stabilized weak form using least-square terms and discontinuity capturing operators, see also [START_REF] Hulbert | Discontinuity-capturing operators for elastodynamics[END_REF]. Least-square terms had previously been developed in another context by one of the authors, see [START_REF] Franca | The galerkin gradient least-squares method[END_REF][START_REF] Loula | Mixed Petrov-Galerkin methods for the Timoshenko beam problem[END_REF][START_REF] Loula | Stability, convergence and accuracy of a new finite element method for the circular arch problem[END_REF]. The Time Discontinuous Galerkin (TDG) method and the Discontinuous Galerkin method (DGM) encountered more success. They have elicited numerous papers presenting theoretical and numerical results, see for instance [START_REF] French | A space-time finite element method for the wave equation[END_REF][START_REF] Li | Implementation and adaptivity of a space-time finite element method for structural dynamics[END_REF][START_REF] Ekevid | Adaptive fea of wave propagation induced by high-speed trains[END_REF][START_REF] Huang | On the use of space-time finite elements in the solution of elasto-dynamic problems with strain discontinuities[END_REF][START_REF] Chien | Three-dimensional transient elastodynamic analysis by a space and time-discontinuous galerkin finite element method[END_REF][START_REF] Costanzo | Proof of unconditional stability for a single-field discontinuous galerkin finite element formulation for linear elasto-dynamics[END_REF][START_REF] Thompson | Adaptive space-time finite element methods for the wave equation on unbounded domains[END_REF][START_REF] Khalmanova | A space-time discontinuous galerkin finite element method for fully coupled linear thermo-elasto-dynamic problems with strain and heat flux discontinuities[END_REF][START_REF] Petersen | A space-time discontinuous galerkin method for the solution of the wave equation in the time domain[END_REF]. TDG methods have also given rise to other original methods such as the TX-FEM (Time Extended Finite Element Method) [START_REF] Réthoré | A combined space-time extended finite element method[END_REF] also known as the enriched space-time FEM (eg. [START_REF] Zilian | The enriched space-time finite element method (est) for simultaneous solution of fluid-structure interaction[END_REF][START_REF] Zilian | Hybridized enriched space-time finite element method for analysis of thin-walled structures immersed in generalized newtonian fluids[END_REF]), the Nitsche XFEM-DG method [START_REF] Lehrenfeld | Analysis of a Nitsche XFEM-DG discretization for a class of two-phase mass transport problems[END_REF], or the combination of the TDG and the PGD (Proper Generalized Decomposition) [START_REF] Boucinha | Space-time proper generalized decompositions for the resolution of transient elastodynamic models[END_REF].

On the other hand, some authors have worked on continuous formulations, such as the Time Continuous Galerkin (TCG) method. These formulations are similar to the original idea of Oden [START_REF] Oden | A general theory of finite elements. II. Applications[END_REF], i.e., no separation of space and time is made and the space-time domain is meshed either with a structured mesh or with an unstructured mesh. In particular, some authors have shown the interest of using the p-version of the finite element method in space and time, [START_REF] Hou | Application of triangular space-time finite elements to problems of wave propagation[END_REF]. Other have provided comparisons between TCG and TDG for elastodynamics problems with and without discontinuity. They have shown that the convergence of TCG is faster than that of TDG for elementary 1D problems [START_REF] Idesman | Solution of linear elastodynamics problems with space-time finite elements on structured and unstructured meshes[END_REF]. An important aspect in the improvement of the numerical efficiency is the use of remeshing strategies, which constitute an open field of research. Applications to 4D problems have also been proposed, e.g. [START_REF] Dumont | 4d Remeshing Using a Space-Time Finite Element Method for Elastodynamics Problems[END_REF][START_REF] Karyofylli | Simplex space-time meshes in two-phase flow simulations[END_REF].

Regardless of the approach used, an essential aspect of the formulation of space-time methods is the choice of the weak form. For elastodynamics problems, a common approach is to use the Hamiltonian principle. It has been used by many authors both for displacement and displacement-velocity formulations. It is also a natural way to develop semi-discrete time integrators [START_REF] Gravouil | Heterogeneous asynchronous time integrators for computational structural dynamics[END_REF]. Clearly, the disadvantage of using a two-field formulation (displacement-velocity) is the increase in the number of unknowns in the resulting system, but such formulations make it possible to impose straightforwardly velocity boundary and initial conditions. Besides the Hamiltonian principle, weak forms can be formulated using arbitrary test functions and/or their derivatives (also by using single or two-field forms). The stability of the different formulations is a current issue, which is discussed in the literature for different types of equations. Nevertheless, as can be seen from previous work, e.g. [START_REF] Hughes | Space-time finite element methods for elastodynamics: Formulations and error estimates[END_REF], space-time formulations can exhibit undesired oscillations and the introduction of dissipative terms using least square terms helps in reducing them.

The Finite Elements Method (FEM) is a highly versatile method but has a major disadvantage in the description of the geometry. It uses simple geometric elements and can therefore induce approximation errors and difficulties in representing complex geometries. Unlike the finite elements method, the Isogeometric1 Analysis (IGA) [START_REF] Hughes | Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement[END_REF][START_REF] Cottrell | Isogeometric analysis: toward integration of CAD and FEA[END_REF], makes it possible to exactly describe complex geometries using NURBS or B-Spline functions in a similar way to what is done by CAD software. Interested readers can find more informations on NURBS and B-Spline representation of geometry in [START_REF] Piegl | The NURBS Book. Monographs in Visual Communication[END_REF] and reference therein. Obviously, alternative methods to IGA that share the idea of an exact description of the geometry can be found in the literature, such as the NURBS-enhanced FEM (e.g. [START_REF] Sevilla | 3d nurbs-enhanced finite element method (nefem)[END_REF]). The IGA also offers great flexibility in the definition of approximation spaces, by allowing the control of the order of continuity of the approximation functions at certain discretization points. Space-Time IGA methods (either continuous or discontinuous) can therefore be an attractive alternative to Space-Time FEM methods. Some authors have reported interesting results obtained with such methods: in [START_REF] Langer | Space-Time Isogeometric Analysis of Parabolic Evolution Equations[END_REF] for a linear parabolic problem using a continuous formulation, in [START_REF] Loli | An efficient solver for space-time isogeometric galerkin methods for parabolic problems[END_REF][START_REF] Montardini | Space-time least-squares isogeometric method and efficient solver for parabolic problems[END_REF] a specific preconditioner is developed to enhance performance of continuous ST-IGA formulation also for linear parabolic problem, in [START_REF] Kuraishi | Space-time computational analysis of tire aerodynamics with actual geometry, road contact, tire deformation, road roughness and fluid film[END_REF] and [START_REF] Kuraishi | Space-time isogeometric flow analysis with built-in reynolds-equation limit[END_REF] and references therein ST-IGA is used for fluid flow or fluid-solid interaction problems or in [START_REF] Hesch | Variational space-time elements for largescale systems[END_REF] for the Kuramoto-Sivashinsky equation and in [START_REF] Bonilla | Maximum-principle preserving space-time isogeometric analysis[END_REF] for a scalar advection-diffusion problem.

In this work, our main objective is to propose a continuous ST-IGA formulation that is sufficiently general and versatile to address a variety of linear and non-linear elastodynamics problems. In the case of compressible elasticity, we adopt a two-field (displacement/velocity) formulation with additional least square terms. These supplementary terms are simpler than the ones found in the literature and allow a straightforward extension to non-linear cases. This form can be interpreted as a penalized formulation with a consistency condition on the acceleration field. It makes sense with the higher continuity order property (compared with piecewise Lagrange polynomials of the FEM) of the NURBS-based approximation. We also propose a comparison of the numerical results obtained using our ST-IGA formulations with the results obtained using continuous ST-FEM and semi-discrete FEM simulations with the HHT scheme for transient dynamic problems (with and without discontinuities). We also show that the continuous form is easily extendable to nearly-incompressible hyperelasticity with a displacement/velocity/pressure formulation.

The paper is organized as follows. In a first section we present the stabilized weak form for small strain elasticity (compressible case). In the second section, the formulation is extended to the case of compressible and quasi-incompressible hyperelasticity. The space-time discretization is detailed in the third section. The fourth and fifth sections are devoted to the presentation of numerical examples in linear and non-linear cases.

Small-strain elastodynamics

Problem statement

We consider the space-time cylinder Q = Ω × [0, T ], where Ω ⊂ R d is a closed domain such that d is the number of space dimensions and T is the final time. The boundary of Ω is denoted by Γ. Let Γ u and Γ σ denote the non-overlapping subregions of Γ such that : Γ = Γ u ∪ Γ σ , Γ u ∩ Γ σ = φ. The conservation of the linear momentum on Q takes the following form in the small strain case:

ρü -div x σ = f ∀(x, t) ∈ Q (1)
where ρ(x) is the material density, u(x, t) is the displacement field, σ(x, t) is the stress tensor and f (x, t) is the applied body load. Assuming linear isotropic elasticity, the stress tensor can be related to the strain field, ε(x, t) = 1 2 (∇ x u + ∇ T x u), by the following constitutive equation

σ = 2µε + λtr(ε)1 (2)
where λ(x, t) and µ(x, t) are Lamé elasticity parameters and 1 is the second order identity tensor. The problem is fully defined by taking into account the boundary and initial conditions:

u = g(x, t) for x ∈ Γ u , t ∈ [0, T ] σ • n = T(x, t) for x ∈ Γ σ , t ∈ [0, T ] u(x, t = 0) = u 0 (x) for x ∈ Ω u(x, t = 0) = v 0 (x) for x ∈ Ω (3)
where n(x) is the normal to the boundary Γ σ .

A two field weak form

We define the Sobolev spaces

H l,k (Q) = {u ∈ L 2 (Q) : ∂ α x u ∈ L 2 (Q) ∀ α with 0 ≤ |α| ≤ l, ∂ i t ∈ L 2 (Q), i = 0, . . . , k} of functions defined in the space-time cylinder Q, where L 2 (Q) denotes the space of square-integrable functions, α = (α 1 , . . . , α d ) is a multi-index with non-negative integers, |α| = α 1 + . . . + α d , ∂ α x u := ∂ α1 x1 ∂ α2 x2 . . . ∂ α d x d u with ∂ αi xi • = ∂ αi • /∂x i αi and ∂ i t u := ∂ i u/∂t i .
In order to obtain a two field weak form, we first have to introduce the velocity field, v(x, t) = u and then transform the second order problem of eq. ( 1) into a first order system, such that:

ρ v -div x σ -f = 0 ρ(v -u) = 0 ∀(x, t) ∈ Q u = g(x, t) for x ∈ Γ u , t ∈ [0, T ] σ • n = T(x, t) for x ∈ Γ σ , t ∈ [0, T ] u(x, t = 0) = u 0 (x) for x ∈ Ω v(x, t = 0) = v 0 (x) for x ∈ Ω (4)
Weighting the first two equations of the system (4) by the time derivatives of test functions defined in appropriate spaces and integrating the equation over the space-time cylinder Q, we obtain the following weak variational formulation:

Find (u(x, t), v(x, t)) ∈ H u × H v , such that ∀ (δu(x, t), δv(x, t)) ∈ H u 0 × H v 0 Q (ρ v -div x σ -f ) δ udQ + Q ρ (v -u) δ vdQ = 0 (5)
where

H u = {u ∈ H 1,1 (Q), u = g on Γ u , u(x, t = 0) = u 0 (x)}, H v = {v ∈ H 0,1 (Q) , v = ġ on Γ u , v(x, t = 0) = v 0 (x)}, H u 0 = {u ∈ H 1,1 (Q), u = 0 on Γ u , u(x, t = 0) = 0}, H v 0 = {v ∈ H 0,1 (Q), v = 0 on Γ u , v(x, t = 0) = 0}.
After integrating by part the divergence operator and using traction boundary conditions, the weak variational formulation of eq [START_REF] Bajer | Notes on the stability of non-rectangular space-time finite elements[END_REF] becomes:

Find (u(x, t), v(x, t)) ∈ H u × H v , such that ∀ (δu(x, t), δv(x, t)) ∈ H u 0 × H v 0 Q ρ v δ udQ + Q σ : ε( δ u)dQ - P T δ udP - Q f δ udQ + Q ρ (v -u) δ vdQ = 0 (6)
where P = Γ σ × [0, T ] is the space time boundary where external forces apply. Remarks 1. The previous system is similar to the formulation adopted in [START_REF] Hughes | Space-time finite element methods for elastodynamics: Formulations and error estimates[END_REF] but a continuous Galerkin formulation is used here and no least-square terms are added (we consider additional terms in the following section). 2. The choice of using the time derivatives of test functions is made here as it has been done previously by other authors, e.g. [START_REF] Hulbert | Space-time finite element methods for second-order hyperbolic equations[END_REF]. The first part of eq. ( 6) can also be obtained by applying the principle of virtual power if we consider δ u to be a virtual velocity field. Other formulations can be found in the literature where the time derivatives of the test functions were not used. By proceeding in this way one obtains a form similar to the principle of virtual work. We can also consider the stationarity of a Hamiltonian form. Consider, for instance, the following Lagrangian L:

L(u, u, v) = K HR ( u, v) -P(u) (7)
where K HR and P are the kinetic energy (modified according to the Hellinger-Reisner mixed form) and the potential energy defined such that:

K HR ( u, v) = Ω ρv u - 1 2 ρvvdΩ P(u) = Ω 1 2 σ : εdΩ - Ω fudΩ - Γσ TudS (8)
Using the stationarity of the Hamiltonian, H(u, u, v) = T 0 L(u, u, v)dt, leads to:

- Q ρ(v -u)δvdQ + Q ρvδ udQ - Q σ : ε(δu)dQ + P TδudP + Q f δudQ = 0 (9)
By integrating by part the second integral and choosing δu such that δu(t = 0) = 0 and δu(t = T ) = 0, we obtain:

- Q ρ(v -u)δvdQ - Q ρ vδudQ - Q σ : ε(δu)dQ + P TδudP + Q f δudQ = 0 ( 10
)
which is the same as eq. ( 6) after replacing δu by δ u and δv by δ v. Note that the variational form of eq. ( 6) leads to a non-symmetric bilinear form while that of eq. ( 10) is symmetric, as expected. For some specific choices for the approximation space of u and v and for structured meshes, classical time integration schemes such as the Crank-Nicholson scheme can be recovered in the case of piecewise linear approximation functions. However, the use of continuous B-Spline approximations for test functions and possibly the use of unstructured meshes allow the construction of richer time integration schemes. 3. One advantage of this two field formulation is that it makes it possible to impose the velocity conditions (initial or boundaries) as Dirichlet conditions.

Stabilization of the discrete form in the presence of discontinuities: GAC and GLS

As shown in previous references (e.g. [START_REF] Franca | The galerkin gradient least-squares method[END_REF]), the discrete forms that can be obtained from the weak form of eq. ( 6) do not make it possible to get rid of spurious numerical oscillations that occur in some typical elastodynamics problems for which time discontinuities appear, such as impact problems (as will be illustrated in section 5). The addition of supplementary terms can be an appropriate and efficient strategy for controlling these oscillations. In this section, continuous Galerkin forms having least square terms are used. Without any restrictions on the choice of the space-time discretization (FE or IGA), u h and v h are taken to be the continuous approximation fields of the displacement and velocity, respectively. The space time cylinder is discretized. We denote by Q h its discrete form, such that Q h = ∪ e Q e where Q e are elemental entities. We denote by F p h a functional space that is build from polynomials of degree p and we introduce the subspaces:

H u h = {u h ∈ (F p h ∩ H 1,1 (Q)) d , u h = g on Γ u h , u h (x, t = 0) = u 0 (x)} (where Γ u h corresponds to a continuous discretization of Γ u ), H v h = {v h ∈ (F p h ∩ H 0,1 (Q)) d , v h = ġ on Γ u h , v h (x, t = 0) = v 0 (x)}, H u h 0 = {u h ∈ (F p h ∩ H 1,1 (Q)) d , u h = 0 on Γ u h , u h (x, t = 0) = 0}, H v h 0 = {v h ∈ (F p h ∩ H 0,1 (Q)) d , v h = 0 on Γ u h , v h (x, t = 0) = 0}.
We first define a standard Galerkin Least Square (GLS) form in a similar manner to what was done by [START_REF] Franca | The galerkin gradient least-squares method[END_REF][START_REF] Hulbert | Space-time finite element methods for second-order hyperbolic equations[END_REF]. The system of eq. ( 6) is modified so that:

Find (u h (x, t), v h (x, t)) ∈ H u h × H v h , such that ∀ (δu h (x, t), δv h (x, t)) ∈ H u h 0 × H v h 0 Q h ρ vh δu h dQ + Q h σ : ε(δu h )dQ - P h T h δu h dP - Q h f h δu h dQ + e Qe (ρü h -div x σ -f h )ρ -1 τ (ρ δu h -div x (C : ε(δu h ))dQ + Q h ρ (v h -uh ) δ v h dQ = 0 ( 11
)
where τ is a stabilization (numerical) parameter and C is the 4th order elasticity tensor. Hence, this form can be easily implemented in the linear case, but not as easily in the non-linear case. We propose a new alternative form based on re-injecting the strong form of the first momentum conservation, ρ v = div x σ + f , into the least square term. We obtain the following form:

Find (u h (x, t), v h (x, t)) ∈ H u h × H v h , such that ∀ (δu h (x, t), δv h (x, t)) ∈ H u h 0 × H v h 0 Q h ρ vh δu h dQ + Q h σ : ε(δu h )dQ - P h T h δu h dP - Q h f h δu h dQ + e Qe ρ(ü h -vh )τ ( δu h -δ v h )dQ + Q h ρ (v h -uh ) δ v h dQ = 0 (12)
This form does not cause any difficulties in the implementation of non-linear cases. We choose to call it Galerkin with Acceleration Consistency (GAC) condition. Both stabilization strategies share the main idea of adding a mesh-dependent perturbation term. The proof of consistency for both formulations is straightforward and does not pose any difficulties. Some evidence of stability for the first formulation can be found in the literature, e.g. [START_REF] Hulbert | Space-time finite element methods for second-order hyperbolic equations[END_REF]. In the following we study the convergence properties of these formulations for both IGA and FE space-time approximations from a numerical point of view.

Finite-strain elastodynamics

Problem statement

In this section, we distinguish between the current domain Ω(t) and the initial domain Ω 0 . A material point in Ω 0 is denoted by X, and by x in Ω(t). The transformation, φ(X, t), makes it possible to relate between both configurations, so that x = φ(X, t). We assume that φ is a bijective function, and so we can always switch from a current description to a description in the initial configuration. In the following we adopt a Lagrangian representation of the problem. The space time cylinder is defined such that Q = Ω 0 × [0, T ]. Therefore, time and space integration commute, so that Q dQ = T 0 Ω0 dQ = Ω0 T 0 dQ. The conservation of the linear momentum on Q is similar to eq. ( 1):

ρ 0 ü -DIV X Π = f (13)
where ρ 0 (X) is the initial material density, Π is the first Piola-Kirchoff stress, DIV X is the Lagrangian divergence operator and f is the body force defined by unit of initial volume. The dot superscript stands for the so called material time derivative (derivative with respect to time holding X constant), such that for a Lagrangian field u(X, t) = (∂u(X, t)/∂t) X .

For the constitutive behavior, we restrict ourselves to the case of isotropic hyperelastic materials and we distinguish between two cases: compressible materials and nearly-incompressible materials. For the first case, the stress/strain relation is simply obtained from the derivative of a specific free energy potential ψ(F):

Π = ρ 0 ∂ψ ∂F ( 14
)
where F is the transformation gradient:

F = ∇ X φ = ∇ X u + 1.
For a nearly incompressible behavior, we adopt a formalism based on a hybrid free energy, as described in [START_REF] Lejeunes | Hybrid free energy approach for nearly incompressible behaviors at finite strain[END_REF]. This formalism is similar to a partial Legendre transformation of the Helmoltz energy formulated with a volumetric part. It makes it possible to introduce the hydrostatic pressure, p, as a thermodynamic state variable and naturally leads to a displacement/pressure formulation of the conservation equations.

The first Piola-Kirchoff stress is therefore defined by:

Π = Π : P F + JpF -T (15)
where J = det(F), P F is a 4th order tensor that corresponds to the deviatoric operator in the Lagrangian configuration, Π = ρ 0 ∂ϕ/∂ F where ϕ is a hybrid free energy potential (see [START_REF] Lejeunes | Hybrid free energy approach for nearly incompressible behaviors at finite strain[END_REF]) and F = J -1/3 F. The additional field p(X, t) stands for the hydrostatic pressure and we need to define a supplementary constitutive equation that relates p to the volume variation J. For the sake of simplicity, we consider the simplest form:

J -1 = g(p) ( 16
)
where g is a compressibility function.

The problem is fully defined by taking into account the same boundary and initial conditions as defined in eqs. (3).

Variational formulations for compressible and nearly-incompressible problems

For the compressible problem, the variational formulation is very similar to the one of eq. ( 5):

Find (u(x, t), v(x, t)) ∈ H u × H v , such that ∀ (δu(x, t), δv(x, t)) ∈ H u 0 × H v 0 Q (ρ 0 v -DIV X Π -f ) δ udQ + Q ρ 0 (v -u) δ vdQ = 0 (17)
By integrating by part the divergence operator and adding a least square term, the stabilized version of the previous equation in its discrete form can be written as:

Find (u h (x, t), v h (x, t)) ∈ H u h × H v h , such that ∀ (δu h (x, t), δv h (x, t)) ∈ H u h 0 × H v h 0 Q h ρ 0 vh δu h dQ + Q h Π : ∇ X δu h dQ -β P h T h δu h dP -β Q h f h δu h dQ + e Qe ρ 0 (ü h -vh )τ ( δu h -δ v h )dQ + Q h ρ 0 (v h -uh ) δ v h dQ = 0 (18)
here T h = ΠN where N is the normal to the initial contour Γ. A load parameter β ∈ [0, 1] is introduced here to control the iterative convergence of the numerical scheme. As for the small strain case, we recover the principle of virtual power in the first part of the system of eq. ( 19).

For the nearly-incompressible formulation, we propose to adopt a three field formulation and we have:

Find (u h (x, t), v h (x, t), p h (x, t)) ∈ H u h × H v h × L 2 (Q), such that ∀ (δu h (x, t), δv h (x, t), δp h (x, t)) ∈ H u h 0 × H v h 0 × L 2 (Q) Q h ρ 0 vh δu h dQ + Q h ( Π : P F + p h JF -T ) : ∇ X δu h dQ -β P h T h δu h dP t -β Q h f h δu h dQ + e Qe ρ 0 (ü h -vh )τ ( δu h -δ vh )dQ + Q h ρ 0 (v h -uh ) δ v h dQ + Q h (J(u h ) -1 -g(p h )) δp h dQ = 0 (19)
4. Space-time discretization

About Isogeometric Analysis

In the following, we briefly recall some fundamental aspects of the Isogeometric Analysis (IGA). For further details we refer the interested reader to the work of [START_REF] Hughes | Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement[END_REF][START_REF] Cottrell | Isogeometric analysis: toward integration of CAD and FEA[END_REF] and references therein. We use the NURBS2 version of the IGA. NURBS are based on B-Splines which are polynomial functions. These functions are defined from a knot vector (a set of points in a parametric space). For instance, in one dimension, if E = {ξ 1 , ξ 2 , ..., ξ n+p+1 } is a knot vector (where p is the polynomial degree and n is the number of basis functions), B-Spline functions are recursively defined by the Cox-de-Boor formula:

For p = 0 N i,0 (ξ) = 1 if ξ i ≤ ξ < ξ i+1 0 otherwise For p > 0 N i,p (ξ) = ξ -ξ i ξ i+p -ξ i N i,p-1 (ξ) + ξ i+p+1 -ξ ξ i+p+1 -ξ i+1 N i+1,p-1 (ξ) (20)
Knots can be repeated at the same position and the number of repetitions of a knot is called the order of multiplicity of this knot. The continuity order of the approximation functions is equal to the order of the functions between knots and is equal to the order minus the order of multiplicity at knots. The NURBS functions are defined by:

R i,p (ξ) = N i,p (ξ)ω i n l=1 N l,p (ξ)ω l (21)
where ω i is the i th weight that can be viewed as a supplementary coordinate that controls the approximation. The NURBS functions allow conic sections to be exactly represented.

The space time cylinder can therefore be described, for simple cases 3 , by a parametric surface or volume for a 1D or a 2D case. For a 1D problem, the space time cylinder is defined by the following parametric surface:

t(ξ, η) = n1 i=1 n2 j=1 R i,p (ξ)R j,p (η)B t i,j x(ξ, η) = n1 i=1 n2 j=1 R i,p (ξ)R j,p (η)B x i,j (22)
where B i,j = {B t i,j , B x i,j } are called control points (similar to nodes' dofs for finite elements), see Figure 1 for a synthetic view of the description of a space time domain within the IGA for a 2D case. All the fields (including virtual fields) are approximated with the same basis as the one used for describing the space-time domain. For instance for a 1D problem in space, the discretized kinematic and velocity fields are:

u h (ξ, η) = n1 i=1 n2 j=1 R i,p (ξ)R j,p (η)d u i,j v h (ξ, η) = m1 i=1 m2 j=1 R i,q (ξ)R j,q (η)d v i,j (23)
where d i,j = {d u i,j , d v i,j } are the values that control the approximation of the fields at control points (which are not interpolant in general). Discretization for 2D and 3D cases can be derived in the same manner. Remarks 1. The main concept of Space-Time IGA (or Space-Time FE) does not assume, a priori, any separations of space and time due to the mapping between the space-time domain and the parametric space. This mapping plays a fundamental role in these methods. 2. In this work, we only consider cases where we use the basis functions of same order in all the parametric directions. The tensorial structure of the approximation makes it possible, for example, to consider a higher polynomial order in a specific parametric direction. This can be relevant to very specific cases (e.g., when the jacobian operator between the parametric space and the space time cylinder is diagonal). 3. For the sake of simplicity we assume in the following that the displacement and the velocity fields are described with the same approximation basis (i.e., p = q and m i = n i , ∀i = 1 . . . d). As for the FE method, the IGA in a multifield context makes it possible to control the meshing (subdiscretization technique) and the order of approximation of each field. We have additional possibilities, we can also control the order of continuity of each field at control points, which we cannot do in the case of FE.

Matrix form of space/time problems

Using the matrix/vector representation of the approximations of eq. ( 23):

u h (ξ, η) = N p d u v h (ξ, η) = N p d v (24)
where N p is a matrix obtained from the approximation functions and d u , d v are the vectors obtained from the unknown values at control points. We can define the following differential operators (that operate on real or test fields):

uh (ξ, η) = B t d u ∇ x u h (ξ, η) = B x d u ε(u h (ξ, η)) = B sym x d u üh (ξ, η) = B tt d u ∇ x uh (ξ, η) = B xt d u ε( uh (ξ, η)) = B sym xt d u (25)
These operators require the computation of the second derivatives of the approximation functions (we therefore need p ≥ 2). We also need to compute the space-time jacobian matrix (for a 1D space problem):

J(ξ, η) =    ∂t ∂ξ ∂x ∂ξ ∂t ∂η ∂x ∂η    (26)
With the operators defined in eq. ( 25) and using the mapping of a space-time element to the parametric space, the previous weak formulations can be rewritten without difficulties into a matrix form. For instance, the formulation of eq. ( 12) can be written as (linear case):

K uu K uv K vu K vv d u d v = f u f v (27)
The matrix and vector introduced in eq. ( 27) are defined by:

K uu = A e Qe B sym T xt DB sym x + ρB tt T τ B tt dQ K uv = A e Qe ρB t T B t -ρB tt T τ B t dQ K vu = A e Qe -ρB t T B t -ρB t T τ B tt dQ K vv = A e Qe ρB t T N p + ρB t T τ B t dQ (28)
and:

f u = A e Qe B t T f h dQ + S e Pe B t T T h dP f v = 0 (29)
The assembly operator on elements is denoted by A e and by S e for assembly on element surfaces (or lines) on the boundary of the space domain (i.e. x = 0 or x = L for a 1D problem). Remarks 1. As mentionned previously, the previous linear system does not exhibit any symmetry properties on its left hand side. 2. Initial and kinematic/velocity boundary conditions are prescribed directly using the values at control nodes located on the boundary (we should pay attention to the non interpolant behavior of B-Spline functions on boundaries when g(x, t) and T(x, t) are not constant functions).

For a non-linear problem, such as the one of eq. ( 18), we obtain a non-linear system to be solved in the vector form:

r u (d, β) = 0 r v (d, β) = 0 ( 30
)
where

d = [d u , d v ]
T and (r u , r v ) are defined by:

r u (d, β) = A e Qe ρ 0 B t T vh + B T xt Π(u h ) + B tt T ρ 0 τ (ü h -vh ) dQ -βf u r v (d, β) = A e Qe ρ 0 B t T ((v h -uh ) + τ (ü h -vh )) dQ (31)
and β ∈ [0, 1] is a load factor. The procedure used to solve this non-linear system is very standard and completely identical to a predictive/corrective Newton-Raphson procedure for a quasi-static problem. Starting from a known solution, (

d k , β k ), at increment k, we search ∆d k = [∆d u k , ∆d v k ] T for a given load increment ∆β k such that r u (d k + ∆d k , β k + ∆β k ) = r v (d k + ∆d k , β k + ∆β k ) = 0.
The linearization of this non-linear problem leads to:

K uu (u k h , v k h ) K uv (u k h , v k h ) K vu (u k h , v k h ) K vv (u k h , v k h ) ∆d u k ∆d v k = -r u (d k , β k ) -f u ∆β k -r v (d k , β k ) (32)
We distinguish the predictive and the corrective phases. For the predictive phase, the previous solution is converged and therefore r u (d k , β k ) = 0, r v (d k , β k ) = 0 and ∆β k = l where l is a prescribed value.

For the corrective phase, the loading is frozen: ∆β k = 0 and r u (d k , β k ) = 0, r v (d k , β k ) = 0. We iterate with (d k+1 , β k+1 ) = (d k + ∆d k , β k + ∆β k ) until convergence to a given tolerance. The matrix terms K vu , K vv , K uv are identical to those defined previously at eq. ( 28), the term K uu is slightly different and is defined as follows:

K uu = A e Qe B T xt D nl B x + ρB tt T τ B tt dQ (33)
where D nl = ∂Π/∂F.

The specific case of nearly-incompressible problems

The three field formulation presented at eq. ( 19) is well adapted for nearly incompressible situations but we need to pay attention to numerical instabilities or oscillations that may occur if we do not carefully choose the approximation spaces for the pressure and displacement. For isogeometric analysis, different propositions can be found in the literature. We choose here the one proposed in [START_REF] Kadapa | Subdivision based mixed methods for isogeometric analysis of linear and nonlinear nearly incompressible materials[END_REF]. It consists in using the same order of approximations for the pressure and the displacement fields but with a grid subdivision for the pressure field (i.e. , a coarser mesh for pressure than the one used for velocity/displacement).

Using the same interpolation order for displacement field and pressure field mostly leads to unstable finite elements, but, because of the subdivision properties of B-Spline, it is not the case for IGA. It was shown in [START_REF] Kadapa | Subdivision based mixed methods for isogeometric analysis of linear and nonlinear nearly incompressible materials[END_REF] that this formulation satisfies the inf-sup (or Ladyzhenskaya-Babuska-Brezzi condition) for small strain formulations and this property seems conserved for finite strain formulations, at least for simple problems. In any case, we did not experience any typical pressure instability or volumetric locking with this formulation.

In this work, for ST-IGA, we have chosen to use a coarser mesh size for the pressure only for space discretization. Figure 2 presents a synthetic view for a 2D space problem. The numerical integration is done on the finest mesh for all fields using a standard Gauss-Legendre scheme. We consider the case of a unit bar (1D space problem) with homogeneous boundary conditions in space and time for the displacement and the velocity, subjected to body load f (x, t) defined by:

f (x, t) = 4π 2 sin(2πx) E log 1 + cos 2 (2πt) 2 -4ρ 1 + 3 cos(4πt) (3 + cos(4πt)) 2 (34)
The exact solution is therefore:

u(x, t) = sin(2πx) log 1 + cos 2 (2πt) 2 v(x, t) = - 4π sin(4πt) sin(2πx) 3 + cos(4πt) (35)
The space and time domain are chosen to be [0, 1] and so the space-time domain is a unit square. Figure 3 illustrates the test considered. The material parameters are: E = 2 and ρ = 1. This test makes it possible to establish a convergence study while considering both space-time FE and space-time IGA. For the finite elements, we use quadratic quadrangle with 9 nodes. We compute the following L 2 errors on the space-time domain Q:

er U = Q u h (x, t) -u(x, t) 2 dQ er V = Q v h (x, t) -v(x, t) 2 dQ (36)
where (u h , v h ) is the approximated solution and (u, v) are given by eq. [START_REF] Kuraishi | Space-time computational analysis of tire aerodynamics with actual geometry, road contact, tire deformation, road roughness and fluid film[END_REF]. We first plot the results without any stabilization (τ = 0). We use a uniform mesh refinement (in space and time) and we plot the L 2 error with respect to the element size. Figure 4 shows that we obtain a monotoneous convergence behavior and by computing the rate of convergence for different polynomial orders (from 2 to 4 for IGA), we can see that we have an optimal rate of convergence for both ST-FE and ST-IGA. If we compare ST-FE and ST-IGA for the same polynomial order (p=2), we can see that there are no differences between the quadratic finite elements and the quadratic NURBS elements for this test.

Figure 5 shows a comparison of the results for different stabilization strategies. It can be noticed that all formulations are convergent. The convergence rate is affected in a similar manner with both stabilization strategies (GAC or GLS) with ST-FE or ST-IGA. For the same value of the stabilization parameter, we obtain a stronger influence on approximations of order 2 while approximations of order 3 do not seem much affected by the additional stabilization terms.

This simple convergence study illustrates the interest of using higher polynomial orders for space-time discretizations. It also confirms that GLS and GAC stabilization strategies behave in a similar manner for both ST-FE and ST-IGA, at least for this simple test.

Impact of an elastic bar

We consider here the case of the impact of a 1D elastic bar against a rigid wall. This test was first investigated with ST-FE in [START_REF] Hughes | Space-time finite element methods for elastodynamics: Formulations and error estimates[END_REF] and it is a typical benchmark test to investigate the behavior of numerical oscillations for a wave propagation problem using these numerical methods. The problem is illustrated in Figure 6(a). A 1D bar is subjected to an initial homogeneous velocity v 0 . The impact against the rigid wall is taken into account with displacement and velocity conditions on one side of the bar (here x = 0). These Using ST-IGA or ST-FE without any stabilization, we obtain the typical stress solution in the space-time domain shown in Figure 6(b). Spurious numerical oscillations are obtained regardless the mesh size or the polynomial order of the approximations. These oscillations propagate over time and solving the problem for a longer period of time will result in increasingly noisy solutions.

Plotting the stress along the bar at t = 0.5 and using unstabilized formulations (i.e., for τ = 0), we clearly observe numerical oscillations in the part of the bar subjected to compression. Figure 7(a) shows a comparison of the results between ST-FE and ST-IGA using the same mesh size, and we observe the same behavior for both methods. However, as can be seen from the plot of the total energy (see Figure 7 Numerical oscillations can be controlled by adding stabilization terms, as can be seen in Figures 8(a) and 8(b). First, we can notice from Figure 8(a) that ST-FE and ST-IGA with GLS stabilization of order 2 behave in the same way. As expected, increasing the order of NURBS while keeping other parameters (mesh size and stabilization parameter) constant leads to a more accurate result in terms of localization of the compression wavefront. Second, if we compare GAC and GLS stabilizations, we can see from Figure 8(b) that the GAC stabilization seems to be less dissipative than GLS for the same order of approximation, both for ST-IGA and ST-FEM. This point is more evident in Figure 8(c) where a comparison of the total energy integrated over the bar upon time is plotted: dotted lines (GAC stablization) are always over solid lines (GLS stabilization). It is also important to notice that by zooming on the compressive part (Figure 8(d)), small oscillations are still present with ST-FEM but not for ST-IGA.

Another interesting aspect is to investigate the behavior of stabilization strategies for different values of the stabilization parameter τ . Figure 9 shows that, as expected, the numerical dissipation is directly related to the value of τ . A good option with ST-IGA of order 2 seems to choose τ approximately in the same order of magnitude of the meshsize for the time because using a larger value leads to a strong numerical dissipation without improving the numerical solution. We can also notice that the GLS stabilization with IGA can lead to a non-monotonic evolution of the total energy for large values of the stabilization parameter.

The impact of the mesh size ratio is also an important point for space-time methods. We define h = ∆t/∆x as the ratio of the time increment to the space discretization parameter. Figure 10 shows the results obtained with ST-IGA for GLS and GAC stabilization parameters using a constant size of the space discretization and different time increments. We obtain a correct description of the compression wave using GAC and GLS only until h = 4. Higher values of the mesh ratio lead to oscillations whatever the value of the stabilization parameter is. Nevertheless, we can obtain satisfactory results for moderate mesh size ratios without being limited by the CF L ≤ 1 condition where CFL is the Courant-Friedrichs-Lewy number (CF L = h here).

This example illustrates some of the principal characteristics of space-time methods that were also found by other authors in the literature. But these results are obviously not general and should be confirmed by other numerical tests.

Two bar impact

This test is very similar to the previous one. We consider the case of an elastic bar with an initial velocity v 0 = 1 that impacts a second elastic bar with different material properties. The boundary conditions and the dimensions of the bars are illustrated in Figure 11. The material properties of the first bar (colored in blue) are: ρ = 1, E = 1, and those of the second bar (colored in red) are: ρ = 1, E = 5.5. The solution just after the impact corresponds to two compressive waves that propagate from the interface (impact zone) with two different velocities. After being reflected on both sides, these waves interact with the bi-material interface (with reflexion and transmission) on the one hand and with each other on the other hand. Therefore, we obtain, after a few interactions, several waves that propagate with constant but different velocities (Figure 11 shows this typical solution obtained with stabilized ST-IGA).

For this test, we compare the results obtained with ST-IGA stabilized using the GAC formulation (for this example the stabilization parameter was fixed such that τ = (∆t) 1/(p-1) where p is the polynomial order of the approximation) for different polynomial orders from 2 to 4 with a standard finite element model on Abaqus. For this model, we used quadratic truss elements with the implicit resolution. The time integration scheme is the Hilber-Hughes-Taylor scheme (see [START_REF] Hilber | Improved numerical dissipation for time integration algorithms in structural dynamics[END_REF]) with the following parameters α = -0.05, β = 0.275, γ = 0.55, which corresponds to a transient fidelity integration. For all models, we used the same space-time discretization with ∆t = ∆x = 0.033. Figure 12 shows the stress along the bar at different times. It can be seen that the finite element models exhibit strong oscillations, unlike what happens in the case of the ST-IGA solutions. Increasing the order of approximation both in space and time for the ST-IGA allows better solutions (more conservative) to be obtained. This can be clearly seen in Figure 13, ST-IGA schemes of order 3 and 4 are more conservative than the ones obtained with the FE model on Abaqus.

(c) t = 2.495 -1.2 -1 -0.8 -0.6 -0.4

Extension of the bar test for truss structures

In this section we show that the previous formulation concerning 1D bar problems can be extended to the case of a bar in space for applications to truss structures. We assume that each truss in a structure is described by a space-time patch. Each patch has a local coordinate system with 2 dimensions in space, (e x , e y ), which is related to a global coordinate system (e x , e y ) with the following rotation matrix:

e x e y = cos(θ) sin(θ) -sin(θ) cos(θ) e x e y = R(θ) e x e y (37)
The axis e x corresponds to the direction of the truss and e y is normal to the truss. We denote by (u x (x, t), u y (x, t)) and (v x (x, t), v y (x, t)) the displacement and velocity of the truss in a patch. We also assume that the normal velocity is constant and the normal displacement is only a function of time, such that:

v y (x, t) = cst u y (x, t) = y(t) (38
∇ x u y = 0 (39)
The uniaxial strain and stress in the truss are defined by:

ε(x, t) = ∇ x u x (x, t) σ(x, t) = Eε(x, t) (40)
The weak form for this problem on a patch is simply obtained from an extension of eq. ( 12). We first assume that the angle θ is constant (not time dependent) and we add two supplementary terms (least-square like terms) that correspond to the constraints defined in eq. ( 39):

Find (u h (x, t), v h (x, t)) ∈ H u h × H v h , such that ∀ (δu h (x, t), δv h (x, t)) ∈ H u h 0 × H v h 0 Q h ρ(R 1 vh )(R 1 δu h )dQ + Q h σ : ε(R 1 δu h )dQ - P h (R 1 T h )(R 1 δu h )dP - Q h (R 1 f h )(R 1 δu h )dQ + e Qe ρ(R 1 üh -R 1 vh )τ (R 1 δu h -R 1 δ vh )dQ + Q h ρ (R 1 v h -R 1 uh ) R 1 δ vh dQ + e Qe k 1 (R 2 vh )(R 2 δ vh )dQ + e Qe k 2 (∇ x R 2 u h )(∇ x R 2 δu h )dQ = 0 (41)
where R 1 and R 2 correspond respectively to the first and second row of R(θ), such that u x = R 1 u h and u y = R 2 u h , k 1 , k 2 are penalty parameters. For a multipatch configuration we simply proceed as in the previous example, two connected patches share a common interface with the same refinement and the same control points.

As an example of application, we consider the case of an assembly of two trusses that are connected with a given angle θ 0 . Each truss is of unit length with the same elastic material whose characteristics are: E = 1, ρ = 1. The trusses have an initial vertical velocity, v 0 = 1, and impact a wall on one side (we assume that the trusses remain fixed onto the wall after impact), see Figure 14.

Figure 15 shows the results obtained for two configurations: θ 0 = 20 o and θ 0 = 90 o . The mesh size is chosen such that ∆t = ∆x = 0.02. The stabilization parameter and the penalty conditions are taken to be , k 2 = 10 9 . It can be observed that in the situation where θ 0 = 90 o , the upper truss moves vertically with a constant velocity but without stress, as expected. The impact wave is not transmitted to the upper truss in this case, unlike what happens in the case where θ 0 = 20 o . This simple example illustrates the potential of space time methods for application to wave propagation problems in structures. Furthermore they could be extended to cases where the rotation matrix R is not constant on time to deal with the case of truss structures with pivot connection between trusses. Obviously, beams, shells and other structural elements can be developed in this framework.

6. Numerical applications at finite-strain 6.1. Nearly-incompressible plane-strain structure subjected to harmonic loading

We consider the case of a non-linear beam like structure (2D plane-strain case) of 5 units of length in height and 1 unit of length in width (see Figure 16(a)). The beam is fixed on one side and subjected to a harmonic homogeneous body load f (t) = f 0 cos(2πF t). The constitutive material is chosen to be hyperelastic and we adopt a Mooney-Rivlin potential:

ρ 0 ϕ = c 10 (I 1 ( C) -3) + c 01 (I 2 ( C) -3)
where C = J -2/3 C = J -2/3 F T F, I 1 ( C) = tr( C), I 2 ( C) = 1 2 ((tr( C)) 2 -tr( C2 )). For the compressibility part, we choose the following compressibility function:

g(p) = exp p k - 1
The material parameters are: c 10 = 0.9, c 01 = 0.1, k = 1000, ρ 0 = 0.001. From a linear vibration analysis, we determine the first three linear modes of this structure. The frequencies of these modes are: F 1 = 0.56, F 2 = 2.92, F 3 = 4.5. In the following we consider a harmonic body load with a frequency chosen between the first and the second modes of the structure: F = 1.3 with an amplitude of f 0 = 0.0184. As for the previous example, we compare the results of ST-IGA with the results obtained on Abaqus using implicit and explicit solvers. For the space domain, we used the following mesh size: ∆x = ∆y = 0.25. For the implicit solver, we consider quadratic hybrid elements (CPE8H) that are adapted to the case of near-incompressible behavior.

For the explicit solver, we are restricted to use linear elements with reduced integration (CPE4R). For the ST-IGA model, we use the GAC form with a stabilization parameter value of τ = ∆t/( E/ρ), where E = 6(c 10 + c 0.1 ) is the infinitesimal Young modulus and ∆t = 0.14 is the mesh size chosen for the time.

A typical response of the ST-IGA model for the time interval, [0, 10], is shown in Figure 16(b). Figure 17 shows a comparison between the different numerical methods for the displacement and the velocity at a node located at the top of the structure (the red dot in Figure 16(a)). As expected, the solutions obtained using ST-IGA and FE with implicit resolution are almost identical, while the explicit resolution exhibits some differences that are certainly due to the fact that only linear elements are used.

Figure 18 illustrates the main advantage of ST-IGA methods for non-linear problems: we can obtain more conservative solutions than with standard implicit integration schemes such as the HHT scheme. Table Figure 19: Impact of a plane strain structure on a wall at finite strain 1 allows the numerical characteristics of the different methods to be compared. For ST-IGA, we obtain the solution on the time slab with only 6 newton iterations (for only one increment for this test) while the implicit resolution requires 250 increments with 1 or 2 newton iterations per increment. Obviously, there is a huge difference between space-time continuous formulations and other methods in terms of the size of the linear system to be solved at each iteration. The purpose of this test is to evaluate the numerical behavior of the ST-IGA method in a context where we have a complexity due to both the dynamics (impact situation and non-linear behavior) and to the geometry. The boundary conditions and the geometry of the structure are given in Figure 19. This plane strain structure has an initial velocity, v 0 = 5 (in the vertical direction), and at t = 0 it impacts a wall at the top of the structure (we prescribe the separation of the top surface of the structure from the wall). The material is assumed to be elastic and isotropic and we adopt the following potential of energy:

ρ 0 ψ = µ 2 (I 1 (C) -3) + λ 2 ln(J) 2 -µ ln(J) (42)
with λ = 0.52, µ = 0.34 and ρ 0 = 0.001. As in the previous example, we compare the results obtained with standard implicit and explicit finite elements on Abaqus with the results of space-time methods (FE or IGA based) stabilized with GAC. On Abaqus the behavior law is defined using the subroutine UHYPER. Therefore we only use the implicit version of Abaqus with transient-fidelity resolution (HHT scheme with α = -0.05, β = 0.275, γ = 0.55). We use the same mesh size for the discretization of the space domain: ∆x = ∆y = 0.036 for all models.

Figure 20 shows a comparison of the methods for solutions concerning the displacement and velocity components of a point located at the bottom of the structure (red point in Figure 19(a)). First, it can be seen that the response obtained with Abaqus exhibits many oscillations of the velocity. The solution obtained using ST-FE is globally in accordance with the one obtained with Abaqus implicit. ST-IGA exhibits a similar response immediately after impact but deviates from FE responses especially on the vertical component (for both displacement and velocity). Obviously, all these models are questionable and we cannot define a reference solution for this problem, as each model has its own bias.

Nevertheless, we can plot the evolution of the total energy (kinetic energy+strain energy integrated over the space domain) upon time. Figure 21 shows that the evolution of this energy with the implicit solver of Abaqus exhibits a non monotonic behavior (increase and decrease in the total energy). This non-monotonous behavior is obviously completely inconsistent: we can expect a loss of energy due to numerical dissipation which is intrinsic to the HHT scheme but we do not expect an energy increase. Space-time methods have a more expected behavior with a monotonic decrease in energy due to the contribution of stabilization terms. It can also be noted that ST-FE exhibits a higher rate of decrease than ST-IGA. Therefore, from an energy point of view, space-time methods behave better than Abaqus Implicit for this test. Table 2 compares the numerical characteristics of each model. As in the previous example, it is obvious that ST methods represent a very interesting alternative for solving non-linear dynamics problems. Furthermore, they offer the possibility of developing of high-performance resolution methods based on space-time domain decomposition strategies and parallelized resolution for time dependent problems. ∆t Number of increments Newton Iterations ST-IGA (p=2,3) 0.01 2 5 for the first increment and 1 for the second ST-FE (p=2) 0.01 2 5 for the first increment and 1 for the second FE Implicit (HHT) 0.001 212 3 to 6 per increment

Discussion and concluding remarks

We have developed new space-time isogeometric formulations for elastodynamics problems. These formulations are based on displacement-velocity (and pressure in case of nearly incompressible behaviors) weak forms. This formalism has made it possible to develop a new stabilization term based on least-squares. This term (GAC) exhibit the main advantage of being independent of the constitutive law and permit both small and finite strains problems to be dealt with. The proposed formulation leads to optimal rates of convergence on trial problems and it is shown that the additional terms play a fundamental role in controlling numerical oscillations that occur in problems with time discontinuities, such as impact problems. Obviously, specific methods that can deal more precisely with time discontinuity problems exist in the literature (some of them are cited in the introduction). However, one main interest of the proposed formulation is that it is general enough to be applied to different problems without the need to develop specific, case-dependent, terms. We have explored different cases: small-strain elasticity, compressible and quasi-incompressible hyperelasticity, and the approach can easily be extended to multi-field and multi-physics problems such as thermo-elasticity.

With regard to the NURBS formulations, we have shown that ST-IGA formulations on the one hand make it possible to obtain more conservative solutions than semi-discrete methods and, on the other hand, offer a richness in the construction of approximation spaces: high-order continuity, sub-division techniques, etc. Obviously the numerical cost of such methods is high, especially for multi-field formulations, but employing high order space-time schemes makes it possible to obtain accurate numerical solutions with coarser mesh in time compared with implicit semi-discrete methods. Therefore, there is several advantages of using NURBS in this space time context. However, in view of considering more complex or more industrial examples several enhancements or developments should be done. The first one could be to deal with multi-pach in space and time with non-conforming meshes using for instance a Nitsche formulation (e.g. [START_REF] Nitsche | Über ein variationsprinzip zur lösung von dirichlet-problemen bei verwendung von teilräumen, die keinen randbedingungen unterworfen sind[END_REF][START_REF] Ruess | Weak coupling for isogeometric analysis of non-matching and trimmed multi-patch geometries[END_REF][START_REF] Nguyen | Nitsche's method for two and three dimensional nurbs patch coupling[END_REF]). The Nitsche formulation has the advantage of allowing a control of the order of continuity between patches for each field of the formulation. This could be interesting for the time dimension in an apriori approach (i.e. having knowledge of the discontinuities in the time domain). Another alternative is to consider one patch with local refinements methods in space and time. This can be done using various methods, such as: T-splines (e.g. [START_REF] Scott | Local refinement of analysis-suitable t-splines[END_REF]), PHT-splines (e.g. [START_REF] Nguyen-Thanh | Rotation free isogeometric thin shell analysis using pht-splines[END_REF]), LRB-splines (e.g. [START_REF] Dokken | Polynomial splines over locally refined box-partitions[END_REF]), the GIFT method (see [START_REF] Atroshchenko | Weakening the tight coupling between geometry and simulation in isogeometric analysis: From sub-and super-geometric analysis to geometry-independent field approximation (gift)[END_REF]), etc. Local refinement methods could be interesting if used in an adaptative process. However, this also requires the development of apriori estimators which have the inconvenient of being problem dependent. The second aspect is related to the computational cost of the quadrature method, some authors have already proposed alternatives to the standard Gauss-Legendre quadrature, and this could be very interesting in a space-time context because we have an extra dimension (e.g. [START_REF] Auricchio | A simple algorithm for obtaining nearly optimal quadrature rules for nurbs-based isogeometric analysis[END_REF]). Space-time formulations are also well adapted for the development of parallel resolution in space and time and references on this topic can be found in the literature (e.g. [START_REF] Hesch | Variational space-time elements for largescale systems[END_REF][START_REF] Langer | Space-Time Isogeometric Analysis of Parabolic Evolution Equations[END_REF]).

Finally, we see an open field of research for ST-IGA methods, that we will explore, for solving non-linear and non-conservative transient problems especially in the case of history dependency of the constitutive behavior (viscosity, damage, viscoplasticity, etc). In these cases, the set of equations leads to a DAE system for which a local-global time integration scheme is the most currently employed strategy, which can cause many difficulties in certain situations (thermo-mechanical coupling, etc).

Figure 1 :

 1 Figure 1: Space-Time cylinder described with B-Spline functions

Figure 2 :Figure 3 :

 23 Figure 2: Mesh subdivision technique for the pressure field (3 field formulation for nearly incompressible problems)

Figure 4 :

 4 Figure 4: Convergence study for the bar test without stabilization

Figure 5 :Figure 6 :

 56 Figure 5: Convergence study for the bar test with stabilization (p = 2)

  (b)), non-stabilised space-time methods are conservative.

Figure 7 :

 7 Figure 7: Results of the bar impact test without stabilization with: ∆t = ∆x = 0.0125

Figure 8 :

 8 Figure 8: Results of the bar impact test with different stabilization strategies with: ∆t = ∆x = 0.0125

Figure 9 :-Figure 10 :

 910 Figure 9: Comparison of stabilization strategies for the bar impact test with ST-IGA of order 2, ∆t = ∆x = 0.0125

Figure 11 :

 11 Figure 11: Impact of 2 elastic bars. The blue bar has an initial velocity v 0 = 1 and it is assumed that the two bars stay perfectly tied after impact. A typical Stress field on the Space-Time domain (ST-IGA with GAC stabilization) is shown.

Figure 12 :

 12 Figure 12: Stress along the two bars for different times after impact. Space and time increments are identical for all models: ∆t = ∆x = 0.033

Figure 13 :

 13 Figure 13: Energy kinectic + strain energy integrated over the two bars upon time

Figure 14 :

 14 Figure 14: Truss impact: Geometry and boundary conditions

  (a) deformed configuration for θ 0 = 20 o (isocolor corresponds to the x component of the displacement field) (b) stress field for θ 0 = 20 o (c) deformed configuration for θ 0 = 90 o (isocolor corresponds to the x component of the displacement field) (d) stress field for θ 0 = 90 o

Figure 15 :Figure 16 :Figure 17 :

 151617 Figure 15: Result of the trusses impact test (ST-IGA with GAC stabilization and p = 2)

Figure 18 :

 18 Figure 18: Energy balance (kinetic energy + strain energy -external power integrated over time) integrated over the space domain

Figure 20 :

 20 Figure 20: Displacement and velocity at the bottom right corner of the structure upon time

Figure 21 :

 21 Figure 21: Energy balance (kinetic energy + strain energy -external power integrated over time) integrated over the space domain

  9

Table 1 :

 1 Comparison of characteristics between incremental and space-time resolutions 6.2. Plane strain structure subjected to an impact

		∆t	Number of increments	Newton Iterations
	ST-IGA	0.14	1	6
	FE implicit (HHT)	0.05	250	1 or 2 per increment
	FE explicit	0.00016	62500	-

Table 2 :

 2 Comparison of characteristics between incremental and space-time resolutions

The concept of isogeometric approximation seems to have been proposed for the first time by[START_REF] Grebennikov | Isogeometric approximation of functions of one variable[END_REF] but for a very different objective from that of the isogeometric analysis.

Non Uniform Rational B-Spline

More complex cases can be treated by using multiple patches (multiple parametric geometries)

Acknowledgments

The project leading to this publication has received funding from Centrale Marseille & Région Sud and the support of Labex MEC (Mechanics and Complexity) an Excellence Initiative of Aix-Marseille University -A*MIDEX, a French "Investissements d'Avenir" program.