# Automatic time stepping algorithms for implicit numerical simulations of blade/casing interactions.

L Noels*, L Stainier*, J-P Ponthot* and J Bonini**

*LTAS - Milieux Continus et Thermomécanique, Université de Liège, Chemin des Chevreuils 1, 4000 Liège, Belgium

**SNECMA Moteurs - Engineering Division, Centre de Villaroche, 77550 Moissy-Cramayel, France

**Abstract** - An automatic time stepping algorithm for non-linear problems, solved by implicit schemes, is presented. The time step computation is based on the estimation of an integration error calculated from the acceleration difference. It is normalised according to the size of the problem and the integration parameters. This time step control algorithm modifies the time step size only if the problem has a long time physical change. Additionally, the Hessian matrix can be kept constant for several iterations, even though the problem is non-linear. A criterion selecting if the Hessian matrix must be calculated or not is developed. Finally, a criterion of iterations divergence is also proposed. It avoids the determination, by the user, of a maximal iteration number. This minimises the total number of iterations, and thus the computation cost. Industrial numerical examples are presented that demonstrate the performances (precision and computational cost) of the algorithms.

## NOTATION

<table><tr><td>x</td><td>vector of nodal positions</td></tr><tr><td>x&#x0307;</td><td>vector of nodal velocities</td></tr><tr><td>x&#x0308;</td><td>vector of nodal accelerations</td></tr><tr><td>M</td><td>mass matrix</td></tr><tr><td>F<sup>int</sup></td><td>vector of internal forces</td></tr><tr><td>F<sup>ext</sup></td><td>vector of external forces</td></tr><tr><td>R</td><td>residual (out-of-balance) vector</td></tr><tr><td>&sigma;</td><td>Cauchy stress tensor</td></tr><tr><td>b</td><td>body force per unit mass</td></tr><tr><td>f</td><td>surface forces</td></tr><tr><td>&rho;</td><td>current mass density of the material</td></tr><tr><td>&rho;<sub>0</sub></td><td>initial mass density of the material</td></tr><tr><td>N</td><td>matrix of FEM shape functions</td></tr><tr><td>B</td><td>matrix of the derivatives of the FEM shape functions</td></tr><tr><td>&beta;</td><td>first Newmark parameter</td></tr><tr><td>&gamma;</td><td>second Newmark parameter</td></tr><tr><td>&alpha;<sub>M</sub></td><td>first free parameter balancing sampling time around [t<sub>n</sub>, t<sub>n+1</sub>] for averaging inertia terms</td></tr><tr><td>&alpha;<sub>f</sub></td><td>second free parameter balancing sampling time around [t<sub>n</sub>, t<sub>n+1</sub>] for averaging forces</td></tr><tr><td>&delta;</td><td>tolerance on the non-dimensional residual</td></tr><tr><td>K<sub>f</sub></td><td>tangent stiffness matrix</td></tr><tr><td>C<sub>T</sub></td><td>tangent damping matrix</td></tr><tr><td>S</td><td>Hessian matrix</td></tr><tr><td>S(t)</td><td>current surface of the body</td></tr><tr><td>V(t)</td><td>current volume of the body</td></tr><tr><td>S<sub>0</sub></td><td>initial surface of the body</td></tr><tr><td>V<sub>0</sub></td><td>initial volume of the body</td></tr><tr><td>&Delta;t</td><td>time increment: &Delta;t = t<sub>n+1</sub> - t<sub>n</sub></td></tr><tr><td>t</td><td>time</td></tr><tr><td>r</td><td>non-dimensional residual of the equations of motion</td></tr><tr><td>&omega;</td><td>pulsation for a one-degree-of-freedom linear system</td></tr><tr><td>&Omega;</td><td>non-dimensional pulsation (&omega;&Delta;t)</td></tr><tr><td>&epsilon;</td><td>integration error for a one-degree-of-freedom linear system</td></tr><tr><td>e</td><td>integration error for the system</td></tr><tr><td>&eta;</td><td>parameter to modify time step size</td></tr><tr><td>CO</td><td>counter limit before dividing time step size</td></tr><tr><td>CT</td><td>counter limit before augmenting time step size</td></tr><tr><td>ERRO</td><td>maximal integration error of last time step before dividing time step size</td></tr><tr><td>ERRT</td><td>maximal integration error of last time step before augmenting time step size</td></tr><tr><td>PRCU</td><td>tolerance on the integration error</td></tr><tr><td>RAPRES</td><td>ratio of two successive non-dimensional residual</td></tr><tr><td>RAT</td><td>factor multiplying the time step size between two steps</td></tr><tr><td>REJL</td><td>maximal integration error accepted</td></tr></table>
# INTRODUCTION

