Thermal lattice Boltzmann method for multiphase flows

TL;DR

This paper introduces a thermal lattice Boltzmann method for simulating heat transfer in multiphase flows, effectively handling phase boundaries and latent heat without interface tracking, with validated invariance and conduction properties.

Contribution

The paper presents a novel LB-based approach that incorporates internal energy via additional distribution functions and suppresses parasitic heat diffusion near boundaries.

Findings

Method maintains Galilean invariance.

Correctly models sound speed dependence on heat capacity.

Low scheme diffusion of internal energy.

Abstract

New method to simulate heat transport in multiphase lattice Boltzmann (LB) method is proposed. The energy transport equation needs to be solved when phase boundaries are present. Internal energy is represented by an additional set of distribution functions, which evolve according to a LB-like equation simulating the transport of a passive scalar. Parasitic heat diffusion near boundaries with large density gradient is supressed by using the interparticle "pseudoforces" which prevent the spreading of energy. The compression work and heat diffusion are calculated by finite differences. The latent heat of a phase transition is released or absorbed in the inner side of a thin transition layer between liquid and vapor. This allows one to avoide the interface tracking. Several tests were carried out concerning all aspects of the processes. It was shown that the Galilean invariance and the scaling of thermal conduction process hold as well as the correct dependence of sound speed on the heat capacity ratio. The method proposed has low scheme diffusion of the internal energy, and it can be applied to modeling a wide range of multiphase flows with heat and mass transfer.

Equations

$$N_k(\mathbf{x}+{\mathbf{c}}_k\Delta t,t+\Delta t)=N_k(\mathbf{x},t)+{\Omega }_k(N)+\Delta N_k,k=0,\dots ,b,$$

$${\Omega }_k=\frac{N_k^{eq}(\rho ,\mathbf{u})-N_k(\mathbf{x},t)}{\tau },$$

$$N_k^{eq}(\rho ,\mathbf{u})=\rho w_k\left(1+\frac{({\mathbf{c}}_k\mathbf{u})}{\theta }+\frac{({\mathbf{c}}_k\mathbf{u}{)}^2}{2{\theta }^2}-\frac{u^2}{2\theta }\right).$$

$$\Delta N_k(\mathbf{x},t)=N_k^{eq}(\rho ,\mathbf{u}+\Delta \mathbf{u})-N_k^{eq}(\rho ,\mathbf{u}),$$

$$\rho =\sum _{k=0}^b{N}_k,\rho \mathbf{u}=\sum _{k=1}^b{\mathbf{c}}_k{N}_k.$$

$$\rho {\mathbf{u}}^{\ast }=\sum _{k=1}^b{\mathbf{c}}_k{N}_k+\mathbf{F}\Delta t/2.$$

$$\mathbf{F}=-\nabla U=-\nabla (p(\rho ,T)-\rho \theta ).$$

$${\mathbf{F}}_N=2A\nabla ({\Phi }^2)+(1-2A)2\Phi \nabla \Phi ,$$

$$\mathbf{F}(\mathbf{x})=\frac{1}{\alpha h}\left[A\sum _{k-1}^b{G}_k{\Phi }^2(\mathbf{x}+{\mathbf{e}}_k){\mathbf{e}}_k+(1-2A)\Phi (\mathbf{x})\sum _{k-1}^b{G}_k\Phi (\mathbf{x}+{\mathbf{e}}_k){\mathbf{e}}_k\right],$$

$$\tilde{p}=\frac{8\tilde{\rho }\tilde{T}}{3-\tilde{\rho }}-3{\tilde{\rho }}^2.$$

$${\left(\frac{\partial p}{\partial \rho }\right)}_T\le {\left(\frac{h}{\Delta t}\right)}^2+\theta .$$

$$\frac{\partial E}{\partial t}+\nabla \cdot (\mathbf{u}E)=\frac{p}{\rho }\frac{d\rho }{dt}+\nabla \cdot (\lambda \nabla T)+\hat{\sigma }:\nabla \mathbf{u},$$

$$\frac{p}{\rho }\frac{d\rho }{dt}=-p\mathrm{div}({\mathbf{u}}^{\ast }).$$

