The standard LBE with the Bhatnagar-Gross-Krook (BGK) estimation can be presented as below Equation 2 [27]:

(2)

In Equation 2, x shows the position in space, t is time, δt is the time step, f_α is the particle distribution function (PDF) along the direction of αth, f_α^eq is its corresponding equilibrium PDF, e_α is the particle velocity in the direction of αth, τ is the single relaxation time and finally N represents the number of discrete particle velocities.

For 2-D flow, the nine-velocity LBE model on the 2-D square lattice can be applied which is, represented as the D2Q9 model. This model has been widely used in the literature. The equilibrium distribution for the D2Q9 model can be identified as below Equation 3 [27]:

(3)

In Equation 3, ω represents the weighting factor, c=δx/δt is the lattice speed, δx is the lattice step and u shows the macroscopic velocity.

Below correlations are the weighting factor and discrete velocity for the D2Q9 model shown in Equation 4 [27]:

(4)

In most LBE simulations, Equation (2) can be solved with two steps. The first step is collision and the second step is streaming. In the collision step, the PDFs for each direction are relaxed toward quasi-equilibrium distributions then the distributions move to the neighboring nodes at the streaming step. The local mass density ρ and the local momentum density ρu are shown as below Equations (5, 6) [27]:

(5)

(6)

The Navier-Stokes (N-S) equations are recovered with LBE model for fluid flow near the incompressible limit and then the viscosity in the N-S equations is obtained as below [27]:

(7)

In Equation 7, c_s shows the lattice sound speed. So applying Equation 7 makes the LBGK scheme a second-order method to solve the incompressible flow [28].

In spite of the reasonable achievements on isothermal flow simulations, the thermal flow simulations have a great potential to be investigated. Generally, the existing thermal LBE models consist of two main categories such as the multispeed and the multi-distribution function approaches. The multispeed approach considers energy conservation by adding the higher-velocity terms in the equilibrium distribution but from theoretical aspect, the multispeed approach is instable numerically [29, 30]. The viscous heat dissipation and compression work with the pressure are considered negligible in the multi-distribution function approach so the temperature field is passively advocated with the fluid flow and then it can be solved with considering an independent temperature distribution function. This term is also called the passive-scalar approach which has better numerical stability and accuracy [31].

Energy equation can be solved in the LBE model framework. In this study, the D2Q9 model has been chosen as the velocity field. Consequently, the discrete Boltzmann equation for temperature field can be updated in Equation 8 [32]:

(8)

where τ_T shows the dimensionless relaxation time of temperature, g_α is the temperature distribution function and finally, g_α^eq is the temperature distribution function at the equilibrium state and it can be calculated with below correlation in Equation (9):

(9)

where u represents whole fluid velocity. The temperature of the whole fluid T is considered with below formula in Equation 10:

(10)

In thermal model, the thermal diffusivity (X) is the key parameter in the simulation of a flow with heat transfer term and it can be expressed as below Equation 11 [32]:

(11)

The Prandtl number (Pr) for component, which related the kinematic viscosity (υ) and the thermal diffusivity (X) can be expressed as below in Equation 12 [32]:

(12)

Microscopically, the reason for fluid system segregation into different phases is related to the inter-particle forces. In the single component multi-phase (SCMP), LBE model has been proposed by Shan and Chen. LBE model is a simple interaction potential model which is defined to describe the fluid/fluid interaction. Then the model incorporates the forcing with shifting the velocity in the equilibrium distribution. The fluid/fluid interaction force is defined as below Equation 13 [31]:

(13)

