Phase-field lattice Boltzmann modeling of boiling using a sharp-interface energy solver

The main objective of this paper is to extend an isothermal incompressible two-phase lattice Boltzmann equation method to model liquid-vapor phase change problems using a sharp-interface energy solver. Two discrete particle distribution functions, one for the continuity equation and the other for the pressure evolution and momentum equations, are considered in the current model. The sharp-interface macroscopic internal energy equation is discretized with an isotropic finite difference method to find temperature distribution in the system. The mass flow generated at liquid-vapor phase interface is embedded in the pressure evolution equation. The sharp-interface treatment of internal energy equation helps to find the interfacial mass flow rate accurately where no free parameter is needed in the calculations. The proposed model is verified against available theoretical solutions of the two-phase Stefan problem and the two-phase sucking interface problem, with which our simulation results are in good agreement. The liquid droplet evaporation in a superheated vapor, the vapor bubble growth in a superheated liquid, and the vapor bubble rising in a superheated liquid are analyzed and underlying physical characteristics are discussed in detail. The model is successfully tested for the liquid-vapor phase change with large density ratio up to 1000.

Figure 1

A schematic of a liquid-vapor system with the phase interface.

Figure 2

INB nodes (solid blue circles) surrounding the phase interface (open circles: vapor bulk nodes, solid red circles: liquid bulk nodes). For interpretation of the references to color in this figure caption, the reader is referred to the web version of this article.

Figure 3

Illustration of a phase interface curve in the regular lattice structure with normal vectors to phase interface (solid blue circles: INB nodes, open circles: vapor bulk nodes, solid red circles: liquid bulk nodes, open squares: phase interface points; see text for notations and definitions.). For interpretation of the references to color in this figure caption, the reader is referred to the web version of this article.

Figure 4

Temperature gradient calculation at phase interface point $P_I$ on the vapor side using the temperatures of the phase interface ($T_I$), the neighboring lattice nodes ($T_1^v$ and $T_2^v$) in the vapor phase, and the ghost node ($T_0^v$).

Figure 5

The mass flux at the phase interface point $P_I$ distributed by a linear weighting function on lattice nodes.

Figure 6

Illustration of different phase interface points on different characteristic directions starting from point $N^v$ in the vapor phase.

Figure 7

A schematic of the steady state heat conduction in a liquid-vapor system along with initial temperature ($T$) and heat flux ($q$) distribution.

Figure 8

(a) Comparison of the temperature distribution in the steady state heat conduction problem between the current result (symbol) and the theory (line) along with initial temperature profile (square symbol) and (b) the sharp profile of the temperature across the phase interface.

Figure 9

(a) Comparison of the temperature distribution at $t^{*}=1$ along with initial temperature profile (square symbol) and (b) comparison of the heat flux jump $〈q〉=q_v-q_l$ in the unsteady heat conduction problem between the current result (symbol) and the theory (line).

Figure 10

A schematic of the two-phase Stefan problem with initial temperature profile.

Figure 11

Comparison of the temperature distribution in the two-phase Stefan problem with small density effect (${\rho }_l/{\rho }_v=1.5$) between the LBE result (symbol) and the theory (line) for different thermal diffusivity ratios and (b) the sharp profile of the temperature across the phase interface.

Figure 12

Comparison of the phase interface location in the two-phase Stefan problem with small density effect (${\rho }_l/{\rho }_v=1.5$) for different thermal diffusivity ratios.

Figure 13

Comparison of (a) the phase interface location and (b) the liquid velocity ($u_l$) in the two-phase Stefan problem between the LBE result (symbols) and the theoretical solution (lines) [24] for different density ratios.

Figure 14

Comparison of (a) the non-dimensional temperature profiles at $t^{*}=0.5$ and (b) the temporal evolution of the heat flux jump in the two-phase Stefan problem between the LBE result (symbols) and the theoretical solution (lines) [24] for different density ratios.

Figure 15

Comparison of the phase interface location in the two-phase Stefan problem for different phase interface thicknesses at ${\rho }_l/{\rho }_v=100$.

Figure 16

Convergence test toward sharp-interface limit for the interface location based on (a) different initial vapor layer thicknesses $X\left(0\right)=10,20,30,$ and 40 and (b) different thermal diffusivities $\chi =0.001,0.002,0.004$ in the two-phase Stefan problem. For interpretation of the references to color in this figure caption, the reader is referred to the web version of this article.

Figure 17

A schematic of the two-phase sucking interface problem and the initial temperature profile.

Figure 18

Comparison of the phase interface location in the two-phase sucking interface problem between the LBE result (symbols) and the theoretical solution (lines) [13] for different density ratios.

Figure 19

The temperature profile at different times in the two-phase sucking interface problem at ${\rho }_l/{\rho }_v=100$.

Figure 20

The scaling relation between the calculated nondimensional phase interface displacement and Stefan number after a fixed number of iterations in the two-phase sucking interface problem.

Figure 21

(a) A schematic diagram of the droplet evaporation. (b) Transient radius of the droplet at $t^{*}=0$, 0.29, 0.67, and 1 and velocity vectors around the evaporating droplet.

Figure 22

(a) Mesh study in the droplet evaporation problem. (b) Transient radius of the droplet for five different initial radii.

Figure 23

Comparison of the Young-Laplace law for the evaporating droplet between the LBE result (symbol) and the theory (line) for two different rates of evaporation corresponding to ${\chi }_v=0.01$ and 0.05.

Figure 24

The investigation of the $d^2$ law in droplet evaporation problem for four different Ja numbers of $\text{Ja}=0.25$, 0.5, 0.75, and 1.0.

Figure 25

Effect of (a) surface tension, (b) relaxation parameter corresponding to three different Laplace numbers ($\text{La}=\sigma {\rho }_l{D}_0/{\mu }_l^2$) on transient radius of the droplet during evaporation.

Figure 26

A schematic of the vapor bubble growth problem in a superheated liquid and the temperature distribution in radial direction where the donut shape with light blue color represents the thermal boundary layer around the vapor bubble.

Figure 27

Transient shape of the vapor bubble at different simulation times (a) $t^{*}=0$, (b) $t^{*}=0.5$, and (c) $t^{*}=1.0$.

Figure 28

Comparison of the transient radius in the vapor bubble growth problem between the LBE result (symbols) and the theoretical solution (lines) [40].

Figure 29

Investigation of a mesh independent solution in a two-dimensional vapor bubble growth in a superheated liquid.

Figure 30

Radial distributions of temperature, velocity magnitude ($U$), density ($\rho$), and interfacial source term across the phase interface at $t^{*}=0.5$ in the vapor bubble growth problem.

Figure 31

A schematic of the three-dimensional vapor bubble rising in a superheated liquid problem.

Figure 32

Snapshots of the vapor bubble shapes at different simulation times $t^{*}=0$, 0.16, 0.28, 0.40, 0.52, 0.64, 0.76, 0.88, and 1 in vapor bubble rising in superheated liquid problem.

Figure 33

Temperature contours at the midplane in $z$ direction along with the vapor bubble transient shape at three different times $t^{*}=0.28$, 0.64, and 1 in vapor bubble rising in a superheated liquid problem. For interpretation of the references to color in this figure caption, the reader is referred to the web version of this article.

Figure 34

The interfacial source term contours at the midplane in $z$ direction at three different times $t^{*}=0$, 0.64, and 1 in vapor bubble rising in a superheated liquid problem. For interpretation of the references to color in this figure caption, the reader is referred to the web version of this article.