Code_Aster ®
Version
7.4
Titrate:
Quasi-static nonlinear algorithm


Date:
06/07/05
Author (S):
P. BADEL, J. LAVERNE, NR. TARDIEU Key
:
R5.03.01-D Page
: 1/28

Organization (S): EDF-R & D/AMA
Handbook of Référence
R5.03 booklet: Nonlinear mechanics
Document: R5.03.01

Quasi-static nonlinear algorithm
(operator STAT_NON_LINE)

Summary:

Operator STAT_NON_LINE [U4.51.03] of Code_Aster allows a quasi-static stress in the case of
to integrate various types of non-linearities coming from the behavior of material or of large
geometrical transformations. One describes the total algorithm of resolution here employed. The integration of the relations
of behavior themselves is described in other documents, like [R5.03.02] for the élasto-
plasticity, to which one will be able to refer for examples. For calculations in great transformations
geometrical, one will be able to consult for example the document [R5.03.20] on nonlinear elasticity into large
displacements, or the document [R5.03.21] on the thermoelastoplasticity with isotropic work hardening.

Handbook of Référence
R5.03 booklet: Nonlinear mechanics
HT-66/05/002/A

Code_Aster ®
Version
7.4
Titrate:
Quasi-static nonlinear algorithm


Date:
06/07/05
Author (S):
P. BADEL, J. LAVERNE, NR. TARDIEU Key
:
R5.03.01-D Page
: 2/28

Count

matters

1 Presentation ........................................................................................................................................... 3
1.1 General information on the relations of behavior ............................................................................... 3
1.2 Position of the nonlinear quasi-static problem ............................................................................ 4
2 Method of resolution ........................................................................................................................... 6
2.1 Principle of the method of NEWTON .............................................................................................. 6
2.2 Adaptation of the method of NEWTON to the problem arising ........................................................... 7
2.2.1 Resolution without boundary conditions dualized .................................................................. 7
2.2.2 Resolution with boundary conditions dualized ........................................................... 9
2.2.2.1 Phase of prediction ................................................................................................... 9
2.2.2.2 Loop on the iterations of NEWTON ..................................................................... 11
2.2.2.3 RIGI_MECA_TANG and FULL_MECA .......................................................................... 16
2.2.3 Case of the following loadings ............................................................................................ 16
2.3 Linear search ......................................................................................................................... 17
2.3.1 Principle ................................................................................................................................. 17
2.3.2 Minimization of a functional calculus .......................................................................................... 18
2.3.3 Method of minimization ..................................................................................................... 18
2.3.4 Application to the minimization of energy ............................................................................ 20
2.3.5 Determination of the step of advance .................................................................................... 21
2.4 Algorithm of STAT_NON_LINE ................................................................................................... 23
3 Control ................................................................................................................................................ 24
4 Great deformations ......................................................................................................................... 24
4.1 Objective ........................................................................................................................................... 24
4.2 Great plastic deformations .................................................................................................. 25
4.3 Functionality PETIT_REAC .......................................................................................................... 26
5 Bibliography ........................................................................................................................................ 28
Handbook of Référence
R5.03 booklet: Nonlinear mechanics
HT-66/05/002/A

Code_Aster ®
Version
7.4
Titrate:
Quasi-static nonlinear algorithm


Date:
06/07/05
Author (S):
P. BADEL, J. LAVERNE, NR. TARDIEU Key
:
R5.03.01-D Page
: 3/28

1 Presentation

1.1
General information on the relations of behavior

STAT_NON_LINE is the operator of Code_Aster making it possible to carry out mechanical calculations not
linear when the effects of inertia are neglected.

Calculation relates a priori only to the variables mechanical (displacements, constraints, variables
interns) by excluding any coupling with other physical phenomena (thermics,…). By
consequent, associated fields influencing the mechanical behavior (thermal fields,
metallurgical) are calculated as a preliminary by other operators (THER_LINEAIRE [U4.33.01],
THER_NON_LINE [U4.33.02]), even by other codes (SYRTHES, SATURNE,…).

This assumption suffers an exception with regard to modeling
thermo hydro-mechanical (modeling known as “THM”) for which milked STAT_NON_LINE
the whole of the coupled problem of the equations of diffusion of thermics, the pressure of
fluid and of mechanical balance [R7.01.10].

In STAT_NON_LINE, two families of behaviors are available:

·
that which corresponds to the key word factor COMP_ELAS (COMPortement ELAStique) led to
through equilibrium equation to a nonlinear system depending explicitly on the field
displacements U compared to the configuration of reference, and parameterized by the moment of
calculation (through inter alia the thermal evolution). For more details, one will be able
to defer, for example, with the document [R5.03.20] concerning elasticity into large
transformations (hyper-elasticity), or the document [R5.03.21] on the thermoelastoplasticity with
isotropic work hardening,

·
the other family, which corresponds to the key word factor COMP_INCR (COMPortement INCRémental),
is associated relations of behavior expressed by a differential equation
implicit (for example elastoplasticity, viscoplasticity, hypo-elasticity). In this case,
relation of behavior is integrated as presented for example in [R5.03.02]: in
connecting an increment of displacement U calculated starting from a mechanical state given (the state
mechanics being represented by a field of displacements U, a stress field
and a field of variables intern) with the stress field at the moment T of calculation.
The equilibrium equation thus leads to a nonlinear system out of U, but which is also
parameterized by the moment of calculation through the facts of the case (variation of the loading
mechanics and thermal evolution for example).

In both cases, one calculates the solution gradually. It is theoretically not
essential in the nonlinear elastic case, but it may be that nonthe linearity of the solution
sought either too strong for the algorithm of resolution used, and that it is essential, for
numerical reasons, to operate step by step.

It is necessary nevertheless to have for the spirit the fundamental difference between the two approaches. The elastic case
suppose the existence of a state of reference per report/ratio to which the elastic strain is written: this
state corresponds in a state without forced deformation nor. It is the “absolute” value of the loading which
create the deformation. The incremental case is based on the state previously calculated and “forgets” all
reference to the former states except that given by the internal variables. In this case, it is
variation of the loading which modifies the state of the system: in particular, one needs a variation of the field of
temperature to create thermal deformations.
Handbook of Référence
R5.03 booklet: Nonlinear mechanics
HT-66/05/002/A

Code_Aster ®
Version
7.4
Titrate:
Quasi-static nonlinear algorithm


Date:
06/07/05
Author (S):
P. BADEL, J. LAVERNE, NR. TARDIEU Key
:
R5.03.01-D Page
: 4/28

1.2
Position of the nonlinear quasi-static problem

Consequently [§1.1], it is seen that it is legitimate to consider that the nonlinear problem that one
draft has like unknown factor a displacement and is parameterized by time. That is to say thus the problem not
linear quasi-static according to:


vTR (U, T) = vTL (T)
v such as Bv = 0


Drunk = ud (T)



where:

·
T represents the variable of moment,
·
U the field of displacement taken starting from a configuration of reference,
·
L the mechanical loading to which the structure is subjected (pressure, imposed force,…),
·
the relation Drunk = ud (T) corresponds to the boundary conditions imposed in displacements
(imposed displacements, connections between degrees of freedom,…) : B is a linear operator of
the space of the fields of displacements on a space of functions defined on part of
edge of the structure, ud is a function given on this part,
·
R represents the internal forces of the problem of quasi-static mechanics nonlinear
(in the linear case, there are R (U, T) = Ku, where K is the matrix of rigidity of the structure).
The notation R (U, T) is a short cut which one will use in the continuation.

