A novel DEM-based pore-scale thermal-hydro-mechanical model for fractured non-saturated porous materials

For fracture propagation, a novel DEM-based pore-scale thermal-hydro-mechanical model of two-phase fluid flow with heat transfer in non-saturated porous materials with low porosity was developed. Numerical computations were performed for bonded granular specimens, using a DEM fully coupled with CFD (based on a fluid flow network) and heat transfer, which integrated discrete mechanics with fluid mechanics and heat transfer at the meso-scale. Both the fluid (diffusion and advection) and bonded particles (conduction) were involved in heat transfer. The numerical findings of the coupled thermal-hydraulic-mechanical (THM) model were first compared to the analytical solution of the classic 1D heat transport problem. The numerical and analytical outcomes were in perfect agreement. Advection's impacts on the cooling of a bonded particle assembly were next numerically demonstrated for low and high Peclet numbers. Finally, the THM model's utility was proved in a thermal contraction test employing a bonded particle assembly during cooling, which resulted in the creation of a macro-crack. The effects of a macro-crack on the distribution of fluid pressure, density, velocity, and temperature were studied.


Introduction
The majority of physical phenomena in engineering challenges happen in non-isothermal environments. Furthermore, even if the physical system is initially in thermodynamic equilibrium, physical processes or chemical reactions may cause local temperature fluctuations, resulting in heat transfer. Only in a few circumstances may the engineering problem under investigation be reduced to isothermal conditions. As a result, numerous scientific disciplines and technical applications, such as environmental science, chemical and food processing, powder metallurgy, energy management, geomechanics, and geological engineering, are interested in understanding heat transfer in particle systems. In the analysis of many multi-field issues in porous and fractured materials, the requirement to address the effect of heat transfer becomes crucial. Diffusion, advection, and radiation are examples of heat transmission mechanisms. Many nonlinear processes affect complex coupled thermal-hydraulic-mechanical (THM) systems, which include heat transfer, fluid flow, and material deformation.
The most popular technique in THM models is a continuum approach, which is based on a mathematical framework linking sets of differential equations to explain thermodynamic, solid, and fluid mechanics laws, such as finite element [19,26,31,36,39,52], or finite difference implementations [37]. Even though they are appealing for macro-scale applications, continuum modelling approaches based on the finite element method (FEM) or the finite volume method (FVM) have significant computational and continuity limitations when applied to discontinuous and highly deformable media like packed or fluidized beds and fractured granular porous materials. Classical approaches have tremendous numerical challenges generating suitably fine meshes in porous media with low porosity (less than around 15%), such as concrete or rocks [1]. When modelling the crack initiation and propagation in porous materials with low porosity, the problem is amplified dramatically (e.g. concrete or rocks). Discrete techniques, such as the discrete element method (DEM) [13] or the finite-discrete element method (FDEM) [49,50], on the other hand, have proven to be successful in modelling the behaviour of particulate systems. The capacity to directly capture and predict fracture and micro-crack evolution through particulate materials is a benefit of DEM. Yan and Zheng [50] created a 2D thermo-mechanical model based on the FDEM (finite-discrete element technique) for modelling thermal rock cracking, while Yan et al. [49] created a 2D coupled thermal-hydro-mechanical model for simulating rock hydraulic fracturing. The ability to anticipate coupled THM processes at the meso-scale has recently been expanded because of DEM's strength in realistic simulations of the behaviour of fragmented particle systems. Various methods were adopted to combine DEM with fluid flow and heat transfer methodologies, to connect TH processes with DEM, direct numerical simulations (DNS) could be used. Different numerical approaches (for example, FEM and FVM) were applied to solve governing equations. Deen et al. [14] suggested a method for implementing submerged boundaries that did not require the use of an effective diameter. THM processes in dense fluidparticle systems were the subject of the approach. The proposed approach was confined to invariant geometries, topologies, and porosities of relatively high values (porosities greater than those of concrete or rock). DNS-DEM models, in reality, are limited to systems with fewer particles than CFD models. The lattice Boltzmann methods (LBMs) [10,18,51] rely heavily on the precise depiction of solid-fluid interfaces and have the same numerical and computational restrictions as DNS-DEM models.
To decrease computing costs and overcome numerical restrictions when utilizing DEM to THM processes in very dense fluid-particle systems with very low porosity (even below 5%), the fluid flow and heat transfer models must be reduced. Recently, most DEM-based THM models separated fluid flow within reservoirs (pores, macro-pores, preexisting cracks, etc.) from the flow between reservoirs. According to the assumption, the fluid in the reservoirs is compressible, but the fluid flowing between the reservoirs is incompressible. Cundall [11], Hazzard et al. [16], and Al-Busaidi et al. [2] were the first to introduce and explore the concept of simplification. In addition to highlighting fluid characteristics, fluid flow regimes are separated. To estimate the mass flow rates at the reservoirs' boundaries, the fluid flow regime in the reservoirs is stationary or close to stagnation, while the fluid flow between the reservoirs is laminar. In most cases, a Poiseuille flow model in pipes or between two parallel plates is utilized to predict mass flow rates. The pressure in the pore is calculated using either the assumed equation of state [2,16] or the Stokes equation [7,32]. All models assume a pure liquid or mixture flowing in a single-phase barotropic fluid flow (without tracking phase fractions). Different heat transfer models are combined with coupled DEM-CFD models in this technique. For each reservoir (cell) volume, Tomac and Gutierrez [46] computed the energy conservation equation. The selected energy conservation equation was related to energy transport in an incompressible fluid's laminar flow. Caulk et al. [8] proposed a more advanced 3D DEM-based THM model based on the framework of the pore-scale finite volume (PFV) scheme, which was first proposed by Chareyre et al. [9] and later extended by Scholtès et al. [38] for up-scaling compressible viscous flow and oriented towards saturated dense grain packing applications in geomechanics.
The purpose of this work is to show a new DEM-based pore-scale hydro-mechanical model of two-phase fluid flow in non-saturated porous materials with very low porosity that is enhanced by heat transfer (e.g. for concrete or rocks). Both the fluid (diffusion and advection) and the bonded particles (conduction) were involved in heat transfer. THM calculations were carried out numerically using a 3D DEM in conjunction with 2D CFD (based on a fluid flow network made up of channels) and 2D heat transfer, which connected discrete mechanics, fluid mechanics, and heat transfer at the meso-scale. Previously, a coupled DEM/CFD model based on the fluid flow network without heat transfer (formulated by the authors [24,25]) was used to describe hydraulic fracturing in unsaturated rock masses with one-or two-phase laminar viscous two-phase fluid flow containing a liquid and gas. This paper's DEM-based THM mesoscopic technique for modelling fluid flow and heat transfer has substantial advantages over other existing ones in the literature. Some of the advantages are as follows: 1. the precise tracking of water/gas fractions in pores, taking into account their variable geometry, size, and position, 2. an effective algorithm for automatically meshing and remeshing particle and fluid domains to account for changes in their geometry and topology, 3. the use of a coarse mesh of solid and fluid domains to generate a virtual fluid flow network (VPN) and to solve the energy conservation equation, 4. using FVM to solve the energy conservation equation in both domains on a very coarse mesh of cells, 5. to examine supercritical fluid flow, the corrected Peng-Robinson equation of state [34] was adopted for both fluid phases (necessary e.g. for studying THM processes in hydraulic fracturing problems), 6. calculating varying heat fields in pores and particles, 7. tracking virtual thermal deformation of discrete elements to precisely compute fluid volume changes over time, and 8. two-phase flow.
There are no coupled DEM-based thermal-hydro-mechanical models of multi-phase supercritical fluid flow in rocks that can be used to mimic THM processes during hydraulic fracturing. The numerical findings in the current paper are first compared to the analytical solution of the classic 1D heat transfer issue to validate the innovative fully coupled THM model. The numerical and analytical results are found to be perfectly in agreement. A numerical representation of the effects of advection on the cooling of a bonded granular bar specimen is also provided. Finally, the THM model's utility is proved in a thermal contraction test employing a bonded particle assembly during cooling, which results in the creation of a macro-crack. The effects of a macro-crack on fluid pressure, density, velocity, and temperature are investigated. The findings of this study can be directly applied to a variety of geothermal applications.
The structure of the current study is outlined below. Following Sect. 1, a mathematical model of the DEMbased linked thermal-hydro-mechanical method is provided in Sect. 2. The calibration of the model is detailed in Sect. 3. The validation of the THM model is described in Sect. 4. The influence of advection on the cooling of a bonded granular assembly is examined in Sect. 5. A damage process in a bonded granular assembly during a thermal contraction test is discussed in Sect. 6. Finally, in Sect. 7, some closing remarks are made. The THM model has been incorporated by the authors into the YADE DEM code, which is an opensource software [20,41].

