Next Article in Journal
Forty Years’ Experience in Teaching Fluid Mechanics at Strasbourg University
Previous Article in Journal
A New Wall Model for Large Eddy Simulation of Separated Flows
Previous Article in Special Issue
Progress in Phenomenological Modeling of Turbulence Damping around a Two-Phase Interface
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Dynamics of Thin Film Under a Volatile Solvent Source Driven by a Constant Pressure Gradient Flow

1
Mathematics Department, University of Mauritius, Réduit 80837, Mauritius
2
Mechanical Engineering Department, University of Canterbury, Private Bag 4800, Christchurch 8140, New Zealand
3
Electrical and Computer Engineering Department, University of Canterbury, Private Bag 4800, Christchurch 8140, New Zealand
*
Author to whom correspondence should be addressed.
These authors contributed equally to this work.
Submission received: 20 June 2019 / Accepted: 19 November 2019 / Published: 29 November 2019
(This article belongs to the Special Issue Free surface flows)

Abstract

:
The evolution of a thin liquid film subject to a volatile solvent source and an air-blow effect which modifies locally the surface tension and leads to Marangoni-induced flow is shown to be governed by a degenerate fourth order nonlinear parabolic h-evolution equation of the type given by t h = div x M 1 h x 3 h + M 2 h x h + M 3 h , where the mobility terms M 1 h and M 2 h result from the presence of the source and M 3 h results from the air-blow effect. Various authors assume M 2 h 0 and exclude the air-blow effect into M 3 h . In this paper, the authors show that such assumption is not necessarily correct, and the inclusion of such effect does disturb the dynamics of the thin film. These emphasize the importance of the full definition t · grad γ = grad x γ + x h grad y γ of the surface tension gradient at the free surface in contrast to the truncated expression t · grad γ grad x γ employed by those authors and the effect of the air-blow flowing over the surface.

1. Introduction

In this work, we study the influence of a volatile solvent on the evolution of a thin liquid film on a solid surface. As the volatile solvent diffuses through the atmosphere, a non-uniform solvent concentration arises which induces surface tension gradient and Marangoni-driven flow [1,2,3]. This effect is well-known and was described in the pioneering works of Marangoni [4] and Thomson [5].
This physical phenomena is exploited in the industrial process, and is known as Marangoni drying whereby a jet of ethanol vapor (or other solvent) is blown onto a wet surface. As the ethanol interacts with the liquid free surface, a Marangoni flow is generated which helps drying the surface. This process has motivated a number of studies [6,7,8,9]. On experimental side, Marangoni drying was first studied by Leenaars et al. [6]. They coined the name Marangoni drying to this process. Leenaars et al. [6] have shown that this drying technique is more effective than spin drying which is usually carried out by centrifugation: Marangoni drying yields to extremely clean surface whereas spin drying can contaminate the surface. On the theoretical side, Matar and Craster [7] have given a mathematical description based on this drying process. The mathematical framework of Matar and Craster [7] permits the quantitative evaluation of the thinning process due to Marangoni drying. The model also makes it possible through parametric study to examine the response of the liquid film to optimize the thinning process. O’Brien [8] has considered the dynamics of a liquid film subjected to the action of alcohol-vapor-induced Marangoni effects using asymptotic methods. Taking into account only terms of leading orders into the generic equations, O’Brien reduces the mathematical description to a single non-linear partial differential equation whose characteristics are similar to the so-called Korteweg-de Vries and Burgers equations.
Preluding these studies, Carles and Cazabat had already demonstrated that the surrounding atmosphere plays an important role in the dynamics of a droplet spreading [10,11]. More specifically, these authors observe that the usual laws of spreading no longer hold. For some droplets, the spreading was strongly accelerated and accompanied by instabilities at the contact line, whereas for others, the process was somewhat reversed: retraction of droplets after a fast initial spreading. In particular, when a non-volatile droplet was investigated under an inhomogeneous surrounding of volatile solvent, they observe both the Marangoni effect and instabilities at its edge (see Reference [10]). In the second paper of Carles and Cazabat (see Reference [11]), the spreading dynamics of droplets under an inhomogeneous vapor phase of its volatile solvent were further investigated. The authors have investigated the spreading of an oil droplet under an atmosphere contaminated with a volatile solvent. Still, the observed behaviors were completely different from the normal ones. For, in some situations the surface tension gradient counter-accelerates the spreading mechanism while in others it accelerates instead, both along its directions. After a short time, when the interfacial effect was no longer operative, it was observed that the droplet either resumes spreading or retracts, again depending on the surface tension of the solvent (see Reference [11]).
Much research has been devoted to the understanding of Marangoni-driven thin liquid films. An exhaustive review of such studies is beyond the scope of this introduction but the interested reader is referred to the reviews by Oron et al. [12] and Craster and Matar [13]. Most of these study have considered a leading order expression for the surface tension gradient at the free surface. One aim of this paper is to explore the effect of these terms which have typically been disregarded by others.
Another important aspect in the study of Marangoni drying is the two way-coupling between the liquid film and the surrounding atmosphere. Few papers appear to have investigated this effect. One notable exception is the work of Sultan et al. [14] who studied the stability of a two-dimensional bi-layered liquid-vapor system over a solid substrate. In this paper, the authors generalize the one-sided study of Burelbach et al. [15] on the evaporation of thin film by including the diffusion of the vapor phase region. Their results are described in terms of both interfacial and mass transport phenomena. In contrast to earlier work, we consider here diffusion driven by an air-blow effect in the surrounding atmosphere and a two-way coupling whereby the air-blow effect in the vapor phase affects the liquid film.
This paper is organized as follows: In Section 2.1, the mathematical descriptions is exposed in dimensional form; the spatial description of a volatile solvent and the character of the interfacial force it induces are argued. Thereafter, Section 2.2 and Section 2.3 are devoted to perturbation theory so as to render the mathematical descriptions of earlier sections dimensionless and at leading order. Section 2.4 concerns the numerical side, yet some leading features are exposed in Appendix A. Finally, the results are discussed through a series of case studies under Section 3 and concluded on this basis under Section 4.

2. Methods

2.1. Mathematical Formulations

2.1.1. Modeling the Inhomogeneous Vapor Phase Region

Following the works of Carles et al. (see References [10,11]), we describe the presence of a volatile solvent in the atmosphere, region Ω v R 2 , by a volatile source function S x of the form
S x = μ s M s π λ v T exp x x 0 2 4 λ v T for   all   x = x , y Ω v .
In Equation (1), x 0 = x s , y s Ω v is the point at which the source is located. Further, the scalar quantities M s , μ s , λ v > 0 define the rate of increase in moles of solvent per unit length, the strength of the source and the diffusion constant of the volatile solvent vapor. The time scale T > 0 serves to establish a diffusion length l T = 2 λ v T . In the following, we will consider the transport of volatile solvent from a point source given by Equation (1).

2.1.2. On the Surface Tension and Its Derivative

For some real t < , let the set I t = 0 , t denote the period of time over which the study is carried out. Suppose the scalar field function c l t , x be the spatio-temporal distribution of the chemical concentration in the liquid domain Ω l , where t , x I t × Ω l . Then, according to the experimental work of Carles et al. [10], the closure relationship relating the scalar field function c l t , x to the surface tension γ t , x at every t , x I t × Γ 4 is linear and hence,
γ t , x = κ 1 + κ 2 c l t , x in t , x I t × Γ 4 ,
where κ 1 > 0 and κ 2 0 are property-dependant constants. To this end, the directional derivative of the scalar field function γ t , x along the unit vector e x yields to
grad x γ = κ 2 x c l + x h y c l for   all   t , x I t × Γ 4 ,
because c l = c l x , y = h , t and h = h x , t on Γ 4 . To the best of our knowledge, the effect of the term x h y c l was disregarded in the past. In the sequel, we will show that it has a non-negligible effect. Before moving on to the next section, we deduce the constants κ 1 , κ 2 of the above expression. Let γ 0 , l > 0 be the surface tensions of the liquid in its uncontaminated state. Then, it is evident that one of these constants describes the surface tension of the liquid film when the liquid film is chemically stable, and the other equates the surface tension gradient of the liquid whenever that stable condition is being disturbed chemically. Thus, these equalities follow:
κ 1 = γ t , x = γ 0 , l for   all   t , x 0 × Γ 4 , κ 2 = D c l γ t , x for   all   t , x I t × Γ 4 ,
respectively, where D η · denotes the ordinary derivative operator with respect to η ; that η = c l in the above description is understood but when η = t , D η · defines itself as the Lagrangian derivative operator. Consequently, it results by substitution that,
γ t , x = γ 0 , l + D c l γ c l t , x for   all   t , x I t × Γ 4 .