In fact, more precisely, R (U, T) is connected to the stress field by the operator of
work of virtual deformations QT [§ 2.2.1] according to the relation:

R (U, T) QT
=
(U):
éq
1.2-1

In small displacements, QT is independent of displacements; for the large ones
displacements, to refer to [R5.03.20]. Subsequently, one will place oneself on the assumption
small displacements and small deformations.

The stress field I at the moment Ti is written (U, T, T, H)
I
I I
i-1, if one notes Ti the field of
temperatures and Hi-1 last history of the structure. For the elastic behaviors, the history
does not intervene: the Hi-1 unit is thus empty. For the incrémentaux behaviors, the history is
the whole of the states (fields of displacements, constraints and variables intern) at the moment
precedent: H
=
1
{U
I
I 1, I 1, I
, T
-
-
-
- 1 I}
1.

In the general case, the dependence of the operator R is, as we saw in [§1.1],
implicit compared to time: it results from the integration of the relation of behavior in
time (for the problems of elastoplasticity for example). The dependence clarifies compared to
time is very rare: it in particular in the case of appears relations of fascinating behavior in
count a phenomenon of work hardening by time known as time-hardening.
Handbook of Référence
R5.03 booklet: Nonlinear mechanics
HT-66/05/002/A

Code_Aster ®
Version
7.4
Titrate:
Quasi-static nonlinear algorithm


Date:
06/07/05
Author (S):
P. BADEL, J. LAVERNE, NR. TARDIEU Key
:
R5.03.01-D Page
: 5/28

The dualisation of the boundary conditions of Drunk Dirichlet = ud (T) led to the problem
according to [R3.03.01]:
R

(U, T) + BT = L (T)


éq
1.2-2
Drunk = ud (T)



The unknown factors are now at any moment T the couple (U,), where represents the “multipliers
of Lagrange " of the boundary conditions of Dirichlet [R3.03.01]. Vector BT is interpreted like
opposite of the reactions of support to the corresponding nodes.

The formulation of the quasi-static problem consists in expressing the balance of the structure (forces
external interns = forces) for a succession of moments of calculation {Ti}
who parameterize it
1i I
loading:
R

(U, T) + BT
= L (T)
I I
I
I

éq
1.2-3
Drunk
= ud (T)


I
I

what amounts cancelling in (U, T)
I
I I the vector F (U, T) defined by:

R (U, T) + BT - L (T)
F (U, T) =


Drunk - ud (T)


The state of the structure in t0 is supposed to be known. One carries out I increments (or not) of load definite
as follows:

not
charge n°
1
2
I
moment
t0 T1 t2
….
Ti-1 Ti

The unknown factors are calculated in an incremental way by the total algorithm of resolution (even for
elastic behaviors). From (ui-1, i-1), solution satisfying Ti-1 balance, one
determine ui and I, which will make it possible to obtain the Ti solution:

U

= U
+ U
I
I

- 1
I



I = I +
- 1
I

The increments ui and I are initially estimated by linearizing the problem compared to time
around (ui, T
- 1
i-1 i-1) (phase known as of prediction or Euler [bib 3]). Then a method is used
NEWTON or one of its alternatives to solve the equation [éq 1.2-3] in an iterative way (one
calculate a continuation (one, N
I
I)).

For the incremental relations of behavior, one needs to know out of Ti-1 it
stress field i-1 and the field of variables intern i-1 (Cf. [R5.03.02] for one
example).
Handbook of Référence
R5.03 booklet: Nonlinear mechanics
HT-66/05/002/A

Code_Aster ®
Version
7.4
Titrate:
Quasi-static nonlinear algorithm


Date:
06/07/05
Author (S):
P. BADEL, J. LAVERNE, NR. TARDIEU Key
:
R5.03.01-D Page
: 6/28

2
Method of resolution

2.1
Principle of the method of NEWTON

The method of NEWTON is a traditional method of resolution of the equations of the type:

F (X) = 0

where F is a vector function (nonlinear) vector X.

It consists in building a vector series {xn} converging towards solution X. To find
N
the new one reiterated xn+1, one approaches F (xn+1) by a development limited to command 1 around xn, and
one expresses that F (xn+1) must be null:

0
N 1
+
N
'N
N 1
= (
) () +
+
F X
F X
F (X) (X
- xn),

that is to say:

- 1
F 'xn xn+1 - xn
= - F xn
or xn+1 = xn - [F '
() (
)
()
(xn)] F (xn)

Recall:

F '(xn) is the tangent linear application associated the function F. The derivative as in point X
in the direction H is defined like:

'
F (X + H) - F (X)
F (X) H = lim
.
0


The matrix of F '(X) in the bases chosen for the vector spaces concerned
be called the matrix jacobienne F as in point X.

When F is related to an Euclidean vector space with actual values, F '(X) is one
linear form, and one can show that there is a vector (single), noted F (X) and called it
gradient of F, such as:

F '(X) H = HT F
(X) (produced scalar of H and the gradient).

When one is close to the solution, the convergence of the method of NEWTON is quadratic (it
numbers of 0 after the comma in the double error with each iteration: 0.19 - 0.036 - 0.0013 -
0.0000017 for example). But this method (using the true tangent) has several disadvantages:

·
it requires the calculation of the tangent to each iteration, which is expensive (especially in
case where a direct solvor is used),
·
if the increment is large, the tangent (known as coherent or consistent) can lead to
divergences of the algorithm,
·
it can not be symmetrical, which obliges to use particular solveurs.

For this reason one can use other matrices in the place of the tangent matrix:
stamp elastic, a tangent matrix obtained before, the symmetrized tangent matrix,…
[§2.2.1].
Handbook of Référence
R5.03 booklet: Nonlinear mechanics
HT-66/05/002/A

Code_Aster ®
Version
7.4
Titrate:
Quasi-static nonlinear algorithm


Date:
06/07/05
Author (S):
P. BADEL, J. LAVERNE, NR. TARDIEU Key
:
R5.03.01-D Page
: 7/28

2.2
Adaptation of the method of NEWTON to the problem arising

2.2.1 Resolution without boundary conditions dualized

If one initially forgets the boundary conditions of Dirichlet, one must solve one
system of the form:

R (U, T) = Lméca
I I
I
,

where Lméca
I
indicate the mechanical loading at the moment Ti.

Using the notations of [§ 2.1], that reverts cancelling the vector function F defined by:

F (U, T) = R (U, T) - Lméca
I I
I I
I


The nodal forces R can symbolically be noted QT, where QT is the matrix associated with
the operator divergence (expression of the agricultural work of virtual deformations), with
(U, T, T, H). One a:
I
I
I
i-1

(QT) =
(U):(W) D


K
K


(L) =
F W D +
T W D




I K
I
K
I
K


where
·
wk indicates the basic function associated with the kth degree of freedom with the structure,
·
F indicates the voluminal forces applying at the moment T to,
I
I
·
T indicates the surface forces applying at the moment T to the border of.
I
I

The matrix Q depends on displacements U in great displacements [R5.03.20].

The application of the method of NEWTON results in solving a linear succession of problems of the type
(N is the number of the iteration of NEWTON, I that of the variable of moment):

KN un+1 = Lméca - RN
I
I
I
I,

