A Potential Field Description for Gravity-Driven Film Flow over PieceWise Planar Topography

Models based on a potential field description and corresponding first integral formulation, embodying a reduction of the associated dynamic boundary condition at a free surface to one of a standard Dirichlet-Neumann type, are used to explore the problem of continuous gravity-driven film flow down an inclined piece-wise planar substrate in the absence of inertia. Numerical solutions of the first integral equations are compared with analytical ones from a linearised form of a reduced equation set resulting from application of the long-wave approximation. The results obtained are shown to: (i) be in very close agreement with existing, comparable experimental data and complementary numerical predictions for isolated step-like topography available in the open literature; (ii) exhibit the same qualitative behaviour for a range of Capillary numbers and step heights/depths, becoming quantitively similar when both are small. A novel outcome of the formulation adopted is identification of an analytic criteria enabling a simple classification procedure for specifying the characteristic nature of the free surface disturbance formed; leading subsequently to the generation of a related, practically relevant, characteristic parameter map in terms of the substrate inclination angle and the Capillary number of the associated flow.


Introduction
Film flow on rigid substrates that are either locally non-planar, due to the presence of isolated topographical features, or exhibit a periodically repeating pattern of the same, has been the topic of considerable research interest over the past four decades or so, motivated by such occurrences in numerous industrially relevant fields.These include, for example, the manufacture of micro-scale sensors and devices [1], thin-film transistors [2], organic light-emitting diode (OLED) displays [3], printed circuits [4], and the deposition of functional coatings using slot-die and slide/cascade arrangements [5].Cumulatively, this has led to a better understanding of important and interrelated phenomena, using gravity-driven film flow in particular as a an investigative platform: ranging from the unwanted free surface disturbances generated by randomly occurring or prescribed micro-scale substrate irregularities, which in turn can become a feature of the final coated product once drying has taken place [6], to the susceptibility of such flows to loss of stability and whether the presence of topography can lead to a delayed onset compared to film flow on completely planar substrate.
Significant contributions to the field experimentally are the early investigations of [7][8][9] and more recently those of [10][11][12] covering both steady films and routes to instability.Much of the associated theoretical and predictive contributions have used models assuming the existence of a steady flow employing the long-wave approximation [13], leading to a lubrication or integral boundary layer approximation, with the corresponding equation sets solved numerically [14].Of note is the combined work of Kalliadasis et al. [15], Mazouchi and Homsy [16] and the comprehensive numerical solutions of Gaskell et al. [17], the results of which being in very good agreement with the experimentally captured free surface data of Decré and Baret [10].In addition, [17] provide a very detailed quantitative parameter study and error analysis of the effect of Reynolds number and flow rate with regard to topography height/depth.The same flow problems were considered a decade later by Veremieiev et al. [18] using a depth averaged form of the NS equations, confirming their findings.Note also that Hayes et al. [19] used a lubrication approximation to model the case of thin film flow over a small particle arrested at a wall, represented as a Dirac function, far exceeding the strict requirements of the underlying assumptions; nevertheless, their results were similarly shown to be realistic when compared against experiment [10].
Few related numerical solutions, of a finite or boundary element type, have appeared focusing in particular on the behaviour of the accompanying 2D flow structure within the film [20] and in the case of 3D flow [21][22][23]; other, related investigations, have used a volume of fluid approach to explore pattern formation and mixing in film flow [24], surface wetting [25], together with the evolution of a single rivulet and its breakup [26,27].
Complementary stability analyses have become increasingly available [28][29][30], revealing interesting features concerning the effect of periodically repeating topography on the stability of films, showing the critical Reynolds number signalling the onset of instability, compared to the corresponding value predicted for film flow on totally planar substrate, is increased.Resulting even in the re-establishment of islands of stability as the Reynolds number is increased further still-in this respect the reader is directed to the following recent contributions [31][32][33][34][35] In connection with the issue of minimising/controlling the free surface disturbance associated with film flow over topography, others have considered the inverse problem of substrate design based on information concerning the free surface, [36][37][38], non-Newtonian effects [39] and the the influence of electric fields [40,41].
The remainder of the manuscript proceeds as follows.The modelling approaches adopted and the corresponding equation sets are developed in Section 2, with solutions obtained in one of two ways, as found in Section 3. Numerically in the case of the first integral equations; analytically by invoking the long-wave approximation which for zero Reynolds number reduces the entire problem to a set of two nonlinear ODEs that are shown to recover the well known degenerate Nusselt solution for film flow on completely planar substrate.Linearisation of these equations to first order, in terms of a sufficiently small value of the ratio of topography height to asymptotic film thickness, as performed in Section 4, leads to analytical solutions for flow over piece-wise features which: (i) compare well to their numerical counterparts and ones in the open literature; (ii) facilitate the construction of a practically useful characteristic parameter map describing the expected nature of the accompanying free surface disturbance.Conclusions are drawn and potential future applications of the method discussed in Section 5.

