Coupled Thermo-Hydro-Mechanical Simulation of Geological Porous Media A Subtitle that can be Very Much Longer if Necessary Master’s thesis in Applied Mechanics FANNY BYSTRÖM Department of Applied Mechanics CHALMERS UNIVERSITY OF TECHNOLOGY Gothenburg, Sweden 2017 Master’s thesis 2017:48 Coupled Thermo-Hydro-Mechanical Simulation of Geological Porous Media Fanny Byström Department of Applied Mechanics Division of Fluid Dynamics Chalmers University of Technology Gothenburg, Sweden 2017 Coupled Thermo-Hydro-Mechanical Simulation of Geological Porous Media FANNY BYSTRÖM © FANNY BYSTRÖM, 2017. Supervisor: Christoffer Källerman, ÅF Examiner: Lars Davidson, Applied Mechanics division of Fluid Dynamics Master’s Thesis 2017:48 ISSN 1652-8557 Department of Applied Mechanics Division of Fluid Dynamics Chalmers University of Technology SE-412 96 Gothenburg Telephone +46 (0)31 772 1000 Cover: Temperature field extracted from STAR-CCM+ Typeset in LATEX Gothenburg, Sweden 2017 iv Coupled Thermo-Hydro-Mechanical Simulation of Geological Porous Media FANNY BYSTRÖM Department of Applied Mechanics Chalmers University of Technology Abstract Countries around the globe have battled with the issue of nuclear waste for decades and the solution appears to be underground repositories. To assure its safety, both experimental and numerical testing is essential. In this master’s thesis, a numerical model has been developed in the software STAR-CCM+ with the intention of cap- turing the difficult thermo-hydro-mechanical (THM) coupling that occurs in porous rocks as a result of thermal loads. Results show that the respective mechanisms have been captured differently well. The thermal field is highly similar to that in the reference experiment while the pressure and displacement deviate to a higher extent. It was experienced difficult to identify "the villain of the piece", but the problems could possibly have originated from the numerical set-up not being ideal (highly sensitive simulation), assumptions and simplifications or weaknesses in the software that has a relatively newly intro- duced stress model. Furthermore, it is worth mentioning that, though the pressure and displacement are higher than expected, they still follow the correct behaviour throughout the simulation. Additionally, a sensitivity study was made that revealed the case’s strong dependence on the heat source and permeability of the rock. Concluding remarks include that THM coupling is highly complex and numerically sensitive. While a model was achieved and deemed relatively satisfactory, it is not yet completely reliable and accurate. To possibly accomplish this, further work is required in the form of general troubleshooting and a reduction of simplifications including for example the use of isotropy and constant values for parameters such as porosity and permeability. It is also advised to consider validating/trouble-shooting the model by using a different software, especially for the stress model. Keywords: Thermo-Hydro-Mechanical Coupling, Porous medium, STAR-CCM+, CFD, FEM, Nuclear waste, Nuclear repository v Acknowledgements The thesis work presented in this report was carried out between January 2017 and June 2017 at the consultant company ÅF in Gothenburg, Sweden. To them I owe a big thank you for taking me in and providing me this opportunity. A special thanks to my supervisor Christoffer Källerman for his support and expertise. I would also like to express my appreciation to Christian Windisch at the STAR-CCM+ support. Last, but not least, I would like to thank Lars Davidson at Chalmers University of Technology for accepting this work and acting as the examiner. Fanny Byström, Gothenburg, June 2017 vii Contents List of Figures xiii List of Tables xv 1 Introduction 1 1.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.3 Delimitations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 2 Theory 3 2.1 Nuclear Repositories and Related Studies . . . . . . . . . . . . . . . . 3 2.1.1 Repositories . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 2.1.2 Related Studies . . . . . . . . . . . . . . . . . . . . . . . . . . 3 2.2 Continuum Mechanics . . . . . . . . . . . . . . . . . . . . . . . . . . 4 2.3 Porous Media . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.3.1 Hydro-Mechanical Coupling . . . . . . . . . . . . . . . . . . . 6 2.4 Thermo-Hydro-Mechanical Coupling . . . . . . . . . . . . . . . . . . 7 2.4.1 Governing Equations . . . . . . . . . . . . . . . . . . . . . . . 8 2.4.1.1 Fluid Process . . . . . . . . . . . . . . . . . . . . . . 8 2.4.1.2 Deformation process . . . . . . . . . . . . . . . . . . 8 2.4.1.3 Heat transport . . . . . . . . . . . . . . . . . . . . . 8 2.4.2 Discretisation Method . . . . . . . . . . . . . . . . . . . . . . 9 2.4.3 Thermo-Hydro-Mechanical Coupling in STAR-CCM+ . . . . . 10 3 Methods and Numerical Set-Up 11 3.1 Computational Domain . . . . . . . . . . . . . . . . . . . . . . . . . . 12 3.2 Meshing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 3.3 Numerical Models and Solver Set-Up . . . . . . . . . . . . . . . . . . 13 3.3.1 User Defined Equation of State . . . . . . . . . . . . . . . . . 14 3.3.2 Properties of the Porous Medium . . . . . . . . . . . . . . . . 14 3.4 Boundary Conditions and Initial Conditions . . . . . . . . . . . . . . 15 3.4.1 Porous Domain . . . . . . . . . . . . . . . . . . . . . . . . . . 15 3.4.2 Solid Domain . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 3.5 Coupling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 3.6 Simulation Steps . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 3.7 Post-Processing and Validation . . . . . . . . . . . . . . . . . . . . . 17 ix Contents 4 Results and Discussion 19 4.1 Thermal Field (THM) . . . . . . . . . . . . . . . . . . . . . . . . . . 19 4.1.1 Comparison to HE-D in-situ test and Previous Simulations . . 20 4.2 Pore Pressure (THM) . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 4.2.1 Comparison to HE-D in-situ test and Previous Simulations . . 23 4.3 Displacement (THM) . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 4.3.1 Comparison to HE-D in-situ test and Previous Simulations . . 26 4.4 Effect and Strength of the Coupling . . . . . . . . . . . . . . . . . . . 27 4.5 Sensitivity Study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 4.6 Sources of Error . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 5 Conclusion 31 5.1 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 Bibliography 33 x Nomenclature Latin bw Water compressibility [1/Pa] b Body load [N/m3] c Specific heat [J/(kgK)] E Young’s modulus [Pa] I Identity tensor [-] k Thermal conductivity [W/(mK)] p Pressure [Pa] P Resistance tensor [kg/(m3s)] q Darcy flux/velocity [m/s] t Time [s] T Temperature [◦C] u Displacement [m] v Fluid Velocity [m/s] V Volume [m3] Greek β Thermal expansion coefficient [1/K] ε Strain [-] κ Permeability [m2] µ Dynamic Viscosity [Pa/s] ν Poisson’s ratio [-] ρ Density [kg/m3] σ Stress [Pa] τ Viscous stress [Pa] φ Porosity [-] Acronyms CFD Computational Fluid Dynamics CSM Computational Solid Mechanics EOS Equation of State FEA Finite Element Analysis FEM Finite Element Method FVM Finite Volume Method GUI Graphical User Interface HM Hydro-Mechanical THM Thermo-Hydro-Mechanical xi Contents xii List of Figures 2.1 THM Coupling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 3.1 Numerical Methodology . . . . . . . . . . . . . . . . . . . . . . . . . 11 3.2 Computational Domain . . . . . . . . . . . . . . . . . . . . . . . . . . 12 3.3 Finite Volume Mesh. Left: full domain, right: zoom close to heater . 13 3.4 Finite Element Mesh. Left: full domain, right: zoom close to heater . 13 3.5 Data Mapper Formulation . . . . . . . . . . . . . . . . . . . . . . . . 16 4.1 Temperature as function of time (left) and temperature as function of heater distance (right) . . . . . . . . . . . . . . . . . . . . . . . . 20 4.2 Temperature field at t= 90, 338 and 518 days, respectively . . . . . . 20 4.3 Temperature field at t = 260 days for reference simulation (left) [21] and here simulated fields (right) . . . . . . . . . . . . . . . . . . . . . 21 4.4 Temperature comparison between this study (red) and HE-D in-situ testing (symbol points) as well as previous simulation (black) [20] . . 21 4.5 Pore Pressure development over time . . . . . . . . . . . . . . . . . . 22 4.6 Pressure field at t = 90, 338 and 518 days, respectively . . . . . . . . 22 4.7 Velocity field at t = 93 days . . . . . . . . . . . . . . . . . . . . . . . 23 4.8 Pressure field at t = 260 days for reference simulation (left) [21] and here simulated field (right) . . . . . . . . . . . . . . . . . . . . . . . . 24 4.9 Pore pressure comparison between this study (red) and HE-D in-situ testing (symbol points) as well as previous simulation (black) [20] . . 24 4.10 Displacement (z-direction) as function of time (left) and distance to heater (right) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 4.11 Strain(z,z) at t = 93 days . . . . . . . . . . . . . . . . . . . . . . . . 25 4.12 Displacement (z-direction) at t = 90, 338 and 518 days, respectively . 25 4.13 Displacement (z direction) from Reference 2 (left) and here simulated in STAR-CCM+ (right) . . . . . . . . . . . . . . . . . . . . . . . . . 26 4.14 Displacement (z-direction) as function of distance from Reference 2 (left) and here simulated in STAR-CCM+ (right) . . . . . . . . . . . 26 4.15 Pressure and displacement as functions of time for TH/TM coupling compared to THM coupling . . . . . . . . . . . . . . . . . . . . . . . 27 4.16 Temperature and pore pressure as function of time for different powers 28 4.17 Temperature and pore pressure as function of time for different values of κ (TH coupling) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 xiii List of Figures xiv List of Tables 3.1 Meshing Set-Up . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 3.2 Solver Set-Up . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 3.3 Porous Properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 3.4 Heat Source Variation . . . . . . . . . . . . . . . . . . . . . . . . . . 16 xv List of Tables xvi 1 Introduction 1.1 Background Countries with nuclear power have been struggling with the issue of nuclear waste for decades. The combination of high-level and long-lived waste could, potentially, harm humans for tens of thousands of years, why the disposal has been heavily debated and evoked ethical questions. Today, there is what can only be described as an international consensus on the advisability of storing nuclear waste in deep geological repositories. While this may solve the main concern of many – avoiding radiation exposure to living creatures – questions regarding its geological impact arises. The radioactive decay of the inserted waste will release heat, applying a thermal load to the porous rock and affecting it in both a structural and fluid flow manner. Within a porous medium, there also exists a relation between pore pressure and mechanical stress, making this a thermo-hydro-mechanically (THM) coupled process. In order to assure the safety, the coupled effects must be predicted and even though THM models have been available since the 1980s, it is still considered a major challenge. The software STAR-CCM+ has, much due to its favourable GUI and powerful post-processing, become highly popular. In recent years FEA was introduced into the software, meaning that an all-in-one solution to the thermo- hydro-mechanical coupling described above could/should be possible. Investigating the possibility to develop such a desired numerical model within the software was the underlying reason and motivation for this study. 1.2 Objectives The general goal of this master’s thesis was to develop a numerical model that can predict the thermo-hydro-mechanically coupled processes that occur in geological porous media surrounding a heat source such as nuclear waste. Particularly, the goal was to develop such a model in the software STAR-CCM+. The study has considered the following questions: • How can coupled THM behaviour in geological porous media be modelled and simulated in STAR-CCM+? – How can the temperature development within the medium numerically be captured? – How can the pore pressure and flow of groundwater within the medium numerically be captured? 1 1. Introduction – How can the deformation and stresses within the medium numerically be captured? – Which mechanisms and parameters are dominant in the coupling process? 1.3 Delimitations The boundaries applied to the project are: • The time of the project is limited to 20 weeks (full time). • Simulations are to be performed in STAR-CCM+. • Solely one type of geological media is considered. • Chemical coupling is not considered. • Heat is given strictly from the heat source. The main assumptions applied to the project are: • The porous medium is homogeneous and isotropic. • The solid skeleton of the porous medium is linear-elastic. • The porous medium is fully saturated. • The flow is laminar. • Local thermal equilibrium throughout the domain. 2 2 Theory In the following chapter, a brief introduction and background to nuclear repositories will be presented together with relevant theory to understand the coupled processes and how they can be captured numerically. 2.1 Nuclear Repositories and Related Studies 2.1.1 Repositories Deep geological repositories are official policy for nuclear disposal in many countries with high-level and long-lived waste. A country is, however, yet to dispose their fuel, but many have come a long way in the planning and research stage with for example sites and strategy decided [1]. The typical repository is a multibarrier system that would be located >300 meters below ground in mined repositories or in boreholes. The purpose is to provide a long- term and safe isolation as well as minimising the requirement for future maintenance. The multibarrier system consists of both a natural barrier in the form of the host rock as well as an engineering barrier system including buffer material and metal waste canisters. Typical host rocks that are considered are graphite and clay due to their low permeability, low diffusion coefficients and a degree of self-healing when fractured [2]. 2.1.2 Related Studies Over the past thirty years both experimental and numerical testing have been per- formed in order to assure the safety of the rock and constrain the transport of radionuclides by obtaining an understanding of the difficult THM coupling [3]. In situ Experiments The presumable use of geological nuclear disposals created the need for underground laboratories. Many such were constructed around the world, such as Äspö Hard Rock Laboratory in Sweden, FEBEX (Full-scale Engineered Barriers Experiment) and the Mont Terri rock laboratory in Switzerland. Experiments related to THM coupling are generally performed by inserting heaters into rocks, with or without engineering barriers, at depths of several hundred meters and sensors are used to monitor parameters such as temperature, pressure and strains. Furthermore, exper- iments are also used to determine material properties of the barriers. For further 3 2. Theory information, see for example the Mont Terri Project’s website [4], where information about the rock laboratory and several experiments are presented. Numerical There exist numerous codes for THM-modelling of porous media. A common ap- proach is to use open source codes based on the finite element method, such as OpenGeoSys [5], ROCMAS [6], THAMES [7], FRACON [8] or FRT-THM [9] to mention but a few. Other existing possibilities are, for example, the commercial finite-difference code FLAC [10] or commercial all-in-one software such as COM- SOL Multiphysics, ABAQUS and ANSYS. Many of the mentioned modelling meth- ods were created in the international project DECOVALEX, short for DEvelopment of COupled models and their VALidation against EXperiments. This is an ongo- ing project, founded in 1992 and consist of nuclear waste management organisations and regulatory authorities from 13 countries, who work on developing and validating models that can predict the behaviour of rocks subjected to thermal loads. Ulti- mately, no information about previous THM coupled simulations in STAR-CCM+ has been found. 2.2 Continuum Mechanics In order to understand the physics of a porous medium, it is essential to under- stand the physics of its constituents - solids and fluids. Their mechanics, i.e. their behaviour when subjected to forces/displacements and their effect on their sur- roundings, is often modelled using the continuum assumption. This means that the assumption of continuous mass is made, or in other words that the body is entirely filled with material matter. This macroscopic view is accurate on length scales much greater than that of interatomic distances. [11] Continuum mechanics is generally modelled using the balance of mass, momentum and energy together with constitutive relations. Below, the general form of the three equations is given for any continuum [12]. ∂ρ ∂t +∇ · (ρv) = Qm (2.1) ∂ρv ∂t +∇ · (ρvv) = ∇ · σ + b (2.2) ∂ρcpT ∂t +∇ · qT = QT (2.3) In the above, v is the fluid/solid velocity, Qm the source/sink term representing added/removed mass, σ the stress tensor, b the body force per unit volume, cp the specific heat (under constant pressure), qT the heat flux and QT the heat sink/- source term. Constitutive relations will vary whether the body is a fluid or solid, giving solids and fluids different versions of the above general equations. Just like the above equation, the fluid and solid specific equations can take many forms. Be- low follow some general constitutive relations and governing equations specified for 4 2. Theory fluids and solids, respectively. Fluid Mechanics The most common approach to describe the behaviour of fluids is to follow the equa- tions developed by Navier and Stokes, which are written similarly as the equations above. Normally, the stress tensor is re-written in terms of pressure and viscous stress (τ ) which are considered the two contributors to forces acting on the fluid. This division yields the following ∂ρv ∂t +∇ · (ρvv) = −∇p+∇ · τ + b (2.4) Assuming a Newtonian fluid, the viscous stress is linearly proportional to the velocity gradient and the above becomes [12] ∂ρv ∂t +∇ · (ρvv) = −∇p+∇ · ( µ(∇v +∇vT ) ) +∇ (−2µ 3 ∇ · v ) + b (2.5) where µ is the dynamic viscosity of the fluid. Solid Mechanics In solid mechanics, it is common to express the balance equations using a Lagrangian approach. This means that the convective term in equation 2.2 vanishes and the time derivative of the velocity reduces to the partial second derivative of the displacement, resulting in the so called equation of motion, or equilibrium equation, shown below [13] ρü = ∇ · σ + b (2.6) where u is the displacement vector and ü is hence the acceleration of the solid. Under the assumption of small deformations the above equation becomes ∇ · σ + b = 0 (2.7) To solve the equation one needs to establish constitutive relations. This can be done by utilising Hooke’s law and thereby assuming a linear-elastic material, meaning that the stress-strain relation is linear and that the material deforms reversibly under stress. Hooke’s law is defined as σ = Cε (2.8) where ε is the strain rate tensor defined as ε = 1 2[∇u+ (∇u)T ] (2.9) and C is the fourth-order stiffness (or elasticity) tensor which, for an isotropic ma- terial, is defined as: C = E (1 + ν)(1− 2ν)  1− ν ν ν 0 0 0 ν 1− ν ν 0 0 0 ν ν 1− ν 0 0 0 0 0 0 1−2ν 2 0 0 0 0 0 0 1−2ν 2 0 0 0 0 0 0 1−2ν 2  (2.10) 5 2. Theory where E is the Young’s modulus and ν Poisson’s ratio. The assumption of isotropy means that the properties of the material is independent of direction. 2.3 Porous Media A porous medium is a material with a solid skeleton containing fluid filled pores. Porous media is typically associated with natural objects, such as soils and rocks. Important parameters that characterise a porous material include porosity, perme- ability and saturation. Porosity is the composition of fluid in contrast to the volume of the medium. It is defined as φ = Vv Vt , (2.11) where Vv is the volume of the void and Vt is the total volume of the medium. The porosity is for example used to define material properties by proportionally blending values from the solid and fluid as well as in the time-derivative term in the scalar transport equation for unsteady flows. Permeability, often denoted κ, measures a material’s ability to transmit fluids, where a high permeability means that the fluid/s can pass rapidly through the medium. The intrinsic permeability depends on porosity and when choosing a host rock for a repository, a rock with low permeability is desirable, since it would limit the move- ment of fluid/s and consequently the radionuclide. The saturation, often denoted S, describes the water, or moisture, content of the material. It ranges from 0, meaning fully dry pores, to 1, meaning that all pores are filled with water. If fully saturated conditions are applicable/assumed, the material consists solely of the solid skeleton and liquid. The density then becomes ρ = φρw + (1− φ)ρs, (2.12) where the subscripts w and s represent water and solid, respectively. 2.3.1 Hydro-Mechanical Coupling Porous media is generally described using poromechanics, or the theory of poroelas- ticity if the assumption of an elastic solid and viscous fluid is applicable. It describes the behaviour of, and dependence between, fluid motion and solid deformation and is central to the phenomena of groundwater flow. In the beginning of 1920s Karl von Terzaghi was the first to express parts of this behaviour in the Terzaghi’s principle, or Terzaghi’s theory of one-dimensional consolidation. The theory states that a soil subjected to a stress will be opposed by the pore pressure. This is the definition of the so called effective stress, which is expressed as [14] σ′ = σ − p (2.13) 6 2. Theory where σ′ is the effective stress, σ is the total solid stress and p is the pore pressure. Since the development of Terzaghi’s principle, many have expanded on the descrip- tion of porous media. A commonly used theory and model for isothermal poroelas- ticity was developed by M.A Biot and uses the coupling of Darcy’s law and the law describing linear elasticity together with the Terzaghi’s principle [15]. Darcy’s law describes the discharge rate, also called Darcy flux or Darcy velocity, through a saturated porous media and is defined as q = −κ µ ∇p (2.14) or in terms of the fluid velocity: v = − κ φµ ∇p (2.15) Darcy’s law is valid for Stokes flows with negligible inertia and thus laminar flow prevails. For flows with higher velocities, or for unsaturated media, the equation requires modifications. Following Biot’s theory, the general stress-strain relationship of solid materials (Sec- tion 2.2), is incorporated with the principle of effective stress, meaning that the stress can be expressed by the strain of the solid skeleton and the pore pressure originating from the fluid. This leads to the following expression: σ′ = σ − pI = E 1 + ν ( ε+ ν 1− 2ν εvI ) − pI (2.16) where εv is the volumetric strain. In the above, the material is once again assumed to be linear-elastic as well as isotropic, meaning that the pore pressure does not cause shearing and it affects the normal strains equally, independent of direction. 2.4 Thermo-Hydro-Mechanical Coupling When a porous medium is experiencing a temperature gradient, an additional cou- pling to the above problem is added, namely the thermal coupling. When the porous medium is heated, thermal expansion will occur and thermal stresses will be induced. The expansion will also result in a rise of pore pressure and this hydraulic gradient can trigger the water to flow away from the heat source. If the porous medium is sufficiently permeable, water is pumped out and consolidation might occur which would cause the pore pressure to dissipate. Furthermore, a change in pore pres- sure will alter the effective stress causing movement of the solid. Ultimately, the expansion and general movement of the solid changes the volume relations and consequently the porosity and permeability. In Figure 2.1 a schematic of the full thermo-hydro-mechanical coupling is illustrated. 7 2. Theory Figure 2.1: THM Coupling 2.4.1 Governing Equations In the following section the governing equations of THM coupling are presented. They describe the non-isothermal behaviour of a saturated porous medium by ex- pressing the main physics of fluid flow, solid deformation and heat transfer. 2.4.1.1 Fluid Process The flow through a porous medium is governed by Darcy’s law, given in equation 2.14. Inserting this into the mass conservation (2.1) and accounting for the impact from the solid movement yields ∂(φρw) ∂t +∇ · (ρwqw) +∇ · (ρwu̇) = Qw (2.17) 2.4.1.2 Deformation process By introducing Biot’s theory into the equation of motion (balance equation of linear momentum), see equation 2.6, and adding the impact from the thermal stresses yields the following [17] ∇ · (σ − pI − αE∆TI) + ρb = 0 (2.18) where α is the thermal expansion coefficient and ∆T is the temperature difference. 2.4.1.3 Heat transport Heat transport occurs through conduction, convection and/or radiation. In a porous medium, conductive and advective transport is normally considered, but in porous geological media, however, the rocks generally have such low permeability that the convective fluxes can be neglected as well. Hence, conduction is the dominant trans- port mechanism and the flux in equation 2.3 is described by Fourier’s law [16] qT = −∇ · (k∇T ) (2.19) 8 2. Theory resulting in the following transport equation ρc ∂T ∂t −∇ · (k∇T ) = QT (2.20) where c is the specific heat, k the thermal conductivity and QT the heat source. The heat capacity of the porous medium, ρc, is defined as: ρc = φcwρw + (1− φ)csρs (2.21) and the thermal conductivity in a similar manner as k = φkw + (1− φ)ks (2.22) where w and s once again represent water and solid, respectively. 2.4.2 Discretisation Method The governing equations presented above rarely have analytic solutions, but need to be transformed into a set of algebraic equations by discretising them in both space and time. Spatial Discretisation By utilising a spatial discretisation method, one can find approximate solutions to the continuous partial differential equations by transforming the equations into dis- crete algebraic equations. The most common methods include the finite difference method (FDM), finite volume method (FVM) and finite element method (FEM). Commonly, FVM is used in fluid dynamics/CFD while FEM is used in solid me- chanics/CSM. In the finite volume method the domain is divided into a finite number of control vol- umes/cells, where the variable of interest is located at the centroid. The governing equations are thereafter integrated over the control volume. Volume integrals that contain a divergence term are converted to surface integrals, using the divergence theorem. These are then evaluated as fluxes at the surfaces of each control volume, where the flux entering a cell must be identical to that leaving the adjacent. The discrete solution of a control volume is thereafter interpolated to surrounding cells, generally by assuming a piece-wise linear variation, in an iterative process, assuring that the conservation is satisfied for each and every cell. [12] In the finite element method, the domain is similarly divided into a finite number of elements. The governing equations are integrated over each element after having been multiplied with a weight function. The integrals are evaluated numerically using integration rules, such as Gauss integration rule. The distribution of the dependent variables is interpolated by the use of so called shape functions (same form as weight function). Additionally, the boundary conditions, which in FVM are added after, are in FEM part of the discretised equations. For more information, see for example [18] 9 2. Theory Temporal Discretisation In transient simulations the additional discretisation in time is necessary. This is achieved by integrating the discrete equations, term by term, over a chosen time- step. It can be executed using a variety of techniques, where the most common ones include the explicit and implicit method. The integral form of their general equation can be written as: ϕn+1 − ϕn ∆t = f · F (ϕn+1) + (1− f) · F (ϕn) (2.23) where ϕ is a scalar quantity and f is a weight between 0 and 1. f=0 → Explicit method, meaning that F(ϕ) is evaluated at current time-step. f=1 → Implicit method, meaning that F(ϕ) is evaluated at future time-step. 2.4.3 Thermo-Hydro-Mechanical Coupling in STAR-CCM+ The assigned software for this thesis work is, as mentioned, STAR-CCM+. STAR- CCM+ was, originally, a pure CFD software, but in 2015 a solid stress model was added [19]. The model, however, is not compatible with a porous medium, in- stead the stress analysis requires an entirely solid region. The mechanical coupling (TM/HM) within the medium is therefore not integrated, but needs to be added manually by the user in order to capture the full coupling. A porous medium is modelled by the addition of a momentum source term in the general momentum equation (2.2). The source term is defined as fp = −P · vs (2.24) where P is the resistance tensor and vs is the superficial velocity, which is the same as the Darcy velocity in fully saturated flows. The resistance tensor is defined as P = Pv + Pi|vs| (2.25) where Pv is the viscous resistance tensor and Pi is the inertia resistance tensor. The viscous resistance relates to the viscosity of the fluid and the structure of the skeleton while the inertia resistance takes into account losses due to expansion, constriction and bending of the pores. In low-velocity simulations, the inertia resistance can be neglected and the flow can be considered laminar. The model then reduces to Darcy’s Law, with the viscous resistance defined as Pv = µ κ (2.26) 10 3 Methods and Numerical Set-Up In the following chapter the methodology and numerical set-up is explained. It ex- plains what was done prior to, during and after simulations as well as the principal steps involved. Since the aim was to develop a simulation model that can capture the difficult thermo-hydro-mechanical coupling, rather than predicting specific values that would exist in a real nuclear repository, strive was to avoid unnecessary complexity. The desire for simplicity combined with the need for a reference to use as guideline and validation, motivated the choice to resemble an in-situ test at Mont Terri rock labo- ratory called HE-D. The experiment was one without any engineering barrier, i.e. it only consisted of heaters inserted in the rock Opalinus clay. Simulations performed on this experiment will also be used as references and numerical guidelines. For more information about the reference experiment/simulations, the reader is referred to [20] and [21]. In Figure 3.1 below, the general outline of the methodology is illustrated. Figure 3.1: Numerical Methodology 11 3. Methods and Numerical Set-Up 3.1 Computational Domain The geometry was constructed using the 3D-CAD modeller in STAR-CCM+ and is visualised below in Figure 3.2. As can be seen, the domain consists of a plane cube with a horizontal cylinder placed in the centre to represent the waste canister/heat source. The cube has the dimensions of 16×20×16 m and the cylinder has a radius of 0.15 m and is 5 m long. Since no interest exists in the heater itself, it was cut out of the domain and hence acts as a boundary to the model. Figure 3.2: Computational Domain Due to the assumed isotropy, the box shown above was cut in two (at y=10m) and symmetry was applied, leading to fewer cells and consequently shorter simulation time. The created geometry part was duplicated in order to create two different regions, continuums and meshes for the flow and stress solvers, respectively. 3.2 Meshing The spatial discretisation was made by using a combination of the finite volume and finite element method. FVM was used for the mesh that solves the flow and energy equations while FEM was used for the mesh that solves the solid stress equations. Both meshes were created using the built-in automated mesh-creator. In Table 3.1 below, the meshing settings are illustrated. As can be seen, the core volume mesher is polyhedral and tetrahedral for the FV and FE discretisation, respectively. In general, polyhedral elements are to prefer over tetrahedral, since the mesh would require less memory and less computing time to reach the same accuracy [22]. The automated mesher in STAR-CCM+ can, however, only create tetrahedral elements when using the solid stress-model and the finite elements it requires. 12 3. Methods and Numerical Set-Up Table 3.1: Meshing Set-Up Finite Volume Finite Element Element type Polyhedral Tetrahedral Base size [m] 0.5 0.5 Minimum Surface Size [%] 10 10 Number of cells 30831 164708 Surface Remesher Automatic Surface Repair Mid-side Vertex Option - The settings expressed above resulted in the final meshes shown below in Figure 3.3 and Figure 3.4. In the figures to the right, a close up on the mesh adjacent to the heater is shown. This is the area where the gradients are expected to be the highest, why the meshes were refined in this zone (10% of base size, as described in the table). Figure 3.3: Finite Volume Mesh. Left: full domain, right: zoom close to heater Figure 3.4: Finite Element Mesh. Left: full domain, right: zoom close to heater 3.3 Numerical Models and Solver Set-Up In Table 3.2 below, the models chosen for the fluid and solid continuums are dis- played together with relevant specifications. 13 3. Methods and Numerical Set-Up Table 3.2: Solver Set-Up Model/Solver Specification Both Regions Implicit Unsteady Time-Step: 1 h Temporal Discretisation: 1st order Three Dimensional Cell-Quality Remediation Porous Region Laminar Liquid Water Gradients Hybrid Gauss-LSQ Segregated Flow ↪→ velocity Under-relaxation factor: 0.1 ↪→ pressure Under-relaxation factor: 0.3 Segregated Energy Under-relaxation factor fluid: 0.9 Under-relaxation factor solid: 0.99 User Defined EOS Compressible Solid Region Linear Isotropic Elastic Solid Solid Stress Interpolation Method: Backward Euler Specified Temperature Load 3.3.1 User Defined Equation of State Water is commonly assumed incompressible due to its high bulk modulus. How- ever, the low permeability makes the clay close to impermeable which increases the pressure remarkably and together with the increase of temperature, the den- sity variation should no longer be ignored. Considering the compressibility of the groundwater means that a user-defined field function had to be created in order to determine the dependence of temperature and pressure on density, ρ = f(p, T ). This can be retrieved from steam tables, or, as an alternative, simplified empirical relationships can be used. The following relation is frequently used and was chosen for this task [23] ρw ρw0 = 1 + βw(p− p0)− bw(T − T0) (3.1) where βw and bw are the water compressibility and expansion coefficients, respec- tively and the subscript 0 represents a reference state. The two coefficients were given constant values of 4× 10−10 1/Pa and 0.00021 1/K. Furthermore a temperature dependent dynamic viscosity was defined as [24] µ(T ) = 2.414× 10−5 × 10247.8/(T−140) (3.2) 3.3.2 Properties of the Porous Medium Properties of the rock were chosen to replicate previous studies and experiments performed on Opalinus clay. In Table 3.3 below, the properties of the clay are 14 3. Methods and Numerical Set-Up shown [20][21][25]. Table 3.3: Porous Properties Parameter Value Petrophysical properties Density, ρ 2450 kg/m3 Porosity, φ 0.16 Hydraulic properties Intrinsic Permeability, κ 5 × 10−20 m2 Mechanical properties Young’s Modulus, E 7.5 GPa Poisson’s Ratio, ν 0.285 Thermal and thermo- Thermal Conductivity, k 2.5 W/(m·K) mechanical properties Heat Capacity, c 800 J/(kg·K) Coefficient of thermal expansion, β 1.4 × 10−5 1/K A porous region in STAR-CCM+ requires a fluid continuum, in this case liquid. Due to this, a number of properties of the solid phase needed to be specified for the solver to obtain the correct porous properties with the help of the specified porosity. For the liquid continuum, water was selected and the only values altered with were density and dynamic viscosity due to the assumed compressibility, see Section 3.3.1. The solid density, solid heat capacity and solid thermal conductivity were obtained from the relations given in equations 2.12, 2.21 and 2.22. Since little variation occurs in the solid density is was assumed constant. Material properties in the solid continuum were in a similar manner, with the help of porosity, altered with in order to account for the presence of water. 3.4 Boundary Conditions and Initial Conditions The two different regions and continuums belonging to the porous and solid domains required their respective boundary conditions and initial conditions, which will now be presented. 3.4.1 Porous Domain The domain was given initial conditions according to field values, namely an initial temperature of 15◦C, initial water pressure of 0.9 MPa and zero velocities. All outer boundaries of the box (except for the symmetry boundary) were set to pressure outlets in order to allow for water to be squeezed out of the porous domain with the Dirichlet boundary conditions T=15◦C and p=0.9 MPa. Furthermore, the inner cylinder boundary was set to a heat source that varied according to the HE-D in-situ experiment at the Mont Terri rock laboratory, as shown below in Table 3.4. When the heat source reached a new stage of heating/cooling, it was increased/decreased linearly during the first day of the new stage. 15 3. Methods and Numerical Set-Up Table 3.4: Heat Source Variation Time [Days] Heat Source [W] 0-90 650 91-338 1950 339-518 0 3.4.2 Solid Domain The solid continuum was given the initial condition of 15◦C, zero displacement and zero velocity. The outer boundaries of the box (except for symmetry boundary) were set to fixed in the normal direction but were allowed lateral displacement, while the cylinder boundary was set to fixed. 3.5 Coupling Since a porous medium in STAR-CCM+ is not compatible with a stress solver the mechanical processes are not accounting for the impact of pore pressure nor the thermal expansion caused by the heat source. In a similar manner, the impact of the solid displacement on the pore-water is not included, why these need to be added manually. This was achieved by using user defined field functions together with so called "Data Mappers", that transfer information from a source to a target mesh. In this study, volume mappers have been used with the formulation "Least Squares Interpolation". In Figure 3.5 below, an illustration of a set-up that transfer data from the centroid of one mesh (blue) to the centroid of a second mesh (red), using the a least square scheme is shown. Face 0 is the closest source to the target Face a and the neighbours of Face 0 are defined as any face that shares at least one vertex with Face 0. As can be seen this includes Faces 2-6 and 8, which together with 0 are included as part of the interpolation stencil for the value assigned to a. Figure 3.5: Data Mapper Formulation 16 3. Methods and Numerical Set-Up FE to FV The porous region and the solvers it uses, do not solve stresses/strains/displace- ments. Looking at equation 2.17, the last term on the left hand side is therefore not a part of the governing equations in STAR-CCM+, why it has to be added to achieve the full coupling. In order to obtain the time derivative of the displacement vector a field sum monitor was first created, where the two just completed (time-step based) displacements are saved. Obtained from the field sum monitor is therefore u(t-1)+u(t). Thereafter a field function was created in order to obtain the time derivative with the following expression: (2*u - field sum)/timestep. For the first iterations, the field sum monitor has not got any values, why a condition to the field function was added (given zero displacement here). When the time derivative was obtained, a final field function was created with the full expression of the term. The function was mapped from the solid region to the porous region where it was used as a source term in the mass balance equation, hence the sign of the field function becomes negative (sink). The mapping was performed at each timestep by interpo- lating data from the cells of the FE mesh to the cells of the FV mesh. FV to FE The solid stress model does not include temperature nor pressure, which are required as seen from equation 2.18. The temperature-term was added by simply "mapping" the temperature from the porous region onto the solid region, where the option to specify the temperature was made available by the model "Specified Temperature Load", see Table 3.2. Furthermore, the pressure-term was added by creating a field function with the definition ∇p and mapping this from the porous region to the solid region where it works as a body load in the equilibrium equation. Both mappings were performed each timestep by interpolating data from the cells of the FV mesh to the vertices of the FE mesh. 3.6 Simulation Steps The simulations were performed in advancing complexity by adding the coupling mechanisms step by step. The main steps were performed in the following order 1. Thermo-hydraulic coupling. 2. Thermo-mechanical coupling is added to the above. 3. Thermo-hydro-mechanical coupling. The thermal load is the driving mechanism and it is of great importance to capture it correctly, why this was made first. The above list describe the main steps, but naturally, a more detailed description could be done, involving for example steady simulations, simulations on a solid only, continuous increase of permeability, simu- lation to test pressure difference in a porous cube etc. 3.7 Post-Processing and Validation The main parameters to show was decided to be temperature, pressure and displace- ment. Thereby, all three couplings are presented and, at the same time, these are 17 3. Methods and Numerical Set-Up the variables that act as inputs into the other respective solver. Additionally, some discussion and illustration of the velocity and strain will be had and shown. Results will be presented in three ways - instantaneous contour plots, functions of time and functions of distance to the heat source. The contour plots are hardcopied directly from STAR-CCM+ at times according to [20] (before and after changes of the heat source), while the time- and space dependent values are exported and thereafter imported into MATLAB where they are plotted. By extracting similar/equal fig- ures and plots as previous simulations it becomes more straightforward to validate the results and draw conclusion about how well the coupling mechanisms have been captured. Two different sources will be used as references. The first one, Reference 1, is taken from [21] and is the one primarily used to resemble. However, the report does not show many results, why also results from [20] are used for comparison. This, further on referred to as Reference 2, simulates the same experiment and also shows results together with experimental values. Due to the assumed isotropy together with equal boundary conditions around the model, only one plane and one direction were considered necessary for illustrating the results. The contours plots will be shown in the x-z plane while z will be the direction to show for parameters such as displacement. Furthermore, 5 points were created along a line, starting from the heater-boundary and going all the way to the top surface, in order to illustrate the difference within the domain as a function of time. The points are at a distance of 0, 1, 3, 5 and 7 meters from the heater boundary. 18 4 Results and Discussion In the following chapter, the results are presented, analysed and compared - to each other as well as to in-situ experiments and previous simulations. Additionally, the impact of the different coupling mechanisms are investigated together with critical and sensitive parameters. 4.1 Thermal Field (THM) In the left figure of Figure 4.2, the temperature development during the 518 days of simulation, at five different points in the domain, is shown. The different stages of the heat source are easily identified and the reaction to changes in power is, as can be seen, instant and strong. The majority of the temperature increase/decrease, of each stage, is achieved during the first couple of days. In fact, 80 % of the increase of the first and second heating stage (at the heater boundary) has occurred after only approximately 9 and 13 days, respectively. Furthermore, the maximum temperature in the domain is 108.5◦C which is reached just before the heat source is turned off (338 days). It can also be seen that, due to the non-existing convective heat transfer, the spread in the domain is relatively slow and that, after the 180 days of cooling, the maximum temperature in the domain is 15.5◦C, i.e. it has almost returned to its original 15◦C. In the right figure, it is shown how temperature varies as you move away from the heater. It is plotted at a variety of times, representing crucial moments (90 days = end of heating stage 1, 93 days = beginning of heating stage 2, 338 days = end of heating stage 2, 346 days = beginning of cooling stage and 518 days: = end of cooling stage). Since the temperature at the edge of the domain is relatively constant, the temperature difference in the domain is largest when the heater temperature has its maximum just before the heater is turned off. It also becomes clear(er) that the difference within the domain is close to zero in the end of the simulation. 19 4. Results and Discussion Figure 4.1: Temperature as function of time (left) and temperature as function of heater distance (right) In Figure 4.2, the instantaneous temperature field in a cross section is shown at three different times, namely 90 days (end of heating stage 1), 338 days (end of heating stage 2) and 518 days (end of cooling stage/simulation). The three figures have been scaled equally in order to get a better view of the difference over time. Figure 4.2: Temperature field at t= 90, 338 and 518 days, respectively 4.1.1 Comparison to HE-D in-situ test and Previous Simu- lations Based on the HE-D experiments and the simulations performed in previously men- tioned reports, the thermal field has been captured very well. In fact, the instanta- neous temperature field in Reference 1 ([21]) has the exact same maximum temper- ature at the heater boundary (105◦C) at t = 260 days, which can be seen below in Figure 4.3. 20 4. Results and Discussion Figure 4.3: Temperature field at t = 260 days for reference simulation (left) [21] and here simulated fields (right) Furthermore, plots from Reference 2 ([20]) were used to compare the temperature variation over time at the heater. In Figure 4.4 below, the compared results are shown in black (where the straight line represent their simulated temperature and the symbol points are the measured temperatures at different sensors) and the here, in STAR-CCM+, simulated temperature is shown in red. This emphasises the state- ment that the temperature is well captured, however slightly faster in the reaction to changes. Figure 4.4: Temperature comparison between this study (red) and HE-D in-situ testing (symbol points) as well as previous simulation (black) [20] 4.2 Pore Pressure (THM) In the left figure of Figure 4.6, the pore pressure at the heater boundary is shown as a function of time together with temperature. Just as temperature reacts in- stantaneously to the heat source, the pressure similarly reacts immediately to the temperature gradient and the higher the temperature increase/decrease, the greater is the pressure jump/drop. However, pressure only follows the evolution of temper- ature to a certain extent, since it also can be seen that, while temperature has is maximum at the end of heating stage 2, pore pressure has its extreme value after only approximately 97 days, when it reaches 4.8 MPa. The reason for pressure not 21 4. Results and Discussion continuing to increase even though temperature does is that the pore pressure dis- sipation due to groundwater flow and consolidation overcomes the thermal effects. When the heat source is turned of and the temperature close to the heater experi- ences a large temperature drop, the pressure drop too becomes massive and in fact reaches negative values. In the right figure below, pressure is shown over time at the five selected points in the domain. It is clearly seen that as you move away from the heater, the less visible are the different stages of the heater. However, compared to the temperature plot displayed in the previous section, the pressure has a much stronger variation further out in the domain with high pressures not only close to the heater. This indicates that even small temperature gradients create a relatively large pressure increase and, of course, the fact that the pore pressure is also affected by the mechanical processes. Furthermore, unlike temperature, pressure is not fully recovered at the end of the simulation. Figure 4.5: Pore Pressure development over time In Figure 4.6 below, the instantaneous pressure field is shown at t=90, 338 and 518 days. The magnitudes are once again scaled, this time with the top value from 338 days and the bottom value from 518 days, due to the negative pressure experienced at this time. Figure 4.6: Pressure field at t = 90, 338 and 518 days, respectively Opalinus clay is almost impermeable with its intrinsic permeability of approximately 5 × 10−20 m2, meaning that even high pressure and temperature gradients should 22 4. Results and Discussion not create any significant ground water motion. As can be seen in Figure 4.7, which shows the velocity field after 93 days (when large gradients are experienced), this is also the case. The water is propagating away from the heater and out of the domain, but with a velocity that, if kept at this value which is to be considered as an extreme value, would require several hundred years to transport the water only 1 meter. The radionuclides that in reality might be released from the waste canisters would therefore not migrate fast in the clay rock. Figure 4.7: Velocity field at t = 93 days 4.2.1 Comparison to HE-D in-situ test and Previous Simu- lations Below in Figure 4.8, pressure contours are shown at t = 260 days. At this point, compared to Reference 1, the magnitudes are highly similar (3 MPa in left figure). The spread, however, seems to be further in the STAR-CCM+ simulation. In Figure 4.9, the pressure over time is shown for pore pressure simulate in this study (red) on top of a graph from Reference 2 (black). As can be seen, the pressure has been captured above that of the reference(s) where the maximum value experi- enced by the rock is approximately 3.6 MPa, hence slightly more than 1 MPa below the ones obtained in this study. It is, however, visible that the simulation follows the correct behaviour. The difference either indicates that the pressure reacts too strongly against temperature changes or the pore pressure dissipation has not been captured to its fullest. Furthermore, it can be seen that the experiment (symbols) also reaches negative values when the heater is shut off (approximately -0.5 MPa), why this is not a numerical issue, but in fact must mean that the "pressure" is pulling instead of pushing. By proper definition, of course, this is no longer a pressure. 23 4. Results and Discussion Figure 4.8: Pressure field at t = 260 days for reference simulation (left) [21] and here simulated field (right) Figure 4.9: Pore pressure comparison between this study (red) and HE-D in-situ testing (symbol points) as well as previous simulation (black) [20] 4.3 Displacement (THM) In Figure 4.10 the z-displacement is shown as a function of time and distance to heat source in the left and right figure, respectively. Once again, the three stages of the heat source are clearly distinguishable. The maximum displacement is slightly below 2 mm and with the size of the domain (16x20x16 m), the assumption of small deformations can therefore be assumed satisfied. The reason behind these small displacements lie in the stiff nature of the clay together with its modest thermal expansion coefficient. Noticeable is also that the displacement has its maximum at a certain distance away from the heater, meaning that the thermal expansion of the rock, which naturally is highest closest to the heater, pushes the clay out, which in turn creates compression in the outer zones of the domain as illustrated by the negative strain in Figure 4.11. 24 4. Results and Discussion Figure 4.10: Displacement (z-direction) as function of time (left) and distance to heater (right) Figure 4.11: Strain(z,z) at t = 93 days Below, the instantaneous displacement is shown at three different times. The figures have been scaled equally in order to clarify the difference. All though barely visible due to the small difference within the domain, the displacement has, like pressure, inverted in the end of the simulation. Figure 4.12: Displacement (z-direction) at t = 90, 338 and 518 days, respectively 25 4. Results and Discussion 4.3.1 Comparison to HE-D in-situ test and Previous Simu- lations From the contour plots in 4.13, which show the z-displacement at t=260 days, it can be observed that the spread in the domain follows a similar behaviour to the simulation in Reference 1, all though with larger magnitudes (just above 1 mm compared to just below 2 mm). Figure 4.13: Displacement (z direction) from Reference 2 (left) and here simulated in STAR-CCM+ (right) In Figure 4.14, a more thorough comparison is possible, where displacement as func- tion of distance to the heater is shown at different times for Reference 2 (left) and STAR-CCM+ simulation (right). They do experience somewhat similar displace- ments at the different times and locations, but it can be noted that the displace- ments are too high in STAR-CCM+. Possible explanations for this include that properties (Young’s modulus and thermal expansion) of the solid continuum were manually altered to adjust for the presence of water and the increase in pressure experienced here compared to the reference, which causes the displacement to in- crease further. The biggest difference is, once again, in the cooling stage, where the simulated displacement has become inverted. This is most likely due to the larger negative pressure that occurs at this time together with the, almost, non-existent temperature gradient. Figure 4.14: Displacement (z-direction) as function of distance from Reference 2 (left) and here simulated in STAR-CCM+ (right) 26 4. Results and Discussion 4.4 Effect and Strength of the Coupling In order to understand how much the different coupling mechanisms affect each other, a simulation was run with the hydro-mechanical coupling removed. Hence, the thermo-hydraulic (TH) and thermo-mechanical (TM) coupling were tested in- dividually. The time-dependent results are shown for pressure and displacement, respectively, in the two figures of Figure 4.15. As seen, the pore pressure increase/de- crease is stronger without the impact of the solid displacement, which is correct since the coupling term acts as a sink term in the mass balance equation. Comparing the two curves, it becomes clear that the impact from the displacement affect mostly the peaks of pressure, i.e. when the gradient of the displacement is largest. It is important to remember that the differences between the simulations most likely is smaller, since the displacement as mentioned is larger than in reality and the larger the displacement, the stronger will it affect the pressure. Finally, the displacement decreases without the impact of pore pressure. The increase due to the pressure gradient in the THM coupled simulation does not seem to be instant and hence stems from pore pressure dissipation. The impact is, however, to consider moderate why the displacements due to thermal expansion clearly are dominant. Comparing the two curves it can also be determined that the reason for the displacements being overestimated is not due to the high pressure, since the values experienced in the reference cases are surpassed even in the decoupled simulation. Figure 4.15: Pressure and displacement as functions of time for TH/TM coupling compared to THM coupling 4.5 Sensitivity Study A numerical simulation of a thermo-hydro-mechanically coupled problem has many sensitive parameters. In this study, two presumably particularly critical parameters were chosen which are the power of the heat source and the permeability of the rock. The heat source is the instigating load in this simulation, why it was decided to see how the results would differ if the power is increased respectively reduced by a factor of two. Due to the time limitation of the project, this was only tested in a 27 4. Results and Discussion thermo-hydraulic coupled simulation (FEA most time consuming here), hence the reaction of the solid skeleton in terms of displacement/stresses due to a change of the heat source was not investigated. In the two figures in Figure 4.16, temperature and pressure are shown as function of time and power of the heat source. When examining the rise and reduction of each stage for temperature and pressure it was found that the relation to the change in power is linear, i.e. twice the power results in twice the rise/reduction of each heating stage. Figure 4.16: Temperature and pore pressure as function of time for different powers The main source of the pore pressure increment lies in the low permeability of the clay, causing the water to expand further than the solid skeleton. Since the perme- ability is not a fixed material parameter, it was thought interesting to investigate how a change would affect the results. Acceptable values of permeability for this specific clay normally ranges between 1 ·10−19 - 1 ·10−20, but due to the weight of the parameter, a wider range was tested (1 ·10−17 - 1 ·10−22). Once again, solely thermo- hydraulic coupling was simulated and the results are shown below in Figure 4.17. As should be, the temperature is (close to) unaffected while pressure experiences an ex- treme difference depending on permeability. For the high permeability, κ = 1 ·10−17, there is an almost non-existing pore pressure increment, while κ = 1 · 10−22 lead to pressures above 30 MPa. Hence, the simulation is, as expected, highly sensitive to the permeability, even within the reasonable range where the difference in maximum pressure is well above a factor two. 28 4. Results and Discussion Figure 4.17: Temperature and pore pressure as function of time for different values of κ (TH coupling) 4.6 Sources of Error CFD and CSM are often used to, in a time and money efficient way, predict reality. While the objective, of course, is to come as close to reality as possible, trade-offs that, due to time/knowledge/computer resources, lead the results to differ from can be hard to avoid. This thesis study has been no exception and below follow a num- ber of possible sources of error. Assumptions and Simplifications: Due to the time limitation and the stated difficulties of THM-modelling many as- sumptions and simplifications were made to reduce the complexity of the simulation, such as the assumption of isotropy and fully saturated conditions as well as the fact that various properties were assigned constant values, e.g. porosity and permeabil- ity. While all the above may be part of the reason why the simulations in this study vary from reality it should not cause any major differences from the reference simulations. Both reference cases used fully saturated conditions, constant perme- ability and had some axisymmetric simulations that applied isotropy. Hence, these simplifications are not assumed the main culprit of the errors. However, it is still considered something to be changed/added in the future. STAR-CCM+: The assigned software does not have a stress model that is compatible with a porous medium. This means that the stress analysis was made on an entirely solid domain which, even though properties were assigned values with the intention of account- ing for the water, might be the reason behind the increased displacements. Solely Young’s modulus and the thermal expansion were available to adjust and it is not impossible that more is required. General numerical set-up and problems Throughout the thesis work, many numerical problems were encountered and dealt with. While the root for the most critical ones causing divergence or completely un- 29 4. Results and Discussion realistic values, have been found it is not to be excluded that there still exist settings that are not ideal. In a coupled problem like this, a larges amount of parameters and settings may affect multiple aspects of the simulation. Deviation from reference cases While this is not an error per se, there is a possibility that differences between the simulation performed here and the reference case may have lead to a bigger deviation in pressure and displacement than expected. Differences include that the two heaters in the experiment were joined to one as well as some variation in initial conditions and boundary conditions (for example varying temperature at one boundary). 30 5 Conclusion The ambition throughout the thesis work was to create a numerical model within the software STAR-CCM+ that can simulate thermo-hydro-mechanical coupling and obtain reasonable results. This was, to some extent, accomplished. A model was built where the thermal field is well captured and properties such as pressure and displacement follow the same behaviour, but show less agreement with the reference experiment and simulations in terms of magnitudes. However, the fact that numer- ous research groups world wide have been working on this problem for decades in order to validate its accuracy is an indication of the difficulties of the THM coupling, why the results are considered fairly satisfactory. Both pressure and displacement have, independently, been overestimated which was proven by the decoupling of the hydro-mechanical interaction. The displacement at this point still was too large and the pressure would increase even further if the displacement was corrected. The reason behind the deviating results has on the other hand not been determined, but only speculated about. Possible reasons in- clude un-ideal parameters and settings, pore pressure dissipation not being fully captured as well as the possibility of STAR-CCM+ not being optimal for THM coupled problems. While it has been determined that a THM model is possible to create, aspects such as the fairly new solid stress model leaves much to be desired as it for examples is not compatible with a porous medium. Additionally the software is not ideal for these low-velocity simulations since it is, according to the software support, originally created to solve flows with velocities above 1 m/s. If this in fact causes the differences is difficult to say, but it is a possible cause, why perhaps a software originally aimed for solid mechanics is preferable (low porosity and almost stagnant problem), at least as a mean of excluding possible errors. Ultimately, it has been understood that trouble-shooting a three-way coupled model is difficult, since it can be problematic to identify the parameter(s) causing the errors when they are affected by each other. 5.1 Future Work Based on the above discussion about STAR-CCM+, a decision should first be made if to continue with the software. An alternative could be the software Abaqus FEA, as it also is a commercial all-in-one solution that has shown accurate results for this purpose previously. If decided to continue in STAR-CCM+, one needs to get to the 31 5. Conclusion bottom of the high displacements and pressures. Several mechanisms/behaviour could, or possibly should, be included into the model, in whatever software is used, to obtain more accurate and interesting (three last items in list below) results. Below follows a list of possible additions to the model: • Anisotropy • Possibility for parameters such as porosity, permeability, Yong’s modulus and thermal conductivity to vary • Effect of moisture change and drainage • Effect of heater orientation • Introduction of engineering barriers 32 Bibliography [1] World Nuclear Association ,"Storage and Disposal of Radioactive Wastes", London, United Kingdom, May 2017. [Online]. Available: http://www.world- nuclear.org/information-library/nuclear-fuel-cycle/nuclear-wastes/storage- and-disposal-of-radioactive-wastes.aspx [2] B. Grambow, "Geological Disposal of Radioactive Waste in Clay", ELEMENTS, vol. 12, no.4, pp. 239–245, August 2016. [3] J. Rutqvist, L. Börjesson, M. Chijimatsuc, A. Kobayashic, L. Jingd, T.S. Nguyene, J. Noorishada and C.-F. Tsanga, "Thermohydromechanics of partially saturated geological media: governing equations and formulation of four finite element models", International Journal of Rock Mechanics & Mining Sciences, vol. 38, no. 1, pp. 105-127, January 2001. [4] Mont Terri Project, Switzerland. [Online]. Available: https://www.mont- terri.ch [5] O. Kolditz, U-J. Görke, H. Shao and W. Wang, Thermo-Hydro-Mechanical- Chemical Processes in Porous Media: Benchmarks and Examples., Heidelberg, Springer, 2012. [6] J. Noorishad, C.-F. Tsang and P.A. Witherspoon, "Coupled ther- mal–hydraulic–mechanical phenomena in saturated fractured porous rocks: nu- merical approach", Journal of Geophysical Research, vol. 89, pp. 10365–1037, 1984. [7] Y. Ohnishi and A. Kobayashi, "THAMES", in: O Stephansson, L. Jing, C.-F. Tsang (Eds.) "Coupled Thermo-Hydro-Mechanical Processes of Fractured Me- dia", Developments in Geotechnical Engineering, vol. 79, pp. 545-549, Elsevier, 1996. [8] T.S. Nguyen, "Description of the computer code FRACON" in: O Stephansson, L. Jing, C.-F. Tsang, "Coupled Thermo-Hydro-Mechanical Processes of Frac- tured Media", Developments in Geotechnical Engineering, vol. 79, pp. 539-544, Elsevier, 1996. [9] Q. Liu, C. Zhang, X. Liu. "A practical method for coupled THM simula- tions of the Yucca Mountain and FEBEX case samples for task D of the DECOVALEX-THMC Project", in Proc. GEOPROC2006 Int. Symp.: Second International Conference on Coupled Thermo-hydro-mechanical-chemical pro- cesses in Geosystems and Engineering, HoHai University, Nanjing, China, May 22-25 2006, pp. 220-225. [10] JI. Israelsson, "Short description of FLAC version 3.2", in: Stephansson O, Jing L, Tsang C-F, editors. "Coupled thermo-hydro-mechanical processes of frac- 33 Bibliography tured media", Developments in Geotechnical Engineering, vol. 79, pp. 513±22, Elsevier, 1996. [11] D. S. Chandrasekharaiah and L. Debnath, Continuum Mechanics, Academic Press Boston, 1994. [12] STAR-CCM+ Theory and User Guide, Version 12.02. CD-Adapco, New York, 2017. [13] B. Sundström, "Jämviktsekvationer - rörelseekvationer" in Handbok och formel- samling i Hållfasthetslära, Institutionen för hållfasthetslära, KTH, Stockholm, 1999. [14] K. Terzaghi, "Erdbaumechanik auf Bodenphysikalischer Grundlage", Franz Deuticke, Liepzig-Vienna, 1925 [15] M.A. Biot, "General Theory of Three-Dimensional Consolidation", Journal of Applied Physics, vol. 12, no. 2, pp. 155-164, February 1941. [16] O. Kolditz and W. Wang, "Object-oriented finite element analysis of thermo- hydro-mechanical (THM) problems in porous media", International Journal for Numerical Methods in Engineering, vol. 69, no. 1, pp. 162-201, January 2007. [17] W. Wang, G. Kosakowski and O. Kolditz, "A parallel finite element scheme for thermo-hydro-mechanical(THM) coupled problems in porous media", Comput- ers and Geosciences, vol. 35, no. 8, pp. 1631-1641, August 2009. [18] O.C. Zienkiewicz and R.L. Taylor, The Finite Element Method for Solid and Structural Mechanics, Elesevier Butterworth-Heinemann, Oxford, 6th edition, 2006. [19] S. Wasserman, "Fluid Software STAR-CCM+ Adds Solid Mechanics", June 2015, [Online]. Available: http://www.engineering.com/DesignSoftware/DesignSoftwareArticles/ ArticleID/10333/Fluid-Software-STAR-CCM-Adds-Solid-Mechanics.aspx [20] A. Gens, J. Vaunat, B. Garitte and Y. Wileveau, "In situ behaviour of a stiff lay- ered clay subject to thermal loading: observations and interpretation", Geotech- nique, vol. 57, no. 2, March 2007. [21] B. Garitte, B.J. Graupner, C. Lee, K. Maekawa, C. Manepally, P. Pan, J. Rutqvist and W. Wang, "The Mont Terri HE-D Experiment as a Benchmark for the Simulation of Coupled THM Processes", In. Proc. Intern. workshop on geomechanics and energy-deep geological disposal of radioactive waste, Lau- sanne, 26-28 November 2013. [22] M. Peric, "Flow simulation using control volumes of arbitrary polyhedral shape", ERCOFTAC Bulletin, no. 62, pp. 25-29, September 2004. [23] A. Gens, K. Wieczorek, I. Gaus, B. Garitte, J.C. Mayor, K. Schuster, G. Ar- mand, J.L. García-Siñeriz and T. Trick, "Performance of the Opalinus Clay under thermal loading: experimental results from Mont Terri rock laboratory (Switzerland)", Swiss J Geosci, no. 110, pp. 269–286, February 2017 [24] S. Moaveni, Engineering Fundamentals: An Introduction to Engineering, 4th ed., Indian Institute of Technology Kanpur, Global Engineering, 2011. [25] P. Bossart and M. Thury, "Characteristics of the Opalinus Clay at Mont Terri", Swiss Geological Survey, Wabern, Tech. Rep., 2008 34 List of Figures List of Tables Introduction Background Objectives Delimitations Theory Nuclear Repositories and Related Studies Repositories Related Studies Continuum Mechanics Porous Media Hydro-Mechanical Coupling Thermo-Hydro-Mechanical Coupling Governing Equations Fluid Process Deformation process Heat transport Discretisation Method Thermo-Hydro-Mechanical Coupling in STAR-CCM+ Methods and Numerical Set-Up Computational Domain Meshing Numerical Models and Solver Set-Up User Defined Equation of State Properties of the Porous Medium Boundary Conditions and Initial Conditions Porous Domain Solid Domain Coupling Simulation Steps Post-Processing and Validation Results and Discussion Thermal Field (THM) Comparison to HE-D in-situ test and Previous Simulations Pore Pressure (THM) Comparison to HE-D in-situ test and Previous Simulations Displacement (THM) Comparison to HE-D in-situ test and Previous Simulations Effect and Strength of the Coupling Sensitivity Study Sources of Error Conclusion Future Work Bibliography