2.1.3. Constitutive Equations

Consider a thin liquid film initially sitting on a solid substrate Γ 2 , having a thickness h 0 > 0 and a length L > 0 (with ord h 0 / L 10 3 , see References [7,8]). For mathematical treatment this liquid film is designated by Ω l R 2 and bounded by Ω l = i = 1 4 Γ i . Henceforth, subscripts v and l indicate the vapor and liquid phases, respectively. The vapor phase surrounding is defined by Ω v R 2 and bounded by Ω v = i = 4 7 Γ i . The volatile solvent source is placed at x 0 = x s , y s Ω v ; Figure 1 supports such geometrical descriptions.
The presence of the source S x gives rise to advection-convection phenomena. Therefore, if c k t , x designates the chemical concentration in Ω k , where t , x I t × Ω k and k l , v , then the following advection-diffusion equations result in:
D t c k t , x = λ k div c k + δ k v S x for   all   t , x I t × Ω k ,
where the bi-subscripts δ k v 0 , 1 is the Kronecker delta (with k l , v ). The operator D t · = t · + u · grad · designates the Lagrangian derivative operator. Evidently, λ l > 0 is the diffusion constant of the volatile solvent into the liquid film, and u = u t , x for all t , x I t × Ω k (with k l , v ); u t , x is called the velocity vector field distribution within either Ω l or Ω v . For every k l , v , the flow field u t , x within Ω k is described by the incompressible Navier–Stokes equations of constant fluids properties in a gravitational field, which in vector invariant forms read:
D t u t , x = 1 ρ div T + f for   all   t , x I t × Ω k , div u = 0 for   all   t , x I t × Ω k ,
where the tensor field T p , u , which stands for the total stress tensor, is given by
T p , u = p I + 2 μ E u , E u = 1 2 grad u + grad u T .
In Equation (4), tensors I and E u , respectively, are the identity and strain-rate tensors, and fields p t , x , u t , x = u x , u y T and f t , x = f x , f y T , respectively, denote the pressure, velocity, and body force (per unit mass). The property constants ρ > 0 and μ > 0 , respectively, define the density and dynamic viscosity of the k th fluid (with k = l for liquid and k = v for vapor phases). The pressure difference sustained across the interface Γ 4 is modeled with the Young-Laplace equation. Thus,
p l t , x = p v γ div n for   all   t , x I t × Γ 4 .
(We note here that grad p l = grad γ · n if the atmosphere has constant pressure; in the sequel, an alternative form for grad p v is established such that the latter operates as an air-blow effect.) In addition, it is understood that by the vector field u t , x defined in Ω l , we implicitly mean u l t , x , and so forth for other property constants and variables. The vector fields n t , x , t t , x are the outer unit normal and tangent vectors of their boundary. In particular, on Γ 4 (see Figure 1), we obtain
n t , x = x h , 1 T 1 + x h 2 and t t , x = 1 , x h T 1 + x h 2 for   all   t , x I t × Γ 4 .

2.1.4. Initial and Boundary Conditions

At time t = 0 , the domains Ω l and Ω v are initially motionless and free from or low in volatile substances. Thus, these initial conditions follow:
u k t , x = 0 for   every   t , x , k 0 × Γ 4 × l , v , c k t , x = c 0 k for   every   t , x , k 0 × Γ 4 × l , v .
In these conditions, c 0 l , c 0 v 0 are arbitrary constants which will be used in the nondimensionalization process so as to obtain a zero concentration in dimensionless forms. On boundaries Γ 5 , Γ 7 , we impose a constant pressure gradient in the atmosphere. If K v > 0 stands for such constant and · denotes the Euclidean norm in R 2 , then
grad p v = K v for   all   t , x I t × Ω v .
On the boundaries Γ 1 , Γ 3 of the domain Ω l , we prescribe a so-called normal stress-free boundary conditions. Thus, for every i 1 , 3 ,
n · T n = p 0 l for   all   t , x I t × Γ i ,
where the scalar field p 0 l > 0 is an arbitrary constant and n t , x is the outer unit normal vector of the boundary in question. For the chemical concentrations, we impose a flux continuity boundary condition on the interface Γ 4 ,
λ l n · grad c l = λ v n · grad c v for   all   t , x I t × Γ 4 .
Further, we assume the mass transfer process to vanish on Γ 6 and to result in a zero-flux boundary conditions on boundaries Γ 1 to Γ 3 and Γ 5 to Γ 7 . Thus,
λ l n · grad c l = 0 for   all   t , x I t × i = 1 , 2 , 3 Γ i , λ v n · grad c v = 0 for   all   t , x I t × i = 5 , 6 , 7 Γ i , c v = c 0 v for   all   t , x I t × Γ 6 .
On the surface Γ 2 , we consider a Navier slip-with-friction conditions,
t · E u n + β u = 0 for   all   t , x I t × Γ 2 , n · u = 0 for   all   t , x I t × Γ 2 ,
where β > 0 is the Navier-slip coefficient. For, the movement of the liquid film entails a singularity of stress at the contact line if the usual no-slip boundary condition, i.e., t · E u n = 0 on Γ 2 , is imposed between the liquid film and the solid substrate (see Reference [16] and references therein). On the interface Γ 4 , we require that the substantial derivative of y h x , t vanishes. Therefore,
D t y h = 0 for   all   t , x I t × Γ 4 .
Finally, to ensure continuity upon tangential shear stress on Γ 4 , we consider
t · T n = t · D c l γ grad c l for   all   t , x I t × Γ 4 .
A remark deserves attention at this point. It is very true that Figure 1 does not indicate the presence of a contact line and, therefore, the discussion of the Navier slip-with-friction condition seems odd. The authors agreed that the Navier-slip-with-friction conditions on Γ 2 should instead be the classical no penetration boundary conditions; that is, the normal and tangential velocity vector fields vanish for all t , x I t × Γ 2 , n · u = 0 and t · u = 0 . However, we wish to derive a degenerate fourth order nonlinear parabolic h-evolution equation that will apply not only in the case considered here (Navier-slip-with-friction coefficient β > 0 set to zero) but also to the case in which a contact line (or trijunction) exists and a liquid on a solid substrate spreads and displaces the surrounding fluid (say, vapor phase region Ω v ).

2.2. Scaling Analysis

2.2.1. Asymptotic Approximations

