D Remeshing Using a Space-Time Finite Element Method for Elastodynamics Problems

In this article, a Space-Time Finite Element Method (STFEM) is proposed for the resolution of mechanical problems involving three dimensions in space and one in time. Special attention will be paid to the non-separation of the space and time variables because this kind of interpolation is well suited to mesh adaptation. For that purpose, we have developed a technique of 4D mesh generation adapted to space-time remeshing. A difficulty arose in the representation of 4D finite elements and meshes. This original technique does not require coarse-to-fine and fine-to-coarse mesh-to-mesh transfer operators and does not increase the size of the linear systems to be solved, compared to traditional finite element methods. Space-time meshes are composed of simplex finite elements. Computations are carried out in the context of the continuous Galerkin method. We have tested the method on a linearized elastodynamics problem. Our technique of mesh adaptation was validated on elementary examples and applied to a problem of mobile loading. The convergence and stability of the method are studied and compared with existing methods. This work is a first implementation of 4D space-time remeshing. A stability criterion for the method is established, as well as a convergence rate of about two. Using simplex elements, it is possible to develop a technique of mesh adaptation able to follow a mobile loading zone.


Introduction
The STFEM (Space-Time Finite Element Method) can be regarded as an extension of the classical finite element method, applied to a boundary problem resulting from a non-stationary problem.Currently, several approaches exist.One can quote for example the Large Time INcrement method (LATIN [1]), the discontinuous Galerkin method [2][3][4] and our method, which is a continuous Galerkin method [5,6].In most publications on the discontinuous Galerkin method, like in [7], the interpolation functions are assumed to be a product of functions of space variables and functions of time variables.We will see in this paper that special attention will be paid to the non-separation of the space and time variables.The reason for this choice is not motivated by the accuracy of the numerical results, but rather by what constitutes the aim of our study: remeshing.We will see that this kind of interpolation is well suited to mesh adaptation.The space-time mesh adaptation we developed is based on a method of mesh generation not structured in space and time.The construction of 4D meshes collides with the limits of representation.To overcome this difficulty, we propose an automatic method of construction inspired by what can be achieved in 2D and 3D.Our technique of mesh adaptation was applied to a problem of mobile load like contact forces.Our approach makes possible the building of an evolutionary mesh able to follow the clamping zone.
Moreover, this technique does not require a mesh-to-mesh transfer operator and allows the preservation of the exact sizes of the linear systems on each space-time slab.
Let us note that one of the drawbacks of the STFEM such as defined in the works of [8,9] is the size of the linear systems to be solved, since it is necessary to solve the full 4D problem at once.The use of a laminated mesh allows us to avoid the assembling of the total matrix of the problem and permits us to consider only submatrices.This drastically reduces the size of the systems to be solved.The size of these linear systems is exactly the same as that obtained in the case of approaches coupling an incremental method of finite differences type to solve time integration, with the "classical" finite element method being used to solve the space problem.
Another large group of methods is based on semi-discretisation, whereby finite elements are used in space and finite differences are used in time.Even if this well-known technique is simpler to use in a classical framework, the remeshing required is expensive due to the necessity of the construction of interpolation/restriction operators between the grids.
The paper is organized as follows.In Section 2, the elastodynamics problem is formulated, and the space-time finite element method is developed.The 4D mesh generation is presented in Section 3, and a paragraph is specially devoted to adaptive mesh refinement.Numerical results are presented and discussed in Section 4.