where un+1 is noted
un+
=
1 - one
I
I
I the increment of displacement between two iterations of NEWTON
successive. The matrix K N
N
N
I is the matrix of tangent rigidity in ui and vector IH represents them
forces intern with the nth iteration of NEWTON of I ème not of load. The Lméca quantity - R N
I
I
represent the not balanced forces, which one can also call the “residue”.

The matrix K nor is the matrix of the tangent linear application of the function F, it is thus worth:

F
R
Lméca
KN
I
I =
=
-

.
U


N
U
N
U
(U T,)
(U T
N
I I
I, I)
(U T
I, I)
Handbook of Référence
R5.03 booklet: Nonlinear mechanics
HT-66/05/002/A

Code_Aster ®
Version
7.4
Titrate:
Quasi-static nonlinear algorithm


Date:
06/07/05
Author (S):
P. BADEL, J. LAVERNE, NR. TARDIEU Key
:
R5.03.01-D Page
: 8/28

In the absence of following forces [§2.2.3], the second term is null. The matrix K nor is thus the derivative
at the plain point of the nodal forces (or interns) compared to displacements:

R
Kni =

U (one T
I, I)

A small error in the evaluation of the residue can have serious consequences, because it is its calculation
exact which guarantees, if one converges, that it will be towards the sought solution. On the other hand, it is not
always necessary to use the true tangent matrix, whose calculation and factorization are expensive.
For example, an alternative of the method uses the elastic matrix K 0.

Method using the true tangent matrix K nor (known as also coherent or consistent matrix)
be called the method of NEWTON; methods using of other matrices (such as for example
stamp elastic K 0) are called methods of modified NEWTON. The choice enters a matrix
tangent (the last obtained or a preceding matrix) is carried out in Code_Aster by
the intermediary of key word MATRICE: “TANGENTE” or “ELASTIQUE” of the key word factor NEWTON. Moreover,
it is possible to use a matrix of discharge, i.e. of a matrix with internal variables
constants (the evolution of nonthe linearities is thus not taken into account in this matrix), in
below of a certain step of time, for certain laws of behavior. One will refer to
documentation [U4.51.03] for the use of this functionality.

The method of NEWTON with consistent tangent matrix can be illustrated simply using
diagram of the figure [Figure 2.2.1-a].

R
Li
L-R2
2
R L-R1
1
R
R0
0
1
U
U
2
U
U
I
I
I
I
U
1
U
2
U
I
I

Appear 2.2.1-a
Handbook of Référence
R5.03 booklet: Nonlinear mechanics
HT-66/05/002/A

Code_Aster ®
Version
7.4
Titrate:
Quasi-static nonlinear algorithm


Date:
06/07/05
Author (S):
P. BADEL, J. LAVERNE, NR. TARDIEU Key
:
R5.03.01-D Page
: 9/28

2.2.2 Resolution with boundary conditions dualized

When one takes into account the conditions of imposed displacements, the system to be solved is written:

R

(U, T) + BT
= Lméca
I I
I
I


Drunk
= ud


I
I

One will use the symbol to note the increments since the preceding balance (out of Ti-1) of
displacements, of the parameters of Lagrange, the loading and imposed displacements:

U = U - U
I
I
i-1
= -

I
I
i-1


Lméca = Lméca - Lméca

I
I
i-1
ud = ud - ud
I
I
i-1

Knowing (U
,)
i-1
1
I
, the couple (U,)
I
I will be determined by the resolution of the system:

R

(U
+ U
, T) + BT (
+) = Lméca + Lméca
i-1
I
I
i-1
I
i-1
I

éq
2.2.2-1
B (U
+ U
) = ud + ud


i-1
I
i-1
I

One will use a method of NEWTON to solve this system. However, the experiment shows
that the convergence of the method of NEWTON is strongly dependant on a judicious choice on
the initial estimate: “more the initial estimate is close to the solution, plus the algorithm converges quickly”.
To start the iterative process of the method, it is thus useful to determine “a good” increment
initial (u0, 0)
I
I. For that, one linearizes compared to time the continuous problem: it is what one
call the phase of prediction (or of initialization). One connects with the loop of the iterations of
NEWTON which makes it possible, with convergence, to obtain the values of (U,)
I
I, and thus those of
(U
)
I, I.

2.2.2.1 Phase of prediction

One linearizes the system [éq 2.2.2-1] compared to time around (U
,)
i-1
1
I
; by taking account of
the balance obtained at the moment Ti-1, one obtains the equations making it possible to calculate values
predictive (u0, 0)
I
I:

R
R
u0 BT0
+
= Lméca -
T


U
I
I
I

I
T

U

I
T
- 1
i-1
Drunk = ud

I
I

R
= ():(T T (T))

indicate the total differential of R
Q U
,
compared to T. This notation
T Ti-1

T
particular be to draw the attention to the fact that
=
+

as one sees it in
T
T
T T
continuation.
Handbook of Référence
R5.03 booklet: Nonlinear mechanics
HT-66/05/002/A

Code_Aster ®
Version
7.4
Titrate:
Quasi-static nonlinear algorithm


Date:
06/07/05
Author (S):
P. BADEL, J. LAVERNE, NR. TARDIEU Key
:
R5.03.01-D Page
: 10/28

Therefore, while using [éq 1.2-1] and dependence from report/ratio with time and compared to
temperature:



T

T
Q


0
0


Q:
+
: U + BT
= Lméca - QT:
T
- QT:
T



U
U
I
I
I
I
I

T
T


ui-1
U

T
T
I
I
I
-
-
-
1
1
1
B
U
= ud

I
I

If the constraints depend on other variables, such as for example of proportions Z of
various metallurgy i.e. phases (T, T (T), Z (T)) , one sees appearing to the second member

the term corresponding: QT:
Z
. It is supposed thereafter that the effects
I
Z Ti-1
metallurgical [R4.04.02] are integrated into the second member representative of the effects
thermics. Moreover, one currently does not hold account for the phase of prediction of
dependence clarifies constraints compared to time. Lastly, dependence of
stamp QT compared to displacements is neglected on the assumption of small
QT
displacements: the term
, known as term of geometrical rigidity, thus disappears from
U ui-1
the preceding equation.

With obvious notations, the system of equations obtained is also written:

K


u0

+ BT
0

= Lméca + Lther
i-1
I
I
I
I


éq
2.2.2.1-1
B u0

=
ud


I
I

Let us note that one can replace the Ki-1 matrix formally, derived from R compared to U in ui-1
by the elastic matrix K 0.

For the developers, let us specify that the calculation of the tangent matrix at the time of the phase of prediction
is carried out via the option of calculation RIGI_MECA_TANG [§ 2.2.2.3].

The effective increment of loading Li appearing to the second member includes/understands two terms:

·
the mechanical increment of loading Lméca
I
, composed of the dead loads
(independent of the geometry, like gravity) and of the following loads (dependant
geometry, like the following pressure [R3.03.04]). Actually, there is cases (it
first increment of load, for example) where méca
L
is unknown: this Lméca fact is
I 1
-
I
rather calculated by the relation
méca
L
= méca
L
- T
Q
- T
B, which returns exactly to
I
I
I 1
-
I 1
-
even, in the case running, by taking account of balance with the increment (i-1). One
will notice whereas the expression utilized the multipliers of Lagrange to the increment
(i-1), which is sometimes unknown (with the first increment of load, for example). But linearity
boundary conditions in imposed displacements, which results in the fact that the matrix
B is constant, allows in this case an abuse language: one poses = 0 and
I
0
truth
0 truth
=
+
I
I 1
-
(the true exhibitor corresponds to the variables appearing indeed
I
in the system of equations above) and one solves:

0
T
0

K U + B = méca
L
+
ther
L
- T
Q
I 1
-
I
I
I
I
I 1

-,
0

Drunk = D
U - U
B
I
I
I 1
-
Handbook of Référence
R5.03 booklet: Nonlinear mechanics
HT-66/05/002/A

Code_Aster ®
Version
7.4
Titrate:
Quasi-static nonlinear algorithm


Date:
06/07/05
Author (S):
P. BADEL, J. LAVERNE, NR. TARDIEU Key
:
R5.03.01-D Page
: 11/28

·
the increment of loading says “thermal” (and metallurgical [R4.04.02]) Lther
I
, resulting from
derivation of the nodal forces compared to the temperature and which is an estimate of the effect
of a variation in temperature. Indeed, if one notes K the module of compression
hydrostatic and the thermal dilation coefficient, the increment of loading
“thermal” is written (since R = QT and ther = - 3K T
I where I is the matrix
I D
D
identity):


Lther = - QT:
T = QT:3 K T
I,
I
I
I
T
D
Ti-1

In the elastic case, they are the nodal forces associated a thermal dilation (it
is not strictly speaking a loading, that is assimilated rather to the effect of a deformation
initial). This estimate is used in the phase of prediction and the criterion of stop. If
thermal dilations make leave the structure of the elastic range (plasticity by
example), this estimate will be corrected at the time of the iterations of NEWTON.

Note:

A particular case relates to the use of an excitation of the type TYPE_CHARGE: Meaning “DIDI”
Differential Dirichlet, i.e. compared to the initial state. That consists, for the conditions with
D
limits of the blockings type, to impose, not B U = U, but B (U - U) = ud
0
. In this case, it
system to be solved in the phase of prediction for the new increment of load is:

K U0 + T
B 0 = méca
L
+
ther
L
- T
Q
i-1
I
I
I
I
i-1

Bu0 = D
U
I
I

2.2.2.2 Loop on the iterations of NEWTON

One must find the values (U,)
I
I of the increments of displacements and parameters of Lagrange
since the values (U
,)
i-1
1
I
obtained with preceding balance (urgent Ti-1). One takes as
initial values (u0, 0)
I
I obtained at the end of the phase of prediction, before beginning them
iterations of the method of NEWTON.

With each iteration, one must solve a system allowing to determine (un+1 n+1
,
), increments
I
I
displacements and parameters of Lagrange since the result (one
N
,)
I
I of the iteration
the preceding one:

KN


un+1 + BT
N
+1 = Lméca - R (one) - BT N
I
I
I
I
I
I


éq
2.2.2.2-1
B un+


1
I
= 0

with R (one) = QT N
N
N
I
I, constraints I being calculated starting from ui displacements by
the intermediary of the relation of behavior of the material [§1.1]. In fact, in the case of them
incrémentaux behaviors, nor is calculated starting from (i-1, i-1) and of the increment of
deformation (un+1) induced by the increment of un+1 displacement = one + un+1 - U
I
I
I
I
i-1.
Handbook of Référence
R5.03 booklet: Nonlinear mechanics
HT-66/05/002/A

Code_Aster ®
Version
7.4
Titrate:
Quasi-static nonlinear algorithm


Date:
06/07/05
Author (S):
P. BADEL, J. LAVERNE, NR. TARDIEU Key
:
R5.03.01-D Page
: 12/28

Stamp “tangent”

As in the phase of prediction, one is not obliged to use the true tangent matrix K nor. In
private individual, operator STAT_NON_LINE authorizes the use of the elastic matrix K 0, or
reactualization of the tangent matrix all i0 not of time (key word REAC_INCR) or all n0
iterations of NEWTON (key word
N
REAC_ITER). Thus, the matrix K I perhaps replaced by one
stamp K
m
J, with J I, or a Ki matrix, with m N.

Caution: a “stiff” matrix too does not pose problems of stability but can produce
a very slow convergence; a “flexible” matrix too can lead to divergence, it is
advised in this case to make linear search [§ 2.3].

It is difficult to give a rule making it possible to know when one must reactualize the tangent matrix:
that strongly depends on the degree of nonlinearity of the problem and the size of the increments of load.

In discharge, it is recommended either to use the elastic matrix for the phase of prediction
and of resolution, the elastic matrix for the phase of prediction then the matrix is to use
tangent for the resolution.

The figures hereafter illustrate the various possibilities of reactualization of the tangent matrix:
stamp elastic K 0 used with each iteration [Figure 2.2.2.2-a], tangent matrix reactualized with
each iteration and with each step of time [Figure 2.2.2.2-b], reactualized tangent matrix all them
i0 not of time (i0 =1 here) [Figure 2.2.2.2-c], and stamps tangent reactualized all the n0 iterations
NEWTON (n0 =2 here) [Figure 2.2.2.2-d].

R
Li
U
Appear 2.2.2.2-a: use of the elastic matrix
Handbook of Référence
R5.03 booklet: Nonlinear mechanics
HT-66/05/002/A

Code_Aster ®
Version
7.4
Titrate:
Quasi-static nonlinear algorithm


Date:
06/07/05
Author (S):
P. BADEL, J. LAVERNE, NR. TARDIEU Key
:
R5.03.01-D Page
: 13/28

R
Li
U
Appear 2.2.2.2-b: use of the true revalued tangent matrix
with each iteration

R
Li+1
L I
U

Appear 2.2.2.2-c: use of the tangent matrix
revalued with each step of time
Handbook of Référence
R5.03 booklet: Nonlinear mechanics
HT-66/05/002/A

Code_Aster ®
Version
7.4
Titrate:
Quasi-static nonlinear algorithm


Date:
06/07/05
Author (S):
P. BADEL, J. LAVERNE, NR. TARDIEU Key
:
R5.03.01-D Page
: 14/28

R
L I
U

Appear 2.2.2.2-d: use of the revalued tangent matrix
all 2 iterations of NEWTON

Method of modified NEWTON (using another matrix that the consistent tangent matrix)
converge less quickly than the method of traditional NEWTON, but is less expensive. It is
all the more economic as the number of degrees of freedom of the system is high. This is why one
can advise, when the system is of big size, to keep the same tangent matrix during
some iterations. Lastly, let us not forget to announce that in certain cases, it is calculation with
stamp elastic which is fastest in terms of calculating time, even if the iteration count
carried out is much more important (they are cheap iterations since the matrix is not
calculated and factorized that only once); moreover, it is the elastic matrix which it is recommended
to use during the discharges.

As the equation [éq 2.2.2.2-1] shows it, it is necessary to carry out with each iteration of NEWTON calculation
possible of the new tangent matrix K N
N
T N
I and of the “nodal forces” R + B
I
I: for
developers, let us specify that these operations are carried out by the option of calculation FULL_MECA
(RAPH_MECA if the tangent matrix is not recomputed).

The tangent matrix obtained by option RIGI_MECA_TANG and the tangent matrix obtained
by option FULL_MECA are in general calculated differently [§ 2.2.2.3].

Updating of the unknown factors