Set ϵ = h 0 / L 1 and, by hypothesis, assume ord h 0 / L 10 3 ; obviously, ϵ denotes the asymptotic parameter. Geometrical properties are scaled based on the initial film thickness h 0 = ϵ L and physical properties are on the contrary scaled based on interfacial and viscous forces. On this basis, it is correct to consider ord x , y T L , ϵ L T within the domain Ω l . Let u 0 l > 0 characterize the flow field in Ω l . Then, by Equation (4), it follows that ord u u 0 l , ϵ u 0 l T in Ω l . Accordingly, ord t L / u 0 l and ord p μ u 0 l / L . On the interface Γ 4 , dynamic equilibrium requires that ord p ϵ 1 γ 0 , l γ 0 , v / L ; hence, ord u 0 l ϵ 1 γ 0 , l γ 0 , v / μ . On the whole, the set of characteristic scales in Ω l reads
Fr l = def ϵ 2 γ 0 , l γ 0 , v 2 μ 2 g L , Pe l = def ϵ γ 0 , l γ 0 , v L μ λ l , Re l = def ϵ γ 0 , l γ 0 , v L μ ν , t ϵ 1 μ L γ 0 , l γ 0 , v t , c l c 0 l + c l c 0 l c l , x , y T L x , ϵ y T , p ϵ 1 γ 0 , l γ 0 , v L p , u x , u y T = ϵ γ 0 , l γ 0 , v μ u x , ϵ u y T .
By q α q is meant that the quantity q is dimensional if it precedes ↦ and dimensionless if it follows ↦; the quantity α is obviously dimensional. The parameters Fr l , Pe l , Re l > 0 are, respectively, the Froude, Péclet, and Reynolds numbers in Ω l . By substitution of the above scalings into the dimensional equations, one obtains a system of dimensionless variables equations. In Ω v , it is obvious that ord x , y T L , L T . Consequently, ord u u 0 v , u 0 v T in Ω v , where the scalar u 0 v > 0 designates a characteristic flow field in Ω v . In addition, ord S λ v c v c 0 v / L 2 because the spreading of the volatile solvent is by diffusion process. On the whole, the following scalings hold reasonable in the domain Ω v :
Fr v = def u 0 v 2 g L , Pe v = def u 0 v L λ v , Re v = def u 0 v L ν , t L u 0 v t , c v c 0 v + c v c 0 v c v , S λ v L 2 ( c v c 0 v ) S , x , y T L x , y T , u x , u y T u 0 v u x , u y T , p ρ u 0 v 2 p .
Substituting the above scalings into Equation (3), there results in
D t c v = 1 Pe v ( div c v + S x ) for   all   t , x I t × Ω v , D t c l = ϵ 2 Pe l ϵ 2 x 2 c l + y 2 c l for   all   t , x I t × Ω l .
In a like manner, the dimensionless forms of Equation (4) resolve into
div u = 0 for   all   t , x I t × Ω l , D t u x = ϵ 2 Re l x p + ϵ 2 x 2 u x + y 2 u x for   all   t , x I t × Ω l , D t u y = ϵ 4 Re l y p + ϵ 4 x 2 u y + ϵ 2 y 2 u y ϵ 1 Fr l for   all   t , x I t × Ω l ;
and
div u = 0 for   all   t , x I t × Ω v , D t u x = x p + 1 Re v x 2 u x + y 2 u y for   all   t , x I t × Ω v , D t u y = y p + 1 Re v x 2 u y + y 2 u y 1 Fr v for   all   t , x I t × Ω v ,
respectively. The surface tension of the above dimensionless scaling set is termed the spreading pressure, π s = γ 0 , l γ 0 , v [7], where γ 0 , k > 0 is the surface tension of the k th fluid (with k = l , v ), as already pointed out. The quantities c 0 k , c k are prescribed values characterizing the initial and final values of the concentrations of the k th fluids (with k = l , v ). The control parameters Fr v , Pe v , Re v > 0 , respectively, are obviously the Froude, Péclet, and Reynolds numbers in Ω v .
It now remains to nondimensionalize the set of initial and boundary conditions. Beginning with the initial conditions given in Section 2.1.4 for the field functions u k t , x , c k t , x (with k = l , v), respectively, we obtain, for every k l , v ,
u k t , x = 0 for   all   t , x 0 × Ω k , c k t , x = 0 for   all   t , x 0 × Ω k ,
where the property constants c 0 k , c k (with k = l , v ) have now disappeared. Likewise, nondimensionalizing the normal stress-free boundary condition given in Section 2.1.4 yields to p + x u x = P l . Imposing on boundaries Γ 1 , Γ 3 , respectively, a zero viscous stress along with a Dirichlet condition on p x , t gives
p t , x = p 0 l for   all   t , x I t × i = 1 , 3 Γ i .
To nondimensionalize the flux-continuity boundary condition, one should be careful. In fact, away from Γ 4 , ord L 1 , L 1 T in Ω v , whereas immediately above Γ 4 , ord L 1 , ϵ 1 L 1 T . Therefore, letting λ l / λ v = c v c 0 v / c l c 0 l , the dimensionless form of the flux-continuity boundary condition writes
ϵ 2 x h x c l + y c l = ϵ 2 x h x c v + y c v for   all   t , x I t × Γ 4 .
In like manner, after nondimensionalizing the other boundary conditions given in Section 2.1.4, one deduces
x c l = 0 for   all   t , x I t × i = 1 , 2 , 3 Γ i , x c v = 0 for   all   t , x I t × i = 5 , 6 , 7 Γ i , c v = c 0 v for   all   t , x I t × Γ 6 .
For the dimensionless form of the so-called Navier slip-with-friction conditions, we proceed in two steps. In the first place, we nondimensionalize the conditions given there. Then, we replace the Navier-slip coefficient β by 2 ϵ β L 1 . This gives
u x = β y u x + ϵ 2 x u y for   all   t , x I t × Γ 2 , u y = 0 for   all   t , x I t × Γ 2 .
Likewise, after nondimensionalizing the kinematic boundary conditions given in Section 2.1.4, we obtain
u y = t h + u x x h for   all   t , x I t × Γ 4 , y u x = x c l + x h y c l D c l γ + ϵ 2 2 x h x u x 2 x h y u y x u y + x h 2 y u x + ϵ 4 x h 2 x u y for   all   t , x I t × Γ 4 .
Note that the dimensionless unit normal vector and tangential vectors n t , x and t t , x , respectively, defined in Section 2.1.3 are approximated as follows, ord n ϵ x h , 1 T , ord t 1 , ϵ x h T , respectively. We now wish to combine Equation (2) with Equation (5), given in Section 2.1.3. To do so, define two dimensionless numbers, namely, the Capillary (Ca) and Marangoni (Ma) numbers as
Ca = γ 0 , l γ 0 , v L μ λ l and Ma = D c l γ μ λ l c l c 0 l L ,
respectively. The construction of these dimensionless numbers Ca, Ma > 0 is based on the following line of thought. Since ord u 0 l ϵ γ 0 , l γ 0 , v / μ and ord γ ϵ μ λ l / L , in order to establish the Capillary number Ca, it suffices to substitute ϵ γ 0 , l γ 0 , v / μ and ϵ μ λ l / L for u 0 l and γ , respectively, into Ca = μ u 0 l / γ . On the one hand, the Marangoni number Ma being regarded as proportional to (thermal-)surface tension forces divided by viscous forces, i.e., Ma = T γ T T 0 L / μ α , where T and α , respectively, designate the temperature and thermal diffusivity, so to adapt this dimensionless number with the problem under consideration, it suffices to replace the field T by c l , and the coefficient α by λ l . Thence, the dimensionless numbers so far defined. Using those dimensionless numbers into Equation (2), the resulting expression in dimensional form reads,
γ c l = γ 0 , l + Ma Ca γ 0 , l γ 0 , v c l c 0 l c l t , x for   all   t , x I t × Γ 4 .
Combining the above expression with Equation (5), the resulting dimensional expression writes
p l = p v + γ 0 , l + Ma Ca γ 0 , l γ 0 , v c l c 0 l c l x 2 h for   all   t , x I t × Γ 4 .
Let us now investigate into p v t , x . In the immediate vicinity of every point x = x , h t , x Γ 4 , it is clear that
ord δ p = ord p l p v ϵ 1 γ 0 , l γ 0 , v L for   all   t , x I t × Γ 4 ,
which, in other words, suggests to scale δ p as such, δ p ϵ 1 γ 0 , l γ 0 , v δ p / L ; obviously, the differential quantity δ p is dimensional if it precedes ↦ and dimensionless otherwise. Hence, taking into account the characteristic scales given in Section 2.2.1, setting γ 0 , l γ 0 , l γ 0 , l γ 0 , v ( γ is dimensional if it precedes ↦ and dimensionless otherwise) and then rescaling the dimensionless constant property γ 0 , l by γ 0 , l Ma / Ca c 0 l / c l c 0 l , the dimensionless form of the pressure field p l t , x resolves into
p l p v + ϵ 2 γ 0 , l + Ma Ca c l x 2 h for   all   t , x I t × Γ 4 .
The dimensionless form of the (chemical)-surface tension is now given by the expression γ c l = γ 0 , l + Ma / Ca c l .

2.3. Leading Order Model

