# Automatic energy conserving space-time refinement for linear dynamic structural problems

Pauline Cavin, Anthony Gravouil, Antonius Lubrecht, Alain Combescure

► To cite this version:

Pauline Cavin, Anthony Gravouil, Antonius Lubrecht, Alain Combescure. Automatic energy conserving space-time refinement for linear dynamic structural problems. International Journal for Numerical Methods in Engineering, 2005, 64 (3), pp.304-321. 10.1002/nme.1366. hal-00373833

HAL Id: hal-00373833

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

Submitted on 23 Jun 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
Automatic energy conserving space-time refinement for linear
dynamic structural problems

P. Cavin, A. Gravouil, A. A. Lubrecht and A. Combescure*,†

Laboratoire de Mécanique des Contacts et des Solides, UMR CNRS 5514 INSA-LYON, France

In this paper a local space-time automatic refinement method (STAR method) is developed to efficiently solve time-dependent problems using FEM techniques. The automatic process is driven by an energy or a displacement error indicator which controls the precision of the result. The STAR method solves the numerical problem on grids with different mesh size. For the Newmark schemes, a general demonstration, using the energy method, gives the interface conditions between two successive grids which is compatible with the stability of the scheme. Finally, using a linear one-dimensional example, the convergence of the method and the precision of the results are discussed.

KEY WORDS: local mesh refinement; controlled precision; error indicator and energy method

# 1. INTRODUCTION

Ever since computers have entered the world of engineering design, their capacity but also the engineer's expectations have been ever increasing. Today, the most common technique for solving engineering problems is through the use of finite element methods (FEM). As the physical problems treated, become more and more complex, the numerical models become larger and larger. Currently, a large problem can require a million degrees of freedom (dof), and include dynamical and non-linear aspects. The memory and run time requirements of such models are often prohibitive using classical numerical techniques and thus only a small part of the full problem is actually modelled. Hence the prediction of the behaviour of a complete and complex structure can often not be realized.

Furthermore, it is difficult for the engineer performing the calculations to judge the precision of the results obtained. However, precision is an important condition for safe design or for evaluating maximal mechanical variables (as stress, displacement, or strain).

<footnote>*Correspondence to: A. Combescure, INSA de Lyon, LaMCoS, UMR 5514, 18-20 rue des Sciences, F-69621, Villeurbanne, Cedex, France.</footnote>

<footnote>†E-mail: alain.combescure@insa-lyon.fr</footnote>
An efficient mesh structure will assure that the precision in different parts of the total structure is comparable, in order not to waste calculational effort. In order to obtain this constant precision, the mesh size in space and time has to vary over the calculational domain.

Different approaches have been explored to obtain an efficient solution of a dynamic problem. Gravouil et al. [1,2] chose to divide the domain into fixed sub-domains with a given spatial and temporal discretization. The problem on each sub-domain can be solved either explicitly or implicitly. The work shows that imposing velocity continuity at the interface between the sub-domains assures the scheme's stability, but adds virtual damping to the system. Recently, Prakash and Hjelmstad [3] determined the complementary interface conditions which do not modify the numerical properties of the system. In this approach, additional boundary conditions are imposed on the displacement and the acceleration at the interface. The main drawback of this approach is that it requires user experience, as the sub-domains have to be defined before the calculation starts. Wiberg et al. [4], using a discontinuous Galerkin scheme (DGS) coupled to an h-adaptive procedure, propose a method to solve time-dependent problems. The h-refinement method, which refines in space over the entire domain, is controlled by a Zienkewicz [5,6] error estimator. Combe et al. [7] developed another error estimator for transient dynamic simulations. Also p-refinement in time, studied by Campbell et al. [8], can be mentioned.

The presented method is related to subcycling methods [9-13] as different times steps are used over the structure. However, for the STAR method, the domain where the time step varies, changes with time to adhere to the physical phenomena.

In another approach specific space-time finite elements were developed by Hughes and Hulbert [14,15]. The method proposed in this paper is inspired by the multigrid method pioneered by Brandt [16,17]. Another way of adapting multigrid methods to FEM is described by Adams [18,19].

This paper proposes an efficient adaptive FEM to solve, as a first step, a linear dynamic problem. With two initial grids in space and in time over the entire calculational domain, an error indicator generates automatic local space-time grid refinements until the required precision is reached. For this reason the method presented is called STAR (space-time automatic refinement) method. A stability analysis demonstrates which time boundary conditions to apply in case of local refinement. To solve time-dependent problems a Newmark scheme is used.

The FEM package used in this study is Castem 2000 which was developed by Verpeaux et al. [20,21].

## 2. THEORY

### 2.1. Reference problem

In this work only linear structural dynamic problems are studied. To solve these problems classical FEM with a displacement formulation are used. The FEM formulation leads to the following matrix equation (1):

$$M\ddot{u} + Ku = F \qquad (1)$$
where M represents the mass matrix, K the stiffness matrix, F the external force vector, u the unknown displacement vector which depends on the time t, $\ddot{u}$ the unknown acceleration vector which depends on the time t.

M and K are supposed to be symmetric, M is positive definite and K is semi-positive definite.

The problem consists of finding {$\underline{u}, \dot{\underline{u}}, \ddot{\underline{u}}$} satisfying (1) and the initial conditions:

$$ \underline{u}(t=0) = \underline{U}_0 $$

$$ \dot{\underline{u}}(t=0) = \underline{V}_0 $$

where $\underline{U}_0$ and $\underline{V}_0$ are given vectors.

A Newmark algorithm is chosen to solve for $\underline{u}$. The general scheme is defined by two parameters $\beta$ and $\gamma$ and the following equations:

$$ \begin{aligned} M\ddot{\underline{u}}_{n+1} + K\underline{u}_{n+1} &= \underline{F}_{n+1} \\ \underline{u}_{n+1} &= \underline{u}p_n + \beta\Delta t^2\ddot{\underline{u}}_{n+1} \\ \dot{\underline{u}}_{n+1} &= \dot{\underline{u}}p_n + \gamma\Delta t\ddot{\underline{u}}_{n+1} \end{aligned} \quad (2) $$