$$g_k(\mathbf{x}+{\mathbf{c}}_k\Delta t,t+\Delta t)=g_k(\mathbf{x},t)+\frac{g_k^{eq}(E,\mathbf{u})-g_k(\mathbf{x},t)}{{\tau }_E}+\Delta g_k(\mathbf{x},t).$$

$$\Delta g_k^{(1)}(\mathbf{x},t)=g_k(\mathbf{x},t)\frac{\Delta E}{E}.$$

$$\tilde{p}=\frac{8\tilde{\rho }{\tilde{T}}_0}{3-\tilde{\rho }}-3{\tilde{\rho }}^2.$$

$${\left(\frac{\partial E}{\partial V}\right)}_T=P{\left(\frac{\partial P}{\partial T}\right)}_V-P=0.$$

$$\Delta g_k^{(2)}(\mathbf{x},t)=g_k^{eq}(E,\mathbf{u}+\Delta \mathbf{u})-g_k^{eq}(E,\mathbf{u}).$$

$${\lambda }_{liq}{\frac{\partial T}{\partial x}|}_{x=\xi -0}-{\lambda }_{vap}{\frac{\partial T}{\partial x}|}_{x=\xi +0}={\rho }_{liq}(T)Q(T)\frac{d\xi }{dt},$$

$$\frac{dE}{dt}=\frac{{\rho }_{liq}Q(T)}{{\rho }_2-{\rho }_1}\frac{d\rho }{dt}=-\frac{{\rho }_{liq}Q(T)}{{\rho }_2-{\rho }_1}\rho \mathrm{div}({\mathbf{u}}^{\ast }).$$

$$c_T={\left(\frac{\partial \tilde{p}}{\partial \tilde{\rho }}\right)}_T=\frac{24{\tilde{T}}_0}{(3-\tilde{\rho }{)}^2}-6\tilde{\rho },c_S={\left(\frac{\partial \tilde{p}}{\partial \tilde{\rho }}\right)}_S=\frac{24\gamma {\tilde{T}}_0}{(3-\tilde{\rho }{)}^2}-6\tilde{\rho }.$$

$$\Delta E_i^n=-p_i^n\frac{(u^{\ast }{)}_{i+1}^n-(u^{\ast }{)}_{i-1}^n}{2h}\Delta t.$$

$$T_{vap}=T_0-\frac{p}{C_V}\left(\frac{1}{{\rho }_{vap}}-\frac{1}{{\rho }_0}\right).$$

Full text

1 Introduction

Simulation of fluid flows with phase transitions between liquid and vapor is difficult because new phase boundaries can apper in the bulk during calculations, and the existing boundaries can disappear or change their topology. Therefore, the application of interface tracking methods is difficult if not impossible. Moreover, the density ratio of liquid and vapor phases is usually high (can reach tens and hundreds of thousands) leading to noticeable numerical diffusion and/or dispesion near the boundaries when using common finite-difference methods.

The lattice Boltzmann method (LBM) [1, 2] is based on the solution of a kinetic equation for pseudoparticles. It was widely applied for simulating flows of single-phase and multiphase media [3, 4, 5, 6, 7, 8, 9]. Moreover, the method is easily parallelizable on graphic accelerators using the GUDA technology [8, 9, 10, 11, 12, 13].

Three essentially different approaches were proposed to simulate the heat transport in the LBM: the model with the extended velocity set [14, 15, 16, 17], the use of the finite-difference equation for energy [18], and the method of a second set of distribution functions (method of passive scalar) [19, 20].

The first approach has rather narrow range of simulated temperature where simulations are stable. Moreover, the size of data increases significantly. This method was applied only to single-phase flows. When the energy equation is solved by finite difference methods, the large numeric diffusion arises in the case of moving fluid that restricts significantly the possibilty of modeling.

