Analytical solution for dynamic evaporation of liquid in isothermal condition

Abstract

Report issue for preceding element

An analytical solution based on a diffuse interface model is presented for an isothermal evaporation problem under sub-saturation pressure. The macroscopic equations are derived from the free-energy method, widely recognized in the lattice Boltzmann literature, distinguishing our approach from conventional evaporation models that rely on jump conditions or pure kinetic theory. The interface behavior is fully described by differential equations, eliminating the need for assumptions such as local equilibrium at the interface. We derive an exact analytical solution for the inviscid case and propose an approximate solution when viscosity effects are considered. Our model unveils a novel relationship between evaporation rate and viscosity, providing new insights that have not been thoroughly explored in the literature. The analytical results are validated through numerical simulations using the open-source parallel library OpenLB, demonstrating excellent agreement in predicting the physical behavior of the evaporation phenomena within the framework of diffuse interface methods.

Report issue for preceding element

I Introduction

Report issue for preceding element

Understanding and predicting the evaporation phenomena is crucial for engineering applications, e.g. food dryingFathi et al. (2022), industrial dehydrationAl-Jammali et al. (2023), distillation processSkiborowski (2023), vacuum evaporationAbdelrahman et al. (2025), paint dryingDi Tullio et al. (2023), evaporative coolingAbdullah et al. (2023), steam power generationPolski et al. (2024), and others. There are different ways to study this phenomenon: analyticallyHertz (1882); Knudsen (1909); Hołyst et al. (2015); Persad and Ward (2016); Schrage (1953); Stefan (1891); Scriven (1995); Cossali and Tonini (2023); Álvarez-Hostos et al. (2024), numericallyMunicchi et al. (2022); Wang et al. (2023); Municchi et al. (2024), and experimentallyManova et al. (2022); Baptistella et al. (2023); Marchetto et al. (2024). For highly complex problems, numerical and experimental approaches are essential. However, the analytical method provides physical insights that lead to a deeper understanding of the phenomenon, making it highly valuable.

Report issue for preceding element

There is extensive literature on evaporation modeling. An important result is the Hertz–Knudsen relation Hertz (1882); Knudsen (1909); Hołyst et al. (2015); Persad and Ward (2016). Another significant contribution involving kinetic theory was made by Schrage Schrage (1953), who calculated evaporation fluxes of pure substances and multi-component systems. Notable results employing macroscopic phase change modeling include the one-dimensional Stefan problem Stefan (1891), originally developed for solid–liquid phase transition and Scriven’s solution Scriven (1995) for droplet evaporation. Recent articles extend these approaches to more complex situations Cossali and Tonini (2023); Álvarez-Hostos et al. (2024).

Report issue for preceding element

Another approach to modeling multiphase systems that has gained popularity in numerical studies is the diffuse interface model. This model treats the interface as a thin region in which fluid properties continuously transition between the vapor and liquid phases. The phase-field model Jacqmin (1999); Lamorgese et al. (2011); Kim (2012), which usually employs the Cahn–Hilliard Novick-Cohen (2008); Lee et al. (2014) or the Allen–Cahn Feng and Prohl (2003) equation, is a diffuse interface model used in combination with several numerical methods, for example: lattice Boltzmann methods (LBM) Wang et al. (2016), finite element methods (FEM) Feng et al. (2007), finite volume (FV) Qiao et al. (2015) or immersed boundary methods (IBM) Hua et al. (2014). Other examples of diffuse interface model-based methods are the free-energy Swift et al. (1996); Semprebon et al. (2016); Simonis et al. (2024) and the pseudo-potential Shan and Chen (1993, 1994) LBM, which are single-component methods capable of modeling liquid-vapor phase transition without any special treatment of the interface. Due to good parallelizability and robustness, LBM are generally well-suited for complex flow applications Ito et al. (2024); Bukreev et al. (2023a, 2024); Mink et al. (2022), even at large scales Kummerländer et al. (2024). Notably, many of the diffuse interface models listed above (or derivates thereof) are implemented in open-source libraries such as OpenLB Krause et al. (2020) and thus readily available for parallel execution on personal devices or clusters of any size.

Report issue for preceding element

While diffuse interface models such as Cahn–Hilliard are well-established in numerical and analytical studies of phase separationFoard and Wagner (2009), their use in analytical descriptions of liquid evaporation processes remains limited. However, these models also hold great potential for the development of analytical solutions since the fluid behavior is fully described by differential equations. This eliminates the need for assumptions such as local equilibrium at the interface, which is commonly used in sharp interface modeling. This represents a significant advancement in the analysis of non-equilibrium thermodynamics.

Report issue for preceding element

In this work, we employ a diffuse interface model to study the evaporation process. To simplify the problem, we consider an isothermal system subject to a sub–saturation pressure condition. A similar system was previously studied by Jamet Jamet (2004), but using a sharp interface approach. With our innovative methodology, we observed a dependence of the evaporation rate on the fluid viscosity – a behavior that has not been captured by any existing model in the literature. Finally, we derived an analytical solution for the evaporation rate in the inviscid case and an approximate solution for the viscous case. The predictions from our solutions are confirmed by simulations using an LBM implemented in the open source parallel library OpenLB Krause et al. (2020).

Report issue for preceding element

The practical applications of our analytical solution are as follows: Firstly, it provides predictions that can be tested in experiments to validate diffuse interface methods or suggest future corrections to this approach. Secondly, our solution can also be used to construct benchmark cases to test numerical approaches for diffuse interface models. Thirdly: predicting evaporation phenomena under sub-saturation pressure is crucial in, e.g., fuel injection systems, especially within low-pressure chambers. This process has significant implications for combustion enginesZeng et al. (2012). In future works, our approach of constructing an analytical solution can be extended to non–isothermal and multi–component problems.

Report issue for preceding element

The article is organized as follows: In Section II, the physical problem of isothermal evaporation under sub-saturation pressure is introduced. In Section III, the governing equations are presented. Section IV details the analytical solution procedure, including the derivation of the interface velocity and its dependence on physical parameters. In Section VI, the analytical solution is evidenced through comparisons with LBM simulations, and the influence of viscosity and pressure differences on the interface dynamics is analyzed. Finally, in Section VII, we summarize and discuss the main findings and propose directions for future work.

Report issue for preceding element

II Problem description

Report issue for preceding element

The problem of interest in this work is a dynamic isothermal liquid-vapor evaporation process similar to the one described by Jamet Jamet (2004). We consider a system with two phases (one liquid and one vapor) separated by a flat interface, which reduces to a one-dimensional case. The system has a fixed temperature $T$. In equilibrium, this system would have a homogeneous and constant pressure called saturation pressure $p_{\mathrm{sat}}(T)$. Our interest is to study what happens in this system when a pressure smaller than the saturation pressure is imposed in the vapor phase.

Report issue for preceding element

