4. Digital implementation

This chapter will discuss the choices to be preferred in terms of numerical methods for DYNA_NON_LINE [bib2].

In general, it is recommended to follow the following logic:

1. once the possibilities of almost static resolution through STAT_NON_LINE have been exhausted (including linear search and control),

2. start dynamic approaches with an implicit resolution,

3. in case of failure of the transient implicit strategies (including by combining different damping like Rayleigh and a HHT [bib2] and [bib3] scheme), we switch to explicit (having previously verified that the global damping does not depend on the stiffness matrix).

For each transitory approach, it is necessary to start with as little damping as possible, therefore, in particular, with non-dissipative time patterns (such as Newmark or centered differences [bib4]).

4.1. Discretization in time

Unlike the almost static, time has a physical meaning. Its discretization is all the more sensitive.

A few rules can be set out:

• the evolution of the imposed loads must be sampled in a sufficiently fine manner (between 5 to 10 time steps per shortest period of the signals in question),

• the modal behavior of the structure must be well represented (as above, we must have between 5 and 10 time steps per period the weakest of the modes considered).

Given the low frequency of the problems that we want to address here, these two rules are generally not very penalizing. However, compared to the time steps of the previous quasistatic calculations, the dynamic time steps can be quite significantly lower.

Explicitly, it is also necessary to respect the Current condition (CFL cf. [bib5], [bib2] and R05.05.05) under penalty of numerical discrepancy. For an integration diagram of the centered differences type, the critical time step is \(2/\omega\) with \(\omega\) which is the highest natural pulsation of the system.

We can calculate this pulsation with CALC_MODES by choosing the option “PLUS_GRANDE”.

For more details, the reader can refer to the U2.06.13 documentation.

For most structures, the Current condition is very penalizing: the speed of the waves is often of the order of a few thousand m/s, we arrive at time steps of less than 10-5 s.

4.2. Choice of integration schemes in time

Implicit schemes can be classified into two categories (we deliberately set aside first-order and/or speed schemes that are more specifically adapted to very irregular problems):

• average acceleration (NEWMARK [bib4]) of order 2 and which does not provide digital dissipation: to be used first,

• HHT complete (MODI_EQUI = “OUI” [bib3]) which remains of order 2, unlike the case of the modified average acceleration (MODI_EQUI = “NON”, option by default). This scheme is specifically developed to introduce high frequency digital damping and therefore not to disturb the low frequency physical response. Damping is directly controlled by parameter ALPHA in the diagram.

If we observe high frequency oscillations in the numerical solution (basically, oscillations whose period is of the order of a few steps of time), we can choose the complete HHT diagram, to start with a value of the order of —0.1 for the parameter ALPHA. A value of —0.3 is a high limit that can still be used.

If more damping is desired at medium frequency, then the modified mean acceleration scheme can be used.

Explicitly, we have two diagrams:

• centered differences (DIFF_CENT [bib4]) which is non-dissipative,

• Tchamwa-Wielgosz (TCHAMWA [bib6]) which is dissipative, in a manner comparable to HHT.

Here again, it is recommended to start by using a non-dissipative scheme.

4.3. Depreciation models

The order of introduction and use of dissipation in the discretized model is as follows:

1. intrinsic dissipation linked to non-linear behavioral relationships, to bonds (friction),

2. global dissipation such as structural damping (Rayleigh or modal),

3. Digital dissipation of the diagram in time.

Ideally, the first category should be sufficient, but in practice, for reasons of simplifying the model, it is often essential to add structural damping, schema damping being the last resort.

We will only discuss here the use of structural damping, in the Rayleigh sense, and the one linked to the diagram.

Let’s just remember that the more sources of dissipation we multiply, the more difficult will be their control and their physical interpretation.

4.3.1. Rayleigh damping

This model makes it possible to define the global damping matrix C as being a linear combination of stiffness and mass matrices (to have a diagonal damping matrix based on the usual dynamic modes):

\(C=\alpha K+\beta M\)

The U2.06.13 documentation explains it in detail.

The Rayleigh damping coefficients are defined, in terms of the characteristics of the material (command DEFI_MATERIAU), by the parameters AMOR_ALPHA and AMOR_BETA. The values to be imposed in order to obtain the desired damping x **** in the natural frequency interval f1 and f2 are deduced from the following equations:

Equation 1:

Equation 2:

Where f1 and f2 are the two natural frequencies limiting the study interval in question. In the context of this document, low-frequency solutions are sought, so the frequencies f 1 and f 2 will be associated with the first frequencies of the model, whose modes are consistent with the imposed load.

4.3.2. Depreciation due to the schema in time [bib2]

The documentation R05.05.05 and U2.06.13 present this aspect. Here we will limit ourselves to recalling the main trends.

To summarize, we can recall that, among the implicit schemes:

• the average acceleration pattern does not disappear,

• only the complete HHT diagram does not disturb the low frequency domain,

• for the same value of the parameter ALPHA the modified mean acceleration introduces much more dissipation than the diagram HHT.

In order to highlight the influence of the high frequency damping of implicit schemes, the U2.06.13 documentation presents examples of applications.

Finally, as far as explicit schemes are concerned, the rate of damping in the Tchamwa scheme is qualitatively close to that of the modified mean acceleration.

4.4. Adaptation of explicit methods

The conditional stability of explicit schemes makes them very unsuitable for simulating slow phenomena. Explicit resolution methods are not used here to capture fast phenomena such as wave propagation, but their use should be seen as a particular solver that can be adapted for slow problems.

In order to be able to increase the critical time step [bib2], we can increase the density of the structure (which decreases the speed of the waves in proportion to its square root):

However, it must be done gradually.

Indeed, two risks exist:

• if the time step becomes too large, the calculated solution may miss certain phenomena such as the appearance of shear bands and go so far as to branch off into a branch that is very different from the expected response,

• the increase in density may be limited by poor conditioning of the mass matrix.

As an indication, the maximum explicit time step (and therefore the maximum density) may be of the same order of magnitude as the time step required for the implicit transitory calculation, in any case, it must remain less than the quasistatic time step. These are crude trends, which do not exempt a parametric study on the explicit time step.

If the model has strong stiffness heterogeneities (definition of several materials), it may be appropriate to modify the densities separately, so as to have a relatively homogeneous Current condition between areas with different materials.

Important note

If we impose boundary conditions while traveling that change over time, we must take into account the fact that these conditions are in fact imposed in acceleration explicitly (because it is the primal unknown). This means that one must enter in DYNA_NON_LINEla second derivative of the moving signal that one wants to impose. This evolution of the imposed displacement must therefore be differentiable at least twice in time…

Finally, it is recommended to use a diagonal mass matrix (lumped), which is obtained by the keyword MASS_DIAG = “OUI” from DYNA_NON_LINE. As this calculation option is not available for all finite elements, the user may be forced to use the consistent mass, if necessary, as if implicitly.