FLAC3D Theory and Background • Fluid-Mechanical Interaction

Numerical Formulation

Substitution of the fluid mass balance equation in the constitutive equation for the pore fluid yields the expression for the fluid continuity equation:

(1)1Mpt+nsst=1s(qi,i+qv)αϵt+βTt

In the FLAC3D numerical approach, the flow domain is discretized into brick-shaped zones defined by eight nodes. (The same discretization is used for mechanical and thermal calculations, when applicable.) Both pore pressure and saturation are assumed to be nodal variables. Internally, each zone is subdivided into tetrahedra (2 overlays are used by default), in which the pore pressure and the saturation are assumed to vary linearly.

The numerical scheme relies on a finite volume nodal formulation of the fluid continuity equation. The formulation can be paralleled to the mechanical constant stress formulation (presented in Finite Volume Approximation to Space Derivatives) that leads to the nodal form of Newton’s law. It is obtained by substituting the pore pressure, specific discharge vector and pore-pressure gradient for velocity vector, stress tensor and strain-rate tensors, respectively. The resulting system of ordinary differential equations is solved using an explicit mode of discretization in time (default mode). An implicit fluid flow calculation scheme is also available. However, it is only applicable to saturated flow simulations (flow in which the unit saturation is not allowed to change).

The principal results are summarized below. The formulation, as it applies to saturated flow, is described first. This formulation is then extended to apply to saturated/unsaturated flow simulations. For now we note that the saturated/unsaturated algorithm is based on variable substitution and upstream weighting techniques and conserves fluid flow.

Starting from a state of mechanical equilibrium, a coupled hydromechanical static simulation in FLAC3D involves a series of steps. Each step includes one or more flow steps (flow loop) followed by enough mechanical steps (mechanical loop) to maintain quasi-static equilibrium.

The increment of pore pressure due to fluid flow is evaluated in the flow loop; the contribution from volumetric strain is evaluated in the mechanical loop as a zone value which is then distributed to the nodes.

For the effective stress calculation, the total stress increment due to pore-pressure change arising from mechanical volume strain is evaluated in the mechanical loop, and that arising from fluid flow is evaluated in the flow loop.

Note that in FLAC3D, all material models are expressed in terms of Terzaghi effective stresses (i.e., effective stresses are used to detect failure in plastic materials). In this context, the pore-pressure field may originate from different sources: a fluid flow analysis; a coupled fluid/mechanical simulation; or an initialization with the zone gridpoint initialize pore-pressure or zone water command.

Saturated Fluid Flow

In this section, as a first step towards the derivation of the more general FLAC3D fluid flow algorithm, we consider flow in a medium which remains fully saturated at all times.

Finite-Volume Approximation to Space Derivatives

By convention, the tetrahedron nodes are referred to locally by a number from 1 to 4, face n is opposite node n, and the superscript (f) relates to the value of the associated variable on face f.

A linear pore-pressure variation and constant fluid density are assumed within a tetrahedron. The pressure head gradient, expressed in terms of nodal values of the pore pressure by application of the Gauss divergence theorem, may be written as

(2)(pρfxigi),j=13V4l=1(plρfxligi) n(l)j S(l)

where n(l) is the exterior unit vector normal to face l, S is the face surface area and V is the tetrahedron volume.

For numerical accuracy, the quantity xix1i, where x1i corresponds to the coordinates of one of the tetrahedron’s corners, may be substituted for xi in the expression of the pressure head in Equation (2) without affecting the value of the pressure head gradient. Accordingly, Equation (2) takes the form

(3)(pρfxigi),j=13V4l=1pl n(l)j S(l)

where the nodal quantity pl is defined as

(4)pl=plρf(xlix1i)gi

Nodal Formulation of the Mass Balance Equation

For fully saturated flow (s is constant and fixed to one), the fluid continuity (Equation (1)) may be expressed as

(5)qi,i+b=0

where

(6)b=1Mptqv

is the equivalent of the instantaneous “body forces,” ρbi, used in the mechanical node formulation, and

(7)qv=qvαϵt+βTt