The updating of displacements and their increments is done as follows, that of the parameters of
Lagrange is carried out in an identical way:

un+


1
= one + un+1
I
I
I
un+1

= one + un+


1
I
I
I

Handbook of Référence
R5.03 booklet: Nonlinear mechanics
HT-66/05/002/A

Code_Aster ®
Version
7.4
Titrate:
Quasi-static nonlinear algorithm


Date:
06/07/05
Author (S):
P. BADEL, J. LAVERNE, NR. TARDIEU Key
:
R5.03.01-D Page
: 15/28

Criteria of convergence

There are three criteria of convergence.

Criterion RESI_GLOB_MAXI consists in checking that the infinite standard of the residue
T
N
Q + T N
B - méca
L
, is lower than the value specified by the user. He is not advised
I
I
I

to use this criterion alone, because one cannot easily have an idea of the absolute orders of magnitude
acceptable.

The criterion of convergence chosen by defect amounts checking that the residue is sufficiently small,
as previously, and this relative with a quantity representative of the loading (it is it
criterion RESI_GLOB_RELA):

QT N
+ BT N
- Lméca
I
I
I

,
éq
2.2.2.2-2
Lméca + Lther - BT N

I
I
I


being desired relative precision given by the user (or the default value of 10-6) and


the standard of the maximum.

One can notice that, in the case of use of RESI_GLOB_RELA, the criterion can become
singular if the external loading Lméca
Lther
BT N
+
-
becomes null. This can arrive in
I
I
I
case of total discharge of the structure. If such a case of figure arises (i.e loading
10-6 time smaller than the smallest loading observed until now the increment), the code
use then criterion RESI_GLOB_MAXI with as value that observed with convergence
preceding increment. When the loading becomes again not no one, one returns to the criterion
initial.

The third criterion is criterion RESI_REFE_RELA: the idea of this criterion is to build a force
nodal of reference, which will be used to estimate term in the long term, the nullity (approximate) of the residue:

J
{
}
ddls
(T N T N méca
Q + B - L
F
I
I
I
)
ref.
J
J

More precisely, the nodal force of reference is built starting from the data of an amplitude of
constraint of reference (in mechanics; in the case of the THM, it is necessary to give one
reference to each physical phenomenon entering the generalized constraint)
ref.
. From
this amplitude of constraint of reference, one defines the tensor test
: it is null for all these
J
components, except j-ième which is worth ref.
. One defines then, for each node of each element
following nodal force (the goal being to give an idea of the importance of a component in a point
of Gauss of the constraint on the nodal force):

NR
M
~
1
E
R =
B



I

test
I, J
J
NR =1 j=1

with NR the number of points of Gauss of the element, M the component count of the tensor of
constraint; the exhibitor being used to note the evalutation of quantity at the point of Gauss, are them
weight of the points of Gauss.
The nodal force of reference is then defined by:

ref.
~e
F
= min R
I
I
E
I
where is the whole of the elements connected to node I.
I
Handbook of Référence
R5.03 booklet: Nonlinear mechanics
HT-66/05/002/A

Code_Aster ®
Version
7.4
Titrate:
Quasi-static nonlinear algorithm


Date:
06/07/05
Author (S):
P. BADEL, J. LAVERNE, NR. TARDIEU Key
:
R5.03.01-D Page
: 16/28

Convergence is issued carried out when all the criteria specified by the user are checked
at the same time. By defect, one makes a test on the relative total residue (RESI_GLOB_RELA) and numbers it
maximum of iterations of NEWTON (ITER_GLOB_MAXI).

2.2.2.3 RIGI_MECA_TANG and FULL_MECA

It is important to stress that the tangent matrix resulting from option RIGI_MECA_TANG and the matrix
tangent resulting from option FULL_MECA are in general not identical.

Let us suppose that one reached convergence for the moment Ti-1 and that one seeks now with
to obtain balance for the moment following Ti. The matrix resulting from IH GI _MECA_TANG comes from one
linearization of the equilibrium equations compared to time around (U
,)
i-1
1
I
i.e around
balance at the moment Ti-1. It is thus the tangent matrix of the system with convergence at the moment Ti-1.
On the other hand, the matrix resulting from FULL_MECA comes from a linearization from the equilibrium equations by
report/ratio with displacement around (one
N
,) i.e around balance at the moment T.
I
I
I

Moreover, one can interpret the differences between RIGI_MECA_TANG and FULL_MECA in others
terms. One can thus show that the matrix resulting from RIGI_MECA_TANG corresponds to the operator
tangent of the problem continuous in time, known as also problem of speed (and connects the speed of constraint
at the speed of deformation), whereas the matrix resulting from FULL_MECA corresponds to the operator
tangent of the problem discretized in time [bib1]. The document [R5.03.02] gives the expression in
each of the two cases for the relation of elastoplasticity of Von Mises with isotropic work hardening.

It is pointed out that the processing of a relation of behavior [R5.03.02 § 5] consists with:

·
to calculate constraints N
N
I and the variables intern I starting from the initial state (
,)
i-1
i-1
and of the increment of linked displacement,
·
to calculate the nodal forces R N = QT N
N
I
I and reactions of BTi support,
·
to calculate (possibly) the tangent matrix (option RIGI_MECA_TANG for the phase of
prediction, option FULL_MECA for the iterations of NEWTON).

2.2.3 Case of the following loadings

A following loading (in mechanics) is a loading which depends on the geometry of the structure,
as for example the pressure which is exerted in the direction opposed to the normal (or the forces
of inertia in a reference mark not galiléen). Thus, when the structure becomes deformed with the evolution of
charge, the loading, expressed in an absolute reference mark, is transformed. The loads which do not depend
no the geometry of the structure are called dead or fixed loads (for example,
gravity). To indicate that a load must be treated like a following load in
STAT_NON_LINE, one indicates TYPE_CHARGE: “SUIV” under key word EXCIT.

A mechanical loading Lméca (T) comprising following loads is written Lfixe () + Lsuiv
T
(U, T)
(the fixed exhibitor indicates here the died loads, and suiv the following loads). The system of equations
to solve becomes then:

R

(U, T) + BT
= Lfixe (T) + Lsuiv (U, T)
I I
I
I
I I


Drunk
= ud


I
I

Operations of derivation allowing to write the phase of prediction and the iterations of the method
NEWTON thus utilize the derivative of Lsuiv compared to displacements U.
Handbook of Référence
R5.03 booklet: Nonlinear mechanics
HT-66/05/002/A

Code_Aster ®
Version
7.4
Titrate:
Quasi-static nonlinear algorithm


Date:
06/07/05
Author (S):
P. BADEL, J. LAVERNE, NR. TARDIEU Key
:
R5.03.01-D Page
: 17/28

The phase of prediction becomes:


Lsuiv
(

K
) u0

BT
0

Lfixe Lsuiv Lther

I -
+
=
+
+
1
U
I
I
I
I
I

U

I

- 1
B

u0

=
ud

I
I

and the iterations of NEWTON consist in solving the system:


Lsuiv
(

KN -
) un+1 + BT
n+
1 = Lfixe + Lsuiv (one) - R (one) - BT N

I

U
I
I
I
I
I
I

one

I
B
un+

1
I
= 0

Thus, with the beginning of each step of load (prediction) and with each iteration of NEWTON, one must
to calculate a matrix of “rigidity” (- Lsuiv/U) and a Lsuiv vector one
(
) related to the loadings
one
I
I
follower.

