Lattice Solid/Boltzmann Microscopic Model to Simulate Solid/Fluid Systems—A Tool to Study Creation of Fluid Flow Networks for Viable Deep Geothermal Energy

0. INTRODUCTION

In the forthcoming years, modelling of hydraulic fracture and flow of fluid and heat in fragmented material will play a vital role in the exploitation of geothermal energy. The deep geothermal energy extraction process involves drilling deep boreholes into hot rocks and injecting water under high pressure which decreases the effective friction on closed (impermeable) inter-meshing joint surfaces, causing them to "fracture" and slip in response to the existing tectonic pre-stress of the region. The fracture and slip of these joints creates high permeability pathways for fluid flow. The thermal energy trapped in the rocks within this region can then be extracted by circulating water between the injection and production wells through this permeable fracture system. Due to the high expense of drilling and inability to access underground wells, there are huge risks and uncertainties associated with this process. One of the challenges in geothermal energy extraction is how to stimulate and sustain the flow of water through the geothermal field and how to generate an efficient hydraulic subsurface heat exchange (fracture) system.

A fully developed thermo-hydro-mechanical coupling model and code which includes the most important physical mechanisms would provide means to study the phenomena and meet this challenge with minimal cost and risk.

Many approaches have been proposed for simulating fragmentation of solids using continuum methodology. While continuum based models can simulate discontinuities to an extent either replacing the discontinuities with material of a different rheology, or through special treatments of the discontinuity nodes, they cannot be used to study emergent behavior and probe the evolution of fracture systems which is a consequence of microscopic processes. Particle based models such as the discrete element method (DEM) and lattice solid model (LSM) naturally overcome such difficulties since displacements and detachment of solid fragments can be simulated.

When modeling fluid flow, the classical continuum approach is based on the numerical solution of Navier-Stokes (N-S) equations (i.e., computational fluid dynamics). As in the case of solids, in addition to the continuum approach based on numerical solution of the N-S equations, there are also other microscopic and mesoscopic approaches, such as the molecular dynamics (MD) method, the lattice Boltzmann method (LBM), and the smoothed particle hydrodynamics (SPH) method (Gingold and Monaghan, 1977). The LBM method is based on kinetic gas theory and simulates fluid flows by tracking the evolution of the single fluid particle distributions. The LBM is a semi-microscopic approach that models particle distributions rather than individual particles. Some advantages of LBM over the classical N-S approach include its ease in implementation and parallelization, and its ability to handle boundary conditions of complicated geometries (Chen and Doolen, 1998).

This paper introduces the coupled LSM/LBM, a new fully microscopic model based on coupling the LSM with the LBM thereby allowing simulations including solid-fluid coupling, fracture, and thermal fluid flow, and hence, to investigate the creation, evolution and efficiency of geothermal reservoir fracture systems as emergent phenomena.

We believe that the coupled LSM/LBM is a natural and elegant approach from a conceptual and physical point of view for simulating complicated solid-fluid systems aimed at studies of emergent behaviour than any approach which involves coupling microscopic and macroscopic methods for solid dynamics and fluid flow.

An interesting alternative to the coupled LSM/LBM is the coupled LSM/SPH (Komoróczi et al., 2013). While it has less conceptual and physical appeal in our view due to involving a coupling between a microscopic and macroscopic method, it is attractive for other reasons such as ease of implementation of the SPH within the LSM by making use of the LSM particle implementation as a basis for SPH nodal points, and possibly greater computational efficiency for some physical problems (although this is speculation at this stage).

1. HISTORICAL REVIEW OF LSM DEVELOPMENT

