[1]

[cor1]Corresponding author

# Pseudopotential Lattice Boltzmann Method for boiling heat transfer: a mesh
refinement procedure

Alfredo Jaramillo<EMAIL_ADDRESS>Vinícius Pessoa Mapelli
<EMAIL_ADDRESS>Department of Mechanical Engineering, Heat Transfer
Research Group, São Carlos School of Engineering (EESC), University of São
Paulo (USP), São Carlos 13566-590, Brazil  Luben Cabezas-Gómez
<EMAIL_ADDRESS>
###### Abstract

Boiling is a complex phenomenon where different non-linear physical
interactions take place and for which the quantitative modeling of the
mechanism involved is not fully developed yet. In the last years, many works
have been published focusing on the numerical analysis of this problem.
However, a lack of numerical works assessing quantitatively the sensitivity of
these numerical simulations to grid parameters can be identified, especially
for the Lattice Boltzmann method (LBM). The main goal of this work is to
propose a mesh refinement methodology for simulating phase-change heat
transfer problems by means of the pseudopotential LBM. This methodology was
based on relating the physical parameters to their lattice counterparts for an
arbitrary mesh under the viscous regime (where $\Delta t\propto\Delta x^{2}$).
A suitable modification of the EOS parameters and the adjusting of
thermodynamic consistency and surface tension for a certain $\Delta x$ were
the main steps of the proposed methodology. A first ensemble of simple
simulations including the droplet vaporization and the Stefan problems was
performed to validate the proposed method and to assess the influence of some
physical mechanisms. Global norms in space and time were used to evaluate the
variations of both the density and temperature fields for pool boiling
simulations when the lattice discretization is refined. It was observed that
the proposed methodology provides convergent results for all the problems
considered, and the convergence orders depend on the complexity of the
simulated phenomena.

###### keywords:

Pseudopotential LBMBoiling Heat transfer Lattice Boltzmann Equation Mesh
refinement

comment

##  1 Introduction

Boiling heat transfer is a phenomenon widely employed on common daily tasks as
long as on industrial applications, such as nuclear reactors, heat pipes,
industrial boilers, and so on (Carey, 2008). The reason for its widespread
adoption is its high efficiency in transferring heat, especially on nucleate
boiling regime, when vapor bubbles cycle through well-known phases of
formation, growth, and departure. The mathematical modeling and numerical
analysis for this phenomenon represent a challenging problem with inherent
non-linear behavior, and physical phenomena evolving along with different
scales (Son et al., 1999, 2002). Moreover, solid-liquid-vapor interparticle
interaction may take place, which has been proved to have an important
influence in quantities of interest like the heat transfer coefficient (Wang
and Dhir, 1993). For that reason, several researchers have dedicated
themselves to better understanding these phenomena, either experimentally,
theoretically, or numerically (Lienhard et al., 1973; Haramura and Katto,
1983; Kandlikar, 2001; Guan et al., 2011; Gong and Cheng, 2015a; Kim et al.,
2016; Liang and Mudawar, 2018; Fariñas Alvariño et al., 2019; Guzella et al.,
2020).

Recently, a numerical method that has drawn attention is the Lattice Boltzmann
Method (LBM) due to its potential to simulate a variety of complex phenomena
such as flow in porous medium (Guo and Zhao, 2002; Liu et al., 2016),
turbulence (Yu et al., 2005; Chen et al., 2004), particle suspensions
(McCullough et al., 2016; Safaei et al., 2016) and others (Gong and Cheng,
2012, 2013, 2015b, 2015a; Guzella et al., 2020). One of the LBM advantages for
simulating interfacial dynamics lies in the fact that it does not involve
front-tracking, as is the case on more traditional methods such as VOF and
level-set. LBM can be seen as the discretization of the Boltzmann transport
equation so that its relation to kinetic theory is also considered to be an
advantage, since it allows to include small and mesoscale particle physical
interactions into the mathematical modeling, mimicking the macroscopic flow
behavior.

Several LBM models have been developed for multiphase flow simulations
(Gunstensen et al., 1991; Swift et al., 1995; Luo, 1998; Shan and Chen, 1993).
Among them, pseudopotential-based models have been the most successful
methodologies for simulating phase-change heat transfer phenomena. The
original pseudopotential model was proposed by Shan and Chen (1993), where an
interparticle force is responsible for maintaining liquid-vapor phases in
equilibrium. Since then, many research groups have contributed to its
development. Yuan and Schaefer (2006) proposed a pseudopotential method in
which a non-ideal equation of states (EOS) could be employed, such as the Van
Der Waals equation and the Peng-Robinson equation. Sbragaglia et al. (2007)
was one of the first to propose a model that allowed to tune the surface
tension, a modeling feature that has a key relevance in this work. Kupershtokh
et al. (2009) proposed a modified interparticle force by combining the
original Shan and Chen method with the one proposed by Zhang and Chen (2003),
which resulted in an enhancement of thermodynamic consistency. In the works of
Li and Luo (2013) and Li et al. (2013), the authors proposed a consistent
approach to adjust surface tension and vapor-liquid coexistence curve. Later,
Lycett-Brown and Luo (2015) performed a third-order analysis and argued that
latter adjustments, along with an interface thickness treatment, are required
to achieve a consistent grid refinement when the pseudopotential method is
employed. More recently, Zhai et al. (2017); Kharmiani et al. (2019) and
Czelusniak et al. (2020) also proposed consistent models to tune these three
aspects of the pseudopotential method.

The pseudopotential method has also been employed in simulations involving
phase-change and heat transfer. Gong and Cheng (2012) proposed a similar
interparticle force to Kupershtokh et al. (2009) looking to improve
thermodynamical consistency. In a further study, Gong and Cheng (2013)
employed a second distribution function to simulate a temperature passive-
scalar equation along the aforementioned interparticle force to study
nucleation on a microheater under constant wall temperature and constant heat
flux conditions. Bubble growth and other aspects such as the effects of
contact angle, superheat, and gravity were investigated. A reasonable
agreement between the simulations results and well-known engineering
correlations was observed. Later, Gong and Cheng (2015a) used that same
methodology to investigate pool boiling on mixed wettability surfaces. In
their study, it was shown that surfaces with hydrophobic and hydrophilic spots
can achieve better heat transfer efficiency by avoiding bubble coalescence.
This improvement is also determined by geometry, such as the size and pitch of
these spots. In the same year, Li et al. (2015) employed a model they proposed
before (Li and Luo, 2013) along with a finite difference method for simulating
the temperature field evolution (called _the hybrid approach_ , firstly
performed by Zhang and Chen (2003)). By increasing the wall superheating, the
authors were able to simulate the regimes of nucleation, transition, and film
boiling. Later, Gong and Cheng (2017) included some solid nodes on their
computational grid to study the heater thermal response under different
boiling regimes. A similar configuration has been employed by Dong et al.
(2018) to investigate film boiling behavior. Mixed wettability surfaces have
also been investigated by Ma and Cheng (2019), with a three-dimensional
extension of the model proposed by Gong and Cheng (2012). More recently, Zhao
et al. (2021) employed a hybrid approach based on Li and Luo (2013) to study
several aspects of pool boiling on confined and free outlet geometries.

Despite its advantages when applied to phase-change simulations, LBM
pseudopotential-based methods could still be considered to be in their early
research stage. A great number of works present qualitative results and/or
values in the so-called lattice units, which are relative to mesh spacing and
sound speed. Even though the conversion to physical units is immediate, their
relationship may lead to confusion when a non-ideal EOS is employed. Also,
studies showing fundamental aspects of the pseudopotential method are still
scarce. For instance, in the work of Li et al. (2018a), the authors clarify a
misconception on the energy equation employed by Gong et al. (2018), which has
been extensively employed in previous studies. Also, Lycett-Brown and Luo
(2015) and Zheng et al. (2019) show some spurious terms on the macroscopic
equations when the forcing scheme proposed by Kupershtokh et al. (2009) is
used, which can also influence the coexistence curve.

Grid refinement when using the pseudopotential method is not a trivial task
(Lycett-Brown and Luo, 2015), and studies of the grid influence on numerical
pool boiling results are scarce. The present work focuses on fundamental
numerical aspects of the pseudopotential method when a hybrid approach is
used. Specifically, the main goal is to provide a mesh refinement procedure
for simulating liquid/vapor phase-change when the pseudopotential-based LBM is
employed. This kind of study is essential to perform LBM simulations of
multiphase problems for fluids of engineering interest. To attain this goal,
the physical parameters were related to their lattice counterparts for an
arbitrary mesh under the viscous regime (where $\Delta t\propto\Delta x^{2}$)
in order to conserve the hydrodynamic and thermodynamic behaviors. The non-
trivial part of this task is the suitable modification of the EOS and the
adjusting of the surface tension for a certain $\Delta x$. Once these
modifications were understood, some simple numerical experiments (the droplet
vaporization problem and the Stefan problem) were analyzed numerically to
validate the proposed methodology. Lastly, typical pool boiling simulations
are analyzed and discussed qualitatively and quantitatively.

The structure of this article is as follows: Section 2 describes the
pseudopotential LBM multiphase model, the contact angle models, and the energy
conservation equation employed. In Section 3, a mesh refinement procedure to
maintain a consistent physical behavior of the fluid between meshes is
presented. The grid refinement procedure is validated for the droplet cooling
and the Stefan problems in Section 4. Pool boiling numerical findings are also
shown and discussed in that section. Finally, some concluding remarks are
drawn in Section 5.

##  2 Mathematical modelling

The LBM is based on a discretized form of the Boltzmann transport equation
known as the lattice Boltzmann equation (LBE) (Chen and Doolen, 1998; Shan and
He, 1998). The LBE is employed to compute the fluid behavior calculating the
evolution of a series of discrete distribution functions $\\{f_{i}\\}$ given
by:

$f_{i}(\mathbf{x}+\,\mathbf{c}_{i},t^{n+1})=f_{i}(\mathbf{x},t^{n})+\Omega_{i}(\mathbf{f},\mathbf{f}^{\text{eq}})(\mathbf{x},t^{n})+\mathbf{R}(\mathbf{x},t^{n})$
(1)

where ${f}_{\alpha}$ are the particle distribution functions related to
velocity $\mathbf{c}_{i}$ and $f^{\text{eq}}_{i}$ are the local equilibrium
distribution functions, $\mathbf{f}$, $\mathbf{f}^{\text{eq}}$ correspond to
vectors whose components are $[\mathbf{f}]_{i}$ = $f_{i}$,
$[\mathbf{f}^{\text{eq}}]_{i}=f^{\text{eq}}_{i}$, and $t$ and $\mathbf{x}$ are
time and space coordinates, respectively. The term $\mathbf{R}$ is a general
representation of a source term in the moment space, defined by the forcing
scheme employed. Through this term, effects of a general external force field
$\mathbf{F}$ are taken into account. The term expressed by
$\Omega_{i}(\mathbf{f},\mathbf{f}^{\text{eq}})$ is the so-called collision
operator, and it can generally be defined as a quantity dependent on
distribution function $\mathbf{f}$ and its equilibrium values
$\mathbf{f}^{\text{eq}}$. For the two-dimensional D2Q9 lattice scheme, a
finite set of velocities vector $\mathbf{c}_{i}$ is formulated as (Qian et
al., 1992):

