Code_Aster ®
Version
5.0
Titrate:
Modelings THHM. General information and algorithms
Date:
22/06/01
Author (S):
C. CHAVANT
Key:
R7.01.10-A
Page:
1/36
Organization (S): EDF/MTI/MN
Handbook of Référence
R7.01 booklet:
Document: R7.01.10
Modelings THHM. General information and algorithms
Summary
Modules THM of Code_Aster are those which treat the equations of the mechanics of the continuous mediums in
using the theory of the porous environments possibly unsaturated and by considering that phenomena
mechanics, thermics and hydraulics are completely coupled. We present the equations here
of balance, or conservation equations solved by these modules. We give a definition of
generalized constraints and generalized deformations, allowing to define way rather general what is
law of behavior THM - at least what the modules considered consider thus - and allowing
to treat the nonlinear equations exhibées within the framework of the algorithms of operator STAT_NON_LINE.
Handbook of Référence
R7.01 booklet:
HI-75/01/001/A

Code_Aster ®
Version
5.0
Titrate:
Modelings THHM. General information and algorithms
Date:
22/06/01
Author (S):
C. CHAVANT
Key:
R7.01.10-A
Page:
2/36
Contents
1 Introduction ............................................................................................................................................ 4
2 Presentation of the problem: Assumptions, Notations ............................................................................... 4
2.1 Description of the porous environment ........................................................................................................... 4
2.2 Notations ......................................................................................................................................... 4
2.2.1 Descriptive variables of the medium ............................................................................................. 5
2.2.1.1 Geometrical variables ............................................................................................. 5
2.2.1.2 Thermodynamic variables of state ........................................................................... 5
2.2.1.3 Descriptive fields of the medium ..................................................................................... 5
2.2.2 Sizes .............................................................................................................................. 6
2.2.2.1 sizes characteristic of the heterogeneous medium ...................................................... 6
2.2.2.2 mechanical magnitudes .............................................................................................. 7
2.2.2.3 hydraulic sizes .............................................................................................. 7
2.2.2.4 Thermal quantities ............................................................................................... 8
2.2.3 External data .............................................................................................................. 9
2.3 Derivative particulate, densities voluminal and mass ............................................................. 9
3 continuous Equations ............................................................................................................................. 10
3.1 Mechanics: conservation of the momentum .............................................................. 10
3.2 Hydraulics: conservation of the mass ........................................................................................ 11
3.3 Equation of energy ..................................................................................................................... 11
3.3.1 The first principle ............................................................................................................... 12
3.3.2 The second principle ........................................................................................................... 12
3.3.3 Equation of energy ............................................................................................................ 13
4 variational Writing of the equilibrium equations .................................................................................. 14
4.1 Mechanics ..................................................................................................................................... 14
4.2 Hydraulics .................................................................................................................................... 14
4.3 Thermics ..................................................................................................................................... 14
5 Discretization in time ........................................................................................................................ 15
5.1 Mechanics ..................................................................................................................................... 15
5.2 Hydraulics .................................................................................................................................... 15
5.3 Thermics ..................................................................................................................................... 16
6 Principle of virtual work, strains and stresses generalized, laws of behavior .......... 16
6.1 Generalized constraints and deformations ..................................................................................... 16
6.2 Principle of virtual work .......................................................................................................... 17
6.3 Laws of behavior ................................................................................................................... 18
6.3.1 Mechanical law of behavior ......................................................................................... 19
6.3.1.1 General writing ...................................................................................................... 19
6.3.1.2 Case of the effective constraints ................................................................................. 19
6.3.1.3 Choice of the constraints .............................................................................................. 19
Handbook of Référence
R7.01 booklet:
HI-75/01/001/A

Code_Aster ®
Version
5.0
Titrate:
Modelings THHM. General information and algorithms
Date:
22/06/01
Author (S):
C. CHAVANT
Key:
R7.01.10-A
Page:
3/36
6.3.2 Hydraulics ........................................................................................................................... 20
6.3.3 Thermal law of behavior ........................................................................................... 20
6.3.4 Homogenized density ......................................................................................... 20
7 Algorithm of resolution ....................................................................................................................... 21
7.1 Nonlinear algorithm of resolution of the equilibrium equations .................................................... 21
7.2 Passage of the nodal values to the values at the points of Gauss ................................................... 21
7.3 Vectors and matrices according to options' .......................................................................................... 23
7.3.1 Residue or nodal force: options RAPH_MECA and FULL_MECA ............................................. 24
7.3.2 Tangent operator: options FULL_MECA, RIGI_MECA_TANG ............................................. 25
7.4 Total algorithm ........................................................................................................................... 29
8 option FORC_NODA ............................................................................................................................. 30
9 space Discretization ........................................................................................................................... 31
10 Bibliography ....................................................................................................................................... 32
Appendix 1 Problème mono dimensional P1P1 ....................................................................................... 33
Handbook of Référence
R7.01 booklet:
HI-75/01/001/A

Code_Aster ®
Version
5.0
Titrate:
Modelings THHM. General information and algorithms
Date:
22/06/01
Author (S):
C. CHAVANT
Key:
R7.01.10-A
Page:
4/36
1 Introduction
Modules THM of Code_Aster are those which treat the equations of the mechanics of the mediums
continuous by using the theory of the porous environments possibly unsaturated and by considering that them
mechanical, thermal and hydraulic phenomena are completely coupled.
We present here the equilibrium equations, or conservation equations solved by these modules.
We give a definition of the generalized constraints and deformations generalized, allowing
to define way rather general what of behavior THM - at least is a law what modules
considered consider thus - and allowing to treat the nonlinear equations exhibées in
tally of the algorithms of operator STAT_NON_LINE.
The laws of behavior THM strictly speaking are not developed in this document, but
in the document [R7.01.11].
Phenomena chemical (transformations of the components, reactions producing of the components
etc…), just as the radiological phenomena are not taken into account at this stage of
development of Code_Aster. The mechanical, hydraulic and thermal phenomena are taken in
count or not according to the behavior called upon by the user in command STAT_NON_LINE,
according to the following nomenclature:
Modeling
Phenomena taken into account
KIT_HM
Mechanics, hydraulics with an unknown pressure
KIT_HHM
Mechanics, hydraulics with two unknown pressures
KIT_THH
Thermics, hydraulics with two unknown pressures
KIT_THM
Thermics, mechanics, hydraulics with an unknown pressure
KIT_THHM
Thermics, mechanics, hydraulics with two unknown pressures
The document present describes the laws of conservation for the most general case said THHM. Cases more
simple are obtained starting from the general case by simply cancelling the quantity absent.
2
Presentation of the problem: Assumptions, Notations
In this chapter, one mainly endeavors to show the porous environment and his characteristics.
2.1
Description of the porous environment
The porous environment considered is a volume made up of a more or less homogeneous solid matrix, more
or less coherent (very coherent in the case of concrete, little in the case of sand). Between
solid elements, pores are found. One distinguishes the closed pores which do not exchange anything with theirs
neighbors and the connected pores in which the exchanges are numerous. When one speaks about
porosity, it is many connected pores about which one speaks. Inside these pores are at the maximum
two components present at more under two phases. The system is regarded as closed.
2.2 Notations
The sizes associated with a component C present under a phase p are noted
p
X. The index of
C
component C can vary from 1 to 2 and that of the phase also.
The porous environment at the current moment is noted, its border
, and it is noted,
at the initial moment.
0
0
N
N


D
indicate the normal in a point of
, image of the normal with
. We will note (
)
0
0
(resp D (

) the element of surface of
(resp
).
0)
0
Handbook of Référence
R7.01 booklet:
HI-75/01/001/A

Code_Aster ®
Version
5.0
Titrate:
Modelings THHM. General information and algorithms
Date:
22/06/01
Author (S):
C. CHAVANT
Key:
R7.01.10-A
Page:
5/36
The medium is defined by:
· parameters (vector position X, time T),
· variables (displacements, pressures, temperature),
· intrinsic sizes (forced and mass deformations, contributions, heat,
hydraulic enthali, flows, thermics…).
For the solid phase, one makes the assumption of small displacements.
The various notations are clarified hereafter.
2.2.1 Descriptive variables of the medium
These are the variables whose knowledge according to time and of the place make it possible to know
completely the state of the medium. These variables break up into two categories:
· variables
geometrical,
· variables of thermodynamic state.
2.2.1.1 Variables
geometrical
In all that follows, one adopts a Lagrangian representation compared to the skeleton (within the meaning of
[bib1]) and the co-ordinates X = X (T are those of a material point attached to the skeleton. All them
S
)
space operators of derivation are defined compared to these co-ordinates.
U

X
Displacements of the skeleton are noted U (X, T) = U.
y

