Simulation of Boiling Heat Transfer at Different Reduced Temperatures with an Improved Pseudopotential Lattice Boltzmann Method

Abstract

The pseudopotential Lattice Boltzmann Method has attracted much attention in the recent years for the simulation of boiling heat transfer. Many studies have been published recently for the simulation of the bubble cycle (nucleation, growth and departure from a heated surface). This paper puts forward two-dimensional simulations of bubble nucleation, growth and departure using an improved pseudopotential Lattice Boltzmann Model from the literature at different reduced temperatures, $T_r=0.76$ and $T_r=0.86$. Two different models using the Bhatnagar–Gross–Krook (BGK) and the Multiple-Relaxation-Time (MRT) collision operators with appropriate forcing schemes are used. The results for pool boiling show that the bubbles exhibit axial symmetry during growth and departure. Numerical results of departure diameter and release period for pool boiling are compared against empirical correlations from the literature by varying the gravitational acceleration. Reasonable agreement is observed. Nucleate boiling trends with heat flux are also captured by the simulations. Numerical results of flow boiling simulations are compared by varying the Reynolds number for both reduced temperatures with the MRT model. It was found that the departure diamenter and release period decreases with the increase of the Reynolds number. These results are a direct effect of the drag force. Proper conclusions are commented at the end of the paper.

1. Introduction

In the recent years, the pseudopotential Lattice Boltzmann Method (LBM) has emerged as a powerful CFD tool, showing great potential and capabilities to handle complex two-phase flow phenomena [1,2,3,4,5]. Among them, boiling heat transfer is very important in many industrial applications. In regard to boiling heat transfer modes, nucleate boiling has been recognized as one of the most effective, and has been used in different cooling applications, such as nuclear reactors [6], computer chips [7] and electronic devices [8]. Regarding the bubble cycle, the bubble nucleation happens when the liquid near the surface is heated until a certain superheating degree is achieved and a bubble nucleus is formed. After the nucleation, the bubble starts the growth stage, due to inertial and heat transfer effects. Finally, the bubble departs due to the buoyancy effects. Regarding these phenomena, a detailed comparison of the results of bubble departure diameter and release period at reduced temperatures below 0.8 for pool boiling and flow boiling applications are still scarce. Recent publications on boiling heat transfer also focused on fundamental aspects, such as the effect of the surface roughness [9]. From the theoretical point of view, nucleate boiling is extremely complex, involving processes of nucleation, growth, departure and coalescence of vapor bubbles. The concept of symmetry is clearly seen during the bubble cycle, since the bubble exhibit axial symmetry during growth and departure.

Several papers have been published, applying the pseudopotential LBM for the simulation of the bubble cycle (nucleation, growth and departure from a heated surface), some of them with limitations of the numerical model. For example, Zhang and Chen [10] developed a hybrid Lattice Boltzmann Method based on the pseudopotential model to simulate boiling heat transfer. The LBM was applied to solve for the density and velocity fields and the extended Lax-Wendroff scheme for the solution of the temperature field. However, the authors disregarded the influence of pressure variation across the interface. Later, Hazi and Markus [11] developed a method based on the pseudopotential method to simulate nucleate boiling. The authors employed two particle distribution functions, one for the density and velocity fields and another for the temperature field. However, the authors adopted an artificial equation of state to couple the hydrodynamic and thermal models. This assumption made by Hazi and Markus [11] might be not reliable to simulate a real physical system.

The study carried out by Gong and Cheng [12] proposed an improved phase-change LBM based on a modification of the original Shan-Chen multiphase model. This new method directly incorporates a non-ideal equation of state into the model. The thermodynamic consistency is ensured by a modification of the calculation of the interparticle force. The authors solved the density, velocity and temperature fields using the LBM, assuming two distribution functions. A validation of the numerical scheme was performed by comparing the bubble departure diameter with the correlation proposed by Fritz [13]. A reasonable agreement was observed. Later, the authors published studies using this phase-change LBM to evaluate surface wettability, [14,15]. In a subsequent study, the authors used the model to analyzed the influence of single and multi-cavities during the nucleation of a bubble [16,17]. However, in all these studies, the energy equation did not take into account the density variation across the interface in the diffusive term [18]. Regarding the density variation across the interface, Gong and Cheng [19] considered this effect. In their study, they evaluated the effect of the thermal response of the heater over the simulation results. In all these studies, the Bhatnagar–Gross–Krook (BGK) collision operator was used.

Li et al. [20] were pioneers to perform simulations of different boiling regimes at different surface wettability conditions using the pseudopotential LBM. In their study, they developed an improved pseudopotential Lattice Boltzmann Method for the simulation of boiling heat transfer. The Multiple-Relaxation-Time (MRT) collision operator was used to simulate the density and velocity fields and the temperature field was obtained using a Runge–Kutta scheme. An improved forcing scheme was applied in order to achieve thermodynamic consistency. The numerical results were presented for different wall superheat and different wettability conditions and agreed reasonably well with the theoretical results.

Based on the previous developed research studies showed before, recent works focused on more applications at different conditions. For example, Hu and Liu [21] studied boiling heat transfer using the improved pseudopotential LBM proposed by Li et al. [20] with an innovative nucleation site treatment on the wall. Then, the authors simulated the different boiling stages and were able to reproduce the boiling curve and well as all stages of the bubble cycle were captured.

Chang et al. [22] performed numerical simulations pool boiling heat transfer enhancement on structured surfaces with columns using the pseudopotential LBM proposed by Gong and Cheng [12]. The authors analyzed the effects of geometrical parameters of the structured surfaces over the boiling curves. The main conclusion drawn by the authors was that, in general, two factors affect the enhanced heat transfer effect of pool boiling on the structured surfaces: the surface area and convection near the bubble.

Zhang et al. [23] studied the effects of wall superheat and surface wettability on nucleation site interactions using the pseudopotential LBM proposed by Li et al. [20]. According to the authors, most boiling heat transfer correlations assume that the nucleation and growth of bubbles at adjacent nucleation sites are independent. However, the experimental results show some interactions between adjacent nucleation sites, while, the experimental results have some divergence due to the complexity of boiling [24,25,26]. Using numerical simulations, they were able to show that the wall superheat could influence the nucleation site interactions by changing the relatively intensity of thermal interaction and hydrodynamic interaction.