Principle of the Method
We consider the motion of an elastic body within the small perturbations hypothesis.Let Ω be the set taken up by the body and [0, T] a time interval.The body is submitted to volume force density f d , boundary force density F d on its boundary part ∂ 1 Ω and imposed displacements u d on its boundary part The dynamic problem is: seek the displacement u and the Cauchy stress tensor σ such that: where ρ is the specific mass, ü is the second derivative of the displacement with respect to time, u 0 is the initial displacement, u0 is the initial velocity, a is the Hooke tensor and ε is the infinitesimal linear strain tensor.The aim of this study is to use a finite element method.Then, the previous dynamic problem has to be considered as a boundary problem on the time interval [0, T].For that purpose, as in the cases of the discontinuous Galerkin method [2,10,11] and the LArge Time INcrement (LATIN) method [1,12], the variational formulation is written on the whole space-time domain Ω × [0, T].The variational formulation of the previous boundary problem can be written as follows: Find u ∈ U ad such that: where U ad is the set of displacements, regular enough, which verifies the boundary kinematic conditions and the initial conditions, v is the virtual displacement and U 0 ad is the set of virtual displacements, regular enough, which verify boundary kinematic conditions only.
The first term on the left-hand side of Equation ( 2) is integrated by parts in time in order to determine the first derivative of u and the initial velocity.It gives: where uT (x) is the velocity at time t = T.
The space-time finite element method (STFEM) was firstly proposed in [5,6].Their discretisation used structured space-time meshes obtained as the Cartesian product of spatial elements and a time interval, which is not generally suitable for space-time mesh adaptations.Since then, numerous papers on STFEM have been published.Most of them, like [11,13], deal with the discontinuous Galerkin method in time, but the discretisation also uses structured meshes obtained as the Cartesian product of spatial elements and a time interval.However, the STFEM proposed in [2] has been developed on unstructured meshes.It also employs the discontinuous Galerkin method in time and incorporates stabilizing terms of the least squares type.The space and time discontinuities of all variables are taken into account.In our study, we use a continuous Galerkin method.Classical Lagrange polynomials are used.The finite elements are isoparametric.On a space-time finite element E e (Figure 1), the displacement verifies: where n e is the total number of nodes of the element E e , ϕ e i are the interpolation functions and u e i the nodal displacements.Using matrix notation, one has: u(x, t) = N e (x, t)U e where U e = (u e 1 , ..., u e n e ) T and N e (x, t) = (ϕ e 1 (x, t), ..., ϕ e n e (x, t)) 18 submitted to MDPI the velocity at time t = T.
-time finite element method (STFEM) has been firstly proposed in se structured space-time meshes obtained as the Cartesian product of sp erval, which is not generally suitable for space-time mesh adaptation ers on STFEM have been published.Most of them, like [17,20], Galerkin method in time, but the discretisation also uses structured me n product of spatial elements and a time interval.However, the STFEM eveloped on unstructured meshes.It also employs the discontinuous Ga orporates stabilising terms of least squares type.The space and time d are taken into account.In our study, we use a continuous Galerkin met omials are used.The finite elements are isoparametrics.On a space-time  The same interpolation is used for the virtual displacement v. Then: v(x, t) = N e (x, t)V e where V e = (v e 1 , ..., v e n e ) T Let p be the total number of space-time elements; the previous discretization gives: V T e M e U e (7) where: is the elementary matrix relative to the inertia forces.One can notice that M e is symmetric.Concerning the discretization of the initial and final impulses' contributions, one has: where Λ e is the elementary vector relative to the initial and final impulses.It is defined by: where Ω 0 is the domain at time t = 0 and Ω T is the domain at time t = T.
Similarly, let B e be the matrix such that: the virtual works of internal and external forces are respectively discretized by: and: where the elementary matrix K e relative to internal forces is: B T e aB e dxdt (14) and the elementary vector F e relative to external forces is: N T e F d dsdt.
This space-time discretization leads to the following linear system: where [ M u ] is the assembled matrix relative to the inertia forces, [ K u ] is the assembled matrix relative to the internal forces, {F u } is the nodal vector of external forces, Λ is the nodal vector of impulses and {U} is the nodal vector of displacements.One can note that the matrices [ M u ] and [ K u ] are symmetric.In order to have band matrices, and because we have in mind making computations incrementally in time, the meshes are built to be stratified in time, as in Figure 2.Moreover, the node numbering is conducted in such a way that all nodes in a same stratum have close numbers, then the left-hand side of the system (16) verifies.
metric.In order to have band matrices, and because we have in mind making comput entally in time, the meshes are built to be stratified in time, as in figure 2.Moreover, the ring is conduced in such a way that all nodes in a same stratum have close numbers, th d side of the system (16) verifies Using the space-time mesh described in Figure 2, one has: With this numbering, the total matrix [T] and the sub-matrices [T ij ] are band matrices. Comments: • Choosing a Lagrange interpolation for displacements implies that displacements are continuous, but the velocities are discontinuous.As a consequence, integration by parts as in ( 3) is not totally rigorous, and it could be necessary to use the discontinuous Galerkin formulation, which amounts to writing the derivative of velocity within the theory of distributions.We will preserve the formulation in (3), knowing that the error here is of the same order as in the case of traditional finite elements in space.Indeed, with a Lagrange interpolation local displacements are continuous, whereas the global deformation is discontinuous.