To illustrate this phenomenon, we use a simulation performed with the free-energy LBM formulation proposed by Wagner Wagner (2006). Figure 1 shows a representation of the simulated system. The top plot shows the density profile $\rho^{\ast}$ at two instants of time ($t^{\ast}=200$ and $t^{\ast}=600$). On the left and right boundaries (vapor phase) a pressure lower than the saturation pressure $p_{\text{v}}=0.99p_{\text{sat}}$ is imposed. The superscript $\cdot^{\ast}$ indicates dimensionless quantities such as $t^{\ast}$. The definition of the dimensionless quantities is presented in Section III. The simulation details are presented in Section V.

Report issue for preceding element

Figure 1: Top: density $\rho^{\ast}$ along domain size $x^{\ast}$ at two instants $t^{\ast}$. Bottom: velocity $u^{\ast}$ for the same conditions. Dimensionless quantities defined in (5). Simulations conducted with the model of Wagner Wagner (2006) implemented in OpenLB Krause et al. (2020) (see Section V for details).

Report issue for preceding element

We observe in the density profile of Figure 1 that the liquid region shrank over time, which means that the imposed boundary condition led to an evaporation process. The lower plot of Figure 1 shows the velocity profile. In the liquid region, the velocity is zero due to the problem symmetry. However, the velocity profile indicates that mass is leaving the domain through both boundaries.

Report issue for preceding element

This simulation was performed for a uniform dimensionless viscosity $\nu^{\ast}=1$. We repeated the simulation for different viscosity values and recorded the interface position for several other time instants. The results are summarized in Figure 2. Two very interesting behaviors are observed. First, for sufficiently long simulation times, the "position vs. time" curve is a straight line, which implies a constant interface velocity. The second surprising result is that the interface velocity depends on the fluid viscosity. This dependence of the evaporation rate on fluid viscosity is not present in classical evaporation modelsHertz (1882); Knudsen (1909); Schrage (1953) nor in analytical solutions to evaporation problemsStefan (1891); Scriven (1995). Showing, to the best of our knowledge, a behavior that has not yet been described in the literature.

Report issue for preceding element

Figure 2: Interface position $x_{i}^{\ast}$ over time $t^{\ast}$ for several viscosities $\nu^{\ast}$. Dimensionless quantities described in (5). Simulations conducted with the model of Wagner Wagner (2006) implemented in OpenLB Krause et al. (2020) (see Section V for details).

Report issue for preceding element

In the next sections, we focus on obtaining an analytical solution that describes this physical phenomenon.

Report issue for preceding element

III Governing equations

Report issue for preceding element

This problem is modeled by the one-dimensional (1D) mass and momentum conservation equations:

Report issue for preceding element

$$\frac{\partial\rho}{\partial t}+\frac{\partial(\rho u)}{\partial x}=0,$$ (1a)

$$\rho\frac{\partial u}{\partial t}+\rho u\frac{\partial u}{\partial x}=-\frac{\partial p}{\partial x}+\frac{\partial\tau_{xx}}{\partial x},$$ (1b)

where $\rho$ denotes the density, $u$ is the velocity, $p$ is the pressure, and $\tau_{xx}$ is the viscous stress as presented in Krüger et al.(Krüger et al., 2017, p. 6):

Report issue for preceding element

$$\tau_{xx}=\left(\frac{4}{3}\nu+\nu_{\text{b}}\right)\rho\frac{\partial u}{\partial x},$$ (2)

and, unless stated otherwise, all variables depend on time $t$ and space $x$. The symbols $\nu>0$ and $\nu_{\text{b}}$ in (2) represent the shear and bulk kinematic viscosity coefficients, respectively. For dilute monatomic gases, kinetic theory Vincenti et al. (1966) predicts $\nu_{\text{b}}=0$ which coincides with Stokes’ hypothesis Stokes (1845). Notably, the bulk viscosity for non-ideal fluids is a complex topic Graves and Argrow (1999); Kosuge and Aoki (2022). In addition, the stress-tensor behavior at the interface is barely described in the literature.

Report issue for preceding element

In particular, in this work, we compare the analytical solution of a problem governed by (1a) and (1b) with the numerical simulations performed with LBM based on an isothermal flow assumption. Thus, we use a different stress tensor definition than (2), i.e.

Report issue for preceding element

$$\tau_{xx}=2\nu\rho\frac{\partial u}{\partial x}.$$ (3)

The alternative form of the stress tensor given in Eq.(3), compared to the complete expression in (2), arises from the fact that the bulk viscosity predicted by the Bhatnagar–Gross–Krook Bhatnagar et al. (1954)–based lattice Boltzmann model is incompatible with the assumption of isothermal flow. Unless stated otherwise, we continue with the stress tensor defined in (3). The implications of this consideration can be tested, and if necessary, future corrections can be made.

Report issue for preceding element

Besides, we use a Korteweg-type pressureKorteweg (1901)

Report issue for preceding element

$$p=p_{\mathrm{EOS}}+\frac{\kappa}{2}\left(\frac{\partial\rho}{\partial x}\right)^{2}-\kappa\rho\frac{\partial^{2}\rho}{\partial x^{2}},$$ (4)

where $p_{\text{EOS}}$ denotes the equation of state (EOS) and $\kappa>0$ is the surface tension coefficient.

Report issue for preceding element

A fixed pressure $p_{\text{v}}$ is imposed on the boundary. Since the temperature is fixed, this will imply a fixed vapor density at the boundary $\rho_{\text{v}}$. Then, we define the following dimensionless quantities:

Report issue for preceding element

$$x^{\ast} \coloneqq\frac{x}{\rho_{\text{v}}}\sqrt{\frac{p_{\text{v}}}{\kappa}},~{}~{}~{}t^{\ast}\coloneqq t\frac{p_{\text{v}}}{\sqrt{\rho_{\text{v}}^{3}\kappa}},~{}~{}~{}u^{\ast}\coloneqq u\sqrt{\frac{\rho_{\text{v}}}{p_{\text{v}}}},$$ (5)

$$\nu^{\ast} \coloneqq\frac{\nu}{\sqrt{\rho_{\text{v}}\kappa}},~{}~{}~{}~{}~{}\rho^{\ast}\coloneqq\frac{\rho}{\rho_{\text{v}}},~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}p^{\ast}\coloneqq\frac{p}{p_{\text{v}}}.$$

IV Ansatz for a model solution

Report issue for preceding element

When performing initial simulations, we observed that the interface was moving with a constant speed (cf. Figure 2). Thus, in the ansatz used further below to derive an analytical solution, we make explicit use of the following assumptions:

Report issue for preceding element

1. After a transient time, the solution reaches a state where the interface propagates with constant velocity $u_{\text{int}}$.

   Report issue for preceding element

2. The center of the drop will have a zero velocity by symmetry and will retain a constant density.

   Report issue for preceding element

3. The previous assumptions suggest that we can pick a frame of reference ($x^{\prime},t^{\prime}$) for half of the system that moves with the interface, where any field $\phi$ is then independent of time, i.e. $\partial\phi/\partial t^{\prime}=0$.

   Report issue for preceding element