Despite all the aforementioned publications that focused on pseudopotential LBM application for the simulation boiling heat transfer, to the best of authors’ knowledge, none of them evaluated the influence of reduced temperature variation on bubble departure diameter and release period results for pool and flow boiling conditions. Most publications only simulate conditions with reduced temperature close to the critical point, which raises a question regarding the simulation of real boiling process for medium to lower reduced temperatures ($T_r\le 0.8$), commonly used in several experimental set-ups, as in Alvarino et al. [9], for example. Hence, the present paper focuses on the use of pseudopotential LBM for simulating both, pool and flow, boiling heat transfer problems investigating the influence of reduced temperature on boiling dynamics, presenting results that include the bubble departure diameter and release period and the average heat flux from the heated surface. Particularly, two thermodynamic conditions, namely $T_r=0.76$ and $T_r=0.86$, were employed in order to investigate pool and flow boiling under influence of distinct values of gravitational and inertial effects. It should be noted that the majority of the simulation results presented in the open literature corresponds to reduced temperatures above $T_r=0.80$. The simulations are performed using the pseudopotential LBM with two improved forcing schemes from the literature to ensure thermodynamic consistency, namely Li et al. [27] and Li et al. [28], along with the D2Q9 velocity scheme. It is worth mentioning that both models obtained from the literature were developed using the BGK and MRT collision operators, respectively. In order to differentiate these two models, they are mentioned in the text as the BGK and MRT models.

Regarding the numerical results, at first, pool boiling results are presented. The bubble cycle is studied and the results of departure diameter and release period are compared with empirical correlations from the literature by varying the gravitational acceleration. The average heat flux is shown for both reduced temperatures. Later, bubble cycle simulations under the influence of a fully-developed flow are performed. In this case, besides the reduced temperature influence, inertial effects on the departure diameter and release period are also investigated by varying the Reynolds number, presenting also the average heat flux temporal variation.

2. Hydrodynamic Model: Improved Pseudopotential Lattice Boltzmann Method

The Boltzmann Equation describes the evolution of a distribution function, $f(\mathbf{u},\mathbf{x},t)$, of a collection of particles. This function represents the density of particles at a given time t, in a position given by $\mathbf{x}=[x_1,x_2,x_3]$, with particle velocity $\mathbf{u}=[u_1,u_2,u_3]$, under the effect of an external force $\mathbf{F}$:

$$\frac{\partial f(\mathbf{u},\mathbf{x},t)}{\partial t}+\mathbf{u}\frac{\partial f(\mathbf{u},\mathbf{x},t)}{\partial \mathbf{x}}+\frac{\mathbf{F}}{m}\frac{\partial f(\mathbf{u},\mathbf{x},t)}{\partial \mathbf{u}}=\Omega (f,f^{eq})$$

(1)

In Equation (1), $f^{eq}$ is the equilibrium distribution function, which describes the equilibrium condition of the collection of particles. The right-hand side is the collision operator, which represents the binary collision between particles. The Lattice Boltzmann Equation is a discrete form of the Boltzmann equation for the density distribution function. The discretization is first based on the introduction of a set of discrete velocities, namely ${\mathbf{c}}_i$, for the simulated particles, resulting in velocity schemes, namely lattices, which are similar to a computational mesh. Later, the equation is discretized in the space and time. This last discretization step can be performed using either the method of characteristics or the finite difference method [29]. In either case, the Lattice Boltzmann Equation is obtained:

$$f_i(\mathbf{x}+{\mathbf{c}}_i\Delta t,t+\Delta t)=f_i(\mathbf{x},t)+\Delta t\Omega (f_i,f_i^{eq})+\Delta t{F}_i$$

(2)

where the subscript i stands for the values related to each of the discretized value of velocity set ${\mathbf{c}}_i$, and the lattice spacing, $\Delta x=|{\mathbf{c}}_i|\Delta t$, and the timestep $\Delta t$, are usually set to 1.

The most difficult step to solving the Lattice Boltzmann Equation is the collision operator. The most popular approach is to consider the BGK collision operator [30], which considers the collision of two simulated particles to be related only by one relaxation parameter, $\frac{\Delta t}{{\tau }_{\nu }}$. In this case, the discretized evolution equation with the BGK collision operator for the density distribution function is represented by:

$$f_i(\mathbf{x}+{\mathbf{c}}_i\Delta t,t+\Delta t)=f_i(\mathbf{x},t)-\frac{\Delta t}{{\tau }_{\nu }}(f_i(\mathbf{x},t)-f_i^{eq})+\Delta t{F}_i.$$

(3)

The discretized equilibrium distribution function, $f_i^{eq}$, is obtained by equalizing its continuum form with the Hermite’s quadrature associated with some weighting coefficients, $w_i$. It should be mentioned that these weighting coefficients should satisfy specific conservation relations [29]. The discretized form of the equilibrium distribution function can be expressed as:

$$f_i^{eq}=w_i\rho \left(1+\frac{{\mathbf{c}}_i.\mathbf{v}}{{c_s}^2}+\frac{{({\mathbf{c}}_i.\mathbf{v})}^2}{{c_s}^2}-\frac{{\left|\mathbf{v}\right|}^2}{2{c_s}^2}\right),$$

(4)

where $c_s$ is a constant—namely, the speed of sound—which is a function of the chosen velocity scheme, $\mathbf{v}$ is the fluid macroscopic velocity and $\rho$ is the fluid density.

In this paper, the D2Q9 velocity scheme is employed to perform the two-dimensional numerical simulations. A schematic representation of this scheme is shown in Figure 1:

For this velocity scheme, the speed of sound is ${c_s}^2=\frac{1}{3}$ and the weighting coefficients, $w_i$ and the discrete velocities, ${\mathbf{c}}_i$, are given by, respectively [29]:

$$w_i=\left\{\begin{array}{cc}\frac{1}{9},\hfill & \mathrm{for}i=1:4\hfill \\ \frac{1}{36},\hfill & \mathrm{for}i=5:8\hfill \\ \frac{4}{9},\hfill & \mathrm{for}i=9\hfill \end{array}\right.$$

(5)

$${\mathbf{c}}_i=\left\{\begin{array}{cc}(\pm 1,0),(0,\pm 1)\hfill & \mathrm{for}i=1:4\hfill \\ (\pm 1,1)\hfill & \mathrm{for}i=5:8\hfill \\ (0,0)\hfill & \mathrm{for}i=9\hfill \end{array}\right.$$

(6)

During boiling heat transfer simulations in the mesoscopic scale, the following forces must be included: the interparticle interaction force, ${\mathbf{F}}_{int}$, which accounts for phase separation and the gravitational force, ${\mathbf{F}}_g$. The interaction force between the liquid and solid wall, which allows a tunable wettability condition, is not considered in the present work.

The interparticle interaction force model proposed by Shan and Chen [31], which is computed based on the pseudopotential, $\psi$, is expressed as:

$${\mathbf{F}}_{int}=-G\psi \left(\mathbf{x}\right)\sum _i{w}_i^{\ast }\psi (\mathbf{x}+{\mathbf{c}}_i\Delta t){\mathbf{c}}_i\Delta t,$$

(7)

where $w_i^{\ast }$ are the redefined weighting coefficients [21] and G is the interaction strength. Following Yuan and Schaefer [32], the pseudopotential is computed from an Equation of State (EOS) for real gases as follows:

$$\psi \left(\mathbf{x}\right)=\sqrt{\frac{2\left(p-\frac{\rho }{3}\right)}{|{\mathbf{c}}_i|G}},$$

(8)

where p is the thermodynamic pressure, which can computed from an Equation of State (EOS) for real gases. By using the previous definition for the pseudopotential, the parameter G no longer represents the interaction strength, and it is used only to ensure that the square root is positive. In all simulations performed in this paper, $G=-1$ is applied.

In this study, the Peng–Robinson EOS is adopted to compute the pressure, since it can handle density ratio as high as 1000 and allows the simulation at conditions at lower reduced temperatures ($T_r\sim 0.59$):

$$p=\frac{\rho RT}{1-b\rho }-\frac{a{\rho }^2{\left[1+(0.37464+1.54226\omega -0.2699{\omega }^2)\left(1-\sqrt{\frac{T}{T_c}}\right)\right]}^2}{1+2b\rho -b^2{\rho }^2}.$$

(9)

In Equation (9), $T_c$ is the critical temperature and $\omega$ is the acentric factor. Parameters $a=\frac{0.45724R^2{T}_c^2}{p_c}$ and $b=\frac{0.0778R{T}_c}{p_c}$, in which $p_c$ is the pressure at the critical point and R is the universal gas constant. The parameters of the Peng-Robinson EOS were chosen as proposed by Li et al. [20]: $a=3/49$, $b=2/21$, $R=1$ and $\omega =0.344$. Using these values, the critical temperature is given by $T_c=0.1094$. It should be mentioned that all quantities are presented in lattice units, as commonly performed in the literature.

The gravitational force effect is included as follows:

$${\mathbf{F}}_g=\mathbf{g}(\rho -\overline{\rho }),$$

(10)

where $\overline{\rho }$ is the average density computed in the computational domain, and $\mathbf{g}$ is the gravitational acceleration.

Focusing on the model using the BGK collision operator, the previous forces are incorporated into the numerical model using two different forcing schemes. For the interparticle force, the forcing scheme proposed by Li et al. [27] is applied, in order to ensure the thermodynamic consistency:

$$F_{i,1}=\left(1-\frac{1}{2{\tau }_{\nu }}\right)w_i\left(\frac{{\mathbf{c}}_i-{\mathbf{v}}_{new}}{c_s^2}+\frac{{\mathbf{c}}_i·{\mathbf{v}}_{new}}{c_s^4}{\mathbf{c}}_i\right)·{\mathbf{F}}_{int}$$

(11)

where the modified velocity, ${\mathbf{v}}_{new}$, is given by:

$${\mathbf{v}}_{new}=\mathbf{v}+\frac{{\sigma }_{BGK}{\mathbf{F}}_{int}}{({\tau }_{\nu }-0.5){\psi }^2}$$

(12)

In this case, the pressure tensor is changed and its new form is expressed as [27]:

$$\mathbf{P}={\mathbf{P}}_{original}+2G^2{\left|{\mathbf{c}}_i\right|}^4{\sigma }_{BGK}\nabla \psi \nabla \psi$$

(13)

In Equation (13), ${\mathbf{P}}_{original}$ is the original pressure tensor and ${\sigma }_{BGK}$ is a parameter to adjust the mechanical stability condition in order to ensure the thermodynamic consistency. It should be noticed that when the forcing scheme proposed by Li et al. [27] is applied, only the anisotropic part of the pressure tensor is changed.

For the gravitational force, the traditional forcing scheme proposed by Guo et al. [33] is applied:

$$F_{i,2}=\left(1-\frac{1}{2{\tau }_{\nu }}\right)w_i\left(\frac{{\mathbf{c}}_i-\mathbf{v}}{c_s^2}+\frac{{\mathbf{c}}_i·\mathbf{v}}{c_s^4}{\mathbf{c}}_i\right)·{\mathbf{F}}_g.$$

(14)

The total external force to be implemented into the evolution equation for the particle distribution function (see Equation (3)) is then $F_i=F_{i,1}+F_{i,2}$. It should be noticed that the forcing scheme [27] will reduce to the forcing scheme [33] when ${\sigma }_{BGK}=0$.

From the Chapman–Enskog expansion for the LBM, the relation between the kinematic viscosity, $\nu$, and the relaxation parameter, ${\tau }_{\nu }$, can be established [19] as follows:

$$\nu =c_s^2({\tau }_{\nu }-0.5)\Delta t$$

(15)

The BGK collision operator has been considered in the literature to simulate boiling heat transfer [12,14,15,16,17,19], especially due to its simplicity. However, this collision operator suffers from numerical instabilities as the Reynolds number increases. This collision operator is not suitable for simulating lower reduced temperatures, since the spurious current near the interface may create strong instabilities. In this sense, the MRT collision operator can be considered.

The MRT collision operator can be seen as an extension of the BGK collision operator, being superior over the BGK formulation in terms of numerical stability. This is associated with possibility to adjust the relaxation parameters individually, in order to achieve optimal stability conditions.

Using the standard MRT collision operator [34,35], the discretized evolution equation for the density distribution function is given by:

$$f_i(\mathbf{x}+{\mathbf{c}}_i\Delta t,t+\Delta t)=f_i(\mathbf{x},t)+\Delta t({\mathbf{M}}^{-\mathbf{1}}\Lambda \mathbf{M})(f_j(\mathbf{x},t)-f_j^{eq})+\Delta t{F}_i,$$

(16)

where $\mathbf{M}={\mathrm{M}}_{ij}$ is an orthogonal transformation matrix and $\Lambda$ is a diagonal matrix, which includes all the relaxation times.

The transformation matrix M for the D2Q9 model is given by:

$$\mathbf{M}=\left[\begin{array}{ccccccccc}1& 1& 1& 1& 1& 1& 1& 1& 1\\ -4& -1& -1& -1& -1& 2& 2& 2& 2\\ 4& -2& -2& -2& -2& 2& 2& 2& 2\\ 0& 1& 0& -1& 0& 1& -1& -1& 1\\ 0& -2& 0& 2& 0& 1& -1& -1& 1\\ 0& 0& 1& 0& -1& 1& 1& -1& -1\\ 0& 0& -2& 0& 2& 1& 1& -1& -1\\ 0& 1& -1& 1& -1& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 1& -1& 1& -1\end{array}\right]$$

(17)

Using this transformation matrix, the collision step in the moment space, which represents the right-hand side of Equation (16) is performed as [36]:

$${\mathbf{m}}^{\ast }=\mathbf{m}-\Lambda (\mathbf{m}-{\mathbf{m}}^{eq})+\left(\mathbf{I}-\frac{\Lambda }{2}\right)\mathbf{S},$$

(18)

where $\mathbf{m}=\mathbf{M}f$, ${\mathbf{m}}^{eq}=\mathbf{M}{f}^{eq}$, $\mathbf{I}$ is the unitary tensor, $\mathbf{S}$ is the forcing term in the moment space with $(\mathbf{I}-0.5\Lambda )\mathbf{S}=\mathbf{M}\mathbf{F}$.

The diagonal matrix $\Lambda$ is given by Lallemand and Luo [37]:

$$\Lambda =\mathrm{diag}({\tau }_e^{-1},{\tau }_{\zeta }^{-1},{\tau }_j^{-1},{\tau }_q^{-1},{\tau }_j^{-1},{\tau }_q^{-1},{\tau }_{\nu }^{-1},{\tau }_{\nu }^{-1},{\tau }_{\rho }^{-1}),$$

(19)

where ${\tau }_{\rho }$ and ${\tau }_j$ are the relaxation parameters of conserved moments and are set to 1, ${\tau }_{\nu }$ determines the dynamic viscosity (Equation (15)), ${\tau }_e$ is associated with the bulk viscosity and is set to 1.1, ${\tau }_{\zeta }^{-1}$ is related to energy square and is set to 1.1 and ${\tau }_q^{-1}$ is related to the energy flux and is set to 1.1.

From ${\mathbf{m}}^{eq}=\mathbf{M}{f}^{eq}$, the equilibrium momentum ${\mathbf{m}}^{eq}$ can be expressed as:

$${\mathbf{m}}^{eq}={\rho (1,-2+3|\mathbf{v}|}^2,1-{3\left|\mathbf{v}\right|}^2,v_x,-v_x,v_y,-v_y,v_x^2-v_y^2,v_x,v_y{)}^T,$$

(20)

where ${\left|\mathbf{v}\right|}^2=v_x^2+v_y^2$.

The forcing scheme proposed by Li et al. [28] is applied in order to ensure the thermodynamic consistency, by adjusting the mechanical stability condition. In this scheme, the forcing term $\mathbf{S}$ is given by:

$$\mathbf{S}=\left[\begin{array}{c}0\\ 6\mathbf{v}·\mathbf{F}+\frac{{\sigma }_{MRT}{\left|{\mathbf{F}}_{int}\right|}^2}{{\psi }^2({\tau }_e-0.5)}\\ -6\mathbf{v}·\mathbf{F}+\frac{{\sigma }_{MRT}{\left|{\mathbf{F}}_{int}\right|}^2}{{\psi }^2({\tau }_e-0.5)}\\ F_x\\ -F_x\\ F_y\\ -F_y\\ 2(v_x{F}_x-v_y{F}_y)\\ v_x{F}_y+v_y{F}_x\end{array}\right]$$

(21)

where ${\sigma }_{MRT}$ is a parameter used to tune the mechanical stability condition, $|{\mathbf{F}}_{int}{|}^2=F_{int,x}^2+F_{int,y}^2$, $\mathbf{F}={\mathbf{F}}_{int}+{\mathbf{F}}_g$ is the total force, ${\left|\mathbf{F}\right|}^2=F_x^2+F_y^2$, where $F_x$ and $F_y$ are the total force components, $v_x$ and $v_y$ are the fluid velocity components.

Using this forcing scheme, the new form of the pressure tensor is given by:

$$\mathbf{P}={\mathbf{P}}_{original}+2G^2|{\mathbf{c}}_i{|}^4{\sigma }_{MRT}{|\nabla \psi |}^2\mathbf{I}.$$

(22)

Similar to the forcing scheme proposed by Li et al. [27] for the model using the BGK collision operator, in the scheme proposed by Li et al. [28], the mechanical stability condition is slightly changed so that the thermodynamic consistency can be achieved. It should be noticed that when this forcing scheme is applied, only the isotropic part of the pressure tensor is changed.

However, a major drawback of the traditional pseudopotential method is the impossibility to adjust the surface tension independently of the density ratio. In order to overcome this limitation, the new method proposed by Li and Luo [38] is employed to ensure a tunable surface tension while it ensures that the density ratio is unchanged.

This approach is based on the addition of a new source term, $\mathbf{C}$, to the collision step (Equation (18)) of the model using the MRT collision operator. According to the authors, the collision step in moment space is modified as follows [38]:

$${\mathbf{m}}^{\ast }=\mathbf{m}-\Lambda (\mathbf{m}-{\mathbf{m}}^{eq})+\left(\mathbf{I}-\frac{\Lambda }{2}\right)\mathbf{S}+\Delta t\mathbf{C}.$$

(23)

In the model proposed by Li and Luo [38], the additional source term $\mathbf{C}$ has the following form:

$$\mathbf{C}=\left[\begin{array}{c}1.5{\tau }_e^{-1}(Q_{xx}+Q_{yy})\\ 1.5{\tau }_{\zeta }^{-1}(Q_{xx}+Q_{yy})\\ 0\\ 0\\ 0\\ 0\\ {\tau }_{\nu }^{-1}(Q_{xx}-Q_{yy})\\ {\tau }_{\zeta }^{-1}{Q}_{xy}\\ 0\end{array}\right]$$

(24)

where the terms $Q_{xx}$, $Q_{yy}$ and $Q_{xy}$ are computed as described in Equation (25):

$$\mathbf{Q}=\kappa \frac{G}{2}\psi \left(\mathbf{x}\right)\sum _i{w}_i^{\ast }[\psi (\mathbf{x}+{\mathbf{c}}_i\Delta t)-\psi \left(\mathbf{x}\right)]{\mathbf{c}}_i{\mathbf{c}}_i,$$

(25)

where $\kappa$ is a parameter that allows a tunable surface tension.

The solution of the Lattice Boltzmann Equation allows the computation of the macroscopic fields. In this case, the density and velocity can be obtained directly from the density distribution function using the following equations, respectively:

$$\rho =\sum _i{f}_i$$

(26)

$$\mathbf{v}=\frac{{\sum }_i{\mathbf{c}}_i{f}_i}{\rho }+\frac{{\mathbf{F}}_{int}+{\mathbf{F}}_g}{2\rho }.$$

(27)

3. Modeling the Energy Conservation Equation

The energy equation is obtained considering the diffuse-interface model [39] and can be expressed as follows:

$$\rho \left(\frac{\partial e}{\partial t}+\mathbf{v}·\nabla e\right)=\nabla ·(k\nabla T)-p\nabla ·\mathbf{v},$$

(28)

where e is the internal energy, k is the thermal conductivity and T is the fluid temperature. It is worth mentioning that viscous dissipation as well as the kinetic and potential energies were not taken into account.

Using the well-known relation for a pure substance: $de=c_{v}dT+\frac{1}{{\rho }^2}\left[p-T{\left(\frac{\partial p}{\partial T}\right)}_{\rho }\right]d\rho$, an evolution equation for the temperature field can be obtained:

$$\frac{\partial T}{\partial t}=-\mathbf{v}·\nabla T+\frac{\nabla ·(k\nabla T)}{\rho c_v}-\frac{T}{\rho c_v}{\left(\frac{\partial p}{\partial T}\right)}_{\rho }\nabla ·\mathbf{v}.$$

(29)

It should be noticed that in Equation (29), the coupling between the hydrodynamic and thermal models is presented in the last term, where $c_v$ is the constant volume specific heat.

The energy equation is solved using the fourth-order Runge–Kutta method. The resulting discrete evolution equation for the temperature is given by:

$$T^{t+\Delta t}=T+\frac{\Delta t}{6}(h_1+2h_2+2h_3+h_4),$$

(30)

where the functions $h_1$, $h_2$, $h_3$ and $h_4$ are given by:

$$h_1={\left[-\mathbf{v}·\nabla T+\frac{\nabla ·(k\nabla T)}{\rho c_v}-\frac{T}{\rho c_v}{\left(\frac{\partial p}{\partial T}\right)}_{\rho }\nabla ·\mathbf{v}\right]}_{T^t}$$

(31)

$$h_2={\left[-\mathbf{v}·\nabla T+\frac{\nabla ·(k\nabla T)}{\rho c_v}-\frac{T}{\rho c_v}{\left(\frac{\partial p}{\partial T}\right)}_{\rho }\nabla ·\mathbf{v}\right]}_{T^t+\Delta t\frac{h_1}{2}}$$

(32)

$$h_3={\left[-\mathbf{v}·\nabla T+\frac{\nabla ·(k\nabla T)}{\rho c_v}-\frac{T}{\rho c_v}{\left(\frac{\partial p}{\partial T}\right)}_{\rho }\nabla ·\mathbf{v}\right]}_{T^t+\Delta t\frac{h_2}{2}}$$

(33)

$$h_4={\left[-\mathbf{v}·\nabla T+\frac{\nabla ·(k\nabla T)}{\rho c_v}-\frac{T}{\rho c_v}{\left(\frac{\partial p}{\partial T}\right)}_{\rho }\nabla ·\mathbf{v}\right]}_{T^t+\Delta t{h}_3},$$

(34)

in which $T^t$ is the temperature field computed at instant t. In order to compute the spatial discretization of the specific quantities, the isotropic schemes proposed by Lee and Lin [40] are employed.

4. Results and Discussion

In this section, the simulation results are presented. In a first moment, the simulation results obtained for the reduced temperature equal to $T_r=T_{sat}/T_c=0.76$ are presented. After that, a detailed discussion and comparison of with those obtained for $T_r=0.86$. The differences between the results obtained using the BGK and MRT collision operators are carefully addressed. The parameters to ensure the thermodynamic consistency for the models using the BGK and MRT collision operators are chosen as ${\sigma }_{BGK}=0.105$ [27] and ${\sigma }_{MRT}=0.1$ [28], respectively. Further details about the forcing schemes can be seen in the mentioned references.

For all simulations, the initial velocity field is considered equal to zero. Using the density and velocity fields, the initialization of the particle distribution function is performed using the equilibrium distribution function (see Equation (4)).

4.1. Computation of the Surface Tension: Young–Laplace Test

In order to obtain the surface tension in the LBM, the Young–Laplace test is performed. In this case, the pressure difference, $\Delta p$, across the interface of a droplet with a radius $R_d$ in equilibrium within the vapor phase is related to the surface tension, $\gamma$, as:

$$\gamma =\frac{\Delta p}{R_d}.$$

(35)

At the defined reduced temperature, $T_r$, the equilibrium densities are obtained from the chosen EOS using the Maxwell Construction Rule [29]. In this case, the thermodynamic coexistence is ensured according to the following relation:

$$\underset{{\rho }_v}{\overset{{\rho }_l}{\int }}(p_0-p)\frac{1}{{\rho }^2}d\rho =0.$$

(36)

where p is the pressure computed by the non-ideal EOS and $p_0=p\left({\rho }_l\right)=p\left({\rho }_v\right)$ [28]. It should be highlighted that all the parameters are set in lattice units, as well as the computed equilibrium densities.

At $T_r=0.86$, the liquid and vapor equilibrium densities are ${\rho }_l=6.50$ and ${\rho }_v=0.38$, respectively. The kinematic viscosity of the liquid is taken as ${\nu }_l=0.1$ and used to compute the relaxation time, following Li et al. [18]. In this case, the computed relaxation time (see Equation (15)) assumes a constant value over the whole computational domain: ${\tau }_{\nu }=0.8$.

A similar procedure is applied at $T_r=0.76$. In this case, the equilibrium densities are obtained: ${\rho }_l=7.59$ and ${\rho }_v=0.12$. However, a more refined set of parameters is applied for each phase for the simulations at this temperature: ${\nu }_l=0.2$ and ${\nu }_v=0.01$, respectively. Using these values, a distribution of the kinematic viscosity inside the computational domain can be computed using the following interpolation scheme [12]:

$$\nu =\frac{\rho -{\rho }_v}{{\rho }_l-{\rho }_v}{\nu }_l+\frac{{\rho }_l-\rho }{{\rho }_l-{\rho }_v}{\nu }_v.$$

(37)

In this case, the computed relaxation time will vary in the computational domain according to the kinematic viscosity, as shown in Equation (15). The computed values for the relaxation time for the liquid and vapor are equal to 1.1 and 0.53, respectively. This is the range of the relaxation time for the simulations performed at $T_r=0.76$.

The static droplet is obtained using the following density distribution [41]:

$$\rho (x,y)=\frac{{\rho }_l+{\rho }_v}{2}+\left(\frac{{\rho }_l-{\rho }_v}{2}\right)tanh\left[\frac{2\left(R(x,y)-R_d\right)}{W}\right],$$

(38)

where $R(x,y)=\sqrt{{(x-x_0)}^2+{(y-y_0)}^2}$, $(x_0,y_0)$ is the coordinate of the center of the droplet, W is the interface width and $R_d$ is the droplet radius. In the present work, $W=5$ and $R_d=50$. The simulation is performed until the density field achieves equilibrium, considering a relative tolerance of ${10}^{-10}$ between two consecutive iterations. The computation domain is 200 × 200. Periodic boundary conditions are set at all boundaries.

Figure 2 presents the results of the density field obtained from simulations with the models using the BGK and MRT collision operators at $T_r=0.76$. In this case, the computed surface tension obtained are equal to 0.382 and 0.177, respectively, using the models with the BGK and MRT collision operators. At $T_r=0.86$, the density field was obtained by the same procedure. The computed surface tension are equal to 0.1881 and 0.0858, respectively.

These surface tension results reveal an interesting feature of the improved forcing schemes proposed by Li et al. [27] and Li et al. [28], regarding the models that are based on the BGK and MRT collision operators, respectively. Since the authors proposed different methodologies to approximately achieve thermodynamic consistency by modifying the pressure tensor, the computed surface tension is not the same for the same reduced temperature. Hence, the method proposed by Li and Luo [38] is used to adjust the surface tension computed by the MRT based model, in order to obtain the same value obtained by the BGK based model.

By performing numerical simulations considering different values of the parameter $\kappa$ (see Equation (25)) for the MRT model, relations between this parameter and the surface tension are obtained considering the different reduced temperatures. These relations are presented in details in Figure 3. Using these results, the surface tension of the MRT collision operator is adjusted to match the one obtained by the BGK model. At $T_r=0.86$, the parameter is set to $\kappa =-1.165$ and the computed surface tension for the MRT model is 0.1881. Similarly, at $T_r=0.76$, when $\kappa =-1.126$, the computed surface tension is 0.382. At this point, a consistent comparison between the numerical results obtained by the BGK and MRT models can be performed.

The surface tension is higher for $T_r=0.76$, due to a lower saturation pressure. These differences in the surface tension lead to bubbles with a higher departure diameter.

4.2. Single Bubble Nucleation: Bubble Departure Diameter and Release Period versus Gravitational Acceleration

For the simulation of the bubble cycle, the computational domain is initialized with a liquid-vapor distribution, considering the same occupied area for each phase [18], as shown in Figure 4.

The upper and lower surfaces are modeled as walls using the scheme proposed by Zou and He [42]. This scheme allows the computation of the unknown density distribution functions at the solid walls. After the streaming step at the upper surface, the unknown density distribution functions are $f_4$, $f_7$ and $f_8$. For the lower surface, the unknown density distribution functions are $f_2$, $f_5$ and $f_6$. The forces acting at these boundary surfaces are included following a methodology introduced by Krüger et al. [29]. For the upper boundary, the following equations comprise the no-slip boundary condition:

$$f_4=f_2+\frac{F_y}{6}$$

(39)

$$f_7=f_5+\frac{1}{2}(f_1-f_3)+\frac{F_x}{4}+\frac{F_y}{6}$$

(40)

$$f_8=f_6-\frac{1}{2}(f_1-f_3)-\frac{F_x}{4}+\frac{F_y}{6}$$

(41)

Similarly, for the lower boundary, the following equations are applied:

$$f_2=f_4-\frac{F_y}{6}$$

(42)

$$f_5=f_7-\frac{1}{2}(f_1-f_3)-\frac{F_x}{4}-\frac{F_y}{6}$$

(43)

$$f_6=f_8+\frac{1}{2}(f_1-f_3)+\frac{F_x}{4}-\frac{F_y}{6}.$$

(44)

Periodic boundary conditions are specified at the lateral boundaries for the density distribution function. The implementation of this boundary condition is straightforward in the LBM [29]. For the temperature field, periodic boundary conditions are also specified at the lateral boundaries. As for the heated surface, constant temperature is used for the pool boiling simulations.

The departure diameter and the release period play an important role on boiling heat transfer. These quantities are chosen to assess the results of the numerical model. The numerical results are compared with the empirical correlations proposed by Fritz [13] and Zuber [43], which are given, respectively, by:

$$D_b=0.0208\theta {\left[\frac{\gamma }{({\rho }_l-{\rho }_v)g}\right]}^{0.5}$$

(45)

$$T_b=\frac{D_b}{0.59}{\left[\frac{\gamma ({\rho }_l-{\rho }_v)g}{{\rho }_l^2}\right]}^{-0.25},$$

(46)

where $\theta$ is the contact angle and g is the value of the gravitational acceleration. From Equations (45) and (46), it can be noticed that $D_b\propto g^{-0.5}$ and $T_b\propto g^{-0.75}$, respectively.

For the initial conditions we had the velocity field, $\mathbf{v}$, with zero value, and the density field was initialized as liquid density ${\rho }_l$ for the lower half of the computational domain, and gas density ${\rho }_v$ for the upper half. For the simulations performed at $T_r=0.86$, the thermal parameters were taken from Li et al. [18]. In this case, $c_v=5.0$ and the thermal conductivity is related to the density field according to $k=\alpha \rho c_v$. The thermal diffusivity is taken as a constant value: $\alpha =0.06$. For the simulations performed at $T_r=0.76$, the following parameters were chosen: $c_{v,l}=6.98$, $c_{v,v}=5.26$, ${\alpha }_l=0.05$ and ${\alpha }_v=0.1$. Using these values, distributions of $c_v$ and $\alpha$ were obtained using Equation (37). The computational mesh was chosen as 300 × 151. The heated surface is considered to be located at the center node of the lower surface and at its immediate two neighbor nodes. At these nodes, the temperature was set to $T_b=1.25T_c$. Hence, the wall superheat was $\Delta T=T_b-T_{sat}=0.0427$ at $T_r=0.86$ and $\Delta T=0.053$ at $T_r=0.76$.

The numerical results of the bubble departure diameter and the release period versus the gravitational acceleration are presented in Figure 5. It should be mentioned that these quantities are presented only for the first bubble. Following Krüger et al. [29], the interface was defined as the region which has the average density computed using liquid and vapor densities. At both reduced temperatures, it can be seen that the bubble departure diameter predicted by the BGK model agreed very well with the ones predicted by the MRT model. From the release period results shown in Figure 5, a reasonable agreement is also observed. The highest discrepancy was noticed at $T_r=0.76$ at the lowest gravitational acceleration, namely $g=3.125\times {10}^{-5}$. In this case, it can be noticed that the release period predicted by the BGK model is smaller than the one obtained by the MRT model. At $T_r=0.86$, a reasonable agreement is observed for the release period, despite the discrepancy observed also at $g=3.125\times {10}^{-5}$. In this case, the BGK model predicted a slightly higher release period. The difference between the results might suggest an influence of the forcing scheme, as showed by Li et al. [44]. The forcing scheme introduced by Gong and Cheng [12] introduces a modified velocity $v_{new}$ (see Equation (12)), in which the additional terms serve as an approximation for $\psi \nabla \psi$. Hence, at large density ratios, the coexistence densities will be affected by the kinematic viscosity and consequently by the single relaxation time. For the MRT model operator, the influence of the viscosity at large density ratios are reduced significantly, since only ${\tau }_v$ changes with viscosity (see Equations (17) and (18)). Hence, the BGK model with the forcing scheme proposed by Li et al. [27] may be not suitable for simulations at lower reduced temperatures. It should be noticed that the simulations with this operator were performed considering a variable relaxation time, since the kinematic viscosity varies in the computational domain.

Considering the correlations given by Equations (45) and (46), fitted curves are also shown for each model at a given temperature. For the departure diameter at $T_r=0.86$, the exponent matches exactly with the theoretical results from literature [13]. On the other hand, at $T_r=0.76$, the exponent is quite close to the one predicted by the correlation. For the release period, the obtained exponents are in reasonable agreement with the ones given by the empirical correlations proposed by Fritz [13] and Zuber [43], respectively.

Now, comparing the results between the two simulated reduced temperatures, it can be noticed that the bubble diameter is higher at small reduced temperature for all simulations. This is an expected behavior, because at lower saturation temperatures, the surface tension is higher. This also leads to higher release periods, because the required values of the buoyancy force are higher for the bubble detachment.

Some snapshots of bubble nucleation, growth and departure from the heated surface for different time-steps are presented for both reduced temperatures in Figure 6 ($T_r=0.76$) and Figure 7 ($T_r=0.86$). The gravitational acceleration was set to $g=1\times {10}^{-5}$ and $g=2.5\times {10}^{-5}$, for $T_r=0.76$ and $T_r=0.86$, respectively. At this point, it should be mentioned that the gravitational acceleration plays an important role on the stability of the simulation. A higher value of this quantity (e.g., $2.5\times {10}^{-5}$) at $T_r=0.76$ resulted in instabilities in the simulation, even for the MRT model. The effect of the gravitational acceleration over the stability of the simulation is not the focus of the present study and will be addressed in details in a future research work. The results shown in the Figure 6 clearly show that, at $T_r=0.76$, the BGK model predicted a smaller release period when compared with the MRT model. In Figure 7, it can be noticed that the release period is very similar for both models, as shown in Figure 5b.

Following the presented snapshots of the bubble cycle at $T_r=0.76$ and $T_r=0.86$, now the dynamics of the bubble diameter and the stages associated with the vapor bubble growth, namely expansion and rewetting, are investigated in detail. In Figure 8, the dynamics of the bubble diameter and the space-averaged heat flux are presented for both reduced temperatures. The space-averaged heat flux is computed on the fluid side, considering the length of the heater (see Figure 4).

These results are presented at the same conditions to the ones shown in Figure 7 using the MRT model. The stages can be identified by analyzing the local heat flux and the space-averaged heat flux, which can be computed by, respectively:

$$q^{{}^{″}}=-k{\left(\frac{\partial T}{\partial y}\right)}_{y=0}$$

(47)

$$\overline{q^{{}^{″}}}=\frac{{\int }_0^{L_x}{q}^{{}^{″}}dx}{L_x},$$

(48)

where $q^{{}^{″}}$ and $\overline{q^{{}^{″}}}$ are the local and space-averaged heat flux, respectively, $L_x$ is the length of the channel, and k is thermal conductivity.

As shown in Figure 8 for both reduced temperatures, the heat flux reduces during the expansion stage. In the rewetting stage, the liquid surrounding the bubble moves toward the heated surface, increasing the heat flux. During both stages, the bubble diameter increases at approximately the same rate. The peak heat flux is obtained at the time of the bubble departure. As expected the surface averaged, and consequently the local, heat fluxes are higher for the lower reduced temperature. This is related to the latent heat of vaporization, as well as the thermal conductivity of the liquid, which are higher at $T_r=0.76$.

4.3. Single Bubble Nucleation under Forced Convection

In this sub-section, simulations of the bubble cycle under forced convection conditions are performed at different Reynolds number. In this case, the bubble cycle is affected by a constant pressure gradient driven flow. This kind of problem is very important, since it is similar to the flow boiling inside a channel. Hence, the results can be used, at least in the qualitative point of view, to understand the phenomena that are associated with the real flow boiling problem. A schematic view of the simulated problem is shown in Figure 9. It should be mentioned that the problem was simulated using the same set-up as the pool boiling, with the assumption of fully-developed flow. This means that the domain is initialized with a liquid–vapor distribution with period boundary conditions at the inlet and outlet of the channel. At the walls, constant temperatures were considered.

First, the isothermal flow is simulated and compared with the dimensionless form of the well-known analytical solution for fully developed flow [45]. For simplicity, the relaxation time is taken as 1. The simulation was performed at $Re=100$. The computational mesh was 20 × 20. The Reynolds number was computed using the maximum flow velocity. The chosen convergence criteria was that velocity field remains unchanged under a ${10}^{-6}$ tolerance between two consecutive iterations. Both collision operators were considered in the simulations. The results of the dimensionless velocity field at the outlet are presented in Figure 10. It can be noticed that the LBM simulation with the models with the BGK and MRT collision operators correctly predicts this flow, exactly matching with the analytical solution.

Now, a small heater was considered at the position $x_{source}=50$. The length of the small heater was $L_{source}=3$, expressed in lattice units. A constant temperature was imposed to represent the heater (see Figure 9). The same parameters used in the simulations which provided the results shown in Figure 6 and Figure 7 are considered. In order to simulate a geometry more similar to a channel, the following computational mesh was considered: 300 × 151. The Reynolds number, defined as $Re=\frac{u_{max}{H}_c}{{\nu }_l}$, was used to characterize the flow, where $H_c$ is the channel height. Before the study of the bubble cycle under forced convection conditions, the isothermal flow was simulated until the equilibrium is reached. Then, the effect of the small heater was included in the simulation. The same parameters and thermophysical properties used in the previous section were considered in the simulations.

The results are presented only for the model using the MRT collision operator. The results for both reduced temperatures are presented together in order to better demonstrate the differences related to this property. The influence of the Reynolds number over the departure diameter and release period was explored.

In Figure 11, numerical results of the density field for the bubble cycle at $T_r=0.76$ and $T_r=0.86$ are presented for $Re=100$. It can be seen that the advection effects promote an earlier departure of the bubble from the heated surface. This is associated with the effect of the drag force. The temperature fields are shown in Figure 12, for both reduced temperatures. The temperature field is also affected by the advection influence of the main flow. In Figure 13, the results of the density field are presented for a higher Reynolds number, namely $Re=500$. As the Reynolds number is increased, the drag force resulting from the flow advection effects promote the bubble departure from the heated surface. Hence, a smaller departure diameter is noticed when the Reynolds number is increased. This is observed for both reduced temperatures. In Figure 14, the temperature fields are shown for the same ondition. In this case, the higher advection effects resulted in a lower growth rate of the bubble.

The results of the departure diameter and the release period for the first bubble are presented in Figure 15. As shown before, the departure diameter decreases with the increase of the Reynolds number. It can be noticed that the decrease rate is more pronounced at $T_r=0.76$. In fact, the departure diameter at $T_r=0.76$ is lower than the one obtained at $T_r=0.86$, when $Re=500$. From the results shown Figure 15a, it can be noticed that the bubble is actually removed earlier from the heated surface due to higher advection effects. The release period for the first bubble follow also reduces when the Reynolds number is increased. This is also due to the effect of the drag force. It should be noticed that, due to this effect, the released period is lower at $T_r=0.76$ than $T_r=0.86$. From the experimental point of view, the effect of the Reynolds number over the departure diameter and release period can be seen in Tibiriçá and Ribatski [46]. In their study, the authors presented an experimental analysis of the departure bubble diameters and release frequencies of $R134a$ and $R245fa$ refrigerants for flow boiling in microscale channels. The obtained experimental results show a decrease of the departure bubble diameter with the $Re$ augmentation for both refrigerants. This behavior is qualitatively the same shown in Figure 15a, for both reduced temperatures, confirming the observed trends of the present simulation results.

Similarly to the analysis of the results performed for the pool boiling simulations, the space-averaged heat flux is computed for each Reynolds number. The space-averaged heat flux is computed at each timestep, considering the full length of the heater (see Figure 9). The displayed results in Figure 16 confirm the results presented in Figure 15b. The heat flux peak is associated with the bubble release period. Hence, the release period of the subsequent bubbles can also be determined using these results. At $T_r=0.76$, it can be noticed that the heat flux peak occurs first for $Re=500$ than $Re=300$, and $Re=100$ as expected. The same behavior is noticed at $T_r=0.86$.

Another interesting observation is related to the time-averaged heat flux over the heated surface domain obtained for the two reduced temperatures as a function of Reynolds number, which are presented in Figure 17. The time-averaged heat flux is computed for the time interval show in Figure 16. At $T_r=0.86$, it can be noticed that the heat flux actually increases slightly as the Reynolds number increases. This can be explained by the effect of the drag force, with removes the bubble from the surface, allowing the cooled liquid to reach the heated surface. At $T_r=0.76$, the heat flux increased from $Re=100$ to $Re=300$ and decreased from $Re=300$ to $Re=500$. The last can be explained by the strong effect of the advection effects over the nucleation process, as also shown in Figure 16a. The heat flux peak value at $Re=500$ is lower than at $Re=300$. Hence, the computed time-averaged heat flux is lower, for the considered time interval.

5. Conclusions

In this paper, the bubble cycle was simulated using an improved pseudopotential Lattice Boltzmann Method from the literature at different saturation temperatures, namely $T_r=0.76$ and $T_r=0.86$, using two models which used the BGK and MRT collision operators. The implementation of these models considers the application of appropriate forcing schemes for each model in order to ensure thermodynamic consistency [27,28]. A detailed comparison of the numerical results under different conditions was performed and the following conclusions can be drawn:

• The comparison between the results for different reduced temperatures revealed that the decrease of the reduced temperature results in bubbles with higher departure diameter and higher release period for the pool boiling case. This behavior is in accordance with boiling heat transfer theory. In the case of flow boiling, the departure diameter, reduces as the Reynolds number is increased, indicating that this quantity is a strong function of the drag force.
• In the pool boiling simulations, at both reduced temperatures, $T_r=0.86$ and $T_r=0.76$, a reasonable agreement is observed for the departure diameter regarding the results obtained with both simulation models. The release period also showed good agreement, presenting some discrepancies for the lowest gravitational acceleration, namely $g=3.125\times {10}^{-5}$. This might suggest an influence of the forcing scheme used for the model with the BGK collision operator over the numerical results. This observation is based on the fact that models with the BGK collision operator use a single relaxation time related to the assumed kinematic viscosity. For the considered forcing scheme, this affects in more extend the coexistence densities at large densities ratios. It is a fact that even in single-phase flow simulations models with the BGK collision operator present less numerical stability than those with the models with the MRT collision operator [29].
• For pool boiling and flow boiling simulations, the space-averaged heat flux is an important quantity. for the pool boiling simulations, the expansion and rewetting stages can be identified. For the flow boiling, the release period of the departure bubbles can be identified by the heat flux peak. For flow boiling simulations at both reduced temperature, the periodic behavior of the heat flux is observed for all simulated Reynolds number. At $T_r=0.76$ and $T_r=0.86$, the effect of the Reynolds number is to anticipate the bubble departure.
• The pool boiling results showed that the gravitational acceleration plays a very important role on the LBM numerical stability for medium to lower reduced temperatures. The same behavior can be attributed to the flow boiling problem. This behavior is present even when the MRT model is used. In order to avoid these numerical instabilities, small values of gravitational acceleration should be used. This introduces a question regarding the simulation of real boiling process for medium to lower reduced temperatures ($T_r\le 0.8$).The use of very small values of gravitational acceleration can produce difficulties in fitting the LBM simulations into the desired limits of dimensionless numbers necessary to simulate a particular heat transfer phase-change process, i.e., the Grashof and Bond numbers. Further investigations are under development to address this issue.