The passive scalar approach is realized in the LBM by introducing the additional set of distribution functions. This model has much lower numeric diffusion than the finite-difference method. Usually the temperature $T$ is used as a passive scalar which can be done only in a case of almost constant density of the fluid. When a phase transition liquid-vapor is present, the change of density is not small, and the tarnsport of internal energy should be considered instead of temperature. Such approach was used in the works [20, 21], and it was combined with extended velocity set in the work [22]. However, only single-phase flows were considered in these works.

2 Lattice Boltzmann method

The lattice Boltzmann method (LBM) is based on the solution of a kinetic equation for pseudoparticles. Only a limited set of particle velocities ${\bf c}_{k}$ is possible such that the vectors ${\bf e}_{k}={\bf c}_{k}\Delta t$ correspond to vectors to neighbor nodes of a regular spatial lattice [23]. The usual choice of velocity set is three vectors in one-dimensional case $|{\bf c}_{k}|=0,h/\Delta t$, D1Q3 model, nine velocities in 2D case (D2Q9, $|{\bf c}_{k}|=0,h/\Delta t,\sqrt{2}h/\Delta t$), and nineteen velocities in 3D case (D3Q19, $|{\bf c}_{k}|=0,h/\Delta t,\sqrt{2}h/\Delta t$). Here, $h$ is the lattice space, $\Delta t$ is time step.

One-particle distribution functions $N_{k}$ are used as main variables, they have the meaning of the parts of fluid density. Evolution equation for $N_{k}$ has the form

[TABLE]

where $\Omega_{k}$ is the collision operator, and $\Delta N_{k}$ is the change of distribution functions under the action of volume forces (both external and internal).

The collision operator is mostly chosen in the form of a relaxation to local equilibrium with one (BGK, Bhatnagar-Gross-Krook model [23]) or several relaxation times (MRT, Multi-Relaxation-Time model [24]). For the BGK model, the collision operator is

[TABLE]

where $\tau=t_{rel}/\Delta t$ is the non-dimensional relaxation time. Equilibrium distribution functions are usually taken as truncated Maxwellians up to the second order in fluid velocity $u$ [25]

[TABLE]

Coefficients $w_{k}$ depend on the lattice geometry [23]. For the one-dimensional model D1Q3 they are $w_{0}=2/3,w_{1-2}=1/6$, for two-dimensional model D2Q9 $w_{0}=4/9,w_{1-4}=1/9,w_{5-9}=1/36$, and for three-dimensional model D3Q19 $w_{0}=1/3,w_{21-6}=1/18,w_{7-18}=1/36$. The kinetic temperature of pseudoparticles in LB models listed is $\theta=(h/\Delta t)^{2}/3$, and the kinematic viscosity is defined by the relaxation time $\nu=\theta(\tau-1/2)/\Delta t$. Change of distribution functions at a node is calculated using the Exact Difference Method (EDM) [26]

[TABLE]

where $\Delta{\bf u}={\bf F}\Delta t/\rho$ is the change of fluid velocity in one time step, and F is the total force acting on the fluid at a node.

Density $\rho$ and velocity u of fluid are calculated as

[TABLE]

Under the action of volume forces, the physical fluid velocity ${\bf u}^{*}$ is defined at half time step [27]

[TABLE]

3 Phase transitions

In order to obtain phase transitions in a fluid, it is necessary to model in the LBM the attractive part of the “intermolecular” interaction. This was done in the work [5] by the introduction of attractive forces acting on the fluid at a node from neighbor nodes. Later, the total force F acting on the fluid at a node was introduced as a gradient of pseudopotential $U$ defined from the equation of state for the fluid [6, 18]

[TABLE]

LBM with such attractive forces represents phase boundary as a thin transition layer between liquid and vapor where density changes smoothly across several lattice nodes (interface capturing). The surface tension arises at phase boundaries.

In the work [28], a new function was introduced $\Phi=\sqrt{-U}$. The equation (6) can be rewritten in the equivalent form

[TABLE]

where $A$ is a free parameter which allows one to obtain correct phase densities at the coexistence curve. The sufficiently isotropic approximation of the formula (7) is

[TABLE]

