# Implicit-explicit time integration algorithms for the numerical simulation of blade/casing interactions

Ludovic Noels, Laurent Stainier, Jean-Philippe Ponthot, Jérôme Bonini

► To cite this version:

Ludovic Noels, Laurent Stainier, Jean-Philippe Ponthot, Jérôme Bonini. Implicit-explicit time integration algorithms for the numerical simulation of blade/casing interactions. Structures under Shock and Impact VII, 2002, Montréal, Canada. 10.2495/SU020171. hal-01008292

HAL Id: hal-01008292

https://hal.science/hal-01008292

Submitted on 5 Nov 2021

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
# Implicit-explicit time integration algorithms for the numerical simulation of blade/casing interactions

L. Noels¹, L. Stainier¹, J.-P. Ponthot¹ & J. Bonini²

¹LTAS-Thermomechanics,

Liège University, Belgium.

²SNECMA-Moteurs, France.

## Abstract

In order to solve fast dynamic problems, an explicit method is the most adapted. But for slower dynamics, an implicit method is more stable. Typical industrial problems are governed by high frequency, e.g. impact, during short time intervals and slower dynamics, e.g. spring-back, during other time intervals. The optimal solution is then to have both implicit algorithm and explicit methods readily available in the same code and to be able to switch automatically from one to the other. Criteria that decide when to shift from a method to another have been developed here. Implicit balanced restarting conditions that annihilate numerical oscillations resulting for an explicit calculation are also proposed.

## 1 Introduction

Many industrial problems of today need to be simulated with non-linear models. This is especially true in the aeronautical domain where the constructor must ensure that components stand up to impact. Adequate non-linear numerical finite elements simulations can avoid a lot of expensive experiments. The choice of a time integration algorithm is an essential criterion to ensure efficiency and robustness of the numerical simulations. Difficulty in this choice resides in being able to combine robustness, accuracy and stability of the algorithm. Implicit ones need to be solved iteratively on each time increment (time step), contrarily to explicit ones. But, for stability reasons, explicit methods use smaller time steps than implicit ones. Explicit methods, avoiding iterations and convergence problems, are there-
fore generally used for fast dynamics simulation for which small time step size is always necessary to capture high frequency solutions. It is especially true for problems with many degrees of freedom, for which iterations are very expensive and convergence problems are frequent [1]. On the other hand, for slower dynamics problems, implicit algorithms allow to work with greater time step size, resulting in more numerical stability and accuracy [1, 2, 3]. Most industrial problems have some time intervals governed by fast dynamics and others governed by slower dynamics. Then, one can take advantage from a solution method that combines both families of integration algorithms.

A solution is to integrate over some time intervals with an implicit method and other time intervals with an explicit one. Few works have been developed with this latter combination and they were all developed for sheet metal forming analysis. Jung and Yang [4] have simulated a stamping simulation that begins with an implicit scheme and shifts to an explicit one when a problem of convergence appears. No return to implicit scheme is actually planned. Another method, developed by Finn et al. [5] and by Narkeeran and Lovell [6], simulates stamping (as a fast dynamics problem) with an explicit scheme and springback phase (slow dynamics) is subsequently analyzed with an implicit one. The time of transition is fixed by the user and initial conditions for the implicit phase, such as velocities and accelerations, are set to zero. This method has been generalized in this work and automatic criteria that decide to shift from a family to another have been developed. They depend on an integration error [7, 8] that allows to determine implicit time step size and they also depend on a ratio between the computational time (or CPU) needed to solve an implicit time step and the CPU needed to solve an explicit time step. Initial conditions, when shifting from explicit to implicit scheme occurs, are also defined to avoid lack of stability and convergence.

This paper will be organized into three sections. First, time integration algorithms will be briefly explained. Second, the mentioned criteria and initial restarting conditions will be detailed. Third, numerical simulations will be exposed to validate the methodologies.

# 2 Numerical integration of transient problems

## 2.1 Equations of motion

FEM (Finite Element Model) semi-discretization of the equations of motion of a nonlinear structure leads to the following coupled set of second order nonlinear differential equations [9, 10, 11]:

$$R = M\ddot{x} + F^{\text{int}}(x, \dot{x}) - F^{\text{ext}}(x, \dot{x}) = 0 \quad (1)$$

where *R* is the residual vector, *x* the vector of the nodal positions at current time, $\dot{x}$ the vector of nodal velocities, $\ddot{x}$ the vector of nodal accelerations. *M* is the mass matrix, $F^{\text{int}}$ the vector of internal forces resulting from body's deformation and $F^{\text{ext}}$ the vector of external forces. Both vectors are non-linear in *x* and in $\dot{x}$ due
to the coupled phenomena of contact, plastic deformations or geometrical non-linearities.