$\mathbf{c}_{i}=\left\\{\begin{array}[]{ll}(0,0)\,c,i=0~{},\\\
(1,0)\,c,(0,1)\,c,(-1,0)\,c,(0,-1)\,c,&i=1,\ldots,4~{},\\\
(1,1)\,c,(-1,1)\,c,(-1,-1)\,c,(1,-1)\,c~{},&i=5,\ldots,8{}\end{array}\right.$
(2)

where $c={\Delta x}{}/{\Delta t}$.

One of the approaches most commonly employed for the collision operator
$\Omega_{i}(\mathbf{f},\mathbf{f}^{\text{eq}})$ was proposed by Bhatnagar,
Gross and Krook (1954), and can be written as:

$\Omega_{i}(\mathbf{f},\mathbf{f}^{\text{eq}})=-\frac{1}{\tau}(f_{i}-f_{i}^{\text{eq}})~{},$
(3)

where the parameter $\tau$ is the relaxation time. The BGK collision operator
is also known as the single-relaxation time operator. A more stable version of
this collision operator allows for multiple-relaxation times (MRT) (Higuera
and Jiménez, 1989; d’Humières, 1992), reading:

$\Omega_{i}(\mathbf{f},\mathbf{f}^{\text{eq}})=-\left[\mathbf{M}^{-1}\mathbf{\Lambda}\mathbf{M}\right]_{ij}(f_{j}-f_{j}^{\text{eq}})~{},$
(4)

$\mathbf{\Lambda}$ being the relaxation matrix

$\mathbf{\Lambda}=\text{diag}\left(\tau_{\rho}^{-1},\tau_{e}^{-1},\tau_{\zeta}^{-1},\tau_{j}^{-1},\tau_{q}^{-1},\tau_{j}^{-1},\tau_{q}^{-1},\tau_{\eta}^{-1},\tau_{\eta}^{-1}\right)~{},$

where $\tau_{\rho}$ and $\tau_{j}$ being relaxation parameters related to
density and momentum conserved moments, $\tau_{\eta}$ defines dynamic
viscosity through $\eta=c_{s}^{2}(\tau_{\eta}-0.5){\Delta t}$, $\tau_{e}$ is
associated with the bulk viscosity, $\tau_{\zeta}$ is related to the energy
square quantity, and $\tau_{q}$ is related to the energy flux. The matrix
$\mathbf{M}$ is the matrix that converts $f$ into the momentum space, defined
as (Lallemand and Luo, 2000):

$\mathbf{M}=\left(\begin{array}[]{rrrrrrrrr}1&1&1&1&1&1&1&1&1\\\
-4&-1&-1&-1&-1&2&2&2&2\\\ 4&-2&-2&-2&-2&1&1&1&1\\\ 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)\,.$

The equilibrium distribution function values $\mathbf{f}^{\text{eq}}$ can be
computed through its relation to the macroscopic density ($\rho$) and velocity
($\mathbf{u}=u_{\alpha}$) fields as follows:

$f_{i}^{eq}=w_{i}\bigg{(}\rho+\frac{c_{i\alpha}}{c_{s}^{2}}\rho
u_{\alpha}+\frac{(c_{i\alpha}c_{i\beta}-c_{s}^{2}\delta_{\alpha\beta})}{2c_{s}^{4}}\rho
u_{\alpha}u_{\beta}\bigg{)},$ (5)

where the weights $w_{i}$ are related to each discrete velocity
$\mathbf{c}_{i}$, and $c_{s}$ is the lattice sound speed. For the D2Q9
velocity set, one has $w_{0}=4{}/9$, $w_{1}=w_{2}=w_{3}=w_{4}=1{}/6$ and
$w_{5}=w_{6}=w_{7}=w_{8}=1{}/36$. In that way, the macroscopic density and
velocity may be recovered through following relations:

$\rho=\sum_{i}f_{i},\quad\rho\mathbf{u}=\sum_{i}f_{i}\mathbf{c}_{i}+\frac{\mathbf{F}}{2}~{}.$
(6)

It is worth mentioning that each $f_{i}$ represents the number of particles
placed between $\mathbf{x}$ and $\mathbf{x}+d\mathbf{x}$ at time $t$, and
moving with velocity $\mathbf{c}_{i}$. In the LBM, the spacial ($\Delta x$)
and temporal $(\Delta t)$ discrete quantities are commonly fixed to the unity
(e.g., 1 m and 1 s when SI units are employed). In Eq. (1) that assumption was
considered, implying that $t^{n+1}=t^{n}+1$ and
$d\mathbf{x}=\mathbf{c}_{i}\Delta x=\mathbf{c}_{i}$.

Therefore, once one has obtained numerical results with a given number of
timesteps and lattices distribution, the refinement of these results, to
assess mesh sensibility, requires recomputing the non-dimensional parameters
that represent the particular physical setting. This issue is further
developed in Section 3.

###  2.1 The pseudopotential method

The last term from Eq. (1), $\mathbf{R}$, is the source term, which is
described by the forcing scheme employed. This source term is used for adding
effects of a general external force field, $\mathbf{F}$. In this work, the
enhanced Li and Luo (2013) forcing scheme, based on original scheme from Guo
and Zhao (2002), along with Li et al. (2013) surface tension correction term
are employed.

In this work, the external force field reads
$\mathbf{F}=\mathbf{F}^{m}+\mathbf{F}^{\text{ads}}+\mathbf{F}^{b}$, where
$\mathbf{F}^{m}$ is the intermolecular force, $\mathbf{F}^{\text{ads}}$ the
adherence force, and $\mathbf{F}^{b}$ the body force. The intermolecular force
field models the bulk two-phase physics through a pseudopotential function,
which depends on $\rho$. For single-component fluids, the interaction force
may be written (Chen and Doolen, 1998):

$\mathbf{F}_{i}^{m}(\mathbf{x},t)=-G\,\psi(\mathbf{x},t)\sum_{i}\text{w}_{i}\psi(\mathbf{x}+\mathbf{c}_{i},t)\,\mathbf{c}_{i},$
(7)

where $\psi$ is the pseudopotential field and $G$ is the interaction strength.
For the intermolecular force, the weights $\text{w}_{i}$ are equal to
$\text{w}_{i}=1{}/3$ for $|\mathbf{c}_{i}|^{2}=1$, and $\text{w}_{i}=1{}/12$
when $|\mathbf{c}_{i}|^{2}=2$ (Sbragaglia et al., 2007). A prescribed non-
ideal EOS can be imposed by using the following relation for the
pseudopotential:

$\psi(\mathbf{x},t)=\sqrt{\frac{2(p_{\text{\scriptsize
EOS}}(\mathbf{x},t)-\rho(\mathbf{x},t)/3)}{G}}~{},$ (8)

where $p_{\text{\scriptsize EOS}}$ is the thermodynamic pressure. When using
this definition, the parameter $G$ no longer represents intermolecular
strength, and it is used only to assure a positive value inside the square
root. The value $G=-1$ is assumed throughout this work.

In this work, the Peng-Robinson EOS is used, which can be written:

$p_{\text{\scriptsize EOS}}=\frac{\rho
RT}{1-b\rho}-\frac{a\rho^{2}\left[1+\left(0.37464+1.54226\,\omega-0.2699\,\omega^{2}\right)\left(1-\sqrt{T/T_{c}}\right)\right]^{2}}{1+2b\rho-b^{2}\rho^{2}}~{},$
(9)

where $T=T(\mathbf{x},t)$ is the temperature field, $T_{c}$ is the critical
temperature, $\omega$ is the acentric factor (fixed to $9.52381\times
10^{-2}$), $R$ is the universal gas constant (fixed to 1),
$a=0.45724\,R^{2}T_{c}^{2}/p_{c}$, $b=0.0778\,RT_{c}/p_{c}$, and $p_{c}$ is
the pressure at the critical point.

Gravitational force effects are incorporated by taking:

$\mathbf{F}^{b}=\mathbf{g}\,(\rho-\overline{\rho})~{},$ (10)

where $\mathbf{g}$ is the gravitational acceleration, and $\overline{\rho}$
the spatial mean density computed throughout entire computational domain at
every timestep.

The adherence force $\mathbf{F}^{\text{ads}}$ may be interpreted as a boundary
treatment of the intermolecular force $\mathbf{F}^{m}$. The contact angle
models considered are better discussed on Section 2.2. The source term
$\mathbf{R}$ proposed by Li et al. (2013) can be defined as:

$R_{ij}=\left[\mathbf{M}^{-1}\left(\mathbf{I}-\frac{\mathbf{\Lambda}}{2}\right)\right]_{ij}S_{j}+\left[\mathbf{M}^{-1}\right]_{ij}C_{j},$
(11)

where the forcing term $\mathbf{S}$ is computed by:

$\mathbf{S}=\begin{pmatrix}0\\\
6(u_{x}F_{x}+u_{y}F_{y})+\frac{12\sigma\,(\mathbf{F}^{m}\cdot\mathbf{F}^{m})^{2}}{\psi^{2}(\tau_{e}-0.5)}\\\
-6(u_{x}F_{x}+u_{y}F_{y})-\frac{12\sigma\,(\mathbf{F}^{m}\cdot\mathbf{F}^{m})^{2}}{\psi^{2}(\tau_{\zeta}-0.5)}\\\
F_{x}\\\ -F_{x}\\\ F_{y}\\\ -F_{y}\\\ 2(u_{x}F_{x}-u_{y}F_{y})\\\
u_{x}F_{y}+u_{y}F_{x}\end{pmatrix}~{}.$ (12)

In Eq. (12), the parameter $\sigma$ is used to adjust the thermodynamic
consistency, which is specially needed when employing the definition given on
Eq. (8) (Lycett-Brown and Luo, 2015). If $\sigma=0$, first term of Eq. (11)
simplifies to traditional Guo et al. (2002) forcing scheme.

As for the additional term $\mathbf{C}$, its definition is as follows:

$\mathbf{C}=\left[\begin{array}[]{c}0\\\
1.5\tau_{e}^{-1}\left(Q_{xx}+Q_{yy}\right)\\\
-1.5\tau_{\zeta}^{-1}\left(Q_{xx}+Q_{yy}\right)\\\ 0\\\ 0\\\ 0\\\ 0\\\
-\tau_{v}^{-1}\left(Q_{xx}-Q_{yy}\right)\\\
-\tau_{v}^{-1}Q_{xy}\end{array}\right]~{},$ (13)

with

$\qquad\mathbf{Q}=\kappa\,\frac{G}{2}\psi(\mathbf{x},t)\sum_{i}w\left(|\mathbf{c}_{i}|^{2}\right)\left[\psi(\mathbf{x}+\mathbf{c}_{i},t)-\psi(\mathbf{x},t)\right]\mathbf{c}_{i}\mathbf{c}_{i}~{}.$
(14)

$\mathbf{C}$ is a correction term used to tune the desired surface tension
through the adjustment of the parameter $\kappa$ present in $\mathbf{Q}$ (Li
and Luo, 2013).

###  2.2 Contact angle modeling and surface tension at the solid walls

The fluid-solid interaction is given by the Zou-He scheme (Zou and He, 1997)
and two contact angle models:

The $\psi$-based model:

In this case, the adherence forces are given by (Benzi et al., 2006; Huang et
al., 2009)

$\mathbf{F}^{\text{ads}}(\mathbf{x},t)=-G\,\psi(\rho(\mathbf{x},t),T(\mathbf{x},t))\,\psi(\rho_{w}(\mathbf{x}),T^{*}_{w}(\mathbf{x},t))\sum_{i}w_{i}\,s(\mathbf{x}+\mathbf{c}_{i})\,\mathbf{c}_{i}~{},$
(15)

where $s$ is an indicator function that is equal to the unity only when
$\mathbf{x}+\mathbf{c}_{\alpha}$ corresponds to a solid lattice, and
$\rho_{w}(\mathbf{x})$ is an auxiliary function used to model the local solid-
fluid interaction.

The modified-$\psi$-based model:

Under this scheme one imposes (Li et al., 2014):

$\mathbf{F}^{\text{ads}}(\mathbf{x},t)=-G_{w}(\mathbf{x})\,\psi^{2}(\rho(\mathbf{x},t),T(\mathbf{x},t))\sum_{i}w_{i}\,s(\mathbf{x}+\mathbf{c}_{i})\,\mathbf{c}_{i}~{}.$
(16)

When using the $\psi$-based scheme, the fluid-solid interaction is modeled
through the imposition of a pair of functions $(\rho_{w},T^{*}_{w})$ that may
depend on both space and time coordinates.

Having selected one of these models for the fluid-solid interaction links, at
the wall lattices the Eq. (7) is replaced by

$\mathbf{F}_{i}^{m}(\mathbf{x},t)=-G\,\psi(\mathbf{x},t)\left(\sum_{j\in
I^{\text{fluid}}}w_{j}\psi(\mathbf{x}+\mathbf{c}_{j},t)\,c_{ji}\right)+\mathbf{F}_{i}^{\text{ads}}(\mathbf{x},t)~{},$
(17)

where $I^{\text{fluid}}$ denotes the set of indices for which
$\mathbf{x}+\mathbf{c}_{j}$ is a fluid lattice.

####  2.2.1 Surface tension at the boundary walls

The correction term defined by Eqs. (13)-(14) corresponds to discrete
approximations of derivatives of the scalar field $\psi$. Thus, its evaluation
at wall lattices requires some attention and it is not clear how this
calculation have been performed in the literature. In this work, we have
calculated this correction term at the boundary walls following the contact
angle model computational procedure. More precisely, when the $\psi$-based
model is imposed, the term $\mathbf{Q}$ is computed as:

$\displaystyle\mathbf{Q}$
$\displaystyle=\kappa\,\frac{G}{2}\psi(\mathbf{x},t)\sum_{i\in
I^{\text{fluid}}}w\left(|\mathbf{c}_{i}|^{2}\right)\left[\psi(\mathbf{x}+\mathbf{c}_{i},t)-\psi(\mathbf{x},t)\right]\mathbf{c}_{i}\mathbf{c}_{i}$
$\displaystyle+\kappa\,\frac{G}{2}\psi(\mathbf{x},t)\sum_{i\in
I^{\text{solid}}}w\left(|\mathbf{c}_{i}|^{2}\right)\left[\psi(\rho_{w}(\mathbf{x},t),T^{*}_{w}(\mathbf{x},t))-\psi(\mathbf{x},t)\right]\mathbf{c}_{i}\mathbf{c}_{i}~{},$
(18)

analogously, when the modified-$\psi$-based model is used, $Q$ reads

$\displaystyle\mathbf{Q}$
$\displaystyle=\kappa\,\frac{G}{2}\psi(\mathbf{x},t)\sum_{i\in
I^{\text{fluid}}}w\left(|\mathbf{c}_{i}|^{2}\right)\left[\psi(\mathbf{x}+\mathbf{c}_{i},t)-\psi(\mathbf{x},t)\right]\mathbf{c}_{i}\mathbf{c}_{i}\,,$
(19)

which is equivalent to substitute $\psi(\mathbf{x}+\mathbf{c}_{i},t)$ by
$\psi(\mathbf{x},t)$ in Eq. (14) when $\mathbf{x}+\mathbf{c}_{i}$ corresponds
to a solid lattice.

###  2.3 The energy conservation equation

The energy conservation equation is solved for computing the scalar
temperature field. In the present work this equation is solved by the Finite
Difference (FD) method, following the hybrid approach (Zhang and Chen, 2003;
Li et al., 2015). The evolution of the temperature field, derived from the
energy conservation equation, is governed by (Li et al., 2020):

$\partial_{t}T+\mathbf{u}\cdot\nabla T=\frac{1}{\rho c_{\text{\scriptsize
V}}}\nabla\cdot(\lambda\nabla T)-\frac{1}{\rho c_{\text{\scriptsize
V}}}T\left(\frac{\partial p_{\text{\scriptsize EOS}}}{\partial
T}\right)_{\rho}\nabla\cdot\mathbf{u}~{},$ (20)

with $c_{\text{\scriptsize V}}$ being the constant volume specific heat, and
$\lambda=\lambda(\mathbf{x},t)$ is the thermal conductivity, given by

$\lambda(\mathbf{x},t)=\frac{\lambda_{\ell}(\rho(\mathbf{x},t)-\rho_{\text{\scriptsize
v}})+\lambda_{\text{\scriptsize
v}}(\rho_{\ell}-\rho(\mathbf{x},t))}{\rho_{\ell}-\rho_{\text{\scriptsize
v}}},$

where $\rho_{\ell}$ ($\lambda_{\ell}$) and $\rho_{\text{\scriptsize v}}$
($\lambda_{\text{\scriptsize v}}$) are the density (thermal conductivity) of
the liquid and vapor phase respectively. As usual, the viscous dissipation
term is assumed to be negligible. For all the simulations performed in this
work, we have taken the values of $c_{\text{\scriptsize V}}$ of both phases as
being equal and constant.

Equation (20) is discretized through a FD method along the lattices discrete
values and weights as described in Li et al. (2018b); Guzella et al. (2020)
regardless of the type of boundary condition for $T$ (Dirichlet or heat flux).
No special treatment is needed at the boundaries since the _wet-node_ approach
is used for the LBM (Zou and He, 1997). While the time derivatives are
computed according to a second-order Runge-Kutta method adapted from LeVeque
(2007).

##  3 Proposed mesh refinement procedure

In this section, we present the proposed mesh refinement procedure for
simulating liquid/vapor phase-change problems with heat transfer. This
procedure will be applied to the solution of the three typical problems
presented in Section 4. However, the proposed approach can be applied for
simulating other phase-change problems with LBM.

To refine the grid, it is necessary to choose the set of non-dimensional
parameters that define the physical problem. The definition of these non-
dimensional parameters is based on working at the viscous regime, considering
$\Delta\hat{t}\propto\Delta\hat{x}^{2}$. $\Delta\hat{t}$ and $\Delta\hat{x}$
are the physical time step and the physical distance between lattices,
respectively. Through this work we use the hat symbol “^" to refer to physical
quantities in SI units, otherwise, the quantities are assumed to be expressed
in lattice units.

###  3.1 Physical dependence between meshes

The proposed procedure is addressed to establish the necessary conditions for
performing a mesh refinement study to simulate liquid/vapor phase-change
problems with heat transfer. The main problem is related to pool boiling, so
the gravitational effects are not negligible. For this reason, the chosen
physical characteristic length, velocity, and time are those proposed by Son
et al. (1999), respectively:

$\hat{\ell}_{0}=\sqrt{\frac{\hat{\gamma}}{\hat{g}\left(\hat{\rho}_{\ell}-\hat{\rho}_{v}\right)}}\,,\qquad\hat{u}_{0}=\sqrt{\hat{g}\,\hat{\ell}_{0}}\,,\qquad\hat{t}_{0}=\hat{\ell}_{0}/\hat{u}_{0}~{},$
(21)

In Eq. (21) $\hat{\ell}_{0}$ is the capillary length, which represents the
ratio between the surface tension and the buoyancy force effects,
$\hat{\gamma}$ is the surface tension coefficient, $\hat{g}$ is the
gravitational acceleration, and $\hat{\rho}_{\ell}$ ($\hat{\rho}_{v}$) is the
liquid (vapor) phase density. In lattice units (l.u), the corresponding non-
dimensional parameters are given by:

$\ell_{0}=\hat{\ell}_{0}/\Delta\hat{x}\,,\qquad
u_{0}=\sqrt{{g}\,\ell_{0}}\,,\qquad t_{0}={\ell}_{0}/{u}_{0}~{},$ (22)

where the gravitational acceleration in lattice units ($g$) is given by
$g=\frac{\gamma}{\ell_{0}^{2}\left({\rho}_{\ell}-{\rho}_{v}\right)}$.

The fluid kinematic viscosity ($\hat{\eta}$) is assumed to be a function of
each phase viscosity, liquid $\hat{\eta}_{\ell}$ and vapor $\hat{\eta}_{v}$,
according to:

$\hat{\eta}(\mathbf{x},t)=\frac{\hat{\rho}(\mathbf{x},t)-\hat{\rho}_{v}}{\hat{\rho}_{\ell}-\hat{\rho}_{v}}\,\hat{\eta}_{\ell}+\frac{\hat{\rho}_{\ell}-\hat{\rho}(\mathbf{x},t)}{\hat{\rho}_{\ell}-\hat{\rho}_{v}}\,\hat{\eta}_{v}~{},$
(23)

where $\hat{\eta}_{\alpha}=\hat{U}_{0}\hat{\ell}_{0}/\text{Re}_{\alpha}$,
$\alpha\in\\{\ell,v\\}$. The symbol $\text{Re}_{\alpha}$ represents the
respective Reynolds number of the $\alpha$ phase. Knowing that the physical
time step $\Delta\hat{t}$ is related to $\Delta\hat{x}$ by
$\Delta\hat{t}=\Delta\hat{x}^{2}\eta/\hat{\eta}$ (Krüger et al., 2017), the
following relations may be computed:

$\Delta\hat{t}=\Delta\hat{x}^{2}\frac{\eta}{\hat{\eta}}=\Delta\hat{x}^{2}\frac{u_{0}\ell_{0}}{\hat{u}_{0}\hat{\ell}_{0}}~{},$
(24)

Therefore, to remain in the viscous regime when changing $\Delta\hat{x}$ the
quantity ${u}_{0}{\ell}_{0}$ must not depend on $\Delta\hat{x}$. On the other
hand, it is possible to compute the relation:

${u}_{0}{\ell}_{0}=\sqrt{\frac{\hat{\ell}_{0}}{\hat{\rho}_{\ell}-\hat{\rho}_{v}}\frac{\gamma}{\Delta\hat{x}}}~{},$

and thus one concludes that it is necessary for $\gamma$ to scale as
$\mathcal{O}\left(\Delta\hat{x}\right)$. As observed by Czelusniak et al.
(2020), the only way to obtain that behavior is by adjusting the parameter
$\kappa$ in Eq. (14) for each $\Delta\hat{x}$, as the surface tension
(computed from the planar interface problem) is proportional to $1-\kappa$.

Before discussing how this is achieved in this work, let us summarize the
physical quantities we want to keep constant among the different meshes. These
quantities are chosen in order to ensure the required consistency from
similarity laws, mainly considering the hydrodynamic and thermodynamic
similarities.

  * •

Hydrodynamic consistency: The Reynolds number $\text{Re}_{\alpha}$ of each
phase, as well as, the surface tension are keep constant. Particularly, the
last condition means that the Laplace equation

$\Delta P_{\Delta\hat{x}}=\frac{\gamma}{R_{\Delta\hat{x}}}~{},$

must remain valid for any $\Delta\hat{x}$, with $\Delta P_{\Delta\hat{x}}$
being the difference in the pressure field (given by an EOS) inside/outside
the bubble, and $R_{\Delta\hat{x}}$ the bubble radius. As we want
$\gamma=\mathcal{O}\left(\Delta\hat{x}\right)$ and $R_{\Delta\hat{x}}$ must
scale as $\mathcal{O}\left(\Delta\hat{x}^{-1}\right)$ this implies that one
must have $\Delta
P_{\Delta\hat{x}}=\mathcal{O}\left(\Delta\hat{x}^{2}\right)$.

  * •

Thermodynamic consistency: The density ratio is kept constant for a particular
imposed saturation temperature ($T_{s}$) and all meshes. The coexistence curve
is adjusted through the parameter $\sigma$ to maintain the same coexistence
densities for the $T_{s}$ imposed. This procedure is not exact, in the sense
that the coexistence curve $(T_{s},\rho)$ is not exactly the same when
changing $\Delta\hat{x}$ (Li et al., 2013). Due to this reason, the surface
tension adjustment through the parameter $\kappa$ needs to be performed for
each mesh as a final step of the computation of the parameters.

itemize

Hydrodynamic: by keeping constant the Reynolds number $\text{Re}_{\alpha}$ of
each phase, and conserving the physics related to surface tension.
Particularly, the last means that the Laplace equation

$\Delta P_{\Delta\hat{x}}=\frac{\gamma}{R_{\Delta\hat{x}}}~{},$

must remain valid for any $\Delta\hat{x}$, with $\Delta P_{\Delta\hat{x}}$
being the difference in the pressure field (given by an Equation of State)
inside/outside the bubble, and $R_{\Delta\hat{x}}$ the bubble’s radius. As we
want $\gamma=\mathcal{O}\left(\Delta\hat{x}\right)$ and $R_{\Delta\hat{x}}$
must scale as $\mathcal{O}\left(\Delta\hat{x}^{-1}\right)$ this implies that
one must have $\Delta
P_{\Delta\hat{x}}=\mathcal{O}\left(\Delta\hat{x}^{2}\right)$.

Thermodynamics: The saturation temperature ($T_{s}$) is kept constant among
meshes. The coexistence curve is adjusted through the parameter $\sigma$ in
order to maintain the same coexistence densities for the $T_{s}$ imposed. This
procedure is not exact, in the sense that the coexistence curve $(T_{s},\rho)$
is not exactly the same when changing $\Delta\hat{x}$ (Li et al., 2013).
Because of this, the surface tension adjustment through the parameter $\kappa$
needs to be performed for each mesh as a final step of the parameters
computations.

The conversion factor for some arbitrary field $\theta=\theta(\mathbf{x},t)$
is the constant $C_{\theta}$ such that the physical field
$\hat{\theta}(\mathbf{x},t)$ (e.g., expressed in SI units) can be obtained
from the lattice units term $\theta(\mathbf{x},t)$ as
$\hat{\theta}(\mathbf{x},t)=C_{\theta}\,\theta(\mathbf{x},t)$ for any pair
$(\mathbf{x},t)$. The conversion factors $C_{\rho}$, $C_{p}$ and $C_{T}$
correspond to the critical values $\rho_{c}$, $p_{c}$ and $T_{c}$ (critical
density, pressure and temperature, respectively), and they are determined by
the EOS in use, which in this work corresponds to the Peng-Robinson EOS, see
Eq. (9). It can be proved that for this EOS $\rho_{c}=\rho_{c}(b)$,
$p_{c}\propto a$ and $T_{c}\propto a$. Due to these relations, in order to
make the pressure field variation proportional to $\Delta x^{2}$, when
changing the mesh parameter $\Delta x$ we choose to fix the value of $b$ and
to vary $a$ as:

$a(\Delta x)=a_{\text{base}}\,\frac{\Delta x^{2}}{\Delta x_{0}^{2}}~{},$ (25)

where $a_{\text{base}}$ is a reference value such that $a(\Delta
x_{0})=a_{\text{base}}$ at the coarsest mesh ($\Delta x_{0}$) set for each
numerical test. In Table 1 we summarize the scaling of each field of interest
with $\Delta x$ (or, equivalently, with $\Delta\hat{x}$).

Symbol | Description | Conversion factor | Field Scaling (lu)
---|---|---|---
| Length | $C_{\ell}=\hat{\ell}_{0}/\ell_{0}$ | $\mathcal{O}\left(\Delta x^{-1}\right)$
| Time | $C_{t}=C_{\ell}\sqrt{g\ell_{0}}/\hat{u}_{0}$ | $\mathcal{O}\left(\Delta x^{-2}\right)$
$\rho$ | Density | $C_{\rho}=\hat{\rho}_{c}/\rho_{c}$ | $\mathcal{O}\left(1\right)$
$p$ | Pressure | $C_{p}=\hat{p}_{c}/p_{c}$ | $\mathcal{O}\left(\Delta x^{2}\right)$
$T$ | Temperature | $C_{T}=\hat{T}_{c}/T_{c}$ | $\mathcal{O}\left(\Delta x^{2}\right)$
$\lambda$ | Thermal conductivity | $\left(C_{\rho}C_{\ell}^{4}\right)/\left(C_{T}C_{t}^{3}\right)$ | $\mathcal{O}\left(1\right)$
$C_{v}$ | Specific heat at constant pressure | $C_{\ell}^{2}/\left(C_{t}^{2}\,C_{T}\right)$ | $\mathcal{O}\left(1\right)$
$\eta$ | Viscosity | $C_{\ell}^{2}/C_{t}$ | $\mathcal{O}\left(1\right)$
$\mathbf{g}$ | Gravitational acceleration | $C_{\ell}/C_{t}^{2}$ | $\mathcal{O}\left(\Delta x^{3}\right)$
Table 1: Scaling of the fields when changing the refinement parameter.

In summary, hydrodynamic and thermodynamic similarities imply in varying three
parameters in order to keep the physical consistency:
$a\propto\Delta\hat{x}^{2}$, $\gamma=\gamma(\Delta\hat{x})$ and
$\kappa=\kappa(\Delta\hat{x})$. Although the relation of $a$ parameter to grid
refinement is clear, $\sigma$ and $\kappa$ parameters dependence on $\Delta x$
is not trivial, and for that reason, they should be computed through numerical
procedures.

###  3.2 Mesh refinement example

For a better exposition of the employed methodology, we take as an example the
solution of a generic one-dimensional problem using a coarse base grid of
$L=128$ lattices, setting $a=6.12\times 10^{-2}$ and $T_{r}=0.86\,T_{c}$.
Also, assume one is interested in a fixed mesh ratio equal to $1.5$ for three
finer meshes. This results in meshes with $L=192$, $288$ and $432$. Let us
remark that these meshes were employed for the Stefan-problem study in Sec.
4.1.

The first step in the grid refinement procedure is to adjust the Peng-Robinson
parameter $a$ for each mesh. This parameter is computed through Eq. (25),
which results in values of $a$ equal to $6.12\times 10^{-2}$, $2.72\times
10^{-2}$, $1.21\times 10^{-2}$ and $5.37\times 10^{-3}$ for, respectively,
$L=128$, $192$, $288$ and $432$. This parameter influences the hydrodynamic
similarity through the Laplace equation, and also acts on the thermodynamic
consistency (Lycett-Brown and Luo, 2015; Czelusniak et al., 2020), specially
on the vapor branch.

The thermodynamic consistency may be ensured by solving the planar coexistence
curve equation, given by the following expression (Li and Luo, 2013; Li et
al., 2013):

$\int_{\rho_{g}}^{\rho_{l}}\left(p_{sat}-p_{\text{\scriptsize
EOS}}\right)\frac{\dot{\psi}}{\psi^{1-16G\sigma}}d\rho=0,$ (26)

where, for a given value of $a$, the parameter $\sigma$ (Li and Luo, 2013) is
adjusted for this equation to be satisfied. This is the second step of the
proposed refinement procedure.

In Fig. 1(a), coexistence curves for meshes from $L=128$ to $L=432$ and their
respective values of $a$ are exhibited for $\sigma$ equal to $0.103$, which is
a well-known value that ensures a better thermodynamic consistency (Li and
Luo, 2013; Li et al., 2013). In this case, it is possible to see how planar
coexistence vapor densities deviate from the Maxwell rule values (Clerk-
Maxwell, 1875), as Peng-Robinson $a$ parameter decreases. To counteract this
effect, $\sigma$ is adjusted for each considered grid. The correct value is
found by employing the secant root method to solve Eq. (26) with an initial
guess of $-16\,G\sigma_{0}=2$ (Lycett-Brown and Luo, 2015). Coexistence curves
obtained are shown in Fig. 1(b), where a better agreement between the
simulated coexistence curves and the Maxwell rule values is observed, which
thus warranties the thermodynamic consistency for each refinement.

(a)

(b)

Figure 1: Planar interface coexistence curves for distinct Peng-Robinson $a$
parameter values with fixed $\sigma=0.103$ parameter (a) and adjusting
$\sigma$ parameter for keeping thermodynamic consistency (b).

The third and last step of the grid refinement procedure consists in ensuring
a linear decrease in the surface tension $\gamma=\mathcal{O}(\Delta\hat{x})$.
For the Li et al. (2013) surface tension calculation method, defined by Eqs.
(13)-(14), it is possible to show that the surface tension for the planar
interface problem is given by:

$\gamma=\frac{Gc^{4}}{6}(1-\kappa)\int\psi\frac{d^{2}\psi}{dx^{2}}\,dx.$ (27)

From Eq. (27), it is possible to verify that the surface tension value for the
planar interface problem is dependent on $\psi$ profile across vapor-liquid
interface, and consequently, on the density profile as well. In Fig. 2, one
can verify how the parameters $\sigma$ and $a$ for $L=128$ to $L=432$
influence the density profile. As aforementioned, the parameter $\sigma$
ensures same vapor and liquid densities, while EOS $a$ parameter ensures
$\Delta P_{\Delta\hat{x}}=\mathcal{O}\left(\Delta\hat{x}^{2}\right)$. However,
the variation of these parameters result in a thicker interface in lattice
units, so that the surface tension $\gamma$ is influenced as well. So, the
procedure to ensure a linear surface tension decrease with $(\Delta\hat{x})$
consists in finding a $\kappa$ parameter so that:

$\gamma_{\text{1}}(\kappa)=\gamma_{\text{0}}\frac{\Delta x}{\Delta
x_{\text{0}}},$ (28)

where subscripts $0$ and $1$ refer to, respectively, a fine and a coarse mesh.
The target finer mesh value $\gamma_{1}$ is computed from Eq. (28), and then
the parameter $\kappa$ is determined from Eq. (27) by substituting
$\gamma=\gamma_{1}$ and solving the right hand side integral for the planar
interface problem (Li and Luo, 2013; Czelusniak et al., 2020). This procedure
results in $\kappa$ equal to 0, $-5.17\times 10^{-2}$, $-7.44\times 10^{-2}$
and $-8.45\times 10^{-2}$ for, $L=128$, $192$, $288$ and $432$, respectively.

Figure 2: Planar interface density profile for Peng-Robinson EOS with $a$
parameter equal to $6.12\times 10^{-2}$, $2.72\times 10^{-2}$, $1.21\times
10^{-2}$ and $5.37\times 10^{-3}$

.

##  4 Numerical results

Liquid/vapor phase-change problems involve many physical phenomena including:
hydrodynamic effects, heat transport (convection and diffusion) between phases
and with external medium, thermodynamic properties influence on phase-change,
surface tension action of liquid/vapor and liquid/vapor-solid interfaces,
solid-fluid interaction and contact angle influence, among others. In this
section the proposed mesh refinement procedure, presented in Sec. 3 is applied
for simulating three problems related with liquid/vapor phase-change phenomena
with heat transfer.

First, two simpler problems are analyzed in order to develop some intuition
about the application of the developed mesh refinement procedure. In Sec. 4.1
a two-dimensional droplet vaporization in superheated vapor is presented. The
droplet is initially put into contact with a superheated vapor resulting in
its vaporization due to the heat transfer in the liquid-vapor interface. Then
in Sec. 4.2 a one-dimensional liquid/vapor Stefan problem, mimicking a one
dimensional boiling like problem, is solved for a configuration similar to the
one-dimensional Stefan problem presented in Safari et al. (2013). Finally, in
Sec. 4.3 the pool boiling problem is assessed considering two models for the
contact angle. In the solution of each problem the procedure explained in Sec.
3.1 is employed for refining the mesh and the resulting physical parameters
(in l.u.) are reported for each analysis. Unless otherwise stated, the results
are re-scaled to the coarsest mesh in order to perform the comparisons.

In the Appendix A, a brief description of the computational tool used for
obtaining all the numerical results of this work is given. Also, the
characteristics of the supercomputers employed are provided together with the
execution times of some pool boiling simulations when varying the mesh
refinement.

###  4.1 Droplet vaporization

The scheme of this problem is shown in Fig. 3. Initially, the liquid phase
temperature is set to $T_{\text{sat}}=0.86\,T_{c}$, while the vapor
temperature is set to $T_{v}=T_{c}$. The hydrodynamics setting is periodic on
both axis while the temperature at the boundary is kept constant and equal to
$T_{v}$. Due to the warmer condition of the vapor (superheated at the same
pressure), the liquid evaporates with time and the droplet radius diminishes
up-to the collapse.

Figure 3: Scheme of the 2D droplet vaporization problem.

The droplet is initialized at the center of the domain. At $t=0$, a Heaviside-
like function is used for both the temperature and density fields, i.e., for
an arbitrary field $\theta$, its shape at $t=0$ is written (where the sub-
indices indicate the respective phase)

$\theta(\mathbf{x})=\frac{1}{2}\left\\{\theta_{v}+\theta_{\ell}+\left(\theta_{\ell}-\theta_{v}\right)\tanh\left(\frac{4.6}{5.0}(d(\mathbf{x},\mathbf{x}_{c})-R_{0})\right)\right\\}~{},$
(29)

where $d(\mathbf{x},\mathbf{x}_{c})$ is the Cartesian distance from
$\mathbf{x}$ to the center of the domain ($\mathbf{x}_{c}$) and $R_{0}\equiv
R(t=0)$ is the initial radius.

For these simulations we fixed the relaxation parameters to
$\tau_{\rho}^{-1}=\tau_{e}^{-1}=\tau_{\zeta}^{-1}=\tau_{j}^{-1}=\tau_{q}^{-1}=1.25$.

####  4.1.1 Grid refinement

The results obtained for four different meshes are reported next. In Fig. 4 is
shown how the radius of the droplet $R(t)$ changes with time, while the
respective parameters in l.u. are given in Table 2. The following parameters
were assumed constant for all meshes: $\eta_{\ell}=\eta_{v}=0.1$,
$\lambda_{\ell}=\lambda_{v}=2/3$, $C_{\text{V};\ell}=C_{\text{V};v}=5.0$ and
$\mathbf{g}=0$. To assure stability at the beginning of the simulation, the
energy equation is not solved from $t=0$ to $t=10.000$. Thus, the temperature
field is constant along that time interval. It can be observed that $R(t)$
first oscillates and the mechanical equilibrium is reached afterward around
$t\approx 3.000$.

Figure 4: Evolution of the droplet radius in time, $R(t)$, for the four different meshes. L | steps $\times 10^{6}$ | $a$ | $\sigma$ | $\kappa$ | $\Delta^{*}R(t)$
---|---|---|---|---|---
128 | 0.40 | $6.12\times 10^{-3}$ | 0.103 | 0 | -
192 | 0.90 | $2.72\times 10^{-3}$ | 0.115 | $-5.17\times 10^{-2}$ | $2.83\times 10^{-2}$
288 | 2.03 | $1.21\times 10^{-3}$ | 0.121 | $-7.44\times 10^{-2}$ | $6.64\times 10^{-3}$
432 | 4.56 | $5.37\times 10^{-4}$ | 0.123 | $-8.45\times 10^{-2}$ | $1.35\times 10^{-3}$
Table 2: Lattice units parameters for the refinement of the droplet
vaporization.

The curve $R(t)$ displayed in Fig. 4 show a fast convergence as a function of
mesh size. This behavior is confirmed by the data presented in Table 2 for a
relative quantity $\Delta^{*}R(t)$, defined as:

$\Delta^{*}R(t)=\frac{\sum_{t=0}^{T}|R(t)-R^{*}(t)|}{\sum_{t=0}^{T}|R^{*}(t)|}~{},$
(30)

In Eq. (31) relation $R(t)$ is evaluated at a certain mesh and $R^{*}(t)$ is
evaluated at its immediate refined mesh. The values of $\Delta^{*}R(t)$ are
shown in the last column of Table 2. The value of the second row ($2.83\times
10^{-2}$) corresponds to the relative difference between the curves for
$L=128$ and $L=192$, while the value at the third row corresponds to the
difference between $L=192$ and $L=288$, and that presented in the last row to
the difference between $L=288$ and $L=432$, respectively. It can be observed
that the results present an asymptotic regime of convergence (related to the
different meshes) because the relative difference decreases at a ratio
approximately constant. In fact:

$\frac{\Delta^{*}R(t)_{128,192}}{\Delta^{*}R(t)_{192,288}}=\frac{2.83\times
10^{-2}}{6.64\times 10^{-3}}\approx
4.27\,,\qquad\frac{\Delta^{*}R(t)_{192,288}}{\Delta^{*}R(t)_{288,432}}=\frac{6.64\times
10^{-3}}{1.35\times 10^{-3}}\approx 4.91~{},$

which gives convergence orders of

$\log(4.27)/\log(1.5)\approx 3.58,\qquad\log(4.91)/\log(1.5)\approx 3.92~{},$

respectively. Let us recall that this convergence order is related to both
spatial and time refinement as we are imposing the viscous regime, where
$\Delta t\propto\Delta x^{2}$.

or this problem, a fast convergence is observed for the curve $R(t)$. In fact,
defining $\Delta^{*}R(t)$ as the relative quantity

$\Delta^{*}R(t)=\frac{\sum_{t=0}^{T}|R(t)-R^{*}(t)|}{\sum_{t=0}^{T}|R^{*}(t)|}~{},$
(31)

where $R(t)$ is evaluated at a certain mesh and $R^{*}(t)$ is its immediate
refinement, $\Delta^{*}R(t)$ is shown in the last column of Table 2. The value
of the second row ($2.83\times 10^{-2}$) corresponds to the relative
difference between the curves for $L=128$ and $L=192$; the value at the third
row corresponds to the difference between $L=192$ and $L=288$, and so on. One
can observe that the results seem to be at an asymptotic regime (related to
the different meshes) as that relative difference decreases at a ratio
approximately constant, in fact:

$\frac{\Delta^{*}R(t)_{128,192}}{\Delta^{*}R(t)_{192,288}}=\frac{2.83\times
10^{-2}}{6.64\times 10^{-3}}\approx
4.27\,,\qquad\frac{\Delta^{*}R(t)_{192,288}}{\Delta^{*}R(t)_{288,432}}=\frac{6.64\times
10^{-3}}{1.35\times 10^{-3}}\approx 4.91~{},$

which gives convergence orders of

$\log(4.27)/\log(1.5)\approx 3.58,\qquad\log(4.91)/\log(1.5)\approx 3.92~{},$

respectively. Let us recall that this convergence order is related to both
spatial and time refinement as we are imposing the viscous regime, where
$\Delta t\propto\Delta x^{2}$.

This high convergence order was the highest one observed during this work.
Also, it was the only case where phase-change and, at the same time, evidence
of an asymptotic convergence regime were found. It is worth noticing that on
this configuration no effect related to wall nodes (where fluid-solid
interaction takes place) is involved.

###  4.2 One-dimensional boiling-like

The second problem analyzed consists of a one-dimensional configuration where
the fluid liquid phase is initially in equilibrium at the saturation
temperature $T_{\text{sat}}=0.86\,T_{c}$ in the whole computational domain.
For $t>0$, the temperature at the left wall ($T_{w}$) assumes a value higher
than $T_{\text{sat}}$ and a phase-change from liquid to gas is triggered,
forming a liquid/vapor interface as shown in Fig. 5. This phase-change process
produces the interface movement from left to right, reaching the right end of
the domain at a finite time, where null Neumann boundary conditions (Lou et
al., 2013) are imposed.

Figure 5: Scheme of the one-dimensional boiling-like problem.

In the idealized incompressible limit, the temperature on the liquid phase
region is homogeneous and constant in time, being equal to the initial
saturation temperature of the fluid $T_{\ell}\equiv T_{\text{sat}}$. On the
other hand, in the vapor region, the temperature decays linearly from $T_{w}$
to the temperature of the liquid. This is the so-called one-dimensional Stefan
problem, see for instance Safari et al. (2013). The solution to this problem
(citar livros de transferência de calor clásicos) is obtained by solving the
energy conservation equation for each phase and applying the boundary
conditions at domain boundaries and the interface. This last condition assumes
a continuity of both mass and energy fluxes at the interface. The analytical
solution results in the following expression for the interface position in
time:

$s(t)=2\beta\sqrt{\alpha_{v}\,t}~{},$ (32)

where
$\alpha_{\theta}\equiv\lambda_{\theta}/(c_{\text{v},\theta}\,\rho_{\theta})$,
(for $\theta=\ell,v$) denotes the thermal diffusivity of each phase and
$\beta$ is the solution of the following transcendental equation:

$\beta\exp(\beta^{2})\,\text{erf}(\beta)=\frac{c_{\text{v},v}(T_{w}-T_{\ell})}{h_{fg}\sqrt{\pi}}~{},$
(33)

being $c_{\text{v},v}$ the specific heat at constant volume of the vapor and
$h_{fg}$ the latent heat of vaporization. This last property is computed
following (Gong and Cheng, 2013) from the Peng-Robinson EOS.

####  4.2.1 Grid refinement and compressibility effects

Two set of parameters were tested: For the first set, the physical parameters
are fixed as in Section 4.1 excepting the thermal conductivities, that are set
such that the thermal diffusivities correspond to
$\alpha_{\ell}=\alpha_{v}=0.3$, and the wall temperature is set to
$T_{w}=0.9\,T_{c}$. For the second set of parameters, it was imposed that
$c_{\text{v},\ell}=c_{\text{v},v}=10$, $\lambda_{\ell}=15$ and $\lambda_{v}=1$
in order to keep the diffusivities relatively constant ($\alpha_{\ell}=0.23$
and $\alpha_{v}=0.26$) while diminishing the effect of the compression term in
Eq. (20). On the left side, using the $\psi$-based contact angle model (see
Sec. 2.2), the wall is modeled as being hydrophobic by setting $\rho_{w}$ to
the vapor density $\rho_{v}$. For these simulations we fixed the relaxation
parameters to $\tau_{\rho}^{-1}=1.0$, $\tau_{e}^{-1}=\tau_{\zeta}^{-1}=0.8$,
$\tau_{j}^{-1}=1.0$ and $\tau_{q}^{-1}=1.1$.

The resulting curves $s(t)$ considering the first set of parameters are shown
in Fig. 6 a) for $L$ varying from 160 to 2560. The analytical solution (Eq.
(32)) is also included. Let us remark that this analytical solution is used
here as a reference and, as the energy equation (Eq. (20)) incorporates
compressibility effects, the numerical solutions obtained by the LBM do not
need to converge to this reference. In fact, one can observe that as the
number of lattices augments there appear oscillations of $s(t)$ around the
reference solution. Moreover, for this first set of parameters, there exist
some time intervals where the interface speed $s^{\prime}(t)$ is negative.