Non-linear dynamics problems can be solved with two kind of time stepping algorithms: explicit or implicit. For an explicit algorithm, the solution at time $t_{n+1}$ depends only on the solution at time $t_n$, while for an implicit algorithm, it also implicitly depends on the solution at time $t_{n-1}$ itself. The problem must then be solved iteratively. Stability (i.e. positive damping of initial perturbations) imposes different restrictions on those two families of algorithms and, with a proper choice of parameters, the time step size can be much larger for an implicit algorithm than for an explicit algorithm. The total number of time steps will thus generally be smaller, and even though the cost of a time step is higher as a consequence of the need for computing and inverting a Hessian matrix, the total computation time for an implicit scheme is often more interesting than for an explicit scheme. In this context, if the time step size is chosen too small, the calculation cost is very expensive, while if it is chosen too large, the integration is not accurate enough or the iterations diverge (when solving the balance equations). Therefore, time step size should be carefully evaluated. Since the problem evolves with time, time step size should be continuously adapted to this evolution. An automatic time stepping algorithm is then the only solution to accurately solve the problem in a reasonably short computation time.

For an industrial problem that has a large number of degrees of freedom, the most expensive operation of an implicit code is the inversion of the Hessian Matrix. For non-linear problems, the Hessian matrix normally changes every iteration, but the Newton-Raphson iterations can sometimes converge while using the old inverted matrix. Still, this inverted matrix must be regularly recomputed to avoid divergence. In a classical strategy, this inversion occurs at the beginning of each time step and for some iterations selected by the user. But if the Hessian matrix is not inverted for too many iterations, the problem diverges, while if the inversion occurs too frequently, the problem becomes too expensive. According to the evolution of the problem with time, an algorithm automatically selecting if the inverted Hessian matrix must be recalculated or not can reduce the total computation cost.

Assuming the inverse Hessian matrix is updated at an acceptable frequency, the Newton-Raphson iterations can still diverge. The time step is then rejected and the time step size is reduced. A problem is to determine when iterations diverge. Usually, a maximum number of iterations is defined. If this number is too small, a time step can be rejected while the problem slowly converges. If this number is chosen too large, some iterations are needlessly computed when the divergence actually occurs. It is then interesting to determine if divergence occurs on the basis of the evolution of the residual. The maximum number of iterations is more difficult to be correctly determined when the inverted matrix is not computed at each iteration. Indeed, this number depends on how frequently the inverted matrix is computed.

This paper proposes an automatic time step control algorithm based on the measure of the integration error. This algorithm modifies the time step size only if durable physical changes occur in the problem evolution. An algorithm choosing if the Hessian matrix is to be recomputed is also proposed. This determination is based on residual evolution with iterations. Finally, a divergence criterion based on this residual evolution is implemented. Industrial numerical examples are then presented to illustrate these new algorithms.

# NUMERICAL INTEGRATION OF TRANSIENT PROBLEMS

## Equations of motion

FEM (space) semi-discretization of the equations of motion of a nonlinear structure leads to the coupled set of second order nonlinear differential equations [1,2,3,4,5]:
$$R = M \ddot{x} + F^{\text{int}}(x, \dot{x}) - F^{\text{ext}}(x, \dot{x}) = 0 \quad [1]$$

This set is completed by two sets of given initial conditions, for both

$$x_0 = x(t_0) \text{ and } \dot{x}_0 = \dot{x}(t_0) \quad [2]$$

Internal and external forces can be written:

$$F^{\text{int}}(x, \dot{x}) = \int_{V(t)} [B]^T \{\sigma\} dV \quad [3]$$

$$F^{\text{ext}}(x, \dot{x}) = \int_{V(t)} \rho [N]^T \{b\} dV + \int_{S(t)} [N]^T \{f\} dS \quad [4]$$

Note that expression (4) collects all types of loading (applied through local or distributed actions, in a follow-up way or not, reactions to imposed displacements and contact situations) and that the consistent mass matrix reads

$$M = \int_{V(t)} \rho [N]^T [N] dV = \int_{V_0} \rho_0 [N]^T [N] dV_0 \quad [5]$$

## Implicit schemes

The most general scheme for implicit integration of (1) is a generalized trapezoidal scheme [6] where updating of positions and velocities are based on “averaged” accelerations stemming from associated values between $t_n$ and $t_{n+1}$. It reads for instance

$$\dot{x}_{n+1} = \dot{x}_n + (1-\gamma) \Delta t \ddot{x}_n + \gamma \Delta t \ddot{x}_{n+1} \quad [6]$$

$$x_{n+1} = x_n + \Delta t \dot{x}_n + \left(\frac{1}{2} - \beta\right) \Delta t^2 \ddot{x}_n + \beta \Delta t^2 \ddot{x}_{n+1} \quad [7]$$

or equivalently

$$\ddot{x}_{n+1} = \frac{1}{\beta \Delta t^2} \left[ x_{n+1} - x_n - \Delta t \dot{x}_n - \left(\frac{1}{2} - \beta\right) \Delta t^2 \ddot{x}_n \right] \quad [8]$$

$$\dot{x}_{n+1} = \frac{\gamma}{\beta \Delta t} \left[ x_{n+1} - x_n + \left(\frac{\beta}{\gamma} - 1\right) \Delta t \dot{x}_n + \left(\frac{\beta}{\gamma} - \frac{1}{2}\right) \Delta t^2 \ddot{x}_n \right] \quad [9]$$

The discretized equations of motion (1) become:

$$R_{n,n+1} = \frac{1 - \alpha_M}{1 - \alpha_F} M \ddot{x}_{n+1} + \frac{\alpha_M}{1 - \alpha_F} M \ddot{x}_n + (F_{n+1}^{\text{int}} - F_{n+1}^{\text{ext}}) + \frac{\alpha_F}{1 - \alpha_F} (F_n^{\text{int}} - F_n^{\text{ext}}) = 0 \quad [10]$$

This general form was introduced, for linear problems, by Chung and Hulbert [6]. Particular choices of parameters lead to well-known [1,2,6] schemes such as

- $\alpha_M = \alpha_F = 0$ for Newmark scheme

- $\alpha_M = 0$ for Hilber-Hughes-Taylor scheme

- $\alpha_F = 0$ for Wood-Bossak-Zienkiewicz scheme
Unconditional stability and second order accuracy of the scheme for linear problems [1,2,6] require that

$$
\begin{align*}
\gamma &\ge \frac{1}{2} - \alpha_M + \alpha_F \\
\alpha_M &\le \alpha_F \le \frac{1}{2} \\
\beta &\ge \frac{1}{4} + (\alpha_F - \alpha_M)/2
\end{align*}
\tag{11}
$$

Associated rules for Newmark scheme are

$$
\begin{align*}
\gamma &\ge \frac{1}{2} \\
\beta &\ge \frac{1}{4} \left(\gamma + \frac{1}{2}\right)^2
\end{align*}
\tag{12}
$$

It's worth pointing out that though classical schemes require $0 < \alpha_f < 1/2$ (i.e. sampling the force in the second half of $[t_n, t_{n+1}]$), here, no such rule is followed for $\alpha_M$, the sampling parameter for the inertia terms: it might be negative for instance, thus leading to an interpolation at $t_{n+1}$ instead of an extrapolation.

**Iterative solution**

Iterative solution of the nonlinear system (10) requires first the elimination of acceleration and velocity at time $t_{n+1}$ with the help of (8) and (9) and, secondly, the writing of the Hessian matrix of the system, i.e.