where coefficients $G_{k}>0$ differ for different lattice directions. For neighbor nodes, they are $G_{k}=1$. For next-neighbor nodes, the vaues of coefficient ensuring isotropy are $G_{5-8}=1/4$ in the two-dimensional model D2Q9 and $G_{7-18}=1/2$ in the three-dimensional model D3Q19. The coefficient $\alpha$ is equal to 1, $3/2$, and 3 for the models D1Q3, D2Q9, and D3Q19, correspondingly.

In present work, we used Van der Waals equation of state which is written in reduced variables as

[TABLE]

Here and below, pressure, density and temperature are scaled by their values at the critical point, $\tilde{p}=p/p_{c},\tilde{\rho}=\rho/\rho_{c},\tilde{T}=T/T_{c}$. For this equation, the approximation (7) gives the best agreement with the phase coexistence curve at $A=-0.152$ (the deviation of density in simulations from the theoretical value is less than 0.2% in the range of temperature from the critical one $\tilde{T}=1$ down to $\tilde{T}=0.4$ [28]). More complex equations of state including tabulated ones for real fluids were considered in [29, 30].

The stability of the LBM with equation of state in the form $p=p(\rho,T)$ is defined by the criterion [31]

[TABLE]

4 Heat transport

Evolution equation for the internal energy per unit volume $E$ is

[TABLE]

where the first term in the right hand side corresponds to the pressure work, the second one represents the heat conduction, and the last one is the viscous heating. Here, $\lambda=\rho C_{V}\chi$ is the heat conductivity, $C_{V}$ is the specific heat at constant volume, $\chi$ is the thermal diffusivity, and $\hat{\sigma}$ is the viscous stress tensor. It is more convenient to express the pressure work from the velocity divergence using the continuity equation

[TABLE]

Calculation of the left hand side of Eq. (8) is the most complicated, i.e., the advection of the internal energy by the fluid flow with velocity calculated from Eq. (5). The viscous heating is usually small, and we will neglect it in the following.

This was done in the works [14, 15, 16, 17] by introducing the extended set of velocity vectors ${\bf c}_{k}$ and the increasing the expension order of euilibrium distribution functions (terms up to fourth order in u were used). The drawbacks of this approach is relatively narrow range of simulated temperature where simulations are stable, and the significant increase of the size of data.

In the work [18], the advection of energy in equation (8) was calculated by a finite-difference method using the values of fluid density and velocity obtained from LBE. However, this method produces large numerical diffusion and dispersion of energy near phase boundaries in simulations of moving fluids which significantly complicates modeling.

Third approach to simulating the advection of energy in LBM is based on the use of passive scalar (additional set of distribution functions $g_{k}$) [19] which has much lower scheme diffusion comparing to finite-difference methods. Earlier, this approch was used for simulating the flows with almost constant fluid density and specific heat when temperature can be used as passive scalar. At phase transitions between liquid and vapor, the change of density is however not small, and the advection of internal energy should be considered instead of temperature.

When using the passive scalar approach $E=\sum\limits_{k=0}^{b}g_{k}$ for the internal energy density $E=\rho c_{V}T$, evolution equations for distribution functions $g_{k}$ can be written in the form analogous to Eq. (1)

[TABLE]

Here, $\Delta g_{k}=\Delta g_{k}^{(1)}+\Delta g_{k}^{(2)}$ is the total change of distribution functions, $\tau_{E}$ is the non-dimensional relaxation time for the energy. Equilibrium distribution functions $g_{k}^{eq}(E,{\bf u})$ have the same form as $N_{k}^{eq}(E,{\bf u})$ (2).

The change of energy at a node $\Delta E$ due to the pressure work and the heat conduction is calculated by usual finite-difference formulas. Corresponding changes of energy distribution functions $\Delta g_{k}^{(1)}$ are proportional to the change of energy

[TABLE]

The main problem in this approach stems from the jump of specific heat at phase boundaries. This leads to parasitic diffusion (spreading) of internal energy from dense phase (liquid) to rarefied one (vapor) even if pressure and temperature are uniform. This effect is readily observed for stationary droplet in the case of barotropic equation of state (pressure depends only on density). Since there is no feedback on temperature, waves of pressure and density does not arise. To demonstrate the parasitic diffusion of energy, the barotropic van der Waals equation of state was used with a constant $\tilde{T}_{0}$ instead of the temperature of fluid

