# Subgrid moving contact line model for direct numerical simulations of bubble dynamics in pool boiling of pure fluids

Vadim S Nikolayev, Linkai Wei, Guillaume Bois

► To cite this version:

Vadim S Nikolayev, Linkai Wei, Guillaume Bois. Subgrid moving contact line model for direct numerical simulations of bubble dynamics in pool boiling of pure fluids. Journal of Physics: Conference Series, 2024, 9th European Thermal Sciences Conference (Eurotherm 2024), 2766 (1), pp.012123. 10.1088/1742-6596/2766/1/012123. cea-04602109

HAL Id: cea-04602109

https://cea.hal.science/cea-04602109

Submitted on 5 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Subgrid moving contact line model for direct numerical simulations of bubble dynamics in pool boiling of pure fluids

Vadim S Nikolayev, Linkai Wei, Guillaume Bois

► To cite this version:

Vadim S Nikolayev, Linkai Wei, Guillaume Bois. Subgrid moving contact line model for direct numerical simulations of bubble dynamics in pool boiling of pure fluids. Journal of Physics: Conference Series, 2024, 9th European Thermal Sciences Conference (Eurotherm 2024), 2766 (1), pp.012123. 10.1088/1742-6596/2766/1/012123. cea-04602109

HAL Id: cea-04602109

https://cea.hal.science/cea-04602109v1

Submitted on 5 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Subgrid moving contact line model for direct numerical simulations of bubble dynamics in pool boiling of pure fluids

Vadim S. Nikolayev<fnref n="1" />, Linkai Wei<fnref n="2" /> and Guillaume Bois<fnref n="2" />

E-mail: vadim.nikolayev@cea.fr

**Abstract.** This contact line vicinity model is conceived as a subgrid model for the DNS of bubble growth in boiling. The model is based on the hydrodynamic multiscale theory and is suitable for the partial wetting case. On the smallest length scale (distance from the contact line) ~100 nm, the interface slope is controlled by the Voinov angle. It is the static apparent contact angle (ACA) that forms due to evaporation, similarly to previous models neglecting the contact line motion. The calculation of the Voinov angle is performed with the generalized lubrication approximation and includes several nanoscale effects like those of Kelvin and Marangoni, vapor recoil, hydrodynamic slip length and interfacial kinetic resistance. It provides the finite values of the heat flux, pressure and temperature at the contact line. The dynamic ACA is obtained with the Cox-Voinov formula. The microscopic length of the Cox-Voinov formula (Voinov length) is controlled mainly by the hydrodynamic slip. The integral heat flux passing through the contact line vicinity is almost independent of the nanoscale phenomena, with the exception of the interfacial kinetic resistance and is mostly defined by the dynamic ACA. Both the dynamic ACA and the integral heat flux are the main output parameters of the subgrid model, while the local superheating and the microscopic contact angle are the main input parameters. The model is suitable for the grid sizes > 1 µm.

## 1. Introduction

In nucleate boiling, the vapor bubble nucleation and growth on the heating wall favors high heat fluxes from the surface to the liquid thus providing a high rate of overall heat exchange. The bubble growth modelling is a longstanding challenge and is still debated in the literature. In the modeling of bubble growth, the line of triple vapor-liquid-solid contact (CL) presents a separate issue because of its geometrical singularity that has two main consequences for its vicinity loosely called microregion. First, high rates of the heat and mass exchange in the microregion need to be included into the model to reproduce correctly the wall cooling at the CL [1] and bubble growth rate [2]. Second, strong hydrodynamic flows caused by the mass exchange cause strong variation of interface slope at the bubble base (called apparent contact angle $\theta_{app}$) [3]. It needs to be correctly modeled as it presents the geometrical boundary condition that affects the bubble shape.

Several versions of microregion models have been developed since the pioneering works [4, 5], most of them considering the static CL and the complete wetting case where the long-range

<footnote n="1">¹ Université Paris-Saclay, CEA, SPEC, CNRS, 91191 Gif-sur-Yvette Cedex, France</footnote>

<footnote n="2">² Université Paris-Saclay, CEA, Service de Thermo-hydraulique et de Mécanique des Fluides, 91191 Gif-sur-Yvette Cedex, France</footnote>
Subgrid moving contact line model for direct
numerical simulations of bubble dynamics in pool
boiling of pure fluids

Vadim S. Nikolayev¹, Linkai Wei² and Guillaume Bois²

