Abstract.
In this article, we present a model of heat transfer occurring through a liquid film flowing down a vertical wall. This new model is formally derived using the method of asymptotic expansions by introducing appropriately chosen dimensionless variables. In our study the small parameter, known as the film parameter, is chosen as the ratio of the flow depth to the characteristic wavelength. A new Nusselt solution should be explained, taking into account the hydrodynamic free surface variations and the contributions of the higher order terms coming from temperature variation effects. Comparisons are made with numerical solutions of the full Fourier equations in a steady state frame. The flow and heat transfer are coupled through Marangoni and temperature dependent viscosity effects. Even if these effects have been considered separately before, here a fully coupled model is proposed. Another novelty consists in the asymptotic approach in contrast to the weighted residual approach which have been formerly applied to these problems.
Key words and phrases: Heat transfer; Thin liquid film; Asymptotic modelling; Long waves; Thermal dependency properties; Marangoni effect
MSC:
2010 Mathematics Subject Classification:
76D33 (primary), 76B25, 76B15 (secondary)[]
[]
[]
[]
[]
\contentsmargin0.5em
\titlecontentssection[]
\contentslabel[\thecontentslabel]
\thecontentspage
[]
\titlecontentssubsection[]
\contentslabel[\thecontentslabel]
\thecontentspage
[]
\titlecontents*subsubsection[]
\thecontentspage
[
]
Marx Chhay
LOCIE, Polytech Annecy–Chambéry, France
Denys Dutykh
CNRS, Université Savoie Mont Blanc, France
Marguerite Gisclon
CNRS, Université Savoie Mont Blanc, France
Christian RuyerQuil
LOCIE, Polytech Annecy–Chambéry, France
New asymptotic heat transfer model in thin liquid films
arXiv.org / hal
Contents
1. Introduction
Liquid films have many significant applications in chemical engineering because of their reduced resistance to heat and mass transfers. This thermal resistance may be further reduced by inducing surface deformations and wave patterns as surface instabilities may be triggered by the dependence of surface tension to temperature (Marangoni effect [Joo1991]) or by an hydrodynamic mechanism when the fluid is set into motion either by gravity (Kapitza instability mode [Kalliadasis2012]) or by centrifugal acceleration [Matar2005]. As an example, the Marangoni effect may be coupled to wall topography to induce thermocapillary convection [Alexeev2005] and to promote heat transfer.
In this paper, we focus on falling liquid films. These flows are generally encountered whenever the pressure drop is critical, e.g. in absorption machines, or whenever a low thermal driving force is required, for instance in the separation of multicomponent mixtures that are temperaturedependent. The dynamics of such flows have attracted a considerable interest as it presents a wavy regime organized around largeamplitude teardrop like solitary waves whose interactions intensify transfers. This wavy regime is triggered by a longwave instability mode corresponding to a zero critical wavenumber. For this reason, the waves are long compared to the film thickness, they emerge at relatively long distances from the liquid inlet and they are slow to interact one with another. As a result, direct numerical simulations (DNSs) of such flows are hindered by the large domain that is necessary to account for their natural evolution, which explains that DNSs are generally restricted to twodimensional, i.e. spanwise independent, situations or to the construction of periodic waves.
Mathematical modeling offers a useful reduction of the numerical cost and a welcome framework for the understanding of the disordered dynamics of such flows with the development of coherentstructure theories. Indeed, the large aspect ratio of the waves enables to introduce a small parameter , or film parameter, which compares the typical length of the wave to the thickness of the film. In this framework, the streamwise () and spanwise () coordinates as well as the time () are slow variables, i.e. , whereas the crossstream coordinate is a fast variable (). It is thus possible to eliminate the fast variable and to obtain a reduced set of equations which describes the slow evolution of the film in a spatial domain whose dimension is reduced from 3D to 2D or from 2D to 1D if spanwise independent solutions are looked after. Following Kapitza’s initial work [Kapitza1948], an important amount of work has been produced in order to derive such reduced set of equations or lowdimensional models (see for instance the review by [Kalliadasis2012]). Benney [Benney1966] thus showed that a series expansion of the flow variables with respect to the film parameter leads to a solution that is fully characterized by the film thickness and its gradients, the film dynamics being governed by a single evolution equation for . Unfortunately, Benney’s equation admits nonphysical singularities in finite time at moderate Reynolds number [Pumir1983] as a result of a too strict slaving of the velocity field to the free surface elevation. A cure to this shortcoming is offered within the SaintVenant framework after averaging the primitive equation across the film depth. This idea dates back to the original work of Kapitza [Kapitza1948] and was successfully applied by Shkadov [Shkadov1967] who derived a set of two evolution equations for the local thickness and the local flow rate . Shkadov’s averaging approach requires a closure hypothesis in the form of a polynomial ansatz for the velocity distribution, which corresponds to the Nusselt parabolic profile in Shkadov’s classical work. More sophisticated ansatz have been proposed e.g. in [Yu1995, Bohr1997]. However, consistent averaging of the primitive equations has been introduced by Roberts [Roberts1996] and RuyerQuil & Manneville [RuyerQuil2000] using different approaches. Currently, one of the widely used approaches is the Weighted Residual Method (WRM), which has been successfully applied to derive approximate equations [RuyerQuil2000, Kalliadasis2012]. However, the main drawback of the WRM is that the mathematical structure of the resulting averaged equations is unclear. Consequently, the models derived in this way may be difficult to justify mathematically. This is the main reason why we opt for the method of asymptotic expansions in the present study. Other references on this topic include [Luchini2010, Rojas2010, Benilov2014]. Note that the approach by Luchini and Charru [Luchini2010], based on the averaging of the kinetic energy equation of the flow, leads to a result that is identical to the WRM method at firstorder of the film parameter. In essence, the Inertial Lubrication Theory introduced by [Rojas2010], is similar to the WRM method, leading to very similar equations, coefficients being almost identical for most terms. Benilov [Benilov2014] derived depthaveraged equations using the Galerkin method, which is one particular weighted residual method.
Benney’s original work has been extended in [Joo1991, Joo1991a] to deal with the conduction of heat across the film and the coupling of the hydrodynamics to the transfer offered by the dependence of surface tension on temperature (Marangoni effects). The classical Benney expansion followed for instance by [GambaryanRoisman2010, Leslie2011, Leslie2012] to deal with the heat transfer across an horizontal film is not appropriate to deal with a falling film. Indeed, this approach requires that heat diffusion overcomes convection, and is thus limited to only small values of the Péclet number. To overcome this limitation, Scheid et al. followed the Weighted Residual technique initiated by RuyerQuil & Manneville [RuyerQuil2000] and derived several models of various accuracy [RuyerQuil2005, Scheid2005]. Though enabling to accurately decipher the complex interplay between the Kapitza hydrodynamic instability and the longwave Marangoni thermocapillary instability [Scheid2008], these models are only valid at still relatively low values of the Péclet number. Indeed, as the Péclet number is raised, these models may predict nonphysical values of the temperature. This behaviour is related to the onset of sharp temperature gradients at the free surface due to flow orientation in largeamplitude solitary waves. Though a cure to this limitation has been proposed by [Trevelyan2007], available lowdimensional models still fail to capture correctly the temperature distribution at large Péclet number.
This article aims at deriving a new conservative formulation for heated falling film wich retains consistency with the classical longwave expansion. This new formulation can be seen as a lowdimensional modelling of heated falling film flows following a derivation procedure that has been proposed by Vila and coworkers [Boutounet2008, FernandezNieto2010, Boutounet2013, Noble2007]. This procedure, which will be referred to hereinafter as the SaintVenant consistent approach, is based on the classical SaintVenant equations that are obtained by indepth averaging of the primitive equations with a uniform weight. However, contrary to the Kapitza–Shkadov approach which assumes the velocity field to be strictly parabolic [Shkadov1967], Vila proposes a closure that is compatible with Benney’s longwave asymptotics and enables to accurately recover the threshold of the Kapitza instability. Our derivation takes into account the thermocapillary effect and the dependence of viscosity with respect to temperature. These are the principal coupling mechanisms between the hydrodynamics of the film and the heat transfer induced by the dependence of the thermophysical properties of the fluid with temperature [DAlessio2014, Alexeev2005, GambaryanRoisman2010].
The structure of the paper is as follows. In the next Section, governing equations are recalled. Then some physical behaviour of the heated falling film are highlighted, using a basic modeling in order to introduce the dynamic of the system. In Section 3 the asymptotic model is derived. In Section 4, some numerical experiments illustrate in one hand the realistic behaviour of the model and, on the other hand, a comparison with an existing model in the literature is performed. Some discussions about the formal derivation conclude this work in Section 5.
2. Problem formulation
2.1. Governing equations
We consider an anisotherm liquid film flowing down a heated vertical plate. The flow is supposed to be twodimensional in space, the axis corresponding to the streamwise direction and the axis to the crossstream direction. The liquid domain is delimited by a vertical wall at and the free surface boundary located at the height . The sketch of the fluid domain is depicted in Figure 1.
The motion of the liquid is governed by the incompressible Navier–Stokes equations
where , and represent the velocity and pressure fields and the gravity acceleration vector. The physical parameters , , correspond to the density, the kinematic viscosity and the dynamic viscosity. The heat transfer and the hydrodynamics of the film are coupled by the dependence of the physical properties with respect to the temperature. In this paper we focus on the principal sources of coupling and assume that only the surface tension and the dynamic viscosity of the liquid depend on the temperature . As a further simplification, we assume linear laws,
and
where and are positive constants as surface tension and viscosity generally decrease with the temperature.
The heat transfer occurring through the liquid domain is modelled by the advectiondiffusion Fourier equation
where corresponds to the temperature field and is the thermal diffusion coefficient, which is assumed to remain constant.
At the wall, a noslip condition
(2.1) 
and a constant wall temperature are imposed while at the free surface, the kinematic condition governing the evolution of the fluid elevation reads
(2.2) 
The continuity of the fluid stresses at the free surface gives
where .
Heat transfer at the free surface is modelled by a thermal exchange coefficient that is assumed to remain constant so that temperature field verifies the Newton’s law of cooling at the free surface
where and denote the thermal conductivity and the unit exterior normal
2.2. Scaled equations
The specific geometry of the falling film is characterized by the typical length scales in both the streamwise direction and the crossstream direction. The evolution of the hydrodynamic instabilities and the thermal diffusion process can also be described through these typical lengths. Introducing the dimensional quantities