•
Even if it is not absolutely necessary, the advantage of using a laminated mesh such as defined here is that it becomes possible, rather than assembling the total matrix [T], to only assemble the sub-matrices [T ij ].This considerably reduces the size of the systems to be solved.More precisely, the size of these linear systems is exactly the same as that obtained in the case of approaches based on the coupling of finite incremental differences in time with finite elements in space.Moreover, the method is not limited to simplex elements, and the spatial position of each set of nodes can vary from one time plane to the other.It is one of the main advantages of the method.

•
We specify that the nodal vector relating the boundary conditions with velocity {Λ} is written as: where {Λ 0 } is given starting from conditions of initial velocity while {Λ n } is unknown.Consequently, the resolution of System ( 16) is the following: The first system of equations, provides {U 1 }; the system of equations: provides the displacements {U i }; and the last system of equations, gives {Λ n }. • Finally, the matrices of resolution [T i/i+1 ] are generally non-symmetric, even if the total matrix [T] is symmetric.Thus, for the algorithm presented above, a non-symmetrical solver should be used.This can appear penalizing in terms of computing time.However, since the final objective is to use this approach to deal with problems of contact with friction and since the nonlinear resolution we developed in [14] is of the Gauss-Seidel nonlinear type, asymmetries do not affect computing time.

4D Mesh and Remeshing
In order to propose a remeshing technique, it is firstly necessary to be able to build 4D meshes.Obtaining only one 4D finite element does not pose real problems, even if some difficulties in graphic representation arise (see Figure 3).On the contrary, building a 4D mesh, even the most elementary is far from being commonplace, except in the case of regular meshes formed by finite elements of multiplexing type (functions of interpolation obtained as the products of functions of space by functions of time).However, in the general case and in particular with the problem of remeshing, which is what is of interest here, the meshing remains an issue.

4D Mesh Generation
Figure 3 identifies parts of elementary 2D, 3D and 4D meshes with their node numbering.One denotes by n 0 the total number of nodes at time t = 0 of the entire space mesh, and we assume this total number is the same at time t = h.
For the 2D mesh, the connectivities are i, j, n 0 + i j, n 0 + i, n 0 + j.
For the 3D mesh, the connectivities are i, j, k, n Using the previous building of connectivities, obtained by circular permutations, we propose the following generalization of connectivities for the 4D mesh.
i, j, k, l, n is far from being commonplace, except in the case of regular meshes formed by finite elements of multiplexing type (functions of interpolation obtained as the products of functions of space by functions of time).But in the general case, and in particular with the problem of remeshing which is what of interest here, the meshing remains an issue.