The only loads which can be treated like following loads in the current state of
operator STAT_NON_LINE are:

·
pressure for modelings 3D, 3d_SI, D_PLAN, D_PLAN_SI, AXIS, AXIS_SI,
C_PLAN, C_PLAN_SI [R3.03.04],

·
the loading of gravity for elements CABLE_POULIE [R3.08.05], elements with three
nodes comprising a pulley and two strands of cables: the force of gravity being exerted on
the element depends on the respective lengths of the two strands,

·
the centrifugal force in great displacements, which for a disk speed is given
by:

[OM]
D
=
[(OM +

)
;
0

D


U]


·
the loading of gravity for all modelings THM of the porous environments not
saturated [R7.01.10]: indeed, the density depends on the nodal variables U, p and T
to take account of the relations of behavior of the géomatériaux one.

2.3 Seek
linear

Linear search here exposed relates to linear search in the absence of control. For
description of linear search in the presence of control, one will refer to documentation
[R5.03.80].

2.3.1 Principle

The introduction of linear search into operator STAT_NON_LINE results from a report:
method of NEWTON with consistent matrix does not converge in all the cases of figure,
in particular when one leaves too much far from the solution. In addition, the use of matrices other than
stamp tangent consistent can, when they are too “flexible”, lead to divergence.
linear search makes it possible to be guarded against such divergences.
Handbook of Référence
R5.03 booklet: Nonlinear mechanics
HT-66/05/002/A

Code_Aster ®
Version
7.4
Titrate:
Quasi-static nonlinear algorithm


Date:
06/07/05
Author (S):
P. BADEL, J. LAVERNE, NR. TARDIEU Key
:
R5.03.01-D Page
: 18/28

It consists in considering (un+1, n+1
I
I
), either like the increment of displacements and of
parameters of Lagrange, but as a direction of search in which one will seek with
to minimize a functional calculus (the energy of the structure). One will find a step of advance in this
direction, and the updating of the unknown factors will consist in making:


un+1 = one + un+1
I
I
I
n+1 = N + n+


1
I
I
I


In the absence of linear search (by defect) the scalar is of course equal to 1.

2.3.2 Minimization of a functional calculus

In order to be better convinced of the founded good of linear search, one can interpret the method of
NEWTON like a method of minimization of a functional calculus (if matrices
tangents are symmetrical). We insist on the fact that the equations obtained are
rigorously those of the method of NEWTON exposed in [§2.2] and that only the way of y
to arrive is different.

“The talk the dualisation of the boundary conditions of Dirichlet and we place Forget” to simplify
on the assumption of the small deformations. The functional calculus is considered:

J: V
IR

U J (U) = W ((U)) D F .u D t.u D




where the density of free energy W makes it possible to connect the tensor of the constraints to the tensor of
W
deformations linearized by the relation =.

The functional calculus J being convex, to find the minimum of J is equivalent to cancel its gradient, that is to say:
J (U) .v = 0 v V,

what is exactly Principe of Travaux Virtuels since:
J (U) .v = (U):(v) D - F .v D t.v D





Thus, to solve the equations resulting from Principe of Travaux Virtuels (bases problem formulated in
[§1.2]) is equivalent to minimize the functional calculus J which represents the energy of the structure (energy
decreased intern of the work of the external forces F and T).

2.3.3 Method of minimization

Minimization is made in an iterative way, classically in two times with each iteration:

·
calculation of a direction of search along which one will seek reiterated according to,
·
calculation of the “best” step of advance in this direction: un+1 = one +.
Handbook of Référence
R5.03 booklet: Nonlinear mechanics
HT-66/05/002/A

Code_Aster ®
Version
7.4
Titrate:
Quasi-static nonlinear algorithm


Date:
06/07/05
Author (S):
P. BADEL, J. LAVERNE, NR. TARDIEU Key
:
R5.03.01-D Page
: 19/28

In a problem of minimization, the natural idea is to advance in the direction opposed to the gradient
functional calculus, which is locally the best direction of descent since this direction is
normal with the lines of isovaleurs and directed in the direction of the decreasing values [Figure 2.3.3-a].

U N
- J (U N)
J (U) = cste

Appear 2.3.3-a

However, it is possible to improve the choice of the direction of descent by using this method
of gradient in a metric news. It is what will enable us to find the equations
traditional of the method of NEWTON [§2.3.4].

Let us take the simple example of a problem with two variables X and there for which the functional calculus with the form
of an ellipse whose minimum is in (/has,/b):

1
1
J (X, y) =
ax2 + by2 - X - y
2
2

While choosing like direction of descent the reverse of the gradient of J, one passes from one reiterated to the following
(only let us reason on X) by:

xn+1 = xn - (axn -),

who does not point towards the solution since the normal in a point of an ellipse does not pass in general
not by the center of the ellipse [Figure 2.3.3-b].

y

B

X
has

Appear 2.3.3-b
Handbook of Référence
R5.03 booklet: Nonlinear mechanics
HT-66/05/002/A

Code_Aster ®
Version
7.4
Titrate:
Quasi-static nonlinear algorithm


Date:
06/07/05
Author (S):
P. BADEL, J. LAVERNE, NR. TARDIEU Key
:
R5.03.01-D Page
: 20/28

On the other hand, if a change of variables is carried out so that the isovaleurs of J become
circles [Figure 2.3.3-c]:


X = has X


y = B y



1


X y =
X 2 + y2 -
X -
y

J (,)
(
)
2
has
B

y

B

X
has

Appear 2.3.3-c

The use of the opposite direction of the gradient of J then makes it possible to obtain the solution in one
iteration:




X n+1
X N
X N
xn+
=
-
-
=

1
(
)
=
has
has
has

Thus, the use of the method of gradient in the metric news allows a convergence
immediate. In a more complicated case (functional calculus convex but different from an ellipse),
convergence is not instantaneous but the change of variables makes it possible to reduce appreciably
the iteration count necessary.

2.3.4 Application to the minimization of energy

To simplify, one will place oneself in the discretized linear case where the functional calculus J is worth:

1
J (U) =
uTK U - uTL
2

O one notes K the matrix of rigidity of the structure, and L the vector of the imposed loadings.
Handbook of Référence
R5.03 booklet: Nonlinear mechanics
HT-66/05/002/A

Code_Aster ®
Version
7.4
Titrate:
Quasi-static nonlinear algorithm


Date:
06/07/05
Author (S):
P. BADEL, J. LAVERNE, NR. TARDIEU Key
:
R5.03.01-D Page
: 21/28

To minimize J, we will use the same method of descent as previously while making with
precondition a completely similar change of variables. The matrix K being symmetrical defined
positive, its eigenvalues are real positive: one can thus define the “square root” of K that
K will be noted (also symmetrical). One poses U =
K U, the minimization of J is then
equivalent to that of:

1
J (U) =
U U - U
K - 1
T
T
L
2


By taking as direction of descent the opposite direction of the gradient of J, one obtains:

un+1
one
one
K -
=
-
-
1
(
L)

Maybe, while returning to the initial variables:

un+1
one
K -
=
-
1 (K one - L)

Or:

K (un+1 - one) =
L - K one

One finds the equations of the method of NEWTON: the matrix K is Hessienne of
functional calculus J (matrix of the derived second) and the second member is the difference of the loading
and of the nodal forces, also called residue.

