# Compression of a Poroelastic Sample – Mandel’s Problem (FLAC2D)

## Problem Statement

Note

The project file for this example is available to be viewed/run in FLAC2D.[1] The main data files used are shown at the end of this example. The remaining data files can be found in the project.

During compression of a poroelastic specimen under constant boundary conditions, the pore pressure will display a non-monotonic variation with consolidation time. At initial consolidation times, an increase in the pore pressure will be induced near the center of the sample when it is subjected to a constant vertical load and drained laterally. Subsequently, the pore pressure falls. This effect pointed out by Mandel (1953). It was also predicted by Cryer (1963), and thus is also known as the Mandel-Cryer effect, and was demonstrated experimentally by Verruijt (1965).

In Mandel’s problem, a sample of saturated poroelastic material is loaded under plane-strain conditions by a constant compressive force applied on rigid impervious platens (see Figure 1). The width of the sample is 2a its height is 2b, and the force intensity is 2F. The application of the load is instantaneous, the platens are impervious and the sample is free to drain laterally.

The short-term response of the material corresponds to a uniform vertical stress across the sample. As lateral drainage takes place, the nonuniform dissipation of induced pore pressure causes an apparent softening of the material near the edges of the sample. The resulting stress concentration in the stiffer (still undrained) core is then responsible for an additional increase in pore pressure in the middle of the sample. As drainage proceeds and the pore pressure gradient decreases, the vertical load is again transmitted in a uniform manner. The non-monotonic variation of pore pressure with time observed in Mandel’s problem serves to illustrate the main difference in prediction between the Biot and Terzaghi theories.

## Analytical Solution

The solution to Mandel’s problem, generalized for the case of compressible constituents, is given by Cheng and Detournay (1988). The expression for the pore pressure is

$p = {{2FB(1+\nu_u)}\over{3a}} \sum_{i=1}^{\infty} {{\sin \alpha_i}\over{\alpha_i - \sin \alpha_i \cos \alpha_i}} \left[ \cos {{\alpha_i x}\over a} - \cos \alpha_i \right] \exp (- \alpha_i^2 ct/a^2)$

where $$B$$ is Skempton’s pore pressure coefficient (the ratio of induced pore pressure to variation of confining pressure under undrained conditions), $$ν$$ and $$ν_u$$ are drained and undrained Poisson’s ratio, $$c$$ is the true diffusivity (or generalized consolidation coefficient), $$t$$ is time and $$α_i, i = 1, ∞$$, are the roots of the equation

$\tan \alpha_i = {{1 - \nu}\over{\nu_u - \nu}} \alpha_i$

If the solid grains that form the material are assumed to be incompressible, as is the case in FLAC2D, the following relations between material constants apply.

\eqalignno{ B &= {{K_w} \over {K_w + n K}} &\cr K_u &= K + {K_w \over n} \cr c &= k \ / \left({n \over K_w} + {1 \over {K + 4 G /3}}\right) & \cr }

where $$K_w$$ is fluid bulk modulus, $$n$$ is porosity, $$K$$ and $$K_u$$ are drained and undrained bulk moduli, $$G$$ is shear modulus and $$k$$ is mobility coefficient. The drained and undrained Poisson’s ratios are related to the moduli $$G$$, $$K$$ and $$K_u$$ as

\eqalignno{ \nu &= {{3 K - 2 G} \over {6 K + 2 G}} &\cr \ \ \ & &\ \ \cr \nu_u &= {{3 K_u - 2 G} \over {6 K_u + 2 G}} & \cr }

The formulae for horizontal displacement at $$x = a$$ and vertical displacement at $$y = b$$ are