To begin, some simple analysis is now in order. With L 10 4 [m], u 0 v 10 4 [m·s−1] and λ v 10 5 [m2·s−1], one obtains ord Pe v 10 3 ; hence, ord Pe v ϵ . (These numerical values are extracted from References [7,8].) On the other hand, we may assert that ord Pe l 1 ϵ 2 , since the spreading pressure π s is a predominant factor. At this point, we are now well-prepared for the search of a simplified model.
In the limit that the quantities ϵ 2 and Pe v , respectively, tend to zero and to unity, the leading order forms of the dimensionless equations for the fields c k t , x (with k = l , v) reduce to
div c v + S x = 0 for   all   t , x I t × Ω v , D t c l y 2 c l = 0 for   all   t , x I t × Ω l .
Some analysis is now in order to support the assumption made upon Pe l . Suppose that the flow induced within the region Ω l is an ϵ k -order of the Péclet number, i.e., ord Pe l ϵ k (with k 0 ). Then, if we express Fr l , Pe l and Re l in terms of the quantity ϵ γ 0 , l γ 0 , v / μ it can be shown that
ord ϵ k L / λ l ord Re l L / ν ord Fr l 1 / g L ( w i t h k 0 ) .
And, consequently,
ord ϵ 2 k Re l ord ϵ k Re l ord ν λ l ord ν g ϵ 2 k Fr l ord ϵ k 2 Fr l ,
where ν ϵ k g (with k 0 ). The above approximation enables us to obtain some information as to the nature of the order of magnitube of one dimensionless number with respect to the other. Moreover, for the case of water, we found k 2 , based on the data of O’Brien [8]; hence, the assertion made in so far. With this view in mind, taking k 2 there follows that
ord ϵ 4 Re l ord ϵ 2 Re l ord ϵ 1 Fr l ,
which simply asserts a small Bond number ( Bo l ), since Bo l = ϵ Re l / Fr l 1 . Consequently, the system of Equation (7) reduces to
x u x + y u y = 0 for   all   t , x I t × Ω l , y 2 u x x p = 0 for   all   t , x I t × Ω l , y p = 0 for   all   t , x I t × Ω l .
On the other hand, paying little attention to the effects arising from body and viscous forces in Ω l , i.e., Fr v 1 , Re v 1 1 , the system of Equation (8) results in
x u x + y u y = 0 for   all   t , x I t × Ω v , D t u x + x p = 0 for   all   t , x I t × Ω v , D t u y + y p = 0 for   all   t , x I t × Ω v .

2.3.1. Alternative Form for the Pressure Field in the Vapor Phase

We postulate that the pressure field in Ω v varies according to
p v t , x ; ϵ p 0 v + K v x + ϵ j p j v t , x for   all   t , x I t × Ω v ,
where j 1 . Consequently, the leading order of the flow equations in Ω v in vector form writes
D t u v t , x ; ϵ K v e x + ϵ j p j v t , x for   all   t , x I t × Ω v ,
where j 1 . Therefore, D t u v + K v e x = 0 at leading order; that is, we force the acceleration of every fluid parcel carrying the volatile solvent to depart from its respective trajectory; thence, a translational movement of the volatile substances towards the x-direction.

2.3.2. Method of Solution

Computing the system Equation (10) using the leading order Neumann and Dirichlet kinematic boundary conditions given in Section 2.2.1, one obtains
u x t , x = x p 2 2 h y + β y 2 + Ma Ca L h c l y + β , u y t , x = x 2 p 6 y 3 3 h y + 2 β y + 1 2 Ma Ca L h 2 c l y + 2 β y + x p 2 x h y + 2 β y for   all   t , x I t × Ω l ,
where L h · = x + x h y · , and
L h 2 · x + x h y 2 · = x 2 + x h 2 y 2 + 2 x h x y + x 2 h y · .
To settle the spatial derivative x p appearing in u x , t , the pressure field model p v t , x ; ϵ is invoked into the pressure expression p l t , x , giving
x p l ϵ 2 x γ 0 , l + Ma Ca c l x 2 h + K v + ϵ k x p j v .
where j 1 . To augment the interfacial hydrodynamics by an ϵ 2 -order, we assume that ord γ 0 , l ϵ 2 . Consequently, this assertion suggests to rescale the surface tension of the liquid film by substituting ϵ 2 γ 0 , l for γ 0 , l . Thence,
x p l t , x K v + γ 0 , l x 3 h t , x for   all   t , x I t × Γ 4 .
A similar approach was suggested by Oron et al. [12]. In the following section, the h-evolution equation is established.

2.3.3. Evolution Equation for the Interface

Using the general form of the Leibniz integral rule, one obtains
t h = div x 0 h u x d y for   all   t , x I t × Γ 2 .
When all these are considered, we deduce after integration the h-evolution equation, which writes
t h = div x h 3 3 1 + 3 β h K v + γ 0 , l x 3 h div x h 2 2 1 + 2 β h Ma Ca L h c l for   all   t , x I t × Γ 2 .
In particular, it may well be pointed out here that when the hitherto neglected term x h y c l , the induced pressure gradient x p v = K v , and the Navier-Slip coefficient β > 0 are disregarded into Equation (13), we deduce the well-known extended thin film model, written as
t h = div h 2 2 μ γ h 3 3 μ p for   all   t , x I t × Γ 2 ,
in its general and dimensional form.

2.4. Numerical Method

2.4.1. Preliminaries

Define three mobility terms, namely, M 1 h , M 2 h , and M 3 h , respectively, as follows
M 1 h = def h 3 3 1 + 3 β h γ 0 , l , M 2 h = def h 2 2 1 + 2 β h Ma Ca y c l , M 3 h = def h 2 2 1 + 2 β h Ma Ca x c l h 3 3 1 + 3 β h K v ,
so that Equation (13) may be rewritten as such:
t h = div x M 1 h x 3 h + M 2 h x h + M 3 h for   all   t , x I t × Γ 2 .
Having settled the h-evolution equation in condensed form, we now set forth the latter into a weak differential formulation by re-defining the functions h, x 2 h as thus, h , x 2 h = h , h ; the idea of doing so is to help numerical implementation. This yields to the following h , h -system:
t h = x M 1 h x h + M 2 h x h + M 3 h ,
0 = x 2 h + h for   all   t , x I t × Γ 2 ,
which is now a degenerate second order nonlinear parabolic system.
To adapt the situation, we construct at every point x of the surface Γ 2 and time t a depth average flow field u ^ = u ^ x , u ^ y T with the velocity field u t , x as follows
u ^ t , x , h = 1 h 0 h u t , x , y d y for   all   t , x I t × Γ 2 .
Substituting the vector field u x , t , established in Section 2.3.2, into the above integral yields
u ^ x t , x , h = 1 2 Ma Ca h + 2 β L h c l h 2 h + 4 β x p ,
u ^ y t , x , h = 1 6 Ma Ca h 2 h + 3 β L h 2 c l + h 2 6 h + 3 β x h x p h 3 8 h + 4 β x 2 p for   all   t , x I t × Γ 2 .
Having laid enough foundations, we are now well-prepared to address the problem numerically; this is the duty of what follows.

2.4.2. Computing with the COMSOL Multiphysics Software

In summary, the flow of interest takes place in the vapor domain Ω v and the liquid domain Ω l with a strong coupling through the free surface Γ 4 among the fields c l t , x , c v t , x , and h t , x . Let their respective systems of equations be designated by S c v , S c v , and S h , respectively. Then, for the evolution of the chemical concentration function c l t , x into Ω l , the corresponding system S c l writes
S c l : D ^ t c l y 2 c l = 0 for   all   t , x I t × Ω l , y c l y c v = 0 for   all   t , x I t × Γ 4 , y c v = 0 for   all   t , x I t × i = 1 , 2 , 3 Γ i ,
where the hatted operator D ^ t · = t · + u ^ · grad · designates the Lagrangian derivative operator subjected to the previously established depth averaged flow field u ^ = u ^ x , u ^ y T . For the evolution of the chemical concentration function c v t , x in Ω v , the corersponding system S c v writes
S c v : div c v + S = 0 for   all   t , x I t × Ω v , y c v y c l = 0 for   all   t , x I t × Γ 4 , y c v = 0 for   all   t , x I t × i = 5 , 7 Γ i , c v = 0 for   all   t , x I t × Γ 6 ,
and, finally, for the evolution of the free-surface function h t , x , the corresponding system writes
S h : t h + div x F = 0 for   all   t , x I t × Γ 2 , x h = 0 for   all   t , x I t × 0 , L ,
where the scalar function F h , h , designating the flux functions of the thin film equation, is defined as
F h , h = def M 1 h x h + M 2 h x h + M 3 h for   all   t , x I t × Γ 2 .
Appendix A gives an account on Modeling Methodology using the COMSOL Multiphysics interfaces. With all these features settled conveniently into appropriate COMSOL Multiphysics interfaces, we ran several numerical simulations, whose solutions are discussed in the next section.