The lattice solid model (LSM) was first developed by Peter Mora in 1992 at Institute de Physique du Globe (Mora and Place, 1993; Mora, 1992). The LSM was motivated by molecular dynamics principles with either particle interactions via Lennard-Jones potential functions or elastic-brittle bonds. The LSM shares similarities with the Discrete Element Mode (DEM) but involves a different computational approach that is based on the use of large-scale supercomputer simulation via domain decomposition and parallel programming. The model was initially developed to study fracture, friction, rock rheology, tectonic processes and earthquake dynamics (Mora and Place, 1994, 1993; Mora, 1992). As more physics was incorporated including heat flow, thermo-porous coupling, thermo-mechanical feedback, and additional degrees of freedom in 3D (tension, shear, torsion and bending) to ensure realistic fracturing behaviour, it became more powerful and was calibrated with laboratory experiments and used to study outstanding scientific research questions and the complex nonlinear behaviour of discontinuous solids including the heat flow paradox in geophysics (Place and Mora, 2000; Mora and Place, 1998). The LSM has since been developed as a parallel open-source C++ software system called ESyS_Particle. The initial team working on the LSM consisted of Peter Mora and David Place, but rapidly grew into a major particle research group at Earth Systems Science Computational Centre (ESSCC) that was directed by Peter Mora at the University of Queensland. Major contributors included Steffen Abe, Yucang Wang, Fernando Alonso-Marroquin, and Shane Latham, among others. The development of the ESyS-particle software infrastructure became funded by the Australian Computational Earth Systems Simulator Major National Research Facility (ACcESS MNRF) established by Peter Mora and a national group of solid earth systems simulation scientists. The ESyS-Particle software is available at: https://launchpad.net/esys-particle/ESyS-particle, which is a site maintained by Dion Weatherley with software developments contributed from around the world and by present and past ESSCC researchers.

The initial particle simulation program of 1992 was written by Peter Mora in Connection Machine Fortran (CMF was a precursor to HPF and F90) at the IPG in Paris and was called LSMearth. David Place worked on an accurate way to simulate inter-particle friction (Place and Mora, 1999), and Peter Mora then added frictional heating and heat flow. Subsequently, Peter Mora, David Place and Steffen Abe added random particle sizes (Place and Mora, 2001), thermo-porous and thermo-mechanical coupling effects (Abe et al., 2000), and proposed an efficient scheme to model arbitrary shaped angular particles (Abe and Mora, 2003). LSMearth was redeveloped in C++ by David Place, and subsequently parallelised by Steffen Abe et al. (Abe et al., 2004) and Shane Latham et al. (Latham et al., 2005) using MPI into a highly efficient and usable tool, with various new micro-physics being added by a group of researchers and computational scientists under the direction of Peter Mora as part of the ACcESS MNRF. Around 20 researchers worldwide are currently contributing to ESyS-particle. Some of the most critical recent developments have been made by Yucang Wang and Fernando-Alonso Marroquin respectively to model the full set of degrees of freedom (six kinds of independent relative movements are transmitted between two 3-D interacting particles) needed to accurately simulate observed fracture patterns (Wang, 2009; Wang and Alonso-Marroquin, 2009; Wang and Mora, 2009, 2008b; Wang et al., 2006), and how to efficiently simulate complex particle shapes (Alonso-Marroquin and Wang, 2009; Alonso-Marroquin et al., 2007; Abe and Mora, 2003). One of the major advantages of the discrete (LSM) model over the continuum models is that the large deformation and dynamics process associated with fracture phenomena can be easily modelled. It has been successfully utilised to the study of physical process such as rock fracture (Wang and Alonso-Marroquin, 2009; Wang and Mora, 2008b; Place et al., 2002; Place and Mora, 2001; Mora and Place, 1993; Mora, 1992), the stick-slip frictional instability and earthquake dynamics (Mora and Place, 1999, 1993), the heat-flow paradox (Alonso-Marroquin et al., 2006; Mora and Place, 1999, 1998; Mora et al., 1997), accelerating acoustic emissions/energy release, load-unload-response-ratio and critical sensitivity prior to catastrophic failure and simulated earthquakes (Mora et al., 2002, 2000), earthquake physics and predictability (Mora and Place, 2002; Mora et al., 2000), localisation phenomena (Place and Mora, 2000; Mora and Place, 1999, 1998), fault gouge evolution (Mair and Abe, 2008; Alonso-Marroquin et al., 2007, 2006; Mora et al., 2000; Mora and Place, 1999) and comminution in shear cells (Mair and Abe, 2008).

2. FRAMEWORK TO SIMULATE DEEP GEOTHERMAL RESERVOIR SYSTEMS

