6. Diagram of centered differences at a constant step#

6.1. Principle#

The explicit pattern of centered differences at a constant step is written as:

\(\begin{array}{c}{\dot{X}}_{n\text{+}\frac{1}{2}}={\dot{X}}_{n\text{-}\frac{1}{2}}+\Delta t{\ddot{X}}_{n}({t}_{n},{X}_{n},{\dot{X}}_{n})+o(\Delta {t}^{2})\\ {X}_{n\text{+}1}={X}_{n}+\Delta t{\dot{X}}_{n\text{+}\frac{1}{2}}({t}_{n},{X}_{n},{\dot{X}}_{n})+o(\Delta {t}^{2})\end{array}\)

in acceleration parametrization (on does not offer displacement parametrization), with the following notations:

_images/Object_110.svg

.

The speed is expressed at half-integer indices of the time discretization while the displacements and accelerations are expressed at integer indices.

Written this way, the diagram is of order two. It can also be written in the following way with the values over three successive steps:

\(\begin{array}{c}{\dot{\mathrm{X}}}_{n\text{+}\frac{1}{2}}={\dot{\mathrm{X}}}_{n\text{-}\frac{1}{2}}+\mathrm{\Delta }t{\ddot{\mathrm{X}}}_{n}\left({t}_{n},{\mathrm{X}}_{n},{\dot{\mathrm{X}}}_{n}\right)+o\left(\mathrm{\Delta }{t}^{2}\right)\\ {\mathrm{X}}_{n+1}={\mathrm{X}}_{n}+\mathrm{\Delta }t{\dot{\mathrm{X}}}_{n\text{+}\frac{1}{2}}\left({t}_{n},{\mathrm{X}}_{n},{\dot{\mathrm{X}}}_{n}\right)+o\left(\mathrm{\Delta }{t}^{2}\right)\end{array}\phantom{\rule{4em}{0ex}}\Rightarrow \phantom{\rule{4em}{0ex}}\begin{array}{c}{\dot{\mathrm{X}}}_{n+1}={\dot{\mathrm{X}}}_{n}+\frac{1}{2}\mathrm{\Delta }t\left({\ddot{\mathrm{X}}}_{n}+{\ddot{\mathrm{X}}}_{n+1}\right)\\ {\mathrm{X}}_{n+1}={\mathrm{X}}_{n}+\mathrm{\Delta }t\left({\dot{\mathrm{X}}}_{n}+\frac{1}{2}\mathrm{\Delta }t{\ddot{\mathrm{X}}}_{n}\right)\end{array}\)

The acceleration in \({t}_{n}\) is not immediately calculable because the speed is only known at the previous half-time (in \({t}_{n-1/2}\)), which poses a problem for evaluating the damping terms. To get around this difficulty, we calculate the acceleration in \({t}_{n}\) using the following approximation to the equation of motion:

\({\ddot{X}}_{n}({t}_{n},{X}_{n},{\dot{X}}_{n})\stackrel{~}{=}{\ddot{X}}_{n}({t}_{n},{X}_{n},{\dot{X}}_{n\text{-}})={M}^{-1}(F({t}_{n})-\mathrm{K.}{X}_{n}-C\mathrm{.}{\dot{X}}_{n\text{-}})\)

which is a valid approximation if the depreciation is low enough (\({\dot{X}}_{n}={\dot{X}}_{n\text{-}1/2}+o(1)\)). The diagram loses its second-order accuracy if the damping of the structure is significant.

Other methods of approximating the acceleration may be considered. The one chosen proved to be a good compromise between simplicity and stability, like the study described in reference [bib4] on the precision and stability of several methods.

The fields are archived at moments \({t}_{n},{t}_{;\text{+}1},\mathrm{...}\), the speed being approximated at these times by the following formula:

\({\dot{X}}_{n\text{+}1}={\dot{X}}_{n\text{+}1/2}+\frac{\Delta t}{2}{\ddot{X}}_{n\text{+}1}({t}_{n\text{+}1},{X}_{n\text{+}1},{\dot{X}}_{n\text{+}1/2})\)

The pattern of centered differences « shortens » the periods (frequency distortion). It maintains angular momentum. The total energy of the system oscillates over time steps.

6.2. Stability#

