FLAC3D Theory and Background • Fluid-Mechanical Interaction

One-Dimensional 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⁵ 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 one-dimensional consolidation problem may be expressed in the framework of Biot theory (see Detournay and Cheng 1993).

Figure 1: One-dimensional consolidation.

The diffusion equation for the pore pressure, \(p\), has the form

(1)\[{{\partial p}\over{\partial t}} - c {{{\partial}^2p}\over{\partial z^2}} = - {{\alpha}\over{\alpha_1 S}} {{d {\sigma}_{zz}}\over{dt}}\]

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

(2)\[{{\partial p}\over{\partial t}} - c {{{\partial}^2p}\over{\partial z^2}} = 0\]

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 one-dimensional mechanical constitutive law (\({\sigma}_{zz} = \alpha_1 {\epsilon}_{zz} - \alpha p\)) to express strain in terms of stresses. It is given by

(3)\[p_0 = {{\alpha}\over{\alpha_1}} {{1}\over{S}} {p_z}\]

The solution for the pore pressure is

(4)\[\hat{p} = 2{{p_0}\over{p_z}} \sum_{m=0}^{\infty} {{\sin (a_m \hat{z})}\over{a_m}} e^{-a_m^2 \hat{t}}\]

where:	\(\hat{p}\)	=	\(p/p_z\);
	\(a_m\)	=	\((2 m +1)\pi/2\);
	\(\hat{z}\)	=	\((H-z)/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,

(5)\[{\hat{u}}_z = 2 \alpha {{p_0}\over{p_z}} \left [ \sum_{m = 0}^{\infty} {{\cos (a_m \hat{z})}\over{a_m^2}} e^{-a_m^2 \hat{t}} \right ] + {\hat{y} - 1}\]

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 × 108 Pa
dry shear modulus, \(G\)	2 × 108 Pa
Biot modulus, \(M\)	4 × 109 Pa
Biot coefficient, \(\alpha\)	1.0
permeability, \(k\)	10-10 m2/Pa-sec

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 mid-height are evaluated using FISH functions, and compared to the numerical solution as time proceeds. The results are plotted versus fluid-flow 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 fluid-flow time are presented. The FLAC3D data file is listed in “1DConsolidation-Coupled.dat”.

break

Figure 2: FLAC3D grid for one-dimensional consolidation.

Figure 3: Comparison between analytical and numerical values of pore pressure in a 1D consolidation test — \(M\) = 4 × 10⁹ Pa.

Figure 4: Comparison between analytical and numerical values for vertical displacement in a 1D consolidation test — \(M\) = 4 × 10⁹ Pa.

Figure 5: Evolution of pore pressure, total and effective stresses in a 1D consolidation test — \(M\) = 4 × 10⁹ Pa.

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 pore-pressure 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 Flow-Only and Coupled-Flow Problems). To demonstrate the validity of this numerical approach, the analytical solutions for pore pressure and displacement corresponding to \(M\) = 1.5 × 10¹⁰ (\(R_k \simeq\) 20) and \(M\) = 1.5 × 10¹⁵ (\(R_k \simeq\) 20 × 10⁵) 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

Figure 6: Comparison between analytical pore-pressure solutions for two large values of M in a 1D consolidation test.

Figure 7: Comparison between analytical vertical displacement solutions for two large values of M in a 1D consolidation test.

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¹⁰ Pa (10 times larger than in the former numerical calculation, see 1DConsolidation-Coupled10M.dat). In the former calculation with \(M\) = 4 × 10⁹ 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¹⁰ 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¹⁰ Pa are shown in Figure 8 and Figure 9, respectively. The maximum relative error is less than 4%.

break

Figure 8: Comparison between analytical and numerical values of pore pressure in a 1D consolidation test - \(M\) = 4 × 10¹⁰ Pa.

Figure 9: Comparison between analytical and numerical values for vertical displacement in a 1D consolidation test - \(M\) = 4 × 10¹⁰ Pa.

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 (1DConsolidation-Uncoupled.dat), the numerical simulation is repeated for a Biot modulus of 4 × 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 pore-pressure 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

Figure 10: Comparison between analytical and numerical values of pore pressure in a 1D consolidation test - \(M\) = 4 × 10¹⁰ Pa.

Figure 11: Comparison between analytical and numerical values for vertical displacement in a 1D consolidation test - \(M\) = 4 × 10¹⁰ Pa.

Figure 12: Evolution of pore pressure, total and effective stresses in a 1D consolidation test - \(M\) = 4 × 10¹⁰ Pa.

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 time-total [c_age]
            model fluid active off mech on
            zone gridpoint initialize biot 0.0
            model solve ratio 2e-5
        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 2e-5
        end_command
    end_loop
end

The results using the above simplified uncoupled method are presented in Figure 13 to Figure 15.

break

Figure 13: Comparison between analytical and numerical values of pore pressure in a 1D consolidation test, simplified method - \(M\) = 4 × 10¹⁰ Pa.

Figure 14: Comparison between analytical and numerical values for vertical displacement in a 1D consolidation test, simplified method - \(M\) = 4 × 10¹⁰ Pa.

Figure 15: Evolution of pore pressure, total and effective stresses in a 1D consolidation test, simplified method - \(M\) = 4 × 10¹⁰ Pa.

This particular example involves fully saturated coupled flow. In such cases, solution times can be reduced by using the saturated fast-flow algorithm described in Fully Saturated Fast Flow. The algorithm is selected by using the zone fluid fastflow on command. If 1DConsolidation-Coupled.dat with \(M\) = 4 × 10¹⁰ Pa (\(R_k\) = 52) is repeated with this command, the runtime is now approximately 12 seconds (see 1DConsolidation-FastFlow.dat). As \(R_k\) increases, the use of the fast-flow algorithm becomes more beneficial.

