FLAC3D Theory and Background • FluidMechanical Interaction
Solving FlowOnly and CoupledFlow Problems
FLAC3D has the ability to perform both flowonly and coupled fluidmechanical analyses. Coupled analysis may be performed with any of the builtin mechanical models in FLAC3D.
Several modeling strategies are available to approach the coupled processes.
One of these assumes that pore pressures once assigned to the grid do not change;
this approach does not require any extra memory to be reserved for flow calculation.
The commands associated with this mode are discussed
in Grid Not Configured for Fluid Flow.
Modeling strategies involving flow of fluid require that the model configure fluidflow
command be issued to configure the grid for fluid analysis,
and the model configure fluidflow
command be issued for all zones in which fluid flow can occur.
The commands associated with model configure fluidflow
mode are discussed
in Grid Configured for Fluid Flow.
The different modeling strategies for coupled analysis will be illustrated in the following sections, the more elaborate requiring more computer memory and time. As a general rule, the simplest possible option should be used, consistent with the reproduction of the physical processes that are important to the problem at hand. Recommended guidelines for selecting an approach based on time scales are given in Selection of a Modeling Approach for a Fully Coupled Analysis, and various modeling strategies are described in Fixed Pore Pressure to Coupled Flow and Mechanical Calculations.
Time Scales
When planning a simulation involving fluid flow or coupled flow calculations with FLAC3D, it is often useful to estimate the time scales associated with the different processes involved. Knowledge of the problem time scales and diffusivity help in the assessment of maximum grid extent, minimum zone size, timestep magnitude and general feasibility. Also, if the time scales of the different processes are very different, it may be possible to analyze the problem using a simplified (uncoupled) approach. (This approach for fully coupled analyses is discussed in detail in Selection of a Modeling Approach for a Fully Coupled Analysis.)
Time scales may be appreciated using the definitions of characteristic time given below. These definitions, derived from dimensional analysis, enter the expression of analytical continuous source solutions. They can be used to derive approximate time scales for the FLAC3D analysis.
Characteristic time of the mechanical process
where \(K_u\) is undrained bulk modulus, \(G\) is shear modulus, \(\rho\) is mass density, and \(L_c\) is characteristic length (i.e., the average dimension of the medium).
Characteristic time of the diffusion process
where \(L_c\) is the characteristic length (i.e., the average length of the flow path through the medium), and \(c\) is the diffusivity defined as mobility coefficient \(k\) divided by storativity \(S\):
There are different forms of storativity that apply in FLAC3D depending on the controlling process:
fluid storage
phreatic storage
elastic storage
where \(M\) is Biot modulus, \(\alpha\) is Biot coefficient (\(M\) = \(K_f/n\) and \(\alpha\) = 1 for incompressible grains), \(K\) is drained bulk modulus, \(G\) is shear modulus, \(\rho_w\) is fluid density, \(g\) is gravity and \(L_p\) is characteristic storage length (i.e., the average height of the medium available for fluid storage).
For saturated flowonly calculations (rigid material), \(S\) is fluid storage and \(c\) is the fluid (Equation (3) and (4)):
For unsaturated flow calculations, \(S\) is phreatic storage and the diffusivity (from Equation (3) and (5)) is
For a coupled, saturated, deformationdiffusion analysis with FLAC3D, \(S\) is elastic storage and \(c\) is the true or generalized coefficient of consolidation, defined from Equation (3) and (6) as
There are some noteworthy properties based on the preceding definitions:
Because explicit timesteps in FLAC3D correspond to the shortest time needed for the information to propagate from one gridpoint to the next, the magnitude of the timestep can be estimated using the smallest zone size for \(L_c\) in the formula for characteristic time. It is important to note that the explicit fluid flow timestep in FLAC3D is calculated using the fluid diffusivity (Equation (7)), even in a coupled simulation. The timestep magnitude may thus be estimated from the formula obtained after substitution of Equation (7) in Equation (2), and using the smallest zone size \(L_z\) for \(L_c\):
(10)\[\Delta t = \min \left( {{L_z^2 } \over{k M}} \right)\]In a saturated fluidflow problem, a reduced bulk modulus (i.e., \(M\) or \(K_f\)) leads not only to an increased timestep, but also to an increased time to reach steady state, so that the total number of steps, \(n_t\), stays the same. This number may be estimated by taking the ratio of the characteristic times for the model, \(t_c\), to the critical timestep, \(\Delta t\), using Equation (2), (7) and (10), which gives
(11)\[n_t = {\left( {L_c \over L_z} \right)}^2\]where \(L_c\) and \(L_z\) are characteristic length for the model and the smallest zone.
In a partially saturated fluidflow problem, adjustments can be made to the fluid modulus (\(M\) or \(K_f\)) to speed convergence to steady state. Be careful not to reduce \(M\) (or \(K_f\)) so much that numerical instability results. A condition for stability may be derived from the requirement that the fluid storage (used in the critical timestep evaluation) must remain smaller than the phreatic storage over one zone height, \(L_z\):
(12)\[M > a L_z \rho_w g/n\]where \(a\) is an adjustment factor chosen equal to 0.3.
Using Equation (2), we see that, to avoid any boundary effects in diffusion problems, the characteristic length, \(L_c\), of the model must be larger than the dimension \(\sqrt{c t_s}\), where \(t_s\) is the maximum simulation time and \(c\) is the controlling diffusivity. In turn, the minimum simulation time is controlled by the relation \(t_{min} > L_z^2/c\).
In a coupled flow problem, the true diffusivity is controlled by the stiffness ratio \(R_k\) (i.e., the stiffness of the fluid versus the stiffness of the matrix):
(13)\[R_k = {{\alpha^2 M} \over {K + {4 \over 3} G}}\]With this definition for \(R_k\), Equation (14) may be expressed in two forms:
(14)\[c = k M {1 \over{1 + R_k}}\]and
(15)\[c = {k \over{\alpha^2}} \left( K + {4 \over 3} G \right) {1 \over{1 + {1 \over R_k}}}\]If \(R_k\) is small (compared to 1), Equation (9) shows that FLAC3D’s standard explicit timestep (see Equation (10)) can be considered as representative of the system diffusivity. An orderofmagnitude estimate for the number of steps needed to reach full consolidation under saturated conditions, for instance, can be calculated using Equation (11).
If \(R_k\) is large (i.e., \(M\) is large compared to \((K + 4G/3)/\alpha^2\) (or \(K_f >>> (K + 4G/3)n\))), FLAC3D’s explicit timestep will be very small, and the problem diffusivity will be controlled by the matrix (see Equation (10) and (14)). The value for \(M\) (or \(K_f\)) can be reduced in order to increase the timestep and reach steadystate computationally faster. Equation (14) indicates that if \(M\) (or \(K_f\)) is reduced such that \(R_k\) = 20, then the diffusivity should be within 5% of the diffusivity for infinite \(M\) (or \(K_f\)) (see 1D Consolidation for an example). The time scale is expected to be respected within the same accuracy. Note that, in any case, \(K_f\) should not be made higher than the physical value of the fluid (2 × 10^{9} Pa for water).
Selection of a Modeling Approach for a Fully Coupled Analysis
A fully coupled quasistatic hydromechanical analysis with FLAC3D is often timeconsuming, and even sometimes unnecessary. There are numerous situations in which some level of uncoupling can be performed to simplify the analysis and speed the calculation. The following examples illustrate the implementation of FLAC3D modeling approaches corresponding to different levels of fluid/mechanical coupling. Three main factors can help in the selection of a particular approach:
the ratio between simulation time scale and characteristic time of the diffusion process;
the nature of the imposed perturbation (fluid or mechanical) to the coupled process; and
the ratio of the fluid to solid stiffness.
The expressions for characteristic time, \(t_c^{f}\) (in Equation (2)), diffusivity, \(c\) (in Equation (9)), and the stiffness ratio, \(R_k\) (in Equation (13)) can be used to quantify these factors. These factors are considered in detail below, and a recommended procedure to select a modeling approach based on these factors is given in Recommended Procedure to Select a Modeling Approach.
Time Scale
We first consider the time scale factor by measuring time from the initiation of a perturbation. We define \(t_s\) as the required time scale of the analysis, and \(t_c\) as the characteristic time of the coupled diffusion process (estimate of time to reach steady state, defined using Equation (2) and Equation (9)).
Shortterm behavior
If \(t_s\) is very short compared to the characteristic time, \(t_c\), of the coupled diffusion process, the influence of fluid flow on the simulation results will probably be negligible, and an undrained simulation can be performed with FLAC3D (using either a “dry” or “wet” simulation — see No Flow. The wet simulation is carried out using
model configure fluidflow
andzone fluid active
off or themodel solvestatic
command). No real time will be involved in the numerical simulation (i.e., \(t_s <<< t_c\)), but the pore pressure will change due to volumetric straining if the fluid modulus (\(M\) or \(K_f\)) is given a realistic value. The footing load simulation is an example of this approach.
Longterm behavior
If \(t_s >>> t_c\), and drained behavior prevails at \(t = t_s\), then the porepressure field can be uncoupled from the mechanical field. The steadystate porepressure field can be determined using a fluidonly simulation (
model solvefluid
) (the diffusivity will not be representative), and the mechanical field can be determined next by cycling the model to equilibrium in mechanical mode with theporepressuregeneration
property at the default value of off (model solvestatic
). (Strictly speaking, this engineering approach is only valid for an elastic material, because a plastic material is pathdependent.)
Another way to describe the time scale is in terms of undrained or drained response. Undrained strictly means that there is no exchange of fluid with the outside world (where outside world means outside the sample in a lab test, and other elements in a numerical simulation or a field situation). Drained means that there is a full exchange of fluid with the outside world, which implies that the fluid pressure is able to dissipate everywhere. The words are typically associated with shortterm and longterm, respectively, because an undrained test usually can be done quickly, while a drained test requires a long time for excess fluid pressures to dissipate. In the field, shortterm means that there is insignificant migration of fluid, and longterm means that the pressure has stabilized (which needs a long time).
Note that in simulations conducted outside of the fluid configuration, or within the fluid configuration but with fluid turned off, the total stress adjustments due to an imposed porepressure change (such as those resulting from the lowering of the water table) will not be done internally by the code. The porepressure increments may, however, be monitored (using a FISH function, for instance) and used to decrement the total normal stresses before cycling to mechanical equilibrium. The saturated and unsaturated mass densities will also need to be adjusted if the water table has been moved within the grid, and the simulation is conducted outside of the fluid configuration. (See Verification of Concepts, and Modeling Techniques for Specific Applications for additional information on this topic.)
Nature of Imposed Perturbation to the Coupled Process
The imposed perturbation to a fully coupled hydromechanical system can be due to changes in either the fluid flow boundary condition or the mechanical boundary condition. For example, transient fluid flow to a well located within a confined aquifer is driven by the change in pore pressures at the well. The consolidation of a saturated foundation as a result of the construction of a highway embankment is controlled by the mechanical load applied by the weight of the embankment. If the perturbation is due to change in pore pressures, it is likely that the fluid flow process can be uncoupled from the mechanical process. This is described in more detail below, and illustrated in Transient Fluid Flow to a Well in a Shallow Confined Aquifer. If the perturbation is mechanically driven, the level of uncoupling depends on the fluid versus solid stiffness ratio, as described below.
Stiffness Ratio
The relative stiffness, \(R_k\) (see Equation (13)), has an important influence on the modeling approach used to solve a hydromechanical problem:
Relatively stiff matrix (\(R_k <<< 1\))
If the matrix is very stiff (or the fluid highly compressible) and \(R_k\) is very small, the diffusion equation for the pore pressure can be uncoupled, since the diffusivity is controlled by the fluid (Detournay and Cheng 1993). The modeling technique will depend on the driving mechanism (fluid or mechanical perturbation):
In mechanically driven simulations, the pore pressure may be assumed to remain constant. In an elastic simulation, the solid behaves as if there were no fluid; in a plastic analysis, the presence of the pore pressure may affect failure. This modeling approach is adopted in slope stability analyses (see Fixed Pore Pressure).
In pore pressuredriven elastic simulations (e.g., settlement caused by fluid extraction), volumetric strains will not significantly affect the porepressure field, and the flow calculation can be performed independently (
model solvefluid
). (In this case, the diffusivity will be accurate because for \(R_k <<< 1\), the generalized consolidation coefficient in Equation (9) is comparable to the fluid diffusivity in Equation (7).) In general, the porepressure changes will affect the strains, and this effect can be studied by subsequently cycling the model to equilibrium in mechanical mode (model solvestatic
). The fluid modulus (\(M\) or \(K_f\)) must be set to zero during mechanical cycling, to prevent additional generation of pore pressure.
Relatively soft matrix (\(R_k >>> 1\))
If the matrix is very soft (or the fluid incompressible) and \(R_k\) is very large, then the system is coupled, with a diffusivity governed by the matrix. The modeling approach will also depend on the driving mechanism.
In mechanically driven simulations, calculations can be timeconsuming. As discussed in note 5 of Time Scale, it may be possible to reduce the value for \(M\) (or \(K_f\)), such that \(R_k\) = 20, and not significantly affect the response. This can be done automatically by using the
zone fluid modulusautomatic
command.In most practical cases of pore pressuredriven systems, experience shows that the coupling between pore pressure and mechanical fields is weak. If the medium is elastic, the numerical simulation can be performed with the flow calculation in flowonly mode (
model solvefluid
) and then in mechanicalonly mode (model solvestatic
) to bring the model to equilibrium. This can be done easily in multiple stages using themodel solvefluiddecoupled
) command.It is important to note that, in order to preserve the true diffusivity (and hence the characteristic time scale) of the system, the fluid modulus \(M\) (or \(K_f\)) must be adjusted to the value
or
during the flow calculation (see Equation (9)), and to make sure the porepressuregeneration
property is off to prevent further adjustments by volumetric strains (Berchenko 1998).
This scaling can be done most easily using the zone fluid modulusscale
command.
Recommended Procedure to Select a Modeling Approach
It is recommended that the selection of a modeling approach for a fully coupled analysis follow the procedure outlined in Table 1. First, determine the characteristic time of the diffusion process for the specific problem conditions and properties (see Time Scale), and compare this time to the actual time scale of interest. Second, consider whether the perturbation to the system is primarily porepressure driven or mechanically driven. Third, determine the ratio of the stiffness of the fluid to the stiffness of the solid matrix. Table 1 indicates the appropriate modeling approach based on the evaluation of these three factors. The table also indicates the required adjustment to the fluid modulus, \(M_a\) (or \(K_f^a\)), for each case. Finally, the table lists several examples from the manual that illustrate each modeling approach.
Time Scale 
Imposed Process Perturbation 
Fluid vs Solid Stiffness 
Modeling Approach & Main Calculation Commands 
Adjust Fluid Bulk Modulus (\(M^a\) or \(K_f^a\)) 
Examples (6) 