Problem Specification and Model Formulations
A schematic of the problem of interest is provided in Figure 1; it shows the associated coordinate system and geometry, for the case of film flow, asymptotic thickness H 0 , down a piece-wise planar substrate containing, for illustration purposes only, a sharp-sided trench topography of depth A and length L. For sufficiently large L the problem reduces to one of flow over a step-down followed by a step-up [16,17].This general problem definition similarly applies to the situation of flow over an equal but opposite peak topography and to smoothly varying topographical features.
Schematic of a continuous film of liquid, asymptotic thickness H 0 , flowing down a rigid piece-wise planar substrate, b(X), inclined at angle α to the horizontal and containing, for illustration purposes, a sharp trench topography, height A and length L; shown also is the orientation of the corresponding (X, Y) co-ordinate system, with f (X) denoting the shape of the free surface disturbance experienced by the film.

Standard Formulation
The liquid film is assumed to be Newtonian and incompressible with the flow governed by the corresponding steady-state NS and Continuity equations: ∂U ∂X where U, V and g cos(α), g sin(α) denote the the velocity and body force components parallel and perpendicular to the substrate, respectively; P is the pressure, the density and η the dynamic viscosity.
The problem constitutes a simply connected domain containing a free surface given by X = X(S); that is, the free surface is assumed to be parametrised with respect to its arclength S. In two dimensions its shape is determined by one kinematic ( U • n = 0) and two dynamic boundary conditions: with normal vector n, tangential vector t and σ denoting the surface tension.At the piece-wise planar substrate, whose shape is characterised by Y = b(X), the velocity field must fulfil the no-slip/no-penetration conditions U = 0, V = 0; while, sufficiently far away from the topographical feature of interest, both upstream and downstream, the Nusselt velocity profile for flow over a planar surface is assumed.While the above formulation represents the usual starting point from which to theoretically investigate the flow of interest, as discussed in the Introduction, the thrust of the present contribution is based entirely on a related approach that employs a potential field description leading to a first integral formulation.

First Integral Formulation
In recent years, different first integral forms of the 2D-NS equations have appeared [42][43][44]; the present work employs the real-valued form given in Marner et al. [45]: where Φ enters the problem via the relationship P + U 2 /2 + E = ∆Φ, and E(X, Y) = g [sin(α)X − cos(α)Y] + E 0 (with E 0 a constant that can be chosen arbitrarily) is the potential energy density of the conservative force.The Continuity equation is fulfilled identically by introducing a streamfunction Ψ in accordance with While the substrate boundary conditions are implemented in the usual way, there is an obvious benefit when it comes to applying those at the free surface via a first integral approach: the dynamic boundary conditions there can be reformulated more conveniently, Scholle et al. [44], by eliminating the pressure; leading to an integral formulation and revealing (4) to be a pure condition on the potential gradient, namely: where the entries in the tensor ε are determined by the Levi-Civita symbol, i.e., ε 12 = −ε 21 = 1 and ε 11 = ε 22 = 0, and the integral limit S 0 can be chosen arbitrarily.Furthermore, knowledge of the potential gradient at the free surface allows the derivation of Dirichlet and Neumann boundary conditions for the potential variable in the simply connected domain if Φ is specified at a single point [45].The kinematic boundary condition can be formulated as: the unknown constant Ψ s indicating the flow rate and the free surface shape characterised in terms of the as yet unknown function f (X) via Y = f (X).After decomposition into component parts the first integral (7) of the dynamic boundary condition takes the form: with E defined in terms of the integrals appearing in (9,10) as: where the constants of integration E 1 , E 2 can be chosen arbitrarily in the same way as E 0 .
For convenience E 0 = − g cos αH 0 is chosen; then E 1 = 0 and E 2 = g cos αH 2 0 /2.Subsequently, use is made of the following scalings.Taking H 0 as a characteristic length scale defines the corresponding non-dimensional quantities: while a characteristic velocity scale based on a value equal to twice the corresponding free surface Nusselt velocity, denoted by U N = gH 2 0 sin(α)/η, leads to the following scalings of the velocities (U, V), potential field Φ, the streamfuction Ψ and the potential energy density E(X, Y): The Reynolds and Capillary number is given by Re = U N H 0 /η and Ca = ηU N /σ, respectively.