First consider a single tetrahedron. Using this analogy, the nodal discharge, Qne [m3/s], n = 1,4, equivalent to the tetrahedron-specific discharge and the volumetric source intensity, b, may be expressed as (see Finite Volume Approximation to Space Derivatives)

(8)Qne=QntqvV4+mn dpdtn

where

(9)Qnt=qin(n)iS(n)3

and

(10)mn=V4Mn

In principle, the nodal form of the mass balance equation is established by requiring that, at each global node, the sum of equivalent nodal discharges (Qne) from tetrahedra meeting at the node (averaged over two overlays) and nodal contribution (Qnw) of applied boundary fluxes and sources be zero.

The components of the tetrahedron-specific discharge vector in Equation (9) are related to the pressure-head gradient by means of the transport law (see equation). In turn, the components of the pressure-head gradient can be expressed in terms of the tetrahedron nodal pore pressures, using Equation (3).

In order to save computer time, local matrices are assembled. A zone formulation is adopted, whereby the sum Qnz of contributions Equation (9) from all tetrahedra belonging to the zone and meeting at a node, n, is formed and averaged over two overlays. Local zone matrices [M] that relate the eight nodal values {Qz} to the eight nodal pressure heads {p} are assembled. Because these matrices are symmetrical, 36 components are computed; these are saved at the beginning of the computation and updated every ten steps in large-strain mode. By definition of zone matrices, we have

(11)Qnz=Mnjpj

where {p} is the local vector of nodal pressure heads for the zone.

In turn, global nodal values QnT are obtained by superposition of zone contributions. Taking some liberty with the notation, we write for each global node n,

(12)QnT=Cnjpj

where [C] is the global matrix and {p} is, in this context, the global vector of nodal pressure heads.

Returning to our previous consideration, we write

(13)Qne+Qnw=0

where, for simplicity of notation, the sign is used to represent summation of the contributions at global node n of all zones meeting at that node. (A zone contribution consists of contributions of tetrahedra involved in the zone, averaged over two overlays.) Using Equation (8) for Qne and Equation (8) for qv in Equation (13), we obtain, after some manipulation,

(14)dpndt=1mn[QnT+Qnapp+Qnthm]

where QnT is a function of the nodal pore pressures defined in Equation (12) together with Equation (4), Qnapp is the known contribution of applied volume sources, boundary fluxes and point sources, having the form

(15)Qnapp=[qvV4+Qw]n

and Qnthm is the thermal-mechanical nodal discharge contribution defined as

(16)Qnthm=ddt[(αϵV4)n(βTV4)n]

Equation (14) is the nodal form of the mass-balance equation at node n; the right-hand side term QnT+Qnapp+Qnthm is referred to as the out-of-balance discharge. This term is made up of two contributions: the out-of-balance flow discharge, QnT+Qnapp; and the out-of-balance thermal-mechanical discharge, Qnthm. When the fluid is at rest (the fluid-flow calculation is turned off), the out-of-balance flow discharge vanishes and pore-pressure changes are caused by mechanical and/or thermal deformations only. For flow calculation only (no thermal-mechanical coupling), the out-of-balance thermal-mechanical discharge is equal to zero.

In FLAC3D, the Biot modulus is a nodal property, and we may write, using Equation (10),

(17)1mn=MnVn

where

(18)Vn=(V4)n

After substitution of Equation (16) in Equation (14), we write, rearranging terms,

(19)ddt[pnpnv]=MnVn[QnT+Qnapp]

where

(20)pnv=MnVn[(αϵV4)n(βTV4)n]

An expression such as Equation (19) holds at each global node involved in the discretization. Together, these expressions form a system of ordinary differential equations that is solved in FLAC3D, for given dpnv/dt, using either explicit or implicit (saturated flow only) finite-volume schemes. The domain of application of each scheme is discussed below.

Explicit Finite-Volume Formulation

In the explicit formulation, the value of the quantity ppv at a node is assumed to vary linearly over a time interval Δt. In principle, the derivative in the left-hand side of Equation (19) is expressed using forward finite differences, and the out-of-balance flow-discharge is evaluated at time t. Starting with an initial pore-pressure field, nodal pore pressures at incremental time values are updated, provided the pore pressure is not fixed, using the expression