\eqalignno{ u_x(a,t) &= {{F \nu} \over {2G}} + {{F (1 - \nu_u)} \over {G}} \sum_{i=1}^{\infty} {{\sin \alpha_i \cos \alpha_i}\over{\alpha_i - \sin \alpha_i \cos \alpha_i}} \exp (- \alpha_i^2 ct/a^2) &\cr \ \ \ & &\ \ \cr u_y(b,t) &= - {{F(1- \nu)b} \over {2Ga}} + {{F (1 - \nu_u)b} \over {G a}} \sum_{i=1}^{\infty} {{\sin \alpha_i \cos \alpha_i}\over{\alpha_i - \sin \alpha_i \cos \alpha_i}} \exp (- \alpha_i^2 ct/a^2) & \cr }

According to the exact solution, the initial (instantaneous) and final vertical displacements of the upper platen are

\eqalignno{ u_y(b,0) &= -Fb {{1-\nu_u} \over{2Ga}} &\cr \ \ \ & &\ \ \cr u_y(b,\infty) &= -Fb {{1-\nu} \over{2Ga}} & \cr }

The lateral boundaries of the sample stay plane during drainage; they first move outward and then inward. The initial (instantaneous) and final lateral displacements of the right boundary are

\eqalignno{ u_x(a,0) &= F {{\nu_u} \over{2G}} &\cr \ \ \ & &\ \ \cr u_x(a,\infty) &= F {{\nu} \over{2G}} & \cr }

The degrees of consolidation in the horizontal and vertical directions, $$D_h$$ and $$D_v$$, are defined as

\eqalignno{ D_h &= {{u_x(a,t) - u_x(a,0)} \over{u_x(a,\infty) - u_x(a,0)}} &\cr \ \ \ & &\ \ \cr D_v &= {{u_y(b,t) - u_y(b,0)} \over{u_y(b,\infty) - u_y(b,0)}} & \cr }

These degrees are found to be identical, and the analytic expression for $$D = D_h = D_v$$ is

$D = 1 - {{4(1-\nu_u)} \over {1-2\nu}} \sum_{i=1}^{\infty} {{\sin \alpha_i \cos \alpha_i}\over{\alpha_i - \sin \alpha_i \cos \alpha_i}} \exp (- \alpha_i^2 ct/a^2)$

Because the discharge has only a horizontal component in Mandel’s problem, the pore pressure, stress and strain solutions are independent of the $$y$$-coordinate.

## FLAC Model

By symmetry, only a quarter of the sample is considered in the numerical model. The grid has 20 zones in the x-direction and 2 in the y-direction. Mandel’s problem is solved for the case

 drained Poisson’s ratio ($$ν$$) 0.2 Skempton coefficient ($$B$$) 0.9

The FLAC2D simulation is carried out to produce results in terms of normalized pore pressure, $$\hat{p}$$, distance, $$\hat{x}$$, and time, $$\hat{t}$$, defined as

\eqalignno{ \hat p &= {{ap} \over{F}} &\cr \hat x &= {{x} \over{a}} &\cr \hat t &= {{ct} \over{a^2}} & \cr }

The normalized results are not affected by the absolute magnitude of material properties used, as long as their combination yields the values specified above for $$ν_u$$ and $$B$$. The FLAC2D grid dimensions, applied force and property values used in the simulation may be viewed as scaled for the purpose of producing the normalized results. They are selected as follows.

 model width ($$a$$) 1 model height ($$b$$) 0.1 applied force ($$F$$) 1 drained bulk modulus ($$K$$) 1 shear modulus ($$G$$) 0.75 fluid bulk modulus ($$K_w$$) 9 porosity ($$n$$) 0.5

With the above values, the true diffusivity in the FLAC2D model is unity, and the model time is t ̂. The model mechanical boundary conditions correspond to roller boundaries along the x- and y-axes of symmetry. The sample is initially stress-free, and the pore pressure is equal to zero. Undrained conditions are established first by applying a constant unit mechanical pressure at the top boundary of the model. In the second part of the simulation, the pore pressure is fixed at zero on the right side of the model to allow drainage to occur. The rigid plate condition is enforced by applying a vertical velocity at the top of the model. The velocity magnitude is derived from the exact displacement solution . As a verification to the numerical solution, the reaction force on the top platen is monitored to check whether it remains constant and equal to one.