Lubrication Model
The first integral Equations ( 5) and ( 6) can be simplified based on a lubrication approximation, noting that the same applied to the NS equations, Stillwagon and Larson [8], leads to a single equation for the local film thickness; the velocity field to leading order being locally parabolic.The requirement underpinning its applicability in the case of topography exhibiting step-like changes is that the commensurate free surface disturbance is slowly varying which is indeed the case; film flow over step-like topographical features have been investigated successfully by many authors within a lubrication framework, as highlighted in the Introduction. Neglecting and omitting the inertial term, 2 UV, in Equation ( 6), results in the following simplified equation: which can be integrated with respect to Y to give: the function F 1 (X) appearing as a consequence.By inserting (13) into (5) and neglecting, as above, the inertial term, U 2 − V 2 , leads to the following relationship which when integrated with respect to Y, gives: with appearance of the function F 2 (X).Taking the derivative of ( 13) with respect to Y and that of (14) with respect to X, and combining the two, ∂ Y (13)−∂ X (14), eliminates the potential Φ, leading to: Furthermore, the assumption that ∂ 2 X ψ ∂ 2 Y ψ allows integration of the above expression, resulting in the locally parabolic velocity profile: and when integrated for a second time the cubic form: of the streamfunction; introducing two additional functions F 3 (X) and F 4 (X).Inserting the solution ( 16) into ( 13) and ( 14) gives: for the gradient of the potential.
The above lubrication model (16)(17)(18) in terms of Ψ and Φ contains four unknown functions F 1 (X), • • • , F 4 (X); additionally the shape of the free surface, given by the function f (X), remains an unknown.The five boundary conditions required to close the problem are dealt with in detail in Appendix A.1; the resulting set of five ODEs can be reduced by means of successive elimination of unknowns, as shown in Appendix A.2, resulting finally in just two nonlinear ODEs (A10, A11) for the functions F 1 (X) and f (X).
Making use of the scalings identified earlier, such that b(x) = b(X)/H 0 , together with the following substitutions: which, motivated by the degenerate, A = 0, Nusselt solution (A12, A13) for Ψ s and F 1 (provided for completeness in Appendix A.3 and in which case both ψ * and F(x) vanish), can reasonably be considered small for film flow over a piece-wise planar substrate.Finally, the non-dimensional free surface profile and local film thickness are given by: Considering all of the above, the following non-dimensional form of the equation set (A10, A11) is obtained: Consequently, the entire problem has been reduced to a set of just two nonlinear ODEs (23,24), which can be solved either numerically as they are or, as is the method of choice derived in Section 3 following further approximation, analytically.

Methods of Solution
Two methods of solution are adopted: the first a general numerical approximation based on a FE formation of Equations ( 5) and ( 6); the second analytical, based on the above lubrication analysis, providing simple and elegant asymptotics.

