Three-Dimensional Long-Wave Instability of an Evaporation/Condensation Film

Abstract

This paper explores the stability and dynamics of a three-dimensional evaporating/condensing film while falling down a heated/cooled incline. Instead of using the Hertz–Knudsen–Langmuir relation, a more comprehensive phase-change boundary condition is employed. A nonlinear differential equation is derived based on the Benny-type equation, which takes into account gravity, energy transport, vapor recoil, effective pressure, and evaporation. The impact of effective pressure and vapor recoil on instability is studied using a linear stability analysis. The results show that spanwise perturbations can amplify the destabilizing effects of vapor recoil, leading to instability. Energy transport along the interface has almost no effect on the stability of the system, but it does influence the linear wave speed. Nonlinear evolution demonstrates that, in contrast to the vapor recoil effect, effective pressure can improve stability and delay film rupture. The self-similar solution demonstrates that the minimal film thickness decreases as ${(t_r-t)}^{1/2}$ and ${(t_r-t)}^{1/3}$ under the dominance of evaporation and vapor recoil, respectively.

1. Introduction

The evaporation/condensation of liquid films is ubiquitous in nature and technological systems with numerous applications, such as ice formation on aircraft wings, liquid films in mammalian lungs, heat exchanger design, and seawater desalination cooling systems [1,2,3,4]. As there exists a complex interaction of mechanisms, including vapor recoil, thermocapillarity, and van der Waals forces, several studies have revealed the importance of these factors’ effects on the stability of films [5,6,7,8]. A detailed review of long-scale models for evaporative thin films has been examined by Oron et al. [9], as well as by Craster and Matar [10] and Chattopadhyay et al. [11].

A number of theoretical studies targeted at describing these films and their potential to destabilize have been conducted to examine the issue of flat thin-film stability. The earliest stability studies focused on non-evaporating films falling down an inclined plane [12,13,14]. Yih [15,16] and Benjamin [17] demonstrated linear stability analyses of the falling thin-film model. Their findings suggested that this basic state is unstable in long-wavelength perturbations. Based on their works, Benney [18] derived an evolution equation for the film thickness, known as the Benny equation, and expanded the analysis of falling films into the nonlinear regime. After that, Williams and Davis [19] created a nonlinear stability model for which they carried out a numerical analysis for a finite perturbation. Comparing to the linear model, the results of their nonlinear model showed that the time to rupture was shortened and the needed perturbation wavelength was changed to one with a smaller value. These conclusions are supported by Sharma and Ruckenstein [20], who added Marangoni effects, and the stability of the film was found to be improved by the nonlinear effects associated with Marangoni motion. In recent studies, the understanding of interfacial phenomena at the microscopic level has been greatly enhanced by advanced experimental techniques. For instance, high-energy X-ray reflectivity measurements have provided significant insights into the hydrophobic gap at the water–octadecyl-trichlorosilane interface, showing the existence of a molecular-scale density deficit at the interface [21]. Similarly, investigations into the capillary wave dynamics of thin liquid polymer films have elucidated the behavior of liquid surfaces at nanometric thicknesses, contributing to our knowledge of surface tension and molecular interactions in confined systems [22].

In the context of liquid films, the Marangoni effect can play an important role in determining stability and flow behavior. Davis [6] examined Marangoni effects on thin-film stability, which mostly focused on thermocapillary effects. This analysis showed that hydrodynamic instabilities, which appear as 2D and 3D surface waves, are powered by heat gradients. The investigation of evaporating and condensing thin liquid films on a hot inclined plane by Bankhoff [23] demonstrated the impact of phase shift on thin-film stability. The findings of the study found that evaporation from the surface has a destabilizing effect, while condensation has the opposite influence, which was recently proved by [24]. Thus, thermocapillarity potentially becomes one of the most researched phenomena in fluid mechanics, which leads to variations in surface tension and can amplify the runaway effect as warmer liquid from the depression is dragged towards the cooler areas of the interface.

Sreenivasan and Lin [25] used the linear theory to study the onset of surface tension-driven stationary circulation in a thin-film flow down a heated incline. Their results were supported by Kelly et al. [26], who investigated the impact of a temperature-dependent surface tension with long-wavelength disturbances. While many studies have adopted a uniform heated plane [27,28,29], Scheid et al. [30] were more concerned with the instabilities of thin-film flow over a localized heater to determine the influence of the temperature from the inclined plane. Similarly, Gambaryan-Roisman [31] investigated how the non-uniformity of substrate thermal conductivity affected the hydrodynamics and heat transport in thin liquid films while taking surface tension, thermocapillarity, and evaporation into account. Recently, a two-dimensional flow of a thin film over a uniformly heated/cooled slippery inclined substrate was analyzed by Ding et al. [32] to capture the Marangoni effect. They employed a weighted-residual model and a Benney-type model for comparison, and they discovered that the weighted-residual model was more logical than the Benney-type model. To express a thin liquid film flow on a non-uniformly heated, slippery inclined plane, Ref. [33] developed a mathematical model within the framework of the long-wave approximation method.

As volatile processes significantly impact the thermal response of the system in actual circumstances, the Hertz–Knudsen–Langmuir relation was adopted by some researchers to describe the phase-change occurrence at the liquid–vapor interface [5,7,34,35,36,37,38,39]. This relationship is predicated on the idea that mass momentum and energy transfer are unrelated to evaporation and only affect the states of the liquid and vapor phases [40]. Burelbach et al. [5] provided a thorough explanation of how the pertinent equations were derived. As a result of the assumption that the ratios between the densities, viscosities, and thermal conductivities of gases and liquids are very tiny, there is significant simplification and the creation of so-called one-sided models [5]. Instead of applying the traditional Hertz–Knudsen–Langmuir relation, Shklyaev and Fried [41] revisited this problem and imposed the more general equations that Fried et al. [42,43] proposed for expressing the balance of configurational momentum at the liquid–vapor interface. They developed a nonlinear evolution equation determining film thickness, and two new terms were added to their model. One was known as the effective pressure accounting for the vapor recoil, $p-E^2{D}^{-1}{J}^2$, and the other was the energy transmission within the liquid–vapor interface, $N(-{\nabla }_s\left((\mathit{u}·\mathit{t})\mathit{t}\right))$. The two terms were stabilized according to the results of a linear stability analysis.