The transformation between the old and new frame of reference is:

Report issue for preceding element

$$\begin{aligned} x^{\prime}&=x-u_{\text{int}}t,\\ t^{\prime}&=t,\end{aligned}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\begin{aligned} x&=x^{\prime}+u_{\text{int}}t^{\prime},\\ t&=t^{\prime}.\end{aligned}$$ (6)

From the definitions in (5) and the consideration of Assumption 3, it follows that:

Report issue for preceding element

$$\frac{\partial\phi}{\partial t} =\left(\frac{\partial\phi}{\partial t^{\prime}}\right)\frac{\partial t^{\prime}}{\partial t}+\left(\frac{\partial\phi}{\partial x^{\prime}}\right)\frac{\partial x^{\prime}}{\partial t}\quad\Rightarrow\quad\frac{\partial\phi}{\partial t}=-u_{\text{int}}\frac{\partial\phi}{\partial x^{\prime}},$$ (7)

$$\frac{\partial\phi}{\partial x} =\left(\frac{\partial\phi}{\partial t^{\prime}}\right)\frac{\partial t^{\prime}}{\partial x}+\left(\frac{\partial\phi}{\partial x^{\prime}}\right)\frac{\partial x^{\prime}}{\partial x}\quad\Rightarrow\quad\frac{\partial\phi}{\partial x}=\frac{\partial\phi}{\partial x^{\prime}}.$$

IV.1 Mass conservation

Report issue for preceding element

Now, the mass conservation equation (1a) is rewritten in terms of $x^{\prime}$. Since the solution is independent of $t^{\prime}$, we write $\partial\phi/\partial x^{\prime}=d\phi/dx^{\prime}$ which gives

Report issue for preceding element

$$-u_{\text{int}}\frac{d\rho}{dx^{\prime}}+\frac{d(\rho u)}{dx^{\prime}}=0\quad\Rightarrow\quad(u_{\text{int}}-u)\frac{d\rho}{dx^{\prime}}=\rho\frac{du}{dx^{\prime}}.$$ (8)

We integrate (8) and obtain a relation between $\rho$ and $u$ up to a constant $C_{1}$:

Report issue for preceding element

$$\int\frac{d\rho}{\rho}=\int\frac{du}{u_{\text{int}}-u}+C_{1}\quad\Rightarrow\quad\rho(u_{\text{int}}-u)=\mathrm{e}^{C_{1}}=\text{const}.$$ (9)

According to Assumpion 2, in the liquid region $\rho=\rho_{\text{l}}$ and $u=0$, respectively. Then, from (9) we obtain

Report issue for preceding element

$$\rho(u_{\text{int}}-u)=\rho_{\text{l}}u_{\text{int}}.$$ (10)

Thus, our main finding from the mass conservation equation (1a) is (10) that leads to the useful relation:

Report issue for preceding element

$$\frac{du}{dx^{\prime}}=\frac{(u_{\text{int}}-u)}{\rho}\frac{d\rho}{dx^{\prime}}\quad\Rightarrow\quad\frac{du}{dx^{\prime}}=\frac{\rho_{\text{l}}u_{\text{int}}}{\rho^{2}}\frac{d\rho}{dx^{\prime}}$$ (11)

We make explicit use of (11) below.

Report issue for preceding element

IV.2 Momentum conservation

Report issue for preceding element

Considering the momentum equation (1b), we start by replacing the derivatives in terms of $t$ and $x$ with derivatives in terms of $x^{\prime}$, which gives

Report issue for preceding element

$$-\rho u_{\text{int}}\frac{du}{dx^{\prime}}+\rho u\frac{du}{dx^{\prime}}=-\frac{dp}{dx^{\prime}}+\frac{d}{dx^{\prime}}\left(2\nu\rho\frac{du}{dx^{\prime}}\right).$$ (12)

Next, we rewrite the left-hand side of (12) using (10) and (11) to obtain

Report issue for preceding element

$$\rho(u-u_{\text{int}})\frac{du}{dx^{\prime}}=-\frac{\rho_{\text{l}}^{2}u_{\text{int}}^{2}}{\rho^{2}}\frac{d\rho}{dx^{\prime}}.$$ (13)

Using (11), we also rewrite the viscous term as

Report issue for preceding element

$$2\nu\rho\frac{du}{dx^{\prime}}=2\rho_{\text{l}}u_{\text{int}}\frac{\nu}{\rho}\frac{d\rho}{dx^{\prime}}.$$ (14)

Now, all terms of (12) that were dependent on $u$ are replaced by terms that depend on $\rho$. Hence, we have reached a closed-form equation that depends on one variable only. All terms of the modified momentum conservation equation based on (12), (13), and (14) can be grouped in the form

Report issue for preceding element

$$\frac{d}{dx^{\prime}}\left(p+\frac{(\rho_{\text{l}}u_{\text{int}})^{2}}{\rho}-2\rho_{\text{l}}u_{\text{int}}\frac{\nu}{\rho}\frac{d\rho}{dx^{\prime}}\right)=0,$$ (15)

which gives

Report issue for preceding element

$$p+\frac{(\rho_{\text{l}}u_{\text{int}})^{2}}{\rho}-2\rho_{\text{l}}u_{\text{int}}\frac{\nu}{\rho}\frac{d\rho}{dx^{\prime}}=C_{2}=\text{const}.$$ (16)

Considering that density gradients are equal to zero in the bulk regions, we can evaluate the constant $C_{2}$ in (16) and obtain

Report issue for preceding element

$$C_{2}=p_{\text{v}}+\frac{\rho_{\text{l}}^{2}u_{\text{int}}^{2}}{\rho_{\text{v}}},\quad\text{and}~{}~{}~{}C_{2}=p_{\text{l}}+\frac{\rho_{\text{l}}^{2}u_{\text{int}}^{2}}{\rho_{\text{l}}}.$$ (17)

Next, we proceed by manipulating the differential equation in a more convenient form to solve.

Report issue for preceding element

IV.3 Final ordinary differential equation

Report issue for preceding element

From (4), (16) and (17) we have the following ordinary differential equation (ODE) in terms of $x^{\prime}$:

Report issue for preceding element

$$\kappa\rho\frac{d^{2}\rho}{dx^{\prime 2}}-\frac{\kappa}{2}\left(\frac{d\rho}{dx^{\prime}}\right)^{2}+2\rho_{\text{l}}u_{\text{int}}\frac{\nu}{\rho}\frac{d\rho}{dx^{\prime}}=f(\rho),$$ (18a)

where

$$f(\rho)=p_{\text{EOS}}-p_{\text{v}}-u_{\text{int}}^{2}\rho_{\text{l}}^{2}\left(\frac{1}{\rho_{\text{v}}}-\frac{1}{\rho}\right).$$ (18b)

The ODE (18a) can be converted into a more convenient form by applying a simplification based on the assumption that the density function $\rho=\rho(x^{\prime})$ is monotonically growing, i.e. $d\rho/dx^{\prime}>0$ for all $x^{\prime}$. Then, the dependency on $x^{\prime}$ can be transformed into a dependency on $\rho$.