FE Analysis
The least squares FEM developed in Marner et al. [45], Marner [46] is employed, based on the first integral formulation of Section 2 and facilitating the use of simple equal order elements together with highly efficient multigrid solvers; for a comprehensive review of least squares methods, including a special treatment of the Stokes and NS equations, the reader is referred to [47].
Practicality suggests rewriting Equations ( 5) and ( 6) in terms of velocity variables u 1 = ∂ψ/∂y, u 2 = −∂ψ/∂x and first order derivatives of the potential field φ, i.e., φ 1 = ∂φ/∂x and φ 2 = ∂φ/∂y, leading to a system of four equations involving first order derivatives only, which is consistent with a first order system least squares methodology.This allows for the use of standard piece-wise polynomial H 1 -elements [47].The introduction of new field variables requires a corresponding prescription of the standard Continuity equation and a further condition for the potential fields, namely: In terms of the scalings defined earlier, Equations ( 5), ( 6) and ( 25) can be rewritten in the following non-dimensional and condensed first-order form: the associated differential operators of which are defined as: with the solution vector given by: The operator L embodies the terms associated with the linear Stokes flow case with N [u] containing terms related to inertial effects, while writing the vector of field variables in the form (28) facilitates the use of the compact notation (26).
The free surface shape is obtained by iterating over the kinematic condition while solving a sequence of flow problems with prescribed dynamic conditions in a fixed domain; the free surface approximation is updated after each iteration step.Such a procedure avoids the necessity of incorporating the corresponding geometric degrees of freedom into the system matrix, as in the case of the spine method [48], which would result in a more complicated implementation process and destroy the favourable structure of the matrix system.The linear systems in each iteration step are constructed via standard FE concepts and solved efficiently via a customised multigrid procedure [46].
Using the above notation and the formula for the dynamic boundary condition ( 7), the non-dimensionalised film flow boundary value problem at each iteration step has the form (i, j = 1, 2): with a semi-parabolic Nusselt velocity profile prescribed (at an appropriate distance from the topographical feature) at the left-and right-hand sides of the flow domain.Note that, via this arrangement Dirichlet boundary conditions are prescribed everywhere, for either both velocity components or both potential fields.For convenience the Dirichlet conditions can be written in the condensed form Bu = g u [46].Since the FE results generated in Section 4 are merely for comparison and verification purposes against corresponding lubrication solutions, Re → 0, subsequent description of the methodology is confined to the case of Stokes flow only.The linear boundary-value problem was solved via the following minimisation procedure: with • 0 denoting the standard L 2 -norm and involving the function spaces: where Ω and ∂Ω denote the flow domain and its entire boundary, respectively.The equations comprising the differential system are multiplied by weighting factors, hence the presence of the weighting operator W, which is actually a diagonal matrix with positive entries √ w i and w 1 , . . . ,w 4 ∈ R >0 .This provides the option to give prominence to selected equations, for example the Continuity equation with the objective of improving mass conservation; typically Equations ( 25) are weighted by a factor of 10 4 against 1 for the other two equations.It is important to note that such a weighting scenario does not affect the principal convergence properties of the method which are shown to be optimal, see [46].The minimisation problem, (30), is equivalent to the weak formulation: and a conformal discrete version is simply obtained by replacing the function spaces U and V by finite dimensional subspaces U h ⊂ U and V h ⊂ V.In the present work, piece-wise quadratic Lagrange polynomials on a triangulated domain were used for all test and solution spaces; typically resulting in about 50,000 basis functions per unknown.The system matrix was constructed following an iso-parametric FE-concept allowing triangles at a boundary with a curved side, with the associated integrals evaluated via standard Gaussian quadrature schemes.The free surface iteration process was begun from an initial shape obtained from a corresponding analytical solution-the details of which are provided in the following subsection-and updated iteratively via the kinematic boundary condition.Given that an approximate solution u (n) (x, y) together with an approximate free surface shape f (n) (x) depending on x = x(s) is available, successive free surface shapes, f (n+1) (x), are found by solving a differential equation of the form: by a high-order Runge-Kutta scheme, in which case the streamline starting from (x 0 , y 0 ) is traced.
Since the current free surface shape f (n) and corresponding velocity field u are approximations only, it may happen that a streamline f (n+1) starting from (x 0 , h 0 ) lies partly outside of the current mesh (flow domain); this is the reason why h 0 /2 < y 0 < h 0 is chosen in a way that ensures the streamline f (n+1) remains inside the meshed domain and reaches the outflow boundary.If this is achieved, a new approximation of the free surface shape is obtained by shifting the calculated streamline f (n+1) by an amount h 0 − y 0 to the correct position.Using f (n+1) as the new free surface mesh boundary, the solution u (n+1) is derived from the weak formulation (33) and once again a free surface approximation is traced, and so on.
The starting value y 0 for the streamline algorithm is adaptive.If the initial guess for the free surface is poor, it might be necessary to choose y 0 significantly smaller than h 0 ; however, in the course of the iteration process the difference h 0 − y 0 is continually decreased until it converges to zero.As a convergence criterion for the whole iteration process the fractional difference between the flow in the normal direction at the free surface and that in the tangential direction is integrated along the free surface and the value required to have changed by no more than 0.01% between iterations.Note that, an appropriate damping methodology, depending on the value of the Capillary number and the whole iteration history, is essential in terms of the stability and speed of convergence of the iterative process [45,46].

Linear Approximation
The approach developed in Section 2.2.1 allows for a simple treatment of the flows of interest, in that since the functions F (x) and ϕ (x) = 0 can be assumed small, for b (x) small, then in the spirit of a lubrication approximation products of their derivatives can be considered negligible compared to all terms in Equations (23,24) containing the product F (x)ϕ (x).Taking the linear combination 2d(x)×(24)−3× (23) and solving for F(x) gives: facilitating the elimination of F(x) from ( 23), requiring the second derivative of ( 35) to be found by making use of the following relationship: leading to: Finally, by writing d(x) = 1 + ϕ(x) − b(x) according to (22), including inserting (36) into ( 23), the following nonlinear third order ODE: is obtained for the surface shape function ϕ(x).Since b(x) is zero beyond the immediate vicinity of the topography and ϕ(x) vanishes for x → ±∞, it follows from (37) that ψ * = 0.The nonlinear ODE (37) applies to any film for arbitrary a = A/H 0 and free surface disturbances away from H 0 .If the film is thin compared to A, both b(x) and ϕ(x) can reach large values requiring a numerical solution; however, if A and the associated free surface disturbances are small compared to H 0 , then b(x) and ϕ(x) both remain small and Equation (37) can be linearised by neglecting quadratic and higher order terms with respect to ϕ(x) and b(x), as employed in the analysis of Hayes et al. [19], resulting in the following linearised ODE: Equation ( 38) has constant coefficients, in contrast to Equation (37), and can be solved analytically for any locally defined piece-wise planar topographical feature given by b(x).