## Results

The numerical simulation is carried out for a total value of normalized time equal to 4, with intermediate results at $$\hat{t} = 0.01, 0.1, 0.5, 1$$ and $$2$$. (Results at $$\hat{t} = 0$$ correspond to the undrained response.) The pore pressure profiles at those times are checked against exact solutions in Figure 2. At early times, the pore pressure at $$x = 0$$ (center of sample) is seen to rise above the undrained value, before decreasing as drainage evolves. Figure 2 This is also shown in Figure 3, which compares FLAC2D to the analytical solution for the actual pore pressure versus consolidation time at the center of the sample. As may be seen in Figure 4, the reaction force stays equal to unity throughout the simulation. (The approach taken to apply the force boundary condition is thus equivalent to a servo-controlled velocity.) Numerical values for the horizontal and vertical degrees of consolidation are plotted versus the analytical solution values in Figure 5.

## References

Cheng, A. H. D., and E. Detournay. “A Direct Boundary Element Method for Plane Strain Poroelasticity,” Int. J. Num. Methods and Anal. in Geomechanics, 12, 551-572 (1988).

Cryer, C.W. “A Comparison of the Three-Dimensional Consolidation Theories of Biot and Terzaghi,” Quart. J. Mech. and Appl. Math., XVI, 4, 401-412 (1963).

Mandel, J. “Consolidation des sols (étude mathématique),” Géotechnique, 3, 287-299 (1953).

Verruijt, A. “Discussion,” Proc. 6th Int. Conf. Soil Mechanics and Foundation Engineering, 3, Montréal, 401-402 (1965).

Data Files

mandel.dat

;-----------------------------------------------------------
; Compression of a Poroelastic Sample – Mandel’s Problem
;-----------------------------------------------------------
model new
model large-strain off
fish automatic-create off
model title "Compression of a Poroelastic Sample– Mandel’s Problem   "
model config fluid
program call 'fishfunctions.dat'
[setup]
; --- model geometry ---
zone create quad size 20 2 p 1=([c_a],0.0) p 2=(0.0,[c_b]) p 3=([c_a],[c_b])
zone face skin
zone cmodel assign elastic
zone property bulk [c_bu] shear [c_sh] density 1.0
;fluid properties
zone fluid cmodel assign isotropic
zone fluid property permeability [c_wk] porosity [c_po]
zone gridpoint initialize fluid-modulus [c_wb]
; --- boundary conditions ---
zone gridpoint fix velocity-x range group 'West'
zone gridpoint fix velocity-y range group 'Bottom
zone face apply stress-yy -1.0 range group 'Top'
; --- undrained response ---
model fluid active off
model solve ratio 1e-5

[check_ppu(10)]
model save 'man0'

; --- drained conditions ---
zone face apply-remove range group 'Top'
zone face apply pore-pressure 0  range group 'East'
zone gridpoint initial velocity-y 0 range group 'Top'
zone gridpoint fix velocity-y range group 'Top'
model fluid active on

fish history ana_pp
fish history num_pp
model history fluid time-total
; --- drained test ---

model fluid timestep fix 8.0e-5
[big_step(125)]
[check_pp(11)]
model save 'man1'

[big_step(1125)]
[check_pp(12)]
model save 'man2'
[big_step(5000)]
[check_pp(13)]
model save 'man3'

[big_step(6250)]
[check_pp(14)]
model save 'man4'

[big_step(12500)]
[check_pp(15)]
model save 'man5'

[big_step(25000)]
[check_pp(16)]
model save 'man6'

[deg_cons(7)]
model save 'man7'


fishfunctions.dat