Therefore, in order to discuss the behavior of an evaporation or condensation film under general boundary conditions in three dimensions, we draw on the research conducted by Shklyaev and Fried [41]. We aim to specifically look into how the stability of a system is impacted by spanwise perturbations, effective pressure, evaporation, and vapor recoil. The structure of this paper as follows. Utilizing the long-wave assumption, we develop the mathematical model in Section 2. The linear stability theory is employed to reveal the influences of effective pressure and energy transport along the interface on the stability of the system in Section 3. In Section 4, we investigate the nonlinear evolution of the evaporation/condensation film. A conclusion is made in Section 5.

2. Mathematical Formulation

We examine the behavior of a three-dimensional incompressible fluid. It falls down along an inclined plane which is either heated or cooled, as depicted in Figure 1. The angle of inclination is denoted by $\vartheta$. The physical properties of the fluid are characterized by its density, kinematic viscosity, and thermal conductivity, denoted by $(\rho ,\nu ,\kappa )$ for the liquid phase and $({\rho }_v,{\nu }_v,{\kappa }_v)$ for the vapor phase. We assume that the physical properties of the vapor phase are significantly smaller than those of the liquid phase, which can be expressed mathematically as

$$\frac{{\rho }_v}{\rho }\ll 1,\frac{{\nu }_v}{\nu }\ll 1,\frac{{\kappa }_v}{\kappa }\ll 1.$$

(1)

The dynamics of the liquid film are governed by the continuity equation, the momentum equation, and the energy equation, which can be written as

$$\nabla ·\mathit{u}=0,$$

(2)

$$\rho [\frac{\partial \mathit{u}}{\partial t}+(\mathit{u}·\nabla )\mathit{u}]=-\nabla p+\mu {\nabla }^2\mathit{u}+\rho \mathit{g},$$

(3)

$$\frac{\partial T}{\partial t}+(\mathit{u}·\nabla )T=k_{th}{\nabla }^{2}T,$$

(4)

where $\mathit{u}=(u,v,w)$ is the velocity field and $\nabla =({\partial }_x,{\partial }_y,{\partial }_z)$ is the Hamilton operator. $\rho$ is the fluid density, p is the pressure, $\mu$ is the dynamic viscosity of fluid, $\mathit{g}=(g\cos \vartheta ,0,g\sin \vartheta )$ is the gravitational acceleration, and t is time. The temperature and the thermal diffusivity are denoted by T and $k_{th}$, respectively. These equations form the basis of numerical simulations and enable us to study the behavior of the liquid film under various conditions.

Here, we consider no slip, no penetration, and a constant temperature condition at the solid wall $z=0$, which can be written as follows:

$$\mathit{u}=0,T=T_0.$$

(5)

At the liquid–vapor phase interface $z=h(x,y,t)$, several conditions need to be satisfied, including mass flux conservation, stress balance, energy balance, and kinetic conditions. The mass flux conservation can be expressed as

$$J=\rho (\mathit{u}-{\mathit{u}}_i)·\mathit{n}={\rho }_v({\mathit{u}}_v-{\mathit{u}}_i)·\mathit{n},$$

(6)

where ${\mathit{u}}_i$ represents the velocity at the interface, ${\mathit{u}}_v$ is the velocity of the vapor, and $\mathit{n}$ is the unit normal vector to the interface.

For the purposes of this study, we assume that the surface tension $\sigma$ is linearly dependent on the temperature T. Specifically, we use the following expression to model the temperature dependence of the surface tension:

$$\sigma ={\psi }_s-{\eta }_s(T-T_s),$$

(7)

where ${\psi }_s$ and ${\eta }_s$ are constants of the interfacial free-energy and entropy densities arising from $T=T_s$; $T_s$ represents the saturation temperature where the liquid and vapor phases are in thermodynamic equilibrium.

The stress balance of the interface conditions is

$$\Sigma ·\mathit{n}+\frac{J^2}{{\rho }_v}·\mathit{n}={\nabla }_s\sigma -\sigma ({\nabla }_s·\mathit{n})\mathit{n},$$

(8)

where $\Sigma =-p\mathit{I}+\mu [\nabla \mathit{u}+{(\nabla \mathit{u})}^T]$ is the stress tensor, $\mathit{I}$ is the identity tensor, ${\nabla }_s=\nabla -\mathit{n}(\mathit{n}·\nabla )$ is the surface gradient operator, and $J^2/{\rho }_v$ is the vapor recoil term, which describes the contribution of the momentum of the evaporating or condensing vapor to the overall momentum balance of the liquid–vapor interface. When the liquid is evaporating, it releases molecules of vapor into the surrounding gas. These vapor molecules have a certain momentum due to their thermal motion, and this momentum can cause a recoil effect on the liquid–vapor interface. Similarly, when vapor condenses back into the liquid, the momentum of the vapor molecules can cause a recoil effect on the interface in the opposite direction.

The energy balance at the liquid–vapor interface is expressed as

$$T_s{\eta }_s(K{\mathit{u}}_i·\mathit{n}-{\nabla }_s·\left(({\mathit{u}}_i·\mathit{t})\mathit{t}\right))=\kappa \nabla T·\mathit{n}+\rho L(\mathit{u}-{\mathit{u}}_i)·\mathit{n},$$