Figure 6: Interface position for the one-dimensional boiling-like problem for
different meshes. a) First set of parameters, b) Second set of parameters.

The relative change between the curves is computed analogously to Eq. (31) and
is shown in Table 3 (row 1st set of parameters). It can be observed that the
convergence ratio for this set of parameters is slow. Moreover, the
convergence order varies between 0.89 for the coarsest meshes and 0.24 for the
finest ones, which indicates that these results do not correspond to an
asymptotic convergence regime (which requires the convergence order to be
invariant among meshes (De Bortoli et al., 2015)).

The oscillations of the liquid/vapor interface seem to be a consequence of the
fluid compressibility. Numerical evidence supporting this assertion was first
obtained using the second set of parameters, for which the effects of the
compression term in Eq. (20) are reduced. The relative change between the
curves when refining the mesh for this second set of parameters is shown in
Table 3 (row 2nd set of parameters) and the corresponding curves $s(t)$ are
shown in Fig. 6 b). This time the oscillations are less pronounced and a
faster convergence is observed, with a convergence order that goes from 3.28
for the coarsest meshes to 1.59 for the finest ones. These results are a sign
of these interface position oscillations being not numerical artifacts but a
consequence of the physics constituting the current LBM simulation model.

| 320 | 640 | 1280 | 2560
---|---|---|---|---
$\Delta^{*}_{1}s(t)$ (1st set) | $8.86\times 10^{-2}$ | $4.77\times 10^{-2}$ | $4.57\times 10^{-2}$ | $3.89\times 10^{-2}$
$\Delta^{*}_{2}s(t)$ (2nd set) | $1.09\times 10^{-1}$ | $3.31\times 10^{-2}$ | $1.47\times 10^{-2}$ | $9.28\times 10^{-3}$
Table 3: Relative change of the curves $s(t)$ between meshes (from $L=320$ to
$L=2560$). Figure 7: Interface position for the one-dimensional boiling-like
problem for different meshes. a) First set of parameters, b) Second set of
parameters.