¹ Université Paris-Saclay, CEA, SPEC, CNRS, 91191 Gif-sur-Yvette Cedex, France

² Université Paris-Saclay, CEA, Service de Thermo-hydraulique et de Mécanique des Fluides, 91191 Gif-sur-Yvette Cedex, France

E-mail: vadim.nikolayev@cea.fr

**Abstract.** This contact line vicinity model is conceived as a subgrid model for the DNS of bubble growth in boiling. The model is based on the hydrodynamic multiscale theory and is suitable for the partial wetting case. On the smallest length scale (distance from the contact line) ~ 100 nm, the interface slope is controlled by the Voinov angle. It is the static apparent contact angle (ACA) that forms due to evaporation, similarly to previous models neglecting the contact line motion. The calculation of the Voinov angle is performed with the generalized lubrication approximation and includes several nanoscale effects like those of Kelvin and Marangoni, vapor recoil, hydrodynamic slip length and interfacial kinetic resistance. It provides the finite values of the heat flux, pressure and temperature at the contact line. The dynamic ACA is obtained with the Cox-Voinov formula. The microscopic length of the Cox-Voinov formula (Voinov length) is controlled mainly by the hydrodynamic slip. The integral heat flux passing through the contact line vicinity is almost independent of the nanoscale phenomena, with the exception of the interfacial kinetic resistance and is mostly defined by the dynamic ACA. Both the dynamic ACA and the integral heat flux are the main output parameters of the subgrid model, while the local superheating and the microscopic contact angle are the main input parameters. The model is suitable for the grid sizes > 1 µm.

# 1. Introduction

In nucleate boiling, the vapor bubble nucleation and growth on the heating wall favors high heat fluxes from the surface to the liquid thus providing a high rate of overall heat exchange. The bubble growth modelling is a longstanding challenge and is still debated in the literature. In the modeling of bubble growth, the line of triple vapor-liquid-solid contact (CL) presents a separate issue because of its geometrical singularity that has two main consequences for its vicinity loosely called microregion. First, high rates of the heat and mass exchange in the microregion need to be included into the model to reproduce correctly the wall cooling at the CL [1] and bubble growth rate [2]. Second, strong hydrodynamic flows caused by the mass exchange cause strong variation of interface slope at the bubble base (called apparent contact angle $\theta_{app}$) [3]. It needs to be correctly modeled as it presents the geometrical boundary condition that affects the bubble shape.

Several versions of microregion models have been developed since the pioneering works [4, 5], most of them considering the static CL and the complete wetting case where the long-range
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Schematics of the curved interface in the contact line (CL) vicinity in the reference of the immobile CL (the wall is moving with the speed $U_{cl}$).</figcaption></figure>

molecular interaction are implemented through the disjoining pressure. However, in most practical cases of boiling, the wetting is partial and CL moves. The microregion model for the partial wetting case has been developed by Janeček and Nikolayev [6, 7]. It includes several nanoscale effects like those of Kelvin and Marangoni, vapor recoil, hydrodynamic slip length and interfacial kinetic resistance. Since such a model provides the finiteness of temperature, pressure and heat flux, its numerical implementation is robust. Its version suitable for description of large $\theta_{app}$ has been developed based on the extended lubrication approach [8–10]. The analysis with the moving CL [10–12] shows that, while the interface shape is defined by evaporation-induced liquid flow at the distance $\lesssim 100$ nm from CL, the flow caused by the CL motion is dominant at larger length scales. Thanks to such a scale separation, the effects of evaporation and CL motion can be considered separately. Elaboration of the subgrid model for bubble growth simulations based on this idea is the objective of this study.

## 2. Generalized lubrication approach

Since the microregion size is typically much smaller than the CL radius, the subgrid problem can be solved in 2D, even for DNS of the 3D or axisymmetric bubble. Instead of operating with the Cartesian coordinates $(x, h)$ of a point A belonging to the interface (Fig. 1), one uses a curvilinear coordinate $s$ that runs along the interface, with $s = 0$ at the CL. Therefore, the following geometrical relations hold:

$$ 
\begin{aligned}
\frac{dh}{ds} &= \sin \phi, \\
\frac{dx}{ds} &= \cos \phi,
\end{aligned}
\qquad (1) $$

where $\phi$ is the local interface slope. For the point A, one introduces a straight wedge formed by the tangent to the interface and the wall and the corresponding circular arc passing through A with the length

$$ \zeta = h\phi/\sin\phi. \qquad (2) $$