;fish functions
fish def setup
global c_F  = 1.        ; applied force
global c_a  = 1.        ; sample half width
global c_b  = 0.1       ; sample half height
global c_wb = 4.5       ; water bulk modulus
global c_wk = 11./18.   ; mobility coefficient
global c_po = 0.5       ; porosity
global c_bu = 1.        ; drained bulk modulus
global c_sh = 0.75      ; shear modulus
global c_ub = c_bu + c_wb/c_po              ; undrained bulk modulus
global c_sk = c_wb/(c_wb + c_po*c_bu)       ; Skempton coefficient
global val  = c_bu/c_sh
global c_nu = (3.*val-2.)/(6.*val+2.)       ; drained Poisson ratio
global val  = c_ub/c_sh
global c_un = (3.*val-2.)/(6.*val+2.)       ; undrained Poisson ratio
global stor = c_po/c_wb+1./(c_bu+4.*c_sh/3.); storativity
global diff = c_wk/stor                     ; diffusivity
global coe  = diff/(c_a*c_a)
global v_co = -(1.-c_un)*c_F*c_b/(c_sh*c_a) ; coeff. for applied vel
global v_co = v_co*coe
global c_nm = 100                    ; number of mech steps per flow step
end
[setup]
;[pnt = list(gp.list)(gp.isgroup(::gp.list,"Top"))]
fish def check_ppu(taba)  ; undrained pore pressure profile
local pcoe = c_F*c_sk*(1.+c_un)/(3.*c_a)
global tabn = taba + 10
local ip = 1
loop foreach local pnt list(gp.list)(gp.isgroup(::gp.list,"Top"))
local xval  = gp.pos.x(pnt)
local pp    = gp.pp(pnt)
table.x(taba,ip) = xval
table.y(taba,ip) = pcoe
table.x(tabn,ip) = xval
table.y(tabn,ip) = pp
ip=ip+1
end_loop
end

;
fish def ini_root
global tabroot = 100
global nroot   = 50
global tab1  = tabroot + 1
global tab2  = tabroot + 2
global tab3  = tabroot + 3
global tab4  = tabroot + 4
global tab5  = tabroot + 5
global c_coe = (1. - c_nu)/(c_un - c_nu)
global pcoe  = 2.*c_F*c_sk*(1.+c_un)/(3.*c_a)
global tbi   = 0
end
[ini_root]

; calculate and store roots
fish def func(c_x0)
global func = math.tan(c_x0)/(c_coe*c_x0) - 1.
end
;Name: froot
;Purpose: Finds the root of a function FUNC=f(x)
;Diagram:
;Input: c_x1/float/0.99/interval lower bound
;Input: c_x2/float/1.01/interval upper bound
;Input: tol/float/1e-4/tolerance
;Code: FLAC
;Diagram:
;Note: Using Brent's method, find the root of a function FUNC=f(x)
; bracketed in the interval [c_x1,c_x2]
; (i.e. such that f(c_x1)*f(c_x2) < 0)
; input:      c_x1 c_x2      interval bounds
;             func           fish function with input c_x0
;             tol
; output:     froot          root value
;
; Numerical Recipes: function zbrent
;
fish def froot(c_x1,c_x2,tol)
local itmax = 100
local c_eps = 3.e-8
local c_x0 = c_x1
local fa = func(c_x0)
c_x0 = c_x2
local fb = func(c_x0)
local c_a = c_x1
local c_b = c_x2
if fa*fb > 0. then
local toto=io.out(' root must be bracketed for froot')
exit
end_if
local fc = fb
loop local iter (1,itmax)
if fb*fc > 0. then
local c_c = c_a
fc = fa
local c_d = c_b - c_a
local c_e = c_d
end_if
if math.abs(fc) < math.abs(fb) then
c_a = c_b
c_b = c_c
c_c = c_a
fa = fb
fb = fc
fc = fa
end_if
local tol1 = 2.*c_eps*math.abs(c_b)+0.5*tol
local xm = 0.5*(c_c-c_b)
if math.abs(xm) <= tol1 then
froot = c_b
exit
end_if
if fb = 0. then
froot = c_b
exit
end_if
if math.abs(c_e) >= tol1 then
if math.abs(fa) > math.abs(fb) then
local c_s = fb/fa
if c_a = c_c then
local c_p = 2.*xm*c_s
local c_q = 1. - c_s
else
c_q = fa/fc
local c_r = fb/fc
c_p = c_s*(2.*xm*c_q*(c_q-c_r)-(c_b-c_a)*(c_r-1.))
c_q = (c_q-1.)*(c_r-1.)*(c_s-1.)
end_if
if c_p > 0. then
c_q = -c_q
end_if
c_p=math.abs(c_p)
if 2.*c_p<math.min(3.*xm*c_q-math.abs(tol1*c_q),math.abs(c_e*c_q)) then
c_e = c_d
c_d = c_p/c_q
else
c_d = xm
c_e = c_d
end_if
else
c_d = xm
c_e = c_d
end_if
else
c_d = xm
c_e = c_d
end_if
c_a = c_b
fa = fb
if math.abs(c_d) > tol1 then
c_b = c_b + c_d
else
c_b = c_b + math.abs(tol1)*math.sgn(xm)
end_if
c_x0 = c_b
fb = func(c_x0)
end_loop
toto = io.out(' froot exeeding maximum iteration')
froot = c_b
end

