CB Manual 2018-03 v7.7
CB Manual 2018-03 v7.7
CB Manual 2018-03 v7.7
USER’S GUIDE
Version 8.6
December 2018
TABLE OF CONTENTS
• I. CODE_BRIGHT. FOREWORD
I. 1. Introduction
I. 2. System basics
I. 3. Using this manual
• II. CODE_BRIGHT. PRE-PROCESS. PROBLEM DATA.
II. 1. Problem type
II. 2. CODE_BRIGHT interface
II. 2.1. Problem data
II. 2.2. Materials
II. 2.3. Conditions
II. 2.4. Intervals data
• III. CODE_BRIGHT. PROCESS.
III. 1. Calculate
III. 2. Data Files
III. 3. General information file ROOT_GEN.DAT
III. 4. Geometrical description file ROOT_GRI.DAT
III. 5. Summary-list of cards
• IV. CODE_BRIGHT. POSTPROCESS.
IV. 1. Facilities description
IV. 2. Read Post-processing
IV. 3. Post process files format
• V. CODE_BRIGHT. THEORETICAL ASPECTS
V. 1. Basic formulation features
V. 2. Governing equations
V. 2.1 Balance equations
V. 2.2 Constitutive equations and restrictions
V. 2.3 Boundary conditions
V. 2.4 Summary of governing equations
V. 3. Numerical Approach
V. 3.1 Introduction
V. 3.2 Treatment of different terms
V. 4. Theoretical approach summary
V. 5. Features of CODE_BRIGHT
V. 6. Parallel version of CODE_BRIGHT
V. 6.1 Matrix storage mode in CODE_BRIGHT
V. 6.2 Iterative solver for nonsymmetrical linear systems of
equations
V. 6.3 Parallel version of CODE_BRIGHT
V. 7. Appendix 1. Thermo-hydro-mechanical interactions
2
• VI. CODE_BRIGHT. CONSTITUTIVE LAWS
VI. a. Hydraulic and thermal constitutive laws. Phase properties
VI. b. Mechanical constitutive laws: Elastic and visco-plastic models
VI. c. Mechanical constitutive laws: Damage-Elastoplastic model for
argillaceous rocks
VI. d. Mechanical constitutive laws: Thermo-elastoplastic model
VI. e. Mechanical constitutive laws: Barcelona Expansive model
VI. f. Mechanical constitutive laws: CASM’s family models
VI. g. Excavation/construction process
VI. h. THM discontinuities
• CODE_BRIGHT. REFERENCES
___________________________________
3
I. CODE_BRIGHT. FOREWORD
I.1. INTRODUCTION
The program described here is a tool designed to handle coupled problems in geological media.
The computer code, originally, was developed on the basis of a new general theory for saline
media. Then the program has been generalised for modelling thermo-hydro-mechanical (THM)
processes in a coupled way in geological media. Basically, the code couples mechanical,
hydraulic and thermal problems in geological media.
The theoretical approach consists in a set of governing equations, a set of constitutive laws and
a special computational approach. The code is written in FORTRAN and it is composed by
several subroutines. The program does not use external libraries.
CODE_BRIGHT uses GiD system for preprocessing and post-processing. GiD is developed by
the International Center for Numerical Methods in Engineering (CIMNE). GiD is an interactive
graphical user interface that is used for the definition, preparation and visualisation of all the
data related to numerical simulations. This data includes the definition of the geometry,
materials, conditions, solution information and other parameters. The program can also
generate the finite element mesh and write the information for a numerical simulation program
in its adequate format for CODE_BRIGHT. It is also possible to run the numerical simulation
directly from the system and to visualize the resulting information without transfer of files.
For geometry definition, the program works quite like a CAD (Computer Aided Design)
system. The most important difference is that the geometry is developed in a hierarchical mode.
This means that an entity of higher level (e.g. a volume) is constructed over entities of lower
level (e.g. a surface); two adjacent entities (e.g. two volumes) will then share the same lower
level entity (e.g. a surface).
All materials, conditions and solution parameters can also be defined on the geometry without
the user having any knowledge of the mesh. The meshing is performed once the problem has
been fully defined. The advantages of doing this are that, using associative data structures,
modifications can be made on the geometry and all other information will be updated
automatically.
Full graphic visualisation of the geometry, mesh and conditions is available for comprehensive
checking of the model before the analysis run is started. More comprehensive graphic
visualisation features are provided to evaluate the solution results after the analysis has been
performed. This post-processing user interface is also customisable depending on the analysis
type and the results provided.
A query window appears for some confirmations or selections. This feature is also extended to
the end of a session, when the system prompts the user to save the changes, even when the
normal ending has been superseded by closing the main window from the Window Manager,
or in most cases with incorrect exits.
4
the geometry while maintaining the attributes and conditions definitions. Alterations to the
attributes or conditions can simultaneously be made without the need of reassigning to the
geometry. New meshes or small modifications on the obtained mesh can also be generated if
necessary and all the information will be automatically assigned correctly.
The system does provide the option for defining attributes and conditions directly on the mesh
once this has been generated. However, if the mesh is regenerated, it is not possible to maintain
these definitions and therefore all attributes and conditions must be redefined. In general, the
complete solution process can be described as:
1. Define geometry - points, lines, surfaces, volumes.
• Use other facilities.
• Import from CAD.
2. Define attributes and conditions.
3. Generate mesh.
4. Carry out simulation.
5. View results.
Depending upon the results in step (5) it may be necessary to return to one of the steps (1), (2)
or (3) to make alterations and rerun the simulations.
Building a geometrical domain in GiD is based on the 4 geometrical levels of entities: points,
lines, surfaces and volumes. Entities of higher level are constructed over entities of lower level;
two adjacent entities can therefore share the same level entity.
All domains are considered in 3-dimensional space but if there is no variation in the third
coordinate (into the screen) the geometry is assumed to be 2-dimensional for analysis and
results visualisation purposes. Thus, to build a geometry, the user must first define points, join
these to form lines, create closed surfaces from the lines and define closed volumes from the
surfaces. Many other facilities are available for creating the geometrical domain; these include:
copying, moving, automatic surface creation, etc.
The geometrical domain can be created in a series of layers where each one is a separate part
of the geometry. Any geometrical entity (points, lines, surfaces or volumes) can belong to a
particular layer. It is then possible to view and manipulate some layers and not others. The main
purpose of the use of layers is to offer a visualisation and selection tool, but they are not used
in the analysis.
The system has the option of importing a geometry or mesh that has been created by a CAD
program outside GiD; at present, this can be done via a DXF, IGES or NASTRAN interface.
Once the geometry and attributes have been defined, the mesh can be generated using the mesh
generation tools supplied within the system. Structured and unstructured meshes containing
triangular and quadrilateral surface meshes or tetrahedral and hexahedral volume meshes may
be generated. The automatic mesh generation facility utilizes a background mesh concept for
which the users are required to supply a minimum number of parameters.
Simulations are carried out by using the calculate menu. The final stage of graphic
visualisation is flexible in order to allow the users to critically evaluate the results quickly and
easily. The menu items are generally determined by the results supplied by the solver module:
this not only reduces the amount of information stored but also allows a certain degree of user
customisation. The post solver interface may be included fully into the system so that it runs
automatically once the simulation run has terminated.
5
I.3. USING THIS MANUAL
This User Manual has been split into several differentiated parts. The part, THEORETICAL
ASPECTS, contains the theoretical basis of CODE_BRIGHT, and the numerical solution. In
CODE_BRIGHT. PREPROCESS. PROBLEM DATA, it is described how to enter the data
of the problem, i. e. general data, constitutive laws, boundary conditions, initial conditions and
interval data. The referred as CODE_BRIGHT. PROCESS is related to the calculation
process. This part also contains the description of input files. The part, CODE_BRIGHT.
CONSTITUTIVE LAWS contains a description of hydraulic, thermal and mechanical
constitutive laws and phase properties. Finally, CODE_BRIGHT. TUTORIAL, introduces
guided examples for a fast and easy familiarization with the system.
___________________________________
6
II. CODE_BRIGHT. PRE-PROCESS. PROBLEM DATA.
Problem data include all the parameters, conditions (see section Conditions), materials
properties (see section Materials), problem data (see section Problem Data) and intervals data
(see section Interval Data) that define the project. Conditions and materials should be assigned
to geometrical entities.
In order to build the data files ROOT_GEN.DAT and ROOT_GRI.DAT, data is introduced into
several window statements associated with these concepts (interface inputs).
Once the geometry has been prepared, it is necessary to go through the different Interface steps,
i.e. PROBLEM DATA, MATERIALS, CONDITIONS, and INTERVAL DATA. See the
tutorials for a guided introduction to the interface between GiD and CODE_BRIGHT.
7
GENERAL DATA
Title of the problem Interface Default: Coupled problem in geological media
Execution Only data file generation: ROOT_gen.dat and ROOT_gri.dat are
built
Full execution: Calculation with the finite element program
CODE_BRIGHT is performed (default option)
Backup No Backup
(IMBACKUP in Save Last: Allows restart the calculation from the last time step computed.
root_gen.dat) Information is saved in file root_save.dat
Save All: Allows restart the calculation from any time step computed.
Information is saved in different files root_tnum_save.dat for each interval
data computed (tnum). To restart the calculation, rename the file
root_tnum_save.dat to root_save.dat
Axisymetry No, Around y-axis
(IAXISYM in In 2-D axisymetry the principal stresses are: σ r (radial), σ y (axial), σ θ
root_gen.dat) (circumferential)
Gravity X (2-D and 3-D) component Interface Default: 0.0
Gravity Y (3-D) component Interface Default: 0.0
Gravity Y (2-D) or Z (3-D) component Interface Default: -9.81
EQUATIONS SOLVED
Stress equilibrium (unknown Yes, No
displacement u) (IOPTDISPL
in root_gen.dat)
Updated lagrangian Yes, No. Updated lagrangian method, i.e. coordinates are modified after
method (IUPDC in each time increment is solved. If deformations are very large, some
root_gen.dat) elements may distort. If distortion is very large the volume of an element
may become negative and the execution will terminate immediately.
Mass balance of water (unknown liquid Yes, No
pressure Pl) (IOPTPL in root_gen.dat)
Constant Pl Constant liquid phase pressure for problems not including the
mass balance of water equation
Mass balance of air (unknown liquid Yes, No
pressure Pg) (IOPTPG in root_gen.dat)
Constant Pg Constant gas phase pressure for problems that do not include
the equation of mass balance of air. Usually equal to 0.1 MPa.
Dissolved air into liquid phase Allowed, Not allowed
(IOPTXAL in root_gen.dat)
Energy balance (unknown temperature) Yes, No
(IOPTTEMP in root_gen.dat)
Vapour into gas phase (IOPTXWG in Allowed, Not Allowed
root_gen.dat)
Constant Temp Constant temperature for problems that do not include the
equation of energy balance.
Mass balance of conservative species Yes, No
(unknown concentration) (IOPTXWS in
root_gen.dat)
8
Combinations of solving options are described below:
Pl Pg T Variable
1 0 0 Compressible water flow, one phase, one species, air is not considered.
0 1 0 Compressible air flow, one phase, one species.
0 0 1 Heat flow (only conduction).
1 1 0 Two phase flow (liquid + gas), air dissolved permitted, vapour not permitted.
1 0 1 Water two phase non-isothermal flow, vapour allowed, gas phase at constant
pressure.
0 1 1 Compressible non-isothermal gas flow, one phase, one species.
1 1 1 Non-isothermal two phase (liquid + gas) flow, vapour and air dissolved are allowed.
SOLUTION STRATEGY
Position of intermediate time tk+ε for matrix evaluation, i.e.
Epsilon (intermediate time for the point where the non-linear functions are computed.
nonlinear functions) (usual values: 0.5, 1). See details on Numerical Method.
Default: 1.0
Theta (intermediate time for Position of intermediate time tk+θ for vector evaluation, i.e.
implicit solution) the point where the equation is accomplished. Default: 1.0
0-4: Time step control based on N-R iterations:
0: no time step prediction is performed.
1: predicts time stepping according to a limit of 4 iterations.
2: predicts time stepping according to a limit of 3 iterations.
3: predicts time stepping according to a limit of 2 iterations.
4: predicts time stepping according to a limit of 1 iteration.
Time step control 6-9: Time step control based on error estimation:
(ITIME in root_gen.dat) 6: controls time stepping by means of a prediction based on
Default: 1 the relative error deviation in the variables (relative error
lower than 0.01).
7: same as 6 but with a tolerance equal to 0.001.
8: same as 6 but with a tolerance equal to 0.0001.
9: same as 6 but with a tolerance equal to 0.00001.
Note: a time step control = 1 will always be considered for
negative time.
Max. number of iterations per Maximum number of Newton Raphson iterations per time
time step (ITERMAX in step. If the prescribed value is reached, time step is reduced.
root_gen.dat) Default: 10
Solver type Direct: LU + BACK
(ISOLVE in root_gen.dat) Iterative: Sparse + CGS
Max number of solver iterations Default: 5000
Solver type = Max abs solver error variable Default: 1.e-9
Iterative: Sparse +
CGS Max abs solver error residual Default: 0
Max rel solver error residual Default: 0
9
Elemental relative Elemental suction (consistent approach)
permeability Average nodal degrees of saturation (default)
computed from: Average nodal relative permeabilities
(IOPTPC in Average nodal relative permeabilities (applies also for derivatives)
root_gen.dat) Maximal nodal relative permeability
Max Abs Displacement (m) Maximum (absolute) displacement
(DELMXU in root_gen.dat) error tolerance (m). When correction
of displacements (displacement
difference between two iterations) is
lower than this value, convergence has
been achieved. Default: 1e-6
Stress equilibrium
Max Nod Bal Forces (MN) Maximum nodal force balance error
(unknown tolerance (MN). If the residual of
displacement u) = yes (DELFMX in root_gen.dat) forces in all nodes are lower than this
value, convergence has been achieved.
Default: 1e-10
Displacement Iter Corr (m) Maximum displacement correction per
(DUMX in root_gen.dat) iteration (m) (time increment is
reduced if necessary). Default: 1e-1
Max Abs Pl (MPa) (DELMXPL Maximum (absolute) liquid pressure
in root_gen.dat) error tolerance (MPa). Default: 1e-3
Mass balance of
Max Nod Bal Forces (MN) Maximum nodal water mass balance
water (unknown
error tolerance (kg/s). Default: 1e-10
liquid pressure Pl) = (DELQWMX in root_gen.dat)
yes Pl Iter Corr (MPa) (DPLMX in Maximum liquid pressure correction
root_gen.dat) per iteration (MPa) (time increment is
reduced if necessary). Default: 1e-1
Max Abs Pg (MPa) Maximum (absolute) gas pressure
(DELMXPG in root_gen.dat) error tolerance (MPa). Default: 1e-3
Mass balance of air Max Nod Air Mass (kg/s) Maximum nodal air mass balance
(unknown liquid (DELQAMX in root_gen.dat) error tolerance (kg/s). Default: 1e-10
pressure Pg) = yes Pg Iter Corr (MPa) (DPGMX in Maximum gas pressure correction per
root_gen.dat) iteration (MPa) (time increment is
reduced if necessary). Default: 1e-1
Max Abs Temp (C) (DELMXT Maximum (absolute) temperature
in root_gen.dat) error tolerance (C). Default: 1e-3
Energy balance Max Nod Energy (J/s) Maximum nodal energy balance error
(unknown (DELQMX in root_gen.dat) tolerance (J/s). Default: 1e-10
temperature) = yes Temp Iter Corr (C) (DTMX in Maximum temperature correction per
root_gen.dat) iteration (C) (time increment is
reduced if necessary). Default: 1e-1
Max Abs Solute (DELMXI in Maximum (absolute) concentration
root_gen.dat) error tolerance. Default: 1e-3
Mass balance of Max Nod Solute mass balance Maximum nodal solute mass balance
conservative species (DELIMX in root_gen.dat) error tolerance. Default: 1e-10
(unknown Solute Iter Corr Maximum solute concentration
concentration) = yes correction per iteration (time
DIMX in root_gen.dat)
increment is reduced if necessary).
Default: 1e-1
10
Comments regarding the use of tolerances
In order to illustrate the use of tolerances the thermal problem is considered with the following
tolerances:
Max Abs Temp (C) T1
Max Nod Energy (J/s) T2
Temp Iter Corr (C) T3
Convergence can be achieved in two ways: the one when δT < T1 for all nodes (condition A)
and the second when (q h < T2) also for all nodes (q h represents here the energy balance or
residual at a node) (condition B).
It is to be mentioned that convergence in terms of δT and convergence in terms of q h should be
reached simultaneously because the Newton - Raphson method is used. For this reason, the
program stops the iteration process when one of the two conditions (A or B) is achieved.
When more than one degrees of freedom are solved per node and one of the recomended options
is used (convergence by variable OR residual), convergence in terms of variable or residual
should be achieved by all the variables simultaneously. In other words, it is not possible that
the mechanical problem converges by residual and the thermal problem converges by the
variable.
Finally, if (δT > T3), time increment will be reduced. This parameter controls the accuracy of
the solution in terms of how large time increments can be. A low value of T3 will force to use
small time increments when large variations of temperature take place.
OUTPUT
Write Iteration information is written in file ROOT_GEN.OUT according to:
numerical NONE: no information about convergence is written. This option should
process be used if the user is very confident with the time discretization and not
information interested in details at every time step or problems with time increment
reductions. Usually this happens when previous runs have shown that
(IOWIT in convergence and time discretization work very well.
root_gen.dat) PARTIAL: partial information is written. Time intervals and time-values,
number of iterations, CPU-time values, etc. are written. Convergence
information (e.g. residuals) is only written if time increment reductions
take place.
ALL: all iteration information is written. Convergence information is
written for all iterations and all time increments. This option may result
in a very large file ROOT_GEN.OUT
11
Writing Writing results frequency in output files according to the number of time steps
frequency (positive integer value) or according to a given time increment (negative integer
(INTER in value).
root_gen.dat) If it is positive, e.g. is set to 20, results for the complete mesh will be written
only every 20 calculated time increments.
If it is negative, then we can obtain the output values in a specified time: e.g.
setting a value of -10 will produce output for 0, 10, 20, 30, … units of time.
Note that you may need to set a suitable maximum time step in the interval data
in order for this implementation to work well (the maximum time step should
be around one order of magnitude lower than the writing time frequency). See
Figures II.2.1a, b and c.
Writing every
20 time steps
Writing every
100 time steps
12
Writing every 10 days
(input -10 as writting frequency
when using days as time units)
Figure II.2.1c. Writing every 10 days (Writing frequency = -10 and days selected in the
interval data).
OUTPUT
(continuation)
Write piezometric head Yes, No
Write boundary flow No (Defaul option)
rates in additional file Use writing frequency
Write all
Write boundary No (Defaul option)
reactions in additional Use writing frequency
file Write all
Output points
Nodes
(IOWCONTOURS in
root_gen.dat) Gauss points: (Default option)
Write all information Yes (default option), No
(IWRALL in If No is selected, the following option appears:
root_gen.dat) Separated output files (IPOLYFILES in root_gen.dat) : Yes, No
and user go to Select output window.
13
SELECT OUTPUT
14
II.2.2. MATERIALS
All materials must be defined from a generic material. The following steps show how to
assign materials and do modifications:
- Creating new materials: In order to create new materials, one should write a material
name and complet the necessary constitutive laws and do an Accept Data to validate the
data entered. It is necessary to create a material before assigning it on the geometry.
- Assignment must respect hierarchical structure of entities (i.e. cannot assign a material
on a line belonging to a surface that have just been identified with another material). This
type of error may create conflicts.
- Posterior modifications on the parameters of assigned materials do not require a re-
meshing process.
- Material names: When introducing a name for a material, it is strongly recommended to
avoid spaces or underscores (e.g. use mat1 instead mat_1 or mat 1). The use of spaces or
underscores ( _ ) might create conflicts when the material is read.
ITYCL P1 P2 P3 P4 P5 P6 P7 P8 P9 P10
A number indicates the intervals where the law will be defined. This number fixes the number
of lines for VALUES to be entered. Every Interval line assumes parameters of INTERVAL
DATA according to the same order.
15
The following constitutive laws are available:
Assign material
With this instruction, the material is assigned to the selected entities. If assigning from a
window, every time the assigned material changes, the button Assign must be pressed again.
The user must select the entity on which to assign the materials, i.e.: line, surface or
volume when working in geometry mode or directly over the elements when working in
mesh mode. It is recommended to assign the materials on the geometry entities rather than on
the elements.
If assigning from the command line, option UnAssignMat erases all the assignments of this
particular material.
When a mesh has been already generated, and changes in the assigned materials are required,
then it is necessary to re-mesh again or assign the materials directly on the mesh.
Draw material
Draws a color indicating the selected material for all the entities that have the required material
assigned. It is possible to draw just one or draw all materials. To select some of them the users
should use a:b and all material numbers that lie between a and b will be drawn.
When drawing materials in 3 dimensions, it may be necessary to change the viewing mode to
polygons or render (see section Render) to diferenciate the front and back of the objects.
Unassign material
Command Unassign unassigns all the materials from all the entities. For only one material,
use UnAssignMat (see section Assign material).
16
New material
When the command NewMaterial is used, a new material is created taking an existing
one as a base material. Base material means that the new one will have the same fields as the
base one. Then, all the new values for the fields can be entered in the command line. It is
possible to redefine an existing material.
To create a new material or redefine an existing one in the materials window, write a new name
or the same one and change some of the properties. Then push the command Accept.
These types of elements are assigned by the interface between GiD and CODE_BRIGHT.
Note that linear triangular elements or linear tetrahedrons, which have been proven to be very
adequate for flow problems, should be avoided for mechanical problems. This is because if the
medium is nearly-incompressible (creep of rocks takes place with very small volumetric
deformation), locking takes place (not all displacements are permitted due to element
restrictions).
17
II.2.3. CONDITIONS
Conditions are all the properties of a problem, excluding materials, that can be assigned to an
entity. In this concept several types of conditions have been included: Force/Disp conditions,
Flux conditions, Initial unknowns, Porosity (and other variables), Initial stress, Joint element
width, time evolution location, etc. The condition window permits to choose entities to assign
on (Point, Line, Surface or Volume in geometry display mode and Node or Element in mesh
display mode) and select different types of conditions. It must be taken into account that
conditions assigned in mesh display mode will be unassigned in every new meshing process.
The following points should be taken into account for condition construction:
• Force/Disp conditions add up all conditions assigned at every node, except for variables
Index (takes last value encountered) and Multiplier (takes the biggest).
• Flux conditions, Initial unknowns, Porosity (and other variables), Initial stress and Joint
element width are assigned with entities priority in the following order: Points, Lines,
Surfaces and Volumes (i.e. the node takes a Flux_Point_B.C. refusing a Line_Flux_B.C.
assigned previously).
If a mesh has already been generated, for any change in the condition assignments, it is
necessary to re-mesh again to transfer these new conditions to the mesh.
Conditions description
II.2.3.1 Force/displacement conditions
The mechanical boundary conditions only exist if the mechanical problem is solved (Solve
displacement). For each time period only the types that undergo changes should be read.
18
= 1 means that displacement rate will be
X direction prescribed prescribed in the X direction. The value
is given above.
= 1 means that displacement rate will be
Y direction prescribed prescribed in the Y direction. The value
is given above.
= 1 means that displacement rate will be
Z direction prescribed prescribed in the Z direction. The value
is given above.
The units of this parameter depend on
whether force or stress is applied:
MN s
γ (multiplier) When applying a force:
m
MPa s
When applying a stress:
m
∆f x o obtained as ramp
loading during the current
interval.
∆f y o obtained as ramp
loading during the current
interval.
∆f z o obtained as ramp
loading during the current
interval.
Points (2-D or 3-D) Lines (usually 2-D) Surfaces (usually 3-D) Volumes (3-D)
Forces Forces Forces Forces
Boundary stresses Boundary stresses
19
II.2.3.2 Flux Boundary Condition
Mass or heat transport problems. These conditions only exist if any balance (water, air, energy
flow) problem is solved. For each time period only the types that undergo changes need to be
read.
The boundary condition is incorporated by adding a flux or flow rate. The mass flux or flow
rate of species i = w as a component of phase α = g (i.e. the inflow or outflow of vapour) is
calculated as:
( ) ( ) ( )
j gw = ω wg j g0 + ω wg γ g Pg0 − Pg + β g ρ g ω wg − ρ g ω wg ( ) ( )
0 0 0
where the superscript 0 stands for the prescribed values, ω is mass fraction, ρ is density, P g is
gas pressure, j g 0 is a prescribed gas flow and γ g and β g are two parameters of the boundary
condition. Particular cases of this boundary condition are obtained for instance in the following
way:
Description (ωg w)0 jg0 γg Pg0 (ρ g ) βg
A prescribed mass flow rate of gas with 0.02 0.02 1e-5
kg/kg of vapor and 0.98 kg/kg of air is injected kg/s
If Pg < Pg0 = 0.1 => a variable mass flow rate 0.02 10 0.1
of gas with 0.02 kg/kg of vapor and 0.98 kg/kg
of air is injected.
If Pg < Pg0 = 0.1 => a variable mass flow rate
of gas with variable composition outflows.
Humidity in the boundary is prescribed to 0.01 1.12 10
0.0112 kg/m3. This is equivalent to a relative
humidity of 0.0112/0.0255 = 0.44 = 44%
jlw = (ω )w 0
l j + ( ω ) γ ( P − P ) + β ((ρ ω ) − (ρ ω ))
0
l
w 0
l l l
0
l l l
w 0
l l
w
l
(ω ) w 0
= 1 − ( ωla )
0
l
20
Positive values of mass flow rate indicate injection into the medium.
For energy, the boundary condition has the general form:
( )
je = je0 + γ e T 0 − T + E gw ( j gw ) + ...
In other words, a von Newman type term plus a Cauchy type term and a series of terms that
represent the energy transfer caused by mass inflow and outflow through the boundary.
The set of parameters that are required for these equations are (note that the symbols used by
Gerard et al., 2009, are not the same):
21
For a positive value of δ a parabolic curve is used; for a negative value an exponentially
decaying curve is used. δ is the distance from the reference pressure to the point of change.
Index → +1.0 means that all flow rates are nodal values. For instance,
(auxiliary a pumping well boundary condition.
index) → -1.0 means that all flow rates are per unit volume (3-D), area
(2-D) or length (1-D) of medium (internal source or
sink). For instance, a recharge due to rain in a 2-D case.
→ +2.0 means that all flow rates are per unit area (3-D) or
length (2-D) (lateral fluxes). For instance, lateral fluxes
from neighbour aquifers.
Prescribed gas, liquid and heat flows must be given in terms of flow units depending on the
way these flows are considered, i.e., depending on the kind of element they pass through and
on the problem dimension. The required units for each case are graphically specified below:
INDEX PROBLEM
ILLUSTRATION FLOW UNITS
PARAMETER DIMENSION
Index = 1.0 3-D kg
Mass:
s
J
Heat:
s
2-D kg
Mass:
s
J
Heat:
s
1-D kg
Mass:
s
J
Heat:
s
Index = -1.0 3-D kg
Mass:
m3 s
J
Heat:
m3 s
2-D kg
Mass:
m2 s
J
Heat:
m2 s
1-D kg
Mass:
ms
J
Heat:
ms
22
Index = 2.0 3-D kg
Mass:
m2 s
J
Heat:
m2 s
2-D kg
Mass:
ms
J
Heat:
ms
Required units
Problem Gas Liquid Heat
index Parameters Parameters Parameter
dimension flow flow flow
γg and γl βg and βl γe
rate jg rate jl rate je
kg kg kg m3 J J
1.0 ---
s s s MPa s s sC
kg kg kg m3 m2 J J
1D =
ms ms m s MPa ms s ms msC
kg kg kg m 3
m J J
-1.0 2D =
2
m s 2
m s 2
m s MPa 2
m s s m2 s m2 s C
kg kg kg m 3
1 J J
3D 3 = 3
3
m s 3
m s m s MPa m3 s s 3
m s m sC
kg kg kg m 3
m2 J J
2D =
ms ms m s MPa ms s ms msC
2.0
kg kg kg m 3
m J J
3D 2 2 =
m s m s 2
m s MPa 2
m s s m2 s m2 s C
The fact that units are different for 3D, 2D and 1D is due to the reduction of one dimension in
2D and two dimensions in 1D. However, if a 2D model is considered to have 1 m associated
thickness, then units would be identical as in 3D. Similarly, if a 1D model is considered to have
1 m2 associated surface then units would be identical as in 3D.
The above boundary conditions are rather general. They incorporate terms of von Newman type
and Cauchy type. The equation includes three terms. The first one is the mass inflow or outflow
that takes place when a flow rate is prescribed at a node. The second term is the mass inflow or
outflow that takes place when a phase pressure is prescribed at a node. The coefficient γ is a
leakage coefficient. This variable allows prescribing a pressure with more or less strength. If γ
is very large, pressure will tend to reach the prescribed value (see Figures II.2.2 and II.2.3).
However, an extremely large value can produce matrix ill conditioning and a lower one can
produce inaccuracy in prescribing the pressure. However, it is not difficult to guess adequate
23
values for a given problem simply by trial. The third term is the mass inflow or outflow that
takes place when species mass fraction is prescribed at a node.
A surface where seepage (only outflow for liquid phase is permitted) is a case that may be of
interest. To indicate that only outflow is permitted γl is entered with negative sign. This negative
sign only indicates that nodes with this kind of boundary condition allow seepage (i.e. only
outflow).
If there is inflow of gas or liquid phase, it is very important to give values of the following
variables: (ωgw)o, (ωl a)o, (ρl)o, (ρg)o and To. Otherwise, they are assumed to be zero, which is
not correct because they will be too far from equilibrium. If outflow takes place, this is not
relevant because the values of the medium are used instead of the prescribed ones.
qi
γi>0
γi
1.0
inflow Poi
Pi
outflow
Figure II.2.2
qi
γi<0
inflow δ Poi
Pi
outflow γi
1.0
Figure II.2.3
24
Boundary conditions variable with time.
a) One boundary condition variable with time.
It is possible to assign a flux condition varying with time using an auxiliary file. In this case,
the user must assign a value of -999 to the particular variable for which a given variation wants
to be assigned and include an ASCII file called “root_bcf.dat” in the GiD project folder with
the variation with time of the given flux variable. The structure of the “root_bcf.dat” file is
illustrated in Table II.2.2. Note that with root we mean the root name of the project: for example,
if our project is DAM.gid then the file should be named DAM_bcf.dat.
Table II.2.2. Illustration of the format of the root_bcf.dat file
Number of data (N)
Number of Flux Id. variable (1) Id. variable (2) … Id. variable
variables (NF) (NF)
Time (1) Value Value … Value
Time (2) Value Value … Value
… … … … …
Time (N) Value Value … Value
It should be noted that the first line of the root_bcf.dat file must contain the number of data (N)
that has to be read. The first column of the second line refers to the number of flux variables
(NF) for which a given variation with time want to be assigned. The other columns of the second
line contain a special flag (or indicator) of the flux variables to be changed. This indicator is
shown in Table II.2.3. The following lines (from third to N) contain the time in the same unit
considered at the interval data (in the first column) and the values of the flux variables assigned
to the specific time (in the other columns).
Table II.2.3. Identification number (Id.) of the flux variables
Id. Flux variable
1 ωgw prescribed mass fraction (kg/kg)
2 jg prescribed gas flow rate (units in Table II.2.1)
3 Pg prescribed gas pressure (MPa)
4 γg parameter for gas pressure term (units in Table II.2.1)
5 βg parameter for humidity term (units in Table II.2.1)
6 ρg prescribed gas density (kg/m3)
7 ωl w prescribed mass fraction of solute (kg/kg)
8 ωl a prescribed mass fraction of air (kg/kg)
9 jl prescribed liquid flow rate (units in Table II.2.1)
10 Pl prescribed liquid pressure (MPa)
11 γl parameter needed to be ≠ 0 when Pl is prescribed (units in Table II.2.1)
12 βl parameter needed only when mass transport problem is considered (Table II.2.1)
13 ρl prescribed liquid density (kg/m3)
14 je prescribed heat flow rate (units in Table II.2.1)
15 T prescribed temperature (C)
16 γe parameter needed to be ≠ 0 when T is prescribed (units in Table II.2.1)
25
b) Multiple boundary conditions variable with time.
For a value in a flux boundary condition that varies with time according to a function of time
(list of values) the following actions have to be taken:
• Set de value of the variable equal to -99901.
• If you have more values for the same boundary condition or in another boundary condition,
use the values -99902, -99903 up to -99920 if necessary.
• Prepare a file (root_bcf.dat) with the following structure:
- First row: m, number of data values to be read (lines).
- Second row: –n 1 2 3 … n where n ≤ 20 number of variables that change with time.
- m rows with time_1 value_1, time_2, value_2, … time_n, value_n.
1
The first implementation of this module is due to Maarten Saaltink.
26
The calls to atmosferic_boundary_condition subroutine appear in the following
subroutines:
• atm_boundary_conditions (bcond_flow.f) itself called by newton_raphson
(nr.f)
• write_boundary_flows (write.f) itself called by main_calculate
(code_bright_main.f)
Running atmosferic_boundary_condition needs index=5 2 as boundary condition
type, this latter being passed in FLUX(20). FLUX vector is read from file root_gen.dat by
read_boundary_conditions (read_general.f). The index of this file is iin1 (the
concerned card is numbered 20).
It should be noted that flow rates computed by the atmosferic_boundary_condition
subroutine are then treated as if index = 2.0 (used in classic flux boundary conditions)
were set i.e. as if flow rates are per unit area.
BEGIN
call get_atm_data
energy flux
(*)
write to output per radiation call sun
per convection
2
This variable is locally called ICON.
27
Input data
Problem data definition file (root_gen.dat)
When atmospheric boundary conditions are considered, the parameters presented in Table II.2.4
have to be entered in CardGroup 20 in file root_gen.dat. These conditions are activated via
FL20, which should be set to 5. In Table II.2.4, latitude, time when autumn begins, time at
noon, dry and wet albedos are used for calculating radiation when radiation type is lower than
3.
Roughness length, screen height and stability factor are used for evaporation estimation and for
estimation of the advective energy flux.
Table II.2.5 presents ranges of roughness lengths for different types of surfaces from which
evaporation has to be calculated.
Table II.2.5: Roughness lengths for different types of surfaces after (Chow et al. 1988).
Type of surface Height main roughness (m)
Ice, mud flats 1.10-5
Water 1.10-4 – 6.10-4
Grass (up to 10 cm high) 1.10-3 – 2.10-2
Grass (10 to 50 cm high) 2.10-2 – 5.10-2
Vegetation (1–2 m) 0.2
Tress (10 – 15 m) 0.4 – 0.7
28
Atmospheric data input file (root_atm.dat)
General parameters are to be entered in the problem data file (see CardGroup 20 description
here above) but time varying atmospheric data which are required to compute mass and heat
fluxes should be entered within an ASCII file called root_atm.dat 3. Data that can be read is
summarized in Table II.2.6.
For each variable, the pair of columns containing available data is organised is the way
schematically presented in Table II.2.7 4. More details about this file and time varying
atmospheric data are given is section 0, dedicated to get_atm_data subroutine.
3
The file root_atm.dat is read with a free format. A dedicated tool developed by J.M. Pereira
(atmdata.exe) can be used to check its general format.
4
In this table, light grey and bold grey cells respectively identify measured data and unused
data (the former being constituted of time (ti) and corresponding values (xi) pairs for each
quantity).
29
Table II.2.6: Time varying atmospheric data to be provided in file root_atm.dat.
Data Unit
Atmospheric temperature, Ta °C
Atmospheric gas pressure, Pga MPa
Relative humidity, Hr -
Radiation 5, Rn J.m-2.s-1
Cloud index 6, In -
Rainfall, P kg.m-2.s-1
Wind velocity, va m.s-1
Table II.2.7: Illustration of the format of root_atm.dat file (excluding first line).
Ta Pga Hr Rn In P va
Flag 0 0 0 0 0 0 0 0 0 0 0 0 0 0
Annual mean 0 x_am 0 x_am 0 x_am 0 x_am 0 x_am 0 x_am 0 x_am
Annual ampl 0 x_aa 0 x_aa 0 x_aa 0 x_aa 0 x_aa 0 x_aa 0 x_aa
Annual gap (s) 0 x_ag 0 x_ag 0 x_ag 0 x_ag 0 x_ag 0 x_ag 0 x_ag
Daily ampl 0 x_da 0 x_da 0 x_da 0 x_da 0 x_da 0 x_da 0 x_da
Daily gap (s) 0 x_dg 0 x_dg 0 x_dg 0 x_dg 0 x_dg 0 x_dg 0 x_dg
Unused 0 0 0 0 0 0 0 0 0 0 0 0 0 0
Measures… ti xi ti xi ti xi ti xi ti xi ti xi ti xi
Measures… ti xi ti xi ti xi ti xi ti xi ti xi ti xi
Measures… ti xi ti xi ti xi ti xi ti xi ti xi ti xi
Measures… ti xi ti xi ti xi ti xi ti xi ti xi ti xi
Measures… ti xi ti xi ti xi ti xi ti xi ti xi ti xi
Measures… ti xi ti xi ti xi ti xi ti xi ti xi ti xi
Measures… … … … … … … … … … … … … … …
It should be noted that the first line of the file root_atm.dat must contain the number of lines
(excluding this one) and number of columns that has to be read. The second line refers to
interpolation or simulation option: it corresponds to a special flag allowing the user to simulate
the atmospheric data on the base of annual and daily characteristics that are furnished (see eq.
(20)). This simulation will be processed if, for a given quantity, its flag is set to ‘0’. The values
necessary to proceed to this simulation are provided in the 5 following lines and correspond to
annual mean, amplitude and gap and daily amplitude and gap. On the contrary, if this flag is set
to ‘1’, Code_Bright will use the measured data provided in the rest of the lines of the data file
and process to linear interpolations in order to obtain the value of a quantity for a given
calculation time.
5
Radiation data will be used only if “Radiation type” is set to ‘2’ or ‘3’ in the boundary
conditions parameters.
6
Cloud index allows to account for a cloudy sky in the radiation computation (In = 1 for a
clear sky and In = 0 for a completely cloudy sky).
30
Subroutines description
Subroutine atmosferic_boundary_condition
General description
This subroutine is the core of the atmospheric boundary condition module. It computes
atmospheric boundary conditions, including evaporation, rain, radiation, advective and
convective energy fluxes. Those are expressed in terms of fluxes of water, air and energy as
functions of the state variables (liquid pressure, gas pressure and temperature of the soil).
Moreover, the subroutine calculates the derivatives of these three fluxes with respect to the state
variables. Positive values always mean entering the system. Negative values always mean
leaving the system e.g. evaporation is negative.
It first calls get_atm_data to read atmospheric data which is stored in a matrix named
atmosferic 7 (see subroutine get_atm_data for more information on the format of this
file). Some general parameters (like for instance dry and wet albedos) are read from FLUX, an
argument passed to the main subroutine (which corresponds to CardGroup 20 of the problem
data file root_gen.dat).
It then computes water flux (through gas and liquid phases, due to evaporation and rain), air
flux and energy flux (radiation, advective and convective energy fluxes). Optionally, it can
write these fluxes to files depending upon the presence of surveyed nodes or not. The general
equations for calculating these fluxes are now presented.
Fluxes of mass
- Flux of gas: The flux of the gas phase qg is given by the following equation, in which Pga is
the atmospheric pressure and γg is a leakage coefficient:
qg = γ g (Pg − Pga ) (1)
- Flux of air: For the flux of air ja only the advective part is considered:
(
ja = ω ga qg = 1 − ω gw qg ) (2)
- Flux of water: Evaporation E is given by an aerodynamic diffusion relation:
k 2 vaϕ
=E 2 ( ρva − ρv ) (3)
za
ln
z0
where ρva and ρv respectively are the absolute humidity (mass of vapour per volume of gas,
which can be calculated from relative humidity Hr and temperature) of the atmosphere and at
the node of the boundary condition, k is the von Karman’s constant (often taken as 0.4), φ is a
stability factor, va the wind velocity, z0 is the roughness length, za is the screen height at which
va and ρva are measured. In theory, ρv must be the value at roughness length (z0). Instead, it is
calculated from the state variables at the node of the boundary condition. Hence, a constant
profile for ρv is assumed between this node and height z0.
7
Matrix atmosferic is stored in bb(n74).
31
The advective flux of vapour by the gas phase jgw is given by:
jgw = ω gw qg if Pg > Pga
jw = ρva (4)
qg if Pg ≤ Pga
g ρ ga
where ρga is the atmospheric gas density and qg is the flux of the gas phase given by equation
(1).
Surface runoff jsr (which corresponds to the flow rate of water through the liquid phase jlw) is
written as:
j sr = γ w (Pl − Pga ) if Pl > Pga
(5)
j sr = 0 if Pl ≤ Pga
where γw is another leakage coefficient. It must be said that ponding is not explicitly simulated,
that is, Code_Bright does not have a special element representing storage of water in a pond.
When one assumes no ponding, a very high value for γw can be used (but not to high to avoid
numerical instabilities). Then, if the soil is saturated (Pl > Pga) all rainfall that cannot infiltrate
will runoff.
The flux of water jw is the sum of rainfall P, evaporation E and advective flux of vapour gas
phase jgw and of surface runoff jsr:
j w = k rain P + k evap E + j gw + j sr (6)
where coefficients krain and kevap are input data passed through FLUX and may be used to disable
their respective flux.
Flux of energy
- Radiation
Several options are available to evaluate radiation (using ISUN which is passed to the
subroutine by FLUX(19)): 0 for horizontal plane, 1 for vertical cylinder and 2 for measured
atmospheric radiation (short + long wavelength radiations). SUN subroutine is called only 8 if
ISUN is lower or equal to 1.
The radiation Rn can be given as a measured data or it can be calculated, depending on the value
of ISUN:
Rn = (1 − Al ) R g + εRa − εσT 4 if ISUN ≤ 1
Rn = (1 − Al ) Rm − εσT if ISUN = 2
4
(7)
R = R else
n m
where Rg is the direct solar short wave radiation, Ra is the long wave atmospheric radiation, Al
is the albedo, ε is the emissivity, σ is the Stefan-Boltzman constant (5.67×10-8 J s-1 m-2 K-4) and
8
For other cases (ISUN = 2 or 3), RAD_DIR is set to values read from root_atm.dat by
subroutine get_atm_data.
32
Rm represents the values of measured radiation and read from file root_atm.dat by subroutine
get_atm_data.
Both the albedo and emissivity are considered function of the liquid saturation Sl:
(
Al = Ad + ( Ad − Aw ) Sl2 − 2Sl ) (8)
The calculation of the solar radiation Rg depends on the value of ISUN. Only the case of
ISUN=0 will be presented here. It takes into account the time of the day and the year according
to:
πR (t − t m + 0.5d s )π
R g = G sin if t m − 0.5d s ≤ t ≤ t m + 0.5d s
2d s ds (11)
R = 0 otherwise
g
where ds is the time span between sunrise and sunset and tm is the time at noon and RG is the
daily solar radiation calculated by an empirical relation:
t m + 0.5 d s
RG = ∫ R g dt = R A (0.29 cos λ + 0.52 I n ) (12)
t m − 0.5 d s
where λ is the latitude and In is the cloud index (= 1 for a clear sky, = 0 for a completely clouded
sky) and RA is the daily solar radiation in absence of atmosphere:
d πd
RA = S0 rs d cos λ cos δ sin s + d s sin λ sin δ (13)
π dd
where S0 is the solar constant (= 1367 J·m-2·s-1), rs is the relation between the average distance
between the earth and the sun and that of a given moment, dd is the duration of a day (= 86400 s)
and δ is the declination of the sun.
The values of ds, rs and δ are calculated as follows:
dd
ds = acos(− tan λ tan δ ) (14)
π
t − t0 t − t0
rs = 1.00011 + 0.034221cos 2π + 0.00128 sin 2π
da d a
(15)
t − t0 t − t0
+ 0.000719 cos 4π + 0.000077 sin 4π
d a d a
33
t − ts
δ = −δ max sin 2π (16)
d a
where da is the duration of a year (= 365.241 days = 3.15568×107 s), t0 is the time of January
1st, ts is the time when autumn starts (September 21st for the northern hemisphere) and δmax is
the maximum declination of the sun (= 0.4119 rad = 23.26°).
Subequations 2 and 3 in Eq. (7) use data read from root_atm.dat file. ISUN=2 suppose that
measured data corresponds to global radiation (short + long wavelength i.e. solar + atmospheric
radiations) whereas ISUN=3 assumes that measured data already are net radiation.
Results output
If values are surveyed at a given node, the following variables are written to files
200+nodout: t+∆t, P, E, jwg, jwl, jw, ja, Rn, Hs, Hc, je.
34
Subroutine get_atm_data
This subroutine computes atmospheric data at time t+dt either by simulation or by interpolation
of input data. In both cases, returned values are summarized in Table II.2.8. This table also
mentions implemented names for these variables and columns concerned in the matrix
atmosferic, where all data needed for simulations or interpolations are stored.
Two options can be used to compute atmospheric data: interpolation and simulation.
Interpolation uses a simple linear interpolation of the specified parameters versus time.
Simulation uses the following sinusoidal expression:
t − ta t − td
x(t ) = x m + x a sin 2π + x d sin 2π (20)
da dd
where x is the value of the parameter, xm is its mean value, aa is its annual amplitude, ad its daily
amplitude, ta is the start of the annual variation, td is the start of the daily variation, da is the
duration of a year and dd is the duration of a day (= 86400 s).
Table II.2.8: Atmospheric data taken into account in the boundary conditions module.
Atmospheric variable Unit Implemented name Columns used
Atmospheric temperature, Ta °C TEMP_ATM 1–2
Atmospheric gas pressure, Pga MPa PG_ATM 3–4
Relative humidity, Hr – RELHUM 5–6
Solar radiation, Rn J/m²/s RAD_DIR 7–8
Insolation fractions, In – FRACINS 9 – 10
Rain, P kg/m²/s RAIN 11 – 12
Wind velocity, va m/s WIND 13 – 14
Matrix atmosferic is read from file root_atm.dat (or root_atm.inp) which index is iin3
=103. This file is opened by read_assign_files subroutine (read_grid.f) and read by
read_atm_bc (read_general.f) subroutine which assigns atmosferic its values and is
called by main_initialize (code_bright_main.f).
Note that atmosferic dimensions also are read from iin3 and that this instruction is present
in read_assign_dim_opt_2 (code_bright_main.f).
Data simulation
If a simulation is performed, a sinus shape function with annual and daily variations is used.
The daily variation is only taken into account if the time increment is lower than one day. Input
data (for each variable) needed for each variable is (Unit represents the unit given in Table
II.2.8):
• annual mean (Unit), xm,
• annual amplitude (Unit), xa,
• annual gap (s), ta,
• daily amplitude (Unit), xd,
• daily gap (s), td.
35
For a given variable, the simulated value a time t+dt is obtained according to the following
relation 9:
1 2π (t + 0.5 dt − t a ) π dt
x(t + dt ) = xm + xa d a sin sin
π dt da da
(21)
2π (t + 0.5 dt − t d ) π dt
+ xd d d sin sin
dd d d
Figure II.2.6 and Figure II.2.7 show the simulation of the annual variation of an atmospheric
variable (case of temperature).
Data interpolation
Interpolation concerns all available data for discrete times ti satisfying t < ti+1 and t + dt > ti
until condition t + dt < ti+1 is satisfied. If code time t is lower (resp. bigger) than lowest (resp.
highest) discrete time, the subroutine is stopped.
Annual amplitude
20
Temperature (°C)
Annual gap
15
Annual mean
10
0
0 50 100 150 200 250 300 350 400
time (days)
Figure II.2.6: Simulation of annual variation of an atmospheric variable (case of temperature).
9
This expression cancels out daily variations if the time increment dt is higher than dd,
duration of one day.
36
Annual variation with and without daily variations
29
with
28
without
27
26
25
Temperature (°C)
24
23
22
21
20
19
18
170 175 180 185 190 195 200
time (days)
Figure II.2.7: Simulation of annual variation of an atmospheric variable (case of temperature;
close view).
Subroutine SUN
This routine calculates the direct solar radiation. Calculation type (ISUN) is passed through
FLUX vector in atmosferic_boundary_condition. The different values for ISUN
that imply a call to SUN subroutine are:
• 0 - horizontal plane,
• 1 - vertical cylinder.
However, ISUN=0 is the unique option proposed in the manual of Retraso.
This calculation takes into account sun distance and declination (function of date from 1st
January), solar day duration (function of latitude and declination). All this data together with
insolation fraction allow computing daily radiation (equations where presented in section 0).
If time increment is bigger than a day, direct radiation directly uses daily radiation. Otherwise,
time with respect to night is taken into account and the direct radiation is calculated.
Figure II.2.8 shows the simulated daily radiation versus time in both cases.
Case ISUN=1
ds
R A = S 0 rs 1 + 1 − (cos λ cos δ + sin λ sin δ )
2
(22)
2
37
30
25
Daily Radiation (MJ/m²)
20
15
10
0
0 50 100 150 200 250 300 350 400
Time (days)
Reference
Chow, V. T., Maidment, D. R., and Mays, L. W. (1988). Applied Hydrology, McGraw-Hill.
Initial values of the unknowns can be assigned on surfaces/volumes on the geometry. A constant
or linear distribution is available.
If distribution is linear, information about unknowns’ values at final point and the coordinates
of the initial and final points are required.
In case of nodes with multiple initial conditions assigned, the ones assigned into entities of
higher levels prevail. It is recommended to assign the materials on the geometry entities, but it
is also possible to assign them directly into mesh elements in case is needed.
38
II.2.3.5 Initial porosity
A constant initial value of porosity can be assigned on surfaces/volumes on the geometry.
Porosity value should be less than 1.
In chapter VI, the description of history variables required for elastoplastic and viscoplastic
models is included.
If distribution is linear, information about stresses and history variables values at final point
and the coordinates of the initial and final points are required.
39
y’=y’’ y
z x’
α
x
x’’
β
y
x’
x
Cross-anisotropy plane
β
z=z’ z’’
Figure II.2.9: Convention of reference axis for transverse isotropic material.
40
Assigning priorities
Conditions assigned on the geometry are distributed over the mesh with priorities. In general
points have priority over lines, lines over surfaces and surfaces over volumes. At mesh level,
nodes have priority over elements.
Mechanical boundary conditions on high entities are superimposed when they are applied to a
lower entity. This means, for instance in a 2-D case, that a point that belongs to two lines will
have the combination of boundary conditions coming from these two lines.
Assign condition
A condition is assigned to the entities with the given field values. If assigning from the
command AssignCond, the option Change allows the definition of the field values. Do not
forget to change these values before assigning. Option DeleteAll erases all the assigned
entities of this particular condition. Conditions can be assigned both on the geometry and on
the mesh but it is convenient to assign them on the geometry and the conditions will then be
transferred to the mesh. If conditions are assigned on the mesh, any remeshing will cause the
conditions to be lost.
Conditions that are to be attached to the boundary of the elements, are assigned to the elements
and GiD searches the boundaries of the elements that are boundaries of the total mesh. Option
Unassign inside AssignCond, permits to unassign this condition. It is also possible to
unassign from only certain entities.
If a mesh has already been generated, for any change in the condition assignments, it is
necessary to re-mesh again.
Draw condition
Option Draw all draws all the conditions assigned to all the entities in the graphical window.
This means to draw a graphical symbol or condition number over every entity that has this
condition. If one particular condition is selected, it is possible to choose between Draw and one
of the fields. Draw is like Draw all but only for one particular condition. If one field is
chosen, the value of this field is written over all the entities that have this condition assigned.
When the condition has any field referred to the type of axes, the latters can be visualized by
means of Draw local axes.
Unassign condition
In window mode, command UnAssign lets the user to choose between unassigning this
condition from the entities that owe it or unassigning all the conditions or select some entities
to unassign. In command mode UnAssing, do it for all the conditions. For only one condition,
use command Delete All (see section Assign condition).
Entities
Create an information window with all entities assigned including values at every one.
41
II.2.4. INTERVALS DATA
Intervals are a way to change some conditions and, eventually, material properties.
Properties for materials can vary at each interval or remain constant. For a problem with several
intervals, a window with the following fields will appear for each constitutive law:
ITYCL P1 P2 P3 P4 P5 P6 P7 P8 P9 P10
In which a different line for each interval will appear. Usually, only the first line should be
filled. The following lines may be left in blank if material properties are unchanged.
Boundary conditions may vary at each interval or remain constant. For conditions the
correct way to proceed is to define all the invariable conditions first (i.e. those that
remain unchanged during all intervals). Then, it is allowed to define as many intervals
as desired with the command NewInterval or update the conditions in different
intervals using the command ChangeInterval. It is possible to define as many
conditions as necessary into particular intervals. The conditions which have not been
duplicated when creating new intervals, are only considered for their interval.
Interval data parameters, decribe temporal limits and time steps for each interval. They can be
entered with the command IntervalData or in the intervals data window. If entered in a
window, the data is not accepted until button Accept is pressed. This data can be entered
before or after meshing.
Description
INTERVAL DATA
Units of time discretization Time units for defined interval. Options:
Seconds, Minutes, Hours, Days, Weeks,
Months, Years
Initial Time (start period) Initial time for defined interval.
(TIMEI in root_gen.dat)
Initial Time Step Initial time step for this time interval.
(DTIME in root_gen.dat)
Final Time (end period) Final time for defined interval.
(TIMEF in root_gen.dat)
Partial Time Time from which time increment is kept constant.
(TIME1 in root_gen.dat)
Partial Time Step Constant time step value.
(DTIMEC in root_gen.dat)
Put displacements to 0 Yes / No
This option put displacement to zero at the
beginning of the next time step.
42
In 'Writing results frequency', the intermittence for writing results is defined, i.e. only after
a given number of time steps the results will be written. This may cause inconveniences
if the user desires the results at precisely fixed times (for instance: 6 months, 1 year, 2
year, etc.). Moreover, if something changes between two runs (e.g. boundary conditions)
and any time increment should be modified, the value of the times in which results are
output will not be identical between these two runs. In this case, it would be difficult to
make a comparison of the two analyses because the output results correspond to
different times.
However, it is possible to decide the values of the times for output using a sequence of
consecutive intervals. In this way, the results will be output for all 'Final time' defined,
and if the user is only interested in these fixed times a very large value may be used for
'Writing results frequency' to avoid output at other times.
_________________________________
43
III.CODE_BRIGHT. PROCESS
III.1. CALCULATE
This part deals with the stage of the process that solves the numerical problem. The system
would allow calling the Finite Element program without necessity of leaving from the work
environment. Pressing Calculate the user can see a Process window, and clicking on Start
the solver module runs.
The option Calculate from the beginning is useful in case that you need to calculate
the project from the beginning (without using the saved files created in a previous calculate
process). WARNING: Calculate from the beginning option will erase all the save
type files from the working directory.
44
Card 3. Dimensions and options
Variables: NZ1, NZ2, MFRONTH, NDF, MNVAL, ISOLVE
Format: (10I5). It is not required if free format is used
NZ1: =MXDIFN: maximum difference between connected nodes, this variable
is read for dimensioning purposes. The node numeration of the grid is
assumed to have been optimised in order to reduce the matrix band width.
If θ=ε=0 are used in a non-mechanical problem, then MXDIFN can be 0
because a quasi-explicit approximation will be used, i.e. only a NDF-
diagonal matrix is solved which contains derivatives of the storage terms.
(See below for NZ=NZ1*NZ2).
NZ2: =MBANDT: total band width (geometrical for 1 variable), (MBANDT =
2(MXDIFN+1)-1, the user should provide a value but the code checks this
value. So this entry is redundant.
NZ= Used only for ISOLVE=5. It is the number of nonzero-blocks in the
NZ1*NZ2: jacobian (i.e. the number of nonzeros for NDF=1). This variable is
computed as NZ=NZ1*NZ2. Since this variable is checked internally, if
the number of nonzeros is not known a priori, a guess can be used and the
code automatically checks its validity. Otherwise, the required value is
output.
MFRONTH: void
NDF: Number of degrees of freedom per node. For instance a 2-dimension
thermomechanical analysis requires NDF=3.
MNVAL: Maximum number of integration points in an element (default=1). For a
two-dimensional analysis with some (not necessarily all) quadrilateral
elements, MNVAL=4. For a three-dimensional analysis with some (not
necessarily all) quadrilateral prism elements, MNVAL=8.
ISOLVE: Solve the system of equations according to different algorithms.
ISOLVE=3: LU decomposition + backsubstitution (NAG
subroutines, fonts available). (recommended option for direct
solution).
ISOLVE=5: Sparse storage + CGS (conjugate gradients squared)
46
IOPTPC =-3, krl-element and krg-element are computed by averaging nodal
values. Derivatives of relative permeabilities are also averaged.
IOPTPC =-4, krl-element and krg-element are set equal to the maximum
nodal value.
IOPTPC = 1: capillary pressure is used (Pc = Pg-Pl) as state variable instead
of Pl. If IOPTPC is = 1 then it is necessary to use IOPTXAL=1 and
IOPTXWG=1, and IOPTDISPL=0 and IOPTTEMP=0. That is,
IOPTPC=1 is only available for two phase immiscible fluids.
IOPTHYS: = 1: option for hysteretic behaviour of retention curve
IUPDC: = 1: updated lagrangian method, i.e., co-ordinates are modified after each
time increment is solved. If deformations are very large, some elements
may distort. If distortion is very large the volume of an element may
become negative and the execution would terminate immediately.
Remarks: vapour and air dissolved are considered automatically depending on options in Card
5. However, if for any reason they want not to be considered, then the auxiliary indexes
IOPTXWG=1 or IOPTXAL=1 can be used.
Card 8. Constants.
Variables: EPSILON, THETA, PGCONS, TCONS, PLCONS
Format: (6F10.0). It is not required if free format is used
EPSILON: Position of intermediate time tk+ε for matrix evaluation, i.e. the point
where the non-linear functions are computed. (frequent values: 0.5, 1)
THETA: Position of intermediate time tk+θ for vector evaluation, i.e. the point
where the equation is accomplished. (frequent values: 0.5, 1)
PGCONS: Constant gas phase pressure for solving with IOPTPG=0, otherwise
ignored. (frequent value: 0.1 MPa = atmospheric pressure).
TCONS: Constant temperature for solving with IOPTTEMP=0, otherwise this
value is ignored.
PLCONS: Constant liquid phase pressure for solving with IOPTPL=0, otherwise
ignored. (if PLCONS is greater than -1.0 x 1010 then wet conditions are
assumed for computing viscous coefficients in creep laws. (Otherwise the
medium is considered dry.)
47
Card 9. Void.
This line should be left blank.
48
ITIME (see table): 0No time step prediction
1Time step prediction according to a limit of 4 iterations.
2Time step prediction according to a limit of 3 iterations.
3Time step prediction according to a limit of 2 iterations.
4Time step prediction according to a limit of 1 iteration.
6A new time step is predicted from the relative error in variables
of the previous time increment. If the relative error is less than
dtol=0.01, time increment is reduced according to error deviation.
7 The same as 6, but with dtol=0.001.
8 The same as 6, but with dtol = 0.0001.
9 The same as 6, but with dtol = 0.00001.
TIME STEP FACTOR IS 1.4 ALWAYS
ITIME = 0
DTIMEC is the upper bound of time step
NUMBER OF NR ITERATIONS AS A MEASURE OF ERROR
0.25
4
ITIME = =
1 f ≥ 0.5
iter
0.25
3
ITIME = =
2 f ≥ 0.5
iter
0.25
2
=
ITIME =3 f ≥ 0.5
iter
0.25
1
=
ITIME =4 f 1.05 ≥ 0.5
iter
RELATIVE ERROR CONTROL
0.5
First Order Approach = DTOL
0.1 ≤ f 0.8 ≤ pdfmax
error
ITIME = 6 DTOL = 0.01
ITIME = 7 DTOL = 0.001
ITIME = 8 DTOL = 0.0001
ITIME = 9 DTOL = 0.00001
In the upper Table, f is the factor for time step reduction and pdfmax
is set to 1.4.
IMBACKUP: 0 No Backup.
1 Backup only for the last time step.
2 Backup for all time steps.
IWRALL: 1 Write all information for output.
0 Write partial information for output.
IPLOYFILES: 1 Write in separated output files (select variables in output
window). Two files are generated by each variable selected:
ROOT_variable.post.msh and ROOT_variable.post.res
0 If IWRALL=1
49
CardGroup 11. Convergence parameters
Variables:
Displacements: DELMXU, FACU, DELFMX, DUMX
(Omit this line if IOPTDISPL=0)
Liquid pressure: DELMXPL, FACPL, DELQWMX, DPLMX
(Omit this line if IOPTPL=0)
Gas pressure: DELMXPG, FACPG, DELQAMX, DPGMX
(Omit this line if IOPTPG=0)
Temperature: DELMXT, FACT, DELQMX, DTMX
(Omit this line if IOPTTEMP=0)
Inclusions conc.: DELMXI, FACI, DELIMX, DIMX
(Omit this line if IOPTXWS=0)
Format: (5F10.0). It is not required if free format is used
Each computed unknown requires a line with its associated parameters. In this way each
equation has different tolerances. If IOPTDISPL=1, only one line with DELMXU, FACU,
DELFMX, DUMX should be read regardless whether the problem is one, two or three
dimensional.
50
FACPG: Maximum (relative) gas pressure error tolerance (-)
DELQAMX: Maximum nodal air mass balance error tolerance (kg/s)
DPGMX: Maximum gas pressure correction per iteration (MPa) (time increment is
reduced if necessary).
DELMXT: Maximum (absolute) temperature error tolerance (oC)
FACT: Maximum (relative) temperature error tolerance (-)
DELEMX: Maximum nodal energy balance error tolerance (J/s)
DTMX: Maximum temperature correction per iteration (oC) (time increment is
reduced if necessary).
DELMXI: Maximum (absolute) water in inclusion mass fraction error tolerance (-)
FACI: Maximum (relative) water in inclusion mass fraction error tolerance (-)
DELIMX: Maximum nodal inclusions balance error tolerance (kg/s)
DIMX: Maximum mass fraction in solid correction per iteration (-) (time
increment is reduced if necessary)
Relative error is defined as the ratio between variable correction (δx) and variable increment
(∆x).
Convergence criteria are as follows (only convergence on the equation of energy balance is
illustrated, but the same applies for the other equations):
If (δT < DELMXT + FACT.T) for all nodes, then convergence has been achieved (condition
A). T is the value of the variable temperature.
If (qh < DELEMX) for all nodes (qh represents here the energy balance or residual in a node),
then convergence has been achieved (condition B).
It should be mentioned that convergence in terms of δT and convergence in terms of qh should
be reached simultaneously because the Newton - Raphson is used. For this reason the program
stops the iteration process and looks for another time step when one of the two conditions (A
or B) is achieved. For instance if the user decides that convergence should be imposed because
the residual has reached a tolerance then, DELMXT and FACT should be set to very low values.
When more than one degree of freedom is solved per node and the last option is used,
convergence in terms of variable or residual should be achieved by all the variables
simultaneously. In other words, it is not possible that the mechanical problem converges by
residual and the thermal converges by the variable.
If (δT > DTMX), time increment will be reduced. This parameter controls the accuracy of the
solution in terms of how large can be the time increments. A low value of DTMX will force to
small time increments when large variations of temperature take place.
Usually, it is difficult to guess the values of the tolerances that should be used in a problem.
The convergence criterion in terms of absolute terms is linked with the unknowns, and hence it
also depends on the range of variation of the variable. In this case the user decides the degree
of accuracy that is needed for each variable. The tolerances in relative terms are usually larger
than the values for absolute.
51
Finally, the tolerance values for residual convergence are more difficult to guess because 'a
priori' it is difficult to know the values of forces or flows equilibrating at nodes. Again the user
should reach a compromise between a very strict value or a less severe condition.
Convergence parameters for Conjugate Gradients Squared method of solution (Omit this
CARD if ISOLVE is not equal to 5).
Variables: DXS,DRS,DRSREL
Format: (5F10.0). It is not required if free format is used
This Card is only required for ISOLVE=5.
DXS: Maximum abs. correction for solver (usually a very low value)
DRS: Maximum abs. residual for solver (< min(DELFMX, DELQWMX, DELQAMX,
DTMX, DELIMX)) assuming all them > 0
DRSREL: Maximum relative residual for solver. The solver residual is normalised with the
RHS of the system of equations to be solved.
This group of Cards ends with ' -1' (forma I5).
The following group of Cards, beginning with time period definition can be repeated several
times to define periods or steps with different material properties and boundary conditions. For
the first step all information should be read and for the subsequent steps only modifications are
required.
52
Time increments during time step [TIMEI1, TIMEF] are adapted by the code according to flag
control ITIME (see Card 10). This may cause inconveniences if the user desires the results at
precisely fixed times (for instance: 6 months, 1 year, 2 year, etc.). Moreover, if something
changes between two runs (e.g. boundary conditions) and any time increment should be
modified, the value of the times in which results are output will not be identical between the
two runs. In this case, it would be difficult to make a comparison of the two analyses because
we would not have the same times for output.
A first way to overcome this inconvenience is to prescribe an upper bound for the time
increment, reflected in the variable DTIMEC. If convergence requires time increments smaller
than DTIMEC, time increment is reduced. But, if convergence is easy and the current time
increment becomes higher than DTIMEC, it is fixed to DTIMEC. Variable TIME1 allows for
setting an intermediate time between TIMEI and TIMEF from which the upper bound for the
increment becomes active, as represented below:
Another way to set fixed times for output results is to use a sequence of Cards number 13
separated by two (only flow or only mechanical problem) or three (flow and mechanical
problem) lines with '-1' (format I5) indicating that nothing changes in the new time period,
except the time discretization. In this way, results will be output for all TIMEF's, and if the user
is only interested in these fixed times a very large value may be used for INTER (see Card 10)
to avoid output at other times.
Example:
-1 indicates no change in material properties
-1 indicates no change in mechanical boundary conditions
-1 indicates no change in flow boundary conditions
350000. 0.0 0.0 10000. 360000. 86400.
-1 indicates no change in material properties
-1 indicates no change in mechanical boundary conditions
-1 indicates no change in flow boundary conditions
360000. 0.0 0.0 10000. 370000. 86400.
-1 indicates no change in material properties
-1 indicates no change in mechanical boundary conditions
-1 indicates no change in flow boundary conditions
in this case for the times 350000, 360000 and 370000 the results would be written. Time step
in this case would be lesser or equal than 10000.
It is is possible to define at the beginning of the calculation a step for equilibration of the initial
stress state. This is done by defining a time step starting from a negative value (TIMEI <0) and
ending at 0 (TIMEF = 0). During this step, gravity is applied as a ramp. Greater is time step
(TIMEF – TIMEMAX), smoother is the gravity ramp.
53
Card 14. Number of material
Variables: IMAT
Format: (I5). It is not required if free format is used
IMAT: index of material (<= NUMMAT)
(if '-1' (format I5) is read, no more materials are read, and hence, parameters will be zero (or
default values when defined) or the value read in a former time period)
54
Card 17. Type of Boundary Condition (Mechanical Problem)
Variables: IF
Format: (I5) . It is not required if free format is used
(Omit Card 17 if IOPTDISPL=0)
IF: index of boundary condition (IF <= NFDTYPE)
(if ' -1' (format I5) is read, no more boundary condition types are expected)
TIT for FD1 Value for FD1 TIT for FD$ Value for FD$
TIT for FD2 Value for FD2 TIT for FD$ Value for FD$
TIT for FD3 Value for FD3 TIT for FD$ Value for FD$
TIT for FD4 Value for FD4 TIT for FD$ Value for FD$
... etc. according to NPFD=2*(NDIM*(NDIM+2)+1).
A20 F10.0 A20 F10.0
FD1: x direction force applied fxo FD17 ∆fxo obtained as ramp loading
during TIMEI and TIMEF
FD2 y direction force applied fyo FD18 ∆fyo obtained as ramp loading
during TIMEI and TIMEF
FD3 z direction force applied fzo FD19 ∆fzo obtained as ramp loading
during TIMEI and TIMEF
FD4 displacement rate, first direction u1o FD20
FD5 displacement rate, second direction u2o FD21
FD6 displacement rate, third direction u3o FD22
FD7 cos(α1), first direction FD23
55
FD8 cos(β1), first direction FD24
FD9 cos(γ1), first direction FD25
FD10 cos(α2), second direction FD26
FD11 cos(β2), second direction FD27
FD12 cos(γ2), second direction FD28
FD13 cos(α3), third direction FD29
FD14 cos(β3), third direction FD30
FD15 cos(γ3), third direction FD31
FD16 γ FD32 index
For a one dimensional problem the general boundary condition is applied by means a force
computed as:
f x = f xo + γ cos(α 1 )(u10 − u x )∆t
where u1 is the computed displacement along the first direction. Obviously, for a one-
dimensional problem cosα1 can only be equal to zero or one.
For a two dimensional problem the general boundary condition is applied by means a force
computed as:
f x = f xo + γ cos(α 1 )(u10 − u1 )∆t + γ cos(α 2 )(u 20 − u 2 )∆t
f y = f yo + γ cos(β 1 )(u10 − u1 )∆t + γ cos(β 2 )(u 20 − u 2 )∆t
where:
u1 = u x cos(α 1 ) + u y cos(β 1 )
u 2 = u x cos(α 2 ) + u y cos(β 2 )
A very large value of γ can be used to impose a fixed displacement rate. If γ is insufficiently
large, the prescription of the displacement will be inaccurate. On the contrary, extremely large
values can cause matrix ill conditioning. Each specific problem requires an adjusted value.
If index is equal to 0, the values of forces calculated above are directly incorporated at the nodal
force balance. If index is equal to 1 then, the forces are considered stresses on the boundary,
and therefore the forces to be applied at nodes are internally obtained by the product with the
lateral areas of elements.
For three dimensional problems, for instance, it is possible to prescribe displacement rate for
three different directions, without any other restriction. In this way, any kind of displacement
boundary condition (ex: displacement zero along a direction 45 degrees with respect to the
vertical, etc) can be imposed. For a constant force applied on the boundary, the three
components along x,y,z axes should be given.
This is a loop for IF=1, NFDTYPE. For each IF, I=1, NPFD. This variable is (NPFD=5
NDIM+1) the number of parameters for mechanical boundary condition.
The last Card of this group must be always ' -1' (format I5) regardless of the number of types
read.
56
This group of Cards (Card 17 and CardGroup 18) (mechanical boundary conditions) only
exists if the mechanical problem is solved. For each time period only the types that change need
to be read.
TIT for FL1 Value for FL1 TIT for FL1 Value for FL1
TIT for FL2 Value for FL2 TIT for FL2 Value for FL2
... ...
TIT for FL20 Value for FL20 TIT for FL20 Value for FL20
A20 F10.0 A20 F10.0
57
FL7 ωl w prescribed mass fraction of solute FL28
(kg/kg)
FL8 ωl a prescribed mass fraction of air FL29
(kg/kg)
FL9 jl prescribed liquid flow rate (units in FL30 increment of jl during
Table II.2.1) time step (units in
Table II.2.1)
FL10 Pl prescribed liquid pressure (MPa) FL31 increment of Pl during
time step (MPa)
FL11 γl (see comments for negative value; FL32
units in Table II.2.1)
FL12 βl (units in Table II.2.1) FL33
FL13 ρl prescribed liquid density (kg/m3) FL34
FL14 je prescribed heat flow rate (units in FL35 increment of je during
Table II.2.1) time step (units in
Table II.2.1)
FL15 T prescribed temperature (C) FL36 increment of T during
time step (C)
FL16 γe (units in Table II.2.1) FL37
FL17 λe (positive values): FL38
[je = je*exp(-abs(λe) t)] (units: 1/s)
λe (negative values): [je = jet-abs(λe)]
FL18 FL39
FL19 δ: parameter for smoothing curve the FL40
seepage (outflow of water only)
boundary condition. For a positive
value a parabolic curve is used; for a
negative value an exponentially
decaying curve is used. δ is the
distance from the reference pressure
to the point of change
FL20 index: auxiliary index. FL40
index = +1.0 means that all flow rates
are nodal values
index = -1.0 means that all flow rates
are per unit volume (3-D), area(2-D)
or length (1-D) of medium (internal
source or sink)
index = +2.0 means that all flow rates
are per unit area (3-D) or length (2-
D) (lateral fluxes).
58
The boundary condition is incorporated by adding a flux. The mass flux of species i=w as a
component of phase α=g (i.e. the inflow or outflow of vapour) is:
( ) 0 dt
( ) dt
(
j gw = ω gw j g0 + ∆j g0 + ω gw γ g Pg0 + ∆Pg0 − Pg + β g ρ g ω gw
∆t
0
∆t
) − (ρ ω )
0
g
w
g
where the superscript ()o stands for the prescribed values, dt is the current time increment and
∆t the current time step. Terms ∆(.) dt/∆t allow for imposing a linear variation of the variable
(.) during the time step. Mass fraction and density prescribed are only required for inflow
because for outflow the values in the medium are automatically considered.
Positive values of mass flow rate indicate injection to the medium.
This general form of boundary condition, includes three terms. The first one is the mass inflow
or outflow that takes place when a flow rate is prescribed at a node. The second term is the mass
inflow or outflow that takes place when a phase pressure is prescribed at a node. The coefficient
γ is a leakage coefficient, that is, a constant that allows to prescribe a pressure with more or less
strength. If γ is large pressure will tend to reach the prescribed value. However, an extremely
large value can produce matrix ill conditioning and a lower one can produce inaccuracy in
prescribing the pressure. However, it is not difficult to guess adequate values for a given
problem simply by trial. The third term is the mass inflow or outflow that takes place when
species mass fraction is prescribed at a node.
A surface where seepage (only outflow for liquid phase is permitted) is possible has a boundary
condition of prescribed liquid pressure. However, only liquid outflow is permitted. To
recognize this fact, γl must be negative. This negative sign only indicates that nodes with this
kind of boundary condition allow seepage.
Another situation occurs when an internal source or sink should be imposed. In this case it is
preferable to use index = -1.0 and the program automatically considers that the nodal flows are
per unit volume and will be multiplied by the volume associated to the cell centered in the node.
If there is inflow of gas or liquid phase, it is very important to give values of the following
variables: ωlh, ωgw, ωla, ρl, ρg and T. Otherwise they are assumed zero which is not correct
because they will be far from the equilibrium. If outflow takes place, this is not relevant because
the values of the medium are used.
For energy the boundary condition has the general form:
dt dt
je = je0 + ∆je0 + γ e T 0 + ∆T 0 − T + E gw ( j gw ) + ...
∆t ∆t
in other words, the last terms imply that mass inflow and outflow through the boundary induces
energy transfer.
In general, this is a loop for IF=1, NFLUXTYPE. For each IF, I=1,NPFLUX.
The last Card of this group must be always ' -1' (format I5).
This group of Cards (Card 19 and CardGroup 20) only exists if any balance (water, air,
energy flow) problem is solved. For each time period only the types that change need to be
read.
59
III.4. GEOMETRICAL DESCRIPTION FILE: ROOT_GRI.DAT
Card 1. Grid writing index
Variables: IOWGRI, IOFILE, IFMT
Format: (5I5). It is not required if free format is used
IOWGRI: =1, a ROOTMSH.DAT file is created on output
IOFILE: =1, four (4) names are read in the following four lines for files containing, respectively,
FILE1: nodes (Card 2),
FILE2: connectivities (Card 3),
FILE3: initial conditions (Card 4) and
FILE4: element-wise variables (stresses and porosities) (Card 5).
IFMT: =1, to read connectivities according to old format
60
MNNEL is the maximum number of nodes that may have a possible element in the finite
element grid that is used in a problem. With the elements that are implemented at present the
following values are internally assigned to MNNEL: for NDIM=1 is MNNEL=2, for NDIM=2
is MNNEL=6, and for NDIM=3 is MNNEL=8.
61
Format: (16I5). It is not required if free format is used
NOUTOT: Number of nodes for which time evolution is required
IVOU: Variable required at these nodes. IVOU can range from 1 (first
unknown) to NDF (last unknown), and from NDF+1 (first nodal
dependent variable (DEPVARN vector)) to NDF+NDVN (last nodal
dependent variable).
INTERNODE: Frequency for output (=1 implies all time steps).
62
IELVOUT: Variable required at these elements. IELVOUT can range from 1 to 2
(DEPVARE vector, i.e. degree of saturation and/or porosity), or, from -
1 to -6 (stress vector) and from -7 to (-7 - nhv/2) (history variables).
INTERELEM: Frequency for output (=1 implies all time steps).
63
III.5. SUMMARY-LIST OF CARDS
This section contains the list of Cards with the variables that are read in each one.
File ROOT_GEN.DAT
Card 1. Problem HEAD
Card 2. Dimensions and NUMNP, NUMEL, NDIM, IAXISYM, NUMMAT, NHV
options
Card 3. Dimensions and MXDIFN, MBANDT, MFRONTH, NDF, MNVAL,
options ISOLVE
Card 4. Dimension NFDTYPE, NFLUXTYPE
boundary conditions
Card 5. Options. IOPTDISPL, IOPTPL, IOPTPG, IOPTTEMP, IOPTXWS
Unknowns to be calculated
Card 6. Other options IOPTXHL, IUPDPOR, IOPTXWG, IOPTXAL, IOPTPC,
IOPTHYS, IUPDC
Card 7. Flags. Auxiliary IFLAG1, IFLAG2, IFLAG3, IFLAG4, IFLAG5
options
Card 8. Constants EPSILON, THETA, PGCONS, TCONS, PLCONS
Card 9. Void
Card 10. Options IOWIT, INTER, ITERMAX, IOWCONTOURS,
ITERMAXS, ITIME, IMBACKUP, IWRALL,
IPOLYFILES
CardGroup 11. DELMXU, FACU, DELFMX, DUMX (Omit if
Convergence parameters IOPTDISPL=0)
DELMXPL, FACPL, DELQWMX, DPLMX (Omit if
IOPTPL=0)
DELMXPG, FACPG, DELQAMX, DPGMX (Omit if
IOPTPG=0)
DELMXT, FACT, DELQMX, DTMX (Omit if
IOPTTEMP=0)
DELMXI, FACI, DELIMX, DIMX (Omit if IOPTXWS=0)
DXS,DRS,DRSREL (Omit if ISOLVE not equal 5)
This group ends with -1
Card 12. Gravity GRAVITY(1), ..., GRAVITY(NDIM)
Card 13. Period time TIMEI, DTIME, TIME1, DTIMEC, TIMEF, FACTTIME
variables
Card 14. Number of IMAT
material
64
Card 15. Number and name ICL, TIT, ITYCL
of constitutive law
CardGroup 16. Parameters TIT, PARCL(1,ICL,IMAT) TIT, PARCL(6,ICL,IMAT)
constitutive law TIT, PARCL(2,ICL,IMAT) TIT, PARCL(7,ICL,IMAT)
TIT, PARCL(3,ICL,IMAT) TIT, PARCL(8,ICL,IMAT)
TIT, PARCL(4,ICL,IMAT) TIT, PARCL(9,ICL,IMAT)
TIT, PARCL(5,ICL,IMAT) TIT, PARCL(10,ICL,IMAT)
(group of Cards from IMAT=1 to NUMMAT and for every IMAT value from ICL=1 to
NCL (not all ICL are required) )
This group ends with -1 (ICL loop)
This group ends with -1 (IMAT loop)
Card 17. Type of boundary IF
condition (Mechanical
problem)
CardGroup 18. TIT, FORDISP(1,IF)
Force/displacement TIT, FORDISP(2,IF)
prescribed TIT, FORDISP(...,IF)
(group of Cards from IF=1 to NFDTYPE ) (Omit if IOPTDISPL=0)
This group ends with -1
Card 19. Type of boundary IF
condition. Mass or heat
transport problems
CardGroup 20. Flux TIT, FLUX(1,IF), TIT, FLUX(21,IF)
problem boundary condition TIT, FLUX(2,IF), TIT, FLUX(22,IF)
TIT, FLUX(...,IF, TIT, FLUX(...,IF))
TIT, FLUX(20,IF), TIT, FLUX(40,IF)
(group of Cards from IF=1 to NFLUXTYPE ) (Omit if IOPTPL + IOPTPG +IOPTTEMP =
0)
This group ends with –1
The group of Cards from 13 to 20 can be repeated in order to make a simulation with several
time periods in which the boundary conditions and material properties are not the same. If any
parameter is not read, the value in the previous interval is used. If a '-1' is read with IMAT, ICL
and IF, then no change takes place in material properties and boundary conditions.
65
File ROOT_GRI.DAT
Card 1. Grid writing index IOWGRI, IOFILE, IFMT
CardGroup 2. Node co- N, COORD(1, N), ..., COORD(NDIM, N), IFORDISP(1,N),
ordinates and boundary IFORDISP(2,N), IFORDISP(3,N), IFLUXTYPE(1,N),
condition type IFLUXTYPE(2,N), IFLUXTYPE(3,N), WIDTH(N)
(group of Cards from N=1 to NUMNP)
CardGroup 3. Node L, MTYPE, LTYPE, KXX(1,L),..., KXX(MNNEL,L)
connectivities, material,
element type,...
(group of Cards with L=1 to NUMEL)
CardGroup 4. Initial N, XOLD(1,N), ..., XOLD(NDF,N)
values of unknowns
(group of Cards with N=1 to NUMNP)
CardGroup 5. Initial L, STRESSOLD(1, 1, L), ..., STRESSOLD(NSTREC, 1, L)
values of stresses
(group of Cards with L=1 to NUMEL) (Omit if IOPTDISPL=0)
CardGroup 6. Other L, POROSITY(L), (FK(I, L), I=1,NDIM), ANISOTPER(1,
element wise properties L), ..., ANISOTPER(NISOT, L), THICKNESS (L),
(FK(I,L), I=NDIM+1, NDIM+3)
(group of Cards with L=1 to NUMEL)
Card 7. Time evolution of NOUTOT, IVOU(1), ..., IVOU(10), INTERNODE
state or dependent variables
at nodes
Card 8. Nodes for time NODOUT(1), ..., NODOUT(NOUTOT)
evolution
Card 9. Piezometric head IWHEAD, NWHEAD
map
Card 10. Nodal flows IWNFLOW
Card 11. Time evolution of LOUT, IELVOUT(1), ..., IELVOUT(10),
dependent variables at INTERELEMENT
elements
Card 12. Element numbers NELOUT(1), ..., NELOUT(LOUT)
for time evolution of
element-wise variables
_________________________________
66
IV. CODE_BRIGHT. POSTPROCESS.
IV.1. POST-PROCESSING TOOLS: SHORT DESCRIPTION
The different results that the system allows to display are the following ones:
• Geometry: GiD displays the whole volumetric mesh, surface sets and boundary surfaces. It
can also cut and divide them, in its original state as well as also in its deformed state, and
switch them on and off. GiD displays how the meshes/sets will be deformed according to a
certain vectorial variable. GiD provides two representations: on the first one all the results
will then be drawn on these original or deformed meshes (Main Geometry); on the second
one, a superimposed representation (Show Geometry) is provided, that can be also deformed,
but the results will still be drawn on the main representation. The user can compare the main
original/deformed meshes/sets with the second representation, which can also be deformed
with the same or another vector deformation, or with no deformation at all. The scaling of
all the displays can be modified interactively.
• Show minimum and maximum: The minimum and maximum values of the variable for
the currently viewed meshes/sets can appear, pointing all the nodes where these limits are
computed, in dark blue the minimum values and in red the maximum ones.
• Vectors: GiD presents a vector distribution according to the vectorial or matrix variables on
each node, showing their magnitudes and directions. The scaling of the vectors can be
modified interactively.
• Contour fields: GiD represents the variables through isosurfaces or contours that comprise
all the values between two given values. GiD takes advantage of the graphical capabilities
of the machine, allowing a smoothing of the results when a high number of colors is used.
• Contour lines: This representation is quite similar to the last one, but the uniform bands are
substituted here by isolines, where each one ties several points with the same value.
67
The file root.post.msh contains information about the mesh, like the coordinates of the nodes,
and the coordinates of the elements, while the file root.post.res presents the information of the
results and the variation in time of different variables in the nodes or Gauss points. To create
the program it is important to take into account the format of these files, specially the
root.pos.res file. Note that with root we mean the root name of the project: for example, if our
project is DAM.gid then the post process files will be DAM.post.res and DAM.post.msh.
The first lines of the file root.pos.res are the header, which contains information about the Gauss
points. Based on the type of element employed, this header may have 5 or 10 lines: 5 lines for
triangular elements and 10 lines for quadrilateral elements. See Figure IV.3.1.
The results of the process appear in the next lines. The results are divided by time, which means
that are presented for variable 1 in time 1, then variable 2 in time 1, until the last variable, and
then start again for all the variables in time 2 and so on.
The line 6 −or 11, depending on the element employed− presents the header for the results of
variable 1 in time 1, which has the following format:
Result “Variable” Isochrones “time” “type” “set”
Variable: Ex. Temperature, liquid pressure, stress
Type: Vector, scalar or matrix.
Set: Nodes (onNodes) or Gauss points (GP) depending of the variable.
For example: Result Temperature Isochrones 0.549584E-06 Scalar onNodes
The next line has the text “Values”. In the following lines 2 columns for scalar type variables
appear, 3 columns for variable type vector or 6 columns for matrix type, as it is shown in Figure
IV.3.2. The first column presents the number of the node or Gauss point, whose coordinates are
presented in the file root.post.msh. After the last value, the line “end values” appears and in the
succeeding line the header for the next variable appears.
68
Result format for a scalar variable.
If the user needs to create another file type root.post.res to visualize the results in the GID
interface, it is important to consider the format presented before, and take into account that the
header of the values has the following format:
‘Result’,1x,a15,1x,’Isochrones `,e12.6,1x,a15,1x,a30
69
V. CODE_BRIGHT. THEORETICAL ASPECTS
In porous media, subjected to thermal, hydraulic, and mechanical conditions, relevant thermo-
hydro-mechanical (THM) phenomena takes place. In fact, there exist a number of mutual
interactions that must be taken simultaneously into account in analyses. For instance, strains
due to thermal loading will induce stress variations and changes in mass storage terms and
hydraulic conductivity. The thermal expansion of the water in the pores itself causes changes
in the degree of saturation or, if the material is saturated or quasi-saturated, increases of water
pressure. Thermal induced vapor diffusion and the dependence of water viscosity on
temperature also affect significantly the water transfer process.
On the other hand, changes in hydraulic conditions influence the temperature field via variations
of thermal conductivity and affect the stress/strain field due to pore water pressure and pore gas
pressure changes. Gas pressure is affected by the increase in vapour pressure with temperature.
This may lead to further changes in the pattern of gas and water flow. Finally, porosity changes
due to volumetric strain influence pore pressure distributions because of associated variations
in storage terms and hydraulic conductivity. The effect on temperature is less important as the
variations of thermal conductivity with porosity are relatively small. In Appendix 1 the most
significant interactions between the various phenomena are presented in a systematic manner.
An unavoidable consequence of all those phenomena interacting simultaneously is the need to
carry out coupled THM analysis in which all the main aspects of the problem can be considered
in an integrated way. Such a formulation and the numerical approach adopted to solve the
governing equations are presented in the following sections.
Solid phase
Liquid phase:
water +
70
The three phases are:
• solid phase (s) : mineral
• liquid phase (l) : water + air dissolved
• gas phase (g) : mixture of dry air and water vapour
The three species are:
• solid (-) : the mineral is coincident with solid phase
• water (w) : as liquid or evaporated in the gas phase
• air (a) : dry air, as gas or dissolved in the liquid phase
The following assumptions and aspects are taken into account in the formulation of the problem:
• Dry air is considered a single species and it is the main component of the gaseous
phase. Henry's law is used to express equilibrium of dissolved air.
• Thermal equilibrium between phases is assumed. This means that the three phases are
at the same temperature
• Vapour concentration is in equilibrium with the liquid phase. Psychrometric law
expresses its concentration.
• State variables (also called unknowns) are: solid displacements, u (three spatial
directions); liquid pressure, Pl; gas pressure, Pg; and temperature, T.
• Balance of momentum for the medium as a whole is reduced to the equation of stress
equilibrium together with a mechanical constitutive model to relate stresses with
strains. Strains are defined in terms of displacements.
• Small strains and small strain rates are assumed for solid deformation. Advective terms
due to solid displacement are neglected after the formulation is transformed in terms
of material derivatives (in fact, material derivatives are approximated as eulerian time
derivatives). In this way, volumetric strain is properly considered.
• Balance of momentum for dissolved species and for fluid phases are reduced to
constitutive equations (Fick's law and Darcy's law).
• Physical parameters in constitutive laws are function of pressure and temperature. For
example: concentration of vapour under planar surface (in psychrometric law), surface
tension (in retention curve), dynamic viscosity (in Darcy's law), strongly depend on
temperature.
71
V.2. GOVERNING EQUATIONS
The governing equations for non-isothermal multiphase flow of water and gas through porous
deformable saline media have been presented by Olivella et al. (1994). A detailed derivation is
given there, and only a brief description is included here.
The equations that govern this problem can be categorised into four main groups. These are:
balance equations, constitutive equations, equilibrium relationships and definition constraints.
Equations for mass balance were established following the compositional approach. That is,
mass balance is performed for water, air and salt species instead of using solid, liquid and gas
phases. Equation for balance of energy is established for the medium as a whole. The equation
of momentum balance for the porous medium is reduced to that of stress equilibrium.
72
Mass balance of solid
Mass balance of solid present in the medium is written as:
𝜕𝜕(𝜌𝜌𝑠𝑠 (1 − 𝜙𝜙)) (1)
+ ∇ ∙ (𝐣𝐣𝑠𝑠 ) = 0
𝜕𝜕𝜕𝜕
where θs is the mass of solid per unit volume of solid and js is the flux of solid. From this
equation, an expression for porosity variation was obtained as:
𝐷𝐷𝑠𝑠 𝜙𝜙 (1 − 𝜙𝜙) 𝐷𝐷𝑠𝑠 𝜌𝜌𝑠𝑠 𝑑𝑑𝐮𝐮 (2)
= + (1 − 𝜙𝜙)∇ ∙
𝐷𝐷𝐷𝐷 𝜌𝜌𝑠𝑠 𝐷𝐷𝐷𝐷 𝑑𝑑𝑑𝑑
The material derivative with respect to the solid has been used and its definition is:
𝐷𝐷𝑠𝑠 (∎) 𝜕𝜕(∎) 𝑑𝑑𝐮𝐮 (3)
= + ∙ ∇(∎)
𝐷𝐷𝐷𝐷 𝜕𝜕𝜕𝜕 𝑑𝑑𝑑𝑑
Equation (2) expresses the variation of porosity caused by volumetric deformation and solid
density variation.
73
Mass balance of air
Once the other mass balance equations have been written it is straightforward to obtain the mass
balance of air taking into account that air is the main component of the gas phase and that it
may be also present as air dissolved in the liquid phase.
𝜕𝜕�𝑒𝑒𝑔𝑔 𝜌𝜌𝑠𝑠 (1 − 𝜙𝜙) + 𝑒𝑒𝑔𝑔 𝜌𝜌𝑔𝑔 𝑆𝑆𝑔𝑔 𝜙𝜙 + 𝑒𝑒𝑙𝑙 𝜌𝜌𝑙𝑙 𝑆𝑆𝑙𝑙 𝜙𝜙� 𝜙𝜙𝜙𝜙𝑔𝑔 𝑝𝑝𝑔𝑔 𝜕𝜕𝜌𝜌𝑔𝑔
− +
𝜕𝜕𝜕𝜕 𝜌𝜌𝑔𝑔 𝜕𝜕𝜕𝜕
+𝛻𝛻 ∙ �𝐢𝐢𝑐𝑐 + 𝐣𝐣𝑒𝑒𝑒𝑒 + 𝐣𝐣𝑒𝑒𝑒𝑒 + 𝐣𝐣𝑒𝑒𝑒𝑒 � = 𝑓𝑓 𝑄𝑄
(8)
𝜕𝜕�ℎ𝑠𝑠 𝜌𝜌𝑠𝑠 (1 − 𝜙𝜙) + ℎ𝑔𝑔 𝜌𝜌𝑔𝑔 𝑆𝑆𝑔𝑔 𝜙𝜙 + ℎ𝑙𝑙 𝜌𝜌𝑙𝑙 𝑆𝑆𝑙𝑙 𝜙𝜙� 𝜕𝜕𝑝𝑝𝑔𝑔
− 𝜙𝜙𝜙𝜙𝑔𝑔 +
𝜕𝜕𝜕𝜕 𝜕𝜕𝜕𝜕
+𝛻𝛻 ∙ �𝐢𝐢𝑐𝑐 + 𝐣𝐣ℎ𝑠𝑠 + 𝐣𝐣ℎ𝑔𝑔 + 𝐣𝐣ℎ𝑙𝑙 � = 𝑓𝑓 𝑄𝑄
where 𝐢𝐢𝑐𝑐 is energy flux due to conduction through the porous medium, the other fluxes (𝐣𝐣𝑒𝑒𝑒𝑒 , 𝐣𝐣𝑒𝑒𝑒𝑒 ,
𝐣𝐣𝑒𝑒𝑒𝑒 , 𝐣𝐣ℎ𝑠𝑠 , 𝐣𝐣ℎ𝑙𝑙 , 𝐣𝐣ℎ𝑔𝑔 ) are advective fluxes of energy or enthalpy caused by mass motions and 𝑓𝑓 𝑄𝑄 is
an internal/external energy supply (heat supply). In this case this term accounts, for instance,
energy dissipation due to medium deformation which is not explicit because it is negligible in
most cases. The use of the material derivative allows obtaining an equation formally similar to
the mass balance of water. The reason for the similarity is that both water and internal energy,
are considered present in the three phases.
Hence, only one equation is required which expresses the balance of internal energy in the
porous medium as a whole. The enthalpy equation is obtained using the definition of enthalpy
as ℎ = 𝑒𝑒 + 𝑝𝑝𝑝𝑝. The density time derivative or the pressure time derivative terms in energy or
enthalpy balance are sometimes neglected.
The fluxes in the divergence term include conduction of heat and advection of heat caused by
the motion of every species in the medium. A non-advective mass flux causes an advective heat
74
flux because a species inside a phase moves and transports energy. Contrary to what happens
with the movement of a contaminant in a groundwater system, the diffusive term for heat
transport (conduction of heat) is much larger than the term concerning hydromechanical
dispersion (non-advective flux caused by the velocity of fluids). For this reason, this term is
sometimes neglected.
Constitutive equations
Darcy's law liquid and gas advective flux 𝐪𝐪𝑙𝑙 , 𝐪𝐪𝑔𝑔
Equilibrium restrictions
Henry's law Air dissolved mass fraction 𝜔𝜔𝑙𝑙𝑎𝑎
The constitutive equations establish the link between the independent variables (or unknowns)
and the dependent variables. There are several categories of dependent variables depending on
the complexity with which they are related to the unknowns. The governing equations are
finally written in terms of the unknowns when the constitutive equations are substituted in the
balance equations.
Another type of relationships that relate dependent variables with unknowns are the equilibrium
restrictions. They are obtained assuming chemical equilibrium for dissolution of the different
species (air and vapour) in phases (liquid, gas). This assumption is sufficiently adequate
because these chemical processes are fast compared to the transport processes that take place
in porous media and, for this reason, they are not rate controlling.
75
V.2.3. Boundary conditions
Application of the Green's theorem to the divergence term (both in the balance or equilibrium
of stresses equations) produces terms which represent fluxes or stresses across or on the
boundaries. These terms are substituted by nodal flow rates or forces in the discretized form of
the equations. For the mechanical problem, the classical approach is followed to impose
external forces. Imposing displacements is made by means of a Cauchy type boundary
condition, i.e. a force computed as the stiffness of a spring times the displacement increment.
The boundary conditions for balance equations are incorporated by means the simple addition
of nodal flow rates. For instance, the mass flow rate of water as a component of gas phase (i.e.
vapour) is:
( ) ( ) γ (P ) (
− Pg + β g ρ g ω wg ) − (ρ ω ) (9)
0 0 0
j gw = ω wg j g0 + ω wg g g
0
g
w
g
where the superscript ()0 stands for prescribed values. This general form of boundary condition,
includes three terms. The first one is the mass inflow or outflow that takes place when a flow
rate of gas (jg0) is prescribed. The second term is the mass inflow or outflow that takes place
when the gas phase pressure (Pg0) is prescribed at a node. The coefficient γg is a leakage
coefficient, i.e., a parameter that allows a boundary condition of the Cauchy type. The third
term is the mass inflow or outflow that takes place when vapour mass fraction is prescribed at
the boundary. This term naturally comes from the nonadvective flux (Fick's law). Mass fraction
and density prescribed values are only required when inflow takes place. For outflow, the values
in the medium are considered. For the energy balance equation, the boundary condition has a
similar form.
76
V.3. NUMERICAL APPROACH
V.3.1. Introduction
The system of PDE's (Partial Differential Equations) is solved numerically. The numerical
approach can be viewed as divided into two parts: spatial and temporal discretizations. Finite
element method is used for the spatial discretization while finite differences are used for the
temporal discretization. The discretization in time is linear and the implicit scheme uses two
intermediate points, tk+ε and tk+θ between the initial tk and final tk+1 times. Finally, since the
problem presented here is non-linear, the Newton-Raphson method was adopted to find an
iterative scheme.
k
j
em
en
i
Once the solid balance is substituted in the other balance equations, computation of porosity at
an intermediate point is not necessary because its variation is expected to occur at slow rates.
For this reason, porosity is integrated explicitly, that is, the values at tk are used. Since the
variation of porosity is expressed by the solid mass balance equation, this assumption leads also
to some advantages for the iterative scheme. After the spatial discretization of the partial
differential equations, the residuals that are obtained can be written (for one finite element) as:
ru d u a u b u 0 (10)
rPl d Pl a Pl b Pl 0
d
rP = d P + a P + b P = 0
g dt g g g
rT d T a T b T 0
where r are the residuals, dd/dt are the storage or accumulation terms, a are the conductance
terms, and b are the sink/source terms and boundary conditions. After time discretization a more
compact form can read as:
k +1 d k +1 − d k (11)
r( X )= + A ( X k + ε ) X k + θ + b( X k + θ ) = 0
∆t k
77
where k is the time step index, : X=[(ux,uy,uz,Pl,Pg,T)(1), ..., (ux,uy,uz,Pl,Pg,T)(n)], is the vector of
unknowns (i.e. a maximum of seven degrees of freedom per node), A represents the
conductance matrix. The Newton-Raphson scheme of solution for this non-linear system of
AE's is:
∂r ( X k +1 ) k +1,l +1 (12)
k +1 ( X − X k +1,l ) = −R ( X k +1,l )
∂X
where l indicates iteration.
In the present approach, the standard Galerkin method is used with some variations in order to
facilitate computations. General aspects related to numerical solution of hydrogeological
problems can be found in Huyakorn and Pinder (1983). As shown in the preceding section, in
the mass and energy balance equations the following terms may be distinguished:
• Storage terms. These terms represent the variation of mass or energy content and therefore,
they are calculated by means of variables such as degree of saturation, density, porosity,
mass fraction and specific energy.
• Advective fluxes. The advective fluxes caused by motion of fluids computed using Darcy's
law and, except for the coefficients, they are explicit in terms of pressure gradients.
• Nonadvective fluxes. These terms, computed through Fick's law, are proportional to
gradients of mass fractions which do not belong to the set of unknowns. Fourier's law is used
for the conductive heat flux and it expresses proportionality to temperature gradients.
• Volumetric strain terms. In fact, these terms are also storage terms. They are proportional to
∇⋅du/dt which is equivalent to the volumetric strain rate.
• Sink/source terms.
Each of these terms requires specific treatment. This is described in detail in Olivella et al
(1996).
In order to explain the treatment of the different terms and equations the following notation is
introduced:
• node i: node in a finite element mesh
• e1,e2,...,em: elements that contain node i, i.e. a cell centered in node i is composed by a
fraction of these elements. m is variable from node to node and it is not related to the number
of nodes per element.
• nem: number of nodes in element em. For example, nem=3 for triangles, nem=4 for
quadrilaterals, nem=4 for tetrahedrons, etc.
• (⋅)k: the quantity is computed at time tk of the temporal discretization. The same for tk+1, tk+ε
or tk+θ
• (⋅)em: the quantity is computed in element em. This means at the center of the element or, in
other words, using the average of nodal unknowns.
• (⋅)i: the quantity is computed in node i as a function of the unknowns in that node. \item--}
• (⋅)i,em: the quantity is computed in node i but with the material properties corresponding to
element em.
• Vem: volume of element em.
• ξi: shape function for node i.
78
V.3.2 Treatment of different terms
Treatment of storage terms
In this sub-section we refer to terms not related to volumetric strain or porosity variation. The
storage or accumulation terms are computed in a mass conservative approach (Allen and
Murphy, 1986; Celia et al., 1990; Milly, 1984). The conservative approach discretizes directly
the accumulation terms while the capacitative approach uses the chain rule to transform time
derivatives in terms of the unknowns. Milly (1984) proposes modifications of the capacitative
approach in order to conserve mass. It seems reasonable that the mass conservative approach
should give a more accurate solution than the capacitative approach.
Mass conservation in time is achieved if the time derivatives are directly approximated by a
finite difference in time. Finite element method for the space discretization conserves mass
(Milly, 1984).
A typical storage term is (from Eq. 5) the variation of water in the gas phase:
𝐷𝐷𝑠𝑠 �𝜔𝜔𝑔𝑔𝑤𝑤 𝜌𝜌𝑔𝑔 𝑆𝑆𝑔𝑔 � 𝜕𝜕�𝜔𝜔𝑔𝑔𝑤𝑤 𝜌𝜌𝑔𝑔 𝑆𝑆𝑔𝑔 � (13)
𝜙𝜙 ≅ 𝜙𝜙
𝐷𝐷𝐷𝐷 𝜕𝜕𝜕𝜕
where the material derivative with respect to the solid is approximated as an eulerian derivative
because the small strain rate assumption. The weighted residual method is applied to the
governing equations and, for node i, (13) is transformed into:
At this point of the development we assume that porosity is defined element-wise. An element-
wise variable (Voss, 1984) is space-constant over every element, but different from element to
element. We will use φemk for porosity in element em at time tk. Similarly, a cell-wise variable
(Voss, 1984) is space constant over the cell centered in the node. It would be very easy to
compute (14) if the time derivative could be computed in a cell-wise way, because one value
would be sufficient for node i and (14) would be transformed into a very simplified form.
However, the degree of saturation is not only a function of nodal unknowns but also of material
properties such as porosity or retention parameters. To overcome this difficulty, the time
derivative in (14) is computed from nodal unknowns but with material properties of every
element in contact with the node. Hence m values are necessary in node i. Obviously if part of
this time derivative is not material dependent (density and concentration are only function of
temperature and pressure) then the corresponding variables are only computed in the node. This
leads to a kind of modified cell-wise variables.
Making use of these approximations, we finally obtain, for example for any integral in (14):
𝑘𝑘+1 𝑘𝑘
𝜕𝜕�𝜔𝜔𝑔𝑔𝑤𝑤 𝜌𝜌𝑔𝑔 𝑆𝑆𝑔𝑔 � �𝜔𝜔𝑔𝑔𝑤𝑤 𝜌𝜌𝑔𝑔 𝑆𝑆𝑔𝑔 � − �𝜔𝜔𝑔𝑔𝑤𝑤 𝜌𝜌𝑔𝑔 𝑆𝑆𝑔𝑔 � (15)
𝑘𝑘 𝑖𝑖,𝑒𝑒𝑚𝑚 𝑖𝑖,𝑒𝑒𝑚𝑚
� 𝑁𝑁𝑖𝑖 𝜙𝜙 𝑑𝑑𝑑𝑑 ≅ 𝜙𝜙𝑒𝑒𝑚𝑚 � 𝑘𝑘+1 𝑘𝑘
� � 𝑁𝑁𝑖𝑖 𝑑𝑑𝑑𝑑
𝑒𝑒𝑚𝑚 𝜕𝜕𝜕𝜕 𝑡𝑡 − 𝑡𝑡 𝑒𝑒𝑚𝑚
where a simple finite difference is used for the time discretization. This approximation allows
us to make the space integration independently of the physical variables. Therefore,
computation of geometrical coefficients is necessary only once for a given finite element mesh.
The integral of the shape function over an element is equal to Vem/nem for the case of linear
shape functions. These geometrical coefficients are also called influence coefficients. Without
loss of generality, they can be computed either analytically or numerically. Finally, it should be
79
pointed out that this formulation gives rise to a concentrated scheme, which means that the
storage term in node i is only a function of unknowns in node i. This is clearly advantageous
from a computational point of view (Huyakorn and Pinder, 1983).
𝐤𝐤𝑘𝑘𝑟𝑟𝑟𝑟 (16)
− � (𝛻𝛻𝑡𝑡 𝑁𝑁𝑖𝑖 )𝜔𝜔𝑔𝑔𝑤𝑤 𝜌𝜌𝑔𝑔 𝐪𝐪𝑔𝑔 𝑑𝑑𝑑𝑑 = �� (𝛻𝛻𝑡𝑡 𝑁𝑁𝑖𝑖 )𝜔𝜔𝑔𝑔𝑤𝑤 𝜌𝜌𝑔𝑔 �𝛻𝛻𝑁𝑁𝑗𝑗 �𝑑𝑑𝑑𝑑� �𝑝𝑝𝑔𝑔 � −
𝑣𝑣 𝑣𝑣 𝜇𝜇𝑔𝑔 𝑗𝑗
𝐤𝐤𝑘𝑘𝑟𝑟𝑟𝑟
− �� (𝛻𝛻𝑡𝑡 𝑁𝑁𝑖𝑖 )𝜔𝜔𝑔𝑔𝑤𝑤 𝜌𝜌𝑔𝑔 𝜌𝜌 𝐠𝐠𝑑𝑑𝑑𝑑�
𝑣𝑣 𝜇𝜇𝑔𝑔 𝑔𝑔
where the subscript j indicates summation over element nodes. Pg is a node-wise (Voss, 1984)
variable, which means that it is defined by its nodal values and interpolated on the elements
using the shape functions. Generalised Darcy's law has been used to compute the flux of the
gas phase:
(∇P − ρgg)
kkrg (17)
qg = −
µg
g
where k is the tensor of intrinsic permeability, krg is the relative permeability of the gas phase,
µg is the dynamic viscosity of gas and g is a vector of gravity forces. For node i the volume v
over which the integrals in (16) have to be performed is composed by the elements e1, e2, ...,
em. In this way, the advective terms (16) represent the lateral mass fluxes to cell associated to
node i from contiguous cells. The pressure term is considered first. The contribution of element
em to the total lateral flux towards node i is approximated as:
𝐤𝐤𝑘𝑘𝑟𝑟𝑟𝑟 (18)
�� (𝛻𝛻𝑡𝑡 𝑁𝑁𝑖𝑖 )𝜔𝜔𝑔𝑔𝑤𝑤 𝜌𝜌𝑔𝑔 �𝛻𝛻𝑁𝑁𝑗𝑗 �𝑑𝑑𝑑𝑑� �𝑝𝑝𝑔𝑔 � ≈
𝑒𝑒𝑚𝑚 𝜇𝜇𝑔𝑔 𝑗𝑗
𝑘𝑘+𝜀𝜀
𝑘𝑘𝑟𝑟𝑟𝑟 𝑘𝑘+𝜃𝜃
≈ �𝜔𝜔𝑔𝑔𝑤𝑤 𝜌𝜌𝑔𝑔 � �� (𝛻𝛻𝑡𝑡 𝑁𝑁𝑖𝑖 )𝐤𝐤 𝑘𝑘𝑒𝑒𝑚𝑚 �𝛻𝛻𝑁𝑁𝑗𝑗 �𝑑𝑑𝑑𝑑 � �𝑝𝑝𝑔𝑔 �
𝜇𝜇𝑔𝑔 𝑒𝑒𝑚𝑚
𝑗𝑗
𝑒𝑒𝑚𝑚
where three different intermediate points may be used, one for the pressure (tk+θ), another for
the intrinsic permeability (tk) and yet another for the remaining coefficients (tk+ε) including the
relative permeability. The intrinsic permeability remains in the integral because it is a tensorial
quantity, but if its product with the shape function gradients is split, then its coefficients can be
taken off from the integral. It should be noticed that intrinsic permeability is handled explicitly
(i.e. evaluated at time tk) because it is a function of porosity structure, which we assume to vary
slowly. Since all physical variables can appear outside the integral because they are considered
element-wise, the integrals of products of shape function gradients are also considered influence
coefficients (Huyakorn et al., 1986). They have to be computed for each element, but only once
for a given mesh.
80
A similar approximation is used for the gravity term in (16). Evaluation of density element-
wise is convenient in order to balance correctly pressure gradients with gravity forces at element
level.
Treatment of nonadvective terms (diffusive/dispersive)
In the balance equation of node i we find, typically, the following diffusive term:
(19)
− � (𝛻𝛻𝑡𝑡 𝑁𝑁𝑖𝑖 )𝐢𝐢𝑔𝑔 𝑑𝑑𝑑𝑑 = �� (𝛻𝛻𝑡𝑡 𝑁𝑁𝑖𝑖 )𝜙𝜙𝜙𝜙𝜌𝜌𝑔𝑔 𝑆𝑆𝑔𝑔 𝐷𝐷𝑔𝑔𝑤𝑤 𝐈𝐈�𝛻𝛻𝑁𝑁𝑗𝑗 �𝑑𝑑𝑑𝑑� �𝜔𝜔𝑔𝑔𝑤𝑤 �
𝑗𝑗
𝑣𝑣 𝑣𝑣
where the subscript j indicates summation over the nodes. ωgw is considered a node-wise
variable. Fick's law has been used to compute the diffusive flux:
i g = −φτρ g S g DgwI∇ω gw (20)
where τ is a tortuosity parameter, Dgw is the molecular diffusion coefficient which is a function
of temperature and gas pressure and I is the identity matrix. The contribution of element em to
the total lateral diffusive flux towards node i is approximated as:
(21)
�� (𝛻𝛻𝑡𝑡 𝑁𝑁𝑖𝑖 )𝜙𝜙𝜙𝜙𝜌𝜌𝑔𝑔 𝑆𝑆𝑔𝑔 𝐷𝐷𝑔𝑔𝑤𝑤 𝐈𝐈�𝛻𝛻𝑁𝑁𝑗𝑗 �𝑑𝑑𝑑𝑑 � �𝜔𝜔𝑔𝑔𝑤𝑤 � ≈
𝑗𝑗
𝑒𝑒𝑚𝑚
𝑘𝑘+𝜀𝜀 𝑘𝑘+𝜃𝜃
≈ (𝜙𝜙𝜙𝜙)𝑘𝑘𝑒𝑒𝑚𝑚 �𝜌𝜌𝑔𝑔 𝑆𝑆𝑔𝑔 𝐷𝐷𝑔𝑔𝑤𝑤 � �� (𝛻𝛻𝑡𝑡 𝑁𝑁𝑖𝑖 )𝐈𝐈�𝛻𝛻𝑁𝑁𝑗𝑗 �𝑑𝑑𝑑𝑑 � �𝜔𝜔𝑔𝑔𝑤𝑤 �
𝑒𝑒𝑚𝑚 𝑗𝑗
𝑒𝑒𝑚𝑚
where various time intermediate points have been used similarly to what was explained for the
advective terms. The treatment of these diffusive terms also takes advantage of the fact that the
Newton-Raphson method is used to obtain the iterative scheme. We directly interpolate mass
fractions (e.g. ωgw) and compute gradients.
The dispersive term is treated in a similar way as the diffusive. In this case dispersivities are
element-wise dependent variables. In principle, the liquid and gas fluxes, used to compute the
dispersion tensor, are also computed element-wise.
81
(23)
𝑑𝑑𝐮𝐮
𝑡𝑡
𝐮𝐮𝑗𝑗𝑘𝑘+1 − 𝐮𝐮𝑗𝑗𝑘𝑘
� 𝑁𝑁𝑖𝑖 𝛼𝛼𝐦𝐦 𝐁𝐁 𝑑𝑑𝑑𝑑 ≈ 𝛼𝛼𝑒𝑒𝑘𝑘+𝜀𝜀 �� 𝑁𝑁𝑖𝑖 𝐦𝐦𝑡𝑡
𝐁𝐁𝑗𝑗 𝑑𝑑𝑑𝑑 � � �
𝑒𝑒𝑚𝑚 𝜕𝜕𝜕𝜕 𝑚𝑚
𝑒𝑒𝑚𝑚 𝑡𝑡𝑘𝑘+1 − 𝑡𝑡𝑘𝑘
where j indicates summation over element nodes, u is the vector of nodal displacements and Bj
is the j-submatrix of B.
( )
r σ k +1 = ∫ B tσ k +1dv − f k +1 = 0 (24)
v
where r(σk+1) represents the residual corresponding to the mechanical problem and σk+1 is the
stress vector. Matrix B (composed by gradients of shape functions) is defined in such a way
that stress is a vector and not a tensor. The body force terms and the boundary traction terms
are represented together by fk+1. The constitutive model relates stresses with strains, with fluid
pressures and with temperatures at a point in the medium. Only if elasticity is included, the total
strain rate is decomposed in the following way:
𝐡𝐡𝑘𝑘+1 = −𝐃𝐃𝑒𝑒 𝐁𝐁(𝐮𝐮𝑘𝑘+1 − 𝐮𝐮𝑘𝑘 ) + (𝛔𝛔𝑘𝑘+1 − 𝛔𝛔𝑘𝑘 ) + 𝐈𝐈�𝑝𝑝𝑔𝑔𝑘𝑘+1 − 𝑝𝑝𝑔𝑔𝑘𝑘 � + 𝐃𝐃𝑒𝑒 𝐈𝐈𝑎𝑎𝑠𝑠 (𝑠𝑠 𝑘𝑘+1 − 𝑠𝑠 𝑘𝑘 ) (26)
+ 𝐃𝐃𝑒𝑒 𝐈𝐈𝑎𝑎 𝑇𝑇 (𝑇𝑇 𝑘𝑘+1 − 𝑇𝑇 𝑘𝑘 )
where h is the residual of stresses at every point. If stress can be obtained in an explicit way
from (26), it is simply substituted in (24). However, when nonlinear models are introduced a
substitution of the differential or incremental forms is necessary.
For the mechanical problem, the approximations that should be made are different from the
ones used for flow problems. According to the numerical approximations proposed for the flow
problems (hydraulic and thermal), we would tend to use element-wise matrices. However, the
mechanical problem has some peculiarities which do not allow this kind of simplified treatment.
First, linear triangular elements (the simplest element in two dimensional analyses), which have
been proven to be very adequate for flow problems, should be avoided for mechanical problems.
This is because if the medium is nearly-incompressible (creep of rocks takes place with very
small volumetric deformation), locking takes place (not all displacements are permitted due to
element restrictions). Second, linear quadrilateral elements with element-wise variables (this is
equivalent to one integration point) lead to hour-glassing (uncontrolled displacement modes
appear).
82
In order to overcome these difficulties, the selective integration (B-bar) method is used. It
consists in using a modified form of matrix B which implies that the volumetric part of
deformation and the deviatoric part are integrated with different order of numerical integration
(Hughes, 1980). For linear quadrilateral elements, four integration points are used to integrate
the deviatoric part while one is used for volumetric strain terms. Although this approximation
is different from what is proposed for the flow problem, element-wise variables or parameters
are maintained (porosity, saturation,...). Stress is not element-wise and it must be computed at
the integration points.
83
V.4. THEORETICAL APPROACH SUMMARY
The governing equations include: stress equilibrium equations (1, 2 or 3 according to the
dimensions of the problem), mass balance equations (different species) and internal energy
balance equation for the medium as a whole (thermal equilibrium is assumed).
The stress equilibrium equations are a simplified form of the balance of momentum for the
porous medium. Mass balance of water, solid and air are established. Since the assumption of
equilibrium is made, the mass of each species as present in any phase (solid, liquid or gas) is
balanced for the porous medium as a whole. In this way, one equation for each species is
obtained. The equilibrium assumption implies that partition functions are required to compute
the fraction of each species in each phase.
Each partial differential equation is naturally associated to an unknown. These unknowns can
be solved in a coupled way, i.e., allowing all possible cross coupling processes that have been
implemented, or, on the contrary, any uncoupled problem to obtain a single unknown can be
solved.
The balance equations that CODE_BRIGHT solves are compiled here:
0 Unknown:
tensor of vector of displacements,
divergence + = vector 0
total stress body forces u=(ux,uy,uz)
0
∇⋅σ +b = 0
84
∂ internal enery in solid, total fluxes external supply Unknown:
+ divergence =
∂t liquid and gas phase of energy of heat tempera-
ture, T (o)
𝜕𝜕�𝑒𝑒𝑔𝑔 𝜌𝜌𝑠𝑠 (1 − 𝜙𝜙) + 𝑒𝑒𝑔𝑔 𝜌𝜌𝑔𝑔 𝑆𝑆𝑔𝑔 𝜙𝜙 + 𝑒𝑒𝑙𝑙 𝜌𝜌𝑙𝑙 𝑆𝑆𝑙𝑙 𝜙𝜙� 𝜙𝜙𝜙𝜙𝑔𝑔 𝑝𝑝𝑔𝑔 𝜕𝜕𝜌𝜌𝑔𝑔
− +
𝜕𝜕𝜕𝜕 𝜌𝜌𝑔𝑔 𝜕𝜕𝜕𝜕
+𝛻𝛻 ∙ �𝐢𝐢𝑐𝑐 + 𝐣𝐣𝑒𝑒𝑒𝑒 + 𝐣𝐣𝑒𝑒𝑒𝑒 + 𝐣𝐣𝑒𝑒𝑒𝑒 � = 𝑓𝑓 𝑄𝑄
The definition of a problem (which of the above described equations should be solved) is
achieved by means of a set of general options (IOPTDISPL, IOPTPL, IOPTPG, IOPTTEMP).
These general options indicate whether one equation is included or not. For instance a
mechanical problem would require IOPTDISPL=1 and the other indexes equal to 0. Other
secondary options allows to include or not any of the possible processes. Specific indexes are
used to decide if the solid is soluble (i.e. the medium is saline), if the air solubility in liquid
phase is taken into account and if vapour is considered. Vapour transfer can only be considered
if the thermal problem is solved.
85
V.5. FEATURES OF CODE-BRIGHT
The implementation of a coupled non-linear approach requires some especific developments
and approximations. In this section the main aspects of the numerical approximation are
reviewed. The program CODE_BRIGHT uses the finite element method to solve the coupled
equations presented above. The main features of the numerical approach are:
• Linear interpolation functions on segments, triangles, quadrilaterals, tetrahedrons,
triangular prisms and quadrilateral prisms (regular). Analytical integration is used
for segments, triangles and tetrahedrons. Numerical integration is used for
quadrilateral, arbitrary triangular prisms (6 points) and quadrilateral prisms (8
points). For the mechanical problem selective integration is used for quadrilateral
and quadrilateral prisms (this means that the volumetric part is integrated with a
reduced quadrature of 1 point). Finally, for all elements the flow equations are
solved using element-wise and cell-wise approximations. This approximation is
independent of the type of integration performed.
• Finite differences and implicit scheme are used for time integration. Two
intermediate points are defined between the two ends of the time interval (tk, tk+1).
One represents the point where the equation will be accomplished (tk+θ) and the other
is the point where the non-linear functions are computed (tk+ε). For instance ε=0 and
θ=1 states for a linearised problem with a fully implicit scheme of integration.
• Newton-Raphson method for solution of the non-linear system of algebraic
equations that results once the space and time discretizations are applied.
• LU decomposition and backsubstitution (non-symmetric matrix) or conjugate
gradients squared to solve the system of linear equations that result from the
Newton-Raphson application.
• Automatic discretization of time. Increase or reduction of time increment according
to convergence conditions or output requirements. Reduction of time increment may
be caused by: excessive variation of unknowns per iteration, excessive number of
iterations to reach convergence and correction larger than in the previous iteration.
The main features of the program CODE_BRIGHT are:
• Options that allow to solve uncoupled and coupled problems. For instance: Hydro-
mechanical, Thermo-mechanical, Hydro-thermal problems can be solved if the
physical situation requires one of these approaches.
• Types of analysis: One dimension (uni-axial confined strain and axi-symmetric).
Two dimension (plane strain and axi-symmetric). Three dimensions.
• Several element types.
• Constitutive laws: each law defined as a set of parameters. Different types of
relationships can be chosen in some cases.
• Boundary conditions:
• Mechanical problem: forces and displacement rate in any spatial
direction and at any node
• Hydraulic problem: mass flow rate of water and air prescribed and
liquid/gas pressure prescribed at any node
• Thermal problem: heat flow rate prescribed and temperature
prescribed at any node
86
• Convergence criteria: Tolerances for absolute and relative error independent for
each unknown. Tolerance for residual convergence of each problem (mechanical,
hydraulic, etc). The node under worse conditions is used to verify the convergence
condition.
• Forces/flows = ε→0
• Absolute variable correction = δx→0
• Variable correction / variable increment = δx/∆x→0
• Output options: Time evolution of variables in nodes or elements. The user should
decide 'a priori' the nodal or element variables that will be output at all times
(absolutely all computed times will be output for a few variables). Contour maps in
the solution domain. Nodal or element variables can be used to draw contour maps.
However, in the second case it is required to perform an interpolation that may be
difficult due to the lack of continuity of the element variables.
87
V.6. SOLUTION OF SYSTEM OF EQUATIONS IN CODE_BRIGHT
As described in the preceding chapters, CODE_BRIGHT is a FE-three dimensional code that
solves the following equations:
• equilibrium of stresses (displacements)
• mass balance of water (liquid pressure)
• mass balance of air (gas pressure)
• balance of energy (temperature)
Associated to every equation there is a variable, i.e. the unknown that is obtained by solving
the corresponding equation. The code can be used to solve problems that need only some of the
equations in the list. The equations are solved together in a monolithic way. Usually problems
require variable time stepping and an iterative method for the solution of non-linearities.
Usually the resulting system of equations is non-symmetric.
88
of non-zero values if NDF was equal to 1. With this storage mode, the vectors JA and IA can
maintain the same size as would be required for NDF=1.
89
x: = x 0 (intial guess);
r: = b − Ax;
r is an arbitrary vector such that (r , ~
~ r ) ≠ 0 (e.g. ~
r = r );
ρ = (r , ~
0 r ); β = ρ ; p = q = 0;
0 0
K is a preconditioner matrix;
for i = 0,1,2,...
u: = r + β i q;
p: = u + β i (q + β i p);
Solve q from Kq = p;
v: = Aq;
α i = ρi / ( ~
r , v );
q: = u − α i v;
Solve v from Kv = u + q;
u: = Av;
x: = x + α i v;
r: = r − α i u;
if x close enough to A −1b then quit;
ρ = (r , ~
i +1 r );
β i+1 = ρi +1 / ρi ;
end i
where it can be seen that two matrix - vector products should be performed per iteration, two
vector-vector products (indicated by (,)), and the remaining operations are vector updates.
90
V.7. APPENDIX 1. THERMO-HYDRO-MECHANICAL INTERACTIONS
In this Appendix the main interactions between the various thermo-hydro-mechanical processes
are presented in a synthetic and systematic manner.
THERMAL PHENOMENA
Heat storage
Effects from:
• Thermal phenomena
− Heat storage proportional to temperature
• Hydraulic phenomena
− Liquid flow modifies the amount of water and air present
− Gas flow modifies the amount of air and water present
− Phase changes modifies heat storage through the latent heat of vapour
• Mechanical phenomena
− Porosity changes modify the amount of space left for fluids
Heat conduction
Effects from:
• Thermal phenomena
− Heat conduction driven by temperature gradients (Fourier’s law)
• Hydraulic phenomena
− Liquid flow modifies thermal conductivity
− Gas flow modifies thermal conductivity
• Mechanical phenomena
− Porosity changes modifies thermal conductivity
Heat advection by liquid flow
Effects from:
• Hydraulic phenomena
− Heat transport by liquid flow
Heat advection by air flow
Effects from:
• Hydraulic phenomena
− Heat transport by gas flow
Heat advection by vapour flow
Effects from:
• Hydraulic phenomena
− Heat transport by vapour diffusion
− Heat transport by gas flow
Phase changes
Effects from:
• Thermal phenomena
− Vapour pressure affected by temperature (water phase diagram and
psychrometric law)
• Hydraulic phenomena
− Vapour pressure affected by liquid flow through suction changes
(psychrometric law)
− Vapour pressure affected by gas flow through suction changes (psychrometric
law)
91
HYDRAULIC PHENOMENA
Water storage
Effects from:
• Thermal phenomena
− Liquid density changes with temperature
− Vapour density changes with temperature
− Phase change modifies the amount of water in liquid and gas phases
• Hydraulic phenomena
− Liquid density changes with liquid pressure
− Vapour density changes with suction and gas pressure
• Mechanical phenomena
− Porosity changes affect the space available for liquid and gas
Air storage
Effects from:
• Thermal phenomena
− Gas density changes with temperature
− Amount of dissolved air changes with temperature
• Hydraulic phenomena
− Gas density changes with gas pressure
− Amount of dissolved air depends on gas pressure
• Mechanical phenomena
− Porosity changes affect the space available for liquid and gas
92
Gaseous air transfer
Effects from:
• Thermal phenomena
− Hydraulic conductivity affected by gas viscosity that increases with
temperature.
− Degree of saturation varies with temperature (thermal expansion and phase
changes)
− Temperature variations influence gas density
• Hydraulic phenomena
− Gas flow controlled by gas pressure gradients (Darcy’s law)
− Hydraulic conductivity affected by degree of saturation, in turn controlled by
the value of suction (retention curve)
• Mechanical phenomena
− Porosity changes affect the value of hydraulic conductivity
− Porosity changes vary the pore space volume available for gas
93
MECHANICAL PHENOMENA
Stress/strain field
Effects from:
• Thermal phenomena
− Thermal expansion of materials
− Dependence of constitutive laws on temperature
• Hydraulic phenomena
− Dependence of constitutive laws on suction
• Mechanical phenomena
− Stress/strain constitutive laws
_________________________________
94
VI. CODE_BRIGHT. CONSTITUTIVE LAWS
This chapter contains the different models available and the corresponding parameters required
by each model. The following constitutive laws are available
HYDRAULIC AND MECHANICAL CONSTITUTIVE MODELS
THERMAL CONSTITUTIVE
MODELS (a) ELASTICITY (b)
RETENTION CURVE NONLINEAR ELASTICITY (b)
INTRINSIC PERMEABILITY VISCOPLASTICITY FOR SALINE MATERIALS
LIQUID PHASE RELATIVE (b)
PERMEABILITY VISCOPLASTICITY FOR GRANULAR
GAS PHASE RELATIVE MATERIALS (b)
PERMEABILITY VISCOPLASTICITY - GENERAL (b)
DIFFUSIVE FLUXES OF MASS DAMAGE-ELASTOPLASTIC MODEL FOR
DISPERSIVE FLUXES OF ARGILLACEOUS ROCKS (c)
MASS AND ENERGY THERMOELASTOPLASTIC MODEL FOR SOILS
CONDUCTIVE FLUX OF HEAT (d)
BARCELONA EXPANSIVE MODEL FOR SOILS
(e)
CASM’s FAMILY MODELS (f)
PHASE PROPERTIES (a) EXCAVATION PROCESS (g)
SOLID PHASE PROPERTIES
LIQUID PHASE PROPERTIES
GAS PHASE PROPERTIES
95
HYDRAULIC AND THERMAL LAWS
ICL NAME ITYCL DESCRIPTION
6 Retention curve 1 Van Genuchten model
2 Linear model
4 Square law
9 Van Genuchten model with asymptotic branch to
negative capillary pressures
12 Van Genuchten model modified for FEBEX project
18 Van Genuchten model modified for freezing model
7 Intrinsic permeability 1 Kozeny's model
2 Exponential law
4 Kozeny's model for matrix + cubic law for discontinuity
(normal strain is used to calculate aperture)
5 Kozeny's model for matrix + cubic law for discontinuity
(volumetric strain is used to calculate aperture)
15 Same as 5 but with different relative permeability for
matrix and discontinuity
16 Barton's law – Joint element
14 Liquid phase relative 1 Van Genuchten model
permeability 5 Liquid perfectly mobile
6 Generalized power
8 Power with initial cut off
12 Van Genuchten model for freezing model
19 Gas phase relative 1 Default law
permeability 5 Gas perfectly mobile
6 Generalized power
11 Diffusive flux of vapor 1 Molecular diffusion of vapour or air
2 Molecular diffusion of vapour or air + tortuosity is
variable with gas pressure (Pg)
12 Diffusive fluxes of 1 Molecular diffusion of dissolved salt and dissolved air
dissolved salt and air
8 Dispersive fluxes of 1 Fick's law (mass flux) and Fourier's law (heat flux)
mass and energy
9 Conductive flux of heat 1 Thermal conductivity dependence on porosity
(1) (geometric weighted mean)
2 Thermal conductivity dependence on porosity
(weighted arithmetic mean)
3 Thermal conductivity dependence on porosity
(nonlinear function of porosity)
20 Conductive flux of heat 1 Dependence on degree of saturation
(2) 2 Dependence on degree of saturation
4 Dependence on degree of saturation
5 Dependence on water/ice content
6 Dependence on degree of saturation (Chen & Ledesma,
2009)
96
PHASE PROPERTIES
ICL NAME ITYCL DESCRIPTION
10 Solid phase properties 1 Solid specific heat, density and expansion coefficient
2 Solid specific heat, density and expansion coeff.
(variation of solid phase specific heat)
15 Liquid phase properties: 1 Liquid density (exponential variation)
density 2 Liquid density (linear variation)
4 Liquid density (CO2)
5 Liquid density (linear dependency of the thermal
expansion coefficient with temperature)
6 Liquid density adjustment for a wide temperature
range (-30ºC—300ºC)
16 Liquid phase properties: 1 Liquid viscosity
viscosity
17 Gas phase properties: 1 Dry air density. Law of ideal gases and Henry's law for
density dry air
2 Usually used to consider a second liquid phase instead
of the gas phase
3 Like ITYCL=1 but with user defined values for gas
molecular mass and Henry's constant
4 Gas density CO2
18 Gas phase properties: 1 Gas viscosity
viscosity 2 Gas viscosity (exponential law)
4 Gas viscosity CO2
97
RETENTION CURVE
98
EQUATIONS ITYCL=18: Freezing model:
-λ
1
Sl − S rl Pi − Pl 1-λ
S= = 1 + ;
Sls − S rl P
e
ρi T
=
Pi Pl − ρil Clausius Clayperon
ρl T + 273
where Pi is the ice pressure, ρi =
ρl 0.91, =
ρil 306 MPa,
T temperature in K
Ice volumetric fraction: Si = 1 − Sl
99
RETENTION CURVE (ICL=6). PARAMETERS FOR ITYCL=4 (square law):
P1 Po MPa Measured P at certain temperature
P2 Void
P3 Void
P4 Srl Residual saturation
P5 Sls Maximum saturation
RETENTION CURVE (ICL=6). PARAMETERS FOR ITYCL=9 (Van Genuchten model with
asymptotic branch that goes to negative capillary pressures):
P1 Po MPa Measured P at certain temperature
Surface tension at temperature in which Po was measured
P2 σo N m-1
(usually σo=0.072 N/m at 20ºC)
P3 λ Shape function for retention curve
P4 Srl Residual saturation
P5 Sls Maximum saturation
P6 f Used for the asymptotic branch
100
RETENTION CURVE (ICL=6). PARAMETERS FOR ITYCL=18 (freezing model):
P1 Po MPa Measured P at certain temperature
Surface tension at temperature in which Po was measured
P2 σo N m-1
(usually σo=0.072 N/m at 20ºC)
P3 λ Shape function for retention curve
P4 Srl Residual saturation
P5 Sls Maximum saturation
Parameter for porosity influence on retention curve:
P6 a
Po(φ)=Po exp(a(φο−φ)
Parameter for porosity influence on retention curve:
P7 b
λ(φ)= λ exp(b(φο−φ))
P8 Void
Reference porosity for porosity influence on retention
P9 φο
curve
Flag to indicate stress concept for use in the mechanical
model:
P10 i_stress
i_stress=0 : Net stress (σn = σ - Pi)
i_stress=-1: Bishop’s stress (σb = σ - Pi + Sr (Pi-Pl))
Srl and Sls are lower and upper bounds of saturation. Effective saturation Se is defined in such a
way that ranges between 0 and 1.
101
INTRINSIC PERMEABILITY
CODES in ICL=7 ITYCL=1,2,4,5,6,15,16
ROOT_gen.dat
DESCRIPTION Intrinsic permeability
EQUATIONS ITYCL=1: For a continuum medium (Kozeny’s model):
ϕ3 (1 − ϕo ) 2
=k ko ϕo : reference porosity
(1 − ϕ) 2 ϕo 3
k o : intrinsic permeability for matrix ϕo
(∇Pα − ρα g)
kk rα
which is used in Darcy’s law: qα = −
µα
where viscosity, density and relative permeability are defined in other laws.
ITYCL=2: Exponential law: =k k o exp {b ( φ − φo )}
𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚 𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑
𝐤𝐤𝑘𝑘𝑟𝑟𝑟𝑟 = 𝐤𝐤 𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚 𝑘𝑘𝑟𝑟𝑟𝑟 + 𝐤𝐤𝑓𝑓𝑓𝑓𝑓𝑓𝑓𝑓𝑓𝑓𝑓𝑓𝑓𝑓𝑓𝑓 𝑘𝑘𝑟𝑟𝑟𝑟
𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚 𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑
𝑘𝑘𝑟𝑟𝑟𝑟 is calculated normally, and 𝑘𝑘𝑟𝑟𝑟𝑟 = 𝑆𝑆𝑔𝑔 = (1 − 𝑆𝑆𝑙𝑙 )
102
ITYCL=16: Zero thickness element
Barton’s law: Longitudinal intrinsic permeability
2
e2 a 2 1 a : Opening of the joint
k=
l = 2.5 JRC : Joint Roughness Coefficient
12 JRC 12
Transversal intrinsic permeability, kt , is considered to be equal to the
continuum media.
For the retention curve (ICL=6), air entry value depends on joint aperture,
kl 0 σ
as: P = P0
kl σ0
b3
P8 a m Spacing of the fractures:=k kmatrix +
12a
Permeability of the matrix is obtained as usual in porous media.
Reference strain to calculate aperture variations:
P9 εo - ∆b = a∆ε = a (ε − ε o ) for ε > ε o
103
INTRINSIC PERMEABILITY (ICL=7). PARAMETERS FOR ITYCL=6
P1-P5 The same as in ITYCL=1
P6 γ MPa-1 Dilatability
P7 εB - Maximum volumetric strain
P9 Kc - Permeability at zero minor principal stress
104
LIQUID PHASE RELATIVE PERMEABILITY
(
k rl = S e 1 − (1 − S e1/λ ) )
λ 2
Liquid phase relative permeability (ICL=14). PARAMETERS FOR ITYCL=1 (Van Genuchten
model):
P1 Void
P2 Void
P3 λ Power
Residual saturation (default = same value as for
P4 Srl
retention curve)
Maximum saturation (default = same value as for
P5 Sls
retention curve)
Liquid phase relative permeability (ICL=14). PARAMETERS FOR ITYCL=5 (liquid perfectly
mobile): None.
Srl and Sls are lower and upper bounds of saturation. Effective saturation Se is defined in such a
way that ranges between 0 and 1. In principle, the same values Srl and Sls should be defined for
liquid and gas relative permeability and for retention curve. However, different values can be
used to define a maximum saturation of liquid with possibility of remnant gas flow or vice-
versa.
106
GAS PHASE RELATIVE PERMEABILITY
107
DIFFUSIVE FLUXES OF VAPOUR
108
DIFFUSIVE FLUXES OF DISSOLVED SALT AND AIR
109
DISPERSIVE FLUXES OF MASS AND ENERGY
110
CONDUCTIVE FLUX OF HEAT (1)
111
CONDUCTIVE FLUX OF HEAT 1 (ICL=9). PARAMETERS FOR ITYCL=1 (geometric
weighted mean):
P1 λdry W m-1 K-1 Thermal conductivity of the dry porous medium
Thermal conductivity of the water saturated porous
P2 λsat W m-1 K-1
medium
P3 (λsolid)o W m-1 K-1 Solid phase thermal conductivity (ignored if λdry, λsat >0)
P4 λgas W m-1 K-1 Gas phase thermal conductivity (ignored if λdry, λsat >0)
P5 λliq W m-1 K-1 Liquid phase thermal conductivity (ignored if λdry, λsat >0)
P6 a1 Ignored if λdry, λsat > 0
P7 a2 Ignored if λdry, λsat > 0
P8 a3 Ignored if λdry, λsat > 0
112
CONDUCTIVE FLUX OF HEAT 1 (ICL=9). PARAMETERS FOR ITYCL=3 (nonlinear
function of porosity):
P1 Void
P2 Void
P3 (λsolid)o W m-1 K-1 Solid phase thermal conductivity
P4 (λdry)o W m-1 K-1 Dry thermal conductivity for reference porosity
P5 (λsat) o W m-1 K-1 Saturated thermal conductivity for reference porosity
P6 Void
P7 Void
P8 Void
P9 φo Reference porosity
P10 n Power of porosity
Heat dispersion is defined in the constitutive law ICL=8, ITYCL=1 (Dispersive fluxes of mass
and energy).
113
CONDUCTIVE FLUX OF HEAT (2)
where A1 represents the value of λ for 𝑆𝑆𝑆𝑆 = 0, A2 the value of λ for 𝑆𝑆𝑆𝑆 = 1,
𝑆𝑆𝑆𝑆 ∗ the degree of saturation for which thermal conductivity is the average of
the two extreme values and b is a parameter. In some cases, A1 and A2 are
equivalent to λsat and λdry, respectively. A1 and A2 are introduced as λsat and
λdry in ICL=9, ITYCL=1: Conductive flux of heat (1).
114
CONDUCTIVE FLUX OF HEAT 2 (ICL=20). PARAMETERS FOR ITYCL=7 (S-shaped
function):
Note this option does not allow using porosity, degree of saturation or temperature
dependencies.
115
PHASE PROPERTIES
EQUATIONS
Thermal dilatation coefficient for grains should be equal to the bulk value if thermal expansion
of the porous medium does not produce porosity variations. Specific heat for water, air and salt
(not solid) are internal values.
116
LIQUID PHASE PROPERTIES. LIQUID DENSITY.
CODES in ICL=15 ITYCL=1,2,3,4,5,6
ROOT_gen.dat
DESCRIPTION Liquid density
EQUATIONS ITYCL=1: Exponential variation:
ρl = ρl 0 exp (β ( Pl − Plo ) + α T + γ ωlh )
ITYCL=2: Linear variation:
ρl = ρl 0 (1 + β ( Pl − Plo ) + α T + γ ωlh )
ITYCL=4: CO2
ρl = ρl 0 exp ( β ( Pl − Pl 0 ) + αT + γωlh )(1 + δωCO
l
2
)
Vϕ
δ = 1 − ρl
M CO2
Vϕ = ( 37.51 − 9.585 ×10 −2
T + 8.740 × 10−4 T 2 − 5.044 × 10−7 T 3 ) × 10−6 m3 /mol
(Garcia, 2003)
ITYCL=5: Linear dependency of the thermal expansion coefficient
with temperature:
ρl = ρl 0 exp( β( Pl − Plo ) + A (T ) T + γωlh )
A (T ) =−3.5 × 10 −6 T − 7.49 × 10 −5
117
LIQUID PHASE PROPERTIES (ICL=15). PARAMETERS FOR ITYCL=1, 2:
P1 ρlo kg m-3 Reference density (default = 1002.6 kg·m-3)
P2 β MPa-1 Compressibility (default = 4.5×10-4)
Volumetric thermal expansion coefficient for water
P3 α o -1
C
(default = -3.4×10-4)
P4 γ Solute variation (default = 0.6923)
P5 Plo MPa Reference pressure (default = 0.1)
118
LIQUID PHASE PROPERTIES. LIQUID VISCOSITY.
EQUATIONS B
ITYCL=1: µ l = A exp
27315
. + T
Remark: liquid and gas density and viscosity are not material dependents. For this reason,
values should be prescribed only once. If these are multiplied defined, the code will use the
values it reads first.
119
GAS PHASE PROPERTIES. GAS DENSITY.
EQUATIONS ITYCL=1: law of ideal gases and Henry's law for dry air (as
ITYCL=3 with Ma = 0.02895 and H = 10000 MPa)
ITYCL=2: ρ g =
a
(ρ )
a
g o exp( β ( Pg − Pgo ) + αT )
Usually used to consider a second liquid phase instead of the
gas phase (in that case do not consider vapour in gas phase)
ITYCL=3: law of ideal gases and Henry’s law for any dry gas
species (as ITYCL=1, but with user defined values for gas
molecular mass and Henry’s constant:
P M dgs
ωldgs = dgs
H Mw
where Pdgs is dry gas species pressure (air pressure in the
formulation), Mw is molecular mass of water and Mdgs is
molecular mass of dry gas species.
120
GAS PHASE PROPERTIES (ICL=17). PARAMETERS FOR ITYCL=3 (gases law with
modified molecular mass and Henry’s constant):
P1 M kg mol-1 Molecular mass
P2 H MPa Henry’s constant
GAS PHASE PROPERTIES (ICL=17). PARAMETERS FOR ITYCL=4 (law for CO2 with
the values adjusted by Spycher et al., 2003):
P1 MCO2 kg mol-1 Molecular mass of CO2 (0.044 kg mol-1)
P2 H MPa Henry’s constant
Alternative CO2 density function at
P8 T=320 K adjusted from Span and Wagner
(1996) tables. To use it write: -777
121
GAS PHASE PROPERTIES. GAS VISCOSITY.
EQUATIONS ITYCL=1:
A 273 + T 1
µg =
B bk
1 + 1+
273 + T Pg
bk = C − Dk
(k : intrinsic permeability)
B
ITYCL=2: µ g = A exp
273.15 + T
Tc ρc
a10 = 0.248566120 a11 = 0.004894942
a20 = −0.373300660 a21 = 1.22753488
a30 = 0.363854523 a31 = −0.774229021
a40 = −0.0639070755 a41 = 0.142507049
122
GAS PHASE VISCOSITY (ICL=18). PARAMETERS FOR ITYCL=2:
P1 A MPa s Pre-exponential parameter
o Exponential parameter (only used if
P2 B C
A=B=0, but not used if A>0 and B=0)
Gas phase properties can be used to consider a second liquid in the case of a two immiscible
phase flow problem in a porous medium. In this case, water vapour and air dissolved must not
be considered, hence, VAPOUR NOT PERMITTED and DISSOLVED AIR NOT
PERMITTED should be used to avoid the species to be mixed.
123
CODE_BRIGHT. CONSTITUTIVE LAWS
This chapter contains the different models available and the corresponding parameters required
by each model. The following constitutive laws are available:
124
LINEAR ELASTICITY (Mechanical data 1)
CODES in ICL=1 ITYCL=1, 2, 16
ROOT_gen.dat
DESCRIPTION Elastic parameters (linear elasticity model). Linear elasticity with
parameters E and ν. Young modulus can be variable.
EQUATIONS dE
ITYCL = 1: E= Eo + (φ − φo ) ≥ Emin E
dφ
ITYCL = 2: Bi-linear elasticity model → Ec
εv limit εv
PARAMETERS
LINEAR ELASTICITY (Mechanical data 1). ITYCL=1
P1 E MPa Young Modulus
P2 dE/dφ MPa Variation of Young modulus with porosity
P3 ν - Poisson’s Ratio
P4 φo - Reference porosity
P5 φ min - Minimum porosity
P6 Emin MPa Minimum elastic modulus
-
P7 b Biot coefficient (default value = 1)
P3 ν - Poisson’s Ratio
Volumetric strain limit to change elastic modulus
P4 εv limit -
(positive value for an open gap)
125
LINEAR ELASTICITY (Mechanical data 1). ITYCL=16 (Elasticity – Zero thickness element)
P1 m MPa Model parameter
P2 Ks MPa/m Shear stiffness
P3 E MPa/m Out of plane stiffness
P4 amin m Minimum aperture
P5 a0 m Initial aperture of the joint
Note: Joint element should be created in GiD using the option: Geometry/create/contact
surface.
Contact surfaces are defined as being between two lines that are physically in the same place,
but which have different line and point entities. Choose the Contact surface option from the
menu, and then select some lines on both bodies. Using contact surface entities is like a meshing
specification. In this way, equal meshes will be generated for the two lines, ensuring a one-to-
one relationship between nodes. You can also select no mesh for the contact entity. This makes
it possible to have exactly the same mesh for both lines but without any additional element
between them.
PARAMETERS
LINEAR ELASTICITY – TEMPERATURE AND SUCTION. ITYCL=1
P1 as MPa-1 Swelling coefficient for changes in suction
P2 void
o -1
P3 bs C Linear thermal expansion coefficient for the medium
P4 Void
P5 Void
126
NONLINEAR ELASTICITY (Mechanical data 1)
127
PARAMETERS
Suction micro is used as history variable and therefore the model cannot be combined with
viscoplastic models.
128
NONLINEAR ELASTICITY. ITYCL=5 (with two independent coupling terms)
129
ANALOGY WITH BBM
where:
κi ( s ) = κio (1 + α i s ) (
κ s ( p ', s ) = κ so 1 + α sp ln ( p ' pref ) ) exp ( α s )
ss
For a3 = 0, the model (ITYCL=1) coincides with the elastic part of BBM for constant
coefficients:
∆e s + 0.1 −κi 0 −κ s 0 s + 0.1
=a1∆ ln ( − p ') + a2 ∆ ln = ∆ ln ( − p ') + ∆ ln
1+ e 0.1 1 + e 1+ e 0.1
For a3 different from zero the equation ( ITYCL=1) can be expanded in the following way.
∆e s + 0.1 s + 0.1
= a1∆ ln ( − p ') + a2 ∆ ln + a3 ∆ ln ( − p ') ln =
1+ e 0.1 0.1
s + 0.1 s + 0.1
= a1 + a3 ln ∆ ln ( − p ') + a2 + a3 ln ( − p ') ∆ ln =
0.1 0.1
a s + 0.1 a3 s + 0.1
= a1 1 + 3 ln ∆ ln ( − p ') + a2 1 + ln ( − p ') ∆ ln
a1 0.1 a2 0.1
Depending on the values of the parameters, negative compressibility can be obtained. This
can be limited with the Kmin indicated above.
For a3 and a4 different from zero the equation (ITYCL=5) can be transformed in the following
way.
∆e s + 0.1 s + 0.1
= a1∆ ln ( p ') + a2 ∆ ln + a3 ln ( p '/ pref ) ∆ ln + a4 s∆ ln ( p ') =
1+ e 0.1 0.1
s + 0.1
=[ a1 + a4 s ] ∆ ln ( p ') + a2 + a3 ln ( p ') ∆ ln
0.1
a a s + 0.1
= a1 1 + 4 s ∆ ln ( p ') + a2 1 + 3 ln ( p '/ pref ) ∆ ln
a1 a2 0.1
κ κ
a1 =
− i0 a2 =
− s0 a4 =
a1α i a3 =
a2 α sp
1+ e 1+ e
_______________________________________
130
VISCOELASTICITY FOR SALINE MATERIALS (Mechanical data 1)
CODES in ICL=2 ITYCL=1
ROOT_gen.dat
DESCRIPTION Parameters for linear viscous deformation model. The deformation
mechanism fluid assisted diffusional transfer (FADT) was applied
to develop an equation for creep of salt under wet conditions.
EQUATIONS Strain rate for a linear viscoelasticity is computed as:
FADT
dε 1 1
= d ( σ ′ - p ′I ) + v p′ I
dt 2 ηFADT 3 ηFADT
where σ ’ is the effective stress tensor (σ ’ = σ + Pf, where Pf =
max(Pg, Pl)), p’ is the mean effective stress (p' = p + Pf), I is the
identity tensor.
1 16 B(T ) S l
= d
g FADT ( e)
2η d
FADT d 3
0
where gFADTd(e) v
and gFADT (e) are internal nonlinear functions of
void ratio (e), and Sl is degree of saturation.
AB − QB
B(T ) = exp
RT RT
3 g 2 e 3/ 2 g2
v
g FADT ( e) = d
g FADT ( e) =
(1 + e) (1 + e)
1 2e
g= f =
(1 − f ) 2 3(1 − e )(1 + e)
3/ 2
PARAMETERS
P1 do m Grain size
P2 AB s-1MPa-1m3 Pre-exponential parameter
P3 QB J mol-1 Activation energy
P4 Void
P5 Void
If the pre-exponential parameter is set to zero (AB = 0.0) the viscous counterpart of the model
does not work. In this way the parameter acts as option because the value of these pre-
exponential parameter is checked to decide if this mechanism is considered.
This viscoelastic model (corresponding to FADT mechanism of deformation) requires that the
liquid pressures are computed or, alternatively, a value of PLCONS greather that -10-12 MPa.
Otherwise liquid is considered inexistent and the mechanism FADT remains inactive.
131
VISCOPLASTICITY FOR SALINE MATERIALS (Mechanical data 1)
CODES in ICL=3 ITYCL=1
ROOT_gen.dat
DESCRIPTION The deformation mechanism referred as dislocation creep (DC)
has been applied to develop an equation for creep of porous salt
aggregates. This mechanism leads to nonlinear dependences on
stresses.
EQUATIONS Strain rate for a nonlinear viscoelasticity is computed as:
DC
dε 1 ∂G
= d Φ (F)
dt ηDC ∂σ ′
where G is a flow rule, F is a stress function and Φ is a scalar
function. These functions are defined as:
2 1
−p ηv n+1
F =G = q 2 + Φ ( F ) =F n α p = dDC
α
p ηDC
1 1
= A(T ) g DC
v
(e) = A(T ) g DC
d
(e)
η DC
v
η d
DC
where n is the power that comes from the rock power law and
gDCd(e) and gDCv(e) are internal nonlinear functions of void ratio
(e) defined as follows:
v
g DC (=
e) 3( g − 1) n f
n −1
1+ g + g 2 2g +1 1
d
g DC =
(e)
f +
3 3 g
where f and g are functions of the void ratio (see expression in
the constitutive law Viscoelasticity for saline materials).
The temperature dependence is considered as:
− QA
A(T ) = A A exp
RT
PARAMETERS
P1 AA s-1 MPa-n Pre-exponential parameter
P2 QA J m-1 Activation energy
P3 n - Stress power
P4 – P10 Void
If the pre-exponential parameter is set to zero the viscous counterpart of the model does not
work. In this way the parameter acts as option because the value of this pre-exponential
parameter is checked to decide if this part of the model is operating.
132
VISCOPLASTICITY FOR GRANULAR MATERIALS (Mechanical data 1)
CODES in ICL=33 ITYCL=1,2 or 3, or 4,5, or 9,11
ROOT_gen.dat
DESCRIPTION Viscoplasticity (general model for soils and rocks)
dε ∂G
EQUATIONS Viscoplastic constitutive model: = Γ Φ( F )
dt ∂σ
and the stress function adopted is: Φ( F ) = F m
where the yield function and the flow rule are defined as:
ITYCL=1: G = F = q 2 − δ( p') n + δ( p') n +1 / po
ITYCL=2: G = F = q 2 − δpo ( p') n + δ( p') n +1
ITYCL=3 (Cam-Clay): G = F = q 2 − δ 2 ( po p'− p' 2 )
ITYCL=4, 5 (Drucker-Prager –based on Mohr-Coulomb
parameters):
G = F = q − δp '− cβ
6 sin φ '
δ
= M=
3 − sin φ '
6 cos φ '
β=
3 − sin φ '
For this model, equations are written assuming p>0 compression, but the
program uses the standard sign criteria for continuum mechanics).
The invers of viscosity can be written as a function of temperature as:
− Q
Γ = Γo exp
RT
The hardening laws are expressed in general as:
1+ e 1+ e
dpo= Dd ( εlv ) + po d εv= Dl εlv−1 + po d εv
λ−κ χ
where D ≠ 0 or (λ-κ) ≠ 0 permit to use each one of the two
possibilities.
d δ= Ed ( εlv )= El εlv−1d εv
where po and δ are parameters in the yield surface and flow rule.
Invariants used in the models are defined as:
p = σ oct =
1 1
(
I = σ + σ y + σz
3 1 3 x
)
3 1
( ) ( ) + (σ − σ x ) + 6( τ 2xy + τ 2yz + τ 2zx )
2 2 2
q= τ oct = σx − σ y + σ y − σz z
2 2
133
VISCOPLASTICITY FOR GRANULAR MATERIALS (Mechanical data 1)
α= a5 + a6Wd + a7 Wd − Wd 0 dWd= qd ε d
2
I1 = (σ x + σ y + σ z )
1 1
p=
3 3
1 1
( σ x − σ y ) + (σ y − σ z ) + (σ z − σ x ) + 6(τ xy2 + τ yz2 + τ zx2 )
2 2
J= =
2
2 Sij S ij
2 6
1 −1 3 3 J 3
θ
= sin (− 3
), −30 ≤ θ ≤ 30
3 2 J 2
134
1
1 3 + tan θT ⋅ tan 3θT + 3 sign (θ ) ⋅
A= ⋅ cos θT
3 ( tan 3θT − 3 tan θT ) sin ϕ
1 1
B=
− sign (θ ) sin θT + sin ϕ ⋅ cos θT
3 ⋅ cos 3θT 3
+1, θ ≥ 0
sign (θ ) =
−1, θ < 0
135
PARAMETERS FOR ITYCL=1 and 2 (Mech. data 1 ─viscoplasticity for granular materials)
Power of the stress function (integer value)
P1 m
(typical value = 3)
1/Viscosity (for plasticity use a sufficiently large
P2 Γ o =1/η s-1 MPa-m
value)
Activation energy (=0 for temperature
P3 Q J mol-1
independent model)
P4 D MPa Constant in hardening law
P5 po MPa Initial value of po (positive)
P6 l - Power of hardening laws
P7 χ - Parameter in hardening law
P8 n - Power in F and G
P9 E - Constant for hardening law
P10 δ - Initial value of δ
P4 b - Cohesion: c ' = ( a + bs ) g ( φ ) ; g ( φ ) = ( f + f ) / 2
P5 a MPa Cohesion: c ' = ( a + bs ) g ( φ ) ; g ( φ ) = ( f + f ) / 2
Porosity function: f ( φ ) = 1 − ( φ φoo )
n
P6 n
Parameter to reduce dilatancy (ranges between 0 and 1)
P7 α
(default = 1)
P8 φoo - Reference porosity
P9 Void
P10 δ - Equivalent to M
136
PARAMETERS FOR ITYCL=5 (Mech. data 1 ─viscoplasticity for granular materials)
P1 – P3 Same as in ITYCL = 1
P4 Void
P5 c MPa Cohesion
P6 Void
P7 α Parameter to reduce dilatancy (ranges between 0 and 1) (default = 1)
P8 – P9 Void
P10 δ - Equivalent to M
137
HISTORY VARIABLES:
The Visco-plasticity model (ICL =33) requires four history variables:
Hist_var 1 Po MPa Evolution of preconsolidation mean stress
Hist_var 2 δ0 Evolution of δ (parameter of the flow rule)
Hist_var 3 EDP Plastic deviatoric strain
Hist_var 4 EVP Plastic volumetric strain
The first two variables can be assigned as initial conditions on surfaces/volumes if an initial
particular distribution on the geometry is required. The procedure is the same as followed by
initial stresses as was described in chapter II. PREPROCESS, PROBLEM DATA, section
II.2.3.5.
If no value is assigned for the first two variables in conditions, internally, the program sets the
input parameters P5 (for Po*) and P10 (for δ0) of the ICL=33 (ITYCL=1, 2 or 3), as initial
values.
The evolution of the four history variables can be visualized as an output in Post-process GID
interface.
138
VISCOPLASTICITY - GENERAL
CODES in ICL=34, 35 and 36 ITYCL=1, 16
ROOT_gen.dat
DESCRIPTION Viscoplasticity (general model for unsaturated soils based on Desai and
Perzyna theory)
EQUATIONS ITYCL = 1: Viscoplastic constitutive model
N
dε ∂G F
=Γ φ( F ) φ( F ) =
dt ∂σ ' Fo
where the yield function is defined as: F ( J1 , J 2 D , J 3 D=
, s ) aJ 2 D − µ 2 Fb Fs
with the following additional functions:
γ − ( J10 ( s ) + k2 s + k4 ) ( J1 + k1s + k4 ) + ( J1 + k1s + k4 ) − k3 sJ10 ( s )
2− n
Fb =
n 2
27 ( −3/ 2 )
Fs= (1 − βs S ) S= J 3D ( J 2 D )
m
2
The viscoplastic potential is defined similarly as:
G ( J1 , J 2 D , J 3 D , =
s ) aJ 2 D − bµ 2 Fb Fs
where b is a non-associativity parameter.
Hardening is described with the following function, which is equivalent to
the BBM model:
λ (0) −κ
J10* λ ( s ) −κ
=J ( s ) 3=
p c 0
1
c
po ( s ) J10 ( s ) / 3
3p
λ ( s ) = λ ( 0 ) (1 − r ) exp ( −β s ) + r
Suction and net stress are defined as:
= (
s max ( Pg − Pl ) ,0 ) σ n =σtotal
n − max( Pg , Pl )
And the invariants are:
p '=
1
3
( σ 'x + σ ' y + σ 'z )= p − max( pg , pl )= J1 / 3 − max( pg , pl )
1 1
J2D = trace(s : s) = q 2 s = σ '− p ' I
2 3
Hardening depends on viscoplastic volumetric strains according to:
1+ e 1+ e dpo* 1 + e vp
dJ=
0*
J1o*d εvvp ⇔ dp
= *
po* d εvp
v ⇔ = d εv
1
(λ ( 0) − κ) o
( λ ( 0) − κ) po* χ ( 0 )
which is equivalent to BBM model as shown.
Note that, using k1=3k, k2=3k, k3=0, k4=0, and Fs=1:
1
p, s ) a q 2 − µ 2 γ32 −( po ( s ) + ks ) 2− n ( p + ks ) n + ( p + ks ) 2
F ( q, =
3
In the same way the viscoplastic potential is described as:
1
, s ) a q 2 − bµ 2 γ32 −( po ( s ) + ks ) 2− n ( p + ks ) n + ( p + ks ) 2
G (q, p=
3
which incorporates a parameter to allow for non-associativity conditions.
Strength can be considered also a function of suction:
s
µ
µ (s) =µ dry − ( µ dry − µ sat ) sat (µ < µ dry )
µ dry
sat
139
PARAMETERS FOR ITYCL=1
140
This model has the following case that can be used:
Cam-Clay model Cap models
a 3 3
n 1 3,5,7,9
γ -1/9 +1/9
k1 0 0
k2 -3k ?
k3 3k ?
k4 3ps0 ?
µ M M
Parameters k, ps0 (see BBM). M: slope of critical state line
141
VISCOPLASTICITY - GENERAL
CODES in ICL=34,35 and 36 (continued)
ROOT_gen.dat
DESCRIPTION Viscoplasticity (general model for unsaturated soils based on Desai and
Perzyna theory).
EQUATIONS ICL=36 (Viscoplasticity – General Parameters 3), ITYCL=2:
The LC curve is defined in the following way:
3 p y λ d ( s ) + J1o* ( λ i − κ ) λ(0) − κ o*
J1 ( s ) = =
3J y + ( J1 − 3 J y )
o
λ + λ (s) − κ
i d
λ( s) − κ
and the form of the compressibility is: λ ( s ) = λ ( 0 ) (1 − r ) exp ( −β s ) + r
Hardening depends on viscoplastic volumetric strains according to:
1 vp
=dJ10* 3 d εv
λ (0) − kio
ICL=36 (Vicoplasticity – General Parameters 3), ITYCL=3:
The LC curve is defined in the following way:
3 p y λ d ( s ) + J10* ( λ i − κ ) λi − κ 0*
J (s) = =
3 p y + i ( J1 − 3 p y ) and the
o
λ + λ (s) − κ λ + λ (s) − κ
1 i d d
s + 0.1
form of the compressibility is: λ ( s ) = λ i + λ d ( s ); λ d ( s ) = λ d − α s ln
0.1
Hardening depends on viscoplastic volumetric strains according to:
1 vp
=dJ10* 3 d εv
λ i − kio
For more details see Oldecop & Alonso, 2001.
142
ICL=36 (Mechanical data 1 Vicoplasticity – General Parameters 3) ITYCL=3
void -
P2 λi-κ - Viscoplastic compression parameter
P3 λd - Viscoplastic compression parameter
P4 αs - Parameter λ(s) curve
P5 py MPa Parameter in LC curve
P6 k1 -
P7 k2 -
P8 k3 -
P9 k4 -
143
VISCOPLASTICITY - GENERAL
CODES in ICL=34,35 and 36 (continued)
ROOT_gen.dat
DESCRIPTION Viscoplasticity – Zero thickness element
EQUATIONS ICL=34,35,36 ITYCL=16
Viscoplastic constitutive model:
N
dε ∂G F
=Γ φ( F ) φ( F ) =
dt ∂σ ' Fo
- The yield function is defined as: F = τ2 − ( c '− σ ' tan φ ') , where τ is the
2
shear stress; c ' is the effective cohesion, σ ' is the net normal stress, tan φ '
is the tangent or internal friction angle and χ is the tensile strength.
- Evolution of strength parameters with suction:
( c0 + c1ψ ) + ( b0 + b1ψ ) (1 − e−b tan α )
c0' ( ψ ,αa ) = 2 a
parameter that controls the shape of the cohesion- α a curve. tan φ0( ψ ,α ) is the
'
a
tangent of the internal friction effective initial angle that depends of ψ and
α a (Figure VIb.3); tan φ0 is the value of tan φ'0 for ψ =0 and α a = 0º ; t1 is
the slope of the tan φ'0 - ψ fit line for α a =0º (Figure VIb.4); d0 and d1 are
model parameters that control the increment of tan φ'0 with suction for
αa = 5º ,15º ,30º , 45º and tan α a is the geometric tangent of the asperity
roughness angle.
- The softening is defined as:
u vp u vp
=c ' c0' 1 − s* tan φ=' tan φ'0 − ( tan φ0' − tan φ'res ) s*
uc uφ
where usvp is the visco-plastic shear displacement; uc* is the critical value
for shear displacement when c’=0; tan φ'res is the tangent of internal friction
effective residual angle and uφ* is the critical value of shear displacement
when tan φ=' tan φ'res .
- The viscoplastic potential is defined as:
∂G
= 2 tan φ ' ( c '− σ ' tan φ ' ) f σdil f cdil , 2τ
T
∂σ
σ' σ' c'
f σdil = χ d tan α a 1 − exp −βd f cdil =
qu qu c0'
where, qu is the compression strength at which dilatancy vanishes; χ d and
βd are model parameters.
144
Figure VIb.1. Effective cohesion vs. αa (Zandarin, 2010)
Figure VIb.3. Effective tangent of the internal friction angle vs. α (Zandarin, 2010)
a
Figure VIb.4. Effective tangent of the internal friction angle vs.suction (Zandarin, 2010)
145
PARAMETERS FOR ITYCL=16
146
HISTORY VARIABLES:
The Viscoplasticity-general model (ICL =34,35,36) requires four history variables:
Hist_var 1 (J1o*)F MPa Evolution of size of F (note that this is the 1s invariant)
Hist_var 2 (J1o*)G MPa Evolution of size of G (note that this is the 1st invariant)
Hist_var 3 EDP Plastic deviatoric strain
Hist_var 4 EVP Plastic volumetric strain
The first two variables can be assigned as initial conditions on surfaces/volumes if an initial
particular distribution on the geometry is required. The procedure is the same followed by initial
stresses as was described in chapter II. PREPROCESS, PROBLEM DATA, section II.2.3.5.
If no value is assigned for the first two variables in conditions, internally, the program sets the
input parameters P7 (for (J1o*)F) and P8 (for (J1o*)G) of the ICL=35 (ITYCL=1), as initial
values.
The evolution of the four history variables can be visualized as an output in Post-process GID
interface.
Note: Effective stresses plotted in the Post-process GID interface correspond with net stresses
for unsaturated conditions and Terzaghi's effective stresses for saturated conditions. Stress and
strain invariants follow the soil mechanics notation (positive for compression).
147
APPENDIX 1. ZERO THICKNESS
MECHANICAL PROBLEM
The mechanical behaviour of the joint elements is defined by the relation between stress and
relative displacements of the joint element (Figure VIb.5) calculated on the mid-plane. The
mid-plane relative displacements are interpolated using the nodal displacements and the shape
functions.
σ ut
3 4
mp1 τ mp2 a0 un
a
1 2 dl
Figure VIb.5. Joint element with double nodes. a) Stress state on the mid-plane of the joint
element. b) Relative displacement defined at mid plane.
The normal and shear displacement increment calculated on the midplane is defined as:
un
=
w mp = r N ump [ −I 4 I4 ] u j
us mp
where un and us are the normal and tangential relative displacements, r is a rotation matrix, Nmpu
is a matrix of shape functions, I4 is an identity matrix of 4th order, uj is the vector of nodal
displacements.
The stress tensor on the mid-plane is calculated as a function of displacement components,
normal and shear:
σ '
σ 'mp
= = D w mp
τ
mp
where σ’mp is the net effective stress on the mid-plane of the element and it is defined as σ’mp
= σmp- max(Pg; Pl) mp; (σ is total mean stress; Pg is the gas pressure and Pl is the liquid
pressure, both interpolated to the mid-plane of the element); τ is the tangential stress on the
mid-plane; D is the stiffness matrix which relate relative displacements with the stress state.
148
ELASTIC BEHAVIOUR
The elastic behaviour of the joint is established as a relationship between the normal-tangencial
effective (σ’,τ) and the normal-tangential (un, us) relative displacement of the joint element.
This is established using a normal (Kn) and a tangential stiffness (Ks) coefficients. The normal
stiffness depends on the opening of the joint.
σ ' K n 0 un m
= Kn = , K s = constant
τ 0 K s us a − amin
where m is a parameter of the model; a is the opening or aperture of the element and amin is
the minimum opening or aperture of the element.
σ'
amin a
Figure VIb.6. Elastic constitutive law of the joint element. Normal stiffness depends on joint
opening.
VISCO-PLASTIC BEHAVIOUR
The constitutive behaviour for the mechanical of rough rock joint was developed based on the
formulations proposed by Gens, et al. (1985) and Carol, et al, (1997). According to these
theories, it is necessary to define a yield surface, a plastic potential and a softening law to
mathematical model the shear behaviour of a joint.
The visco-plastic displacements occur when the stress state of the joint reaches a failure
condition. The failure surface can be defined linearly (the one implemented in code_bright):
F ≡ τ 2 − ( c '− σ ' tan φ ')
2
where τ is the shear stress; c’ is the effective cohesion; σ’is the net normal stress and tanφ’ is
the tangent of effective angle of internal friction. χ is a parameter.
τ F>0
F=0
F<0
Φ c
σ
Figure VIb.7. Hyperbolic (continuous) and linear (dashed) failure surface and strength
parameters
149
SOFTENING LAW
The strain-softening behaviour of the joint subjected to shear stress is modelled considering the
degradation of the strength parameters tanφ’ and c’.
The degradation of the parameters tanφ’ and c’ is considered dependent on the accumulated
viscoplastic shear displacement. This is based on the slip weakening model introduced by
Palmer & Rice (1973). In this way the tangent of the friction angle decays from the peak (intact
material) to the residual value and the cohesion from the value c to zero. Two different values
u* permit to define the decrease of cohesion (u*c’) and friction angle (u*tanФ’). The
mathematical expressions are:
u vp
=
c ' c '0 1 − s*
uc
where c’ is the effective cohesion that corresponds to the visco-plastic shear displacement usvp;
c’0 is the initial value of the effective cohesion; u*c is a parameter.
u vp
tan φ ' =tan φ '0 − ( tan φ '0 − tan φ 'res ) s*
uφ
tanφ’ is the tangent of effective angle of internal friction that corresponds to the visco-plastic
shear displacement usvp; tanφ’0, tanφ’res and u*φ are parameters.
1
τ c tanΦ
2 1
c 1 tanΦ 0
Φ0
0
170°
Φres
σ 2
tanΦres
2
u* c u vp
s u * tanΦ u vp
s
Figure VIb.8 a) Evolution of the failure surface during softening. b) Softening law of cohesion.
c) Softening law of tanφ.
VISCO-PLASTIC DISPLACEMENTS
A viscoplastic yield surface implies that when F < 0 the stress state of the element falls inside
of the elastic region. In contrast, if F >= 0 the displacements of the element undergo a visco-
plastic component. The viscoplastic displacements are calculated as:
vp
dw F ∂G
= Γ <φ >
dt F0 ∂σ
where Γ is a fluidity parameter. The visco-plastic displacement rate is given by a power of
law considered for the function φ:
∂G ∂G
∆unvp =
ΓFN ∆t ∆usvp =
ΓFN ∆t
∂σ ∂τ
150
PLASTIC POTENTIAL SURFACE AND DILATANCY
To calculate the direction of displacements it is necessary to define the derivatives of G with
respect to stresses:
∂G
= 2 tan φ ' ( c '− σ ' tan φ ') fσdil f cdil
T
, 2τ
∂σ
This is a non-associated flow rule, because of the inclusion of both parameters fσdil and fcdil
which consider the dilatant behaviour of the joint with shear stresses (Lopez, et al. 1999). The
amount of dilatancy depends on the level of the normal stress and on the degradation of the
joint surface.
The following expressions describe these effects:
∂∆unvp
ΓN F N −1 2τ fσdil f cdil 2 tan φ ' ( c '− σ ' tan φ ') ∆t
=
∂τ
∂∆usvp ∂∆usvp
ΓN F N −1 2τ 2 tan φ ' ( c '− σ ' tan φ ') ∆t
= = ΓN F N −1 2τ ∆t + Γ F N 2∆t
∂σ ∂τ
Finally, the elasto-viscoplastic mechanical model of the joint is expressed by the tangent
stiffness matrix:
−1
D=
evp
Ce + Cvp
More information about this joint element can be found in Zandarin (2010) –see References.
151
APPENDIX 2. EXAMPLE OF USE
The use of the BBM model can be achieved by combination of the non-linear elasticity and the
viscoplasticity for unsaturated soils.
εtotal = ε elastic + ε viscoplastic
The elastic part can be linear or nonlinear and may depend on suction and temperature. The
viscoplastic part can also be a function of suction.
The following parameters are required and example values are given:
The parameter P7 limits the stiffness coefficient, so it cannot go beyond a certain value. A value
in the range of 1 MPa to 20 MPa is normally required. The parameter P6 limits the possibility
of tensions. The model may use Poisson (P5) or shear stiffness (P4) alternatively. The use of
one of them implies the other is variable as they are related.
The so-called state surface is an old model for unsaturated soils is based on reversibility. The
volumetric strain is calculated in a reversible way according to:
∆e s + 0.1 s + 0.1
= a1∆ ln ( p ') + a2 ∆ ln + a3 ∆ ln ( p ') ln
1+ e 0.1 0.1
where p' is mean effective stress (mean stress plus maximum of liquid and gas pressure) and s
is suction (gas pressure minus liquid pressure). Shear strain is linearly elastic with modulus G
or, alternatively, a constant value of the Poisson’s ratio can be used.
For a3 = 0, the model coincides exactly with the elastic part of BBM for constant coefficients:
∆e s + 0.1 −κi 0 −κ s 0 s + 0.1
a1∆ ln ( p ') + a2 ∆ ln
= = ∆ ln ( p ') + ∆ ln
1+ e 0.1 1 + e 1+ e 0.1
For a3 different from zero the equation can be expanded in the following way:
152
∆e s + 0.1 s + 0.1
= a1∆ ln ( p ') + a2 ∆ ln + a3 ∆ ln ( p ') ln =
1+ e 0.1 0.1
s + 0.1 s + 0.1
= a1 + a3 ln ∆ ln ( p ') + a2 + a3 ln ( p ') ∆ ln =
0.1 0.1
a s + 0.1 a3 s + 0.1
= a1 1 + 3 ln ∆ ln ( p ') + a2 1 + ln ( p ') ∆ ln
a1 0.1 a2 0.1
Depending on the values of the parameters, negative compressibility can be obtained. This can
be limited with the Kmin indicated above.
153
VISCOPLASTICITY - GENERAL
154
CODE_BRIGHT. CONSTITUTIVE LAWS
This chapter contains the different models available and the corresponding parameters required
by each model. The following constitutive laws are available:
HYDRAULIC AND MECHANICAL CONSTITUTIVE MODELS
THERMAL CONSTITUTIVE
MODELS (a) ELASTICITY (b)
RETENTION CURVE NONLINEAR ELASTICITY (b)
INTRINSIC PERMEABILITY VISCOPLASTICITY FOR SALINE MATERIALS
LIQUID PHASE RELATIVE (b)
PERMEABILITY VISCOPLASTICITY FOR GRANULAR
GAS PHASE RELATIVE MATERIALS (b)
PERMEABILITY VISCOPLASTICITY - GENERAL (b)
DIFFUSIVE FLUXES OF MASS DAMAGE-ELASTOPLASTIC MODEL FOR
DISPERSIVE FLUXES OF ARGILLACEOUS ROCKS (c)
MASS AND ENERGY THERMO-ELASTOPLASTIC MODEL FOR SOILS
CONDUCTIVE FLUX OF HEAT (d)
BARCELONA EXPANSIVE MODEL FOR SOILS
(e)
CASM’s FAMILY MODELS (f)
PHASE PROPERTIES (a) EXCAVATION PROCESS (g)
SOLID PHASE PROPERTIES
LIQUID PHASE PROPERTIES
GAS PHASE PROPERTIES
155
DAMAGE-ELASTOPLASTIC MODEL FOR ARGILLACEOUS ROCKS
CODES ICL=70 to 79 ITYCL=see below
EQUATIONS For this model, equations are written assuming Soils Mechanics
convention (p > 0, εv > 0, compression). p is the mean effective stress,
J the square root of the second invariant of deviatoric stress tensor, θ the
Lode’ angle (-30º in triaxial compression, +30º in triaxial extension)
(see ICL 21 to 27 for their definition).
Elastic law:
eM ds p
dσ ijM= Dijkl d ε kl − δ kl M − d ε kl
M
Ks
eM
σ ijM are the stresses prevailing at clay particles contact, Dijkl is the
mechanical elastic stiffness matrix of the clay, d ε klM are the strains
corresponding to the clay matrix deformation (equal to the external
strains), K sM is the bulk modulus against suction changes (if any), d ε klp
are the plastic strains of the clay matrix.
156
Value of EMh is computed as:
EMh = max(EMh0, EMh1 pM + EMh2)
vhere EMh0, EMh1 and EMh2 are parameters of the model
K M
=
(1 + e )σ hM
κM
h
Kv =
M (1 + e )σ vM
κM
K sM =
(1 + e )( s + patm )
κ sM
3
is clay matrix tensile strength, c’M clay matrix cohesion,
ptM = c 'M cot φ 'M
φ’M clay matrix friction angle.
Cohesion depends on suction following the law:
'M ( s ) c 'M ( 0 ) + s tan φbM
c=
157
Hoek & Brown (1980) (ITYCL=2):
π
4sin 2 θ M −
6 M 2 2 m sin θ
M M
=F p J − J M − m M ( p M + ptM ) ≥ 0
M
Rc 3
is clay matrix tensile strength, RcM clay matrix uniaxial
ptM = RcM / m M
compressive strength, mM a parameter defining the shape of the
parabolic yield criterion.
where ptM is the clay matrix tensile strength, p0M the clay matrix isotropic
yield locus and M the slope of the critical state line in the pM –vqM
diagram. The following dependencies on suction are considered:
ptM = k M s
λ M ( 0 ) −κ M
p*M λ M ( s )−κ M
= ( s ) λ M ( 0 ) (1 − r M ) + re− β s
with λ M=
M
M M
p p 0M
0 c
pc
Mohr Coulomb and Hoek & Brown yield criteria present corners in the
deviatoric plane. They are smoothed using Sloan & Booker (1986)
procedure. Lode’s angle θt at which smoothing starts must be defined
(see ICL = 74).
158
Mohr-Coulomb (ITYCL=1):
1
cos θ +
Gp = sin θ M sin φ 'M J M − ω M sin φ '( p M + ptM )
M
3
ptM, c’M and φ’M are parameters defining the yield criterion. ωΜ is a
parameter defining the non associativity of the flow. It takes a value
equal to 1 when associated and equal to 0 for null dilatancy.
Hardening law:
Mohr-Coulomb (ITYCL=1):
A softening law is introduced through the following dependency of the
tensile strength on the plastic strain:
ptM = c 'M cotanφ 'M
with
if ε eq < ξ rc 0
pM M
c 'M =c'M
peak
c'Mpeak (1 − α M )
(ε − ξ rcM0 ) if ξ rc 0 ≤ ε eq ≤ ξ rc
M pM M
c 'M =c'Mpeak + pM
(ξ M
rc −ξ M
rc 0 ) eq
if ξ rc < ε eq
M pM
c 'M =α M c'M
peak
and
if ε eq < ξ rφ 0
pM M
φ 'M =φ 'M
peak
φ 'Mpeak (1 − β M )
(ε − ξ rMφ 0 ) if ξ rφ 0 ≤ ε eq ≤ ξ rφ
M pM M
φ 'M =φ 'Mpeak + pM
(ξ M
rφ −ξ M
rφ 0 ) eq
if ξ rφ < ε eq
M pM
φ 'M =β M φ 'M
peak
159
αΜ = c’resM / c’peakM is a brittleness parameter for the cohesive
component, βΜ = φ’resM / φ’peakM is a brittleness parameter for the
frictional component, εeqpM is the equivalent plastic strain, ξrc0M the
accumulated equivalent plastic strain at which the cohesion starts to
degrade, ξrcM the accumulated equivalent plastic strain at which the
residual cohesion cres’M is reached, ξrφ0M the accumulated equivalent
plastic strain at which the friction angle starts to degrade, ξrφM the
accumulated equivalent plastic strain at which the residual friction angle
φ’res’M is reached. αΜ = βΜ = 1 means perfect plasticity, αΜ = βΜ = 0,
total strength degradation (residual cohesion and friction angle equal to
0).
Bond behaviour
Elastic law:
dσ ijb Dijkl
= eb
( dε klb − dε kld )
eb
Dijkl is the secant damaged elastic matrix. It is related to the secant
−L
undamaged elastic tensor Dijkl by Dijkl = e Dijkl . L is the damage
eb 0 eb eb 0
variable, related to the ratio of bond mickocraks area over the whole
eb 0
bond area. Dijkl is defined by the undamaged bond Young’s modulus Eb
and bond Poisson’s ratio νb through the classical linear orthotropic
elasticity.
160
Damage locus: Damage locus is defined as an energy threshold
1 b b
=
Fd σ ij ε ij − r b ( s )
2
rb is the value of energy threshold that depends on suction following:
r b ( s=
) r b + r0bs s
is a parameter which controls the change of bond damage locus with
r0bs
suction.
Rate dependency: Rate dependency is introduced as a delayed
microcracking and use the visco-damage formalism. Damage variable
is expressed as a function of the distance between the current bond stress
point and the infinitely slow damage locus:
dt
dL = Fd
ηb
where dt is the time increment, ηb is the bond viscosity and 〈〉 are the
Macauley brackets. Infinitly low damage locus takes the form:
ηb
Fd =
Fd − dL ≤ 0
dt
Damage rule: Damage rule gives the evolution of damage strain d ε kld
with damage variable L. This relation is constrained by bond elastic
moduli evolution and must take the form: d ε kld = ε klb dL
Damage evolution law: It gives the evolution of damage locus rb with
damage variable L. Three different expressions may be considered:
a) linear: r= r0 + r1 L
b b b
b
a) exponential: r b = r0b exp r1 L
a) logarithm: r= r0 + r1 ln L
b b b
r0 is the damage of the intact material and r1 a parameter giving the rate
of evolution (higher value of r1 gives lower damage rate). r1 is taken
function of suction following: r1= b
r10b + r1bs s
r10b is a parameter which controls the damage evolution rate for the
saturated bond material and r1bs is a parameter which controls the change
of damage evolution rate with suction.
Coupling behaviour: Coupling comes from the restrictions that local
strain ε ij and ε ij must be compatible with the external strain ε ij and
M b
161
Since this model requires a substantial number of parameters, several ICL's are included in
Mechanical Data 3:
Argillite Bonding (ICL=70) contains elastic parameters for the bonding
Argillite Matrix (ICL=71) contains elastic parameters for the soil matrix
Argillite – Coupling (ICL=72) contains parameters for the coupling between bonding
and soil matrix
Argillite – Yield vol (ICL=73) contains parameters for the shape of the yield function
in the p-q plane (soil matrix)
Argillite – Yield dev (ICL=74) contains parameters for the shape of the yield function
in the deviatoric plane (soil matrix)
Argillite – Plastic vol (ICL=75) contains parameters for the shape of the plastic
potential in the p-q plane (soil matrix)
Argillite – Plastic dev (ICL=76) contains parameters for the shape of the plastic
potential in the deviatoric plane (soil matrix)
Argillite – Hardening (ICL=77) contains parameters for the hardening law (soil matrix)
Argillite – Visco (ICL=78) contains parameters for the viscoplastic model (bonding and
soil matrix)
Argillite – Control parameters (ICL=79) contains parameters to control the
integration of the constitutive law
162
PARAMETERS FOR ARGILLITE BONDING ICL=70 (Damage orthotropic elastic model for
the bond material). ITYCL=1: Linear damage evolution law. ITYCL=2: exponential damage
evolution law. ITYCL=3: logarithm damage evolution law.
Young’s modulus for the bond material in the plane
P1 Ebh MPa orthogonal to the direction of orthotropy
If P6 = 0 or 1 (no anisotropy), isotropic Young’s modulus
Poisson’s ratio for the bond material in the plane orthogonal to
P2 νbh - the direction of orthotropy
If P6 = 0 or 1 (no anisotropy), isotropic Poisson’s ratio
P3 rb10 MPa Damage evolution rate for the saturated bond material
b
P4 r 0s MPa Change of bond damage locus with suction
b
P5 r 1s MPa Change of bond damage evolution rate with suction
Ratio of anisotropy (ratio between the value of Young’s
modulus for bond material in the direction perpendicular and
Ebh / Ebv =
P6 - parallel to the orthotropy axis)
νbhv / νbvh
If P6 = 0, a default value equal to 1 (no anisotropy) is assigned
to this parameters
P7 νbhv - Cross Poisson’s ratio of the bond material (νbhv = dεbh / dεbv)
Shear modulus of the bond material along the axis of
P8 Gv MPa
orthotropy (ratio dεh / dεv)
Bulk modulus against suction changes for the bond material
P9 Kb s MPa
(considered isotropic)
163
PARAMETERS FOR ARGILLITE MATRIX ICL=71, ITYCL = 1 or 2 (Linear orthotropic
elastic model for the soil matrix)
Young’s modulus for the saturated matrix in the plane
M
P1 E h* MPa orthogonal to the direction of orthotropy
If P6 = 0 or 1 (no anisotropy), isotropic Young’s modulus
Poisson’s ratio in the plane orthogonal to the direction of
P2 νMh - orthotropy
If P6 = 0 or 1 (no anisotropy), isotropic Poisson’s ratio
P3 K Ms MPa Bulk modulus against suction changes (considered isotropic)
Coefficient setting the change of E with suction according to
P4 EMs -
the equation EMh = EMh* + EMs s
P5 - - Void
Ratio of anisotropy (ratio between the value of Young’s
modulus in the direction perpendicular and parallel to the
EMh / EMv =
P6 - orthotropy axis)
νMhv / νΜvh
If P6 = 0, a default value equal to 1 (no anisotropy) is
assigned to this parameters
P7 νMhv - Cross Poisson’s ratio (νMhv = dεMh / dεMv)
P8 GMvh = GMhv MPa Cross shear modulus (GMvh = dτvh / dγvh)
Coefficient giving the change of the Young’s modulus with
P9 EM1h - the stress in the plane orthogonal to the direction of
orthotropy
Value of Young modulus in the plane orthogonal to the
P10 EM2h MPa
direction of orthotropy at null mean (net or effective) stress.
Elastic thermal deformation will be taken into account in the model in a coupled way if the
elasticity law ICL=5 ITYCL=1 is used. Parameter P3 corresponds to the thermal expansion
coefficient.
164
PARAMETERS FOR ARGILLITE MATRIX ICL=71 (Elastic model for the soil matrix),
ITYCL = 3 (Orthotropic Camclay type elastic model for the soil matrix)
165
PARAMETERS FOR ARGILLITE – YIELD VOL ICL=73, ITYCL=1 (Mohr Coulomb
criterion for the soil matrix – shape in pM-qM diagram).
P1 φ’M º Friction angle
P2 c’ M MPa Cohesion
Coefficient setting the change in cohesion with suction
P3 φ bM º
following: c’M = s tan(φbM)
P4 β M MPa-1 Coefficient setting the change in friction angle with suction
PARAMETERS FOR ARGILLITE – YIELD VOL ICL=73, ITYCL=2 (Hoek & Brown
criterion for the soil matrix – shape in pM-qM diagram).
Ratio of uniaxial compressive strength divided by tensile
P1 mM -
strength
Uniaxial compressive strength at the reference temperature
P2 Rc M MPa
(P8 field)
P3 rM - Coefficient setting the change in cohesion with suction
P4 βM MPa-1 Coefficient setting the change in cohesion with suction
P5 - - void
P6 - - void
Coefficient setting the decrease of uniaxial compressive
P7 kTM
strength with temperature
P8 T0 Reference temperature
166
PARAMETERS FOR ARGILLITE – PLASTIC VOL ICL=75, ITYCL=1,2,3 (Mohr Coulomb,
Hoek & Brown or Basic Barcelona model plastic potential for the soil matrix – shape in pM-
qM diagram).
Coefficient of non associativity (0: no volumetric plastic strain,
P1 ωΜ -
1: full volumetric plastic strain –associative plasticity)
P4 βΜ - =
φ 'Mpeak
167
PARAMETERS FOR ARGILLITE – HARDENING ICL=77, ITYCL=4 (Camclay model plus
Gallipoli et al. (2003) proposal).
P1 λ M(0) - Slope of the virgin loading line in the e-ln(pM) diagram
P2 a - Parameter for the expression of the state parameter e/eSS
P3 b - Parameter for the expression of the state parameter e/eSS
P4 R - Mean pore radius
168
Important notes:
1. Damage only, elastoplasticity only or coupled damage-elastoplasticity can be defined
depending on the combinations of ITYCL used. They are:
Damage
ITYCL = 1 ITYCL = 0 ITYCL = 0 ITYCL = 0
only
HISTORY VARIABLES
Argillite model requires in total 35 history variables. They are listed in the following table.
Variable 1 is the history variable of the elastoplastic model for the matrix. In the most general
case, value for this variable need to be specified in the input file (see format in the Part III).
In the case of Mohr Coulomb model (ITYCL(73)=1), variable 1 corresponds to the tensile
strength of the material pt. If a null value is specified in the input file for this variable, the
program will compute the value of the tensile strength for each element from the cohesion and
the friction angle of the material assigned to the element. The expression used to compute the
default value is: pt = c’ cotanφ’.
In the case of Hoek & Brown model (ITYCL(73)=2), variable 1 corresponds to the tensile
strength of the material pt. If a null value is specified in the input file from this variable, the
program will compute the value of the tensile strength for each element from the uniaxial
compressive strength Rc and the parameter m of the material assigned to the element. The
expression used to compute the default value is: pt = Rc / m.
In the case of Cam clay model (ITYCL(73)=3), variable 1 corresponds to the preconsolidation
pressure of the material p0*. IF A NULL VALUE IS SPECIFIED IN THE INPUT FILE FOR
THIS VALUE, THE PROGRAM WILL ABORT DUE TO A MATH ERROR (this ellipse of
Cam clay model degenerates into a point).
169
Variable 2 is the history variable of the damage model for the bond. In the most general case,
value for this variable need to be specified in the input file (see format in the Part III). However,
in many cases, this value can be set to 0. It means that damage will start for any stress changes.
Variable 3 is used internally to keep memory of the type of constitutive matrix (elastic, tangent
elastoplastic, tangent damage or tangent damage elastoplastic) used to build the global tangent
stiffness. No input from the user is associated to this variable.
Variable 4 is the plastic multiplier of the matrix. This variable measures the amount of plastic
strain. It is usually set to 0 at the beginning of the computation and updated within the argillite
subroutine. This variable is used only by the user for output visualization purposes (it provides
the spatial distribution of plastic strain intensity within the mesh at any output time).
Variable 5 is the damage multiplier of the bonds. This variable measures the amount of
damage. It is usually set to 0 at the beginning of the computation and updated within the argillite
subroutine. This variable is used by the user for output visualization purposes only (it provides
the spatial distribution of damage intensity within the mesh at any output time).
Variables 6 to 11 are the current stresses within the bonds. Variables 12 to 17 are the stresses
within the matrix at time of bond formation. In the most general case, all these variables have
to be specified in the input file*. However, from a practical point of view, most of the problems
can be run setting these values to 0 at the beginning of the computation (in this case, the damage
locus is centered on the initial stress state at the beginning of the computation).
Variables 18 to 35 contain the total, plastic and damage strain tensor components. They are
usually set to 0 at the beginning of the computation and updated within the argillite subroutine.
These variables are used by the user for output visualization purposes only (they provides the
spatial distribution of the total, plastic and damage strain within the mesh at any output time).
A summary of the different options for history variables input is summarized in the following
table.
170
Type of
input of Characteristics of the problem
history
variables
Mohr Coulomb Hoek & Brown Cam clay model
ITYCL(73)=1 ITYCL(73)=2 ITYCL(73)=3
Tensile strength is Tensile strength is NOT POSSIBLE
constant within each constant within each
material and material and computed
All variables
computed from from uniaxial
=0
cohesion and friction compressive strength and
angle parameter m
Damage locus is centred on the initial stress state of the problem. Material
is not damaged at time 0 of the computation.
From t=0, damage occurs for any stress changes
Mohr Coulomb Hoek & Brown Cam clay model
ITYCL(73)=1 ITYCL(73)=2 ITYCL(73)=3
Variable 1 ≠ Tensile strength can Tensile strength can vary Preconsolidation
0 vary from element to from element to element pressure can vary from
Variable 2 = element and is equal and is equal to the element to element and
0 to the specified value specified value for is equal to the specified
Variable 6 to for Variable 1 Variable 1 value for Variable 1
17 = 0 Damage locus is centred on the initial stress state of the problem. Material
is not damaged at time 0 of the computation.
Damage occurs for any stress changes
Mohr Coulomb Hoek & Brown Cam clay model
ITYCL(73)=1 ITYCL(73)=2 ITYCL(73)=3
Tensile strength can Tensile strength can vary Preconsolidation
Variable 1 ≠ vary from element to from element to element pressure can vary from
0 element and is equal and is equal to the element to element and
Variable 2 ≠ to the specified value specified value for is equal to the specified
0 for Variable 1 Variable 1 value for Variable 1
Variable 6 to Damage locus is centred on the stress state at time of bond stress state of
17 ≠ 0 the problem. It reflects the history of damage previous to time 0 of the
computations.
Damage occurs only when the energy input to the material is equal to the
specified value for Variable 2
Variables 3, 4 and 18 to 35 can be specified in the input file when the history of plastic and
damage strain previous to the current computations is wanted to be reproduced (particularly in
case of restart).
171
CODE_BRIGHT. CONSTITUTIVE LAWS
This chapter contains the different models available and the corresponding parameters required
by each model. The following constitutive laws are available:
22 THERMOELASTOPLASTIC MODEL FOR UNSATURATED SOILS 1 Contain parameters for the thermal terms
172
THERMO-ELASTOPLASTIC MODEL FOR SOILS
EQUATIONS For this model, equations are written assuming Soils Mechancis
compression (p > 0, εv > 0, compression).
The mechanical constitutive equation takes the incremental general form:
dσ ' = Ddε + hds
∂G
This equation is derived from: d=
ε d ε e + d ε=
p
( D e ) −1 d σ '+ αΙds + Λ
∂σ '
where an elasto-plastic constitutive law has been selected that is based on
a generalized yield surface that depends not only on stresses but on suction
as well: F = F ( σ ', εvp , s )
[
λ ( s) = λ ( o) (1 − r ) exp ( − β s) + r ]
ps= ps 0 + k s exp( − ρ∆T ) , ∆T = T − Tref
Hardening depends on plastic volumetric strain according to:
1+ e
= dpo* po*d εvp
λ ( 0 ) − kio
173
EQUATIONS The plastic potential is taken as:
(continuation) 3J 2
G= α 2 − L2p ( p '+ ps )( po − p ' )
gp
where gp is a function of the Lode angle and
Lp = M / g p θ= − π / 6
α is a non-asociativity parameter.
k i ( s) dp' k s ( p' , s) ds
dε ev = + + (α o + 2α 2 ∆T )dT
1 + e p' 1 + e s + 01 .
where:
s + 0.1
ki ( s ) = kio 1 + α i s + α il ln
0.1
(
', s ) k so (1 + α sp ln p ' pref ) exp α s s s
ks ( p = )
For deviatoric elastic strains, a constant Poisson’s ratio is used.
0.1 MPa (default = 0) and γ 0.7 is the shear strain at which the modulus G
is 0.722·G0 (approximately 0.7·G0).
ref
In order to “activate” this model, then G0 (ICL=22; ITYCL=1; P7) and
γ 0.7 (ICL=22; ITYCL=1; P8) must be greater than 0.
Since this model requires a substantial number of parameters, several ICL's are included:
TEP – Elastic parameters (ICL=21) contains elastic parameters (ITYCL = 1).
TEP – Thermal and other parameters (ICL=22) contains other parameters (ITYCL = 1).
TEP – Plastic parameters 1 (ICL=23) contains plastic parameters (ITYCL = 1).
TEP – Plastic parameters 2 (ICL=24) contains parameters for different aspects (ITYCL = 1).
TEP – Parameters Shape Yield Surf. (ICL=25) contains parameters for the function gy
(ITYCL = 1,2,3).
TEP – Parameters Shape Plastic Pot. (ICL=26) contains parameters for the function gp
(ITYCL = 1,2,3).
174
TEP – Integration Control Parameters (ICL=27) contains parameters for the integration of the
model (ITYCL = 1).
PARAMETERS FOR ICL=21 (TEP Elastic Parameters), ITYCL=1
P1 κ io - Initial (zero suction) elastic slope for specific volume-mean stress.
P2 κ so - Initial (zero suction) elastic slope for specific volume-suction.
P3 Kmin MPa Minimum bulk module.
P4 - -
P5 ν - Poisson’s ratio (-1 < ν < 0.5).
P6 α ss - Parameter for κs (only for expansive material).
P7 α il - Parameter for κi (only for expansive material).
P8 αi - Parameter for κi (only for expansive material).
P9 α sp - Parameter for κs (only for expansive material).
P10 pref MPa Reference mean stress (only for expansive material).
175
P5 k Parameter that takes into account increase of tensile strength due to suction
P6 ps0 MPa Tensile strength in saturated conditions
176
PARAMETERS FOR ICL=27 (TEP Integration Control Parameters), ITYCL=1
P1 Tole1 Yield surface tolerance (typically 1.e-8).
P2 Tole2 Elastic integration tolerance (typically between 1.e-4 and 1.e-6).
P3 Tole3 Plastic integration tolerance (typically between 1.e-4 and 1.e-2).
P4 µ Integration weight (ranges from 0 to 1) (typically 1).
P5 Index -1 elastoplastic matrix (typical value).
+1 elastic matrix.
P6 Itermaxc Maximum allowed subincrementations. When this value is reached, the
execution continues even if stress errors are large. Default=100.
Note that if Itermaxc > Itermaxs this condition is not used.
P7 Itermaxs Maximum allowed subincrementations. When this value is reached, the
stress and stiffness matrices calculation stops and the time step is reduced. A
message about the type of problem encountered appears. Default=10.
HISTORY VARIABLES:
The Thermo-elastoplastic (ICL=21 to 27) model requires two history variables:
Hist_var 1 Po* MPa Evolution of preconsolidation mean stress for saturated soil
Hist_var 2 eo Evolution of void ratio
Note: Effective stresses plotted in the Post-process GID interface correspond with net stresses
for unsaturated conditions and Terzaghi's effective stresses for saturated conditions. Stress and
strain invariants follow the soil mechanics notation (positive for compression).
______
177
CODE_BRIGHT. CONSTITUTIVE LAWS
This chapter contains the different models available and the corresponding parameters required
by each model. The following constitutive laws are available:
178
BARCELONA EXPANSIVE MODEL FOR SOILS (BExM)
DESCRIPTION Elastoplastic constitutive law for expansive soils (BExM by Alonso et al.,
1999).
EQUATIONS For this model, equations are written assuming Soils Mechanics convention
(p > 0, εv > 0, compression). p is the mean effective stress, J the square root
of the second invariant of deviatoric stress tensor, θ the Lode’ angle (-30º
in triaxial compression, +30º in triaxial extension) (see ICL 21 to 27 for
their definition).
Two levels of soil structure are defined: macrostructural level (Macro) and
microstructural level (micro) (see Appendix at the end of section VIe):
VpM VpM
Vp
VT Vpm
Vag
Vs Vs
where:
VT = Total volume
V pM = Volume of macropores
V pm = Volume of micropores
=V p Total volume of pores (=V pM + V pm )
Vs = Volume of solid particles
=Vag Volume of aggregates (=Vs + V pm )
179
m
Micro void ratio: e m = V p
Vs
= (
nm n m 1 − nM ) (
n =n M + n m 1 − n M )
The following additive equation for the total strain rate holds:
d ε kl d ε klMacro + d ε klmicro
=
where
d ε kl is the total strain rate: tr ( ε kl ) = de ;
1+ e
dε Macro
kl
is the macrostructural strain rate, associated to the
de M
macroskeleton: tr ( ε klMacro ) = ;
1 + eM
d ε klmicro is the microstructural strain rate (only volumetric behaviour is
dem
considered at this level): tr ( ε klMicro ) =
1 + em
Elastic law:
Macro ds Macro
=dσ ij Dijkl d ε kl
Macro
− δ kl Macro
− d ε klLC − d ε klSD − d ε klSI
3K s
micro
micro
ds
=dσ ij Dijkl d ε kl − δ kl
micro
micro
3K s
Macro
Dijkl is the mechanical elastic stiffness matrix which relates stress and
elastic deformation at macro level ( d ε klMacro );
micro
Dijkl is the mechanical elastic stiffness which relates stress and elastic
deformation at micro level ( d ε klmicro ).
K sMacro is the bulk modulus against macro suction changes and K smicro the bulk
modulus against micro suction changes.
d ε klLC are the macro plastic strains if LC is activated.
d ε klSD and d ε klSI are the macro plastic strains if SD and SI are activated,
respectively.
180
micro
Dijkl is defined by the bulk modulus Kmicro.
Macro bulk modulus against net stress changes is linearly dependent of the
logarithmic mean net stress (p) following the relationship:
K Macro =
(1 + e ) pM
κ Macro
Micro bulk modulus against effective stress changes is linearly dependent
of the logarithmic effective net stress (p+s) following the relationship:
micro
K
(1 + em ) ( p + s )
=
κ micro
where κ Macro and κ micro the slopes of the unloading/reloading lines in the plan
eMacro-lnp and emicro-ln(p+s), respectively. Note that the behaviour of the
microstructure is formalized by means of an effective stress concept
generalized for unsaturated conditions (effective stress is recovered).
Bulk modulus against suction changes is linearly dependent of the
logarithmic suction following the relationship:
K sMacro =
(1 + e ) ( s + p )
M
atm
κs
where κ s is the slope of the drying/wetting line in the plane eM-ln(s+patm).
patm is the atmospheric pressure and is taken equal to 0.1 MPa by default.
Yield function
M2
F LC = J 2 −
3
( p + p )( p t 0 − p) ≤ 0
where pt is the clay tensile strength, p0 the clay matrix isotropic yield
locus and M the slope of the critical state line in the p –q diagram. The
following dependencies on suction are considered:
pt = k s s
λ ( 0 ) −κ Macro
p* λ ( s )−κ Macro
p0 = pc 0 with λ ( s )= λ ( 0 ) (1 − r ) e − β s + r
pc
F SD =γSD − p − s F SI = p + s − γ SI
Rate dependency
Rate dependency is introduced as a visco-plastic mechanism. Plastic
multiplier λp is expressed as a function of the distance between the current
stress point and the inviscid plastic locus:
dt
dλ p = F
η
where dt is the time increment, η is the clay viscosity and 〈〉 are the
Macauley brackets. Inviscid plastic locus takes the form:
η
F=
F− dλ p ≤ 0
dt
181
where F can be either the LC, SD or SI yield criterion
Plastic potential
M2
J2 ω
G p =−
3
( p + p )( p − p )
t 0
Hardening law:
The hardening/softening law is introduced through the following
dependency of the saturated isotropic yield locus on the plastic strain:
dp0* (1 + e)
=
p0* λ − κ
( dε volLC + dε volSD + dε volSI )
K micro SD K micro SI
d=
γ SD d εvol + d εvol
f SD f SI
K micro SD K micro SI
d=
γ SI d εvol + d εvol
f SD f SI
SD
where f and f SI are the micro-macro interaction functions defined as
follow:
nSD nSI
p and p
f SD =f SD 0 + f SD1 1 − f=
SI
f SI 0 + f SI 1
p0 p0
182
PARAMETERS FOR BEXM – ELASTIC ICL=81, ITYCL=1 (Elastic model)
Matrix elastic stiffness parameter at macro level for changes in
P1 κ Macro - mean stress ( p )
Matrix elastic stiffness parameter at micro level for changes in
P2 κ micro - mean effective stress ( p + s micro )
P3 void -
P4 void -
Elastic macro stiffness parameter for changes in macro suction
P5 κs -
( s Macro )
P6 νM - Poisson’s ratio
Macro
P7 K min MPa Minimum bulk modulus at macro level
micro
P8 K min MPa Minimum bulk modulus at micro level
183
PARAMETERS FOR BEXM – INACTIVE2 ICL=84, ITYCL= 1 (INACTIVE)
HISTORY VARIABLES:
The Barcelona Expansive Model (BExM) (ICL=80 to 88) model requires the following
history variables:
Hist_var 1 P0* MPa Initial preconsolidation mean stress for saturated soil
Hist_var 2 γSD MPa Initial value of the parameter that defines the position
of SD yield surface
Hist_var 3 γSI MPa Initial value of the parameter that defines the position
of SI yield surface
Hist_var 4 nm Initial microstructural porosity
(micropore volume Vpm / total volume VT)
184
Initial values of P0* defining the initial position of LC surface are required.
By default, initial position of the surfaces γSD and γSI is fixed using the initial stress state
involving no elastic region between SD and SI surface. To do this, first γSI is computed as γSI =
p+smicro, then, the position of γSD is computed using an initial separation between SI and SD
surfaces which is prescribed as an input parameter in P4 of the ICL=88. If user wants to control
the initial position of γSD and γSI, then, specific values should be introduced as initial conditions
on surfaces/volumes, otherwise, these values should be keep as zero.
The evolution of history variables can be visualized as an output in Post-process GID interface.
Note: Effective stresses plotted in the Post-process GID interface correspond with net stresses
for unsaturated conditions and Terzaghi's effective stresses for saturated conditions. Stress and
strain invariants follow the soil mechanics notation (positive for compression).
185
APPENDIX: ABOUT THE NOTIONS OF MICRO- AND MACRO-STRUCTURES IN
THE DOUBLE-STRUCTURE MODEL IMPLEMENTED IN CODE_BRIGHT
The double structure models have been historically developed in order to reproduce the
behavior of unsaturated expansive clays (Gens & Alonso, 1992; Alonso et al., 1999). In this
seminal model, two levels of structure are considered (see Figure VIe.1):
1. The microstructure, corresponding to the clay particles made of active minerals take
place. As such, the microstructure is provided with a reversible strain-stress
relationship derived from considerations about double-layer theory:
dϵv = 𝛽𝛽𝑚𝑚 𝑒𝑒 𝛼𝛼𝑚𝑚 𝑝𝑝 𝑑𝑑𝑑𝑑
2. Macrostructure: Responsible of the structural rearrangements. At this level the
response of collapse and loading occur. The relation of stress and strains is defined by
the BBM model (Alonso et al. 1990).
A clear picture of what is the micro and macro-structure can be illustrated by looking at the
pore size distribution obtained in FEBEX bentonite (Figure VIe.2). Macropores corresponds in
this case to voids with entrance radii close to 30 µm while the microstructure has radii around
70 Å.
Even so, this model does not refer to unique process, or size of pores. The model is able to
reproduce other types of problems. The mathematical formulation for double structure soils
presented by Sánchez et al. 2005, is referred to the FEBEX bentonite, which present a clear
double structure evidenced by the pore size distribution test presented in Figure VIe.2.
186
Figure VIe.2 Pore Size Distribution of FEBEX Bentonite (Sánchez et al. 2005)
Figure VIe.3 Pore Size Distribution of compacted Jossigny silt (Casini et al. 2012)
Formulation can be also used to model types of materials, provided two main pores families
can be observed in the pore size distribution curve. An example is for example provided by
the work on compacted silty clay. In this case, the macrostructure refers to the arrangement of
silt particles and the large pores between them with entrance radii of order of 10 µm. The
microstructure refers to the clay particles and with entrance radii close to 1 mm.
Double-structure model can be used in this case to reproduce the fact that the stress dependency
of microstructure is much lower than that of macrostructure. In this case, as the scale of the
microstructure corresponds to the size of clay particles, the phenomenological law considered
is taken from classical expression for soils:
187
(1 + 𝑒𝑒𝑚𝑚 )𝑑𝑑𝑑𝑑
dϵv =
𝜅𝜅𝑚𝑚
Model can be also used to model three-level structures materials as, for example, a mixture of
bentonite powder with high density bentonite pellets. Three main pore families have been
detected in these types of materials, as shown in Figure VIe.4.
Figure VIe.4. Pore Size Distribution of mixture Bentonite powder/high density bentonite
pellets (Alonso et al., 2011)
In this case, the use of the double structure model needs an arbitrary split of the pore size
distribution into two pore family, according to the requirements of the modelling
(compressibility, permeability changes, …).
As a conclusion, concepts of micro and macro-structure define essentially two different levels
of scale within the material, not related to absolute value of pore size. For this reason, this type
of models can also be used to model fissured materials when the upper scale corresponds to
fissures and the lower scale to matrix.
188
CODE_BRIGHT. CONSTITUTIVE LAWS
This chapter contains the different models available and the corresponding parameters required
by each model. The following constitutive laws are available:
189
CASM's FAMILY MODELS
CODES in ICL=90, 91 (Mechanical data 3) ITYCL=1, 2 ,3, 4
ROOT_gen.dat
DESCRIPTION Clay and Sand model for soils (CASM)
EQUATIONS ICL=90 (CASM – General) ITYCL=1: Saturated CASM model
Yield function:
The yield function for CASM (Yu, 1998) expressed for a general stress state
n
3J 1 p'
is: = f + ln '
Mθ p ' ln r P0
where: p’ is the mean effective stress
1/ 2
1
=
2
(
)
J trace σ ij − p ' δ ij (for axi-symmetric conditions: q = 3 J )
'
P0 is the preconsolidation pressure
M θ is the slope of the critical state line
n is a constant used to specify the shape of the yield surface
r is a spacing ratio introduced to control the location of the intersection of
the critical state line with the yield surface.
where M is the slope of the CSL under triaxial compression ( θ = −30º ). The
parameter α controls the difference of the strengths between the triaxial
3 − sin φcs
compression and extension. Often, α is taken as: α = , where
3 + sin φcs
φcs denotes the friction angle of the soil at the critical state.
190
EQUATIONS ICL=90 (CASM – General) ITYCL=1: Saturated CASM model
(continued)
Hardening parameter:
Similar to the Cam-clay models, in CASM the change in size of the yield
surface ( dP0' ) is assumed to be related to the incremental plastic volumetric
by:
P0'
'
dP = *
0 d ε vp
λ −κ *
λ κ
where: λ* = and κ* = are the
1 + e0 1 + e0
modified compression index and modified elastic swelling index.
Plastic potential:
The plastic potential for CASM follows the stress-dilatancy relation of
Rowe (1962).
p' 2 3J 3J
=g 3M θ ln + ( 3 + 2 M θ ) ln + 3 − ( 3 − M θ ) ln 3 −
ζ p' p'
Where ζ is a size parameter, which can be determinate easily for any given
stress state (p’,J) by solving the above equation and it is internally
computed.
Elastic behaviour:
The elastic behaviour of this critical state model is the same as in the Cam
–clay models with the tangent modulus (K) and shear modulus (G) being
defined by the following expressions (a constant Poisson’s ratio (ν) is
assumed):
p' 3 (1 − 2υ ) K
K= * G=
κ 2 (1 + υ )
191
EQUATIONS ICL=91 (CASM – Specific) ITYCL=2: Unsaturated CASM model
This model require P1 to P9 parameters of the ICL=90 and ICL=91.
Model follows the Barcelona Basic Model (BBM) characteristics. The
model is formulated using two alternative options for the constitutive
stresses variables: Net stress and Bishop’s or average stress, which can be
selected using an indicator parameter.
- Bishop’s stress: σij =σij − Pg δij + S r sδij
'
Yield function:
In terms of Bishop’s stress:
n
3J 1 p'
= f + ln
M θ p ' ln r Pc
In terms of Net stress:
n
3J 1 p + ps
= f + ln
M ( p + p ) ln r P + p
θ s c s
where, r* and β are material parameters. The first is related to the maximum
stiffness of the soil (for an infinite suction), =r * λs ( s → ∞ )/ λo , and the
second controls the rate of increase of soil stiffness with suction.
192
Net stress formulation requires an explicit variation of apparent cohesion
with suction. The increase in cohesion follows a linear relationship with
suction,
ps = k s s
where ks is a scalar variable.
Plastic potential:
Three options are available for the flow rule. An scalar variable is
introduced to select a particular option as follow,
2 M θ ( 6 − M θ ) − ( 3M θ )
n n
λ*
with α= n −1
3 ( 6 − M θ )( 3M θ ) λ − κ*
*
Hardening parameter:
Isotropic hardening is controlled by the plastic volumetric strains ( dε vp )
through,
Po (1 + e )
dPo = d ε vp
λo − κ
Elastic behaviour:
Elastic behaviour is the same as defined in ICL=90 ITYCL=1. However,
net stress formulation requires an explicit relation to consider the effect of
suction on volumetric elastic strains, through the incorporation of the
elastic compressibility parameter for changes in suction, κs ,
κs 1
= dε ije,s = ds δij ds δij
(1 + e ) ( s + patm ) 3K s
Where, K s , is the Bulk modulus for changes in suction. If Bishop’s stress
is configured, the effect of suction on volumetric elastic strains is accounted
for via the variation in Bishop’s stress with suction. Elastic strain
increments are computed using the following equations for Net stress
approach and Bishop’s stress approach, respectively.
2G - 3 K 1 1
dε ije = dpδ ij + dσ ij + dsδ ij
6 KG 2G 3K s
2G - 3 K 1
dε ije = d ( p + Sr ⋅ s ) δ ij + dσ 'ij
6 KG 2G
193
CODES in ICL=90, 91 (Mechanical data 3) ITYCL=1, 2 ,3, 4
ROOT_gen.dat
DESCRIPTION Clay and Sand model for soils (CASM)
EQUATIONS ICL=91(CASM – Specific) ITYCL=3: Cemented CASM model
This model require P1 to P9 of ICL=90 and P1 to P7 of ICL=91
Cemented CASM model is based on the formulation for cemented soils
proposed by Gens & Nova (1993), in which a new state variable denoted
as ‘bonding’ b is incorporated.
Yield function:
3J 1
n
( p ' + pt ) Pc' Po' (1 + b )
=
=f + ln
M θ ( p '+ pt ) ln r ( Pc' + pt ) pt' = Po' (α t b )
Where, b is a non-dimensional variable that represents the degree of
' '
bonding. Po is the preconsolidation pressure of the unbonded material. Pc
'
controls the yielding of the bonded soil in isotropic compression and pt is
related to the cohesion and tensile strength of the material; αt is a parameter.
'
Both Pc and pt increase with the magnitude of bonding. The unbonded
behaviour is recovered when b goes to zero.
The function defining the reduction of bonding (b) with increased
degradation and the relationship controlling the evolution of degradation in
response to plastic strains, are:
b = b0 e −( h − h0 )
=dh h1 d ε vp + h2 d ε qp
The above expressions ensure that degradation increases monotonically,
independently of the sign of the plastic strains. h1 y h2 are material
parameters (greater than zero) defining the rate of degradation. d ε v and
p
d ε qp are the plastic volumetric strain increment and plastic deviatoric strain
increment, respectively.
The three options for the plastic potential used in the ICL=91 ITYCL=2
(Unsaturated model) are also available for this model.
Hardening parameter:
Either a volumetric hardening law or a combined volumetric and shear
hardening law can be used for this propose. If a combined hardening is
dP0' 1
adopted then,: = d ε p + ω d ε qp
P0'
λ −κ* v
*
'
hardening parameter ( P0 ).
Elastic behaviour is the same as defined in ICL=90 ITYCL=1.
194
CODES in ICL=90, 91 (Mechanical data 3) ITYCL=1, 2 ,3, 4
ROOT_gen.dat
EQUATIONS
ICL=91 (CASM – Specific) ITYCL=4: Double hardening soil model
(DHSM)
This model require P1 to P9 of ICL=90 and P1 to P8 of ICL=91
The model involves two plastic mechanisms: a plastic volumetric-driven
mechanism with isotropic hardening, by means of the use of the CASM
yield surface (ICL=90), and a plastic shear-driven mechanism that gives
rise to a nonlinear stress-strain relationship of a hyperbolic type, by the use
of a shear yield surface based on the Hardening Soil Model (HSM) (Schanz
et al., 1999).
Yield functions:
- Volumetric yield surface: CASM surface (see ICL=1, ITYCL=1)
- Shear yield surface:
1 J 2J
=fs − −γ p
E50 J Eur
1 −
Ja
with:
Jf ( p '+ c ' cot φ ' ) g(θ ) sin φ '
=
Ja = and θ)
g (=
Rf Rf sin θ sin φ '
cosθ +
3
J a is the asymptotic value of the shear strength, Rf stands for the relation
between the ultimate deviatoric stress, J f , and the asymptotic stress, a
standard setting is Rf = 0.9. γ p is the hardening parameter of the shear yield
surface and corresponds to the plastic deviatoric strain ( ε d ), as is described
p
P'o
p'
195
E50 is a secant stiffness modulus for primary loading and Eur is a secant
stiffness modulus for elastic unloading and reloading, these modulus are
computed as,
m m
p '+ c ' cot φ ' p '+ c ' cot φ '
=E50 E=
ref and Eur Eurref ref
ref
50
p + c ' cot φ ' p + c ' cot φ '
σ1 - σ3
qa asymptote
qf failure line
E50
Deviatoric stress
Eur
Axial strain ε1
Plastic potentials:
A non-associative flow rule is used for the shear yield surface which has
the following form,
d ε vp = g (ψ m ) d γ p
where, g (ψ m ) , is a function of mobilized dilatancy angle, ψ m , and
Lode’s angle, θ , as,
sinψ m
g (ψ m ) =
sin θ sinψ m
cosθ +
3
The mobilized dilatancy angle, ψ m , is defined as,
sin φm − sin φcs
sinψ m =
1 − sin φm sin φcs
196
φcs is the friction angle at the critical state and φm is the mobilized
friction angle. Material contracts for small stress ratios φm < φcs , while
dilatancy occurs for high stress ratios φm > φcs . Considering dense
materials contraction is excluded by taking ψ m = 0 for a mobilized friction
angle φm < φcs .
The volumetric surface uses a non-associated flow rule, the plastic potential
is defined in the ICL=90 ITYCL=1.
Hardening parameters:
d ε dp = 2 dJ 2ε
{ }
1/2
2
( ) ( )
+ ( d ε xp − d ε zp ) + d ε yp − d ε zp + d γ xy2 + d γ xz2 + d γ yz2
2 2 2
d ε=p
d ε xp − d ε yp
d
3
ε
where, J 2 is the second strain invariant of deviatoric strain tensor, γ p is the
integral of plastic deviatoric strain rates.
The hardening parameter of the volumetric yield surface is the so-called
'
preconsolidation pressure ( P0 ) defined in the ICL=90, ITYCL=1.
The two surfaces can be activated simultaneously and move together as in
a multi-surface plasticity problem, if the stress path reaches the intersection
between both surfaces.
Elastic behaviour of the DHSM is the same as in basic CASM model
defined in the ICL=90, ITYCL=1.
where, Aφ and Ac are two input parameters. It these parameters are zero
the basic DHSM is recovered.
197
PARAMETERS FOR CASM – GENERAL ICL = 90 ITYCL=1 (General: required for ICL=91,
ITYCL=2,3,4)
P1 ν - Poisson ratio
P2 κ - Slope of unload/reload compression curve
P3 λ - Slope of the normal compression curve
P4 r - Spacing ratio
P5 n - Shape parameter
P6 M - Slope of Critical State Line
Friction angle at CS (computed as a function of M). If φcs=0,
P7 φcs º
the shape of the YS plots as a circle in the deviatoric plane)
P8-P9 - -
P10 Su MPa Undrained shear strength (optional). By default = 0
198
PARAMETERS FOR CASM – SPECIFIC ICL = 91 ITYCL=3 (Cemented)
P1 b0 - Initial bonding. Evolution of bonding is stored as a history variable
P2 h0 - Degradation threshold
P3 h1 - Degradation rate for compression
P4 h2 - Degradation rate for shear
P5 αt - Parameter for tensile strength
P6 ω - Contribution of the plastic deviatoric strain to hardening parameter
Flow rule indicator:
1: Non-associated: Rowe
P7 Flow_rule -
2: Non-associated: K0 fit. Parameter α is required
3: Associated
P8 α - Non-associated parameter in case of flow_rule option equal to 2.
199
HISTORY VARIABLES:
The CASM’s family of constitutive models (ICL =90, 91) use the following history variables.
Output variables can be visualized in the Post-process interface of GID.
Model Hist_var Description Type
ICL = 90 1 P0 (MPa): Evolution of preconsolidation pressure Input/Output
(ITYCL 1) 2 F: Value of the yield function Internal
General
CASM 3 e: Void ratio Internal
model 4 Id_F1: Indicator of plasticity Internal
1 P0 (MPa): Evolution of preconsolidation pressure Input/Output
2 s: Suction Output
ICL = 91 3 e: Void ratio Internal
(ITYCL 2) Pc (MPa): Evolution of preconsolidation
4 Output
Unsaturated pressure due to suction
CASM ps (MPa): Evolution of tensile strength due to
model 5 Output
suction
5 F: Value of the yield function Internal
6 Id_F1: Indicator of plasticity Internal
1 P0 (MPa): Evolution of preconsolidation pressure Input/Output
Pc (MPa): Evolution of preconsolidation
2 Internal
pressure due to bond
ICL = 91 pt (MPa): Evolution of tensile strength due to
3 Output
(ITYCL 3) bond
Cemented 4 b: Evolution of bonding Input/Output
CASM
model 5 h: Evolution of degradation rate Output
6 F: Value of the yield function Internal
7 e: Void ratio Internal
8 Id_F1: Indicator of plasticity Internal
P0 (MPa): Evolution of preconsolidation pressure
1 Input/Output
(volumetric surface)
γp : Evolution of plastic shear strain
2 Output
(deviatoric surface)
ICL = 91 3 F1: Value of the volumetric yield function Internal
(ITYCL 4)
Double 4 F2: Value of the shear yield function Internal
hardening 5 c: Evolution of cohesion with plastic strain Internal
CASM 6 φ: Evolution of friction angle with plastic strain Internal
model
7 e: Void ratio Internal
Id_F1: Indicator of plasticity (volumetric
8 Internal
surface)
9 Id_F2: Indicator of plasticity (deviatoric surface) Internal
200
The input variable (P0) is introduced as initial condition on surfaces/volumes in the conditions
window of GiD. Void ratio (e) is computed internally as a function of porosity.
The evolution of some history variables can be visualized as an output in Post-process GID
interface.
Note: Effective stresses plotted in the Post-process GID interface correspond with Terzaghi's
effective stresses for saturated conditions (ITYCL=1,3,4). For ITYCL =2 (Unsaturated CASM
model) if the indicator iunsat (P6 of ICL=91 ITYCL=2) is equal to 0, net stress are plotted, if
iunsat=1, Bishop’s effective stresses are plotted. Stress and strain invariants follow the soil
mechanics notation (positive for compression).
_________
201
CODE_BRIGHT. CONSTITUTIVE LAWS
This chapter contains the different models available and the corresponding parameters required
by each model. The following constitutive laws are available:
202
VI.g. EXCAVATION/CONSTRUCTION PROCESS
EXCAVATION/CONSTRUCTION
EQUATIONS
203
Figure VIg.1. Illustration of the two main steps followed during excavation problem
204
VI.h. THM DISCONTINUITIES
205
Conditions Equation for plastic parameters
206
CODE_BRIGHT. REFERENCES
Formulation and numerical methods
Allen, M.B. and Murphy C.L. (1986): A Finite-Element Collocation Method for Variably
Saturated Flow in Two Space Dimensions, Water Resources Research, 22, no. 11:
1537:1542.
Altunin, VV & Sakhabetdinov, MA (1972). Viscosity of liquid and gaseous carbon dioxide at
temperatures 220-1300 K and pressure up to 1200 bar. Teploenergetika, 8:85-89.
Celia, M.A., Boulotas, E.T. and Zarba, R. (1990): A General Mass Conservative Numerical
Solution for the Unsaturated Flow Equation, Water Resources Research, 26, no. 7:
1483:1496.
Chen, G., Ledesma, A. (2009). Coupled Thermohydromechanical Modeling of the Full-Scale
In Situ Test "Prototype Repository. Journal of geotechnical and geoenvironmental
engineering, vol. 135, núm. 1, p. 121-132.
Garcia, J.E. (2003). Fluid Dynamics of Carbon Dioxide Disposal into Saline Aquifers. PhD
thesis, University of California, Berkeley.
Gerard, P., Léonard, A., Masekanya, J.-P., Charlier, R. and Collin, F. (2010), Study of the soil–
atmosphere moisture exchanges through convective drying tests in non-isothermal
conditions. Int. J. Numer. Anal. Meth. Geomech., 34: 1297–1320. doi:10.1002/nag.866
Hughes, T. J. R. (1980): Generalisation of Selective Integration Procedures to Anisotropic and
Nonlinear Media, Int. J.Num. Meth. Eng. 15, 1413-1418.
Huyakorn, P.S. and Pinder, G.F. (1983): Computational Methods in Subsurface Flow,
Academic Press, Inc. ISBN 0-12-363480-6.
Huyakorn, P.S., Springer, E.P., Guvanasen, V. and Wadsworth, T.D. (1986): A Three-
dimensional Finite-Element Model for Simulating Water Flow in Variably Saturated
Porous Media, Water Resources Research, 22, no. 13 : 1790:1808}.
McCutcheon, S.C., Martin, J.L, Barnwell, T.O. Jr. 1993. Water Quality in Maidment, D.R.
(Editor). Handbood of Hydrology, McGraw-Hill, New York, NY (p. 11.3).
Milly, P.C.D. (1984): A mass-conservative procedure for time stepping in models of
unsaturated flow, in Proceedings Fifth International Conference on Finite Elements in
Water Resources, edited by J. P. Laible et al., Springer-Verlag, New York: 103-112.
Olivella, S., J. Carrera, A. Gens, E. E. Alonso, (1994). Non-isothermal Multiphase Flow of
Brine and Gas through Saline media. Transport in Porous Media, 15, 271:293
Olivella, S., A. Gens, J. Carrera, E. E. Alonso, (1996), Numerical Formulation for a Simulator
(CODE_BRIGHT) for the Coupled Analysis of Saline Media, Engineering
Computations, Vol. 13, No 7, pp: 87-112.
Olivella, S., J. Carrera, A. Gens, E. E. Alonso. (1996) Porosity Variations in Saline Media
Caused by Temperature Gradients Coupled to Multiphase Flow and
Dissolution/Precipitation. Transport in Porous Media, 25:1-25.
Pini, G. and Gambolatti, G. (1990), Is a simple diagonal scaling the best pre-conditioner for
conjugate gradients on super-computers? Advances in Water Resources.
Sonnelveld, P., (1989), CGS, A fast Lanczos-type solver for non-symmetric linear systems,
SIAM J. Sci. Stat. Comput., 10, 36-52.
207
Spycher, N, Pruess, K & Ennis-king, J (2003). CO2-H2O Mixtures in the Geological
Sequestration of CO2. I. Assessment and calculation of mutual solubilities from 12 to
100°C and up to 600 bar. Geochim. Cosmochim. Acta, 67:3015-3031.
Voss, C.I. (1984): SUTRA Users's Guide, U.S. Geological Survey, Water Resources
investigations, Report 48-4369.
Van der Vorst, Henk (1990), Iterative methods for the solution of large systems of equations
on super-computers, Advances in Water Resources, Vol. 13, No. 3, 137-146.
Viscous/creep models
Olivella, S.; A. Gens, J. Carrera, and E. E. Alonso (1993): Behaviour of Porous Salt Aggregates.
Constitutive and Field Equations for a Coupled Deformation, Brine, Gas and Heat
Transport Model. The Mechanical Behaviour of Salt III, Trans Tech Publications, pp.
269:283. ISBN 0-87849-100-7.
Olivella, S., A. Gens, J. Carrera, E. E. Alonso (1996) Analysis of Creep Deformation of
Galleries Backfilled with Porous Salt Aggregates. The Mechanical Behaviour of Salt
IV, Trans Tech Publications, 379:386.
208
Wood, D.M., 1990. Soil behaviour and critical state soil mechanics. Cambridge University
Press, Cambridge.
Other
Gens, A. , A.J. Garcia-Molina, S. Olivella, E. E. Alonso, F. Huertas. (1998) Analysis of Full
Scale In-situ Heating Test Simulating Repository Conditions. Int. Journal for Numerical and
Analitical Methods in Geomechanics. 22:515-548.
___________________________________
210