Report issue for preceding element

Considering $d\rho/dx^{\prime}$ as a function of $\rho$, we have that

Report issue for preceding element

$$\frac{d\rho}{dx^{\prime}} =z(\rho),$$ (19)

$$\frac{d^{2}\rho}{dx^{\prime 2}} =\frac{dz}{dx^{\prime}}\quad\Rightarrow\quad\frac{d^{2}\rho}{dx^{\prime 2}}=\frac{dz}{d\rho}\frac{d\rho}{dx^{\prime}}\quad\Rightarrow\quad\frac{d^{2}\rho}{dx^{\prime 2}}=\frac{1}{2}\frac{dz^{2}}{d\rho}.$$

Also, the first two terms of the left-hand side of (18a) can be grouped such that

Report issue for preceding element

$$\kappa\rho\frac{d^{2}\rho}{dx^{\prime}}-\frac{\kappa}{2}\left(\frac{d\rho}{dx^{\prime}}\right)^{2}=\frac{\kappa\rho}{2}\frac{dz^{2}}{d\rho}-\frac{\kappa}{2}z^{2}$$ (20)

$$\quad\Rightarrow\quad\kappa\rho\frac{d^{2}\rho}{dx^{\prime}}-\frac{\kappa}{2}\left(\frac{d\rho}{dx^{\prime}}\right)^{2}=\frac{\rho^{2}}{2}\frac{d}{d\rho}\frac{z^{2}}{\rho}.$$

By applying the transformations of (20) to (18a), we obtain a non-linear ODE:

Report issue for preceding element

$$\frac{d}{d\rho}\frac{z^{2}}{\rho}+4\frac{\rho_{\text{l}}u_{\text{int}}}{\kappa}\frac{\nu}{\rho^{3}}z=\frac{2}{\kappa}\frac{f(\rho)}{\rho^{2}}.$$ (21)

The ODE (21) is an Abel’s equation of the second kind. This equation has the general formBougoffa (2010):

Report issue for preceding element

$$[g_{0}(\rho)+g_{1}(\rho)z]\frac{dz}{d\rho}=f_{2}(\rho)z^{2}+f_{1}(\rho)z+f_{0}(\rho).$$ (22)

Equation (21) is a particular case of (22) with coefficients:

Report issue for preceding element

$$g_{0}(\rho)=0,{} ~{}g_{1}(\rho)=\frac{2}{\rho}$$ (23)

$$f_{0}(\rho)=\frac{2}{\kappa}\frac{f(\rho)}{\rho^{2}},~{}~{}f_{1}(\rho)= -4\frac{\rho_{\text{l}}u_{\text{int}}}{\kappa}\frac{\nu}{\rho^{3}},~{}~{}f_{2}(\rho)=\frac{1}{\rho^{2}}.$$

Looking for a general solution to (21) may not be worth the effort, as the procedure could be too complex for practical applications. Therefore, we will proceed with an approximate solution.

Report issue for preceding element

IV.4 Exact solution of inviscid case

Report issue for preceding element

The simplest case occurs when there is no viscosity, $\nu=0$. Thus, (21) can be simplified

Report issue for preceding element

$$\frac{d}{d\rho}\frac{z^{2}}{\rho}=\frac{2}{\kappa}\frac{f(\rho)}{\rho^{2}}.$$ (24)

Next, we integrate the whole expression

Report issue for preceding element

$$\frac{z^{2}(\rho)}{\rho}-\frac{z_{\text{v}}^{2}}{\rho_{\text{v}}}=\int_{\rho_{\text{v}}}^{\rho}\frac{2}{\kappa}\frac{f(\rho)}{\rho^{2}}d\rho,$$ (25)

where $z_{\text{v}}=z(\rho_{\text{v}})$. The boundary condition $z_{\text{v}}=0$ (no density gradients in the bulk phases) is applied to obtain a solution for $z(\rho)$:

Report issue for preceding element

$$z^{(0)}(\rho)=\sqrt{\rho\int_{\rho_{\text{v}}}^{\rho}\frac{2}{\kappa}\frac{f(\rho)}{\rho^{2}}d\rho}.$$ (26)

The symbol $z^{(0)}$ represents the solution of $z$ for the inviscid case. The selection of the positive sign in the solution (26) means that the interface density grows with respect to the $x$ direction. In Figure 1 this represents the interface in the left-hand side.

Report issue for preceding element

Then, we replace $f(\rho)$ in (26) by (18b), and apply the boundary conditions $z_{L}=z(\rho_{L})=0$ (no density gradients in the bulk phases). After this step, we obtain the following expression:

Report issue for preceding element

$$0=\int_{\rho_{\text{v}}}^{\rho_{\text{l}}}\frac{p_{\text{EOS}}-p_{\text{v}}}{\rho^{2}}d\rho-\int_{\rho_{\text{v}}}^{\rho_{\text{l}}}\frac{u_{\text{int}}^{2}\rho_{\text{l}}^{2}}{\rho^{2}}\left(\frac{1}{\rho_{\text{v}}}-\frac{1}{\rho}\right)d\rho.$$ (27)

Finally, we perform the integral in the second term of the right-hand side of (27) and isolate $u_{\text{int}}$:

Report issue for preceding element

$$u_{\text{int}}^{(0)}=\frac{\sqrt{\int_{\rho_{\text{v}}}^{\rho_{\text{l}}}\frac{p_{\text{EOS}}-p_{\text{v}}}{\rho^{2}}d\rho}}{\sqrt{\frac{\rho_{\text{l}}}{\rho_{\text{v}}}\left(\frac{1}{2}\frac{\rho_{\text{l}}}{\rho_{\text{v}}}-1\right)+\frac{1}{2}}},$$ (28)

where $u_{\text{int}}^{(0)}$ represents the interface velocity for the inviscid case. The positive velocity is compatible with the orientation of the interface (26) (see Figure 1). The inviscid interface velocity is completely independent on $\kappa$. An important consideration regarding Eq.(28) is the dependence of the liquid density $\rho_{\text{l}}$ on the interface velocity. As a result, determining the correct liquid velocity requires solving a coupled system comprising Eq.(28) and the following two additional equations:

Report issue for preceding element

$$p_{\text{l}}=p_{\text{v}}+\rho_{\text{l}}u_{\text{int}}^{2}\left(\frac{1}{\rho_{\text{v}}}-\frac{1}{\rho_{\text{l}}}\right);\quad p_{l}=p_{EOS}(\rho_{\text{l}},T).$$ (29)

Next, we consider the case with viscous effects.

Report issue for preceding element

IV.5 Approximate solution of viscid case

Report issue for preceding element

The equivalent of (25) for the viscid case is:

Report issue for preceding element