(9)

where $K=-\nabla ·\mathit{n}$ is the curvature of the interface, L is the latent heat of phase change, and $\mathit{t}$ is the unit tangent vector on the interface. The left-hand side of the equation represents the energy flux at the interface due to interfacial shear and convection. Note that the term $K\mathit{u}·\mathit{n}$ is much smaller than the term ${\nabla }_s·\left(({\mathit{u}}_i·\mathit{t})\mathit{t}\right)$ in the assumption of the long-wave model; thus, this term can be neglected in this paper [41]. On the right-hand side, the term $\kappa \nabla T·\mathit{n}$ represents the heat flux as conduction across the interface, while the term $\rho L(\mathit{u}-{\mathit{u}}_i)·\mathit{n}$ denotes the energy released or absorbed since phase change at the interface.

In order to describe the phase-change phenomenon, we utilize the equation proposed by Fried et al. [43]:

$$\frac{1}{\rho }{\beta }_{s}J=\rho L(T/T_s-1)+(p-J^2/{\rho }_v)-\frac{1}{2}\rho {\left|u\right|}^2,$$

(10)

where ${\beta }_s$ is a modulus associated with the kinetics of attachment and detachment at the interface [41,43]. The term $(p-J^2/{\rho }_v)-\frac{1}{2}\rho {\left|u\right|}^2$ on the right-hand side of Equation (10) is defined as the effective pressure term. The kinetic energy term $\frac{1}{2}\rho {\left|u\right|}^2$ is much smaller than the other two terms, and is therefore negligible [41]. Wayner took a different approach and discussed the influence of pressure on the evaporating meniscus, where he did not consider the vapor recoil’s influence on the mass flux [44]. However, in the present work we consider the vapor recoil effect and use the model as described in Equation (10).

2.1. Nondimensionalization

We start by introducing the following scales: the initial mean thickness of the film, $h_0$, $\frac{\nu }{h_0}$, $\frac{h_0^2}{\nu }$, $\frac{\rho {\nu }^2}{h_0^2}$, $T_0-T_s$, $\frac{\kappa (T_0-T_s)}{L{h}_0}$ as the length, velocity, time, pressure, temperature, and mass flux scales, respectively. Using these scales, we can nondimensionalize the system of Equations (2)–(10). The resulting nondimensional governing equations are

$$\nabla ·\mathit{u}=0,$$

(11)

$$\frac{\partial \mathit{u}}{\partial t}+(\mathit{u}·\nabla )\mathit{u}=-\nabla p+{\nabla }^2\mathit{u}+G{\mathit{e}}_3,$$

(12)

$$Pr[\frac{\partial \theta }{\partial t}+(\mathit{u}·\nabla )\theta ]={\nabla }^2\theta ,$$

(13)

where $G=\frac{g{h}_0^3}{{\nu }^2}$, ${\mathit{e}}_3$ represents the direction of gravity, the Prandtl number is denoted as $Pr=\frac{\kappa }{\nu }$, $G\sin \vartheta$ is defined as the Reynolds number, and $\theta$ is the dimensionless temperature.

The dimensionless boundary condition at solid wall $z=0$ is

$$u=v=w=0,\theta =1.$$

(14)

At the liquid–vapor interface, the boundary conditions can be written as

$$\begin{array}{c}\hfill (\mathit{u}-{\mathit{u}}_i)·\mathit{n}=D({\mathit{u}}_v-{\mathit{u}}_i)·\mathit{n}=EJ,\end{array}$$

(15)

$$\begin{array}{c}\hfill N(-{\nabla }_s\left((\mathit{u}·\mathit{t})\mathit{t}\right))-(\nabla \theta ·\mathit{n}+J)=0,\end{array}$$

(16)

$$\begin{array}{c}\hfill \Sigma ·\mathit{n}·\mathit{n}+E^2{D}^{-1}{J}^2=KS(1-C\theta ),\end{array}$$

(17)

$$\begin{array}{c}\hfill \Sigma ·\mathit{n}·\mathit{t}+Ma·P{r}^{-1}(\mathit{t}·\nabla )\theta =0.\end{array}$$

(18)

where $D=\frac{{\rho }_v}{\rho }$ is the density ratio; $E=\frac{k_{th}▵T}{\rho \nu L}$ is the evaporation number. $S=\frac{{\psi }_s{h}_0}{\rho {\nu }^2}$ is the dimensionless surface tension. $C=\frac{{\eta }_s▵T}{{\psi }_s}$ is the Crispation number, $N=\frac{T_s\nu {\eta }_s}{k_{th}▵T{h}_0}$, which accounts for the thermal diffusion along the interface, and $Ma=\frac{{\eta }_s▵T{h}_0}{\rho \nu \kappa }$ is the Marangoni number.

Thus, recall that the dimensionless phase-change condition at the interface becomes

$$A_{1}J=\theta +A_2(p-E^2{D}^{-1}{J}^2),$$

(19)

where $A_1=\frac{k_{th}{\beta }_s{T}_s}{{\rho }^2{L}^2{h}_0}$; $A_2=\frac{T_s{\nu }^2}{L▵T{h}_0^2}$. $A_1$ characterizes how far the system is from thermodynamic equilibrium, and $A_2$ expresses the magnitude of the effective pressure.

The dimensionless form of the kinematic condition of the moving interface is

$$EJ=\frac{w-h_t-u{h}_x-v{h}_y}{{(1+h_x+h_y)}^{1/2}}.$$

(20)

where $h_t$, $h_x$, and $h_y$ are the partial derivations of liquid film thickness h.