uz
2.2.1.2 Thermodynamic variables of state
The thermodynamic variables are:
· pressures of the components: since we consider that there is to more both
components, there will be with more the two conservation equations of the mass, and thus by duality
with more the two variables of pressure,
· the temperature of the medium T (,
X T).
2.2.1.3 Descriptive fields of the medium
The principal unknown factors, which are also the nodal unknown factors (noted U (,
X T) in this document)
are:
· 2 or 3 (according to the dimension of space) displacements U (,
X T), U (,
X T), U (,
X T) for
X
y
Z

modelings KIT_HM, KIT_HHM, KIT_THM, KIT_THHM,
·
temperature
T (,
X T) for modelings KIT_THH, KIT_THM, KIT_THHM,
· two
pressures
p (,
X T), p (,
X T) for modelings
1
2
KIT_HHM, KIT_THH, KIT_THHM,
· one
pressure
p (,
X T) for modelings
1
KIT_HM, KIT_THM.
Handbook of Référence
R7.01 booklet:
HI-75/01/001/A

Code_Aster ®
Version
5.0
Titrate:
Modelings THHM. General information and algorithms
Date:
22/06/01
Author (S):
C. CHAVANT
Key:
R7.01.10-A
Page:
6/36
2.2.2 Sizes
The equilibrium equations are:
· conservation of the momentum for mechanics,
· conservation of the masses of fluid for hydraulics,
· conservation of energy for thermics.
The equilibrium equations utilize directly the generalized constraints. Constraints
generalized are connected to the deformations generalized by the laws of behavior.
generalized deformations are calculated directly starting from the variables of state and their gradients
space temporal.
The laws of behavior can use additional quantities, often arranged in the variables
interns. These quantities are not described in this document which does not treat laws of
behavior strictly speaking.
2.2.2.1 sizes characteristic of the heterogeneous medium
· Porosity eulérienne: .

If one notes the part of volume occupied by the vacuums in the current configuration,

one a:


=

éq 2.2.2.1-1



The definition of porosity is thus that of porosity eulérienne.
· The saturation of the phase p: p
S.

If one notes
p
the total volume occupied by the phase p, in the current configuration, one has by
definition:

p

S p =
éq 2.2.2.1-2



This saturation is thus finally a proportion varying between 0 and 1. In the equations of
assessment, it is clear that it is the product of porosity by saturation
p
S
who will intervene. One
can thus legitimately wonder why it is not this quantity which is taken as
unknown factor. The answer comes from what it is saturation p
S which intervenes more simply
in the laws of behavior.
· Density eulérienne of the component C in the phase p: p
.
C

If one notes
p
M mass of the phase p of the component C, in volume of the skeleton in
C
the current configuration, one has by definition:

p
Mc = p
D p
p
C
=
PS p
C

D
=
PS p
C

D







éq 2.2.2.1-3
Handbook of Référence
R7.01 booklet:
HI-75/01/001/A

Code_Aster ®
Version
5.0
Titrate:
Modelings THHM. General information and algorithms
Date:
22/06/01
Author (S):
C. CHAVANT
Key:
R7.01.10-A
Page:
7/36


The density of the phase p is simply the sum of the densities of its
components:
p = PC

C

· Lagrangian homogenized density: r.
At the current moment, the mass of volume, M is given by:
M =
rd


0
0
éq 2.2.2.1-4
2.2.2.2 sizes
mechanics
1
· The tensor of the deformations U
() (X,) = (u+T
T
U),
2
· The tensor of the constraints which are exerted on the porous environment: .
This tensor breaks up into a tensor of the effective constraints plus a tensor of constraints of
pressure =
+
'
1.

and
'

are components of the generalized constraints. This cutting
p
p
is finally rather arbitrary, but corresponds all the same to an assumption rather commonly
allowed, at least for the mediums saturated with liquid.
2.2.2.3 sizes
hydraulics
· Mass contributions in components p
m (unit: kilogram per cubic meter). They
C
represent the mass of fluid brought between the initial and current moments. They form part
generalized constraints.

p
p
p
p
p
m = J S
- S

0
C
C
C 0 0
éq 2.2.2.3-1


The mass contributions make it possible to define the total density seen compared to
configuration of reference: R = R
, where R indicates the density
0 + m
+ m + m
lq
vp
have
0
homogenized in an initial state.

· Hydraulic flows:

p
W (unit: kilogram/second/square meter) of representation eulérienne
C

p
M (unit: kilogram/second/square meter) of Lagrangian representation
C


P is noted
v the speed of the component C in the phase p, J Jacobien of the material transformation
C
and the v speed of the skeleton. p
,
,
p
S indicate the densities, porosity and them
S
0
C 0
0
saturations at the initial moment. By definition:

p
p
p
W = S

v - v
C
C
(PC S)

éq 2.2.2.3-2
Handbook of Référence
R7.01 booklet:
HI-75/01/001/A

Code_Aster ®
Version
5.0
Titrate:
Modelings THHM. General information and algorithms
Date:
22/06/01
Author (S):
C. CHAVANT
Key:
R7.01.10-A
Page:
8/36


The Lagrangian form of
p
W noted
p
M is obtained while writing:
C
C

p
M N
. D
W N
.
C
0 (
0) = p D
C
(
)

éq 2.2.2.3-3


The variables m, M and m, M refer each one to a component of conservative mass.
1
1
2
2

One poses by principle:

1
2
1
2
m = m + m
;
M = M + M
1
1
1
1
1
1
1
2
1
2
m = m + m
;
M = M + M
2
2
2
2
2
2


What we will write:

m
m
component =
phase
component
Nb phasedu
component
M
M
component =
phase
component
Nb phasedu
component


In the applications, one could for example have:


2 components: air and water

2 phases for water

1 phase for the air


One would have then: 1
1
m and
M: contribution of mass and liquid water flow
1
1
2
2
m and
M: contribution of mass and vapor flow
1
1
1
1
m and
M: contribution of mass and flow of dry air
2
2
2
2
m and M: non-existent
2
2

· Pressures:
Since we consider that there can be two components other than the solid, there are two
conservation equations of the mass, and thus two multipliers associated, i.e. two
pressures p and p. Aucune assumption is not made on what these two pressures p and p. mean.
1
2
1
2
That will depend on the laws of behavior and the way of writing them. For example one can choose:
p =
(p (gas)

thin cable
pressure
-
)
p (liquid)
1
p pressure
=
ga

of
Z
air)