(21)pn<t+Δt>=pn<t>+Δpnv<t>+Δpn<t>

where

(22)Δpn<t>=χn[QnT<t>+Qnapp<t>]
(23)χn=MnVnΔt

and

(24)Δpnv<t>=MnVn[(αΔϵV4)n(βΔTV4)n]<t>

Numerical stability of the explicit scheme can only be achieved if the timestep remains below a limiting value.

Stability Criterion

To derive the stability criterion for the flow calculation, we consider the situation in which a node, n, in an assembly of zones is given a pore-pressure perturbation, p0, from a zero initial state. Using Equation (12), we obtain

(25)QnT=Cnn p0

If node n belongs to a leaky boundary, we have

(26)Qnapp=Dnn p0

where Dnn is used to represent the pressure coefficient in the global leakage term at node n.

After one fluid-flow timestep, the new pore pressure at node n is (see Equation (21) through Equation (23))

(27)pn<Δt>=p0[1MnVn(Cnn+Dnn)Δt]

To prevent alternating signs of pore pressures as the computation is repeated for successive Δt, the coefficient of p0 in the preceding relation must be positive. Such a requirement implies that

(28)Δt<VnMn1Cnn+Dnn

The value of the poro-mechanical timestep used in FLAC3D is the minimum nodal value of the right-hand side in Equation (28), multiplied by the safety factor of 0.8.

To assess the influence of the parameters involved, it is useful to keep in mind the following representation of the critical timestep. If Lc is the smallest tetrahedron characteristic length, we may write an expression of the form

(29)Δtcr=1a[cL2c+hMLc]1

where Δtcr is the critical timestep, c is the fluid diffusivity, and a is a constant, larger than unity, that depends on the medium geometrical discretization (e.g., a = 6 for a regular discretization in cubes; see Karlekar and Desmond (1982) for a thermal analogy).

The critical timestep in Equation (29) corresponds to a measure of the characteristic time needed for the diffusion “front” to propagate throughout the tetrahedron. To estimate the time needed for the front to propagate throughout a particular domain, a similar expression can be used, provided Lc is interpreted as the characteristic length of the domain under consideration. Consider the case where no leakage occurs and the properties are homogeneous. By taking the ratio of characteristic times for the domain and the tetrahedron, we see that the number of steps needed to model the propagation of the diffusion process throughout that domain is proportional to the ratio of square characteristic lengths for the domain and the tetrahedron. That number may be so large that the use of the explicit method alone could become prohibitive. However, in many groundwater problems, the advantage of this first-order method is that the calculated timestep is usually small enough to follow nodal pore-pressure fluctuations accurately.

Implicit Finite-Volume Formulation

The requirement that Δt should be restricted in size to ensure stability sometimes results in an extremely small timestep (of the order of a fraction of a second), especially when transient flow in layers of contrasting properties is involved. The implicit formulation eliminates this restriction, but it involves solving simultaneous equations at each timestep and, besides, it only applies to saturated fluid-flow simulations.

The implicit formulation in FLAC3D uses the Crank-Nicolson method, in which the pore pressure ppv at a node is assumed to vary quadratically over the time interval Δt. In this method, the derivative d(ppv)/dt in Equation (19) is expressed using a central difference formulation corresponding to a half timestep, while the out-of-balance flow discharge is evaluated by taking average values at t and t+Δt. In this formulation we have

(30)pn<t+Δt>=pn<t>+Δpnv<t>+Δpn<t>

and

(31)Δpn<t>=χn[12(QnT<t+Δt>+QnT<t>)+¯Qnapp]

where

(32)χn=MnΔtVn

Δpnv<t> is given by Equation (24), and

(33)¯Qnapp=12(Qnapp<t+Δt>+Qnapp<t>)

From Equation (12), we may write

(34)QnT<t>=Cnjpj<t>

and

(35)QnT<t+Δt>=Cnjpj<t+Δt>

For no variation of the gravity term during the time interval Δt, Equation (30) can be written using Equation (4).

(36)pn<t+Δt>=pn<t>+Δpnv<t>+Δpn<t>