3. Results and Discussion

Herein, the results are discussed, emphasizing the hitherto neglected term x h y γ and the hitherto overlooked term grad p v = K v (air-blow effect; K v , a prescribed constant).

3.1. Preliminaries

By an interfacial surface tension gradient is meant the vector gradient Γ 4 γ , where the operator grad Γ 4 · = I n n grad · is termed the Γ 4 -interface gradient operator on Γ 4 ; ⊗ is called the tensor product operator: for given vector fields v = v x , v y , w = w x , w y R 2 , their tensor product results in a 3 × 3 matrix with the entries v w = v i w j i , j x , y 2 . Let grad Γ 4 be its characteristic scale. Then, grad Γ 4 · grad Γ 4 grad Γ 4 · . (Note: grad Γ 4 · is dimensional if it precedes ↦ and dimensionless otherwise.) After nondimensionalizing the latter, this results in
grad Γ 4 grad Γ 4 · 1 L e x x + e x x h y + ϵ L e y x + ϵ 2 L e y x h 2 x .
Further, asserting that ord grad Γ 4 L 1 —which is reasonable—yields
grad Γ 4 · e x x + x h y · ( = e x L h · ) ,
at leading order. With these features beforehand together with the following definition and approximation,
e x · grad · = · , · grad Γ 4 · · ,
respectively, there follows that grad Γ 4 γ γ + γ . Upon perusal of earlier works, one finds that the term grad Γ 4 γ was approximated as follows: grad Γ 4 γ γ . In other words, their authors considered thin film equations of the form
t h γ = 0 = H π ; h , γ for   all   t , x I t × Γ 2 .
For example, in the paper of O’Brien [8], the explicit form of H π ; h , γ at leading order—which is the ϵ 0 -order treated in the present work—reads
H π ; h , γ = h 2 2 γ for   all   t , x I t × Γ 2 ,
whereas, in the present work the following extended model is proposed, namely,
t h γ 0 = H π ; h , γ + δ H π ; h , γ , p v for   all   t , x I t × Γ 2 ,
where the property π stands as the set of control parameters, for instance, ours reads π = def x 0 , μ s , Ca , Ma ; and the differential δ H goes to zero whenever γ and p v simultaneously go to zero. Thus, to distinguish our results with respect to others, it suffices to demonstrate through numerical experiments the following inequality:
t h γ 0 t h γ = 0 = δ H π ; h , γ , γ , p v 0 ,
Graphically, the effects/intepretations of the terms γ = e x x γ , γ = e x x h y γ and p v = K v e x , and δ H can be discussed by the aid of two new concepts, namely, coarse, fine approximations.
Physically speaking, if an orthogonal projection is considered the vector quantity Γ 4 γ approximates to Γ 4 γ γ . This is the case of a coarse approximation of Γ 4 γ . The works of Matar et al. [7] and O’Brien [8] are, among others, two examples in which are found coarse approximation. Contrarily, if a rotational projection is adopted instead, this would result to Γ 4 γ γ + γ . Thus, the h-evolution proposed here yields a better approximation, termed fine approximation. The left-hand side of Figure 2 illustrates those two concepts. In addition, when the free surface of the liquid film is subjected to an air-blow effect, p v = K v ( K v constant), the right-hand side of Figure 2 follows.
Throughout what follows, by the statement steady-state regime/profile is meant that there exists a property t > 0 such that the property h t , x at every point x Γ 2 is invariant for all t t .

3.2. Case 1. The h-Evolution Equation versus the Term γ

3.2.1. On Spatial Inhomogeneity

We describe here the h-evolution equations h γ = 0 and h γ 0 to stress on the effects of an inhomogeneous vapor phase and the consequence entails when γ 0 . The plots of h γ = 0 and h γ 0 are depicted in Figure 3.
It will be remarked in Figure 3 that in both cases the observed film thickness distribution is somewhat very similar, differing only in magnitudes. Consequently, this put forward the following fact:
| h γ 0 h γ = 0 | 0 ,
which demonstrates that the term γ t , x does affect the solution of the thin film proplem. It is, therefore, of interest to examine in which situations the term γ t , x is of negligible consequence, which is argued in what follows.

3.2.2. Of the Strength of the Source

What happens when we augment the strength μ s 0 of the source is here discussed. We carried out two parametric studies, one including and the other excluding the term γ t , x , with respect to μ s ranging from μ s = 10 to μ s = 41 . In Figure 4 are exposed the graphical results which we shall now discuss.
Glancing at Figure 4, one sees that the evolution of min h μ s with respect to μ s is linear only when γ x , t = 0 . So to speak, when γ x , t 0 the quantity min h t , x decays nonlinearly. On this basis, if the nonlinear effect arises when γ t , x is included in the mathematical model, obviously the rate at which the liquid film will dry or coat its substrate will be faster. Consequently, the higher be the effect of the volatile solvent, i.e., μ s 1 , the more effective will be the effect of the term γ t , x ; thus, this shows that its presence does affect the solution of the problem.

3.2.3. A Liquid Film of Infinitesimal Thickness

It be interesting to see how fast the quantity min h t , x goes to zero. In this section, upon the following substitution, 0.15 h 0 for h 0 , we consider an extreme case, a liquid film of infinitesimal thickness, to describe these situations. With these features beforehand, two simulations were ran, of which one takes account of the term γ t , x and the other disregards it. The resulting superimposed curves are shown in Figure 5.
Of these curves, the lower curve—the one such that min h t , x = 0 (i.e., h γ 0 )—describes the case for which γ t , x 0 ; the upper curve—the other one such that min h t , x 0 (i.e., h γ = 0 )—describes the case for which γ t , x = 0 . Physically speaking, it would appear that if the liquid film is of infinitesimal thickness, the higher the magnitude of the diffusive flux, the larger will be the deformation of the liquid film.
For that matter of comparison, cases for which γ 0 and γ = 0 are exposed in Figure 6. As can be easily seen in Figure 6, the curve h γ 0 goes to zero faster than h γ = 0 .
Numerically, we found in the linear regime that
min h γ 0 0.08 t + 0.14 , min h γ = 0 0.012 t + 0.03 , ( only   in   the   l i n e a r   regime ) .
As a result, both expressions can, therefore, be described by a linear (dimensionless) relationship of the form h t , x = L / 2 = π 1 t + π 2 , where π 1 , π 2 > 0 are property dependent constants.
At this point, we thought it desirable to ascertain how long does min h t , x take to goes to zero. Based on the foregoing statements, it would also appear that the time at which min h γ 0 0 is t 2.5 and that at which min h γ 0 0 is t 1.75 - so to speak, a discrepancy of about 30%. Recall the characteristic scales of t and h t , x , namely, t ϵ 1 μ L / γ 0 , l γ 0 , v and h ϵ L , respectively, the dimensional form of h t , x = L / 2 = π 1 · t + π 2 writes
h t , x = L / 2 = ϵ 2 π 1 μ γ 0 , l γ 0 , v t + ϵ L π 2 .
In particular, if τ ref > 0 designates a reference time at which min h t , x vanishes, it follows that
τ ref γ 0 π 2 π 1 ϵ 1 μ L γ 0 , l γ 0 , v γ 0 < π 2 π 1 ϵ 1 μ L γ 0 , l γ 0 , v γ = 0 τ ref γ = 0 ;
thence, knowing π 1 ( = min t h ), π 2 ( α h 0 ; α 1 , a correcting factor fitting the above linear law), one can deduce τ ref γ 0 , τ ref γ = 0 , and conversely. Using the following numerical values, γ 0 , l = 7.26 × 10 2 [N·m−1], γ 0 , v = 3.26 × 10 2 [N·m−1], μ 10 3 [kg·m−1·s−1], and L = 5 × 10 2 [m], respectively, we obtain
τ ref γ 0 2.18   [ s ] and τ ref γ = 0 3.13   [ s ] ,
respectively. Thus, an error of one second in the time at which min h x , t 0 will account for the term γ . (Of these numerical values, the values of the physical property constants γ 0 , l , γ 0 , v , and μ are those found in the paper of [8].)

3.2.4. On Capillary and Marangoni Effects