Thus method of NEWTON perhaps interpreted like resulting from the minimization of energy
structure via a method of gradient applied after a change of metric.

2.3.5 Determination of the step of advance

Let us return to the real problem, that of the resolution of R (U, T) = L
I I
I. This problem can be interpreted
like the minimization of:

W
T
(ui) - ui Li,

where W (ui) corresponds to the discretization, on the basis of function of form, of energy interns
structure
((U))
W D.


One calculated by the method of NEWTON an increment of un+1 displacement
I
who, in the problem
of minimization, is interpreted like a direction of search. One will calculate the step of advance
in this direction allowing to minimize the value of the functional calculus:

Min {W N
(U +
n+
U 1) -
N
L (U +
n+1
I
I
I
I
ui)}
R
Handbook of Référence
R5.03 booklet: Nonlinear mechanics
HT-66/05/002/A

Code_Aster ®
Version
7.4
Titrate:
Quasi-static nonlinear algorithm


Date:
06/07/05
Author (S):
P. BADEL, J. LAVERNE, NR. TARDIEU Key
:
R5.03.01-D Page
: 22/28

To find the minimum of this scalar function of that one will note F (), one seeks the point where
its derivative is cancelled (that amounts making orthogonal the final residue and the direction of search):

F '() = [N T
+1
T
N
n+
U
1
I
] [Q (ui +ui) - Li] = 0