The generalized lubrication equation, written in the reference where CL is immobile but the wall moves with the CL velocity $U_{cl}$ in the opposite direction, reads [10]

$$ \frac{d}{ds} \left\{ \frac{1}{\mu G(\phi)} \left[ \frac{\zeta}{2}(\zeta + 2l_s) \frac{d\sigma}{ds} + \frac{\zeta^2}{3}(\zeta + 3l_s) \frac{d\Delta p}{ds} \right] - U_{cl} \zeta \frac{F(\phi)}{G(\phi)} \right\} = -\frac{J}{\rho_l}, \qquad (3) $$

where $\Delta p$ is the interfacial vapor-liquid pressure jump, $J$, the evaporation mass flux, $l_s$, the hydrodynamic slip, $\sigma$, the surface tension, $\mu$, liquid shear viscosity, $\rho_l$, the liquid density, and the coefficients

$$ F(\phi) = \frac{2\phi^2}{3} \frac{\sin \phi}{\phi - \sin \phi \cos \phi}, \qquad (4) $$

$$ G(\phi) = \frac{\phi^3}{3} \frac{4}{\sin \phi \cos \phi - \phi \cos 2\phi} \qquad (5) $$
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Schematics of the curved interface in the contact line (CL) vicinity in the reference of the immobile CL (the wall is moving with the speed $U_{cl}$).</figcaption></figure>

molecular interaction are implemented through the disjoining pressure. However, in most practical cases of boiling, the wetting is partial and CL moves. The microregion model for the partial wetting case has been developed by Janeček and Nikolayev [6, 7]. It includes several nanoscale effects like those of Kelvin and Marangoni, vapor recoil, hydrodynamic slip length and interfacial kinetic resistance. Since such a model provides the finiteness of temperature, pressure and heat flux, its numerical implementation is robust. Its version suitable for description of large $\theta_{app}$ has been developed based on the extended lubrication approach [8–10]. The analysis with the moving CL [10–12] shows that, while the interface shape is defined by evaporation-induced liquid flow at the distance $\lesssim 100$ nm from CL, the flow caused by the CL motion is dominant at larger length scales. Thanks to such a scale separation, the effects of evaporation and CL motion can be considered separately. Elaboration of the subgrid model for bubble growth simulations based on this idea is the objective of this study.

## 2. Generalized lubrication approach

Since the microregion size is typically much smaller than the CL radius, the subgrid problem can be solved in 2D, even for DNS of the 3D or axisymmetric bubble. Instead of operating with the Cartesian coordinates $(x, h)$ of a point A belonging to the interface (Fig. 1), one uses a curvilinear coordinate $s$ that runs along the interface, with $s = 0$ at the CL. Therefore, the following geometrical relations hold:

$$ 
\begin{aligned}
\frac{dh}{ds} &= \sin \phi, \\
\frac{dx}{ds} &= \cos \phi,
\end{aligned}
\qquad (1) $$

where $\phi$ is the local interface slope. For the point A, one introduces a straight wedge formed by the tangent to the interface and the wall and the corresponding circular arc passing through A with the length

$$ \zeta = h\phi/\sin\phi. \qquad (2) $$

The generalized lubrication equation, written in the reference where CL is immobile but the wall moves with the CL velocity $U_{cl}$ in the opposite direction, reads [10]

$$ \frac{d}{ds} \left\{ \frac{1}{\mu G(\phi)} \left[ \frac{\zeta}{2}(\zeta + 2l_s) \frac{d\sigma}{ds} + \frac{\zeta^2}{3}(\zeta + 3l_s) \frac{d\Delta p}{ds} \right] - U_{cl} \zeta \frac{F(\phi)}{G(\phi)} \right\} = -\frac{J}{\rho_l}, \qquad (3) $$

where $\Delta p$ is the interfacial vapor-liquid pressure jump, $J$, the evaporation mass flux, $l_s$, the hydrodynamic slip, $\sigma$, the surface tension, $\mu$, liquid shear viscosity, $\rho_l$, the liquid density, and the coefficients

$$ F(\phi) = \frac{2\phi^2}{3} \frac{\sin \phi}{\phi - \sin \phi \cos \phi}, \qquad (4) $$

$$ G(\phi) = \frac{\phi^3}{3} \frac{4}{\sin \phi \cos \phi - \phi \cos 2\phi} \qquad (5) $$
are the corrections to the conventional lubrication theory satisfying $F(\phi \to 0) = G(\phi \to 0) = 1$. The surface tension gradient is

