iMechanica Journal Club Article Posted Online

Published online at http://imechanica.org/node/17096

Numerical modeling of thermo-hydro-mechanical coupling processes in porous media

Thermo-hydro-mechanics (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 multi-physical 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 pore-fluid 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 thermo-hydro-mechanical processes. This journal club article will focus on the derivation and implementation of the thermo-hydro-mechanical model for numerical simulations. Any feedback or critique are greatly appreciated.

1. Monolithic vs. operator-splitting schemes

The governing equation of the thermo-hydro-mechanical 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 pore-fluid diffusion but not vice versa. Li, Liu and Lewis, introduce a co-rotational FEM formulation and incorporate plasticity into THM model to model the non-isothermal 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 open-source project OpenGeoSys, which takes advantage of an object-orient framework and provides software engineering tools such as platform-independent 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 TOUGH-FLAC, which links flow simulator TOUGH2 with a small strain finite difference code FLAC (Rutqvist,2011). Klar et al. 2013 employ an explicit time-marching scheme to simulate the fully coupled thermal, and pore-fluid 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 operator-splitting 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 reactive-diffusion 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. Taylor-Hood, Raviart-Thomas)  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 inf-sup 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 non-physical results are due to the kernel (null space) of the discrete gradient operator being non-trivial. This non-trivial 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 inf-sup 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 Taylor-Hood finite element (quadratic basis functions for displacement and linear basis functions for pore pressure) to satisfy inf-sup 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 lowest-order Raviart-Thomas space, and piecewise constant pore pressure may also lead to stable solutions for isothermal poromechanics problems. Nevertheless, inf-sup 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 post-processing 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 equal-order 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 inf-sup stable finite elements (e.g. quadratic-displacement/linear-pore-pressure 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 over-diffusing 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 over-diffusion (Sun 2013a and Sun et al. 2013b). The stabilization term is derived by solving a small strain one-dimensional deformation-diffusion 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 three-field thermo-hydro-mechanical problem, the situation is more complicated. Since the thermal convection-diffusion and pore fluid diffusion may reach steady-state in profoundly different rates, it is unclear how to stabilize both the pore pressure and temperature fields without over-diffusing one or both of them.

3.   Length scale and bifurcations 

Unlike single-phase 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 one-dimensional 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 post-bifurcation 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 sub-scale 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 post-bifurcation regime even though the macroscopic responses are inferred from grain-scale simulations. 




Reference  

1.Abellan, M. A., & De Borst, R. (2006). Wave propagation and localisation in a softening two-phase medium. Computer Methods in Applied Mechanics and Engineering, 195(37), 5011-5019.

2.Armero F (1999) Formulation and finite element implementation of a multiplicative model of cou- pled poro-plasticity 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 inf-sup condition and its evaluation for mixed finite element methods. Com- 
puters and Structures 79(2):243 – 252, DOI 10.1016/S0045-7949(00)00123-1 

8.Bazant, Z. P. (2010). Can multiscale-multiphysics 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 low-order 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. Springer-Verlag 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) Assumed-deformation 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) Patch-averaged 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 Navier-Stokes 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), 789-818.

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/s11440-007-0033-0, URL 
http://dx.doi.org/10.1007/s11440-007-0033-0 

30.Jing L, Tsang CF, Stephansson O (1995) Decovalexan international co-operative 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, Regenauer-Lieb 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 open-source initiative for numerical simulation of thermo-hydro-mechanical/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 thermo-hydro-mechanical process in poro-elasto-plastic 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 convection-dominated 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 fluid-saturated 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 thermal-mechanical-hydrological 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 multi-scale simulations. Acta Geophysica, 62(5), 1109-1126.

40.Notz PK, Pawlowski RP, Sutherland JC (2012) Graph-based 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: Template-based 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: Template-based 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 multi-phase thermo-poromechanical 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 fluid-saturated elastic porous media 
with compressible constituents. Rev Geophys 14(2):227–241, DOI 10.1029/RG014i002p00227 

52.Rutqvist J (2011) Status of the tough-flac 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 component-based partial 
differential equation code built on trilinos. ACM Transactions on Mathematical Software 

54.Schrefler B (1995) Fe in environmental engineering: coupled thermo-hydro-mechanical 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- 
hydro-mechanical problems. Computer methods in applied mechanics and engineering 144(1):93–109 

56.Selvadurai A, Nguyen T (1997) Scoping analyses of the coupled thermal-hydrological-mechanical 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 poro-elastic and poro-elasto-plastic 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 elasto-plasticity. 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 thermo-hydro-mechanical 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 shear-enhanced compaction bands using multi- scale computations. Geophysical Research Letters 38(10):L10,302 

64.Sun W, Chen Q, Ostien J (2013a) Modeling the hydro-mechanical responses of strip and circular punch loadings on water-saturated collapsible geomaterials. Acta Geotechnica pp 1–32, DOI 10.1007/s11440- 013-0276-x 

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):65-75.  

66.Terzaghi Kv, Rendulic L (1934) Die wirksame flachenporositat des betons. Z O ́ st Ing-u 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 two-phase 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 low-order finite elements for coupled solid-deformation/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 Cohesive-frictional 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