In this paper, $A_2$, N, and $E^2{D}^{-1}$ are the three important dimensionless parameters that we will focus on to study the influences of the effective pressure, the interfacial energy transport, and the vapor recoil effect on the dynamics of the film.

2.2. Long-Wave Model

The long-wave model is a common simplification used to study the dynamics of thin liquid films in which the characteristic length scales of the system are much larger than the mean film thickness. In this model, the governing equations are derived by assuming that the perturbations in the film thickness, velocity, and temperature vary slowly in the streamwise and spanwise directions.

To derive the long-wave model for our system, we introduce the small parameter $ϵ=h_0/l$, where l is the characteristic length scale of the system. In this paper, we fix $ϵ=0.2$, which means that the characteristic length l is much larger than the mean film thickness $h_0$, and this is a typical assumption for thin liquid films. Based on the long-wave assumption, we make the following transformation:

$$({\partial }_t,{\partial }_x,{\partial }_y)=ϵ({\partial }_{\tilde{t}},{\partial }_{\tilde{x}},{\partial }_{\tilde{y}}){\partial }_{xx}={ϵ}^2{\partial }_{\tilde{x}\tilde{x}},$$

(21)

For ease of writing, we will omit the wave number symbol in the following statements.

We use the perturbation method to study each physical quantity of the system, assuming that the solution of the parameter follows a power law with respect to $ϵ$. We expand the variables in powers of $ϵ$ as shown below:

$$\begin{array}{c}\hfill u=u_0+ϵu_1+\dots ,\end{array}$$

(22)

$$\begin{array}{c}\hfill v=v_0+ϵv_1+\dots ,\end{array}$$

(23)

$$\begin{array}{c}\hfill w=ϵw_0+{ϵ}^2{w}_1+\dots ,\end{array}$$

(24)

$$\begin{array}{c}\hfill p=p_0+ϵp_1+\dots ,\end{array}$$

(25)

$$\begin{array}{c}\hfill \theta ={\theta }_0+ϵ{\theta }_1+\dots \end{array}$$

(26)

$$\begin{array}{c}\hfill J=J_0+ϵJ_1+\dots ,\end{array}$$

(27)

Substituting Equations (22)–(27) into governing Equations (11)–(13) together with using boundary conditions (14)–(20), we can derive a set of long-wave equations that describe the dynamics of the liquid film with the leading order and first order in term $ϵ$. These equations involve the dimensionless parameters $A_1$, $A_2$, N, E, and S, which were defined in the previous section.

The leading orders of the system are

$$\begin{array}{c}\hfill u_{0,x}+v_{0,y}+w_{0,z}=0,\end{array}$$

(28)

$$G\sin \vartheta +u_{zz}=0,$$

(29)

$$v_{0,zz}=0,$$

(30)

$$p_{0,z}+G\cos \vartheta =0,$$

(31)

$${\theta }_{0,zz}=0.$$

(32)

With the velocity and temperature boundary conditions at $z=0$,

$$u_0=v_0=w_0=0,{\theta }_0=1,$$

(33)

and at $z=h(x,y,t)$,

$$\begin{array}{c}\hfill p_0=E^2{D}^{-1}{J}_0^2-\overline{S}(h_{xx}+h_{yy}),\end{array}$$

(34)

$$u_{0,z}=v_{0,z}=0,$$

(35)

$$\overline{N}(u_{0,x}+v_{0,y})+J_0+{\theta }_{0,z}=0,$$

(36)

$${\theta }_0=A_1{J}_0-A_2(p_0-E^2{D}^{-1}{J}_0^2),$$

(37)

$$\overline{E}{J}_0=-h_t-u_0{h}_x-v_0{h}_y+w_0.$$

(38)

It is easy to obtain the solution of the following parameters:

$$u_0=G\sin \vartheta (hz-\frac{z^2}{2}),v_0=0,w_0=\frac{-G\sin \vartheta h_x{z}^2}{2},$$

(39)

$$p_0=G\cos \vartheta (h-z)+E^2{D}^{-1}{J}_0^2-\overline{S}(h_{xx}+h_{yy}),$$

(40)

$${\theta }_0=1+\left[A_1{J}_0+A_2\overline{S}(h_{xx}+h_{yy})-1\right]\frac{z}{h},$$

(41)

$$J_0=\frac{-G\overline{N}\sin \vartheta h^2{h}_x+1-A_2\overline{S}(h_{xx}+h_{yy})}{A_1+h}.$$

(42)

The first orders of the system are

$$\begin{array}{c}\hfill u_{1,x}+v_{1,y}+w_{1,z}=0,\end{array}$$

(43)

$$u_{0,t}+u_0{u}_{0,x}+w_0{u}_{0,z}=-p_{0,x}+u_{1,zz}=0;$$

(44)

$$v_{1,zz}=p_{0,y},$$

(45)

$$w_{0,t}+p_{1,z}=w_{0,zz},$$

(46)

$$Pr({\theta }_{0,t}+u_0{\theta }_{0,x}+w_0{\theta }_{0,z})={\theta }_{1,zz}.$$

(47)

And the first-order boundary conditions at $z=0$ are

$$u_1=v_1=w_1=0,{\theta }_1=0,$$

(48)

At the liquid–vapor interface, the first-order boundary conditions are

$$-p_1+2(-u_{0,z}{h}_x-v_{0,z}{h}_y+w_{0,z})=0,$$

(49)

$$u_{1,z}+MaP{r}^{-1}({\theta }_{0,x}+h_x{\theta }_{1,z})=0,$$

(50)

$${\theta }_1=A_1{J}_1-A_2{p}_1$$

(51)

$$\overline{E}{J}_1=-h_{1,t}-u_0{h}_{1,x}+w_1$$

(52)

Integrating twice with respect to $u_{1,zz}$ in Equation (44), while incorporating boundary conditions Equations (50) and (48), yields