## 2.2 Implicit schemes

The most general scheme for implicit integration of eqn (1) is a generalized trapezoidal scheme [9, 10, 12] where updating of positions and velocities is 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 (2) $$

$$ 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 (3) $$

The discretized motion eqn (1) can be rewritten under the form proposed by Chung and Hulbert [12]:

$$ \begin{aligned} R_{n,n+1} = & \frac{1-\alpha_M}{1-\alpha_F} M \ddot{x}_{n+1} + \frac{\alpha_M}{1-\alpha_F} M \dot{x}_n + (F_{n+1}^{int} - F_{n+1}^{ext}) \\ & + \frac{\alpha_F}{1-\alpha_F} (F_n^{int} - F_n^{ext}) = 0 \end{aligned} \quad (4) $$

where $R_{n,n+1}$ is the residual vector of time step $n$ to $n+1$.

Iterative solution of the nonlinear eqn (4) requires 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 (5) $$

where $K_T$, $C_T$ are respectively the tangent stiffness and damping matrices. Using eqns (2) to (5) and a Newton-Raphson technique, the iterative solution of the problem can be written as:

$$ S \Delta x = -R \quad (6) $$

## 2.3 Explicit Scheme

Chung and Hulbert have extended their implicit scheme to an explicit one, taking $\alpha_F = 1$ in eqn (4) [13]. Its principal advantage is its numerical dissipation property. Time integration is then:

$$ \ddot{x}_{n+1} = \frac{M^{-1} (F_n^{ext} - F_n^{int}) - \alpha_M \ddot{x}_n}{1 - \alpha_M} \quad (7) $$

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

$$ x_{n+1} = x_n + \Delta t \dot{x}_n + \Delta t^2 \left[ \left( \frac{1}{2} - \beta \right) \ddot{x}_n + \beta \ddot{x}_{n+1} \right] \quad (9) $$
This scheme is conditionally stable and time step size is limited, depending on maximal model frequency $\omega_{max}$, but also depending on spectral radius ($\rho_b$):

$$ \Delta t = \gamma_s \Delta t_{crit} = \gamma_s \frac{\Omega_s(\rho_b)}{\omega_{max}} \quad (10) $$

with

$$ \Omega_s(\rho_b) = \sqrt{\frac{12(1+\rho_b)^3(2-\rho_b)}{10 + 15\rho_b - \rho_b^2 + \rho_b^3 - \rho_b^4}} \quad (11) $$

In eqn (10), $\gamma_s$ is a safety factor (< 1) that accounts for the destabilizing effects of non-linearities.

## 2.4 Implicit time step size control

The implicit time step size control is the one proposed by Geradin [7], extended to highly non-linear problems by Noels et al. [8]. This scheme continuously adapts time step size to physical modes evolution and keeps time step size constant during long time intervals. To estimate time step size an integration error is computed.

The integration error is deduced from truncated terms of eqn (2) and eqn (3). This error is of the third order: $O(\frac{1}{6}\Delta t^3 \ddot{x}) \approx O(\frac{1}{6}\Delta t^2 \Delta \ddot{x})$. To have a problem independent error, it is made non-dimensional, using $x_0$ (the initial position vector) and a reference error $\epsilon$ [7, 8]. To take into account the rotation, the integration error is then rewritten by taking the variation of the nodal acceleration modulus ($N$ is the number of nodes)[8]:

$$ e_{int} = \frac{\Delta t^2 \sum_{i=1}^{N} \Delta \|\ddot{x}_i\|}{6 \epsilon \|x_0\|} \quad (12) $$

Time step size is deduced from the integration error defined in eqn (12) and from a tolerance $PRCU$ fixed by the user. The relation to be verified is:

$$ e_{int} < PRCU \quad (13) $$

The new time step size $\Delta t_{new}$ to reach a reference integration error (half of the tolerance $PRCU$) is deduced from the current time step size ($\Delta t_{cur}$) and from the current integration error ($e_{int,cur}$), using the following relation developed by Geradin [7]:

$$ \left( \frac{\Delta t_{new}}{\Delta t_{cur}} \right)^n = \frac{PRCU}{2e_{int,cur}} \quad (14) $$

with $\eta \in [2, 3]$ a user specified parameter [7, 8]. The time step size management, based on eqns (13) and (14), is the one developed by Noels et al. [8].
# 3 Shifts from an algorithm family to another

## 3.1 Shift from an implicit algorithm to an explicit algorithm