The impact of Capillary (Ca) and Marangoni (Ma) numbers on the behavior of the thin film equation is now discussed. The resulting graphical results are shown in Figure 7, Figure 8 and Figure 9, respectively. On inspecting those figures, it is easily seen that, on augmenting the number Ca the liquid film thins (Figure 8), while on the contrary, on incrementing the number Ma, it thickens (Figure 9). Thus, out of these situations, it follows that, by careful variations of the control parameters Ca, Ma, respectively, one can regulate the thinning rate.
As can been easily seen, Figure 7 shows that the evolution of the quantity min h Ma thins almost linearly with respect to Ma when γ t , x = 0 . But, since this is not the case when such term is taken into account, there obviously follows that the higher be the quantity Ma the more effective be consequence arising from γ t , x . Clearly, in the context of so-called Marangoni Drying, it acts as a additive effect, and, therefore, should not be left out of account.
Since the Capillary number Ca is a measure of the relative effect of viscous effects with respect to surface tension forces acting across the interface Γ 4 , an interesting way of interpreting this number is as follows: Let F μ , F γ denote the viscous and surface tension forces. Write Ca ref = F μ / F γ to designate a reference Capillary number. Then, variations of these forces by δ F μ , δ F γ results in the following strict inequalities
F μ F γ + δ F γ < Ca ref < F μ + δ F μ F γ .
The resulting curves are exposed in Figure 8 for Ca = Ca ref ± 0.5 , where Ca ref = 1.0 .
On the other hand, the Marangoni number Ma may be regarded as that dimensionless quantity which describes the relative effect of chemical surface tension forces with respect to viscous forces acting across the interface Γ 4 . Thus, if the quantity F γ c l designates the chemical surface tension force, and δ F γ c l an increment of it, the following strict inequalities hold
F γ c l F μ + δ F μ < Ma ref < F γ c l + δ F γ c l F μ .
In the above condition, the quantity Ma ref = F γ c l / F μ designate a reference Marangoni number.
The resulting curves are depicted in Figure 9 for Ma = Ma ref ± 10 , where Ma ref = 80 .

3.3. Case 2. The h-Evolution Equation Versus the Term p v

3.3.1. Effect of an Air Blow with γ = 0

When volatile solvent issues under a constant pressure gradient, i.e., p v K v , from the source S x in the vapor phase region Ω v , the interface Γ 4 undergoes remarkable transformations. For example, an air blow over the liquid film is an instance which reflects this situation and is discussed here. In Figure 10 is displayed such a situation where p v K v = 0.1 . It shows how the shapes of the function h t , x describe themselves for subsequent instants of time when subjected to such a constant pressure gradient.
From these, it is observed that they all retain their sinusoidal forms and, furthermore, resolve themselves in a more or less regular manner into a deformed travelling waves, as can be easily seen in Figure 10. It is interesting to see how the interface h t , x evolves when we increase the quantity K v > 0 by 0.1 , the resulting quantity by the same amount, and so forth. In other words, we wish to see what happens when we augment the strength of the air blow. Physically, the following statements do hold. Increasing the quantity p v will both weaken and deviate the diffusive fluxes hitting the interface Γ 4 . Consequently, min h t , x h 0 with a constant pressure gradient amounting to K v = 0.1 would be less than that which would amount to K v = 1.0 , and thence the graphical results of Figure 10. (Note that the placement of the curved arrow appearing there is merely to aid theoretical understanding.)

3.3.2. Effect of an Air Blow with γ 0

We saw in earlier sections that including the term p v > 0 results in translating the contaminated zone to the east, whereas including the term γ results in thinning the liquid film faster. In this section, both terms are operative. Clearly, under this condition, the dynamics of the interface h t , x would be a combined translational movement accompanied by a fast thinning process. Indeed, glancing at Figure 11, one observes that this is clearly the case. However, while in Figure 10 the curved arrow is towards the north-east direction, that of Figure 11 is towards the south-east. We found a simple explanation to explain this mechanism. In fact, if for a prescribed value of p v = K v = 1.0 we additionally take into account the term γ = x h y γ , the resulting effect would be, more diffusive fluxes will strike the interface Γ 4 while undergoing motion. Consequently, when time evolves, that region of the interface Γ 4 which is contaminated will move to the east while thinning the thickness of the liquid film to the south; thence, the curved arrow shown in Figure 11 is towards the south-east.
If the gradient p is disregarded in Equation (14), and the latter expressed in dimensionless form, a γ -dominated h-evolution equation is obtained, whose mathematical description writes
t h + d i v x h 2 2 Ma Ca grad x C l = 0 for   all   t , x I t × Γ 2 ,
This equation is similar to that leading order model proposed by O’Brien [8]; it is identical when C l = C O Brien l , where
C O Brien l t , x = 2 Λ 1 t exp Λ 2 x x 0 2 4 t for   all   t , x I t × Γ 4 .
If variable field h O Brien t , x stands for that solution when Equation (20) is computed using Equation (21), then to validate the h-evolution equation of the present work, it suffices to prove through numerical experiments that Ord h h O Brien , Ord C l C O Brien l ϵ . Indeed, a simple calculation, after neglecting terms of orders h h O Brien k and C l C O Brien l k (with k 2 ), shows that
t h h O Brien + div x h 2 2 Ma Ca grad x C l C O Brien l Ord h h O Brien , C l C O Brien l for   all   t , x I t × Γ 2 .
When those two orders of magnitudes are proved negligible, one says that the solutions of the variable fields in question give agreement with respect to themselves. Consequently, in order to show such agreement, it only necessitates to display by means of numerical experiments graphical results elucidating the relative departure of C l t , x with respect to C O Brien l t , x and of h t , x with respect to h O Brien t , x .
In Figure 12 is portrayed the graph of h h O Brien x , t , and in Figure 13 and Figure 14, respectively, are portrayed the superposition of the graphs C l t , x and C O Brien l t , x , and of h t , x and h O Brien t , x ; these variable field functions are plotted against the substrate coordinate x L , + L , where L = 5 / 2 [1].
Exploring the discussion a step further, one deduces from Figure 12 the following inequalities:
3 ϵ / 10 h h O Brien ϵ for   all   t , x t × Γ 2 .
Figure 15 illustrates the distribution of chemical concentration c l t , x into the liquid film Ω l when the γ = 0 (left-hand side), and when γ 0 (right-hand side).
The Marangoni flow induced in Ω l by the presence of the volatile solvent is illustrated in Figure 16 for the case γ = 0 . In the next section, we conclude the work.

4. Conclusions

We have studied the movement of a liquid film in an endeavor to examine in which situation does the hitherto neglected term γ = x h y γ have an effect, and a constant pressure-gradient-driven flow p v = K v promotes the equilibrium thickness of the liquid film. On trial, the means of several numerical experiments gave graphs which do not have the same shapes. When γ 0 , our results showed that the thinning process produced by the presence of the volatile source is much increased by the inclusion of the term—a phenomenon evidently occasioned in Marangoni drying. Thus, if a liquid film is used to realise a specific functional property relative to its contact surface, asserting that the term γ is negligibly small could in some cases have a profound effect upon that property, and, therefore, can led to inaccurate predictions. When p v = K v is considered, the liquid film thins less rapidly but undergoes translational movement in the direction of the pressure-gradient-driven flow so as to deviate the diffusive fluxes from hitting the free surface of the liquid film horizontally. Although there are applications in coating and drying processes when thin film flow instability resulting from a pressure-gradient-driven flow is undesirable, it is to be noted that the latter is desirable insofar that other transport phenomena come into play as, for instance, when heat transfer is concerned. Furthermore, if the effect of an air blow is too stiff and the thickness of the liquid film quite thin, instabilities result. On the whole, we have expounded conclusively at the following points:
  • The discrepancy between the several results is owing to the effect of the hitherto neglected term γ on the dyamics of the liquid film.
  • The inclusion of a constant pressure-gradient-driven flow p v = K v in the vapor phase domain might be of some advantage in supporting the thinning process.

Author Contributions

Conceptualization, M.I.K. and M.S.; methodology, M.I.K.; software, M.I.K.; validation, M.I.K., M.S. and V.N.; formal analysis, M.I.K.; investigation, M.I.K.; resources, M.I.K.; data curation, M.I.K.; writing—original draft preparation, M.I.K.; writing—review and editing, M.I.K. and M.S.; visualization, M.I.K. and M.S.; supervision, M.S. and V.N.; project administration, M.S.; funding acquisition, M.S.