$$\begin{array}{cc}\hfill u_1=& \left[G\cos \vartheta h_x+E^2{D}^{-1}{\left({J_0}^2\right)}_x-\overline{S}(h_{xxx}+h_{xyy})\right](\frac{z^2}{2}-hz)\hfill \\ \hfill & +\frac{1}{6}{G}^2{\sin }^2\vartheta h{h}_x(\frac{z^4}{4}-h^{3}z)-\frac{G\sin \vartheta }{2}(\frac{z^3}{3}-h^2\mathrm{z})(\overline{E}{J}_0+G\sin \vartheta h^2{h}_x)\hfill \\ \hfill & -MaP{r}^{-1}{\left[\frac{A_1-A_1\overline{N}G\sin \vartheta h^2{h}_x+A_2\overline{S}h(h_{xx}+h_{yy})}{A_1+h}\right]}_{x}z,\hfill \end{array}$$

(53)

We can also obtain that the velocity component in the y direction is

$$\begin{array}{ll}{v}_1& =\left[G\cos \vartheta h_y+E^2{D}^{-1}{\left({J_0}^2\right)}_y-\overline{S}(h_{xxy}+h_{yyy})\right](\frac{z^2}{2}-hz)\\ & +Ma{\mathrm{Pr}}^{-1}{\left[\frac{h+A_1\overline{N}G\sin \vartheta h^2{h}_x-A_2\overline{S}h(h_{xx}+h_{yy})}{A_1+h}\right]}_{y}z,\end{array}$$

(54)

Note that to investigate the effects of these parameters we rescale the parameters as $(\overline{E},\overline{N},\overline{S})=({ϵ}^{-1}E,ϵN,{ϵ}^{2}S)$, which is shown in the above equations. The term ${ϵ}^{2}C$ was neglected due to its much smaller magnitude compared to the unity.

With the rescaling parameters, the integral form of kinematic condition Equation (20) can be written as

$$h_t+{\partial }_x{\int }_0^{h}udz+{\partial }_y{\int }_0^{h}vdz+\overline{E}J=0,$$

(55)

Substituting the velocity into the kinematic equation, we obtain the evolution equation of the film thickness $h(x,y,t)$:

$$h_t+Q+ϵ(R+\overline{E}{J}_1)+O\left({ϵ}^2\right)=0.$$

(56)

where

$$Q=G{h}^2{h}_x\sin \vartheta +\overline{E}{J}_0,$$

(57)

$$\begin{array}{cc}\hfill R={\left\{\frac{-h^3{\varphi }_1}{3}-\frac{3}{40}{G}^2{\sin }^2\vartheta h^6{h}_x+\frac{5G\sin \vartheta }{24}{h}^{4}Q+\frac{h^2}{2}MaP{r}^{-1}{\left(\frac{h+Fh}{A_1+h}\right)}_x\right\}}_x& \\ \hfill +{\left\{\frac{-h^3{\varphi }_2}{3}+\frac{Ma{\mathrm{Pr}}^{-1}}{2}{h}^2{\left(\frac{h+Fh}{A_1+h}\right)}_y\right\}}_y& ,\hfill \end{array}$$

(58)

and $J_1=H_1+H_2+H_3+H_4$, with

$$\begin{array}{cc}\hfill & H_1=\frac{\mathrm{Pr}{h}^3}{3{(A_1+h)}^3}\left\{(A_1+h)\left[-A_1\overline{N}G\sin \vartheta {\left(hQ\right)}_x+A_2\overline{S}(Q_{xx}+Q_{yy})\right]+(1+F)Q\right\},\hfill \end{array}$$

(59)

$$\begin{array}{cc}\hfill & H_2=\frac{Pr}{A_1+h}\left\{-\frac{3h^5}{20}G\sin \vartheta {\left(\frac{F+1}{A_1+h}\right)}_x+\frac{h^4{h}_{x}G\sin \vartheta }{8}\frac{F+1}{A_1+h}\right\}-\frac{2A_{2}G\sin \vartheta h{h}_x}{A_1+h},\hfill \end{array}$$

(60)

$$\begin{array}{cc}\hfill & H_3=\frac{\overline{N}h}{A_1+h}{\left\{\frac{{\varphi }_1{h}^2}{2}+\frac{1}{8}{G}^2{\sin }^2\vartheta h^5{h}_x-\frac{h^{3}G\sin \vartheta }{3}Q-MaP{r}^{-1}{\left(\frac{h+Fh}{A_1+h}\right)}_{x}h\right\}}_x,\hfill \end{array}$$

(61)

$$\begin{array}{cc}\hfill & H_4=\frac{\overline{N}h}{A_1+h}{\left\{\frac{{\varphi }_2{h}^2}{2}-MaP{r}^{-1}{\left(\frac{h+Fh}{A_1+h}\right)}_{y}h\right\}}_y,\hfill \end{array}$$

(62)

where

$$J_0=\frac{-G\overline{N}\sin \vartheta h^2{h}_x+1-A_2\overline{S}(h_{xx}+h_{yy})}{A_1+h},$$

(63)

$${\varphi }_1=G\cos \vartheta h_x+E^2{D}^{-1}{\left({J_0}^2\right)}_x-\overline{S}(h_{xxx}+h_{xyy}),$$

(64)

$${\varphi }_2=G\cos \vartheta h_y+E^2{D}^{-1}{\left({J_0}^2\right)}_y-\overline{S}(h_{xxy}+h_{yyy}),$$

(65)

$$F=A_1\overline{N}G\sin \vartheta h{h}_x-A_2\overline{S}(h_{xx}+h_{yy}).$$

(66)