4D mesh generation
Figure 3 identifies parts of elementary 2D, 3D and 4D meshes with their node numbering.One denotes by n 0 the total number of nodes at time t = 0 of the entire space mesh, and we assume this total number is the same at time t = h.For the 2D mesh, the connectivities are i, j, n 0 + i j, n 0 + i, n 0 + j.
For the 3D mesh, the connectivities are Using the previous building of connectivities, obtained by circular permutations, we propose the following generalization of connectivities for the 4D mesh.
i, j, k, l, n 0 + i j, k, l, n 0 + i, n 0 + j k, l, n 0 + i, n 0 + j, n 0 + k l, n 0 + i, n 0 + j, n 0 + k, n 0 + l.This 4D mesh is constituted by four hypertetrahedrons (The hypertetrahedron is the four-dimensional tetrahedron.Other names of this element type are simplex or pentatope).We propose to build 4D space-time meshes, resulting from unspecified 3D space meshes, by applying the building technique developed above, for each 3D finite element of the 3D space mesh.However, in this case, it must be checked that the total space-time volume is covered by the 4D mesh.For this purpose, we computed the sum of the volumes of the hypertetrahedrons of the 4D mesh and compared it with the total volume generated by the 3D object multiplied by the time interval.
For a 3D mesh made up of tetrahedrons, the interfaces between the elements are triangles.In 4D, these interfaces are tetrahedrons (see the diagram at the bottom of Figure 4).Therefore, we must thus check that for our technique of mesh generation, all couples of adjoining 4D finite elements have a common tetrahedron.As we use a building technique containing circular permutations, it is necessary to respect a particular order in the numbering of the nodes of each 3D finite element.A way of doing this is to arrange the nodes of each 3D element in ascending order.Table 1 gives an example of a table of connectivities for an elementary 4D mesh, resulting from the 3D space mesh represented by the left-hand diagram of Figure 4. Let us note that this 4D mesh contains eight finite elements, against two for the 3D mesh source, and that n 0 = 5.It is observed that the connectivities are arranged in ascending order.In this case, it is checked that the tetrahedra filling the space-time interface (diagram at the right-hand of Figure 4) are common to the adjoining elements.Indeed, Elements 1 and 5 contain the tetrahedron (1; 2; 4; 6), and Elements 2 and 6 contain the tetrahedron (2; 4; 6; 7).Lastly, Elements 4 and 7 contain the tetrahedron (4; 6; 7; 9).Let us notice that the tetrahedron (1; 2; 3; 4; 5) at time t and the related tetrahedron (6; 7; 8; 9; 10) (not represented) at time t + ∆t could also have different shapes and could be localized at different places (see for example Figure 5).In this case, the triangles (1; 2; 4) and (6; 7; 9) in the right part of Figure 4 can be different.contain the tetrahedron (4; 6; 7; 9).Let us notice that the tetrahedron (1; 2; 3; 4; 5) at time t and the 164 related tetrahedron (6; 7; 8; 9; 10) (not represented) at time t + ∆t could also have different shapes and 165 could be localized at different places (see for exemple in figure 5).In this case, the triangles (1;2;4) and 166 (6;7;9) in the right part of figure 4 can be different. .3D initial mesh (scheme at the left); part of the 4D mesh generated by the triangle (1;2;4) common to the two finite elements of the 3D initial mesh (scheme at the right)

Remeshing technique 168
In this paragraph we present our technique of space-time mesh adaptations.In the literature, 169 many articles on mesh adaptations, [7], [9], [12,13], [23,24], [8], [10] and [22], can be found.Among 170 these papers a large number deals with space-time mesh adaptations.They use the discontinuous 171 Galerkin method.In most of them, the approach is incremental, i.e. remeshing is carried out at given 172 steps of time.Generally the values of the unknown of the new mesh are obtained by approximation or 173 interpolation from those of the old mesh, which we will call "mesh-to-mesh transfer".Moreover, the 174 The hypertetrahedron is the four-dimensional tetrahedron.Other names of this element type are simplex or pentatope.
ain the tetrahedron (4; 6; 7; 9).Let us notice that the tetrahedron (1; 2; 3; 4; 5) at time t and the ed tetrahedron (6; 7; 8; 9; 10) (not represented) at time t + ∆t could also have different shapes and d be localized at different places (see for exemple in figure 5).In this case, the triangles (1;2;4) and ) in the right part of figure 4 can be different.In this paragraph we present our technique of space-time mesh adaptations.In the literature, y articles on mesh adaptations, [7], [9], [12,13], [23,24], [8], [10] and [22], can be found.Among e papers a large number deals with space-time mesh adaptations.They use the discontinuous rkin method.In most of them, the approach is incremental, i.e. remeshing is carried out at given of time.Generally the values of the unknown of the new mesh are obtained by approximation or polation from those of the old mesh, which we will call "mesh-to-mesh transfer".Moreover, the he hypertetrahedron is the four-dimensional tetrahedron.Other names of this element type are simplex or pentatope.

