7. Adaptive time step diagram#

7.1. Principle#

Explicit calculation methods are particularly suitable for the simulation of rapid phenomena, such as the propagation of waves in solids. On the other hand, they are less suitable for slower phenomena since the condition of stability of the diagram imposes a time step of the order of the smallest natural period of the system.

The adaptive schema, based on the schema of centered differences, was developed to allow the calculation of transient responses in which « fast » and « slow » phenomena. For example, during an impact, at first, high-frequency waves propagate and dissipate in the structure. Then, the structure only responds to its low frequency modes, the high frequencies having been dampened. The idea is therefore to adapt the time step gradually according to the phenomena involved, by setting a precision criterion on the solution.

7.2. Diagram#

The explicit diagram of centered differences with varying steps is written as:

\(\begin{array}{c}{\dot{X}}_{n\text{+}\frac{1}{2}}={\dot{X}}_{n\text{-}\frac{1}{2}}+\frac{\Delta {t}_{n\text{-}1}+\Delta {t}_{n}}{2}{\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}\)

with the following notations:

It can be seen that the time step varies. It is indexed: \(\Delta {t}_{n}\).

This means that the schema is no longer strictly of the second order, since it is no longer « centered ». The more different \(\Delta {t}_{n-1}\) and \(\Delta {t}_{n}\) are, the closer the order of the diagram is to one. Strong variations in the time step are therefore accompanied by a decrease in precision. The speed formula used leads to good results when the time step decreases but lowers the stability limit when the time step increases. This is why it is forced to increase only very gradually.

Finally, the same approximations are used as for the centered differences in terms of calculating accelerations and speeds at « integer » time steps:

acceleration is estimated by \({\ddot{X}}_{n}({t}_{n},{X}_{n},{\dot{X}}_{n})\stackrel{~}{=}{\ddot{X}}_{n}({t}_{n},{X}_{n},{\dot{X}}_{n-})\) and \({\ddot{X}}_{n}({t}_{n},{X}_{n},{\dot{X}}_{n\text{-}})={M}^{\text{-}1}\mathrm{.}(F({t}_{n})-\mathrm{K.}{X}_{n}-C{\dot{X}}_{n\text{-}})\);
and the stored speed is evaluated by \({\dot{X}}_{n\text{+}1}={\dot{X}}_{n\text{+}}+\frac{\Delta t }{2}{\ddot{X}}_{n\text{+}1}\).

As with the centered differences scheme, from which it is inspired, the adaptive step pattern requires the mass matrix to be inverted. This is why we require the diagonalization of the mass matrix as well as the same restrictions on the Lagrange degrees of freedom as for the centered differences diagram.

7.3. Estimation of the time step according to the required precision#

To define a criterion on the time step according to the precision required on the solution, we introduce the concept of apparent disturbed frequency [bib4]:

\({f}_{{\text{AP}}_{n}}=\frac{1}{2\pi }\sqrt{\mid \frac{{\ddot{X}}_{x}-{\ddot{X}}_{x\text{-}1}}{{X}_{x}-{X}_{x\text{-}1}}\mid }\)

This quantity can be interpreted as the « instantaneous frequency » of the system. It is in fact an approximation of the local slope of the force/displacement curve. It is related to the error due to truncation in the limited expansions of the diagram. It also makes it possible to take account of external forces and their frequency fluctuations.

In the case of a system with several degrees of freedom, it is necessary to calculate an apparent frequency for each degree of freedom. The maximum over all the frequencies calculated is then used to determine the time step.

If the denominator approaches zero, the apparent frequency may become very large and lose its physical meaning. An unjustified refinement of the time step is then obtained when the speed is cancelled out. In the case of sinusoidal oscillations, this is the case twice per period. The criterion is then modified by introducing the following condition:

\(\frac{\mid {X}_{x}-{X}_{x\text{-}1}\mid }{\Delta t }\le {\dot{X}}_{\text{min}}\Rightarrow {f}_{{\text{AP}}_{n}}=\frac{1}{2\pi }\sqrt{\mid \frac{{\ddot{X}}_{x}-{\ddot{X}}_{x\text{-}1}}{{\dot{X}}_{\text{min}}\Delta t }\mid }\)

An intermediary between the disturbed apparent frequency and the truncation error is obtained. The value of \({\dot{X}}_{\text{min}}\) is not easy to determine*a priori* and a poorly chosen value can lead to an artificial moderation of the apparent frequency.

Two methods are proposed.

7.3.1. Influences of neighboring nodes#

In the case of a system with several degrees of freedom, we can use the information given by the \(1\le j\le \text{nv}\) nodes neighboring node \(i\):