The pattern of centered differences is conditionally stable. In the case of a system without damping [bib2], the pattern is stable for a time step verifying \(\Delta t <\frac{2}{{\omega }_{\text{max}}}\) where \({\omega }_{\text{max}}\) is the largest natural pulsation of the system, i.e. \(\Delta t <\frac{{T}_{\text{min}}}{\pi }\). It takes a minimum of time steps to describe the smallest period of system \({T}_{\text{min}}\).

The limit value for the time step decreases slowly as the damping increases [bib4]. The stability criterion then requires:

\(\mathrm{\Delta }t<\frac{2}{{\omega }_{\text{max}}}\left(\sqrt{1+{\xi }^{2}}-\xi \right)\)

where \(\xi\) refers to reduced depreciation.

For example, for a depreciation of 0.5%, the condition becomes \(\Delta t <\frac{{T}_{\text{min}}}{5}\).

6.3. Algorithm#

In summary, the schema as introduced in the*Code_Aster* looks like this:

0

initialization:

\(\Delta t,{X}_{0},{\dot{X}}_{0}\) data \({\ddot{X}}_{0}={M}^{-1}(F(t=0))-K\mathrm{.}{X}_{0}-C\mathrm{.}{\dot{X}}_{0}\) \({\dot{X}}_{\text{-}1/2}={\dot{X}}_{0}-\frac{\Delta t}{2}{\ddot{X}}_{0}\)

1

At every step of time \({X}_{n},{\dot{X}}_{n\text{-}1/2},{\ddot{X}}_{n}\) known \(\begin{array}{c}{\dot{X}}_{n\text{+}1/2}={\dot{X}}_{n\text{+}1/2}+\Delta t{\ddot{X}}_{n}({t}_{n},{X}_{n},{\dot{X}}_{n\text{+}1/2})\\ {X}_{n\text{+}1}={X}_{n}+\Delta t{\dot{X}}_{n\text{+}1/2}\\ {\ddot{X}}_{n\text{+}1}={M}^{\text{-}1}(F({t}_{n})-\mathrm{K.}{X}_{n\text{+}1}-C\mathrm{.}{\dot{X}}_{n\text{-}1/2})\\ {\dot{X}}_{n\text{+}1}={\dot{X}}_{n\text{+}1/2}+\frac{\Delta t}{2}{\ddot{X}}_{n\text{+}1}\end{array}\)

2

possible archiving of \({X}_{n\text{+}1},{\dot{X}}_{n\text{+}1},{\ddot{X}}_{n\text{+}1}\) then go back to step 1) for the next step.

6.4. Diagonal mass matrix#

Calculating the acceleration requires the mass matrix to be inverted. This explicit scheme becomes more efficient if a concentrated mass lumping matrix is used in such a way that it is diagonal. The reversal then no longer requires factorization and is immediate.

This is why the centered difference scheme is only legal with mass matrices constructed diagonally, by the “MASS_MECA_DIAG” option of the CALC_MATR_ELEM operator.

6.5. Checking the time step#

We saw that the pattern of centered differences is stable provided that the time step, in the absence of damping, is less than a limit value, equal to \(\Delta t <\frac{2}{{\omega }_{\text{max}}}\). In practice, a time step is used that is equal to 5% to 20% of the critical time step. A time step test was therefore introduced which verifies that:

\(\Delta t <\mathrm{0,}\text{05}\frac{2\pi }{\underset{1\le i\le \text{nddl}}{\text{max}}(\sqrt{\frac{{k}_{\text{ii}}}{{m}_{\text{ii}}}})}\)

where \({k}_{\text{ii}}\) and \({m}_{\text{ii}}\) are the diagonal terms for the stiffness and mass matrices.

If this condition is not met, the user is stopped with a message telling him the maximum time step that he can use.

6.6. Acceleration calculation#

The acceleration is calculated as follows:

for each degree of freedom, we test whether the diagonal term of the corresponding mass matrix is zero.

  • if it is not zero, the acceleration term is calculated using the formula: \({\ddot{X}}_{n\text{+}1}={M}^{\text{-}1}(F({t}_{n})-K\mathrm{.}{X}_{n\text{+}1}-\mathrm{C.}{\dot{X}}_{n\text{-}})\)

  • if it is zero, the acceleration term is not calculated. This is the case for so-called Lagrange degrees of freedom. If they correspond to blocked degrees of freedom, it is lawful not to take into account the line in question and not to calculate its acceleration. In the case where the Lagrange degree of freedom was introduced to define a connection between two degrees of freedom, this no longer makes sense. The schema is therefore unusable and a test stops the execution with an explicit message.