Using this last expression and Equation (34), Equation (35) becomes

(37)QnT<t+Δt>=QnT<t>+CnjΔpj<t>+CnjΔpjv<t>

After substitution of Equation (37) into Equation (31), we obtain, using Equation (35),

(38)Δpn<t>=χn[Cnj(pj<t>+12Δpjv<t>)+12CnjΔpj<t>+¯Qnapp]

Finally, regrouping terms:

(39)[δnjχn2Cnj]Δpj<t>=χn[Cnj(pj<t>+12Δpjv<t>)+¯Qnapp]

For simplicity of notation, we define the known matrix [A] and vector [b<t>] as

(40)Anj=δnjχn2Cnj

and

(41)bn<t>=χn[Cnj(p<t>j+12Δpjv<t>)+¯Qnapp]

With this convention, we may write Equation (39) in the form

(42)AnjΔpj<t>=bn<t>

One equation holds for each of the nodes involved in the grid, and the resulting system is solved in FLAC3D using Jacobi’s iterative method. In this approach, pore-pressure increments at iteration r+1 and node n are calculated using the recurrence relation

(43)Δpn<r+1><t>=Δpn<r><t>+1Ann[AnjΔpj<r><t>+bn<t>]

where Einstein’s summation convention applies to index j only. Using the Equation (40) of [A], Equation (43) takes the form

(44)Δpn<r+1><t>=Δpn<r><t>+11χn2Cnn[χn2CnjΔpj<r><t>Δpn<r><t>+bn<t>]

The initial approximation is chosen such that

(45)Δpn<0><t>=0

(Note that in FLAC3D, pressure-dependent boundary conditions (contained in bn) are updated in the implicit iterative procedure.)

Convergence Criterion

In FLAC3D a minimum of 3 and a maximum of 500 iterations are considered, and the criterion for detection of convergence has the form

(46)maxn  |  Δpn<r+1><t>Δpn<r><t>  |  <102 (maxn  |  Δpn<r><t>  |)

The magnitude of the timestep must be selected in relation to both convergence and accuracy of the implicit scheme. Although the Crank-Nicolson method is stable for all positive values of Δt (for no leakage), the convergence of Jacobi’s method is not unconditionally guaranteed unless the matrix [A] is strictly diagonally dominant:

(47)nj=1jk  |  Akj  |  <  |  Akk  |

for 1kn (see Dahlquist and Bjorck 1974). According to the definitions of Anj ( Equation (40)) and χn ( Equation (32)), this sufficient condition is always fulfilled for sufficiently small values of Δt. If convergence of Jacobi’s method is not achieved, an error message is issued. It is then necessary to either reduce the magnitude of the implicit timestep, or use the explicit method. (The explicit timestep may be used as a lower bound.)

Also, although the implicit method is second-order accurate, some insight may be needed in order to select the appropriate timestep. Indeed, its value must remain small compared to the wavelength of any nodal pore-pressure fluctuation. Typically, the explicit method is used earlier in the run or in its perturbed stages, while the implicit method is preferred for the remainder of the simulation. (Alternatively, the implicit method could be used with the explicit timestep value for extra accuracy.)

Remember that the implicit approach should not be selected for simulations in which changes of saturation are to be expected. In addition, computation time and computer memory are two factors that must be taken into consideration when selecting the implicit approach in FLAC3D. In the implicit method, a set of equations requiring a minimum of three iterations must be solved at each timestep. The amount of calculation required for one iterative step is approximately equal to that needed for one timestep in the explicit scheme. Also, intermediate values must be stored in the iterative procedure, requiring extra memory to be allocated in comparison with the explicit scheme. Those inconveniences, however, can be more than offset by the much larger timestep generally permitted by the implicit method, or by the gain in accuracy.

Saturated/Unsaturated Flow

For saturated fluid flow, the nodal volumetric flow rates in a zone {Qz} are related to the nodal pore pressures {p} by Equation (11), which may be expressed in matrix notation as

(48){Qz}=[M]{pρfxigi}