Free Surface Disturbance Types
Before generating solutions for film flow over isolated piece-wise, step-like topography the above linear approximation can be used to identify the type of free surface disturbance that can arise: damped oscillatory behaviour synonymous with a pronounced capillary ridge/trough, or simply a smooth monotonic transition.For this purpose it is sufficient to focus on the homogeneous solution of (38), i.e., for b(x) = 0, requiring the zeros of the associated characteristic polynomial to be determined: Based on the detailed analysis provided in Appendix A.4 three modes are discernible, depending on the sign of the discriminant D defined according to (A18).These are: Monotonic transition If D < 0, three real-valued roots (A19-A21) result with κ 1,2 > 0 and κ 3 < 0, giving rise to the general solution: as a superposition of pure exponential functions, indicating an exponential like spatial decay of the free surface disturbance, with decay length H 0 /|κ 3 | and H 0 / min{κ 1 , κ 2 } downstream and upstream of the topography, respectively.Damped oscillations If D > 0 , two of the roots, (A19, A20), are complex-valued with positive real part with one the complex conjugate of the other, i.e., κ 1,2 = κ 0 ± ik with κ 0 > 0, while the third root, (A21), remains a negative real number, i.e., κ 3 < 0. Although by analytic continuation the general solution ( 40) is also valid in this case, it is useful to rearrange it using Euler's formula exp(κ 1,2 x) = exp(κ 0 x) [cos (kx) ± i sin (kx)], to give: where ϕ ± are given by: Unlike the previous mode, spatial oscillations occur upstream of the topography that are exponentially damped with decay length H 0 /κ 0 ; while a monotonic exponential transition downstream of the topography persists.Aperiodic limit The limit case indicated by a vanishing discriminant, D = 0, requires separate treatment, which is elegantly achieved by applying the limit k → 0 to (41), giving: implying again a monotonic transition of the free surface disturbance spanning the upstream and downstream regions of the topography.

Film Flow over Isolated Step-Like Topographical Features
Taking the case of a piece-wise step-down topography as an exemplar, defined as: the corresponding solution of equation ( 40) for a free surface disturbance which decays monotonically is given by: fulfilling the asymptotic conditions automatically.The remaining three constants ϕ 1 , ϕ 2 , ϕ 3 have to be determined by the matching condition that ϕ, ϕ and ϕ have to be continuous, leading to the following system of linear equations: the solution of which reads: The same solution is also valid for the case of damped spatial oscillations but can be more conveniently expressed according to (40) as: where κ 0 = (κ 1 + κ 2 )/2, k = i(κ 2 − κ 1 )/2 and The aperiodic borderline is obtained by applying the limit k → 0 to (50-53), leading to the non-oscillatory solution: with coefficients: (57) The above procedure can also be applied to a step-up topography, defined as: the solution for which is obtained from Equations ( 45) and ( 46) for a step-down by the inversion operation a → −a.Similarly solutions can be formulated for linear superpositions of both solutions to form a trench or peak feature.

Green's Function and Linear Solution for Arbitrary Topography
For linear problems a Green's function methodology provides an elegant way to determine the explicit solution of the respective differential equation for arbitrary substrate topography via a convolution integral.In this resect, consider the solution ϕ θ (x) obtained in Section 3.2.2 for a step topography taking a = −1, with corresponding substrate profile b(x) = θ(x) where θ(x) denotes the Heaviside step function.Depending on the the parameter regime, ϕ θ (x) takes the form ( 46), ( 50) or (54); its derivative, G(x) := ϕ θ (x), is the solution for b(x) = θ (x) = δ(x) where δ(x) denotes the Dirac delta function.Since any arbitrary topography b(x) fulfils the identity b the free surface shape corresponding to b is given by the convolution integral: This elegant method only works for the linearised theory; the nonlinear ODE (37) would need to be solved for different topography b(x) individually.Taking, as an example, film flow over a trench topography defined as: where as above a = A/H 0 and in addition = L/H 0 denoting the aspect ratio of the trench, the solution for the free surface shape, (59), is given by: a linear combination of the solution, −ϕ θ , for a step-down at x = 0 with one for a step-up topography, +ϕ θ , at x = .