$$
S = \left[ \frac{1}{\beta \Delta t^2} \left( \frac{1 - \alpha_M}{1 - \alpha_F} \right) M + \frac{\gamma}{\beta \Delta t} C_T + K_T \right] \quad [13]
$$

where $K_T, C_T$ are respectively the tangent stiffness and damping matrices defined by

$$
K_T = \frac{\partial}{\partial x} (F^{\text{int}} - F^{\text{ext}}) \quad [14]
$$

$$
C_T = \frac{\partial}{\partial \dot{x}} (F^{\text{int}} - F^{\text{ext}}) \quad [15]
$$

The residual for iteration number $i+1$ is defined by:

$$
R = \frac{1 - \alpha_M}{1 - \alpha_F} M \ddot{x}_{n+1} + \frac{\alpha_M}{1 - \alpha_F} M \ddot{x}_n + \left( F^{\text{int}}(x_{n+1}^t, \dot{x}_{n+1}^t) - F^{\text{ext}}(x_{n+1}^t, \dot{x}_{n+1}^t) \right) + \frac{\alpha_F}{1 - \alpha_F} (F_n^{\text{int}} - F_n^{\text{ext}}) \quad [16]
$$

Using equations (13) to (16), the iterative solution of system (8,9,10) can be written as:

$$
S \cdot \Delta x = -R \tag{17}
$$

Iterations stop when the non-dimensional residual *r* becomes lower than the accuracy tolerance *δ*.
Therefore, the following relation is verified:

$$
r = \frac{|R|}{|F^{\text{ext}}| + |F^{\text{int}}|} < \delta \qquad [18]
$$
# AUTOMATIC TIME STEP SIZE CONTROL

## Introduction

A relatively simple method, proposed by Ponthot [3], aims at an optimal number of iterations. If the number of iterations exceeds this optimal number, time step size is reduced, while, if the number of iterations is lower than the optimal number, time step size is augmented. Givoli and Henisberg [7] propose to modify the time step size to keep the displacements difference, between two successive times, lower than a given limit. Géradin [8] estimates the integration error from the accelerations and the inertial forces difference between two successive times multiplied by the square of time step size. This error is divided by a constant depending on the initial positions and by a constant that is the average error for a one-degree-of-freedom linear system. The error must be lower than a given tolerance. If the error is higher, the time step is rejected and its size is divided by two. If the error is lower than the tolerance but higher than half the tolerance, the time step is divided by the ratio between the error and the half tolerance to the power one third. If the error is lower than the tolerance divided by sixteen, the time step size is multiplied by two. For Cassano and Cardona [9], the time step control is the same than for Géradin, but the error is calculated only from the accelerations difference and is not divided by a constant depending on the initial positions but by a term that evolves with positions. Hulbert and Jang [10] estimate the error from the accelerations difference multiplied by the square of time step size. This error is divided by a term that depends on the positions difference. Time step control is characterised by two tolerances (TOL1 and TOL2) and by a counter of maximal index LCOUNT. If the error is higher than TOL2 then the step is rejected and time step size is reduced. If error is lower than TOL2 and higher than TOL1, the time step is accepted and time step size is kept constant. If the error is lower than TOL1 then time step is accepted. If it occurs successively LCOUNT times, then the time step size is augmented. The counter is introduced to avoid undesirable change in time step size due to periodic nature of the local error. Dutta and Ramakrishnan [11] also calculate the error from the accelerations difference multiplied by the square of the time step size. It is made non-dimensional by dividing it by the maximum norm of the positions vector, for the previous time step. The time interval is divided in sub-domains. In each sub-domain there are a certain number of time steps of constant size. Once the time marching scheme has gone through a whole sub-domain, an average error is calculated. A time step size for the next sub-domain is then computed from this average error.

In this paper, the automatic control scheme is based on the algorithm proposed by Géradin [8]. Nevertheless, due to the non-linear characteristics of the problems we are interested in, we will assume that the time step size reacts only on evolution in physical mode and not on numerical mode. Changes in time step size will also occur only if the new time step size can be kept constant for several steps.

## Error estimator

