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.

Figure 1: Mandel’s problem.

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.

Figure 2: Pore pressure profile comparison : \(\hat{p}\) vs. \(\hat{x}\).

Figure 3: Pore pressure versus consolidation time at the center of the sample.

Figure 4: History of y-reaction force on top platen.

Figure 5: Degree of consolidation versus log time.

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

[1]

To view this project in FLAC2D, use the program menu.

Help ▼ Examples…

  ⮡   FLAC
    ⮡   VerificationProblems
      ⮡   MandelsProblem
        ⮡   Mandels-Problem.prj