\({f}_{{\text{AP}}_{n}}=\underset{\text{DX},\text{DY},\text{DZ},\text{DRX},\text{DRY},\text{DRZ}}{\text{max}}(\underset{1\le i\le \text{nb}\text{noeud}}{\text{max}}\frac{1}{2\pi }\sqrt{\mid \frac{{\ddot{X}}_{x}^{i}-{\ddot{X}}_{x\text{-}1}^{i}}{{b}_{n}^{i}}\mid })\)

Where \({b}_{n}^{i}=\Delta {t}_{n}\text{max}({\text{10}}^{-\text{15}}{\text{ms}}^{-1}\text{,}{\dot{X}}_{n\text{+}}^{i}\text{,}\frac{1}{\text{100}}\underset{1\le j\le \text{nv}}{\text{max}}({\dot{X}}_{n\text{+}}^{j}))\)

This method requires the identification of neighboring nodes and the estimation of « speeds » according to each type of degree of freedom (translation DX, DY, DZ, and possibly rotation DRX, DRY and DRZ) for these neighboring nodes.

The programmed method simplifies this formula somewhat and consists, for a given degree of freedom, \(i\), in doing an upward search and a downward search on the degrees of freedom in their numbering order defined by NUME_DDL from this position. The first two degrees of freedom, \(k\) and \(l\), of the same nature, found respectively before and after the degree of freedom \(i\), are considered to be the « neighbors ». To limit the cost of this technique, the search is done once and for all at the beginning of the transitory calculation and the « neighbors » are recorded in two integer tables.

The use of this method is started with the keyword VITE_MIN =” NORM “.

7.3.2. Use of information in the previous time#

We can also rely on the information provided by the previous time steps to estimate the minimum speed. It is then estimated by the following formula:

\({\dot{X}}_{\text{min}}=\underset{k<n}{\text{max}}(\frac{\mid {\dot{X}}_{k}^{i}\mid }{\text{100}},{\text{10}}^{-\text{15}}{\text{ms}}^{\text{-}1})\)

We then have:

with \({b}_{n}^{i}=\Delta {t}_{n}\text{max}({\text{10}}^{\text{-15}}{\text{ms}}^{\text{-1}},{\dot{X}}_{n\text{+}}^{i},\frac{1}{\text{100}}\underset{k<n}{\text{max}}({\dot{X}}_{k}^{j}))\)

This method is triggered by the keyword VITE_MIN =” MAXI “.

This method cannot be used if the speed varies too much during the calculation, because in this case we would have at each step:

\(\frac{\mid {X}_{n}-{X}_{n\text{-}1}\mid }{\Delta t }\le {\dot{X}}_{\text{min}}^{i}\)

7.4. Choice of the number of steps per apparent period, \(N\)#

The error calculations and stability criteria established for a system with only one degree of freedom (see [bib4]) made it possible to estimate the number of \(N\) steps required per apparent period to obtain good precision. These tests showed that a minimum of 20 steps per period is necessary. This number can be set by the user in the command file using the “NB_POIN_PERIODE” keyword. Its default value, set at 50, leads to a precision on the temporal integration of the order of 1 to 2%.

The initial time step is used as the absolute maximum time step: \(\mathrm{\Delta }{t}_{\text{mac}}=\mathrm{\Delta }{t}_{\text{initial}}\). Weighted by a coefficient that can be set by PAS_LIMI_RELA, it is used as a minimum time step: \(\mathrm{\Delta }{t}_{\text{min}}=\text{PLR}\ast \mathrm{\Delta }{t}_{\text{initial}}\)

7.5. Heuristic for the evolution of the time step#

We define an indicator, called an « error », on the choice of the time step:

\(\text{erreur}=\mathrm{\Delta }{t}_{n}{\mathit{Nf}}_{{\text{AP}}_{n}}\)

This indicator must be less than 1 in order to hope to guarantee a good temporal integration of the smallest specific period. However, the adaptive schema must simultaneously avoid using too small a time step, which would then cause additional calculation costs, or even the appearance of parasitic « noises ».

Depending on the indicator, the algorithm will increase or decrease the time step. For this purpose, two coefficients are defined, \(\text{CDP}\), the time step refinement coefficient (keyword “COEF_DIVI_PAS”, default value: 1.334) and \(\text{CMP}\), the time step amplification coefficient (keyword “COEF_MULT_PAS”, key word “”, default value: 1.1).

During this search for the optimal time step, a maximum number of iterations for reducing the time step is defined, \({\text{iter}}_{\text{max}}\), in order to avoid the time step from evolving too abruptly, which is harmful to the order of the schema, and in order not to launch an optimization that is too expensive.

