FLAC3D Theory and Background • FluidMechanical Interaction
OneDimensional Consolidation
Note
The project file for this example is available to be viewed/run in FLAC3D.[1] 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 fluid
. 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). 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:
fish define my_solve
loop local ii (1,10)
global c_age = 500.0*ii
command
zone property bulk [c_bulk]
zone gridpoint initialize biot [c_biotma]
model fluid active on mech off
model solve fluid timetotal [c_age]
model fluid active off mech on
zone gridpoint initialize biot 0.0
model solve ratio 2e5
end_command
end_loop
end
A simplified method is to use the zone consolidation
command:
zone property bulk [c_bulk]
zone gridpoint initialize biot [c_biotm]
fish define my_solve
loop local ii (1,10)
global c_age = 500.0*ii
command
zone consolidation [c_age] ratio 2e5
end_command
end_loop
end
The results using the above simplified uncoupled method are presented in Figure 13 to Figure 15.
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 16 and Figure 17. 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 largestrain off
model title "Onedimensional consolidation (coupled)"
fish automaticcreate off
model configure fluid
program call 'fishFunctions'
[setup(4.0e9)]
;  model geometry 
zone create brick size 1 1 [int(hh)]
zone face skin
;  mechanical model 
zone cmodel assign elastic
zone fluid biot on
zone property bulk [c_bulk] sh [c_shear]
zone face apply velocityx 0.0 range group 'West' or 'East'
zone face apply velocityy 0.0 range group 'North' or 'South'
zone face apply velocityz 0.0 range group 'Bottom'
zone initialize stresszz 0.
zone face apply stresszz [sig0] range group 'Top'
;  fluid flow model 
zone fluid cmodel assign isotropic
zone fluid property permeability [c_perm] biot [c_biotc]
zone gridpoint initialize biot [c_biotm]
zone gridpoint initialize porepressure 0
;  first establish undrained response 
model fluid active off
model solve ratio 1e4
;  histories 
history interval 200
fish history pp10
fish history c_szz
fish history c_eszz
fish history c_uz
fish history ft
;  drained response 
zone face apply porepressure 0 range group 'Top'
model fluid active on
model fluid substep 1
model mechanical substep 1
model mechanical follower on
model solve fluid timetotal 500 or mechanical ratio 1e4
model fluid substep 100
model mechanical substep 10
model mechanical follower on
model solve fluid timetotal 5000 or mechanical ratio 1e4
history export '1' vs '5' table '1'
history export '4' vs '5' table '2'
history export '2' vs '5' table '3'
history export '3' vs '5' table '4'
table '1' label 'Pore Pressure'
table '2' label 'Vertical Displacement'
table '3' label 'Total Stress'
table '4' label 'Effective Stress'
model save '1dcons_coup'
program return
1DConsolidationCoupled10M.dat
model new
model largestrain off
model title "Onedimensional consolidation (coupled)"
fish automaticcreate off
model configure fluid
program call 'fishFunctions'
[setup(4.0e10)]
;  model geometry 
zone create brick size 1 1 [int(hh)]
zone face skin
;  mechanical model 
zone cmodel assign elastic
zone fluid biot on
zone property bulk [c_bulk] sh [c_shear]
zone face apply velocityx 0.0 range group 'West' or 'East'
zone face apply velocityy 0.0 range group 'North' or 'South'
zone face apply velocityz 0.0 range group 'Bottom'
zone initialize stresszz 0.
zone face apply stresszz [sig0] range group 'Top'
;  fluid flow model 
zone fluid cmodel assign isotropic
zone fluid property permeability [c_perm] biot [c_biotc]
zone gridpoint initialize biot [c_biotm]
zone gridpoint initialize porepressure 0
;  first establish undrained response 
model fluid active off
model solve ratio 2e5
;  histories 
history interval 1000
fish history pp10
fish history c_szz
fish history c_eszz
fish history c_uz
fish history ft
;  drained response 
zone face apply porepressure 0 range group 'Top'
model fluid active on
model fluid substep 1
model mechanical substep 1
model mechanical follower on
model solve fluid timetotal 500 or mechanical ratio 2e5
model fluid substep 100
model mechanical substep 10
model mechanical follower on
model solve fluid timetotal 5000 or mechanical ratio 2e5
history export '1' vs '5' table '1'
history export '4' vs '5' table '2'
history export '2' vs '5' table '3'
history export '3' vs '5' table '4'
table '1' label 'Pore Pressure'
table '2' label 'Vertical Displacement'
table '3' label 'Total Stress'
table '4' label 'Effective Stress'
model save '1dcons_coup10M'
program return
1DConsolidationUncoupled.dat
model new
model largestrain off
model title "Onedimensional consolidation (uncoupled)"
fish automaticcreate off
model configure fluid
program call 'fishFunctions'
[setup(4.0e10)]
;  model geometry 
zone create brick size 1 1 [int(hh)]
zone face skin
;  mechanical model 
zone cmodel assign elastic
zone property bulk [c_bulk] shear [c_shear]
zone face apply velocityx 0.0 range group 'West' or 'East'
zone face apply velocityy 0.0 range group 'North' or 'South'
zone face apply velocityz 0.0 range group 'Bottom'
zone initialize stresszz 0.
zone face apply stresszz [sig0] range group 'Top'
;  fluid flow model 
zone fluid cmodel assign isotropic
zone fluid property permeability [c_perm] biot [c_biotc]
zone fluid biot on
zone gridpoint initialize biot [c_biotm]
zone gridpoint initialize porepressure 0
;  first establish undrained response 
model fluid active off mech on
zone property bulk [c_bulku]
zone gridpoint initialize biot 0.0
model solve
zone gridpoint initialize porepressure [p0]
;  drained response 
zone face apply porepressure 0 range group 'Top'
fish define my_solve
loop local ii (1,10)
global c_age = 500.0*ii
command
zone property bulk [c_bulk]
zone gridpoint initialize biot [c_biotma]
model fluid active on mech off
model solve fluid timetotal [c_age]
model fluid active off mech on
zone gridpoint initialize biot 0.0
model solve ratio 2e5
end_command
table.y('1',ii) = pp10
table.x('1',ii) = ft
table.y('2',ii) = c_uz
table.x('2',ii) = ft
table.y('3',ii) = c_szz
table.x('3',ii) = ft
table.y('4',ii) = c_eszz
table.x('4',ii) = ft
end_loop
table.label('1') = 'Pore Pressure'
table.label('2') = 'Vertical Displacement'
table.label('3') = 'Total Stress'
table.label('4') = 'Effective Stress'
end
[my_solve]
model save '1dcons_uncoup'
program return
1DConsolidationSimplified.dat
model new
model largestrain off
model title "Onedimensional consolidation simplified"
fish automaticcreate off
model configure fluid
program call 'fishFunctions'
[setup(4.0e10)]
;  model geometry 
zone create brick size 1 1 [int(hh)]
zone face skin
;  mechanical model 
zone cmodel assign elastic
zone property bulk [c_bulk] shear [c_shear]
zone face apply velocityx 0.0 range group 'West' or 'East'
zone face apply velocityy 0.0 range group 'North' or 'South'
zone face apply velocityz 0.0 range group 'Bottom'
zone initialize stresszz 0.
zone face apply stresszz [sig0] range group 'Top'
;  fluid flow model 
zone fluid cmodel assign isotropic
zone fluid property permeability [c_perm] biot [c_biotc]
zone fluid biot on
zone gridpoint initialize biot [c_biotm]
zone gridpoint initialize porepressure 0
;  first establish undrained response 
model fluid active off mech on
zone property bulk [c_bulku]
zone gridpoint initialize biot 0.0
model solve
zone gridpoint initialize porepressure [p0]
;  drained response 
zone face apply porepressure 0 range group 'Top'
;
zone property bulk [c_bulk] ; drained bulk modulus here!
zone gridpoint initialize biot [c_biotm] ; either biotmodulus
;zone gridpoint initialize fluidmodulus 2.0e10 ; or fluidmodulus
;
fish define my_solve
loop local ii (1,10)
global c_age = 500.0*ii
command
zone consolidation [c_age] ratio 2e5
end_command
table.y('1',ii) = pp10
table.x('1',ii) = ft
table.y('2',ii) = c_uz
table.x('2',ii) = ft
table.y('3',ii) = c_szz
table.x('3',ii) = ft
table.y('4',ii) = c_eszz
table.x('4',ii) = ft
end_loop
table.label('1') = 'Pore Pressure'
table.label('2') = 'Vertical Displacement'
table.label('3') = 'Total Stress'
table.label('4') = 'Effective Stress'
end
[my_solve]
model save '1dcons_simplified'
program return
1DConsolidationFastFlow.dat
model new
model largestrain off
model title "Onedimensional consolidation (coupled  fast flow)"
fish automaticcreate off
model configure fluid
program call 'fishFunctions'
[setup(4.0e10)]
;  model geometry 
zone create brick size 1 1 [int(hh)]
zone face skin
;  mechanical model 
zone cmodel assign elastic
; set fluid biot on
zone property bulk [c_bulk] sh [c_shear]
zone face apply velocityx 0.0 range group 'West' or 'East'
zone face apply velocityy 0.0 range group 'North' or 'South'
zone face apply velocityz 0.0 range group 'Bottom'
zone initialize stresszz 0.
zone face apply stresszz [sig0] range group 'Top'
;  fluid flow model 
zone fluid cmodel assign isotropic
zone fluid property permeability [c_perm]
zone gridpoint initialize porepressure 0
;  fast flow setting 
zone fluid fastflow on
;  first establish undrained response 
model fluid active off
model solve ratio 1e5
;  histories 
fish history pp10
fish history c_szz
fish history c_eszz
fish history c_uz
fish history ft
;  drained response 
zone gridpoint fix porepressure 0 range positionz 20
model fluid active on
model mechanical substep 100
model mechanical follower on
model fluid substep 10
model step 0
model fluid timestep fix 1
model solve fluid timetotal 200 or mechanical ratio 1e7
model fluid timestep auto
model solve fluid timetotal 5000 or mechanical ratio 1e7
history export '1' vs '5' table '1'
history export '4' vs '5' table '2'
history export '2' vs '5' table '3'
history export '3' vs '5' table '4'
table '1' label 'Pore Pressure'
table '2' label 'Vertical Displacement'
table '3' label 'Total Stress'
table '4' label 'Effective Stress'
model save '1dcons_ff'
program return
Endnotes
Was this helpful? ...  Itasca Software © 2022, Itasca  Updated: Mar 09, 2023 