Following a procedure similar to the one exposed by Welch and Wilson (2000),
one can find that continuity of the mass and energy fluxes, at the idealized
incompressible limit, imply the following equation for the interface velocity,
${ds}/{dt}$:

$h_{fg}\rho_{v}\frac{ds}{dt}=-\lambda_{v}\left.\frac{\partial T}{\partial
x}\right|_{x=s^{-}}+h_{fg}\rho_{v}\left.u\right|_{x=s^{-}}+\lambda_{\ell}\left.\frac{\partial
T}{\partial x}\right|_{x=s^{+}}~{},$ (34)

where the superscripts “$+,-$" indicate the right or left limits respectively
and $u$ denotes the fluid velocity (liquid or vapor phases). The last two
terms of this equation are null at the incompressible limit. In fact,
disregarding these terms and integrating the resulting equation one recovers
the formula (32) (Welch and Wilson, 2000). Nevertheless, Eq. (34) can be used
to obtain some insight in our LBM simulations where compressible effects are
not negligible.

The density, temperature, and velocity profiles of the fluid, for the first
set of parameters and a series of timesteps, are shown in Fig. 7 a) around the
instant of time where the first interface oscillation takes place. Contrary to
what happens in the incompressible limit, the interface temperature oscillates
in time around $T_{\text{sat}}$. Let us remark from the last term in Eq. (20)
that the vaporization process tends to cool the interface, while the
liquefaction process has the opposite effect. One can notice that this
temperature different from $T_{\text{sat}}$ at the interface produces non-null
gradients of both the fields $T$ and $\rho$ (and thus non-null velocities due
to pressure gradients) that make the LBM solution differ from the analytical
one. For instance, at $t=50$ the derivative $\frac{\partial T}{\partial x}$ at
the right side of the interface is positive while the velocity of the fluid at
the left side (vapor region) is negative. Considering Eq. (34) and these two
facts, one could expect the velocity of the interface to decrease (relative to
the incompressible regime) or even to make it reverse its movement, as
observed around that instant of time. This relation between the oscillation of
the interface movement and these non-null gradients, depending on whether
$\left.T\right|_{x=s(t)}$ is greater or lower than $T_{\text{sat}}$ was
observed along with the whole simulation.