$$\frac{z^{2}(\rho)}{\rho}-\frac{z_{\text{v}}^{2}}{\rho_{\text{v}}}=\int_{\rho_{\text{v}}}^{\rho}\left[\frac{2}{\kappa}\frac{f(\rho)}{\rho^{2}}-4\frac{\rho_{\text{l}}u_{\text{int}}}{\kappa}\frac{\nu}{\rho^{3}}z(\rho)\right]d\rho.$$ (30)

The problem with (30) is that we cannot use it to explicitly solve for $z$. To overcome this challenge and obtain an approximate solution, we consider that the viscosity $\nu$ has a small influence on the interface profile, which depends on $z$. At least, this approximation will be valid for sufficiently small $\nu$. Thus, we assume that $z$ in the viscous case can be approximated by $z^{(0)}$ from the inviscid case. We use this approximation to compute the integral of the second term on the right-hand side of (30).

Report issue for preceding element

Applying the boundary condition $z_{\text{l}}=0$ and the previous assumption to (30) yields

Report issue for preceding element

$$0=\int_{\rho_{\text{v}}}^{\rho_{\text{l}}}\left[\frac{2}{\kappa}\frac{f(\rho)}{\rho^{2}}-4\frac{\rho_{\text{l}}u_{\text{int}}}{\kappa}\frac{\nu}{\rho^{3}}z^{(0)}(\rho)\right]d\rho.$$ (31)

Replacing $f(\rho)$ in (26) by (18b) leads to a final equation that describes the interface velocity:

Report issue for preceding element

$$Au_{\text{int}}^{2}+Bu_{\text{int}}+C=0,$$ (32a)

$$A=\frac{\rho_{\text{l}}}{\rho_{\text{v}}}\left(1-\frac{1}{2}\frac{\rho_{\text{l}}}{\rho_{\text{v}}}\right)-\frac{1}{2},$$ (32b)

$$B=-2\rho_{\text{l}}\int_{\rho_{\text{v}}}^{\rho_{\text{l}}}\frac{\nu}{\rho^{3}}z^{(0)}d\rho,$$ (32c)

$$C=\int_{\rho_{\text{v}}}^{\rho_{\text{l}}}\frac{p_{\text{EOS}}-p_{\text{v}}}{\rho^{2}}d\rho.$$ (32d)

Notably, the solution approach above is valid even if $\nu$ depends on $\rho$. The interface velocity is obtained by simply solving the quadratic equation, (32a) toward

Report issue for preceding element

$$u_{\text{int}}=\frac{-B-\sqrt{B^{2}-4AC}}{2A},$$ (33)

the choice of sign in (33) is motivated by the following reasons:

Report issue for preceding element

• A is always negative ($A<0$) if $\rho_{\text{v}}<\rho_{\text{l}}$,

  Report issue for preceding element

• B is always negative ($B<0$) if $\rho_{\text{v}}<\rho_{\text{l}}$,

  Report issue for preceding element

• C is always positive ($C>0$) for $p_{\text{v}}<p_{\text{sat}}(T)$.

  Report issue for preceding element

Consequently, we conclude that a positive solution for $u_{i}$ is only possible for (33). The positive solution was selected due to the interface orientation choice in (26). The interface velocity must be solved together with the liquid density by coupling (33) with (29). The validity of the solution (33) is examined in Section VI.

Report issue for preceding element

V Lattice Boltzmann method and OpenLB implementation

Report issue for preceding element

In this work, we compare the analytical results given by (28) and (32a) with LBM simulations performed using OpenLB Krause et al. (2020). This is a powerful open-source LBM library applied in various fluid dynamics applications, including sub-grid particulate flows Bukreev et al. (2023b), fully resolved particle flows Hafen et al. (2023); Marquardt et al. (2024), turbulence simulations Siodlaczek et al. (2021); Simonis et al. (2022), optimization Reinke et al. (2022); Jeßberger et al. (2022), sub-grid multiphase flows Bukreev et al. (2023a), and fully resolved multiphase flows Simonis et al. (2024). The codes used in this work are implemented in OpenLB release 1.8 Kummerländer et al. (2025) and are available at https://www.openlb.net/download/.

Report issue for preceding element

The basic equation is named the lattice Boltzmann equation Krüger et al. (2017):

Report issue for preceding element

$$f_{i}(t+\Delta t,\bm{x}+\bm{c}_{i}\Delta t)-f_{i}(t,\bm{x})=\frac{\Delta t}{\tau}(f^{\mathrm{eq}}_{i}-f_{i})+\Delta tF_{i},$$ (34)

where $f_{i}$ and $f_{i}^{\mathrm{eq}}$ are the distribution function and its equilibrium counterpart. The indexes $i$ indicate the lattice velocity $\bm{c}_{i}$ in which $f_{i}$ is evaluated. The relaxation time $\tau$ is related to the fluid kinematic viscosity $\nu=c_{s}^{2}(\tau-0.5\Delta t)$. The parameter $c_{s}^{2}=(1/3)(\Delta x/\Delta t)^{2}$ is called lattice sound speed Krüger et al. (2017). The equilibrium distribution function is a function of the fluid density $\rho$ and equilibrium velocity $u_{\alpha}^{\mathrm{eq}}$Li et al. (2012):

Report issue for preceding element

$$f_{i}^{\mathrm{eq}}=w_{i}\rho\left(1+\frac{c_{i\alpha}}{c_{s}^{2}}u_{\alpha}^{\mathrm{eq}}+\frac{c_{i\alpha}c_{i\beta}-c_{s}^{2}\delta_{\alpha\beta}}{2c_{s}^{4}}\rho u_{\alpha}^{\mathrm{eq}}u_{\beta}^{\mathrm{eq}}\right),$$ (35)

where $w_{i}$ are the lattice weights for each lattice direction $i$. The equilibrium velocity is a quantity used to define $f_{i}^{\mathrm{eq}}$ and its relation with the real fluid velocity is introduced later.

Report issue for preceding element

The term $F_{i}$ is called forcing scheme and is responsible for the addition of an external force $F_{\alpha}$ to the LBM. Wagner Wagner (2006) proposed the following forcing scheme for the free-energy LBM:

Report issue for preceding element

$$F_{i}=w_{i}\left(\frac{c_{i\alpha}}{c_{s}^{2}}F_{\alpha}+\frac{c_{i\alpha}c_{i\beta}-c_{s}^{2}\delta_{\alpha\beta}}{2c_{s}^{4}}(F_{\alpha}u_{\beta}^{\mathrm{eq}}+F_{\beta}u_{\alpha}^{\mathrm{eq}}+\psi\delta_{\alpha\beta})\right),$$ (36)

with

Report issue for preceding element

$$\tau\psi=\left(\tau-\frac{1}{4}\right)\frac{F_{\alpha}F_{\alpha}}{\rho}+\frac{1}{12}\nabla^{2}\rho.$$ (37)

The Wagner forcing was validated only for 1D cases. Meaning that an extension for 2D cases was not provided yet. The thermodynamic force is related to the gradient of the chemical potential:

Report issue for preceding element

