Abstract
An Enhanced Geothermal System (EGS) is an artificially created geothermal reservoir formed by hydrofracturing hot dry rock. Thermal shock occurs when the cold water contacts the hot rock near the injection borehole, creating a network of small, disorganized, closely spaced micro cracks. As the cold-water injection continues, the hot rock cools down and the micro cracks coalesce, becoming a better-defined network of thermal fractures. Thermal fractures in an EGS reservoir are believed to improve reservoir performance by increasing the surface area for heat exchange and lowering flow impedance; however, it is difficult to precisely predict how they grow and affect the permeability of the reservoir. The goal of this paper is to provide an insight into the transport mechanisms within the thin, permeable, thermally shocked region of an EGS reservoir. COMSOL Multiphysics® is used to set up an indealized porous region with identical geometrical features at different domain scales to show the scale dependence of heat and mass transport in the initial microscale crack network and in the later coalesced thermal fractures. This research shows the importance of EGS maturity in determining how heat and mass are transferred and how to select appropriate analytical tools for different stages of development.
1. Introduction
Global demand for electricity generation from alternative energy sources is increasing. More countries are evaluating the potential of geothermal energy as their alternative energy source for electricity generation [1]. An Enhanced Geothermal System (EGS) has the potential to take geothermal energy production to a new level of utility-scale energy production. For the last several decades, starting with the Fenton Hill project, several EGS projects have been developed with the hope of understanding the complex nature of man-made geothermal reservoirs. This type of man-made geothermal reservoir requires a cold-water injection to a hot but dry granitic rock at several kilometers deep to create an artificial reservoir and recover the injected water as hot steam through production wells. Hot dry granitic rocks are easier to find at that depth. One of the main engineering problems has been the inability to connect the injection well with production well that are several hundred meters apart. This is partially due to the lack of our understanding of how the state of stresses from the weight of the overlaying strata and locked-in stresses of the tectonic origin combine to interact with the induced thermal changes. The authors assume that the initial development of thermal fractures caused by the cold-water injection controls the patterns of fracture propagation, as this pattern development will initially define the connecting paths between injection and production wells. Thus, it is important to investigate how heat and mass are transferred in the thermally shocked region and how they will influence the development of reservoir permeability.
Thermal fractures are believed to improve reservoir performance [2–4]. In an EGS reservoir, as the cooling process continues, thermal fractures are widened and penetrate deeper into the hot rock. This stimulation increases the reservoir’s ability to transport heat and mass by lowering the flow impedance and by increasing the heat exchange surface areas. Thus far, thermal fractures have not been extensively studied or incorporated into reservoir simulations because they add another layer of complexity to the modeling; however, their inclusion is vital to understanding the thermal and hydrologic effects of cold-water injection.
Several authors have studied how thermal stresses are developed in brittle materials due to cold-water injection. Chen and Marovelli [5] conducted an experiment to analyze thermal stresses in a rock disk subjected to an external thermal shock. Perkins and Gonzalez [6] as well as Kocabas [7] proposed analytical models to investigate the state of stresses induced by cold fluid injection. Ghassemi et al. [8] developed an integral equation to calculate thermally induced stresses associated with the cooling of a planar fracture in a hot rock. The findings from these studies are relevant to the induced tensile stresses present in the cooled region of the EGS reservoir that stimulate thermal fractures. Another key finding from these studies is the direct correlation between the difference in the applied temperature and the crack density.
A series of studies have also been conducted to understand the mechanical and thermal behavior of thermally induced fractures and their interactions in brittle-elastic materials [9–12]. Nemat-Nasser et al. [13] investigated stability of growing thermal fractures. Barr [14] investigated the branching of thermal fractures by examining the potential of crack propagation away from the surface of the primary fracture in a hot dry rock (HDR) system. Bažant and Ohtsubo [15] as well as Murphy [3] and Barr [14] show that, during the heat extraction process, a network of closely spaced thermal fractures, which resembles a “waffle” grid of grooves [2], is formed adjacent to the primary fractures. One of the common observations of these studies is that the growth patterns of thermal fractures are complex, particularly under the large thermal strain close to the surfaces of primary fractures.
Although these studies present a comprehensive view of the impact of thermal fracturing, thermal fractures have not clearly been integrated into EGS reservoir modeling. In fact, a common practice of reservoir modeling neglects the positive effects of thermal fractures by excluding them from the core structure of simulation, and this exclusion hinders progress towards a better understanding of the contribution of thermal fractures in reservoir performance.
Only a small number of studies have been conducted to investigate the effect of thermal fractures on reservoir performance. Harlow and Pracht [16] were probably the first to integrate thermal fractures into a reservoir model as an additional porosity. Stephenes and Voight [17] show that the presence of thermal stresses lowers the pressure requirement to initiate hydraulic fractures in the reservoir. Tran [18] also indicates that the presence of thermal fractures in a reservoir could change both the shape and aperture of the primary hydraulic fractures, and as a consequence, reservoir performance is affected. Huang et al. [12] developed a quasi-static discrete element model to simulate fracture propagation induced by thermal stresses. Huang et al.’s model of thermal fractures has a significant influence on heat conduction within the rock matrix. These studies highlight the importance of thermal fractures in an EGS reservoir model; however, the heat and mass transport mechanisms within the thermally shocked region of an EGS reservoir have not yet been clearly investigated.
Motivated by these observations, this paper investigates the transport mechanisms within a thermally shocked region to provide a better understanding of the overall heat and mass transport in an EGS reservoir. In response to the cold-water injection, a thin thermally shocked region is formed adjacent to the hydraulically induced fracture. The thermally shocked region initially acts as a transition zone between the hydraulically induced fracture and thermal fractures. The structure of the thermally shocked region and its transport properties can affect overall reservoir performance in two ways: first by controlling the growth behavior of the thermal fractures and second by altering heat transfer from the rock mass to the working fluid. In this study, a porous medium consisting of a solid skeleton and fluid pathways is considered as a conceptual model of the thermally shocked region. The pore-scale approach is adopted, and mathematical formulation and numerical examples are developed to investigate the transport mechanisms within the region of interest.
2. Conceptual Model
Figure 1(a) shows a conceptual model of an EGS reservoir comprising a series of parallel hydraulically induced fractures connected to define a doublet system of an injection and a production well embedded in hot dry crystalline rock. This system is commonly used to characterize an idealized EGS reservoir [2, 20, 21], although the geometrical configuration of the model does not fully portray a realistic EGS reservoir.
As mentioned earlier, the occurrence of the thermal fractures in an EGS reservoir has been proven theoretically and experimentally [3, 5, 13–15]. Their presence has also been acknowledged in induced seismicity studies in several operating geothermal sites [22]. However, thermal fractures are excluded from the geometrical representation of the conceptual model of an EGS reservoir. A representation of the thermally shocked region adopted from earlier studies [3, 10] is shown in Figure 1(b).
This paper focuses on the transport behavior of a thin thermally shocked region adjacent to a hydraulically induced lenticular fracture. The goal of simulation of the thermally shocked region is to generate randomized flow paths at different length scales to gain a better insight into how heat and mass are transported. The thermally shocked region with micro cracks is idealized as a granular porous medium. Following this assumption, MATLAB® was used to generate a random 2D porous region composed of randomly placed, nonoverlapping circles. To this end, the authors are less concerned with the specific shape of particles, as long as the nonoverlapping packing of circular particles can create similar tortuous paths that are created by a collection of small fractures. Various arrangements of clustered and differently sized circles approximate irregular shaped matrix elements. The resulting clusters of randomly sized particles create preferential flow paths that can be used to effectively measure the scale effect on mass transport and temperature distribution.
The degree of fragmentation due to thermal shock depends on several factors, such as the severity of the temperature gradient, the depth of reservoir, the temperature mismatch between the rock and the water at the fracture surfaces, the thermal expansion coefficient, Young’s modulus of the rock, and the rate of water circulation [2]. Since we cannot determine the exact density of the fractures within the thermally shocked region, defining a Representative Elementary Volume (REV) with different length scales can be a practical option to characterize the density of the thermal fractures within this region. In this work, we utilize REVs with different side lengths, from 1 mm to 500 mm. The REV with the side length of 1 mm can yield a value representative of the thermally shocked region with the highest density of fractures. Therefore, the REV with the length scale of 1 mm is considered a severely thermally shocked region. Proportionally, a model with the length scale of 100 mm can be considered a moderately thermally shocked region.
The 1 mm scale model represents the initial stage of the thermally shocked region, and the 100 mm represents the latter stage once the cracks have coalesced or the initial region that can be developed when the thermal shock was not severe enough so that fragmentation was coarse. Porous media with length scales lower than 1 mm and larger than 500 mm are beyond the scope of this paper. The former length range may cause nanoscale pore spaces, and the latter may cause an unwanted turbulent flow regime.
3. Contributing Physics
In order to describe the transport mechanisms within the thermally shocked porous region of an EGS reservoir, a set of equations that describe the fluid flow, heat, and mass transport are required. Furthermore, these equations must be coupled to capture the multiphysics nature of the described transport phenomena. In this work, a sequential coupling approach is taken. Initially, fluid flow and mass transport are considered to be the only contributing transport processes without heat transfer. Thus, the connected pore network alone is responsible for fluid and mass transport and they are coupled in an isothermal manner. Later, the solid rock matrix is added in conjunction with the pore network, and the equation of heat transfer is added to the mass transport model to build the sequential coupling for heat and mass transport. This sequentially but fully coupled system is capable of analyzing the nonisothermal transport behavior of a thermally shocked region. The physics of fluid flow, heat, and mass transport in porous media exhibit complex phenomena that stem from different length scales. The transport phenomena in the severely thermally shocked region with many small fractures can be best described by coupling the Navier-Stokes and advection-diffusion equations together with conjugate heat transport. In order to investigate the distinct transport phenomena of different scales, we create computational domains that are geometrically identical but different in scale. In each scale, the small identical Reynolds number assures that the flow is laminar.
4. Pore-Scale Modeling
In this section, we will describe our sequential coupling method to analyze heat and mass transport at the pore-scale.
4.1. Mass Transport (Isothermal Coupling)
Following the sequential coupling approach, the physics of fluid flow and mass transport are first coupled. There is no thermal interaction between the fluid in the pore and the solid rock matrix. The flow is assumed to be isothermal single-phase laminar flow, and the Navier-Stokes equation governs the fluid flow. Figure 2 shows the computational domain. Relative Concentrations of mass and temperature arrivals are monitored at discrete locations (A, B, C, and D) of this simulation domain.where denotes the density of the fluid (kg/m3), the fluid velocity (m/s), the pressure (Pa), a body force term (N/m3), and the identity matrix. The low-pressure gradient of Pascals is applied on the boundaries to keep the flow regime laminar. The initial and boundary conditions can be expressed asA no-slip boundary condition is applied at the interior walls (solid skeleton surfaces), as well as on the upper and the lower walls. The transient advection-diffusion equation governs transport of nonreactive species of mass in a laminar flow regime.where denotes mass concentration (mol/m3), the diffusion coefficient (m2/s), the velocity (m/s), and the source or sink term (mol/m3s).
The complex geometry of the pore network in conjunction with the coupling of the contributing physics makes the simulation computationally expensive. To maintain numerical accuracy with a manageable simulation time, a Gaussian profile is given to express the initial condition of the mass transport. For the initial condition, the concentration of the mass within the domain is defined as follows:According to (5), within the pore-network domain very close to the left boundary (adjacent to the inlet, ), the mass concentration is equal to . However, within the rest of the pore-network domain the mass concentration is almost zero (see Figure 13).
For the boundary conditions, a constant concentration of is applied to the inlet boundary, and the species are transported out of the domain by fluid motion to the right side of the model. where is the normal vector, used to identify the component of the mass flux perpendicular to the boundary. All other boundaries, including solid skeleton surfaces and the upper and lower boundaries, are assumed to be insulating boundaries where no mass flows in or out.For a complex flow field geometry such as the one considered here, the Stokes equations do not have simple analytical solutions; however, the transport of species in one dimension has been analytically solved [23–27] and it can be expressed aswhere is the Relative Concentration (RC) defined as the ratio of the species concentration to its initial value at a distance from the inlet (m), erfc is the complementary error function, is the flow speed (m/s), is the diffusion coefficient (m2/s), and is the time (s).
Using (8), the Mass Breakthrough Time (MBT) is calculated against velocity for a given RC value. This RC specific MBT is calculated at the outlet D, that is, the time to travel 1 mm in distance. Two major trends can be observed in Figure 3. First, when the velocity becomes small enough, that is, less than m/s, then MBT for each RC value becomes constant. This is an indication that mass is transported by diffusion, and the MBT depends only on the concentration gradient. On the other hand, when the velocity becomes large enough, MBT values for each RC value merge to a single value that can be calculated as the time that takes the fluid velocity to travel 1 mm in distance. This is the indication that the mass is now transported along the movement of fluid. More specifically, for example, for RC of 0.1, the MBT due to pure diffusion can be found by reading off the value for the point where the plotted red line intersects with the smallest velocity value. On the other hand, the MBT due to pure advection can be found as the point where the red line intersects with the largest velocity value. Any MBT between these two extremes are due to the combined effect of diffusion and advection. The precise contribution of each effect is unknown. However, it is clear from Figure 3 that, for slower velocities, the MBT is influenced more by the diffusional mass transfer. Here, only the case with the length scale of 1 mm is considered for validation. It is noted that the coefficient of diffusion is assumed to be equal to that of self-diffusion (or tracer-diffusion) of water, m2/s.
In Figure 3, based on our theoretical analysis, we show which transport mechanisms are at work for a given velocity and Relative Concentration. Now, the model is validated by measuring the RC values as the flow passes through the outlet D (see Figure 2). Figure 4 shows that, for a given RC value, the MBT is the shortest when both diffusion and advection participate in transfer mechanisms and the longest when pure advection controls the dispersion. It is also noted that, for the advection dominated flow, the mass passes the outlet D with a large concentration gradient.
4.2. Heat and Mass Coupled Transport (Nonisothermal)
A fully coupled nonisothermal, transient model is needed to capture the transport phenomena within the thermally shocked porous region of the EGS. Unlike the isothermal coupling model where only the fluid phase in the pore network was responsible for mass transport, for the nonisothermal coupled case, both the pore fluid network and solid rock matrix participate in heat and mass transfer. The conjugate heat transport concept [28] is used to couple the physics of heat transfer in both domains.
It is assumed that heat is transferred within the solid matrix only due to conduction; therefore, Fourier’s law can define the temperature field in this domain.where is the effective volumetric heat capacity of the solid phase at constant pressure (J/m3°C) and is the thermal conductivity of the solid (J/ms°C). In the pore fluid network, the energy equation is modified to reflect the influence due to viscous effect and the temperature-dependent pressure work. Accordingly, the following equations are used to govern the energy transport in the system:where indicates the effective volumetric heat capacity of the fluid at constant pressure (J/m3°C), the thermal conductivity of fluid (J/ms°C), and u the fluid velocity. Initial and boundary conditions are applied as follows:A constant temperature of is prescribed at the inlet. It is also assumed that there is no heat flux across the upper and lower boundaries.By accommodating these coupled equations in the numerical model, COMSOL Multiphysics was used to numerically study coupled heat and mass transport phenomena in a complex fractured rock. COMSOL Multiphysics is a Finite Element (FE) simulation software, well suited for simulating mechanical-thermal-hydrological coupled problems. The computational domain of interest is discretized explicitly for both pore-network and solid rock domains. Time is also discretized by utilizing the Backward Differentiation Formula (BDF) time stepping method. The BDF is a family of implicit, linear multistep methods used to numerically integrate ordinary differential equations, and its transient solver is both stable and versatile and provides extra robustness required for transport applications. Parameters used in the simulation are listed in Table 1. For more details on the numerical model and the initial and boundary conditions see Appendix.
5. Results and Discussion
Fully coupled models were simulated at the domain length scales of 1 and 100 mm. A small pressure difference of 0.1 Pa between the inlet and outlet quickly led to a steady-state condition. Due to the Venturi effect, fluid velocity locally increases as the fluid passes through constrictions and evolves to an uneven flow field. The fluid velocity in the preferential paths is higher and rapidly drops to almost zero as the fluid diverts from the preferred paths to locally stagnant regions. The flow patterns of the preferred flow paths for both length scales are identical. Although the flow patterns are identical to both domain length scales because the local velocity gradients define flow paths, the values are different by two orders of magnitude proportional to the scale of the domain, as can be seen in Figure 5. In the following discussion of the results, the mass that is dispersed is interpreted as a tracer.
(a)
(b)
Figure 6 shows the transport of a tracer (as an example of mass transport) within the pore network for both domain length scales. The RC value of 0.5 measured at the center line B is used as a criterion to define this specific MBT for both length scales. For the model with the domain length scale of 1 mm, it takes about 300 seconds for the tracer to reach the RC of 0.5 at B, whereas, for the model with the length scale of 100 mm, it takes about 500 seconds. It is important to note that there is a significant difference in the way the tracer is spatially dispersed. For the 1 mm domain length scale, the tracer transport is distributed uniformly across the domain with a flat propagating concentration front, while, for the 100 mm scale model, unevenly propagating tracer front with preferential paths can clearly be observed. This difference clearly shows that, at the smaller domain scale, the tracer is dispersed by diffusion due to the concentration gradient, and at the larger scale, mass is dispersed by advection. It is reminded that similar flow patterns do not assume similar mass transport. As shown in Table 2, the Reynolds numbers for flows at 1 and 100 mm scales are and 1.40 m/s, respectively. In both scales, the flow is extremely slow and laminar.
(a)
(b)
(c)
(d)
Figure 7 shows how temperature is transported at both domain scales. While the fluid with lower temperature sweeps the rock domain with a higher initial temperature, thermal energy is transported through conduction in the rock matrix and through both conduction and convection in the pore network. The Relative Temperature (RT) is defined as the ratio of the current temperature to its initial value, and the RT value of 0.5 at B is used to define the Temperature Breakthrough Time. Results in Figure 7 indicate that, for the model with the length scale of 1 mm, the average temperature at B achieves the criterion in less than 1 (0.26 sec) second of the cold fluid injection. However, for the model with domain length scales of 100 mm, it takes 471 seconds for the RT to reach the value of 0.5 at B. It is noted that, for the smaller geometry, heat is transported uniformly across the model with a straight temperature propagation front, while, for the larger geometry, an uneven heat front can clearly be observed. This difference can also be explained based on (10). In the case with the length scale of 100 mm, the energy transported by the fluid motion is dominant, whereas, in the case with 1 mm scale length, the dominant heat transport mechanism is conduction.
(a)
(b)
Table 2 summarizes which mechanisms are dominant for heat and mass transport at each domain scale. It is noted that the small Reynolds numbers indicate that the flow remains laminar at all scales, and furthermore for the domain length scale up to a few millimeters, the given pore structure produces a Stokes flow with creeping motion. The Reynolds number defines the relative influence between inertial and viscous resistance on flow behavior. The Péclet number (Pe) was also calculated to define the ratio of the rate of advection to the rate of diffusion of the same physical quantity driven by an appropriate gradient. In the context of mass transfer, the Péclet number is calculated as the product of the Reynolds number and the Schmidt number. In the context of heat transfer, the Péclet number is equivalent to the product of the Reynolds number and the Prandtl number. As can be seen in Figure 8, mass transfer is dominated by diffusion for the length scale of 1 mm but as the length scale increases, mass transfer due to advection quickly becomes more dominant. On the other hand, heat transfer remains diffusion dominated up to the scale of 10 mm, whereas advection and diffusion are equally important for the length scales between 10 and 30 mm, and beyond the length scale of 50 mm, heat transfer will be dominated by advection (Figure 9). In summary, for the flow considered here at different domain length scales, advection has a greater influence on how mass is transported for length scales for a few millimeters or greater. However, the way heat is transported goes through three different combinations of mechanisms, and they are diffusion up to around 10 mm, advection-diffusion between 10 and 30 mm, and advection greater than 50 mm of domain length scales.
5.1. Thermal versus Mass Breakthrough Time
Quantitative criteria to assess EGS reservoir performance are divided into two groups. The criteria are used to judge the hydraulic performance and those that evaluate the thermal performance of the EGS reservoir. Mass Breakthrough Time (MBT) and Thermal Breakthrough Time (TBT) are usually used to identify the distribution of fluid residence times within the reservoir and temperature decline in the recovery as a function of time. The MBT and TBT have been used in the previous sections of this paper to evaluate the thermal and transport properties in the different scale studies.
The scale of the thermally shocked region and its transport properties can affect overall reservoir performance. Comparing MBT and TBT at different length scales can give a better understanding of the transport mechanism within a specific region and its effect on overall EGS performance.
Figure 10 shows the mass and heat breakthrough times at the outlet for the length scales of 1 and 100 mm. Results show that, for the case with the length scale of 1 mm (Figure 10(a)), heat transport is much faster than mass transport. For example, it takes 193 seconds for the mass to reach RC of 0.1 at the outlet; however, for the heat, it takes less than one second to reach a RT of 0.9 at the outlet. For both cases, the time for 10% of the original concentration of mass (tracer) or temperature to reach at the outlet was measured. Three orders of magnitude difference in breakthrough times of heat and mass indicate the critical role of diffusion at the model with the length scale of 1 mm. For the case with the length scale of 100 mm, the mass transport is faster than the heat; however, unlike the smaller scale, heat and mass transport operate on the same order of magnitude (Figure 10(b)). For example, it takes 529 seconds for mass to reach the RC of 0.1 at the outlet of the model. The time for the RT of 0.9 at the outlet is about 844 seconds. Unlike the smaller case, the role of advection is pronounced in the model with the length scale of 100 mm.
(a)
(b)
6. Conclusion and Recommendation
The mechanisms of the heat and mass transport within a thermally shocked region of an EGS reservoir were studied. In this numerical study, the thermally shocked region adjacent to the primary hydraulic fractures was assumed to be a porous medium. A pore-scale simulation approach was considered to study its applicability and limitations in capturing coupled transport physics of heat and mass. For the pore-scale modeling, mass and heat transfer at two domain length scales of 1 and 100 mm were considered to represent a thermally shocked region with two different degrees of fragmentation or fracture coalescence. Using the COMSOL platform, contributing physical laws were sequentially coupled to capture multiphysics features of the problems investigated. The results showed that, for the severely thermally shocked region represented by the domain length scale of 1 mm, diffusion is responsible for both heat and mass transfer. The TBT was found to be three orders of magnitude faster than the MBT. In other words, the ratio of thermal diffusivity to molecular diffusivity was also found to be 103 for this model. This ratio can determine the relationship between heat and mass transfer in a diffusion-dominated small system where the flow speed is insignificant. The pore-scale simulation results also indicate that, for a moderately thermally shocked region represented by the domain length scale of 100 mm, mass transfer is mainly accomplished by advection, but heat is transferred via both fluid motion and conduction. Thermal and Mass Breakthrough Times were found to be in the same orders of magnitude.
Although a simple model has been presented in this paper, the results of pore-scale analyses confirm that the transport behavior of the thermally shocked region of an EGS strongly depends on the degree of fragmentation. In actual field cases, one can identify the following appropriate regions of an EGS reservoir where the presented pore-scale analysis may be beneficial:(i)Near the injection point of water: this is where intense cooling on the surface of the primary hydraulic fracture produces a thin thermally shocked region with a high degree of fragmentation (i.e., higher fracture density).(ii)Near the production point of hot water: this is where the cooling effect becomes less intense, and the adjacent thermally shocked region is less fragmented.(iii)Early stage of EGS fluid circulation: this is when the thermal front advances very rapidly along the primary hydraulic fracture and isotherms in the rock matrix are parallel to the main hydraulic fracture.(iv)Later stage of EGS fluid circulation: this is when the cold front penetrates more into the rock matrix and isotherms in the rock matrix become perpendicular to the main hydraulic fracture.Differences in the transport behavior of the thermally shocked region due to the variation in its fragmentation may cause temporally and spatially sensitive behavior in the reservoir as speculated above. The exact fragment geometry of the inside of an EGS reservoir will never be known, but incorporating representative elements into a reservoir simulation may help more accurately depict and predict the heat and mass transport processes in the EGS circulation system.
Appendix
Numerical Model Description
The computational domain in COMSOL Multiphysics was created in 2D space dimension. In the direct pore-scale approach, the geometry is discretized explicitly for both pore network and solid phases. As shown in Figure 11 the domain is discretized spatially with triangular elements for both pore network and solid skeleton subdomains. Table 3 shows the statistics of the elements. Time is also discretized by utilizing the Backward Differentiation Formula (BDF) time stepping method. The initial time step size is set to 0.001 and as the simulation approaches the convergence, the time steps are increased incrementally to the value of 0.1 s and step size of 0.1 seconds is kept to the end of simulation.
The pore network was assumed to be fully saturated with water and the granite with predefined properties was designated as the solid domain. The properties of the assigned fluid, such as heat capacity, density, and thermal conductivity, were assumed to be temperature dependent (Figure 12) while the properties of the solid phases were considered constant. The values of the basic properties of the solid phases are = 2600 kg/m3, = 850 J/(kg °C), and = 2.9 W/(m °C). An unstructured triangular mesh was created for both solid and pore-network domains. The mesh is adequately sized and is refined at the region adjacent to the inlet where both the species and temperature gradients are higher. The element size of the pore network and solids is calibrated according to the contributing physics, and the average element quality of the mesh for both 1 and 100 mm models is computed higher than 0.94. The porosity of the generated domain is calculated through dividing the integration of the pore-network region by the length and width of the structure:where and denote dimensions of the computational domain in and directions. The porosity results in a value of 0.431.
(a)
(b)
As mentioned, to avoid high numerical dispersion, a gradual profile was given to express the initial condition for both the physics of heat and mass transport (Figure 13). The initial mass concentration everywhere in the domain is zero except in the vicinity of the inlet where the concentration is equal to the inlet concentration. The initial temperature was also assumed to be 85°C everywhere, except adjacent to the inlet where the temperature is assigned to be the same as the inlet temperature. Table 4 describes the initial and boundary conditions for the fluid flow, mass, and heat transport.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.