Remeshing Technique
In this section, we present our technique of space-time mesh adaptations.In the literature, many articles on mesh adaptations [2,3,15-21] can be found.Among these papers, a large number deals with space-time mesh adaptations.They use the discontinuous Galerkin method.In most of them, the approach is incremental, i.e., remeshing is carried out at given steps of time.Generally, the values of the unknown of the new mesh are obtained by approximation or interpolation from those of the old mesh, which we will call "mesh-to-mesh transfer".Moreover, the interpolation used is of a multiplexing type; the function of interpolation is defined by the product of a function of space by a function of time.
In [22], we proposed an incremental technique for mesh adaptation, which does not require mesh-to-mesh transfer.It was coupled with problems of rubbing contact (see [14]).In addition, we developed a non-incremental technique of mesh adaptation, based on non-structured space-time meshes.
Some teams have already worked on this problem.One can quote the works of Hugues and Hulbert [2,3], Tezduyar et al. [20,21] and Idesman et al. [9,23].They use the continuous or discontinuous Galerkin method.In these approaches, calculations are carried out on the whole space-time domain Ω × [0, T].Thus, for a field Ω of dimension d and a total number N of nodes of the space-time mesh, the dimension of the linear problem to be solved is d × N, which quickly becomes large when d = 2 or d = 3.A solution to decrease the computational time is to use parallel computations.This is the option chosen in [8,9].
In the context of the continuous Galerkin method, we suggested, in [22,24], a non-incremental solution, which substitutes the concept of a step of time by that of a "space-time front".Erickson et al. [16] have also proposed an advancing-front mesh generation, in the context of the discontinuous Galerkin method.This technique was successfully used by Miller et al. [25] in their multi-field space-time discontinuous Galerkin method, for d = 1 and 2 in linearized elastodynamics applications.The advantage of this frontal resolution is that it decreases the size of the linear systems to be solved.Due to technical difficulties, the frontal resolution in the case of d = 3 has not yet been implemented.Nevertheless, we propose a particular incremental remeshing technique based on the construction of 4D space-time meshes that are able to follow an evolutionary loaded zone.This technique of mesh generation uses simplex finite elements.Figure 5 gives an illustration of the technique.The principle is to maintain the same number of nodes during the simulation, but to locate a sufficiently large number of them under the loaded area.In this particular case, it is possible to preserve the matrices T i/i−1 , T i/i , T i/i+1 identical for all i, which involves a reduction of the computational time.An example of mechanical application is provided in the following section.This technique is aimed at applications in simulating problems of wear between two bodies in contact.
Version May 15, 2018  Preliminary results on the stability of the method have been established in [3] and compared with the Newmark integration scheme.Here we sumarise the main results.Let δ and θ be two real parameters, the Newmark integration scheme reads where ∆t is the time step of integration, { Ui+1 } and { Üi+1 } are respectively the assembled vector of nodal velocities and accelerations at time (i + 1)∆t.We showed in [3] that: • For 1D space-time elastodynamic applications, the use of the STFEM method with linear simplex elements is similar to the use of the implicit Newmark integration scheme with δ = 1/2 and θ = 1/3.The method is then unconditionally stable.
• For 2D space-time elastodynamic applications, the use of the STFEM method with linear simplex elements is similar to the use of the explicit Newmark integration scheme with δ = 1/2 and θ = 0.The method is then conditionally stable.Classically, the time step has to verify the CFL condition : ∆t ≤ min J 2 ω J , where each ω J is the frequency of a normal mode of vibration.
• For higher dimensions (3D and 4D) no direct relationship between the STFEM and the Newmark method has been established.Nevertheless, we noted that our method required a sufficiently small space-time slabs, of the same order of the discretization time step necessary with explicit methods of integration.
• Furthermore, the use of the STFEM method with multiplex elements is similar to the use of the implicit Newmark integration scheme with δ = 1/2 and θ = 1/3, for 1D, 2D, 3D and 4D