$$ \frac{d\sigma}{ds} \approx -\gamma \frac{dT^i}{ds}, \qquad (6) $$

where $T^i$ is the interfacial temperature and $\gamma = -d\sigma/dT$ is generally positive for pure fluids. The interface shape is defined with the equation

$$ \frac{\sigma}{\cos \phi} \frac{d^2h}{ds^2} - \Delta p = J^2 \left( \frac{1}{\rho_v} - \frac{1}{\rho_l} \right), \qquad (7) $$

that involves the vapor recoil; $\rho_v$ is the vapor density. The heat flux is determined as

$$ q = k(T_w - T^i)/\zeta, \qquad (8) $$

where $k$ is the liquid thermal conductivity. The mass flux is $J = q/\mathcal{L}$, where $\mathcal{L}$ is the latent heat. The wall temperature $T_w = T_{sat} + \Delta T$ differs from the saturation temperature $T_{sat}$ for the given system pressure by the superheating $\Delta T$. The vapor recoil also affects the interfacial temperature

$$ T^i = T_{sat} \left[ 1 + \frac{\Delta p}{\mathcal{L}\rho_l} + \frac{J^2}{2\mathcal{L}} \left( \frac{1}{\rho_v^2} - \frac{1}{\rho_l^2} \right) \right] + R^i q. \qquad (9) $$

The interfacial thermal resistance (unity accommodation coefficient is assumed) reads

$$ R^i = \frac{T_{sat} \sqrt{2\pi R_v T_{sat}} (\rho_l - \rho_v)}{2\mathcal{L}^2 \rho_l \rho_v}, \qquad (10) $$

where $R_v$ is the specific vapor constant. The boundary conditions for these equations are defined by the geometry at the contact line. At the CL ($s = 0$), the geometry implies

$$ \begin{aligned} h &= 0, \\ \phi &= \theta_{micro}, \end{aligned} \qquad (11) $$

where $\theta_{micro}$ is the static contact angle measured under adiabatic conditions. More precisely, it should be static receding angle for $U_{cl} > 0$, and static advancing angle in the opposite case. Another boundary condition of zero interface curvature

$$ \frac{\partial \phi}{\partial s} = 0 \qquad (12) $$

is given at the point M that is the right end of the integration interval $s = s_M$, which is, typically, commensurate with the DNS mesh size, 5–20 µm. The set of lubrication equations (3, 7) is of fourth order, so one more boundary condition is required. The temperature continuity along the interfaces at the CL $T^i = T_w$ is necessary so $J(s \to 0)$ is bounded. From (9), one obtains the condition

$$ \Delta p = \frac{\mathcal{L}\rho_l}{T_{sat}} (\Delta T - R^i J \mathcal{L}) - \frac{J^2 \rho_l}{2} \left( \frac{1}{\rho_v^2} - \frac{1}{\rho_l^2} \right), \qquad (13) $$

valid as $s \to 0$. The asymptotic analysis [10] of the above model yields

$$ J(s \to 0) = \frac{U_{cl} F(\theta_{micro})}{\frac{G(\theta_{micro})}{\theta_{micro} \rho_l} + \frac{l_s \mathcal{L} \theta_{micro}}{\mu k} \gamma}. \qquad (14) $$
are the corrections to the conventional lubrication theory satisfying $F(\phi \to 0) = G(\phi \to 0) = 1$. The surface tension gradient is

$$ \frac{d\sigma}{ds} \approx -\gamma \frac{dT^i}{ds}, \qquad (6) $$

where $T^i$ is the interfacial temperature and $\gamma = -d\sigma/dT$ is generally positive for pure fluids. The interface shape is defined with the equation

$$ \frac{\sigma}{\cos \phi} \frac{d^2h}{ds^2} - \Delta p = J^2 \left( \frac{1}{\rho_v} - \frac{1}{\rho_l} \right), \qquad (7) $$

that involves the vapor recoil; $\rho_v$ is the vapor density. The heat flux is determined as

$$ q = k(T_w - T^i)/\zeta, \qquad (8) $$

where $k$ is the liquid thermal conductivity. The mass flux is $J = q/\mathcal{L}$, where $\mathcal{L}$ is the latent heat. The wall temperature $T_w = T_{sat} + \Delta T$ differs from the saturation temperature $T_{sat}$ for the given system pressure by the superheating $\Delta T$. The vapor recoil also affects the interfacial temperature