$$F_{\alpha}=-\rho\partial_{\alpha}\mu.$$ (38)

The spatial derivatives are computed using second order central finite difference stencils.

Report issue for preceding element

The velocity set is the standard two-dimensional nine velocities scheme (D2Q9):

Report issue for preceding element

$$\bm{c}_{i}=\begin{cases}(0,0),~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}i=0,\\ (c,0),(0,c),(-c,0),(0,-c),~{}~{}~{}~{}~{}~{}i=1,...,4,\\ (c,c),(-c,c),(-c,-c),(c,-c),~{}i=5,...,8.\\ \end{cases}$$ (39)

The macroscopic variables are computed from the moments of the distribution function:

Report issue for preceding element

$$\rho=\sum_{i}f_{i},$$ (40a)

$$\rho u_{\alpha}^{\mathrm{eq}}=\sum_{i}f_{i}c_{i\alpha}.$$ (40b)

The real fluid velocity $u_{\alpha}$ depends on $u_{\alpha}^{\mathrm{eq}}$ and the force $F_{\alpha}$:

Report issue for preceding element

$$u_{\alpha}=u_{\alpha}^{\mathrm{eq}}+\frac{\Delta t}{2}\frac{F_{\alpha}}{\rho}.$$ (41)

In this work, we use an equation of state based on the Landau free energy functional Briant et al. (2004):

Report issue for preceding element

$$p_{\mathrm{EOS}}=p_{c}(\nu_{\rho}+1)^{2}(3\nu_{\rho}^{2}-2\nu_{\rho}+1-2\tau_{w}),$$ (42a)

$$\nu_{\rho}=\frac{\rho-\rho_{c}}{\rho_{c}},~{}~{}~{}\tau_{w}=\frac{T_{c}-T}{T_{c}},$$ (42b)

where $\rho_{c}$, $p_{c}$ and $T_{c}$ are the critical density, pressure, and temperature. The chemical potential is described by:

Report issue for preceding element

$$\mu=\frac{4p_{c}}{\rho_{c}}\nu_{\rho}(\nu_{\rho}^{2}-\tau_{w})-\kappa\nabla^{2}\rho.$$ (43)

For a planar interface in equilibrium, the bulk densities ($\rho_{\text{v}}^{\text{sat}}$ and $\rho_{\text{l}}^{\text{sat}}$), interface thickness $\xi$, and surface tension $\gamma$ are given by:

Report issue for preceding element

$$\rho_{\text{v}}^{\text{sat}} =\rho_{c}(1-\tau_{w});~{}~{}\rho_{\text{l}}^{\text{sat}}=\rho_{c}(1+\tau_{w});$$ (44)

$$\xi =\sqrt{\frac{\kappa\rho_{c}^{2}}{4\tau_{w}p_{c}}};~{}~{}~{}~{}\gamma=\frac{4}{3}\sqrt{2\kappa p_{c}}(\tau_{w})^{3/2}\rho_{c}.$$

The superscripts in $\rho_{\text{v}}^{\text{sat}}$ and $\rho_{\text{l}}^{\text{sat}}$ mean that these densities are defined under saturation condition. We define the density ratio under saturation condition:

Report issue for preceding element

$$r_{\rho}^{\text{sat}}=\frac{\rho_{\text{l}}^{\text{sat}}}{\rho_{\text{v}}^{\text{sat}}}\quad\Rightarrow\quad r_{\rho}^{\text{sat}}=\frac{1+\tau_{w}}{1-\tau_{w}}.$$ (45)

The fluid temperature is related to $r_{\rho}^{\text{sat}}$. Thus, throughout this work, we opt to inform the reader $r_{\rho}^{\text{sat}}$ instead of the fluid temperature. For instance, our problem has only three independent parameters, which are dimensionless quantities: $p_{\mathrm{v}}/p_{\mathrm{sat}}$, $r_{\rho}^{\text{sat}}$ and $\nu^{\ast}$.

Report issue for preceding element

Considering a system of length $L^{\ast}$, we first initialize the density profile using the hyperbolic tangent function:

Report issue for preceding element

$$\rho(x^{\ast})=\rho_{\mathrm{v}}+\frac{\rho_{\mathrm{l}}-\rho_{\mathrm{v}}}{2}\left[\mathrm{tanh}\left(\frac{x^{\ast}_{1}}{\sqrt{2}\xi^{\ast}}\right)-\mathrm{tanh}\left(\frac{x^{\ast}_{2}}{\sqrt{2}\xi^{\ast}}\right)\right]$$ (46)

with $x_{1}^{\ast}=x^{\ast}-0.25L^{\ast}$, $x_{2}^{\ast}=x^{\ast}-0.75L^{\ast}$ and $\xi^{\ast}=\xi\sqrt{p_{\mathrm{v}}}/(\rho_{\mathrm{v}}\sqrt{\kappa})$.

Report issue for preceding element

The initialization of the velocity profile is discussed later. After defining the macroscopic quantities, the distribution function is initialized to be equal to its equilibrium counterpart (35). Periodic boundary conditions are applied in the bottom and top to reproduce the 1D case. In the side boundaries, the chemical potential and the density are fixed for the vapor value.

Report issue for preceding element

There is a key point to discuss regarding viscosity. We have observed that in some simulations, particularly at low viscosity, oscillations in the density and velocity profiles can arise during the initial transient regime. To ensure that these oscillations are quickly damped, we assign a relaxation time of 1 in the bulk region while keeping the desired relaxation time at the interface. The higher viscosity in the bulk rapidly suppresses these oscillations, allowing the system to converge more quickly to a steady interface velocity. Importantly, this approach has no impact on the final result.

Report issue for preceding element

The final velocity of the interface should not depend on the initial conditions of the system but rather on the boundary conditions and fluid properties. To verify this, we simulated Cases 1 and 2, as shown in Figure 3. In Case 1, the system is initialized with the velocity profile obtained from the analytical solution (33) and (10). In Case 2, the system is initialized with zero velocity. It is observed that, although the transient period differs between the two cases, the final interface velocity remains entirely independent of the initial condition. Therefore, we chose to initialize the system using the analytical solution to achieve a shorter transient period.

Report issue for preceding element

Figure 3: Interface velocity $u_{\mathrm{int}}^{\ast}$ along time $t^{\ast}$ for LBM simulations in OpenLB Krause et al. (2020). All cases runned with $r_{\rho}^{\text{sat}}=32$, $p_{\mathrm{v}}=0.95p_{\mathrm{sat}}$ and $\nu^{\ast}=1$. Case 1: $L^{\ast}=83$, initial velocity equal to analytical solution (33). Case 2: $L^{\ast}=83$, initial velocity equal to zero. Case 3: $L^{\ast}=166$, initial velocity equal to zero.

Report issue for preceding element

Furthermore, even when initializing with the analytical solution, a slight oscillation in the interface velocity can be observed, as recorded in Figure 3. This occurs because $f_{i}$ is initialized equal to $f_{i}^{\mathrm{eq}}$, which is only an approximation rather than the exact condition.