with the predictors $\underline{u}p_n$ and $\dot{\underline{u}}p_n$:

$$ \begin{aligned} \underline{u}p_n &= \underline{u}_n + \Delta t \dot{\underline{u}}_n + (\frac{1}{2} - \beta) \Delta t^2 \ddot{\underline{u}}_n \\ \dot{\underline{u}}p_n &= \dot{\underline{u}}_n + (1 - \gamma) \Delta t \ddot{\underline{u}}_n \end{aligned} \quad (3) $$

and the initial conditions:

$$ \begin{aligned} \underline{u}_0 &= \underline{U}_0 \\ \dot{\underline{u}}_0 &= \underline{V}_0 \end{aligned} \quad (4) $$

where $\Delta t$ represents the time step, $\underline{F}_n = \underline{F}(t_n)$, $\underline{u}_n$, $\dot{\underline{u}}_n$ and $\ddot{\underline{u}}_n$ are the approximations of $\underline{u}(t_n)$, $\dot{\underline{u}}(t_n)$ and $\ddot{\underline{u}}(t_n)$. For simplicity a state vector $\underline{X} = \{\underline{u}, \dot{\underline{u}}, \ddot{\underline{u}}\}$ is defined. Note that Equations (2) can be rewritten as a function of $\underline{X}$.

## 2.2. Local refinement in space and time

The purpose of this work is to efficiently solve the reference problem with an automatic local refinement technique. The STAR method links the refinement in space and in time. To generate this coupling and to simplify the programming, a space-time element (ste) $\Omega e$ is created. For example on level j, this element $\Omega e^j$ is composed of two finite elements in space $\omega e^j$ and two time steps $t_n^j$ and $t_{n+1}^j$. These elements are illustrated in Figure 1(a). $t_n^j$ represents the time step on level j at the instant n. For the variables, the superscript indicates the level and the subscript the instant.

The horizontal axis $X$ represents space and the vertical axis $T$ represents time. For every space-time element the time is defined by $t_{n+1}^j = t_n^j + \Delta t$. The space element $\omega e_j$ can be either one, two or three-dimensional. In this paper one-dimensional elements are used for simplicity reasons. The space-time element $\Omega e^j$ is defined by a space boundary and a time boundary as:
$\Omega e^j = [\omega e^j] * [n, n+1]$. The space-time domain of level $j$, $\Omega^j$, consists of all the space-time elements on the same level $\Omega e^j$, thus the following relation $\Omega^j = \sum_{\text{ste}} \Omega e^j$ is obtained.

The spatial refinement is the same as presented in a previous paper by the authors [22]. In this study, the refinement is simultaneously applied in the spatial and temporal direction.

Initially, the problem is defined with a first discretization in time and in space which represents the first level on space-time domain $\Omega^1$. Then $\underline{X}^1$ can be calculated on this first space-time domain. To initialize the process a second level is created. Each space-time element is divided in $p$ in each direction. Thus the second space-time domain $\Omega^2$ is defined by $p$ times as many time steps and $p^d$ times as many (in $d$ dimension) space elements for each time step. Then $\underline{X}^2$ can be calculated for level 2. Next the difference between the solutions on level 1 and 2, called the error, is computed. The error, evaluated by an appropriate error indicator, is calculated only at the coinciding time steps. Whenever the error exceeds the precision required by the user, the current space-time element of level 2 is subdivided into $p^{d+1}$ space-time elements, forming the level 3 and the space-time domain $\Omega^3$. Level 3 can consist of the entire space-time domain, or of one (or several) space-time sub-domain(s). Appropriate boundary conditions are added to level 3 equations, such that $\underline{X}^3$ can be calculated. The new error is computed by comparing the solutions on level 2 and 3. Whenever this error exceeds the required precision, the STAR process is repeated creating finer and finer levels.

The solution, the corresponding 'optimal mesh' of the reference problem and the number $k$ of levels required to reach the given precision are known only at the end of the automatic process. The final solution depends on the solution at each level and on the precision required.

For solving the reference problem, note that the level index varies from $j = 1$ to $k$ and the time index varies from $n = 0$ to $r$.

Finally, this process can be summarized by Equations (5)–(7).

$$ \begin{bmatrix} [M^1] \\ [M^2] \\ \vdots \\ [M^k] \end{bmatrix} \left\{ \begin{matrix} \underline{\ddot{u}^1} \\ \underline{\ddot{u}^2} \\ \vdots \\ \underline{\ddot{u}^k} \end{matrix} \right\} + \begin{bmatrix} [K^1] \\ [K^2] \\ \vdots \\ [K^k] \end{bmatrix} \left\{ \begin{matrix} \underline{u^1} \\ \underline{u^2} \\ \vdots \\ \underline{u^k} \end{matrix} \right\} = \left\{ \begin{matrix} \underline{F^1} \\ \underline{F^2} \\ \vdots \\ \underline{F^k} \end{matrix} \right\} \quad (5) $$

<figure><img src="image_2.png" /><figcaption>Figure 1. Space-time element $\Omega e^j$ for a coarse level $j$ and the corresponding space-time elements on the successive level $j+1$. The $i$ index varies from 2 to $p-2$.</figcaption></figure>
Note that $[M^j]$ and $[K^j]$ for $j=1...k$ represent space-time matrices. Hence, they are defined for each time $t_n$ for $n=0...r$. For example the matrix $[K^j]$ is given by

$$[K^j] = \begin{bmatrix} [K_0^j] \\ [K_1^j] \\ \vdots \\ [K_r^j] \end{bmatrix}$$

The final state vector $\underline{X}_{\text{final}}$ is obtained as follow:

$$\underline{X}_{\text{final}} = \underline{X}^k \cup \underline{X}^{k-1} \cup \cdots \cup \underline{X}^{j*} \cup \cdots \cup \underline{X}^{2*} \quad (6)$$

Where the solution $\underline{X}^{j*}$ is the part of $\underline{X}^j$ on the space-time domain $\Omega^j - \Omega^{j+1}$. Please note that $\underline{X}^{k*} = \underline{X}^k$ and the sum ends at 2 as $\Omega^2 - \Omega^1 = \emptyset$.