In order to simulate the creation, evolution, and long-term energy output of a deep geothermal reservoir system, one must be able to simulate complex nonlinear dynamical processes in complicated solid-fluid systems. An ideal simulation may involve first simulating the joint and fracture patterns in the region of interest by modelling the tectonic stress and thermal history of a region. Subsequently, one may simulate the drilling of boreholes into the rock, injection of water under high pressure, the ensuing hydraulic fracturing process, the subsequent long-term evolution of the fluid-flow pathways, and the overall energy output of the region before either the heat of the region is fully depleted, or no further energy can be produced due to blocking of flow pathways. Figure 1 depicts this simulation strategy.

Figure  1.  Illustration of simulation based framework to study creation of fluid flow networks for viable deep geothermal energy. Upper left. creation of a pre-existing joint/fracture system via microscopic LS/LB simulation of tectonic stress and thermal history acting on a rock volume (left. real fracture systems; right. simulated fracture systems. The upper and lower images are at different scales). Upper right: Drilling, high pressure water injection, and hydraulic fracturing via microscopic LS/LB simulation. Lower left. microscopic LS/LB simulation to study the evolution of fractures and their effective permeability as they slip, grains rub off, fluids flow, and chemical processes occur over a long time scale. Lower right. macroscopic simulation to calculate the long-term energy output of the geothermal reservoir/fracture system predicted by the micro-scale simulation using the permeability evolution of fractures derived from micro-scale simulation (Xing and Mora, 2006).

DownLoad: Full-Size Img PowerPoint

During the first stage of creating a fracture system, water is injected at the "injection wells" and this causes a decrease in the effective normal stress on the initially closed joint surfaces around these wells. This then allows these joints to slip and partially release the tectonic pre-stress in the region. This induced seismicity is used in the field to help map the fracture system which provides a pathway (or pathways) for fluid flow. However, it would be costly to vary the injection characteristics by trial and error in the field until a suitable fracture system with multiple fluid flow pathways could be developed. Multiple pathways are crucial for any given field to be fully produced, that is, to fully deplete the anomalously high heat for a given region such that the investment of the injection and production wells is has been fully covered and the field is economically viable.

Simulations may provide a means to study how to vary injection characteristics such as depths and pressures, in order to create multiple fractures and hence, flow pathways between injection wells and production wells. However, such simulations would need to accurately simulate fracture and other emergent behavior in complicated solid-fluid models. Here, we propose use of the LSM for simulation of the solid dynamics due to its proven ability to accurately reproduce realistic fracture behavior as illustrated in the following sections.

3. REALISTIC FRACTURING BEHAVIOR USING THE LSM

Figure 2 shows the 3D complex fracturing behaviour of a leftwards travelling impacting ball in a laboratory experiment (Khanal et al., 2008) compared to a numerical simulation using the LSM by Yucang Wang. Figure 3 shows a 3D simulation result of axial breakage compared to a laboratory experiment of axial breakage (see Wang and Alonso-Marroquin, 2009, 2008).

Figure  2.  Laboratory result of an impacting ball by Khanal et al. (2008) (left) compared with numerical result using the LSM by Wang and Alonso-Marroquin (2008) (right).

DownLoad: Full-Size Img PowerPoint

Figure  3.  Numerical simulation of axial breakage (top) compared to laboratory experiment (bottom).

DownLoad: Full-Size Img PowerPoint

The additional examples shown in Fig. 4 of the simulation of the development of brittle fracture systems also exhibit a close match to the patterns observed in laboratory experiments (Wang and Mora, 2009, 2008b), and provides a further illustration that the lattice solid model and software system, ESyS-Particle, is capable of matching fracture patterns seen in laboratory experiments. Hence, implementation of the lattice solid model including thermo-porous flow, frictional heat generation and heat flow, thermo-mechanical coupling effects and realistic solid-fluid dynamics such as hydraulic fracturing, provides a unique platform for scientific breakthroughs on how to create and sustain the multiple fluid flow pathways required for proof of concept and development of the new geothermal green energy industry.

Figure  4.  Brittle failure experiments—laboratory result versus simulation result.

DownLoad: Full-Size Img PowerPoint

4. TENSION, SHEAR, TORSION AND BENDING DEGREES OF FREEDOM

