A CONSERVATIVE SAINT-VENANT TYPE MODEL TO DESCRIBE THE DYNAMICS OF THIN PARTIALLY WETTING FILMS WITH REGULARIZED FORCES AT THE CONTACT LINE

This paper deals with the numerical simulation of thin liquid films flowing on partially wetting solid substrates. A 2D Saint-Venant like model is proposed. Its originality lies in the conservative formulation of the capillary forces and in the model used for the disjoining pressure that accounts for the contact line capillary forces. A finite volume scheme is proposed for the resolution of the system and various numerical examples are presented and discussed. In particular, when the mesh resolution is fine enough, the model is proved to be able to predict correctly the spreading of a film with the exact contact angle in the vicinity of the contact line. When the mesh size is larger than the film thickness (which could be the case for many industrial applications), it is of course no longer possible to recover the contact angle. However, the model is proved to correctly predict the spreading of the film. This important feature is related to the thermodynamic consistency of the model in the sense that the latter ensures by construction the decrease of the film total free energy in the absence of external driving forces.


Introduction
The motion and stability of liquid thin films and droplets which wet a solid substrate are present in a lot of natural and industrial processes and have been the object of a lot of research studies for several decades [8,13,33]. For instance, in the context of in-flight aircraft deicing, when a thermal protection system is activated, the supercooled water droplets impacting onto an aircraft surface do not freeze instantaneously and can coalesce and form a thin liquid film as a result of aerodynamic forces. Experimental studies show that this liquid film is not always stable and can split into rivulets that may refreeze on unprotected surfaces [38]. The modeling of rivulet flows and the accurate prediction of wet and dry surfaces is important since it has a direct influence on the exchanges between the film, the wall and the airflow boundary layer.
A classical approach for the modeling of thin liquid films is based on the averaging of the Navier-Stokes (N-S) equations over the film thickness h under the assumption that the ratio between h and any in-plane length-scale is small (long-wave approximation). The reduction of the N-S equations can lead to the single non-linear fourth-order partial differential lubrication equation for h [2], or to the Saint-Venant third-order system composed of two equations for h and the flow mass rate q [32]. The basic formulation of these models is limited to fully wetting films which spread indefinitely without contact line. Film rupture to rivulets or isolated droplets through the destabilizing molecular forces is not included [12,20,27,30,34,35].
To take into account the long-range forces in the vicinity of the contact line with possible film breaking into rivulets and partial wetting, an additional pressure term introduced by Frumkin and Derjaguin [4,5,11,15], called the disjoining pressure, is included into the lubrication equation. The disjoining pressure Π d (h) is an additional term arising from the interaction between two interfaces separated by a distance h. The first interface is the one between the liquid film and the solid substrate. The second interface is the one between the liquid film and the gaseous boundary layer above it. The disjoining pressure Π d (h) can be classified into several interactions like the van der Waals interactions, the electrical double layer or the polar interaction. However, to investigate the behavior of a macroscopic film, Π d (h) only has to verify the augmented Young Laplace equation [21]. Following this interpretation, alternative phenomenological models for the disjoining pressure have appeared in the literature like a diffuse interface model by Thiele & Pomeau that is based on exponential functions [26,36], power-law expressions [16,19,29], or combinations of power-law and exponential expressions [1,31]. These models are based on a disjoining pressure defined between two unbounded parallel interfaces. In reality, and especially near the contact line, the interfaces are generally not flat. Dai [7] and Hocking [17] incorporated the gradient and curvature of the liquid film into the expression of the disjoining pressure.
Another difficulty is the hydrodynamic contact line singularity in the configuration of partial wetting [18,24]. To prevent viscous dissipation divergence for vanishing film thicknesses, a first solution is the add of a slip length that allows the fluid to slip at the boundary [3,14,28]. Some authors rather use a disjoining pressure combining long-range repulsive and short-range attractive forces leading to an energetically stable precursor film of thickness h prec h [1,16,19,29,31]. The latter, which is present everywhere on the substrate actually acts as a slip length.
This paper is the first of a series of two articles dedicated to the 3D macroscopic numerical simulation of the transition from a continuous liquid film to rivulets. Having in mind future industrial applications, in particular in the field of ice protection system (IPS) numerical simulation, this involves the use of complex general geometries with any unstructured surface grids. To our knowledge, only the lubrication equation has been generalized to partial wetting by adding a disjoining pressure to the pressure term [4,5,11,15] with a 4 th order capillary term. The objective of this paper is to propose a general formulation for partially wetting films adapted to a classical 2 nd order finite volume space discretization on any grids. This involves a conservative formulation of the equations with the writing of an augmented system [25]. Particular attention is paid to obtaining the most generic model possible, verifying the existence of a mathematical entropy having here the meaning of a free energy. A second requirement which is related to the industrial context is the need to propose a simplified model for the capillary effects. Indeed, a description at molecular scales is not possible. The model must account for the macroscopic properties of a partially wetting film without solving the smallest scales. This paper is focused on the general equations in 2D with the presentation of the model and its numerical discretization. A forthcoming paper in 3D, based on the same approach, will complement the results already presented in [22].
The model is established in the first part. The equivalent Saint-Venant type system is derived as well as its associated energy equation. The second part is dedicated to the numerical method. The last part of the paper is dedicated to numerical applications of the model for some academic but generic test cases. The influence of the regularizing length scale involved in the definition of the model chosen for the disjoing pressure is discussed and best practises are derived depending on the mesh resolution.
1. Model derivation based on the lubrication equation

