Code_Aster ®
Version
6.4
Titrate:
Algorithm of integration of DYNA_LINE_TRAN
Date:
28/02/03
Author (S):
E. BOYERE, G. JACQUART, LIGHT A.C.,
Key:
R5.05.02-C Page
: 1/20
Organization (S): EDF-R & D/AMA, EDF-Pôle Industrie/CNPE of Tricastin, EDF-R & D/TESE
Handbook of Référence
R5.05 booklet: Transitory or harmonic dynamics
Document: R5.05.02
Algorithms of direct integration of the operator
DYNA_LINE_TRAN
Summary:
This document describes the diagrams of temporal integration which are used to solve in a direct way of
problems of dynamics in transitory linear mechanics. The diagrams of NEWMARK and WILSON are
detailed, as well as the diagrams “centered differences with constant step” and “not of adaptive time”.
Handbook of Référence
R5.05 booklet: Transitory or harmonic dynamics
HT-66/03/005/A
Code_Aster ®
Version
6.4
Titrate:
Algorithm of integration of DYNA_LINE_TRAN
Date:
28/02/03
Author (S):
E. BOYERE, G. JACQUART, LIGHT A.C.,
Key:
R5.05.02-C Page
: 2/20
Count
matters
1 Introduction ............................................................................................................................................ 3
2 Methods of implicit direct temporal integration and clarifies of a dynamic problem .................. 4
3 the diagram WILSON [bib1] ................................................................................................................. 6
3.1 Presentation of the diagram .................................................................................................................. 6
3.2 Complete algorithm of the WILSON method:.............................................................................. 7
3.3 Stability condition of the diagram WILSON .................................................................................. 7
4 the diagram of NEWMARK [bib1], [bib2] ................................................................................................ 8
4.1 Presentation of the diagram .................................................................................................................. 8
4.2 Complete algorithm of the method of NEWMARK ......................................................................... 8
4.3 Stability conditions of the diagram of NEWMARK:........................................................................ 9
5 numerical Damping of the implicit schemes ............................................................................ 10
6 Diagram of the centered differences with constant step ............................................................................... 12
6.1 Principle .......................................................................................................................................... 12
6.2 Stability .......................................................................................................................................... 13
6.3 Algorithm ...................................................................................................................................... 13
6.4 Stamp of diagonal mass .......................................................................................................... 13
6.5 Checking of the step of time ......................................................................................................... 14
6.6 Calculation of acceleration .................................................................................................................. 14
7 Diagram with step of adaptive time ....................................................................................................... 14
7.1 Principle .......................................................................................................................................... 14
7.2 Diagram .......................................................................................................................................... 15
7.3 Estimate of the step of time according to the precision required ................................................... 16
7.3.1 influences of the close nodes .............................................................................................. 16
7.3.2 use of information at previous time .................................................................... 17
7.4 Choice of the number of steps per apparent period, NR ...................................................................... 17
7.5 Heuristics of evolution of the step of time T ............................................................................... 17
N
7.6 Algorithm ...................................................................................................................................... 18
8 Conclusion ........................................................................................................................................... 19
9 Bibliography ........................................................................................................................................ 20
Handbook of Référence
R5.05 booklet: Transitory or harmonic dynamics
HT-66/03/005/A
Code_Aster ®
Version
6.4
Titrate:
Algorithm of integration of DYNA_LINE_TRAN
Date:
28/02/03
Author (S):
E. BOYERE, G. JACQUART, LIGHT A.C.,
Key:
R5.05.02-C Page
: 3/20
1 Introduction
The goal of the transitory dynamic analysis is to determine according to time the response of one
structure, being given an external loading, or boundary conditions functions of time
when the effects of inertia cannot be neglected.
In a certain number of physical configurations, one cannot do without a transitory analysis
while being satisfied with a modal or harmonic analysis:
· if the history of the phenomenon has an importance in the study,
· if the external loading is complex (seism, excitations multi-components…),
· if the system is nonlinear (plasticity, shocks, friction…).
The methods of transitory analysis dynamic which can then be used are divided into two
main categories:
· methods known as of direct integration,
· methods of RITZ including/understanding inter alia the recombination of modal projection.
The methods of direct integration are thus called because of the fact that no transformation
is not carried out on the dynamic system after the discretization in finite elements.
We will make a presentation of the algorithms of direct integration used to solve one
dynamic problem in mechanics for linear structures. These algorithms are employed
in operator DYNA_LINE_TRAN of Code_Aster.
The methods of RITZ, on the contrary, proceed to a transformation of the initial dynamic system,
very often a projection on a subspace of the space of starting discretization. The resolution
dynamics is done then on a modified system, which gives access only one approximation of
response of the initial system. They are presented in another document [R5.06.01].
Handbook of Référence
R5.05 booklet: Transitory or harmonic dynamics
HT-66/03/005/A
Code_Aster ®
Version
6.4
Titrate:
Algorithm of integration of DYNA_LINE_TRAN
Date:
28/02/03
Author (S):
E. BOYERE, G. JACQUART, LIGHT A.C.,
Key:
R5.05.02-C Page
: 4/20
2 Methods of implicit direct temporal integration and
clarify of a dynamic problem
It is supposed that the studied structure has a linear behavior and that equations governing sound
balance dynamic were discretized by finished differences or finite elements. One is obtained
discrete system of differential equations of the second command which it is a question of integrating in time.
In a general way, these equations take the following form:
Mr. X
& + C.X +
T
&
K.X = R
T
T
T
·
M is the matrix of mass of the system,
·
C is the matrix of viscous damping of the system,
·
K is the elastic matrix of rigidity of the system,
·
R is the vector of the external forces applied to the viscous system.
The system is of the second command.
Two classes of methods of integration can be distinguished to integrate them step by step
equilibrium equations: they are the methods of explicit and implicit integration.
Let us see what distinguishes them by examining temporal integration from the following linear system:
Mr. X
& + C.X +
T
&
K.X = R
T
T
T
This differential connection of the second command can be brought back to a first order system:
A.U & = B.u + F
éq 2-1
where:
X
X&
I 0
0
I
0
U = U = A =
B =
F =
X
&
X
&
0 M
- K - C
R
To integrate this differential equation, one thus uses a discretization Ti of the interval of study
that a formula of differences finished to express the derivative &u.
One calls methods of integration clarifies the methods where, in [éq 2-1] written at time Ti, only
the derivative &u utilizes the variable U at time ti+1. In this way, determination of
sizes sought at the moment ti+1 does not result from an inversion of system utilizing
the operator K. If moreover, one carries out a “farmhouse-lumping” in order to return the matrix M diagonal,
determination of ui+1 is particularly simple. They are the principal characteristics there of
methods of explicit integration.
The implicit or semi-implicit methods utilize the discretization of U in [éq 2-1] to one
posterior moment with Ti, generally ti+1, in order to determine the variables of the problem with ti+1.
Their determination thus passes by the resolution of a system utilizing the operator K.
Two concepts concerning the diagrams of integration are important: consistency and stability.
Handbook of Référence
R5.05 booklet: Transitory or harmonic dynamics
HT-66/03/005/A
Code_Aster ®
Version
6.4
Titrate:
Algorithm of integration of DYNA_LINE_TRAN
Date:
28/02/03
Author (S):
E. BOYERE, G. JACQUART, LIGHT A.C.,
Key:
R5.05.02-C Page
: 5/20
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 derived 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 (cf [bib2]).
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 the diagram.
Certain implicit algorithms have the characteristic to be unconditionally stable, according to the choice
certain parameters, which makes their interest and makes it possible to integrate the dynamic phenomenon with one
no arbitrarily large time.
The diagram of WILSON and the diagram of NEWMARK can be explicit for certain choices of
their parameters. In Code_Aster, they are employed for their properties of stability
unconditional, clean with the implicit schemes. They will thus be classified here in the category of
implicit schemes and one will see under which conditions they give the properties of stability
wished.
Two explicit diagrams of integration were also introduced into Code_Aster. It is about
diagrams DIFF_CENTER and ADAPT which are both based on the method of the centered differences.
They are conditionally stable and requires to be powerful a matrix of mass
diagonalized. Conditional stability leads to a control of the step of time which, exploited in the case
diagram ADAPT, allows an adaptation of the step of time according to the speed of
modelled phenomena.
Handbook of Référence
R5.05 booklet: Transitory or harmonic dynamics
HT-66/03/005/A
Code_Aster ®
Version
6.4
Titrate:
Algorithm of integration of DYNA_LINE_TRAN
Date:
28/02/03
Author (S):
E. BOYERE, G. JACQUART, LIGHT A.C.,
Key:
R5.05.02-C Page
: 6/20
3
The diagram WILSON [bib1]
3.1
Presentation of the diagram
It will be supposed in what follows that the solids are elastic linear. This method leaves
the assumption that acceleration is linear between T and T +. T
:
&X
= &X +
.
+
(&X
-
+.
&X
T
T
T
T
T)
. T
éq
3.1-1
&
&
X
X
T +.t
T
T
+
&Xt
T
T + T T +. T
While integrating [éq 3.1-1] according to the variable, one obtains:
2
&X
= &X + &X +
.
+
(&X
-
+.
&X
T
T
T
T
T
T)
2. T
éq
3.1-2
2
3
X = X + X & +
X
+
& +
(X
.
&
.
- X
T
T
T
t+
T
& T)
2
6 T
éq
3.1-3
One writes the equilibrium equations at time T +. T
with 1:
Mr. X
&
+ C.X
t+ T
& +
K.X =
R
t+ T
t+ T
t+ T
éq
3.1-4
while expressing &
Xt+. T and &Xt+. T according to &Xt+. T and of Xt, &Xt and &Xt by the system
[éq 3.1-2], [éq 3.1-3], and while replacing in [éq 3.1-4], it comes:
~
~
K. X
=
+.
R
T
T
where
~
3
6
K = K + (
. C +
. M
.t)
(.t) 2
~R = R t+.(R
t+ T
- R T) + Mr. (A.X
0
T + A.
2 &
X
T + 2. &Xt) + C. (A.X
1
T + 2. &
X
T + A.
3 &X
T)
6
3
. T
a0 = (
has =
=
has
2.a
=
has
. T
) 2
1
(. T
)
2
1
3
2
Handbook of Référence
R5.05 booklet: Transitory or harmonic dynamics
HT-66/03/005/A
Code_Aster ®
Version
6.4
Titrate:
Algorithm of integration of DYNA_LINE_TRAN
Date:
28/02/03
Author (S):
E. BOYERE, G. JACQUART, LIGHT A.C.,
Key:
R5.05.02-C Page
: 7/20
One goes back with displacements, speeds and accelerations to the step T + T
by the relations:
&X =
a.
t+ T
4 (X - X
+
t+. T
T)
a.
5 &
X + A.
T
6 &Xt
&X
=
t+ T
&X + A.
T
7 (&X
+
t+ T
&Xt)
X
= X + T
.
t+ T
T
&X + A.
T
8 (&X
+ 2.
t+ T
&Xt)
has
- has
3
T
t2
has
0
=
has
2
=
= 1 has -
has =
has
4
5
6
7
8 =
2
6
3.2
Complete algorithm of the WILSON method:
) initialization has:
1) initial conditions 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 rigidity K and mass M
~
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 +.(Rt+ T - Rt) + Mr. (A. Xt + A.X&t + 2.X&t) + C. (A. Xt + 2.X&t + A.X&t
0
2
1
3
)
~
~
2) to solve K. X
=
+.
R
T
T
3) to calculate displacements at time T + T
&X =
a.
t+ T
4 (X - X
+
t+. T
T)
a.
5 &
X + A.
T
6 &Xt
&X
=
t+ T
&X + A.
T
7 (&X
+
t+ T
&Xt)
X
= X + T
.
t+ T
T
&X + A.
T
8 (&X
+ 2.
t+ T
&Xt)
4) calculation of the step of next time: return to the beginning
3.3
Stability condition of the WILSON diagram
The method is unconditionally stable for WILSON > 1.37, a value usually employed
for being 1.4. Moreover, the method presents numerical dissipation for > 1, all the more
important what increases.
The key word WILSON factor: (THETA: HT) makes it possible to specify the use of this algorithm and the choice
value of. By defect, the value of is taken to 1.4.
Handbook of Référence
R5.05 booklet: Transitory or harmonic dynamics
HT-66/03/005/A
Code_Aster ®
Version
6.4
Titrate:
Algorithm of integration of DYNA_LINE_TRAN
Date:
28/02/03
Author (S):
E. BOYERE, G. JACQUART, LIGHT A.C.,
Key:
R5.05.02-C Page
: 8/20
4
The diagram of NEWMARK [bib1], [bib2]
4.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
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 =
- 1
has =
- 2
has
T
. 1
.
4
5
6
(
) has
T
7
2
=
-
=
4.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 and mass M
~
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:
~
to calculate the effective loading R
2
~
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}
~
~
to solve K. X
=
+.
R
T
T
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
5
calculation of the step of next time: return to the beginning
Handbook of Référence
R5.05 booklet: Transitory or harmonic dynamics
HT-66/03/005/A
Code_Aster ®
Version
6.4
Titrate:
Algorithm of integration of DYNA_LINE_TRAN
Date:
28/02/03
Author (S):
E. BOYERE, G. JACQUART, LIGHT A.C.,
Key:
R5.05.02-C Page
: 9/20
4.3
Stability conditions of the diagram of NEWMARK:
The method of NEWMARK is unconditionally stable if:
(
2 +) 2
1
>.
0 5 and >
4
One introduces a positive numerical damping if >.
0 5 and negative if <.
0 5. When
=.
0 5 and = 0, the formula of NEWMARK are reduced to the diagram centered differences. One
1
combination very often employed is =.
0 5 and =, because it leads to a diagram of command
4
2, unconditionally stable without numerical damping.
This diagram of integration is 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. It is integrated in Code_Aster in operator DYNA_LINE_TRAN.
key word factor NEWMARK: (ALPHA: Al, DELTA: of) allows to specify the employment of this
algorithm and the choice of the value of and. By defect, the value of is taken to 0.25 and the value
of is taken to 0.5.
Handbook of Référence
R5.05 booklet: Transitory or harmonic dynamics
HT-66/03/005/A
Code_Aster ®
Version
6.4
Titrate:
Algorithm of integration of DYNA_LINE_TRAN
Date:
28/02/03
Author (S):
E. BOYERE, G. JACQUART, LIGHT A.C.,
Key:
R5.05.02-C Page
: 10/20
5
Numerical damping 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 substanciellement large compared to the smallest clean period of the system without
to be likely to cause an instability of the results.
However, if the contents of the answer reside in a whole of clean modes, of which highest
an Eigen frequency has Fmax, one will have to still respect a criterion on the step of time of the form:
1
1
T < (
with
10 * F
) (100 * F
max
max)
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
these high modes.
One can see on the graph hereafter the reduction in amplitude of a system to a degree of freedom, without
damping, when one integrates it by various methods (WILSON and NEWMARK
1
1
=,
=):
2
4
NEWMARK method
1
1
Percentage amplitude decay (AD X 100%)
=, =
2
4
It is checked here that the algorithm of NEWMARK with these parameters does not present any damping
numerical.
Handbook of Référence
R5.05 booklet: Transitory or harmonic dynamics
HT-66/03/005/A
Code_Aster ®
Version
6.4
Titrate:
Algorithm of integration of DYNA_LINE_TRAN
Date:
28/02/03
Author (S):
E. BOYERE, G. JACQUART, LIGHT A.C.,
Key:
R5.05.02-C Page
: 11/20
On the other hand, the implicit algorithms also have a rather significant effect of elongation of the periods
clean contained in the response of the structure which leads to a dephasing of the calculated solution.
The graph below presents percentages of elongation of the clean period of a system at one
ddl without damping.
Percentage period elongation (EP/T X 100%)
T/T
On these 2 graphs, one notes that to guarantee a precision on the amplitude and the phase of
calculated displacements, it is necessary to respect a criterion close to:
.
01
.
0 01
T <
with
F
F
max
max
where Fmax is the high frequency of the movement which one wishes to correctly capture in
numerical analysis.
Handbook of Référence
R5.05 booklet: Transitory or harmonic dynamics
HT-66/03/005/A
Code_Aster ®
Version
6.4
Titrate:
Algorithm of integration of DYNA_LINE_TRAN
Date:
28/02/03
Author (S):
E. BOYERE, G. JACQUART, LIGHT A.C.,
Key:
R5.05.02-C Page
: 12/20
6
Diagram of the centered differences with constant step
6.1 Principle
The diagram clarifies centered differences with constant step is written:
&X
2
1 = &
X 1 + T &X T
+ O T
N (
, X,
N
N
&Xn) (
)
n+
N
2
2
X
2
1 = X
+ T
T
1
+ O T
n+
N
&X (, X,
N
N
&Xn) (
)
n+ 2
with the following notations:
X
X
N
&x
- 1
X
&x
n-12
N
n+12
n+1
tn-1
tn-1/2
tn
tn+1/2
tn+1
T
T
.
The speed is expressed with indices half-entireties of the discretization in time whereas them
displacements and accelerations are expressed with the whole indices.
Written this way, the diagram is of command 2.
Acceleration in T is not immediately calculable because speed is known only with the half-not of
N
previous time (in T
), which poses problem to evaluate the terms of damping. For
N 12
to circumvent this difficulty, one calculates acceleration in T per L `following approximation:
N
&X T
T
-
T
1
-
-
N (
, X,
N
N
&Xn) &X
, X,
N
N
&X
M 1 F
N
1
2
=
(N)
K X
C
N
&X
N
N 2
-
what constitutes a valid approximation if damping is sufficiently weak
(&
X =
+ O 1). The diagram loses its precision of command 2 if the damping of the structure is
N
&X
()
N 12
important.
Other methods of approximation of acceleration can be considered. That selected is
revealed a good compromise enters simplicity and stability, like the study described in the reference
[bib4] on the precision and the stability of several methods.
The fields are filed at the moment T, T,
N
;+1 L, speed being approximate at these moments by
following formula:
&
T
X
= &X
+
&X T, X,
1
1
1
&X
n+
N
N
N
N
+
+
+
+
n+ 1
2
2
1
12
Handbook of Référence
R5.05 booklet: Transitory or harmonic dynamics
HT-66/03/005/A
Code_Aster ®
Version
6.4
Titrate:
Algorithm of integration of DYNA_LINE_TRAN
Date:
28/02/03
Author (S):
E. BOYERE, G. JACQUART, LIGHT A.C.,
Key:
R5.05.02-C Page
: 13/20
6.2 Stability
The diagram of the centered differences is conditionally stable. In the case of a system without
damping [bib2], the diagram is stable for a step of time checking T < 2/
where
is
max
max
the greatest own pulsation of the system, is T < T
/. One needs a minimum of step of time
min
to describe the smallest period of the system T
.
min
The limiting value for the step of time decreases slowly when damping increases [bib4]. By
example, for a damping of 0,5%, the condition becomes T < T
/5.
min
6.3 Algorithm
In short, the diagram such as it is introduced into Code_Aster arises in the following way:
0 inialisation
:
T, X, &X given
0
0
&X = M 1
T = 0 -
-
0
(F (
)) K X C
0
&X0
T
&X
1
= &X0 -
&X
-
0
2
2
1
With each step of time X, X
&
, X
& known
N
N
N
- 12
&X
1
= &X 1 + T &X T, X,
N
&X
n+
N
N
N
+
1
2
2
n+ 2
X
1 = X
+ T
n+
N
&Xn+12
&X
= M-1 F T -
1 -
n+1
(N) KX
C
N
&X
+
N 12
T
&X 1 =
1
+
N
&X
&X
+
N
N
+
+1
2
2
2
possible filing of X
, X
&
, X
&
n+1
n+1
n+1
then return at stage 1) for the following step.
6.4
Stamp of diagonal mass
The calculation of acceleration requires the inversion of the matrix of mass. This explicit diagram becomes
more powerful if a matrix of concentrated mass is used (`MASS_LUMPING `) so that it
that is to say diagonal. The inversion then does not require any more factorization and is immediate.
This is why in Code_Aster, the diagram of centered differences is licit only with matrices
of mass built in a diagonal way, by option “MASS_MECA_DIAG” of the operator
CALC_MATR_ELEM.
Handbook of Référence
R5.05 booklet: Transitory or harmonic dynamics
HT-66/03/005/A
Code_Aster ®
Version
6.4
Titrate:
Algorithm of integration of DYNA_LINE_TRAN
Date:
28/02/03
Author (S):
E. BOYERE, G. JACQUART, LIGHT A.C.,
Key:
R5.05.02-C Page
: 14/20
6.5
Checking of the step of time
It was seen that the diagram of centered differences is stable provided that the step of time, in
the absence of damping, is lower than a limiting value, equal to T < 2/
. In practice
max
one employs a step of time which is worth from 5% to 20% of the step of critical time. It was thus introduced one
test on the step of time which checks that:
2
T < 0 0
, 5
where K and m are the diagonal terms of the matrices of stiffness and of
K
II
II
II
max
1 inddl mii
mass.
If this condition is not checked, the user is stopped with a message indicating the step to him of
maximum time which it can use.
6.6
Calculation of acceleration
The calculation of acceleration is done as follows:
for each degree of freedom, one tests if the diagonal term of the matrix of mass corresponding is
no one.
· if it is not null, the term of acceleration is calculated according to the formula
:
&X
= M-1 F T
N
-
N 1 -
+
() KX
C &X
N 1
+
N 12
· if it is null, the term of acceleration is not calculated. It is the case for degrees of freedom
said of Lagrange. If they correspond to blocked degrees of freedom, it is licit of not
to take account of the line in question and not to calculate its acceleration. If it
degree of freedom of Lagrange was introduced to define a connection between two freedom degrees,
that does not have any more a smell. The diagram is thus then unusable and a test stops the execution with
an explicit message.
7
Diagram with step of adaptive time
7.1 Principle
The methods of calculation clarifies are particularly indicated in simulation of phenomena
rapids, such as the wave propagation in the solids. On the other hand, they agree less
well with slower phenomena since the stability condition of the diagram imposes a step of
time of about a smallest clean period of the system.
The adaptive diagram, based on the diagram of centered differences, was developed to allow it
calculation of transitory answers in which fast” and “slow” phenomena “. For example
at the time of an impact, initially of the high frequency waves are propagated and dissipate themselves
in the structure. Then, the structure does not answer any more but on its modes low frequency, them
high frequencies being deadened. The idea is thus to adapt the step of time progressively in
function of the phenomena brought into play, by fixing a criterion of precision on the solution.
Handbook of Référence
R5.05 booklet: Transitory or harmonic dynamics
HT-66/03/005/A
Code_Aster ®
Version
6.4
Titrate:
Algorithm of integration of DYNA_LINE_TRAN
Date:
28/02/03
Author (S):
E. BOYERE, G. JACQUART, LIGHT A.C.,
Key:
R5.05.02-C Page
: 15/20
7.2 Diagram
The diagram clarifies centered differences with variable steps is written:
T
1 +
T
&X
2
1 = &
X
N
N
1 +
&X T
+ O T
N (
, X,
N
N
&Xn) (
)
n+
N
2
2
2
X
2
1 = X
+ T
T
1
+ O T
n+
N
&X (, X,
N
N
&Xn) (
)
n+ 2
with the following notations:
T1 +
-
T
N
N
2
X
X
&x
N
&x
- 1
n-1
N
N
+1
X
2
2
N + 1
T
T
n+1/2
tn+1
n-1
tn-1/2
tn
T
T
n-1
N
.
It is noted that the step of time varies. It is subscripted: T.
N
That has as a consequence which the diagram is not rigorously any more of the second command, since it is not
“centered more”. More T
and T are different, plus the command of the diagram is close to 1. The strong ones
n-1
N
variations of the step of time are thus accompanied by a fall of precision. The formula speed
employee leads to good results when the step of time decreases but cause a drop in the limit of
stability when the step of time increases. This is why one it constrained to increase only very
gradually.
Lastly, one uses the same approximations as for the differences centered with regard to
calculation of accelerations and speeds to the steps of “whole” times:
· acceleration is estimated by &X T
T
and
N (
, X,
N
N
&Xn) &X
, X,
N
N
&X
N
N
- 12
&X T, X,
-
T -
-
;
N
N
&X
M 1 F
N
=
(N)
K X
C
N
&X
N
N 12
12
T
· and stored speed is evaluated by &X
.
1 = &
X 1 +
&X
n+
N
N
+
+1
2
2
As for the diagram of centered differences, of which it is inspired, the diagram with adaptive step
require the inversion of the matrix of mass. This is why the diagonalisation of the matrix is required
of mass as well as the same restrictions on the degrees of freedom of Lagrange as for the diagram
with centered differences.
Handbook of Référence
R5.05 booklet: Transitory or harmonic dynamics
HT-66/03/005/A
Code_Aster ®
Version
6.4
Titrate:
Algorithm of integration of DYNA_LINE_TRAN
Date:
28/02/03
Author (S):
E. BOYERE, G. JACQUART, LIGHT A.C.,
Key:
R5.05.02-C Page
: 16/20
7.3
Estimate of the step of time according to the precision required
To define a criterion on the step of time according to the precision required on the solution, one
introduced the concept of frequency connect disturbed [bib4]:
1
X & - X&
F
X
X
=
- 1
AP N
2
X - X
X
x-1
This size can be interpreted like the “instantaneous frequency” of the system. It is indeed one
approximation of the local slope of the curve forces/displacement. It is related to the error due to
truncation in the limited developments of the diagram. It makes it possible moreover to take account of
forces external and of their fluctuations in frequency.
In the case of a system with several degrees of freedom, it is necessary to calculate a frequency
connect for each degree of freedom. One then employs the maximum on all the frequencies
calculated to determine the step of time.
If the denominator tends towards zero, the apparent frequency can become very large and lose its
physical significance. One then obtains an unjustified refinement of the step of time when speed
cancel yourself. In the case of sinusoidal oscillations, it is the case twice per period. It then is modified
criterion by introducing the following condition:
X - X -
&
&
1
1
-
X
X
X
X
X&
F
X
X
min
=
- 1
T
AP N
2 p
X&
T
min
One obtains an intermediary between the frequency connect disturbed and the truncation error. The value
of &
X
is not easy to determine a priori and a badly selected value can lead to one
min
artificial moderation of the apparent frequency.
Two methods are proposed.
7.3.1 influences of the close nodes
In the case of a system with several degrees of freedom, one can be useful of the information given by
the 1 J nv nodes close to node I:
1
X I
I
& X - X&
F
X
=
- 1
max
max
AP N
DX, DY, DZ, DRX, DRY, DRZ
1 inb node 2
Bi
N
1
- 15
- 1
where Bi
I
J
N =
T
max 10 ms, X & 1,
max &
1
N
n+
N
+
2
100
1 jnv (X
2)
This method requires the census of the nodes close and the estimate “speeds” according to
each type of degree of freedom (translation “DX”, “DY”, “DZ”, and possibly rotation “DRX”, DRY' and
“DRZ”) for these close nodes.
The method programmed in Code_Aster simplifies this formula somewhat and consists, for one
degree of freedom given, I, to make starting from this position an ascending search and one
downward search on the degrees of freedom in their order of classification defined by
NUME_DDL. The first two degrees of freedom, K and L, of comparable nature respectively found
before and after the degree of freedom I are regarded as the “neighbors”. To limit the cost of this
technique, search is made once for all at the beginning of transitory calculation and the “neighbors”
are recorded in two tables of entireties.
The use of this method is started by key word VITE_MIN: “NORM”.
Handbook of Référence
R5.05 booklet: Transitory or harmonic dynamics
HT-66/03/005/A
Code_Aster ®
Version
6.4
Titrate:
Algorithm of integration of DYNA_LINE_TRAN
Date:
28/02/03
Author (S):
E. BOYERE, G. JACQUART, LIGHT A.C.,
Key:
R5.05.02-C Page
: 17/20
7.3.2 use of information at previous time
One can be also based on the information brought by the steps of previous times to estimate
minimal speed. One then estimates it by the following formula:
&X I
&
K
X
= max
,
-
10 15 ms-1
min
K <n 100
One has then:
1
X I
I
& X - X&
F
X
=
- 1
max
max
AP N
DX, DY, DZ, DRX, DRY, DRZ
1 inb node 2
Bi
N
- 15
-
1
1
with Bi
I
J
N =
T
max 10 ms, X & 1,
max &
N
n+ 2
(X K)
K
100 <n
This method is engaged by key word VITE_MIN: “MAXI”.
This method cannot be employed if speed varies too much during calculation, bus in this case
one would have with each step:
X - X
N
n-1 X I
&
T
min
7.4
Choice of the number of steps per apparent period, NR
Error analyzes and criteria of stability established for a system with only one degree of freedom (see
[bib4]) allowed to estimate the number of steps NR necessary per period connect to obtain one
good precision. These tests showed that a minimum of 20 steps per period is necessary. It
a number is skeletal by L `user in the command file thanks to the key word
“NB_POINT_PERIODE”. Its default value established to 50 leads to a precision on integration
temporal of about 1 to 2%.
The step of initial time is useful like step of maximum time in the absolute: T
= T
. Balanced
mac
initial
by a skeletal coefficient by “PAS_LIMI_RELA”, it is used as step of minimal time:
T
= PLR *
min
tinitial
7.5
Heuristics of evolution of the step of time T
N
One defines an indicator, known as “error”, on the choice of the step of time:
error = T
NR F
N
AP N
It is necessary that this indicator is lower than 1 to hope to guarantee a good temporal integration of
smaller clean period. However the adaptive diagram must concomitantly avoid the use of one
no the too small time, which would cause a overcost of calculation then, even appearance of “noises”
parasites.
According to the indicator, the algorithm will increase or decrease the step of time. One defines for that
two coefficients, CDP, the coefficient of refinement of the step of time (word key “COEF_DIV_PAS”,
default value: 1,334) and CMP, the coefficient D `amplification of the step time (word key
“COEF_MULT_PAS”, default value: 1,1).
Handbook of Référence
R5.05 booklet: Transitory or harmonic dynamics
HT-66/03/005/A
Code_Aster ®
Version
6.4
Titrate:
Algorithm of integration of DYNA_LINE_TRAN
Date:
28/02/03
Author (S):
E. BOYERE, G. JACQUART, LIGHT A.C.,
Key:
R5.05.02-C Page
: 18/20
At the time of this search of the step of optimal time, one defines a maximum iteration count of reduction of
no time, iter
, to avoid with the step of time to evolve/move in a too brutal way, which is
max
prejudicial with the command of the diagram, and not to launch a too expensive optimization.
· if the indicator of error is higher than its limiting value, that one did not exceed the limiting number
of refinement for a step of time and that the step of time remains larger than its value
minimal fixed a priori, the step of time is refined
:
1
T
T >
, iter < iter
and
,
max
T
N
>
min
N
T
T
N
NR F
CDP
N
AP N
· if the indicator shows that since five consecutive steps the step of time appears too fine, i.e.
0 7
, 5
T <
, then semi (
N T, CMP T
N
T
min
)
N
NR F
N
AP N
7.6 Algorithm
the algorithm was programmed in Code_Aster according to the following flow chart:
0 Initialization
:
X, X & given
0
0
&X = M 1
T = 0 -
-
0
(F (
)) K X C
0
&X0
T
&X
1
= &X0 -
&X
-
0
2
2
recovery of the parameters of integration:
T
initial
CMP coefficient of performance of the step of time
CDP coefficient of reduction of the step of time
PLR limits to refinement such as Dt PLR Dt
initial
NR numbers of steps of time per apparent period
itermax a maximum number of reductions of the step of time
1
with each step of time:
X, X&
, X
& known
N
N
N
- 12
T
1 = T
+ T
n+
N
N
1.0
iter=0
1.1 : temporal integration
T
1 +
T
&X 1 = &X
N
N
1 +
&Xn
n+
N
2
2
2
X
&
1 = X
+ T X
n+
N
N
1
n+ 2
&X
1 = M 1
- T1 -
1 -
n+
(
F n+)
K X
C
N
&X
+
1
n+
2
T
&X
1 = &
X
N
1 +
&X
n+
n+1
n+
2
2
Handbook of Référence
R5.05 booklet: Transitory or harmonic dynamics
HT-66/03/005/A
Code_Aster ®
Version
6.4
Titrate:
Algorithm of integration of DYNA_LINE_TRAN
Date:
28/02/03
Author (S):
E. BOYERE, G. JACQUART, LIGHT A.C.,
Key:
R5.05.02-C Page
: 19/20
1.2 calculation of the apparent frequency and the error on the step of time
1
X I
I
& X - X&
F
X
=
- 1
max
max
AP N
DX, DY, DZ, DRX, DRY, DRZ
1 inb node 2
Bi
N
error = T
NR F
N
AP N
1.2 test on the relevance of the step of time
· if error > 1et iter < iter
max
then T/CDP T
N
N
but if T < T
stop of calculation with error message
N
min
iter + 1 iter and return in 1.1
· if error > 1et iter > iter
max
then emission of an alarm and passage as in point 2.
· if error < 1 passage as in point 2
with if error < 0 7
, 5 since 5 consecutive steps:
amplification of the step of time T = semi (
N T
, CMP
max
T
N
N)
2
acceptance of the solution: possible filing of X
, X &, X
&
n+1
n+1
n+1
then N + 1 N: return in 1 for the step of next time
8 Conclusion
operator DYNA_LINE_TRAN allows the choice between several methods of temporal integration. In
their parameter setting by defect, the diagrams of WILSON and NEWMARK are implicit schemes
unconditionally stable. They thus require a linear inversion of system to each step of
times but n the other hand offer a choice of the step time which is restricted only by the smoothness with
which one wishes to describe the temporal evolution of the modelled phenomena.
Diagrams DIFF_CENTER and ADAPT are explicit what avoids to them, in the case of a matrix of
mass diagonal, an inversion of expensive matrix. But the conditional stability of this type of
diagram generally leads to the use of small steps of times, conditioned by smallest
clean period of the system. It is thus not guaranteed that the explicit diagrams are
systematically faster. That depends on the simulated phenomena. If the physics of these
phenomena requires a fine temporal discretization, the step of time employed is naturally
in the interval of stability. In the contrary case, the constraints of numerical stability involves one
inflation in the number of steps of time very expensive.
Diagram ADAPT makes profitable information on the frequential contents of the answer to adapt it
no time. The discretization of time is not thus imposed any more by the smallest clean period of
system but by its answer. That can be an advantage when the frequency of the answer evolves/moves
in time, like in the case of the impacts.
Handbook of Référence
R5.05 booklet: Transitory or harmonic dynamics
HT-66/03/005/A
Code_Aster ®
Version
6.4
Titrate:
Algorithm of integration of DYNA_LINE_TRAN
Date:
28/02/03
Author (S):
E. BOYERE, G. JACQUART, LIGHT A.C.,
Key:
R5.05.02-C Page
: 20/20
9 Bibliography
[1]
K.J. BATHE, E.L. WILSON: “Numerical Methods in Finite Element Analysis” - PRENTICE
HALL INC.
[2]
T.J.R. HUGHES: “The Finite Element Method, Linear Static and Dynamic Finite Element
Analysis " - 1987 - PRENTICE HALL Inc.
[3]
S. COURTIER-ARNOUX and coll: “Notes of the course of function: Dynamics of Structures”
- note HP-61/94/189
[4]
G. JACQUART, S. GARREAU: “Algorithm of integration to step of adaptive time in
Code_Aster - note HP-61/95/023
[5]
A.C. LIGHT: “Introduction of the explicit diagrams “centered differences” and “not of
adaptive time “in operator DYNA_LINE_TRAN of Code_Aster.
Handbook of Référence
R5.05 booklet: Transitory or harmonic dynamics
HT-66/03/005/A
Outline document