Funding

This research received no external funding.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A. Modeling Methodology with COMSOL

The modeling methodology is herein explained. For illustration, Figure A1 is considered.
Figure A1. An overview of numerical method. The COMSOL File includes two models, Mod-k (with k = 1 , 2 ). Of these, Mod-1 is a 2 D -Model containing two Physics Interfaces for the scalar fields c j (with j = 1 , 2 ); Mod-2 is a 1 D -Model comprising only one Physics Interface, that computes the h-evolution equation. Additionally, each models includes a coupling operator Opt-k(·) (with k = 1 , 2 ) which passes data from one model.
Figure A1. An overview of numerical method. The COMSOL File includes two models, Mod-k (with k = 1 , 2 ). Of these, Mod-1 is a 2 D -Model containing two Physics Interfaces for the scalar fields c j (with j = 1 , 2 ); Mod-2 is a 1 D -Model comprising only one Physics Interface, that computes the h-evolution equation. Additionally, each models includes a coupling operator Opt-k(·) (with k = 1 , 2 ) which passes data from one model.
Fluids 04 00198 g0a1
Let those models be Mod-1, Mod-2, respectively. With this view in mind the methods of solutions are based on the following sequences of operations. In the first place, add two models to the COMSOL model builder, one two dimensional model (i.e., Mod-1) to solve problem statements S c l , S c v , respectively, and one one dimensional model (i.e., Mod-2) to solve problem statement S h . In the second place, add two Coefficient Form PDE interfaces to the two dimensional model and to the one dimensional model one Coefficient Form PDE interface. To construct the so-called 1 D , 2 D computational geometries, we proceed as follows. Under Mod-1 construct a large rectangular box whose length from x , y = 0 , 0 spans L = 1 to the east and h 0 + L = 2 to the north. To describe the free-surface of the liquid film, a horizontal line defined by the equation y x = h 0 is inserted. (we note that the respective scalings L, h 0 are all dimensionless quantities, but, however, they are here employed to reflect the nature of the thin film.) In like manner, under Mod-2 construct a line segment whose length from x = 0 spans L = 1 to the east. We pointed out that the fields c l t , x , c v t , x and h t , x are coupled. Consequently, to ensure these couplings from Mod-1 to Mod-2, we add a coupling operator Opt-1 to Mod-1 and another coupling operator Opt-2 to Mod-2. Logically, the inclusion of these operators Opt-1, Opt-2 is merely to pass data from a (source) domain to a destination domain. For the operator Opt-1, we let the interface Γ 2 be its source domain and the line segment defined under the model Mod-2 be its destination domain. For the operator Opt-2, we reverse the situation, letting the interface Γ 2 be its destination domain and the line segment defined under the model Mod-2 be its source domain. Thus, to invoke the functions c l t , x , 0 of Mod-1 into Mod-2 and h t , x of Mod-2 into Mod-1, we shall use
Opt 1 c l t , x , h 0 , Opt 2 h t , x ,
respectively. (More precisely, we used two linear extrusion operators, since, we wish to force out (that is, extrude) the source domain of definition of its functions onto an equivalent destination domain; see the COMSOL Reference Guide for further details.) To implement S h , we switch the coefficient form PDE interface of model Mod-2 to a two-dimensional system h , h -system so to speak we set Ξ = h , h T . Note that for all the problem statements S c l , S c v , and S h , the discretization method for computing the solution is quadratic. In Table 1 is exposed all the scalar quantities used, of which, some are partly extracted from the papers of O’Brien, Matar et al. and Buckingham et al. [7,8,16], while others are based on the underlying philosophy of the problem at hand. These are entered in the parameters table, found under the COMSOL global definitions interface. We made no mention of the imposed initial and boundary conditions. Their implementation are immediate, for, it suffices to consider the default and added nodes to implement them.

References

  1. Kim, S.; Kim, J.; Kim, H.-H. Dewetting of Liquid Film via Vapour-Mediated Marangoni Effect. J. Fluid Mech. 2019, 872, 100–114. [Google Scholar] [CrossRef]
  2. Li, C.; Zhao, D.; Wen, J.; Cheng, J.; Lu, X. Evolution of Entrained Water Film Thickness and Dynamics of Marangoni Flow in Marangoni Drying. RSC Adv. 2018, 8, 4995–5004. [Google Scholar] [CrossRef]
  3. Hérnandez-Sánchez, J.H.; Eddi, A.; Snoeijer, J.H. Marangoni Spreading due to a Localized Alcohol Supply on a Thin Water Film. Phys. Fluids 2015, 27, 032003. [Google Scholar] [CrossRef]
  4. Marangoni, C. On the Expansion of Liquid Floating on the Surface of Another Liquid; Tipographia dei Fratelli Fusi: Pavia, Italy, 1865. [Google Scholar]
  5. Thomson, J. On Certain Curious Motions Observable at the Surfaces of Wine and Other Alcoholic Liquors. Philos. Mag. 1855, 10, 330–333. [Google Scholar] [CrossRef]
  6. Leenaars, A.F.M.; Huethorst, J.A.M.; Van Oekel, J.J. Marangoni Drying: A New Extremely Clean Drying Process. Am. Chem. Soc. 1990, 6, 1701–1703. [Google Scholar] [CrossRef]
  7. Matar, O.K.; Craster, R.V. Models for Marangoni Drying. Phys. Fluids 2001, 13, 1869–1883. [Google Scholar] [CrossRef]
  8. O’Brien, S.B.G.M. On Marangoni Drying: Nonlinear Kinematic Waves in a Thin Film. J. Fluid Mech. 1993, 254, 649–670. [Google Scholar] [CrossRef]
  9. O’Brien, S.B.G.M.; Schwartz, L.M. Theory and Modeling of Thin Film Flows. Encycl. Surf. Colloid Sci. 2002, 63, 52–83. [Google Scholar]
  10. Carles, P.; Cazabat, A.M. Spreading of Oil Drops under a Solvent Vapor: Influence of Marangoni Effect. Progr. Colloid Polym. Sci. 1990, 82, 76–81. [Google Scholar]
  11. Carles, P.; Cazabat, A.M. Spreading Involving the Marangoni Effect: Some Preliminary Results. Colloids Surf. 1989, 41, 97–105. [Google Scholar] [CrossRef]
  12. Oron, A.; Davis, S.H.; Bankoff, S.G. Long-Scale Evolution of Thin Films. Rev. Mod. Phys. 1997, 69, 931–980. [Google Scholar] [CrossRef]
  13. Craster, R.V.; Matar, O.K. Dynamics and Stability of Thin Liquid Films. Rev. Mod. Phys. 2009, 81, 1131–1198. [Google Scholar] [CrossRef]
  14. Sultan, E.; Boudaoud, A.; Amar, M.B. Evaporation of a Thin Film: Diffusion of the Vapour and Marangoni Instability. J. Fluid Mech. 2005, 543, 183–202. [Google Scholar] [CrossRef]
  15. Burelbach, J.P.; Bankoff, S.G.; Davis, S.H. Nonlinear Stability of Evaporating-Condensing Liquid Film. J. Fluid Mech. 1988, 195, 463–494. [Google Scholar] [CrossRef]
  16. Bertozzi, A.; Shearer, M.; Buckingham, R. Thin Film Traveling Waves and the Navier Slip Condition. SIAM J. Appl. Math. 2003, 63, 722–744. [Google Scholar]