Lubrication equation for partially wetting films
We consider an incompressible liquid film moving on a planar solid surface. The substrate curvature effects will not be addressed in this paper. The x and z axis are respectively parallel and perpendicular to the substrate. The interface between the air gas and the liquid surface corresponds to z = h(x, t), where h is the liquid film thickness. The velocity field in the liquid film is written u (x, z, t). The liquid film thickness is assumed to be small (ε = ho Lo 1) with h o and L o respectively the characteristic length scales related to the liquid film thickness and spreading. The inertial forces inside the liquid film are supposed to be small or at most comparable to the viscous ones leading to a Reynolds number lower than 1 (Re 1). Depending on the targeted applications, the Froude number is supposed to be large for air-sheared liquid film (F r 1) and of the order of 1 for falling films (F r ∼ 1). Both Reynolds and Froude numbers are computed from h 0 and a characteristic velocity scale u 0 of the film motion. Under these assumptions, the liquid film motion can be described by the lubrication equation [2]: where the density of the liquid film ρ is a positive constant that does not depend on space and u h is a function of the unknown h defined by: µ is the dynamic viscosity of the liquid film, P a stands for gas pressure above the film, ρg n h is the hydrostatic pressure, γ lg κ is the Laplace pressure, with κ being the curvature of the gas/film interface and γ lg the liquid/gas surface tension coefficient. κ is given by: g n and g t are respectively the components of the gravity acceleration g projected onto the z axis and the x axis. τ a is the shear stress of the outer air (x component). The term Π d (h) represents an additional pressure term that models the interactions between the molecules of the fluid and the substrate when the liquid film thickness h tends toward zero [11]. It is the cause of partial wetting effects (dewetting and film rupture phenomena). The so-called disjoining pressure Π d (h) can be derived from the disjoining energy e d including long-range contributions [11] by: The expression of e d (h) will be given later in Sec. 1.4. In Eq. (2), u h is the z-average of the velocity profile of a 2D Poiseuille flow u h reading: By construction, u h satisfies: and: b is a slip length ( Fig. 1) preventing the divergence of the viscous dissipation rate for vanishing film thicknesses [18]. From [9], b can be physically related to a characteristic molecular length scale. Typically, it can be assumed that b ∼ 10 −9 . This is this value that has been used for all the simulations presented in Sec. 3. It is worth mentioning that for such a small value of b, no influence is expected on the numerical results since it is too (Eq. (7)).
small compared to the film thickness even in the cells located in the vicinity of the contact line. Equation (7) may be rewritten:

Equivalent Saint-Venant model
Back to Eqs. (1) and (2), the term corresponding to the capillary forces (γ lg κ) is of the fourth order in space. To avoid numerical difficulties due to the spatial discretization of such a high order term (especially regarding the 3D extension of the model with unstructured grids on 2D surfaces), it may be interesting to reduce it by two orders in space. A first step to do this is to replace Eqs. (1) and (2) by the following third-order Saint-Venant system, where the unknowns are now h and q = hu: Similarly to the lubrication equation (Eq. (1)), the closure formulation for the wall shear stress τ w is obtained by assuming a Poiseuille flow inside the liquid film (Eq. (8)): In Eq. (9b), the closure hypothesis u 2 = u 2 has been done for the convective term. Although this assumption is not consistent with a Poiseuille profile inside the film, it has no consequences on the validity of our model under the assumptions ε = ho Lo 1 and Re 1. Indeed, Eq. (9b) is equivalent to: where u h is defined by Eq.(2) and T ν = ρh(h + b)/3µ is the viscous time scale. In our study, the characteristic time scale at which physical phenomena are observed is the convective time scale T conv = L o /u o . From the hypothesis ε = ho Lo 1 and Re 1, the ratio between T ν and T conv is given by: Hence, it can be assumed that the solution of Eq. (11) satisfies u ≈ u h . This proves the formal equivalence between the proposed Saint-Venant system (9) and the lubrication equation (1)+(2) at the relevant time scale for our targeted applications. For the rest of the manuscript, u will be written simply u to simplify the notations. In Eq. (9b), the pressure terms (except P a ) can be re-written in a conservative way as follows: This leads to an equivalent conservative formulation of the Saint-Venant system (9) given by: where F is given by : The term ∂F ∂x in Eq. (14b) is still of third order in space. The spatial order of the system can be reduced yet by introducing the new independant unknown p = (∂h/∂x) [25]. By spatially deriving Eq. (14a), the following equation for p is obtained: Adding this equation to system (14) results in the following second-order in space system: where the unknowns are now ρh, p and ρq.

Energy equation associated to the Saint-Venant model
From physical arguments, the film free energy e per unit surface area is given by: where e k , e p and e c , respectively the kinetic, potential and capillary energy densities, are defined as: The first (resp. second) term in the definition of e c stands for the liquid/gas (resp. liquid/substrate) interfacial energy. The third one is detailed in Sec. 1.4. It stands for the disjoining energy associated to partially wetting effects.
Defining the following variable ϕ by: the flux function F (Eq. (15)) can be rewritten: with L the Legendre transformation of e defined as: where U = (ρh, p, ρq) T . System (17) can then be rewritten: where S = h ρg t − ∂Pa ∂x − τ w + τ a represents the source terms and A is the following skew symmetric matrix: Multiplying Eq. (23) by ϕ and using Eq. (21) and the skew symmetric property of A, we get, after some algebra: This equation describes the decrease of the free energy of the film in the absence of external driving forces. The term Φ = u · (e + L) − ϕ · A ∂ϕ ∂x corresponds to the energy flux. The dissipation term reads: System (17) is thus thermodynamically consistent. It is worth mentioning that if the disjoining energy e d is a convex function of h (which is not the case for a partially wetting film as shown in Sec. 1.4), e can be proved to be a convex function of the conservative variables ρh, p, ρq, and thus plays the role of a Lax entropy [23] for Eqs. (17). This implies in particular that in this case system Eqs. (17) has a hyperbolic structure compatible with the second-order terms. However this mathematical structure is no longer present in the case where the e d function is concave.

Model for the disjoining energy e d
The objective is to model the competition among the following surface energies, each of them being associated with an interface between two phases: • The surface energy between the solid substrate and the liquid film γ sl , • The surface energy between the liquid film and the gas flow above it γ lg , • The surface energy between the solid substrate and the gas flow above it γ sg .
This balance drives the behavior at the contact line, in particular the wetting and unwetting properties of the film. To ensure that e d (h) is able to address partial wetting situations with both wet and dry areas, it is necessary for the film energy formulation (Eq. (18)) to degenerate into the one of a dry substrate when h = 0 and into the one of a wet substrate when h R where R is the characteristic length scale of the long-range molecular forces. In this paper, only macroscopic scales will be considered. It is then recommended to define e d as simple and regular as possible so that the following conditions are verified: • For dry conditions, the only surface energy involved is the one at the substrate/gas interface γ sg . This corresponds to x = 0 in Fig. 2 where h(0) = 0 and p(0) = dh dx x=0 = 0. Therefore, from Eqs. (18) and (19): which leads to: where S is the so-called spreading parameter (here negative since the substrate is supposed to be partially wetting) and θ s is the static contact angle derived from γ sl , γ sg , γ lg and the Young-Dupré equation [8].
• For wet conditions such as h R: As it will be shown in Sec. 2, the two conditions written for (h, p) = (0, 0) and h R are relevant to address partial wetting at a macroscopic scale. On the other hand, these conditions alone are not accurate enough if physics has to be described at the very microscopic level.
For S < 0 and for partially wetting configurations, the disjoining energy e d must be a concave function [10]. In the present work, we propose to use the following simple expression for e d that meets the previous boundary conditions (Eqs. (28) and (29)): As will be shown, this simple smooth model for the disjoining energy allows to recover the macroscopic behavior of the liquid film in the vicinity of the contact line. Indeed, provided a concave shape for e d , only its asymptotic behavior matters but not its own expression. Long-range forces located at the contact line are simply regularized over the characteristic length h which will play the role of a "numerical" scale for the long-range molecular forces that does not need to be as small as the real physical scale to get the correct macroscopic behavior.
Recommended best practices regarding the choice of h are proposed in Sec. 3.