This algorithm has been implemented in the module MECANO of the general-purpose finite element code SAMCEF [13]. The integration error used is the one developed by Géradin [8], which was already implemented in MECANO. This error is given by:

$$e = \frac{\Delta t^2}{6 \epsilon(\Omega_k) [x_0^T M x_0]^{1/2}} [\Delta \ddot{x}^T \Delta(M\ddot{x})]^{1/2} \quad [19]$$

In this expression, $x_0$ is the initial vector of positions and $\epsilon(\Omega)$ is the average error for a one-degree-of-freedom linear system with pulsation $\omega$, time step $\Delta t$ and non-dimensional pulsation $\Omega=\omega\Delta t$. It is evaluated for $\Omega_k=0.6$. This value gives correct results for a one-degree-of-freedom linear system. Solving (8,9,10) for a one-degree-of-freedom linear system, $\epsilon$ takes the following expression:
$$ \epsilon(\Omega) = \frac{(1 - \alpha_F) \Omega^3 \sqrt{1 + \frac{\Omega^2}{4}}}{3 \pi [1 - \alpha_M + (1 - \alpha_F) \Omega^2 \beta]} \quad [20] $$

### Remarks:

- Expression (19) come from the truncature error of equations (6,7). Indeed, the error is of the third order: $e = O(\Delta t^3 \ddot{x}) \approx O(\Delta t^2 \ddot{x})$. Géradin demonstrates that, for a linear multi-degree-of-freedom system, this expression can be written as (19). However, for non-linear systems, no advantage is gained in replacing the acceleration difference by a term depending on the accelerations and the inertial forces difference as in (19).

- The error is divided by expression (20) evaluated at $\Omega_k$ to have an estimation of the error that does not depend on integration parameters.

## Time step size control

The computed error must be of the order of a user-defined tolerance $PRCU$. A low $PRCU$ gives a good accuracy but a longer computation time. A higher $PRCU$ gives a shorter computation time but a lower accuracy. If $PRCU$ is too high, the time step size can result in an error lower than $PRCU$ but can be not small enough to allow for the iterations to converge. Therefore, if a problem of convergence appears, the algorithm reduces $PRCU$.

If the iterations converge, the algorithm tries to adjust the time step size to have an error equal to one half of $PRCU$. There exist three possibilities:

- The error is larger than $PRCU/2$. It is considered to be too high. To ensure a good accuracy, time step size must be reduced.

- The error is smaller than a limit $TRHLD$. It is considered to be too small. To ensure a reduced computation cost, time step size must be augmented.

- The error is in the interval $[TRHLD, PRCU/2]$. The time step size ensures a good accuracy with a relatively low computation cost. Time step size is kept constant

Let us first examine the problem of too high an error. The time step size must therefore be reduced. But to avoid needless changes of time step, we will assume that the variation of the integration error is due to a durable and physical evolution in the problem. The time step is then reduced only if there are a number (CO) of successive time steps for which integration error is larger than $PRCU/2$. This number CO can be taken equal to three. The factor by which the time step size is reduced depends on the maximum error (ERRO) of CO successive time steps. Géradin demonstrates that for a linear one-degree-of-freedom system, the factor by which the time step size needs to be multiplied to reduce the error from e to $PRCU/2$ can be written:

$$ RAT = \left[ \frac{PRCU}{2e} \right]^{1/\eta}, \quad \eta \in [2, 3] \quad [21] $$

For non-linear systems $\eta$ can be out of this interval. To ensure that the time step size is sufficiently reduced, $\eta$ is taken smaller than two. The factor that finally multiplies the time step size is $RAT=[0.5 PRCU/ERRO]^{2.3}$. But if there is a rapid change in the physical problem (impact...), the time step is not immediately adapted. Therefore, if the error e is larger than $PRCU$, the time step size is immediately multiplied by $RAT=[0.5 PRCU/e]^{2.3}$. If the error e is larger than a limit $REJL$, the time step size is multiplied by $RAT=[0.5 PRCU/e]^{2.3}$, but the time step is rejected. $REJL$ can be taken equal to $1.5 PRCU$.

