40 6 Numerical results ‘The equations governing poroelasticity form a compli- cated coupled system of partial differential equations, ‘and hence, there are very few analytical solutions avai able, However, there are some known solutions that provide an ideal opportunity to examine the numerical ‘accuracy of the algorithms presented herein, In this Paper, We turn to the well-known case of Mandel’s problem. This problem gives us a chance to test our algorithm in a two-dimensional setting, and it also pro- vides a nice example of the dynamics involved in a fuid-solid interaction. 6.1 Mandel’s problem Mande'’s problem is important because it admits an an- alytical solution in two dimensions on a finite domain. ‘The original paper by Mandel [19] presented only an ‘analytical form for the pressure; later, Abousleiman [1] extended the results to include analytical expressions for the displacement and stress. It is therefore an ex- cellent model to verily the accuracy of a poroelasticity algorithm, Mandel’s problem is also fascinating from the point of view that it clearly illustrates that sotid— fluid interactions can lead to unexpected behavior, and thus highlights the need for poroelasticity theory in practice. For a precise formulation of the problem, we con- sider a poroelastic slab of extent 2a in the x direction, 2b in the y direction, and infinitely long in the z di- rection. The slab is sandwiched in between two rigid plates as shown in Fig, 1. At time = 0, a downward surface force of magnitude 2F is applied to the top uy =0 Fig. Mandel’ problem D sping Comput Geosei (2007) 11131144 plate, and an equal but upward force is applied to the bottom plate. Because the plates are considered rigid, we must impose the additional constraint that the vertical displacements at the top and bottom, respec- tively, are uniform. This requirement ensures that the slab remains in contaet with the plates. We incorporate this constraint into our methods by use of Lagrange ‘multipliers. In addition, at all times, the slab at x = +a remains drained, thus leading to the boundary condi- tion p(cka, y,1) = 0. ‘The application of a load causes an instantaneous and uniform pressure increase throughout the do- ‘main; the theory predicts this to be p* = lim p(x.0) = 4,B( + ¥4)F. Also, the x displacement depends only con x and f, and the y displacement depends on y and ¢ We can thus define the following x and y settlements: S:() = Use = 4,1) and sy =uj(y=b,0). It is. pre dicted that there will be instantaneous and asymptotic settlements, So = limse() my Fv a" =FU= wb Sy His (9 = Fw 50 -- AL this point, it is useful to remark that the problem is symmetric about the x and the y axes, and is indepen- dent of the z direction. Therefore, we may reduce the computational domain to only the upper-right quadrant of the xy plane, as illustrated in Fig. 1. Thus, we add to the governing equations 2.3a and 2.3b the following ‘boundary and initial conditions on the domain (0, a) x (0,5): n=O, ¥=0; ty =0, y=0; tyv=-2F, y=b, =Ootherwise: p=0, x=0 BERR PHO x20 y=0. y u(x, y =) = constant (Rigid Plate Constraint) p=u=0, 1=0; We also remark that the constant in the rigid plate constraint is time-dependent.Comput Geosc (2007) 1113-144 1a We refer now to the notation in Appendix Al as wwe list the analytical series solutions for the pressure, displacement, and stress as found in Abousleiman [1]: — 2H S x (cos “2 — cosa) exp(—azeit/a*) x exp(-afedt/a") x FR cosay FO) , FU= v9) snag cose ane pe xexp(-azeat/a")]y, SiN. C086 ooo ott fay ‘We note that all other components of the stress tensor ake 710, C44 = Ory = 0, We note that, for the above 2. Mandel’s pressure solution a boundary layer ast 0 +f tin S ept-ageit/a%y, 8 jg tpt? d= (11) tet Fig. 3 Mandel’ problem. The pressure integrals ||p.(0)\f and Pex( 5 equations, aj, represents the positive solutions to the nonlinear equation ny and must be solved for numerically, 6.11 Results Recall that the optimal error estimates presented herein stipulate a sufficiently smooth analytical solu- tion, However, full-order convergence for even the lowest-order methods requires that p « L*(#°). One can see why by examining the proof of the error es- timates given earlier. Assuming suflicient regularity implied that the following integrals are finite: f JEp(eyl lade = where kis the order of the mixed space. However, even fork =0, if p ¢ LAUP), then J IIa(o)|R'de = 00. Indeed, for Mandel’s problem, Fig. 2 shows the true pressure solution for a typical choice of physical para- meters. As {> 0, a very large gradient occurs. Natu- rally, this leads one to suspect that the regularity of the pressure might be limited. This is confirmed by looking at Fig. 3, which shows the regularity of the pressure by computing [lp (OIlj and lips (ll. Its discovered that f leCOlibyade <2, lipeOlig = OC) (61) [pels OW") (62) D springer1 Comput Geosei (2007) 11131144 ‘The order approximation (6.1) implies that pe 12(H!), but Eq. 6.2 shows that p ¢ L?(H° ‘Thus, it is reasonable to assume that p € L2(I1"*) for some s¢ [0, 1). Furthermore, it is precisely this value of s that limits the best rate of convergence that ‘one might expect from the CG/Mixed algorithm for Mandel’s problem. For numerical verification of our algorithms, we use a force F = 2,000, Standard continuous linear ele- ments are used to approximate the displacement, and the lowest Raviart-Thomas space is used for the flow variables, The backward Euler scheme is used. The top image in Fig. 4 shows how the computed pressure error, lIp ~ Bllz~(22), varies with element size /. The bottom image shows the displacement error, {ju — @lli~u- ‘To minimize the effects of the error produced by the time discretization, a small time step of Ar= le-6 ate 10% 0.500117, 4.0000 hh element size ate 10% uw “ o.so0202 e = ~olment sae Fig. 4 LCL?) pressure error (Wop) and L™(H') displacement error (bottom) in Mandel's problem D sping is chosen. The convergence rate is determined by the slope of the logarithm of the error (= 0.500117). ‘The convergence rate (~ 0.5) for both the pressure and displacement conforms to the expectation that the lack of regularity in the pressure solution would degrade performance. From the convergence rate and the theorem on error estimates, one might be led to conclude that p € L2(H!**) for some small 0 < ¢<<1. For a larger time step, the finite element solution might exhibit a better convergence rate because the large pressure gradient subsides over time, However the time-diseretization error can become problematic for those larger time steps. See Phillips and Wheeler [22] for the analysis of the discrete-in-time case and an example using Terzaghi’s consolidation problem, which has a pressure solution profile similar to that for Mandel’s problem. ‘An interesting aspect of Mandet’s problem is shown in Fig. 5 — the increase in pressure (above the initial increase) near the center of the medium. This is known as the Mandel-Crever effect and is a phenomenon par- ticular to poroelastic material. The increase occurs be- cause the deformation and rigid plate condition create an effect similar to that ofa source term in the pressure equation. Indeed, Coussy [6] shows that one ean reduce the diffusion equation to one involving the normalized pressure and a source as expl-a3i). (63) ‘The source term is time-dependent only (indepen- dent of x), and it can be quite large at early time: The interesting behavior has been confirmed exper imentally in the work of Gibson et al. [11] and Verruijt [32} Appendix A1: Additional notation and terms in the literature For completeness, we present here some additional terms and notations commonly found in poroelastic literature. Though not essential for the development of our model and the subsequent error estimates, fami iarity with this additional nomenclature will be helpful ‘when discussing numerical results. The first of the additional terms can be written in terms of the Lamé coefficients and 3. and are common in the field of linear elasticity. They are, respectively,Comput Geosc (2007) 1113-144 Big.§ The Mandel-Creyer effect is captured by the approximate solution the skeleton bulk modulus K, Young's modulus E, and Poisson's coefficient v _ ae marews 2GK+u) There are also “undrained” versions of the above coef ficients and are denoted by the subscript u: Ky, Ey. Ys Additionally, Ky has the following relation, Kaka With the above, two useful terms have been defined: the fluid diffusivity coefficient c, and Skempton's coeff cient B: 1 Kee TG Kot ia a ak Clearly, ccis usefully defined only for the ease where the permeability « is a constant, Skempton’s cocificient can be used to reveal a useful relation between Poisson’s coeflicient and its undrained version, References 1, Abousleiman, ¥., Cheng, A-D., Cui, L., Detournay, E. Roepiers, J: Mandel’s problem revisited. Géotechnique 46, 187-195 (1996) 2, Biot, M: General theory of three-dimensional cons J. Appl Phys. 2, 155-164 (1941) 3, Biot, Mz Theory of elasticity and consolidation for @ porous anisotropic media. J. Appl. Phys. 26(2), 182-185 (1955) 4, Boyee, W., DiPrima, Ru Elementary Differential Equations, ‘ith edn. Wiley, New York (2001) 5, Brenner, S., Scott, L: The Mathematical Theory of Finite Element Methods. Springer, Berlin Heidelberg New York (4394) 6, Coussy, 0. Poromechanies. Wiley, New York (2008) 7. Detoutnay, B,, Cheng, A.-D.: Poroelastie response ofa bore hole in non-hydrostati tress field Int J. Rock Mech, Mining Sei 28, 171-182 (1988) 8, Evans, L. Partial Dlferential Equations. American Mathe- ‘matical Society, Providence (1998) 9, Gai, X: A Coupled Geomechanics and Reservoir Flow Model on Parallel Computers. Pa. thesis, University of “Texas at Austin (2008) 10, Garagash, D., Detoumay. E: An analysis ofthe inuence of, the pressurization rate on the borchole breakdown pressure 1 Solids Struct. $4(2), 3099-3118 (1997) D springerwae Comput Geose (2007) 11131144 11, Gibson, R., Knight, K, Taylor, P= A critical experiment to examine theories of three-cimensional consolidation. Eu Cont Soil Mech. 1, 69-73 (1963) 12, Hudson, J., Stephansson, O,, Andersson, J, Tsang, C-F Ling, L! Coupled T-H-M issues related to radioactive waste repository design and performance. Int, J. Rock Mech, Mining Se. 38, 13-161 (2001) 13, Kim, J-M, Parizek, R: Numerical simulation of the noord bergum effect resulting from groundwater pumping ina lay ered aquifer system, J, Hydrol 202, 251-243 (1997) 14, Langtord, T! Northwest Houston Sinking Faster than Coastal Areas. ttp:iwww [Link] (1997) 15, Lewis, RIW., Schrefler, B.A: The Finite Element Method in the Deformation and Consolidation of Porous Media, Wiley ‘& Sons, Chichester (1987) 16, Lin, T: An Interview with Fan Lin, PRD. hitp/www. cs-topic¢.comlearthquakes/interviews/TanLin htm] (2003) 17, Lis, R: Discontinuous Galerkin finite element solution for oromechanies, PhD. thesis, University of Texas at Austin (2004) 18, Lubick, N Modeling complex, multiphase porous media sys- tems, SIAM News $(3) (2002) 19, Mandel, 1: Consolidation des sols (Gude mathématique), Géotechnique 30, 287-209 (1953) 20, Nedelee, 1: Mixed finite elements in R*. Numer. Math, 38, 315-341 (1980) 21, Philips, PJ Finite clement methods for linear poroelas- Licty: theoretical and computational results, Ph.D. thesis, University of Texas at Austin (2005) 22, Philips, PJ, Wheeler, MI. : A. coupling of mixed and continuous Galerkin nite element. methods for poroclasticity Il: the discrete-in-time case. Comput. Geoec {oi:10 10071810596.007- 9044-2 23, Rajapakse, Rs Stress analysis of borehole in poroelastic redhuen. J Eng. Mech, 19(6), 1205-1227 (1993) 24, Raviart, RA, Thomas, JW.: A mixed finite element method for 2nd order elliptic problems. In: Mathematieal Aspects of the Finite Element Method. Lecture Notes in Mathematics, vol. 606, pp. 292-315, Springer, Berlin Heidelberg New York 977) D springer 26, 2. 23, 28, 30. 31 33 34 3 36, a. 38. Rivire, B., Wheeler, M: Optimal error estimates applied to linear elasticity, Technieal report, ICES Report. ICES, Austin (2000) Roose, T, Nett, P., Munn, L, Boucher, Y., Jai, R. Solid stress generated by epheroid growth estimated using 3 lineax poroelastic model. Microvase: Res. 66, 204-212 (2003) Rutgvst, J, Tsang, C-F: Analysis of thermal-hydrologic~ ‘mechanical behavior near an emplacement drift at Yueea mountain I. Contam. Hydrol. 62-68, 637-852 (2003) Showalter. RE: Diffusion in poro-elastic media. J. Math, Anal. Appl. 281, 310-340 (2000) Smilli, A., Sobey, L, Molnar, Z.: A bydro-clastic model of hydrocephalus. Technical report, Oxford University Comput ing Laboratory: Numerical Analysis Group, Oxford (2004) Swan, C, Lakes, R., Brang, R, Stewart, Ks Micromechan: ically based poroclastic modeling of Muid flow in Maversi bone. Biomech, Eng, 125(1), 25-87 (2003) Terzaghi, K: Theoretical Soil Mechanics. Wiley, New York (sas) Verruijt. Az Discussion. Proc. 6th Int. Cont. Soil Mech. 3, 401-40 (1965) Wang, ILE: Theory of Linear Poroclasticity with Ap- plications to Geomechanies and Hydrogeology. Princeton University Press, Princeton (2000) ‘Wang, ¥,, Dusseault, M: A coupled conduetive-convective thermo-poroelastic colution and implications for wellbore stability J. Petrol. Sei. Eng, 38, 187-198 (2003) ‘Weiser, A., Wheeler, M: On convergence of block-centered Maite differences and elliptic equations. SIAM J. Numer. Anal. 251(2), 351-375 (1988) Wheeler, MIF A priori L? error estimates for Galerkin approximations to parabolic partial differential equations. SIAM. Numer, Anal 10, 723-759 (1973) Wheeler, ME, Philips, PJ: A coupling of mixed and Galerkin methods for pore-clasticty. In: Proceedings of the Second MIT Conference on Computational Fluid and Solid ‘Mechanics, MIT, Cambridge, 17-20 June 2003 Zenisek, A: The existence and uniqueness theorem in Biot’s consolidation theory. Aplik. Matem. 29, 194-210 (cases)