A 2D/3D FINITE ELEMENT SOFTWARE

The ADELI package is free of charge for an academic, non-profit use. A new user must complete and sign a licence agreement that specify the conditions for a legal use of the software, and send it to the first address on the front page. Please acknowledge the use of ADELI software in any published work by mentioning this document.

2. SUMMARY

ADELI is a FORTRAN 77 Finite Element software developed to model the thermo-mechanical behaviour of the crust and the lithosphere at geological time scales for 2D and 3D quasi-static problems. The analysis is performed in large strain using the concept of objective derivative. While the space is discretized using linear elements (triangles in 2d and tetrahedrons in 3d), the time approximation is done using an explicit finite method based on the Dynamic Relaxation Method, and more specifically on the algorithm proposed by Cundall (1988).

The main capabilities of the program are:

1. Meshes of various sizes and shapes can be automatically generated starting from the definition of their boundaries, and an arbitrary number of materials can be defined.

2. The rheology can be chosen elastic (linear compressible), elastoplastic (Von Mises or Drucker-Prager), viscoelastic (linear or non-linear Maxwell body), or a combination of two anelastic behaviour.

3. Thermal properties can be used in order to compute a transient or steady state thermal solution.

4. Body forces corresponding to a constant gravity field can be included.

5. Boundary conditions on the sides are given in term of velocities and/or stress on the mesh border, and in term of temperature and heat flow for the associated thermal problem.

6. Contact problems between bodies are treated using Coulomb friction via an implicit algorithm for steep contact and dry friction.

7. Surface erosion by diffusion and transport due to water circulation can be computed using a finite difference formulation coupled to mesh update.

8. Initial conditions can be adjusted for internal stress and temperature.

The output consists in nodal and elementary values in a single file that can be visualized at different time steps (X-windows / Postscript), or like a movie for the 2d version.

Syntax highlights

Unix/ADELI commands lines are in italic

NUMERICAL input parameters are in upper case

KEYWORDS input parameters are in bold upper case

file names are underscored

input / output listing are in Courier font 10

notes and legends are in Helvetica font 10

recent changes in the user’s guide are in red

less recent changes in the user’s guide are in blue

3. INSTALLATION

Although the program should be able to run under various platforms, we have only used it under Unix and Linux and we describe all the operations under these operating systems.

The steps for installation are the following:

1. Create a directory (named X) where you want to install the code;

2. Get the software package via ftp:

ftp ftp.dstu.univ-montp2.fr (login = anonymous and password = your e-mail

cd pub/GTS/jean/adeli

get tarfile_adeli3p4.z

get tarfile_adeli3d2.z

3. Unshrink the code to get a tarfile

gunzip tarfile_adeli3p4.z

4. Extract files

tar xvf tarfile_adeli3p4

5. Add the paths X/bin and X/com in the \PATH variable in the .cshrc file.

The software should then be ready to use for Solaris Sun systems. If the binary programs are not compatible with your operating system, a makefile can be used for SUN, Linux, IBM and SGI UNIX/AIX operating systems in order to rebuilt the binary code.

4. STARTING AN EXPERIMENT

First create a work directory where the experiment will be processed. Although it is not required, each experiment should be done in a specific directory. Also avoid using one of the subdirectory of X, all the files produced by the software will be in the work directory.

Initializing an experiment is done via the command "ia".

Try ia for self information.

Type ia X 2d or ia3d X 3d on the work directory.

The command copy the two input files jessai and iessai into the work directory, as well as the files aadeli.ini and kcnusc.pal (2d), and the files p2x_ini (3d). It also makes a symbolic link with the file X/doc/ADELI2d.MESSAGES or X/doc/ADELI3d.MESSAGES.

In this document, the input files jessai and iessai will be called by their generic names j-file and i-file, where -file is a five digit name.

The experiment can then be started using the command “ea”.

Try ea for self information and examples.

The command use four arguments:

1. the version name 3.4b (2d) or 3d.2 (3d);

2. a code for the operating system;

3. f (foreground) or b (background);

4. the name -file of the experiment.

For example, to run the 2d code in foreground on a SUN using the input files jessai and iessai, try:

ea 3.4b sun f essai

You should see on the screen first:

ADELI VERSION 3.4b (> 16 Mars 1998 )

NOM DE L'ESSAI (5 LETTRES) ?

LECTURE DES COORDONNEES

LECTURE DES CONTOURS

LECTURE DES PROPRIETES DES MATERIAUX

then later something like :

CALCUL DES MASSES NODALES

ITERATION \% TEMPS IREDUT \% PLAST \% INERTD \% INERT \% ERGLOB S/CONT S/STRESS

1: 1 0.00 0 0.00 1.22 1.21 0.00 2 0/ 0

2: 1 0.01 0 0.00 1.18 1.16 0.00 3 0/ 0

3: 1 0.01 0 0.00 0.99 0.97 0.00 3 0/ 0

and at the end:

SOIT : 0 HEURES

0 MINUTES

27 SECONDES

iflag = 0

ADELI_3.4b : TERMINAISON NORMALE

Such a sequence indicates that the code is technically working. The program has then created different files, among them:

pessai ( p-file ) which contains the output for plotting;

oessai ( o-file ) which contains the output for erosion and friction;

eessai ( e-file ) which may contain debugging messages.

If the program is used in background (recommended for large runs), the standard output is redirected in the file -file , here essai. Visualising the results is done with reading the p-file with aadeli or xadeli2d or xadeli3d

5. PARAMETRIZING AN EXPERIMENT

5.1 Introduction

Most of the physical information about the experiment is contained in the i-file , and changes have usually to be made only is this file. The j-file contains numerical parameters that modify the behaviour of some algorithms and also the memory allocation. The input files for 2d and 3d experiments are quite similar, as well as the meaning of the parameters. We describe in the following the use of the 2d files, and note the differences with the use 3d files when needed.

Note : see also the benchmark section and the experiment section on the web site to see various examples of j-file and i-file.

The structure of the j-file is quite simple, as always the same number of parameters are read. A typical file is (3d3 version) :

VERSION

3d.3

nvmax nprmax nfrmax nelmax npmax nmamax nthmax nflmax nvalmax

2000 1000 7000 9000 8000 10 1000 100 100

ndime3 npmax3 nnod3 nelmax3 nfacex nfrmax3

3 40000 4 100000 20000 10000

ncouchx nnode ndime ntype

3 3 2 3

fract redut tredut

4. 1.0 0.1

toler redtol nerrspa

0.88 0.9 0

nmixe nmixs idual nbifu

0 0 0 0

iobj icinem

1 1

inter nliss

1 2

inorm inormdev ciner cisup ciinf civar convth

2 1 .05e5 1.5 0.5 0.01 1.e-4

niter naugt

1 500

imdf damp ctraine

1 0.8 0.0

ioptmsh interp xmesh

0 0 1.05

xplotp xplott xplotc

0.5 0.5 0.5

npxt npyt xmint ymint xmaxt ymaxt

0 0 0.001 250 200 450

nfacep lfacep

1 5

gradmax devmin devmax redelm aigumin

0.66 5.0 30.0 0.01 8.0

The 2*I+1 lines of the j-file correspond to comments recalling the variable names below, and the 2*I lines correspond to variables. For a classical use, let the variables to their current values. However the following variables can be changed:

nvmax : maximum number of pressure boundary nodes;

nprmax : maximum number of velocity boundary edges;

nfrmax : maximum number of contour nodes;

nelmax : maximum number of elements;

npmax : maximum number of nodes;

nmamax : maximum number of materials;

nthmax : maximum number of thermal boundary nodes;

nflmax : maximum number of heat flow boundary edges;

nvalmax : maximum number of some scalar arrays;

npmax3 : maximum number of velocity boundary nodes;

nnod3 : maximum number of 3D mesh nodes;

nelmax3 : maximum number of 3D mesh elements;

nfacex : maximum number of 3D surface faces;

nfrmax3 : maximum number of 3D mesh contour nodes.

These variables control the dynamic memory allocation of numerous arrays (~100), which is done within a large array in the main program. This array has a maximum value of WMAX0 setup as a parameter value. Increase it and recompile the code if more memory is needed.

The other variables of the j-file that can be frequently changed are xplotp xplott xplotc. They control output in the p-file and o-file for global output (xplotp), topographic output (xplott) and contact output (xplotc). They correspond to the normalized time interval for which you want to have an output. To be clear, 0.5 will provide 3 output at 0*TFIN, 0.5*TFIN and 1.0*TFIN.

The structure of the i-file corresponds to some integer parameters at the beginning, which controls the use of different blocks below (mesh, boundary conditions, ...). The text above each parameter is a comment that reminds the name of the variable in the program (as it is in the j-file). Do not remove it. All the descriptions refer to these names. All the integer and real variables are read in free format. Each block is limited by a keyword that is identified by the program. The end of file also correspond to a keyword (i.e. COORDONNE) of 9 characters. Any misreading will cause the program to abort with a error message (see Appendix ) and a non zero value for iflag at the end of the standard output.

The first parameter of the i-file is IECHO. This parameter is used for debugging, and the output is written in the e-file. However the use of the value 102 writes an echo of the input files in the o-file. We recommend to use this value. IMECA controls the mechanical solution, which will be computed if IMECA = 1. ITHERM controls the thermal solution, which will be computed if ITHERM = 1.

Note : The thermal solution is not computed in 3d.

5.2 Mesh generation

The mesh generation is done starting from the contour of the future mesh. This contour corresponds to a 1D mesh of the material(s) you want to define.

Figure 5.2.1 : Example of edges definition of two materials I and II. Large numbers refer to points numbers, circled numbers to edges.

The linear element size is computed using the variable NELEM which defines the approximate number of elements of the future mesh. The definition of the contour edges corresponds to the position of NFRON0 initial points defined 2 lines after the keyword COORDONNE. Each line contains 1 integer (the point number), ranging from 1 to NFRON0, and the x and y position as real numbers in a Cartesian system (see figure 5.2.1 and table 5.2.1):

Table 5.2.1 : Edges definition in the i-file

COEF

1000.

COORDONNEES DE NFRON0 POINTS DU FRONT

numP X Y

1 0 0

2 20 0

3 20 0

4 30 0

5 0 -15

6 30 -15

7 10 -10

8 0 -20

9 30 -20

x and y units are defined by the coefficient COEF which multiply x and y values to produce values in S.I. (meters). For example, using 1000 for COEF leads to have x and y in kilometres in the file, that is convenient for crustal scale problems.

Note : The coordinates contour is the only field using non S.I. values as an input.

Each material is described by a list of edges corresponding to successive points. The idea is to define the contour of each material by writing the points numbers in a trigonometric direction (counter clockwise rotation). This list corresponds to NFACE0 lines defined 2 lines after the keyword CONTOURS_ ( _ is a blank). Each line contains an integer (the edge number) increasing from , 1 to NFACE0, another integer corresponding to the first point of the edge, and three other numbers that can be used to define curved edges (set the first one to 0 to have linear edges). The next line describe the next point of the contour, etc, until the point numbered -999 is reached. This means that the current material is described, and the last non-negative point (the line above the -999 line) will connect to the first line of the sequence. The next material is described after the -999 line using the same convention. However, a rule of thumb is that the first edge of the next material should be attached to a previously defined material if you intend to connect to it. The contour of figure 5.2.1 corresponds to the following sequence of table 5.2.2.

Table 5.2.2 : Contour edges definition in the i-file

CONTOURS DE MATERIAUX (nface0 valeurs)

NUMF POINT JCERC TGIP1 TGIP2

1 2 0 0 0

2 1 0 0 0

3 5 0 0 0

4 6 0 0 0

5 4 0 0 0

6 3 1 0.1 0

7 7 1 -0.1 0

8 -999 0 0 0

9 6 0 0 0

10 5 0 0 0

11 8 0 0 0

12 9 0 0 0

13 -999 0 0 0

The sequence of the list of edges will end with a line containing the number NFACE0 and -999 (plus three other number). The number of materials will therefore correspond to the number of -999 encountered in the list of edges.

Some curvature may be imposed to the edges with coding the third value JCERC of the corresponding line to values greater than 0. The following possibilities exist: JCERC = 1 creates a sinus arc with the value of TGIP1 is the ratio between the sinus arc arrow and the distance between the starting point and the ending point of the sinus arc. Positive values of TGIP1 create curvature inward the mesh, negative values outward the mesh. JCERC = 2 creates a circle arc with the same convention. JCERC = 3 creates a fourth order polynomial defined by the orientation of the tangents at the starting and ending points. TGIP1 is the orientation (in degrees) of the tangent at the starting point with respect to the orientation of the straight line between these two points. Positive values of TGIP1 correspond to the trigonometric direction. The same convention apply to the ending point using TGIP2. TGIP1 = TGIP2 = 0 creates a straight line. JCERC = 4 creates an ellipse with TGIP1 = 0.5 * small diameter / large diameter.

Note: that the use of large curvature will affect the linear element size. This can lead to large element size, specially with JCERC = 4.

Running the code with the values of tables 1 and 2 and a value of 1000 for NELEM should produce the contour and mesh of figure 2. The elements of the mesh are generated using an algorithm that create new points within the contour of a material, forms the associate elements, and reduces the current inner contour according to the newly defined elements.

Figure 5.2.2: Mesh generated according to tables 5.2.1 and 5.2.2.

The mesh of the material is completed when the contour corresponds to an empty sequence of nodes. The mesh generation is automatic and is only affected by the initial contour, NELEM, NMIXE and COEF_MESH. The purpose of NMIXE is to create a mesh with pairs of triangular elements in order to use an averaged value for the volumetric strain rate. This option is activated when NMIXE = 3, and should be used in conjunction with elasto-plastic rheologies like Von Mises or Drucker Prager plasticity. The mesh of figure 5.2.2 is created in such a way. Employ NMIXE = 0 not to use this option. COEF_MESH is found in the rheological block for each material, and allows to multiply the current element size by COEF_MESH. This way it is possible to generate different element sizes.

3d mesh generation

The mesh generation at 3d is based on the expansion of a 2d mesh in the direction orthogonal of the 2d plane. This expansion corresponds to the creation of one or more layers of variable thickness, that contain tetrahedral elements. We assume that the x-y plane is horizontal, and that z increases when elevation increases. The primary 2d mesh generation follows the same rules as above. However, two options exist, controlled by the parameter IMAP3: If IMAP3 = 0, the primary 2d plane is assumed vertical (x-z plane), and the direction of expansion is horizontal along y. This option may be used if you want to start from a cross-section. If IMAP3 = 1, the primary 2d plane is assumed horizontal (x-y plane), and the direction of expansion is vertical along z.

The controlling block for the mesh size appears in table 5.2.3

Table 5.2.3 : 3D-mesh element size control in the i-file

imapg npmapg

0 30

tmapgmin telem0 (echelle utilisateur)

25 50

dmapcrit (echelle utilisateur)

150

liste des npmapg points de controle (x,y, echelle utilisateur)

1900.000 700.000

1840.100 721.329

1779.499 740.575

1718.269 757.715

1656.483 772.729

1594.215 785.598

1531.539 796.307

1468.531 804.844

1405.265 811.197

1341.818 815.360

1278.264 817.328

1214.681 817.098

1151.143 814.671

1087.728 810.049

1024.510 803.238

961.565 794.245

898.968 783.083

836.795 769.764

775.119 754.304

714.014 736.721

653.554 717.037

593.810 695.275

534.854 671.461

476.756 645.624

419.586 617.794

363.413 588.005

308.302 556.292

254.320 522.693

201.532 487.248

150.000 450.000

The linear element size is directly controlled using the parameter TELEM0. Two options are useful:

1) If imapg = 1 the code will use an element size that will vary between TMAPGMIN (around NMAPG control points, with a damping parameter DMAPCRIT) and TELEM0 far from these points.