Macroscopic consistency of the model
The objective here is to highlight the model capability to reproduce the typical physical behaviors observed in the vicinity of the contact line in both static and dynamic regimes. A volume of liquid film near the contact line is considered (in light blue in Fig. 2). Its borders are defined by: • The solid substrate, • The gas phase, The point x where the thickness h(x) and the angle θ(x) are defined. Using Eq. (17a), Eq. (17c) can be rewritten: in the absence of external driving forces and where D· Dt stands for the total derivative. The volume of liquid film is assumed to advance with the constant velocity of the contact line u CL [8] so that the previous equation may be rewritten: Equation (32) is integrated from the positions ξ = 0 to ξ = x, leading to: After some developments, some terms of the previous equation can be simplified as follows: • From Eq. (30), it is obtained: where • The boundary condition (Eq. (28)) for the disjoining energy gives: where dη = tan(θ)dξ, the angle θ being defined in Fig. 2. Equation (33) can be rewritten in a dimensionless form: where C a = µu CL γ lg is the capillary number, l c = γ lg ρg n is the capillary length and the function G is defined as: Equation (38) is now written for the central region where R h(x) l c (Fig. 2). In that region, the dimension of the liquid film is small when compared to the capillary length l c so that the gravitational forces can then be neglected [8]. Figure 3 shows a steep decrease of the function G and the associated long-range can then be neglected in the central region. Regarding the curvature effects, the hypothesis of a simple wedge shape in the central region is done [9], leading to κ(x) ∼ 0. Therefore, in the central region, Eq. (38) can be rewritten as: Regarding the integral contribution and provided the simple wedge shape hypothesis in the central region [9], the angle θ is constant in the central region. Therefore, one obtains: On the other hand, near the contact line, as a consequence of the boundary conditions h(0) = 0 and dh dx x=0 = tan(θ(0)) = 0, one has tan(θ) < tan(θ(x)) for 0 < h(ξ) < R, leading to: So, We can thus introduce a macroscopic length scale H(x) associated to the film thickness in the central region and defined as follows: By definition, because of Eq. (43), one has H(x) > h(x). Equation (40) can then be rewritten: leading finally to: The logarithmic behavior of the right-hand side can been found in all dynamic contact line models and is in particular in line with [9] for the case of small angles (θ, θ s 1). For a receding contact line (C a > 0 in Fig. 2), cos (θ s ) < cos (θ(x)) and θ < θ s whereas for an advancing contact line (C a < 0 in Fig. 2), cos (θ s ) > cos (θ(x)) and θ > θ s . This is consistent with the experimental observations. If the film is considered to be at rest (C a = 0), Eq. (46) leads to: cos (θ(x)) = cos (θ s ) (47) which is consistent with the Young-Dupré equation.

Numerical methods
This section is devoted to the presentation of the finite volume numerical scheme used to perform the simulations of Sec. 3 and to some theoretical considerations regarding its consistency with the energy balance equation (Eq. (25)). The proposed scheme is only of first order in space and time but could be easliy extended to second order by using the classical MUSCL approach.

Discretization scheme
For the sake of simplicity, let us first consider the spatial discretization of system (23) only and let us restrict ourselves to a uniform mesh of size ∆x. The finite volume method, which consists of integrating system (23) over a given cell i, leads to the following semi-discrete system: where U i and S i are respectively the averaged values of U and (0, 0, S) T over the cell i, and F i+1/2 and B i+1/2 are respectively the numerical fluxes at interface i + 1/2 between cells i and i + 1 that respectively correspond to the terms u U + (0, 0, L) T and A ∂ϕ ∂x . Regarding first F i+1/2 , we have chosen the following scheme: with u + i+1/2 (resp. u − i+1/2 ) being respectively the positive (resp. negative) part of u i+1/2 defined as u i+1/2 = (u i + u i+1 ) /2 and: This corresponds to an upwind scheme for the convective part of the flux and a centered scheme for the force term L. Regarding now the term B i+1/2 , the scheme used is a centered scheme that reads: For the time discretization, we used the implicit Euler method in order to avoid any resctrictive stability condition on the time step (at least from heuristic considerations). This finally leads to the following numerical scheme:

Semi-discrete energy equation
Let us now consider the problem of deriving a discrete counterpart of the continuous energy equation (Eq. (25)). Taking the scalar product of Eq. (48) with ϕ i (t) leads to: Using the expression of F i+1/2 and the definition of ϕ i we get: Then, using that L i = ϕ i · U i − e i (by definition of L), we get after some algebra: where: The skew-symmetry of A i+1/2 leads to: In the absence of external driving forces: Gathering all the intermediate results, we finally obtain the following local semi-discrete energy inequality (in the absence of external driving forces): where: Summing Eq. (61) over all cells, we finally get: with E = i∈Z e i ∆x. The decrease of the total energy at the discrete level is therefore guaranteed if all the terms D i are non-negative. By definition of D i , this property is verified if and only if the film energy density e is a convex function of U (i.e. of h, p and ρq). It can be easily verified that e k , e p and e c − e d are convex functions of U . So, for fully wetting films (e d ≡ 0), the proposed semi-discrete scheme ensures the decrease of the film total energy. But, for partially wetting films, the disjoining energy e d being intrinsically a concave function of h, D i can not be proved to be non-negative for all conditions. Thus, for partially wetting films, unlike the case of the continuous formulation, the semi-discrete scheme (Eq. (48)) cannot be proved to be consistent with the decrease of the film energy. This question will have to be further investigated in the future but is out of the scope of the present paper. Finally, it is worth mentionning that in practice (using an implicit euler scheme for the time discretization), the monotonous decrease of the film energy has always been verified for all the applications presented in the next section.