The results of the fast-flow calculation are shown in Figure 16 and Figure 17. The maximum relative error is less than 2%.

break

Figure 16: Comparison between analytical and numerical values of pore pressure in a 1D consolidation test.

Figure 17: Comparison between analytical and numerical values for vertical displacement in a 1D consolidation test.

Reference

Detournay, E., and A. H. D. Cheng. “Fundamentals of Poroelasticity,” in Comprehensive Rock Engineering, Vol. 2, pp. 113-171. J. Hudson et al., eds. London: Pergamon Press (1993).

Data File

1DConsolidation-Coupled.dat

model new
model large-strain off
model title "One-dimensional consolidation (coupled)"
fish automatic-create 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] shear [c_shear]
zone face apply velocity-x 0.0 range group 'West' or 'East'
zone face apply velocity-y 0.0 range group 'North' or 'South'
zone face apply velocity-z 0.0 range group 'Bottom'
zone initialize stress-zz 0.
zone face apply stress-zz [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 pore-pressure 0
; --- first establish undrained response ---
model fluid active off
model solve ratio 1e-4
; --- 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 pore-pressure 0 range group 'Top'
model fluid active on
model fluid substep 1
model mechanical substep 1
model mechanical follower on
model solve fluid time-total 500 or mechanical ratio 1e-4
model fluid substep 100
model mechanical substep 10 
model mechanical follower on
model solve fluid time-total 5000 or mechanical ratio 1e-4
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

1DConsolidation-Coupled10M.dat

model new
model large-strain off
model title "One-dimensional consolidation (coupled)"
fish automatic-create 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] shear [c_shear]
zone face apply velocity-x 0.0 range group 'West' or 'East'
zone face apply velocity-y 0.0 range group 'North' or 'South'
zone face apply velocity-z 0.0 range group 'Bottom'
zone initialize stress-zz 0.
zone face apply stress-zz [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 pore-pressure 0
; --- first establish undrained response ---
model fluid active off
model solve ratio 2e-5
; --- 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 pore-pressure 0 range group 'Top'
model fluid active on
model fluid substep 1
model mechanical substep 1
model mechanical follower on
model solve fluid time-total 500 or mechanical ratio 2e-5
model fluid substep 100
model mechanical substep 10 
model mechanical follower on
model solve fluid time-total 5000 or mechanical ratio 2e-5
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

1DConsolidation-Uncoupled.dat

model new
model large-strain off
model title "One-dimensional consolidation (uncoupled)"
fish automatic-create 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 velocity-x 0.0 range group 'West' or 'East'
zone face apply velocity-y 0.0 range group 'North' or 'South'
zone face apply velocity-z 0.0 range group 'Bottom'
zone initialize stress-zz 0.
zone face apply stress-zz [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 pore-pressure 0
; --- first establish undrained response ---
model fluid active off mechanical on
zone property bulk [c_bulku]
zone gridpoint initialize biot 0.0
model solve
zone gridpoint initialize pore-pressure [p0]
; --- drained response ---
zone face apply pore-pressure 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 mechanical off
            model solve fluid time-total [c_age]
            model fluid active off mechanical on
            zone gridpoint initialize biot 0.0
            model solve ratio 2e-5
        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

1DConsolidation-Simplified.dat

model new
model large-strain off
model title "One-dimensional consolidation simplified"
fish automatic-create 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 velocity-x 0.0 range group 'West' or 'East'
zone face apply velocity-y 0.0 range group 'North' or 'South'
zone face apply velocity-z 0.0 range group 'Bottom'
zone initialize stress-zz 0.
zone face apply stress-zz [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 pore-pressure 0
; --- first establish undrained response ---
model fluid active off mechanical on
zone property bulk [c_bulku]
zone gridpoint initialize biot 0.0
model solve
zone gridpoint initialize pore-pressure [p0]
; --- drained response ---
zone face apply pore-pressure 0 range group 'Top'
;
zone property bulk [c_bulk] ; drained bulk modulus here!
zone gridpoint initialize biot [c_biotm] ; either biot-modulus 
;zone gridpoint initialize fluid-modulus 2.0e10 ; or fluid-modulus 
;
fish define my_solve
    loop local ii (1,10)
        global c_age = 500.0*ii
        command
            zone consolidation [c_age] ratio 2e-5 
        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

1DConsolidation-FastFlow.dat

model new
model large-strain off
model title "One-dimensional consolidation (coupled - fast flow)"
fish automatic-create 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] shear [c_shear]
zone face apply velocity-x 0.0 range group 'West' or 'East'
zone face apply velocity-y 0.0 range group 'North' or 'South'
zone face apply velocity-z 0.0 range group 'Bottom'
zone initialize stress-zz 0.
zone face apply stress-zz [sig0] range group 'Top'
; --- fluid flow model ---
zone fluid cmodel assign isotropic
zone fluid property permeability [c_perm]
zone gridpoint initialize pore-pressure 0
; --- fast flow setting ---
zone fluid fastflow on
; --- first establish undrained response ---
model fluid active off
model solve ratio 1e-5
; --- histories ---
fish history pp10
fish history c_szz
fish history c_eszz
fish history c_uz
fish history ft
; --- drained response ---
zone gridpoint fix pore-pressure 0 range position-z 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 time-total 200 or mechanical ratio 1e-7
model fluid timestep automatic
model solve fluid time-total 5000 or mechanical ratio 1e-7
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