Step-Up and Step-Down Topography
Initially the linear approximation to the lubrication model derived in Section 3.2 was validated against existing experimental data and related numerical predictions available in the open literature, namely those of Decré and Baret [10] and Gaskell et al. [17], respectively, for film flow over step-up and step-down topography.The former still remain pretty much unique in the context of step-like topography while the latter have become a benchmark for such flows.
Step-down and step-up topography is considered for the same set of parameter values and fluid properties: α = 30°, a = 0.2 and Ca = 7 × 10 −4 .Figure 2 shows the free-surface disturbance predicted via equation (41) to be in excellent agreement for both configurations: the figures in the left column compare the linear solution (blue curve) against their experimental counterpart (dashed curve), while those in the right column compare the linear solution against both experiment and the numerical prediction (black curve) of [17].Given the related experimental error involved [10] all three profiles can be said to agree extremely well.Figure 3 shows the predicted free surface disturbance resulting from film flow down substrate containing a step-down topography, inclined at angle α = 60°to the horizontal, for step heights a = 1/2, 1/4 and 1/32 for different values of Ca (12.0, 0.774 and 0.1).The entire flow regime, obtained using the linear model, is shown in the left column for a = 1/2 and 1/4; while in the right column only the predicted free surface profile (both analytic and FE) is presented for all three step heights.Using A as a alternative scaling length in the vertical direction, the analytically predicted profile becomes independent of a due to the linearity of the analytical approach, while the FEM solution of the full nonlinear problem clearly reveals a dependence of the surface profile on a, manifesting itself as a shift of the free surface in the upstream direction with decreasing a; such that for very small a the numerical solution approaches the analytical one.In such cases, agreement between the analytical and numerical solution is better still the smaller the Capillary number, since due to the stronger influence of surface tension the free surface disturbance is more gradual (smoother); which is a key requirement of the lubrication assumption.
For small surface tension (Ca = 12), case (a), the free surface height transitions monotonically from upstream to downstream.In contrast, for a Capillary number below a critical value for the chosen inclination angle (cases (b) and (c)), on the upstream side of the step the free surface experiences a maximum value, above that of the Nusselt thickness, before recovering to attain this asymptotic thickness at a distance upstream of the step dependent on the value of Ca.As expected, this capillary ridge is more pronounced for smaller values of Ca and is shifted further in the upstream direction.Film flow over a step-up feature exhibits qualitatively similar behaviour but in this case, and in contrast, a capillary trough forms on the upstream side of the step.Equivalent plots to those for a step-down are shown in Figure 4.A trough is clearly visible for sufficiently small Capillary numbers and comparison of the analytical solution with the accompanying numerical results reveals once again a systematic shift of the film surface in one direction, this time downstream with decreasing a; the correspondence between the two becoming better for smaller a, and as observed in the case of a step-down the agreement between the two solutions is better still the smaller the Capillary number.Otherwise, overall good agreement between the analytical and numerical predictions is observed.

Trench Topography
In the left column of Figure 5 resulting free surface predictions, (61), for film flow over a trench topography, (60), are shown for the case α = 80°, with a = 1/2, Ca = 0.3 and for different trench widths, .These are compared to corresponding FE solutions in the right column, the latter revealing also the internal flow structure in the form of streamlines.
In all cases, in the region spanning the side walls bounding the trench a distinct depression in the free surface is exhibited, the minimum of which shifts towards the step-up side of the trench with increasing , while on the upstream of the step-down a capillary ridge forms that overshoots the Nusselt film thickness before relaxing to this value further upstream.In the region 0 ≤ x ≤ , for = 2 the local free surface position remains everywhere above a value equal to the unperturbed asymptotic film thickness, denoted by the horizontal dashed line; while for the wide trench, = 6, the minimum film thickness undershoots the value of H 0 in the middle of the tench.
The complementary FE results, shown in the right column of Figure 5, confirm these observations qualitatively revealing slight quantitative differences; the reason for the latter is arguably related to the vortical structures, or eddies, present in the corner regions of the trench formed by the side walls.It can be seen that the eddies are larger and extend further up the side of the trench, the smaller the value of .Since they occupy the space available in the trench proportionally, they shift the overlying flow and hence the free surface position in the vicinity of the trench, a feature that is not captured, for obvious reasons, using a lubrication approach.The analytically predicted local film thickness, depicted by the red curve, shows a broadening of the the free surface depression compared to the corresponding FE free surface streamline (the upper most blue curve).Nevertheless the expectation, based on the earlier investigations involving individual step-down and step-up features, is that agreement between the two will improve for smaller values of the depth parameter a.This is demonstrated quantitatively by a comparative study in terms of the predicted free surface profiles presented in Figure 6 based on, as used previously, an alternative scaling with respect to A and in which case the analytically obtained curves are independent of a.
Taking the third case presented in Figure 5 with a trench width = 6, capillary number Ca = 0.3, inclination angle α = 80°and a = 1/2, as a starting point, reveals the broader analytically calculated free surface depression compared to the numerically computed one, which leads also to an overestimation of the maximum depression in the middle of the trench.However, it can be seen that as a is reduced, first to a = 1/4 and then a = 1/32, the quantitative difference between the predicted free surface profiles diminishes considerably, leading to very good agreement between two.