$$ T^i = T_{sat} \left[ 1 + \frac{\Delta p}{\mathcal{L}\rho_l} + \frac{J^2}{2\mathcal{L}} \left( \frac{1}{\rho_v^2} - \frac{1}{\rho_l^2} \right) \right] + R^i q. \qquad (9) $$

The interfacial thermal resistance (unity accommodation coefficient is assumed) reads

$$ R^i = \frac{T_{sat} \sqrt{2\pi R_v T_{sat}} (\rho_l - \rho_v)}{2\mathcal{L}^2 \rho_l \rho_v}, \qquad (10) $$

where $R_v$ is the specific vapor constant. The boundary conditions for these equations are defined by the geometry at the contact line. At the CL ($s = 0$), the geometry implies

$$ \begin{aligned} h &= 0, \\ \phi &= \theta_{micro}, \end{aligned} \qquad (11) $$

where $\theta_{micro}$ is the static contact angle measured under adiabatic conditions. More precisely, it should be static receding angle for $U_{cl} > 0$, and static advancing angle in the opposite case. Another boundary condition of zero interface curvature

$$ \frac{\partial \phi}{\partial s} = 0 \qquad (12) $$

is given at the point M that is the right end of the integration interval $s = s_M$, which is, typically, commensurate with the DNS mesh size, 5–20 µm. The set of lubrication equations (3, 7) is of fourth order, so one more boundary condition is required. The temperature continuity along the interfaces at the CL $T^i = T_w$ is necessary so $J(s \to 0)$ is bounded. From (9), one obtains the condition

$$ \Delta p = \frac{\mathcal{L}\rho_l}{T_{sat}} (\Delta T - R^i J \mathcal{L}) - \frac{J^2 \rho_l}{2} \left( \frac{1}{\rho_v^2} - \frac{1}{\rho_l^2} \right), \qquad (13) $$

valid as $s \to 0$. The asymptotic analysis [10] of the above model yields

$$ J(s \to 0) = \frac{U_{cl} F(\theta_{micro})}{\frac{G(\theta_{micro})}{\theta_{micro} \rho_l} + \frac{l_s \mathcal{L} \theta_{micro}}{\mu k} \gamma}. \qquad (14) $$
The pressure value at the CL is easily obtained by substituting (14) into (13). One can see that the pressure is finite at the CL.

Within such an approach, the power per unit CL length supplied to the microregion is

$$Q_{\text{micro}} = \int_{0}^{s_{M}} qds. \qquad (15)$$

This value is one of four output parameters of the subgrid model. The second parameter is $\theta_{app} = \phi(s = s_M)$. The other two parameters are $x_M$ and $h_M$ (typically, of the order of the grid size), i.e. the Cartesian coordinates of the point M obtained as the solutions of Eqs. (1) for $s = s_M$. As for the input parameters, in addition to the fluid parameters, $l_s$ and $\theta_{micro}$ (which are typically all constant for a given simulation), they are two: the superheating $\Delta T$ and the CL velocity $U_{cl}$. The strategy of the integration of the subgrid model into the DNS is discussed in [13].

<figure><img src="image_3.png" /><figcaption>**Figure 2.** The subgrid problem results computed for water at atmospheric pressure; $s_M = 10$ µm, $l_s = 10$ nm, and $\theta_{micro} = 40^{\circ}$.</figcaption></figure>

The problem posed above is discretized by using the finite volume discretization scheme. The grid is chosen to be progressively denser when approaching CL. The system of equations is solved as linear, the nonlinear terms are treated with iterations. All the fluid parameters mentioned above are assumed to be constant. An example of $\theta_{app}$ variation resulting from such a problem is given in Fig. 2a.

### 3. Multiscale semi-analytical approach

The above calculation is quite delicate numerically because of several reasons. First, it is the necessity to treat a large range of the length scales at the same time (Fig. 3). Second, the convergence of the above algorithm is quite slow at small $\theta_{micro}$. A simpler method can be used to obtain the same output parameters. As mentioned in the introduction, there is a scale
The pressure value at the CL is easily obtained by substituting (14) into (13). One can see that the pressure is finite at the CL.

Within such an approach, the power per unit CL length supplied to the microregion is

$$Q_{\text{micro}} = \int_{0}^{s_{M}} q \, ds. \qquad (15)$$

