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)\[{1 \over M} {{\partial p} \over {\partial t}} + {n \over s} {{\partial s} \over {\partial t}}= {1 \over s} (-q_{i,i} + q_v) - \alpha {{\partial \epsilon} \over {\partial t}} + \beta {{\partial T} \over {\partial t}}\]

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 - {\rho}_f x_i g_i)_{,j} = - {{1}\over{3V}} \sum_{l=1}^4 (p^l - {\rho}_f x_i^l g_i) \ n_j^{(l)}\ 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 \(x_i - x_i^1\), where \(x_i^1\) corresponds to the coordinates of one of the tetrahedron’s corners, may be substituted for \(x_i\) 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 - {\rho}_f x_i g_i)_{,j} = - {{1}\over{3V}} \sum_{l=1}^4 p^{*l} \ n_j^{(l)}\ S^{(l)}\]

where the nodal quantity \(p^{*l}\) is defined as

(4)\[p^{*l} = p^l - {\rho}_f (x_i^l - x_i^1) g_i\]

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)\[q_{i,i} + b^* = 0\]

where

(6)\[b^* = {1 \over M} {{\partial p}\over{\partial t}} - q_v^*\]

is the equivalent of the instantaneous “body forces,” \(\rho b_i\), used in the mechanical node formulation, and

(7)\[q_v^* = q_v - \alpha {{\partial \epsilon} \over {\partial t}} + \beta {{\partial T} \over {\partial t}}\]

First consider a single tetrahedron. Using this analogy, the nodal discharge, \(Q_e^n\) [m³/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)\[Q_e^n = Q_t^n - {{q_v^* V}\over{4}} + m^n\ {{dp}\over{dt}}^n\]

where

(9)\[Q_t^n = {{q_i n_i^{(n)} S^{(n)}}\over 3}\]

and

(10)\[m^n = {{V}\over{4M^n}}\]

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 (\(-Q_e^n\)) from tetrahedra meeting at the node (averaged over two overlays) and nodal contribution (\(Q_w^n\)) 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 \(Q_z^n\) 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 \(\{ Q_z \}\) 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)\[Q_z^n = M_{nj} p^{*j}\]

where \(\{ p^* \}\) is the local vector of nodal pressure heads for the zone.

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

(12)\[Q_T^n = C_{nj} p^{*j}\]

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)\[- \sum Q_e^n + \sum Q_w^n = 0\]

where, for simplicity of notation, the \(\sum\) 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 \(Q_e^n\) and Equation (8) for \(q_v^*\) in Equation (13), we obtain, after some manipulation,

(14)\[{{d p^n}\over{d t}} = - {{1} \over {\sum m^n}} \left[Q_T^n + \sum Q_{app}^n+ \sum Q_{thm}^n \right]\]

where \(Q_T^n\) is a function of the nodal pore pressures defined in Equation (12) together with Equation (4), \(\sum Q_{app}^n\) is the known contribution of applied volume sources, boundary fluxes and point sources, having the form

(15)\[\sum Q_{app}^n = - \sum {\left[q_v {{V}\over{4}} + Q_w\right]}^n\]

and \(\sum Q_{thm}^n\) is the thermal-mechanical nodal discharge contribution defined as

(16)\[\sum Q_{thm}^n = {{d} \over {d t}} \left[ \sum {\left( \alpha \epsilon {{V}\over{4}} \right)}^n - \sum {\left( \beta T {{V}\over{4}} \right)}^n \right]\]

Equation (14) is the nodal form of the mass-balance equation at node \(n\); the right-hand side term \(Q_T^n + \sum Q_{app}^n+ \sum Q_{thm}^n\) is referred to as the out-of-balance discharge. This term is made up of two contributions: the out-of-balance flow discharge, \(Q_T^n + \sum Q_{app}^n\); and the out-of-balance thermal-mechanical discharge, \(\sum Q_{thm}^n\). 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)\[{{1} \over {\sum m^n}} = {{M^n} \over{V^n}}\]

where

(18)\[V^n = \sum {\left( {V \over 4} \right)}^n\]

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

(19)\[{{d}\over{d t}} \left[ p^n - p_v^n \right] = - {{M^n} \over {V^n}} \left[Q_T^n + \sum Q_{app}^n \right]\]

where