Thermo-hydro-mechanical model
The THM model's original concept is based on the notion that in a physical system, two domains coexist: the 3D discrete (solid) domain and the 2D continuous (fluid) domain. The solid domain is made up of a single layer of discrete 3D elements (spheres), whereas the fluid domain was two-dimensional (Fig. 1a). Discrete elements are placed so that their centres of gravity are in the middle of the plane (2D surface). They are projected next onto the plane to form circles (Fig. 1b). As a result, discrete element equations of motion are solved in the 3D discrete domain, whereas fluid flow equations are solved in the 2D fluid continuous domain (red colour in Fig. 1), and heat transfer equations are solved in the 2D fluid and solid continuous domains (red and grey colours in Fig. 1b). Despite the use of a 3D DEM solver, the entire DEM-based THM model is two-dimensional.

DEM for cohesive-frictional materials
DEM calculations were performed with the 3D spherical explicit discrete element open code YADE [20,41] that allows for a small overlap between two contacted bodies (soft-particle model). Utilizing an explicit time-stepping scheme, particles in DEM interact with one another during translational and rotational motions using a contact law and Newton's 2nd law of motion [13]. A cohesive bond at the grain contact is postulated in the model, with brittle failure beneath the critical normal tensile force. Under normal compression, shear cohesion failure causes contact slip and sliding, which follows the Coulomb friction law. The mechanical response of DEM is presented in Fig. 2. The DEM equations are below F s À F s max À F n Â tan l c 0 ðbefore contact breakageÞ; ð4Þ where F ! n -the normal contact force, U-the overlap between discrete elements, N ! -the unit normal vector at contact point, F ! s -the tangential contact force, F ! s;prev -the tangential contact force in the previous iteration, X ! s -the relative tangential displacement increment, K n -the normal contact stiffness, K s -the tangential contact stiffness, E c -the elastic modulus of the particle contact, v c -the Poisson's ratio of particle contact, R-the particle radius, R A and R B contacting particle radii, l c -the Coulomb inter-particle friction angle, F s max -the critical cohesive contact force, F n min -the minimum tensile force, C-the cohesive contact stress (maximum shear stress at pressure equal to zero), and T-the tensile normal contact stress, F ! k damp -the dampened contact force, F ! k and v ! k p -kth-the components of the residual force and translational particle velocity v p and a d -the positive damping coefficient smaller than 1 (sgn(Á) that returns the sign of the kth component of velocity) [12]. Non-viscous damping was applied [12] to accelerate convergence in quasi-static simulations (Eq. 7).
The following material constants: E c , t c , l c , C and T are required for DEM simulations. In addition, R, q, (mass density) and a d are needed. The C/T-ratio is crucial to adequately simulate the failure type of specimens (brittle or quasi-brittle), the distribution of shear and tensile cracks, and the ratio between the uniaxial compressive and tensile strength. The material constants are usually identified by running a series of DEM simulations and comparing them with experimental results of simple tests (e.g. uniaxial compression, triaxial compression, simple shear) [28]. The damping parameter is always a d = 0.08. For this value, the loading velocity v does not affect the results [28]. Damage occurs if a cohesive joint between spheres (Eq. 6) disappears after reaching a critical threshold. If any contact between spheres after failure re-appears, the cohesion does not appear more (Eq. 5). Note that material softening is not considered in the DEM model. Although bonds can break by shear, the essential micro-scale mechanism for damage in the pre-failure regime is bond damage in tension. An arbitrary micro-porosity might be achieved in DEM since the particles may overlap. The fracture is not allowed to propagate through aggregate, i.e. the particle breakage is not taken into account. The model was successfully used by the authors for describing the behaviour of different engineering materials with a granular structure (mainly of granular materials [21][22][23]48] and concrete materials [29,30,40,[42][43][44] by taking shear localization and fracture into account).