where x and x^'(=x + e_iδt) show the position in space G(x, x^') and denotes intermolecular forces between particles. The term G(x, x^') ≤ = 0 represents attractive forces between particles. In our study, the interactions of nearest and next-nearest neighbors are considered, therefore, for a D2Q9 lattice model, this leads to below correlation in Equation 14 [31]:

(14)

Phase separation between components occurs automatically whenever the interaction potentials are chosen appropriately. The term ψ(x) is a function of the local density. In the single component model, the function of ψ(x) can be varied and different choices could provide various equations of state. Applying the Chapmen-Enskog method which is a successive approximation method, can obtain the macroscopic fluid equation of the LBE model. Meanwhile, the EOS can be stated as below term in Equation 15 [31]:

(15)

c ₀ is a constant term and it is equal to 3.0 for the considered D2Q9 model. In this study, ψ(x) is considered as below in Equation 16:

(16)

which provides a non-monotonic pressure-density relationship. The body forces can be stated as below in Equation 17:

(17)

where a is the body force acceleration. The influence of the forces has been incorporated into the LBE model through the momentum variation incorporation in the dynamics of the distribution functions. That means that the velocity ‘u’ in Equation 3 can be replaced with the below in Equation 18 [32]:

(18)

where F_total = F1 + F2. Then, with consideration of averaging of the moment before and after the collision, the whole fluid velocity U can be written as below in Equation 19 [33]:

(19)

The geometry of the present problem is shown in Fig. (1). It was assumed that fracture domain consists of a two-dimensional enclosure of height H and width W. For the present case, Aspect-Ratio (AR) is defined as the ratio of the width of the enclosure to the height of the enclosure (AR=W/H). The temperatures of the top and the bottom horizontal walls of the domain are maintained at T_c and T_H respectively. The two side walls have been considered to be adiabatic i.e., non-conducting. The fracture domain is filled with a Newtonian, incompressible fluid and the flow is laminar.

Fig. (1). Geometry of the present study [2].

The Boussinesq approximation is used in the natural convection simulation. With this approximation, all material properties of flows are assumed to be constant except the density in the term of gravity. After absorbing the constant part of the gravity into the pressure, the effective external force can be written as Equation 20 [33]:

(20)

where g_y is the acceleration due to gravity, β is the thermal expansion coefficient, T ₀ is the reference temperature, which is set as T_c (temperatures of the low temperature wall) in this simulation. ρ ₀ is the density of the fluid at T₀. Therefore, the macroscopic dimensionless counterpart of the meso-scale governing equations for the natural convection with Boussinesq approximation for incompressible flow is as follows [5, 26]:

Continuity Equation:

(21)

Momentum conservation:

(22)

Energy conservation:

(23)

Equations (21-23) were obtained using the characteristic length H, velocity scale time scale t ₀ = H/u ₀ and pressure scale . The dimensionless temperature T* is defined in terms of the wall temperature difference and a reference temperature in Equation 24:

(24)

Nusselt number, Nu is one of the most important dimensionless parameters in describing the natural convection heat transfer. The local Nusselt numbers at the hot and cold walls are calculated as by Equation 25:

(25)

Note that the velocity scale in natural convection is proportional to u _o. In order to ensure the problem is in the incompressible regime, this value is less than about 0.1. Therefore, appropriate gβ values must be selected [32].

In practical applications, boundary conditions play important roles in lattice Boltzmann methods in that they will influence the accuracy and stability of the LBE method.

For flow equation, bounce-back boundary conditions were applied on all solid boundaries, which means that incoming boundary populations equal to out-going populations after the collision. For instance, for south boundary by using D2Q9 (Fig. 2) the unknowns are f2, f5 and f6, which are evaluated by using Equation 26 [27]:

(26)

Fig. (2). Domain boundaries and direction of streaming velocities for D2Q9 model [34].

where n is the lattice local on the boundary. For temperature equation, bounce-back boundary condition (adiabatic) is used on the east and west of the boundaries. For instance, at the east boundary, the following conditions are imposed (Fig. 2) as shown in Equation 27:

(27)

The temperature at the north and west walls are known. Since we used D2Q9, in the south wall (T_H), the unknowns are g2, g5, and g6, which are evaluated by using Equation 28:

(28)

Table 1. Comparison of results between CFD simulation [26] and LBE model.

Ra	H		umax		vmax
	CFD	LBM	CFD	LBM	CFD	LBM
103	1.0004	1.0011	3.6464×10-6	2.5534×10-6	3.9684×10-6	4.0012×10-6
104	2.1581	2.1646	0.25228	0.25531	0.26369	0.26687
105	3.9103	4.0013	0.34434	0.341231	0.37569	0.37901
106	6.3092	6.3901	0.37088	0.36977	0.40600	0.41004

Grid independence test was carried out for a case of AR = 1 in Re = 105 and Pr =0.71. The results for the surface average Nusselt number ( ) at the hot wall in different grid numbers (from 40 to 120 grids) on the lower side of the block surface are plotted in Fig. (3). It is mentioned that the reported number of grids is the number of nodes is defined on each side. As it is observed from the figure, the average Nusselt number values change slightly with the further increase of the nodes from 80 to 120 (less than 1%). The test results (Table 2) reveal the good accuracy of all resolutions greater than 80 nodes per square side. However, to ensure the desired accuracy, the 100-grid case is selected for this study.

Table 2. Variation of average Nusselt with grids and percent of its variation.

Grids	40	60	80	100	120
	3.75	3.90	3.98	4.00	4.01
Percent of variation	-	4%	2%	0.50	0.25

Fig. (3). Evaluation of average Nusselt number ( ) on hot wall in different grid numbers.

Fig. (4). Comparison of variations of dimensionless horizontal velocity between LBE model and CFD results at the middle width of the domain for Pr=0.71 and Ra=105.

The present numerical code was validated against the results of benchmark CFD [26]. The results are compared in Table 1. As one can see from this table, there are good agreements between two methods. Furthermore, for Pr = 0.71 and Ra = 105, which are reasonable values for petroleum fluid, the dimensionless horizontal velocity (u*) at mid-width (x_D= 0.5) and dimensionless vertical velocity (v*) at mid-height (y_D= 0.5) of domain are plotted in Figs. (4 and 5) respectively. These figures also show a good agreement between LBE model and CFD results [26]. However, the slight variations in the values are generally due to the employed numerical methods.

Fig. (5). Comparison of variations of dimensionless vertical velocity between LBE model and CFD results at the middle height of domain for Pr=0.71 and Ra=105.

Fig. (6). Streamline contours for different Ra numbers at Pr=0.71.

Figs. (6 and 7) show streamlines and isothermal contours for Pr = 0.71 at different Ra numbers respectively. As could be seen from Fig. (6) for Ra =104, the flow is symmetrical and has circulating motion in the core region and by increasing Ra, some eddies are then observed at the corners. Fig. (7) also illustrates the isothermal contours, which are symmetrical at the onset of convective motion. However, the contours are indeed more distorted for higher Ra numbers. Moreover, Fig. (8) shows local Nusselt number variations on the hot wall (Nu_H) at different Ra numbers. As could be easily deduced from this figure, the Nusselt number increases with the increase of Ra number, which indicates the rising heat transfer by natural convection mechanism. Additionally, to gain more quantitative insight into the problem, a summary of average values of local Nusselt numbers and the location with maximum local Nusselt number is presented in Table 3. Values show that both the location of maximum Nusselt number and the average Nusselt of the wall increase with Ra. It should also be mentioned here that average Nusselt numbers are found by surface numerical integration of the local values on the wall.

Fig. (7). Isothermal contours for different Ra numbers at Pr=0.71.

Fig. (8). Local Nusselt number variations on the hot wall for Pr=0.71 at different Ra numbers.

Table 3. Average Nusselt number of the wall and the location with maximum local Nusselt number.

Ra	104	105	106
(XD)max	0.74	0.70	0.65
	1.8	3.22	5.42

Figs. (9 and 10) illustrate streamline and isothermal contours respectively for Pr=0.71 and Ra=105 at different aspect-ratios. In a low AR, temperature field is progressively more linear and thus indicates a proportional increase in the importance of the purely conductive heat transfer mechanism. Because of the no-slip wall boundary condition, a vigorous fluid motion cannot be propagated in the domain in this case. On the other hand, by increasing AR, several convection cells propagated in domain and temperature field became more unsymmetrical.

Fig. (9). Streamline contours at different aspect-ratios (AR) at Pr=0.71 and Ra=105.

Fig. (10). Isothermal contours for different aspect-ratios (AR) at Pr=0.71 and Ra=105.

In petroleum reservoir, one phase may exist in some operation conditions. Because of that, the results of the fluid natural convection with two different phases in a fractured domain are presented in this section which shows heat transfer in the multi-phases model. Other external forces such as drag force, solid/fluid interaction forces, etc have been considered neglected in this simulation.

For the two-phase flow validation analysis, the dynamics of the adherence of two circular droplets have been tested. For instance, initially, two identical droplets with the radius of Ri = 20 were situated very close to each other. The distance of their centers at the centerline of computation domain (x_D) is considered 0.5. The higher density section (liquid phase) is inside the droplet and lower density section (vapor phase) is situated outside. The computational resolution is considered 100×100 with gravity free and boundary conditions have been employed periodically in all directions. Once the simulation starts, the droplets coalesce instantly.Fig. (11) represents a snapshot of the droplet shapes at four different lattice times. Moreover, the droplet’s radius at y_D = 0.5 as a function of time has been shown in Fig. (12). The droplet’s radius at y_D = 0.5 is the distance from the interface node at y_D = 0.5 to the center of the domain. The fluctuation of the droplet shape is clear in Fig. (12) and that is due to the surface tension. The relationship of has been maintained where R_f represents the final radius of the droplet, and R_i shows the initial radius of two droplets. The result is in a reasonable agreement with other previous investigations [35, 36].

Fig. (11). Snapshots of the two droplets cohesion at different lattice times.

Fig. (12). The variation of droplet’s radius as a function of lattice time.

The existing LBE model can lead to spurious currents which shows the deviation from the reality in the simulation of multiphase flow. These unphysical velocities can be reached to their maximum value at the interface region. The big spurious currents can introduce error in the temperature field and even produce instability in the simulation. The maximum spurious current in this simulation as a function of g_f and its comparison with a previously published work [37] is depicted in Fig. (13). Additionally, fluid/fluid interaction strength, g_f reduces under some critical value g_fC. This shows that a system separates from a single phase to two phases or in other words from heavier liquid phase and a lighter vapor phase. Fig. (14) also represents the density ratio variation versus g_f and its comparison with the results in [37]. It is obvious that the LBE method is valid only for the incompressible limit u/c_s → 0. This means that u is smaller than 0.13 [31, 33]. In our simulation, to maintain a higher accuracy, maximum velocity needs to be smaller than 0.6 with consideration of the velocity field influence on the field temperature. Therefore, in our simulation, we restrict (corresponding to the highest density ratio of 40) which is adequate and valid for many liquid-vapor systems. Obviously, g_f cannot be selected a very small value in order to separate phases.

Fig. (13). Maximum spurious currents in LBE simulation as a function of gf.

Fig. (14). Density ratio in LBE simulation as a function of gf.

To determine the surface tension numerically, the Laplace law has been introduced in the simulation [38]:

(29)

where p_i and p_o are the pressure inside and outside the droplet respectively and r is the droplet radius. Therefore, the surface tension of the SCMP model is numerically determined with producing circular droplets with different radii in a gravity free periodic boundary domain and then the slope of the pressure difference versus the inverse of radius relation has been evaluated. Fig. (15) represents the pressure difference versus the inverse of droplet’s radius for , and the concluded linear relation is pretty well satisfied. The surface tension has been calculated with slope and equation 29 calculations. It was found to be 0.0062 in lattice units in the higher correlation coefficient, which was reasonable for petroleum two-phase flow. However, since σ is related to g_f, the appropriate surface tension is established with this g_f value. Therefore, is applied to the simulation of the thermal multi-phase model.

Fig. (15). Pressure difference versus the inverse of droplet’s radius for (the correlation coefficient of determination is 0.992).

To prevent the droplet size expansion or shrinking, we set the initial droplet density of inside/outside close to the maximum/ minimum density value which is obtained from special g_f value in Fig. (14). Otherwise, the droplet’s radius is not well-behaved and sometimes the droplet may even expand to the boundaries.

The 100×100 lattices are applied and the boundary conditions set as the previous section. The radius of the droplet which was situated at the center of the domain was 20 lattices. Streamline and isothermal contours at different Rayleigh numbers are brought in Figs. (16) and 17) respectively. As it is obvious in Fig. (16), a central vortex appears as a typical feature of the flow for low Ra values and for the Ra value increment, the vortex tends to become elliptical and streamlines tend to break down. At the same time, as Ra increases the isothermal lines tend to become vertical in the center of the domain and horizontal only in the thin boundary layers near the sidewalls (see Fig. (17). The variation of isotherms indicates the change of the dominant heat transfer mechanism from conduction to convection. This outcome is in agreement with the results of natural convection in a cavity domain which was investigated in previous studies [33].

Fig. (16). Streamline contours in SCMP LBE model at different Ra numbers.

Fig. (17). Isothermal contours in SCMP LBE model at different Ra numbers.

Fig. (18) displays the average bulk Nusselt number variations with the time step at various Rayleigh numbers. The results show that the convection cannot occur in Ra=103 and higher Ra is needed in order to the thermally induced buoyancy can overcome the effects of flow resistance. However, the bulk average Nusselt number rises as the Rayleigh number increases. Moreover, the magnitude of the fluctuations also increases with the Rayleigh number increment. Therefore, the stability of model decreases with increasing the Rayleigh number in the SCMP model.

Fig. (18). The average bulk Nusselt number and the magnitude of its fluctuations versus time step.