Report issue for preceding element

As previously mentioned, the analytical solution depends solely on the boundary conditions and fluid properties. Consequently, the final interface velocity should be entirely independent of the domain size $L^{\ast}$. This behavior can be observed in the comparison between Cases 2 and 3 in Figure 3. Case 3 shares the same simulation conditions as Case 2 but with a domain size $L^{\ast}$ that is twice as large. Despite differences in the transient regime, the final interface velocity is identical in both cases.

Report issue for preceding element

Figure 4: Interface velocity $u_{\mathrm{int}}^{\ast}$ dependency on interface resolution $\xi/\Delta x$ for LBM simulations in OpenLB Krause et al. (2020). Results at $t^{\ast}=200$ with $L^{\ast}=83$, $r_{\rho}^{\mathrm{sat}}=32$, $p_{\mathrm{v}}=0.95p_{\mathrm{sat}}$ and $\nu^{\ast}=1$.

Report issue for preceding element

Finally, before presenting the results, we conducted a mesh study to determine the appropriate resolution for the simulations. To this end, we ran a case with $r_{\rho}^{\mathrm{sat}}=32$, $p_{\mathrm{v}}=0.99p_{\mathrm{sat}}$ and $\nu^{\ast}=1$. Since the domain length is not relevant, we used the interface resolution, given by $\xi/\Delta x$, as a measure of the simulation resolution. This value provides an estimate of the number of nodes that occupy the interface. The results are presented in Figure 4. The difference in $u_{\mathrm{int}}$ between resolutions 10 and 20 was only 0.027%, indicating that a resolution of 10 is appropriate for simulations.

Report issue for preceding element

VI Results

Report issue for preceding element

VI.1 Verifying modeling validity

Report issue for preceding element

Before testing the inviscid case analytical solution (28) and the viscid case approximate solution (33), we validate the physical equations that lead to our solution. The validations are based on comparisons with numerical simulations using the free-energy LBM proposed by WagnerWagner (2006) and implemented in OpenLBKrause et al. (2020). Based on the results, we show that the effects described here are actually induced by the physical equations.

Report issue for preceding element

Equation (16) is crucial to our work, since it defines the changes of the pressure, given in (4), due to momentum exchange and friction forces acting at the interface. To evaluate the effect of each of those terms, we adopt the following definitions:

Report issue for preceding element

$$p_{\text{non-ideal}}\coloneqq p_{\mathrm{EOS}}+\frac{\kappa}{2}\left(\frac{\partial\rho}{\partial x^{\prime}}\right)^{2}-\kappa\rho\frac{\partial^{2}\rho}{\partial x^{\prime 2}},$$ (47a)

$$p_{\text{momentum}}\coloneqq p_{\mathrm{v}}+\rho_{\mathrm{l}}^{2}u_{\mathrm{int}}^{2}\left(\frac{1}{\rho_{\mathrm{v}}}-\frac{1}{\rho}\right),$$ (47b)

$$p_{\text{friction}}\coloneqq p_{\mathrm{v}}+2\rho_{\mathrm{l}}u_{\mathrm{int}}{\rho}\frac{d\rho}{dx^{\prime}},$$ (47c)

$$p_{\text{total}}\coloneqq p_{\mathrm{v}}+\rho_{\mathrm{l}}^{2}u_{\mathrm{int}}^{2}\left(\frac{1}{\rho_{\mathrm{v}}}-\frac{1}{\rho}\right)+2\rho_{\mathrm{l}}u_{\mathrm{int}}\frac{\nu}{\rho}\frac{d\rho}{dx^{\prime}}.$$ (47d)

We perform an LBM simulation using $r_{\rho}^{\mathrm{sat}}=2$, $p_{\mathrm{v}}=0.99p_{\mathrm{sat}}$, and $\nu^{\ast}=0.5$, running the simulation until a steady-state velocity solution is reached. From the resulting density and velocity fields, we compute (47a), (47b), (47c), and (47d). The first- and second-order derivatives in these equations are approximated using central finite differences with second-order accuracy.

Report issue for preceding element

Figure 5: Pressure $p^{\ast}$ (left vertical axis) and density $\rho^{\ast}$ (right vertical axis) along domain length $x^{\ast}$. Dimensionless quantities are defined in (5). Pressure definitions are given in (47). Simulation results with $r_{\rho}^{\mathrm{sat}}=2$, $p_{\mathrm{v}}=0.99p_{\mathrm{sat}}$ and $\nu^{\ast}=1$ obtained from LBM implemented in OpenLBKrause et al. (2020).

Report issue for preceding element

In Figure 5 we show the density profile for this simulation and how $p_{\text{non-ideal}}$ and $p_{\text{total}}$ change along the interface region. Both pressures are equivalent which is consistent with (16). We observe that going from the vapor region and crossing the interface, $p_{\text{non-ideal}}$ grows, reaching a peak inside the interface region and then decreasing towards the liquid region. The value of $p_{\text{non-ideal}}$ is larger inside the liquid region in comparison with the vapor region.

Report issue for preceding element

In Figure 5 we observe that $p_{\text{momentum}}$ grows continuously from the vapor region to the liquid region due to momentum exchange across the interface. In the bulk regions $p_{\text{non-ideal}}$ is equal to $p_{\text{momentum}}$ indicating that the momentum exchange is the only responsible for the pressure difference between bulk phases.

Report issue for preceding element

Finally, we also plot $p_{\text{friction}}$ in Figure 5. The value of this quantity grows inside the interface due to the velocity gradient across the interface. However, the friction does not contribute to change the pressure between the two bulk phases. The numerical results support our physical model. Any possible physical inconsistency of this solution will be related to the diffuse interface model and not with any further consideration.

Report issue for preceding element

VI.2 Comparison with lattice Boltzmann method

Report issue for preceding element

We compare the analytical approximation given by (33) against results obtained from LBM simulations. To assess the model’s robustness, we consider three representative density ratios, $r_{\rho}^{\text{sat}}=2$, $8$, and $32$. We then analyze how the interface velocity $u_{\mathrm{int}}^{\ast}$ varies as a function of the kinematic viscosity $\nu^{\ast}$. Simulations are performed across a range of pressure ratios defined as $f_{\text{p}}=p_{\text{v}}/p_{\text{sat}}$.

Report issue for preceding element

It is important to note that the viscid solution given by (33) is an approximation and may deviate from the exact solution of the underlying differential equation. However, the free-energy LBM is a consistent numerical method for solving the mass and momentum conservation equations, and its results converge to the exact solution as the spatial and temporal resolution increases. Therefore, we use high-resolution LBM simulations as a reference to assess the accuracy of our analytical approximation.

Report issue for preceding element

Figure 6: Dependency of interface velocity $u^{\ast}_{\mathrm{int}}$ on fluid kinematic viscosity $\nu^{\ast}$ for fixed $r_{\rho}^{\text{sat}}$ and different values of $f_{\text{p}}=p_{\mathrm{v}}/p_{\text{sat}}$. Dimensionless quantities are defined in (5). Points represent simulation results using LBM implemented in OpenLBKrause et al. (2020). Lines represent approximate solution (33).

