Published online at http://imechanica.org/node/17096
Numerical modeling of thermohydromechanical coupling processes in porous media Thermohydromechanics (THM) is a branch of mechanics aimed to predict how deformable porous media behave, while heat transfer and fluid transport simultaneously occur in the pores filled by liquid and/or gas. Understanding these multiphysical responses is important for a wide spectrum of modern engineering applications, such as tissue scaffolding, geothermal heating, mineral exploration and mining, hydraulic fracture, energy piles, tunneling with frozen soil and nuclear waste storage and management. The major difficulty encountered for modeling this coupling physical processes is that the thermal and porefluid diffusions, and the deformation of solid skeleton may take place in profoundly different spatial (micron vs. kilometer) and temporal scales (second vs decades). Bundling all the physical processes in simulations while maintaining accuracy and robustness is therefore a difficult task.In the last three decades, a considerable progress has been made for deriving mathematical theories and implementing computer models to replicate the fully coupled thermohydromechanical processes. This journal club article will focus on the derivation and implementation of the thermohydromechanical model for numerical simulations. Any feedback or critique are greatly appreciated. 1. Monolithic vs. operatorsplitting schemes The governing equation of the thermohydromechanical problem can be obtained from the balance principles. The displacement, pore pressure, Darcy’s velocity and temperature are usually the primary solution field stored at the nodal points. If the whole set of governing equations are solved simultaneously, then the solver is called monolithic. The monolithic solver is particularly important for strongly coupling problems in which passing parameters among fluid and solid solvers may not be sufficient (Preisig & Prevost, 2011).For instance, a monolithic small strain finite element code, FRACON, has been developed by (Nguyen & Selvadurai,1995). In this code, the balance of linear momentum and mass are fully coupled, while thermal transport may affect the solid deformation and porefluid diffusion but not vice versa. Li, Liu and Lewis, introduce a corotational FEM formulation and incorporate plasticity into THM model to model the nonisothermal elastoplastic responses of porous media at large strain in (Li,2005). Recent work by Preisig and Prevost employed a fully coupled THM simulator to compare the numerical solutions against the field data in a case study for carbon dioxide injection at In Salah, Algeria (Preisig & Previst, 2011). Kolditz et al., 2012 introduces an opensource project OpenGeoSys, which takes advantage of an objectorient framework and provides software engineering tools such as platformindependent compiling and automated benchmarking for developers. In addition to the monolithic finite element scheme, attempts have been made to sequentially coupled multiphase flow and geomechanical simulators by establishing proper feedback and information exchange mechanisms. One such example is TOUGHFLAC, which links flow simulator TOUGH2 with a small strain finite difference code FLAC (Rutqvist,2011). Klar et al. 2013 employ an explicit timemarching scheme to simulate the fully coupled thermal, and porefluid diffusions and the path dependent responses of the solid skeleton. This sequential coupling approach is an attractive alternative to the monolithic approach, as it is easier to implement and maintain flow and solid simulators separately. In addition, the operatorsplitting approach, if used properly, can offer tremendous saving on memory and CPU times. The operator splitting approach also enables one to assign different time steps to different physical processes. For instance, in reactivediffusion simulations, it is common to have the reaction component updating with a several order finer time step than the diffusion counterpart. However, proper communication must be established to ensure the correctness and numerical stability of numerical solutions (Jha & Juanes, 2007, Preisig, 2011, Klar et al. 2013, Sun et al., 2013a, Sun et al. 2013b). As noted in (Jha & Juanes,2007), the sequential coupling scheme used to link the fluid and solid simulators may have profound impact on the efficiency, stability and accuracy of the numerical solutions (Schrefler et al.,1995, Schrefler et al.,1997). If the fluid and solid simulators use different grids or meshes (eg. finite volume for fluid and finite element for solid), then a proper data projection scheme is required to transfer information from Gauss points and nodes of the solid mesh to the fluid mesh and vice versa (Goumiri & Previst, 2011). For large scale parallel simulations, the sequential couplings must be carefully designed to avoid causing bottleneck due to the difference in solver speed. This can be a significant problem if either the solid or fluid solver runs only in serial. 2. Mixed finite elements (e.g. TaylorHood, RaviartThomas) vs. stabilized procedures As noted in Zienkiewicz et al., 1999, Wan, 2002, Mira et al., 2003, Truty & Zimmermann, 2006, Simoni et al., 2008, White & Borja, 2008, Preisig & Prevost , 2011, Sun et al, 2013a, Sun et al. 2013b Borja, 2013, numerical stability is often a major challenge for poromechanics models. Due to the lack of infsup condition (Babuvska, 1973, Brezzi et.al, 1985, Bathe, 2001, Bochev et al. 2006), pore pressure and temperature fields may exhibit spurious oscillation patterns and/or checkerboard modes if the displacement, pore pressure and temperature are spanned by the same set of basis function. While these spurious oscillations are less severe at the drained/isothermal limit, they may intensify when small time step is used or when materials are near undrained/adiabatic limit. From a mathematical viewpoint, these nonphysical results are due to the kernel (null space) of the discrete gradient operator being nontrivial. This nontrivial kernel makes it possible to have certain spatially oscillating pore pressure and temperature fields that have no impact on the solid deformation in a numerical simulation. To cure the numerical instability due to the lack of infsup condition, previous researches have established a number of techniques that employ different basis functions to interpolate displacement and pore pressure and obtain stable solutions. For instance, Zienkiewicz and coworkers Zienkiewicz et al. 1999, and Borja & Alarcon 1995 used TaylorHood finite element (quadratic basis functions for displacement and linear basis functions for pore pressure) to satisfy infsup condition and maintain numerical stability for isothermal hydromechanics problems. On the other hand, Jha & Juanes 2007 have shown that linear displacement combined with pore fluid velocity in the lowestorder RaviartThomas space, and piecewise constant pore pressure may also lead to stable solutions for isothermal poromechanics problems. Nevertheless, infsup stable mixed finite element models require multiple meshes for displacement, pore pressure and/or fluid velocity. As a result, they require additional programming effort to pre and postprocessing data and maintain the more complex data structure. To avoid the complications of using multiple meshes for each solution field, an alternative is to use one single equalorder finite element mesh with stabilization procedures. Many stabilization procedures have been proven to be able to eliminate the spurious oscillation modes without introducing extra diffusion for small strain isothermal poromechanics problems. For instance, White & Borja 2008 employed a polynomial projection scheme originated from Dohrmann & Bochev, 2004 to simulate slip weakening of a fault segment. A numerical example in this paper shows that spurious oscillation may persist in the infsup stable finite elements (e.g. quadraticdisplacement/linearporepressure pair) if the diffusivity is very low for a given time step size. Perhaps the major drawback for the stabilized finite element approach is that it is often necessary to find the stabilization parameter that just give the right amount of stabilization effect without overdiffusing the solution. Recently, my collaborators and I have attempted to address this issue for isothermal poromechanics problem in the large deformation regime, where the stabilization term is adaptively adjusted to avoid overdiffusion (Sun 2013a and Sun et al. 2013b). The stabilization term is derived by solving a small strain onedimensional deformationdiffusion problem. By determining the critical value of the diffusivity that leads to complex valued growth/decay rate, the optimal value for the stabilization parameter can be estimated. For the threefield thermohydromechanical problem, the situation is more complicated. Since the thermal convectiondiffusion and pore fluid diffusion may reach steadystate in profoundly different rates, it is unclear how to stabilize both the pore pressure and temperature fields without overdiffusing one or both of them. 3. Length scale and bifurcations Unlike singlephase materials, porous media are multiphase in which solid, liquid(s) and gas(s) can all occupy fractions of volume of the representative elementary volume. As a result, the deformation of the solid skeleton are influenced by the heat transfer and pore fluid diffusions. One direct consequence is that the diffusion will introduce rate dependence to the mechanical responses. This rate dependence is sometime considered by some as a localization limiter which regularizes the governing equations when deformation band formed in numerical simulations. Zhang, et al (1999) analyzed the characteristic equation of a onedimensional wave propagation problem and found that a length scale may be derived for compressive wave even when softening occurs for solid skeleton, but this is not the case for shear wave. Abellan & de Borst (2005) conducted a similar dispersion analysis and found that the physical length scale disappear in short wavelength limit. This result is supported by numerical simulations in which the strain of a softening bar composed of saturated porous media is found to be mesh dependence. In other words, using diffusion alone as a mean to circumvent mesh dependence is not productive, according to both Zhang et al (1999) and Abellan & de Borst (2005). Nevertheless, a mesh independent postbifurcation response may still be obtained in numerical simulations, if a length scale can be introduced through other means (e.g. grain rotation, gradient dependence, nonlocal plasticity or damage model, rate dependent solid constituent…etc). A similar observation is made by Bazant (2010) for models that couple multiple spatial scales. The author observed that incorporating subscale simulations to macroscopic simulations alone does not introduce length scale for the softening materials. This treatment merely move the strain localization phenomenon one scale down. Recent multiscale discrete element/finite element models proposed by Guo & Zhao (2014) and by Nguyen et al. (2014) both confirm that the macroscopic finite element responses are mesh dependence in the postbifurcation regime even though the macroscopic responses are inferred from grainscale simulations. Reference 1.Abellan, M. A., & De Borst, R. (2006). Wave propagation and localisation in a softening twophase medium. Computer Methods in Applied Mechanics and Engineering, 195(37), 50115019. 2.Armero F (1999) Formulation and finite element implementation of a multiplicative model of cou pled poroplasticity at finite strains under fully saturated conditions. Computer methods in applied mechanics and engineering 171(3):205–241 3.Athy LF (1930) Density, porosity, and compaction of sedimentary rocks. AAPG Bulletin 14(1):1–24 4.Auricchio F, Beir ̃ao da Veiga L, Lovadina C, Reali A (2005) A stability study of some mixed finite elements for large deformation elasticity problems. Computer methods in applied mechanics and engineering 194(9):1075–1092 5.Auricchio F, da Veiga LB, Lovadina C, Reali A, Taylor RL, Wriggers P (2013) Approximation of incompressible large deformation elastic problems: some unresolved issues. Computational Mechanics pp 1–15 6.Babuˇska I (1973) The finite element method with lagrangian multipliers. Numerische Mathematik 20(3):179–192 7.Bathe KJ (2001) The infsup condition and its evaluation for mixed finite element methods. Com puters and Structures 79(2):243 – 252, DOI 10.1016/S00457949(00)001231 8.Bazant, Z. P. (2010). Can multiscalemultiphysics methods predict softening damage and structural failure?. International Journal for Multiscale Computational Engineering, 8(1). 9.Bear J (1972) Dynamics of fluids in porous media. Elsevier Publishing Company, New York, NY 10.Biot M (1941) General theory of three dimensional consolidation. Journal of Applied Physics 12(2):155 –164 11.Bochev PB, Dohrmann C, Gunzburger M (2006) Stabilization of loworder mixed finite elements for the stokes equations. SIAM J Numer Anal 44(1):82–101 12.Borja RI (2013) Plasticity Modeling and Computation. Springer, Berlin 13.Borja RI, Alarco ́n E (1995) A mathematical framework for finite strain elastoplastic consolidation part 1: Balance laws, variational formulation, and linearization. Computer Methods in Applied Me chanics and Engineering 122(1):145–171 14.Borja RI, Tamagnini C, Alarc ́on E (1998) Elastoplastic consolidation at finite strain part 2: Finite element implementation and numerical examples. Computer Methods in Applied Mechanics and Engineering 159(1):103–122 15.Brezzi F, Fortin M (1991) Mixed and hybrid finite element methods. SpringerVerlag New York, Inc. 16.Brezzi F, Douglas J, Marini LD (1985) Two families of mixed finite elements for second order elliptic problems. Numerische Mathematik 47:217–235 15.Broccardo M, Micheloni M, Krysl P (2009) Assumeddeformation gradient finite elements with nodal integration for nearly incompressible large deformation analysis. International journal for numerical methods in engineering 78(9):1113–1134 16.Callari C, Armero F (2004) Analysis and numerical simulation of strong discontinuities in finite strain poroplasticity. Computer Methods in Applied Mechanics and Engineering 193(27):2941–2986 17.Castellazzi G, Krysl P (2012) Patchaveraged assumed strain finite elements for stress analysis. International Journal for Numerical Methods in Engineering 90(13):1618–1635 18.Coussy O (2004) Poromehcanics 19.Coussy O (2007) Revisiting the constitutive equations of unsaturated porous solids using a lagrangian saturation concept. International Journal for Numerical and Analytical Methods in Geomechanics 31(15):1675–1694 20.Cowin SC, Doty SB (2009) Tissue mechanics. Springer 21.Dohrmann CR, Bochev PB (2004) A stabilized finite element method for the stokes problem based on polynomial pressure projections. International Journal for Numerical Methods in Fluids 46(2):183– 201 22.Girault V, Raviart PA (1986) Finite element methods for NavierStokes equations, Theory and algorithms, volume 5 of Springer Series in Computational Mathematics. Springer, Berlin 23.Goumiri IR, Prevost JH (2011) Cell to node projections: An assessment of error. International Journal for Numerical and Analytical Methods in Geomechanics 35(7):837–845, DOI 10.1002/nag.927, URL http://dx.doi.org/10.1002/nag.927 24.Guo, N., & Zhao, J. (2014). A coupled FEM/DEM approach for hierarchical multiscale modelling of granular media. International Journal for Numerical Methods in Engineering, 99(11), 789818. 25.Heroux MA, Bartlett RA, Howle VE, Hoekstra RJ, Hu JJ, Kolda TG, Lehoucq RB, Long KR, Pawlowski RP, Phipps ET, et al (2005) An overview of the trilinos project. ACM Transactions on Mathematical Software (TOMS) 31(3):397–423 26.Hiroshi H, Minoru T (1986) Equivalent inclusion method for steady state heat conduction in com posites. International Journal of Engineering Science 24(7):1159–1172 27.Holzapfel GA (2000) Nonlinear solid mechanics: a continuum approach for engineering 28.Howell JS, Walkington NJ (2011) Inf–sup conditions for twofold saddle point problems. Numerische Mathematik 118(4):663–693 29.Jha B, Juanes R (2007) A locally conservative finite element framework for the simulation of coupled flow and reservoir geomechanics. Acta Geotechnica 2:139–153, DOI 10.1007/s1144000700330, URL http://dx.doi.org/10.1007/s1144000700330 30.Jing L, Tsang CF, Stephansson O (1995) Decovalexan international cooperative research project on mathematical models of coupled thm processes for safety analysis of radioactive waste repositories 32(5):389–398 31.Karrech A, Poulet T, RegenauerLieb K (2012) Poromechanics of saturated media based on the logarithmic finite strain. Mechanics of Materials 51(0):118 – 136 32.Kolditz O, Bauer S, Bilke L, Bo ̈ttcher N, Delfs J, Fischer T, Go ̈rke U, Kalbacher T, Kosakowski G, McDermott C, et al (2012) Opengeosys: an opensource initiative for numerical simulation of thermohydromechanical/chemical (thm/c) processes in porous media. Environmental Earth Sci ences 67(2):589–599 33.Lewis R, Majorana C, Schrefler B (1986) A coupled finite element model for the consolidation of nonisothermal elastoplastic porous media. Transport in Porous Media 1(2):155–178 34.Li X, Liu Z, Lewis RW (2005) Mixed finite element method for coupled thermohydromechanical process in poroelastoplastic media at large strains. International Journal for Numerical Methods in Engineering 64(5):667–708 35.Liu R, Wheeler M, Dawson C, Dean R (2009) Modeling of convectiondominated thermoporome chanics problems using incomplete interior penalty galerkin method. Computer Methods in Applied Mechanics and Engineering 198(9):912–919 36.McTigue D (1986) Thermoelastic response of fluidsaturated porous rock. Journal of Geophysical Research 91(B9):9533–9542 37.Mira P, Pastor M, Li T, Liu X (2003) A new stabilized enhanced strain element with equal order of in terpolation for soil consolidation problems. Computer methods in applied mechanics and engineering 192(37):4257–4277 37.Moran B, Ortiz M, Shih C (1990) Formulation of implicit finite element methods for multiplicative finite deformation plasticity. International Journal for Numerical Methods in Engineering 29(3):483– 514 38.Nguyen T, Selvadurai A (1995) Coupled thermalmechanicalhydrological behaviour of sparsely frac tured rock: Implications for nuclear fuel waste disposal. International Journal of Rock Mechanics and Mining Sciences and Geomechanics Abstracts 32(5):465 – 479 39.Nguyen, T. K., Combe, G., Caillerie, D., & Desrues, J. (2014). FEM× DEM modelling of cohesive granular materials: Numerical homogenisation and multiscale simulations. Acta Geophysica, 62(5), 11091126. 40.Notz PK, Pawlowski RP, Sutherland JC (2012) Graphbased software design for managing complexity and enabling concurrency in multiphysics pde software. ACM Transactions on Mathematical Software (TOMS) 39(1):1 41.Nur A, Byerlee J (1971) An exact effective stress law for elastic deformation of rock with fluids. Journal of Geophysical Research 76(26):6414–6419 42.Pantuso D, Bathe KJ (1997) On the stability of mixed finite elements in large strain analysis of incompressible solids. Finite elements in analysis and design 28(2):83–104 43.Pawlowski RP, Phipps ET, Salinger AG (2012) Automating embedded analysis capabilities and man aging software complexity in multiphysics simulation, part i: Templatebased generic programming. Scientific Programming 20(2):197–219 44.Pawlowski RP, Phipps ET, Salinger AG (2012) Automating embedded analysis capabilities and man aging software complexity in multiphysics simulation, part i: Templatebased generic programming. Scientific Programming 20(2):197–219 45.Pawlowski RP, Phipps ET, Salinger AG, Owen SJ, Siefert CM, Staten ML (2012) Automating em bedded analysis capabilities and managing software complexity in multiphysics simulation, part ii: Application to partial differential equations. Scientific Programming 20(3):327–345 46.Postelnicu A (2004) Influence of a magnetic field on heat and mass transfer by natural convection from vertical surfaces in porous media considering soret and dufour effects. International Journal of Heat and Mass Transfer 47(6):1467–1472 47.Preisig M, Pr ́evost JH (2011) Coupled multiphase thermoporomechanical effects. case study: Co2 injection at in salah, algeria. International Journal of Greenhouse Gas Control 5(4):1055 – 1064, DOI 10.1016/j.ijggc.2010.12.006 48.Prevost JH (1982) Nonlinear transient phenomena in saturated porous media. Computer Methods in Applied Mechanics and Engineering 30(1):3 – 18 49.Radovitzky R, Ortiz M (1999) Error estimation and adaptive meshing in strongly nonlinear dynamic problems. Computer Methods in Applied Mechanics and Engineering 172(1):203–240 50.Rajagopal KR, Tao L (1995) Mechanics of mixtures, vol 754. World scientific Singapore 51.Rice J, Cleary M (1976) Some basic stress diffusion solutions for fluidsaturated elastic porous media with compressible constituents. Rev Geophys 14(2):227–241, DOI 10.1029/RG014i002p00227 52.Rutqvist J (2011) Status of the toughflac simulator and recent applications related to coupled fluid flow and crustal deformations. Computers and Geosciences 37(6):739 – 750 53.Salinger AG, Pawlowski RP, Phipps ET, Bartlett RA, Hansen GA, Kalashnikova I, Ostien JT, Sun W, Chen Q, Mota A, Muller RP, Nielsen E, Gao X (2013) Albany: A componentbased partial differential equation code built on trilinos. ACM Transactions on Mathematical Software 54.Schrefler B (1995) Fe in environmental engineering: coupled thermohydromechanical processes in porous media including pollutant transport. Archives of Computational Methods in Engineering 2(3):1–54 55.Schrefler B, Simoni L, Turska E (1997) Standard staggered and staggered newton schemes in thermo hydromechanical problems. Computer methods in applied mechanics and engineering 144(1):93–109 56.Selvadurai A, Nguyen T (1997) Scoping analyses of the coupled thermalhydrologicalmechanical behaviour of the rock mass around a nuclear fuel waste repository. Engineering Geology 47(4):379– 400 57.Selvadurai A, Suvorov A (2012) Boundary heating of poroelastic and poroelastoplastic spheres. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Science 468(2145):2779– 2806 58.Simo J, Miehe C (1992) Associative coupled thermoplasticity at finite strains: Formulation, numerical analysis and implementation. Computer Methods in Applied Mechanics and Engineering 98(1):41– 104 58.Simo J, Taylor R, Pister K (1985) Variational and projection methods for the volume constraint in finite deformation elastoplasticity. Computer Methods in Applied Mechanics and Engineering 51(1):177–208 59.Simoni L, Secchi S, Schrefler B (2008) Numerical difficulties and computational procedures for thermohydromechanical coupled problems of saturated porous media. Computational Mechanics 43(1):179–189 60.Skempton A (1984) Effective stress in soils, concrete, and rocks. Selected papers on soil mechanics pp 4–16 61.de Souza Neto EA, Peri ́c D, Owen DRJ (2008) Computational Methods for Plasticity. John Wiley & Sons, Ltd 62.Sun W, Andrade J, Rudnicki J (2011) Multiscale method for characterization of porous microstruc tures and their impact on macroscopic effective permeability. International Journal for Numerical Methods in Engineering 88(12):1260–1279 63.Sun W, Andrade JE, Rudnicki JW, Eichhubl P (2011) Connecting microstructural attributes and permeability from 3d tomographic images of in situ shearenhanced compaction bands using multi scale computations. Geophysical Research Letters 38(10):L10,302 64.Sun W, Chen Q, Ostien J (2013a) Modeling the hydromechanical responses of strip and circular punch loadings on watersaturated collapsible geomaterials. Acta Geotechnica pp 1–32, DOI 10.1007/s11440 0130276x 65.Sun W, Ostien JT, Salinger AG (2013b) A stabilized assumed deformation gradient finite ele ment formulation for strongly coupled poromechanical simulations at finite strain. International Journal for Numerical and Analytical Methods in Geomechanics, 37(16):6575. 66.Terzaghi Kv, Rendulic L (1934) Die wirksame flachenporositat des betons. Z O ́ st Ingu ArchitVer 86(1):2 67.Tezduyar TE, Osawa Y (2000) Finite element stabilization parameters computed from element ma trices and vectors. Computer Methods in Applied Mechanics and Engineering 190(3):411–430 68.Truty A, Zimmermann T (2006) Stabilized mixed finite element formulations for materially nonlin ear partially saturated twophase media. Computer methods in applied mechanics and engineering 195(13):1517–1546 69.Wan J (2002) Stabilized finite element methods for coupled geomechanics and multiphase flow. PhD thesis, University of Texas at Austin. 70.White JA, Borja RI (2008) Stabilized loworder finite elements for coupled soliddeformation/fluid diffusion and their application to fault zone transients. Computer Methods in Applied Mechanics and Engineering 197(49):4353–4366 71.Wriggers P, Reese S (1996) A note on enhanced strain methods for large deformations. Computer Methods in Applied Mechanics and Engineering 135(3):201–209 72.Zhang H, Sanavia L, Schrefler B (1999) An interal length scale in dynamic strain localization of multiphase porous media. Mechanics of Cohesivefrictional Materials 4(5):443–460 73.Zienkiewicz OC, Chan A, Pastor M, Schrefler B, Shiomi T (1999) Computational geomechanics with special reference to earthquake engineering. Wiley Chichester, UK
0 Comments

Group NewsNews about Computational Poromechanics lab at Columbia University. Categories
All
Archives
December 2021