Fluid flow model
The general concept of a fluid flow algorithm using DEM and a channel network is adopted from [2,11,16]. In this concept, fluid flow is simulated by assuming that each particle contact is an artificial flow channel (between two parallel plates in 2D or along a duct in 3D) and those artificial channels connect real reservoirs in the particulate medium (pores, fractures, and pre-existing cracks) that store fluid pressures. Thus, the pressure in reservoirs depends both on the mass transported along channels from/ to other reservoirs and the volume changes of reservoirs. Since the volume of reservoirs changes due to the material deformation (described by discrete elements in DEM), the fluid density must also change (the fluid in reservoirs is compressible). Thus, the fluid moves in channels while the reservoirs solely store pressure. The artificial channels create a fluid flow network. The fluid flow in those artificial channels is characterized by a simplified laminar flow of an The fluid model in the current paper significantly differs from the general concept [2,11,16]. The reservoirs (pores, cracks, pre-existing cracks, etc.) store now not only pressures but also phase fractions, fluids densities, energy, and temperature. The continuity equation is employed to compute the density of fluid phases stored in reservoirs. The fluid-phase fractions in reservoirs are computed by applying the equation of state for each phase assuming that fluid phases share the same pressure (as in the Euler model of multi-phase flow). The mass flow rate in artificial channels of a fluid flow network is now estimated by solving continuity and momentum equations for the laminar flow of incompressible fluid. The gravity centres of the 3D spheres are located on the XOY plane. The 3D spherical particles are projected onto the 2D mid-plane and next discretized into the 2D polygons [24]. To get a more realistic distribution of the unknown variables (pressure, fluid-phase fractions, and densities), a remeshing procedure discretizes the overlapping circles, determines the contact lines, and deletes the overlapping areas [24]. As a result, each reservoir is discretized into a number of triangles (in 2D problems). Each triangle in the fluid domain is called the Virtual Pore (VP) (Fig. 3). The artificial channels connect the gravity centres of triangles (VPs) to create a fluid flow network called the Virtual Pore Network (VPN). The displacements of discrete elements in the perpendicular direction OZ and the rotations around the axes OX and OY are fixed. All forces converted from pressure and shear stress are applied to the spheres. The forces are computed, based on the pressure and shear stress for the specimen thickness (in the OZ direction) equal to the largest diameter of spheres. All numerical parameters might be used without the necessity of re-calibration. VPs accumulate pressure and store both fluid-phase fractions and densities, energy, and temperature. The mass change in VPs is related to the density change in a fluid phase that results in pressure variations. The equation of momentum conservation is neglected in triangles but the mass is still conserved in the This algorithm is repeated for each VP in VPN and each solid cell (stage 'd') using an explicit formulation. According to the above algorithm, incompressible laminar two-phase fluid (liquid/gas) flow under non-isothermal conditions is assumed in the channels of the fluid flow network. The liquid and gas initially exist in the material matrix and pre-existing discontinuities. Two types of channels are introduced in the Virtual Pore Network [24] ( Fig. 3): (A) the channels between discrete elements of the material matrix in contact (called the 'S2S' channels) and (B) the channels connecting grid triangles in pores that touch each other by a common edge (called the 'T2T' channels). The channel length is assumed to be equal to the distance between the gravity centres of adjacent grid triangles. In real 3D problems, the fluid flows around the  [24] spheres in contact. However, in 2D problems, there is no free space for fluid flow. Therefore, the concept of virtual S2S channels is introduced [24]. The hydraulic aperture h of the channels 'S2S' is related to the normal stress at the particle contact by a modified empirical formula of Hökmark et al. [35]: where h inf is the hydraulic aperture for the infinite normal stress, h 0 -the hydraulic aperture for the zero normal stress, r n -the effective normal stress at the particle contact and bthe aperture coefficient. The hydraulic aperture of the channel type 'T2T' is directly related to the geometry of the adjacent triangles as where e-the edge length between two adjacent triangles, xthe angle between the edge with the length e and the centre line of the channel that connects two adjacent triangles and c-the reduction factor, necessary to fit the fluid flow intensity to real complex fluid flow conditions in materials.
The reduction factor c is determined in parametric studies to keep the maximum Reynolds number R e along the main flow path always lower than the critical one for laminar flow.