Now, let us examine the problem of too small an error. The time step size could be augmented without degrading the solution. To avoid needless time step size changes, another counter is introduced. If $CT$ successive time steps give an error lower then the limit $TRHLD$, the time step size is then augmented.
ERRT is the maximal error of those CT steps. To ensure the time step size is not augmented too much, η from equation (21) is taken larger than 3. The factor multiplying the time step size is finally RAT=[0.5 PRCU/ERRT]$^{1.5}$. A problem due to the introduction of a counter occurs when the solution becomes smoother (external forces diminish...). Indeed, TRHLD must be taken small (e.g. PRCU/16) and CT relatively large (e.g. 5) to ensure a good accuracy. In those conditions, time step size augments slowly. To reduce the computation cost, TRHLD can be increased and CT can be decreased when the time step size is augmented. TRHLD can be multiplied by 1.3 while CT is reduced to 4 first and to 2 next. Once a time step size is reduced, TRHLD and CT are set back to their respective initial values PRCU/16 and 1.3. In some problems (translation at constant velocities), the error becomes nil. To avoid a division by zero, ERRT is limited by TRHLD PRCU/10.

## SELECTIVE UPDATING OF THE INVERSE HESSIAN MATRIX

For non-linear problems, if the Hessian matrix is not recomputed and inverted, the convergence of Newton-Raphson iterations is slower than if the Hessian matrix were recomputed and inverted at each iteration. For some step, divergence could also occur. Therefore, the criterion must consider two facts:

- Convergence of the iterations must be ensured.

- Not updating the Hessian matrix must reduce the total computation cost. Indeed, a problem with a small number of degrees of freedom and with strong non-linearities can converge in a few iterations when the Hessian matrix is updated at each iteration, but converge with more iterations when the Hessian matrix is not updated. When the number of degrees of freedom is reduced, an iteration without recalculation is not much less expensive. The total cost is then reduced when the Hessian matrix is often recalculated. On the other hand, if the problem has a large number of degrees of freedom and only a few non-linear elements, not updating the Hessian matrix can then reduce the computation cost.

The evolution of the non-dimensional residual *r* (18) could indicate if the problem converges or not. While *r* decreases, iterations converge even if the Hessian matrix is not recalculated and not inverted. An indication of how it could be interesting not to recalculate the Hessian matrix is the ratio *VALRF* between the time needed for an iteration with recalculation and an iteration without recalculation. This ratio indicates how much an iteration without recalculation could advantageously replace an iteration with recalculation. This value *VALRF* is an integer in the range 2 to 15.

The proposed algorithm is the following:

- The Hessian matrix is recalculated at the first iteration if the time step size has changed. Indeed, *S* depends on *Δt* (13). No convergence is possible if *S* is not recalculated.

- If the number of the iterations is greater than *VALRF*, the next iteration is made with recalculation of the Hessian matrix. Then, iterations occur without recalculation only if it is less expensive.

- If the number of the iterations is lower than *VALRF*, the Hessian matrix is recalculated only if the non-dimensional residual *r* has not been reduced by a ratio chosen equal to *RAPRES* = *VALRF*/10 ∈ [0.75, 0.95].

- If the non-dimensional residual has not been divided by *RAPRES*, the next iteration then needs recalculation of the Hessian matrix. But ideally, this iteration does not take as initial values (*x*, *x*̇, *x*̈) the values at the end of the previous iteration, but the value at the end of the last iteration which has converged. Some divergences of the iterations are then avoided. For practical reasons, implementation of this last remark in MECANO was not possible. Thus, the following solution has been adopted. If the non-dimensional residual has not been divided by *RAPRES*, the next iteration occurs with recalculation and the initial values are the prediction values.

- If an iteration, with a number lower than *VALRF* and larger than one, with recalculation of the Hessian matrix occurs, all the subsequent iterations of this step will occur with recalculation. Convergence is then ensured.
- If the last iteration of the previous time step has needed recalculation of the Hessian matrix, the first iteration of the present step occurs with recalculation.

