FLAC3D Theory and Background • Fluid-Mechanical Interaction

Properties and Units for Fluid-Flow Analysis

The properties that relate to fluid flow in FLAC3D are:

the mobility coefficient, \(k\), or the hydraulic-conductivity;

the fluid mass density, \({\rho}_f\);

the Biot coefficient, \(\alpha\), and Biot modulus, \(M\) (for flow through a material with compressible grains), or the fluid bulk modulus, \(K_f\), and porosity, \(n\) (for flow through a material with incompressible grains only), or the consolidation coefficient.

Properties used to specify fluid-mechanical interaction behavior are density-saturated, effective-cutoff, and pore-pressure-generation.

The thermal-coupling parameters are the coefficient of drained linear thermal expansion, \({\alpha}_t\) (thermal-mechanical), and the undrained thermal coefficient, \(\beta\) (poro-thermal).

These properties are examined below.

All thermal-poro-mechanical quantities must be given in a consistent set of units. No conversions are performed by the program.

Mobility Coefficient

The mobility coefficient, \(k\) (e.g., m²/(Pa-sec) in SI units) is used in FLAC3D. This term is used to avoid confusion about different ways the term \(permeability\) can be interpreted in different fields. It is the coefficient of the pressure term in Darcy’s law and is related to the \(hydraulic\) \(conductivity\), \(k_h\) (e.g., m/s), by the expression

(1)\[k = {k_h \over {{\rho_f} g}}\]

where \(g\) is the gravitational acceleration.

The \(intrinsic\) \(permeability\), \(\kappa\) (e.g., in m²), is related to \(k\):

(2)\[k = {\kappa \over \mu}\]

where \(\mu\) is the dynamic viscosity (e.g., units of N-sec/m²).

Equation (1) or (2) may be used to derive \(k\) (required by FLAC3D) from either \(k_h\), in velocity units, or \(\kappa\), in [length]² units (remember that \(k\) must end up with units [L³T/M] – e.g., in SI units this would be m³ sec/kg (or m²/Pa-sec))). Using the values \(\mu\) = 1.011 × 10^-3 kg/(m-sec) for water at 20° C, \(\rho_w g\) = 9.81 × 10³ Pa/m and 1 Darcy = 10^-8 cm², conversions may be derived to calculate \(k\) in SI units for water in FLAC3D:

\(k\) (in SI units) \(\equiv \kappa\) (in cm²) × 9.9 × 10^-2

\(k\) (in SI units) \(\equiv k_h\) (in cm/sec) × 1.02 × 10^-6

\(k\) (in SI units) \(\equiv\) permeability in millidarcies × 9.9 × 10^-13

If there is a variation of mobility-coefficient across the grid, the timestep will be controlled by the largest value (see this equation). For problems in which steady state (but not transient behavior) is required, it may be beneficial to limit the variations in mobility-coefficient to improve convergence speed. For example, there will probably be little difference in the final state between systems where there is a 20:1 variation in mobility, compared to a 200:1 variation, when the flow is moving entirely from one region to the other.

The mobility-coefficient (scalar or tensor) is a zone property specified using the zone fluid property command (see Grid Configured for Fluid Flow).

Mass Density

Three different mass densities may be given as input to FLAC3D in different circumstances: the dry density of the solid matrix, \(\rho_d\); the saturated density of the solid matrix, \(\rho_s\); and the density of the fluid, \(\rho_f\). Note that densities are only required if gravitational loading is specified.

If FLAC3D is configured for fluid flow (model configure fluid-flow), then the dry density of the solid material must be used. FLAC3D will compute the saturated density of each element using the known density of the fluid, the porosity, \(n\), and the saturation, \(s\): \(\rho_s = \rho_d + n s \rho_f\).

The only case when the saturated density is given as input is for an effective stress calculation (static pore-pressure distribution) not carried out in model configure fluid-flow mode. The zone gridpoint pore-pressure command specifies the location of the water table. The dry density is specified for zones above the water table, and the saturated density for zones below.

The solid density (dry or saturated) is given using the zone property density command. The fluid density can be imposed globally using the zone fluid density command, or it can be allowed to vary with position by using the zone fluid property fluid-density command. However, this capability of assigning spatial variation of fluid density must only be used with extreme caution because fluid mass density is not advected in FLAC3D. All densities are zone variables in FLAC3D, and are mass densities (e.g., kg/m³ in SI units).

Fluid Moduli

Biot Coefficient and Biot Modulus

The Biot coefficient, \(\alpha\), is defined as the ratio of the fluid volume gained (or lost) in a material element to the volume change of that element when the pore pressure is changed. It can be determined in the same drained test as that used to determine the drained bulk modulus, \(K\), of the material. Its range of variation is between \(3n/(2+n)\) and 1, where \(n\) is the porosity. In the particular case of an incompressible solid constituent, \(\alpha\) = 1. This value is the default value adopted by FLAC3D.

For an ideal porous material, the Biot coefficient is related to the bulk modulus of the solid component \(K_s\):

(3)\[\alpha = 1 - {{K}\over{K_s}}\]