$$M_{\text{final}} \ddot{\underline{u}}_{\text{final}} + K_{\text{final}} \underline{u}_{\text{final}} = \underline{F}_{\text{final}} \quad (7)$$

The error indicator can be defined using different mechanical variables such as displacement, velocity, acceleration, stress or energy (strain and kinetic) etc. The study of the error indicator in Reference [22] shows that the energy indicator is the most appropriate for solving engineering problems. In this work, only the displacement and the energy indicator are used. The error is defined as a relative indicator. In order to create level $j+2$, the error between the solution on level $j$ and $j+1$ is evaluated. The error is calculated for each space-time element of level $j$, $\Omega e^j$, considering the $p^{d+1}$ corresponding space-time elements of level $j+1$. In order to obtain the error, it is necessary to restrict mechanical quantities from the $p^{d+1}$ space-time elements $\Omega e^{j+1}$ to the corresponding space-time element $\Omega e^j$. Two transfer operators, in space $Ts$ and in time $Tt$, defined in Equation (8), are needed to transfer quantities from level $j+1$ to level $j$ and from the time interval $[pn, pn+p]$ to time interval $[n, n+1]$.

$$\tilde{q}^j = T s_{j+1}^j T t_{[pn, pn+p]}^{[n, n+1]} q^{j+1} \quad (8)$$

The restriction operators used here are the simplest ones. $Ts$ consists of 'keeping' the variables defined on the nodes of the $p^d$ space elements on level $j+1$ only on the nodes of the corresponding space element on level $j$. $Tt$ consists of 'keeping' the variables defined on the time interval $[pn, pn+p]$ on the level $j+1$ only on the time interval $[n, n+1]$ on level $j$.

The error indicator using the displacement is defined as

$$\| \tilde{\underline{u}}^k - \underline{u}^k \| / \max \| \underline{u}^1 \| \quad (9)$$

The energy indicator as:

$$\sqrt{|\tilde{e}^k - e^k|} / \max \sqrt{e_{\text{tot}}^1} \quad (10)$$

With

$$e^k = \frac{1}{2} \dot{\underline{u}}^{kT} M \dot{\underline{u}}^k + \frac{1}{2} \underline{u}^{kT} K \underline{u}^k$$
and

$$ \tilde{e}^k = \frac{1}{2} \tilde{\underline{u}}^{kT} M \tilde{\underline{u}}^k + \frac{1}{2} \tilde{\underline{u}}^{kT} K \tilde{\underline{u}}^k $$

$\max \sqrt{e_{\text{tot}}^1}$ represents the maximum energy of the entire beam on the level 1, and the superscript T indicates the transpose.

For each space-time element on level j the displacement error indicator gives four error values (one for each node of the space-time element $\Omega e^j$) and the energy error indicator gives two error values (one for each space element $\omega e^j$). Only the maximum of the different error values is kept for each space-time element $\Omega e^j$.

This method is efficient because local space-time refinement, using a relative error indicator, occurs only where it is necessary. Furthermore, each space finite element, $\omega e$, which compounds space-time elements, is created only once and can be used many times.

When local refinement is used, the portion of the space-time domain not subdivided has to be taken into account through the boundary conditions. These boundary conditions need to be chosen carefully, to ensure the stability of the algorithm. The boundary condition contains the time boundary conditions as well as the space boundary conditions. For the one-dimensional case the space boundary conditions do not present any difficulty. But for the two or three-dimensional case this aspect has to be detailed. The stability study, presented in the next section, allows one to choose the most appropriate time boundary conditions.

## 2.3. Stability

To study the stability of the STAR, the energy method is used. This method is presented by Hughes and Liu [23, 24], and extended by Gravouil [1] then by Prakash and Hjelmstad [3] to non-matching time interfaces.

As the Newmark algorithm is chosen in this work, the energy method is presented only for this particular case.

Firstly, the energy method is presented in the case without the space-time refinement. The energy method works between two time steps n and n + 1. The following notations are introduced:

$$ \begin{aligned} \langle q \rangle_n^{n+1} &= \frac{1}{2}(q_{n+1} + q_n) \\ [q]_n^{n+1} &= q_{n+1} - q_n \end{aligned} \qquad (11) $$

Using these notations it is possible to rewrite Newmark's equations (2b) and (2c) as (12).

$$ \begin{aligned} [\underline{u}]_n^{n+1} &= \Delta t \langle \dot{\underline{u}} \rangle_n^{n+1} + \frac{\Delta t^2}{2} (2\beta - \gamma) [\ddot{\underline{u}}]_n^{n+1} \\ [\dot{\underline{u}}]_n^{n+1} &= \Delta t \langle \ddot{\underline{u}} \rangle_n^{n+1} + \Delta t (\gamma - \frac{1}{2}) [\ddot{\underline{u}}]_n^{n+1} \end{aligned} \qquad (12) $$

Then the discrete equilibrium equations at the time step n and n + 1 are considered.

$$ \begin{aligned} M \ddot{\underline{u}}_n + K \underline{u}_n &= \underline{F}_n \\ M \ddot{\underline{u}}_{n+1} + K \underline{u}_{n+1} &= \underline{F}_{n+1} \end{aligned} \qquad (13) $$
Premultiplying the two equations (13) by $[\dot{\underline{u}}]_n^{n+1}$, taking the difference between the two new equations and using (12) one obtains:

$$ \frac{1}{2} [\dot{\underline{u}}^T A \underline{u} + \dot{\underline{u}}^T K \dot{\underline{u}}]_n^{n+1} = -\left(\gamma - \frac{1}{2}\right) [\dot{\underline{u}}^T]_n^{n+1} A [\underline{u}]_n^{n+1} + \frac{1}{\Delta t} [\dot{\underline{u}}^T]_n^{n+1} [\underline{F}]_n^{n+1} \quad (14) $$

with

$$ A = M + \frac{\Delta t^2}{2} (2\beta - \gamma)K $$