Numerical results and recommended best practices regarding the choice of h
In the definition of the disjoining energy e d (Eq. (30) On the other hand, if the mesh size can be chosen such that ∆x h ref , then it will be shown that the model allows to compute both the contact angle and the spreading of the film provided that h satisfies ∆x h h ref . In the test cases presented bellow, we focus on a film moving on an horizontal flat plate (g t = 0) in the absence of outer air flow ( ∂Pa ∂x and τ a = 0). Two configurations are studied: the steady state of the stabilization of a puddle subject to gravity, and the contact line dynamics of a spreading film in the absence of gravity.

Static test-case
The equilibrium shape of a 2D heavy droplet deposited on a horizontal plate and subject to gravity is studied. The initial radius of the droplet is supposed to be large enough compared to the capillary length l c for the equilibrium shape to be that of a puddle [10]. The thickness profile h i of the initial droplet is given by: where R i = 0.025m. The initial contact angle θ i between the liquid phase and the solid wall is 15 • . The properties of the liquid droplet are similar to those of silicone oil (Tab. 1).
The equilibrium (Fig. 4) is characterized by the puddle thickness h 0 far from the contact line and given theoretically by [10]: Table 1. Properties of the liquid droplet.   Table 2. Geometrical parameters of the final steady puddle (Fig. 4)

initialized from Eq. (65)
The volume V of the droplet is obtained from the integration of Eq. (65): The spreading W 0 is then defined by: The numerical values for l c , h 0 and V are given in Tab. 2. An analysis is proposed at two different space scales. Firstly, provided that both scales h 0 and h are accurately resolved by the mesh size ∆x, the influence of h on the solution is studied. Both global (like the computed puddle thickness) and local (like the angle θ at the contact line) quantities will be compared with the theoretical values h 0 and θ s respectively. Then, and this is the configuration encountered in practice with a film whose spreading is large compared to its thickness, the influence of h will be studied when h and even h 0 can no longer be resolved spatially. The ability to estimate macroscopic quantities such as spreading will then be investigated.
As an example, time evolution from the initial droplet (Eq. (65)) to the final puddle (Fig. 4) is represented in Fig. 5 with h /h 0 = 1/100 and ∆x/h = 0.25. The equilibrium shape of a static puddle is reached at t = 1s. To ensure stationarity, all the computations have been conducted up to at least 10s.   (Fig. 4). h /h 0 = 1/100 and ∆x/h = 0.25.
value h 0 (Eq. (66)) by computing: The influence of h on ε h0 is shown in Fig. 6. We can observe the existence of an optimal value for the ratio h /∆x. As expected for very large h /∆x (and therefore large h /h 0 ), the numerical error increases. Indeed, increasing h leads to an increase in the thickness over which the long-range forces are regularized around the contact line. For large values of h , the local nature of the long-range forces is then poorly addressed by the model. More surprinsingly, the numerical error increases as well for small h /∆x. This is due to the fact that for h ∆x, the long-range forces are under-resolved, which leads to a non-exact prediction of the static contact angle and in turn of the puddle thickness. This is confirmed in Fig. 7 where a focus on the contact line is proposed. The grey dash-dot line is the reference slope at the contact line derived from the theoretical static contact angle θ s . The dashed lines stand for a poor spatial discretization of the long-range forces (h /∆x 1) whereas the solid lines represent an accurate spatial discretization of the long-range forces (h /∆x 1). Once again, only the optimal configurations (lowest h /h 0 ratio while h /∆x 1) allow an accurate computation of the static contact angle in the vicinity of the contact line (here h /h 0 = 1/50 and h /h 0 = 1/100 in Fig. 7). The most accurate solution h /h 0 = 1/100 (black solid line in Fig. 7) is now chosen as the reference solution. A mesh convergence study is proposed for the most accurate solution h /h 0 = 1/100 in Fig. 8. To do so, ∆x ranges between h /8 and 2h . The following errors are defined: where h # α and p # α refer to the numerical solutions of h and p computed on the grid defined by ∆x = αh with α ∈ [1/8, . . . , 2]. Regarding ε h , h # 1/8 stands for the solution on the finest grid (∆x = h /8). First order is obtained for both ε h and ε p , which is line with the first-order accurate numerical schemes used for space discretization. Figure 8b also shows the consistency of the model for p. Indeed, although p is considered as an independent variable, the relation p = dh dx is well verified at discrete level.  The challenge here is to compute accurately the spreading W 0 of the liquid film, without necessarily accurately capturing the contact angle. The thickness of the initial droplet is given by Eq. (65) with R i = 0.25m and θ i = 15 • . The properties of the liquid droplet are given in Tab. 1. Its volume is given by V =110.1cm 2 . The puddle thickness remains unchanged with h 0 =1.749mm (Tab. 2). This leads to a spreading W 0 =6.3m. Therefore, a clear space scale separation is observed between the thickness and the spreading of the puddle. The computational domain is discretized into 200 cells with L 0 =8m leading to a space step ∆x = 0.04m. The ratios ∆x/W 0 ≈ 0.006 and ∆x/h 0 ≈ 23 indicate an accurate space resolution for W 0 and a poor one for h 0 .
The influence of h on the final shape of the puddle is represented in Fig. 10.     lines), the long-range forces are regularized over a thickness of the order of h 0 (and even larger), leading to a poor estimation of the puddle spreading. On the other hand, when h /h 0 1 (solid lines), the spreading is accurately computed from x = −W 0 /2 to x = W 0 /2. The smaller the ratio h /h 0 , the better the prediction on the spreading. However, the lack of space resolution (∆x h 0 ) prevents any accurate reproduction of the static contact angle near the contact line. Care must be taken when h becomes very small, in the order of the thickness h num of a very thin and non-physical residual film that may remain on the wall due to the artificial diffusion of the first-order numerical scheme. The numerical scale of the long-range forces being h , that may lead to a film that spuriously spreads infinitely. Therefore, h must remain larger than h num .
Two main conclusions can be drawn from this first test case based on the computation of a static puddle. The first point is that, whatever the space resolution, the condition h /h 0 1 must be verified to accurately compute macroscopic properties like puddle spreading (h /h 0 < 1/10 seems to be sufficient in practice). However, and this is the second point, for an accurate description of the local contact angle in the vicinity of the contact line, the additional condition ∆x h is necessary.