First the ratio $r^*$ between the CPU needed for an implicit time step computation and the CPU needed for an explicit time step computation, is evaluated. In this paper, this ratio is averaged for each step to be able to shift from a method to another for non-linear simulation. Shift to explicit method occurs if:

$$ \mu \Delta t_{impl} < r^* \Delta t_{expl} \quad (15) $$

where $\Delta t_{expl}$ is evaluated with eqn (10). The factor $\mu$ is taken greater than unity (typical value is discussed in section 4) to avoid shifting from a method to another too frequently. This methodology allows to take into account the number of degrees of freedom, the algorithms efficiency, the residual tolerance required and the non-linearities evolution.

## 3.2 Shift from an explicit algorithm to an implicit algorithm

While the method used is an implicit one, the explicit time step size could always be easily computed from eqn (10). When the current method is explicit, the implicit time step size, which correctly integrates the problem, does not remain directly accessible. Using developments of section (2.4), nodal acceleration variations can provide us with this implicit time step size. Using eqn (14), acceleration variation is proportional to $\Delta t^n$. Inverting eqn (12) the implicit time step size is (with $N$ the number of nodes):

$$ \Delta t_{impl} = \left[ 6 \frac{\frac{\mu R C U}{2} \varepsilon \|x_\bullet\| (\Delta t_{expl})^{n-2}}{\sum_{i=1}^{N} \Delta \|\ddot{x}_i\|} \right]^{\frac{1}{n}} \quad (16) $$

Therefore the explicit to implicit shift criterion is similar to eqn (15). It yields:

$$ \Delta t_{impl} > \mu r^* \Delta t_{expl} \quad (17) $$

with $\Delta t_{expl}$ the current explicit time step size.

## 3.3 Initial conditions when shifting from an explicit scheme to an implicit scheme

Classical explicit schemes such as the central difference method [9] are well known to generate oscillatory (though stable) solutions. Two solutions are here provided to stabilize and balance the Gauss points values and the nodal values.

First, numerical oscillations of the Gauss points values and of the nodal values are annihilated thanks to the numerical dissipation property of the used $\alpha$-generalized explicit scheme. Indeed, when eqn (17) is satisfied, thus resulting in the choice to switch to implicit, at step number $n$ (at time $t_n$), $r^*$ explicit steps
<figure><img src="image_2.png" /><figcaption>Figure 1: Transition scheme from an explicit scheme to an implicit one.</figcaption></figure>

occur with a spectral radius $\rho_b$ (section 2.3) set equal to zero ($\rho_b$ is a user parameter). Thus, numerical oscillations have been greatly reduced at time $t_{n+r^*}$ (Figure 1).

The second step in the algorithm will be to determine a balanced configuration at time $t_{n+r^*+r^{**}}$. Therefore, we act twofold. First an explicit solution using $r^{**}$ ($r^{**}$ will be defined on next paragraph) explicit steps is computed. This solution results in $x_{n+r^*+r^{**}}^{\text{expl}}$, $\dot{x}_{n+r^*+r^{**}}^{\text{expl}}$ and in $\ddot{x}_{n+r^*+r^{**}}^{\text{expl}}$, which in turn is used as a predictor value for an implicit solution in one time step between time $t_{n+r^*}$ (where numerical oscillations have been reduced) and time $t_{n+r^*+r^{**}}$. This procedure proved to be very effective in order to restart an implicit solution based on explicit unbalanced solution. Therefore, a balanced step of size equal to the implicit time step size is reached. The methodology is illustrated on Figure (1). This balanced solution is reached considering an implicit time step size equal to $\Delta t_{r^{**}} = r^{**}\Delta t_{\text{expl}}$. In general, the iterative process necessary to reach this equilibrium quickly converges and this allows to begin the implicit method with a balanced solution at time $t_{n+r^*+r^{**}}$. Anyway $r^{**}$ must be defined. It is always lower or equal to $\mu r^*$. It is lower if $r^*$ is too large to lead to convergence of the first truly implicit step after time $t_{n+r^*+r^{**}}$. In this work $r^{**}$ is limited to 100. But if $r^{**}$ is lower than $\mu r^*$, time step size is increased (multiplied by 2) each two steps to reach $\Delta t = \mu r^*\Delta t_{\text{expl}}$.

# 4 Numerical examples

Numerical examples will be computed with three methods. The first method uses a purely implicit algorithm, the second method a purely explicit algorithm, and the third method can shift from an algorithm family to another. The accuracy of the numerical results and the computation costs are compared. The parameters $\eta$ of eqn (16) and $\mu$ of eqns (15, 17) are respectively taken equal to 2.5 and 1.5. Decreasing $\eta$ or $\mu$ will result in more shifts from a method to another and thus will degrade the efficiency of the algorithm. Since a return to an implicit scheme leads to some iterations (section 3.3), computation costs can increase.
<figure><img src="image_3.png" /><figcaption>Figure 2: Central point displacement comparison for Taylor bar impact.</figcaption></figure>