### Remark

According to Hughes, Liu and Gravouil [1, 23, 24] the external forces do not influence the stability of the algorithm. So Equation (27) can be considered without the terms containing the external forces.

If the following theorems are verified the scheme is stable.

### Theorem 2.1

If $\gamma \ge \frac{1}{2}$ and A is positive definite, $\underline{u}$ and $\dot{\underline{u}}$ are bounded.

### Theorem 2.2

If $K^{-1}$ exists, $\underline{u}$ is also bounded.

In order to satisfy these conditions, the eigenvalue $\omega$ of A needs to verify

$$ 1 + (\beta - \frac{\gamma}{2})(\omega\Delta t)^2 \ge 0 $$

With this method the well known stability conditions of the Newmark algorithm are found.

$$ \begin{gathered} \frac{1}{2} \le \gamma \le 2\beta \quad \text{always stable} \tag{15} \\ \frac{1}{2} \le \gamma \quad \text{and} \quad 2\beta \le \gamma \quad \text{stable if} \quad \Delta t \ge 1/(\omega_{\max}\sqrt{\gamma/2 - \beta}) \end{gathered} $$

When space-time refinement occurs, it is possible to create a finer level only on one part of the initial space-time domain. The influence of the other part of the space-time domain is replaced by an interface force $\underline{F}_i$.

Now it is supposed that one time step on level $j$ is divided in $p$ time steps for level $j+1$. In order to study the stability of the scheme using the STAR method, it is necessary to verify if the energy term $E_{if}^{j+1}$ of the interface force on level $j+1$ between the time step $pn$ and $pn+p$ is the same as the energy term $E_{if}^j$ on level $j$ between the time step $n$ and $n+1$. The two energy terms are defined by Equations (16) and (17). The following demonstration is inspired by Reference [3]:

$$ E_{if}^j = \frac{1}{\Delta t} [\dot{\underline{u}}^T]_n^{n+1} [\underline{F}_i^j]_n^{n+1} \quad (16) $$

$$ E_{if}^{j+1} = \sum_{i=1}^{p} \frac{p}{\Delta t} [\dot{\underline{u}}^{j+1T}]_{pn+i-1}^{pn+i} [\underline{F}_i^{j+1}]_{pn+i-1}^{pn+i} \quad (17) $$

Please note that $\underline{F}_i$ is a non-zero vector only at the interface of a space-time domain.
Using different hypotheses, given in the following derivation, Equation (17) can be simplified.

In order to express $[\dot{\underline{u}}^{j+1T}]_{pn+i-1}^{pn+i}$, a linear interpolation of the velocity field between the time step $pn$ and $pn + p$ is used (see Equation (18)).

$$ \dot{\underline{u}}_{pn+i}^{j+1} = \frac{i}{p} \dot{\underline{u}}_{pn+p}^{j+1} + \left(1 - \frac{i}{p}\right) \dot{\underline{u}}_{pn}^{j+1} \quad \text{for } i = 0, 1, \dots, p \qquad (18) $$

Using (18) the term $[\dot{\underline{u}}^{j+1}]_{pn+i-1}^{pn+i}$ can be written as below:

$$ \begin{aligned} [\dot{\underline{u}}^{j+1}]_{pn+i-1}^{pn+i} &= (\dot{\underline{u}}_{pn+i}^{j+1} - \dot{\underline{u}}_{pn+i-1}^{j+1}) \\ &= (i/p)\dot{\underline{u}}_{pn+p}^{j+1} + (1 - (i/p))\dot{\underline{u}}_{pn}^{j+1} + ((i-1)/p)\dot{\underline{u}}_{pn+p}^{j+1} + (1 - ((i-1)/p))\dot{\underline{u}}_{pn}^{j+1} \\ &= (1/p)(\dot{\underline{u}}_{pn+p}^{j+1} - \dot{\underline{u}}_{pn}^{j+1}) \end{aligned} \qquad (19) $$

In order to develop the second term $[\underline{\underline{F}}_{i}^{j+1}]_{pn+i-1}^{pn+i}$ of (17), equilibrium equation (20) on level $j+1$ at time steps $pn + i$ for $i = 0, 1, \dots, p$ is used.

$$ M^{j+1} \ddot{\underline{u}}_{pn+i}^{j+1} + K^{j+1} \underline{u}_{pn+i}^{j+1} = \underline{\underline{F}}_{i}^{j+1} \quad \text{for } i = 0, 1, \dots, p \qquad (20) $$

Please note that $M^{j+1}$ and $K^{j+1}$ are constant on the time interval $[pn, pn + p]$. Substituting (20) in $[\underline{\underline{F}}_{i}^{j+1}]_{pn+i-1}^{pn+i}$ one obtains

$$ \begin{aligned} [\underline{\underline{F}}_{i}^{j+1}]_{pn+i-1}^{pn+i} &= \underline{\underline{F}}_{i}^{j+1}_{pn+i} - \underline{\underline{F}}_{i}^{j+1}_{pn+i-1} \\ &= M^{j+1}(\ddot{\underline{u}}_{pn+i}^{j+1} - \ddot{\underline{u}}_{pn+i-1}^{j+1}) + K^{j+1}(\underline{u}_{pn+i}^{j+1} - \underline{u}_{pn+i-1}^{j+1}) \end{aligned} \qquad (21) $$

Linear interpolation of the displacement and the acceleration field between the time step $pn$ and $pn + p$ are now introduced (Equation (22)).

$$ \underline{u}_{pn+i}^{j+1} = \frac{i}{p} \underline{u}_{pn+p}^{j+1} + \left(1 - \frac{i}{p}\right) \underline{u}_{pn}^{j+1} \quad \text{for } i = 0, 1, \dots, p \qquad (22) $$

$$ \ddot{\underline{u}}_{pn+i}^{j+1} = \frac{i}{p} \ddot{\underline{u}}_{pn+p}^{j+1} + \left(1 - \frac{i}{p}\right) \ddot{\underline{u}}_{pn}^{j+1} \quad \text{for } i = 0, 1, \dots, p $$