Numerical Analysis
Our space-time finite elements method was programmed using MATLAB software and was validated on elementary examples.

Stability
Preliminary results on the stability of the method have been established in [22] and compared with the Newmark integration scheme.Here, we summarize the main results.Let δ and θ be two real parameters; the Newmark integration scheme reads: where ∆t is the time step of integration, { Ui+1 } and { Üi+1 } are respectively the assembled vector of nodal velocities and accelerations at time (i + 1)∆t.We showed in [22] that: • For 1D space-time elastodynamic applications, the use of the STFEM method with linear simplex elements is similar to the use of the implicit Newmark integration scheme with δ = 1/2 and θ = 1/3.The method is then unconditionally stable.

•
For 2D space-time elastodynamic applications, the use of the STFEM method with linear simplex elements is similar to the use of the explicit Newmark integration scheme with δ = 1/2 and θ = 0.
The method is then conditionally stable.Classically, the time step has to verify the CFLcondition: ∆t ≤ min J 2 ω J , where each ω J is the frequency of a normal mode of vibration.

•
For higher dimensions (3D and 4D), no direct relationship between the STFEM and the Newmark method has been established.Nevertheless, we noted that our method required sufficiently small space-time slabs, of the same order of the discretization time step necessary with explicit methods of integration.• Furthermore, the use of the STFEM method with multiplex elements is similar to the use of the implicit Newmark integration scheme with δ = 1/2 and θ = 1/3, for 1D, 2D, 3D and 4D space-time applications.In this case, the method is unconditionally stable.
In the present study, a specific numerical investigation has been carried out to estimate the stability conditions for the STFEM method with linear simplex elements for 4D space-time elastodynamic applications.The stability was tested on a beam of length L = 0.1 m and a square section of 0.01 × 0.01 m 2 (see Figure 6).The Young modulus E was equal to 1000 Pa; the Poisson's ratio ν was equal to 0.3; and the density ρ was equal to 680 kg/m 3 .
Version May 15, 2018 submitted to MDPI 10 of 18 m 2 (see figure 6).The Young modulus E was equal to 1,000 Pa, the Poisson's ratio ν was equal to 0.3 229 and the density ρ was equal to 680 kg/m 3 .

Conditions aux limites imposées
-Sur ⌃ c , on peut imposer des CL de Dirichlet : De plus : ge espace-temps de la figure 2.2, on a  We built boundaries conditions in order to obtain the following analytic solution: where c is the velocity of the wave propagation.On faces (1) and ( 2) null Neumann conditions were imposed.On the other 4 faces (Σ c ) Dirichlet conditions were imposed in order to satisfy the analytic solution (20).The volume external force density f d was assumed to be vanishing.This impose that the wave velocity c must verifies Results plotted in figure 7 give the dependence of the time discretization (size of the space-time 231 slab) ∆t to the average size h of the 3D finite elements, to ensure the stability of the method.

232
Version May 15, 2018 submitted to MDPI 10 of 18 m 2 (see figure 6).The Young modulus E was equal to 1,000 Pa, the Poisson's ratio ν was equal to 0.3 -Sur ⌃ c , on peut imposer des CL de Dirichlet : De plus : llage espace-temps de la figure 2.2, on a  We built boundaries conditions in order to obtain the following analytic solution: where c is the velocity of the wave propagation.On faces (1) and ( 2) null Neumann conditions were imposed.On the other 4 faces (Σ c ) Dirichlet conditions were imposed in order to satisfy the analytic solution (20).The volume external force density f d was assumed to be vanishing.This impose that the wave velocity c must verifies Results plotted in figure 7 give the dependence of the time discretization (size of the space-time 231 slab) ∆t to the average size h of the 3D finite elements, to ensure the stability of the method.