where the matrix [M] is a known function of the zone geometry and the saturated value of the mobility coefficient. This relation was derived by application of the Gauss divergence theorem for pore-pressure gradient in the tetrahedron, taking into account energy considerations, and using superposition. It is extended to the case of unsaturated flow in coarse soils (constant air pressure, no capillary pressure) by application of modifications:

  1. The gravity term ρfxigi is multiplied by the average zone saturation, ˉs, to account for partial zone filling.

  2. The nodal flow rates are multiplied by the relative mobility, ˆk (see equation), which is a function of the average saturation at the inflow nodes for the zone, ˆsin. This upstream weighting technique is applied to prevent unaccounted flow in the filling process of an initially dry medium (numerical dispersion caused by the nodal definition of saturation). The relative permeability function has the same form as this equation:

    (49)ˆk(ˆsin)=ˆsin2(32ˆsin)
  3. Nodal inflow rates are scaled according to local saturation.

For saturated/unsaturated flow, the explicit finite volume nodal formulation of the fluid continuity equation (1) takes the form (see Equation (21) through Equation (23))

(50)ΔpnM+Δss=1snV[(QT+Qapp)Δt+s(αΔϵVT4βΔTVT4)]

where the subscript <t> and superscript n have been omitted for simplicity of notation, and VT stands for tetrahedron volume. At a saturated node, s = 1 and Equation (50) becomes

(51)Δp=MV[(QT+Qapp)Δt+(αΔϵVT4βΔTVT4)]

At a partially saturated node, p = 0 and Equation (50) may be expressed as

(52)Δs=1nV[(QT+Qapp)Δt+s(αΔϵVT4βΔTVT4)]

The approach leads to a technique of variable substitution whereby the saturation variable is updated using Equation (52) at an unsaturated node (where the pore pressure is equal to zero), and Equation (51) is used to increment the pore-pressure variable at a saturated node (where saturation is constant). Those relations are applied in such a way that the fluid volume balance is preserved (i.e., Equation (50) is satisfied).

The explicit fluid timestep in FLAC3D is selected on the basis of a stability criterion derived using the saturated fluid diffusivity c=kM (k is the largest principal value of saturated mobility coefficient tensor, M is Biot modulus). (See Equation (28) and Equation (29).) The time scale associated with fluid storage is usually much smaller than the characteristic time scale for phreatic storage. While an estimate of the first one (two-dimensional flow) can be calculated using tL2/c, where L is to be interpreted as average flow path length, a rough estimate of the second may be obtained from the relation tL2/(kˉp/n), where L is the average height of the medium associated with phreatic storage, and ˉp is the related mean pore-pressure change.

Mechanical Timestep and Numerical Stability

The presence of the fluid will increase the apparent mechanical bulk modulus of the medium which, in turn, affects the magnitude of the nodal mass in the density scaling scheme used for numerical stability in FLAC3D. An upper bound for this apparent modulus may be found by considering undrained, saturated isothermal poro-mechanical conditions for which ζt = 0, st = 0 and Tt = 0. The incremental expression of the fluid constitutive law then takes the form

(53)Δp=αMΔϵ

Adopting the incremental form of the elastic law to describe the volumetric part of the stress-strain constitutive relation in a time interval, we can write

(54)13Δσii+αΔp=KΔϵ

Substitution of Equation (53) in Equation (54) gives

(55)13Δσii=(K+α2M)Δϵ

This apparent tangent bulk modulus, K+α2M, is used in place of K when calculating nodal inertial mass (see this equation in the Theory and Background section).

Also note that in FLAC3D, the dry density, ρd, of the medium is considered in the input. The saturated density, ρ=ρd+nsρf, is calculated internally using the input value for fluid density, ρf, porosity, n, and saturation, s. Body forces are then adjusted accordingly in Newton’s equations of motion.

Total Stress Correction

In the FLAC3D formulation, nodal pore pressures are calculated first. The zone pore pressures are then derived by volume averaging of the tetrahedra values (obtained as arithmetic averages of nodal values).

The total stresses are then adjusted, consistent with the definition of Biot effective stress, by adding an increment Δσij to account for the contribution from pore pressure changes (see this equation). The stress increment Δσij is expressed as