Substituting (22) in (21) Equation (23) is obtained.

$$ [\underline{\underline{F}}_{i}^{j+1}]_{pn+i-1}^{pn+i} = (1/p)(\underline{\underline{F}}_{i}^{j+1}_{pn+p} - \underline{\underline{F}}_{i}^{j+1}_{pn}) \qquad (23) $$
Then (19) and (23) are introduced in (17).

$$
\begin{align}
E_{if}^{j+1} &= (p/\Delta t) \sum_{i=1}^{p} [\underline{\dot{u}}^{j+1T}]_{pn+i-1}^{pn+i} [\underline{Fi}^{j+1}]_{pn+i-1}^{pn+i} \nonumber \\
&= (1/(p\Delta t)) \left\{ \sum_{i=1}^{p} (\underline{\dot{u}}_{pn+p}^{j+1T} - \underline{\dot{u}}_{pn}^{j+1T}) (\underline{Fi}_{pn+p}^{j+1} - \underline{Fi}_{pn}^{j+1}) \right\} \nonumber \\
&= (1/\Delta t)(\underline{\dot{u}}_{pn+p}^{j+1T} - \underline{\dot{u}}_{pn}^{j+1T})(\underline{Fi}_{pn+p}^{j+1} - \underline{Fi}_{pn}^{j+1}) \nonumber \\
&= (1/\Delta t)[\underline{\dot{u}}^{j+1T}]_{pn}^{pn+p} [\underline{Fi}^{j+1}]_{pn}^{pn+p} \tag{24}
\end{align}
$$

At the coinciding time step on level j and j + 1, the state vector is projected in space using
a Ps operator. So Equation (25) gives the relation between X^j and X^{j+1}.

$$
\begin{align}
\underline{X}_n^j &= P_s \quad \underline{X}_{pn}^{j+1} \notag \\
\underline{X}_{n+1}^j &= P_s \quad \underline{X}_{pn+p}^{j+1} \tag{25}
\end{align}
$$

With an appropriate projection operator *Ps*, see the following remark, relation (26) is verified.

$$
\begin{align}
\underline{Fi}_n^j &= P_s \mathrm{T} \quad \underline{Fi}_{pn}^{j+1} \tag{26} \\
\underline{Fi}_{n+1}^j &= P_s \mathrm{T} \quad \underline{Fi}_{pn+p}^{j+1}
\end{align}
$$

Remark

The projection operator *Ps* used in the example has to be compatible with the FE assembly
matrix. Thus the relations *M*<sub>*n*</sub><sup>*j*</sup> = *Ps*<sup>*T*</sup>*M*<sub>*pn*</sub><sup>*j*+1</sup>*Ps*, *M*<sub>*n*+1</sub><sup>*j*</sup> = *Ps*<sup>*T*</sup>*M*<sub>*pn*+*p*</sub><sup>*j*+1</sup>*Ps*, *K*<sub>*n*</sub><sup>*j*</sup> = *Ps*<sup>*T*</sup>*K*<sub>*pn*</sub><sup>*j*+1</sup>*Ps*
and *K*<sub>*n*+1</sub><sup>*j*</sup> = *Ps*<sup>*T*</sup>*K*<sub>*pn*+*p*</sub><sup>*j*+1</sup>*Ps* have to be verified, see Appendix A.

Finally, the energy term of the interface force can be written as (27)

$$
\begin{align}
E_{if}^{j+1} &= (1/\Delta t)[\underline{\dot{u}}^{j+1T}]_{pn}^{pn+p} [\underline{Fi}^{j+1}]_{pn}^{pn+p} \nonumber \\
&= (1/\Delta t)[\underline{\dot{u}}^j]^n_n P_s^T [\underline{Fi}^{j+1}]_{pn}^{pn+p} \nonumber \\
&= (1/\Delta t)[\underline{\dot{u}}^j]^n_n [\underline{Fi}^j]_n^{n+1} \nonumber \\
&= E_{if}^j \tag{27}
\end{align}
$$

So, this section allows to conclude that using the projection of X̄ for the time boundary
conditions on the coinciding time steps from level j to level j + 1 and a linear interpolation
for the others time steps, no energy is injected at the interface between the levels j + 1 and j
in the Newmark scheme. So using these time boundary conditions the STAR method is stable,
and the displacement, velocity and acceleration are continuous at the interface between two
levels.
### Remark 1

This development can be extended if there is a constant damping matrix C.

### Remark 2

For the particular case $p = 2$, which means that each space-time element is subdivided in four (in one dimension, Figure 2), only the linear time interpolation of the velocity is necessary for the non-coinciding time steps.

$$ 
\begin{aligned}
E_{if}^{j+1} &= \sum_{i=1}^{2} \frac{2}{\Delta t} [\dot{\underline{u}}_{2n+i}^{j+1T}]_{2n+i-1}^{2n+i} [\underline{Fi}_{2n+i}^{j+1}]_{2n+i-1}^{2n+i} \\
&= \frac{2}{\Delta t} [(\dot{\underline{u}}_{2n+2}^{j+1} - \dot{\underline{u}}_{2n+1}^{j+1})^T (\underline{Fi}_{2n+2}^{j+1} - \underline{Fi}_{2n+1}^{j+1}) + (\dot{\underline{u}}_{2n+1}^{j+1} - \dot{\underline{u}}_{2n}^{j+1})^T (\underline{Fi}_{2n+1}^{j+1} - \underline{Fi}_{2n}^{j+1})]
\end{aligned} 
\quad (28) $$

Now, only the interpolation of the velocity is used for $\dot{\underline{u}}_{2n+1}^{j+1}$.

$$ \dot{\underline{u}}_{2n+1}^{j+1} = \frac{1}{2} (\dot{\underline{u}}_{2n+2}^{j+1} + \dot{\underline{u}}_{2n}^{j+1}) $$

Developing (28) and using the linear interpolation of the velocity it appears that no further hypotheses are needed because all the terms with $\underline{Fi}_{2n+1}^{j+1}$ vanish. Finally, one obtains:

