Code_Aster ®
Version
6.4
Titrate:
Algorithms of temporal integration of operator DYNA_TRAN_MODAL
Date:
28/02/03
Author (S):
E. BOYERE, LIGHT A.C., G. JACQUART Key
:
R5.06.04-B Page
: 1/18
Organization (S): EDF-R & D/AMA, TESE, EDF-Pôle Industrie/CNPE of Tricastin
Handbook of Référence
R5.06 booklet: Dynamics in modal base
Document: R5.06.04
Algorithms of temporal integration of the operator
DYNA_TRAN_MODAL
Summary
This document describes the diagrams of temporal integration which are used to solve in the space of
modes of the problems of transitory dynamics in linear mechanics, with, for certain diagrams, the catch
in possible account of nonlocalized linearities of shocks type, frictions or fluid blade, and the use
possible of the under-structuring. Diagrams of NEWMARK, EULER, DEVOGELAERE, and a diagram with step of
adaptive time, ADAPT, are presented.
Handbook of Référence
R5.06 booklet: Dynamics in modal base
HT-66/03/005/A
Code_Aster ®
Version
6.4
Titrate:
Algorithms of temporal integration of operator DYNA_TRAN_MODAL
Date:
28/02/03
Author (S):
E. BOYERE, LIGHT A.C., G. JACQUART Key
:
R5.06.04-B Page
: 2/18
Count
matters
1 Introduction ............................................................................................................................................. 3
2 Method of temporal integration of a dynamic problem ................................................................ 3
2.1 Introduction ...................................................................................................................................... 3
2.2 Methods of implicit integration ...................................................................................................... 5
2.2.1 Introduction ............................................................................................................................. 5
2.2.2 Method of NEWMARK [bib1] .......................................................................................... 5
2.2.2.1 Presentation of the diagram ............................................................................................. 5
2.2.2.2 Complete algorithm of the method of NEWMARK .................................................... 6
2.2.2.3 Stability conditions of the diagram of NEWMARK .................................................... 6
2.2.2.4 Employment ......................................................................................................................... 7
2.2.2.5 Numerical damping of the implicit schemes ................................................... 7
2.3 Methods of integrations clarifies .................................................................................................... 7
2.3.1 Introduction ............................................................................................................................. 7
2.3.2 Diagram clarifies of modified Euler of command 1 ............................................................................ 7
2.3.2.1 Presentation ............................................................................................................... 7
2.3.2.2 Command and stability of the diagram ...................................................................................... 8
2.3.3 Method of Devogelaere-Fu .................................................................................................. 9
2.3.3.1 presentation ................................................................................................................ 9
2.3.3.2 Command and stability of the diagram ...................................................................................... 9
2.3.3.3 Employment ....................................................................................................................... 10
2.3.4 Diagram of integration to step of adaptive time ................................................................... 10
2.3.4.1 Introduction: interest of a step of adaptive time .................................................... 10
2.3.4.2 Diagram of the centered differences with constant step ................................................... 10
2.3.4.3 Adaptation of the diagram to the variable step of time .................................................... 11
2.3.4.4 Stability and precision of the diagram .............................................................................. 12
2.3.4.5 Criteria of adaptation of the step of time .................................................................... 13
2.3.4.6 Algorithm of the diagram of the centered differences with adaptive step ............................ 15
2.3.4.7 Comments on the parameters of the algorithm .................................................. 16
2.3.4.8 performance of the algorithm ..................................................................................... 17
3 Conclusion ............................................................................................................................................ 17
4 Bibliography ......................................................................................................................................... 18
Handbook of Référence
R5.06 booklet: Dynamics in modal base
HT-66/03/005/A
Code_Aster ®
Version
6.4
Titrate:
Algorithms of temporal integration of operator DYNA_TRAN_MODAL
Date:
28/02/03
Author (S):
E. BOYERE, LIGHT A.C., G. JACQUART Key
:
R5.06.04-B Page
: 3/18
1 Introduction
The goal of the transitory dynamic analysis is to determine according to time the response of one
structure, being given a loading external or boundary conditions functions of time, in
cases where the effects of inertia cannot be neglected.
In a certain number of physical configurations, one cannot be satisfied with a modal analysis
or harmonic and one must carry out a transitory analysis. It is in particular the case if:
· the history of the phenomenon has an importance in the study,
· if the external loading is complex (seism, excitations multi-components, etc…),
· if the system is nonlinear (plasticity, shocks, frictions, etc…).
The methods of transitory analysis which can then be used divide into two large
categories:
· methods known as of direct integration,
· the methods of Ritz, which include/understand inter alia the recombination of projections
modal.
The methods of direct integration are thus called because no transformation is carried out
on the dynamic system after the discretization by finite elements. They are presented in
document [R5.05.02], algorithms of direct integration of operator DYNA_LINE_TRAN.
The methods of Ritz, on the other hand, proceed to a transformation of the initial dynamic system, by
a projection on a subspace of the space of discretization departure. The resolution is done then on
a modified system, which, if it is reduced, gives access only one approximation of the answer of
real system.
Algorithms of temporal integration on a frame of reference generalized are used
to solve the dynamic problems in mechanics for linear structures, with catch in
account possible of nonthe localized linearities such fluid shocks, frictions or blades.
Certain algorithms allow moreover the under-structuring.
These algorithms are programmed in operator DYNA_TRAN_MODAL of Code_Aster [U4.53.21].
2
Method of temporal integration of a dynamic problem
2.1 Introduction
It is supposed that the equations governing the dynamic balance of the solids were discretized by
finite elements. One obtains a discrete system of equations which it is a question of integrating in time. For
that one chooses a discretization {T, I NR of the interval of time of the study [0, T] and one writes
I
}
the balance of the structure at the moment T.
I
In a general way these equations take the following form:
MR. X
& + CX +
=
T +
T
&
K X
R
T
T
ext. ()
R nl (X, X
T
&, X
T & T)
where
·
MR. X
& + CX +
=
T +
T
&
K X
R
T
T
ext. ()
R nl (X, X
T
&, X
T & T) are the matrix of mass of the system,
·
MR. X
& + CX +
=
T +
T
&
K X
R
T
T
ext. ()
R nl (X, X
T
&, X
T & T) are the matrix of rigidity of the system,
·
MR. X
& + CX +
=
T +
T
&
K X
R
T
T
ext. ()
R nl (X, X
T
&, X
T & T) are the matrix of damping of
system,
·
MR. X
& + CX +
=
T +
is the vector of the external forces,
T
&
K X
R
T
T
ext. ()
R nl (X, X
T
&, X
T & T)
·
MR. X
& + CX +
=
T +
is the vector of the nonlinear forces.
T
&
K X
R
T
T
ext. ()
R nl (X, X
T
&, X
T & T)
Handbook of Référence
R5.06 booklet: Dynamics in modal base
HT-66/03/005/A
Code_Aster ®
Version
6.4
Titrate:
Algorithms of temporal integration of operator DYNA_TRAN_MODAL
Date:
28/02/03
Author (S):
E. BOYERE, LIGHT A.C., G. JACQUART Key
:
R5.06.04-B Page
: 4/18
The matrix of damping C is in general difficult to evaluate because damping is often
function of the frequency. It is however frequent to simplify the catch in depreciation account in
employing the model of damping proportional, or model of Rayleigh: .
The methods of reduction of Rayleigh-Ritz are presented in the document [R5.06.01], Méthodes
of Ritz in linear and nonlinear dynamics.
If the term R (X, X &, X
& is not null, the technique of the pseudo-forces consists with
nl
T
T
T)
to project on the basis of linear system and to maintain the forces nonlinear with the second member.
technique of the pseudo-forces is always associated a diagram of explicit integration. This fact
taken into account of nonthe linearities is available only for explicit diagrams. The addition of not
linearities does not modify the form of the equations.
In the method of Ritz, the field of displacement X is replaced by its projection on the basis
T
modal such as X = where is the vector of the generalized co-ordinates and is the base
T
T
T
modal, generally reduced.
The projected dynamic system takes the following form:
MT
T
T
T
T
& + C & + K
= R
T + R
,
with
p
T
T T ext. () nl (T &, t&t &)
T
When the assumption of Basile does not apply (damping nonproportional), the matrix
of damping projected is not diagonal. The integration of the coupled system is done then
obligatorily with one of the three following diagrams: the implicit scheme NEWMARK, the diagram
clarify EULER or explicit diagram ADAPT.
The equation obtained is same form as the equation in X. So in the continuation of
T
T
document, one will use as well notation X for displacement in generalized co-ordinates
T
that for displacement in physical space. In the case of operator DYNA_TRAN_MODAL, it
acts of generalized co-ordinates.
Two classes of method can be distinguished in integration step by step from the equations
of balance, methods of explicit integration and methods of implicit integration.
That is to say the linear dynamic system according to integrating in time:
MR. X
& + CX +
=
T
T
&
K X
R
T
T
ext. ()
This differential connection of the second command can be brought back to a first order system:
NR U & = H U + F
T
T
T
X
I 0
0
I
0
T
where U =
=
, H =
, F =
T
X &, NR 0 M
- K - C
T
R
T
T
To integrate this differential equation, one uses a discretization {T, I NR, like one
I
}
formulate differences finished to express the derivative U
&.
T
One will call methods of integration clarifies the methods where only the derivative U
& utilizes
T
unknown factors at time T
. In this way determination of the sizes sought at the moment T
I 1
+
I 1
+
do not result from an inversion of system utilizing the operator H. If moreover, one carries out one
“farmhouse-lumping” in order to return the matrix M diagonal, the determination of is
particularly simple. They are there the principal characteristics of the methods of integration
explicit.
Handbook of Référence
R5.06 booklet: Dynamics in modal base
HT-66/03/005/A
Code_Aster ®
Version
6.4
Titrate:
Algorithms of temporal integration of operator DYNA_TRAN_MODAL
Date:
28/02/03
Author (S):
E. BOYERE, LIGHT A.C., G. JACQUART Key
:
R5.06.04-B Page
: 5/18
The implicit or semi-implicit methods utilize the discretization of U at one moment
T
posterior with T, generally T
. The determination of the variables thus passes by the resolution of one
I
I 1
+
system utilizing the operator H.
Two concepts are important: consistency, or the command of the diagram of integration, and stability.
The approximations used to obtain the differential operators define consistency, or
the command of the diagram of integration. One can indeed consider that the approximation with which one
displacement with each step of time obtains is related to the command of approximation of the derivative
first and seconds compared to time.
The study of stability of a diagram consists in analyzing the propagation of the numerical disturbances with
run from time. A stable diagram preserves a finished solution, in spite of the disturbances, whereas one
unstable diagram led to a numerical explosion or divergence of the solution.
To make a study of stability of a diagram of integration, one puts this last in the form of one
linear recursive system and one determine the particular characteristics of this system. If all them
eigenvalues of the operator of recursivity are smaller than 1 modulates some, the diagram is stable.
If not it is unstable.
The diagrams of integration clarifies are generally conditionally stable, which means that
the step of time must be sufficiently small to ensure the stability of temporal integration.
Certain implicit algorithms have the property to be unconditionally stable, which makes them
interest and makes it possible to use a step of arbitrarily large time.
The diagrams retained for operator DYNA_TRAN_MODAL are an implicit scheme, NEWMARK, and
three explicit diagrams, EULER, DEVOGELAERE and ADAPT (with step of adaptive time). The choice is done
by key word METHODE: “EULER”, “DEVOGE”, “NEWMARK”, or “ADAPT”.
2.2
Methods of implicit integration
2.2.1 Introduction
The implicit methods utilize the resolution of a matric system with the operator
previously definite. If the solids are supposed to be elastic linear, that results in the resolution
of a linear system to each step of time.
The advantage of these methods is their unconditional stability, which enables them to integrate the equations
dynamics with a step of relatively important time while representing it correctly
behavior of the modes low in frequency of the structure.
An implicit version of the method of NEWMARK, which was programmed in DYNA_TRAN_MODAL for
linear problems.
2.2.2 Method of NEWMARK [bib1]
2.2.2.1 Presentation of the diagram
NEWMARK introduced two parameters and for the calculation of the positions and speeds to the step
T + T
:
&X
= X +
&
T
.
+
[(1).&X +.&X
T
T
T
T
t+ T
]
2
1
X = X + T
.X & + T
-
.
X
. & .X
t+ T
T
T
T
& t+ T
2
+
Let us consider the equilibrium equations at time T + T
:
Mr. X
&
+ C.X
t+ T
& +
K.X =
R
T + T
t+ T
T +t
Handbook of Référence
R5.06 booklet: Dynamics in modal base
HT-66/03/005/A
Code_Aster ®
Version
6.4
Titrate:
Algorithms of temporal integration of operator DYNA_TRAN_MODAL
Date:
28/02/03
Author (S):
E. BOYERE, LIGHT A.C., G. JACQUART Key
:
R5.06.04-B Page
: 6/18
Let us defer the preceding relations while eliminating &
Xt T
+ and &Xt
T
+, it comes:
~
~
~
K. X
R where: K
K
.M
.C
T
T
has
has
+
=
=
+ 0
+ 1
~
R = R
t+ T
+ C. {A.X
1
T + A.
4 &
X
T + A.
5 &Xt} + Mr. (A.X
0
T + A.
2 &
X
T + A.
3 &Xt)
1
1
1
has =
=
=
=
- 1
0
(
has
has
has
. t2
)
1
(. T
)
2
(. T
)
3
2
with:
T
has =
has =
=
has
T
. 1 -
=
.
4
5
6
(
) has
T
7
- 1
2.
- 2
2.2.2.2 Complete algorithm of the method of NEWMARK
) initialization has
1) conditions
initial
X, X
0
& 0 and &X0
2) choices
of
T and, and calculation of the coefficients a1,… a8 (cf above)
3) to assemble the matrices of stiffness K, mass M and damping C
~
4) to form the matrix of effective rigidity K = K + A.M
0
+ A.C
1
~
5) to factorize
K
b) with each step of time
~
1) to calculate the effective loading R:
~
R = Rt+ T + Mr. (A.X
0
T + A.
2 &
Xt + A.
3 &Xt) + C. {A.X
1
T + A.
4 &
Xt + A.
5 &Xt}
~
~
2) to solve
K.
X
t+t = R
3) to calculate speeds and accelerations at time T + T
&X =
a.
t+ T
0 (X -
X
-
t+ T
T)
a.
2 &
X - A.
T
3 &Xt
&X
=
t+ T
&X + A.
T
6 &X + A.
T
7 &Xt+ T
4) calculation of the step of next time: return in b)1)
2.2.2.3 Stability conditions of the diagram of NEWMARK
Method of NEWMARK used in a rather widespread way in the field of mechanics, because it
allows to choose the command of integration, to introduce or not numerical damping, and has
a very good precision.
(
2 +)
1 2
It is unconditionally stable if: >.
0 5 and >
4
One introduces a numerical damping positive if > 12 and negative if < 12.
When = 12 and = 0, the formula of NEWMARK are reduced to the diagram differences
centered. It is thus then an explicit diagram.
A combination very often employed is = 12 and = 14, because it leads to a diagram
of a nature 2, unconditionally stable without numerical damping. In fact the choice was made in
operator DYNA_TRAN_MODAL. The diagram of Newmark of this operator is thus implicit.
Handbook of Référence
R5.06 booklet: Dynamics in modal base
HT-66/03/005/A
Code_Aster ®
Version
6.4
Titrate:
Algorithms of temporal integration of operator DYNA_TRAN_MODAL
Date:
28/02/03
Author (S):
E. BOYERE, LIGHT A.C., G. JACQUART Key
:
R5.06.04-B Page
: 7/18
2.2.2.4 Employment
In DYNA_TRAN_MODAL, this diagram allows integration only linear problems. In
tally of the dynamic under-structuring, it allows to employ a modal base calculated by under
structuring but it does not support direct calculation on the basis of modal substructure.
2.2.2.5 Damping
numerical of the implicit schemes
The numerical advantage of the direct diagrams of implicit integration lies in the fact that the step of
time can be substantially large compared to the smallest clean period of the system without
to be likely to cause an instability of the results.
For modes of period clean about the step of time or lower than the step of time, them
algorithms of integration introduce a strong damping which contributes to erase the contribution of
high modes (cf [R5.05.02]).
It there not of numerical damping in the particular case of the algorithm of NEWMARK with = 14
and = 12.
On the other hand, implicit algorithms one a significant effect of lengthening of the periods of the answer
structure.
One notes that to guarantee a good precision on the amplitude and the phase of displacements
calculated, it is necessary to respect a criterion close to:
1
1
T < (
with
10 * F
100
max)
(* Fmax)
where F
is the high frequency of the movement which one wishes to capture.
max
2.3
Methods of integrations explicit
2.3.1 Introduction
Three methods of integration clarifies are presented: a diagram of modified Euler of command 1, one
diagram of Devogelaere-Fu of command 4 and one diagram with step of adaptive time ADAPT. These three
methods are available in operator DYNA_TRAN_MODAL. The diagrams are presented in
considering that linear forces. However the taking into account of the nonlinear forces of
easily deduced with the technique from the pseudo-forces.
2.3.2 Diagram clarifies of modified Euler of command 1
2.3.2.1 Presentation
This diagram is commonly called “modified Euler” because it is about a very simple alternative but
conditionally stable of the diagram of Euler of command 1, which is, him, unstable. It is thus a diagram
often employed into explicit for mechanics. In Code_Aster, it is quite simply called
EULER.
This diagram was used in the module TRANSIS of POUX [bib3], code finite elements of beam, and
in code CADYRO [bib4] for the calculation of the lines of trees in rotation.
The diagram uses formulates it of Euler of command 1 to estimate the derivative in time, with a formula
of front Euler for the speed and a formula of Euler postpones for displacement, as follows:
&X
1
1
&X
T M -
=
+
-
-
+
+
(R K X C &X O T
N
N
N
N)
()
X
&
1 = X
+ T X 1 + O T
n+
N
n+
()
Handbook of Référence
R5.06 booklet: Dynamics in modal base
HT-66/03/005/A
Code_Aster ®
Version
6.4
Titrate:
Algorithms of temporal integration of operator DYNA_TRAN_MODAL
Date:
28/02/03
Author (S):
E. BOYERE, LIGHT A.C., G. JACQUART Key
:
R5.06.04-B Page
: 8/18
the algorithm is thus the following:
) initialization has: X, X
& given
0
0
b) with each step of time:
&X
1
1
&X
T M -
=
+
-
-
+
(R K X C &X
N
N
N
N)
X
&
1 = X
+ T X
n+
N
n+1
2.3.2.2 Command and stability of the diagram
the approximations used in obtaining this diagram are of command 1. One can thus consider
that the approximation with which one obtains displacement with each step of time is of command 1. It
acts of the consistency of the diagram.
If one puts the diagram of integration in recursive form by eliminating the terms speed, one obtains
the relation of following recurrence (without external force, nor damping):
X
1
2
2
0
1 + (M - K
T -
+
1 =
+
) X X
N
N
N
The eigenvalue of this diagram are for a system with a degree of freedom:
- - 1
2
T ±
- 1
2
I
T (- -
2
4
1
2
MR. K
MR. K
MR. K T)
2
=
if T <
.
2
- 1
MR. K
The module of the eigenvalues is worth 1. One realizes that one is in a limiting situation
but favorable. There will not be uncontrolled increase in the error. Without damping, one is right
on the terminal of stability. It can be an asset for the diagram: it does not introduce dissipation
numerical.
2
If T >
-
one can show that one of the two eigenvalues has a module larger than the unit
1
MR. K
and thus that the diagram is unstable.
2
The criterion of stability of diagram EULER is thus T <
-
.
1
MR. K
This study can be extended to a system with a finished number of degrees of freedom. In this case, the criterion
of stability becomes:
2
T <.
max
L `analysis can be refined by considering a system with damping [bib13].
Handbook of Référence
R5.06 booklet: Dynamics in modal base
HT-66/03/005/A
Code_Aster ®
Version
6.4
Titrate:
Algorithms of temporal integration of operator DYNA_TRAN_MODAL
Date:
28/02/03
Author (S):
E. BOYERE, LIGHT A.C., G. JACQUART Key
:
R5.06.04-B Page
: 9/18
2.3.3 Method of Devogelaere-Fu
2.3.3.1 Presentation
To present the algorithm of Devogelaere-Fu, shortened in DEVOGE in Code_Aster, it is put
dynamic problem in the form:
MR. X
& + CX =
T
where the matrix C is supposed to be diagonal.
T
& T
(
G, Xt)
Displacements and speeds are calculated as follows:
) initialization has
2
T
T
X
X
X
1
1
1
&
4
T,
T,
&
1
0
0
(Mr. G (X
0
0)
M
(
G
X
0
0)
M -
=
-
+
-
-
CX0)
-
2
8
2
&
-
-
X
1
1
1
1
4 4
T
4
T
T
T
T
1 =
(I m) C
(I +
M -
) C &X0 - G, X
1
1 +
(
G
, X
0
0)
-
-
-
2
2
2
b) with each step of time
T
T 2
X
1
1
1
&
4
,
,
4 &
&
1 = X
+
X +
M
T
-
-
T
-
N
N
(
G
X
N
N)
MR. G
X
MR. C
X
X
1
1 -
-
N
1
n+
2
24
N
N
N
2
2
2
2
&
-
T
X
1
1
1
4 4
1
(I MT) C &X
G T
T
-
=
+
+
+
N
(, X
N
N)
G
, X
MR. C
1
1
&X
-
N
n+
4
n+
n+
2
2
2
T
T 2
X
1
1
1
&
4
,
2
,
&
2 &
1 = X
+
X +
M
-
-
+
+
(
G T X)
MR. G T
X
MR. C X
X
N
N
N
N
N
1
1 -
+
N
1
2
6
N
N
n+
2
2
2
&
T
X
-
-
1
1
- 1
6.6
4
4
1 =
+
+
1
1 +
+
(I MT) C &X
(
G T, X
+
+)
G T
, X
MR. C
1
1
&X 1 &X
N
N
N
N
-
+
N
6
n+
n+
n+
2
2
2
2.3.3.2 Command and stability of the diagram
The diagram is of command 4, the approximations in the writing of the derivative temporal being in
(O t4).
It thus has an excellent aptitude for the integration of regular solutions. Its interest is in
less manifest revenge if the functions to be integrated present discontinuities (shocks, friction,
etc)
One can show that for a linear system not deadened the step of time guaranteeing stability
2 2
is worth: T <.
max
Handbook of Référence
R5.06 booklet: Dynamics in modal base
HT-66/03/005/A
Code_Aster ®
Version
6.4
Titrate:
Algorithms of temporal integration of operator DYNA_TRAN_MODAL
Date:
28/02/03
Author (S):
E. BOYERE, LIGHT A.C., G. JACQUART Key
:
R5.06.04-B Page
: 10/18
2.3.3.3 Employment
This method is expensive in calculating times because it twice requires the evaluation of the vector of
forces intern G, operation particularly heavy. Consequently it is used little in
mechanics for direct integration. On the other hand it is employed by the ECA [bib2] in the case of them
systems projected on modal basis.
This diagram allows the taking into account of nonlocalized linearities of shocks type and frictions.
Within the framework of the dynamic under-structuring, it makes it possible to employ a modal base calculated by
under-structuring but it does not support direct calculation on the basis of modal substructure.
2.3.4 Diagram of integration to step of adaptive time
2.3.4.1 Introduction: interest of a step of adaptive time
To carry out the temporal integration of the transient of a structure in a nonlinear phase poses
always problems as for the choice of the step time. The estimate of the error is seldom accessible
during integration.
The diagrams of explicit integration oblige to respect a step of maximum time for not
to diverge. In the case of nonlinear behavior, this step cannot be a priori given and can
to change with each iteration. When rigidity very strongly varies, a step of constant time and very
end to preserve the stability of the diagram led to a very large iteration count and to a time of
considerable calculation.
An algorithm of integration to step of adaptive time was thus developed for DYNA_TRAN_MODAL
and was named ADAPT. It is based on the diagram of centered differences, of command 2.
One can notice that this type of diagram was also programmed in DYNA_LINE_TRAN
(cf [R5.05.02]).
2.3.4.2 Diagram of the centered differences with constant step
One presents initially the diagram of the differences centered at constant step on which diagram ADAPT
bases itself.
It is written as follows:
&X
2
1 = &
X 1 + T &X T
+ O T
N (
, X,
N
N
&Xn) (
)
n+
N
2
2
X
2
&
1 = X
+ T X 1 + O T
n+
N
()
n+ 2
with the following notations:
xn
&x
- 1
X
&x
n-1
X
2
N
n+12
N +1
T
T1
T
T
n-1
N
N
1
T
2
N + 2
n+1
T
T
Handbook of Référence
R5.06 booklet: Dynamics in modal base
HT-66/03/005/A
Code_Aster ®
Version
6.4
Titrate:
Algorithms of temporal integration of operator DYNA_TRAN_MODAL
Date:
28/02/03
Author (S):
E. BOYERE, LIGHT A.C., G. JACQUART Key
:
R5.06.04-B Page
: 11/18
It is noted that the speed is expressed with indices half entireties of the discretization in time then
that displacements and accelerations are expressed with the whole indices. Writing in this way it
diagram is of command 2. However acceleration is not immediately calculable because speed is not
known that with the half not preceding. To circumvent this difficulty, one can use several
approximations speed to the step of whole time.
· method 1: to suppose that &X
T
T what constitutes one
N (X,
N
&X,
N
N)
&X X,
X
N
&X
,
1
N
N
- 2
approximation validates if damping is sufficiently weak (&
X =
1). If
1 + O
N
&X
()
N 2
damping is important, the diagram loses then its precision of command 2.
· method 2
: to use an approximation of command 1 for speed
:
&
T
X =
what makes it possible to preserve command 2 of the diagram.
1 +
- 1 + O
T
N
&X
&X
()
N
N
2
2
· method 3: to use a diagram of correct the predictor type/
&X p = &X
1 +
T &X
N
N
- 1
predictor:
N
- 2
&X p
p
=
T
N
& (
X
, X,
N
N
&Xn)
&
T
X
p
=
1
1 +
+ -
N
&X
(&Xn () &Xn-1)
corrector:
N
-
2
2
&X =
T
N
& (
X
, X,
N
N
&Xn)
where and are two parameter to be chosen. Park and Underwood [bib6] report that to carry out
additional iterations does not improve in a significant way the stability of the diagram.
2.3.4.3 Adaptation of the diagram to the variable step of time
When the step of time varies, the expressions of the preceding paragraph are not valid any more,
acceleration &X being more necessarily expressed in the center of the interval &
X
,
,
1
&X
N
1
N
n+
2
2
as one sees it on the diagram below:
T1 +
-
T
N
N
2
X
X
n-1
&x
N
&x
n-1
N + 1
X
2
2
N + 1
T
T
T
n-1
1
T
1
T
N
n+
n+1
2
N
2
T
T
n-1
N
Handbook of Référence
R5.06 booklet: Dynamics in modal base
HT-66/03/005/A
Code_Aster ®
Version
6.4
Titrate:
Algorithms of temporal integration of operator DYNA_TRAN_MODAL
Date:
28/02/03
Author (S):
E. BOYERE, LIGHT A.C., G. JACQUART Key
:
R5.06.04-B Page
: 12/18
To take account of this, speed is calculated as follows:
T
1 +
T
&X
.
1 = &
X
N
N
1 +
&Xn
n+
N
2
2
2
Complete diagram ADAPT is written then as follows:
1. estimate of &
X according to methods 1, 2 or 3
X
T
1 +
T
2.
&X
2
1 = &
X
N
N
1 +
&X T
+ O T
N (
, X,
N
N
&Xn) (
)
n+
N
2
2
2
3.
X
2
&
1 = X
+ T X 1 + O T
n+
N
()
n+ 2
The command of the diagram is not rigorously any more equal to 2, the diagram having lost its centered character.
More T and T
are different, more the command of the diagram tends towards 1. Strong variations of the step of
N
n+1
times thus lead to a loss of precision.
It is possible to find expressions more complex, which use speed or acceleration with
the preceding iteration [bib7]. However the formula presented here gives satisfactory results
when the step of time decreases but it cause a drop in the limit of stability when the step of time
increase. The remedy is to control the step so that it increases only slowly.
2.3.4.4 Stability and precision of the diagram
To study the diagram, one was satisfied with the analysis of a system to only one degree of freedom, free and
linear, of own pulsation and reduced damping:
&
X + 2
&x + 2
X = 0.
The approached solution, by using the diagram with step of constant, is obtained by the relation of recurrence
following:
WITH Y + B Y -
0
1 =
N
N
X
& N &
with Y = X
N
& 1
n+
2
X
n+1
With and B are two matrices which depend on the selected method to calculate the contribution of
term of damping.
A solution of the form is sought: Y = Y.
N
n-1
is an eigenvalue of A-1 B and can be written in the following form:
= exp
T
- ± I 1 - 2 where and are the calculated pulsation and reduced damping
C
(C (C
c)
C
C
by the algorithm.
Handbook of Référence
R5.06 booklet: Dynamics in modal base
HT-66/03/005/A
Code_Aster ®
Version
6.4
Titrate:
Algorithms of temporal integration of operator DYNA_TRAN_MODAL
Date:
28/02/03
Author (S):
E. BOYERE, LIGHT A.C., G. JACQUART Key
:
R5.06.04-B Page
: 13/18
One can compare them with the exact solution, = exp
T
- ± I 1 - 2, which allows
E
((
)
-
-
C
C
D `to evaluate the error on the pulsation and the error on damping:
and
.
One studied [bib 8] and [bib10] the properties of the diagram according to the method employed to estimate
speed with the whole steps. It was empirically found that method 3 is at the same time more precise and
more stable than methods 1 and 2. The method, without overcost of calculation, makes it possible to increase the command of
diagram and gives in the majority of the cases a better precision, except in the event of weak damping.
It is however less stable than method 1. It is the method 2 which was finally adopted
in diagram ADAPT. These studies made it possible moreover to estimate the number of points per period
necessary to guarantee a stable integration. 20 is a value which gives a good margin of
security. It is the value chosen by defect.
2.3.4.5 Criteria of adaptation of the step of time
The preceding developments make it possible to quantify the errors introduced during the calculation of one
free and linear system. These criteria do not make it possible however to adapt the step of time. They are
indeed delicate to implement in the nonlinear cases and do not take account of the variations
excitation.
Another approach consists in studying the site error introduced by the diagram using
limited developments.
The exact solution of a system to a degree of freedom checks:
T
T
T
2
T
3
X T +
X (T)
X & (T)
X & (T)
X
& (&t)
(O T3)
2 =
+
+
+
+
2
8
48
T
T
T
2
T3
X T -
X (T)
X & (T)
X & (T)
X
& (&t)
(O T3)
2 =
-
+
-
+
2
8
48
T
T
T
3
X T +
X T
T
X & (T)
X
& (&t)
(O T3)
2 =
-
2 +
+
+
24
The formula of integration of the differences thus leads to a truncation error being worth:
T
3
T
2
E =
X
& (&t)
(X & - X &
N
N
N
n-1)
24
12
One can normalize this error to obtain a relative error:
T
2 X & - X&
E
N
N
=
- 1
X 0
N
12
X
N
N
Park and Underwood [bib7] interpreted this error by defining a “apparent pulsation”:
X&
2
N
=
With N
X N
Handbook of Référence
R5.06 booklet: Dynamics in modal base
HT-66/03/005/A
Code_Aster ®
Version
6.4
Titrate:
Algorithms of temporal integration of operator DYNA_TRAN_MODAL
Date:
28/02/03
Author (S):
E. BOYERE, LIGHT A.C., G. JACQUART Key
:
R5.06.04-B Page
: 14/18
Applied to the diagram with centered difference, this definition makes it possible to interpret the relative error E
N
like a variation of the apparent pulsation:
T
2
E
2
-
2
N
With N
With n-1
12
Many algorithms use a criterion of adaptation of the step of time based on the error of
truncation ([bib9], [bib11]). However in the case of a conditionally stable diagram, this
method neither to ensure itself of the stability of integration, nor to guarantee a precision for the calculation of
transients.
Other methods use an approximation of the instantaneous own pulsation of the system [bib12],
using the matrices of mass and stiffness. They have the defect not to adapt to the forces
external and with their fluctuations in frequency.
It is thus useful to find a criterion which takes account of the two approaches. This is why Park and
Underwood introduced the concept of “frequency connect disturbed”:
1
X & - X&
F
N
N
=
- 1
AP N
2
X - X
N
n-1
This size is interpreted like the “instantaneous” frequency of the system.
In the case of a system with several degrees of freedom, it is necessary to calculate a frequency connect for
each degree of freedom and to take the maximum. The step of time can be then selected to respect
a minimum of points per apparent period.
If the denominator of the expression of the apparent frequency tends towards zero, this one can become very
large and not to have significance more. This leads to an unjustified refinement when speed
cancel yourself. To cure it a criterion of the type is added:
X - X -
&
&
1
1
-
N
N
X
X
< X&
F
N
N
min
=
- 1
T
AP N
2
X&
T
min
It is an intermediary between the disturbed apparent frequency and the truncation error. The value
adequate of &
X
is difficult to choose a priori and an unsuited value involves a reduction
min
artificial of the apparent frequency. In the case of a system with several degrees of freedom one
circumvent this difficulty by employing the “close” degrees of freedom:
1 X I
I
& N - X&
F
N
=
- 1
max
AP N
1 inb DLL 2
Bi
N
Bi
J
J
N = max max X N - X n-1, X
&
T
min
I J lb
-
Handbook of Référence
R5.06 booklet: Dynamics in modal base
HT-66/03/005/A
Code_Aster ®
Version
6.4
Titrate:
Algorithms of temporal integration of operator DYNA_TRAN_MODAL
Date:
28/02/03
Author (S):
E. BOYERE, LIGHT A.C., G. JACQUART Key
:
R5.06.04-B Page
: 15/18
where lb indicates for example the bandwidth of the matrix of stiffness and where &
X
can be selected
min
very small.
This method appears very effective if X I N indicate physical components
(displacements). In the case of a projection on modal basis it is not relevant to employ them
close components to calculate the apparent frequency. In this case it is to better return to
first criterion and to use one of the two following methods, specified by key word VITE_MIN:
&X
·
N
if VITE_MIN: “NORM” then it is a variable parameter equal to
100
(&
X
I
=
X 2). This method gives good results when the number of degrees
N
& N
1 inb ddl
of freedom is large and is inapplicable with the case with only one degree of freedom. It is not any more
indicated if the order of magnitude speed is very different from a degree of freedom to another.
· if VITE_MIN: “NORM” then it is a parameter variable and different for each degree from
&X jm
freedom, &
X
J = max
. This method has the advantage of functioning whatever it
min
1 mn 1001
degree of freedom of the system numbers but it cannot be used if the command of
size speed varies too much during calculation because, in this case, one would obtain
X J
J
N - X n-1
systematically:
X J
&
.
T
min
2.3.4.6 Algorithm of the diagram of the centered differences with adaptive step
The rules mentioned above make it possible to fix a number of steps of time desired per period
response according to the wanted precision, NR. It is adjustable by the key word
1
NB_POINT_PERIODE. The step of time T must then be lower than
. The key word
N
NR F
PAS gives
AP N
at the same time the step of initial time, T, and the step of maximum time not to exceed, T
.
ini
max
Handbook of Référence
R5.06 booklet: Dynamics in modal base
HT-66/03/005/A
Code_Aster ®
Version
6.4
Titrate:
Algorithms of temporal integration of operator DYNA_TRAN_MODAL
Date:
28/02/03
Author (S):
E. BOYERE, LIGHT A.C., G. JACQUART Key
:
R5.06.04-B Page
: 16/18
The algorithm is described schematically below:
0) initialization: X and &
X given
0
0
T
0, T =
and &X = & (
X T, X, &X
0
0
0
0)
0
T
1 =
-
ini
initialization of &
X
min
with each step of time
1) initialization of the search of the step of time: NR
= 0
iter
2) calculation of &
X
then X
:
1
n+1
n+ 2
T
estimate speed (method 2): &
X = &X
1 +
&X
N
n-1
N
2
2
T
1 +
T
speed with the semi step: &
X
1 = &
X
N
N
1 +
&X T
N (
, X,
N
N
&Xn)
n+
N
2
2
2
displacement: X
&
1 = X
+ T X
n+
N
1
n+ 2
3) calculation of acceleration &X
n+1
4) calculation of the apparent frequency:
X - X -
1
&
&
1
-
N
N
X
X
X&
F
N
N
min
=
- 1
T
APn
2
X - X
N
n-1
X - X -
1
&
&
1
-
N
N
X
X
< X&
F
N
N
min
=
- 1
T
APn
2
X&
T
min
5) checking of the adequacy enters the step of time and the apparent frequency:
calculation of the indicator err = T
NR F
N
AP N
· if err 1 and NR
< NR
then reduction of the step of time and new iteration of search
itez
iter max
step: T 0 7
, 5t, NR
NR + 1, return in 2)
N
N
iter
iter
· if err 0 7
, 5 since more than 5 steps of consecutive times, then increase in the step of time
T 11
, T
N
N
6) filing of the solution, possible calculation of &
X
and return in 1) for the following iteration.
min
2.3.4.7 Comments on the parameters of the algorithm
The fact of fixing an upper limit NR
by the key word
iter
NMAX_ITER_PAS with the number of
max
reductions of the step of time makes it possible to be ensured of the convergence of the algorithm in the cases
difficult (for example in the event of discontinuity in the external forces).
When the indicator err is higher than 1, the step of time is multiplied by a factor fixes (0,75 by
defect but it can be modified by the user thanks to operand COEF_DIVI_PAS). It would have been
1
possible to write directly: T
T = NR F
, which more intuitive. But this strategy
N
err
N
AP N
conduit with an excessive refinement, the calculated apparent frequency being often largely higher
at the real frequency, when the error is large. However, in only one step of time, T can be
N
considerably reduced (factor 0 7
, 5Niter).
Handbook of Référence
R5.06 booklet: Dynamics in modal base
HT-66/03/005/A
Code_Aster ®
Version
6.4
Titrate:
Algorithms of temporal integration of operator DYNA_TRAN_MODAL
Date:
28/02/03
Author (S):
E. BOYERE, LIGHT A.C., G. JACQUART Key
:
R5.06.04-B Page
: 17/18
On the other hand the increase in the step of time is always much slower (coefficient of
multiplication owing to lack of 1,1 definable by COEF_MULT_PAS) and takes place only if the indicator is
lower than 1 during five steps of consecutive time. These restrictions are justified by the risks of
loss of stability or precision of the diagram when the step of time varies too quickly. A coefficient of
1, 2 or 1, 3 can allow a faster calculation but exposes sometimes at the risks of error.
In short the default values were validated by many tests and give in general
satisfaction in terms of precision and stability [bib8].
2.3.4.8 Performance of the algorithm
With equal precision, the iteration count carried out by diagram ADAPT is at least five times more
weak that with a constant step in the phenomena which justify the use of a variable step by
the irregular aspect of their evolution (shocks, fluid blade, discontinuous excitations, etc).
The empirical studies showed that the step of adaptive time allows in the successful outcomes of
to gain a factor two or three in calculating times. This diagram makes it possible moreover to control
precision of integration by the method of the control of the number of points per period of the answer.
In the case of the very deadened systems, the gains can be even more important (calculations five with
ten times faster).
On the other hand when the step of “ideal” time is about constant, the use of diagram ADAPT
appears useless.
It of course allows the taking into account of nonthe localized linearities of shocks type or frictions, thus
that fluid blades.
In dynamic under-structuring, it is compatible as well with the transitory analysis on the basis
modal restored on the whole system or transitory calculation on the bases distinct from under
structures.
3 Conclusion
As a conclusion, summarized here various possibilities of temporal integration which offers
the operator:
· Euler (“EULER”) clarifies modified to ensure a conditional stability,
· Diagram of Newmark (“NEWMARK”) parameterized in order to be implicit,
· Diagram of Devogelaere-Fu (“DEVOGE”) of command 4,
· Explicit adaptive diagram (“ADAPT”).
The diagram by defect is EULER but it is not systematically adapted more. The diagram of
NEWMARK available in Code_Aster is implicit and guarantees an unconditional stability but
is valid only for purely linear problems. Diagram DEVOGE is of command 4 and thus is
more precise but it is is expensive in calculating times. Diagram ADPAT is more particularly
indicated for the problems with nonlocated linearities, where the step of “ideal” time is not
constant during the transient. It is thus the experiment of the modeling which makes it possible to choose it
diagram adapted best to the problem according to the report/ratio (calculating time)/precision.
Handbook of Référence
R5.06 booklet: Dynamics in modal base
HT-66/03/005/A
Code_Aster ®
Version
6.4
Titrate:
Algorithms of temporal integration of operator DYNA_TRAN_MODAL
Date:
28/02/03
Author (S):
E. BOYERE, LIGHT A.C., G. JACQUART Key
:
R5.06.04-B Page
: 18/18
4 Bibliography
[1]
K. - J. BATHE, E. - L. WILSON: “Numerical Methods in Finite Element Analysis”, Prentice
Hall Inc.
[2]
J. ANTUNE, F. AXISA, H. BUNG, F. DOVEIL, E. LANGRE: “Methods of analysis in
nonlinear dynamics of the structure “- IPSI
[3]
J. - R. LEVESQUE, P. LABBE and Al: “Module TRANSIS of code POUX in documentation of
reference of code POUX “- internal Rapport EDF
[4]
X. RAUD, P. RICHARD: “Structure and validation of the nonlinear calculation of the bearings
hydrodynamics of code CADYRO “ Note HP-61/92/041
[5]
G. - D. HAHN: “A Modified Euler Method for Dynamic Analysis” - International Journal for
Numerical Methods in Engineering Vol. 32 (1991), p 943-955
[6]
K. - C. PARK, P.-G. UNDERWOOD: “A Variable-Step Central Difference Method for
Structural Dynamic Analysis Part II: Implementation and performance evaluation “-
Methods computer in Applied Mechanics and Engineering Vol. 23 (1980) p 259-279
[7]
K. - C. PARK, P.-G. UNDERWOOD: “A Variable-Step Central Difference Method for
Structural Dynamic Analysis Part I: Theoritical Aspects “- Computer Methods in Applied
Mechanics and Engineering Vol. 22 (1980) p 241-258
[8]
G. JACQUART, S. GARREAU: “Algorithm of integration to step of adaptive time in
Code_Aster “- Note HP-61/95/023
[9]
O.C. ZIENKIEWICZ, Y.M XIE: “A posteriori Local Error Estimation and Adaptative
Time-Stepping Procedure for Dynamic Analysis “Earthquake Engineering and Structural
Dynamics Vol. 20 (1991) pp. 871-887
[10]
A.-C. LIGHT, G. JACQUART: “Algorithms of integration of operator DYNA_TRAN_MODAL
of Code_Aster: reference material “- Note HP-51/96/072
[11]
R. Mr. THOMAS, I. GLADWELL: “Variable Order Variable Step Algorithms for Second Order
systems “- International Journal for Numerical Methods in Engineering Vol. 26 (1998)
pp. 39-53
[12]
P.G. BERGAN, E. MOLLESTAD: “Year Automatic Time-Stepping Algorithm for Dynamic
Problem “- Computer Methods in Applied Mechanics and Engineering Vol. 49 (1985)
pp. 299-318
[13]
NR. GAY, S. GRANGER, T. FRIOU: “Presentation of a method of simulation of the coupling
fluidelastic in nonlinear mode “- Note HT-32/94/015
Handbook of Référence
R5.06 booklet: Dynamics in modal base
HT-66/03/005/A
Outline document