FLAC3D Theory and Background • FluidMechanical Interaction
OneDimensional Consolidation
Note
The project file for this example is available to be viewed/run in FLAC3D. The project’s main data files are shown at the end of this example.
A saturated layer of soil of thickness \(H\) = 20 m (shown in Figure 1) and large horizontal extent rests on a rigid impervious base. A constant surface load, \(p_z\) = 10^{5} Pa, is applied on the layer under undrained conditions. The soil matrix is homogeneous and behaves elastically; the isotropic Darcy’s transport law applies. The applied pressure is initially carried by the fluid, but as time goes on the fluid drains through the layer surface, transferring the load to the soil matrix. The solution to this onedimensional consolidation problem may be expressed in the framework of Biot theory (see Detournay and Cheng 1993).
The diffusion equation for the pore pressure, \(p\), has the form
where: 
\(c\) 
= 
\(k/S\); 
\(S\) 
= 
\(1/M + \alpha^2/\alpha_1\) = the storage coefficient; 

\(\alpha_1\) 
= 
\(K + 4G/3\); 

\(M\) 
= 
the Biot modulus; and 

\(\alpha\) 
= 
the Biot coefficient. 
The boundary conditions are \(\sigma_{zz}\) = \(p_z H(t)\) and \(p\) = 0 at \(z\) = H (\(H(t)\) is the step function), and \(u_z\) = 0 and \({{\partial p}\over{\partial z}}\) = 0 at \(z\) = 0.
Because the stress is constant, Equation (1) reduces to
with boundary conditions \(p\) = 0 at \(z = H\), and \({{\partial p}\over{\partial z}}\) = 0 at \(z\) = 0.
The initial value, \(p_0\), for the pore pressure induced from loading of the layer may be derived from the fluid constitutive law (\(p = M(\zeta  \alpha {\epsilon}_{zz})\)) by considering undrained conditions (\(\zeta\) = 0) and using the onedimensional mechanical constitutive law (\({\sigma}_{zz} = \alpha_1 {\epsilon}_{zz}  \alpha p\)) to express strain in terms of stresses. It is given by
The solution for the pore pressure is
where: 
\(\hat{p}\) 
= 
\(p/p_z\); 
\(a_m\) 
= 
\((2 m +1)\pi/2\); 

\(\hat{z}\) 
= 
\((Hz)/H\); and 