The profiles for $\rho$, $T$, and $u$ for the same timesteps but for the
second set of parameters are shown in Fig. 7 b). As for the case discussed
above, it was observed that $\left.T\right|_{x=s(t)}$ is lower than
$T_{\text{sat}}$ for the first time instants too, but its value is slightly
higher. As a consequence, this time the gradients in $T$ and $\rho$ are lower,
and also the vapor’s velocity near the interface (in absolute value), which
has the overall effect of the LBM solution to be more similar to the analytic
solution by Eq. (32).

###  4.3 Pool boiling

In this section, the refinement procedure described in Sec. 3 is applied for
the pool boiling problem. The numerical tests presented were performed using
the $\psi$-based contact angle model. To justify this choice, let us first
show results of refinement tests considering the physical phenomena related to
the wall/fluid interaction alone and the two contact angle models presented in
Section 2. For the results shown hereafter the timescale will be fixed by
always giving the time values relative to the coarsest mesh considered, so the
timestep value of a certain fine mesh is re-scaled considering that $\Delta
x=\Delta t^{2}$. For these simulations we fixed the relaxation parameters to
$\tau_{\rho}^{-1}=1.0$, $\tau_{e}^{-1}=\tau_{\zeta}^{-1}=0.8$,
$\tau_{j}^{-1}=1.0$ and $\tau_{q}^{-1}=1.1$.