(20)\[p_v^n = - {{M^n} \over {V^n}} \left[ \sum {\left( \alpha \epsilon {{V}\over{4}} \right)}^n - \sum {\left( \beta T {{V}\over{4}} \right)}^n \right]\]

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 \(d p_v^n / 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 \(p - p_v\) at a node is assumed to vary linearly over a time interval \(\Delta 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)\[p_{<t+\Delta t>}^n = p_{<t>}^n + \Delta p_{v<t>}^n + \Delta p_{<t>}^n\]

where

(22)\[\Delta p_{<t>}^n = {\chi}^n \left[ Q_{T <t>}^n + \sum Q_{app<t>}^n \right]\]

(23)\[{\chi}^n = - {{M^n}\over{V^n}} \Delta t\]

and

(24)\[\Delta p_{v<t>}^n = - {{M^n} \over {V^n}} {\left[ \sum {\left( \alpha \Delta \epsilon {{V}\over{4}} \right)}^n - \sum {\left( \beta \Delta T {{V}\over{4}} \right)}^n \right]}_{<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, \(p_0\), from a zero initial state. Using Equation (12), we obtain

(25)\[Q_T^n = C_{nn}\ p_0\]

If node \(n\) belongs to a leaky boundary, we have

(26)\[\sum Q^n_{app} = D_{nn}\ p_0\]

where \(D_{nn}\) 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)\[p_{<\Delta t>}^n = p_0 \left[ 1 - {{M^n}\over{V^n}} (C_{nn} + D_{nn}) \Delta t \right]\]

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

(28)\[\Delta t < {{V^n} \over {M^n}} {{1} \over {C_{nn} + D_{nn}}}\]

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 \(L_c\) is the smallest tetrahedron characteristic length, we may write an expression of the form

(29)\[\Delta t_{cr} = {{1}\over{a}} {\left [ {{c}\over{L_c^2}} + {{h M}\over{L_c}} \right ]}^{-1}\]

where \(\Delta t_{cr}\) 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 \(L_c\) 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 \(\Delta 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 \(p - p_v\) at a node is assumed to vary quadratically over the time interval \(\Delta t\). In this method, the derivative \(d(p - p_v)/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 + \Delta t\). In this formulation we have

(30)\[p^n_{<t + \Delta t>} = p_{<t>}^n + \Delta p_{v<t>}^n + \Delta p_{<t>}^n\]

and

(31)\[\Delta p_{<t>}^n = {\chi}^n \left[{{1}\over{2}} \left(Q_{T <t + \Delta t>}^n + Q_{T <t>}^n \right) + \sum \overline{Q^n_{app}}\right]\]

where

(32)\[{\chi}^n = - {{M^n \Delta t}\over{V^n}}\]

\(\Delta p_{v<t>}^n\) is given by Equation (24), and

(33)\[\sum \overline{Q^n_{app}} = {1 \over 2} \left( \sum Q^n_{app<t + \Delta t>} + \sum Q^n_{app<t>} \right)\]

From Equation (12), we may write

(34)\[Q_{T<t>}^n = C_{nj} p_{<t>}^{*j}\]

and

(35)\[Q_{T<t + \Delta t>}^n = C_{nj} p_{<t + \Delta t>}^{*j}\]

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

(36)\[p^{*n}_{<t + \Delta t>} = p_{<t>}^{*n} + \Delta p_{v<t>}^n + \Delta p_{<t>}^n\]

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

(37)\[Q_{T<t + \Delta t>}^n = Q_{T<t>}^n + C_{nj} \Delta p_{<t>}^j + C_{nj} \Delta p_{v<t>}^j\]

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

(38)\[\Delta p_{<t>}^n = {\chi}^n \left[C_{nj} (p_{<t>}^{*j} + {{1}\over{2}} \Delta p_{v<t>}^j) + {{1}\over{2}} C_{nj} \Delta p_{<t>}^j + \sum \overline{Q_{app}^n}\right]\]

Finally, regrouping terms:

(39)\[\left[\delta_{nj} - {{\chi}^n\over{2}} C_{nj}\right] \Delta p_{<t>}^j = {\chi}^n \left[C_{nj} (p_{<t>}^{*j} + {{1}\over{2}} \Delta p_{v<t>}^j) + \sum \overline{Q_{app}^n}\right]\]

For simplicity of notation, we define the known matrix \([A]\) and vector \([b_{<t>}]\) as

(40)\[A_{nj} = \delta_{nj} - {{\chi}^n\over{2}} C_{nj}\]

and

(41)\[b_{n<t>} = {\chi}^n \left[C_{nj} (p{<t>}^{*j} + {{1}\over{2}} \Delta p_{v<t>}^j) + \sum \overline{Q_{app}^n}\right]\]

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

(42)\[A_{nj} \Delta p_{<t>}^j = b_{n<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)\[\Delta p_{<t>}^{n <r+1>} = \Delta p_{<t>}^{n <r>} + {{1}\over{A_{nn}}} \left[-A_{nj} \Delta p_{<t>}^{j <r>} + b_{n<t>}\right]\]

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

(44)\[\Delta p_{<t>}^{n <r+1>} = \Delta p_{<t>}^{n <r>} + {{1}\over{1 - {{\chi}^n\over{2}} C_{nn}}} \left[{{\chi}^n\over{2}} C_{nj} \Delta p_{<t>}^{j <r>} - \Delta p_{<t>}^{n <r>} + b_{n<t>}\right]\]

The initial approximation is chosen such that

(45)\[\Delta p_{<t>}^{n <0>} = 0\]

(Note that in FLAC3D, pressure-dependent boundary conditions (contained in \(b_n\)) 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)\[\max\limits_n \ \ \big\vert\ \ \Delta p^{n<r+1>}_{<t>} - \Delta p^{n<r>}_{<t>} \ \ \big\vert\ \ < 10^{-3}\ \left( \max\limits_n \ \ \big\vert\ \ \Delta p^{n<r>}_{<t>} \ \ \big\vert \right)\]

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 \(\Delta t\) (for no leakage), the convergence of Jacobi’s method is not unconditionally guaranteed unless the matrix \([A]\) is strictly diagonally dominant:

(47)\[\sum\limits_{{j=1}\atop{j\ne k}}^n \ \ \bigg\vert\ \ A_{kj} \ \ \bigg\vert\ \ < \ \ \bigg\vert\ \ A_{kk} \ \ \bigg\vert\]

for \(1 \le k \le n\) (see Dahlquist and Bjorck 1974). According to the definitions of \(A_{nj}\) ( Equation (40)) and \(\chi^n\) ( Equation (32)), this sufficient condition is always fulfilled for sufficiently small values of \(\Delta 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.)

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.

Advanced solvers used in fluid implicit logic

Two new advanced solvers are available in fluid implicit logic in addition to the iterative Jacobi solver. These are the direct solver and preconditioned conjugate gradient solver. Both solvers are unconditionally stable with regard to the implicit timestep. For more information about the solvers, see Pyatigorets and Russell, 2020.

Saturated/Unsaturated Flow

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

(48)\[\{ Q_z \} = [M] \{ p - \rho_f x_i g_i \}\]

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:

The gravity term \(\rho_f x_ig_i\) is multiplied by the average zone saturation, \(\bar s\), to account for partial zone filling.
The nodal flow rates are multiplied by the relative mobility, \(\hat{k}\) (see equation), which is a function of the average saturation at the inflow nodes for the zone, \(\hat{s}_{in}\). 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)\[\hat{k}(\hat{s}_{in}) = \hat{s}_{in}^2(3 - 2\hat{s}_{in})\]
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)\[{{\Delta p} \over{nM}} + {{\Delta s }\over{s}}= - {{1} \over {snV}} {\left[ \left( Q_{T} + \sum Q_{app} \right) \Delta t + s \left( \sum \alpha \Delta \epsilon {{V_T}\over{4}} - \sum \beta \Delta T {{V_T}\over{4}} \right) \right]}\]

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

(51)\[\Delta p = - {{M} \over {V}} {\left[ \left( Q_{T} + \sum Q_{app} \right) \Delta t + \left( \sum \alpha \Delta \epsilon {{V_T}\over{4}} - \sum \beta \Delta T {{V_T}\over{4}} \right) \right]}\]

At a partially saturated node, it is assumed that \(p\) = 0 and Equation (50) may be expressed as

(52)\[\Delta s = - {{1} \over {nV}} {\left[ \left( Q_{T} + \sum Q_{app} \right) \Delta t + s \left( \sum \alpha \Delta \epsilon {{V_T}\over{4}} - \sum \beta \Delta T {{V_T}\over{4}} \right) \right]}\]

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 = k M\) (\(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 \(t \propto L^2/c\), where \(L\) is to be interpreted as average flow path length, a rough estimate of the second may be obtained from the relation \(t \propto L^2/(k \bar{p}/n)\), where \(L\) is the average height of the medium associated with phreatic storage, and \(\bar 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 \({{\partial\zeta}\over{\partial t}}\) = 0, \({{\partial s}\over{\partial t}}\) = 0 and \({\partial{T}\over{\partial t}}\) = 0. The incremental expression of the fluid constitutive law then takes the form

(53)\[\Delta p = - \alpha M \Delta \epsilon\]

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)\[{{1}\over{3}} \Delta {\sigma}_{ii} + \alpha \Delta p = K \Delta \epsilon\]

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

(55)\[{{1}\over{3}} \Delta {\sigma}_{ii} = (K + {\alpha}^2 M) \Delta \epsilon\]

This apparent tangent bulk modulus, \(K + {\alpha}^2 M\), 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, \({\rho}_d\), of the medium is considered in the input. The saturated density, \(\rho = {\rho}_d + n s{\rho}_f\), is calculated internally using the input value for fluid density, \({\rho}_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 \(\Delta \sigma_{ij}\) to account for the contribution from pore pressure changes (see this equation). The stress increment \(\Delta \sigma_{ij}\) is expressed as

(56)\[\Delta {\sigma}_{ij} = \Delta {\sigma}_{ij}^f + \Delta {\sigma}_{ij}^{th}\]

where

(57)\[\Delta {\sigma}_{ij}^f = - \alpha \overline{\Delta p_{<t>}^n} {\delta}_{ij}\]

is the contribution from flow, and

(58)\[\Delta {\sigma}_{ij}^{th} = \alpha \overline{M} {(\alpha \overline{\Delta \epsilon} - \beta \overline{\Delta T}})_{<t>} {\delta}_{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)\[{1 \over M} {{\partial p} \over {\partial t}} = {{\partial \zeta} \over {\partial t}} - \alpha {{\partial \epsilon} \over {\partial t}}\]

where \(\zeta\) 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, \(\epsilon\) is the volumetric strain, \(M\) is Biot modulus, and \(\alpha\) 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/3 G)}\]

or, if \(\alpha\) = 1, and a combination of high fluid bulk modulus and low porosity is considered,

(61)\[{K_f > 20 (K + 4/3 G)} n\]

\(K_f\) 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):

The stable explicit fluid timestep can be extremely small, because it is inversely proportional to Biot modulus.

Mechanical convergence is very slow, as a consequence of the additional stiffness, \(\alpha^2 M\), due to the effect of the fluid on mass scaling.

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)\[{{\partial \epsilon_{unb}} \over {\partial t}} = 0\]