$$ 
\begin{aligned}
E_{if}^{j+1} &= \frac{1}{\Delta t} [\dot{\underline{u}}_{2n+2}^{j+1T} \underline{Fi}_{2n+2}^{j+1} - \dot{\underline{u}}_{2n}^{j+1T} \underline{Fi}_{2n+2}^{j+1} - \dot{\underline{u}}_{2n+2}^{j+1T} \underline{Fi}_{2n}^{j+1} + \dot{\underline{u}}_{2n}^{j+1T} \underline{Fi}_{2n}^{j+1}] \\
&= \frac{1}{\Delta t} (\dot{\underline{u}}_{2n+2}^{j+1T} - \dot{\underline{u}}_{2n}^{j+1T}) (\underline{Fi}_{2n+2}^{j+1} - \underline{Fi}_{2n}^{j+1})
\end{aligned} 
\quad (29) $$

As previously at the coinciding time step the projection of $\dot{\underline{u}}$ from level j to level j + 1 and relation (26) lead to Equation (30).

$$ E_{if}^{j+1} = \frac{1}{\Delta t} (\dot{\underline{u}}_{n+1}^{jT} - \dot{\underline{u}}_{n}^{jT}) (\underline{Fi}_{n+1}^{j} - \underline{Fi}_{n}^{j}) \quad (30) $$

<figure><img src="image_3.png" /><figcaption>Figure 2. Space-time element $\omega e^j$ for a coarse level j and the four corresponding space-time elements on the successive level j + 1 in the particular case $p = 2$.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Figure 3. Example of a two level (j and j + 1) space-time refinement with the boundary conditions on level j + 1.</figcaption></figure>

## 2.4. Space-time boundary conditions

In the rest of this paper only the particular case where $p = 2$ is treated. So this section presents the space-time boundary conditions for this particular case. Figure 2 illustrated the space-time elements for $p = 2$.

Using the result of the stability analysis, the strategy chosen for the space-time boundary conditions is depicted in Figure 3. In this figure a specific space-time refinement illustrates all the possible space and time boundary conditions. Four different types can be identified:

* the circles represent the nodes where the projection of $\dot{u}$ from level j to level j + 1 is used (coinciding nodes only).

* the crosses represent the nodes where the linear interpolation in time of the velocity is used (nodes without coarse grid information).

* the diamonds represent the nodes where the linear projection in space of $u$ and $\dot{u}$ from level j to level j + 1 is used (nodes where the $Ps$ operator is used). The $\ddot{u}$ unknown are calculated with the equilibrium equation at this instant.

* the squares represent the nodes where $X$ is calculated using the Newmark scheme.

# 3. RESULTS

## 3.1. Reference case

The example presented is a one-dimensional beam. One extremity is embedded and a traction force is applied on the other end. The external force is a 'smooth shock', see Figures 4 and 5. Classical steel parameters are used. The characteristic values are a Young modulus of $E = 2 \times 10^{11}$ Pa, a Poisson ratio of $v = 0.3$ and a specific mass of $\rho = 7800$ Kg m$^{-3}$. The length of the beam is $L = 2$ m, its section is $S = 5 \times 10^{-3}$ m$^2$.
<figure><img src="image_5.png" /><figcaption>Figure 4. One-dimensional beam in traction.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 5. ‘Smooth shock’ load function $F(t)$.</figcaption></figure>

The Newmark characteristic values, used in this problem, are $\beta = 1/4$ and $\gamma = 1/2$, which correspond to the mean acceleration scheme. This implicit scheme is of second order precision. The initial mesh supplied, is uniform and consists of 128 space-time elements. This corresponds to 8 space finite elements for each time step and 32 time steps. The initial time step is equal to $L/(8 * \sqrt{E/\rho})$ which corresponds to the critical time step. The space finite elements are classical rod elements with 1 dof at each node.

## 3.2. Refinement

This section presents the meshes obtained at the end of the automatic process for a displacement and an energy error indicator. Figure 6 illustrates the final mesh at different time steps using the energy error indicator and a relative precision of $10^{-4}$. In this figure the left extremity of the beam is embedded and the traction force is applied on the right extremity. To facilitate viewing only the nodes of the meshes are depicted. This figure shows that the spatial mesh changes with time. Mesh refinements move with the waves propagating through the beam. The mesh captures both wave propagation and reflection.

It is possible to present the final mesh in a more condensed way, see Figure 7. The horizontal axis represents the space discretization and the vertical the time discretization. Each horizontal line corresponds to the spatial mesh for one time step. In this representation all space-time elements, introduced in the previous part, appear. The remarks made for the previous figure are more clearly visible in this configuration.

Before analysing the velocity or the displacement solution, the final space-time mesh allows one to identify the wave propagation in the structure.
<figure><img src="image_7.png" /><figcaption>Figure 6. Final mesh obtained with a relative energy error indicator for $t = 0$ s to $t = 1.54 \times 10^{-3}$ s with a time step of $\Delta t = 9.64 \times 10^{-4}$ s.</figcaption></figure>

Figure 7 highlights the originality of the method. Using the STAR technique on the space-time elements, spatial coarsening occurs without the use of a specific coarsening technique. Physical phenomena can be followed by the evolution of the space-time refinement.

This approach controls the precision of the calculated solution. This topic is treated in the next section. Before discussing precision, the numerical energy conservation was verified for two specific cases. The order of magnitude of the different energy terms is $10^6$ J, the energy balance is checked up to $10^{-9}$ J so computer precision is reached. These results indicate that the STAR method is stable.
<figure><img src="image_8.png" /><figcaption>Figure 7. Final space-time mesh using a displacement error indicator (left, required precision $5 \times 10^{-3}$) and an energy error indicator (right, required precision $10^{-4}$).</figcaption></figure>

### 3.3. Precision

The automatic process is driven by the precision required by the user. Finer levels are created whenever the error exceeds the given precision. Figure 8 illustrates the maximum relative error as a function of the level, using the displacement and the energy error indicator. These curves are obtained, respectively, for a given relative precision of $5 \times 10^{-3}$ and $10^{-4}$ corresponding to Figure 7. In this figure the incremental effect of the Newmark scheme appears. As time progresses finer levels are needed to reach the requested precision. Furthermore, the size of a specific space-time domain increases with time.