+
(vapor

2
2.2.2.4 Sizes
thermics
· Not convectée heat Q (see further) (unit: Joule),
p
· Mass enthali of the components m
H C (unit: Joule/Kelvin/kilogram),
· Heat flow: Q (unit: J/S/square meter).
Handbook of Référence
R7.01 booklet:
HI-75/01/001/A

Code_Aster ®
Version
5.0
Titrate:
Modelings THHM. General information and algorithms
Date:
22/06/01
Author (S):
C. CHAVANT
Key:
R7.01.10-A
Page:
9/36
2.2.3 Data
external
· The mass force m
F (in practice gravity),
· Heat sources,
· Boundary conditions relating either to variables imposed, or on imposed flows.
2.3
Derivative particulate, densities voluminal and mass
Description that we make of the medium is Lagrangian compared to the skeleton. One will find in
[bib1] a definition of the concept of skeleton: “the matrix (left occluded solide+porosity) component
the material part of the skeleton and the connected porous space of elementary volume in question
the material point of the skeleton or particle of the skeleton constitutes “.
Either has an unspecified field on, or X (T the punctual coordinate attached to the skeleton that
S
)
we follow in his movement and is X (T the punctual coordinate attached to the fluid. One notes
fl
)
D its
= dt has
! the temporal derivative in the movement of the skeleton:
D its
(X has +, + - X
,
S (T
T) T
T) has ((T) T)
has
S
=
= lim
dt
T
0
T

!
da
has
! is called particulate and often noted derivative (for example in [bib1]). We prefer
dt
to use a notation which recalls that the configuration used to locate a particle is that of
skeleton by report/ratio to which a particle of fluid has a relative speed. For a particle of fluid
location X (T is unspecified, i.e. that the particle of fluid which occupies position X (T with
S
)
S
)
the moment T is not the same one as that which occupies position X (T at another moment T.
S
)
That is to say then A =
AD a quantity related to a voluminal density has, which density is it even

p
range partly by the solid matter constituents and the fluids. That is to say m
has C the mass density of carried by
the liquid phase p of the component C and is with the voluminal density of bound to the solid matter constituents. All
S
these definitions finally amount writing:


p
With =
D =A + A has =
D has +
D has =
has
p
+
S p
amndt


éq 2.3-1




C D
S
fl
S
fl
S
C




p, C

D fl A
While following [bib1], we note
fl the derivative of A if we follow in the movement of
dt
fl
D S.A.
fluid and
S the derivative of A if we follow in the movement of the skeleton.
dt
S
We define then:
DA
D S.A.
D fl A
S
fl
=
S +
fl = D
D

p
p
p m
D has
S +

Cs has Cd
Dt
dt
dt
dt
dt p, C
éq 2.3-2
Handbook of Référence
R7.01 booklet:
HI-75/01/001/A

Code_Aster ®
Version
5.0
Titrate:
Modelings THHM. General information and algorithms
Date:
22/06/01
Author (S):
C. CHAVANT
Key:
R7.01.10-A
Page:
10/36
p
The density m
C has is transported with a relative speed of (p
v - v compared to the skeleton.
C
S)
D its
Taking into account the definition of A =
p
p
p
p
W = S

v - v, one see easily
dt
how the total derivative of A per R!
, and of the definition C
C
(C S)
contribution at time is written finally:
DA
p

=
m
p

+ Div has (cwc has)

D
Dt
!

p, C

éq 2.3-3
Note:
Insofar as we made the assumption of small displacements of the skeleton,
D its
has

has =
and v can be
dt
! can merge with the partial derivative compared to time T S
regarded as null. In the same way, in the continuation of the note we will confuse them
Lagrangian representations and eulériennes of flows,
p
M and p
W.
C
C
3 Equations
continuous
3.1
Mechanics: conservation of the momentum
We note the tensor of the constraints of Cauchy and S the second tensor (symmetrical) of Piola
Kirchoff.
We note P the gradient of the transformation X = X 0 X X,
0
S ()
(T
S
0
)
X
(X, 0)
P =
T
S
X
0
One a:
- 1
- T
S = det P P
.
P
.
.
The equilibrium equations mechanical are written in the configuration:
0
Div P S +
m
RF =
0 (. )
0
We noted Div the operator of divergence compared to the variables of space X of
0
0
configuration.
0
Insofar as we make the assumption of small displacements and the small deformations of
skeleton, this equation can be approached by:
Div +
m
RF = 0
éq 3.1-1
We will further see we always adopt the decomposition =
+ I, where indicates
p
effective constraint. It is thus with the load of the module of integration of the equilibrium equations to make
summon: =
+ I.
p
Handbook of Référence
R7.01 booklet:
HI-75/01/001/A

Code_Aster ®
Version
5.0
Titrate:
Modelings THHM. General information and algorithms
Date:
22/06/01
Author (S):
C. CHAVANT
Key:
R7.01.10-A
Page:
11/36
3.2
Hydraulics: conservation of the mass
The writing eulérienne of the conservation of the fluid mass for the component C is written:
D fl pcSp
D = 0
dt p
p
One can then apply [éq 2.3-1] while taking: has
and m
C = 1 has
S = 0
and [éq 2.3-3] will give:
S
D p S
C

p + Div (p
WC) = 0
dt
p
p
By using the definition of the mass contributions [éq 2.2.2.3-3], the definition of Lagrangian flows
[éq 2.2.2.3-2] one finds the form Lagrangian of the conservation of the fluid mass:
m Div M
1 +
0 (
1) = 0

m
Div M
2 +
0 (
2) = 0
!
!
éq 3.2-1
3.3
Equation of energy
For the function thermodynamic, we adopt a decomposition of the type systematically
[éq 2.3-1]. That corresponds to the fact that various energies have a whole a part carried by the solid
and a part carried by the fluids. The part carried by the solid is characterized by a density
voluminal whereas the parts carried by the fluid are characterized by mass densities,
as we showed in the paragraph [§2.3].

p
Total internal energy: E = E
éq 3.3.1
S + GCV pemc D





p, C


p
Total entropy: S =
S
éq 3.3.2
S + GCV psmc




D

p, C


p
Total enthalpy H =
H
éq 3.3.3
S + GCV phmc




D

p, C

= E - TS

Free energy:
E
Ts
éq 3.3.4
S =
S -
S

p
p
p
mc = m
E C - m
Ts C
G = H - TS

Free enthalpy: G
H
Ts
éq 3.3.5
S =
S -
S

p
p
p
m
G C = m
H C - m
Ts C
Handbook of Référence
R7.01 booklet:
HI-75/01/001/A

Code_Aster ®
Version
5.0
Titrate:
Modelings THHM. General information and algorithms
Date:
22/06/01
Author (S):
C. CHAVANT
Key:
R7.01.10-A
Page:
12/36
Lastly, by noting Q ()
! the rate of heat received by a volume, one has by definition:
Q () = qn.
D +
D
!

éq 3.3.6
It is pointed out finally that the enthalpy of the fluids is calculated by the formula:
H = + p
E
éq 3.3.7
3.3.1 The first principle
With the definitions given higher, he is written:
- E - Div (H pmcMp: +
c)+
M p.
Fm
C
+ -
Q
Div = 0
!
!
p, C
p, C
éq 3.3.1-1
This writing corresponds to the equation (22) chapter III-2-3 of [bib1], in which we have
neglected the terms of inertia. For the homogeneous mediums, it corresponds to the equation (31) of
paragraph IV-3-2 of [bib3].
3.3.2 The second principle
Its form rather well known is:
Q
S + Div (S p
m
p
C Mc)

+ Div - 0
! p, C
T T
éq 3.3.2-1
By using the traditional thermodynamic considerations [bib1] related to the introduction of the enthalpy
free [éq 3.3.5], it is shown that one must necessarily have:

-
= 0

éq 3.3.2-2
p
m

G C -
= 0
éq 3.3.2-3
p
mc

S +
= 0
T
éq 3.3.2-4
Handbook of Référence
R7.01 booklet:
HI-75/01/001/A

Code_Aster ®
Version
5.0
Titrate:
Modelings THHM. General information and algorithms
Date:
22/06/01
Author (S):
C. CHAVANT
Key:
R7.01.10-A
Page:
13/36
3.3.3 Equation of energy
Enough often, one considers that, the transformations being reversible, the second principle provides
finally an equality. Moreover, one replaces in [éq 3.3.2-1] the unknown temperature T by one
value constant known as temperature of reference. It is finally about a linearization of [éq 3.3.2-1]
justified if the variations in temperatures are “small”. Let us note that the term of transport
Div (pm p
C M M complicates the processing of nonthe linearity due to the presence of the temperature in
c)
p, C
denominator of the other terms of [éq 3.3.2-1].
We work in enthalpy in order to overcome this difficulty. One leaves the equation of the first principle
[éq 3.3.1-1] into which one injects the equations [éq 3.3.2-2], [éq 3.3.2-3], [éq 3.3.2-4], and the definition
free enthalpy [éq 3.3-5] and one obtain:
Ts + (H p
mc m p - Ts p
mc m p = -
+. + -
éq 3.3.3-1
C
c)
Div (H p
mc M PC)
M p F m
Q
Div
C
!
!
!
p, C
p, C
p, C
One poses then:
Q = Ts - T p
m
C M p
mc
p, C
éq 3.3.3-2
The quantity Q with the dimension of an energy per unit of volume. It represents the heat received by
the system in a transformation for which there are no contributions of heat per input of fluid
having a enthalpy. Although Q
is not an exact differential, we take this quantity
like variable of state.
Finally, the equation of energy retained with the following form:
pm p
H C m
Q
Div H M
Divq
Mr. F
C + +
(pm p
C
c)+
- p m
C
=
! !
p, C
p, C
p, C
éq 3.3.3-3
Handbook of Référence
R7.01 booklet:
HI-75/01/001/A

Code_Aster ®
Version
5.0
Titrate:
Modelings THHM. General information and algorithms
Date:
22/06/01
Author (S):
C. CHAVANT
Key:
R7.01.10-A
Page:
14/36
4
Variational writing of the equilibrium equations
4.1 Mechanics
We note U
the whole of the fields of displacement kinematically acceptable, i.e.
AD
elements of (H () 3
1 checking the boundary conditions in displacement on the part of

supporting such conditions [bib3].
The variational form of [éq 3.1-1] is:

= '
+

p I

éq 4.1-1
.(v) D =
m
RF. D
v +
ext.
F
D
v
v




U

AD

4.2 Hydraulics
We note P (resp. P
) the whole of the acceptable fields of pressure, i.e. them
1ad
2ad
elements of
1
H () checking the boundary conditions in pressure P (resp. P) on the part of

1
2
supporting such conditions [bib3]. The variational form of [éq 3.2-1] is:
- m1 m2 D
1
2
M
Mr.
D
(1 +
1) 1 + (
1 +
1) 1 =


! ! 1 2
M
M
. D
P

(
1ext +
1 ext.) 1

1 1ad
éq 4.2-1
1
2
1
2

- m m
D
M
Mr.
D
(2 +
2) 2 + (
2 +
2) 2

! ! = 1 2
M
M
. D
P

(
2
+
ext.
2ext)

1

2 2ad
4.3 Thermics
We note T the whole of the acceptable fields of temperature, i.e. the elements of
AD
1
H () checking the boundary conditions in temperature on the part of
supporting of such
conditions. [bib3]. The variational form of [éq 3.3.3-3] is:



Q' D
+
M Q.

p
p
m
p
m
p
H C m D
-
H C
+
D


=
C
C




!
!
p, C
p, C





p

p
m
m
p
+
Mr. F
D


-
H C
M +q. D






éq 4.3-1
C
C ext.
ext.




p, C
p, C





Tad
Handbook of Référence
R7.01 booklet:
HI-75/01/001/A

Code_Aster ®
Version
5.0
Titrate:
Modelings THHM. General information and algorithms
Date:
22/06/01
Author (S):
C. CHAVANT
Key:
R7.01.10-A
Page:
15/36
Let us note that, contrary to other presentations, and in particular [bib8] we did not inject them
conservation equations of the mass, and we integrated by part the term of transport
Div (pm p
H C Mr. This last point has the advantage of not revealing of derived from command
c)
p, C
superior, and, contrary to revealing naturally boundary conditions relative to
p
the input of heat related to hydraulic flows:
hm
p
C M

. D
cext


.

p, C
One will be able in makes consider that the conditions of heat flux define directly:
p
m
p
Q
~ = H cm +q
ext.
C ext.
ext.
5
Discretization in time
In this chapter, we are satisfied to take again the variational formulations in their
applying a discretization compared to the time of the type téta diagram. It is about a method
general of integration of the differential equations [bib12] and [bib13].
is a numerical parameter ranging between 0 and 1. For the linear differential equations (what
is not our case…) this diagram is unconditionally stable for 1/2, it is of command 1 for
1/2 and of command 2 for = 1/2. Nevertheless, it can be preferable to use a value different from
1/2, and this for parasitic reasons of oscillations [bib12].
The subscripted quantities by + are the quantities at the end of the step of time, and those subscripted by ­ are
those of the beginning of the step of time. One notes:
+
-
T = T - T
has = has
+ + (1 -) A
has

5.1 Mechanics
+
+
+

= '+ I
p

éq 5.1-1

+
. (v)
+ m+
ext+


D = R F v
.
D +
F
v
D
v

U



AD
5.2 Hydraulics
- (+
m1
M
Mr.
1 +
+
m21) D
1 +


+
+
T
1
2
D

(1 + 1) 1 =


- (-
m1
1
M
Mr.
1 +
-
m21)

1
D - (-)
-
-
T
1
2
D

(1 + 1)

1

+ T (1
2
M
M
. D
P
ext.
1
+

ext.
1
)


1

1


1ad

éq 5.2-1
- (+
m1
M
Mr.
2 +
+
m22) D
2 +


+
+
T
1
2
D

(2 + 2)

2 =


- (-
m1
1
M
Mr.
2 +
-
m22)

2
D - (-)
-
-
T
1
2
D

(2 + 2) 2

+ T (1
2
M
M
. D
P
2ext +

2ext)

2

2



2ad

Handbook of Référence
R7.01 booklet:
HI-75/01/001/A

Code_Aster ®
Version
5.0
Titrate:
Modelings THHM. General information and algorithms
Date:
22/06/01
Author (S):
C. CHAVANT
Key:
R7.01.10-A
Page:
16/36
5.3 Thermics
(+
Q' - -
Q')

p +



D -

m
p+
T
H C M
Q
D
C
+ +







p, C



- (1 -)

p



m
p
T
H C M
Q
D
H
m
m
D
C
+ -
+

p
m
p
p
C


+
+
-

(C - c)










p, C

p, C

+ (1 -) p
m
H C (p+
m
m
D
éq 5.3-1
C
- p
C
)


=


+

p +
m
T
Mr. F D
1
T
Mr. F D
C
+ (-)


p -
m
C




p, C
p, C




p



m
p


+ T
D - T
H C M
Q
. D
T




cext +
ext.







AD

p, C



One can again consider that the conditions of heat flux define directly:

p

~
m
p

Q
H
M
Q
ext. =
C
C
+
ext.
ext.
p, C
6
Principle of virtual work, deformations and forced
generalized, laws of behavior
6.1
Generalized constraints and deformations
While referring to the variational formulations [éq 4.1-1], [éq 4.2-1] and [éq 4.3-1], it appears that one
can choose:
For the generalized constraints:

;


p

1
m, M1, 1
hm
m
1;
2
m, m2, 2
H 1;
= 1 1
1
1

éq 6.1-1
1
m, M1, 1
hm
m
2;
2
m, m2, 2
H2;
2
2
2
2

Q' Q
,

For the generalized deformations:
= {U, (U);p, p
1
; p,
1
2 p T
; ,
2
T}
éq 6.1-2
The fact is noticed that the generalized deformations contain displacements. That is due to
term Fm
R
v
. variational formulation of the conservation equation of the quantity of

movement [éq 4.1-1], which term couples finally the generalized constraints and displacements
because of [éq 6.3.4-1]. The generalized deformations contain the pressure and the temperature parce
that the associated equations are parabolic.
Handbook of Référence
R7.01 booklet:
HI-75/01/001/A

Code_Aster ®
Version
5.0
Titrate:
Modelings THHM. General information and algorithms
Date:
22/06/01
Author (S):
C. CHAVANT
Key:
R7.01.10-A
Page:
17/36
6.2
Principle of virtual work
The whole of the nonlinear equations to solve can be put in the form:
R (U)
meca
= L
éq 6.2-1
where U indicates generalized displacements, i.e.: U = {U U
, U
, p, p, T in the case more
X
y
Z
1
2
}
General. The forces intern R express themselves starting from a principle of virtual work generalized. In
the case of the mechanics of the continuous mediums “traditional”, i.e. when there is no other
component that the solid, one is accustomed to defining the forces intern by:
T
W .R = W
.
D W, field of displacement kinematically acceptable.
()
In this formulation, the field of deformation depends only on the field of displacement and on
its derivative space, possibly in a nonlinear way if one takes into account deformations
finished. One writes symbolically:
R = T
Q
The law of behavior connects the constraints to the deformations.
Within the framework of the theory of the porous environments developed here, we try to bring us closer it
more possible of this formulation by introducing generalized constraints and deformations
generalized E; The generalized deformations depend only on the field of displacements
generalized U and of its derivative space. The operator U

E (U) is an operator of derivation
compared to the field of co-ordinates.
The law of behavior makes it possible to calculate according to E.
On the other hand, we cannot write directly
T
W R
. = W
.

D, for the reasons

()
following:
· the equations which we treat are evolutionary equations in times and the derivative by
report/ratio at the time of the quantities intervene,
· the equations are nonlinear because of the terms of transport related to the representation
eulérienne of the fluids: these nonlinear terms appear only in the equation of thermics,
· the choice of the unknown factors makes that the nonlinear terms of transport intervene in
generalized constraints. That is to say a term of transport in the equation [éq 4.3-1],

p
m
p

H C M
Q.
, since one took as principal unknown factors for
C +

D



p, C

hydraulics pressures, quantities
p
M related to speeds of the fluids belong to
C
p
generalized constraints, just as the enthali m
H C, and the term of transport given in
example is linear in deformation generalized and quadratic in generalized constraint. This
fact a difference with the formulation of the theory of the traditional continuous mediums where
terms of great deformation are quadratic quantities of the deformations.
For all these reasons, we introduce a field noted such as:
R QT
=
éq 6.2-2
T
Q is defined by:
WT .R = E (W).

D

éq 6.2-3
W

cinématiqu

generalized
T
déplacemen

of

field
acceptable
ement
Handbook of Référence
R7.01 booklet:
HI-75/01/001/A

Code_Aster ®
Version
5.0
Titrate:
Modelings THHM. General information and algorithms
Date:
22/06/01
Author (S):
C. CHAVANT
Key:
R7.01.10-A
Page:
18/36
It is very seen easily and classically that T
Q is transposed of the operator Q such as:
E = QU
éq 6.1-4
The field is a linear function of
! and nonlinear of:
= (,)
!
éq 6.1-5
After discretization in time, +
will become a nonlinear function of +
and -
:
+
+
= (+ -
,)
éq 6.1-6
Let us note finally that for algorithmic reasons (inter alia), one needs to know the derivative
forces intern compared to generalized displacements:
R

R


E

T

=
= Q
Q
U



E

U



E


It is clear that
depends only on the form of the equilibrium equations.

6.3
Laws of behavior
A law of behavior will be simply defined like an unspecified relation between constraints
generalized and generalized deformations. The variables intern are defined like fields
necessary to the calculation of the constraints, whose evolution is given by the laws of behavior, but
who do not intervene directly in the equilibrium equations.
Moreover, we consider that the laws of behavior are written in incremental form and
that they are local. By noting the internal variables, a law of behavior is thus one
relation:
, E
,

!

! !
After discretization in time, the law of behavior becomes a relation:
-
-
-
+
+
+
, E, E

,
The law of behavior will have to also provide the only term which in the expression
R


T

= Q
Q depends on it, namely
. Finally a relation of behavior is one
U



E

E

relation:
+
-
-
-
+
+
+
, E, E

,
éq 6.3-1
+
E
In the following paragraphs, we specify certain aspects of the laws of behavior in
distinguishing the mechanical, hydraulic and thermal contributions.
Handbook of Référence
R7.01 booklet:
HI-75/01/001/A

Code_Aster ®
Version
5.0
Titrate:
Modelings THHM. General information and algorithms
Date:
22/06/01
Author (S):
C. CHAVANT
Key:
R7.01.10-A
Page:
19/36
6.3.1 Mechanical law of behavior
6.3.1.1 Writing
general
The variables intern are noted. A law of behavior of mechanics, within framework THM
is written:

+
= +
(+ + + + - - - - - -
, p, p, T;
p p T
1
2
,
1
2
,)

éq 6.3.1.1-1

+
= +
(+ + + + - - - - - -
, p, p, T;
p p T
1
2
,
1
2
,)
6.3.1.2 Case of the effective constraints
In the case of the assumption of the effective constraints, one with the decomposition: =
+ I where is
p
the tensor of the effective constraints and is a scalar.
p
The variables intern are separate in two parts: the variables intern mechanical and them
variables intern hydraulic. The mechanical law of behavior is divided then into two laws,
H
whose first can be an already existing law within the usual framework of thermomechanics.

+
= +
(+ + - - - -
, T; , T,
,)

éq 6.3.1.2-1

+

, T; , T,
,
=
+
(+
+
-
-
-
-
)

+
p, p; p, p,
1
2
1
2
p = +p (+
+
-
-
-
H)

éq 6.3.1.2-2

+

p, p; p, p,
1
2
1
2
H =
+
H (+
+
-
-
-
H)
The dependences indicated by the equations [éq 6.3.1.2-1] and [éq 6.3.1.2-2] do not have a justification
theoretical a priori. It is simply a question of showing the most general possible dependences
point of view of the data-processing programming. One notices in this decomposition that
dependence compared to thermics was left in the effective constraints; typically, one
think that the laws on the effective constraints are written as in traditional thermomechanics:
+
+
= (+ + + - - - - -
- T I, - T I,
,)
6.3.1.3 Choice of the constraints
Because of rather frequent use of the assumption of the effective constraints, one decides that it
vector of the constraints for the mechanical part contains in all the cases the tensor of the constraints
effective and the scalar. In the general case where the assumption of the effective constraints is not
p
reserve, one will have simply:
. It is thus with the load of the module of integration of the equations
p = 0
of balance (and not of the laws of behavior) to make the sum: +
+
+
=
+ I.
p
Handbook of Référence
R7.01 booklet:
HI-75/01/001/A

Code_Aster ®
Version
5.0
Titrate:
Modelings THHM. General information and algorithms
Date:
22/06/01
Author (S):
C. CHAVANT
Key:
R7.01.10-A
Page:
20/36
6.3.2 Hydraulics
The hydraulic law of behavior provides the following relations:
+
p
m
m
p p T p p T m M
C
=
+
p
C
(
-
-
+
+
+
+
-
-
-
-
p
p
_
,
,
,
; ,
,
,
,
,
,
1
2
1
2
C
C
H)



+ +
, p, p, p, p, T, T;
C
1 +
+
1
2 +
+
2
+






+
+

p
M
M
p
p p
p
p
C
=
p -
-
,
,
,
,
;
C
1 -
-
1
2 -


2




-
+
-
-
p
_
m




T, T, M;F
C
H





-
+

p p T p p T m
H =
+
H (+
+
+
+
-
-
-
-
p
_
,
,
,
; ,
,
,
,
,
1
2
1
2
C
H)
éq 6.3.2-1
It is noticed that the field of gravity is a data of the hydraulic law of behavior parce
that the evolution of the vector of flow follows relations of the type:
fl
M = - P + F.
H
[
fl
m]
6.3.3 Thermal law of behavior
The laws of behavior give:
Q
'+ = Q'+ (+
, p+, p+, T +; -
1
2
, p, p, T -, Q'-
1
2
)

p+
p+
p
hm
m
+ + + + - - - - m
C
= H C, p, p, T;

1
2
, p, p, T, H
1
2
C
C and p








éq 6.3.3-1
+
+
Q

= Q (+
, p+, p+, T +, T +
1
2
; -
, p, p, T -, T -
1
2
, -
Q)

+
+
+

=


(, p+, p+, T+, T+; -, p, p, T, T, -

1
2

T)
T
T
1
2
Let us note that we introduced possible internal variables related to thermics.
6.3.4 Homogenized density
By definition, the homogenized density, which intervenes in the equilibrium equation of
mechanics [éq 3.1-1] is given by:
+
+
+
+
+
R = R
m
m
m
m
0 +
1
1
+ 21 + 12 + 22
éq 6.3.4-1
This equation is not a law of behavior, but it belongs to the conservation equations.
It is integrated in the module of calculation of the equilibrium equations, the modules of calculation of the laws of
behavior not having to treat it.
Handbook of Référence
R7.01 booklet:
HI-75/01/001/A

Code_Aster ®
Version
5.0
Titrate:
Modelings THHM. General information and algorithms
Date:
22/06/01
Author (S):
C. CHAVANT
Key:
R7.01.10-A
Page:
21/36
7
Algorithm of resolution
7.1
Nonlinear algorithm of resolution of the equilibrium equations
In the general case of modeling (variable coefficients, desaturation, convection) the problem
variational presented above [éq 4.1-1] with [éq 4.3-1] is nonlinear compared to the fields of
displacement, pressure and temperature. After discretization by finite elements, one obtains a system
matric nonlinear. The matrix of the tangent operator contains moreover one term nonsymmetrical treaty
as tel. One uses in all the cases of modeling nonlinear solvor STAT_NON_LINE of
Code_Aster resting on a method of Newton-Raphson, described in [bib5]. One introduces
vectorial functional calculus:
F (U) = R (U)
meca
- L
éq 7.1-1
F

The associated tangent operator is noted: F
D = U

For the modules THM, objects of this note, the operator meca
L
does not depend on displacements
generalized. All the terms depending on generalized displacements were introduced into R, and
it is precisely for this reason that displacements are found in the deformations
generalized. Let us note on this subject the very particular processing of the term Fm
R
v
. equation [éq 4.1-1].

+
+
+
+
According to [éq 6.3.4-1], R m
F v
.
D =

R
m1
m2
m1
m2
m
F. D
v

(0 + 1 + 1 + 2 + 2)


We chose to divide this term into two:
The term R m
F v
.
is a contribution to meca
L
0

D
if the user informed operand PESANTEUR

loading used (defined by command AFFE_CHAR_MECA), whereas the term
(+
m1
, which depends on the generalized constraints is a contribution to
1 +
+
m21 + +
m12 +
+
m22) m
F v
.
D

R.
7.2
Passage of the nodal values to the values at the points of Gauss
As in all the codes of finite elements, the terms are calculated by loop on the elements and
loop on the points of gauss. While noting
el
R and
el
DF values at the point of gauss G of the element el
G
G
nodal forces and of the tangent operator, and el
W the weight of integration related to this point of gauss, one a:
G
(
R U) = el el
W R U
G
G (
)
el
G
éq 7.2-1
DF (U) = el el
W DF U
G
G (
)
el
G
éq 7.2-2
Handbook of Référence
R7.01 booklet:
HI-75/01/001/A

Code_Aster ®
Version
5.0
Titrate:
Modelings THHM. General information and algorithms
Date:
22/06/01
Author (S):
C. CHAVANT
Key:
R7.01.10-A
Page:
22/36
Let us note el then
U the vector of the nodal unknown factors on a finite element el. One can thus have:
U
v
W
1

node

p1
p
2
T
U
v
W
el

example
by
U =
2

node

p1
p
2
T
U
v
W
3

node

p1
p
2
T
.
Handbook of Référence
R7.01 booklet:
HI-75/01/001/A

Code_Aster ®
Version
5.0
Titrate:
Modelings THHM. General information and algorithms
Date:
22/06/01
Author (S):
C. CHAVANT
Key:
R7.01.10-A
Page:
23/36
Let us note also el
E the vector of the deformations generalized at the point of gauss G of the element el and el

G
G
the vector of constraints generalized for the point of Gauss G of the element el. In the case more
complete one has as follows:
el




p
m1
1
M11
el
U
m


H 1
1
(U)



m21


p
2

1
M1


m
p1
H 21
el
el
E =
;
=
G
G
1
p2
m


2
p

2
M1


21
m
T


H2
T



m2
G
2
M22
m
H2
2
Q'


Q G
The functions of form of the finite elements make it possible to calculate the matrix then
el
Q of passage of
G
nodal unknown factors with the deformations generalized at the points of gauss defined by:
el
el
el
E = Q U
.
éq 7.2-3
G
G
7.3
Vectors and matrices according to options'
The presentations of the two following paragraphs are made in the most general case where one has one
equation of mechanics, two equations of hydraulics and an equation of thermics. The indices G and el
from now on are omitted, but it is clear that what is described applies to each point of gauss of
each element.
Handbook of Référence
R7.01 booklet:
HI-75/01/001/A

Code_Aster ®
Version
5.0
Titrate:
Modelings THHM. General information and algorithms
Date:
22/06/01
Author (S):
C. CHAVANT
Key:
R7.01.10-A
Page:
24/36
7.3.1 Residue or nodal force: options RAPH_MECA and FULL_MECA
One distributes the terms of the variational formulation according to the following principle:
el
el
If
*
E
*
E G = (v, v,
calculated with
1

, 1, 2, 2,)
G indicates a virtual field of deformation,
()
el
to start from a vector of virtual nodal displacements

*
U, one can define
:
T
*
E elg .elg (U) = 1v + 2 (v) + 3 + 4
. They then are taken again
1
+ 5
1
+ 6
2

+ 7
2
+ 8
discrete variational formulations [éq 5.1-1], [éq 5.2-1], [éq 5.3-1], and one replaces the integrals there
fd by el el
W F for all integral F. One distinguishes the terms multiplying

G
G
el
G
respectively v, (v), and, and one finds:
1
1
2
2
Index
associated
1
- (+
+
+
+
v
m11 + m21 + m12 + 22) +m
m
F
2
+
+
+ I
(v)
p
3
+
+
-
-
- 1
m
m
m
m
1
1
- 21 + 11 + 21
4

T (+
+
1
M
M
1
T M
M
1
1
+ 21) + (-) (-
-
1
1
+ 21)
5
+
+
-
-
- 1
m
m
m
m
2
2
- 22 + 12 + 22
6

T (+
+
1
M
M
1
T M
M
2
2
+ 22) + (-) (-
-
1
2
+ 22)
7
- Q'+ +Q'-

+
-
+
-
m
- H 1
+
-
+
-



1 + (1 -) m
H 11 (m1 - m1 -
+ -

-
1
1)
m
H 21
(1) m
H 21 (m2 m2
1
1
)




+
-
+
-
m
- H 1
+
-
+
-



2 + (1 -) m
H 12 (m1 - m1 -
+ -

-
2
2)
m
H 22
(1) m
H 22 (m2 m2
2
2)




+ T
(
M1+ + M2+ + M1+ + M2+
-
-
-
-
+ -
+
+
+
1
1
2
2) m
F
.
T (1
) (M1 m2 M1 m2
1
1
2
2) m
F
.
8

+
+
+
+

+

1
+
m
m
m
m
T H 1 M1
H
H
H
1
+ 2
+
1 m2
1
+ 1
+
2 M12 +
2
+
2 M 22 + +
Q

+


-
-
-
-
+ (1 -)

1
-
m
m
m
m
T H 1 M1
H
H
H
1
+ 2
-
1 m2
1
+ 1
-
2 M12 +
2
-
2 M 22 + -
Q



Note:
In first term 1 does not appear the term
m
- R F because it is put in the loading
0
outside meca
L
and calculated by the option of calculation of the loading external of gravity.
By using the definition [éq 7.2-1] of el
R, one a:
G
T
T
el
* el
el
* el
U
R
.
= E
. G, which still gives us:
G
G
T
el
el
Rel = Q
. G
G
G
This last equality is only the local form on the level of a point of gauss of [éq 6.2-2].
Handbook of Référence
R7.01 booklet:
HI-75/01/001/A

Code_Aster ®
Version
5.0
Titrate:
Modelings THHM. General information and algorithms
Date:
22/06/01
Author (S):
C. CHAVANT
Key:
R7.01.10-A
Page:
25/36
7.3.2 Tangent operator: options FULL_MECA, RIGI_MECA_TANG
In what follows, if X indicates a vector of components X I and Y a vector of components Y J,
X
I
X


a matrix will indicate whose element occupying line I and the column J is
.
Y
J
Y

To calculate the tangent operator, the following quantities will be calculated:
[DRDE] =
DR1U
DR1E
DR1P1
DR1GP1
DR1P2
DR1GP2
DR1T
DR1GT
DR2U
DR. 2nd
DR2P1
DR2GP1
DR2P2
DR2GP2
DR2T
DR2GT
DR3U
DR. 3rd
DR3P1
DR3GP1
DR3P2
DR3GP2
DR3T
DR3GT
DR4U
DR. 4th
DR4P1
DR4GP1
DR4P2
DR4GP2
DR4T
DR4GT
DR5U
DR. 5th
DR5P1
DR5GP1
DR5P2
DR5GP2
DR5T
DR5GT
DR6U
DR. 6th
DR6P1
DR6GP1
DR6P2
DR6GP2
DR6T
DR6GT
DR7U
DR. 7th
DR7P1
DR7GP1
DR7P2
DR7GP2
DR7T
DR7GT
DR8U
DR. 8th
DR8P1
DR8GP1
DR8P2
DR8GP2
DR8T
DR8GT
Where one noted:
F
F
DRiU
I
=
DRiGP
I
1 =
U
p

1
F
F
DRiE
I
=
DRiGP
I
2 =


p2
F
F
DRiP
I
1 =
DRiT
I
=
p

T
1

F
F
DRiP
I
2 =
DRiDT
I
=
p

T
2


Handbook of Référence
R7.01 booklet:
HI-75/01/001/A

Code_Aster ®
Version
5.0
Titrate:
Modelings THHM. General information and algorithms
Date:
22/06/01
Author (S):
C. CHAVANT
Key:
R7.01.10-A
Page:
26/36
To make these calculations one considers that the laws of behavior provide, for the options
corresponding, all the derivative following:
'
'
'
'
'
'
'
'


U

p

1

p1 p2

p2
T
T
p
p
p
p
p
p
p
p
U

p

1

p1 p2

p2
T
T


m1
1
1
1
1
1
1
1
1
m1
m1
m1
m1
m1
m1
m1
U

p

1

p1 p2

p2
T
T


M
1
M1
1
1
1
1
1
1
1

M
1
1 M1 M1 M1 M1 M1
U

p

1
1
p
2
p
2
p
T
T

1
1
1
1
1
1
1
1
1
hm
1
hm
1
hm
1
hm
1
hm
1
hm
1
hm
1
hm


U

1
p
1
p
2
p
2
p
T
T
2
2
2
2
2
2
2
2
1
m
1
m
1
m
1
m
1
m
1
m
1
m
1
m


U



1
p
1
p
2
p
2
p
T
T
M2
M2

M2

M2

M2

M2

M2

M2
1
1
1
1
1
1
1
1


U

1
p
1
p
2
p
2
p
T
T



2
2
2
2
2
2
2
2

1
m
1
m
1
m
1
m
1
m
1
m
1
m
1
H
H
H
H
H
H
H
hm


U

p
p
p
p
T
T

D D
E =
1
1
2
2

1
1
1
1
1
1

1
1
2
m
2
m
2
m
2
m
2
m
2
m
m

m
2
2


U

1
p
1
p
2
p
p
T
2

T

1
1
1
1
1
1
1
1
M

2
M2 M2 M2 M2 M2 M2 M2


U

p
1

p
p
1
2

p
T
2

T
1
1
1
1
1
1
1
1


hm2
hm2 hm2 hm2 hm2
hm2
hm2

hm2



U

p


1

p
p
1
2

p
T
2

T


m2

m2

m2

m2
2
2
2
2
2
2
2
2
2
m
2
m
2
m
2
m
U

p


1
1
p
2
p
2
p
T
T

2
2
2
2
2
2
2
2
M2 M2 M2 M2 M2
M2 M2 M2
U



1
p
1
p
2
p
2
p
T
T

2
2
2
2
2
2
2
2
2
m
2
m
2
m
2
m
2
m
2
m
2
m
2
H
H
H
H
H
H
H
hm


U



1
p
1
p
2
p
2
p
T
T
Q'
Q'
Q'
Q'
Q'
Q'
Q'
Q'


U

p1 p1 p2 p2
T
T



Q

Q

Q

Q

Q

Q
Q
Q
U

p


1
p1 p2
p2
T
T

Note:
In these expressions, the derivative compared to U are all null, but we keep
the writing taking into account the definition of the matrices
el
Q which we adopted.
G
Handbook of Référence
R7.01 booklet:
HI-75/01/001/A

Code_Aster ®
Version
5.0
Titrate:
Modelings THHM. General information and algorithms
Date:
22/06/01
Author (S):
C. CHAVANT
Key:
R7.01.10-A
Page:
27/36
The call to the laws of behavior will provide the pieces of the matrix D OF

according to the equations
present:




















[
]
DMECDE =
p
p
p
p
;
[
]
DMECP1 = 1

1;
[DMECP2] = 2

2 [DMECDT]

= T

T
p
p
p
p
p
p
p


p
p
p
p

1

1
2




2
T

T
1

m
m11
m1

1
m11
m1
1

1
m11
m1






1
p
p
p
p



1

1
2

2

1
1
1
1
1
T

T
1
1
[DP11DE] M
M
M
M
M
M
M
1
=
;
[DP11P1] 1
=
1;
[DP11P2] 1
=
1 [DP11DT] 1 1
=

p1

p1
p2

p

2
T

T
1




m
1
1
m1
m1
m1
m1
m
m
H 1
H 1. 1
H 1. 1


H 1. 1





p1

p1
p2

p
T

T

2
2

m
m21
m2

1
m21
m2
1

1
m21
m2






1
p
p
p
p



1

1
2

2

2
2
2
2
2
T

T
2
2
[DP12DE] M
M
M
M
M
M
M
1
=
;
[DP12P1] 1
=
1;
[DP12P2] 1
=
1 [DP12DT] 1 1
=

p1

p1
p2

p

2
T

T
2




m
2
2
m2
m2
m2
m2
m
m
H 1
H 1. 1
H 1. 1


H 1. 1





p1

p1
p2

p
T

T

2
1

m
m12
m1

2
m12
m1
2

2
m12
m1






2
p
p
p
p



1

1
2

2

1
1
1
1
1
T

T
1
1
[DP21DE] M
M
M
M
M
M
M
2
=
;
[
] 2
DP21P1 =
2;
[DP21P2] 2
=
2 [DP21DT] 2
2
=

p1

p1
p2

p

2
T

T
1




m
1
1
m1
m1
m1
m1
m
m

H2
H2 H2
H2 H2


H2 H2 m





p1

p1
p2

p
T

T

2
2

m
m22
m2

2
m22
m2
2

2
m22
m2






2
p
p
p
p



1

1
2

2

2
2
2
2
2
T

T
2
2
[DP22DE] M
M
M
M
M
M
M
2
=
;
[DP22P1] 2
=
2;
[DP22P2] 2
=
2 [DP22DT] 2 2
=

p1

p1
p2

p

2
T

T
2




m
2
2
m2
m2
m2
m2
m
m
H2
H2 H2
H2 H2


H2 H2





p1

p1
p2

p
T

T

2
'
Q
Q
'
Q'
Q'
Q'
Q' Q'
[DTDE]
=
p
p
p
p

;
[
]

DTDP1 = 1

1; [DTDP2]

= 2

2 [DTDT]

= T T
Q
Q
Q
Q
Q
Q
Q

p1

p


1
p2

p

2
T
T
Handbook of Référence
R7.01 booklet:
HI-75/01/001/A

Code_Aster ®
Version
5.0
Titrate:
Modelings THHM. General information and algorithms
Date:
22/06/01
Author (S):
C. CHAVANT
Key:
R7.01.10-A
Page:
28/36
In addition, by deriving the expression from the residue compared to the constraints, one defines:
1 1 1 1
1
1
1
1
1
1
1
1
1
1
1 1


1
1
1
2
2
2
1
1
1
2
2
2



p m
hm
m
hm
m
hm
m
hm
Q'
1
M1 1 1 M1 1 2 m2 2 2 m2 2
Q

2 2 2 2 2 2 2
2
2
2
2
2
2
2
2 2


1
1
1
2
2
2
1
1
1
2
2
2



p m
hm
m
hm
m
hm
m
hm
Q'
1
M1 1 1 M1 1 2 m2 2 2 m2 2
Q

3 3 3 3 3 3 3
3
3
3
3
3
3
3
3 3


1
1
1
2
2
2
1
1
1
2
2
2



p m
hm
m
hm
m
hm
m
hm
Q'
1
M1 1 1 M1 1 2 m2 2 2 m2 2
Q

4 4 4 4 4 4 4
4
4
4
4
4
4
4
4 4


1
1
1
2
2
2
1
1
1
2
2
2



p m
hm
m
hm
m
hm
m
hm
Q'
1
M1 1 1 M1 1 2 m2 2 2 m2 2
Q


D D = 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5

1
1
1
2
2
2
1
1
1
2
2
2



p m
hm
m
hm
m
hm
m
hm
Q'
1
M1 1 1 M1 1 2 m2 2 2 m2 2
Q

6.6 6.6 6.6 6
6
6
6
6
6
6
6
6


6
1
1
1
2
2
2
1
1
1
2
2
2



p m
hm
m
hm
m
hm
m
hm
Q'
1
M1 1 1 M1 1 2 m2 2 2 m2 2
Q

7 7 7 7 7 7 7 7 7 7 7 7 7 7 7


7
1
1
1
2
2
2
1
1
1
2
2
2



p m
hm
m
hm
m
hm
m
hm
Q'
1
M1 1 1 M1 1 2 m2 2 2 m2 2
Q



8 8 8
8
8
8
8
8
8
8
8
8
8
8
8

8
1
1
1
2
2
2
1
1
1
2
2
2




p m
hm
m
hm
m
hm
m
hm
Q'
1
M1 1 1 M1 1 2 m2 2 2 m2 2
Q


All these quantities not being inevitably calculated, one will note, for I from 1 to 8:
[
D I
D] I I
=
,

[
D iDP] I I I
21 =
,
,
1
1
1


m
p
2 m2 m
H2
[
D iDP] I I I
11 =
,
,
2
1
2
[
D iDP22] I I I
=
,
,
2
2
2
m1 M1 m
H 1
m2 M2 m
H2
[
D iDP12] I I I
=
,
,
2
1
2
[
D iDT] I I
=
,

m1 M1 m
H 1
Q' Q
It is then clear that:
D OF
= D D
D
. OF

And the contribution of the point of gauss to the tangent matrix
el
DF is obtained by:
G
T
el
el
el
DF = Q
D
.
OF
Q
.
G
G
G
Handbook of Référence
R7.01 booklet:
HI-75/01/001/A

Code_Aster ®
Version
5.0
Titrate:
Modelings THHM. General information and algorithms
Date:
22/06/01
Author (S):
C. CHAVANT
Key:
R7.01.10-A
Page:
29/36
7.4 Algorithm
total
The algorithm becomes then:
Initializations:
+
Calculation of meca
L
(option CHAR_MECA)
Calculation of
-
F
D
(option RIGI_MECA_TANG)
+
-
Calculation of U
by:
-
F
. U0 = meca
D
L
- meca
L
0
Iterations of balance of Newton N
Loop elements el
Loop points of gauss G
calculation
el
HQ
-
-
+
+
calculation el
E
Q U
.
and el
E
Q U
.
G
= el el
G
= el el
G
G
+
el
+
+
el
el
G
-
-
-
+
calculation of:
N
,
(according to option) starting from el
el
el
el
E, E
G
G
N
+
el
G
G
G
G
E
N
G N
+
el
+
+
+
T
el
calculation of
el
el
el
G
from
; R
Q
.
G
=
N
G
G
N
G
N
N
+
+
+
el
el
el


G
+
+
T
G
el
el
G
calculation of
N starting from el
;
N
N
el
DF
= Q.
.
Q
.
(
+
according to option)
el

G N
G
G
+
+
G
N
el
el

E

G N
G
G
N
N
Calculation of U
by:
N 1
+
+
+
+
DF. U
R
L
N
N 1 = -
N +
meca
+
Updating:
U
U
U
N
= N +
1
+
N 1
+
IF test convergence OK
end Newton: no next time
If not
N = n+1
Handbook of Référence
R7.01 booklet:
HI-75/01/001/A

Code_Aster ®
Version
5.0
Titrate:
Modelings THHM. General information and algorithms
Date:
22/06/01
Author (S):
C. CHAVANT
Key:
R7.01.10-A
Page:
30/36
8 the option
FORC_NODA
To the level of the continuous equations, option FORC_NODA corresponds to the calculation of the operator
T
el
R QT
=
el
el
.Au discrete level, option FORC_NODA amounts calculating the vector R = Q
. G.
G
G
As we already noted that depends not only on, but also of
!, one should not
to be astonished to see appearing the step of time T and the constraints at the same time at time + and time -.
The algorithm of Newton-Raphson of command STAT_NON_LINE uses option FORC_NODA for
the calculation of the prediction at the beginning of each step of time. It is not thus pain-killer to calculate
correctly all terms for this option, including those which depend on the step of time. Us
let us illustrate this question for a simple example corresponding to the only equation of hydraulics.
That is to say a simplified version of the hydraulic equation:
DM
-
p D
* +
p D
* =

F p *
M
dt

ext.


After discretization in time:
-
mp D
*
+

T (+ + (1 -) - p D
*
F p *



M
M


) =

ext.
Revealing
+
-
+
-
M = M - M

and

p = p - p, and writing a law of
behavior: m
= NR p
, one finds:
-
* +
* =
* -


*


M
M


-
NR p D
T
p D
F p
T
p D
ext.



By definition the phase of prediction STAT_NON_LINE is written:
0
1
K U = F - T
Q
ext.

0
0
It is then clear that one must take QT p D
* = - T
M-p D
*

0




Handbook of Référence
R7.01 booklet:
HI-75/01/001/A

Code_Aster ®
Version
5.0
Titrate:
Modelings THHM. General information and algorithms
Date:
22/06/01
Author (S):
C. CHAVANT
Key:
R7.01.10-A
Page:
31/36
9 Discretization
space
Finite elements THM of Aster are mixed elements, in the sense that they have at the same time
unknown factors of displacements, pressures and temperatures. A choice of discretization where
displacements, the pressures and the temperatures are interpolated with the same command of approximation
conduit with oscillations, especially for choices of step of too small times compared to
discretization spaces some One will consult on this subject for example [bib10]. This problem is also in
relationship with the manner of calculating the matrix known as of mass, and one will be able to consult on this subject [bib14].
We give in addition in appendix, to illustrate our matter, the solution for the first step of
time of a problem of mono consolidation dimensional with an interpolation P1P1. It is seen that
for a small step of time, it is very oscillating.
For this reason, quadratic elements THM are elements in P2P1, i.e.
the interpolation of displacements is quadratic and that of the temperatures and pressures is linear.
We nevertheless kept all the unknown factors on all the nodes, including the nodes mediums,
but we imposed in the calculation of the matrices of rigidity that the pressure of a node medium of
segment is equal to the half summons nodes node of the segment to which it belongs.
In addition, in the programming, we took account of the following property:
That is to say S a node node, 1
W its function of form as pertaining to a linear element (by
S
example QUAD4), and 2
W its function of form as pertaining to a quadratic element (by
S
example QUAD8). That is to say Na the number of edges having S like end and 2
W the function of form of
my
the quadratic interpolation attached to a node medium of edge, one has the relation then:
Na
1
2
w2
W
W
S =
+ my
S
ma=1 2
This said, and including in interpolation P2P1, of the conditions of nonoscillation exist on the step of
time. [bib10] the relation gives:
x2

T
>
C
20 v
kE (1 -)
where C is the coefficient of consolidation: C
, K
v =
v

being the permeability measured in
lq (1 +) (1 -

2)
m/s.
Handbook of Référence
R7.01 booklet:
HI-75/01/001/A

Code_Aster ®
Version
5.0
Titrate:
Modelings THHM. General information and algorithms
Date:
22/06/01
Author (S):
C. CHAVANT
Key:
R7.01.10-A
Page:
32/36
10 Bibliography
[1]
O. COUSSY: “Mechanical of the porous environments”. Editions TECHNIP.
[2]
P. GERMAIN: “Course of mechanics of the continuous mediums”.
[3]
G. DUVAUT, J.L. LIONS: “Inequations in mechanics and physics”.
[4]
C
. CHAVANT, P. CHARLES, Th. DUFORESTEL, F. VOLDOIRE
: “
Thermo hydro
mechanics of the porous environments unsaturated in Code_Aster “. Note HI-74/99/011/A.
[5]
I. VAUTIER, P. MIALON, E. LORENTZ: “Quasi static nonlinear Algorithm (operator
STAT_NON_LINE) “Document Aster [R5.03.01] Index B.
[6]
C. CHAVANT: “Model of behavior THHM” Document Aster [R7.01.11] Indice A
[7]
A. GIRAUD: “Adaptation to the nonlinear model poroelastic of Lassabatère-Coussy to
modeling porous environment unsaturated “, (ENSG).
[8]
J. WABINSKI, F. VOLDOIRE: Thermohydromécanique in saturated medium. Note EDF/DER
HI-74/96/010, of September 1996.
[9]
T. LASSABATERE: “Hydraulic Couplings in porous environment unsaturated with
phase shift: application to the withdrawal of desiccation of the concrete “. Thesis ENPC.
[10]
PH. MESTAT, Mr. PRAT: “Works in interaction”. Hermes
[11]
J.F. et al. THILUS: “Poro-mechanics”, Biot Conference 1998.
[12]
CROUZEIX and MIGNOT, MASSON: “Numerical Analysis of the differential equations” 1992
[13]
J.P. LEFEBVRE: “Transitory Algorithm of thermics” Document Aster [R5.02.01].
[14]
J. Mr. PROIX: “Diagonalisation of the thermal matrix of mass” Document Aster
[R3.06.07].
Handbook of Référence
R7.01 booklet:
HI-75/01/001/A

Code_Aster ®
Version
5.0
Titrate:
Modelings THHM. General information and algorithms
Date:
22/06/01
Author (S):
C. CHAVANT
Key:
R7.01.10-A
Page:
33/36
Appendix 1 Problème mono dimensional P1P1
A unidimensional problem of consolidation is considered whose unknown factors vary only according to
only variable of space X.
A rectangular field length L, is filled with a porous material of coefficients of Lamé and µ, of
1
coefficient of biot B and module of biot NR =
and of hydraulic conductivity. Density of
M
H
U

fluid is noted. One notes the constraint, U displacement in direction X, =
deformation,
xx
X

p pressure, m the mass contribution in fluid, M the flow of fluid.
The boundary conditions are:
in X = 0: = 0; p = 0
in X = L: U = 0; p = 1 for T > 0
The initial conditions in T = 0 are = U = p = 0.
X = 0;
= 0;
p = 0
X = L;
U = 0;
p = 1
The setting in equation gives:
Balance mechanical and linear elasticity: = (+ 2µ) - LP = 0
m M
Conservation of the mass:
+
= 0
T
X
p

Law of Darcy: M = -
H
X

m
Couplings:
= Np +
B

We make a discretization in implicit time and we are interested in calculation of the first step of
time:
The system of two equations is obtained:
= (+ µ
2) (U) - LP = 0

Np +
B - htp = 0
The mixed variational formulation of this problem is:


*
*
*
2
(U) U
LP U
0
U
(
[+ µ) () - ()]=




*
*
*
*
Np
B (U) p
T p p
0
p
[
+
+ H] =

éq Year 1-1
Handbook of Référence
R7.01 booklet:
HI-75/01/001/A

Code_Aster ®
Version
5.0
Titrate:
Modelings THHM. General information and algorithms
Date:
22/06/01
Author (S):
C. CHAVANT
Key:
R7.01.10-A
Page:
34/36
Concerning the space discretization, we cut out the mediums out of N finite elements. Nodes of element I
are I and i+1. 1 is noted
U and 1
p the displacement and pressure of the first node of the element E, 2
U and 2
p it
E
E
E
E
displacement and pressure of its second node. One supposes that one uses finite elements P1P1, it is with
to say that displacements as the pressures are linearly interpolated. Discretization of the first
equation gives then:
2
1
*
*
U.E. - U.E.
U
U
p
p
E -
E
E +

(
2
1
1
2
+ 2µ)
-
E
B
= 0
X
X
E


2



From where one deduces:
(2

1
1
+ 2
U - U =

E
E)
B X p
p
E
E
E
+ µ
2
2
éq Year 1-2
The discretization of the second equation gives:
2
1
2
1
*
*
2
1
1
2
*
*
2
1
p E + p E U
U
p
p
p
p p
p
E -
E
E +

E
E
- E E -

B
+ NR
+ T
H
E = 0
X
X
X
E
2


2



E


In there bearing [éq Year 1-2], one finds:
2
1
*
*
2
1
2
p E + p E
B
p
p
T
E +
E


NR +
+
p E p p p
H
- E E - E =




2
(2 1
*
*) (2
1) 0
E
2
+ 2µ
2


X E
éq Year 1-3
So now, we make tighten the step of time towards zero with step of constant space,
T

2
B


, and [éq Year 1-3] is reduced to:
H
<< NR
2

+


X

+ µ
2
(2 1
*
*
p E + p E) (1
2
p
p
E +
E) = 0
E
A total classification of the nodes and unknown factors of pressure are introduced:
1
2
p
p
p
p
J =
;
J
J =
J +1
One can see that this whole of relations gives finally:
p1 = 0

p + 2p + + p + = 0
I
I
I 1
I 2
[
1, N -]
1
p
n+1 = 1
Handbook of Référence
R7.01 booklet:
HI-75/01/001/A

Code_Aster ®
Version
5.0
Titrate:
Modelings THHM. General information and algorithms
Date:
22/06/01
Author (S):
C. CHAVANT
Key:
R7.01.10-A
Page:
35/36
The solution of this continuation is:

J - 1
p2 j-1 =


N


2 J - 1
p2 J = -

2n
What gives the distribution of following pressure:
Pressure with the first step of time
1,5
1
0,5
P
0
- 0,5 0
5
10
15
20
25
- 1
- 1,5
X
As an indication, we give Ci below a comparison between numerical results obtained and elements
P1P1, P2P2 and P2P1.
It is seen that element P2P1 does not remove the oscillation, but attenuates it appreciably.
Handbook of Référence
R7.01 booklet:
HI-75/01/001/A

Code_Aster ®
Version
5.0
Titrate:
Modelings THHM. General information and algorithms
Date:
22/06/01
Author (S):
C. CHAVANT
Key:
R7.01.10-A
Page:
36/36
Intentionally white left page.
Handbook of Référence
R7.01 booklet:
HI-75/01/001/A

Outline document