This value is one of four output parameters of the subgrid model. The second parameter is $\theta_{app} = \phi(s = s_M)$. The other two parameters are $x_M$ and $h_M$ (typically, of the order of the grid size), i.e. the Cartesian coordinates of the point M obtained as the solutions of Eqs. (1) for $s = s_M$. As for the input parameters, in addition to the fluid parameters, $l_s$ and $\theta_{micro}$ (which are typically all constant for a given simulation), they are two: the superheating $\Delta T$ and the CL velocity $U_{cl}$. The strategy of the integration of the subgrid model into the DNS is discussed in [13].

<figure><img src="image_3.png" /><figcaption>**Figure 2.** The subgrid problem results computed for water at atmospheric pressure; $s_M = 10 \, \mu\text{m}$, $l_s = 10 \, \text{nm}$, and $\theta_{micro} = 40^{\circ}$.</figcaption></figure>

The problem posed above is discretized by using the finite volume discretization scheme. The grid is chosen to be progressively denser when approaching CL. The system of equations is solved as linear, the nonlinear terms are treated with iterations. All the fluid parameters mentioned above are assumed to be constant. An example of $\theta_{app}$ variation resulting from such a problem is given in Fig. 2a.

### 3. Multiscale semi-analytical approach

The above calculation is quite delicate numerically because of several reasons. First, it is the necessity to treat a large range of the length scales at the same time (Fig. 3). Second, the convergence of the above algorithm is quite slow at small $\theta_{micro}$. A simpler method can be used to obtain the same output parameters. As mentioned in the introduction, there is a scale
separation between the evaporation-induced liquid flow and the flow caused by the CL motion. While the evaporation defines the interface curvature in a region of a length $ℓ_V \sim 100$ nm (called the Voinov length) from the CL, the CL motion controls the interface shape at a larger (but still micrometric) distance. Note that the interface slope provided by the evaporation only is called the Voinov angle $\theta_V$. It can be obtained within a theory presented in sec. 2 corresponding to $U_{cl} = 0$ (dashed line in Fig. 2b). This quantity is similar to the ACA obtained with all the previous microregion models [2, 5]. The ACA accounting for the CL motion can be obtained [10–12] with the generalized Cox-Voinov equation

<figure><img src="image_4.png" /><figcaption>Figure 3. $q(x)$ variation calculated for the same parameters as Fig. 2b and for $\Delta T = 5$ K.</figcaption></figure>

$$ \theta_{app}^3 = \theta_V^3 - 9Ca \ln \frac{x_M}{\ell_V}, \qquad (16) $$

where $Ca = \mu U_{cl}/\sigma$ is the capillary number. This equation corresponds to the solution of the CL motion problem without evaporation satisfying the boundary condition (12). The Voinov length was shown to be controlled by the slip length:

$$ \ell_V \approx \frac{3l_s}{e\theta_V}, \qquad (17) $$

where $e = \exp(1)$. To obtain $\theta_{app}$, the multiscale approach consists in solving Eq. (3) without the last term in the l.h.s. thus getting $\theta_V$, and then use Eq. (16). The deviation of Eq. (16) from the rigorous solution is as small as $\sim 1.5^\circ$ (Fig. 2b). One mentions a stronger deviation of $\theta_{app}$ from $\theta_V$ (dashed line) that shows the necessity to account for the CL motion.

To propose a simplified approach for the calculation of $Q_{micro}$, one approximates the interface as a straight wedge with the opening angle $\theta_{app}$ so

$$ \begin{aligned} x &\approx s \cos \theta_{app}, \\ h &\approx s \sin \theta_{app}, \end{aligned} \qquad (18) $$

These approximations provide a high enough accuracy and, if used for $s = s_M$, define the output parameters $x_M$, $h_M$ of the subgrid model. One can now compare the heat flux variation resulting from sec. 2 and a “$R^i$ only” formula coming from Eqs. (8, 9) where all the other nanoscale effects are neglected:

$$ q = \frac{k\Delta T}{s\theta_{app} + R^i k}. \qquad (19) $$
separation between the evaporation-induced liquid flow and the flow caused by the CL motion. While the evaporation defines the interface curvature in a region of a length $ℓ_V \sim 100$ nm (called the Voinov length) from the CL, the CL motion controls the interface shape at a larger (but still micrometric) distance. Note that the interface slope provided by the evaporation only is called the Voinov angle $θ_V$. It can be obtained within a theory presented in sec. 2 corresponding to $U_{cl} = 0$ (dashed line in Fig. 2b). This quantity is similar to the ACA obtained with all the previous microregion models [2, 5]. The ACA accounting for the CL motion can be obtained [10–12] with the generalized Cox-Voinov equation