This algorithm avoids some needless recalculations and inversions of the Hessian matrix. For strongly non-linear problems with a small number of degrees of freedom, this algorithm is at worst as expensive as an algorithm with recalculation at each iteration. For problems with more degrees of freedom, this algorithm is less expensive than an algorithm where the user decides of the number of the iterations with recalculation. In fact, this algorithm allows a lot of iterations without recalculation when possible, and recalculates frequently the Hessian matrix when needed.

## CRITERION OF DIVERGENCE

Usually, the user specifies a maximum number of iteration. If upon reaching this number, the non-dimensional residual *r* is not lower than the tolerance *δ*, the time step is rejected and the time step size is divided. But when the residual *r* decrease slowly, the maximum number of iteration is exceeded before *r* becomes lower than *δ*. On the other hand, iterations can diverge after a few iterations. More iterations are then needless. Finally, if we accept the problem to be solved without recalculation of the Hessian matrix, the number of iterations is higher than when frequently recalculation occurs. A solution consists to consider that divergence occurs if, after 5 iterations with recalculation, the non-dimensional residual has not been divided by two.

## NUMERICAL EXEMPLES

The three algorithms (automatic time step size control, selective updating of the inverse Hessian matrix, divergence criterion) have been implemented in the dynamics module MECANO of SAMCEF [13]. In the commercial version of MECANO, time step size is chosen with the scheme of Géradin [8]. The user defines the number of the iterations with recalculation of the Hessian matrix and the maximum number of iterations. With the current algorithm, recalculation occurs at the first iteration of all time steps, according to the fact that time step size can change. Two industrial problems from SNECMA have been computed with the old and the new algorithm. These problems cannot be precisely described for confidentiality reasons. For the same reasons, axes are not graduated. Nevertheless we are authorised to tell that the problems were three dimensional models with thousands of degrees of freedom. Some elements are non-linear, simulating contacts, rupture... Problems are also non-stationary. When comparing the new and the old algorithms, precision parameters (*δ*, PRCU) are taken identical. For the old algorithm, the iterations number with recalculation of the Hessian matrix are chosen to minimise the total computation cost. Some attempts were necessary to define these parameters for having a low cost without divergence of the problems.

### Case 1

This problem is solved with the old algorithm and the new one. In both cases, tolerance $\delta$ on the non-dimensional residual *r* is taken equal at $10^{-3}$. Tolerance PRCU on the integration error is taken equal to $10^{-3}$. The initial time step size is the same. With the old algorithm, there is recalculation of the Hessian matrix for iterations 1, 3, 6, 7, 8, 9... Parameters of integration are the same for both methods ($\alpha_M = 0$, $\alpha_f = 0.1$). Figure 1 is showing the time evolution of a displacement degree of freedom of case 1. New algorithms give a solution nearly equal to the old ones. Figure 2 is the energy balance of case 1. The energy balance is the potential energy plus kinetic energy min external forces work. If this balance is positive, energy appears with times and the time integration is non-stable. If this balance is negative, energy disappears. It could be due to physical dissipation or numerical dissipation. If numerical dissipation is too high, the integration is not accurate. On Figure 2, we see that dissipation with new algorithms is a little lower (0.5%) than with the old algorithms. New algorithms thus give a good accuracy. Computation costs are given in
<figure><img src="image_1.png" /><figcaption>Figure 1 Displacement of case 1</figcaption></figure>

<figure><img src="image_2.png" /><figcaption>Figure 2 Energy balance of case 1</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Figure 3 Calculation costs of case 1</figcaption></figure>
## Case 2

This problem is solved with the old algorithms and the new ones. In both cases, tolerance $\delta$ on the non-dimensional residual $r$ is taken equal to $10^{-4}$. Tolerance $PRCU$ on the integration error is taken equal to $10^{-3}$. The initial time step size is the same. With the old algorithms, there is recalculation of the Hessian matrix for iterations 1, 3, 5, 6, 7, 8... Parameters of integration are the same for both methods ($\alpha_M = 0, \alpha_f = 0.1$). Figure 4 represents force evolution for a degree of freedom of case 2. The new algorithms give the same solutions than the old ones. Figure 5 shows the evolution of time step size. With the new algorithm of time step control, the time step size is constant during longer periods. Costly recalculations of the Hessian matrix because of time step size changes can thus be avoided. Figure 6 gives the calculation time needed for both methods. New algorithms reduce calculation time to 60% of the old ones.