Mass flow rate in channels
Three flow regimes in the VPN channels are distinguished: (a) single gas-phase flow with gas-phase fraction a p = 1, (b) single liquid-phase flow with liquid phase fraction a q-= 1 and (c) two-phase flow (liquid and gas) with 0 \ a q-\ 1. For single-phase flow in channels (flow regime 'a' and 'b'), the fluid moves in channels through a thin film region separated by two closely spaced parallel plates according to a classical lubrication theory [35], based on the Poiseuille flow law [5]. As a result, the mass flow rate of the single-phase flow along channels is where (P i and P j are the pressures in the adjacent VPs). A two-phase flow of two immiscible and incompressible fluids in a channel (Fig. 4) is assumed to simulate a twophase fluid flow (flow regime 'c'), driven by a pressure gradient in adjacent VPs. The liquid-gas interface is parallel to the channel plates and constant along the channel.
Gravity forces are neglected. The interface between the fluids, labelled as j = q, p (q-the lower liquid phase and pthe upper gas phase), is assumed to be flat in the undisturbed flow state. Under this assumption, the model allows for a plane-parallel solution. The interface position is known and is related to fractions of fluid phases in adjacent VPs while the volumetric flow rates of fluid phases are unknown.
Continuity and momentum equations characterize the flow in each phase. The time and pressure are scaled by h p =u i and q p u 2 i (h p -the height of the upper layer and u i -the interfacial velocity). The dimensionless continuity and momentum equations are [3] where u j ¼ u j ; v j À Á and p j are the velocity and pressure of the fluid phase j, q j and l j are the corresponding density and dynamic viscosity. The Reynolds number is Re p ¼ q p u i h p =l p and the density and viscosity ratios are r ¼ q q =q p and m ¼ l q =l p . In the dimensionless formulation, the lower and upper phases occupy the regions Àn d y 0 and 0 y 1, where n d ¼ h q =h p . At the channel walls, the velocities meet the no-slip boundary condition The boundary conditions at the interface y ¼ 0 require the continuity of velocity components and tangential stresses [3] and l q ou q oy The solution details are presented in [25]. Solving Eqs. 11 and 12 with boundary conditions (Eqs. [13][14][15], the mass flow rates M q;x and M p;x for both fluid phases are calculated (as well as the shear stress s f 0 at the channel surfaces for y ¼ Àn d and y ¼ 1).

Fluid flow in virtual pores
VPs, in contrast to the channel flow model, assume that the fluid is compressible. The fluid pressure can reach 70 MPa in specific situations, such as during the hydraulic fracturing process. The gas phase exceeds the critical point and becomes a supercritical fluid under these conditions. For both fluid phases in VPs, the Peng-Robinson equation of state [34] is used to describe fluid behaviour above the critical point where P is the pressure [Pa], R denotes the gas constant (R = 8314, 4598 J/(kmol K)), V q=p is the molar volume of liquid (q) and gas (p) fraction [m 3 /kmol] and T denotes the temperature [K]. The parameters in Eq. 16 are: where T q=p;c is the critical temperature of phase [K], P q=p;c denotes the critical pressure of phase [Pa], x q=p is the acentric factor of phase [-] and T r denotes the reduced temperature T T c . When c 1 = c 2 = c 3 = 0, the original model is obtained. The extra factors help connect vapour pressure data from highly polar liquids like water and methanol. For most substances, Eqs. 17-22 provide a good fit for the vapour pressure, however predicting molar volumes can be very inaccurate. The forecast of saturated liquid molar quantities might deviate by l0-40% [27]. Peneloux and Rauzy [33] proposed an effective correction term where s is the small molar volume correction term that is component dependent; V q is the molar volume predicted by Eq. 16 and V corr q refers to the corrected molar volume. The value of s is negative for higher molecular weight nonpolar and essentially for all polar substances. The molar volume correction term is considered to be 0.0 m 3 /kmol and -0.0034 m 3 /kmol for the gas phase and liquid phase (water), respectively. The Peng-Robinson equation of state has the advantage of being able to describe the behaviour of supercritical fluids at extremely high fluid pressures and temperatures. For each phase, the mass conservation equation is used. The mass transfer between phases and the grid velocity is ignored when there is no internal mass source. The discretized form of the mass conservation equation for the liquid phase is where f is the face (edge) number, U n f denotes the volume flux through the face [m 3 /s], based on the average velocity in the channel, a n q;f is the face value of the fluid-phase volume fraction [-], t is the time step [s], n denotes the time increment and i is the VP number [-]. The explicit formulation is used instead of an iterative solution of the transport equation during each time step since the volume fraction at the current time step is directly computed from known quantities at the previous time step. Similarly, the mass conservation equation for a gas phase is introduced. The product q q U n f a n q;f in Eq. 24 is the mass flow rate M q;f of the liquid phase flowing through the face f (edge of a triangle) of VP i . The density of the liquid phase can be calculated by solving the mass conservation equation for both phases The density of the gas phase can also be computed in the same way. It should be noted that the molar volume V (q/p) is related to the gas density.
and to the liquid density Due to the fact that the fluid phases share the same pressure the fluid-phase fractions are computed. Inserting Eq. 27 for the gas phase and Eq. 28 for the liquid phase into Eq. 29, a polynomial equation is obtained with respect to the liquid fraction a nþ1 q;i . The gas-phase fraction is computed as a nþ1

Heat transfer
Heat must be transferred in both the fluid and solid domains. When it comes to heat transfer in multi-phase fluid flow, the temperature is shared, but enthalpy is transferred. The heat transfer model is simplified in the same way that the fluid flow model is simplified. As a result, a homogeneous heat transfer model is assumed in multi-phase flow (mass transfer between phases is not supposed to be taken into account). The multi-phase fluid is homogenized to a single-phase fluid. The effective fluid properties and velocity are calculated using volume averaging over the phases. The numerical solution uses the same very coarse mesh of both domains (Fig. 1b) that is used to create the fluid flow network to solve the governing equations.

Heat transfer in fluid
In multi-phase fluid flow, we assume homogenous heat transfer. The fluid is incompressible and homogeneous. The viscous dissipation of energy is not taken into account.

The energy conservation equation is shared by all phases in the homogeneous model and is expressed in integral form
The specific heat capacity is assumed to be independent of composition and pressure: Equation 30 is applied to each fluid cell (triangle) in the computational domain. The finite volume method is used to solve Eq. 30. The discretization of Eq. 30 on a given cell produces in 2D, rT f denotes the gradient of T at the face f and V is the cell volume. If the time derivative is discretized using backward differences, the first-order accurate temporal discretization is given by where n ? 1 is the value at the next time step t ? Dt and n is the value in the current time t. The energy conservation equation can be expressed in terms of temperature T by assuming that total energy E equals enthalpy h and applying the enthalpy equation of state to Eq. 38 where S h is related to the internal enthalpy source of the diffusive energy [W/m 3 ] of heat transported by diffusion along the channel ''S2S'' where L k is the length [m] of the channel S2S and A k is the area of the channel cross section [m 2 ]. The total energy can be calculated using the enthalpy equation of state E ¼ c p T À T ref ð Þsince total energy equals enthalpy. The Green-Gauss theorem is utilized to generate a scalar value at cell faces and compute secondary diffusion terms and velocity derivatives. To obtain the gradient of the scalar T at the cell centre c0 (Fig. 5), the following discrete form is stated as where T f represents the value of T at the cell face centroid.
The face value T f in Eq. 41 is derived from the arithmetic average of values at neighbouring cell centres using Green-Gauss Cell-based gradient assessment The cell grid is assumed to be very coarse, therefore when the surface f is at the domain boundary or at the liquid-solid interface, the temperature T f is analytically calculated. Several scenarios may be considered: case I: the temperature T f ;bc t ð Þ on the boundary is known (constant or  , where h e is the convective heat transfer coefficient in the surroundings (T e is the temperature of surroundings), case IV: if the fluid cell face is adjacent to the solid cell, the interface is assumed as transfer coefficient h f is calculated using the Nusselt number Nu k , which may be empirically obtained using the macroscopic porosity of the particle assembly and Reynolds number (0 \ Re \ 102) (according to Tavassoli et al. [45]) After that, the convective heat transfer coefficient h ik is calculated The diffusive flux D f in Eq. 39 across the face f is given by The parameter D f may be approximated as where the first term on the right-hand side represents the primary gradient oriented along the vector e s ! and the second term represents the 'cross' diffusion term. In Eq. 47, the parameter A ! f is the area normal vector of the face f directed from the cell c0 to c1, ds is the distance between cell centroids and e s ! is the unit normal vector in this direction (Fig. 6) and rT is the average of gradients at two adjacent cells. Discrete values of the scalar are stored at the cell centres. However, face values are required for convection terms in Eq. 39 and must be interpolated from cell centre values. This is accomplished through the use of an upwind scheme. As a result, the second-order upwind scheme is used to compute face values. A multidimensional linear reconstruction method is used to calculate quantities at cell faces. Thus, the face value u f is computed using the expression below where u-the cell-centred value in the upstream cell, ruthe cell-centred gradient of the value in the upstream cell, r ! -the displacement vector from the upstream cell centroid to the face centroid (Fig. 5).
The gradient in each cell must be determined for this formula to work. Finally, the gradient must be limited to prevent the introduction of additional maxima or minima. As a result, the concept of gradient limiters [4] is introduced.

Heat transfer in solid
The energy conservation equation has the following integral form if there is no convective energy transfer, no internal heat sources, and constant density in solid regions where E is the total energy, equal to enthalpy h ¼ R T T ref c p dT, q s denotes the solid density [kg/m 3 ], k s is the thermal conductivity of solid [W/(mK)], T ref denotes the reference temperature and c p is the specific heat in constant pressure. Equation 49 is applied to each cell (triangle) in the solid domain. On a given cell, the discretization of Eq. 49 produces The gradients and face values are computed in the same way that the fluid gradients and face values are computed (Sect. 2.3). The total energy is calculated using the enthalpy equation of the state

Coupling scheme
The discretization algorithm is based on the Alfa Shapes theory [6]. The triangular pore mesh deforms extensively and even affects the geometry topology in specific problems (e.g. in hydraulic fracturing when the fracture begins to propagate). Equation 25's volume changes are transmitted to CFD. When the topological features of the grid geometry change, the grid is automatically re-meshed [24]. Even though the mesh is coarse, the re-mesh procedure takes a relatively long computational time. Therefore, an algorithm that tracks configuration changes of discrete elements is developed and implemented in the model. The algorithm tracks displacements, overlaps, and size changes of discrete elements in each iteration (time increment). The remeshing procedure runs when the configuration changes are large enough. The criterion for triggering the re-mesh procedure is arbitrarily defined but is based on the results of a parametric study. As a result, the repeated remeshing procedure is run on average every several hundred or even several thousand iterations and slightly increases the total computation time. However, the total number of times the remeshing procedure is repeated strongly depends on the problem studied. The calculation outputs (e.g. pressures, fluid fractions) are accurately converted from the old mesh to the new mesh in this scenario, providing that the mass is a topological invariant [24]. Equation 29 is used to determine the new fluid-phase fractions by using mass quantities rather than mass flow rates. In the VPs, Eq. 16 is utilized to calculate the new pressure after transformation.
The coupling scheme of DEM with CFD and heat transfer (Fig. 7) involves three sets of discrete equations to be solved: the fluid flow equations for all VPs in the fluid domain, the law of motion in DEM for all discrete elements (e.g. spheres) and heat transfer equations for all grid cells in the fluid and particle (solid) domain. The two-way coupling scheme is based on a transfer of pressures, shear stress forces, and temperatures from CFD to DEM and the time derivative of VP volumes from DEM to CFD. The pressure and shear forces from the fluid cause displacements of discrete elements (e.g. spheres) in DEM that change coordinates and triangle volumes in VPN. Temperature changes the diameters of spheres and, as a result, the volume of the solid and fluid. The VP volume's time derivative, dV/dt, is computed in DEM and then transmitted to CFD (Eq. 25, which includes the term dV/dt, takes into account the volume change). As a result, pressure and temperature changes in the fluid are affected by volume changes in DEM. The forces acting on spheres are then transformed into fluid pressures in VPs and channels. CFD calculates the pressures and shear stresses in the 'S2S' channels as well as the temperature in each time step Dt C . They are transferred from CFD to DEM in each time step Dt D . They are next utilized to compute fluid forces, which are added to the contact forces F before time integration to update discrete element displacements. The fluid pressures in VPs and 'S2S' channels are converted into the forces F ! P;j acting on spheres where n ! is the unit vector normal to the discretized sphere's edge, P i -the pressure in VP, i-the VP number, jthe sphere number and A k -the contact area between the fluid in VP i and sphere where r j is the maximum sphere radius and e k -the sphere edge length. The shear stresses are converted into the forces acting on spheres in 'S2S' channels as where I ! denotes the unit vector parallel to the channel wall and oriented in the fluid flow direction, s f 0;i -the shear stress in the channel, i-the channel number, j-the sphere number and A k -the contact area between the channel fluid and sphere and L k -the channel length. The new sphere radii are computed through a thermal (linear) expansion coefficient where r nþ1 -the new radius [m], r n -the previous value of radius [m], a exp -the thermal (linear) expansion coefficient [1/K], T nþ1 avg -the new volume average sphere temperature [K], T n avg -the previous volume average sphere temperature [K] and n-the time increment [-].
The DEM-based THM model was implemented by the authors into the open-source software package YADE [41] which is already parallelized by computer cluster nodes using the OpenMPI library. Technically, the DEM-based THM model (called VPNengine) is a native YADE plugin. However, VPNengine is parallelized with the use of the OpenMP library (by single computer threads) and only a few loops are parallelized. Overall, VPNengine only uses one CPU core, while several loops can use a certain number of threads. Therefore, the plugin's performance is limited. The calculation time of bar specimen thermal contraction during cooling (Sect. 6, case 'a') is about 1 day using a processor Intel Xeon 2.7 GHz. However, in simulations of short-term processes, e.g. hydraulic fracturing, a specimen of about 1100 spheres is processed in a very similar computer runtime. To overcome this limitation, the authors began developing a 3D DEM-based THM model to Fig. 9 Pure DEM simulations of uniaxial tension for cohesive granular specimen of Fig. 8: a deformed specimen (sphere displacements were 10-times enlarged), b distribution of normal sphere contacts in non-deformed specimen at failure (broken contacts are in colour) and c particle displacement vectors before peak and at failure in non-deformed specimen (displacement vectors were 50-times enlarged) be parallelized by the cluster computer nodes using the OpenMPI library. This will greatly increase the efficiency of the software and enable 2D and 3D simulations to be performed on samples containing thousands of spheres within a reasonable computer runtime.

Calibration of DEM and CFD
DEM and CFD were separately calibrated, based on mechanical tests (DEM) and a permeability test [24] (CFD), respectively. For calibration, a simple bonded granular specimen in the form of a bar was used (Fig. 8). A random distribution of spheres was employed with a 20 9 80 mm 2 bar specimen approximating artificial porous material with a porosity of roughly p = 10%. The pores represented the empty area between the spheres. The sphere diameter ranged from 4 to 8 mm with a mean diameter of d 50 % 6 mm. Along the depth of the specimen, one layer of spheres was applied.

DEM results
Pure DEM was used during uniaxial tension of a bonded granular bar specimen of Fig. 8. To maintain quasi-static conditions during uniaxial tension, the smooth bottom and top of the specimen (Fig. 8) were moved at a constant opposing vertical velocity of 1 mm/s. Based on early calculations, DEM simulations used a time step of 2 9 10 -7 s. In DEM simulations with spheres, the following material constants were assumed: E c = 3.36 GPa and t c = 0.30, C = 150 MPa and T = 30 MPa (C/T = 5), l c = 18°and q = 2 600 kg/m 3 , based on quasi-static uniaxial compression and splitting tests on shale specimens [24]. The DEM results were in agreement with the corresponding experimental results [15]. The greatest compressive normal stress was 47 MPa, while the maximum tensile stress was 8 MPa. Figure 9 depicts the deformed specimen, broken normal contact distribution at the failure, and sphere displacement vectors (before the stress peak and at the failure) during uniaxial tension. The red and blue colours in Fig. 9c denote the displacements [ 0.08 mm and \ -0.08 mm, and the orange and green colours the displacements between 0 and 0.08 mm, and -0.08 mm to 0 mm.
A nearly horizontal macro-crack appeared at the top specimen region due to the neighbourhood of the top boundary (Fig. 9a). In the major macro-crack, six normal contacts were broken, for a total of twelve contacts broken in the specimen (Fig. 9b). The sequence of damage was the following: first red contacts, then orange and finally green ones. Due to boundary conditions assumed, the sphere displacements were mostly vertical (Fig. 9c). The mean The behaviour of the specimen was elastic-brittle due to a low number of spheres and 2D conditions [28]. The specimen's tensile strength was 8.8 MPa for the vertical displacement/strain of 0.22 mm/0.28% and the elastic modulus was 3.5 GPa [15,24]. Before the peak, the mean tensile normal contact force between spheres was 20.58 N with a maximum tensile contact force equal to 69.20 N (= TR 2 , Eq. 6). The coordination number was 6.53 (initially) and 4.17 (at the failure).

CFD results
A basic Darcy permeability test was performed on a bonded granular specimen of Fig. 8 to calibrate the fluid flow model regarding the macroscopic permeability. Simulated was single-phase water flow (the specimen was fully saturated). The top edge was subjected to a constant water pressure of 0.1 MPa, while the bottom edge was subjected to a constant water pressure of 5 MPa. On the side edges, zero mass flow rate requirements were imposed (Fig. 10a). The macroscopic permeability coefficient was determined using Darcy's law, assuming that the volumetric flow rate at horizontal edges was the same at equilibrium where Q is the volumetric flow rate at the equilibrium state [m 3 /s], A is the specimen cross section [m 2 ], L denotes the specimen height and DP is the pressure difference between the bottom and top edges [Pa]. The simulation was carried out until the mass flow rate of the fluid along the lower and the upper boundary stopped changing. The fluid pressure isolines were almost parallel to each other, which confirmed that the fluid flow at the macro-level was one-dimensional (Fig. 10b). The volumetric flow rate at the equilibrium state reached 3.5.31 9 10 -5 m 3 /s. Hence, the cohesive granular bar specimen had an estimated permeability coefficient of j ¼ 2:47 Â 10 À16 m 2 . The basic mechanical, fluid and heat material constants for the solid, water and gas are summarized in Table 1. 4 Validation of THM model for 1D heat transfer problem The THM model was validated by comparing the numerical findings with the analytical solution for the classic 1D heat transfer problem (diffusion) in the cohesive granular bar specimen where a eqv is the effective thermal diffusivity [m 2 /s] and t is the time [s]. The initial and boundary conditions for the analytical solution of the 1D heat equation are as follows where L is the length of the bar. Using the Fourier series, the unsteady solution becomes where The calculations were performed using a bonded granular bar specimen with a random distribution of spheres (Fig. 8). The effective thermal diffusivity a eqv was calculated for the volume-averaged phase properties where Hence, the density depends on temperature and pressure. In this study, it was assumed for the sake of simplicity that other fluid properties were independent of temperature since most of the simulations presented were performed for a maximum temperature difference of 30 K. In that temperature range, the dynamic viscosity of the fluid slightly changes only. Therefore, we neglected the temperature dependence of viscosity in the current paper and the thermal properties of the fluid.
From the phase properties (Table 1), the effective material properties of the equivalent solid were determined using Eq. 63: k eff = 3.357 W/m K, c p,eff = 929.51 J/kg K, and q eff = 2422.74 kh/m 3 . For four different time values (100 s, 200 s, 300 s, and 400 s) and the corresponding Fourier numbers F o , the comparison results in Fig. 12 are shown.
The numerical results agree with the analytical solution. After 400 s of cooling, the largest difference between numerical and analytical values was 0.52 K. Figure 13   shows the temperature distribution and liquid density in the bar after 400 s of cooling. The density ranged from 1000.02 to 1014.28 kg/m 3 after 400 s of cooling. In the estimation of the water density, the Peng-Robinson equation of state (Eq. 16) with a correction (Eq. 23) generated

Effect of advection on cooling of bar specimen
A random distribution of spheres (Fig. 8) was used to evaluate the effect of advection on the cooling of the cohesive granular bar specimen. Two simulation types were carried out: for a low Peclet (Pe) number (case 'a') and a high Peclet Pe number (case 'b'). The assumption was that water flowed in a single phase. Heat transfer in the bar specimen was simulated by diffusion (due to temperature differences) and advection (due to fluid mass movement) using the initial and boundary conditions in Fig. 14.
The bar specimen was cooled down again by 30 K. The two various pressure differences between the two edges of the bar specimen were defined: 0.9 MPa (low Peclet number) and 9.9 MPa (high Peclet number). Table 1  In the case 'a' after t = 400 s of cooling, the temperature difference between diffusion without advection and diffusion with advection reached 1.14 K (Fig. 15). The maximum Peclet number in the fluid domain was 24. The advection helped to speed up the cooling process. The fluid velocity was small; it did not exceed a velocity of 8.92Á10 -5 m/s. The velocity vectors were almost parallel to the vertical sides (Fig. 16c), indicating that the fluid flow in the bar specimen was 1D. From the bottom to the top, the fluid pressure varied approximately linearly along the bar   (Fig. 16a). After 400 s of cooling, the fluid density ranged from 995.7 to 1012.5 kg/m 3 (Fig. 16b). In the density estimation, the Peng-Robinson equation of state (Eq. 16) with a correction (Eq. 23) produced an insignificant error (less than 1.3%) (the fluid density was again slightly overestimated). In comparison with diffusion-only heat transfer, the temperature change along the vertical centreline was not symmetric. As a result of a colder fluid flowing (advection) in the same direction as the temperature shift, a very minor spatial shift of the temperature (e.g. temperature about 298 K or 303 K) was seen in the fluid flow direction (Fig. 16b). The temperature of the bar fluctuated very little along its horizontal crosssection. The largest temperature difference in the bar's middle cross-section was only 0.008 K.
The comparison of the results between cases 'a' and 'b' with different Peclet numbers is presented after t = 300 s of cooling of the bar specimen (Fig. 17). An increase in the pressure difference between the two edges of the bar specimen significantly increased the Peclet number from 24 to 179. The difference in the maximum temperature in the specimen reached 2 K between Pe = 24 and Pe = 179. A higher mass flow rate for Pe = 179 resulted in a significant right shift of the temperature distribution along the vertical centreline of the bar.

Bar specimen thermal contraction during cooling
A thermal contraction test during cooling of the cohesive granular bar specimen with a random distribution of spheres was performed (Fig. 8). To eliminate advection at the start of the cooling process and keep the fluid beyond the phase change conditions, the boundary and initial conditions of Fig. 18 were used. The flow of two-phase fluid was studied. The bar specimen was composed of 80%  water and 20% air. Two cases were investigated: (a) sphere displacements and rotations were fixed only at both ends of the specimen and (b) sphere displacements and rotations were fixed in all spheres of the specimen. As a result, temperature variations had an impact on mechanical properties. In contrast to earlier simulations, the thermal properties of the solid and fluid were changed to accelerate the heat transfer process and its effect on the damage mechanism in the specimen. As a result, the water and air heat transfer coefficients indicated in Table 1 were multiplied by 10, and the solid heat transfer coefficient of 420 W/(m K) (typical for silver) was adopted. Solids were assumed to have a thermal (linear) expansion coefficient of 0.00083 1/K. Table 1 summarizes the material properties of water, air, and solids. The bar specimen's initial temperature (solid and fluid) was rather high (368.16 K), although it was still below the boiling point. The bar specimen was then cooled by setting a constant temperature of 278.16 K at both ends (bottom and top) (the temperature difference was 90 K). Figure 19 depicts the temperature distribution in the bar specimen after 17 s. A deformed bar specimen with macrocrack, broken normal sphere contacts, and sphere displacement vectors after 17 s of cooling is illustrated in Fig. 19. The distribution of grain diameter changes during cooling is presented in Fig. 21. Figures 22 and 23 show the distribution of the fluid pressure and fluid temperature in the macro-crack zone.

Case 'a'
The temperature was the highest (333.44 K) in the specimen near the macro-crack edge (Fig. 19a). A macrocrack occurred in the bar specimen due to tension (Fig. 20). It was created at the bar specimen's weakest zone near the top edge in a slightly different location than during pure uniaxial tension (Fig. 9a). The greatest normal contact tensile force before breakage was 29.313 N. Five normal contacts were broken in the macro-crack zone (34 in the entire bar specimen) (Fig. 20b). There was a noticeable vertical movement of spheres in the specimen mid-region (Fig. 20c). The mean diameter of spheres d 50 reduced from 6.032 to 5.837 mm (by around 3%) with the smallest changes in the bar's mid-region (Fig. 21). In the damaged location, the fluid pressure dropped to 0.0116 MPa (by 0.0884 MPa below the original pressure) (Fig. 22). Once the normal contacts broke and the macro-crack began to expand, the pressure drop was produced by an increase in fluid volume in the damaged area due to tensile strain.
The increasing width of the macro-crack generated temperature variations (Fig. 23). Between the particles on both sides of the macro-crack, the temperature in the macro-crack area was reduced by as much as 20.9 K. The particle mobility and deformation had, thus, a significant impact on fluid pressure and temperature. Heat transfer in the bar specimen caused fluctuations in fluid pressure, temperature, and density, which resulted in fluid flow, even though the pressure differential between the bottom and upper specimen borders was initially 0 kPa. The fluid's velocity, on the other hand, was still very low, measuring only 0.0057 m/s.

Case 'b'
An additional simulation was performed to show a decoupling effect between the deformation of the cohesive granular bar specimen and fluid flow and heat transport. In the DEM-CFD simulations, the same specimen was used ( Fig. 8) with the same initial and boundary conditions (Fig. 18). A decoupling effect was obtained by an assumption that the granular assembly was unmovable (the displacements and rotations of all spheres were fixed). Thus, the bar cooling did not cause a crack. The temperature distribution in the specimen (Fig. 24) was similar to the simulation results in Sects. 4 and 5, and completely different from case 'a' (Fig. 19).
Due to temperature differences along the specimen height, the fluid began to flow. This caused fluid pressure changes (Fig. 25). After 17 s of cooling, the maximum fluid pressure dropped in the specimen to about 73.0 kPa (below the initial conditions of 100 kPa). The maximum specimen temperature decreased to 328.1 K (case 'a'), while in case 'b', the maximum temperature reduced to 333.3 K. The temperature difference was due to fluid volume changes in case 'a', while in case 'b', the DEM model was decoupled from the fluid flow and heat transfer model and thus the changes of sphere diameters did not occur, and also the changes of fluid volumes.

Summary and conclusions
The research proposes a new DEM-based pore-scale thermal-hydro-mechanical model of two-phase fluid flow in unsaturated porous materials with low porosity that was coupled with heat transfer. The numerical findings were compared to the analytical solution of the 1D heat transfer problem in an analogous bar specimen to validate the model. The effect of advection was also investigated. In addition, the thermal contraction of the bar specimen during cooling was carefully studied, associated with the creation of a macro-crack. Based on meso-scale simulations, the following conclusions can be offered: • During the 1D heat transfer problem, the largest temperature difference between numerical and analytical findings was just 0.52-0.64 K. • Advection increased the cooling of the specimen. The highest temperature difference between cooling by diffusion and cooling by diffusion with advection was 2.26 K after 400 s of cooling with a pressure decrease of 0.5 MPa at the specimen height. The increase in the pressure difference between the two edges of the specimen increased the Peclet number from 24 to 179 and the maximum temperature difference from 1 to 2 K and caused a significant right shift of the temperature distribution along the vertical centreline of the specimen. • In the thermal contraction of the bar specimen after cooling, the tensile failure mechanism was comparable to that of purely mechanical uniaxial tension. • During specimen thermal contraction, the fluid pressure in the macro-crack dropped to 0.0116 MPa, which was 88.4% lower than the initial pressure. • During specimen thermal contraction, the propagating macro-crack created temperature differences. The temperature difference between particles on both sides of the macro-crack was 20.9 K. • Fluid pressure was affected more by particle displacement changes than by temperature changes during bar thermal contraction. • Heat transfer in the cooled bar specimen during bar thermal contraction resulted in pressure and fluid density discrepancies, causing the fluid to flow at a very low velocity of less than 0.0057 m/s.

Declarations
Conflict of interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.