####  4.3.1 Contact angle model results

For these tests, a droplet is initialized near a wall where the fluid/solid
interaction takes place. The initial interface is showed in Fig. 8. The
simulations run a fixed number of timesteps large enough for the dynamics to
reach a static regime. The resulting density fields for a mesh of $160\times
160$ lattices are shown in Fig. 9 for both the $\psi$-based and the modified
$\psi$-based models and three different contact angles. This initial mesh was
refined four times, using meshes of $240\times 240$, $360\times 360$,
$540\times 540$, and $810\times 810$ lattices (corresponding to a refinement
ratio of 1.5).

The relative change of the density fields (when reaching the static regime)
are shown in Table 4 for the results obtained along with the $\psi$-based
model. These results indicate a convergent behavior of the $\psi$-based
contact angle model without the need to adjust its defining parameter
$\rho_{w}$, with a convergence rate of around 0.9. On the other hand, the
simulations performed for the modified $\psi$-based model show that the
parameter $G_{w}$ needs to be adjusted to keep the contact angle (and thus the
implied physical behavior) constant when applying the proposed refinement
procedure. This can be observed in Fig. 8 where the static liquid/vapor
interface is shown for three different meshes and $G_{w}=-0.8$. Not knowing an
analytic method to determine this parameter, a numerical procedure should be
needed for each mesh and contact angle. For this reason, and to simplify the
exposition, in the next section we present pool boiling results by using only
the $\psi$-based model.

Figure 8: Resulting static liquid/vapor interfaces for the modified $\psi$-based model for $G_{w}=-0.8$ and different meshes. The initial interface is also displayed by a continuous curve. Figure 9: Resulting droplet density fields for two contact angle models for different surfaces smoothness. | 240 | 360 | 540 | 810
---|---|---|---|---
$\Delta^{*}\rho$ ($\rho_{w}=4.0$) | $10.4\times 10^{-2}$ | $8.7\times 10^{-2}$ | $6.3\times 10^{-2}$ | $4.4\times 10^{-2}$
$\Delta^{*}\rho$ ($\rho_{w}=2.3$) | $12.9\times 10^{-2}$ | $10.1\times 10^{-2}$ | $7.1\times 10^{-2}$ | $4.9\times 10^{-2}$
$\Delta^{*}\rho$ ($\rho_{w}=1.7$) | $12.4\times 10^{-2}$ | $9.1\times 10^{-2}$ | $6.0\times 10^{-2}$ | $4.1\times 10^{-2}$
Table 4: Relative change of the density profiles for the $\psi$-based contact
angle model.

####  4.3.2 Pool boiling results

Figure 10: Scheme of the 2D pool boiling problem, boundary and initial
conditions for the phases.

The two-dimensional configuration of the simulated pool boiling problem is
schematized in Fig. 10. Periodic boundary conditions are set at both the left
and right boundaries ($\mathbf{x}=0$ and $\mathbf{x}=L$) and the wet-node
schemes of Zou and He (1997) are used to prescribe no-slip walls at the bottom
and top positions. A region of $L_{0}$ lattices is identified at the middle of
the bottom wall to impose a particular condition either for temperature or for
the contact angle.

The initial interface vapor/liquid is set at half the domain height. Thus, the
initial conditions at the interior of the domain are:

$\hat{T}(\mathbf{x},t=0)=\hat{T}_{\text{sat}}~{},\qquad\hat{\rho}(\mathbf{x},t=0)=\begin{cases}\hat{\rho}_{\ell}&\text{if}~{}y<L/2\\\
\hat{\rho}_{v}&\text{otherwise}\end{cases}~{},$ (35)

along with a null velocity field.

At the upper wall, the temperature is fixed to the saturation temperature,
which is set to $T_{\text{sat}}=0.86\,T_{c}$. The same temperature is set at
the lower wall, excepting for the heater region (placed at half of the domain
and with a length equal to $L_{0}$) where the temperature, denoted by $T_{h}$,
assumes two possible values $T_{h}\in\\{1.0\,T_{c},1.2\,T_{c}\\}$. The
subscript $h$ is used to denote the parameter value at the heater.

Regarding the contact angle along the surfaces, the upper wall was set as
being hydrophobic by imposing $\rho_{w}=0.38$ there. At the lower surface, for
the wall lattices that are outside the heater, the surface was modeled as
being hydrophilic by setting $\rho_{w}=6.50$. At the heater lattices, two
configurations were tested: one hydrophobic, where
$\rho_{w}\equiv\rho_{wh}=0.38$, and one _intermediate_ , where
$\rho_{w}\equiv\rho_{wh}=3.44$. This gives a total of four possible
configurations.

The baseline set of parameters for the next simulations corresponds to
$160\times 160$ lattices, for which $a=6.12\times 10^{-3}$, $\sigma=0.103$,
$\kappa=0$, $c_{\scriptsize{\text{V}}}=5$, $g=-2.5\times 10^{-5}$,
$\eta_{\ell}=\eta_{v}=0.1$, $\lambda_{\ell}=1.95$ and $\lambda_{v}=0.114$.
Remembering that some of these parameters change when refining the mesh
according to Section 3.1.

#### Qualitative comparisons

Let us first show what happens if one just changes the number of lattices
without adjusting these three parameters ($a$, $\sigma$, and $\kappa$). An
example of this is shown in Fig. 11 for the case $\rho_{wh}=0.38$ and
$T_{h}=1.2\,T_{c}$. In that figure, the spatial and temporal meshes are
refined from left to right. In the first row of results the physical time is
fixed (considering $\Delta t\propto\Delta x^{2}$) at a time instant where the
first bubble is developing; and at the second row of results, the physical
time is set for the instant just before the departure of the bubble in each
simulation. The time values that are shown in that last row correspond to the
time values respective to the coarsest mesh ($160\times 160$) considering
$\Delta t\propto\Delta x^{2}$. From these results, it is clear that the
physical behavior is not consistent when refining the mesh, as the bubble
shape, size, and time departure have no relation among meshes.

Figure 11: Density profiles obtained by refining a simulation of the pool
boiling problem without adjusting the physical parameters among meshes.

Repeating these simulations but this time adjusting the physical parameters as
described in Section 3.1 the results show a coherent relation among meshes.
This is shown qualitatively in Fig. 12 for $\rho_{wh}=0.38$ and
$T_{h}=1.2\,T_{c}$, where the first row of snapshots corresponds to an instant
of time (all the timestep values are re-scaled to the coarsest mesh) in which
the first bubble is developing and its boundary is far from the initial
interface, and the second row corresponds to an instant of time just before
the bubble’s departure takes place. Let us remark two facts that can be
observed from these results: 1) for the density profiles obtained for the
coarsest mesh at both $t=12800$ and $t=23360$, there exists a timestep at each
refined simulation that is similar to the coarsest mesh density profile, and
2) there exists a _delay_) between the simulations that diminishes as $\Delta
t$ goes to 0. In fact, denoting by $\delta$ the value of this delay at an
arbitrary timestep, for the event of the first row (the growing bubbles) one
has $\delta=12800-9280=3520$ between the two first meshes,
$\delta=9280-7680=1600$ between the configurations $360\times 360$ and
$240\times 240$, and $\delta=7680-7360=320$ between the last pair of finer
meshes. A similar behavior can be computed for the time instant where the
first bubble is about to depart, this time $\delta$ assumes the values $4800$,
1920, and 960 from the coarser to the finer pair of consecutive grids. These
qualitative observations are numerically assessed in the next section.

Figure 12: Density profiles obtained by refining a simulation of the pool
boiling problem adjusting the physical parameters among meshes according to
Section 3.1.

#### Time and space differences between meshes

To analyze quantitatively the convergence of the fields of interest let us
introduce the space relative difference for an arbitrary field
$\theta=\theta(\mathbf{x},t)$ between two meshes:

$\text{d}^{N}\theta(t)=\frac{\sum_{i\in\mathcal{S}}\left|\mathcal{I}(\theta^{N})(x_{i},t)-\theta^{N^{*}}(x_{i},t)\right|}{\sum_{i\in\mathcal{S}}\left|\theta^{N^{*}}(x_{i},t)\right|}~{},$
(36)

for arbitrary $N$ and $N^{*}$ denoting the immediate finer mesh (e.g., for
$N=160$ one has $N^{*}=240$), where $\mathcal{S}$ refers to the spatial points
indices of the fine mesh with $N^{*}\times N^{*}$ lattices, $\theta^{N^{*}}$
denotes the discrete function defined on the fine grid, and
$\mathcal{I}(\theta^{N})$ denotes the discrete function obtained interpolating
linearly the discrete function $\theta^{N}$ into the fine mesh. The global
space-time relative difference is computed as:

$\text{D}^{N}\theta=\frac{\sum_{i,k\in\mathcal{S},\mathcal{T}}\left|\mathcal{I}(\theta^{N})(x_{i},t_{k})-\theta^{N^{*}}(x_{i},t_{k})\right|}{\sum_{i,k\in\mathcal{S},\mathcal{T}}\left|\theta^{N^{*}}(x_{i},t_{k})\right|}~{},$
(37)

where $\mathcal{T}$ refers to the temporal points indices.

The simulations run a number of timesteps such that the first bubble
departures and the second bubble have been formed for three of the four
configurations tested here (more precisely, running along 40000 timesteps for
the coarsest mesh of $160\times 160$ lattices). The case $\rho_{wh}=3.44$ and
$T_{h}=1.0\,T_{c}$ corresponds to a configuration where nucleation is not
observed. The relative spatial differences for density and temperature,
computed according to Eq. (36), are shown in Fig. 13 and 14 respectively. As
expected, the spatial relative difference grows in time, and the refinement of
the meshes (augmenting the number of lattices and recomputing the physical
parameters as discussed in Section 3.1) has the effect of decreasing the
values of these differences in time. This happens more clearly in the case
without nucleation.

The overall order of convergence is obtained by computing the global error
given by Eq. (37). The resulting global errors for each configuration and the
respective orders of convergence are shown in Fig. 15 left column. It is worth
noticing that the convergence order depends on the parameters of the
simulation and the considered field. For instance, fixing $\rho_{wh}=3.44$,
one observes a convergence ratio of order 1 for the temperature field,
independently of the heater temperature, while for the same value of
$\rho_{wh}$ the convergence ratio for the density field is 1.2 for
$T_{h}=1.2\,T_{c}$ and 1.8 for the lower temperature $T_{h}=1.0\,T_{c}$. The
case $\rho_{wh}=0.38$ and $T_{h}=1.0\,T_{c}$ is more complex. As it can be
observed, the relative differences between meshes for the density field
presented almost no change for that case, diminishing with a ratio of $\approx
0.06$. However, next, we show that this low convergence order is due to the
delay of the simulations between meshes detected in Section 4.3. Below, we put
this into numerical terms.

#### Time-shifting the simulations

Instead of evaluating the relative difference employing Eq. (37), given two
meshes and two shift parameters $\delta$ and $\delta^{*}$, the shifted
relative difference is computed as:

$\text{D}^{N;\delta,\delta^{*}}\theta=\frac{\sum_{i,k\in\mathcal{S},\mathcal{T}}\left|\mathcal{I}(\theta^{N})\left(x_{i},t_{k}+\frac{\delta}{\tau}t_{k}\right)-\theta^{N^{*}}\left(x_{i},t_{k}+\frac{\delta^{*}}{\tau}t_{k}\right)\right|}{\sum_{i,k\in\mathcal{S},\mathcal{T}}\left|\theta^{N^{*}}\left(x_{i},t_{k}+\frac{\delta^{*}}{\tau}t_{k}\right)\right|}~{},$
(38)

where $\delta$ ($\delta^{*}$) is the shift necessary to _synchronize_ , at a
certain timestep denoted by $\tau$, the mesh $N$ ($N^{*}$) with the finest
mesh, which corresponds to $540\times 540$ lattices. For a given mesh with
$N\times N$ lattices, these $\delta$’s are computed by solving a one-
dimensional minimization problem:

$\delta=\underset{-m\leq\delta\leq
m}{\operatorname*{arg\,min}}~{}\sum_{i\in\mathcal{S}}\left|\mathcal{I}(\theta^{N})(x_{i},\tau+\delta)-\theta^{\text{finest}}(x_{i},\tau)\right|~{},$
(39)

where $m$ defines the searching interval size. For the four configurations
tested the reference time was set to $\tau=5120$ and the size of the search
interval size was set to $m=4800$, obtaining always $|\delta|<m$, which
warranties this $\delta$ is defined in terms of a local minimum. The values of
$\delta$ obtained for each case and mesh (relative to the finest mesh) are
given in Table 5. Notice that, as expected, the values of $\delta$ diminish as
the mesh is refined, which is an indicator of local convergence in time. The
case $\rho_{wh}=3.44$ and $T_{h}=1.0\,T_{c}$ presented no nucleation and for
$t=\tau$ the system already reached a stationary regime, which had as a
consequence that $\delta$ resulted to be null for the coarser meshes.

| $\rho_{wh}=0.38$ | $\rho_{wh}=3.44$
---|---|---
| $T_{h}/T_{c}=1.0$ | $T_{h}/T_{c}=1.2$ | $T_{h}/T_{c}=1.0$ | $T_{h}/T_{c}=1.2$
$N=160$ | 6720 | 3840 | 0 | 4480
$N=240$ | 2880 | 1600 | 0 | 1600
$N=360$ | 960 | 640 | 0 | 320
Table 5: Values of the timeshifts $\delta$ computed by Eq. (39) to synchronize
the simulations of each coarse mesh with the finest one at the time
$\tau=5120$. Figure 13: Relative differences in space per timestep, computed
according to Eq. (36) for the density field along three pairs of refinements
$(N,N^{*})\in\\{(160,240),(240,360),(360,540)\\}$. $\rho_{wh}$ and $T_{h}$
denote the $\psi$-based contact angle parameter and the temperature set at the
heater lattices respectively. Figure 14: Relative differences in space per
timestep, computed according to Eq. (36) for the temperature field along three
pairs of refinements $(N,N^{*})\in\\{(160,240),(240,360),(360,540)\\}$.
$\rho_{wh}$ and $T_{h}$ denote the $\psi$-based contact angle parameter and
the temperature set at the heater lattices respectively.

The resulting global relative differences with these timeshifts are shown at
the right side of Fig. 15. First, it is worth noticing that the trend of the
global differences among meshes is to diminish when the shifting is applied,
where the only exception is the case where no nucleation takes place (a case
for which $\delta$ is 0). Regarding the temperature field, the convergence
order was around the unit for all the cases tested. On the other hand, all the
cases where nucleation takes place exhibit a convergence order around the unit
for the density field, while for the case where nucleation is not observed the
convergence order is $\approx 2$.

Figure 15: Global relative differences for the pool boiling according to Eq.
(37) (left side) and to Eqs. (38)-(39) (right side).

These results are evidence of time and space convergence for the pool boiling
problem. As the relative differences employed here summarize the variations
along with both time and space in norm $L^{1}$, these differences going to
zero also represent evidence of convergence for a series of quantities
depending on $\rho$, like bubbles’ geometry at a certain time, time departure
and contact angle.

A quantity of engineering interest that depends both on $\rho$ and $T$ is the
average heat flux entering the computational domain, defined as

$Q(t)=\frac{1}{L}\int_{0}^{L}\left(\lambda(\rho(x,y,t))\frac{\partial
T}{\partial y}(x,y,t)\right)_{y=0}\,dx~{}.$ (40)

This quantity is shown in Fig. 16 for $T_{h}=1.2$ and two values of the
contact angle parameter at the heater lattices: a) $\rho_{wh}=0.38$, and b)
$\rho_{wh}=3.44$. For $\rho_{wh}=3.44$, taking the relative differences of the
curves in time, resulted in a convergence order of around 1.0. While for
$\rho_{wh}=0.38$ the convergence order was lower, around 0.6. This behavior
indicates that this ratio is very sensitive to the contact angle
configurations. Also, from these curves, one can observe the delay taking
place between the simulations and how they tend to synchronize as the number
of computational lattices is augmented. It should be noted that for these
cases we have changed the mesh refinement ratio to 2.0. For the refinement
ratio of 1.5 the same behavior was observed for both cases, and we choose to
show the more refined results.

With the mesh refinement we can also observe that the number of nucleated
bubbles in time increases for finer meshes, resulting in three bubbles for the
finest mesh when the heater contact angle parameter is $\rho_{wh}=0.38$ (see
Fig. 16 a). In this case the number of nucleated bubbles raised from two to
three, while for $\rho_{wh}=0.38$ the number of nucleated bubbles increase
from one to two (see Fig. 16 b). This fact, which is related to the previously
commented time delay for different mesh sizes, shows the importance of
reaching numerically converged solutions in order to correctly describe the
dynamics of the boiling process. Moreover, one observes that the increase of
the bubbles departure frequency is associated with a higher overall heat flux,
which is a well known physical behavior.

Figure 16: Heat fluxes calculated for four different meshes considering
$T_{h}=1.2$. a) $\rho_{wh}=0.38$ and b) $\rho_{wh}=3.44$.

##  5 Conclusions

The following conclusions can be drawn from the present work:

  * •

An original mesh refinement procedure has been proposed to obtain mesh
independent simulation results with pseudopotential LBM of phase-change heat
transfer problems, specially pool boiling problems, in the viscous regime. The
procedure was based on conserving the hydrodynamic and thermodynamic
similarities between different meshes by suitably adjusting the parameters
$a$, $\sigma$ and $\kappa$.

  * •

It has been shown that the attainment of numerically converged simulation
results of pool boiling phenomena depends on both the choice of the contact
angle model and its set-up parameters ($G_{w}$ and $\rho_{w}$ in this work).
The $\psi$-based model was the method that produced consistent results without
requiring the adjustment of its defining parameter.

  * •

A delay of pool boiling converged simulations results for different meshes has
been observed. This could affect quantities of engineering interest like
bubbles’ frequency departure, spatial average heat flux and boiling curves.

  * •

The pool boiling and the Stefan problems results suggest that the convergence
order is affected by the presence of nucleation and the compressibility of the
fluid. This shows that the incorporation of more complex physics into the
simulations affects the convergence ratio and thus would require more
computational effort to well describe the physical phenomena occurring in the
simulated problems.

  * •

The proposed procedure and consequently the simulation of numerically
converged results required an extensive use of computational resources,
indicating the necessity of using LBM parallel numerical codes for this
purpose.

  * •

The paper obtained results suggest the importance of presenting mesh refined
simulation results in order to correctly assess pool boiling problems with the
LBM.

## Acknowledgments

The authors acknowledge the National Laboratory for Scientific Computing
(LNCC/MCTI, Brazil) for providing HPC resources of the SDumont supercomputer,
which have contributed to the research results reported within this paper.

Part of this research was carried out using the computational resources of the
Center for Mathematical Sciences Applied to Industry (CeMEAI) funded by FAPESP
(grant 2013/07375-0).

The authors acknowledge the support received from the São Paulo Foundation for
Research Support (FAPESP), grants 2016/09509-1, 2020/12919-2, and the National
Council for Scientific and Technological Development (CNPq), grants
431782/2018-0, 140634/2019-3, 305941/2020-8.

## Appendix A

In this appendix some details on the computational machines and the execution
times are provided.

An in-house code was developed in the programming language C along with the
parallelization paradigm of the Message Passing Interface (MPI).

The simulations for the heat flux curves provided in Section 4.3.2 were
executed in the _Euler Supercomputer_ 111https://euler.cemeai.icmc.usp.br,
that is located at the Institute of Mathematical and Computer Sciences (ICMC)
of the University of São Paulo at São Carlos. All the other simulations were
executed in the cluster Santos Dumont, located at the National Laboratory for
Scientific Computing (LNCC) 222https://sdumont.lncc.br. The hardware
specifications of both clusters are listed in Table A1.

| _Cluster Sdumont_ | _Cluster Euler_
---|---|---
Processors | 2x Intel Xeon Cascade Lake Gold 6252 | 2x Intel Xeon `E5-2680v2`
Cores per processor | 24 (48 per node) | 10 (20 per nodes)
RAM per node | 384 Gb | 128 Gb
Communications | InfiniBand | InfiniBand
Table A1: Specifications of the computational nodes were the simulations were executed. Mesh | Execution time | # Cores used | Mlups
---|---|---|---
$160\times 160$ | 0 m 37 s | 40 | 30
$320\times 320$ | 4 m 21 s | 80 | 69
$640\times 640$ | 45 m 35 s | 80 | 105
$1280\times 1280$ | 391 m 36 s | 160 | 196
Table A2: Some efficiency data of the in-house code during the heat flux
simulations.

For completeness, the computational times needed to execute each of the 8
simulations summarized in Fig. 16 are given in Table A2. These computational
runs were performed in the Supercomputer Euler where the Intel Suite Compilers
was employed. Observe that, as it is typical in High Performace Computations,
the more refined the mesh is the more efficient the code execution. In fact,
for the smaller mesh our delivers around 0.75 Mlups per core, while for the
finer mesh the efficiency is of about 1.23 Mlups per core.