where the quantity

(63)\[{{\partial \epsilon_{unb}} \over {\partial t}}= {{\partial \zeta} \over {\partial t}} - \alpha {{\partial \epsilon} \over {\partial 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)\[{{\partial \sigma} \over {\partial t}} + \alpha {{\partial p} \over {\partial t}} = K {{\partial \epsilon} \over {\partial t}}\]

In this scheme, the quantity \({{\partial \epsilon_{unb}} \over {\partial t}}\) 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 {{\partial \epsilon_{unb}} \over {\partial t}}\), is interpreted as a pore pressure increment, \(\alpha {{\partial p} \over {\partial t}}\), which is then used to “correct” the (total) stress (see Equation (64) ), using

(65)\[\sigma^n = \sigma - \alpha {{\partial p} \over {\partial t}}\]

where

(66)\[\alpha {{\partial p} \over {\partial t}} = K {{\partial \epsilon_{unb}} \over {\partial t}}\]

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 \(-\alpha {{\partial p} \over {\partial t}} = -K {{\partial \epsilon_{unb}} \over {\partial t}}\) 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, \(V_{unb}\), 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)\[\epsilon_{unb} = {V_{unb,z} \over V_z} - \alpha \Delta\epsilon\]

where \(V_z\) is zone volume, and the unbalanced volume for the zone, \(V_{unb,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)\[{{\partial p} \over {\partial t}} = F_p^{\star} K \epsilon_{unb}\]

where \(F_p^{\star}\) is the relaxation coefficient. Also, the quantity \(- \alpha \Delta\epsilon 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)\[{{\partial p} \over {\partial t}} = F_p^{\star} K {V_{unb} \over V_n}\]

where \(V_n\) is the gridpoint volume.

The relaxation coefficient, \(F_p^{\star}\), 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 \(F_p^{\star}\) 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, \(\hat{V}_{unb}^a\), and maximum unbalanced volume ratio, \(\hat{V}_{unb}^m\), are monitored to detect convergence of the scheme:

(70)\[\hat{V}_{unb}^a = max \biggl[{{ K_n V_{unb}} \over {|\sigma| V_{n}}} \biggr]\]

(71)\[\hat{V}_{unb}^m = {{\sum K_n \sum V_{unb}} \over {|\sigma| \bigl[\sum V_{n}\bigr]^2}}\]

where the summation and maximum are taken over all gridpoints in the model, \(K_n\) is the average bulk modulus of the zones meeting at the gridpoint, and \(|\sigma|\) is the maximum value of the stress norm in the model. The balloon is considered empty when both \(\hat{V}_{unb}^a\) and \(\hat{V}_{unb}^m\) fall below \(10^{-5}\). 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:

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 \(F_p^{\star}K\) in the saturated fast-flow scheme. In fact, the Biot modulus, \(M\), in Equation (59) is replaced by \(F_p^{\star}K\) in the fast-flow scheme (see (69) and (63)).

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.

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:

unconfined fluid flow, in which a phreatic surface, or water table, exists or develops;

models with interfaces;

models with attached grids; and

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

⇄

⇐ FLAC3D Fluid-Thermal-Mechanical-Formulation — Mathematical Description | Calculation Modes for Fluid-Mechanical Interaction ⇒

Was this helpful? ...

Updated: Sep 26, 2024