Figure 1. Geometrical sketch of a thin liquid film ( Ω l ) over a solid substrate ( Γ 2 ) under a line source S x of volatile surfactant. The functions n t , x , t t , x are the unit normal and tangential vectors; y s = h 0 + L s and y L = h 0 + L .
Figure 1. Geometrical sketch of a thin liquid film ( Ω l ) over a solid substrate ( Γ 2 ) under a line source S x of volatile surfactant. The functions n t , x , t t , x are the unit normal and tangential vectors; y s = h 0 + L s and y L = h 0 + L .
Fluids 04 00198 g001
Figure 2. Geometrical interpretations of surface tension and air-blow effects.
Figure 2. Geometrical interpretations of surface tension and air-blow effects.
Fluids 04 00198 g002
Figure 3. Plots of the h-evolution equation in the steady-state regime. The point x 0 = L / 2 , 5 H / 4 is occupied by the source S x . The term p v is set to zero. (Note that the numerical experiments are carried out based on the data given in Table 1.)
Figure 3. Plots of the h-evolution equation in the steady-state regime. The point x 0 = L / 2 , 5 H / 4 is occupied by the source S x . The term p v is set to zero. (Note that the numerical experiments are carried out based on the data given in Table 1.)
Fluids 04 00198 g003
Figure 4. Plots of the function min h t , x with respect to the strength μ s of the source within the range 10 , 40 in the steady-state regime. These plottings are based on a zero pressure gradient basis, i.e., p v = 0 . (Note that the numerical experiments are carried out based on the data given in Table 1.)
Figure 4. Plots of the function min h t , x with respect to the strength μ s of the source within the range 10 , 40 in the steady-state regime. These plottings are based on a zero pressure gradient basis, i.e., p v = 0 . (Note that the numerical experiments are carried out based on the data given in Table 1.)
Fluids 04 00198 g004
Figure 5. Graphs of h t , x versus x-coordinate in the steady-state regime. (Note that for t ranging from 1.75 to 2.5 the results gives h γ 0 = 0 , while, contrarily, h γ = 0 0 within the same range; also, plottings h t , x = 0 (solid line) and h t , x = h 0 (dotted line) are used to delineate the solid substrate and its free surface initially.)
Figure 5. Graphs of h t , x versus x-coordinate in the steady-state regime. (Note that for t ranging from 1.75 to 2.5 the results gives h γ 0 = 0 , while, contrarily, h γ = 0 0 within the same range; also, plottings h t , x = 0 (solid line) and h t , x = h 0 (dotted line) are used to delineate the solid substrate and its free surface initially.)
Fluids 04 00198 g005
Figure 6. Plots of the function min h t , x versus time when the numerical experiments includes (circular markers, ∘) and excludes (square markers, □) the term γ ; solid and dotted lines are approximations based upon the linear expression h t , x = L / 2 = π 1 t + π 2 . (We note that this expression does not hold when t goes to zero, say, when t is within the range 0 , 0.5 , for min h t = 0 , x h 0 , as can be easily seen.)
Figure 6. Plots of the function min h t , x versus time when the numerical experiments includes (circular markers, ∘) and excludes (square markers, □) the term γ ; solid and dotted lines are approximations based upon the linear expression h t , x = L / 2 = π 1 t + π 2 . (We note that this expression does not hold when t goes to zero, say, when t is within the range 0 , 0.5 , for min h t = 0 , x h 0 , as can be easily seen.)
Fluids 04 00198 g006
Figure 7. Graphs of the functions h t , x with respect to the Marangoni number Ma in the steady-rate regime. The term γ is set to zero in one case. (Note the numerical experiment is carried out based on the data given in Table 1.)
Figure 7. Graphs of the functions h t , x with respect to the Marangoni number Ma in the steady-rate regime. The term γ is set to zero in one case. (Note the numerical experiment is carried out based on the data given in Table 1.)
Fluids 04 00198 g007
Figure 8. The vanishing rate of the quantity min h t , x relative to the Capillary numbers Ca = Ca ref ± 0.5 , where Ca ref = 1.0 designates a so-called reference Capillary number.
Figure 8. The vanishing rate of the quantity min h t , x relative to the Capillary numbers Ca = Ca ref ± 0.5 , where Ca ref = 1.0 designates a so-called reference Capillary number.
Fluids 04 00198 g008
Figure 9. The vanishing rate of the quantity min h t , x relative to the Marangoni numbers Ma = Ma ref ± 10 , where Ma ref = 80 designates a so-called reference Marangoni number.
Figure 9. The vanishing rate of the quantity min h t , x relative to the Marangoni numbers Ma = Ma ref ± 10 , where Ma ref = 80 designates a so-called reference Marangoni number.
Fluids 04 00198 g009
Figure 10. Plots of the function h t , x at time t = 1.0 [1]. Here, p v K v = k / 10 , k = 1 , , 10 . The term γ , p v is set to zero. (Note that the numerical experiments are carried out based on the data given in Table 1.)
Figure 10. Plots of the function h t , x at time t = 1.0 [1]. Here, p v K v = k / 10 , k = 1 , , 10 . The term γ , p v is set to zero. (Note that the numerical experiments are carried out based on the data given in Table 1.)
Fluids 04 00198 g010
Figure 11. Plots of h t , x with respect to time, for t = k / 10 such that k = 1 , , 10 ; and p v K v = 1.0 . (Note that the curved arrow indicates the evolution of the h-evolution equation with respect to time; and the numerical experiments are carried out based on the data given in Table 1).
Figure 11. Plots of h t , x with respect to time, for t = k / 10 such that k = 1 , , 10 ; and p v K v = 1.0 . (Note that the curved arrow indicates the evolution of the h-evolution equation with respect to time; and the numerical experiments are carried out based on the data given in Table 1).
Fluids 04 00198 g011
Figure 12. Graph of the thickness difference h h O Brien t , x against the substrate coordinate x L , + L , where L = 5 / 2 [1].
Figure 12. Graph of the thickness difference h h O Brien t , x against the substrate coordinate x L , + L , where L = 5 / 2 [1].
Fluids 04 00198 g012
Figure 13. Graphs of the chemical concentration functions C l t , x and C O Brien l t , x against the substrate coordinate x L , + L , where L = 5 / 2 [1].
Figure 13. Graphs of the chemical concentration functions C l t , x and C O Brien l t , x against the substrate coordinate x L , + L , where L = 5 / 2 [1].
Fluids 04 00198 g013
Figure 14. Graphs of the thickness functions h t , x and h O Brien t , x against the substrate coordinate x L , + L , where L = 5 / 2 [1].
Figure 14. Graphs of the thickness functions h t , x and h O Brien t , x against the substrate coordinate x L , + L , where L = 5 / 2 [1].
Fluids 04 00198 g014
Figure 15. Distribution of chemical concentration c l t , x into the liquid film Ω l . The arrows designate the directions along which they are distributed. (Note that at leading order the computational domain is a rectangular box h 0 × L .)
Figure 15. Distribution of chemical concentration c l t , x into the liquid film Ω l . The arrows designate the directions along which they are distributed. (Note that at leading order the computational domain is a rectangular box h 0 × L .)
Fluids 04 00198 g015
Figure 16. Marangoni flow in the liquid film Ω l at time t = 1.0 (dimensionless unit). The arrows designate the directions of the flow fields. (Note that at leading order the computational domain is a rectangular box h 0 × L .)
Figure 16. Marangoni flow in the liquid film Ω l at time t = 1.0 (dimensionless unit). The arrows designate the directions of the flow fields. (Note that at leading order the computational domain is a rectangular box h 0 × L .)
Fluids 04 00198 g016
Table 1. Data for case studies 1. through 4. Some data are extracted from the papers of O’Brien, Matar et al., and Buckingham et al. [7,8,16], respectively.
Table 1. Data for case studies 1. through 4. Some data are extracted from the papers of O’Brien, Matar et al., and Buckingham et al. [7,8,16], respectively.
ϵ β γ 0 , l μ s K v x 0 CaMa
10 3 10 2 1.810 10 2 k , 5 ϵ L / 4 0.5500
Reference [7]Reference [16]Reference [8]

Share and Cite

MDPI and ACS Style

Khodabocus, M.I.; Sellier, M.; Nock, V. Dynamics of Thin Film Under a Volatile Solvent Source Driven by a Constant Pressure Gradient Flow. Fluids 2019, 4, 198. https://0-doi-org.brum.beds.ac.uk/10.3390/fluids4040198

AMA Style

Khodabocus MI, Sellier M, Nock V. Dynamics of Thin Film Under a Volatile Solvent Source Driven by a Constant Pressure Gradient Flow. Fluids. 2019; 4(4):198. https://0-doi-org.brum.beds.ac.uk/10.3390/fluids4040198

Chicago/Turabian Style

Khodabocus, Mohammad Irshad, Mathieu Sellier, and Volker Nock. 2019. "Dynamics of Thin Film Under a Volatile Solvent Source Driven by a Constant Pressure Gradient Flow" Fluids 4, no. 4: 198. https://0-doi-org.brum.beds.ac.uk/10.3390/fluids4040198

Article Metrics

Back to TopTop