2) If imapg = 0 the code doesn’t take into account of the parameters described above and will use telem0 as an element’s size.

Note that all distances are controlled by the variables COEF.

The coordinate and the contour blocks (keywords COORDONNE and CONTOURS_ are the same as for the 2d version. A new block corresponding to the creation of layers is controlled by the keyword COUCHES__ and has the structure as shown in table 5.2.4. NCOUCH defines the number of material layers, and is followed by NCOUCH+1 lines that define each layer. A line contains the layer number ICC, the out of plane coordinate ZCOUCH of the top of the layer , and the number of sublayers KCOUCH of the layer ICC. A line with ICC > 1 also defines the bottom layer coordinate of layer ICC-1. The last number (RCOUCH) allows to define a variation of the sublayers thickness according to the expression thickness = RCOUCH**(sublayers number – 1).

Table 5.2.4 : 3D-mesh layers definition in the i-file

COUCHES

NCOUCH / NCOUCH+1 LIGNES / ICC,ZCOUCH(ICC),KCOUCH(ICC),RCOUCH(ICC)

1 0.0 2 1.0

2 4e3 1 1.0

The 3d structure is produced using the table 5.2.4. The faces of the 3d mesh are 2d surfaces that are numbered in a similar way as for edges in 2d as shown in figure 3. The understanding of this faces numbering is needed to setup the boundary conditions. First are numbered the faces out of 2d plane (lateral faces), i.e. the faces normal to the plane xz in the example. A face with a given number contains the same points that the edge of the same number. For example, the face 3 contains the same points that the edge 3 in 2d, that are 5 and 6. As a result, if one layer is used, the number of lateral edges is equal to NFACE0-1 (because the last edge number NFACE0 is just used to close the last material. If other layers are used, the numbering follow the same principle, and the number of lateral edge is given by NCOUCH*NFACE0-1.

Note: the maximum allowed layers number is 99 in this present compiled version.

Then are numbered the faces corresponding to the primary 2d mesh, starting from NCOUCH*NFACE0 + 1 for the first material, until NCOUCH*NFACE0 + NMATS where NMATS is the last material. Finally are numbered the faces at the bottom of the last layer, starting from NCOUCH*NFACE0 + NMATS+1 for the first material, until NCOUCH*NFACE0 + 2*NMATS for the last material.

Figure 5.2.3: 3D contour numbering according to table 5.2.3.

The vertical position of the nodes of a 3D mesh can be changed in order to fit an isostatically compensated topography. The option is used if the value of the variable ISOMESH equals 1. This field must start by the keyword ISOSTASIE, followed by the name (5 characters) of the file (table 5.2.5) containing the topographic values. The first line contains the x-y origin (lower-left corner) for the altimetric values. The second line contains the span in x and y. The third line contains the number of points nxiso and nyiso in x and y. On the fourth line starts the list of the nxiso*nyiso altimetric values. All values are given in meters.

Table 5.2.5 : Description of the topography in a 5-character file ( ISOMESH=1)

600000 400000 xiso yiso

400000 200000 dxiso dyiso

4 4 nxiso nyiso

0 0 0 0

0 1000 1000 0

0 0 0 0

5.3 Constitutive laws

The parameters of the constitutive laws of the different materials are defined in the block starting with the keyword RHEOLOGIE. The 6 following lines are comments that recall the meaning of rheological parameters. Then each material is defined by 6 lines (see table 5.3.1): Line 1 contains material number NUMAT, type of constitutive law ITYP, and the mesh parameter COEF_MESH, that will multiply the current element size. The next 5 lines must contain each five parameters, even if they are not used (nul = 0). Line 2 contains elastic parameters. YOUNG is the Young modulus (Pa), POISS the Poisson coefficient, RHO0 the density (kg/m²).The 2 others parameters are not used. Line 3 contains viscous parameters. GAMM0 is the multiplicative term (Pa^-n s^-1) of the fluidity g. EACTI is the activation energy E_a of the power law (J/mol), EXPOS is the exponent n of the power law, and TKELV is a limit temperature. The resulting fluidity γ is computed according to

where T is a temperature defined by the statement T=min(T,TKELV). Line 4 contains plastic parameters of a Drucker-Prager material. COHES is the material cohesion in Pa. PHINI is the initial friction angle (degrees). PHFIN is the friction angle after a phase of strain softening. KAPPAC is the limit of the equivalent plastic strain after which the friction angle f is set to PHFIN. Set PHFIN=PHINI to avoid f to evolve whatever the plastic strain be. PSI is the dilatancy angle in degres. Line 5 contains plastic parameters of a Von Mises material. SEUIL is the initial yield stress (Pa). HARD is a hardening/softening coefficient that means that the yield stress may evolve between SEUIL and SEUIL + (max(kappa,KAPPAC)/KAPPAC) * HARD*SEUIL. KAPPAC is the limit of the equivalent plastic strain after which the yield stress do not evolve. The 2 others parameters are not used. Line 6 contains thermal parameters, that are used if ITHERM = 1. CSPEF is the specific heat. CONDT is the thermal conductivity. SRCTH is an internal source term (W/m³). DILAT is the thermal expansion coefficient. The last parameter is not used.

Setting ITYP allows to use a combination between elasticity and the three other constitutive equations for rheology:

ITYP = 1 involves only elasticity (line 1);

ITYP = 2 involves Maxwell viscoelasticity (lines 1+2);

ITYP = 3 involves Drucker-Prager elastoplasticity (lines 1+3);

ITYP = 4 involves Von-Mises elastoplasticity (lines 1+4);

ITYP = 5 involves Drucker-Prager elastoplasticity OR Maxwell viscoelasticity according to the current stress state (lines 1+2+3);

ITYP = 6 involves Von-Mises elastoplasticity OR Maxwell viscoelasticity (lines 1+2+4).

Table 5.3.1 : Example of a rheological setting

RHEOLOGIE DE CHAQUE MATERIAU :

NUMAT ITYP COEF_MESH

YOUNG POISS RHO0 NUL NUL

GAMM0 EACTI EXPOS TCUTT NUL

COHES PHINI PHFIN KAPPAC PSI

SEUIL HARD KAPPAC NUL NUL

CSPEF CONDT SRCTH DILAT NUL

1 5 1.0

1.e11 0.25 2.8e3 0.0 0.0

0.84912E-14 110e3 1.0 948. 0.0

1.e6 15. 15. 0.01 0.0

1.e8 -0.1 1e-2 0.0 0.0

1070. 3.0 0.0e-6 0.e-5 0.0

Note : Non linear viscoelasticity is now implemented in 3d (Riad Hassani update).

5.4 Velocity boundary conditions at 2D

Velocity boundary conditions can be imposed on the outer edges of the mesh with setting NVFIX0 to the number of edges you want to use. The corresponding block starts with the keyword VITESSES_, followed by 1 line of comment and NVFIX0 lines of variables. The series of edges that forms the boundaries can be listed in any order if they are not connected each others. If two or more edges are connected, they must be listed in trigonometric order. An example is given in table 5.4.1. The first integer of the line (NUMF) refers to the edge number. On this edge it is possible to constrain the normal and/or the tangential velocity, with setting the second integer (normal code) and the third integer (tangential code) to 0 (free displacement) or 1 (locked displacement). These integers are called locking codes. Note that the code 0 is working only if the edge is parallel to the x or y axis. The code 1 is working for any edge orientation. If the code 1 is used, then the two next real variables are used for respectively the normal and tangential velocity. The convention is here that a positive normal velocities are directed inward the material, and positive tangential velocities are directed in the trigonometric direction.

Time variable intensity

The last integer of the line refers to a possible time variable intensity of the velocity vector. If a value of 1 is given for boundary condition line, the number following “EVOLUTION DES SOLLICITATIONS” is used to read a specified number of lines after “VALEURS NORMALISEES DU TEMPS ET...”. Each of these lines (3 in the example) contains 2 real numbers and 1 integer. A given line i contains a relative time trel_i, a coefficient coef_i and a parameter itype_i.

At a given time t, we define the relative time alpha such as t/TFIN, and we must have two lines i and i+1 such that:

0 < beta = (alpha - trel_i)/(trel_i+1- trel_i) <1

coef(t) is then computed according to the value of itype_i:

If itype_i = 1 : coef(t) = coef_i;

If itype_i = 2 : coef(t) = (1-beta) coef_i + beta coef_i+1

If itype_i = 3 : coef(t) = (1-val) coef_i + val coef_i+1 with val = 0.5 [1 - cos(p beta)]

The velocity values are finally modified according to:

v = coef(t) v

Note : this procedure can also be applied to other scalar or vector variables such as pressure boundary condition, gravity, temperature, heat flow, friction coefficient and pluviometry, using the same convention.

Table 5.4.1 : 2d velocity boundary condition

VITESSES IMPOSEES

NUMF CODES_N-T VITESSE_N VITESSE_T EVOLUTION (O=1,N=0)

2 1 0 -1.0e-10 0. 1

10 1 0 -1.0e-10 0. 1

12 1 0 0.0e-10 0. 0

4 1 0 0.0e-10 0. 0

EVOLUTION DES SOLLICITATIONS (NBE DE POINTS NEVOLV)

VALEURS NORMALISEES DU TEMPS ET DE L'INTENSITE, TYPE DE FONCTION

0.0 0.0 1

0.5 0.5 1

1.0 1.0 1

5.5 Velocity boundary conditions at 3D

Coding the boundary conditions in 3d is almost similar than in 2d (see table 5.5.1). The differences are:

- 3 locking codes (CODES_XYZ) and 3 values have to be initialised for each face;

- These locking codes and values refer to boundary conditions expressed in the global coordinate system x-y-z (see figure 3);

- 6 velocity values allow to have a parametric variation along the face using the code of the 5^th column. A value of 0 for this code will use only v1_x v1_y v1_z for the face. A value code of 1, 2 or 3 will use a linear variation between v1_x v1_y v1_z and v2_x v2_y v2_z. A value of 1 uses x-coordinate as a variable, a value of 2 uses y-coordinate as a variable and a value of 3 uses the distance along the face as a variable.

- There is the possibility to define a spatially variable velocity field for a given edge using the parameter NFILEV in the header.

Table 5.5.1 : 3d velocity boundary condition

VITESSES IMPOSEES

coef. multiplicateur des vitesses imposees.

3.1709792e-11

numF codes_N-T1-T2 liaison v1_x v1_y v1_z v2_x v2_y v2_z evolution (o=1,n=0)

1 1 1 0 3 -16.3 -26.6 0 -11.5 -9.9 0 0

2 1 1 0 3 -11.5 -9.9 0 -9.3 8.5 0 0

3 1 1 0 3 -9.3 8.5 0 -14.0 12.5 0 0

4 1 1 0 3 -14.0 12.5 0 -3.0 5.0 0 0

15 1 1 0 3 -3.0 5.0 0 0 0 0 0

16 1 1 0 3 0 0 0 0 0 0 0

17 1 1 0 3 0 0 0 -15 -25 0 0

18 1 1 0 3 -15 -25 0 -15 -25 0 0

9 1 1 0 3 -15 -25 0 -16.3 -26.6 0 0

nom de fichier de conditions en vitesse (utilisé si nfilev = 1)

fich5

EVOLUTION DES SOLLICITATIONS (NBE DE POINTS NEVOLV)

VALEURS NORMALISEES DU TEMPS ET DE L'INTENSITE, TYPE DE FONCTION

0.0 0.0 1

0.5 1.0 1

1.0 0.0 1

This spatially variable velocity field can be set as follow: Put NFILEV to 1. After the line “nom de fichier de conditions…” must be the name in 5 character of a file (here “fich5”). This file has to be in the working directory and must contain a sequence as in the example of table 8. The first line of the file contains the face number where you want to apply a spatially variable boundary condition ( ivitnum). This face number has to be specified in the list of faces. The second line specify the directions ( ivit1 and ivit2) where the velocities are specified (will evolve) (1 for x, 2 for y, 3 for z). These two directions must correspond to locking codes set to 1 below the keyword CODES_XYZ. For example, the sequence of table 8 should correspond to the sequence “1 0 1” for CODES_XYZ of the face 15 (NUMF). The third line corresponds to the grid origin in the frame defined by ivit1 and ivit2. The fourth line corresponds to grid spacing dxvit and dyvit for the same frame. The fifth line corresponds to data number dxvit and dyvit for the same frame. After follows a sequence of 2* nxvit* nyvit real values which correspond to velocity values for ivit1 and ivit2 directions, starting from the origin and increasing in direction ivit1 then ivit2. For a given node of the face, the program will search which nodes of the grid allow to interpolate the value of velocity for each specified direction.

Note: When this option is in use, the corresponding values V_X, V_Y or V_Z (see table 5.5.1) are not used for the directions ivit1 and ivit2. However a time variable intensity can still be used as mentionned above.

Table 5.5.2 : Spatially variable boundary condition file in 3d

15 ivitnum

1 3 ivit1, ivit2

30.5e3 -14.9e3 xvito, yvito

29.5e3 15.0e3 dxvit, dyvit

2 2 nxvit, nyvit

5 0 5 0

5.6 Pressure boundary conditions at 2D

The pressure boundary conditions are defined in the same way than the velocity boundary conditions. The keyword is PRESSIONS, and the number of boundary conditions is given by NPRES0. The table 5.6.1 corresponds to NPRES0 = 1. Each line is composed by the edge number (integer), the type of pressure called itypp (integer), the two values of normal and tangential pressure presn and prest (real numbers), and a code which controls the time evolution (integer). If itypp = 0, the values presn and prest are used. The sign conventions are identical than for velocity boundary conditions in 2d. If itypp = 1, you may add to the static values defined by presn and prest an hydrostatic pressure P_h(y) that is defined by its density RHYDRO and a lithostatic prestress (see after the keyword PRECONTRA). The applied pressure is computed using the actual depth y (if the view in cross-section was chosen) along of the edge, the hydrostatic free elevation Y_h and the Oy gravity value GRAVY. Its value is given by :

P_h(y) = - RHYDRO*GRAVY*(Y_h - y)

The hydrostatic free elevation Y_h is computed with assuming the hydrostatic pressure at the depth YREFW is equal to the integration of the density layers given under the text ylitho and rlitho. An analytical integration is performed to compute the equivalent forces on the two adjacent nodes. If itypp = 2, the same kind of pressure calculation as for itypp = 1 is applied. However, this flag means that the corresponding face may appear dynamically as a pressure boundary condition because a part of the face n acts as a frictional boundary condition. The part of the face where the pressure boundary condition is active starts on the first node of the face and stops when this face becomes in contact with the next face of the list. The next face of the list must correspond to itypp = 1.

An example should be given in Appendix later....

note : As this option has been tested only on one example, the user should perform some simple test (hydrostatic test) before to get into complex geodynamical model.

Table 5.6.1 : Pressure boundary conditions at 2d

PRESSIONS IMPOSEES

numF itypp pression_N pression_T evolution (o=1,n=0)

11 1 0.0 0.0 0

5.7 Pressure boundary conditions at 3D

The use of pressure boundary condition in 3d is about the same as for 2d. The main difference is that the 3 arguments used to define the stress vector must correspond to x-y-z components ( p_x, p_y, p_z columns, see table 5.7.1). The use of the code itypp is the same as for 2d.

Table 5.7.1 : Pressure boundary conditions at 3d

PRESSIONS IMPOSEES

numF itypp p_x p_y p_z evolution (o=1,n=0)

1 1 0.0 0.0 0.0 0

5.8 Friction law at 2D

The use of friction law corresponds to define two or more zones that can interact between them according to a Coulomb law. This feature is controlled by the parameter NFRIC which defines the number of different zones that may interact. The zone description and the friction law is described below the keyword CONTACT_ . Each zone is defined by a line which includes (see table 5.8.1) the zone number, the number of edges belonging to the zone, and the list of these edges, keeping the convention of trigonometric description. After the text “nbr de visibilites mutuelles” comes the number of associations nvisi between the different zones. If you have only two zones this number is one, but may vary until A_nvisi* A_nvisi more than two zones are used. However, the software may be unable to compute contact reactions associated to some complex geometrical situations. After the text zone_cand_num... you have to define for each association 5 numbers:

- the master zone number;

- the slave zone number;

- a parameter irevers which allows to use a normal and revers

configuration for the master/slave zones association if irevers=1;

- the friction coefficient (must be positive);

- a possible time dependence of the friction coefficient ( ievol=1).

Note that a specific output for these friction zones may be setup using the parameter xplotc in the j-file.

Table 5.8.1: Definition of contact zones

CONTACT (tolerance, nbr max d'iteration)

1.e-3 50

numero, nbe faces, liste des numeros des faces (nfric ligne)

1 1 6

2 1 7

nbr de visibilites mutuelles

zone_cand_num1 zone_antg_num2 irevers friction ievol

1 2 1 0.3 0

5.9 Earthquake génération at 3D

The zone description and the friction law is described below the keyword SEISMES_ as shown in table 5.9.1.

Note : Earthquake génération is still experimental is the 3d3 version. Therefore, this option must be use with care.

A useful option of this block is to allow an output for a arbitrary number of surface points (NPQOUT) selected by their x-y initial coordinates. The frequency of output is monitored by ISPANQ (1 means output for each time step). If this option is used, it is also needed to specify the faces that form the surface with coding the parameters NFSUR and LFSUR in the erosion-sedimentation block.

Table 5.9.1: Definition of friction parameter for stick-slip behavior

SEISMES fricb dfric

0.10 0.03

timecos dtcos0 timepost ntpost

100 100 3.15e10 200

nom de fichier de deplacement cosismique max

depla

nombre, pas d'echantillonage et coordonnées des points de controle en surface

2 npqout

100 ispanq

1273e3 539e3

1600e3 600e3

5.10 Erosion-sedimentation laws

Modelling erosion and sedimentation is possible when the parameter NERO is set to 1 or 2. In these cases the parameters below the keyword EROSION__ are shown in table 12.

Table 12: definition of erosion parameters

EROSION (ibase hbase0 transp pluvio)

0 0.5e3 0.00 0e-8

XLF ICONLIM, DIFF0, PANTC, COENL

20000. 00 1.e-7 1.0 0.0

rhosedi qsed

2.3e3 0.0e-13

nombre et liste des faces (nfsur, lfsur(i),i=1,nfsur)

3 7 8 1

evolution des sollicitations (nbe de points nevolv)

valeurs normalisees du temps et de l'intensite, type de fonction

NIVEAU ET DENSITE DE L'EAU

Y_water rho_water

0.0e3 1.0e3

As for other surface boundary conditions, the edges list must be given following the trigonometric convention. The surface submitted to erosion can be crossed by faults. If NERO = 1, a law which models the small scale diffusion of the topography and the large scale transport by the rivers is used. It is assumed that the origin of mass redistribution is caused by the pluviometry. Therefore, the amount of water circulating on the topography depends on the value of the pluviometry (constant in space and time) and of the distance between the current point and the divide line. The option NERO = 2 should not be used.

Note : Erosion sedimentation laws are not written in 3D. However, the description of the surface with NFSUR and LFSUR is active and may be used for other purposes (i.e., output, initializing isostatic topography).

5.11 Initial stress

Modeling the lithosphere involves large models submitted to gravity. In these cases it can be useful to settup an initial stress that corresponds to the lithostatic stress in the lithosphere. Two possibilities exist and are controlled by the keyword ISOSTRES and the parameters encoded below the keyword PRECONTRA.

Table 5.11.1: definition of initial stress parameters

PRECONTRAINTES LITHOSTATIQUES (nlitho, ylitho et rlitho par couche)

ylitho rlitho

0.5e3 0.0

0.0e3 1.0e3

-5.0e3 2.6e3

-6.0e3 2.6e3

-7.0e3 2.6e3

rhydro yrefw

2.2e3 -5.0e3

If ISOSTRES = 1, then the values defined in the table 5.11.1 are used. NLITHO defines the number of lines after the text YLITHO RLITHO. Each line below this text contains respectively the elevation of an horizontal density interface and the density of the material above. Therefore, the first line should contain the elevation of the surface and 0.0 as the approximate air density. A vertical integration of this density structure is given for each integration point. If one wants to have a body in a perfect initial balance, it is necessary to have a perfect agreement between the layered density structure define in the table and the density distribution defined by the mesh and the material properties table. RHYDRO is the density used for hydrostatic restoring (see pressure boundary conditions)

If ISOSTRES = 2, then the table 13 is not used to compute the initial stress. Rather, a direct integration of the weight of the medium above the considered point is computed using the mesh and density information. Using this option allows to take into account the topographic weight. However it should be noticed that in this case the initial stress at a given depth may be not a constant.

Modeling the initial stress at 3D is done in a similar way. If ISOSTRES = 2, then LFSUR must be settup in the erosion-sedimentation block.

5.12 Initial temperature at 2D

Initial temperature is defined following the keyword TEMPERAT.

Table 5.12.1: definition of initial temperature parameters

TEMPERATURE INITIALE (ncoucht)

ytemp vtemp jtemp

1 273. 2

-50E3 773. 2

xtgaus ytgaus dxgaus dygaus tgaus

100e3 -10e3 100e3 5e3 0

nom du fichier de temperature (utilisé si nfilet = 1)

tempi

Table 5.12.2: definition of initial temperature in a user-defined file

10000 -50000 x0 y0

20000 20000 dx dy

2 2 (x,y)

673 673 673 673

5.13 Initial temperature at 3D

Initial temperature is defined following the keyword TEMPERAT.

Table 5.13.1: definition of initial temperature parameters

TEMPERATURE INITIALE (ncoucht)

ztemp vtemp jtemp

3.0e3 273. 2

-15.0e3 273 2

-15.0e3 673 2

-25.0e3 673 2

-50.0e3 673 2

xtgaus ytgaus ztgaus

10e3 0e3 -12.5e3

dxgaus dygaus dzgaus tgaus

10e3 10e3 10e3 0.0

nom du fichier de temperature (utilisé si nfilet = 1)

tempi

Table 5.13.2: definition of initial temperature in a user-defined file

100000 500000 -10000 x0 y0 z0

200000 200000 -20000 dx dy dz

2 2 2 (x,y,z)

673 673 673 673 673 673 673 673

5.14 Time stepping

Time stepping is an important aspect of the modelling as it directly influences the mass matrix computation used to compute the inertial forces. A requirement for a quasi-static solution is that these ine at each time step in the standard output in foreground, and in the file -file in background. These norms INERTD and INERT correspond to the ratio between the inertial forces norm on unconstrained nodes and the sum of reaction forces norms on nodes included in the boundary conditions:

INERT is computed using the total stress tensor in order to compute vector forces. Particularly the increase of mean stress with depth is represented in this inertia estimate. INERTD is computed using the deviatoric stress tensor in order to remove the effect of mean stress increase with depth. The use of the code that has been done in geodynamics (see Examples and reference list) has been done with values of INERTD between 0 and 2%. A sudden increase of these norms generally means that the mesh becomes too deformed due to localized deformation.

Figure 5.13.1 : Evolution of inertial norm INERTD with time

6. VISUALIZING AND PLOTTING

Using the model results consist in reading the nodal and elementary values written the p-file. A first application corresponds to 2D/3D plotting using the graphic software xadeli. This program is a graphic display of geometry, displacements, strain and stress for a specific time step. Also, this program can generate a script that can be further processed by GMT (public domain program Generic Mapping Tool) in order to provide a postscript file of the graphic display.

6.1 2D plotting

Type the command xadeli2d. A window normally appears. Go into the menu “file” and use the submenu “load” to read the file p-file. the other menu parameters is used to define “axes”, “time step” (the first one is 0), “plotted values”, etc. Note that using the “plotted values” submenu allows to display mesh, nodal vectors, elementary invariants, and principal tensor axis for elements. Use the “add zone” submenu to add other drawing zones below. Finally use the submenu “GMT script” to save the corresponding GMT script in the directory PLOT1. Note that this directory is created by the software and must NOT be present before this action.

6.3 3D plotting

3D plotting can be performed using the same principle as for 2D plotting. However, the procedure needs two steps. First, the program p2x reads the p-file and the file p2x_ini, projects the 3D mesh on a specified 2D plane, and write the result in a q-file. Then using the program xadeli3d allows to visualize and plot the results. Note that some options (such as principal tensor axes) cannot be used here.

The structure of the file p2x_ini is given in table 14. In order to build a q-file using p2x, the following parameters should be set up:

-iproj must be setup to 1;

-itrans can be 0 (no hidden edges) or 1 (hidden edges);

-azimuth defines the azimuth of the projection in degrees (geographical convention, North corresponds to y-axis);

-dip defines the dip angle of the projection with respect to horizontal (90 means a view from the top).

The other parameters are not used by the option iproj = 1. The option iproj = 2 is used to perform a profile of interpolated velocity or displacement at the surface of the 3D mesh. The surface where interpolation is made is defined by the line “nombre et faces”. The profile is defined with azimuth, dip, and the origin of the profile (line “longueur et pas de coupe”). The output is written is the files profil and profilxy. The option iproj = 3 is used to perform a computation of stress, strain or strain rate tensors on a profile within the mesh. The surface information is therefore nor used. In addition to the previous option, setup the value of itens to 1, 2 or 3. the lines defined by ivalp et ivalnum are used to choose the tensor type and the principal values to be plotted. The output is written is the files forage and foragexy. The option iproj = 4 allows the user to make horizontal cross sections through the model.

Table 6.3.1: Definition of parameters used by the program p2x

1 iproj ( 1=peau 2=profil 3=forage 4=coupe horiz. 5=coupe vert.)

0 itrans ( 1=projection des faces cachees )

1 itens (1=stress 2=strain 3=strain\_rate )

1 ivalp (1=vp 0=xyz)

1 ivalnum (1-2-3 ; 0 si pas interesse)

2 5 igmtx et igmty

0.7 rayon : rayon de recherche autour du point du profil en km

90 azimuth : azimuth geographique en degres)

0 dip: angle d'incidence, 0=rasant, 90=vertical

0 0 -25 point origine de la coupe

60 1 longueur et pas de coupe

3 3 7 12 nombre et faces

3 1 9 14 nombre et faces

6.4 2D-3D animation

The program aadeli can be used to generate a sequence of graphic files that can be processed by the public software ximage (NCSA X Image Version 1.2.1) developed at the University of Illinois. aadeli reads the initialization file aadeli.ini, the p-file and produces a sequence of files with the generic name m-file following by numerical extensions 0000, 0001, 0002, etc... The meaning of the most useful parameters in file aadeli.ini are the following:

alpha controls the interpolation between two time steps. 1 means no interpolation. The number of interpolated images is 1/alpha-1. Zero is forbidden;

npx and npy define the pixel size of the graphical image;

ilisse performs a linear interpolation of elementary values if 1 is given. 0 performs no interpolation;

imesh allows to draw the mesh if 1 is given;

rgb allows to write the output using a RGB format for video processing.

Appendix

Table A1: example of a i-file for the version 3p4 (2D)

IECHO

102

IMECA ITHERM

0 1

NFRON0 NFACE0 ICERC

12 30 0

NVFIX0 NPRES0 NGRAV

0 0 1

NTHER0 NFLUX0 nfilet

3 3 0

NFRIC NERO ISOMESH ISOSTRS INIVIT INITEMP IRHOREF

0 0 0 0 0 2 0

NTIME TFIN

10000 3.0E15

NELEM

3000

COEF

1000.

COORDONNEES DE NFRON0 POINTS DU FRONT

numP X Y

1 0 0

2 29. 0

3 31 0

4 60 0

5 0 -25

6 29 -25

7 31 -25

8 60 -25

9 0 -50

10 29 -50

11 31 -50

12 60 -50

CONTOURS DE MATERIAUX (nface0 valeurs)

numF point fleche position type

1 2 0 0 0

2 1 0 0 0

3 5 0 0 0

4 6 0 0 0

5 -999 0 0 0

6 2 0 0 0

7 6 0 0 0

8 7 0 0 0

9 3 0 0 0

10 -999 0 0 0

11 3 0 0 0

12 7 0 0 0

13 8 0 0 0

14 4 0 0 0

15 -999 0 0 0

16 6 0 0 0

17 5 0 0 0

18 9 0 0 0

19 10 0 0 0

20 -999 0 0 0

21 7 0 0 0

22 6 0 0 0

23 10 0 0 0

24 11 0 0 0

25 -999 0 0 0

26 8 0 0 0

27 7 0 0 0

28 11 0 0 0

29 12 0 0 0

30 -999 0 0 0

RHEOLOGIE DE CHAQUE MATERIAU :

numat ityp coef_mesh

young poiss rho0 nul nul

gamm0 eacti expos tkelv nul

cohes phini phfin kappac psi

seuil hard kappac nul nul

cspef condt srcth dilat nul

1 5 1.0

1.e11 0.25 2.8e3 0.0 0.0

0.84912E-14 110e3 1.0 948. 0.0

1.e6 15. 15. 0.01 0.0

1.e8 -0.1 1e-2 0.0 0.0

1070. 3.0 0.0e-6 0.e-5 0.0

2 5 1.0

1.e11 0.25 2.8e3 0.0 0.0

0.84912E-14 110e3 1.0 948. 0.0

1.e6 3. 3. 0.01 0.0

1.e6 0.0 1e-2 0.0 0.0

1070. 3.0 0.0e-6 0.e-5 0.0

3 5 1.0

1.e11 0.25 2.8e3 0.0 0.0

0.84912E-14 110e3 1.0 948. 0.0

1.e6 15. 15. 0.01 0.0

1.e8 -0.1 1e-2 0.0 0.0

1070. 3.0 0.0e-6 0.e-5 0.0

4 5 1.0

1.e11 0.25 3.3e3 0.0 0.0

0.48E-18 83e3 1.0 1623. 0.0

1.e6 15. 15. 0.01 0.0

1.e8 -0.1 1e-2 0.0 0.0

1070. 3.0 0.0e-6 0.e-5 0.0

5 5 1.0

1.e11 0.25 3.3e3 0.0 0.0

0.84912E-14 110e3 1.0 948. 0.0

1.e6 3. 3. 0.01 0.0

1.e6 0.0 1e-2 0.0 0.0

1070. 3.0 0.0e-6 0.e-5 0.0

6 5 1.0

1.e11 0.25 3.3e3 0.0 0.0

0.48E-18 83e3 1.0 1623. 0.0

1.e6 15. 15. 0.01 0.0

1.e8 -0.1 1e-2 0.0 0.0

1070. 3.0 0.0e-6 0.e-5 0.0

VITESSES IMPOSEES

numF codes_N-T vitesse_N vitesse_T evolution (o=1,n=0)

evolution des sollicitations (nbe de points nevolv)

valeurs normalisees du temps et de l'intensite, type de fonction

PRESSIONS IMPOSEES

numF type_de_pression pression_N pression_T evolution (o=1,n=0)

evolution des sollicitations (nbe de points nevolp)

valeurs normalisees du temps et de l'intensite, type de fonction

GRAVITE EN X ET EN Y (si ngrav > 0)

0 -10.0

evolution des sollicitations (nbe de points nevolg)

valeurs normalisees du temps et de l'intensite, type de fonction

TEMPERATURES IMPOSEES

numF valeur evolution

14 273 0

9 273 0

1 273 0

evolution des sollicitations (nbe de points nevolt)

valeurs normalisees du temps et de l'intensite, type de fonction

FLUX THERMIQUE IMPOSE

numF flux_normal evolution

18 -0.02 0

23 -0.02 0

28 -0.02 0

evolution des sollicitations (nbe de points nevolq)

valeurs normalisees du temps et de l'intensite, type de fonction

CONTACT (tolerance, nbr max d'iteration)

1.e-3 50

numero, nbe faces, liste des numeros des faces (nfric ligne)

nbr de visibilites mutuelles

zone_cand_num1 / zone_antg_num2, irevers, friction, ievol

1 2 1 0.3 0

3 5 1 0.6 0

4 5 1 0.6 0

evolution des coefficients de friction

numero de la visibilite

valeurs normalisees du temps et de l'intensite, type de fonction

EROSION (ibase,hbase0,transp,pluvio)

0 0.5e3 0.00 0e-8

XLF ICONLIM, DIFF0, PANTC, COENL

20000. 00 1.e-7 1.0 0.0

rhosedi qsed

2.3e3 3.17e-13

nombre et liste des faces (nfsur, lfsur(i),i=1,nfsur)

3 7 8 1

evolution des sollicitations (nbe de points nevolv)

valeurs normalisees du temps et de l'intensite, type de fonction

NIVEAU ET DENSITE DE L'EAU

y_water rho_water

0.0e3 1.0e3

PRECONTRAINTES LITHOSTATIQUES (nlitho, ylitho et rlitho par couche)

ylitho rlitho

0.5e3 0.0

0.0e3 1.0e3

-5.0e3 2.6e3

-6.0e3 2.6e3

-7.0e3 2.6e3

rhydro yrefw

2.2e3 -5.0e3

VITESSE INITIALE (xvit0 xvit1 vit0 vit1 radvit0 epvit0 angvit0)

0.0 10e-2 0. 0. 100e3 10e3 30

TEMPERATURE INITIALE (ncoucht)

ytemp vtemp jtemp

1 273. 2

-50E3 773. 2

xtgaus ytgaus dxgaus dygaus tgaus

100e3 -10e3 100e3 5e3 0

nom du fichier de temperature (utilisé si nfilet = 1)

tempi

FIN_DU_FICHIER

Table A2: example of a i-file for the version 3d4 (3D)

IECHO

5801

IMECA ITHERM irwmemo

1 0 0

NFRON0 NFACE0 nfond ICERC IMAP3

10 16 0 0 1

NVFIX0 NPRES0 NGRAV ISTRIKE nfilev

8 0 1 0 0

NTHER0 NFLUX0 nfilet

0 0 0

NFRIC NERO ISOMESH ISOSTRS INIVIT INITEMP nfilero iquake

2 0 0 1 0 1 0 0

NTIME COEFTIME TFIN

1000 3.1536e7 1000

COEFCOORD

1e3

COORDONNEES DE NFRON0 POINTS DU FRONT

numP X Y

1 0 0

2 50 0

3 100 0

4 50 20

5 50 150

6 50 150

7 50 280

8 0 300

9 50 300

10 100 300

CONTOURS DE MATERIAUX (nface0 valeurs)

numF point fleche position type

1 9 0 0 0

2 8 0 0 0

3 1 0 0 0

4 2 0 0 0

5 4 0 0 0

6 5 0 0 0

7 7 0 0 0

8 -999 0 0 0

9 9 0 0 0

10 7 0 0 0

11 6 0 0 0

12 4 0 0 0

13 2 0 0 0

14 3 0 0 0

15 10 0 0 0

16 -999 0 0 0

imapg npmapg

0 2

tmapgmin telem0 (echelle utilisateur)

2 8

dmapcrit (echelle utilisateur)

liste des npmapg points de controle (x,y, echelle utilisateur)

0 210

0 290

COUCHES

ncouch / ncouch+1 lignes / icc,zcouch(ic),kcouch(ic),coef

1 0.0 9 1

1 -25e3 1 1

ISOSTASIE

topom

FONDATIONS RIGIDES (5 x nfond lignes)

numF / 4 x (x,y,z)

RHEOLOGIE DE CHAQUE MATERIAU :

numat ityp coef_mesh

young poiss rho0 nul nul

gamm0 eacti expos tkelv nul

cohes phini phfin kappac psi

seuil hard kappac nul nul

cspef condt srcth dilat nul

1 2 1.0

1.e11 0.25 2.8e3 0.0 0.0

0.84912E-14 110e3 1.0 948. 0.0

1.e6 15. 15. 0.01 0.0

1.e8 -0.1 1e-2 0.0 0.0

1070. 3.0 0.0e-6 0.e-5 0.0

2 2 1.0

1.e11 0.25 2.8e3 0.0 0.0

0.84912E-14 110e3 1.0 948. 0.0

1.e6 15. 15. 0.01 0.0

1.e6 0.0 1e-2 0.0 0.0

1070. 3.0 0.0e-6 0.e-5 0.0

VITESSES IMPOSEES

coef. multiplicateur des vitesses imposees.

3.1709792e-11

numF codes_N-T1-T2 liaison v1_x v1_y v1_z v2_x v2_y v2_z evolution (o=1,n=0)

1 1 1 0 1 0 12 0 0 0 0 0

2 1 1 0 0 0 0 0 0 0 0 0

3 1 1 0 1 0 0 0 0 12 0 0

13 1 1 0 1 0 12 0 0 24 0 0

14 1 1 0 0 0 24 0 0 0 0 0

15 1 1 0 1 0 24 0 0 12 0 0

18 0 0 1 0 0 0 0 0 0 0 0

19 0 0 1 0 0 0 0 0 0 0 0

nom de fichier de conditions en vitesse (nfilev > 0)

fich5

evolution des sollicitations (nbe de points nevolv)

valeurs normalisees du temps et de l'intensite, type de fonction

PRESSIONS IMPOSEES

numF type_de_pression p_N p_T1 p_T2 evolution (o=1,n=0)

evolution des sollicitations (nbe de points nevolp)

valeurs normalisees du temps et de l'intensite, type de fonction

STRIKE-SLIP CONDITION (2 lignes (nfaces, liste))

1 5

1 6

GRAVITE EN Z.

-10.0

evolution des sollicitations (nbe de points nevolg)

valeurs normalisees du temps et de l'intensite, type de fonction

TEMPERATURES IMPOSEES

numF valeur evolution

evolution des sollicitations (nbe de points nevolt)

valeurs normalisees du temps et de l'intensite, type de fonction

FLUX THERMIQUE IMPOSE

numF flux_normal evolution

evolution des sollicitations (nbe de points nevolq)

valeurs normalisees du temps et de l'intensite, type de fonction

CONTACT (tolerance, nbr max d'iteration)

1.e-3 50

numero, nbe faces, liste des numeros des faces (nfric ligne)

1 2 5 6

2 2 10 11

nbr de visibilites mutuelles

zone_cand_num1 / zone_antg_num2, irevers, friction, ievol

1 2 0 0.10 0

SEISMES fricb dfric

0.10 0.03

timecos dtcos0 timepost ntpost

100 100 3.15e10 200

nom de fichier de deplacement cosismique max

depla

nombre, pas d'echantillonage et coordonnees des points de controle en surface

12 npqout

1 ispanq

0e3 150e3

10e3 150e3

20e3 150e3

30e3 150e3

40e3 150e3

49.9e3 150e3

51.1e3 150e3

60e3 150e3

70e3 150e3

80e3 150e3

90e3 150e3

100e3 150e3

EROSION (ibase,hbase0,transp,pluvio)

0 -1.0e3 0.1 1.e-8

XLF ICONLIM, DIFF0, PANTC, COENL

20000. 00 1.e-7 1.0 0.0

nombre et liste des faces (nfsur, lfsur(i),i=1,nfsur)

2 16 17

nom du fichier d'erosion (si nfilero = 1)

erode

evolution des sollicitations (nbe de points nevolv)

valeurs normalisees du temps et de l'intensite, type de fonction

NIVEAU ET DENSITE DE L'EAU

y_water rho_water

0.0e3 1.0e3

PRECONTRAINTES LITHOSTATIQUES (nlitho, ylitho et rlitho par couche)

ylitho rlitho

0.0e3 0.0

-25.0e3 2.8e3

rhydro

2.8e3

VITESSE INITIALE (xvit0 xvit1 vit0 vit1 )

0.0 10e-2 .5e-4 -0.5e-4

TEMPERATURE INITIALE (ncoucht)

ztemp vtemp jtemp

0.0 624. 2

-13.0e3 624. 2

-25.0e3 1083. 2

xtgaus ytgaus ztgaus

10e3 0e3 -12.5e3

dxgaus dygaus dzgaus tgaus

10e3 10e3 10e3 0.0

nom du fichier de temperature (utilisé si nfilet = 1)

tempi

FIN_DU_FICHIER