fish def store_root
global val_root = -math.pi
global c_eps2   = 1.e-2
global tol      = 1.e-6
global c_int    = math.pi*0.5
loop local ii (1,nroot)
local val   = float(ii-1)*math.pi
global c_eps1 = c_eps2*(val + c_int - val_root - math.pi)
if ii > 50 then
c_eps1 = c_eps1*10.
end_if
local c_x1  = val - c_eps1
local c_x2  = val + c_int - c_eps1
val_root = froot(c_x1,c_x2,tol)
table.x(tabroot,ii) = ii
table.y(tabroot,ii) = val_root
local sr  = math.sin(val_root)
local cr  = math.cos(val_root)
local den = val_root-sr*cr
local r2  = val_root * val_root
table.x(tab1,ii)    = ii
table.y(tab1,ii)    = sr/den
table.x(tab2,ii)    = ii
table.y(tab2,ii)    = cr/den
table.x(tab3,ii)    = ii
table.y(tab3,ii)    = sr*cr/den
table.x(tab4,ii)    = ii
table.y(tab4,ii)    = r2
table.x(tab5,ii)    = ii
table.y(tab5,ii)    = r2*sr*cr/den
end_loop
end
[store_root]

; reaction force
;[global pnt2 = list(gp.list)(gp.isgroup(::gp.list,"Top"))]

fish def rf
local sum = 0.0
loop foreach local pnt list(gp.list)(gp.isgroup(::gp.list,"Top"))
sum = sum+gp.force.unbal.y(pnt)
endloop
rf = sum
end
; outflow
[global pnt3 = list(gp.list)(gp.isgroup(::gp.list,"East"))]
fish def qout
local outflow = -1.0*list.sum(gp.flow(::pnt3))
qout = outflow
end
; velocity per mechanical step
fish def ana_yv
global velo  = 0.0
global veln  = 0.0
local jjsav = 0
loop local jj (1,nroot)
local al2  = table.y(tab4,jj)
veln = velo + table.y(tab5,jj)*math.exp(-al2*coe*	fluid.time.total)
if  veln = velo then
jjsav = jj
jj    = nroot
end_if
velo = veln
end_loop
ana_yv = veln*v_co*fluid.timestep
end
; pore pressure profile
fish def check_pp(taba)
global pcoe = 2.*c_F*c_sk*(1.+c_un)/(3.*c_a)
local tabn = taba + 10
local id =0
loop foreach local ip gp.list