<figure><img src="image_4.png" /><figcaption>Figure 4 Force for case 2</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 5. Time step size evolution for case 2</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 6 Calculation costs for case 2</figcaption></figure>

CONCLUSIONS

A new time step size control algorithm has been presented. This algorithm is based on the measure of an integration error. By introducing counters, the time step size is modified only for physical and durable variations in dynamical problems. But for a sudden changes such as an impact or a contact, the integration error increases in one time step and the algorithm reduces instantaneously the time step size. By modifying the limit under which the time step could be augmented, if the problem becomes smoother, time step size can increase rapidly. This algorithm thus gives a good accuracy with a low computation time and a constant time step for long period. If problems of convergence occur, tolerance on the integration error is reduced to adapt the time step size. Costly time step rejections are thus avoided. Next, an algorithm deciding if the Hessian matrix must be recalculated has been proposed. This algorithm recalculates the Hessian matrix only if it is necessary for convergence. If not, the old Hessian matrix is used in the iterative process and computation time is reduced. Finally a criterion of divergence was implemented. It considers that the problem does not converge if the non-dimensional residual does not decrease when iterating. A lot of needless iterations are thus avoided.

The proposed algorithms were implemented in MECANO of SAMCEF, and two industrial problems were computed. The solutions with the new algorithms are as accurate as with the existing algorithms, but computation costs were reduced to about 50%.

REFERENCES

1 T. Belytschko and T-J-R. Hughes, editors. *Computational Methods for Transient Analysis*. North Holland, 1983.

2 T.J.R. Hughes. *The Finite Element Method*. Prentice-Hall, 1987.

3 J-P. Ponthot. *Traité unifié de la Mécanique des Milieux Continus solides en grandes transformations par la méthode des éléments finis*. Ph D thesis (in French), Université de Liège, Liège, Belgium, 1995.

4 J-P. Ponthot and M. Hogge. *On relative merits of implicit / explicit algorithms for transient problems in metal forming simulation*. In International Conference on Numerical Methods for Metal Forming in Industry, Vol. 2, pp. 128-148, Baden-Baden, Germany, Sept. 1994.

5 M. Hogge and J-P. Ponthot. *Efficient implicit schemes for transient problems in metal forming simulation*. In NUPHYMAT'96, Numerical and Physical Study of Material Forming Processes, CEMEF - Ecole Nationale supérieure des Mines de Paris Sophia-Antipolis, France, June 5-7, 1996.

6 J. Chung and G-M. Hulbert. *A time integration algorithm for structural dynamics with improved numerical dissipations: the generalized-α method*. Journal of Applied Mechanics, Vol. 60, pp. 371-375, 1993.
7 D. Givoli and I. Henisberg. *A simple time-step control scheme*. Communication in Numerical Methods in Engineering, Vol. 9, pp. 873-881, 1993.

8 M. Géradin. *Analyse, simulation et conception de systèmes polyarticulés et structures déployables*. Cours IPSI Paris, 11-13 mars 1997.

9 A. Cassano and A. Cardona. *A comparison between three variable-step algorithms for the integration of the equations of motion in structural dynamics*. Latin American Research, Vol. 21, pp. 187-197, 1991.

10 G.M. Hulbert and I. Jang. *Automatic time step control algorithms for structural dynamics*. Computer Methods in Applied Mechanics and Engineering, Vol. 126, pp. 155-178, 1995.

11 A. DUTTA and C.V. Ramakrishnan. *Accurate computation of design sensitivities for structures under transient dynamic loads using time marching scheme*. International Journal for Numerical Methods in Engineering, Vol. 41, pp. 977-999, 1998.

12 M. Géradin and D. Rixen. *Théorie des vibrations. Application à la dynamique des structures*. Masson, Paris, 1993.

13 SAMTECH. Manuel d'utilisation de Samcef, v8.0. Liège, 1999.