Capillary Criterion and Characteristic Free Surface Parameter Map
The dominant free surface disturbance (a capillary ridge or trough) present in the predictions of Figures 3-6 for different values of α and Ca, raises the question as to whether a simple criterion can be found for identifying the appearance of such capillary-induced phenomena.As discussed earlier, the roots of the characteristic polynomial (39) are indicative in this respect: capillary ridges/troughs in tandem with a series of smaller strongly decaying oscillations occur if and only if two of its roots are complex, which is equivalent to inequality D(α, Ca) > 0 for the discriminant (A18).Since D depends on α and Ca only, a characteristic parameter map in terms of these variables can be generated, as shown in the upper plot of .Upper plot: characteristic parameter map for α against Ca identifying the regions associated with (i) a monotonic transition of the free surface disturbance (blue region) and (ii) the occurrence of spatial oscillations with a dominant capillary ridge/trough (green region); the two regions being separated by a neutral limit curve along which the discriminant (A18) vanishes.The different step-up/-down cases, (a)-(c), explored in Figures 3 and 4 are identified by the labels a-c, respectively, the trench case considered in Figures 5 and 6 by the label t.Lower plot: Neutral limit curve (blue curve) compared with the corresponding curves obtained using the approximations (A26) and (A27), derived in Appendix A.5, when α is small (red curve) and close to 90°(green curve), respectively.
Accordingly, the regime that persists for a given α and increasing Ca, or vice versa, for the case of film flow over a piece-wise, step-like topography can be identified.The three different flow cases, (a)-(c), shown in Figures 3 and 4 are identified on the map with the corresponding labels a-c; the point with the label t denoting the corresponding solution for the trench topography explored in Figures 5  and 6.For a given α a capillary ridge/trough is only present below a critical Ca, the value of which is implicitly given by D(α, Ca) = 0, becoming larger with increasing α.The free surface profile inserts obtained analytically for the case of film flow over a step-down topography have been included to illustrate the distinction between the green region denoting the expected presence of a capillary ridge/trough, compared to the blue region signifying a monotonic transition; the locus of the curve D(α, Ca) = 0 depicts the neutral, or aperiodic, limit separating the two regions.A striking feature is that even for a vertically aligned film (α = 90°) a finite critical capillary number Ca ⊥ can be identified by solving D(90°, Ca ⊥ ) = 0, giving Ca ⊥ = 9 √ 2 ≈ 12.73.For any capillary number above this value, the free surface can only transition monotonically in the direction of flow.
In Appendix A.5 explicit approximation formulae for the neutral limit curve for small inclination angles and those close to 90°are provided; the curves resulting from the latter are shown in the lower plot of Figure 7 in relation to the corresponding exact neutral limit curve.