232
Linear fitting suggests that the stability criterion is We built boundary conditions in order to obtain the following analytic solution: where c is the velocity of the wave propagation.On Faces (1) and ( 2), null Neumann conditions were imposed.On the other four faces (Σ c ), Dirichlet conditions were imposed in order to satisfy the analytic solution (20).The volume external force density f d was assumed to be vanishing.This imposes that the wave velocity c must verify: Results plotted in Figure 7 give the dependence of the time discretization (size of the space-time slab) ∆t on the average size h of the 3D finite elements, to ensure the stability of the method.analysis of the STFEM in 4D, we used the previous example.The time step of discretization ∆t is scaled with respect to h, using the stability criterion obtained in the previous paragraph.We computed the maximum error over space at the last time step between the analytic solution and the numerical solution for each mesh size h.The results are plotted in figure 8.
We can note that the convergence rate for the method is nearly quadratic.

Numerical results on mesh adaptation
To illustrate our technique of mesh adaptation we consider the example of a brake disc subjected to the clamping of a plate on one of its faces (see source model with coarse mesh in figure 9) and blocked on the opposite face.The disc is made of steel with a Young's modulus equal to 210,000 MPa, a Poisson's ratio equal to 0.3 and a density equal to 7,800 kg/m 3 .The internal radius of the disc is equal to 40 mm, its external radius is equal to 100 mm and its thickness is 10 mm.Two 3D initial meshes have been tested: a coarse mesh which contains 634 nodes and 1,752 elements, refined only in the area of the clamping, and a fine mesh which contains 6,476 nodes and 25,856 elements and has a uniform mesh fineness over all the sample (see figure 13).
The clamping area is modeled by a constant pressure of 100 MPa.This area is moved along the circumference of the disc with a rotational speed equal to that of the propagating wave where R = 70 mm is the average radius of the disc.For the both 3D initial meshes, we built an Linear fitting suggests that the stability criterion is: ∆t ≤ αh with α 1/2c.Indeed, in our example, c = 1.407 m/s.

Convergence
In the case of simplex finite elements, the convergence with the STFEM method is comparable to the convergence with the Newmark scheme, for the 1D and 2D problem.Concerning the convergence analysis of the STFEM in 4D, we used the previous example.The time step of discretization ∆t is scaled with respect to h, using the stability criterion obtained in the previous paragraph.We computed the maximum error over space at the last time step between the analytic solution and the numerical solution for each mesh size h.The results are plotted in Figure 8.   Numerical values are gathered in figure 11 for the results at time t = 2.10 −5 s and in figure 12 for the results at time t = 4.10 −5 s.Each check point is defined by its angle, in polar coordinates.We can note that the convergence rate for the method is nearly quadratic.

Numerical Results on Mesh Adaptation
To illustrate our technique of mesh adaptation, we consider the example of a brake disc subjected to the clamping of a plate on one of its faces (see the source model with the coarse mesh in Figure 9) and blocked on the opposite face.The disc is made of steel with a Young's modulus equal to 210,000 MPa, a Poisson's ratio equal to 0.3 and a density equal to 7800 kg/m 3 .The internal radius of the disc is equal to 40 mm; its external radius is equal to 100 mm; and its thickness is 10 mm.Two 3D initial meshes have been tested: a coarse mesh, which contains 634 nodes and 1752 elements, refined only in the area of the clamping, and a fine mesh, which contains 6476 nodes and 25,856 elements and has a uniform mesh fineness over all the sample.
The clamping area is modeled by a constant pressure of 100 MPa.This area is moved along the circumference of the disc with a rotational speed equal to that of the propagating wave V = 1 2πR E ρ , where R = 70 mm is the average radius of the disc.For both 3D initial meshes, we built an incremental 4D space-time mesh, which preserves the 3D mesh at each time step, by imposing an axial rotation to keep the finest zone of the 3D mesh under the loading area.An illustration is shown schematically in Figure 5.
The results of the calculations presented were obtained for space-time slabs of 10 −7 s (this is equivalent to using a time step equal to 10 −7 s).The vertical displacements obtained with fine and coarse meshes have been compared for points located on a circle of control, of radius R equal to 70 mm (see Figure 10).Numerical values are gathered in Figure 11 for the results at time t = 2.10 −5 s and in Figure 12 for the results at time t = 4.10 −5 s.Each check point is defined by its angle, in polar coordinates.
Numerical comparisons show that the coarse and fine meshes give similar results in the clamping zone.However, apart from this zone, the results are somewhat different.Let us note that nodal displacements cannot be identical because dynamic effects depend on the fineness of the mesh.Figures 13 and 14 show the norm of incremental displacements at time t = 2.10 −5 s and t = 4.10 −5 s for the coarse and the fine mesh, respectively.For the coarse mesh, we can observe that the refined zone really remains under the zone of clamping.The distribution of the norm of node displacements is similar for the two positions of the load.It is important to note that they are incremental displacements and not total displacements.Finally, it must be noted that the computational time is 6.6-times faster using the coarse mesh.