: streamwise typical length scale,

: crossstream typical length scale,

: typical average velocity corresponding to hydrodynamic Nusselt solution
and the following change of variables [Noble2007]:
six dimensionless numbers characterize the problem at hand:

the Reynolds number ,

the Péclet number ,

the Biot number ,

the Weber number ,

the Marangoni number ,

the dimensionless rate of change of the dynamic viscosity
besides the film parameter . In what follows, the reference temperature is further chosen to correspond to the temperature of the air, i.e. . For most liquids, surface tension is high and the Weber number is typically large. We thus introduce with . This classical assumption (see [Kalliadasis2012]) enables to include surface tension effects at an early stage of the asymptotic longwave assumption, as surface tension is the sole physical effect which prevent the breaking of the waves. The coupling of the flow dynamics with respect to the heat transfer is accounted for by the Marangoni number and the dimensionless rate of change of the viscosity that are assumed to be order one quantities. Note that the positivity of viscosity and surface tension demands that and .
This set of parameters is usefully completed with the Kapitza number and , where is the capillary length, and is a viscousgravity length [Kalliadasis2012]. The dimensionless groups and are independent of the film thickness and depend only on the fluid properties. The thin liquid depth is characterized by the ratio i.e. .
The dimensionless incompressible Fourier–Navier–Stokes equations read
(2.3a)  
(2.3b)  
(2.3c)  
(2.3d)  
For convenience, the over bar notation for the dimensionless quantities have been dropped in the above equations. The noslip boundary condition at the wall (2.1) and the kinematic condition at the free surface (2.2) remain formally unmodified. The continuity conditions of the fluid stress across the free surface become  
(2.3e)  
and  
(2.3f)  
where the total derivative is given by  
The dimensionless heat transfer between the heated liquid and the ambient air becomes  
(2.3g)  
whereas the Dirichlettype boundary condition at the wall is  
(2.3h) 
2.3. Fourier full 2D model
In order to illustrate numerically the heat transfer behaviour depending on relevant physical parameters, we present below some numerical simulations of the Fourier equation, the solution to the Navier–Stokes equation being approximated by the lowdimensional model that is presented below.
The numerical solution is looked after in a stationary rectangular domain thanks to the change of variables
The transformed heat field becomes a solution of
(2.4) 
where the differential operators are given by
and corresponds to the velocity vector field shifted by the wave celerity .
The boundary condition at the wall remains
and the Robintype condition at the free surface becomes
The velocity field is computed using the parabolic approximation for the downstream component with and . The crossstream velocity component is computed such that the incompressibility is verified.
The numerical results for the steady temperature profile are obtained using an implicit second order scheme. The isothermal lines plotted in Figure 2 have been computed from equation (2.4) for an analogous configuration as is [Trevelyan2007]. The solitary wave profile and the velocity field under the wave were obtained from Vila’s model [Boutounet2008]. The Reynolds number is fixed at , for various fluid properties ( and ).
When no heat transfer is allowed between the liquid film and the surrounding air, the falling film reaches the uniform temperature given by the heated wall (). When heat exchange between the two fluids phases is maximal (), the liquid film behaves as a conductive medium between the heated wall and the colder air.
The temperature follows a linear distribution across a flat film (Nusselt solution). However, as soon as hydrodynamic instabilities occur, a recirculation zone within largeamplitude waves may appear when the fields are described from the wave moving frame. The heat transfer through the liquid film is locally far from being linear. This corresponds to the physical mechanism of heat enhancement, as used in engineering process. In case of vertical falling film, Benjamin [Benjamin1957] has shown the appearance of such inertial hydrodynamic instabilities. Therefore, when considering the vertical configuration of anisotherm falling film, the intensification of heat transfer by the hydrodynamic instabilities must be taken into account. The hypothesis of flat falling film does not stand anymore.
3. Asymptotic model
In what follows, we present the derivation of a system of averaged equations for the mass, momentum and heat balances (2.3a) – (2.3d) following the classical SaintVenant approach. We consider a formal expansion for the velocity, the pressure and the temperature field with respect to the order parameter :
(3.1) 
Consistency with the longwave asymptotic expansion of the flow variables with respect to the parameter is guaranteed by computing the higher order corrections to the thermal Nusselt solution, corresponding to the leading order term of the asymptotic development for the temperature field.
3.1. Thermal Nusselt solution
The substitution of the formal development (3.1) into the dimensionless Fourier equation (2.3d) yields
by taking the limit . In the following, we set
(3.2) 
Taking into account the associated boundary conditions
the resolution of the second order differential equation gives the expected expression for the main order temperature profile called the thermal Nusselt solution [Burelbach1988]:
(3.3) 
with
(3.4) 
We can remark that the Nusselt thermal gradient is function of and its first derivative with respect to .
Although its computation is straightforward, the expression of the thermal Nusselt solution differs from the linear temperature profile found in the literature [DAlessio2010, RuyerQuil2008]. The thermal gradient involves the liquid film height depending on the . The proposed thermal Nusselt solution gives with more accuracy the influence of the hydrodynamic instabilities than as if its linear factor would just have been constant along the downward direction. It is worth to point out that, even when , travelling waves may appear as solutions of the Shallow Water model. Thus, even in this limit case (), the falling film profile may vary along the downward direction and it is expected for the thermal Nusselt solution to describe this behaviour as well.
Remark that the average Nusselt solution is consistent with the limit cases of heat transfer. Indeed out of the instability neighborhood, the uniform temperature profile is reached (adiabatic case: , ) and whereas no resistance to the transfer at the free surface happens, , the linear heat profile is realized yielding .
3.2. Formal derivation of the velocity
We consider the flatfilm velocity distribution, i.e. . Integration of (2.3a) using the boundary conditions (2.1), (2.2) and (2.3f) gives
The graphical comparison of different profiles is shown in Figure 3 for an already large value of the rate of change of viscosity and a thermal gradient . The velocity of the flow increases with as expected from the lowering of the viscosity and the subsequent reduction of the viscous stresses. A significant departure of the velocity profile from the Nusselt parabolic distribution () can be observed. However, the polynomial velocity distribution obtained from a Taylor expansion is still very close to the exact distribution which involves a logarithmic correction. This suggests that the assumption holds up to already large values of . Since the logarithmic correction to the velocity profile will significantly complicate the algebra and the resulting models, we will hereinafter assume only small values of the rate of change of viscosity in our derivation process with the expectation that the resulting simplified equations will remain valid up to order one values of .
We thus assume and introduce a constant such that
(3.5a)  
Using the free divergence condition and the no slip condition , the transverse component velocity is determined as  
(3.5b)  
3.3. Formal derivation of the pressure
Truncated at , integration of the crosstream momentum balance (2.3b) yields
The boundary condition (2.3e) gives
Where the dependence of the pressure with respect to the surface tension has been made explicit even though it formally appears as a correction. As underlined above, inclusion of surface tension effects is required to capture the onset of capillary waves at the front of the solitary waves and to prevent wave breaking, which justifies the assumption of large Weber numbers . As a consequence, the pressure distribution at leading order reduces to the sole contribution of surface tension
(3.6) 
3.4. Consistent averaged momentum equation
In this Section, we present the derivation of a consistent averaged momentum equation following [Boutounet2013]. The proposed momentum balance includes the coupling with heat transfer resulting from the dependence of surface tension and viscosity with respect to the temperature.
We first integrate the streamwise momemtum balance (2.3a) using the boundary conditions (2.1) and (2.2), to obtain an averaged momentum balance
(3.7) 
We know that with and
where .
We introduce the depthaveraged velocity
(3.8) 
so that
(3.9) 
Equation (2.2) gives
(3.10) 
The boundary condition (2.3f) gives
(3.11) 
whereas the wall shear stress is evaluated from the asymptotic expansion which leads to
(3.12) 
The chosen expressions (3.9) and (3.12) of the momentum flux and of the wall shear stress are not unique as one can play with the asymptotic expression of the averaged speed (3.8). Equation (3.9) introduces the classical momentum flux of the shallowwater equations corrected with a term that is function only of the film thickness. Following [Boutounet2013], where the analogy of the shallowwater equations with the compressible Euler equations is underlined, the film thickness being the analogue of the density, the correction is a barotropic contribution to the pressure. The wall shear stress expression (3.12) is chosen so that the dependency of the viscosity with respect to the temperature appears as a correction to which corresponds to the classical Nusselt parabolic profile.
The asymptotic expression of the pressure distribution (3.6) gives
(3.13) 
From (3.9), (3.11), (3.12) and (3.13), the averaged momentum balance (3.7) reads
(3.14) 
The asymptotic momentum balance (3.14) extends the result obtained by [Boutounet2013] by including the principal coupling terms with the heat transfer, i.e. a reduced wall friction and a driving stress proportional to the gradient of free surface temperature, which accounts for the Marangoni effect.
3.5. Asymptotic heat transfer model
Before turning to the averaging of the heat equation (2.3d), we first compute the secondorder correction to the linear distribution given by (3.3). The substitution of , and gives
which can be easily integrated with the help of the boundary conditions
and the expressions (3.5) of and . We thus obtain
3.6. The conservative formulation of the model
In this section, we derive formally an evolution equation for the freesurface temperature temperature following the consistent SaintVenant approach.
We introduce the free surface temperature
(3.15) 
so that the temperature field may be written
where the correction verifies as required by the definition (3.15) so that .
Integration of the the heat equation (2.3d) across the film height yields
(3.16) 
which involves the average temperature and the mixing temperature . The average temperature can be easily derived from the freesurface temperature through the linear temperature distribution (3.3) by . An expression of the mixing temperature in terms of the variables , and requires a closure which can be provided again from the Nusselt solution . Yet, this closure is not unique as one can play with the asymptotic expressions of the freesurface temperature and flow rate
(3.17a)  
(3.17b)  
The point of view that is adopted below is to express the convective terms at the l.h.s. of (3.16) as a function of . This choice stems from i) the classical form of the shallow water equation and the choice of , anf as the variables of the model, ii) the fact that the hyperbolic structure of the obtained evolution equations is then guaranteed as will be underlined in Section 4. The r.h.s. of (3.16) is written as a function of using