[TABLE]

Parasitic diffusion at phase boundaries is shown in Fig. 1. The pressure work and heat diffusion were switched off for clarity. Internal energy “leaks” from the stationary liquid droplet to the surrounding saturated vapor. This leads to generation of non-physical temperature peaks in vapor and drops in liquid near the boundaries. In thermal simulations, such peaks and drops will lead to instability.

In this work, we modify the passive scalar approach in order to be used for the advection of internal energy. The idea is to introduce special “pseudoforces” for energy scalar which prevent the spreading at phase boundaries. The currently realized variant works in the case of a constant specific heat of fluid $C_{V}$, and the internal energy at given temperature proportional to the fluid density. This is valid for van der Waals and other linear in temperature equations of state since for them

[TABLE]

For water, the liquid and vapor specific heat along the coexistence curve are also close in a certain range of pressure (see Fig. 2). “Pseudoforces” are taken into account in the evolution equation for distribution functions (9) by the EDM similar to Eq. (3)

[TABLE]

Here u is the fluid velocity defined by the main set of distribution functions (4).

5 Latent heat of phase transition

It is known that the latent heat of phase transition should be taken into account at the conditions at a moving phase boundary. Corresponding boundary conditions are

[TABLE]

where $\xi$ is the coordinate of a planar pahse boundary between liquid and vapor, and $\rho_{liq}(T)$ is the liquid density at the phase coexistence curve (see Fig. 3,b). The latent heat of phase transition $Q(T)$ decreases with the increasing temperature, and comes to zero at the critical temperature $T=T_{c}$. Tracking the phase boundaries is difficult in simulations because in many cases new boundaries can appear, and existing ones can disappear or change their topology. The advantage of the LBM is its capturing of interfaces. Phase boundaries in the LBM are represented as transition layers where the density continuously change from the liquid to the vapor values according to the phase coexistence curve. The density of a portion of fluid at a phase transition also changes continuously in time. We propose to take into account the latent heat in the following way. If we do not resolve exactly the inner structure of the transition layer (see Fig. 3a) but take into account the latent heat of phase transition only integrally accross the transition layer, we can assume that the latent heat is released or absorbed continuously inside the transition layer in a certain range of density $\rho_{1}<\rho<\rho_{2}$ (Fig.  3b) according to the equation

[TABLE]

Equilibrium densities of the vapor and the liquid $\rho_{liq}$ at every temperature can be used as $\rho_{1}$ and $\rho_{2}$, correspondingly.

6 Numerical validation

6.1 Galilean invariance and scheme diffusion of energy

The initial state for the tests of Galilean invariance and scheme diffusion is shown in Fig. 4. The temperature and the density were distributed stepwise, and the pressure was constant. Periodic boundary conditions were used. The coefficient of scheme diffusion for energy was $D_{E}=\theta(\tau_{E}-1/2)\Delta t=0.001h^{2}/\Delta t$ with $\tau_{E}=0.503$.

Figure 5 shows distribution of temperature and density at diffrent time in the case of zero flow velocity (a–c), and for uniform flow velocity equal to $u=0.1h/\Delta t$ (d–f). Results are independent on the flow velocity, hence, the Galilean invariance holds.

Figure 6 shows the distribution of temperature in a resting fluid for two different thermal diffusivity and corresponding time. The similarity relation $l\sim\sqrt{\chi t}$ is fulfilled.

Since an explicit numerical scheme was used for calculating the heat diffusion, the stability criterion is $\chi\Delta t/h^{2}<0.5/d$ where $d$ is the number of spatial dimensions. One-dimensional calculations with the flow velocity equal to $u=0.1h/\Delta t$ and $\chi\Delta t/h^{2}=0.49$ were indeed stable.