<figure><img src="image_4.png" /><figcaption>Figure 3. $q(x)$ variation calculated for the same parameters as Fig. 2b and for $\Delta T = 5$ K.</figcaption></figure>

$$ \theta_{app}^3 = \theta_V^3 - 9Ca \ln \frac{x_M}{\ell_V}, \qquad (16) $$

where $Ca = \mu U_{cl}/\sigma$ is the capillary number. This equation corresponds to the solution of the CL motion problem without evaporation satisfying the boundary condition (12). The Voinov length was shown to be controlled by the slip length:

$$ \ell_V \approx \frac{3l_s}{e\theta_V}, \qquad (17) $$

where $e = \exp(1)$. To obtain $θ_{app}$, the multiscale approach consists in solving Eq. (3) without the last term in the l.h.s. thus getting $θ_V$, and then use Eq. (16). The deviation of Eq. (16) from the rigorous solution is as small as $\sim 1.5°$ (Fig. 2b). One mentions a stronger deviation of $θ_{app}$ from $θ_V$ (dashed line) that shows the necessity to account for the CL motion.

To propose a simplified approach for the calculation of $Q_{micro}$, one approximates the interface as a straight wedge with the opening angle $θ_{app}$ so

$$ \begin{aligned} x &\approx s \cos \theta_{app}, \\ h &\approx s \sin \theta_{app}, \end{aligned} \qquad (18) $$

These approximations provide a high enough accuracy and, if used for $s = s_M$, define the output parameters $x_M$, $h_M$ of the subgrid model. One can now compare the heat flux variation resulting from sec. 2 and a “$R^i$ only” formula coming from Eqs. (8, 9) where all the other nanoscale effects are neglected:

$$ q = \frac{k\Delta T}{s\theta_{app} + R^i k}. \qquad (19) $$
The corresponding “R<sup>i</sup> only” curve is displayed in Fig. 3. The difference with the rigorous calculation of the sec. 2 model with Eq. (19) exists only for tiny s. They cause a negligible contribution to the integral (15), which explains an extremely weak difference between the numerical approach (solid red line in Fig. 2b) and the formula (dotted red line)

$$Q_{\text{micro}} = \frac{\Delta T \lambda_l}{\theta_{\text{app}}} \ln \left( \frac{s_M \theta_{\text{app}}}{R^i \lambda_l} + 1 \right) \qquad (20)$$

obtained by integration of the flux (19). This feature has been previously mentioned in [14]. It should be mentioned however that for a smaller $\theta_{\text{micro}}$, the microregion contribution becomes more important because the fluxes are larger and the influence of microscopic effects extends over larger distances. For the sake of comparison, we display in Fig. 3 also the $q(x)$ dependence without the account of nanoscale effects (Eq. (19) with $R^i = 0$, dashed curve) that results in the infinite $Q_{\text{micro}}$.

## 4. Conclusion

We propose a subgrid model for inclusion into the DNS describing the moving contact line problem with evaporation. Two approaches are compared: a rigorous calculation based on the generalized lubrication approximation and a semi-analytical model. The latter is much more economical as it requires a computational study with only one parameter, the local superheating at the contact line. The other parameter, the contact line velocity, is accounted for by using an analytical expression based on the multiscale argument. The semi-analytical approach provides a sufficient accuracy for the evaluation of the output parameters of the subgrid model, which are the apparent contact angle, heat power per unit contact line length and corresponding microregion sizes.

## References

[1] Sodtke C, Kern J, Schweizer N and Stephan P 2006 *Int. J. Heat Mass Transf.* **49** 1100 – 1106

[2] Torres L, Urbano A, Colin C and Tanguy S 2024 *J. Comput. Phys.* **497** 112602

[3] Park J H, Park I S and Jo D 2023 *Nucl. Eng. Technol.* **55** 1269 – 1279

[4] Wayner P C, Kao Y K and LaCroix L V 1976 *Int. J. Heat Mass Transf.* **19** 487 – 492

[5] Stephan P and Hammer J 1994 *Heat Mass Transf.* **30** 119 – 125

[6] Janeček V and Nikolayev V S 2012 *Europhys. Lett.* **100** 14003

[7] Janeček V and Nikolayev V S 2013 *Phys. Rev. E* **87** 012404