Conclusions
It is shown that starting from a first integral of the NS equations, written in terms of the streamfunction Ψ and a potential Φ, is beneficial for investigating, from both a numerical and analytical standpoint, gravity-driven film flow down an inclined substrate featuring piece-wise planar step-like and trench topography.Using a least square FEM, the field equations and related boundary conditions can be solved without restrictions on Reynolds or Capillary number and substrate geometry.In tandem, a lubrication approximation of the first integral equations facilitates the modelling of such flows for zero Reynolds number when the height/depth of the topography is small compared to the overlying Nusselt film thickness.Reducing the problem to one requiring the solution of two nonlinear ODEs for the film thickness, which is shown to reproduce the degenerate Nusselt solution for film flow on a uniformly planar surface.Subsequent linearisation of these equations enables solutions to be found via a standard analytical approach, for the topographies of interest which once derived take a matter of seconds to obtain.For the case of step-down and step-up features they are shown to agree almost identically, for the same geometry and fluid properties, with the detailed experimental results of Decré and Baret [10] and the corresponding numerical predictions of Gaskell et al. [17].
Further detailed comparative studies concerning the free surface disturbance generated by step-like and trench topographies, including for the latter visualisation of the internal flow structure, confirm the closeness between corresponding analytic and FE predictions.As expected, better consistency between the two is reached for small values of topography height/depth compared to the asymptotic film thickness and when the value of the Capillary number is small.In all of the cases considered the analytical approach is found to provide an acceptable prediction of the resulting free surface disturbance.
A key result emerging from the adopted first first integral formulation, is the generation of a practically useful characteristic parameter map in terms of inclination angle of the substrate and the Capillary number, as demonstrated in Section 4.3.It allows instantaneous insight as to the behaviour of the free surface disturbance induced by capillary effects associated with a particular topographical feature: namely whether the feature gives rise to spacial oscillations exhibiting a dominant capillary ridge/trough or experiences a simple monotonic transition.A maximum value for the Capillary number, Ca ⊥ = 9 √ 2 ≈ 12.73, is identified above which the free surface disturbance will always transition monotonically.For Capillary number values below this upper bound the presence of spacial oscillations or otherwise, depends on the inclination angle only, with a clear demarkation between the two possibilities represented by a neutral limit curve.It is shown that approximate neutral limit curves can be derived which are in excellent agreement with the master curve both when Ca is small and tends towards the value Ca ⊥ .In some respects this parameter map plays a similar role to a typical stability map in which a curve of neutral stability separates regions of stable from an unstable film flow [30], and indeed from a wider perceptive an analogy with the problem of a buckling beam suggests itself.
Looking ahead, there is considerable scope for developing the methodology beyond the essentially linear approach used in the present work; numerical solution of the nonlinear ODE, (37), derived for the more general problem, offers the prospect of more accurate predictions of the free surface disturbance associated with film flow over topography other than simple piece-wise planar ones.The flexibility of the approach is such that it can be readily extended to include periodically repeating surface features, sharp or smooth, and to encompass the problem of film flow on curved 2D surfaces.Moreover, the recent 3D generalisation of the first integral approach reported by Scholle et al. [49] opens up the prospect of investigating film flows over 3D topographical features in a similar manner, again both analytically and numerically, and for extension to coating flows involving smooth 3D curved surfaces in general-for example film flow on a hemisphere as investigated by Takagi and Huppert [50].

Figure 2 .
Figure 2. Comparison of the free surface disturbance predicted by (41) (blue curve) against the experimental data (dashed curve) of [10] (left column) as well as the numerical solution (black curve) of Gaskell et al. [17] (right column), both for a step-down (top row) and a step-up (bottom row) topography: the flow is from left to right with α = 30°, and Ca = 7 × 10 −4 and the dotted line shows the associated nature of the topography for comparison purposes.

Figure 3 .
Figure 3. Linearised lubrication model solutions (left column) for film flow over a step-down (a = 1/2, 1/4 and 1/32 for different values of Ca (12.0 (top), 0.7774 (middle), 0.1 (bottom)) and α = 60°) showing the predicted free surface disturbance (blue curve), compared with FE solutions of the first integral equations (right column; using a different scaling for clarity): (a) monotonic transition; (b,c) damped spatial oscillations upstream.The flow is from left to right, the dashed lines in the left-hand column indicating the unperturbed film thickness on the upstream and downstream sides of the step.

Figure 4 . 6 Figure 5 .Figure 6 .
Figure 4. Linearised lubrication model solutions (left column) for film flow over a step-up (a = 1/2, 1/4 and 1/32 for different values of Ca (12.0 (top), 0.7774 (middle), 0.1 (bottom)) and α = 60°) showing the predicted free surface disturbance (blue curve), compared with FE solutions of the first integral equations (right column; using a different scaling for clarity): (a) monotonic transition; (b,c) damped spatial oscillations upstream.The flow is from left to right, the dashed lines in the left-hand column indicating the unperturbed film thickness on the upstream and downstream sides of the step.

Figure 7
Figure 7. Upper plot: characteristic parameter map for α against Ca identifying the regions associated with (i) a monotonic transition of the free surface disturbance (blue region) and (ii) the occurrence of spatial oscillations with a dominant capillary ridge/trough (green region); the two regions being separated by a neutral limit curve along which the discriminant (A18) vanishes.The different step-up/-down cases, (a)-(c), explored in Figures3 and 4are identified by the labels a-c, respectively, the trench case considered in Figures5 and 6by the label t.Lower plot: Neutral limit curve (blue curve) compared with the corresponding curves obtained using the approximations (A26) and (A27), derived in Appendix A.5, when α is small (red curve) and close to 90°(green curve), respectively.