The Biot modulus, \(M\), is defined as

(4)\[M = {{K_u - K}\over{\alpha^2}}\]

where \(K_u\) is the undrained bulk modulus of the material.

For an ideal porous material, the Biot modulus is related to the fluid bulk modulus, \(K_f\):

(5)\[M = {{K_f}\over{n + (\alpha - n)(1 - \alpha)K_f/K}}\]

where \(n\) is the porosity. Thus, for an incompressible solid constituent (\(\alpha\) = 1),

(6)\[M = K_f / n\]

The Biot coefficient and Biot modulus are zone properties specified using the zone fluid property biot-coefficient and zone fluid property biot-modulus command.

Fluid Bulk Modulus

In analyses where the grain compressibility can be neglected, the user has the choice to either use the default value of Biot coefficient (i.e., \(\alpha\) = 1) and assign a value equal to \(K_f/n\) to Biot modulus, or give the fluid bulk modulus \(K_f\) as input.

The bulk modulus, \(K_f\), is defined as

(7)\[K_f={\Delta P\over\Delta V_f/V_f}\]

where \(\Delta P\) is the change in pressure for a volumetric strain of \(\Delta V_f/V_f\).

When the fluid modulus \(K_f\) is given as an input, the Biot modulus is computed internally using Equation (7) for incompressible grains. In this calculation, the porosity (a zone property) is evaluated at the nodes using nodal volume averaging. he Biot coefficient is then set to 1 throughout the flow domain, irrespective of any value given for that property. The fluid calculation for incompressible grains is the default calculation in FLAC3D. The fluid modulus is a gridpoint variable specified using the zone fluid property fluid-modulus command.

The “compressibility” of the fluid \(C_f\) is the reciprocal of \(K_f\) (i.e., \(C_f = 1/K_f\)). For example, for pure water at room temperature, \(K_f\) = 2 × 10⁹ Pa, in SI units. In real soils, pore water may contain some dissolved air or air bubbles, which substantially reduce its apparent bulk modulus. For example, Chaney (1978) reports that the fluid modulus can be reduced by an order of magnitude for an air/water mixture at 99% saturation in compacted sand. For groundwater problems, the bulk modulus of water may be different in different parts of the grid, to account for the varying amounts of air present. A FISH function may be used to change local moduli according to some law (e.g., the modulus may be made proportional to pressure to represent a gas), but the user should be careful not to make nonphysical assumptions.

Consolidation Coefficient

In some cases the property you have available is a consolidation coefficient from a test of the site. FLAC3D allows the user to input the consolidation coefficient as a property to specify the bulk modulus. The conversion from consolidation-coefficient to Biot modulus is handled by the equation:

(8)\[c = {{k} \over {{1 \over M} + {\alpha^2 \over {K+ 4 G/3}} }}\]

where \(M\) is the Biot modulus, \(c\) is the given consolidation coefficient, \(k\) is the current isotropic mobility coefficient of the material, \(\alpha\) is the biot coefficient, \(K\) is the bulk modulus of the material, and \(G\) is the shear modulus of the material. If the Biot coefficient is 1.0, then this is equivalent to a fluid modulus described by

\[c = {{k} \over {{n \over K_f} + {1 \over {K+ 4 G/3}} }}\]

where \(K_f\) is the fluid modulus and \(n\) is the porosity of the material.

Note that all three forms of fluid moduli are kept consistent with each other when one is set, unless the Biot coefficient is not 1.0 in which case the returned fluid modulus will be 0.0. During actual calculations the Biot modulus is what is used.

Fluid Moduli and Convergence

If steady-state, fully saturated flow is required, the modulus \(M\) (or \(K_f\)) is unimportant to the numerical convergence process, because the response time of the system and the timestep are both inversely proportional to \(M\) (or \(K_f\)): the same number of steps are necessary, independent of \(M\) (or \(K_f\)). For systems containing a phreatic surface, however, a low bulk modulus will speed convergence to steady state, because the calculation for saturation change involves \(\Delta t\), not the product \(M \Delta t\) (or \(K_f \Delta t\)) (see this equation and this equation). Systems in which solid/fluid interaction is important are more complicated to assess (some guidelines may be found in Time Scale). However, a high value of \(M\) (or \(K_f\)) compared to the mechanical \(K\) will lead to slowly converging solutions. In any case, from a numerical point of view, it is not necessary to use values of \(M\) (or \(K_f\)) that are larger than 20 times \((K + 4/3 G)/\alpha^2\) (or \((K + 4/3 G)n\)) in the simulation (see Time Scale, this equation, and 1D Consolidation for an example).

The zone fluid modulus-automatic can be used to set or limit the fluid bulk modulus to a given ratio (the default is 20).

Note that if you are using the zone fluid steady-state command, it is not necessary to specify a fluid bulk modulus.

Fluid Moduli for Drained and Undrained Analyses

In FLAC3D, whenever the fluid bulk modulus (or Biot modulus) is selected and model configure fluid-flow is specified, the drained bulk modulus must be specified for the solid matrix. The apparent (undrained) bulk modulus of the solid matrix will then be computed, and will evolve with time.