The equation presented in Equation (56) is the same as the one given by Shklyaev and Fried [41] when ${\partial }_y=0$ and $\beta =0$. Additionally, it is identical to the equation presented by Burelbach et al. [5] when $A_2=N=0$.

The values of $(A_1,A_2,\overline{N},\overline{E},\overline{S})$ were originally proposed by Shklyaev and Fried and Wei et al. [41,45] for water, sodium and 1,2-Ethanediol. In this paper, we assume that the film thickness is macroscopic at about $h_0$ ∼ $0.01$ mm and that the system is primarily governed by hydrodynamic effects and heat transfer rather than intermolecular interactions. Thus, the Van der Waals force can be neglected, together with $\rho \approx {10}^3$ kg m⁻³, ${\eta }_s\approx$10⁻⁴ N m⁻¹K⁻¹, ${\beta }_s\approx$10⁶ kg m⁻²s⁻¹, $\Delta T\approx$(10⁻³∼1) K, ${\mathsf{\Psi }}_s\approx$(0.01∼0.1) Nm⁻¹, $T_s\approx$(10²∼10³) K, $k_{th}\approx$(1∼10) Wm⁻¹ K⁻¹, $L\approx$(10⁵∼10⁶) J kg⁻¹, and $\nu \approx ({10}^{-7}$∼10⁻⁶) m²s⁻¹, where the evaporation and vapor recoil effects are the dominant factors. The value ranges of dimensionless parameters in this paper are as follows: $A_1$ ranges from $\left[0\right]$ to $\left[0.1\right]$, $A_2$ is in the range $[0,0.1]$, $\overline{S}$ is in the range $[0.1,2]$, $\overline{N}$ is in the range $[0,0.06]$, $Ma$ is in the range $[0,0.5]$, the Prandtl number is fixed at $Pr=1$, and the inclined angle is fixed at $\vartheta =\pi /4$. The evaporation number E is in the range of $[-0.15,0.15]$. When E is positive, the system is evaporating. When $E=0$, it is a nonvolatile film; otherwise, it is condensing. We discovered that the thermocapillary effect destabilizes the system.

3. Linear Stability Analysis

First, we conduct a linear stability analysis of the system. We consider a small perturbation of the base state of the system by using normal mode analysis. Specifically, the perturbation is expressed as $h=h_0\left(\tau \right)+\widehat{h}\exp [i(\alpha x+\beta y)+\lambda t]$, where $\alpha$ and $\beta$ represent the streamwise and spanwise wavenumbers, respectively. $\lambda ={\lambda }_r+i{\lambda }_i$ is the complex growth rate; the real part represents the effective growth rate, which characterizes the rate where the perturbation either grows or decays in time; and the imaginary part is related to the wavespeed.

The base state of the system was investigated through the resolution of Equation (56), which is expressed as follows:

$$h=-A_1+\sqrt{{(A_1+1)}^2-2\overline{E}\tau }+O\left(ϵ\right),$$

(67)

It should be noted that $\tau$ is a parameter that is used to define the base state of the system, which can change with time. The value of $\tau$ is usually chosen to be a specific value, while the value of t changes continuously as times progresses, and it is used to track the evolution of perturbation over time. Furthermore, Equation (67) demonstrates that an evaporating film will desiccate within a definite time frame of $\tau =\frac{2A_1+1}{2\overline{E}}$.

For the leading-order problem, the dispersion relationship was obtained:

$$\lambda =-i\alpha G\sin \vartheta (h_0^2-\frac{\overline{E}\overline{N}}{A_1+h_0})+\frac{\overline{E}}{{(A_1+h_0)}^2}-\frac{A_2\overline{E}\overline{S}}{A_1+h_0}({\alpha }^2+{\beta }^2)+O\left(ϵ\right).$$

(68)

We also tested the first-order dispersion relationship, and the growth rate is expressed in by $ϵ$; the detailed expression is shown in Appendix A.

The contour plot in Figure 2 illustrates the variation in the first-order real growth rate with respect to the base-state time $\tau$. It can be observed that the effective growth rate exhibits a saddle shape. As time progresses, the system experiences a sharp increase in growth rate, causing it to become highly unstable. This is attributed to the significant reduction in the liquid film thickness resulting from film evaporation within a limited time. The stability study following film rupture will not be further detailed since it has no practical relevance. In the subsequent analysis, we examine the scenario at the base-state time $\tau =0$, i.e., $h_0=1$.

We first start the analysis by ignoring the spanwise influence, i.e., $\beta =0$. The impact of parameter $A_1$, which represents the distance of the system from thermodynamic equilibrium, is demonstrated in Figure 3a. It can be inferred that parameter $A_1$ plays a crucial role in determining the stability of the liquid film, where the real growth rate decreases as $A_1$ increases. Another important stability parameter, effective pressure $A_2$, is depicted in Figure 3b. Both the growth rate and the critical wave number decrease with the increasing effective pressure, which is consistent with the analytical findings in Equation (68). However, at small wave numbers (i.e., $\alpha \ll 1$), stability is barely affected by effective pressure $A_2$, since $A_2$ is proportional to the wave number in Equation (68). Extensive research has been conducted in previous studies regarding the Marangoni effect [32]. As shown in Figure 3c, it can be observed that within the given range of values, the destabilizing impact of the Marangoni effect is minimal. Thus, we will proceed with setting $Ma=0$ for the forthcoming studies.

