A discussion on the validity domain of the weighted residuals model including the Marangoni effect for a thin film flowing down a uniformly heated plate

ABSTRACT In this paper, the validity of the Weighted Residuals Model (WRM), including the Marangoni effect, is investigated through linear stability analyses and two-dimensional nonlinear numerical simulations. The linear stability analyses with the WRM show that the model's accuracy decreases nearly linearly with the Marangoni number for each Reynolds number and achieves a maximum at approximately for each Marangoni number. This is quite different from the isothermal case, for which the error increases monotonically with the Reynolds number and remains small for small-to-moderate Reynolds numbers. This is very important for application of the WRM but has yet to be reported or investigated. The effects of the Reynolds number and Marangoni number on the nonlinear evolution of film layers are then investigated through numerical simulations. At small Reynolds number, it is found that the error caused by the Marangoni effect in predicting the phase speed can be ignored if the Marangoni number is small () or makes the wave in the spatial numerical simulation considerably out of phase if the Marangoni number is large (). On the other hand, the saturation states can be generated by the WRM no matter whether the Marangoni number is small or large. When the Reynolds number is increased to a moderate value and the Marangoni number is taken as zero, it is found that the saturation wave produced by the WRM is very similar to the experimental one, except for the amplitude of the wave being somewhat larger and the wave speed as well as the wavelength being slightly smaller. Hence, it can be inferred that the WRM predicts the saturation waves well for small-to-moderate Reynolds numbers if the Marangoni numbers are limited to a small range depending on the Reynolds numbers.