An undrained analysis may also be performed without specifying model configure fluid-flow. In this case, the undrained bulk modulus \(K_u\) for the solid matrix should be specified. The undrained bulk modulus is

(9)\[K_u = K + {\alpha}^2 M\]

For an incompressible solid constituent: \(\alpha\) = 1, \(M = K_f/n\), and this formula becomes

(10)\[K_u = K + {{K_f}\over n}\]

Porosity

Porosity, \(n\), is a dimensionless number defined as the ratio of void volume to total volume of an element. It is related to the void ratio, \(e\), by

(11)\[n = {e\over 1+e}\]

The default value of \(n\), if not specified, is 0.5. \(n\) should be given as a positive number between 0 and 1, but small values (e.g., less than 0.2) should be used with great caution because the apparent stiffness of the pore fluid is proportional to \(K_f/n\). For low values of \(n\), the stiffness may become very large in comparison to the stiffness of the solid material, causing the FLAC3D solution to take a very long time to converge. Consider reducing \(K_f\), in this case (see Fluid Moduli and Convergence and Time Scale for further guidance).

Porosity is used by FLAC3D to calculate the saturated density of the medium and evaluate Biot modulus in the case when the fluid bulk modulus is given as an input. FLAC3D does not update porosity during the calculation cycle, since the process is time-consuming and only the slope of the transient response is affected. Porosity is a zone property specified using the zone fluid property porosity command.

Saturation

Saturation, \(s\), is defined as the ratio of pore volume occupied by fluid to total pore volume. In FLAC3D’s formulation not using simple-fluid mode, pore pressure is set to zero if the saturation at any point is less than 1; In FLAC3D’s formulation using simple-fluid mode, the saturation \(s\) is assumed a constant value of 1, pore pressure can be a negative value (or, suction can be a positive value), the relative saturation (or called saturation ratio), \(\hat{s}\) is a function of suction, and the permeability ratio is a function of the relative saturation. The effect of dissolved and trapped air may be allowed by reducing the local fluid modulus while keeping the saturation at 1 (i.e., we imagine that there is an equivalent fluid present throughout the pore space). Although no pore pressures are present in a partially saturated region, the trapped fluid still has weight (i.e., body forces act), and the fluid moves under the action of gravity (at a reduced apparent permeability — see this equation).

The initial saturation may be given by the user, but it is also updated during FLAC3D’s calculation cycle as necessary to preserve the mass balance. Note that in FLAC3D, saturation is not considered as an independent variable; it cannot be fixed at any gridpoint.

Note that in an unsaturated zone with non-zero pore-pressure, FLAC3D uses Bishop effective stress (the value of pore-pressure used to calculate effective stress is multiplied by zone saturation).

Undrained Thermal Coefficient

The undrained thermal coefficient, \(\beta\), is defined as the pore-pressure variation divided by 3\(\alpha M\) per unit temperature change in an undrained constrained test (no deformation).

For an ideal porous material, this coefficient is related to the thermal-expansion coefficients for the grains, \(\beta_g\), and fluid, \(\beta_f\), through the formula

(12)\[\beta = \beta_g(\alpha - n) + \beta_f n\]

The undrained thermal coefficient is a zone property, specified using the zone fluid property undrained-thermal-coefficient command.

Effective Stress Cutoff

For the specific purpose of calculating effective stress, by default the average pore-pressure is cut off at 0.0 and not allowed to go negative. This is because often the existence of negative pore-pressures are an artifact of a simple fully-saturated flow assumption, and are not considered physical. If unsaturated flow is active (see the zone fluid unsaturated command), negative pore-pressures are a valid part of the model, and it is designed to have the negative pore pressures conctribute to effective stress, then this cutoff can be changed. See the zone fluid property effective-cutoff command.

Note that in an unsaturated zone with non-zero pore-pressure, FLAC3D uses Bishop effective stress (the value of pore-pressure used to calculate effective stress is multiplied by zone saturation).

Saturated Densiy

By default FLAC3D does not include saturation in the calculation of wet density – the wet density always assumes full saturation regardless of the actual average saturation of the zone. If the zone fluid property density-saturated command is set to true, then the zone will include saturation in this calculation.

Pore Pressure Generation

By default FLAC3D does not include two-way mechanical coupling. Volumetric deformation of the material does not affect pore-pressures. If the zone fluid property pore-pressure-generation command is used to set the proeprty to true, then full two-way coupling will be activated for those zones.

Fluid Tension Limit

In fine soils, the pore water may be able to sustain a significant tension. In FLAC3D, when using the cutoff unsaturated flow model, negative pore pressures can develop, by default, up to a limit beyond which desaturation will take place. The negative pore pressure limit is set using first optional parameter to the zone fluid unsaturated cutoff command, where the value is a negative number (set to zero by default).

Note that a negative pore pressure is not the same as “tension” due to capillary, electrical or chemical forces. The latter forces are best represented by an increased effective stress within the constitutive model. Negative fluid pressures are unrelated to the fact that a material is composed of grains – the pressures simply arise from the expansion of a volume filled with fluid.