(56)Δσij=Δσfij+Δσthij

where

(57)Δσfij=α¯Δpn<t>δij

is the contribution from flow, and

(58)Δσthij=α¯M(α¯Δϵβ¯ΔT)<t>δij

is the contribution from thermal-mechanical coupling. Also, in those expressions, the overline symbol stands for zone average in the sense mentioned above.

Fully Saturated Fast Flow

In the standard fluid-flow formulation in FLAC3D the fluid continuity equation for fully saturated flow of a compressible fluid is given (based upon this equation for isothermal conditions) as

(59)1Mpt=ζtαϵt

where ζ is the variation of fluid content (with respect to a reference state), defined as the variation of fluid volume per unit volume of porous material, ϵ is the volumetric strain, M is Biot modulus, and α is Biot coefficient. The fluid mass balance for the standard formulation is given in Governing Differential Equations by this equation, the balance of momentum by this equation, and the material constitutive law has the form given by this equation.

The standard calculation procedure is applicable to cases of relatively high material stiffness (compared to the fluid stiffness) and average porosity. The performance of the standard scheme (in terms of solution speed and accuracy) degrades when a large value for Biot modulus (compared to the solid matrix stiffness) is used. A high value is generally characterized by the equation

(60)M>20(K+4/3G)

or, if α = 1, and a combination of high fluid bulk modulus and low porosity is considered,

(61)Kf>20(K+4/3G)n

Kf is fluid bulk modulus, n is porosity, and K and G are drained elastic moduli of the solid matrix.

There are several disadvantages to working with the standard scheme with high values of Biot or fluid bulk modulus (i.e., when Equation (60), or Equation (61), holds):

  1. The stable explicit fluid timestep can be extremely small, because it is inversely proportional to Biot modulus.
  2. Mechanical convergence is very slow, as a consequence of the additional stiffness, α2M, due to the effect of the fluid on mass scaling.
  3. The accuracy of the solution is poor because the standard scheme is not well-adapted to solve the degenerate form of the pore pressure equation, Equation (59).

The poor performance of the standard scheme develops because, when M is very large, Equation (59) degenerates to the expression

(62)ϵunbt=0

where the quantity

(63)ϵunbt=ζtαϵt

is defined as the unbalanced volumetric strain rate. The standard scheme is not well-designed to solve this degenerate form of the pressure-based equation.

An Alternative Fast-Flow Algorithm

An alternative strain-based formulation (saturated fast-flow logic) is better adapted for Equation (59) for high values of Biot modulus.[1]

To illustrate the principle of this scheme, consider, for simplicity, the case of a material constitutive law with elastic volumetric behavior, expressed as

(64)σt+αpt=Kϵt

In this scheme, the quantity ϵunbt is evaluated from Equation (63) after performing a flow step and a mechanical step. In general, the quantity will not be negligible, and an iterative scheme is applied to reduce the value of the quantity and enforce Equation (62) . First, the impact on Biot effective stress of the “unbalanced” volumetric strain rate contribution is evaluated; this quantity, Kϵunbt, is interpreted as a pore pressure increment, αpt, which is then used to “correct” the (total) stress (see Equation (64) ), using

(65)σn=σαpt

where

(66)αpt=Kϵunbt

The procedure is repeated in the course of a series of mechanical steps until the unbalanced volumetric strain increment falls below a negligible value. Gridpoint pore pressures are also updated to ensure proper feedback in the fluid-mechanical calculation.

For stability and convergence of the scheme, the magnitude of the stress correction increment αpt=Kϵunbt cannot be applied in full at each calculation step. Instead, the accumulated quantity is leaked gradually into the grid using a “balloon” technique.[2] In this technique, an unbalanced volume, Vunb, is calculated and stored in virtual containers (balloons) associated with each gridpoint (or node). There are two contributions to the unbalanced volume at a gridpoint: one from change of fluid content and one from volumetric strain increment of zones meeting at the gridpoint. The content of the balloon is updated as a result of volumetric deformation, after each mechanical step, and unbalanced flow, after each fluid-flow step. (The variation of fluid content is evaluated at gridpoints using the fluid mass balance equation, as in the standard scheme.)

