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 permeability coefficient, \(k\), the fluid mass density, \({\rho}_f\), and either 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). 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.

Permeability Coefficient

The isotropic permeability coefficient, \(k\) (e.g., m²/(Pa-sec) in SI units), used in FLAC3D is also referred to in the literature as the mobility coefficient. 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 permeability across the grid, the timestep will be controlled by the largest permeability (see this equation). For problems in which steady state (but not transient behavior) is required, it may be beneficial to limit the variations in permeability 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 permeability, compared to a 200:1 variation, when the flow is moving entirely from one permeability region to the other.

The permeability 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), 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 mode. The zone water command (or zone gridpoint initialize 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. [DR: This should be added in :cmd:`zone.property] The fluid density can be imposed globally using the zone water density command, or it can be allowed to vary with position by using the zone initialize 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 calculation mode for compressible grains is turned on with the zone fluid biot on command. The Biot coefficient is a zone property specified using the zone fluid property biot command. The Biot modulus is a gridpoint variable specified using the zone gridpoint initialize 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 (i.e., zone fluid biot off). The fluid modulus is a gridpoint variable specified using the zone gridpoint initialize 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.

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).

Fluid Moduli for Drained and Undrained Analyses

In FLAC3D, whenever the fluid bulk modulus (or Biot modulus) is selected and model configure fluid 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. In this case, the undrained bulk modulus \(K_u\) for the solid matrix should be specified. The undrained bulk modulus is

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

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

(9)\[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

(10)\[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.

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

(11)\[\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.

Fluid Tension Limit

In fine soils, the pore water may be able to sustain a significant tension. In FLAC3D, 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 the zone gridpoint initialize fluid-tension command, where the value is a negative number (set to -10¹⁵, 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.