The accurate modelling of realistic fracturing behaviour that matches with experimental laboratory data was only possible by including rotational dynamics as new degrees of freedom (see Fig. 5). Details of the contact model are provided in the literature (Wang, 2009; Wang and Mora, 2009, 2008b).

Figure  5.  The 6 degrees of freedom for 3D particle interaction. F_r is force in radial direction, F_s1 and F_s2 are shear forces, M_t is twisting torque, and M_b1, M_b2 are bending torques.

DownLoad: Full-Size Img PowerPoint

Elastic properties in regular lattices have been derived in these cases (Wang and Mora, 2008a). As a test of fracture behaviour, we have accurately reproduced "wing cracks" (see Fig. 6) that are observed in the laboratory (Wang and Mora, 2008b). ESyS-Particle has unique capabilities to model such wing cracks as all six degrees of freedom per particle are modelled to include contact forces, bending and twisting moments (Wang, 2009; Wang and Mora, 2009). The 2D example of wing cracks to the right is from Wang and Mora (2008b).

Figure  6.  Simulation of wing cracks.

DownLoad: Full-Size Img PowerPoint

5. THE THERMAL LATTICE BOLTZMANN METHOD

To model dynamics of fluids in the cracks and pore space, we use the thermal energy distribution type BGK thermal lattice Boltzmann method (Hung and Yang, 2011; Guo et al., 2007; He et al., 1998) which has been shown to yield the Navier-Stokes equations for fluid flow combined with heat flow. Namely, in one time step Δt, the mass density of particles moving in the α-direction of a regular lattice moving with velocity c_α denoted f_α is updated as

where the first term on the right denotes the streaming of particles moving one lattice spacing Δx_α in the α-direction in one time step, and the second term on the right denotes the redistribution of mass density flow due to collisions where τ_f is the dimensionless relaxation time constant for the collision term which is related to kinematic viscosity via ${\nu _f} = ({\tau _f} - 0.5)c_s^2\Delta t$ where c_s is the speed of sound in the fluid. The collision term is calculated using

where Δf_c^α is the velocity change due to collision, f_α^eq is the equilibrium distribution expressed in Eq. 3. In this work, we model the D2Q9 BGK model (the two-dimensional, nine speeds Bhatnagar-Gross-Krook model) so the particle distributions travel with speeds of c_α=0 (α=0), c_α=Δx/Δt (α=1, 2, 3, 4) particles moving in the ±x and ±y directions and ${c_\alpha } = \sqrt 2 \Delta x/\Delta t{\rm{ }}$ (α=5, 6, 7, 8) for particles travelling in the diagonal directions. The equilibrium distribution is calculated using

where the equilibrium distribution weights are ${{w}_{\alpha }}=({}^{4}\!\!\diagup\!\!{}_{9}\;,{}^{1}\!\!\diagup\!\!{}_{9}\;,{}^{1}\!\!\diagup\!\!{}_{9}\;,{}^{1}\!\!\diagup\!\!{}_{9}\;,{}^{1}\!\!\diagup\!\!{}_{9}\;,{}^{1}\!\!\diagup\!\!{}_{36}\;,{}^{1}\!\!\diagup\!\!{}_{36}\;{}^{1}\!\!\diagup\!\!{}_{36}\;{}^{1}\!\!\diagup\!\!{}_{36}\;)$, ρ is density, u is the macroscopic velocity of the fluid. For this case, the RMS velocity, and hence, the speed of sound in the fluid is $c_{s}^{{}}={}^{1}\!\!\diagup\!\!{}_{\sqrt{3}}\;$. The macroscopic density and momentum are calculated using $\rho = \sum\limits_\alpha {{f_\alpha }} $ and $\rho \mathit{\boldsymbol{u}} = \sum\limits_\alpha {{f_\alpha }} {\mathit{\boldsymbol{c}}_\alpha }$. In the thermal lattice Boltzmann method, a second distribution g_α is introduced which relates to kinetic energy within the fluid and hence the heat. This is also modeled in the two steps of streaming and collision

where collision term is $\Delta g_\alpha ^C = g_{eq}^\alpha - {g^\alpha }$ and the equilibrium distribution for g_α is calculated using