## 4.1 Numerical example 1: Taylor bar impact

A cylindrical bar (external diameter $d_e = 6.4mm$, length $l = 32.4mm$, density $\rho = 8930kg/m^3$, Young's modulus $E = 117e^3N/mm^2$, Poisson's ratio $\nu = 0.35$), discretized with 1080 elements (27 in each section by 40 along the axis), is at distance $x_0 = 1mm$ from a rigid plane. The plastic hardening law ($\sigma_v$ depending on effective plastic strain $\epsilon^{pl}$) is $\sigma_v = \sigma_0 + h\epsilon^{pl}$ with the yield stress $\sigma_0 = 400N/mm^2$ and the hardening parameter $h = 100N/mm^2$. It has an initial velocity $\dot{x}_0 = 227m/s$ directed towards the plane. After the wave has propagated across it, the bar loses contact with the wall. Numerical parameters used for the time integration scheme (section 2) are reported in Table 1.

Central point displacements obtained with the three methods (full implicit, full explicit and combined implicit/explicit) are compared in Figure 2. The combined method reacts as expected. It computes the rigid translation motion and rebound with an implicit scheme and the plastic impact with an explicit scheme. Comparing computed displacement, it appears that all the solution are quite similar. Computational costs (Table 2) confirm the marked advantage of the combined algorithm. It is about 30% cheaper than the implicit solution and about 700% cheaper than the explicit solution. When shift to implicit scheme occurs, using the defined initial conditions (section 3.3) does not lead to convergence problem.

<table><caption>Table 1: Numerical properties for the schemes.</caption><thead><tr><th></th><th>γ<sub>s</sub></th><th>ρ<sub>b</sub></th><th>α<sub>M</sub></th><th>α<sub>F</sub></th><th>β</th><th>γ</th></tr></thead><tbody><tr><td>explicit</td><td>0.8</td><td>0.2</td><td>-1.6</td><td>1.</td><td>5.5</td><td>3.1</td></tr><tr><td>implicit</td><td>-</td><td>-</td><td>-0.97</td><td>0.01</td><td>0.9801</td><td>1.48</td></tr></tbody></table>
<table><caption>Table 2: CPU comparison for Taylor bar impact.</caption><thead><tr><th></th><th>Implicit</th><th>Explicit</th><th>Implicit/explicit</th></tr></thead><tbody><tr><td>CPU time (min.)</td><td>11.45</td><td>66.2</td><td>8.3</td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>Figure 3: Squetch of the contact between a blade and a rigid casing.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4: Potential energies comparison for contact between a blade and a rigid casing.</figcaption></figure>

## 4.2 Numerical example 2: contact between a blade and a casing

We consider a blade (Figure 3), with a twisted rectangular cross section (casing diameter $c_d = 1.305m$, external blade diameter $d_e = 1.3m$, internal blade diameter $d_i = 0.1m$, length of section $L = 0.11m$, width of section $l = 0.032m$, twisting angle $\alpha_t = 45$, density $\rho = 4450kg/m^3$, young's modulus $E = 110e^9N/m^2$, Poisson's ratio $\nu = 0.31$), discretized with 204 elements, with an uniform rotation motion in a rigid casing. The rotational velocity is $\Omega_r = 3333.3r.p.m$. (Figure 3). The plastic hardening law ($\sigma_v$ depending on effective plastic strain $\bar{\epsilon}^{pt}$) is $\sigma_v = \sigma_0 + (\sigma_{\infty} - \sigma_0)(1 - e^{-h_e \bar{\epsilon}^{pt}}) + h \bar{\epsilon}^{pt}$ with the initial yield stress $\sigma_0 = 1e^9N/m^2$, saturated yield stress $\sigma_\infty = 1.3e^9N/mm^2$, exponential hardening
parameter $h_e = 100$, hardening parameter $h = 300N/mm^2$.

