Rivulet formulation in the flow of film down a uniformly heated vertical substrate

ABSTRACT Three-dimensional regularized weighted residuals model is developed to study vertical falling film on a uniformly substrate with constant heat flux. The model is verified through linear stability analyses in both the stream and the span directions. It is shown that rivulet formation is due to the dominance of the thermocapillary Marangoni effect in the span direction, i.e., the component in the span direction of the surface temperature gradient is much larger than the component in the stream direction, which results from the different heat transfer modes in the stream and the span directions. Four types of initial disturbance are separately imposed on the film flow. It is shown that the initial disturbance does have an influence on the rivulet formation process, but it has no effect on the eventual film structure. Rivulet formation on the film enhances the heat transfer on a wavy film surface, which is attributed for the most part to the increase in the surface area of the film. On the effect of the Reynolds and the Marangoni numbers on rivulet formation, an increase in the former results in an increase in rivulet width, while an increase in the latter causes an advance in the film rupture.


Configuration of film flow
A surface area of film exposed in ambient air A w substrate area spread with film E x deformation energy in the stream direction E z deformation energy in the span direction E HT enhancement of heat transfer on wavy film surface CONTACT Riqiang Duan duanrq@mail.tsinghua.edu.cn k wavenumber p a ambient pressure Q loss heat loss from substrate to ambient air ( = α w (T w − T a )) Q w heat flux distributed in substrate T surface temperature of wavy film T a ambient temperature T N surface temperature of Nusselt flat film T average surface temperature of wavy film T w temperature in substrate α heat transfer coefficient between film and ambient air α w heat transfer coefficient between substrate and ambient air β inclination of substrate T temperature difference related with heat flux Q w θ dimensionless average surface temperature of wavy film scaled with T N . θ N dimensionless surface temperature of Nusselt flat film scaled with T N rivulet-to-rivulet spacing σ a surface tension at T a ω complex angular frequency Viscous-gravity scaling l υ viscous-gravity length scale(= (υ 2 /gsinβ) 1/3 ) t υ viscous-gravity time scale ( = l υ /t υ ) u υ viscous-gravity velocity scale (= (υgsinβ) 1/3 )

Nusselt flat film
h N Nusselt film thickness h N modified Nusselt film thickness (=h N /l υ ) q N Nusselt flat film flow ratē u N Nusselt average velocity

Free variables
h local film thickness p pressure of film q local flow rate in stream direction q ⊥ local flow rate in span direction θ dimensionless surface temperature scaled with T N .