Two-dimensional simulations were carried out using the D2Q9 model. A round droplet surrounded by a saturated vapor moved with uniform velocity along the diagonal of the simulation region. Periodic boundary conditions were used for both $x$ and $y$ directions. The initial temperature was constant, hence, the density of internal energy was higher inside the droplet. During the simulation time of $t=62000$, the droplet made more than six revolutions together with the flow which corresponds to 27 droplet diameters. The isotropy is preserved (the droplet remained round), and almost no parasitic diffusion of energy was present. Three-dimensional simulations with D3Q19 model gave similar results.

6.2 Pressure work

In order to check the calculations of pressure work, we investigated the dependence of the speed of sound on the specific heat. If pressure work is switched off, one obtains the isothermal speed of sound $c_{T}$, and taking into account the pressure work gives the adiabatic speed of sound $c_{S}$. For Van der Waals equation of state, the reduced values of both speeds at the temperature $\tilde{T}_{0}$ are

[TABLE]

Here, $\gamma=C_{P}/C_{V}$ is the heat capacity ratio. The speed of sound was calculated from the dispersion relation for a standing harmonic wave, $c=\omega L/2\pi$ where $L$ is the wavelength, and $\omega$ is the frequency. Figure 8 shows the dependence of the speed of sound on the inverse specific heat $C_{V}$. The isothermal speed of sound is constant, and the adiabatic speed of sound depends linearly on $1/C_{V}$ with agreement of the theoretical result (10).

Another test was the simulation of a spinodal decomposition (decay of an initially uniform fluid which is in thermodynamic state below the spinodal into a mixture of liquid and vapor). Figure 9 shows the simulation results for the case of neglected pressure work (Fig. 9a), and the case of pressure work taken into account (Fig. 9b). The pressure work was calculated with a finite-difference expression

[TABLE]

Thermal diffusivity was set to zero, and only small scheme diffusion was present. Without the pressure work, temperature remains constant. The pressure work results in the increase of the internal energy in liquid phase (which is compressed) and decrease of the internal energy in gas phase (where rarefaction occurs). Since the compression of liquid is relatively low, and the specific heat is significantly larger than that of the vapor, the temperature of the liquid increases only slightly. In contrast, the temperature of the gas phase decreases significantly. One can estimate the change of the vapor temperature as

[TABLE]

From Fig. 9 one can see that simulations give close value of the vapor temperature (Fig. 9b).

6.3 Latent heat of phase transition

The case of spinodal decomposition was simulated with zero and non-zero latent heat of phase transition. Results are shown in Fig. 10. When the latent heat is non-zero, the temperature of liquid phase is significantly higher than the initial one due to the release of the latent heat at the condensation of vapor.

The two-dimensional spinodal decomposition was simulated taking into account the pressure work and the latent heat of phase transition. The initial uniform fluid density was $\tilde{\rho}_{0}=1$, and the temperature was everywhere $\tilde{T}_{0}=0.8$. Simulation results are shown in Fig. 12. The temperature of vapor decresed to $\tilde{T}_{v}\approx 0.77$ which is lower than the initial one, and the temperature of liquid reached the value $\tilde{T}_{l}\approx 0.83$ which is higher than the initial one due to the release of latent heat at the condensation.

In further evolution of the system, the number of droplets decreases due to the coalescence and the evaporation of smaller droplets and growth of larger ones. The distribution of tempersture tends to the uniform one due to heat conductivity. Figure 12 shows the stage when only two droplets remain. The temperature of smaller droplet is lower due to the evaporation, the larger droplet is heated due to the condensation. The temperature difference is $\Delta\tilde{T}\approx 0.01$. When the process ends, and only one droplet remains, the nonuniformity of the temperature decreases to $\Delta\tilde{T}<0.001$.

7 Conclusion

The method of additional LB component was developed for multiphase thermal flows. The algorithm takes into accont heat conduction, pressure work and latent heat of phase transition. The method is interface-capturing, no tracking of phase boundaries and conditions at them is needed. Numerical tests shows Galelean invariance of the method, low scheme diffusion of energy, isotropy and stability. Calculated behavior of sound speed agrees well with theoretical predictions. The method developed is applicable for simulating flows with heat and mass transport and phase transitions.

Acknowledgements

This work was supported by the Russian Scientific Foundation, Grant N 16-19-10229.