In Figure 8 the reference curves are obtained by comparing the solutions on two grids covering the entire space-time domain. The convergence rates are different for the two error indicators. The displacement error indicator seems to converge with an order between 1 and 2, whereas the energy converges with order 0.5. Similar results have been observed and explained in Reference [22]. For these two cases, 6 levels are necessary before the requested precision is reached. The requested precision for the displacement and the energy error indicator are different ($5 \times 10^{-3}$ and $10^{-4}$) because the order of magnitude of the numerator and the denominator for the two error indicators (Equations (9) and (10)) are different. This difference explains the gap between the two curves on Figure 8.
<figure><img src="image_9.png" /><figcaption>Figure 8. Convergence using the displacement and the energy indicator.</figcaption></figure>

<table><caption>Table I. Correlation between the requested precision and the precision of the solution with the displacement error indicator.</caption><thead><tr><th>Requested precision</th><th>Level max reached</th><th>Max Δu last level</th><th>Max Δu all levels</th></tr></thead><tbody><tr><td>5 × 10<sup>-2</sup></td><td>4</td><td>4.3 × 10<sup>-3</sup></td><td>1.7 × 10<sup>-2</sup></td></tr><tr><td>2 × 10<sup>-2</sup></td><td>5</td><td>1.7 × 10<sup>-3</sup></td><td>1.7 × 10<sup>-2</sup></td></tr><tr><td>10<sup>-2</sup></td><td>5</td><td>1.8 × 10<sup>-3</sup></td><td>6.9 × 10<sup>-3</sup></td></tr></tbody></table>

Table I presents the relation between the required precision and the precision obtained in the solution using the displacement error indicator. In this table the maximum level attained is given. The two other columns give the maximal difference between the displacement field obtained with the local STAR and a reference solution obtained with a uniform space-time mesh corresponding to the maximum level reached. The first case compares only the displacement field obtained on the last level with the reference field. The second case compares the final displacement field obtained with the reference field. The final displacement is defined on different levels, and thus on different mesh sizes. In order to perform the comparison, the final displacement field has to be 'projected' to the finest level. The maximum displacement value corresponding to a coarse grid is interpolated to the equivalent coinciding finest grid. This table shows that the precision obtained in the local refinement solution satisfies the requested precision. Hence, the accuracy of the solution is controlled using the space-time refinement strategy.

## 4. CONCLUSIONS

This paper presents a novel strategy to efficiently solve time-dependent FEM problems and to control the precision of the solution. Specific space-time elements are created in order to link space and time refinement and to simplify the programming. The space time refinement is
driven by an automatic process based on a relative error indicator. A stability analysis using the energy method defines the correct time boundary conditions to apply. It was shown theoretically and numerically that the scheme conserves energy up to machine precision.

Furthermore, performing STAR, spatial ‘coarsening’ occurs automatically. One observes that the spatial mesh moves with the wave propagating in the structure: the numerical solution process follows the physical phenomena occurring. The convergence rate of the method is studied using the displacement and the energy error indicator. Finally, the required precision of the solution has been compared to show quantitatively the precision of the solution obtained. This strategy is illustrated using a one-dimensional, linear example and an implicit Newmark scheme. Future work should extend this strategy to problems in higher dimensions, and including non-linear aspects. For non-linear problems it might be necessary to cycle between the different levels in order to determine the correct boundary of the domain for each level. Furthermore, explicit solution techniques should be explored.

## APPENDIX A: COMPATIBILITY BETWEEN THE PROJECTION MATRIX $Ps$ AND THE FE ASSEMBLY MATRIX

When local space-time refinement occurs, for one-dimensional problem, a rod element from level $j$ is divided into two space elements belonging to the level $j+1$. The following notations are used: $X^j$ represents the state vector on level $j$ and $X^{j+1}$ represents the state vector on level $j+1$. In the configuration of Figure A1, the state vectors are: $X^j = \{X1^j, X2^j\}$ and $X^{j+1} = \{X1^{j+1}, X2^{j+1}, X3^{j+1}\}$.

In the STAR method, the quantities defined on the level $j$ are projected in space on level $j+1$ with a $Ps$ projection operator. This operator is linear and is defined as below:

$$Ps = \begin{bmatrix} 1 & 0 \\ 0.5 & 0.5 \\ 0 & 1 \end{bmatrix} \qquad (A1)$$

Thus the relations between $X^{j+1}$ and $X^j$ are given in Equations (A2) and (A3).

$$X^{j+1} = PsX^j \qquad (A2)$$

<figure><img src="image_10.png" /><figcaption>Figure A1. One rod element (a); and the two corresponding rod elements (b) for $p = 2$.</figcaption></figure>
$$ \begin{bmatrix} X1^{j+1} \\ X2^{j+1} \\ X3^{j+1} \end{bmatrix} = \begin{bmatrix} 1 & 0 \\ 0.5 & 0.5 \\ 0 & 1 \end{bmatrix} \begin{bmatrix} X1^j \\ X2^j \end{bmatrix} \quad (A3) $$

For configuration (a) of Figure A1 the mass matrix $m^j$ and the stiffness matrix $k^j$ are defined by Equation (A4).

$$ m^j = \frac{\rho Sl}{6} \begin{bmatrix} 2 & 1 \\ 1 & 2 \end{bmatrix} \quad k^j = \frac{ES}{l} \begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix} \quad (A4) $$

For configuration (b) of Figure A1 the mass matrix $m^{j+1}$ and the stiffness matrix $k^{j+1}$ are defined by Equation (A5).

$$ m^{j+1} = \frac{\rho Sl}{12} \begin{bmatrix} 2 & 1 & 0 \\ 1 & 4 & 1 \\ 0 & 1 & 2 \end{bmatrix} \quad k^{j+1} = \frac{2ES}{l} \begin{bmatrix} 1 & -1 & 0 \\ -1 & 2 & -1 \\ 0 & -1 & 1 \end{bmatrix} \quad (A5) $$