Introduction
Film falling down a heated substrate exhibits intuitively some typical thermohydrodynamic behaviors detected by easy visual observation, i.e. the transition from laminar flow to turbulent flow or the occurrence of film rupture and dry-out. Meanwhile, this phenomenon is widely applied and occurs in many industrial and natural fields, e.g. it is a key process for passive containment cooling systems in the AP1000 nuclear reactor and in multipleeffect distillation (MED) for sea water desalination due to its high efficiency in heat and mass transfer and low resistance in flow. For a thin, heated film flowing down an incline, the film is subject to the formation of a rivulet pattern aligned with the mean flow that may progress to film rupture under the interaction between hydrodynamic waves and longwave thermocapillary instability (Goussis & Kelly, 1991). Rivulet formation in film flow has two effects. The advantageous effect is that rivulet formation largely increases the effective heat and mass transfer. This effect partly results from the increase in the effective area of heat and mass transfer on the film surface due to its deformation and partly from the enhancement of heat and mass transfer due to a local increase in the Reynolds number in the rivulet and a decrease in the film thickness in the furrow between the arrays of the rivulet. The disadvantageous effect is that furrows between the arrays of the rivulet are susceptible to dry-out, which greatly lowers the device efficiency and leads to failure of equipment construction integrity.
The pioneering study on rivulet structure in film flow was carried out by Ludviksson and Lightfood (1968) who set up a film on a vertical substrate heated with a constant downward directed temperature gradient to study the thermohydrodynamic instability of a Marangoni film. Various numerical models simulating thin-film flow have been built to study the evolution of instability of falling film. There are several representative models such as the Benney equation (hereinafter abbrev. BE) (Benney, 1966), the lubrication approximation model, the regularized weighted residuals model (hereinafter abbrev. RWRM) (Ruyer-Quil,  and the integration method of the three-dimensional conservation equations for mass, momentum and energy (Ramaswamy, Krishnamoorthy, & Joo, 1997). Joo, Davis, and Bankoff (1996) have modeled rivulet structure using the BE model. However, their computer simulations experienced finite-time blowup, which is known to be an intrinsic drawback of the BE model. Holland et al. (Holland, Duffy, & Wilson, 2001) used the lubrication approximation model to undertake numerical analyses of film flow down a uniformly heated-cooled substrate and showed that a rivulet is formed due to a transverse flow driven by the thermocapillary Marangoni effect. Ramaswamy et al. (1997) integrated the full Navier-Stokes and Fourier equations in three dimensions to simulate a heated falling film flow down a substrate maintained at a constant temperature and analyzed the nonlinear dynamics of rivulet formation from onset to rupture. However, their computations were restricted to a small domain size and a small Reynolds number. Scheid, Kalliadasis, Ruyer-Quil, and Colinet (2008) used the RWRM method to explore the three-dimensional heated film flow under constant temperature conditions. They discovered that a rivulet pattern is formed and aligned with the mean flow at small Reynolds numbers. The RWRM method, based on a weighted residuals model, is applicable for small to moderate Reynolds numbers (Liu, Jiang, & Duan, 2017;Ruyer-Quil et al., 2005). Compared with other methods, the RWRM model is simpler than the three-dimensional integration method of the conservation equations and is more accurate than either the BE or lubrication approximation models at moderate Reynolds numbers (Kalliadasis, Ruyer-Quil, Scheid, & Velarde, 2012).
All the abovementioned previous studies have primarily focused on an isothermal falling film or a falling film that is heated at a specified temperature boundary condition (hereinafter abbrev. ST) or with local heating. In this paper, a falling film heated using a uniform heat flux boundary condition (hereinafter abbrev. HF) will be investigated. The HF condition is more realistic than the ST condition because it is much closer to actual situations and thus has very good prospects for practical applications. Compared with the ST thermal boundary condition, the HF boundary condition is relatively complicated since both heat transfer from a heated substrate to the film and heat loss from the substrate to the ambient air need to be taken into account. Therefore, the first objective of this paper is to fill the gaps in the formulation of the RWRM method of film flow with HF thermal boundary conditions.
The second objective of this paper is to elucidate the mechanism of rivulet formation in film flow, wherein two issues are concerned. One issue is the role of the thermocapillary Marangoni effect in film flow rivulet formation. With regard to rivulet formation, it is definitely ascribed to the thermocapillary Marangoni effect due to the temperature heterogeneity on the film surface. In fact, the thermocapillary Marangoni effect has no preferred direction. Goussis and Kelly (1991) thought the emergence of a preferred direction of the rivulet was entirely due to the presence of shear flow. We think that this is still an open question that needs to be further investigated. Another issue is the effect of initial disturbance on the characteristics of the formed rivulet array, i.e. the pattern and density of the rivulet. For isothermal film flow, as is well known, the linear development stage of film flow is a low-pass filter in which the disturbance with the highest growth rate is dominant and other modes are diminished. For heated film flow, does the interaction between hydrodynamic waves and the thermocapillary Marangoni effect change this situation, especially in the stream direction?
The third objective of this paper is to explore the mechanism of enhancement of heat transfer of film flow under film instability. Two questions are answered; one question is how much enhancement in heat transfer is attributed under film instability, and the other is the rupture time of the film, which is of special interest to the engineering design of heat exchangers.

Mathematical formulation of the problem
The configuration of film flow along an inclined flat substrate is shown in Figure 1. The angle of inclination between the substrate and the horizon is β. Film flow is uniformly heated by a heat source embedded in the substrate with a constant heat flux. Considering practical applications, it is assumed that some heat, Q loss , is leaked into the ambient air on the opposite side of the substrate, and the heat transfer coefficient between the substrate wall and the ambient air is α w . The rest of the heat, Q w , is absorbed by the film flow and eventually released into the ambient air as well, and the heat transfer coefficient on the free film surface is α.
The heat flux imposed on the film flow is correlated with the Marangoni number, which measures the surface gradient of the surface tension resulting from the heterogeneous distribution of the surface temperature. For the heat flux boundary conditions on the substrate wall, the temperature difference scale, T, is related to the heat flux imposed on the film flow, as shown in Eq. (1) The surface tension of the film liquid is σ = σ a −γ (T−T a ) where σ a is the surface tension at the ambient temperature T a , and γ = dσ dT is the temperature dependence of the surface tension of the film liquid.
The governing equations for a uniformly heated film flow are the continuity equation, the Navier-Stokes equation and the energy equation with appropriate boundary conditions. Accounting for the large film thickness to characteristic length scale ratio in both the stream and the span directions, the mathematical formulation in this paper is simplified using the gradient expansion strategy (Kalliadasis et al., 2012) in which it is postulated that all hydro-and thermodynamic qualities vary slowly in the stream direction, the span direction and time in such a way as is observed in the longwave expansion. An ordering parameter, , is introduced to estimate the amplitude order of each term in the controlling equations and the corresponding boundary conditions through the following transformations.
is also called the 'film parameter' and is a small number that compares the film thickness to the characteristic length scale in the stream direction.
The governing equations and their corresponding boundary conditions are nondimensionalized based on the Shkadov scaling in which the coordinates in the stream and the span directions are compressed by taking their scale as 1/ε times the scale for the coordinate in the cross-stream direction (Shkadov, 1977). The details are as follows.
where variables with primes are dimensionless. For simplicity, primes are dropped in the following. Compared with both the Nusselt scaling and the viscous-gravity scaling (Kalliadasis et al., 2012), the advantage of the application of the Shkadov scaling is that it locates the transition between the drag-gravity and drag-inertial regimes at δ 1. Moreover, an additional technical advantage of this scaling is that its parameters retain values close to unity in the region for moderate Reynolds numbers and large Weber numbers, which is rather useful for the convergence of numerical schemes (Kalliadasis et al., 2012). By algebraic manipulation, the dimensionless versions of the governing equations and the corresponding boundary conditions up to second-order approximation are derived as follows.
3εRe(∂ t u + u∂ x u + v∂ y u + w∂ z u) 3εRe(∂ t w + u∂ x w + v∂ y w + w∂ z w) The dimensionless boundary conditions on the substrate wall at y = 0 are given as follows and at the free surface, y = h(x, z, t), the kinematic condition is given by The dynamic boundary conditions at the free surface are given by Due to the thinness of the film, the momentum equation in the cross-stream direction is simplified by integrating it across the film, as is done in the boundary layer approximation. One then obtains then, substituting p into Eq. (8) and (10), yields where the velocity component, v, in the cross-stream direction is calculated through the continuity equation by integrating in the cross-stream direction In the following, the Buckingham π-theorem is used to make a similarity analysis. Film flow is governed by the film thickness, velocity distribution and surface temperature, which are determined by the physical properties of the fluid, the configuration of the substrate, the hydroand thermodynamic qualities, and the disturbance waves on the free surface. Their dependence relation can be written as follows.
(23) In this system, there are sixteen variables and four dimensions. According to the Buckingham π-theorem, this is rearranged into a relationship among twelve independent dimensionless groups of the original variables. Therefore, the relationship (23) is rewritten in dimensionless form as follows: In the relationship (24), Pr, are the physical property dimensionless group, and β is related to the configuration of the substrate, Bi and Bi w are related to the thermal boundary conditions of the substrate and free surface, which are herein assumed constant. Therefore, in this system, the only free variables are Re, Ma and the disturbance wavenumber, k.

Numerical method
The spectral Galerkin method is used to numerically solve the problem. The velocity components in the stream and the span directions are expanded as The set of trial functions, f j (y/h(x, z, t)),for the velocity components in the stream and the span directions are taken as the following orthogonal polynomials where a j (x, z, t) and b j (x, z, t) are the associated amplitudes. (29) f 0 (ȳ) is a semiparabolic function corresponding to the Nusselt flat film velocity profile. f 1 (ȳ) and f 2 (ȳ) correct the velocity components in the stream and the span directions up to the first order approximation for free surface deformation. However, considering the convenience of working with an amplitude that is homogeneous in the flow rates q and q ⊥ in the stream and the span directions, respectively, the two velocity components are finally expanded as where s 1 , and s 2 are the first order inertia corrections to the Nusselt flat semiparabolic velocity distribution in the stream direction, and r 1 , and r 2 are the corresponding corrections for the velocity components in the span direction. Their expressions are given in Appendix A.
In addition, considering the convenience of working with a surface temperature, θ , that is homogeneous in the temperature field, the temperature field across the film is expanded in the following form.
T(x, y, z, t) The first term on the right side is the linear, zeroth-order temperature profile across the film in terms of the surface temperature, θ, and F is the effective heat flux, defined as follows, The set of trial functions for the temperature field, g j (ȳ), is taken as the following set of orthogonal polynomials, which can also correct the temperature distribution up to the first order approximation.
where t 1 , t 2 , t 3 , and t 4 are the first order corrections to the Nusselt flat linear temperature distribution; their expressions are also presented in Appendix A. The spectral Galerkin method is a particular case of the Weighted Residual method obtained when the set of test functions, w j (ȳ), is chosen as the set of the trial functions separately for the stream-direction momentum equation, the span-direction momentum equation and the temperature equation. After substituting the velocity and temperature expansions, Eqs. (30, 31, 32), into the governing equations and their corresponding boundary conditions, their projection on the set of the test functions gives the equations for h, q , q ⊥ , θ , s j , r j , t j . The resulting system of differential equations is complicated because it has twelve unknowns and a high degree of nonlinearity. Therefore, a certain simplification should be done. Two measures of simplification are taken in this paper without losing the second-order approximation. First, the equations for s j , r j , t j are truncated at their first order approximations based on their nature. It should be remembered that s j , r j , t j are the first order correction to the Nusselt flat film parabolic velocity profile and the linear temperature profile corresponding to the residuals produced from f 0 (ȳ) and g 0 (ȳ). As justified in Kalliadasis et al. (2012), these terms appear through inertial terms involving their space and time derivatives or through products with derivatives of h, q , q ⊥ , and θ , which are terms of O(ε 2 ). Meanwhile, they are absorbed into the stream-direction and span-direction flow rates and the surface temperature through the boundary conditions on the surface in the viscous and conductive terms, which are already of O(ε 2 ). Therefore, the amplitudes s j , r j , t j can be approximated as functions of h,q , q ⊥ , and θ, and their derivatives can be truncated at O(ε). Second, to reduce the singularity due to the strong nonlinearity, the Pade-like approximation approach is used in the streamdirection momentum equation to reduce the degree of the nonlinear terms originating from the inertial term. By introducing a Pade-like regularization factor, the original nine-degree nonlinear term appearing in the inertial term is reduced to a five-degree term. Eventually, four equations for h, q , q ⊥ , and θ are obtained as follows.
From the definition of these dimensionless numbers, Pr, , Bi, Bi w and β are determined once a gas-liquidsolid system is fixed because they depend entirely on the physical property parameters of the system. Re, Ma and the disturbance wavenumber, k, are free variables and can vary with different flow rates, heat fluxes and disturbance frequencies.

Linear stability analysis
Linear stability analyses of the Nusselt flat film solution are carried out to investigate the validity of the above model. In this paper, a temporal stability analysis is made to seek the solution around the Nusselt flat film solution in the form of normal modes: ⎛ where k x, and k z are the real wavenumber components in the stream and span directions, respectively, and ω is the complex angular frequency. The temporal growth rate is the imaginary part ω i of ω. For fixed gas-liquid-solid film flow systems, the physical properties of the fluid, the inclination angle of the substrate and the thermal boundary conditions to ambient air for both the substrate wall and the film surface are fixed. In this paper, the fluid medium is water at 20°C and 1.0 atm; thus, the Kapitza and Prandtl number are taken to be the values = 3375 and Pr = 7, respectively. The inclination angle of the substrate is set to 90°(between the substrate and the horizon), and thus the inclination number is taken to be Ct = 0. On the Bi and Bi w values taken, Our previous paper (Liu et al., 2017), Kobov experiments (Kabov, Scheid, Sharina, & Legros, 2002), and Kallidasis et al. (2012) have all shown that as Bi and Bi w are less than 1, their effects on the Marangoni effect for film flow instability is small. Due to poor heat transfer characteristics at the liquid-gas interface and solid-gas interface, the actual Bi and Bi w values are very small. In the following, the Biot numbers at the substrate wall and free surface are taken to be Bi w = 0.6 and Bi = 0.12, respectively as done by Scheid, Ruyer-Quil, Kalliadasis, Velarde, and Zeytounian (2005). Therefore, a complete instability analysis of the film flow involves identifying the frequency and growth rate of the disturbance mode by varying the Reynolds number, Re, the Marangoni number, Ma, and the disturbance wavenumbers, k x and k z . For the above model and the BE model, the analyses of linear stability of the governing equations and the boundary conditions degenerates into an algebraic relation, ω = ω(k, Re, Ma), which is referred to as the dispersion relation. In the following, the solution to the Orr-Sommerfeld eigenvalue problem (hereinafter abbrev. O-S) and the dispersion relation of the Benney equation are chosen as references to evaluate the validity and the accuracy of the above model. The BE model, derived by perturbation techniques, is an accurate model of film flow for small Reynolds numbers. For the sake of comparison and clarity, the dispersion relations in the stream and the span directions are derived by substituting Eq. (43) into Eqs. (39-42) and imposing the conditions ∂ z = 0 and ∂ x = 0.

Dispersion relation in the stream direction
The dispersion relation in the stream direction is derived by imposing ∂ z = 0 as Eq. (43) is substituted into Eqs. (39)(40)(41)(42). Figure 2 shows the dependence of the growth rate on the wavenumber in the stream direction at small Reynolds number, Re = 1.25, predicted by the above model. Meanwhile, the corresponding results predicted by the solution to the O-S eigenvalue problem and the BE model are also presented for comparison. On the whole, the growth rate predicted by the above model is slightly smaller than the values predicted by both the O-S and BE models. Table 1 shows the deviation of the two most important wavenumbers calculated by the above model and the BE model, i.e. the cut-off wavenumber and the most unstable wavenumber, from the corresponding values predicted by the solution to the O-S eigenvalue problem, which is defined as where k is the wavenumber calculated by the above model or the BE model, and k O−S is the value predicted by the solution to the O-S eigenvalue problem. The most unstable wavenumber and the cut-off wavenumber predicted by the above model are lower by approximately 3%-4% and 2%, respectively, than that of the O-S. In addition, the most unstable wavenumber and the cut-off wavenumber    predicted by the BE model are also lower than the values of the O-S. However, they more closely approach the values of the solution to the O-S eigenvalue problem. Figure 3 shows the effect of the Reynolds number on the growth rate predicted by the above three methods at Ma = 20 where the wavenumbers in the stream and the span directions are taken as k x = 0.25 and k z = 0.25, respectively, which are less than all the cutoff wavenumbers for the investigated Reynolds numbers. It is demonstrated that the prediction of the above model is in good agreement with that of the O-S up to moderate Reynolds numbers. However, the BE model is only applicable at small Reynolds numbers; when the Reynolds number exceeds approximately 10, it is impractical and approaches infinity, losing its validity. Therefore, the above model can give a sufficiently accurate growth rate of the linear instability under the thermocapillary Marangoni effect at both small and moderate Reynolds numbers. Figure 4 shows the effect of the Marangoni number on the linear growth rate at Re = 1.25; herein the stream-direction and span-direction wavenumbers are taken as k x = 0.25 and k z = 0.25, respectively, as above. The range of the Marangoni numbers spans 30, which covers small to large heat fluxes. It is revealed that the predictions of the three methods are in very good agreement with each other, especially at a small Ma number. With the increase in the Marangoni number, the predictions of the above model and the BE model begin to be lower than the value of the O-S, e.g. at Ma = 30 the growth rates predicted by the above model and the BE model are smaller by 2% and 1% than the value of the O-S, respectively.

Dispersion relation in the span direction
As to the linear development of perturbations in the span direction, the three methods agree well with each other with respect to the dispersion relation and the dependence of the growth rate on the Reynolds and Marangoni numbers, as shown in Figures 5-7. In contrast to the perturbation development in the stream direction, the growth rate in the span direction decreases as the Reynolds number increases. Figure 6 demonstrates the effect of the interaction of the hydrodynamic instability and the thermodynamic instability on the perturbation development in the span direction. At small Reynolds numbers, the thermocapillary Marangoni effect is predominant in the span direction. In addition, with an increase in the Reynolds number, the hydrodynamic instability in the stream direction weakens the Marangoni instability in the span direction. Therefore, the increase in the Reynolds number has a mainly destabilizing effect on the instability development in the span direction.

Numerical results
We conducted a three-dimensional numerical simulation of liquid film flow under the HF condition to investigate the effect of the thermocapillary Marangoni effect and initial imposed disturbance on the evolution of film    instability at small Reynolds numbers. The range of the Marangoni number is taken from 2 to 50 and the corresponding heat flux and surface temperature variation, T, are shown in Table 2. The four types of initial imposed disturbances on the film flow are random disturbance, superposition of random disturbance with the most unstable disturbance in the span direction, the most unstable disturbance in the stream direction, and the two most unstable disturbances in both the stream and the span directions in which the most unstable modes are amplitude-dominant over random disturbance with five times the amplitude of that in the random disturbance. In the last three types of initial imposed disturbance, the reason that the random disturbance is still being superposed is to compare with practical experiments because natural random disturbances are ubiquitous in practice.
In performing the calculation, the discretization schemes of the differential operators in the space and time domains are adopted as follows. The computational domain L x × L z is discretized with M × N grid points whose coordinates are x i = iL x /M, and z j = jL z /N (where i = 0, 1 . . . . . . M, j = 0, 1 . . . . . . N). Spatial derivatives are calculated using the Fourier spectral method (Trefethen, 2000). The wavenumber aliasing is treated by keeping the first 2/3 of the Fourier modes in each direction. Grid independence tests are performed by increasing M and N until the variation in the deformation energies, E x and E z , in both the stream and the span directions, as defined in Eqs. (45,46), is limited to 2%. Eventually, the computational domain is taken as 8 times the wavelength of the most unstable disturbance wave in the stream direction and 6 times the wavelength of the most unstable disturbance wave in the span direction, i.e. L x = 2M π /k x (M = 8) and L z = 2N π /k z (N = 6). The time derivatives are numerically solved using the fourth-order Runge-Kutta scheme, and the time step is determined based on the Courant-Friedrichs-Lewy (CFL) condition, which demands that the time step be less than the propagation time of a disturbance traveling over one mesh. The criterion of film rupture is set as the local film thickness being less than h < 10 −3 at which point the computation is terminated.
An initial perturbation of the base-state film thickness is imposed in the following form, h(x, z, 0) = 1 + A x cos(2π n x x/L x ) + A z cos(2π n z z/L z ) + A noiser (x, z), (44) where A x , A z and A noise represent the disturbance amplitudes in the stream and the span directions and in a random direction, respectively, n x and n z represent the wavenumbers in the stream and the span directions, respectively, andr is a random function with values in the interval [−1, 1].
The deformation energies, E x and E z , are defined as the root mean square value of amplitudes covering all modes and all nodes: where E x , and E z are the deformation energies in the stream and the span directions, respectively, and a i , and b i are written as where i is the imaginary unit. The value of the deformation energy represents the degree of the film deformation from the Nusselt flat film. In the case where E x << E z , the film deformation in the span direction is much larger than the deformation in the stream direction, which represents the emergence of a rivulet structure on the film surface. Furthermore, the effect of the film instability on the heat transfer on the film surface between the film and the ambient air is investigated in this research, which is of special interest for engineering applications. The Nusselt flat film is chosen as a reference to analyze the enhancement of heat transfer under film instability, which is measured by the ratio of the heat transfer on the wavy surface of the film to that on the Nusselt flat film, as shown in Eq. (49).
On the right side of Eq. (49), the numerator is the heat transfer on the wavy surface of the film, and the denominator is the heat transfer corresponding to the Nusselt flat film. Eq. (49) is further reduced to Eq. (50) when the Shkodov scaling is used in the temperature.
where A * = A/A w is the ratio of the surface area of the wavy film to the corresponding substrate area spread with the film the latter of which is equal to the surface area of the Nusselt flat film where θ N is the dimensionless surface temperature of the Nusselt flat film based on Shkodov scaling, as shown in Eq. (51), It is shown in Eq. (50) that the heat transfer on the wavy film is enhanced. This enhancement is possibly due to two factors: one factor is the increase in area due to the film surface deformation, i.e. A* > 1, and the other factor is the increase in the temperature difference on the wavy film surface, i.e.θ/θ N ≥ 1.

Random disturbance
The evolution process of the film flow instability and the corresponding contours of the film surface temperature distribution at Re = 1.25 and Ma = 20 under imposition of a random disturbance on the film's surface are shown in Figures 8 and 9 in which the distribution of the film's surface temperature is exactly out of phase with the one for the film's thickness. The evolution process can be categorized into five stages based on the film evolution characteristics. The first stage, covering a time from approximately t = 0 to 1000, as shown in Figure 8(a), features simultaneous development of instabilities in both the stream and span directions whose deformation energies are linearly increased, as shown in Figure 10. As shown in Figure 8(a), the initial two-dimensional, well-arranged wave fronts aligned with the mean flow on the film surface are rapidly distorted into being three-dimensional due to the development of instability in the span direction. The instability in the span direction is driven by the thermocapillary Marangoni effect. The initial disturbance in the span direction on the film surface induces a surface temperature heterogeneity, as shown in Figure 9, which brings about the gradient of the surface tension on the film surface, viz. the thermocapillary Marangoni effect. For heated film, the wave crest on the film surface has higher surface tension than the wave trough. Liquid in the wave crest pulls more strongly on the surrounding liquid than liquid in the wave trough. Furthermore, the thermocapillary Marangoni effect is a self-acceleration behavior, which means that this effect is enforced by the increase in the thickness difference between crest and trough, i.e.  the wave crest is accelerated and becomes thicker, while the wave trough is accelerated and becomes thinner. This process eventually exaggerates the initial disturbance in the span direction, enhances the film thickness variation along the wave front, and makes bumps on the wave front much steeper and cavities on the wave front deeper. Therefore, the propagation speed of the wave front is highly varied and eventually the initial two-dimensional wave front is sheared into being three-dimensional.
In the second stage, covering the time from approximately t = 1000 to 2000 (shown in Figure 8(b)), the deformation energies in both the stream and the span directions stop increasing and remain almost unchanged, as shown in Figure 10. With the steepness increase of bumps, the wave fronts are torn off under the gravitational force, and then the original continuous wave fronts are broken up into short, caterpillar-shaped pieces, which are almost uniformly spread on the film surface, as shown in Figure 8(b). This process is, in fact, a wave breaking, which is justified from the evolution of the deformation energies in both the stream and the span directions in the corresponding period, as shown in Figure 10, over which they remain nearly constant. This observation means that wave generation is counterbalanced with wave breaking during this period. Meanwhile, it is interesting e that the short, caterpillar-shaped wave pattern formed in this stage will soon be reformulated into spanwise-arranged rivulets on the film surface that are spaced almost equally from each other, as shown in Figure 8(c). This reorganization is supposed to be ascribed to the linear filtering of the disturbance in the span direction under the thermocapillary Marangoni effect, which is inferred from the fact that the space between rivulets approximates the wavelength of the most unstable wave mode in the span direction.
In the third stage, covering the time from t = 2000 to 3500 (shown in Figure 8(c)), rivulets are germinated under the thermocapillary Marangoni effect, and meanwhile the deformation energies in both the stream and the span directions also recommence growing, as shown in Figure 10. In the span direction, although the surface wave undergoes braking, the undulation of the film surface in this stage has apparently been much larger than in the initial surface. Thus, the surface temperature gradient is correspondingly higher in magnitude, which makes the thermocapillary Marangoni effect a dominant factor and promotes the development of film instability in the span direction. In the stream direction, the film surface has been flattened by the above wave breaking process, and thus the drag in the stream direction for the film flow is decreased; a situation that is favorable for the development of hydrodynamic instability in the stream direction. Therefore, hydrodynamic instability in the stream direction is also reestablished under the hydrodynamic force. The above can be proved from the evolution of the deformation energies in both the stream and the span directions in the period between t = 2000 and 3500, as shown in Figure 10.
In the fourth stage, covering times from t = 3500 to 4500 and shown in Figure 8(d,e), the deformation energy in the stream direction is decreased, while the deformation energy in the span direction is promptly increased, as shown in Figure 10. In this stage, the wave pattern on the film surface is further reformulated into almost equally spaced and parallel rivulets. The thermocapillary Marangoni effect acts mainly spanwise, as schematically shown in Figure 11. The surface thickness and temperature distribution in both the stream and the span directions are shown in Figure 12(a,b) in which it is shown that the surface temperature variation in the span direction is much larger than the one in the stream direction with magnitudes in the stream and the span directions of 0.07°C and 0.22°C, respectively. Thus, it is estimated that the span component of the surface temperature gradient is approximately three times that of the stream component. This observation is supposedly due to the different heat transfer modes in the stream and the span directions in the film flow. The heat transfer modes are advection and conduction in the stream and the span directions, respectively, and the former of which has a much larger magnitude than the latter. This finding means that the thermocapillary Marangoni effect in film flow heated from the substrate is anisotropic and chiefly applies to the spanwise direction to pump liquid from furrows up to its adjacent rivulets. Therefore, on the one hand, a rivulet crest grows at a much higher rate than its adjacent trough in the same rivulet, and, consequently, the wave front of the rivulet crest becomes much steeper than that of the trough. On the other hand, the difference in the growth rate between the rivulet crest and trough produces increasingly high mean curvatures on the rivulet, which induces a capillary force directed from crest to trough and thus impedes the growth of waves riding on the rivulet. Meanwhile, the flow rate in the furrows is continually decreased due to liquid drawing up to the rivulets. The local Reynolds number in the furrows is correspondingly decreased; therefore, hydrodynamic waves on the furrow are diminished. Thus, in this stage, the deformation energy in the stream direction tends to be decreased in both furrow and rivulet.
In the fifth stage, covering times from approximately t = 4500 to 5000 and shown in Figure 8(e,f), the deformation energy in the span direction continues to increase as sharply as in the last stage, and meanwhile the deformation energy in the stream direction also increases, as shown in Figure 10. The increase in the deformation energy in the stream direction implies that there are new hydrodynamic waves generated on the film surface. As time goes on, liquid in the furrows is incessantly pulled up to the rivulets until they approach dryout. Eventually, a wave pattern on the film surface evolves into an alternate stripe-like rivulet structure, which is in good agreement with the photograph of FC-72 liquid film flow at small Reynolds number, Re = 5, and low heat flux, q = 0.32 W/cm 2 , taken with an infrared scanner by Chinnov & Kabov, as shown in Figure 4 and in reference (Chinnov & Kabov, 2003). As to the flows in the rivulet, the local Reynolds number becomes larger than the one in the mean film flow. A new hydrodynamic wave in the stream direction, driven by hydrodynamic forces, begins to grow on the rivulet. Finally, a typical solitary-like wave pattern, riding on the rivulet and preceded by many small capillary ripples, is formed. Consequently, the deformation energy in the stream direction is again increased in this stage, while the deformation energy in the span direction is still rapidly increased, as in the last stage.

Random disturbance superposed with amplitude-dominant most unstable disturbance in the span direction
The evolution of film instability under the superposition of a random disturbance with the amplitudedominant (A z > A noise ) most unstable disturbance in the span direction is shown in Figure 13. The deformation energies in both the stream and the span directions are shown in Figure 14. In comparison with the above case under random disturbance, a wave in the span direction driven by the thermocapillary Marangoni effect soon becomes visibly evident in the initial period even though the hydrodynamic wave in the stream direction has not yet developed to a level of being perceptible by the naked-eye, as shown in Figure 13(a). Afterwards, a hydrodynamic wave in the stream direction is developed enough to be visible riding on the rivulet aligned with the mean flow, as shown in Figure 13(b). Meanwhile, there are hydrodynamic waves engendered on the furrow, although, they have smaller amplitude compared to the rivulet wave. This finding is supposed to be due to a locally high Reynolds number in the rivulet and a very small one in the furrow. As time goes on, the hydrodynamic wave in the stream direction riding on both rivulet and furrow is suppressed, as shown in Figure 13(c). The reason for this result is the same as stated above in the fourth stage under a random disturbance. The subsequent evolution process of the film flow becomes almost identical to the case under random disturbance. Finally, well-ordered rivulets are formed on the film surface with solitary-like waves riding on their crests, as shown in Figure 13(e,f).
The deformation energies in both the stream and the span directions are shown in Figure 14. It can be seen that, from the beginning, the span-direction deformation energy exceeded the stream-direction deformation energy. This result implies that the rivulet structure on the film surface was formed in the beginning.

Random disturbance superposed with amplitude-dominant most unstable disturbance in the stream direction
The instability evolution process on the film under the amplitude-dominant most unstable disturbance in the stream direction superposed with random disturbance is shown in Figure 15. Except for in the initial stage, the instability of the film flow evolves in approximately the same way as in the case under random disturbance. In the initial stage, the instability development on the film  surface is chiefly identical to the isothermal film flow in which a series of wave fronts in the stream direction is formed on the film surface. Up to about t = 1500, wave fronts aligned with the mean flow are torn off into being three-dimensional under the impact of the thermocapillary Marangoni effect in the span direction. The deformation energies in both the stream and the span directions are shown in Figure 16, which shows that a rivulet structure on the film surface is formed at approximately t = 3200.

Random disturbance superposed with amplitude-dominant most unstable disturbances in both the stream and the span direction
The evolution of film instability is shown in Figure 17 under a random disturbance superposed with the two amplitude-dominant most unstable disturbances in the stream and the span directions. The two most unstable disturbances have the same amplitude, which is 5 times that of the random disturbance. In such a situation, a very regular checkerboard structure is initially formed on the film surface, as shown in Figure 17(a), which is driven simultaneously by hydrodynamic force and the thermocapillary Marangoni effect independently in the stream and the span directions. Then, the dependence of the phase velocity on the film thickness causes the above pattern to be further transformed into a regular, teardrop-shaped, wavy structure, as shown in Figure 17(b). The subsequent evolution process is similar to the above cases, as shown in Figure 17(c-e). Finally, a rivulet structure on the film surface is formed on which solitary waves proceeded by small ripples ride, as shown in Figure 17(f). A rivulet structure is formed at approximately t = 1000, and the deformation energy in the span  direction exceeds that in the stream direction, as shown in Figure 18.

Heat transfer on wavy surface of film
In practical engineering applications, we are concerned with heat transfer on the film surface and the time duration of the film integrity. For heat transfer on the wavy surface of the film, the evolution of E HT on the film surface is shown in Figure 20 under the above four initial disturbances. It is certain that the rivulet formation on the film enhances the heat transfer on the wavy surface; although, there is a slight difference in the E HT under different initial disturbances (disturbance including a component in the span direction engenders a value of E HT that is slightly higher than without it), as shown in detail in Table 3. As orderly rivulets are developed on the film for all of the four initial disturbances imposed on the film, the E HT begins rapidly surging until the film ruptures, as shown in Figure 20. As to the mechanism of the enhancement of heat transfer on the film surface, it is primarily attributed to the increase in the surface area of the film and in part to the increase in surface temperature, as shown in Table 3. The variation of the area of the film surface along with the E HT under a random disturbance is shown in Figure 19, which manifests itself in the same way as the development of E HT . As to the rupture time of the film, a disturbance with a component in the span direction renders the film susceptible to the occurrence of spot dry-out, as shown in Table 3. The duration for which the film retains integrity for the initial disturbance with a component in the span direction is about half the time for the disturbance without it, which is unfavorable for engineering applications.
Note: point set 1 refers to the present calculation, point set 2 refers to the experimental data (Chinnov & Kabov, 2003) with water at a heater dimension of 150 × 150 mm, point set 3 refers to the experimental data (Chinnov, 2009) with water at a heater dimension of 60 × 120 mm, the curve is fitted using the Kabov group experimental data (Kabov & Chinnov, 2001).
at Re = 1.25, 3 and 4 with the same Marangoni number Ma = 20, respectively. It is shown that the rivulet-to-rivulet spacing increases with the Reynolds number, which is an important parameter for deeper insight into the mechanisms of heat transfer and crisis phenomena during heat transfer to liquid film and is strongly dependent on the Reynolds number of the film flow (Chinnov, 2017). Figure 22 shows a plot of the dimensionless rivulet-to-rivulet spacing, /l σ , versus the Reynolds number in which the data with black solid squares are from the present calculation, the data with red solid circles are from the Chinnov & Kabov experiments with water (Chinnov & Kabov, 2003), the data with blue solid triangles are from the Chinnov experiments with water (Chinnov, 2009), and the solid curve is the fit using the Kabov group experimental data at Re = 0.42-24 with a confidence of 0.78 (Kabov & Chinnov, 2001). It can be seen from Figure 22 that the data from the present calculations fall into a reasonable range compared with the experimental data. Note: = 3375, β = π / 2, Pr = 7, Bi = 0.12, Bi w = 0.6, L x × Lz = 2n x π /k x × 2n z π /k z , n x = 8, n z = 6, k x = 0.3891, k z = 0.2997, A x = 0, A z = 0, A noise = 0.001, M × N = 128 × 128.
Meanwhile, at higher Reynolds numbers, flow behavior and rivulet structure become more complicated than at smaller Reynolds numbers, e.g. at Re = 4 no regular solitary-like structures riding on the rivulet appear and a strong three-dimension irregularity exists on the crest of the rivulet, as shown in Figure 21(c). Figure 23(a-c) show wave patterns at Ma = 20, 30 and 40 with the same Reynolds number, Re = 1.25. Two aspects of the effect of the Ma on the wave pattern are noted: one aspect is that the number of rivulets per width increases with Ma because the wavenumber of the most unstable modes under the thermocapillary Marangoni effect is increased with Ma; the other aspect is that increasing the Ma leads to an advance in the occurrence of film rupture, e.g. at t = 5156 for Ma = 20, and t = 2280 for Ma = 30. Even at Ma = 40, the furrow thickness has approached the critical rupture condition at some points, while regular rivulet patterns on the film surface have not yet completely formed.

Conclusions
A three-dimensional, regularized model with heat flux boundary conditions is developed by extension of the existing two-dimensional regularized model. This model is used to study liquid film flowing down a uniformly heated, vertical substrate with constant heat flux and to explore the mechanism of rivulet formation. The conclusions are as follows: (1) The model is verified by linear stability analysis in comparison to the BE model and the solution to the O-S eigenvalue problem. Linear stability analysis proves that the model is applicable from small to moderate Reynolds numbers with high accuracy in reference to the solution to the O-S eigenvalue problem. The model has better performance than the BE model for moderate Reynolds numbers. Meanwhile, as to the prediction of the dispersion relation in the stream direction, it approaches more closely the solution to the O-S eigenvalue problem than the BE model. As to the linear development of the perturbation in the span direction, the three methods agree well with each other in both dispersion relation and dependence of growth rate on the Reynolds and Marangoni numbers. This finding is because the thermocapillary Marangoni effect is entirely dependent on the disturbance development in the span direction.
(2) A rivulet on the film surface, aligned with the mean flow, is formed for a certain range of the Reynolds and Marangoni numbers. Although the process of rivulet formation is slightly different for different initial disturbances imposed on the film, five stages of evolution can be roughly discerned based on the deformation energies and film structure. The thermocapillary Marangoni effect, which is the main driving force in film surface rivulet formation, mainly acts spanwise on the film surface due to weak heat transfer in the span direction, which makes the surface temperature gradient have a higher component in the span direction than in the stream direction. (3) As to the effect of initial imposed disturbance on the evolution of film instability, four types of initial disturbances are investigated. The initial imposed disturbance certainly affects the evolution process; however, it has little effect on the eventual rivulet structure on the film surface, which is fully controlled by the Reynolds and the Marangoni numbers. (4) As to the heat transfer on the wavy film surface, rivulet formation on the film enhances heat transfer, which is attributed for the most part to the increase in the surface area of the film. As to the rupture time of the film, the disturbance with a component in the span direction renders the film susceptible to the occurrence of spot dry-out. (5) As to the effect of the Reynolds and the Marangoni numbers on rivulet formation. Rivulet width incre ases with Reynolds number, while it decreases with Marangoni number. Flow behavior and rivulet structure become more complicated with the increase in the Reynolds number. In addition, up to a certain value of the Reynolds number, a rivulet may not even be formed and the film evolves in almost the same way as in isothermal film flow. In addition, the increase in Marangoni number causes the film rupture occurrence to advance.
Finally, it should be mentioned that the falling film in this study is assumed to be exposed to an air atmosphere, and thus the heat transfer of the falling film to ambient air is via natural convection. Thus, the Bi and Bi w numbers are set as constants. However, in engineering applications, e.g. film heat exchangers, the falling film is exposed to flowing vapor or other gases, and the Bi and Bi w may vary along the falling film flow. As is known, the surface temperature gradient may be reduced by forced convection. Therefore, it is necessary to consider the effect of variable Bi and Bi w on rivulet formation in film flow for engineering applications. This will be investigated in future works.