The influences of spanwise perturbation are illustrated in Figure 4a. We observe that spanwise perturbations always play a stabilizing role in the absence of gravity and evaporation (solid lines), as the maximum effective growth rate consistently decreases with an increase in $\beta$. In fact, spanwise perturbations are typically considered stabilizing in linear stability analysis of a falling liquid film. This is because spanwise perturbations introduce cross-stream vorticity, which interacts with the streamwise vorticity generated by the primary instability to produce three-dimensional flow structures. These three-dimensional structures can enhance the stabilizing effect of viscosity and delay the onset of the primary instability. However, in the presence of film evaporation, as demonstrated by the dashed lines in Figure 4a, this stabilizing effect of spanwise perturbations does not hold true. Since the evaporation effect causes the film to become thinner, making it more susceptible to destabilizing mechanisms such as vapor recoil, the ${\phi }_2$ term in Equation (A7) can explain this phenomenon. Therefore, the stabilizing effect of spanwise perturbations is diminished in the presence of evaporation.

In the present work, surface energy transport has a negligible effect on stability, as shown in Figure 4b, as the parameter $\overline{N}$ exists in the first-order part of the effective growth rate ${\lambda }_r$. However, the influence of energy transport on wave speed cannot be discounted, and Figure 5 provides insight into this effect. For an evaporating film $(E>0)$, an increment in energy transport leads to an increase in wave speed, as shown in Figure 5a. In the case of a condensing film, $\overline{N}$ rises when the wave speed increases at a small wave number, while it decreases at a large wave number, with the critical wave number being $\alpha \approx 0.9$. These results suggest that surface energy transport may not have a significant impact on the stability of the film, but it can affect the wave speed and, hence, the dynamics of the system. Thus, it is important to consider these factors when studying the behavior of evaporating or condensing films.

Figure 6 illustrates the destabilizing effects of vapor recoil and gravity on the system. Note that when $G<0$, it represents an inverted liquid film. As observed in Figure 6a, increasing the vapor recoil effect can result in a higher maximum growth rate, whereas this instability can be mitigated by increasing the effective pressure coefficient $A_2$. Moreover, the maximum real growth rate is relatively small under microgravity conditions. When the gravity increases, the maximum effective growth rate also increases. Under the same conditions, the inverted liquid film appears to be more unstable, as demonstrated by Figure 6b. It is noteworthy that a larger vapor recoil term leads to a larger maximum growth rate, indicating its destabilizing effect on the system. These findings show that both vapor recoil and gravity are critical factors affecting the stability of the system.

4. Nonlinear Evolution

This section aims to provide insights into the nonlinear dynamical evolution of a uniform interface subjected to small finite disturbances. To perturb the interface, we introduce a harmonic wave:

$$h(x,y,0)=1+0.1\cos \left(\alpha x\right)\cos \left(\beta y\right).$$

(69)

Equation (56) subjected to periodic boundary conditions was solved numerically by using the Fourier spectral method. We assume that the solution of $h(x,y,t)$ can be expanded by the Fourier series:

$$h(x,y,t)=\sum _{n,m}{\widehat{h}}_{n,m}\exp (i{k}_{n}x+i{k}_{m}y),$$

(70)

where ${\widehat{h}}_{n,m}\left(t\right)$ are the time-dependent Fourier coefficients. $k_n=\frac{2n\pi }{L_x}$ and $k_m=\frac{2m\pi }{L_y}$ are the wave numbers in the x and y directions, respectively, where $L_x$ and $L_y$ are the lengths of the computational domain in the x and y directions. We chose the number of Fourier modes in the x direction to be $N_x=64$ and those in the y direction to be $N_y=32$. Thus, there exist a total of $N=2048$ Fourier modes, which is sufficient to ensure numerical accuracy.

To investigate the influence of evaporation and the vapor recoil effect, we first consider the quasi-equilibrium case where $A_1=0$, which means the temperature is equal to its saturation value. We then focus specifically on the impact of evaporation and vapor recoil, by setting $G=A_2=0$. With these assumptions, the evolution equation (Equation (56)) can be simplified as

$$h_t+\frac{E}{h}+ϵ\left\{\frac{2}{3}{E}^2{D}^{-1}(h_{xx}+h_{yy})+\frac{\overline{S}}{3}{\left[h^3(h_{xxx}+h_{xyy})\right]}_x+\frac{\overline{S}}{3}{\left[h^3(h_{xxy}+h_{yyy})\right]}_y\right\}=0.$$

(71)

The spatiotemporal evolution in both 2D and 3D simulations revealed that an increase in the evaporation number and vapor recoil could accelerate the breakdown of the film, as shown in Figure 7a,c. Specifically, when evaporation is the dominant factor (Figure 7a), the system exhibits greater stability in the 3D cases. On the other hand, when vapor recoil is the determining element (Figure 7c), the system becomes more unstable in 3D, which is matched with the predictions of linear stability theory. The results indicate that the 3D disturbances exacerbate the vapor recoil effect.

To gain a better understanding of the breakup process, we utilize a self-similar solution to extract the underlying scaling law:

$$h=\Delta t^{\zeta }F\left(\xi \right),\Delta t=t_r-t,\xi =\frac{z-z_r}{{(t_r-t)}^{\gamma }},$$

(72)

where $t_r$ is the breakup time and $z_r$ is the breakup location. The order of each term in Equation (71) yields

$$\begin{array}{cc}\hfill & h_t\sim \Delta t^{\zeta -1}{h}^{-1}\sim \Delta t^{-\zeta }{h}_{xx}+h_{yy}\sim \Delta t^{\zeta -2\gamma }\hfill \\ & {\left[h^3(h_{xxx}+h_{xyy})\right]}_x+{\left[h^3(h_{xxy}+h_{yyy})\right]}_y\sim \Delta t^{4\zeta -4\gamma }\hfill \end{array}$$

(73)

When the evaporation plays an important role, balancing the terms $h_t$, $h^{-1}$, and ${\left[h^3(h_{xxx}+h_{xyy})\right]}_x+{\left[h^3(h_{xxy}+h_{yyy})\right]}_y$ gives

$$\zeta =1/2,\gamma =5/8.$$