[8] Mathieu B 2003 *Etudes physiques, expérimentale et numérique des mécanismes de base intervenant dans les écoulements diphasiques* Ph.D. thesis Polytech Marseille, Université de Provence

[9] Janeček V and Nikolayev V S 2014 *Interfacial Phenom. Heat Transf.* **2** 363 – 383

[10] Zhang X and Nikolayev V S 2022 *J. Fluid Mech.* **948** A49

[11] Janeček V, Andreotti B, Pražák D, Bárta T and Nikolayev V S 2013 *Phys. Rev. E* **88** 060404(R)

[12] Nikolayev V S 2022 *The Surface Wettability Effect on Phase Change* ed Marengo M and De Coninck J (Springer) chap 6, pp 133 – 187

[13] Wei L, Bois G and Nikolayev V S 2024 *9th EuroTherm*

[14] Huber G, Tanguy S, Sagan M and Colin C 2017 *Int. J. Heat Mass Transf.* **113** 662 – 682
The corresponding “R<sup>i</sup> only” curve is displayed in Fig. 3. The difference with the rigorous calculation of the sec. 2 model with Eq. (19) exists only for tiny s. They cause a negligible contribution to the integral (15), which explains an extremely weak difference between the numerical approach (solid red line in Fig. 2b) and the formula (dotted red line)

$$Q_{\text{micro}} = \frac{\Delta T \lambda_l}{\theta_{\text{app}}} \ln \left( \frac{s_M \theta_{\text{app}}}{R^i \lambda_l} + 1 \right) \quad (20)$$

obtained by integration of the flux (19). This feature has been previously mentioned in [14]. It should be mentioned however that for a smaller $\theta_{\text{micro}}$, the microregion contribution becomes more important because the fluxes are larger and the influence of microscopic effects extends over larger distances. For the sake of comparison, we display in Fig. 3 also the $q(x)$ dependence without the account of nanoscale effects (Eq. (19) with $R^i = 0$, dashed curve) that results in the infinite $Q_{\text{micro}}$.

## 4. Conclusion

We propose a subgrid model for inclusion into the DNS describing the moving contact line problem with evaporation. Two approaches are compared: a rigorous calculation based on the generalized lubrication approximation and a semi-analytical model. The latter is much more economical as it requires a computational study with only one parameter, the local superheating at the contact line. The other parameter, the contact line velocity, is accounted for by using an analytical expression based on the multiscale argument. The semi-analytical approach provides a sufficient accuracy for the evaluation of the output parameters of the subgrid model, which are the apparent contact angle, heat power per unit contact line length and corresponding microregion sizes.

## References

[1] Sodtke C, Kern J, Schweizer N and Stephan P 2006 *Int. J. Heat Mass Transf.* **49** 1100 – 1106

[2] Torres L, Urbano A, Colin C and Tanguy S 2024 *J. Comput. Phys.* **497** 112602

[3] Park J H, Park I S and Jo D 2023 *Nucl. Eng. Technol.* **55** 1269 – 1279

[4] Wayner P C, Kao Y K and LaCroix L V 1976 *Int. J. Heat Mass Transf.* **19** 487 – 492

[5] Stephan P and Hammer J 1994 *Heat Mass Transf.* **30** 119 – 125

[6] Janeček V and Nikolayev V S 2012 *Europhys. Lett.* **100** 14003

[7] Janeček V and Nikolayev V S 2013 *Phys. Rev. E* **87** 012404

[8] Mathieu B 2003 *Etudes physiques, expérimentale et numérique des mécanismes de base intervenant dans les écoulements diphasiques* Ph.D. thesis Polytech Marseille, Université de Provence

[9] Janeček V and Nikolayev V S 2014 *Interfacial Phenom. Heat Transf.* **2** 363 – 383

[10] Zhang X and Nikolayev V S 2022 *J. Fluid Mech.* **948** A49

[11] Janeček V, Andreotti B, Pražák D, Bárta T and Nikolayev V S 2013 *Phys. Rev. E* **88** 060404(R)

[12] Nikolayev V S 2022 *The Surface Wettability Effect on Phase Change* ed Marengo M and De Coninck J (Springer) chap 6, pp 133 – 187

[13] Wei L, Bois G and Nikolayev V S 2024 *9th EuroTherm*

[14] Huber G, Tanguy S, Sagan M and Colin C 2017 *Int. J. Heat Mass Transf.* **113** 662 – 682