Numerical simulation of non-isothermal thin liquid film flow on inclined plane using an implicit finite difference scheme

The classical problem of the stability and dynamics of thin liquid films on solid surfaces has been studied extensively. Particularly, thin liquid films subjected to various physico-chemical effects such as thermocapillarity, solutal-Marangoni and evaporative instabilities at the film surface has been the focus of research for more than two decades. Various flow configurations of thin film such as thin film on plane, inclined, and wavy surfaces has been the subject of recent investigations. An inclined film compared to a horizontal film, also experiences the gravity force which may significantly influence the nonlinear dynamics of the film coupled with other forces. In this research, we attempt to study the stability and dynamics of thin liquid films subjected to thermocapillarity and evaporative instabilities at the free surface besides instability owing to ubiquitous van der Waals attraction, using numerical simulations. For a Newtonian liquid, flow in thin liquid film on a planar support and bounded by a passive gas, is represented by Navier-Stokes equation, equation of continuity and appropriate boundary conditions. The external effects are incorporated in the body force term of the Navier-Stokes equation. These governing equations are simplified using the so called long-wave approximation to arrive at a nonlinear partial differential equation, henceforth called equation of evolution (EOE), which describes the time evolution of the interfacial instability in the film caused by internal and/or external effects. Efficient numerical method is required for the solution of the equation of evolution (EOE) in order to comprehend the nonlinear dynamics of the thin film. Here we present the results of our numerical simulation using Crank-Nicholson implicit finite difference scheme applied to the thin film model incorporating instabilities owing to gravity, evaporation and thermo-capillarity. Comparison of our results with those obtained from Spectral method, show remarkable agreement for most of the cases investigated.


Introduction
The stability and dynamics of thin liquid films on solid substrates has been studied extensively for the past three decades. Particular attention has been focused on the case of thin liquid films subjected to various physical-chemical effects such as thermocapillarity, solutal Marangoni and evaporative instabilities at the surface. Furthermore, if the film is inclined at an angle to the horizontal, it experiences the gravitational force which might influence the nonlinear dynamics of the film coupled to other forces. Heat transfer in thin liquid film flow on an inclined plane occurs in numerous chemical and biochemical processes such as solid coating, food processing and in systems with evaporative cooling as well as high throughput exchange devices.
In the case of Newtonian fluid flow of thin liquid film on a solid support and bounded by a passive gas, the dynamics are represented by Navier Stokes equations and equation of evolution together with appropriate boundary conditions. The external effects are generally incorporated in the body force term of the Navier Stokes equations. These governing equations can then be simplified using the so called long wave approximation to arrive at fourth order nonlinear partial differential equation, henceforth called the equation of evolution which describes the time evolution of the interfacial instability caused by internal and/or external effects. The details of the derivation of the equation of evolution are available in the literature.
The linear stability characteristics can be obtained by the solution of the linearized equation of evolution. However, for the accurate simulations of the nonlinear dynamics and surface morphology of thin-film, we require efficient numerical methods for the solution of the equation of evolution (EOE). The extent of nonlinearity and the stiffness of the resulting differential equation depend upon the nature of various physico-chemical effects incorporated in the thin-film model [1][2][3][4][5][6][7][8][9][10][11]. Numerous works on numerical simulations of various liquid thin-film models have been reported in the literature [e.g., 12,13,14,15]. Burelbach et al. (1988) [14] has solved the problem of evaporating/condensing liquid film on horizontal plane using an implicit finite difference scheme, Crank-Nicholson mid-point rule, by employing centered difference for space and forward difference for time derivatives. The resultant nonlinear difference equations are solved by Newton-Raphson iterations. Ali et al. (2005) [13] have employed similar discretization scheme for their isothermal liquid thin-film on inclined plane, however IMSL subroutines were employed for the solution of resulting difference equations which essentially uses an adoption of Newton-Raphson method. Fourier spectral method was used by Joo et al. (1991) [12] for numerical simulations of heated falling films under the influence of thermocapillarity and evaporative effects at the surface, however, they have ignored van der Waals effects in their formulation. There are many other works that have been reported in the literature.
The present paper illustrates some results from the numerical simulation of the flow of the thin liquid layer on inclined plane and subjected to evaporative and thermocapillarity instabilities. An implicit time averaged Crank Nicholson finite difference scheme formulation of the EOE has been used for numerical simulations. To validate our numerical codes of our implicit Crank Nicholson finite difference scheme we apply them to the solution of the equation of evolution of long-wave instabilities of heated falling film discussed by Joo et al. (1991) [12]. The results generated using our code are compared with the results reported by Joo et al. (1991) [12] which uses Fourier spectral method. The results from our implicit finite difference code closely match with the Fourier spectral thereby proving the viability of the implicit scheme as a reliable alternative to an explicit finite difference method. The validated solution procedure has thereafter been employed to study the time evolution of surface instabilities in a non-isothermal thinfilm by incorporating the thermocapillarity and evaporation effects.

Thin film model
We assume a thin viscous Newtonian liquid film. The film is bounded above by its vapour and below by a uniformly heated rigid plane inclined at an angle θ to the horizontal. The film is sufficiently thick so that a continuum theory of the liquid is valid. The liquid film consists of an incompressible Newtonian fluid with constant material properties and is laterally unbounded. Evaporation occurs at the layer resulting in mass loss, momentum transfer and energy consumption. Plate has a fixed constant temperature . Liquid temperature on the free surface is controlled by losses to the passive gas above. The surface tension σₒ depends on temperature so that the thermocapillarity effects are present.
We consider a 2-D Cartesian coordinate system as shown in figure 1. The vapor=liquid interface is located at z=h(x,t). Liquid is assumed to evaporate in a direction normal to the interface. The Navier-Stokes and continuity equations in a non-dimensional Cartesian co-ordinate system (x, z) defined as in figure 1 are: : Where, the spatial coordinates are scaled by mean film thickness ℎ 0 ; time by viscous time scale ℎ 0 2 . u and w are the components of velocity in x-and z-directions, respectively and p is the pressure. The subscripts denote the derivatives.

Equation of Evolution
The length scale of the initial instability in a thin film is much larger than the mean film thickness ℎ 0 . For flows with a characteristic length in the x-direction typically proportional to the disturbance wavelength and much larger than the thickness of the liquid layer, the long wave analysis of Benney (1966) [15] can be applied to the problem at hand.

Numerical method
We have employed an implicit finite difference numerical discretization scheme discussed below.

Implicit Finite Difference Scheme-Crank Nicholson Mid Point Rule
We use the conservative form of the rescaled equation of evolution. This equation of evolution is discretized using forward difference in time and central difference in space with time averaging, resulting into a set of difference equations for the film thickness h. These equations are then solved as a set of nonlinear coupled algebraic equations by an iterative procedure of the type of Newton's method. Let subscript i denotes space discretization and superscripts n and n + 1 denote values of variables at nth and (n + 1)th times. Thus, Eq. (5) after discretization appears as:  .

Conclusion
The results obtained by our implicit Crank Nicholson finite difference method demonstrate a close conformity to the results obtained by Fourier spectral method of Joo et. al. (1991), with some slight deviation in some cases and significant deviation in one of the cases investigated. This indicates that the implicit finite difference method is a reliable method for investigating the highly nonlinear thin film model to a high degree of accuracy.