\(t_s >>> t_c\) 
mechanical or pore pressure 
any \(R_k\) 
Effective Stress (1) 
no fluid 

(steadystate analysis) 
with no fluid flow 

or 

Effective Stress (2) 

\(t_s <<< t_c\) 
mechanical or pore pressure 
any \(R_k\) 
Pore Pressure Generation (3) 
realistic value for \(M^a\) (or \(K_f^a\)) 

(undrained analysis) 



\(t_s\) in the range of \(t_c\) 
pore pressure 
any \(R_k\) 
Uncoupled FlowMechanical (4) 

step 1 

step 2 

\(t_s\) in the range of \(t_c\) 
mechanical 
any \(R_k\) 
Coupled FlowMechanical (5) 

Notes to Table 1:
The effective stress approach with no fluid flow is discussed in Fixed Pore Pressure. In order to establish the initial conditions for this effective stress analysis, use the
zone gridpoint porepressure
command or a FISH function to establish steadystate pore pressures. Specify the correct wet density to zones below the water table, and dry density to zones above.The effective stress approach with groundwater flow is discussed in FlowOnly Calculation. In order to establish the initial conditions for this effective stress analysis, use the
zone gridpoint initialize porepressure
command or a FISH function to establish steadystate pore pressures, or usezone fluid steadystate
or specifymodel solvefluid
and step to steady state if the location of the phreatic surface is not known. Set \(M^a\) (or \(K_f^a\)) to a small value to speed convergence for a partially saturated system. Note that \(M^a\) (or \(K_f^a\)) should be greater than \(0.3nL_{z}\rho_{w}g\) (or \(0.3L_{z}\rho_{w}g\)) to satisfy numerical stability (see Equation (12)).The porepressure generation approach is discussed in No Flow — Mechanical Generation of Pore Pressure. In order to establish the initial conditions for the porepressure generation analysis, use the
zone gridpoint initialize porepressure
command or a FISH function to establish steadystate flow, or use thezone fluid steadystate
command ormodel solvefluid
and step to steady state if the location of the phreatic surface is not known. Set \(M^a\) (or \(K_f^a\)) to a small value to speed convergence for a partially saturated system. Note that \(M^a\) (or \(K_f^a\)) should be greater than \(0.3nL_{z}\rho_{w}g\) (or \(0.3L_{z}\rho_{w}g\)) to satisfy numerical stability (see Equation (12)).The uncoupled fluidmechanical approach is described in Selection of a Modeling Approach for a Fully Coupled Analysis. This approach is recommended for pore pressuredriven systems, and should be used carefully if \(R_k >>> 1\). Note that the adjusted value for \(M^a\) (or \(K_f^a\)) during the flowonly step should satisfy Equation (16) so that the coupled diffusivity will be correct  for this purpose use the
zone fluid modulusscale
command.The fully coupled approach is discussed in Coupled Flow and Mechanical Calculations. Note that for \(R_k >>> 1\), if \(M^a\) (or \(K_f^a\)) is adjusted to reduce \(R_k\) = 20, easiest done by using the
zone fluid modulusautomatic
command, the time response will be close (typically within 5%) to that for infinite \(M\) (or \(K_f\)).The saturated fastflow option for fully coupled analysis (
zone fluid fastflow
on) is discussed in Fully Saturated Fast Flow. Note that this method can only be applied for fully saturated flow. ther limitations are discussed in Fully Saturated Fast Flow.Examples are given in Verification Examples on Fluid. Example Applications (E.A.) and Verification Problems (V.P.) that demonstrate the various methods are in the Example Problems and Verification Problems.
Fixed Pore Pressure (Used in Effective Stress Calculation)
In some calculations, the porepressure distribution is important only because it is used in the computation of effective stress at all points in the system. For example, in modeling slope stability, we may be given a fixed water table. To represent this system with FLAC3D, it is sufficient to specify a porepressure distribution that is unaffected by mechanical deformations that may occur. Because no change in pore pressures is involved, we do not need to configure the grid for groundwater flow.
The zone gridpoint porepressure geometry
command can be used to specify an initial hydrostatic porepressure distribution below a given fixed phreatic surface.
The water density must be provided (with the fluiddensity
keyword or the zone fluiddensity
command),
and appropriate dry and saturated material densities supplied by the user above and below the water table, respectively.
Alternatively, one of the other options of the zone gridpoint porepressure
command or a FISH function can be used
to generate the required static porepressure distribution.
The porepressure distribution corresponds to an initial state for which there is no strain. It remains constant and is unaffected by mechanical deformation. Fluid flow does not take place. The influence of this porepressure distribution is on failure in material for which yield depends on mean effective stress. For an example application, see Curvature Slope Stability.
FlowOnly Calculation to Establish a PorePressure Distribution
Flowonly calculations may be performed to determine the flow and pressure distribution in some system independent of any mechanical effects. For example, it may be necessary to evaluate the groundwater changes that result from the digging of a drainage ditch or the activation of a pumping well. In other instances, an initial porepressure distribution may be needed for a coupled calculation. In both cases, FLAC3D may be run in the fluidflow mode without any mechanical calculation being done. Mechanical calculations may or may not be done subsequently.
The first step in the command procedure for a flowonly calculation is to issue a model configure fluidflow
command
so that extra memory can be assigned for the fluidflow calculation.
Then a choice must be made between the explicit and implicit fluidflow solution algorithm.
By default, the explicit algorithm will be selected,
but the implicit mode of calculation may be activated (and deactivated) at any stage of the calculation
using the zone fluid implicit
on or zone fluid implicit
off command.
You can also use the zone fluid implicit servo
on command, which will cause FLAC3D to attempt
to automatically switch from explicit to implicit when it is computationally advantageous.
The implicit calculation mode may fail to converge in nonfully saturated problems.
In the explicit mode, the fluidflow timestep will be calculated automatically,
but a smaller timestep can be selected using the model fluid timestep
command.
The magnitude of the timestep must be specified by the user in the implicit mode, unless the implicit servo is activated.
This is done by issuing a model fluid timestep
command.
For saturated flow, it is often more efficient to use the implicit solution mode when contrasting permeabilities exist.
The fluidflow model and properties must be specified for all zones in which fluid flow may occur.
Initial and boundary conditions are assigned to complete the fluidflow problem setup.
The fluidflow domain in a fluidonly or fluidmechanical simulation is defined by the assembly of zones with a nonnull fluid flow model.
Flux boundary conditions, for instance, can be assigned by specifying the zone face apply flux
command
with the range keyword to correspond to the boundaries of that domain.
The model solvefluid
command may be specified to execute a given number of fluidflow steps or to stop the calculation automatically when a particular fluidflow time is reached using the time
, timetotal
, or cycles
keywords.
A steadystate flow condition can be calculated in one step by using the zone fluid steadystate
command.
This may fail to converge for an unsaturated model, in which case a steadystate flow condition can also be calculated by using the model solvefluid
command with the ratioflow keyword to specify the limiting unbalanced fluidflow ratio defining the steady flow state.
If the computed porepressure distribution is then to be used in a mechanical calculation
where pore pressure can be assumed to remain constant,
the model solvestatic
commands should be used.
The porepressuregeneration
property should be left off to prevent extra pore pressures from being generated by mechanical deformation.
3D SteadyState Fluid Flow with a Free Surface and 2D SteadyState Fluid Flow with a Free Surface present example flowonly analyses in which unconfined flow through a vertical embankment is calculated. The analysis produces the steadystate phreatic surface. Note that in this example, the fluid bulk modulus is set to a low value to enable the phreatic surface to develop quickly. This can be done when the time taken to reach steady state is unimportant. Note, however, that there is a lower limit for fluid modulus in order to avoid numerical instability (see Equation (12)).
No Flow — Mechanical Generation of Pore Pressure
The undrained (shortterm) response may be analyzed in FLAC3D using both dry and wet approaches.
In a dry simulation, the generation of pore pressure under volumetric strain is not simulated directly.
Instead, its effect on mechanical deformation is taken into account by assigning the undrained bulk value,
\(K_u = K + \alpha^2 M\), to the material bulk modulus, \(K\), in the FLAC3D simulation.
In this case, two different techniques can be applied to detect failure in a MohrCoulomb material.
In the first one, the constant pore pressures are initialized using the zone gridpoint porepressure
command,
and undrained cohesion and friction are given as input.
In the second, the material is assigned a zero friction and a cohesion equal to the undrained shear strength, \(C_u\).
The first technique applies to problems where changes in pore pressure are small compared to the initial values.
The second technique strictly applies to planestrain problems with Skempton coefficient equal to one (\(M >>> K + 4G/3\)).
Note that a dry simulation can be carried out whether or not the model configure fluidflow
command has been issued.
However, if the command has been used, the porepressuregeneration
property must be set off (the default) to prevent additional generation of pore pressure.
In a wet simulation, the shortterm response of a coupled system is analyzed in the fluid configuration of FLAC3D. In this case, drained values should be used for the material bulk modulus, friction and cohesion. If fluid calculations are inactive, the Biot modulus (or fluid modulus) is given a realistic value, and the ``porepressuregeneration` property is on, then pore pressure will be generated as a result of mechanical deformations. For example, the “instantaneous” pore pressures produced by a footing load can be computed in this way. If the fluid bulk modulus (\(M\) or \(K_f\)) is much greater than the solid bulk modulus, convergence will be slow; it may be possible to reduce the fluid bulk modulus without significantly affecting the behavior. (See Note 5 in Table 1.) The data file FootingInstantPP.dat illustrates pore pressure buildup produced by a footing load on an elastic/plastic material contained in a box. To view this project, use the menu command . This corresponds to the “Footing.prj” project. The left boundary of the box is a line of symmetry, and the pore pressure is fixed at zero along the top surface to prevent porepressure generation there. By default, the porosity is 0.5; permeability is not needed, since flow is not calculated. Note that the pore pressures are fixed at zero at gridpoints along the top of the grid. This is done because at the next stage of this model, a coupled, drained analysis in which drainage will be allowed at the ground surface will be performed (see Coupled Flow and Mechanical Calculations). The zero porepressure condition is set now to provide the compatible porepressure distribution for the second stage.
As a large amount of plastic flow occurs during loading, the normal stress is applied gradually by using the ramp
keyword in the zone apply
command.
Figure 1 shows porepressure contours and vectors representing the applied forces.
It is important to realize that the plastic flow will occur in reality over a very short period of time (on the order of seconds);
the word “flow” here is misleading since, compared to fluid flow, it occurs instantaneously.
Hence, the undrained analysis (using model solvestatic
) is realistic.
The saturated fastflow logic can be used for this example because the material is fully saturated, and the stiffness of the fluid is significantly greater than that of the solid matrix.
It is only necessary to include the zone fluid fastflow
on command in FootingInstantPP.dat to perform a saturated fastflow calculation.
Calculational stepping stops, and equilibrium is achieved, when the unbalanced volumes, \(V_{av}\) and \(V_{max}\), and the unbalanced force ratio are smaller than the predefined limits.
The instantaneous porepressure generation is nearly the same (within approximately 5%) as that obtained by the basic flow scheme, as illustrated in Figure 2.
The runtime for the saturated fastflow scheme is roughly three times faster than that for the basic flow scheme.
FootingInstantPP.dat:
model new
model title 'Instantaneous pore pressures generated under an applied load'
zone create brick size 40 1 20 point 0 (0,0,0) point 1 (20,0, 0) ...
point 2 (0,1,0) point 3 ( 0,0,10)
zone face skin
;  mechanical model 
model largestrain off
zone cmodel assign mohrcoulomb
zone property bulk 5e7 shear 3e7 friction 25 cohesion 1e5 tension 1e10
zone face apply velocitynormal 0 range group 'Top' not
;  apply load slowly 
zone face apply stressnormal 0.3e6 ...
servo ramp lowerbound 2e3 ratio convergence ...
range positionx 0 3 group 'Top'
;  fluid flow model 
model configure fluidflow
zone fluid property porepressuregeneration on
zone fluid property fluidmodulus 9e8
;  pore pressure fixed at zero at the surface 
zone face apply porepressure 0 range group 'Top'
;  histories 
zone history porepressure position (2,1,9)
;  timing 
model solvestatic
model save 'load'
program return
Coupled Flow and Mechanical Calculations
FLAC3D can do a coupled fluidflow and mechanical calculation if the grid is configured for fluid,
and if the Biot modulus, or fluid bulk modulus, and permeability are set to realistic values, and if the porepressuregeneration property
is set to on.
The full fluidmechanical coupling in FLAC3D occurs in two directions: porepressure changes cause volumetric strains that influence the stresses;
in turn, the pore pressure is affected by the straining that takes place.
The relative time scales associated with consolidation and mechanical loading should be appreciated. Mechanical effects occur almost instantaneously (on the order of seconds, or fractions of seconds). However, fluid flow is a long process: the dissipation associated with consolidation takes place over hours, days or weeks.
Relative time scales may be estimated by considering the ratios of characteristic times for the coupled and undrained processes. The characteristic time associated with the undrained mechanical process is found by using the saturated mass density for \(\rho\) and the undrained bulk modulus, \(K_u\), as defined in this equation (or this equation for incompressible grains) for \(K\) in Equation (1). The ratio of the diffusion characteristic time (Equation (2) and (9)) and the undrained mechanical characteristic time is thus
In most cases, \(M\) is approximately 10^{10} Pa, but the mobility coefficient \(k\) may differ by several orders of magnitude; typical values are
10^{19}m^{2}/Pasec for granite;
10^{17}m^{2}/Pasec for limestone;
10^{15}m^{2}/Pasec for sandstone;
10^{13}m^{2}/Pasec for clay; and
10^{7}m^{2}/Pasec for sand.
For rock and soil, \(\rho\) is of the order of 10^{3}kg/m^{3}, while \(K + 4/3 G\) is approximately 10^{10} Pa. Using those orders of magnitude in Equation (18), we see that the ratio of fluidtomechanical time scales may vary between \(L_c\) for sand, 10^{6}\(L_c\) for clay, 10^{8}\(L_c\) for sandstone, 10^{10} \(L_c\) for limestone and 10^{12} \(L_c\) for granite. If we exclude materials with mobility coefficients larger than that of clay, we see that this ratio remains very large, even for small values of \(L_c\).
In practice, mechanical effects can then be assumed to occur instantaneously when compared to diffusion effects. This is also the approach adopted in FLAC3D, where no time is associated with any of the mechanical steps taken together with the fluidflow steps. The use of the dynamic option in FLAC3D may be considered, to study the fluidmechanical interaction in materials such as sand, where mechanical and fluid time scales are comparable.
In most modeling situations, the initial mechanical conditions correspond to a state of equilibrium that must first be achieved before the coupled analysis is started.
Typically, at small fluidflow time (compared to \(t_c\) — see Time Scale),
a certain number of mechanical steps must be taken for each fluid step to reach quasistatic equilibrium. At larger fluidflow time,
if the system approaches steadystate flow,
several fluidflow timesteps may be taken without significantly disturbing the mechanical state of the medium.
A corresponding numerical simulation may be controlled manually by alternating
between flowonly (model solvefluid
) and mechanicalonly (model solvestatic
) modes.
The model solve
command can be used to perform calculations for both the flowonly and mechanicalonly processes.
As an alternative approach, the tedious task of switching between flowonly and mechanicalonly modes may be avoided by using the model solvefluidcoupled
command in combination with appropriate settings;
the time
keyword to specific the incremental fluid time to calculate, and the convergence
keyword to set the convergence limit used after each fluid step to return to quasiequilibrium.
In a third approach, the model step
command may be used while both mechanical and fluidflow modules are on.
In this option, one mechanical step will be taken for each fluidflow step.
Here, fluidflow steps are assumed to be so small that one mechanical step is enough.
Input Instructions should be consulted for a complete list of available command options for a coupled analysis.
To illustrate a fully coupled analysis, we continue the footing simulation
(done in No Flow — Mechanical Generation of Pore Pressure),
using the model solvefluidcoupled
command to calculate the consolidation beneath the footing.
The data file is FootingCoupled.dat:
FootingCoupled.dat:
model restore 'load'
zone gridpoint initialize velocity (0,0,0)
zone gridpoint initialize displacement (0,0,0)
;  histories 
model history fluid timetotal
zone history displacementz position 0,1,10
zone history displacementz position 1,1,10
zone history displacementz position 2,1,10
zone history porepressure position 2,1,9
zone history porepressure position 5,1,5
zone history porepressure position 10,1,7
zone gridpoint initialize ratiotarget 1e3 range positionx 4 21
;  turn on fluid flow model 
zone fluid property mobilitycoefficient 1e12
;  solve to 3,000,000 sec 
model solvefluidcoupled time 3e6
model save 'age_3e6'
program return
The cycling output should be watched during the calculation process.
The mechanical convergence used is 1.0 (the default) for this problem;
enough mechanical steps are taken to keep the convergence ratio below this tolerance.
However, the limit to mechanical steps, defined by model mechanical substep
, is set to 1000 (the default) in this example.
If the actual number of mechanical steps taken is always equal to the set value of model mechanical substep
,
then something must be wrong.
Either the force limit or model mechanical substep
has been set too low, or the system is unstable and cannot reach equilibrium.
The quality of the solution depends on the convergence used:
a small convergence will give a smooth, accurate response,
but the run will be slow;
a large convergence will give a quick answer, but it will be noisy.
The characteristic diffusion time for this coupled analysis is evaluated from Equation (2), using Equation (9) for the diffusivity and a value of \(L_c\) = 10 m corresponding to the model height. Using the property values in this example, \(t_c^{f}\) is estimated at 1,200,000 seconds. Full consolidation is expected to be reached within this time scale; the numerical simulation is carried out for a total of 3,000,000 seconds.
Figure 3 shows the time histories to 3,000,000 seconds of displacements under the footing load. In this simulation, pore pressures remain fixed at zero on the ground surface; hence, the excess fluid escapes upward. The leveling off of the histories indicates that full consolidation has been reached.
The coupled simulation can also be run using the saturated fastflow scheme because the foundation material is fully saturated.
In this case the zone fluid fastflow
on command is invoked.
When the fluid flow is turned on, the coupled calculation will continue using the saturated fastflow scheme.
The results are similar to the results for the basic flow case, as indicated
by Figure 4 (compare to Figure 3).
The unbalanced force ratio tolerance is set to a convergence of 0.1 for this simulation, in order to provide a smoother displacement history.
Even with this lower tolerance, the saturated fastflow run is approximately 10 times faster than the basic flow simulation.
Finally, a comparison is made with an uncoupled drained analysis.
First, a flow only calculation is made, then a mechanical only calculation is performed, with porepressuregeneration
set to off.
The data file for the uncoupled run is FootingUncoupled.dat.
The displacement histories for the uncoupled drained simulation are plotted in Figure 5.
The final settlement match within 3% of the basic flow results.
FootingUncoupled.dat:
model restore 'load'
;  turn on fluid flow model 
zone fluid property mobilitycoefficient 1e12
zone gridpoint initialize velocity (0,0,0)
zone gridpoint initialize displacement (0,0,0)
;  histories 
model history fluid timetotal
zone history displacementz position 0,1,10
zone history displacementz position 1,1,10
zone history displacementz position 2,1,10
zone history porepressure position 2,1,9
zone history porepressure position 5,1,5
zone history porepressure position 10,1,7
;  solve to 3,000,000 sec 
;  uncoupled 
model solvefluid time 3e6
zone fluid property porepressuregeneration off
model solvestatic
model save 'age_3e6_uncoup'
program return
If a sudden change of loading or mechanical boundary condition is applied in a coupled problem,
it is important to allow the undrained (shortterm) response to develop before allowing flow to take place.
In other words, FLAC3D should be run to equilibrium under model solvestatic
following the imposed mechanical change.
The model solvefluidcoupled
command can then be used to compute the subsequent coupled flow/mechanical response.
If changes in fluid boundary conditions occur physically at the same time as mechanical changes,
then the same sequence should be followed (i.e., mechanical changes  equilibrium  fluid changes  coupled solution).
Another example of fully coupled behavior is the timedependent swelling that takes place following the excavation of a trench in saturated soil.
In this case, negative pore pressures build up immediately after the trench is excavated;
the subsequent swelling is caused by the gradual influx of water into the region of negative pressures.
We model the system in two stages:
in the first, we allow mechanical equilibrium to occur, without flow; then we allow flow,
using the model solvefluidcoupled
command to maintain quasistatic equilibrium during the consolidation process.
The fluid tension is initialized to a large negative number to prevent desaturation.
The data file used is Swelling.dat.
To view this project, use the menu command . This corresponds to the “Swelling.prj” project.
The trench is excavated in the lefthand part of a flat soil deposit that is initially fully saturated and in equilibrium under gravity. The material is elastic in this case, but it could equally well have been a cohesive material, such as clay. In this run, we assume impermeable conditions for the free surfaces. Figure 6 shows the displacement vectors that accumulate during the time that flow is occurring; the trench is seen at the lefthand side of the model. Figure 7 shows the time history of pore pressure near the crest of the trench; note that there is an initial negative excursion in pressure arising from the instantaneous expansion of the soil toward the trench. Figure 8 shows histories of horizontal and vertical displacement at the crest. The characteristic time for this problem, evaluated using the model length of 40 m for \(L_c\), is approximately 5 × 10^{8} seconds (based on Equation (2) and (9)); the numerical simulation is carried out to that time.
Swelling.dat:
model new
model title "Maintaining equilibrium under timedependent swelling conditions"
; Geometry
zone create brick size 40 1 8
zone face skin
zone group 'Remove' range positionx 0,2 positionz 2,8
;  mechanical model 
zone cmodel assign elastic
zone property bulk 2e8 shear 1e8 density 1500
;  fluid flow model 
model configure fluidflow
zone fluid property mobilitycoefficient 1e14 fluidmodulus 2e9 ...
fluiddensity 1000
zone fluid property porepressuregeneration on ; 2 way coupling
; initial conditions
model gravity 10
zone initializestresses
; Boundary conditions
zone face apply velocitynormal 0 range group 'North' or 'South' or 'East' or 'West'
zone face apply velocity (0,0,0) range group 'Bottom'
zone gridpoint porepressure head 8 fluiddensity 1000
;;; Null region
zone null range group 'Remove'
;  settings 
model largestrain off
; Histories
model history mechanical convergence
; Solve
model solvestatic
model save 'swell1'
zone gridpoint initialize displacement (0,0,0)
model history fluid timetotal
zone history porepressure position 4.5,0.5,6.5
zone history displacementx position 2,0,8
zone history displacementz position 2,0,8
zone gridpoint fix porepressure range group 'East'
zone fluid implicit servo on
model solvefluidcoupled time 5e8 convergence 20
model save 'swell2'
program return
break
Was this helpful? ...  Itasca Software © 2024, Itasca  Updated: Sep 26, 2024 