Mechanics Input (mech_3D.inp)
(January, 2022)
The file mech_3D.inp contains all flags, input and output file
names, and geomechanics properties for the geomechanical simulation.
File mech_3D.inp is free format file. As for TOUGH2, inputs are
preceeded by a keyword heading. Comments and blank lines placed
at the end of records are ignored, which allows for detailed user
descriptions of input parameters.
Elastic properties may be specified for a general horizontally isotropic
(vertically anisotropic) elastic medium. However, equations for failure
strains assume isotropic material, and any failing material is treated
as if it had isotropic elastic properties prior to failure (given by Ep
and nu_pp, below). For isotropic material, Ez = Ep, nu_zp = nu_pp, and
the Gzp should be Ep / 2 (1+nu_pp) to be consistent with them.
order
record 1: 'x' 'z' 'y'
parameters
record 1: nx ny nz Number of mesh blocks in x, y and z directions.
record 2: ftol stol
ftol - Relative tolerance for solution of mechanics force
balance equations (ftol), e.g., 1e-6 .
stol - Maximum change in porosity (in one time step) before
reducing time step (stol), e.g., 5.e-3. If failure is
not expected, stol can be set to a fairly small value
(1.e-4) and used to control time stepping. If hydraulic
fracturing is expected, a fairly large value (e.g., 5.e-3)
should be used.
record 3: initstressfile
File name for initial stresses. If the file name is given
as "smoothedlocalloading" or "uniformstressgradient",
no file is read, and the initial stresses are calculated
internally. These two cases take additional inputs;
record 3a (smoothedlocalloading): initstressratio(1) initstressratio(2) initstressratio(3) horzstressrotation
Ratios of sigma_x'x', sigma_y'y', and sigma_zz to computed
vertical stress, and rotation of principal stress coordinates
x' and y' (rotation of x' from x towards y, in degrees).
Vertical stress is computed from local loading by rock mass
computed from grain density (flow.inp rocks.1 record,
characters 11-20), porosity, and internally computed fluid
density (at P, T conditions), and from vertical stress
specified at boundary normal to z. Assumes that vertical
(zz) stress is set on one surface normal to z.
Loading stress at each level is smoothed before integration
to next level. When smoothed local loading stresses are
computed, they are subsequently written to a file named
smoothedlocalloading.
record 3a (uniformstressgradient): refstress(1) refstress(2) refstress(3)
record 3b (uniformstressgradient): str_grad(1,1) str_grad(1,2) str_grad(1,3)
record 3c (uniformstressgradient): str_grad(2,1) str_grad(2,2) str_grad(2,3)
record 3d (uniformstressgradient): str_grad(3,1) str_grad(3,2) str_grad(3,3)
Diagonal stresses sigma_xx, sigma_yy, sigma_zz (Pascal) at (0,0,0)
and gradient of diagonal stresses (Pa/m). Indexing is such that
(str_grad(1,1), str_grad(1,2), str_grad(1,3)) =
d (sigma_xx, sigma_yy, sigma_zz)/ dx. Compressive stresses are
considered negative.
record 4: iBC(1,1) iBC(2,1)
Mechanics boundary conditions at x(1), x(nx+1). Choices are
0 no displacement normal to boundary
1 no displacement on boundary
2 normal traction on boundary constant to be read in
3 traction on boundary constant value to be read in
4 normal traction on boundary held at initial normal stress in adjacent cell
5 traction on boundary held at initial stress in adjacent cell
record 4a (iBC(1,1)= 2 case): surftractn(1,1,1)
x traction on x(1) surface.
record 4a (iBC(1,1)= 3 case): surftractn(1,1,1) surftractn(2,1,1) surftractn(3,1,1)
x,y,z tractions on x(1) surface.
Record 4a present when iBC(1,1)= 2 or 3.
record 4b (iBC(2,1)= 2 case): surftractn(1,2,1)
x traction on x(nx+1) surface.
record 4b (iBC(2,1)= 3 case): surftractn(1,2,1) surftractn(2,2,1) surftractn(3,2,1)
x,y,z tractions on x(nx+1) surface.
Record 4b present when iBC(2,1)= 2 or 3.
record 5: iBC(1,2) iBC(2,2)
Mechanics boundary conditions at y(1), y(ny+1).
Choices 0, 2, 3, 4, or 5, as above.
record 5a (iBC(1,2)= 2 case): surftractn(2,1,2)
y traction on y(1) surface.
record 5a (iBC(1,2)= 3 case): surftractn(1,1,2) surftractn(2,1,2) surftractn(3,1,2)
x,y,z tractions on y(1) surface.
Record 5a present when iBC(1,2)= 2 or 3.
record 5b (iBC(2,2)= 2 case): surftractn(2,2,2)
y traction on y(ny+1) surface.
record 5b (iBC(2,2)= 3 case): surftractn(1,2,2) surftractn(2,2,2) surftractn(3,2,2)
x,y,z tractions on y(ny+1) surface.
Record 5b present when iBC(2,2)= 2 or 3.
record 6: iBC(1,3) iBC(2,3)
Mechanics boundary conditions z(1), z(nz+1).
Choices 0, 2, 3, 4, or 5, as above.
record 6a (iBC(1,3)= 2 case): surftractn(3,1,3)
z traction on z(1) surface.
record 6a (iBC(1,3)= 3 case): surftractn(1,1,3) surftractn(2,1,3) surftractn(3,1,3)
x,y,z tractions on z(1) surface.
Record 6a present when iBC(1,3)= 2 or 3.
record 6b (iBC(2,3)= 2 case): surftractn(3,2,3)
z traction on z(nz+1) surface.
record 6b (iBC(2,3)= 3 case): surftractn(1,2,3) surftractn(2,2,3) surftractn(3,2,3)
x,y,z tractions on z(nz+1) surface.
Record 6b present when iBC(2,3)= 2 or 3.
record 7: 0 0
record 8: 0 0
record 9: 0 0
record 10: BCvol_thresh
Threshold for zeroing mechanics porosity changes (e.g.,
1.e50): elements with (MESH file) ELEME field volume greater
than or equal to BCvol_thresh are not given porosity changes
due to mechanics. (Zeroing of porosity changes due to
chemistry depends separately on whether ELEME volume is
greater than 1020 m3.)
petsc_mnp (optional field)
record 1: m n p
number of processes to assign in x, y and z directions.
When program is compiled for MPI parallel processing and
run under mpirun, used to override program’s default domain
sub-division, otherwise ignored. When not ignored, product
m*n*p must match number of MPI processes specified for mpirun.
materials
record 1: nummat
number of materials, i.e., number of rock types. Must be
the same number of rock types as in flow.inp. Mechanical
properties of each material are given on nummaterlines
consecutive lines.
record 2: materlinecontrols
three 0 or 1's. The first controls whether (1) or not (0)
a line of temperature-stress bilinear strain parameters
follows the first four lines of material properties for each
rock type. The second controls whether or not a line of
porosity dependent moduli parameters follows first four
(or five) lines of material properties for each rock type.
The third controls whether or not a line of volumetric
effective stress permeability parameters is read. The number
of lines of material properties for each rock type
nummaterlines is either 4, 5, 6, or 7, depending on the number
of 1's in materlinecontrols (000, 100, 010, 110, 001, 101,
011, or 111).
records 3, 3+ nummaterlines , ..., 3 + i*nummaterlines, ...: mname young_par young_perp shear poisson_par poisson_perp angle
five character material name (rock type) as in flow.inp ROCKS block,
Youngs modulus parallel; Ep (in plane normal to z (anisotropy) direction) (Pascal),
Youngs modulus perpendicular; Ez (in anisotropy direction) (Pascal),
shear modulus; Gzp for shearing between parallel and z directions,
Poisson ratio parallel; nu_pp between x and y directions,
Poisson ratio perpendicular; nu_zp ratio of (-)parallel
displacement to vertical displacement for an applied vertical stress,
(angle of anisotropy from vertical); not used, vertical anisotropic assumed.
records 4, 4+ nummaterlines , ..., 4+i x nummaterlines, ...: biot cohesion tensile_stren MC_friction bulk_dens coeff_plast dilation_ang thermal_dilat
Biot coefficient, e.g., 1 - K_dr/K_grain . (Value = -1 treated specially, for water elements in which thermal stresses from rock matrix expansion are not present.)
Mohr-Coulomb cohesion (Pascal), e.g., 0,
tensile strength (Pascal), e.g., 0,
Mohr-Coulomb friction angle (degrees), e.g., 30,
bulk density (kg/m3) (either this, or flow.inp grain density and porosity is used for mechanics, according to idensitysource below)
coefficient of plastic porosity- not used,
dilation angle (dilation on M.C. failure) (degrees)
thermal dilation coefficient (1/ ?C).
records 5, 5 + nummaterlines , ..., 5 + i*nummaterlines, ...: imatrixpermporotype amatrixpermparam bmatrixpermparam cmatrixpermparam ifracpermporotype afracpermparam bfracpermparam
index of matrix porosity permeability relationship for rock type (0-6).
Types 1-6 as in TOUGHREACT V3.0-OMP. Index 0 defaults to type
specified in chemical.inp Zppr record ippzon value,
when chemistry is run, and to none otherwise.
a_param of matrix porosity permeability relation.
b_param of matrix porosity permeability relation.
c_param of matrix porosity permeability relation.
index of fracture porosity permeability relationship for rock type (0,
7 or 8). 0 neglects permeability changes due to failure porosity.
a_param of fracture porosity permeability relation,
for ifracpermporotype=7, fraction of initial porosity assumed
due to fractures (generally a small number),
for ifracpermporotype=8, assumed fracture spacing (m).
b_param of fracture porosity permeability relation,
for ifracpermporotype=7 or 8, fraction of initial
permeability assumed due to fractures (generally close to 1).
records 6, 6+ records 6, 6+ nummaterlines , ..., 6 + i*nummaterlines, ...: idensitysource idissolv_poro
idensitysource =0 for bulk density from flow.inp grain density
and porosity and initial fluid density;
=1 for bulk density from here;
=2 for bulk density from here and grain density
reset from bulk density and flow.inp porosity,
assuming a fluid density of 1000 kg/m3 ;
=3 for bulk density from here and grain density
reset from bulk density and flow.inp porosity,
assuming a fluid density of 1000 kg/m3
idissolv_poro = 0, no imposed dissolution strains,
= 1, sigma_eff*(T-T_0 )*coef , stress temperature
product strains imposed. (idisolv_poro=1 may
need quite small time steps for stability.)
records 7, 7+ nummaterlines , ..., 6 + i*nummaterlines, ... (when materlinecontrl(1:1)=1): c_bilin T_bilin K_bilin d_bilin
bilinear strain parameters (1-4) for this rock type.
records 8, 8+ nummaterlines , ..., 7 + i*nummaterlines, … (when materlinecontrl(2:2)=1): phi_crit a_k b_k c_k d_k a_sh b_sh c_sh d_sh
porosity dependent moduli parameters for this rock type.
Bulk modulus is given by
K * ( a_k +b_k (phi -phi_crit )^(1/2)
+c_k (phi -phi_crit )+d_k (phi -phi_crit )^2 )
/ ( a_k+b_k (phi_0 -phi_crit )^(1/2)
+c_k (phi -phi_crit ) + d_k (phi_0 -phi_crit ))
and shear moduls by
mu * ( a_sh +b_sh (phi -phi_crit )^(1/2)
+c_sh (phi -phi_crit ) +d_sh (phi-phi_crit )^2 )
/ (a_sh+b_sh (phi_0- phi_crit )^(1/2)
+c_sh (phi -phi_crit ) + d_sh (phi_0- phi_crit )^2 ),
where K and mu are the starting bulk and shear moduli.
records 9, 9+ nummaterlines , ..., 8 +i*nummaterlines, ... (when materlinecontrl(3:3)=1): log10(sigma_k0) dlog10(sigma_efftrans)/dlog10(T) dlog10(k)/dlog10(sigma_vol)|low dlog10(k)/dlog10(sigma_v)|high
volumetric effective stress dependent permeability parameters,
after Watanabe, et al. (2017). Transistion volumetric effective
stress is eigma_efftrans = log10(sigma_k0) + log10(T) * dlog10(sigma_efftrans)/dlog10(T).
For sigma_effvol <= sigma_efftrans have
delta log10(k) = dlog10(k)/dlog10(sigma_vol)|low & delta log10(sigma_effvol).
For sigma_effvol > sigma_efftrans have
delta log10(k) = dlog10(k)/dlog10(sigma_vol)|high delta log10(sigma_effvol).
Watanabe et al.'s values for these parameters are 12.8, -4.3, -0.9, and -3.8.
Here, T is temperature Celsius, k is permeability in m^2, and,
with respect to these parameters, volumetric stressess are in MPa,
with compressive stress positive, both contrary to the rest of
this file. When using this option, for materials for which it
is not to be applied, set values to 0 0 0 0. Generally, when
using this option, one wants imatrixpermporotype=0 for the
materials it is applied to.
bilinearstrain (optional field)
record 1: c_bilin T_bilin K_bilin d_bilin tau_smoo
For use with rock types with idisolv_poro=1 (above).
Sets parameters for temperature -volumetric effective stress
product dependent strain of form:
epsilon_d = sigma_eff [ ( T - T_bilin) c_bilin - 1/(3 K_bilin) ] - epsilon_d0
where
epsilon_d0 = sigma_eff(0) [ ( T(0) - T_bilin) c_bilin - 1/(3 K_bilin) ]
when indicator trace(sigma_eff ) == (T-T_bilin) is increasing,
that is, when volumetric effective stress magnitude trace (sigma_eff)
or temperature T increases. When this indicator decreases,
strain epsilon_d reverses, but with only d_bilin as much change
as indicated above (0 <= d_bilin <= 1). These strains are
smoothed over time scale ?_smoo for stability.
When materials field nmaterlines < 5, all five bilinearstrain
entries here are used. When materials nmaterlines =5, only fifth
bilinearstrain parameter (tau_smoo) is used from here, other
values are taken from fifth line of material properties for
each rock type. bilinearstrain parameters have no effect on rock
types with idissolv_poro = 0.
gridmaterials (optional field)
record 1: ngridmater
number (<=14) of mechanics properties to take from
geomechanprop.inp file instead of from mech_3D.inp
materials block.
record 2: igridmater_1 igridmater_2 ... igridmater_ ngridmater
indices (between 1 and 14) of material properties to take
from geomechanprop.inp instead of from mech_3D.inp materials
block. (igridmater=1: young_par, igridmater=14: thermal_dilat.)
(Properties 15 to 22 are always set by rock type, and
properties 23-26 either by rock type or by the bilinearstrain
field parameters.)
mopm
record 1: mopmech(1) mopmech(2) mopmech(3) mopmech(4) mopmech5) numMINCMe iverbose ForPoroC chm_rel_p chm_rel_T maxiter_me ibiot_eff_stress ifailureporo poromix pelastmix
mopmech(1) = 0 no mechanics; = 1 mechanics, =2 read mechanics input files and stop.
mopmech(2) -- not used
mopmech(3) -- not used
mopmech(4) = 0 write total stress, biot*pressure, temp, total
strain only at elements requested by elgm (below).
(RMECH.out files produced with mopmech(4)=0
are insufficient for use as an initial stresses file.)
= 1 write total stress, biot*pressure, temp, total strain,
average displacement at every element, to file RMECH.out
every time step
= 2 write total stress, biot*pressure, temp, total strain,
average displacement at every element, to file RMECH.out
when save files are updated
= 3 write total stress, biot*pressure, temp, total strain
at every element to file RMECH.out each time step,
overwriting old file
= 4 write total stress, biot*pressure, temp, total strain
at every element to file RMECH.out after first time
step and stop
= 5 write initstressratio(1:3) * smoothed local loading stresses
to file smoothedlocalloading and stop
Cases 0, 3, 4, & 5 are the most practical, as RMECH.out files from them
either grow slowly (case 0), or remain a constant size.
When ibiot_eff_stress = 0 (below), pressure is written to file RMECH.out
in place of biot*pressure.
mopmech(5) not used
numMINCMe number of MINC levels (0 no MINC)
iverbose = 0 quiet
> 0 diagnostic information written to standard output
mopmech(8) =0, default,
=1, on restarts, zero displacements, strains, failure strains,
and dissolution strains, (leaves stresses unchanged)
=2, on restarts, same as 1, plus, restores any cohesion
lost on failure
chm_rel_p -- not used
chm_rel_T -- not used
maxiter_me -- not used
ibiot_eff_stress = 0 use unweighted pressure in effective stress for
Mohr-Coulomb and tensile failure criteria (following
Dempsey et al., 2014, and Jaeger, et al., 2007),
= 1 weight pressure by Biot coefficient in computing
effective stress for Mohr-Coulomb and tensile failure
criteria (Zoback, 2014, undated lecture notes).
ifailureporo = 0 neglect changes in porosity due to material
failure when computing pressure,
= 1 compute pressure taking into account changes
in porosity due to material failure.
poromix fraction of change in failure porosity (e.g., 0.25) to
incorporate immediately in pressure calculations.
poromix < 1 increases stability under hydraulic
fracturing conditions. Failure porosity changes
get averaged in in 1/poromix time steps (in a
`leaky bucket' auto-regressive average).
When hydraulic fracturing is not expected to occur,
poromix = 1 is fine.
pelastmix fraction of change in poroelastic porosity (e.g., 0.25)
to incorporate immediately in pressure calculations.
pelastmix < 1 increases stability under hydraulic
fracturing conditions. Poroelastic porosity changes
get averaged in in 1/pelastmix time steps (in a `leaky bucket'
auto-regressive average). When hydraulic fracturing is not
expected to occur, pelastmix = 1 is fine.
elgm (optional field)
record 1: numprelem
number of mechanics elements to write Monitor_CNNNN files for,
where CNNNN corresponds to flow.inp element name.
record 2: monelem_mechname (1)
ELEME name (in MESH file) of element
to write mechanics monitor information,
...
record 1+numprelem: monelem_mechname (numprelem)
ELEME name (in MESH file) of element
to write mechanics monitor information,
permlimit (optional field)
record 1: ixpermlimit iypermlimit izpermlimit
indices of element to center perm limit function on.
record 2: radiuspermlimit permlim_atradius
radius (m) of uniform maximum permeability limit zone.
Permeability limit increases as r for r > radiuspermlimit.
maximum allowed permeability (m2) in uniform maximum permeability zone.
record 3; nexemptrocktype
number of rock types exempt from permeability limit.
record 4; rtype
name of first rock type to be exempt
(not present for nexemptrocktype=0).
...
record 3+nexemptrocktype; rtype
name of last rock type to be exempt from permeability limit.
vertaverage (optional field)
record 1; iz1vertavg iz2vertavg
z indices of elements to average stresses and strains between
for output file stress_strain_vertavg.tec
done (optional field)
end of mech_3D.inp file.
Mechanics initial stresses file.
The geomechanics initial stresses file name is read from mech_3D.inp
file parameters field record 3. The initial stresses file has
free format. The default format of the file is;
Title
time
1 sigma_xx (1,1) sigma_xx (2,1) sigma_xx (3,1) sigma_xx (4,1) sigma_xx (5,1) sigma_xx (6,1) sigma_xx (7,1) sigma_xx (8,1)
sigma_yy (1,1) sigma_yy (2,1) sigma_yy (3,1) sigma_yy (4,1) sigma_yy (5,1) sigma_yy (6,1) sigma_yy (7,1) sigma_yy (8,1)
sigma_zz (1,1) sigma_zz (2,1) sigma_zz (3,1) sigma_zz (4,1) sigma_zz (5,1) sigma_zz (6,1) sigma_zz (7,1) sigma_zz (8,1)
sigma_yz (1,1) sigma_yz (2,1) sigma_yz (3,1) sigma_yz (4,1) sigma_yz (5,1) sigma_yz (6,1) sigma_yz (7,1) sigma_yz (8,1)
sigma_xz (1,1) sigma_xz (2,1) sigma_xz (3,1) sigma_xz (4,1) sigma_xz (5,1) sigma_xz (6,1) sigma_xz (7,1) sigma_xz (8,1)
sigma_xy (1,1) sigma_xy (2,1) sigma_xy (3,1) sigma_xy (4,1) sigma_xy (5,1) sigma_xy (6,1) sigma_xy (7,1) sigma_xy (8,1)
2 sigma_xx (1,2) sigma_xx (2,2) sigma_xx (3,2) sigma_xx (4,2) sigma_xx (5,2) sigma_xx (6,2) sigma_xx (7,2) sigma_xx (8,2)
sigma_yy (1,2) sigma_yy (2,2) sigma_yy (3,2) sigma_yy (4,2) sigma_yy (5,2) sigma_yy (6,2) sigma_yy (7,2) sigma_yy (8,2)
sigma_zz (1,2) sigma_zz (2,2) sigma_zz (3,2) sigma_zz (4,2) sigma_zz (5,2) sigma_zz (6,2) sigma_zz (7,2) sigma_zz (8,2)
sigma_yz (1,2) sigma_yz (2,2) sigma_yz (3,2) sigma_yz (4,2) sigma_yz (5,2) sigma_yz (6,2) sigma_yz (7,2) sigma_yz (8,2)
sigma_xz (1,2) sigma_xz (2,2) sigma_xz (3,2) sigma_xz (4,2) sigma_xz (5,2) sigma_xz (6,2) sigma_xz (7,2) sigma_xz (8,2)
sigma_xy (1,2) sigma_xy (2,2) sigma_xy (3,2) sigma_xy (4,2) sigma_xy (5,2) sigma_xy (6,2) sigma_xy (7,2) sigma_xy (8,2)
...
nelem sigma_xx (1,nelem) sigma_xx (2,nelem) sigma_xx (3,nelem) sigma_xx (4,nelem) sigma_xx (5,nelem) sigma_xx (6,nelem)
sigma_xx (7,nelem) sigma_xx (8,nelem) sigma_yy (1,nelem) sigma_yy (2,nelem) sigma_yy (3,nelem) sigma_yy (4,nelem)
sigma_yy (5,nelem) sigma_yy (6,nelem) sigma_yy (7,nelem) sigma_yy (8,nelem) sigma_zz (1,nelem) sigma_zz (2,nelem)
sigma_zz (3,nelem) sigma_zz (4,nelem) sigma_zz (5,nelem) sigma_zz (6,nelem) sigma_zz (7,nelem) sigma_zz (8,nelem)
sigma_yz (1,nelem) sigma_yz (2,nelem) sigma_yz (3,nelem) sigma_yz (4,nelem) sigma_yz (5,nelem) sigma_yz (6,nelem)
sigma_yz (7,nelem) sigma_yz (8,nelem) sigma_xz (1,nelem) sigma_xz (2,nelem) sigma_xz (3,nelem) sigma_xz (4,nelem)
sigma_xz (5,nelem) sigma_xz (6,nelem) sigma_xz (7,nelem) sigma_xz (8,nelem) sigma_xy (1,nelem) sigma_xy (2,nelem)
sigma_xy (3,nelem) sigma_xy (4,nelem) sigma_xy (5,nelem) sigma_xy (6,nelem) sigma_xy (7,nelem) sigma_xy (8,nelem)
where each stress component is specified at the eight Gauss points in each element.
There are nelem lines of (48) stresses, and nelem = nx * ny * nz. Time is in seconds and stresses are in Pascals.
Extra items at the ends of lines are ignored. Files with title Geomechanics Restart are treated specially, assuming
all information written to a geomechan.restart file is present.
For an element with corners at (x(i),y(j),z(k)), (x(i+1),y(j),z(k)), (x(i),y(j),z(k+1)), (x(i+1),y(j),z(k+1)),
(x(i),y(j+1),z(k)), (x(i+1),y(j+1),z(k)), (x(i),y(j+1),z(k+1)), (x(i+1),y(j+1),z(k+1)), Gauss points are
1/[2 sqrt(3)] grid spacings out, in each dimension, from the cell center, towards corners (x(i),y(j),z(k)),
(x(i+1),y(j),z(k)), (x(i+1),y(j+1),z(k)), (x(i),y(j+1),z(k)), (x(i),y(j),z(k+1)), (x(i+1),y(j),z(k+1)),
(x(i+1),y(j+1),z(k+1)), (x(i),y(j+1),z(k+1)), respectively.
Stresses at elements are in MESH file element order.
Initial stresses file names smoothedlocalloading and uniformstressgradient in mech_3D.inp are treated specially
and do not require files of those names. Instead, they invoke reading additional lines of input from the
mech_3D.inp file and either calculate stresses internally from local loading stresses, or impose a uniform gradient
diagonal stress field, respectively. More details are given in mech_3D.inp file format above.