In order to verify the compatibility requirement between the $Ps$ operator and the FE assembly matrix, the relations $Ps^T m^{j+1} Ps$ and $Ps^T k^{j+1} Ps$ have to be equal, respectively, to $m^j$ and $k^j$.

$$ Ps^T m^{j+1} Ps = \begin{bmatrix} 1 & 0.5 & 0 \\ 0 & 0.5 & 1 \end{bmatrix} \frac{\rho Sl}{12} \begin{bmatrix} 2 & 1 & 0 \\ 1 & 4 & 1 \\ 0 & 1 & 2 \end{bmatrix} \begin{bmatrix} 1 & 0 \\ 0.5 & 0.5 \\ 0 & 1 \end{bmatrix} = \frac{\rho Sl}{6} \begin{bmatrix} 2 & 1 \\ 1 & 2 \end{bmatrix} = m^j \quad (A6) $$

$$ Ps^T k^{j+1} Ps = \begin{bmatrix} 1 & 0.5 & 0 \\ 0 & 0.5 & 1 \end{bmatrix}^T \frac{2ES}{l} \begin{bmatrix} 1 & -1 & 0 \\ -1 & 2 & -1 \\ 0 & -1 & 1 \end{bmatrix} \begin{bmatrix} 1 & 0 \\ 0.5 & 0.5 \\ 0 & 1 \end{bmatrix} = \frac{ES}{l} \begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix} = k^j \quad (A7) $$

Equations (A6) and (A7) prove that the chosen projection operator $Ps$ is compatible with the FE assembly matrix. So this operator satisfies the stability requirement.

REFERENCES

1. Gravouil A. Méthode multi-échelles en temps et en espace avec décomposition de domaines pour la dynamique non-linéaire des structures. Ph.D Thesis, Ecole Normale Supérieure de Cachan, 2000.

2. Gravouil A, Combescure A. Multi-time-step and two-scale domain decomposition method for non-linear structural dynamics. International Journal for Numerical Methods in Engineering 2003; 58:1545–1569.

3. Prakash A, Hjelmstad KD. A feti based multi-step coupling method for newmark schemes in structural dynamics. International Journal for Numerical Methods in Engineering 2004; 61:2183–2204.

4. Li XD, Wiberg N-E. Implementation and adaptivity of a space-time finite element method for structural dynamics. Computer Methods in Applied Mechanics and Engineering 1998; 156:211–229.
5. Zienkiewicz OC, Zhu JZ. A simple error estimator and adaptative procedure for practical engineering analysis. *International Journal for Numerical Methods in Engineering* 1987; **24**:337–357.

6. Zienkiewicz OC, Taylor RL. *The Finite Element Method*, vol. 1 (4th edn). McGraw-Hill: New York, 1989.

7. Combe JP, Ladevéze P, Pelle JP. Discretization error estimator for transient dynamic simulations. *Advances in Engineering Software* 2002; **33**:553–563.

8. Campbell JS, Robinson JC. *p*-adaptativity in the time domain using mvpq multivalue time-marching integrators. *International Journal for Numerical Methods in Engineering* 1990; **40**:2343–2367.

9. Daniel WJT. A study of the stability of subcycling algorithms in structural dynamics. *Computer Methods in Applied Mechanics and Engineering* 1998; **156**:1–13.

10. Daniel WJT. A partial velocity approach to subcycling structural dynamics. *Computer Methods in Applied Mechanics and Engineering* 2003; **192**:375–394.

11. Smolinski P, Belytschko T, Neal M. Multi-time-step integration using nodal partitioning. *International Journal for Numerical Methods in Engineering* 1988; **26**:349–359.

12. Smolinski P. Subcycling integration with non-integer time steps for structural dynamics problems. *Computers and Structures* 1995; **59**:236–244.

13. Wu YS, Smolinski P. A multi-time step integration algorithm for structural dynamics based on the modified trapezoidal rule. *Computer Methods in Applied Mechanics and Engineering* 2000; **187**:641–660.

14. Hughes TJR, Hulbert GM. Space-time finite element methods for elastodynamics: formulations and error estimates. *Computer Methods in Applied Mechanics and Engineering* 1988; **67**:339–363.

15. Hulbert GM, Hughes TJR. Space-time finite element methods for second-order hyperbolic equations. *Computer Methods in Applied Mechanics and Engineering* 1990; **84**:327–348.

16. Brandt A. Multi-level adaptative technique (mlat) for fast numerical solution to boundary value problems. *Lecture Notes in Physics* 1977; **18**:82–89.

17. Brandt A. *Multigrid Techniques: 1984 Guide with Applications to Fluid Dynamics*. GMD-Studien Editions, 1984.

18. Adams M, Taylor RL. Parallel multigrid solvers for 3d-unstructured large deformation elasticity and plasticity finite element problems. *Finite Element in Analysis and Design* 2000; **36**:197–214.

19. Adams M. Evaluation of three unstructured multigrid methods on 3d finite element problems in solids mechanics. *International Journal for Numerical Methods in Engineering* 2002; **55**:519–534.

20. Verpeaux P, Charras T, Millard A. Castem 2000 une approche moderne du calcul des structures. In *Calcul des Structures et Intelligence Artificielle Pluralis*, Fouet JM, Ladevéze P, Ohayon R (eds). 1988; 261–271.

21. Verpeaux P, Millard A. Castem 2000 guide du développement. *Commissariat à l'Energie Atomique Rapport DMT*, 1991.

22. Cavin P, Gravouil A, Lubrecht AA, Combescure A. Efficient fem calculation with predefined precision through automatic grid refinement. *Finite Element in Analysis and Design* 2005; **41**:1043–1055.

23. Hughes TJR, Liu WK. Implicit-explicit finite element in transient analysis: stability theory. *Journal of Applied Mechanics* 1978; **45**:371–374.

24. Hughes TJR, Liu WK. Implicit-explicit finite elements in transient analysis: implementation and numerical example. *Journal of Applied Mechanics* 1978; **45**:375–378.