global sumo  = 0.
global sumn  = 0.
local jjsav = 0
local xval  =gp.pos.x(ip)
local yval  =gp.pos.y(ip)
if yval= 0.05 then
loop local jjp (1,nroot)
local al2  = table.y(tab4,jjp)
local al   = table.y(tabroot,jjp)
local term = math.cos(al*xval/c_a)-math.cos(al)
sumn =sumo+table.y(tab1,jjp)*term*math.exp(-al2*coe*fluid.time.total)
if  sumn = sumo then
jjsav = jjp
jjp    = nroot
end_if
sumo = sumn
end_loop
id=id+1
table.x(taba,id) = xval
table.y(taba,id) = sumn*pcoe
table.x(tabn,id) = xval
table.y(tabn,id) = gp.pp(ip)
endif
end_loop
end
; pore pressure history
fish def ana_pp
global sumo  = 0.
global sumn  = 0.
local jjsav = 0
local pnt    = gp.near(0.,0.1)
local xval  =  gp.pos.x(pnt)
loop local jjp (1,nroot)
local al2  = table.y(tab4,jjp)
local al   = table.y(tabroot,jjp)
local term = math.cos(al*xval/c_a)-math.cos(al)
sumn = sumo + table.y(tab1,jjp)*term*math.exp(-al2*coe*fluid.time.total)
if  sumn = sumo then
jjsav = jjp
jjp    = nroot
end_if
sumo = sumn
end_loop
ana_pp = sumn*pcoe
global num_pp = gp.pp(pnt)
end
; degrees of consolidation history
fish def deg_cons(taba)
local ccoe  = 2.*(1.-c_un)/(c_un-c_nu)
local  uxa0  = c_F*c_un/(2.*c_sh)
local  duxa  = c_F*(c_nu - c_un)/(2.*c_sh)
local  uyb0  = -c_F*c_b*(1. - c_un)/(2.*c_sh*c_a)
local  duyb  = -c_F*c_b*(c_un - c_nu)/(2.*c_sh*c_a)
local  taba1 = taba+1
local  taba2 = taba+2
loop local ip (1,tbi)
local sumo  = 0.
local sumn  = 0.
local jjsav = 0
local tt = table.x(3,ip)
loop local jjp (1,nroot)
local al2  = table.y(tab4,jjp)
sumn = sumo + table(tab3,jjp)*math.exp(-al2*coe*tt)
if  sumn = sumo then
jjsav = jjp
jjp    = nroot
end_if
sumo = sumn
end_loop
local xval = math.log(tt)
table.x(taba,ip)  = xval
table.y(taba,ip)  = 1.-ccoe*sumn
table.x(taba1,ip) = xval
table.y(taba1,ip) = (table.y(4,ip)-uxa0)/duxa
table.x(taba2,ip) = xval
table.y(taba2,ip) = (table.y(3,ip)-uyb0)/duyb
end_loop
end

fish def big_step(nfstep)
local val_yv = 0.0
local rf1 =0.0
loop local kk (1,nfstep)
command
model fluid active on mech active off
model step 1
model fluid active off mech active on
end_command
val_yv = ana_yv
rf1 =rf
command
zone gridpoint initialize velocity-y [val_yv] range group 'Top'
model step 1
zone gridpoint initialize velocity-y 0.0 range group 'Top'
model solve ratio 1.e-5
end_command
tbi = tbi + 1
local pnt43 = gp.find(43)
local pnt41 = gp.find(41)
table.x(1,tbi) = fluid.time.total
table.y(1,tbi) = rf1              ; reaction force
table.x(2,tbi) = fluid.time.total
table.y(2,tbi) = val_yv  ; velocity increment per step
table.x(3,tbi) = fluid.time.total
table.y(3,tbi) = gp.disp.y(pnt43) ; top y-displacement
table.x(4,tbi) = fluid.time.total
table.y(4,tbi) = gp.disp.x(pnt41) ; right x-displacement
end_loop
end