Introduction
Falling film flow is a flow in which the thickness is small compared to the characteristic length scale. This has intrigued numerous researchers due to its relevance to certain practical applications and fundamental problems of fluid mechanics. Falling film is widely applied in practice, such as cooling towers, condenser tubes, and compact multi-phase heat exchangers. A falling film has a simple configuration but exhibits rich dynamics analogous to other classical hydrodynamic instabilities. It was the first open-flow hydrodynamic instability to be fully deciphered from inception to fully developed spatial-temporal dynamics. It is hoped that a comprehensive examination of this intriguing instability will inspire similar approaches to other hydrodynamic instabilities (Chang & Demekhin, 2002). In 1949, the Soviet Union physicists Pyotr and Sergey Kapitza (Kapitza & Kapitza, 1949) observed that, when a liquid film falls down an inclined plate, a rich wave dynamics arises at the surface of the thin film. Yih (1963) performed a linear stability analysis and showed that a vertical falling film is always unstable to long-wave disturbances at any Reynolds number. These unstable waves are known as interfacial modes, and they are the main instability for small-to-moderate Reynolds numbers. When the Reynolds number is increased further, a shear mode may arise in the interior of the film (Bruin, 1974). However, in most practical applications, the flow rates, or equivalently the Reynolds numbers, are not too large. It is thus adequate to consider the interfacial mode only. Benney (1966), based on an assumption of small Reynolds number and long waves, performed a perturbation expansion in terms of the so-called 'film parameter' (the ratio of the film thickness to the streamwise length scale) and derived an evolution equation for film height, known as the Benney equation (hereinafter referred to as the BE). By limiting the wave amplitude at different magnitude orders, the BE can be simplified into a linear wave equation, the Kuramoto-Sivashinsky equation (weak nonlinearity), the KDV equation (strong nonlinearity), and the Burger equation. Moreover, the linear form of the BE exactly predicts the growth rate and critical condition of instability obtained from the eigenvalue problem of the Orr-Sommerfeld equation, which is derived from the linearized Navier-Stokes/energy equations and is always accurate under all flow conditions. The BE predicts the existence of linear periodic waves, traveling waves, and blow-up. In subsequent studies (Cheng & Chang, 1991;Pumir, Manneville, & Pomeau, 1983;Nakaya, 1989) based on the BE, the periodic traveling waves were constructed numerically and theoretically and found to be in good agreement with observations in experiments. Involving only the film thickness and being accurate at small Reynolds number makes the BE fascinating. Many scholars extended the BE to including more physical effects, such as uniform heating as well as evaporating (Joo, Davis, & Bankoff, 1991) and non-uniform heating (Scheid, Oron, Colinet, Thiele, & Legros, 2002). Pumir et al. (1983) integrated the BE numerically and found a singularity; the BE exhibits an unphysical finitetime blow-up behavior at somewhat large Reynolds number. Several authors extended Benney's work to include more higher-order terms in terms of film parameters but still failed to overcome the singularity. Ooshida Takeshi (Takeshi, 1999) realized that the perturbation expansion leading to the BE is divergent, and he proposed a Padé-like regularization method to replace the divergent asymptotic series. Although the equation that he obtained removes the singularity of the finite-time blowup, it does fail to describe accurately the dynamics of a film at moderate Reynolds numbers because its solitarywave solutions exhibit unrealistically small amplitudes and speeds. However, it should be noted that the BE remains valid if the wave amplitude and Reynolds number remain small.
In 1967, with an assumed self-similar semi-parabolic velocity profile and Karman-Pohlhausen averaging method, Shkadov (1968) derived an integral-boundarylayer model for the equation of motion, which, together with the integral continuity equation, constitutes a set of evolution equations for the interface height h(x, t) and local flow rate q(x, t). This formulation was first introduced by Kapitzas to describe stationary waves. Hereinafter, we refer to this model as the Kapitza-Shkadov (K-S) model. The derivation of the K-S model, based on the assumption of long-wave with a self-similar velocity profile, places no restriction on the Reynolds number, Weber number, or wave amplitude. The K-S model captures the fully nonlinear dynamics in the region of moderate Reynolds number (Lee, & Mei, 1996). However, the K-S model does not predict the growth rate of an instability at small Reynolds numbers or the neutral condition well, except for the vertical case being considered.
The shortcoming of the K-S model has been removed by combining the gradient expansion with a 'weighted residuals' technique using polynomials as test functions for the velocity field (Ruyer-Quil & Manneville, 1998, 2000. The resulting models will be referred to as Weighted Residuals Models (WRMs), which are characterized by different complexities and accuracies according to the approximation equations and test functions adopted. Among them, the 'simplified second-order model' is the most widely used because of its eclectic complexity and accuracy. Hence, we will refer to the 'simplified second-order model' as the 'WRM' for simplicity. The WRM is similar to the K-S model but with different numerical coefficients of the terms. It correctly predicts the linear instability threshold for all inclination angles as well as the growth rates of the instability, and it is also applicable to moderate Reynolds numbers.
Moreover, in the limit of vanishing Reynolds numbers, a perturbation expansion of the WRM can lead back to the BE (Kalliadasis, Ruyer-Quil, Scheid, & Velarde, 2012). Hence, it is clear that the WRM covers the BE and is applicable for small-to-moderate Reynolds numbers. Scheid (2004) improved the WRM by introducing the Marangoni effect. However, its validity domain has not yet been examined. The analyses in the present paper show that the accuracy of the WRM decreases as the Marangoni number is increased if the Marangoni effect is included.
From the derivation of the WRM, it can be seen that the validity of the WRM depends on the distributions of velocity and temperature. Perfect description of the film flow indicates excellent agreement between the assumed velocity (temperature) distribution and the actual one. Thus, the WRM has an advantage over the Navier-Stokes equation in respect of its physical interpretation of film flow. Furthermore, the WRM is also practically useful for its simplicity compared with the Navier-Stokes equations and its satisfactory accuracy in describing the nonlinear dynamics of a film layer. This is more obvious if the WRM is extended to threedimensional cases. A three-dimensional numerical simulation with a three-dimensional WRM (Scheid, Kalliadasis, Ruyer-Quil, & Colinet, 2008) or the Benney type equation (Mayo, Mccue, & Moroney, 2013) is usually much simpler than a corresponding numerical simulation with the Navier-Stokes equation combined with the VOF approach (Xu, Yuan, Repke, & Wozny, 2012). Hence, an examination of the validity of the WRM is important for the physical interpretation of film flow and the practical application of the WRM.
In this paper, the validity of the WRM is investigated through linear stability analyses and nonlinear numerical simulations. The paper is organized as follows. First, we formulate the problem and introduce dimensionless groups. By conducting a perturbation analysis, we show that the WRM is exact if the Reynolds number and modified Marangoni number are sufficiently small. Second, we conduct temporal linear stability analyses with the WRM to investigate its accuracy in predicting the growth rates and phase speeds under larger Reynolds numbers and Marangoni numbers by comparing with the results predicted by the eigenvalue problem of the Orr-Sommerfeld equation. Then, we perform numerical simulations with the WRM to investigate its performance in predicting wave shapes and wave speeds by comparing with the solutions of the BE, which is accurate at sufficiently small Reynolds numbers, and published experimental results with moderate Reynolds numbers. Finally, we summarize our results and present some conclusions.

Formulation
We consider an incompressible viscous fluid of constant density ρ and viscosity μ freely falling under gravity g along a uniformly heated inclined plate at an angle β with respect to the horizontal plate and having a free interface, y = h(x, t), as shown in Figure 1. The waveless flat film has a thickness h N , which is known as the Nusselt film thickness and is related to the mean velocity, or the Nusselt velocity, of the semi-parabolic profile by U(y) = h 2 N g sin β(2ỹ −ỹ 2 )ρ/2μ, whereỹ = y/h N is the reduced coordinate. The film layer has a fixed constant temperature T w at the wall, and the temperature of the liquid on the free surface is T s . The thermal conductivity of the liquid is λ, the heat capacity is c p , the thermal diffusivity is χ = λ/ρc p , and the heat transfer coefficient that reflects the rate of heat transport between the liquid film and ambient air is α; σ = σ a − γ (T s − T a ) is the surface tension, σ a is the surface tension at air temperature T a , and γ = −(dσ/dT s )| T a is positive for most liquids.
The final film thickness h, velocity profile u, and temperature distribution T depend on the physical properties of the fluid, the acceleration due to gravity g, the angle of incline β, the flow rate q N , the temperature difference between the air and the wall T = T w − T a , and the wavenumber k or frequency f of the disturbance wave. This can be written as follows (Tan, 2011): Note that ρ, μ, λ, γ , g are dimensional independent parameters, scaled by an expression (1) that can be written as follows (Tan, 2011): where are the viscous-gravity length and viscous-gravity velocity, respectively. -Pr = ν/χ is the Prandtl number, which compares the momentum diffusivity to the thermal diffusivity; Γ = σ a l ν /ρν 2 is the Kapitza number, which compares the surface tension force to the force of inertia; Bi = αl ν /λ is the Biot number, a dimensionless heat transfer coefficient describing the rate of heat transport from liquid to ambient air; Ma = Tγ l ν /ρν 2 is the Marangoni number, which compares the force induced by the surface tension gradient to the inertia; Re = q N /ν = u N h N /ν (here, u N is the average velocity, which is 2/3 of the waveless surface velocity) is the Reynolds number, which can also be interpreted as the dimensionless flow rate; and ω = fl ν /u ν is the dimensionless disturbance frequency. For fixed air-liquid-solid systems, Γ , Pr, Bi, β are all fixed and independent of the flow; the only values that vary are the Reynolds number Re, the Marangoni number Ma, and the disturbance frequency ω or wavenumber k. If the liquid phase is water at 20 • C, Γ ≈ 3375 and Pr ≈ 7. Air is expected to be a poor conductor, and thus the parameter Bi should be small, e.g. Bi = 0.1. In the absence of experimental values for the physical property parameters, we take Γ = 3375, Pr = 7, Bi = 0.1 for convenience. The angle of incline is fixed at β = π/2; then cot β = 0. The temperature difference T w − T a between the wall and ambient air is considered to be less than 5 • C (Ma ≤ 50). The liquid film is assumed to be sufficiently thin for the buoyancy effect to be neglected if it is heated from below.

The weighted residuals model
The governing equations are the continuity equation, the Navier-Stokes equation, and the energy equation. The boundary conditions at the wall are no-slip, nopenetration for velocities, and a constant value for the temperature. The air is assumed to be mechanically 'passive', so that the influence of the viscous stress of the air on the film layer can be ignored.
The system constituted by the equations and boundary conditions can be written in dimensionless form by taking the Nusselt scaling: where t ν = l ν /u ν is a timescale. The system can also be non-dimensionalized using the Shkadov scaling: The primed variables are dimensionless, and the primes will be dropped in the following for simplicity; ε is a small number called the 'film parameter', which measures the ratio of the film thickness to the streamwise characteristic length scale. The film parameter will be present as coefficients in the equations and boundary conditions. We can estimate the magnitude orders of different terms and make a somewhat significant simplification by omitting the higher-order terms. Shkadov scaling, introduced in 1977 by Shkadov, who realized the significance of balancing gravity and viscous drag with surface tension for large-amplitude solitary waves and took h N /ε (where ε = (We) −1/3 and We is the Weber number) as the streamwise characteristic length scale. This scaling makes apparent the balance between all forces, i.e. inertia, gravity, viscosity, and surface tension, which is necessary to sustain strongly nonlinear waves. Furthermore, this scaling will make h,q,T as well as their derivatives retain values close to unity, which makes it favorable for numerical integration. Shkadov scaling is useful for theoretical analyses and numerical integrations, and the use of Nusselt scaling is convenient for comparing with experimental results. Thus, we will use these two scalings alternately. The results shown in the figures are in terms of Nusselt scaling to make them intuitive. Assuming the velocity distribution and temperature distribution and applying the 'weighted residuals' technique, Scheid et al. (Kalliadasis et al., 2012;Scheid, 2004) derived a set of model equations in a two-dimensional form: where h is the local film thickness, q is the local is a modified dimensionless thickness (which compares the thickness to the viscosity diffusion length and relates the Reynolds number to other parameters M, We, B), δ = 3εRe is the reduced Reynolds number (which can be interpreted as a continuation parameter measuring the transition from the 'drag-gravity regime', where the streamwise component of gravity is mainly balanced by the viscous drag, to the 'drag-inertia regime', where inertia balances viscous drag, as pointed out by Kalliadasis et al., 2012, andScheid et al., 2008), Ct = cot β is required to be less than 56 because β is required to be larger than 1 • (as pointed out by Bruin, 1974), η = ε 2 , and ε is the film parameter. The Reynolds number plays an crucial role in determining the validity of the WRM because the modified film thickness h N = (3Re) 1/3 is relevant to other parameters We, ε, M, B. For given physical property parameters Γ , Pr, Bi and inclined angle β, the validity of the WRM is only affected by the Reynolds number and the Marangoni number. This is very useful for analyses because, once the physical property parameters are fixed, the validity of the WRM is determined entirely by the Reynolds number and Marangoni number. The first equation in (4) is the continuity equation, and the other equations are the momentum and energy equations. This model is derived from the second-order approximation of the equations and boundary conditions and it can be called a 'simplified model'.
It should be noted that the velocity distribution U(y) = (y − y 2 /2) (where y = y/h) has been verified by direct simulation (Salamon, Armstrong, & Brown, 1994) but as yet the temperature distribution has not. Since the waveless film (y ≤ h ≡ 1) has a velocity distribution U(y) = (y − y 2 /2) and a temperature distribution T(y) = 1 − By/(1 + B) (Kalliadasis et al., 2012), it is reasonable to expect that the temperature distribution should take the form of 1 + b 0 y if the temperature difference T w − T a is small enough.

Benney equation
We have mentioned that, if the Marangoni effect is not included, the WRM becomes equivalent to the BE as the Reynolds number tends to zero. As Re → 0, h N = (3Re) 1/3 tends to zero, and We = Γ /h 2 N tends to infinity. When ε 2 We ∼ O(1), ε, 3εRe, ε 2 , ε cot β → 0. Following the method used in Kalliadasis et al. (2012), we can perform a gradient expansion of q, θ in system (4) (1). Substituting q = q 0 + εq 1 and θ = θ 0 + εθ 1 into the momentum and energy equations (4b) and (4c), and collecting terms of the same power in ε to construct a set of equations, one can solve these resulting equations successively.
For the equation with the terms of ε 0 , we obtain Substituting these into the equation with the terms of ε 1 , we obtain Substituting q = q 0 + εq 1 into the evolution equation (4a), one can obtain the BE for non-isothermal cases (Scheid et al., 2005): The equation has a structure similar to the BE for the isothermal case but with some additional terms to account for the Marangoni effect. The BE for the isothermal case is recovered if Ma = 0. For the isothermal case, the BE is exact in the limit of vanishing Reynolds number; however, it breaks down at an O(1) Reynolds number. The same is true for the non-isothermal case, as shown by Kalliadasis, Demekhin, Ruyer-Quil, and Velarde (2003). By conducting a perturbation analysis of the WRM, we have shown that the equivalence between the WRM and the BE requires Re → 0, M ∼ O(1), RePr ∼ O(1), Ct ∼ O(1). Note that the physical property parameter Pr is usually not large and both the WRM and the BE are derived implicitly requiring β > 1 • (Ct < 56) to obtain a semi-parabolic velocity distribution. We will not discuss the range of Pr and Ct for convenience. Because the accuracy of the BE is based on only the assumption of small Reynolds number and long waves, and is independent of and Bi, it can be inferred that the WRM is exact when Re → 0, M ∼ O(1) , no matter what and Bi are. However, we remain uncertain about its validity when the Reynolds number and the Marangoni number are increased. In the following sections, we will examine the validity domain of the WRM in terms of the Reynolds number and Marangoni number in more detail through linear stability analyses and nonlinear numerical simulations.

Linear stability analyses
To investigate the validity range of the Reynolds number and the Marangoni number, we fixed Γ = 3375, Pr = 7, Bi = 0.1, Ct = 0 for convenience.
If the Marangoni effect is not included, i.e. Ma = 0, the energy equation can be removed. Substituting h = 1 + he i(kx−ωt) , q = 1/3 +qe i(kx−ωt) (i = √ −1,ĥ andq are disturbance amplitudes) into the continuity and momentum equations (4a) and (4b), we obtain the dispersion relation in terms of the Shkadov scaling: By replacing ω with ω/ε and k with k/ε, the dispersion relation can be changed into one in terms of the Nusselt scaling: Only the temporal mode is considered in this paper since a simpler dispersion relation can be obtained. The k is taken as real, and ω is imaginary. This results in a quadratic equation with respect to ω, which can be solved with the radical formula. The growth rates versus the wavenumber predicted by the WRM, the BE, and the Orr-Sommerfeld (O-S) equation are shown in Figure 2. The BE and the WRM are very accurate if the Reynolds number is sufficiently low (Figure 2(a)). As the Reynolds number is increased to a moderate value (Re = 20), the BE quickly loses its validity, although the WRM still gives an approximate result (Figure 2(b)). The Error Ratio (ER) of the WRM, which defined as where ω(k) i,WRM is the growth rate predicted by the WRM and ω(k) i,OS is that by O-S, increases with the wavenumber before becoming stable; however, it remains small(< 2%) at the most unstable wavenumber, which plays a key role in the evolution of the film layer. Hence, we can conclude that the WRM gives a sufficiently accurate growth rate of the linear instability at both small and moderate Reynolds numbers if the Marangoni effect is neglected.
If the Marangoni effect is considered, we can obtain a cubic equation for ω by substituting h = 1 +ĥe i(kx−ωt) , q = 1/3 +qe i(kx−ωt) , θ = (1 + B) −1 +θe i(kx−ωt) (ĥ,q, andθ are disturbance amplitudes) into the WRM including the Marangoni effect. Its solutions are substantially more complex than the quadratic equation and will be given in a numerical instead of an analytical form. Figure 3 shows the growth rates and error ratios (ERs) as a function of wavenumber under small Reynolds numbers with different Marangoni numbers. The WRM, BE, and O-S give nearly the same results if the Marangoni number is sufficiently small (M = 2.07 ∼ O(1)), as shown in Figure 3(a). Otherwise (M = 20.7 ∼ O(10)), the growth rate of the WRM clearly deviates from the exact rate, as shown in Figure 3 (b). This is consistent with the conclusion that we have mentioned before: the WRM is accurate if Re → 0, M ∼ O(1). However, when the Reynolds number is increased to Re = 20, the error of the WRM becomes non-ignorable (Figure 4 Because the modified Marangoni number is a function of the Reynolds number and the Marangoni number (M = Ma/(3Re) 2/3 ), the Reynolds number and the Marangoni number will be investigated instead of the modified Marangoni number. To assess the error caused by the Reynolds number and Marangoni number in detail, we introduce a new variable called the Maximum Error Ratio (MER) to measure the error of the WRM in predicting the growth rate of the most unstable mode. The MER is defined as where ω imax,WRM and ω imax,OS are the growth rates of the modes with the most unstable wavenumbers predicted by the WRM and the O-S, respectively. It should be noted that the most unstable wavenumber predicted by the WRM is not necessarily equal to that of the O-S. Figure 5 shows the maximum error ratio (MER) as a function of the Marangoni number. The MER increases nearly linearly with the Marangoni number, and the rate of increase varies with the Reynolds number, indicating a varied validity domain of the WRM in terms of the Marangoni number depending on the Reynolds number. It should be noted that the WRM achieves nearly the same accuracies at both Re = 1 and Re = 20 and presents   the lowest accuracy at Re = 4 among the four cases considered in Figure 5. The MER increases quickly with the Marangoni number at Re = 4. Figure 6 shows the MER as a function of the Reynolds number. When the Marangoni effect is not included, the MER increases with the Reynolds number but remains as a small value (< 1%), indicating the applicability of the WRM for small-to-moderate Reynolds numbers. When the Marangoni effect is included, the MER achieves its maximum at Re ≈ 4, which is nearly independent of the Marangoni number, indicating that the lowest accuracy of the WRM occurs at approximately Re = 4. This is consistent with the finding from Figure 5 and quite different from that without the Marangoni effect.
We have investigated the performance of the WRM in predicting the linear growth rates under different Reynolds numbers and Marangoni numbers. The accuracy of the WRM is found to decrease proportionally to the Marangoni number and achieves a minimum at approximately Re = 4. Because the linear analysis is valid only if the disturbance amplitude remains small, a nonlinear numerical simulation with the WRM is required to complete this study.

Numerical simulations
In this section, we pose a series of boundary-value problems and initial-value problems that probe the nonlinear evolutions of the film flow. Boundary-value problems correspond to the spatial analyses in linear stability theory, and the initial-value problems correspond to the temporal analyses. Boundary-value problems require addressing inlet and outlet boundary conditions and a large computational domain to allow for sufficient evolution of the waves. They are more suitable for comparison with experiments, but they are time consuming. Initial-value problems apply periodic conditions in the spatial direction and require a smaller computational domain. They are more suitable for comparison with temporal stability analyses. The quasi-wavelet method and the Fourier spectral method are used to approximate the spatial derivatives for spatial and temporal numerical simulations, respectively. Time advancing is achieved by a fourth-order Runge-Kutta method.

Small Reynolds number
In this section, the validity of the WRM is examined through comparison with the BE, which is exact at small Reynolds numbers for traveling waves (Scheid et al., 2005).
For Re = 1.25, the most unstable wavenumber k is 0.026 (see Figure 3(a)), and the corresponding frequency ω r is approximately 0.026 at Ma = 5; the most unstable wavenumber k is 0.043 (see Figure 3(b)), and the corresponding frequency ω r is approximately 0.043 at Ma = 50. In these two cases, the reduced Reynolds number δ = 0.3354, lying inside the validity domain of the BE, which is valid up to δ ∼ 1 (Kalliadasis et al., 2012;Scheid et al., 2005).
A spatial numerical simulation is preferred to reproduce experimental results. The solutions (h, q, θ) are approximated by where f (x, t) represents h, q, θ , the value of f (x, t) at x k isf k (t), and 2W + 1 ≥ 35 is the number of nodes used to approximate f k (t). A temporal periodic forced flow rate q = q N (1 + A sin(ωt)) is imposed upstream of the filmwith with an initial waveless height h = 1, where q N = 1/3 is the average flow rate, A is the disturbance amplitude, and ω is the frequency. Due to the convective nature of the instability, the disturbance amplitude grows downstream gradually, as shown in Figure 7. For Ma = 5, the wave predicted by the WRM fits well with that predicted by the BE, as shown in Figure 7(a). For Ma = 50, the WRM predicts a lower growth rate and a lower phase speed, as shown in Figure 7(b). Clearly, the error of the WRM in the growth rate fades away as the wave tends toward saturation, while the error in the phase speed makes the wave become out of phase as the wave travels downstream. Now we consider the temporal evolutions of films. These are initial-value problems and the periodic boundary conditions are applied in the spatial direction. The simulations are started with h(x, 0) = 1 + A sin(kx), q(x, 0) = 1 3 , θ(x, 0) = (1 + B) −1 . When a harmonic disturbance A sin(kx) (where 0 ≤ x ≤ 2π/k) with the most unstable wavenumber k predicted by linear analysis is added to the flat film, the film elevation h(x, 0) = 1 + A sin(kx) will grow over time due to the intrinsic instability. The solutions of the WRM are approximated by the finite Fourier series h(x, t) = 1 + N n=1 a n (t)e inkx + c.c., with N ≥ 80 as the number of Fourier modes, which is determined by increasing N until no quantitative differences in the Fourier coefficients are obtained. As time progresses, higher modes are excited, and the wave shape tends to a saturation state. Figure 8 gives the evolutions of the Fourier spectral coefficients as a function of time.
It can be observed that each wave tends to the saturation state after a sufficiently long time, and only a few modes are adequate to restructure the saturation waves.   The final permanent waveforms predicted by the BE and the WRM are shown in Figure 9, where excellent agreements are achieved. The growth rates and phase speeds do not affect the final saturation states. It can be inferred that the temporal simulations with the WRM still predict the saturation states well, even if the errors of the WRM in predicting the growth rates and phase speeds increase. However, this may not be true if the Reynolds number is increased since the flow regime may change into the drag-inertia regime, which is different from the drag-gravity regime occurring at small Reynolds numbers. An isothermal experiment with moderate Reynolds number will be used to validate the WRM instead, in the absence of non-isothermal experiments. The effect of the Reynolds number on the saturation wave will become clear as the Marangoni number is taken as zero.

Moderate Reynolds number
The previous experimental results (Liu & Gollub, 1994) are used to assess the performance of the WRM at moderate Reynolds number. The configuration and parameters adopted here are identical to those in the experiments: Re = 19.33, We = 35, ν = 6.28 × 10 −6 m 2 /s, angle of incline β = 6.4 • , and disturbance frequency f = 1.5 Hz. The inclined angle still lies inside the validity domain of the WRM for β > 1 • . To compare with the experimental results, spatial evolution should be considered. A temporal periodic disturbance in the form of q(0, t) = q N (1 + A) sin(ωt) (q N = 1/3, A = 0.03) is imposed at the inlet.
The spatial evolution of a long wave predicted by the WRM is shown in the right column of Figure 10. It is very similar to that observed in the experiments, as shown in the left column of Figure 10. Due to the convective nature of the instability, the amplitudes of the waves are amplified as the waves travel downstream. It can be observed from Figure 11 that the disturbance waves attain steady traveling permanent waveforms at x ≈ 100 cm (approximately five wavelengths), which were first called 'solitary waves' by Benjamin to distinguish them from solitons in deep-water waves. The shape of the nearly identical solitary waves is constituted by a main hump with a large amplitude preceded by approximately ten small capillary ripples. The solitary waves have an average film thickness of one wavelength, approximately 0.9645, which is slightly smaller than the waveless state. This is consistent with Kapitzas and Alekseenko's finding in their experiments (Alekseenko, Ye, & Pokusaev, 1985;Kapitza & Kapitza, 1949) that the average thickness decreases as the wave travels downstream.
The nearly identical solitary waves shown in Figures 10 and 11 are in a saturation state, which can also be confirmed in corresponding temporal simulation (not shown here) and will be called a 'saturation wave' later. Clearly, the maximum wave amplitude (h ≈ 1.58) of the saturation wave is higher than the experimental one (h ≈ 1.53). The wavelength of the saturation wave is approximately l ≈ 19.4 cm, and the corresponding wavenumber k in terms of the Shkadov scaling is approximately 0.1355. The wave speed of the saturation wave can be calculated immediately with ω/k = 0.1385/0.1355 ≈ 1.02 (in terms of the Shkadov scaling). This is lower than the experimental wave speed (c ≈ 1.09; see Figure 9 in Liu & Gollub, 1994). We can infer from c = L/T that the dimensional wavelength l in the experiment is approximately 20.7 cm, which is longer than the computational wavelength. The errors of the WRM in predicting the wave amplitude, wavelength, and phase speed may result from the departure of the semi-parabolic velocity distribution assumed in the WRM from the actual velocity distribution.
In general, the saturation wave predicted by the WRM at moderate Reynolds number is very similar to the experimental wave. It is reasonable to expect that the saturation wave predicted by the WRM with Re = 4, Figure 10. Comparison between experiment (left) and simulation (right) with Re = 19.33, We = 35, β = 6.4 • , ν = 6.28 × 10 −6 m 2 /s, and f = 1.5 Hz. Reproduced from Physics of Fluids, Vol. 6, Paper No. 1702; doi:10.1063/1.868232, Figure 3, with the permission of AIP Publishing. Ma = 0 will be closer to the exact wave, as the accuracy of the WRM increases as the Reynolds number is decreased for Ma = 0, as shown in Figure 6. Furthermore, the WRM can be valid at Re = 4 if the Marangoni number is sufficiently small since the accuracy of the WRM increases as the Marangoni number is decreased, as shown in Figure 5. Similar conclusions can be drawn for other Reynolds numbers. Thus, it can inferred that the WRM can be applicable for small-to-moderate Reynolds numbers if the Marangoni number is limited to sufficiently small ranges, which depend on the Reynolds numbers.

Summary and conclusions
The validity of the WRM including the Marangoni effect is investigated through linear stability analyses and nonlinear numerical simulations, by which the growth rate and phase speed of the linearly unstable disturbance can be obtained and the interactions among different harmonics excited by the nonlinearity occurring subsequently following the linear instability can be accounted for, respectively. A spatial simulation, which is more suitable for comparison with experiment, is performed by considering a large computational domain with periodic forcing at the inlet and free boundary conditions at the exit, and a temporal simulation, which is more convenient for theoretical analyses, is performed by considering a periodic domain without inlet or exit conditions.
The linear analyses show that the WRM including the Marangoni effect is no longer as accurate as the isothermal analysis, which is applicable for small-to-moderate Reynolds numbers, and may become invalid at approximately Re = 4. This is very important for application of the WRM and has not yet been reported or investigated. The error of the WRM is found to increase nearly linearly with the Marangoni number for each Reynolds number and achieves a maximum at approximately Re = 4 (the error exceeds 20% if Ma = 50), which is nearly independent of the Marangoni number. In other words, the WRM is valid only if the Marangoni number is sufficiently small or if the Reynolds number is far away from four when the Marangoni effect is included.
The error in the WRM caused by the Marangoni effect at small Reynolds number is then investigated through numerical simulation with different Marangoni numbers. The BE is used as the reference model in the absence of corresponding experiments. The numerical simulations show that the error of the WRM in predicting the growth rate and phase speed can become considerable in the spatial simulation if the Marangoni number is sufficiently large (e.g. Ma = 50). Fortunately, the saturation state can be obtained correctly through temporal simulation with the WRM. One may expect that the growth rate and phase speed do not affect the final saturation state, and the temporal simulations with the WRM always predict the saturation states well. This may not be true if the Reynolds number is increased since the flow regime may change into the drag-inertia regime, which is different from the drag-gravity regime occurring at small Reynolds numbers.
An isothermal experiment with moderate Reynolds number (Re = 19.3) is then used to validate the WRM in the absence of non-isothermal experiments. The effect of the Reynolds number on the saturation wave becomes clear if the Marangoni effect is not included. With the experimental parameters and configuration, the spatial simulation with the WRM generates a solitary-wave structure very similar to the experimental wave, except for the wave amplitude being somewhat larger and the wave speed as well as the wavelength being slightly smaller. It can be inferred that the saturation wave predicted by the WRM with Re = 4, Ma = 0 must be more similar to the exact wave since the WRM without the Marangoni effect is more accurate at smaller Reynolds numbers. Furthermore, the WRM can be valid at Re = 4 if the Marangoni number is sufficiently small since the accuracy of the WRM increases as the Marangoni number is decreased. Similar conclusions can be drawn for other Reynolds numbers. Thus, it can be inferred that the WRM predicts the saturation waves well for small-tomoderate Reynolds numbers if the Marangoni numbers are limited to a small range depending on the Reynolds number.
The WRM is accurate as the Reynolds and Marangoni numbers tend to zero, and its accuracy decreases as the Marangoni number is increased. It seems reasonable to infer that the increase of the Marangoni number reshapes the distributions of velocity and temperature within the film layer, because the validity of the WRM depends on the distributions of velocity and temperature, as we have mentioned in the introduction. To improve the accuracy of the WRM, a more accurate assumption of the distributions of velocity and temperature should be adopted in the derivation of the WRM. This is a subject for further studies.