Conclusions
The method we have presented for space-time mesh generation for 4D domains, using simplex elements, made it possible to develop a technique of mesh adaptation able to follow a mobile loading zone.This original technique has been carried out to ensure a minimal computational time and does

Figure 1 .
Figure 1.2D space-time finite element

Figure 4
Figure 4. 3D initial mesh (scheme at the left); part of the 4D mesh generated by the triangle (1;2;4) common to the two finite elements of the 3D initial mesh (scheme at the right)

Figure 4 .
Figure 4. 3D initial mesh (scheme at the left); part of the 4D mesh generated by the triangle (1;2;4) common to the two finite elements of the 3D initial mesh (scheme at the right)

Figure 4 .
Figure 4. 3D initial mesh (scheme at the left); part of the 4D mesh generated by the triangle (1; 2; 4) common to the two finite elements of the 3D initial mesh (scheme at the right).

submitted to MDPI 9 of 18 tFigure 5 .
Figure 5. Space-time mesh generation by rotation of the loaded area.The loaded area is represented by arrows

Figure 5 .
Figure 5. Space-time mesh generation by rotation of the loaded area.The loaded area is represented by arrows.

Figure 6 .
Figure 6.Geometry and 3D mesh of a beam with h = 6.6 10 −3 m

Figure 6 .
Figure 6.Geometry and 3D mesh of a beam with h = 6.6 10 −3 m

Figure 7 .
Figure 7. Size of the time step of discretization ∆t necessary for stability with respect to the average size h of the 3D finite elements

Figure 7 .
Figure 7. Size of the time step of discretization ∆t necessary for stability with respect to the average size h of the 3D finite elements.

Figure 8 .
Figure 8. Maximum error over the space, at the last time step, between the analytic solution and the numerical solution for each mesh size h, in logarithmic scale

Figure 9 .
Figure 9. 3D initial coarse mesh: Mesh is finer under the clamping zone

Figure 8 .
Figure 8. Maximum error over the space, at the last time step, between the analytic solution and the numerical solution for each mesh size h, in the logarithmic scale.

Figure 9 .
Figure 9. 3D initial coarse mesh: the mesh is finer under the clamping zone.

Figure 10 .
Figure 10.Location zone of checked displacements: circle of radius equal to 70 mm.

Figure 11 .
Figure 11.Comparison of vertical displacements for points situated on the circle of control, expressed in mm, at time t = 2.10 −5 s.Each check point is defined by its angular coordinate, expressed in radians.

Figure 12 .
Figure 12.Comparison of vertical displacements for points situated on the circle of control, expressed in mm, at time t = 4.10 −5 s.Each check point is defined by its angular coordinate, expressed in radians.

Figure 13 .
Figure 13.Isovalues of the norm of nodal displacements, expressed in mm, at t = 2.10 −5 s for the coarse mesh (image at the top) and for the fine mesh (image at the bottom).

Figure 14 .
Figure 14.Isovalues of the norm of nodal displacements, expressed in mm, at t = 4.10 −5 s for the coarse mesh (image on the top) and for the fine mesh (image on the bottom).

Table 1 .
Table of connectivities of the 4D space-time mesh resulting from the elementary 3D space mesh of Figure 4.