Dynamic case: contact line spreading in the absence of gravity
The objective here is to numerically check the validity of the formal calculation of Sec. 1.5 and to study the influence of h on the dynamic behavior of the spreading of the contact line in the absence of gravity and for a substrate of low static contact angle (θ s = 4 • ). Focus is put on Eq. (46) where the constant angle θ(x) in the central region is written θ m . It can be seen as the macroscopic angle observed near the contact line. Equation   Table 4. Simple linear regression between θ 3 s − θ 3 m and C a . Slope of the regression lines: α. Root mean-squared error: ε α . The influence of h is studied. h /∆x > 1.
whereH is the characteristic thickness of the liquid film in the central region and C a = µu CL γ lg is the dimensionless capillary number with u CL the velocity of the contact line. Under the hypothesis that θ s 1 and θ m 1 and given thatH b, the previous equation leads to: which has been introduced by de Gennes [9]. As a reminder, the similar Cox-Voinov equation [6,37] is given by: The numerical test case consists in a droplet of initial thickness h i given by Eq. (65), with R i = 500µm, θ i = 50 • , leading to V = 162mm 3 (Eq. (67)). The properties of the liquid droplet are given in Tab. 3. Unlike a totally wetting film that spreads indefinitely over the solid substrate, a partially wetting film has a parabolic equilibrium profile (Eq. (65)) whose radius R s is given by R s = D s /2 = 1.865mm (obtained from volume conservation, Eq. (67)). The maximum thickness h s is then given by h s = 65.1µm (Eq. (65)). In this section, subscript 's' refers to the asymptotic steady state. The computational domain, set to [−L 0 /2; L 0 /2] with L 0 = 0.0044m, is discretized into 4400 cells leading to a space step ∆x = 1.0 µm, which allows h /∆x > 1 for any value of the investigated ratio h /h s (Tab. 4). As an example, time evolution of the shape of the droplet is represented in Fig. 11 for h /h s = 1/33. At t = 5s, the steady state is reached.
The model capability to reproduce the linear dependence between θ 3 s − θ 3 m and C a is now investigated at the numerical level. To do so, the simple linear regressions θ 3 s − θ 3 m = α · C a are computed where α is the slope of the regression lines. The influence of h /h s on the root mean-squared error ε α of the linear regression is shown in Tab. 4. The linear dependency of θ 3 s − θ 3 m on C a is all the stronger when the ratio h /h s is low. Figure 12  The influence of h /h s on the time evolution of the contact line position x CL (Fig. 13) starts to become significant at late times for the "long-range + capillary" regime. For the "visco-capillary" regime, the influence of h on the dynamics of the contact line is weak. To clarify the concepts of "visco-capillary" and "long-range + capillary" regimes, the different contributions of the force balance in Eq. (31) are integrated in space along Figure 11. Contact line spreading in the absence of gravity. Time evolution from the initial droplet. h /h s = 1/33, h /∆x > 1.
the half droplet between x = −L/2 and x = 0 for h /h s = 1/33 (Fig. 14). The "visco-capillary" regime is characterized by the predominance of the capillary (red curve) and viscous terms (green curve) with the regular decrease of the inertial term (black curve). For the "long-range + capillary" regime, both the capillary and long-range (blue curve) terms are predominant.
In conclusion, film dynamics is accurately accounted for. In the case of low inertial regimes corresponding to predominant contribution of the long-range forces, small h /h s ratios are needed.

Conclusions
In this article, we have proposed a simple Saint-Venant type model to predict the movement, under lubrication regime, of a two-dimensional, partially wetting, thin liquid film on a solid substrate. This model is written in conservative form and does not involve partial differential operators of order greater than two, so it is well adapted to finite volume discretization schemes. Moreover, we have proved that this model is thermodynamically consistent in the sense that it leads to a decrease of the film free energy in the absence of external driving forces. With regard to the long-range molecular forces acting in the vicinity of the contact line, these are taken into account by introducing a disjoining energy whose simplified expression is in our case a function of a single parameter h acting as a regulation length scale. We have shown that if we are only interested in the macroscopic behaviour of the film (spreading and apparent contact angle if the mesh size is fine enough), the model gives very good results with h much larger than the physical scale of the disjoining forces, which allows us to consider in the long term the application of the model for industrial configurations. Ongoing work is focused on taking into account the contact angle hysteresis and on the generalization of the model in the case of three-dimensional films. Encouraging results have already been obtained [21] and have been published in [22]. The main difficulty remains the construction of robust and precise numerical schemes on any mesh, and whose compatibility with the decrease of the free energy can be proven.   31)) for h /h s = 1/33. The black dashed vertical lines (t = 4.0 × 10 −3 s and t = 0.2s) outline the boundaries between the "visco-capillary" and "longrange + capillary" regimes. h /∆x > 1.