To simulate a load unbalance phenomenon, the center of rotation of the blade is moved during the first rotation so that the blade comes into contact with the casing and is subsequently deformed. The displacement of the center of rotation is the following. Its initial position is coincident with the center of the casing and there is no contact. During the first half revolution, the rotation center is moved 3cm to the opposite direction of the initial blade position so that contact interactions between the blade and the casing is generated. These interactions lead to irreversible plastic deformation of the blade. During the second half of the first revolution, the center of rotation is brought back to its initial position. Scheme parameters used (section 2) are given in Table 1. Potential energies (i.e. work done by internal forces) obtained with the three methods (full implicit, full explicit and combined implicit/explicit) are compared in Figure 4. Time intervals treated with implicit or explicit scheme when the combined implicit/explicit algorithm is used are indicated. This combined scheme behaves as expected. It computes rotation with no large plastic strain (before first half revolution and after first revolution) with an implicit scheme and interval with important plastic deformation (between first and second half revolution) with an explicit scheme. Comparing full implicit and full explicit computed energies, it appears that both solution are qualitatively identical but are about 25% different quantitatively. Nevertheless the implicit/explicit solution is bracketed by the two other solutions. It proves that the combined method does not lead to numerical instability. Computational costs (Table 3) confirm the advantage of the combined algorithm. It is about 8% cheaper than the implicit solution and about 50% cheaper than the explicit solution. When shift to implicit scheme occurs, using the defined initial conditions (section 3.3), it does not lead to any convergence problem.

<table><caption>Table 3: CPU comparison for contact between a blade and a rigid casing.</caption><thead><tr><th></th><th>Implicit</th><th>Explicit</th><th>Implicit/explicit</th></tr></thead><tbody><tr><td>CPU time (min.)</td><td>5.4</td><td>11.1</td><td>5.1</td></tr></tbody></table>

# 5 Conclusions

An integration scheme that combines implicit and explicit schemes was presented. This scheme integrates some time intervals with an implicit scheme, and others with an explicit scheme. First, automatic criteria that decide to shift from an algorithm family to another were developed. Next, stable balanced initial conditions have also been proposed when shifting from an explicit algorithm to an implicit algorithm. Finally, numerical examples were proposed that confirm the interest of the combined algorithm.
## Acknowledgements

L. Noels and L. Stainier wish to acknowledge support from the Belgian National Fund for Scientific Research (F.N.R.S.).

## References

[1] D.Y. Yang, D.W. Jung, I.S. Song, D.J. Yoo, and J.H. Lee. Comparative investigation into implicit, explicit, and iterative implicit/explicit schemes for the simulation of sheet-metal forming processes. *Journal of Material Processing Technology*, 50:39–53, 1995.

[2] J.C. Gelin, L. Boulmane, and P. Boisse. Quasi-static implicit and transient explicit analyses of sheet-metal forming using a c three-nodes shell element. *Journal of Material Processing Technology*, 50:54–69, 1995.

[3] J.S. Sun, K.H. Lee, and K.P. Lee. Comparaison of implicit and explicit finite element methods for dynamic problem. *Journal of Material Processing Technology*, 105:110–118, 2000.

[4] D.W. Jung and D.Y. Yang. Step-wise combined implicit-explicit finite-element simulation of autobody stamping process. *Journal of Material Processing Technology*, 83:245–260, 1998.

[5] M.J. Finn, P.C. Galbraith, L. Wu, J.O. Hallquist, L. Lum, and T.-L. Lin. Use of a coupled explicit-implicit solver for calculating spring-back in automotive body panels. *Journal of Material Processing Technology*, 50:395–409, 1995.

[6] N. Narkeeran and M. Lovell. Predicting springback in sheet metal forming: an explicit to implicit sequential solution procedure. *Finite Elements in Analysis and Design*, 33:29–42, 1999.

[7] M. Geradin. *Analyse, simulation et conception de systèmes polyarticulés et déployables*. Cours IPSI, Paris, mars 11-13 1997.

[8] L. Noels, L. Stainier, and J.P. Ponthot. Self-adapting time integration management in crash-worthiness and sheet metal forming computations. Accepted for publication in International Journal of Vehicle Design, 2001.

[9] T. Belytschko and T.J.R. Hughes. *Computational methods for transient analysis*. North Holland, 1983.

[10] T.J.R. Hughes. *The finite element method*. Prentice Hall, 1987.

[11] J.-P. Ponthot and M. Hogge. On relative merits of implicit schemes for transient problems in metal forming simulation. In *International Conference on Numerical Methods for Metal Forming in Industry*, volume 2, pages 128–148, Baden-Baden, Germany, Sept 1994.

[12] J. Chung and J.M. Hulbert. A time integration algorithms for structural dynamics with improved numerical dissipations: the generalized-α method. *Journal of Applied Mechanics*, 60:371–375, 1993.

[13] G.M. Hulbert and J. Chung. Explicit time integration algorithms for structural dynamics with optimal numerical dissipation. *Computer Methods in Applied Mechanics and Engineering*, 137:175–188, 1996.