The macroscopic internal kinetic energy and temperature are calculated using $\rho E = \sum\limits_\alpha {{g_\alpha }} {\rm{ }}$ and $E = DRT/2$ where D=2 is the number of dimensions, R is the gas constant, and T is the temperature. In the above equations, the different relaxation time τ_g allows the thermal diffusivity of the fluid to be controlled.

6. MECHANICAL COUPLING OF THE LATTICE SOLID AND LATTICE BOLTZMANN METHOD

To implement mechanical coupling between the LSM and LBM, the following issues need to be considered: moving boundary conditions for a curved solid-fluid interface; momentum transfer between solid particles and the fluid; and force transfer between fluid nodes and solid particles. Here Yu's moving boundary condition is adopted (Yu et al., 2003). The lattice node on the fluid side of the boundary is denoted as x_f and that on the solid side is denoted as x_b (see Fig. 7). The particle momentum moving from x_f to x_b is e_α and the revised momentum from x_b to x_f is ${e_{\tilde \alpha }} = - {e_\alpha }$. Here, x_w denotes the intersection of the wall with the lattice link. Due to the arbitrary position of the particles and the curved particle surface, the particle surface can intersect the link between two nodes at an arbitrary distance.

Figure  7.  The moving curved wall boundary condition.

DownLoad: Full-Size Img PowerPoint

To accurately capture the position of the particle surface, the fraction of an intersected link in the fluid region can be computed as (Fig. 6)

The reflected distribution function at nodes can be calculated using an interpolation scheme

The fluid force acting on the particle surface can be obtained using

where the first summation is taken over all fluid nodes at x_b adjacent to the particle, and the second is taken over all possible lattice directions pointing towards a particle cell. This force is added to the particle force in LSM code.

7. THERMAL COUPLING OF THE LATTICE SOLID MODEL AND LATTICE BOLTZMANN METHOD

The thermal coupling between the lattice solid model and lattice Boltzmann method has the same issues as the mechanical coupling and can be achieved in the same way. However, because the thermal diffusion time scale is longer than the time scale of solid and fluid dynamical processes (e.g., fracture etc.—solid; turbulence etc.—fluid), it is less sensitive to the precision of the implementation. On the other hand, the thermal LBM can become unstable and care must be taken in the implementation, particularly when the fluid is undergoing rapid dynamical interaction with a solid and hence, may be far from equilibrium during such processes. We have implemented heat transfer within the Lattice Solid using the thermal LBM with a third distribution h_α, and a different relaxation time τ_h which allows the solid to have a different thermal diffusivity than the fluid. This approach of modeling the heat flow in the moving lattice solid fragments (groups of bonded particles) consistently with the LBM leads to a stable approach for modeling heat in the solid-fluid system with the advantage of yielding relatively homogeneous numerical precision for the heat flow calculations within both solid and fluid regions (c.f., more heterogeneous precision had different methods been used for heat flow calculations within solid and fluid regions).

We use the thermal-energy-distribution type BGK thermal lattice Boltzmann method that has been shown to yield the NS equations for fluid flow combined with heat flow. Key issues of the LSM-LBM coupling: moving boundary conditions for a curved solid-fluid interface; momentum transfer between solid particles and the fluid; and force transfer between fluid nodes and solid particles. Moving boundary condition is adopted. Verification of the model has been achieved using simulation of particles in fluids, see Fig. 8. The example simulates dragging of particles by the fluids, convective and diffusive heat flow, and turbulent flow for high Reynolds numbers where particle surfaces act as boundary conditions.

Figure  8.  Snapshots of a coupled LSM-LBM simulation of initially hot particles in a cold fluid moving to the right through a channel. The colours represent temperature (left) and velocity (right). The dots within the flow are tracer particles to assist visualisation of the flow field. The edge of the circular particles are coloured dark blue on the left plot to assist in visualising the locations of the particles in the temperature field.

DownLoad: Full-Size Img PowerPoint

8. HYDRAULIC FRACTURE EXAMPLE