Report issue for preceding element

The results of this comparison are shown in Figure 6. The overall behavior of the solution is qualitatively similar across all cases. For very low viscosities, the solution converges to the inviscid limit given by (28). Conversely, when $\nu^{\ast}$ is sufficiently large, the interface velocity becomes small and the quadratic term $u_{\text{int}}^{2}$ in (32a) can be neglected. In this high-viscosity regime, the interface velocity exhibits an inverse dependence on the viscosity.

Report issue for preceding element

A selection of relative errors between the analytical approximation and the LBM results is compiled in Table 1. Three main trends can be identified from the data: First, the approximation becomes less accurate as the viscosity $\nu^{\ast}$ increases. This is consistent with the construction of our solution, which is derived from the inviscid limit and, therefore, performs best in low-viscosity regimes. Second, the error also increases as the pressure fraction $f_{\text{p}}=p_{\text{v}}/p_{\text{sat}}$ decreases. A third trend is observed with respect to the density ratio. For lower values such as $r_{\rho}^{\text{sat}}=2$, the approximation becomes significantly less accurate. This suggests that, near the critical point—where the density ratio approaches unity—the interface becomes more sensitive to viscous effects. Since our approximation relies on the assumption of proximity to the inviscid case, its accuracy degrades in this regime, where the influence of viscosity becomes more pronounced.

Report issue for preceding element

Table 1: Relative error between the analytical approximation and LBM simulations for different density ratios $r_{\rho}^{\text{sat}}$, pressure fractions $f_{\text{p}}$, and viscosities $\nu^{\ast}$.

$r_{\rho}^{\text{sat}}$	$f_{\text{p}}$	$\nu^{\ast}$	Relative error (%)
32	0.99	0.025	0.02
32	0.99	6.4	0.28
32	0.90	0.025	0.35
32	0.90	6.4	2.91
8	0.99	0.025	0.024
8	0.99	6.4	0.26
8	0.90	0.025	0.10
8	0.90	6.4	3.10
2	0.99	0.025	0.20
2	0.99	6.4	10.02
2	0.97	0.2	3.10
2	0.97	6.4	24.81

Report issue for preceding element

These observations help delineate the validity range of the proposed approximation. The solution consistently yields accurate results in the low-viscosity regime. Moreover, we observe that if the density ratio remains above 8 and the pressure fraction exceeds 0.90, the approximation remains reliable across a broad range of viscosities. For $r_{\rho}^{\text{sat}}=2$, however, some low-viscosity cases could not be simulated due to stability limitations inherent to the BGK collision operator at small relaxation times.

Report issue for preceding element

VI.3 Comparison with sharp interface solution

Report issue for preceding element

JametJamet (2004) modeled a similar problem using a sharp-interface approach. His formulation is based on the framework originally proposed by IshiiIshii (1975), which derives from applying jump conditions across the interface under the assumption of thermodynamic equilibrium. Under this assumption, the following equation is obtained (see Ishii(Ishii, 1975, p.41, Eq.(2-106)); Jamet(Jamet, 2004, Eq. (16))):

Report issue for preceding element

$$p_{\text{v}}=p_{\text{sat}}-\frac{1}{2}\frac{\rho_{v}\rho_{l}}{\rho_{l}-\rho_{v}}(u_{v}-u_{l})^{2}.$$ (48)

Assuming symmetry in our test case, such that $u_{\text{l}}=0$, and using (10), we obtain:

Report issue for preceding element

$$u_{\text{int}}=\sqrt{\frac{p_{\text{sat}}-p_{\text{v}}}{\frac{\rho_{\text{l}}^{2}}{2}\left(\frac{1}{\rho_{\text{v}}}-\frac{1}{\rho_{\text{l}}}\right)}}.$$ (49)

This relation shares similarities with (28). In fact, (49) can be derived from (28) if we assume:

Report issue for preceding element

$$\int_{\rho_{v}}^{\rho_{l}}\frac{p_{\text{EOS}}}{\rho^{2}}\,d\rho=\int_{\rho_{v}}^{\rho_{l}}\frac{p_{\text{sat}}}{\rho^{2}}\,d\rho.$$ (50)

This is essentially the Maxwell construction, which holds under thermodynamic equilibrium. However, under non-equilibrium conditions, this approximation may no longer be valid. Moreover, this macroscopic model does not capture the dependence of the evaporation flux on viscosity.

Report issue for preceding element

Figure 7: Dependency of interface velocity $u^{\ast}_{\mathrm{int}}$ on the density ratio $r_{\rho}^{\text{sat}}$ for different values of $f_{\text{p}}=p_{\mathrm{v}}/p_{\text{sat}}$. Plot shows a comparison between our inviscid solution (28) and Jamet solution (49).

Report issue for preceding element

Next, we compare our inviscid solution (28) with the sharp-interface model given by (49). The results for different values of $r_{\rho}^{\text{sat}}$ and $f_{\text{p}}$ are shown in Figure 7. Overall, the two solutions exhibit excellent agreement. Noticeable deviations are observed only when $f_{\text{p}}$ is significantly reduced (e.g., to 0.7), or near the critical point, where $r_{\rho}^{\text{sat}}$ approaches unity.

Report issue for preceding element

For general purposes, we conclude that the two models yield comparable predictions across a wide range of conditions. In future work, we intend to extend this comparison by incorporating LBM simulations using more stable collision operators. At present, the use of the BGK operator restricts us from accessing sufficiently low viscosity to fully explore the inviscid regime.

Report issue for preceding element

VII Conclusion

Report issue for preceding element

This work presents an analytical approximation to the problem of isothermal evaporation under sub-saturation pressure conditions using a diffuse interface model. For the inviscid case, our solution is exact. By avoiding traditional assumptions such as local thermodynamic equilibrium at the interface, the proposed approach offers a new perspective on the relationship between evaporation rates and viscosity. The derived solution is evidenced against numerical simulation results produced with an LBM implemented in OpenLB, demonstrating excellent agreement and reinforcing its physical consistency.

Report issue for preceding element

The findings highlight the potential of diffuse interface models in studying phase change phenomena, not only as a predictive tool but also as a benchmark for numerical methods. The results provide clear evidence that the evaporation dynamics, including interface velocity, are strongly influenced by fluid viscosity, which is often overlooked in existing literature. Furthermore, this found analytical solution can serve as a reference for experimental validation of diffuse interface methodologies or as a basis for future refinements of these models.

Report issue for preceding element

Future work could extend this approach to non-isothermal conditions and multi-component systems, enabling broader applications in practical scenarios such as fuel atomization and heat transfer processes. By bridging the gap between analytical models and numerical simulations, this study contributes to a more comprehensive and accurate modeling of phase-change phenomena in complex systems.

Report issue for preceding element