The unbalanced volumetric strain expression for the zone is calculated using

(67)ϵunb=Vunb,zVzαΔϵ

where Vz is zone volume, and the unbalanced volume for the zone, Vunb,z, comes from the contributions of the gridpoints associated with the zone. The pore pressure increment used to correct normal stresses in a mechanical step (see (65)) is calculated in proportion to the current unbalanced volumetric strain:

(68)pt=FpKϵunb

where Fp is the relaxation coefficient. Also, the quantity αΔϵV is distributed from one zone to balloons at gridpoints (using weighting functions) after each mechanical step. Finally, the pore pressure at a gridpoint is updated using the nodal form of (68):

(69)pt=FpKVunbVn

where Vn is the gridpoint volume.

The relaxation coefficient, Fp, is a dimensionless number which must be small enough to enforce convergence of the fast-flow scheme. The default value is 0.004. This value should be adequate to enforce numerical stability for most problems. It is possible to reset the value using the zone fluid fastflow-relaxation command. (Additional information and recommendations for Fp are given in Section 1.4.1.1 of the Fluid-Mechanical Interaction volume of the FLAC Version 7 manual.)

In an undrained simulation, the only contribution to unbalanced volumetric strain in the balloon comes from mechanical volumetric changes. The volumetric changes are converted gradually into pore pressure changes during the mechanical iterations. The process is interrupted when the balloon is emptied.

In a fluid-mechanical simulation, the flow contribution to unbalanced volumetric strain is added in the balloon during the flow calculation step, and enough mechanical steps are taken afterwards to dissipate this contribution and that from the mechanical volumetric changes.

During cycling using the saturated fast-flow logic, the average unbalanced volume ratio, ˆVaunb, and maximum unbalanced volume ratio, ˆVmunb, are monitored to detect convergence of the scheme:

(70)ˆVaunb=max[KnVunb|σ|Vn]
(71)ˆVmunb=KnVunb|σ|[Vn]2

where the summation and maximum are taken over all gridpoints in the model, Kn is the average bulk modulus of the zones meeting at the gridpoint, and |σ| is the maximum value of the stress norm in the model. The balloon is considered empty when both ˆVaunb and ˆVmunb fall below 105. This default value can be changed with the model solve unbalanced-average and model solve unbalanced-maximum commands.

The saturated fast-flow scheme is applicable only for the case of an incompressible fluid. Otherwise, the standard scheme should be used (i.e., when (60), or (61), does not apply). Therefore, the two schemes complement each other.

When (60), or (61), does apply, the coupled, incompressible saturated fast-flow scheme has several advantages over the standard scheme:

  1. Large fluid timesteps can be taken. The dependence of the stable explicit timestep on the inverse of Biot modulus (or on porosity and inverse of fluid bulk modulus) is removed in this scheme. The stable explicit timestep is inversely proportional to FpK in the saturated fast-flow scheme. In fact, the Biot modulus, M, in Equation (59) is replaced by FpK in the fast-flow scheme (see (69) and (63)).
  2. Convergence is much faster. There is no need to account for fluid stiffness during the mass scaling calculation in the fast-flow scheme. This allows much faster mechanical convergence.
  3. There is better accuracy. The degenerate form of the coupled fluid-mechanical equations, which apply for large values of Biot modulus, are solved directly in the fast-flow scheme. This provides a more accurate solution.

At present, the saturated fast-flow scheme cannot be used for several conditions:

  1. unconfined fluid flow, in which a phreatic surface, or water table, exists or develops;
  2. models with interfaces;
  3. models with attached grids; and
  4. dynamic analysis.

Two verification problems for the saturated fast-flow logic are provided. See One-Dimensional Consolidation and Consolidation Settlement at the Center of a Strip Load.

Endnotes

[1]This scheme is derived from the saturated fast-flow scheme already available in FLAC. See Section 1.4.1 of the Fluid-Mechanical Interaction volume of the FLAC Version 7 manual for reference.
[2]See Section 1.4.1.1 of the Fluid-Mechanical Interaction volume of the FLAC Version 7 manual for reference.