if the error indicator is greater than its limit value, if we have not exceeded the refinement limit number for a time step and if the time step remains greater than its fixed minimum value*a priori, we refine the time step: \(\mathrm{\Delta }{t}_{n}>\frac{1}{{\mathit{Nf}}_{{\text{AP}}_{n}}},\text{iter}<{\text{iter}}_{\text{max}}\text{et}\mathrm{\Delta }{t}_{n}>\mathrm{\Delta }{t}_{\text{min}}\Rightarrow \frac{\mathrm{\Delta }{t}_{n}}{\text{CDP}}\to \mathrm{\Delta }{t}_{n}\),

if the indicator shows that for five consecutive steps the time step has proved to be too fine, i.e. \(\mathrm{\Delta }{t}_{n}<\frac{0,\text{75}}{{\mathit{Nf}}_{{\text{AP}}_{n}}}\), then \(\text{min}\left(\mathrm{\Delta }{t}_{\text{min}},\text{CMP}\mathrm{\Delta }{t}_{n}\right)\to \mathrm{\Delta }{t}_{n}\)

7.6. Algorithm#

the algorithm was programmed in*Code_Aster* according to the following flowchart:

0	Initialization: \({X}_{0},{\dot{X}}_{0}\) data \({\ddot{X}}_{0}={M}^{\text{-1}}(F(t=0))-K{X}_{0}-C{\dot{X}}_{0}\) \({\dot{X}}_{\text{-}}={\dot{X}}_{0}-\frac{\Delta t }{2}{\ddot{X}}_{0}\) retrieving integration parameters: \(\Delta {t}_{\text{initial}}\) CMP time step amplification coefficient CDP time step reduction coefficient PLR limit to refinement such as \(\Delta t \ge \mathrm{PLR}\Delta {t}_{\text{initial}}\) N number of time steps per apparent period itermax maximum number of time step reductions
1	at each time step: \({X}_{n},{\dot{X}}_{n\text{-}},{\ddot{X}}_{n}\) known \({t}_{n\text{+}1}={t}_{n}+\Delta {t}_{n}\) 1.0 iter=0 1.1: temporal integration \({\dot{X}}_{n\text{+}\frac{1}{2}}={\dot{X}}_{n\text{-}\frac{1}{2}}+\frac{\Delta {t}_{n\text{-}1}+\Delta {t}_{n}}{2}{\ddot{X}}_{n}\) \({X}_{n\text{+}1}={X}_{n}+\Delta {t}_{n}{\dot{X}}_{n\text{+}\frac{1}{2}}\) \({\ddot{X}}_{n\text{+}1}={M}^{\text{-1}}\mathrm{.}(F({t}_{n\text{+}1})-K\mathrm{.}{X}_{n\text{+}1}-\mathrm{C.}{\dot{X}}_{n\text{+}\frac{1}{2}})\) \({\dot{X}}_{n\text{+}1}={\dot{X}}_{n\text{+}\frac{1}{2}}+\frac{\Delta {t}_{n}}{2}{\ddot{X}}_{n\text{+}1}\) 1.2 calculation of the apparent frequency and the error on the time step \({f}_{{\text{AP}}_{n}}=\underset{\text{DX},\text{DY},\text{DZ},\text{DRX},\text{DRY},\text{DRZ}}{\text{max}}(\underset{1\le i\le \text{nb}\text{noeud}}{\text{max}}\frac{1}{2\pi }\sqrt{\mid \frac{{\ddot{X}}_{x}^{i}-{\ddot{X}}_{x\text{-}1}^{i}}{{b}_{n}^{i}}\mid })\) \(\text{erreur}=\Delta {t}_{n}{\mathrm{Nf}}_{{\text{AP}}_{n}}\) 1.2 test on the relevance of the time step * if \(\text{erreur}>1\) and \(\text{iter}<{\text{iter}}_{\text{max}}\) So \(\Delta {t}_{n}/\text{CDP}\to \Delta {t}_{n}\) But if \(\Delta {t}_{n}<\Delta {t}_{\text{min}}\) stop the calculation with error message \(\text{iter}+1\to \text{iter}\) and back to 1.1 * if \(\text{erreur}>1\) and \(\text{iter}>{\text{iter}}_{\text{max}}\) then emission of an alarm and passage to point 2. * if \(\text{erreur}<1\) goes to point 2 with if \(\text{erreur}<0,\text{75}\) for 5 consecutive steps: amplifying the \(\mathrm{\Delta }{t}_{n}=\text{min}\left(\mathrm{\Delta }{t}_{\text{max}},\text{CMP}\mathrm{\Delta }{t}_{n}\right)\) time step
2	acceptance of the solution: possible archiving of \({\mathrm{X}}_{n\text{+}1},{\dot{\mathrm{X}}}_{n\text{+}1},{\ddot{\mathrm{X}}}_{n\text{+}1}\) then \(n+1\to n\): go back to 1 for the next step of time