\(\hat{t}\) 
= 
\(ct/H^2\). 
The vertical displacement, \(u_z\), is found by considering the equilibrium equation \(\partial {\sigma}_{zz} / \partial z\) = 0 together with the mechanical constitutive equation \({\sigma}_{zz} = \alpha_1 {\epsilon}_{zz}  \alpha p\). By expressing \({\epsilon}_{zz}\) as \(\partial u_z / \partial z\), we obtain upon integration, taking Equation (4) and the boundary conditions into account,
here \({\hat{u}}_z\) is defined as \({\hat{u}}_z = {{\alpha_1 u_z}\over{H p_z}}\).
The following properties are prescribed for this example:
dry bulk modulus, \(K\) 
5 × 10^{8} Pa 
dry shear modulus, \(G\) 
2 × 10^{8} Pa 
Biot modulus, \(M\) 
4 × 10^{9} Pa 
Biot coefficient, \(\alpha\) 
1.0 
permeability, \(k\) 
10^{10} m^{2}/Pasec 
The FLAC3D model grid for this problem is a column of 20 zones of unit dimensions lined up along the \(z\)axis. (See Figure 2.) The base of the column is fixed, and lateral displacements are restricted in the \(x\) and \(y\)directions. A mechanical pressure, \(p_z\), is applied at the top of the column.
At first, flow is prevented and the model is stepped to equilibrium to establish the initial undrained conditions. At the end of this stage, the stress \({\sigma}_{zz}\) has the value \(p_z\), in equilibrium with the applied pressure, and the pore pressure has the initial undrained value \(p_0 = \alpha p_z/\alpha_1 S\) (see Figure 3). Drainage is then allowed by setting the pore pressure to zero at the top of the column. The FLAC3D model is cycled to a flow time of 5000 seconds, which is approximately the magnitude of the characteristic time, \(t_c = H^2/c\), for this problem.
The analytical solution for the pore pressure, \(\hat{p}\), and vertical displacement, \({\hat{u}}_z\), at the column midheight are evaluated using FISH functions, and compared to the numerical solution as time proceeds. The results are plotted versus fluidflow time in Figure 3 and Figure 4. The transfer of pore pressure to effective stress is illustrated in Figure 5, where the evolution of normalized total stress, \({\sigma}_{zz} / p_z\), effective stress, \(({\sigma}_{zz} + p) / p_z\), and pore pressure, \(p / p_z\), with fluidflow time are presented. The FLAC3D data file is listed in “1DConsolidationCoupled.dat”.
break
Note that as the stiffness ratio \(R_k = \alpha^2 M/(K + 4G/3)\) increases, the number of zones should
increase to keep the error small. This effect is caused by the numerical technique used to evaluate
porepressure changes caused by volumetric strain in model configure fluidflow
. Also, as the ratio
\(R_k\) increases, the timestep will decrease. However, to keep the error small and the timestep
sufficiently large without significantly affecting the solution, it is acceptable to limit the value of
Biot modulus (or \(K_f\)) used in the simulation to a value of approximately twenty times
\((K + 4G/3)/{\alpha^2}\) (see the fifth item in Solving FlowOnly and CoupledFlow Problems, and the zone fluid modulusautomatic
command). To demonstrate the validity of this numerical approach, the analytical solutions for pore pressure and displacement corresponding to \(M\)
= 1.5 × 10^{10} (\(R_k \simeq\) 20) and \(M\) = 1.5 × 10^{15} (\(R_k \simeq\) 20 ×
10^{5}) are compared in Figure 6 and
Figure 7. As seen in those plots, the responses are indeed similar
(maximum relative error less than 5%).
break
For values of the stiffness ratio, \(R_k\), much larger than 20, the convergence of the basic algorithm for coupled fluid flow will be very slow (i.e., numerous mechanical steps will be necessary to bring the system to mechanical equilibrium after each flow step). This is illustrated by repeating the couple test with Biot modulus \(M\) = 4 × 10^{10} Pa (10 times larger than in the former numerical calculation, see 1DConsolidationCoupled10M.dat). In the former calculation with \(M\) = 4 × 10^{9} Pa (\(R_k\) = 5.2), the coupled flow calculation ran in approximately 4 seconds (on a 1.6 GHz i7 computer). For the calculation with \(M\) = 4 × 10^{10} Pa (\(R_k\) = 52), the calculation time is approximately 40 seconds (ten times slower). The comparison of the pore pressure and displacement history results to the analytical solutions for \(M\) = 4 × 10^{10} Pa are shown in Figure 8 and Figure 9, respectively. The maximum relative error is less than 4%.
break
The solution time can be reduced by using the uncoupling technique described in Stiffness Ratio (the applicability of the uncoupling technique to this problem follows from the irrational character of the solution). In this test (1DConsolidationUncoupled.dat), the numerical simulation is repeated for a Biot modulus of 4 × 10^{10} Pa using the uncoupling technique. In this test, the initial undrained conditions are obtained using the undrained bulk value for the material and setting Biot modulus to zero. The pore pressure, which is not updated in this calculation mode, is then initialized at the value \(p_0\). The rest of the simulation is carried out in a series of ten time increments to enable a recording of the variables history. In each increment, the fluid flow calculation is performed for a time interval of 500 seconds. During that stage, a scaled value for Biot modulus is used in order to preserve the coupled system true diffusivity (see this equation). Next, Biot modulus is set to zero to prevent additional generation of pore pressure, and the system is run to mechanical equilibrium using the drained value of the bulk modulus. Fluid and mechanical calculations are repeated until the total simulation time reaches 5000 sec. The results of this simulation (which does not involve any direct calculation of porepressure change due to volumetric straining) are presented in Figure 10 to Figure 12. The calculation time for this uncoupled simulation is approximately 4 seconds. The maximum relative error is less than 1%.
break
The uncoupled method switches calculations between fluid and matrix mechanical using the model solvefluiddecoupled
command.
The results using the above simplified uncoupled method are presented in Figure Figure #fig1dconsolidationsimplifiedpp
to Figure Figure #fig1dconsolidationsimplifiedevolution
.
break
This particular example involves fully saturated coupled flow.
In such cases, solution times can be reduced by using the saturated fastflow algorithm described in Fully Saturated Fast Flow.
The algorithm is selected by using the zone fluid fastflow
on command.
If 1DConsolidationCoupled.dat with \(M\) = 4 × 10^{10} Pa (\(R_k\) = 52) is repeated with this command,
the runtime is now approximately 12 seconds (see 1DConsolidationFastFlow.dat).
As \(R_k\) increases, the use of the fastflow algorithm becomes more beneficial.
The results of the fastflow calculation are shown in Figure 13 and Figure 14. The maximum relative error is less than 2%.
break
Reference
Detournay, E., and A. H. D. Cheng. “Fundamentals of Poroelasticity,” in Comprehensive Rock Engineering, Vol. 2, pp. 113171. J. Hudson et al., eds. London: Pergamon Press (1993).
Data File
1DConsolidationCoupled.dat
model new
model title "Onedimensional consolidation (coupled)"
;  model geometry 
zone create brick size 1 1 20 ; point 1 (0.5,0,0) point 2 (0,0.5,0) point 3 (0,0,10)
zone face skin
;  mechanical model 
model largestrain off
zone cmodel assign elastic
zone property bulk 5e8 shear 2e8
zone face apply velocitynormal 0.0 range group 'Top' not
zone face apply stressnormal 1e5 range group 'Top'
;  fluid flow properties 
model configure fluidflow
zone fluid property biotmodulus 4e9
zone fluid property porepressuregeneration on ; activate two way coupling.
;  first establish undrained response 
model solvestatic
;  histories 
history interval 10
zone history name="pp" porepressure position (0,0,10)
zone history name="szz" stresszz position (0.5,0.5,10.5)
zone history name="zdisp" displacementz position (0,0,10)
model history name="time" fluid timetotal
;  drained response 
zone face apply porepressure 0 range group 'Top'
; Loosely coupled till time 5000
zone fluid property mobilitycoefficient 1e10
zone fluid implicit servo on
model solvefluidcoupled time 5000 convergence 1
model save '1dconsolidation_coupled'
program return
1DConsolidationCoupled10M.dat
model new
model title "Onedimensional consolidation (coupled)"
;  model geometry 
zone create brick size 1 1 20
zone face skin
;  mechanical model 
model largestrain off
zone cmodel assign elastic
zone property bulk 5e8 shear 2e8
zone face apply velocitynormal 0.0 range group 'Top' not
zone face apply stressnormal 1e5 range group 'Top'
;  fluid flow properties 
model configure fluidflow
zone fluid property biotmodulus 4e10
zone fluid property porepressuregeneration on ; activate two way coupling.
;  first establish undrained response 
model solvestatic
;  histories 
history interval 100
zone history name="pp" porepressure position (0,0,10)
zone history name="szz" stresszz position (0.5,0.5,10.5)
zone history name="zdisp" displacementz position (0,0,10)
model history name="time" fluid timetotal
;  drained response 
zone face apply porepressure 0 range group 'Top'
zone fluid property mobilitycoefficient 1e10
; Loosely coupled till time 5000
zone fluid implicit servo on
model solvefluidcoupled time 5000 convergence 1
model save '1dconsolidation_coupled_10M'
program return
1DConsolidationUncoupled.dat
model new
model title "Onedimensional consolidation (uncoupled)"
zone create brick size 1 1 20
zone face skin
;  mechanical model 
model largestrain off
zone cmodel assign elastic
zone property bulk [5e8+4e10] shear 2e8 ; Add fluid modulus
zone face apply velocitynormal 0.0 range group 'Top' not
zone face apply stressnormal 1e5 range group 'Top'
model solve convergence 1
zone property bulk 5e8
;  fluid flow model 
model configure fluidflow
zone fluid property mobilitycoefficient 1e10 biotmodulus 4e10
zone fluid modulusscale
; This is 1e5/(confinedModulus*storage)
zone gridpoint initialize porepressure 98119.4
zone face apply porepressure 0 range group 'Top'
zone results porepressure modelmechanical modelfluid
;  histories 
zone fluid implicit servo on
model solvefluiddecoupled time 5000 convergence 2 stages 10 ...
stageresult '1dconsolidation_uncoupled_'
model save '1dconsolidation_uncoupled'
program return
1DConsolidationFastFlow.dat
model new
model title "Onedimensional consolidation (coupled  fast flow)"
;  model geometry 
zone create brick size 1 1 20
zone face skin
;  mechanical model 
model largestrain off
zone cmodel assign elastic
zone property bulk 5e8 shear 2e8
zone face apply velocitynormal 0.0 range group 'Top' not
zone face apply stressnormal 1e5 range group 'Top'
;  fluid flow model 
model configure fluidflow
zone fluid property mobilitycoefficient 1e10
;  fast flow setting 
zone fluid fastflow on ; implies explicitsaturation off and implicit off, ppgen on
;  first establish undrained response 
model solvestatic convergence 0.01
;  histories 
zone history name="pp" porepressure position (0,0,10)
zone history name="szz" stresszz position (0.5,0.5,10.5)
zone history name="zdisp" displacementz position (0,0,10)
model history name="time" fluid timetotal
;  drained response 
zone face apply porepressure 0 range group 'Top'
model solvefluidcoupled time 5000 convergence 0.01
model save '1dconsolidation_fastflow'
program return
Was this helpful? ...  Itasca Software © 2024, Itasca  Updated: Sep 26, 2024 