## References

  * Benzi et al. (2006) Benzi, R., Biferale, L., Sbragaglia, M., Succi, S., Toschi, F., 2006.  Mesoscopic modeling of a two-phase flow in the presence of boundaries: The contact angle.  Physical Review E 74, 021509\.
  * Bhatnagar et al. (1954) Bhatnagar, P.L., Gross, E.P., Krook, M., 1954.  A model for collision processes in gases. i. small amplitude processes in charged and neutral one-component systems.  Physical Review 94, 511–525.
  * Carey (2008) Carey, V.P., 2008.  Liquid-vapor phase-change phenomena.  Taylor & Francis.
  * Chen et al. (2004) Chen, H., Orszag, S.A., Staroselsky, I., Succi, S., 2004\.  Expanded analogy between Boltzmann kinetic theory of fluids and turbulence.  Journal of Fluid Mechanics 519, 301–314.
  * Chen and Doolen (1998) Chen, S., Doolen, G.D., 1998\.  Lattice Boltzmann method for fluid flows.  Annual review of fluid mechanics 30, 329–364.
  * Clerk-Maxwell (1875) Clerk-Maxwell, J., 1875.  On the dynamical evidence of the molecular constitution of bodies.  Nature 11, 357–359.
  * Czelusniak et al. (2020) Czelusniak, L.E., Mapelli, V.P., Guzella, M.S., Cabezas-Gómez, L., Wagner, A.J., 2020.  Force approach for the pseudopotential lattice Boltzmann method.  Physical Review E 102, 033307\.
  * De Bortoli et al. (2015) De Bortoli, A.L., Andreis, G.S., Pereira, F.N., 2015.  Modeling and Simulation of Reactive Flows.  Elsevier.
  * d’Humières (1992) d’Humières, D., 1992.  Generalized lattice-Boltzmann equations.  Rarefied gas dynamics .
  * Dong et al. (2018) Dong, L., Gong, S., Cheng, P., 2018.  Direct numerical simulations of film boiling heat transfer by a phase-change lattice Boltzmann method.  International Communications in Heat and Mass Transfer 91, 109–116.  URL: https://doi.org/10.1016/j.icheatmasstransfer.2017.12.013, doi:10.1016/j.icheatmasstransfer.2017.12.013.
  * Fariñas Alvariño et al. (2019) Fariñas Alvariño, P., Sánchez Simón, M.L., dos Santos Guzella, M., Amado Paz, J.M., Sáiz Jabardo, J.M., Cabezas-Gómez, L., 2019.  Experimental investigation of the CHF of HFE-7100 under pool boiling conditions on differently roughened surfaces.  International Journal of Heat and Mass Transfer 139, 269–279.
  * Gong and Cheng (2012) Gong, S., Cheng, P., 2012.  Numerical investigation of droplet motion and coalescence by an improved lattice Boltzmann model for phase transitions and multiphase flows.  Computers & Fluids 53, 93–104.
  * Gong and Cheng (2013) Gong, S., Cheng, P., 2013.  Lattice Boltzmann simulation of periodic bubble nucleation, growth and departure from a heated surface in pool boiling.  International Journal of Heat and Mass Transfer 64, 122–132.
  * Gong and Cheng (2015a) Gong, S., Cheng, P., 2015a.  Lattice Boltzmann simulations for surface wettability effects in saturated pool boiling heat transfer.  International Journal of Heat and Mass Transfer 85, 635–646.
  * Gong and Cheng (2015b) Gong, S., Cheng, P., 2015b.  Numerical simulation of pool boiling heat transfer on smooth surfaces with mixed wettability by lattice Boltzmann method.  International Journal of Heat and Mass Transfer 80, 206–216.
  * Gong and Cheng (2017) Gong, S., Cheng, P., 2017.  Direct numerical simulations of pool boiling curves including heater’s thermal responses and the effect of vapor phase’s thermal conductivity.  International Communications in Heat and Mass Transfer 87, 61–71.
  * Gong et al. (2018) Gong, W., Yan, Y.Y., Chen, S., Wright, E., 2018.  A modified phase change pseudopotential lattice Boltzmann model.  International Journal of Heat and Mass Transfer 125, 323–329.
  * Guan et al. (2011) Guan, C.K., Klausner, J.F., Mei, R., 2011.  A new mechanistic model for pool boiling chf on horizontal surfaces.  International Journal of Heat and Mass Transfer 54, 3960–3969.
  * Gunstensen et al. (1991) Gunstensen, A.K., Rothman, D.H., Zaleski, S., Zanetti, G., 1991\.  Lattice Boltzmann model of immiscible fluids.  Physical Review A 43, 4320–4327.
  * Guo and Zhao (2002) Guo, Z., Zhao, T.S., 2002.  Lattice Boltzmann model for incompressible flows through porous media.  Physical Review E 66, 036304\.
  * Guo et al. (2002) Guo, Z., Zheng, C., Shi, B., 2002.  Discrete lattice effects on the forcing term in the lattice Boltzmann method.  Physical Review E 65, 046308\.
  * Guzella et al. (2020) Guzella, M.d.S., Czelusniak, L.E., Mapelli, V.P., Alvariño, P.F., Ribatski, G., Cabezas-Gómez, L., 2020\.  Simulation of boiling heat transfer at different reduced temperatures with an improved pseudopotential lattice boltzmann method.  Symmetry 12, 1358\.
  * Haramura and Katto (1983) Haramura, Y., Katto, Y., 1983\.  A new hydrodynamic model of critical heat flux, applicable widely to both pool and forced convection boiling on submerged bodies in saturated liquids.  International Journal of Heat and Mass Transfer 26, 389–399.
  * Higuera and Jiménez (1989) Higuera, F.J., Jiménez, J., 1989\.  Boltzmann approach to lattice gas simulations.  Europhysics Letters 9, 663\.
  * Huang et al. (2009) Huang, H., Li, Z., Liu, S., Lu, X.y., 2009.  Shan-and-Chen-type multiphase lattice Boltzmann study of viscous coupling effects for two-phase flow in porous media.  International Journal for Numerical Methods in Fluids 61, 341–354.
  * Kandlikar (2001) Kandlikar, S.G., 2001.  A Theoretical Model to Predict Pool Boiling CHF Incorporating Effects of Contact Angle and Orientation .  Journal of Heat Transfer 123, 1071–1079.
  * Kharmiani et al. (2019) Kharmiani, S.F., Niazmand, H., Passandideh-Fard, M., 2019.  An alternative high-density ratio pseudo-potential lattice Boltzmann model with surface tension adjustment capability.  Journal of Statistical Physics 175, 47–70.
  * Kim et al. (2016) Kim, B.J., Lee, J.H., Kim, K.D., 2016.  Improvements of critical heat flux models for pool boiling on horizontal surfaces using interfacial instabilities of viscous potential flows.  International Journal of Heat and Mass Transfer 93, 200–206.
  * Krüger et al. (2017) Krüger, T., Kusumaatmaja, H., Kuzmin, A., Shardt, O., Silva, G., Viggen, E., 2017\.  The Lattice Boltzmann Method.  Springer International Publishing.
  * Kupershtokh et al. (2009) Kupershtokh, A.L., Medvedev, D.A., Karpov, D.I., 2009.  On equations of state in a lattice Boltzmann method.  Computers & Mathematics with Applications 58, 965–974.
  * Lallemand and Luo (2000) Lallemand, P., Luo, L.S., 2000\.  Theory of the lattice Boltzmann method: Dispersion, dissipation, isotropy, galilean invariance, and stability.  Physical Review E 61, 6546\.
  * LeVeque (2007) LeVeque, R., 2007.  Finite Difference Methods for Ordinary and Partial Differential Equations.  Society for Industrial and Applied Mathematics.
  * Li et al. (2018a) Li, Q., Huang, J.Y., Kang, Q.J., 2018a.  On the temperature equation in a phase change pseudopotential lattice Boltzmann model.  International Journal of Heat and Mass Transfer 127, 1112–1113.
  * Li et al. (2015) Li, Q., Kang, Q.J., Francois, M.M., He, Y.L., Luo, K.H., 2015.  Lattice Boltzmann modeling of boiling heat transfer: The boiling curve and the effects of wettability.  International Journal of Heat and Mass Transfer 85, 787–796.
  * Li and Luo (2013) Li, Q., Luo, K., 2013.  Achieving tunable surface tension in the pseudopotential lattice boltzmann modeling of multiphase flows.  Physical Review E 88.
  * Li et al. (2014) Li, Q., Luo, K., Kang, Q., Chen, Q., 2014.  Contact angles in the pseudopotential lattice boltzmann modeling of wetting.  Physical Review E 90.
  * Li et al. (2013) Li, Q., Luo, K.H., Li, X.J., 2013.  Lattice Boltzmann modeling of multiphase flows at large density ratio with an improved pseudopotential model.  Physical Review E 87, 053301\.
  * Li et al. (2020) Li, Q., Yu, Y., Wen, Z.X., 2020.  How does boiling occur in lattice boltzmann simulations?  Physics of Fluids 32, 093306\.
  * Li et al. (2018b) Li, Q., Yu, Y., Zhou, P., Yan, H., 2018b.  Enhancement of boiling heat transfer using hydrophilic-hydrophobic mixed surfaces: A lattice boltzmann study.  Applied Thermal Engineering 132, 490–499.
  * Liang and Mudawar (2018) Liang, G., Mudawar, I., 2018\.  Pool boiling critical heat flux (chf) – part 1: Review of mechanisms, models, and correlations.  International Journal of Heat and Mass Transfer 117, 1352–1367.
  * Lienhard et al. (1973) Lienhard, J.H., Dhir, V.K., Riherd, D.M., 1973.  Peak Pool Boiling Heat-Flux Measurements on Finite Horizontal Flat Plates.  Journal of Heat Transfer 95, 477–482.
  * Liu et al. (2016) Liu, H., Kang, Q., Leonardi, C.R., Schmieschek, S., Narváez, A., Jones, B.D., Williams, J.R., Valocchi, A.J., Harting, J., 2016.  Multiphase lattice Boltzmann simulations for porous media applications.  Computational Geosciences 20, 777–805.  URL: http://link.springer.com/10.1007/s10596-015-9542-3, doi:10.1007/s10596-015-9542-3.
  * Lou et al. (2013) Lou, Q., Guo, Z., Shi, B., 2013.  Evaluation of outflow boundary conditions for two-phase lattice boltzmann equation.  Physical Review E 87.
  * Luo (1998) Luo, L.S., 1998.  Unified theory of lattice Boltzmann models for nonideal gases.  Physical review letters 81, 1618\.
  * Lycett-Brown and Luo (2015) Lycett-Brown, D., Luo, K.H., 2015\.  Improved forcing scheme in pseudopotential lattice Boltzmann methods for multiphase flow at arbitrarily high density ratios.  Physical Review E 91, 23305\.  doi:10.1103/PhysRevE.91.023305.
  * Ma and Cheng (2019) Ma, X., Cheng, P., 2019.  3d simulations of pool boiling above smooth horizontal heated surfaces by a phase-change lattice Boltzmann method.  International Journal of Heat and Mass Transfer 131, 1095–1108.
  * McCullough et al. (2016) McCullough, J.W.S., Leonardi, C.R., Jones, B.D., Aminossadati, S.M., Williams, J.R., 2016.  Lattice Boltzmann methods for the simulation of heat transfer in particle suspensions.  International Journal of Heat and Fluid Flow 62, 150–165.
  * Qian et al. (1992) Qian, Y.H., d’Humières, D., Lallemand, P., 1992.  Lattice bgk models for navier-stokes equation.  EPL (Europhysics Letters) 17, 479–484.
  * Safaei et al. (2016) Safaei, M.R., Jahanbin, A., Kianifar, A., Gharehkhani, S., Kherbeet, A.S., Goodarzi, M., Dahari, M., 2016.  Mathematical modeling for nanofluids simulation: A review of the latest works, in: Modeling and Simulation in Engineering Sciences. InTech.  URL: https://doi.org/10.5772/64154, doi:10.5772/64154.
  * Safari et al. (2013) Safari, H., Rahimian, M.H., Krafczyk, M., 2013.  Extended lattice Boltzmann method for numerical simulation of thermal phase change in two-phase fluid flow.  Physical Review E - Statistical, Nonlinear, and Soft Matter Physics 88, 1–12.
  * Sbragaglia et al. (2007) Sbragaglia, M., Benzi, R., Biferale, L., Succi, S., Sugiyama, K., Toschi, F., 2007\.  Generalized lattice Boltzmann method with multirange pseudopotential.  Physical Review E 75, 026702\.
  * Shan and Chen (1993) Shan, X., Chen, H., 1993.  Lattice Boltzmann model for simulating flows with multiple phases and components.  Physical Review E 47, 1815–1819.
  * Shan and He (1998) Shan, X., He, X., 1998.  Discretization of the velocity space in the solution of the Boltzmann equation.  Physical Review Letters 80, 65\.
  * Son et al. (1999) Son, G., Dhir, V.K., Ramanujapu, N., 1999.  Dynamics and Heat Transfer Associated With a Single Bubble During Nucleate Boiling on a Horizontal Surface.  Journal of Heat Transfer 121, 623–631.
  * Son et al. (2002) Son, G., Ramanujapu, N., Dhir, V.K., 2002.  Numerical simulation of bubble merger process on a single nucleation site during pool nucleate boiling.  J. Heat Transfer 124, 51–62.
  * Swift et al. (1995) Swift, M.R., Osborn, W.R., Yeomans, J.M., 1995.  Lattice Boltzmann simulation of nonideal fluids.  Physical Review Letters 75, 830–833.
  * Wang and Dhir (1993) Wang, C.H., Dhir, V.K., 1993\.  Effect of Surface Wettability on Active Nucleation Site Density During Pool Boiling of Water on a Vertical Surface.  Journal of Heat Transfer 115, 659–669.
  * Welch and Wilson (2000) Welch, W.J., Wilson, J., 2000\.  A Volume of Fluid Based Method for Fluid Flows with Phase Change.  Journal of Computational Physics 160, 662–682.
  * Yu et al. (2005) Yu, H., Girimaji, S.S., Luo, L.S., 2005.  DNS and LES of decaying isotropic turbulence with and without frame rotation using lattice Boltzmann method.  Journal of Computational Physics 209, 599–616.
  * Yuan and Schaefer (2006) Yuan, P., Schaefer, L., 2006\.  Equations of state in a lattice Boltzmannn model.  Physics of Fluids 18, 042101\.
  * Zhai et al. (2017) Zhai, Q., Zheng, L., Zheng, S., 2017.  Pseudopotential lattice Boltzmann equation method for two-phase flow: A higher-order chapmann-enskog expansion.  Physical Review E 95, 023313\.
  * Zhang and Chen (2003) Zhang, R., Chen, H., 2003.  Lattice boltzmann method for simulations of liquid-vapor thermal flows.  Phys. Rev. E 67, 066711\.
  * Zhao et al. (2021) Zhao, W., Liang, J., Sun, M., Wang, Z., 2021.  Investigation on the effect of convective outflow boundary condition on the bubbles growth, rising and breakup dynamics of nucleate boiling.  International Journal of Thermal Sciences 167, 106877.  doi:10.1016/j.ijthermalsci.2021.106877.
  * Zheng et al. (2019) Zheng, L., Zheng, S., Zhai, Q., 2019.  Analysis of force treatment in lattice Boltzmann equation method.  International Journal of Heat and Mass Transfer 139, 747–750.  URL: https://doi.org/10.1016/j.ijheatmasstransfer.2019.05.059, doi:10.1016/j.ijheatmasstransfer.2019.05.059.
  * Zou and He (1997) Zou, Q., He, X., 1997.  On pressure and velocity boundary conditions for the lattice boltzmann BGK model.  Physics of Fluids 9, 1591–1598.