(F '() is the projection of the residue on the direction of search)

With the notations of [§2.2.2] and by taking of account reactions of support, the scalar equation with
to solve to determine the step of advance, is written:

[T
un+1] [QT (a un+1) BT (N n+
+
+
+
1) - Lméca
I
I
I
I
I
I
] = 0

So that the determination of is not too expensive, one uses a method of secant of which it
numbers maximum iterations is fixed by the user. The method of secant can be interpreted like
a method of NEWTON where the derivative at the point running is approached by the direction joining it
not running and the preceding point:

p - p-1
p-1 p
G - p p
G - 1
p+1 = p
p
-
G
=
p
p-1
p
p
G - G
G - G - 1
,

where G p was noted
F
p
= '().

One leaves 0 = 0 and 1 = 1. The method of secant has a command of convergence of about 1.6
[bib2]. It is represented schematically on the figure [Figure 2.3.5-a].

1
G = F '()
4
0
1

3
2
0

Appear 2.3.5-a

At the end of linear search, one brings up to date displacements and parameters of Lagrange by:


un+1 = one + un+1
I
I
I
n+1 = N + n+


1
I
I
I

Handbook of Référence
R5.03 booklet: Nonlinear mechanics
HT-66/05/002/A

Code_Aster ®
Version
7.4
Titrate:
Quasi-static nonlinear algorithm


Date:
06/07/05
Author (S):
P. BADEL, J. LAVERNE, NR. TARDIEU Key
:
R5.03.01-D Page
: 23/28

The test of convergence carries:

·
on the maximum number of iterations of linear search indicated by the user under
key word ITER_LINE_MAXI of the key word factor NEWTON (the default value 0 inhibits
linear search, and is worth 1 then),
·
on criterion RESI_LINE_RELA given by F (). F ()
0, where is worth by defect 0.1.

Linear search is to some extent a “policy insurance” allowing of
to secure against serious divergences of the method of NEWTON. When direction of
seek (un+1, n+1
I
I
) is “bad” (if the tangent matrix is too flexible, by
example), the linear algorithm of search leads to a low value of, which avoids
of going “in the decoration”. It is not necessary to do many iterations in the method
of secant (2 or 3 are enough to avoid the catastrophes) because each one is rather expensive (it
is necessary to recompute the internal forces) and there is not the ambition to find with each iteration of
NEWTON the really optimal one.

2.4 Algorithm
of
STAT_NON_LINE

One will use like previously index I (like “moment”) to note the number of an increment of
charge and exposing it N (like “newton”) to note the number of the iteration of NEWTON in progress.
The algorithm used in operator STAT_NON_LINE can then be written way schematically
following:

(u0,0) and 0 known

Loop over moments Ti (or increments of load): loading L = L
I
(Ti)

·
(ui-1, i-1) known
·
Prediction: calculation of u0
0
I and I
·
Loop on iterations of NEWTON: calculation of a continuation (one, N
I
I)
-
(N
N
N
I, I) and (U,
I
I) known
-
Calculation of the matrices and vectors associated with the following loads
-
Expression of the relation of behavior
-
calculation of constraints N
N
I and of the variables intern I starting from the values
i-1 and i-1 with preceding balance (Ti-1) and of the increment of displacement
one = one - U
I
I
i-1 since this balance
-
calculation of the “nodal forces”: Qn + BT N
I
I
-
possible calculation of the matrix of tangent stiffness: K N = K (one
I
I)
-
Calculation of the direction of search (un+1, n+1
I
I
) by resolution of a system
linear
-
Iterations of linear search:
-
Updating of the variables and their increments:
un+


1 = one + un+1

un+1

= one + un+1
I
I
I
I
I
I

and


n+1
N
n+1
n+1
N
n+


1
I
= I + I

I
= I + I
-
Test of convergence
Handbook of Référence
R5.03 booklet: Nonlinear mechanics
HT-66/05/002/A

Code_Aster ®
Version
7.4
Titrate:
Quasi-static nonlinear algorithm


Date:
06/07/05
Author (S):
P. BADEL, J. LAVERNE, NR. TARDIEU Key
:
R5.03.01-D Page
: 24/28

·
Filing of the results at the moment Ti

U
= U
+ U
I
i-1
I


I = I +
- 1
I



I


I

It is noticed that there are three overlapping levels of iterations: a loop external on the steps of time,
a loop of iterations (known as total) of NEWTON and subloops possible for
linear search (if it is asked by the user) and certain relations of behavior
requiring iterations (known as interns), for example for elastoplasticity in plane constraints.

3 Control

One will refer to documentation [R5.03.80].

4 Large
deformations

4.1 Objective

Until now, we made the assumption that displacements and deformations remained moderate
so as to respect the assumption of the small disturbances. This assumption becomes null and void in
many cases and the method of resolution of the problem must then integrate the evolution of
geometry of the problem, to handle a particular kinematics and to use an adequate formulation of
the law of behavior.

In practice, the assumption of the small deformations can be applied as long as the square of
modulus of deformation remains lower than the precision of calculations considered. In the same way,
the assumption of small rotations can be applied as long as the product enters the square of
the swing angle and the modulus of deformation remain lower than the precision of calculations
considered.

Various alternatives exist within Code_Aster; our objective is not here to do one of them
detailed presentation and we return to the various documents treating specifically each
problems:

·
Relation of nonlinear elastic behavior in great displacements: [R5.03.20],
·
Beams in great displacements: [R5.03.40],
·
Voluminal elements of hulls into nonlinear geometrical: [R5.07.05],
·
Elasto (visco) plasticity, metallurgy and great deformations: [R4.04.03],
·
Elastoplastic modeling with isotropic work hardening in great deformations: [R5.03.21].

The objective is here to present a general functionality of the code which allows a simple processing of
problems great deformations: argument PETIT_REAC of key word DEFORMATION under
key word factor COMP_INCR.
In the continuation, we will limit ourselves to the case of the great plastic deformations, which allows of good
to include/understand problems PETIT_REAC. We base ourselves for that on [bib4] and on
documentation [R5.03.21].
Handbook of Référence
R5.03 booklet: Nonlinear mechanics
HT-66/05/002/A

Code_Aster ®
Version
7.4
Titrate:
Quasi-static nonlinear algorithm


Date:
06/07/05
Author (S):
P. BADEL, J. LAVERNE, NR. TARDIEU Key
:
R5.03.01-D Page
: 25/28

4.2
Great plastic deformations

Initially let us point out the requirements of a modeling great plastic deformations
in mechanics of the continuous mediums in term of equilibrium equations, of kinematic description and of
relation of behavior so as to determine the limitations of modeling PETIT_REAC well.

Equilibrium equations

If one makes the choice write the equilibrium equations on the current configuration and use it
tensor of constraints of Cauchy, they are summarized with:
: v
D = F .v
D +
.
T v
D
v

Vad


éq
4.2-1




The preceding equation can be written on the configuration of reference; it is besides what is made
for nonlinear elasticity in great deformations [R5.03.20]. In the case of laws of behavior
of incremental nature i.e whose evolution is controlled by the current state, this writing loses of sound
interest. One prefers then the writing [éq 4.2-1], although the current configuration is an unknown factor of
problem.

Kinematic description

For the laws of plastic behavior, one uses commonly like measures deformation it
rate of deformation D:
1
D = (U
T
& + U
2
&)

Moreover, one defines also the rate of rotation W:

1
W = (U
T
& - U
2
&)

Let us stress that in the preceding expressions, the operator gradient is defined on the configuration
current.

Elastoplastic relation of behavior

Before writing the relations of behavior, we will make three simplifying assumptions:

·
the elastic strain are small in front of the plastic deformations (what is checked
generally well in the case of metals). This makes it possible to break up the rate of
deformation in an additive way in a plastic part and an elastic part:

D = Of + Dp

·
the plastic deformation is isochoric. This makes it possible to write that the tensors of constraints of
Kirchhoff (stress measurement adequate to the great deformations) and of Cauchy
(intervening in the equilibrium equations) are identical to the first command.
·
the behavior is isotropic.
Handbook of Référence
R5.03 booklet: Nonlinear mechanics
HT-66/05/002/A

Code_Aster ®
Version
7.4
Titrate:
Quasi-static nonlinear algorithm


Date:
06/07/05
Author (S):
P. BADEL, J. LAVERNE, NR. TARDIEU Key
:
R5.03.01-D Page
: 26/28

Under these assumptions, the laws of behavior in great deformations can be written under one
form near to the small transformations. Thus for the plasticity of Von-Mises with isotropic work hardening
one obtains:


J =
E
p

A: (D - D)

F (, p) = - R (p)




éq 4.2-2
eq

3
~

Dp = p

&
2

eq
where

·
~
is the diverter of,
·
A is the tensor of elasticity,
·
J indicates it derived from Jaumann of, which takes account of terms of transport defined in
assistance of the rate of rotation W and which is given by:

J = & - W. +. W

4.3 Functionality
PETIT_REAC

The principle of modeling PETIT_REAC simply consists in reactualizing the geometry of
problem during iterations of Newton (and not at the end of each step of time). This means
that all the quantities intervening in the equations of the problem are evaluated on the configuration
current. Anything else is not modified compared to the case small disturbances.
We have just seen which ingredients are necessary to a “clean” integration of a law of
elastoplastic behavior in great deformations, with the help of three simplifying assumptions.
Now let us detail the differences which the resolution by PETIT_REAC implies.

Equilibrium equations

Taking into account the geometrical reactualization, the interior efforts, first term of [éq 4.2-1], are
estimated in an exact way, with the help of the approximation by finite elements. The calculation of the efforts
outsides east him independent of the resolution by PETIT_REAC: the dead loads are calculated
on the configuration of reference and the loads following on the current configuration. Equations
with balance are thus dealt in an exact way.

Kinematic description

Kinematic description is the same one as that of the small disturbances. This means that one
increment of deformation is calculated by:

1
= ((U) + T (U))



2

The total deflection is then the sum of each one of these increments of linearized deformation,
calculated on different configurations. It is thus delicate to give him a physical direction and better
is worth to use it like an indicator of the level of deformation reached.

Elastoplastic relation of behavior

The difference between the formulation [eq 4.2-2] and the formulation small disturbances lies in
replacement of derived from Jaumann by the simple derivative in time. The disadvantage of this
last is due to the fact that it is not incrémentalement objective i.e not invariant by rigid rotation
structure (for more details, to see documentation [R5.03.21]).
The use of PETIT_REAC is thus not appropriate to great rotations but it is it with large
deformations.
Handbook of Référence
R5.03 booklet: Nonlinear mechanics
HT-66/05/002/A

Code_Aster ®
Version
7.4
Titrate:
Quasi-static nonlinear algorithm


Date:
06/07/05
Author (S):
P. BADEL, J. LAVERNE, NR. TARDIEU Key
:
R5.03.01-D Page
: 27/28

Integration within Code_Aster

In term of finite elements, the resolution by PETIT_REAC implies with each step of load
resolution of the same nonlinear system as previously:

R

(U, T) + BT
= L (T)
I I
I
I


Drunk
= ud (T)


I
I

With the difference close the forces intern are formally calculated by:

R (U, T)
T
= Q (U):
I
I
I

where the operator Q depends on displacements.

Within this framework, the calculation of the tangent matrix carries out formally to:

R

Q (U)
K N =
= Q (U):
+
:
I
U


N
U
N
U
(U, T)
(U, T)
(one, T
I
I
I
I
I
I)


The first term is the contribution of the behavior, similar to what was presented into small
transformations, with the difference which this contribution is evaluated here in current configuration.
second term is the contribution of the geometry which is not present in small transformations.
Within the framework of resolution PETIT_REAC, this term is not present in the calculation of the matrix
tangent. One thus has:


K N = Q (U):
I
U (one, tii)

One can announce that the absence of the geometrical contribution in the tangent matrix can sometimes
to make convergence difficult. Moreover, the presence of important plastic deformations can
to make the problem quasi-incompressible, because of isochoric nature of these deformations.
problems of numerical blocking can then appear which can be circumvented by employment
elements under integrated
Handbook of Référence
R5.03 booklet: Nonlinear mechanics
HT-66/05/002/A

Code_Aster ®
Version
7.4
Titrate:
Quasi-static nonlinear algorithm


Date:
06/07/05
Author (S):
P. BADEL, J. LAVERNE, NR. TARDIEU Key
:
R5.03.01-D Page
: 28/28

5 Bibliography

[1]
P. MIALON: Elements of analysis and numerical resolution of the relations of elastoplasticity.
EDF - Bulletin of Direction of Etudes and Recherches - Série C - N° 3 1986, p. 57 - 89.
[2]
J.F. MAITRE: Numerical analysis, duplicated course of the ENTPE.
[3]
J. SHI & Mr. A. CRISFIELD: “Combining arc-length control and line searches in path
following ", Communications in Numerical Methods in Engineering, Vol. 11, 793-803 (1995).
[4]
E. LORENTZ: Great plastic deformations, modeling in Aster by PETIT_REAC,
EDF - CR MN 1536/07.

Handbook of Référence
R5.03 booklet: Nonlinear mechanics
HT-66/05/002/A

Outline document