Figure 9 shows a 2D simulation of hydraulic fracturing. In this test, rock is modeled as 1 026 bonded LSM particles with a hole in the middle representing a borehole in which pressurized water is injected. Particle sizes are variable and range from 0.1 to 1 units. Water is simulated using the LBM. The minimum particle size is 2 times of the LBM grid size of the fluid. Fluid pressure in the centre of the hole is increased slowly to model the injection of water. The pressure increase is modelled by addition of a source term on the right-hand side of the standard LBM scheme. The cracks start from the surface of the borehole and propagate inside the rock due to the pressure of the fluid. Tensile fractures are dominant in the beginning of the crack initiation. Our simulations include most of the mechanisms of the hydraulic fracturing process: (i) mechanical deformation and fracturing induced by the fluid pressure; (ii) flow of fluid within the fracture that was generated; and (iii) fracture propagation.

Figure  9.  From top to bottom: snapshots of hydraulic fracture; fracturing event distribution, with blue for tensile fractures and red for shear fracture; velocity vector field of fluid flow during hydraulic fracturing; and fluid pressure during hydraulic fracturing (in this case dark green means higher fluid pressure).

DownLoad: Full-Size Img PowerPoint

9. HEAT FLOW EXAMPLE—COLD FLUID THROUGH HOT ROCK

We digitised a 2D porous sandstone to generate the lattice representing the rock. The digitised grains were shrunk to generate permeability that would be present in 3D, and were then initialised with a uniform high temperature. Space between grains was filled with a cold fluid with a velocity u=(0.1c_s, 0). Left and right boundary conditions imposed continued flow to the right and upper and lower boundaries were non-slip. Figure 10 depicts a snapshot of the temperature and velocity fields at t=700 time-steps of a simulation. The results illustrate the ability of the coupled lattice solid/lattice Boltzmann model to simulate the combined fluid/heat flow within a complex solid-fluid system.

Figure  10.  Temperature T (left) and fluid velocity magnitude |u| (right) at t=700 time steps in a simulation of an initially cold fluid with high thermal diffusivity, through a hot solid porous and permeable 2D sandstone matrix with low thermal diffusivity. Colours are red=hot/fast, blue=cold/slow. Lines are tracer particles that started as vertical columns. These allow fluid flow to be visualised simultaneously with temperature.

DownLoad: Full-Size Img PowerPoint

10. CONCLUSIONS

The lattice solid model particle based method has the ability to provide breakthroughs in the understanding of the mechanical properties of hot rocks, their efficient fracturing, and the ability to produce sustained fluid pathways through the induced fracture network. In situ measurements to assess the fracture of rocks require theoretical interpretation of the energy budget and fracture propagation. These theories can be incorporated into a continuum model by interconnecting the length scales of fragmentation. Here particle-based models can be used to extract data of fracture creation that are subsequently passed to the continuum models (i.e., using up-scaling). This ab-initio approach has the potential to exploit a model of large-scale fragmentation based on processes at the smallest scale which represents an independent source of data more closely related to the complexity of fragmentation.

The examples shown in this paper illustrate the potential, when fully developed, of the coupled lattice solid/lattice Boltzmann model to simulate the creation, dynamics, evolution and energy yield of realistic deep geothermal reservoir systems as emergent phenomena. Additional physical processes that may need to be incorporated include chemical reactions and precipitation since these may affect the viability of flow pathways. Carefully controlled laboratory tests can be used to validate the model. Also, if the temperature data and seismicity data are obtained in the real bores, the model results can be compared to experiments. Triggered seismicity can be modeled, with the potential energy in the bonds corresponding to the seismic energy released, and Richter scale can be defined accordingly by analyzing the frequency of the events. If enough real seismic data are obtained, the simulation results can be used to compare with the real data. Currently the simulations presented in this paper are just qualitative and stay at rudimental stages. More detailed and comprehensive work need to be done to produce more reasonable/realistic results.

Further development and application of the thermo- and hydro-mechanically coupled ESyS_Particle simulation software system will benefit geothermal exploitation, by delivering new multi-scale, numerical capabilities to model fracturing, heat flow, and solid-fluid coupling; and hence, to model the generation and sustainability of the hot-fractured rock geothermal energy fracture systems required to exploit this new green-energy resource. This may prove to be the key to realise sustainable green geothermal energy to contribute to power the world.