(74)

This means the minimal film thickness reduces to ${(t_r-t)}^{1/2}$, which is demonstrated in Figure 7b. When the vapor recoil effect is the dominant fact, balancing the terms of $h_t$, $h_{xx}+h_{yy}$, and ${\left[h^3(h_{xxx}+h_{xyy})\right]}_x+{\left[h^3(h_{xxy}+h_{yyy})\right]}_y$ gives

$$\zeta =1/3,\gamma =1/2.$$

(75)

This indicates that the film will break down within a finite time, and the minimum thickness will thin as ${(t_r-t)}^{1/3}$. The contribution under evaporation has a higher order, ${(t_r-t)}^{-1/3}$, than that under the vapor recoil effect, which gives ${(t_r-t)}^{-2/3}$. This implies that the scaling governing the breakdown process is not affected by the vapor recoil effect, and it is in line with the conclusion obtained by Wei and Duan [45].

Figure 8 illustrates the behavior of a nonvolatile film ($E=0$) in zero gravity. The film’s initial state is a uniform thickness of $h_0=1$, but it will be perturbed with a random amplitude of $0.1$. Over time, the film surface becomes gradually stabilized, with the peaks and valleys decreasing in amplitude. Eventually, the film will reach a fully stable state where the surface is once again flat with $h=1$. Since there are no external forces that could cause the film to become unstable in the absence of gravity and evaporation, this behavior is expected. Instead, any initial perturbations are simply damped out over time, leading to a fully stable state.

To study the impact of evaporation on a liquid film, we set the evaporation number to a fixed value of $\overline{E}=0.1$. Figure 9 demonstrates the gradual thinning of the film due to evaporation, with partial drying occurring around $t=5.4$ s. Additionally, the amplitude of the liquid film’s vibration steadily increases over time, indicating that the evaporation process significantly reduces the stability of the system.

Figure 10a presents the influence of the effective pressure, which illustrates that the film develops slower for larger values of $A_2$. This result is constant with the linear stability analysis that the system is stabilized as $A_2$ increases. However, when gravity is increased, the breakup time is delayed, which is contradicted by the prediction of the linear stability analysis. This different performance can be attributed to the nonlinear interactions between gravity and surface tension. As the film thickness decreases due to evaporation, the surface tension forces become more dominant, resulting in a more pronounced curvature of the free surface. This can amplify the effect of gravity and stabilize the film, which fails to achieve linear capture in the linear stability analysis. The delay in breakup time for larger gravity is consistent with the physical intuition that flow can enhance stability. The influence of energy transport on liquid film evolution is still small, especially in the case of microgravity, but overall, energy transport still has a small stabilizing effect, as shown in Figure 10b.

As shown in Figure 11, we investigate the evolution of a condensing film with $E=-0.1$ while fixing the gravity at $G=5$, which makes the film flow down in the x direction with a constant velocity. Moreover, at $t=2.5$ s, we notice that the film thickness decreases as $A_2$ increases. Figure 12 supports the conclusion by showing how the amplitude of film thickness reduces with an increase in $A_2$, proving that the effect pressure acts as a stabilizing factor.

In many practical situations, such as inkjet printing, the liquid is often deposited under a substrate, and the system is controlled by the Rayleigh–Taylor instability mechanism. Therefore, we investigated the performance of an inverted film as shown in Figure 13. It was observed that the film initially moves along the incline due to the force of gravity. At a certain point, a droplet appears on the film’s surface, which grows rapidly over a short period of time. The formation and growth of the droplet can be attributed to the destabilizing effect of gravity on the inverted film. Due to gravity, the thickness of the film varies along the incline, leading to non-uniform flow and the formation of a droplet. The layer thins as the droplet grows, causing an increase in the curvature of the droplet’s interface and accelerating its growth.

5. Conclusions

This study set out to study a three-dimensional evaporating/condensing liquid film with a nonlinear evolution equation under the consideration of the balance of the configurational momentum at the liquid–vapor interface. We considered effect pressure, vapor recoil, and gravity especially and discussed the impact of complex interactions on both linear stability and fully nonlinear dynamics with a long-wave model. Two terms were introduced in our model. One was the effective pressure, $p-E^2{D}^{-1}{J}^2$, accounting for vapor recoil, $E^2{D}^{-1}{J}^2$, and the other was the energy transport along the interface, $-N{\nabla }_s·{\mathit{u}}_s$ (${\mathit{u}}_s$ is the tangential component velocity at the liquid–vapor interface). The dimensionless parameter N accounts for the magnitude of the energy flux along the interface. The effect of the effective pressure and vapor recoil of the system on instability was studied using linear stability. The results indicate that the system is always stabilized by the effective pressure, whereas the vapour recoil effect acts to destabilize the system. Additionally, in this paper spanwise perturbation was considered as a stabilizing term in a three-dimensional nonvolatile film. However, the results demonstrate that vapour recoil reduces the stabilizing effect of spanwise perturbations in the presence of evaporation. For the influence of energy transport, while it had a very minor effect on the stability of the system, evaporation along the interface had an impact on the linear wave speed.

The nonlinear evolution of the nonvolatile, evaporating and condensing film falling down an inclined plane was examined using a Fourier spectral method, taking into account the influences of evaporation number, effective pressure, and vapor recoil. The nonlinear dynamics of the system reveal that effective pressure can enhance stability, thus delaying the rupture of an evaporating film. The study of the self-similar solution reveals that the minimal film thickness reduces as ${(t_r-t)}^{1/2}$ and ${(t_r-t)}^{1/3}$ under the dominance of evaporation and vapor recoil, respectively. The behavior of nonvolatile films under zero gravity was also discussed, showing signs of stabilization over time and reaching a fully stable state where the film surface is flat again.