P Trontin
email: pierre.trontin@onera.fr

J Lallement
email: julien.lallement@onera.fr

P Villedieu
email: philippe.villedieu@onera.fr

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 [START_REF] De Gennes | Wetting: statics and dynamics[END_REF][START_REF] Snoeijer | Moving contact lines: Scales, regimes and transitions[END_REF][START_REF] Dussan | On the spreading of liquids on solid surfaces: static and dynamic contact lines[END_REF]. 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 [START_REF] Zhang | An experimental investigation on the surface water transport process over an airfoil by using a digital image projection technique[END_REF]. 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 [START_REF] Benney | Long waves on liquid films[END_REF], or to the Saint-Venant third-order system composed of two equations for h and the flow mass rate q [START_REF] Shkadov | Wave flow regimes of a thin layer of viscous fluid subject to gravity[END_REF]. 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 [START_REF] Thiele | Dewetting of an evaporating thin liquid film: Heterogeneous nucleation and surface instability[END_REF][START_REF] Thiele | On the importance of nucleation solutions for the rupture of thin liquid films[END_REF][START_REF] Diez | On the breakup of fluid films of finite and infinite extent[END_REF][START_REF] Kheshgi | Dewetting: Nucleation and growth of dry regions[END_REF][START_REF] Reiter | Dewetting of thin polymer films[END_REF][START_REF] Seemann | Dewetting patterns and molecular forces: A reconciliation[END_REF].

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 [START_REF] Frumkin | On wetting and sticking phenomena[END_REF][START_REF] Derjaguin | Tiyra Kapillyarnoy Kondesatsii and Drugiz Kapillapnvix Yavlenii Uchetom Rasklinivayushchevo Daystviya Polimolekuyarnox Shidi Plenox (phonetic translation of the original Russian title)[END_REF][START_REF] Churaev | Wetting films and wetting[END_REF][START_REF] Churaev | Prediction of contact angles on the basis of the frumkin-derjaguin approach[END_REF], 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 [START_REF] Lallement | Modélisation et simulation numérique d'écoulements de films minces avec effet de mouillage partiel[END_REF]. 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 [START_REF] Pismen | Disjoining potential and spreading of thin liquid layers in the diffuse-interface model coupled to hydrodynamics[END_REF][START_REF] Thiele | Film rupture in the diffuse interface model coupled to hydrodynamics[END_REF], power-law expressions [START_REF] Schwartz | On the shapes of droplets that are sliding on a vertical wall[END_REF][START_REF] Glasner | Spreading of droplets under the influence of intermolecular forces[END_REF][START_REF] Khanna | Pattern formation in spontaneous dewetting of thin apolar films[END_REF], or combinations of power-law and exponential expressions [START_REF] Sharma | Relationship of thin film stability and morphology to macroscopic parameters of wetting in the apolar and polar systems[END_REF][START_REF] Beltrame | Rayleigh and depinning instabilities of forced liquid ridges on heterogeneous substrates[END_REF]. 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 [START_REF] Dai | Disjoining pressure for nonuniform thin films[END_REF] and Hocking [START_REF] Hocking | The influence of intermolecular forces on thin fluid layers[END_REF] 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 [START_REF] Huh | Hydrodynamic model of steady movement of a solid/liquid/fluid contact line[END_REF][START_REF] Moffatt | Viscous and resistive eddies near a sharp corner[END_REF]. 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 [START_REF] Dussan | The moving contact line: the slip boundary condition[END_REF][START_REF] Bertozzi | Linear stability and transient growth in driven contact lines[END_REF][START_REF] Roux | Modèles de fronts pour films minces[END_REF]. 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 [START_REF] Schwartz | On the shapes of droplets that are sliding on a vertical wall[END_REF][START_REF] Glasner | Spreading of droplets under the influence of intermolecular forces[END_REF][START_REF] Khanna | Pattern formation in spontaneous dewetting of thin apolar films[END_REF][START_REF] Sharma | Relationship of thin film stability and morphology to macroscopic parameters of wetting in the apolar and polar systems[END_REF][START_REF] Beltrame | Rayleigh and depinning instabilities of forced liquid ridges on heterogeneous substrates[END_REF]. 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 [START_REF] Frumkin | On wetting and sticking phenomena[END_REF][START_REF] Derjaguin | Tiyra Kapillyarnoy Kondesatsii and Drugiz Kapillapnvix Yavlenii Uchetom Rasklinivayushchevo Daystviya Polimolekuyarnox Shidi Plenox (phonetic translation of the original Russian title)[END_REF][START_REF] Churaev | Wetting films and wetting[END_REF][START_REF] Churaev | Prediction of contact angles on the basis of the frumkin-derjaguin approach[END_REF] 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 [START_REF] Noble | Stability theory for difference approximations of some dispersive shallow water equations and application to thin film flows[END_REF]. 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 [START_REF] Lallement | A shallow water type model to describe the dynamic of thin partially wetting films for the simulation of anti-icing systems[END_REF].

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.

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 [START_REF] Benney | Long waves on liquid films[END_REF]:

∂(ρh) ∂t + ∂ (ρhu h ) ∂x = 0 (1)
where u h is a function of the unknown h defined by:

u h (x) = - h(h + b) 3µ
∂ ∂x P a + ρg n h + γ lg κ -Π d (h) -ρg t - 3τ a 2h h + 2b/3 h + b (2)
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:

κ(x) = - ∂ ∂x   ∂h ∂x 1 + ∂h ∂x 2   (3)
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 [START_REF] Derjaguin | Tiyra Kapillyarnoy Kondesatsii and Drugiz Kapillapnvix Yavlenii Uchetom Rasklinivayushchevo Daystviya Polimolekuyarnox Shidi Plenox (phonetic translation of the original Russian title)[END_REF]. 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 [START_REF] Derjaguin | Tiyra Kapillyarnoy Kondesatsii and Drugiz Kapillapnvix Yavlenii Uchetom Rasklinivayushchevo Daystviya Polimolekuyarnox Shidi Plenox (phonetic translation of the original Russian title)[END_REF] by:

Π d (h) = - de d dh (h) (4)
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:

u h (x, z) = - 1 µ ∂ ∂x P a + ρg n h + γ lg κ -Π d (h) -ρg t h(z + b/3) - z 2 2 + τ a µ (z + b/3) (5)
By construction, u h satisfies:

µ ∂u h ∂z x,z=h = τ a (x) (6)
and:

µ ∂u h ∂z x,z=0 = µ u h (x, 0) b/3 (7)
b is a slip length (Fig. 1) preventing the divergence of the viscous dissipation rate for vanishing film thicknesses [START_REF] Huh | Hydrodynamic model of steady movement of a solid/liquid/fluid contact line[END_REF]. From [START_REF] De Gennes | Dynamics of wetting: local contact angles[END_REF], 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 small compared to the film thickness even in the cells located in the vicinity of the contact line. Equation ( 7) may be rewritten: (Eq. ( 7)).

µ ∂u h ∂z x,z=0 = 3µu h h + b - τ a 2 h h + b (8)

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:

∂(ρh) ∂t + ∂ (ρhu) ∂x = 0 (9a) ∂(ρq) ∂t + ∂ ∂x ρhu 2 + h ∂ ∂x P a + ρg n h + γ lg κ -Π d (h) = ρhg t -τ w + τ a (9b)
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)):

τ w = 3µu h + b - τ a 2 h h + b (10)
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:

du dt = ∂u ∂t + u ∂u ∂x = - u -u h T ν (11)
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:

T ν /T conv = 1 3 h + b h ε Re 1 (12)
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:

h ∂ ∂x (ρg n h) = ∂ ∂x ρg n h 2 2 (13a) h ∂ ∂x (γ lg κ) = ∂ ∂x (hγ lg κ) + γ lg ∂h ∂x ∂ ∂x   ∂h ∂x 1 + ∂h ∂x 2   (13b) = ∂ ∂x (hγ lg κ) -γ lg ∂ ∂x   1 1 + ∂h ∂x 2   (13c) h ∂ ∂x -Π d (h) = ∂ ∂x -hΠ d (h) + Π d (h) ∂h ∂x (13d) = ∂ ∂x h ∂e d ∂h (h) - ∂e d ∂h (h) ∂h ∂x (13e) = ∂ ∂x h ∂e d ∂h (h) - ∂e d ∂x (13f)
This leads to an equivalent conservative formulation of the Saint-Venant system [START_REF] Diez | On the breakup of fluid films of finite and infinite extent[END_REF] given by:

∂ (ρh) ∂t + ∂ (ρuh) ∂x = 0 (14a) ∂(ρq) ∂t + ∂ ∂x (ρuq) + ∂ ∂x F h, ∂h ∂x , ∂ 2 h ∂x 2 = h ρg t - ∂P a ∂x -τ w + τ a ( 14b
)
where F is given by :

F h, ∂h ∂x , ∂ 2 h ∂x 2 = ρg n h 2 2 + h ∂e d ∂h (h) -e d (h) - γ lg 1 + ∂h ∂x 2 + γ lg hκ (15)
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) [START_REF] Noble | Stability theory for difference approximations of some dispersive shallow water equations and application to thin film flows[END_REF]. By spatially deriving Eq. (14a), the following equation for p is obtained:

∂p ∂t + ∂ (u • p) ∂x = - ∂ ∂x h ∂u ∂x (16)
Adding this equation to system [START_REF] Derjaguin | Tiyra Kapillyarnoy Kondesatsii and Drugiz Kapillapnvix Yavlenii Uchetom Rasklinivayushchevo Daystviya Polimolekuyarnox Shidi Plenox (phonetic translation of the original Russian title)[END_REF] results in the following second-order in space system:

∂ (ρh) ∂t + ∂ (ρuh) ∂x = 0 (17a) ∂p ∂t + ∂ (u • p) ∂x = - ∂ ∂x h ∂u ∂x (17b) ∂ (ρq) ∂t + ∂ (ρuq) ∂x + ∂ ∂x F h, p, ∂p ∂x = h ρg t - ∂P a ∂x -τ w + τ a ( 17c
)
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: e (ρh, p, ρq) = e k (ρh, ρq) + e p (ρh) + e c (ρh, p)

where e k , e p and e c , respectively the kinetic, potential and capillary energy densities, are defined as:

e k (ρh, ρq) = ρhu 2 2 = (ρq) 2 2ρh e p (ρh) = ρg n h 2 2 = g n (ρh) 2 2ρ e c (ρh, p) = γ lg 1 + p 2 + γ sl + e d (h) = γ lg 1 + p 2 + γ sl + e d ( ρh ρ ) (19)
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:

ϕ = (ϕ h , ϕ p , ϕ q ) T = ∂e ∂ρh , ∂e ∂p , ∂e ∂ρq T =    g n h -u 2 2 + 1 ρ ∂e d ∂h (h) γ lg p √ 1+p 2 u    (20)
the flux function F (Eq. ( 15)) can be rewritten:

F = L -h ∂ϕ p ∂x (21)
with L the Legendre transformation of e defined as:

L = ϕ • U -e (22)
where U = (ρh, p, ρq) T . System [START_REF] Lallement | Modélisation et simulation numérique d'écoulements de films minces avec effet de mouillage partiel[END_REF] can then be rewritten:

∂U ∂t + ∂ (u U ) ∂x + ∂ ∂x   0 0 L   = ∂ ∂x A ∂ϕ ∂x +   0 0 S   (23)
where S = h ρg t -∂Pa ∂xτ w + τ a represents the source terms and A is the following skew symmetric matrix:

A =   0 0 0 0 0 -h 0 h 0   (24)
Multiplying Eq. ( 23) by ϕ and using Eq. ( 21) and the skew symmetric property of A, we get, after some algebra:

∂e ∂t + ∂ ∂x (u • (e + L)) - ∂ ∂x ϕ • A ∂ϕ ∂x = u • S (25)
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:

u • S = -u • τ w = - 3µu 2 h + b ≤ 0 (26)
System ( 17) is thus thermodynamically consistent and the energy e formally plays the role of a mathematical entropy function (in the Lax sense [START_REF] Lax | Hyperbolic systems of conservation laws ii[END_REF]) for Eqs. [START_REF] Lallement | Modélisation et simulation numérique d'écoulements de films minces avec effet de mouillage partiel[END_REF].

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):

γ sg = e f ilm h = 0 p = 0 = γ sl + γ lg + e d (0) (27)
which leads to:

e d (0) = γ sg -γ sl -γ lg = S = γ lg (cos θ s -1) ( 28
)
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 [START_REF] De Gennes | Wetting: statics and dynamics[END_REF].

• For wet conditions such as h R:

e d (h R) = 0 and e = γ sl + γ lg (29)
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 [START_REF] De Gennes | Capillarity and Wetting Phenomena: Drops, Bubbles, Pearls, Waves[END_REF]. 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)):

e d (h) = S exp (-h/h ) = γ lg (cos θ s -1) exp (-h/h ) (30)
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. 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 = 0 where h(0) = 0 and p(0) = dh dx x=0 = 0, • The point x where the thickness h(x) and the angle θ(x) are defined.

Using Eq. (17a), Eq. (17c) can be rewritten:

ρh Du Dt + ∂ ∂x F h, p, ∂p ∂x = - 3µu h + b (31)
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 [START_REF] De Gennes | Wetting: statics and dynamics[END_REF] so that the previous equation may be rewritten:

∂ ∂x F h, p, ∂p ∂x + 3µu CL h + b = 0 (32)
Equation ( 32) is integrated from the positions ξ = 0 to ξ = x, leading to:

ρg n h(x) 2 2 + h(x) ∂e d ∂h (h(x)) -e d (h(x)) + γ lg h(x)κ(x) - γ lg 1 + p 2 + (e d (0) + γ lg ) + x 0 3µu CL h(ξ) + b dξ = 0 (33)
After some developments, some terms of the previous equation can be simplified as follows:

• γ lg 1 + p(x) 2 = γ lg 1 + tan 2 (θ(x)) = γ lg cos (θ(x)) (34)
• From Eq. ( 30), it is obtained:

h(x) ∂e d ∂h (h(x)) -e d (h(x)) = -S • [1 + h(x)/h ] • exp [-h(x)/h ] (35)
where S = γ lg (cos θ s -1) • The boundary condition (Eq. ( 28)) for the disjoining energy gives:

e d (0) + γ lg = γ lg cos (θ s ) (36
θ(x) h(x) l c 0 h(x) ∼ R h(x) R u CL Figure 2.
Vicinity of the contact line for a liquid film at rest. In a similar way to [START_REF] De Gennes | Dynamics of wetting: local contact angles[END_REF], four regions are defined: the molecular region, not addressed in this manuscript, where the continuum description breaks down, the proximal region (h(x) ∼ R) controlled by long range forces, the central region (h(x) R and h(x) l c ) controlled by capillarity and the distal region (h(x) l c ) controlled by gravitational forces. The liquid film is assumed to advance with the constant velocity of the contact line u CL .

• x 0 1 h(ξ) + b dξ = h(x) 0 1 tan(θ) dη η + b (37)
where dη = tan(θ)dξ, the angle θ being defined in Fig. 2.

Equation ( 33) can be rewritten in a dimensionless form:

1 2 h(x) l c 2 + [1 -cos (θ s )] • G h(x) h + h(x)κ(x) + [cos (θ s ) -cos (θ(x))] + 3C a h(x) 0 1 tan(θ) dη η + b = 0 ( 38
)
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:

G (ψ) = (1 + ψ) • exp (-ψ) (39)
Equation ( 38) is now written for the central region where R h(x) l c (Fig. 2). The gravitational forces can then be neglected. Figure 3 shows a steep decrease of the function G and the associated long-range forces for h(x)/h 1. Assuming R ∼ h , the contribution G h(x) h 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 [START_REF] De Gennes | Dynamics of wetting: local contact angles[END_REF], leading to κ(x) ∼ 0. Therefore, in the central region, Eq. ( 38) can be rewritten as:

cos (θ s ) -cos (θ(x)) + 3C a h(x) 0 1 tan(θ) 1 η + b I(η,b) dη = 0 (40)
Regarding the integral contribution and provided the simple wedge shape hypothesis in the central region [START_REF] De Gennes | Dynamics of wetting: local contact angles[END_REF], the angle θ is constant in the central region. Therefore, one obtains:

I(η, b) = 1 tan(θ(x)) 1 η + b if R < η < h(x) (central region) (41)
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:

I(η, b) > 1 tan(θ(x)) 1 η + b if 0 < η < R (proximal region) (42) So, h(x) 0 1 tan(θ) dη η + b > 1 tan(θ(x)) h(x) 0 dη η + b (43)
We can thus introduce a macroscopic length scale H(x) associated to the film thickness in the central region and defined as follows: 1 tan(θ(x))

H(x) 0 dη η + b = h(x) 0 1 tan(θ) dη η + b (44)
By definition, because of Eq. ( 43), one has H(x) > h(x). Equation (40) can then be rewritten:

cos (θ s ) -cos (θ(x)) + 3C a tan(θ(x)) H(x) 0 dη η + b = 0 (45)
leading finally to:

tan (θ(x)) • [cos (θ s ) -cos (θ(x))] = -3C a ln 1 + H(x) b if R h(x) l c (46)
The logarithmic behavior of the right-hand side can been found in all dynamic contact line models and is in particular in line with [START_REF] De Gennes | Dynamics of wetting: local contact angles[END_REF] 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. One last comment about Eq. ( 46) is that it does not depend on the expression of the disjoining energy (Eq. ( 30)) but only on its asymptotic behavior given by Eqs. ( 28) and ( 29).

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 [START_REF] Sharma | Relationship of thin film stability and morphology to macroscopic parameters of wetting in the apolar and polar systems[END_REF] 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:

dU i dt + 1 ∆x F i+1/2 -F i-1/2 = 1 ∆x B i+1/2 -B i-1/2 + S i (48)
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:

F i+1/2 = u + i+1/2 U i + u - i+1/2 U i+1 + 0, 0, L i+1/2 T ( 49
)
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:

L i+1/2 = 1 2 (L i + L i+1 ) (50)
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:

B i+1/2 = A i+1/2 ϕ i+1 -ϕ i ∆x ( 51
)
where A i+1/2 = 1 2 (A i + A i+1 ). 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:

U n+1 i -U n i + ∆t ∆x F n+1 i+1/2 -F n+1 i-1/2 = ∆t ∆x B n+1 i+1/2 -B n+1 i-1/2 + ∆tS n+1 i (52)

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:

∆x de i dt + ϕ i • F i+1/2 -F i-1/2 = ϕ i • B i+1/2 -B i-1/2 + ϕ i • S i ∆x (53)
Using the expression of F i+1/2 and the definition of ϕ i we get:

ϕ i • F i+1/2 -F i-1/2 = ϕ i • u + i+1/2 U i + u - i+1/2 U i+1 + 1 2 u i (L i + L i+1 ) (54)
Then, using that L i = ϕ i • U ie i (by definition of L), we get after some algebra:

ϕ i • F i+1/2 -F i-1/2 = (ue) i+1/2 + (uL) i+1/2 -(ue) i-1/2 -(uL) i-1/2 + D i (55)
where:

(ue

) i+1/2 = u + i+1/2 e i + u - i+1/2 e i+1 ( 56
) (uL) i+1/2 = 1 2 (u i+1 L i + u i L i+1 ) (57)
D i = -u - i+1/2 (e i+1 -e i -ϕ i • (U i+1 -U i )) + u + i-1/2 (e i-1 -e i -ϕ i • (U i-1 -U i )) (58)
The skew-symmetry of A i+1/2 leads to:

ϕ i • A i+1/2 ϕ i+1 -ϕ i ∆x = ϕ i + ϕ i+1 2 • A i+1/2 ϕ i+1 -ϕ i ∆x (59)
In the absence of external driving forces:

ϕ i • S i = - 3µu 2 i h i + b ≤ 0 (60)
Gathering all the intermediate results, we finally obtain the following local semi-discrete energy inequality (in the absence of external driving forces):

∆x de i dt + G i+1/2 -G i-1/2 + Φ i+1/2 -Φ i-1/2 + D i ≤ 0 ( 61
)
where:

G i+1/2 = u + i+1/2 e i + u - i+1/2 e i+1 + 1 2 (u i+1 L i + u i L i+1 ) (62) Φ i+1/2 = ϕ i + ϕ i+1 2 • A i+1/2 ϕ i+1 -ϕ i ∆x ( 63
)
Summing Eq. ( 61) over all cells, we finally get: 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.

dE dt = - i∈Z D i (

Numerical results and recommended best practices regarding the choice of h

h stands for the regularizing scale of the long-range molecular forces in the vicinity of the contact line. If we are only interested by the macroscopic behavior of the film, it can be expected that h does not need to be as small as R. That is what will be shown in the numerical examples presented in this section. In practice, two different cases must be considered depending on the targeted application.

On the one hand, if the mesh size needs to be chosen such that ∆x > h 0 , it is obvious that the numerical method cannot be expected to capture the correct contact angle whatever the value of the parameter h . In that case, it will be shown that the choice h 0 /100 h h 0 /10 allows to capture the correct spreading of the film. Larger values of h may lead to spurious dewetting phenomena (due to spinodal instabilities). Smaller values may lead to a wrong localization of the contact line forces since h needs to be greater than the residual film thickness downstream of the contact line due to the numerical diffusion of the scheme.

On the other hand, if the mesh size can be chosen such that ∆x h 0 , 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 0 . 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 [START_REF] De Gennes | Capillarity and Wetting Phenomena: Drops, Bubbles, Pearls, Waves[END_REF]. The thickness profile h i of the initial droplet is given by:

h i (x) = max   max R 2 i sin (θ i ) 2 -x 2 , 0 - R i tan (θ i ) , 0   ( 65
)
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 [START_REF] De Gennes | Capillarity and Wetting Phenomena: Drops, Bubbles, Pearls, Waves[END_REF]:

h 0 = 2l c sin θ s 2 (66)
l c h 0 V 1.749mm 1.749mm 110.1mm 2
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):

V = Ri -Ri h i (x) dx = θ i R 2 i sin 2 (θ i ) - R 2 i tan (θ i ) (67)
The spreading W 0 is then defined by:

W 0 := V h 0 (68)
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.

3.1.1. Reference numerical solution. Both h 0 and h are spatially well resolved. The computational domain, set to [-L 0 /2; L 0 /2] with L 0 = 0.08m, is discretized into 14720 cells leading to a space step ∆x = 5.43 µm. Different values for h are studied ranging from h /h 0 = 1.55 × 10 -3 to h /h 0 = 2.0 × 10 -1 , or equivalently, from h /∆x = 0.5 to h /∆x = 64.3. For each configuration, the steady state solution h(x) is computed and the thickness h(0) at the middle of the puddle is compared to the theoretical value h 0 (Eq. (66)) by computing:

ε h0 = |h(0) -h 0 | h 0 (69)
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. But, 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 65)) to the final puddle (Fig. 4). h /h 0 = 1/100 and ∆x/h = 0.25. 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 focus on the reference solution is proposed in the upper left corner of Fig. 7. The three regions (proximal, central and distal) are represented. As shown in Sec. 1.5, the geometric angle θ is equal to the static contact angle θ s in the central region (Young-Dupré equation, Eq. ( 47)).

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: Figure 6. Influence of h on ε h0 . The theoretical puddle thickness h 0 is given by Eq. (66). The X axis is ordered from the highest to the lowest values. ∆x = 5.43 µm.

ε h (α) = 1 L 0 L 0 2 - L 0 2 h # α (x) -h # 1/8 (x) 2 dx 1/2 (70a) ε p (α) = 1 L 0 L 0 2 - L 0 2 dh # α dx (x) -p # α (x) 2 dx 1/2 (70b)
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.

Figure 9 represents the time evolution of the film energy e (Eq. ( 18)) and de dt for h /h 0 = 1/100 and ∆x = h /4. In the absence of external driving forces, the energy e is indeed a decreasing function over time, which illustrates the result presented in Sec. 1.3 and Sec. 2.2.

h and h 0 are poorly spatially resolved

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 Figure 7. Focus on the contact line. Grey dash-dot line: reference slope at the contact line derived from the theoretical static contact angle θ s . Dashed lines: poor spatial discretization of the long-range forces (h /∆x 1). Solid lines: accurate spatial discretization of the long-range forces (h /∆x 1). Upper left corner: focus on the three regions (proximal, central and distal) defined in Sec. 1.5 for the reference solution h /h 0 = 1/100. Note that the boundaries between the three regions are given here approximately. The theoretical puddle thickness h 0 is given by Eq. (66). ∆x = 5.43 µm. θ 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. When h /h 0 1 (dashed 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 (46) may then be rewritten as follows:

tan (θ m ) • [cos (θ s ) -cos (θ m )] = -3C a ln 1 + H b (71)
Under the hypothesis that θ s 1 and θ m 1 and given that H b, the previous equation leads to:

θ m θ 2 s -θ 2 m = 6C a ln H b (72)
which has been introduced by de Gennes [START_REF] De Gennes | Dynamics of wetting: local contact angles[END_REF]. As a reminder, the similar Cox-Voinov equation [START_REF] Voinov | Hydrodynamics of wetting[END_REF][START_REF] Cox | The dynamics of the spreading of liquids on a solid surface. part 1. viscous flow[END_REF] is given by:

θ 3 s -θ 3 m = 9C a ln H b (73) ρ (kg.m -3 ) ν (m 2 .s -1 ) γ lg (N.m) θ s ( • ) 1000 10 -6 0.068 4
Table 3. Properties of the liquid droplet.

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)). 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 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. vs. C a is accurately reproduced for large values of C a . At early times, the deviations with the linear curves are more significant. This may be explained by the fact that Eqs. ( 72) and (73) are valid for θ s , θ m 1 and that the initialization angle θ i is set to 50 • , far from the validity range of Eqs. ( 72) and (73).

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 Figure 13. Time evolution of the contact line position x CL . Influence of h /h s . Highlighting of the "visco-capillary" and "long-range + capillary" regimes (see Fig. 14). The black dashed vertical lines (t = 4.0 × 10 -3 s and t = 0.2s) outline the boundaries between the two regimes. h /∆x > 1 (Tab. 4). 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 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. 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.

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

Figure 1 .

 1 Figure 1. Definition of the slip length b so that µ ∂u h ∂z x,z=0 = µ u h (x,0) b/3

Figure 3 .

 3 Figure 3. G function vs. h/h .

Figure 4 .

 4 Figure 4. Puddle in equilibrium. h 0 : puddle thickness. W 0 : puddle spreading. θ s : static contact angle.

Figure 5 .

 5 Figure 5. Time evolution from the initial droplet (Eq. (65)) to the final puddle (Fig.4). h /h 0 = 1/100 and ∆x/h = 0.25.

  (a) Error ε h (Eq. (70a)).(b) Error εp (Eq. (70b)).

Figure 8 .

 8 Figure 8. Mesh convergence study for the most accurate solution h /h 0 = 1/100. The X axis is ordered from the highest to the lowest values.

Figure 9 .

 9 Figure 9. Time evolution of e (Eq. (18)) and de dt . As shown in Sec. 1.3, the film energy e is a decreasing function over time. h /h 0 = 1/100 and ∆x = h /4.

( a )

 a Only half of the puddle is represented (x ≤ 0). Dashed lines: h /h0 1. Solid lines: h /h0 1. The curves h /h0 = 1/100 and h /h0 = 1/1000 are superimposed. (b) Focus near the contact line.

Figure 10 .

 10 Figure 10. Poor spatial resolution (∆x/h 0 ≈ 23). Influence of h on the spreading. For visualisation purposes, the X to Y ratio is set to 140.

Figure 11 .

 11 Figure 11. Contact line spreading in the absence of gravity. Time evolution from the initial droplet. h /h s = 1/33, h /∆x > 1.

3 m

 3 Figure 12. θ 3 sθ 3 m vs. C a . The linear regression lines are plotted (dashed lines). Comparison between h /h s = 1/5 and h /h s = 1/33. h /∆x > 1.

Figure 14 .

 14 Figure 14. Contact line spreading. Time evolution of the magnitudes of the different contributions in the force balance (Eq. (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.

Table 1 .

 1 -3 ) ν (m 2 .s -1 ) γ lg (N.m) θ s ( Properties of the liquid droplet.

	64)

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 ce d are convex functions ρ (kg.m

Table 4 .

 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.

	h (m) 1.302 × 10 -5 6.51 × 10 -6 3.255 × 10 -6 1/20 h /h 0 h /∆x 1/5 13.02 24.76 α 1/10 6.51 27.39 3.255 27.74 2.0 × 10 -6 1/33 2.0 27.1	ε α ×10 -5 18.6 12.5 9.0 7.8

Time evolution from the initial droplet. h /h s = 1/33, h /∆x > 1.

steady state is reached.

three-dimensional films. Encouraging results have already been obtained [17] and have been published in [33]. 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.


Keywords: Centrifuge Modelling, Inclined Pile Foundation, Seismic response, Soil-structure interaction. 1

Although inclined pile foundations are sometimes used in civil engineering structures, their behavior under seismic loadings is not yet thoroughly understood. This paper presents an experimental work on inclined pile foundations under earthquake loadings. Dynamic centrifuge tests are performed to investigate the influence of inclined piles on the superstructures and to quantify the soil-structure interaction. Some key parameters are analyzed such as the residual bending moment, the total bending moment and the axial forces. It is found that in certain cases inclination plays a beneficial role to the seismic behavior of the pile foundation system and that the presence of inclined piles increases its natural frequency. Comparing to a vertical pile foundation, the residual bending moment envelope curve is different and larger values are observed. This novel experimental work provides significant data to understand the behavior of inclined pile foundations in seismic regions.

INTRODUCTION

Inclined piles, also called batter piles or raked piles are widely introduced in civil engineering constructions where a substantial lateral resistance is required. However nowadays, many building codes or standards like AFPS 1990 [START_REF] Afps | Recommandations AFPS 90[END_REF] and Eurocode-8 [START_REF]Eurocode8-Part5, Eurocode 8: Design of structures for earthquake resistance-Part 5: Foundations, retaining structures and geotechnical aspects[END_REF] do not recommend using inclined piles in seismic regions. To be less conservative, codes like ACI318-05 and ACI318-11, although they do not prohibit the use of inclined piles in earthquake-resistant structures, pay attention to the potential damages at the junction of the inclined piles and the buildings caused by important seismic forces. The main drawbacks of the inclined piles are generally identified as: large forces onto the pile cap, reduction of the bending capacity due to the induced tensile forces, unfavorable rotation on the cap and additional bending moment due to the settlement of the soil [START_REF] Gerolymos | Evidence of beneficial role of inclined piles: observations and summary of numerical analyses[END_REF].

Different field observations, like the wharf failure in the port of Okland in 1989 (Loma Prieta earthquake), and in the port of Los Angeles in 1994 (Nothridge earthquake), reveal the unsatisfactory performance of inclined piles. However, recent studies put in evidence a possible positive performance of inclined piles. As reported by Gazetas and Mylonakis [START_REF] Gazetas | Seismic Soil-Structure Interaction: New Evidence and Emerging Issues[END_REF] inclined piles, if properly designed, can be beneficial rather than detrimental both for the structure they support and the piles themselves. Pender [START_REF] Pender | Aseismic pile foundation design analysis[END_REF] and Berrill et al. [START_REF] Berrill | Case study of lateral spreading forces on a piled foundation[END_REF] suggest also important beneficial effects regarding the use of inclined piles. However, the argument about whether the use of inclined piles is detrimental or beneficial is still unsettled. In 2004, Harn [START_REF] Harn | Have Batter Piles Gotten a Bad Rap in Seismic Zones (Or Everything You Wanted to Know About Batter Piles But Were Afraid to Ask)[END_REF] pointed out that the poor performance of inclined piles in the past earthquakes might be due to the lack of knowledge and of analytical tools. The use of displacement based design and of carefully detailed inclined piles may result in a better behavior and significant project savings.

In addition, Giannakou et al. [START_REF] Giannakou | Seismic Behavior of Inclined Piles[END_REF][START_REF] Giannakou | Seismic Behavior of Batter Piles : Elastic Response[END_REF] studied numerically, in the time domain, the performance of batter piles. In their study, linear constitutive laws were assumed for the soil and the inclined pile groups and five different inclinations were considered. They found that, for seismic loadings and purely kinematic conditions, batter piles tended to confirm their negative reputation. However, when the total response is considered (kinematic and inertial interaction) their influence could be beneficial. The work of Giannakou et al. [START_REF] Giannakou | Seismic Behavior of Inclined Piles[END_REF][START_REF] Giannakou | Seismic Behavior of Batter Piles : Elastic Response[END_REF] has shown that in order to understand the performance of inclined piles the influence of the superstructure and the inertial response of the upper structural system have to be taken into account.

In this paper, dynamic centrifuge tests are performed to investigate the influence of inclined piles on the response of some superstructures and on the soil-structure interaction. The effects of inclined pile are highlighted by comparing the behavior of a 1X2 symmetric configuration of 15 o inclined piles supporting a SDOF (single degree of freedom) superstructure with that of a 1X2 vertical pile group that support the same SDOF superstructure. For both pile groups, short and tall superstructures are considered in order to study the influence of the inertial loading coming from different types of superstructures on the performance of inclined piles. The short superstructure represents a shear force dominated loading case and the tall building an overturning moment dominated loading case. For both configurations floating piles were considered and a stiff pile cap imposes a rigid connection at the pile head. The behavior of the different soilpile-superstructure configurations under dynamic excitations is analyzed in terms of residual bending moment, total bending moment and axial force.

EXPERIMENTAL PROGRAM

Introduction

All the dynamic centrifuge tests have been performed at 40 g level. In the following all the data ware presented at the prototype scale otherwise mentioned.

Centrifuge modelling of inclined pile foundations under seismic actions

Zheng Li 1,2 , Sandra Escoffier 1 , Panagiotis Kotronis 2 1 LUNAM Université, IFSTTAR, GERS, SV, F-44341 Bouguenais, France 2 LUNAM Université, Ecole Centrale de Nantes, Université de Nantes, CNRS, Institut de Recherche en Génie Civil et Mécanique (GeM), 1 Rue de la Noë, F-44321 Nantes, France email: zheng.li@ifsttar.fr, sandra.escoffier@ifsttar.fr, Panagiotis.Kotronis@ec-nantes.fr

The various model configurations tested are shown in Figure 1. In order to investigate the behavior induced by the inclined piles, each time the response is compared with the response of the corresponding vertical pile group configuration.

Experimental set-up

All the test are performed in dry dense N34 Fontainebleau sand (Dr=80%). The properties of this sand reported in Table 1.

Table 1: Fontainebleau sand NE34: Material properties [START_REF] Irina | Caractérisationi dynamique de sables de référence en laboratoire -Application `a la réponse sismique de massifs sableux en centrifugeuse[END_REF] The piles are made of aluminum alloy and have a hollow section. They are instrumented with pairs of strain gauges that enable the determination of the bending moment profile and the axial stress at different depths. The connection between the piles and the pile cap is assumed rigid. The pile heads are thus rigidly inter-connected via a stiff cap that prohibits rotations between the pile heads and the cap [START_REF] Mokwa | Rotational Restraint of Pile Caps during Lateral Loading[END_REF]. Table 2 gives the characteristics of the piles and of the pile cap at the model and the prototype scales.

Table 2: Characteristics of piles and of pile cap

The possible soil-cap interaction and the lateral resistance of the pile cap are not considered hereafter. To limit the soil-cap interaction, piles are designed to stand at a certain distance above the ground surface, see Figure 3. Two types of SDOF superstructures with different heights of gravity center are considered i.e. a short and a tall building. Both buildings are designed such as they have:

 the same resonant frequency under fixed base conditions,  the same top mass,  the same total weight (foundation and superstructure system). The frequencies of the short and tall building under fixed-base conditions are tuned at 2.0 Hz (prototype scale).

Instrumentation layout of the tests is presented in Fig. 4. The movement of soil column is followed by two vertical arrays of accelerometers (From CH03 to CH09 and from CH14 to CH19). In order to highlight potential boundary effects, these accelerations are compared with that recorded on the flexible container (Equivalent Shear Beam box). In addition to measurement of bending moment and axial load in both pile, the kinematic response of the cap and superstructure are followed through accelerometers and laser sensors. The double integration of the accelerometers, combined with the use of pass band filter, enables the determination of the dynamics displacement. 3 and plotted in Figure 5. Table 3 also presents the series of seismic tests that has been performed. All the dynamic centrifuge tests that have been performed using real earthquake excitations are listed in Table 4.

Table 4: List of dynamic centrifuge tests 3 PERFORMANCE OF INCLINED PILES WITH SEISMIC SOIL-PILE-SUPERSTRUCTURE INTERACTION (SSPSI) UNDER EARTHQUAKE EXCITATIONS

Frequency response of the superstructure

For each soil-pile-superstructure configuration the first six small base shaking tests (Table 3) are used to evaluate the frequency response of the superstructures. In this framework, functions are calculated between the accelerations measured on the top of the buildings and those measured near the soil surface. The estimated frequency responses are listed in Table 5. As mentioned before, the frequencies of the buildings (short and tall) under fixed-base conditions are tuned at 2.0 Hz. When the buildings are supported by pile groups embedded in dense sand, the SSPSI induces a decrease of the response frequency of the superstructure.

Table 5: Frequency response of the superstructure (short and tall buildings) on inclined and vertical pile groups

The presence of inclined pile reduces the effect of the SSPSI on the frequency response of the superstructure. For the case of the short and slender superstructures installed on the inclined pile groups, the frequencies are 0.14 Hz and 0.08 Hz higher than when they are supported by the vertical pile groups. In addition, an increase of the height of the gravity center of the superstructure enhances the SSPSI on the frequency response of the superstructure for both pile group configurations (-0.24 Hz and -0.18 Hz for the inclined and the vertical pile groups respectively).

3.2

Stresses in piles

Residual bending moments

Due to the permanent deformation of the soil and the influence of the superstructure, embedded deep pile foundations deviate from their initial position during an earthquake event. As a result, residual bending moments appear along the pile. Figure 6 shows the accumulated residual bending moment profiles for the successive earthquake events. The peak values, such as the peak values normalized by the peak values obtained for the short building supported by the vertical pile group, are given in Figure 7.

The presence of inclined piles changes the shape of the residual bending profiles. Inclined piles induce a "C" shape curve while an "S" shape curve appears for vertical piles. For both configurations, the type of superstructure has an effect on the depth of the maximum residual bending moment. In the case of the inclined pile group, when a short building is supported the depth of the maximum residual bending moment increases from 4 to 6 pile diameter and almost constant for the strong earthquakes. The increase of the height of the gravity center increase the depth of the maximum residual bending moment: it appears at a depth of around 6D for the first two small earthquakes and during and after the strong earthquakes, the depth increases up to 10D. Concerning the vertical pile group, when a slender building is supported it is difficult to identify the depth of the maximum residual bending moment due to the shape of the profile). However, it can be located at about 5 pile diameters when a short building is supported. Concerning the performance of inclined piles, for all configurations (type of superstructure and input), their presence increase the maximum residual bending moment in the piles. However the impact level depends on the type of superstructure, the base shaking and the loading history. It is shown in Figure 7 that for the first two small earthquakes the influence of inclined piles is similar for the two types of superstructures. When strong earthquakes are applied however, the effect of inclined piles is more pronounced when a short super-structure is supported. For the last two small earthquakes, the effect of inclined piles is more pronounced for the slender superstructure. This is certainly due to the accumulation of residual bending moments from the previous earthquake events.

Total bending moments

Figure 8 shows the envelop curve of the maximum total (residual + dynamic) bending moments normalized by the maximum total bending moment obtained in the case of the short building supported by the vertical pile group. Due to the presence of the residual bending moments the curves are no symmetric and it is difficult to define a general shape. Qualitatively speaking, the response of the vertical pile group is more significant.

Figure 8: Normalized envelope curves of the maximum total resultant bending moments Figure 9 shows the ratios between the peak values of the total bending moment. The presence of inclined piles reduces the maximum total bending moment for both short and slender superstructures. However, the good performance of the inclined piles is influenced by the type of superstructure and the base shaking frequency content. A better performance is found for the short superstructure. Contrary to the slender superstructure, this performance is not influenced by the frequency of the base shaking.

Axial forces

Figure 10 shows the normalized total axial force profiles for the different configurations. For the friction piles, a reduction of the axial force with increasing pile embedded depth is observed, although in some cases the axial force profiles are distorted. Due to small number of instrumented compression strain gauges it is difficult to precisely determine the shape of the axial force profiles. According the results shown in Figure 11, for both, short and slender superstructures, inclined piles introduce a smaller or equal axial force to the vertical ones. For both configurations, the increase of the gravity center results in more significant axial forces. 1. Superstructures supported by inclined piles group have higher resonant frequencies. This is probably due to the higher horizontal stiffness provided by the inclined piles. In addition, regardless the pile group configuration, an increase of the gravity center of the superstructure results in a decrease of the resonant frequency. 2. For all the configurations tested, in terms of residual bending moments, inclined piles have a poor performance compared with the vertical piles. Furthermore the performances are influenced by the superstructure, the input and the loading history. 3. For the same type of superstructure (short or slender), less total bending moment and axial force appear in the piles for the inclined pile configuration. 4. For both inclined and vertical pile configurations, an increase of the gravity center induces higher axial forces in the piles. 5. In terms of total bending moments, the presence of inclined pile reduces the maximum total bending moment below and above the soil surface. In addition a higher gravity center of the superstructure may have a beneficial effect on the performance of both inclined and vertical configurations. However this effect is more pronounced for vertical pile group.

Figure 1 :Figure 2 :

 12 Figure 1: Model configurations used in the dynamic centrifuge tests 2.2 Base shaking signals Three real earthquakes are selected. The first two signals, Martinique Jara and the Northridge records, are broad band frequency range earthquakes (Figure 2(a) and (b)). They have been selected to enables the determination of the frequency response of the soil column and of the soil-pile-superstructure system. In order to determine the natural frequency of the superstructure system, these signals are first attenuated to small intensities to avoid significant plasticity in the soil. The third signal is a recording of the Kobe earthquake. It has been selected because most of the energy comes from frequencies below 1.5 Hz, see Figure 2(c), that are lower than the response frequency of the soil column and the soil-pile-superstructure systems. To be in the capacity range of the IFSTTAR shaker [10], the three earthquake signals have been filtered within the frequency range of 20∼350 Hz at the model scale. Furthermore, when large intensity earthquakes are concerned, they should be attenuated to be in the acceleration range of the shaker (-9db and -4db for respectively the Northridge and the Kobe records)

Figure 3 :

 3 Figure 3: Spacing of piles in inclined pile and vertical pile groups (prototype scale).

Figure 4 :

 4 Figure 4: Distribution of sensors in the soil-pilesuperstructure system for the centrifuge tests The adopted shaking sequence of the earthquake signals and the corresponding variation of Arias intensity are listed in Table3and plotted in Figure5. Table3also presents the series of seismic tests that has been performed.

Figure 5 :

 5 Figure 5: Variation of arias intensity

Figure 6 :

 6 Figure 6: Normalized residual bending moment profiles

Figure 7 :

 7 Figure 7: Normalized residual bending moments

Figure 9 :

 9 Figure 9: Normalized total bending moments at pile head

Figure 10 :

 10 Figure 10: Normalized axial forces

Figure 11 :

 11 Figure 11: Normalized maximum axial forces 4 CONCLUSIONS In this study, the performance of inclined piles and the soilpile-superstructure interactions (SSPSI) are analyzed and discussed. Cross comparisons are carried out considering the presence of inclined piles and the position of the gravity center of the superstructure. The main conclusions are the following:1. Superstructures supported by inclined piles group have higher resonant frequencies. This is probably due to the

Table 3 :

 3 Shaking sequence

Proceedings of the 9th International Conference on Structural Dynamics, EURODYN 2014

ACKNOWLEDGMENTS

The financial supports of IFSTTAR (Institut français des sciences et technologies des transports, de l'aménagement et des réseaux) and of the Région Pays de la Loire are acknowledged. The authors would like to thank the technical staff of the IFSTTAR centrifuge team for the valuable support and help.


Azeddine Baalal
email: abaalal@gmail.com

Mohamed Berghout

S Barza

DENSITY PROPERTIES FOR FRACTIONAL SOBOLEV SPACES WITH VARIABLE EXPONENTS

Keywords: Mathematics Subject Classification, Primary 46E35; Secondary 35A15, Keywords, fractional Sobolev spaces with variable exponents, density properties

In this article we show some density properties of smooth and compactly supported functions in fractional Sobolev spaces with variable exponents. The additional difficulty in this nonlocal setting is caused by the fact that the variable exponent Lebesgue spaces are not translation-invariant.

Introduction

Function spaces with variable exponents have been the subject of intensive investigation in recent years. Examples of such spaces include Lebesgue and Sobolev spaces with variable exponents. Introduced by Orlicz [START_REF] Orlicz | Über konjugierte Exponentenfolgen[END_REF] in 1931, their properties were further developed by Nakano [START_REF] Nakano | Modulared Semi-Ordered Linear Spaces[END_REF] as special cases of the theory of modular spaces. In the ensuing decades they were primarily considered as important examples of modular spaces or the class of Musielak-Orlicz spaces. Initially of theoretical interest, by the end of the twentieth century these function spaces moved beyond theory and into the area of variational problems and studies of the p(x)-Laplacian operator, which in turn further fueled the development of this theory. Also stimulating the wide-ranging investigation of these spaces were their application to various problems in applied mathematics, for example, in the areas of nonlinear elasticity theory, fluid mechanics, and mathematical modeling of physical phenomena. (For more details on these spaces, see the monographs [6], [START_REF] Kokilashvili | Integral Operators in Non-Standard Function Spaces[END_REF], [18]; see also [START_REF] Scapellato | Homogeneous Herz spaces with variable exponents and regularity results[END_REF].)

Let Ω be a bounded domain of R n . We fix s ∈ (0, 1), and we consider two variable exponents-that is, q : Ω → [1, +∞) and p : Ω × Ω → [1, +∞)-to be two continuous functions. In the following, we use the symbol := to make the left-hand side equal by definition to the right-hand side.

The variable exponent Lebesgue space L q(•) (Ω) is defined by

L q(•) (Ω) := u : Ω → R measurable, Ω |u(x)| q(x) dx < ∞ .
We define a norm, the so-called Luxembourg norm, in this space by u L q(•) (Ω) := inf λ > 0 :

Ω u(x) λ q(x)
dx ≤ 1 .

We define the fractional Sobolev space with variable exponents via the Gagliardo approach as follows:

W s,q(•),p(•,•) (Ω)

:= u ∈ L q(•) (Ω), Ω Ω
|u(x) -u(y)| p(x,y) λ p(x,y) |x -y| n+sp(x,y) dx dy < ∞, for some λ > 0 .

Let [u] s,p(•,•) (Ω) := inf λ > 0 :

Ω Ω |u(x) -u(y)| p(x,y) λ p(x,y) |x -y| n+sp(x,y) dx dy ≤ 1 be the corresponding variable exponent Gagliardo seminorm. It is easy to see that W s,q(•),p(•,•) (Ω) is a Banach space with the norm

u W s,q(•),p(•,•) (Ω) := u L q(•) (Ω) + [u] s,p(•,•) (Ω).
The definition of the spaces L q(•) (R n ) and W s,q(•),p(•,•) (R n ) is analogous to that of L q(•) (Ω) and W s,q(•),p(•,•) (Ω); one just changes every occurrence of Ω by R n (for more detail, see [START_REF] Baalal | Traces and fractional Sobolev extension domains with variable exponent[END_REF], [START_REF] Bahrouni | On a new fractional Sobolev space and applications to nonlocal variational problems with variable exponent[END_REF], [4], and [16]).

Fractional Sobolev spaces with variable exponents have major applications in variational problems related to a well-known fractional version of the p(x)-Laplace operator, given by div(|∇u| p(x)-2 ∇u), that is associated with the variable exponent Sobolev space W 1,p(•) (Ω). We refer, for instance, to [6] and [START_REF] Harjulehto | Overview of differential equations with non-standard growth[END_REF].

The denseness of regular functions such as C ∞ (Ω), functions which are continuously differentiable arbitrarily many times (smooth functions), or C ∞ 0 (Ω), the subset of C ∞ (Ω) of functions which have compact support, was one of the central tools in the function spaces setting-it might have some additional ramifications in the variable exponent Sobolev spaces. Moreover, the density of regular functions was one of the questions considered early (from 1986) in the context of minimizers of variational integrals. It is well known that the regularity of the Dirichlet energy integral minimizer is related to the density of smooth functions in the corresponding function space (see [14]

, [23]-[26]).
The denseness of smooth functions in the variable exponent Sobolev spaces W k,p(•) (Ω) has proved to be a more difficult problem. Nowadays, denseness is a topic of intense focus and study by many researchers (we refer the reader to We recall that the continuity of the shift operator T h : L q(•) → L q(•) , T h f (x) := f (x + h) may fail (see [5]). In connection with the density problem for smooth function, we introduce the most important condition on the exponent in the study of variable exponent spaces, the well-known log-Hölder continuity condition: a function q : Ω → R is called log-Hölder continuous on Ω if there exists C > 0 such that

q(x) -q(y) ≤ C -log |x -y| , |x -y| ≤ 1 2 . (1.1)
We say that a function p :

Ω × Ω → R satisfies condition (B-B) on Ω × Ω if there exists C > 0 such that p(x, y) -p(x , y ) ≤ C -log(|x -x | + |y -y |) , (1.2)
for all (x, y) and (x , y

) ∈ Ω × Ω such that |x -x | + |y -y | ≤ 1 2 .
Notice that conditions (1.1) and (1.2) ensure that the typical mollifier function approximate function in Lebesgue Spaces with variable exponents (we refer the reader to [5], [21], [START_REF] Bahrouni | On a new fractional Sobolev space and applications to nonlocal variational problems with variable exponent[END_REF] and [START_REF] Dipierro | A density property for fractional weighted Sobolev spaces[END_REF]).

We define the following class of variable exponents P log (Ω) := {q : Ω → R : q is measurable and log-Hölder continuous} and P log (Ω × Ω) := p : Ω × Ω → R : p is measurable and satisfies condition (B-B) .

We set p -:= essinf (x,y)∈Ω×Ω p(x, y), p + := esssup (x,y)∈Ω×Ω p(x, y), q -:= essinf x∈Ω q(x), and q + := esssup x∈Ω q(x).

We assume that

1 < p -≤ p(x, y) ≤ p + < ∞, (1.3) 1 < q -≤ q(x) ≤ q + < ∞, (1.4) p (x, y) -(z, z) = p(x, y), ∀(x, y), (z, z) ∈ Ω × Ω.
(1.5)

The aim of this paper is to prove a number of density properties of smooth and compactly supported functions in fractional Sobolev spaces with variable exponent. The first result is an approximation with a continuous and compactly supported function.

Theorem 1.1. Let u ∈ W s,q(•),p(•,•) (R n ). Then for any fixed δ > 0, there exists a continuous and compactly supported function u δ such that

u -u δ W s,q(•),p(•,•) (R n ) -→ 0 as δ -→ 0.
The second result is an approximation with smooth and compactly supported functions. More precisely, we have the following. Theorem 1.2. Let q ∈ P log (Ω) and p ∈ P log (Ω × Ω). Then the space

C ∞ 0 (R n ) is dense in W s,q(•),p(•,•) (R n ).
We consider Ω ⊂ R n to be a W s,q(•),p(•,•) -extension domain if there exists a continuous linear extension operator

E : W s,q(•),p(•,•) (Ω) -→ W s,q(•),p(•,•) (R n )
such that Eu| Ω = u for each u ∈ W s,q(•),p(•,•) (Ω). From Theorem 1.2, we derive the density of C ∞ (Ω) also in W s,q(•),p(•,•) -extension domains, which include in particular domains with Lipschitz boundary (see [START_REF] Baalal | Traces and fractional Sobolev extension domains with variable exponent[END_REF]).

Theorem 1.3. Let q ∈ P log (Ω), p ∈ P log (Ω × Ω) and suppose that Ω is a W s,q(•),p(•,•) -extension domain. Then the space C ∞ (Ω) is dense in W s,q(•),p(•,•) (Ω).
This paper is organized as follows. In Section 2, we describe some properties of fractional Sobolev spaces with variable exponents that will be useful in our exposition. In Section 3, we prove the density results for fractional Sobolev spaces with variable exponents.

Notation and preliminary results

In this section we give some notation and prove several useful properties of fractional Sobolev spaces with variable exponents.

We use the following notation throughout this article. R n is the n-dimensional Euclidean space, and n ∈ N always stands for the dimension of the space. A domain Ω ⊂ R n is a connected open set equipped with the n-dimensional Lebesgue measure. For constants, we use the letter C whose value may change even within a string of estimates. A ball with radius R and center 0 will be denoted by B R . The closure of a set A is denoted by A. We use the usual convention of identifying two µ-measurable functions on A (almost everywhere (a.e.) in A, for short) if they agree almost everywhere-that is, if they agree up to a set of µ-measure zero. The Lebesgue integral of a Lebesgue measurable function f : Ω -→ R is defined in the standard way and denoted by Ω f (x) dx. By ω n-1 , we denote the (n -1)-dimensional measure of the unit sphere S n-1 . By supp f , we denote the support of f , which is the complement of the biggest open set on which it vanishes; in other words, supp f is the closure of the set {x; f (x) = 0}. We denote by C(Ω) the space of uniformly continuous functions equipped with the supremum norm

f ∞ = sup x∈Ω |f (x)|. By C k (Ω), k ∈ N, we denote the space of all functions f such that ∂ α f := ∂ |α| f ∂ α 1 x 1 •••∂ αn xn ∈ C(Ω) for all multi-indices α = (α 1 , α 2 , . . . , α n ), |α| := α 1 + α 2 + • • • + α n ≤ k. The space is equipped with the norm sup |α|≤k ∂ α f ∞ , C ∞ (Ω) = k C k (Ω).
The set of smooth functions in Ω is denoted by C ∞ (Ω)-it consists of functions in Ω which are continuously differentiable arbitrarily many times. The set C ∞ 0 (Ω) is the subset of C ∞ (Ω) of functions which have compact support.

The important role of manipulating Lebesgue-Sobolev spaces with variable exponents is played by the modular of the L q(•) (Ω) space, which is the mapping

ρ : L q(•) (Ω) -→ R + defined by ρ(u) = Ω |u(x)| q(x) dx.
From [10, Theorems 1.4, 1.8], we obtain the following proposition. Proposition 2.1. Let Ω be an open subset of R n , and let u n , u ∈ L q(•) (Ω) be such that n ∈ N. Then

(1) C ∞ 0 (Ω) is dense in the space (L q(•) (Ω), • L q(•) (Ω) ), (2) lim n-→+∞ ρ(u n -u) -→ 0 ⇐⇒ lim n-→+∞ u n -u L q(•) (Ω) = 0.
Next, we prove several properties of W s,q(•),p(•,•) (R n ) that are needed to obtain our main results.

Lemma 2.2. Let ϕ ∈ C ∞ 0 (R n ).
Then there exists a positive constant C such that

R n R n |ϕ(x) -ϕ(y)| p(x,y) |x -y| n+sp(x,y) dx dy ≤ C. Proof. Let ϕ ∈ C ∞ 0 (R n ) with supp ϕ ⊆ B R . Therefore, if x, y ∈ R n \ B R , then ϕ(x) = ϕ(y) = 0. Thus I = R n R n |ϕ(x) -ϕ(y)| p(x,y)
|x -y| n+sp(x,y) dx dy

= B R R n |ϕ(x) -ϕ(y)| p(x,y) |x -y| n+sp(x,y) dx dy = I 1 + I 2 ,
where

I 1 := B R B 2R |ϕ(x) -ϕ(y)| p(x,y)
|x -y| n+sp(x,y) dx dy and

I 2 := B R R n \B 2R |ϕ(x) -ϕ(y)| p(x,y)
|x -y| n+sp(x,y) dx dy.

We first estimate I 1 . We have

I 1 = B R B 2R |ϕ(x) -ϕ(y)| p(x,y) |x -y| n+sp(x,y) dx dy = B R B 2R |ϕ(x) -ϕ(y)| p(x,y) |x -y| p(x,y) × |x -y| p(x,y) |x -y| n+sp(x,y) dx dy.
Hence

I 1 ≤ C B 2R B 2R |x -y| -n+(1-s)p(x,y) dx dy,
where C depends on the C 1 -norm of ϕ, p + , and p -. We set

d := sup |x -y| : (x, y) ∈ B 2R × B 2R .
Observe that |x -y| -n+(1-s)p(x,y) ≤ max{ -n+(1-s)p + , d -n+(1-s)p -}.

Therefore,

I 1 ≤ C B 2R B 2R dx dy ≤ C.
Now we estimate I 2 . For this, we observe that if

x ∈ B R and y ∈ R n \ B 2R , then |x -y| ≥ |y| -|x| ≥ |y| 2 .
Thus

I 2 ≤ C n, s, p + , p -, ϕ L ∞ (R n ) B R R n \B 2R dy |y| n+sp(x,y) dx. Now, if |y| > 1, then R n \B 2R dy |y| n+sp(x,y) ≤ R n \B 2R dy |y| n+sp -= ω n-1 C(R, s, p -) ≤ C,
where ω n-1 denotes the (n -1)-dimensional measure of the unit sphere S n-1 . Now, we assume that |y| < 1. We claim that

R n \B 2R dy |y| n+sp(x,y) ≤ R n \B 2R dy |y| n+sp + = ω n-1 C(R, s, p + ) ≤ C.
Hence, I is bounded.

As an obvious consequence of Lemma 2.2 and Proposition 2.1, we have

C ∞ 0 (R n ) ⊆ W s,q(•),p(•,•) (R n ). Now we give two approximation results. Lemma 2.3. Let u ∈ L q(•) (R n ).
Then there exists a sequence of functions

u m ∈ L q(•) (R n ) ∩ L ∞ (R n ) such that u -u m L q(•) (R n ) -→ 0 as m -→ +∞.
Proof. We set

u m (x) :=      m if u(x) ≥ m, u(x) if u(x) ∈ (-m, m), -m if u(x) ≤ -m. We have u m -→ u a.e. in R n and |u m (x)| q(x) ≤ |u(x)| q(x) ∈ L 1 (R n ).
By the dominated convergence theorem and Proposition 2.1, the claim follows.

Lemma 2.4. Let u ∈ L p(•,•) (R n × R n ).
Then there exists a sequence of functions

u M ∈ L p(•,•) (R n × R n ) ∩ L ∞ (R n × R n ) such that u -u M L p(•,•) (R n ×R n ) -→ 0 as M -→ +∞.
Proof. We set

u M (x, y) :=      M if u(x, y) ≥ M, u(x, y) if u(x, y) ∈ (-M, M ), -M if u(x, y) ≤ -M.
We have

u M -→ u a.e. in R n × R n and |u M (x, y)| p(x,y) ≤ |u(x, y)| p(x,y) ∈ L 1 (R n × R n ).
The claim follows from the dominated convergence theorem and Proposition 2.1.

The proofs of density properties of smooth and compactly supported functions are mainly based on a basic technique of convolution (which makes functions C ∞ ), joined with a cutoff (which makes their support compact). In the remainder of this article, we describe properties of these operations with respect to the norm in fractional Sobolev spaces with variable exponents.

Let η ∈ C ∞ 0 (R n ) be such that η ≥ 0 in R n and supp η ⊆ B 1 . We also assume that

B 1 η(x) dx = 1.
Let ε > 0, and let η ε be the mollifier defined as

η ε (x) = 1 ε n η x ε , x ∈ R n .
For any u ∈ W s,q(•),p(•,•) (R n ), we will denote by u ε the function defined as the convolution between u and η ε ; that is,

u ε (x) = (u * η ε )(x) = R n u(x -z)η ε (z) dz, x ∈ R n .
Note that, by construction, u ε ∈ C ∞ (R n ). For small ε, the convolution does not change the norm too much, according to the following result (see, e.g., [2, Lemma 3.2]), which will be useful in the rest of this article.

Lemma 2.5. Let u ∈ W s,q(•),p(•,•) (R n ), q ∈ P log (Ω) and p ∈ P log (Ω × Ω). We assume that (1.3), (1.4) and (1.5) hold. Then u -

u ε W s,q(•),p(•,•) (R n ) → 0 as ε → 0.
Proof. Let u ∈ W s,q(•),p(•,•) (R n ). Since u ∈ L q(.) (Ω) and q ∈ P log (Ω), by [21, Corollary], we know that

u -u ε L q(•) (R n ) → 0 as ε → 0.
Hence, from Proposition 2.1, it suffices to prove that Ω×Ω u ε (x) -u(x) -u ε (y) + u(y) p(x,y) K(x, y) dx dy → 0 as ε → 0, where K(x, y) = |x -y| -n-sp(x,y) . By the same way in [2, Lemma 3.2], we can prove that Ω×Ω u ε (x) -u(x) -u ε (y) + u(y) p(x,y) K(x, y) dx dy

≤ |B 1 | p -+p + -1 Ω×Ω×B 1 u(x -εz) -u(y -εz) -u(x) + u(y) p(x,y) × K(x, y) η(z) p + + η(z) p -dx dy dz.
Fix z ∈ B 1 and put w = (z, z) ∈ Ω × Ω. We define the function v : Ω × Ω -→ R by v(x, y) := u(x) -u(y) K(x, y)

1 p(x,y) , ∀(x, y) ∈ Ω × Ω. Then v ∈ L p(•,•) (Ω × Ω). Let ε > 0. Since p ∈ P log (Ω × Ω), by [21, Corollary], there exists g ∈ C ∞ 0 (Ω × Ω) with v -g L p(•,•) (Ω×Ω) ε 3 , hence v(• -εw) -v L p(•,•) (Ω×Ω) v(• -εw) -g(• -εw) L p(•,•) (Ω×Ω) + g(• -εw) -g L p(•,•) (Ω×Ω) + v -g L p(•,•) (Ω×Ω) ε 3 + ε 3 + ε 3 = ε ,
with ε is sufficiently small. Hence lim ε→0 Ω×Ω u(x -εz) -u(y -εz) -u(x) + u(y) p(x,y) × K(x, y) dx dy = 0.

Moreover, for a.e. z ∈ B 1 , there exists a positive constant C such that

η(z) p + + η(z) p - Ω×Ω u(x -εz) -u(y -εz) -u(x) + u(y) p(x,y) × K(x, y) dx dy 2C η(z) p + + η(z) p - Ω×Ω u(x) -u(y) p(x,y) × K(x, y) dx dy ∈ L ∞ (B 1 ),
for any ε > 0. Hence, using the dominated convergence theorem, we deduce that

lim ε→0 B 1 η(z) p (x,y) Ω×Ω u(x -εz) -u(y -εz) -u(x) + u(y) p(x,y)
× K(x, y) dx dy dz = 0.

Consequently Ω×Ω u ε (x) -u(x) -u ε (y) + u(y) p(x,y) K(x, y) dx dy → 0 as ε → 0, which concludes the proof. Now, we will discuss the cutoff technique needed for the density argument. For

any j ∈ N, let τ j ∈ C ∞ (R n ) be such that 0 ≤ τ j (x) ≤ 1, ∀x ∈ R n , τ j (x) = 1 if x ∈ B j , 0 if x ∈ R n \ B j+1 , (2.1)
where B j denotes the ball centered at zero with radius j. We have the following result.

Lemma 2.6.

Let u ∈ W s,q(•),p(•,•) (R n ). Then τ j u ∈ W s,q(•),p(•,•) (R n ). Proof. Let u ∈ W s,q(•),p(•,•) (R n ). It is clear that τ j u ∈ L q(•) (R n ) since |τ j | ≤ 1.
Furthermore, y) |x -y| n+sp(x,y) dx dy

I = R n R n |τ j (x)u(x) -τ j (y)u(y)| p(x,y) λ p(x,y) |x -y| n+sp(x,y) dx dy ≤ 2 p + -1 R n R n |τ j (x)(u(x) -u(y))| p(x,y) λ p(x,
+ 2 p + -1 R n R n |u(y)(τ j (x) -τ j (y))| p(x,y) λ p(x,y) |x -y| n+sp(x,y) dx dy ≤ 2 p + -1 R n R n |u(x) -u(y)| p(x,y) λ p(x,y) |x -y| n+sp(x,y) dx dy + 2 p + -1 R n R n |u(y)(τ j (x) -τ j (y))| p(x,y)
λ p(x,y) |x -y| n+sp(x,y) dx dy,

where the integral

R n R n |u(x) -u(y)| p(x,y) λ p(x,y) |x -y| n+sp(x,y) dx dy is finite since u ∈ W s,q(•),p(•,•) (R n ).
According to Lemma 2.3, we can assume that u ∈ L ∞ (R n ). Hence,

|u(y)| p(x,y) ≤ u p + L ∞ (R n ) + u p - L ∞ (R n ) . Therefore, R n R n |u(y)(τ j (x) -τ j (y))| p(x,y) λ p(x,y) |x -y| n+sp(x,y) dx dy ≤ C α R n R n |(τ j (x) -τ j (y))| p(x,y) |x -y| n+sp(x,y) dx dy,
where α = min{λ p + , λ p -} and the constant C depends on p + , p -, and u L ∞ (R n ) . Finally, using Lemma 2.2, we get

R n R n |u(y)(τ j (x) -τ j (y))| p(x,y) λ p(x,y) |x -y| n+sp(x,y) dx dy < ∞.
This concludes the proof. Now, we will prove that we can approximate a function in W s,q(•),p(•,•) (R n ) with another function with compact support, by keeping the error small. More precisely, we have the following result. Lemma 2.7. Let u ∈ W s,q(•),p(•,•) (R n ). Then supp(τ j u) ⊆ B j+1 ∩ supp u, and

τ j u -u W s,q(•),p(•,•) (R n ) -→ 0 as j -→ +∞.
Proof. By (2.1) and [11, Lemma 9], we get supp(τ j u) ⊆ B j+1 ∩ supp u.

Now, let us show that

τ j u -u W s,q(•),p(•,•) (R n ) -→ 0 as j -→ +∞.
From Proposition 2.1, it suffices to prove that

R n τ j (x)u(x) -u(x) q(x) dx -→ 0 as j -→ +∞ and R n R n |τ j (x)u(x) -u(x) -τ j (y)u(y) + u(y)| p(x,y) |x -y| n+sp(x,y)
dx dy -→ 0 as j -→ +∞.

We observe that

τ j (x)u(x) -u(x) q(x) = τ j (x) -1 q(x) u(x) q(x) ≤ 2 q + -1 τ j (x) q(x) + 1 u(x) q(x) ≤ 2 q + u(x) q(x) ∈ L 1 (R n ).
Moreover, by (2.1)

τ j (x)u(x) -u(x) q(x) -→ 0 as j -→ +∞ a.e. in R n .
Then, by the dominated convergence theorem, we get

R n τ j (x)u(x) -u(x) q(x) -→ 0 as j -→ +∞. Now, let us show that R n R n |τ j (x)u(x) -u(x) -τ j (y)u(y) + u(y)| p(x,y)
|x -y| n+sp(x,y) dx dy -→ 0 as j -→ +∞.

We set η j = 1 -τ j . Then η j u = u -τ j u. Moreover,

τ j (x)u(x) -u(x) -τ j (y)u(y) + u(y) = η j (x) u(x) -u(y) -τ j (y) -τ j (x) u(y) .
Therefore,

R n R n |τ j (x)u(x) -u(x) -τ j (y)u(y) + u(y)| p(x,y) |x -y| n+sp(x,y) dx dy ≤ 2 p + -1 R n R n |τ j (x) -τ j (y)| p(x,y)
|x -y| n+sp(x,y) u(y) p(x,y) dx dy

+ 2 p + -1 R n R n |u(x) -u(y)| p(x,y) |x -y| n+sp(x,y) η p(x,y) j (x) dx dy.
Note that by Lemma 2.3, we can assume that u ∈ L ∞ (R n ). Hence,

|τ j (x) -τ j (y)| p(x,y) |x -y| n+sp(x,y) u(y) p(x,y) ≤ C u L ∞ (R n ) , p + , p -|τ j (x) -τ j (y)| p(x,y) |x -y| n+sp(x,y) .
By Lemma 2.2, we deduce that

|τ j (x) -τ j (y)| p(x,y) |x -y| n+sp(x,y) ∈ L 1 (R n × R n ).
Furthermore,

|τ j (x) -τ j (y)| p(x,y) |x -y| n+sp(x,y) u(y) p(x,y) -→ 0 as j -→ ∞ a.e. in R n × R n .
Hence, by the dominated convergence theorem,

R n R n |τ j (x) -τ j (y)| p(x,y)
|x -y| n+sp(x,y) u(y) p(x,y) dx dy -→ 0 as j -→ ∞.

Also, |u(x) -u(y)| p(x,y) |x -y| n+sp(x,y) η j (x) p(x,y) ≤ η j (x)

p + + η j (x)) p -) |u(x) -u(y)| p(x,y) |x -y| n+sp(x,y)
and

|u(x) -u(y)| p(x,y) |x -y| n+sp(x,y) ∈ L 1 (R n × R n ).
Again by (2.1), |u(x) -u(y)| p(x,y) |x -y| n+sp(x,y) η j (x) p(x,y) -→ 0 as j -→ ∞ a.e. in R n × R n .

Hence, by the dominated convergence theorem,

R n R n |u(x) -u(y)| p(x,y)
|x -y| n+sp(x,y) η j (x) p(x,y) dx dy -→ 0 as j -→ ∞.

This concludes the proof.

Approximation by smooth and compactly supported functions

In this section we give some density properties of smooth and compactly supported functions in fractional Sobolev spaces with variable exponents. The additional difficulty in this nonlocal setting is caused by the fact that the variable exponent Lebesgue spaces are not translation-invariant (see [6, Section 3.6]). To overcome this difficulty, we exploit the techniques developed in Section 2.

Theorem 3.1. Let u ∈ W s,q(•),p(•,•) (R n ). Then for any fixed δ > 0, there exists a continuous and compactly supported function u δ such that

u -u δ W s,q(•),p(•,•) (R n ) -→ 0 as δ -→ 0.
Proof. Let u ∈ W s,q(•),p(•,•) (R n ). According to Lemma 2.2, we can approximate u with a sequence of bounded functions. Consequently, we can also assume that u ∈ L ∞ (R n ).

Let τ j ∈ C ∞ (R n ) be as in Section 2, with τ j (P ) = 1 if |P | ≤ j and τ j (P ) = 0 if |P | ≥ j + 1. Let u j = τ j u. Then u j -→ u a.e. in R n as j -→ ∞,

and u(x) -u j (x) q(x) ≤ 2 q + u(x) q(x) ∈ L 1 (R n ).
As a consequence, by the dominated convergence theorem, R n u(x) -u j (x) q(x) dx -→ 0 as j -→ ∞.

Hence, for any fixed δ > 0, there exists j δ ∈ N such that

R n u(x) -u j δ (x) q(x) dx ≤ δ. (3.1)
Note that u j δ is supported in B j+1 and that µ(A) = A dx is finite over compact sets. Hence, we can use Lusin's theorem (see [START_REF] Folland | Real Analysis: Modern Techniques and Their Applications[END_REF]Theorem 7.10,p. 121] for the definition of the uniform norm). We obtain that there exist a closed set E δ ⊂ R n and a continuous and compactly supported function u δ : R n -→ R such that

u δ = u j δ in R n \ E δ , µ(E δ ) ≤ δ and u j δ L ∞ (R n ) ≤ u L ∞ (R n ) .
In particular, since 0 ≤ τ j δ (x) ≤ 1, we have

u δ L ∞ (R n ) ≤ u L ∞ (R n ) < ∞. Therefore, R n u j δ (x) -u δ (x) q(x) dx = E δ u j δ (x) -u δ (x) q(x) dx ≤ 2 q + -1 E δ | u j δ (x) | q(x) dx + E δ u δ (x) q(x) dx) ≤ C(q + , q -, u L ∞ (R n )µ(E δ ). Hence R n u j δ (x) -u δ (x) q(x) dx ≤ Cδ. (3.2)
On the other hand,

u(x) -u δ (x) q(x) ≤ 2 q + -1 u(x) -u j δ (x) q(x) + u j δ (x) -u δ (x) q(x) . Hence R n u(x) -u δ (x) q(x) dx ≤ 2 q + -1 R n u(x) -u j δ (x) q(x) dx + 2 q + -1 R n u j δ (x) -u δ (x) q(x) dx.
Now, from (3.1) and (3.2), we deduce that

R n u(x) -u δ (x) q(x) dx -→ 0 as δ -→ 0.
Therefore, by Proposition 2.1, we obtain

u -u δ L q(•) (R n ) -→ 0 as δ -→ 0. Moreover, R n R n |(u -u j )(x) -(u -u j )(y)| p(x,y) |x -y| n+sp(x,y) dx dy -→ 0 as j -→ ∞.
Hence, for any fixed δ > 0, there exists j δ ∈ N such that

R n R n |(u -u j δ )(x) -(u -u j δ )(y)| p(x,y)
|x -y| n+sp(x,y) dx dy ≤ δ.

Notice that

v j δ (x, y) = |u j δ (x) -u j δ (y)| p(x,y) |x -y| n+sp(x,y) is supported in {P ∈ R n × R n ; |P | ≤ j δ + 1} and µ(A) = A A
dx dy is finite over compact sets. Hence, we can again use Lusin's theorem. We get that there exist a closed set E δ ⊂ R n × R n and a continuous and compactly supported function

u δ : R n × R n -→ R such that u δ = v j δ in R n × R n \ E δ , µ(E δ ) ≤ δ and u j δ L ∞ (R n ×R n ) ≤ v j δ L ∞ (R n ×R n ) .
In particular, since 0 ≤ τ jδ (x) ≤ 1, we have |x -y| n+sp(x,y) dx dy -→ 0 as δ -→ 0, which concludes the proof.

u δ L ∞ (R n ×R n ) ≤ |u(x) -u(y)| p(x,y) |x -y| n+sp(x,y) L ∞ (R n ×R n ) . Now,
Theorem 3.2. Let q ∈ P log (Ω) and p ∈ P log (Ω × Ω). Then the space C ∞ 0 (R n ) is dense in W s,q(•),p(•,•) (R n ).

Proof. We show that for any u ∈ W s,q(•),p(•,•) (R n ), there exists a sequence ρ ε ∈ C ∞ 0 (R n ) such that ρ ε -u W s,q(•),p(•,•) (R n ) -→ 0 as ε -→ 0.

Let u ∈ W s,q(•),p(•,•) (R n ), and let us fix δ > 0. Let τ j be as in Section 2. By Lemma 2.7, we know that

u -τ j u W s,q(•),p(•,•) (R n ) ≤ δ 2 ,
for j large enough.

For any ε > 0, let us consider

ρ ε := τ j u * η ε ,
where η ε is the mollifier function defined in Section 2. By construction, ρ ε ∈ C ∞ (R n ). Furthermore, by [3, Proposition IV.18],

supp ρ ε ⊆ supp(τ j u) + B ε .

Also, by Lemma 2.7, supp(τ j u) ⊆ B j+1 ∩ supp u. Hence, supp ρ ε ⊆ (B j+1 ∩ supp u) + B ε . As a consequence of this, ρ ε ∈ C ∞ 0 (R n ), for ε small enough. Furthermore, Lemma 2.5 gives ρ ε -τ j u W s,q(•),p(•,•) (R n ) ≤ δ 2 for ε small enough. Therefore,

u -ρ ε W s,q(•),p(•,•) (R n ) ≤ u -τ j u W s,q(•),p(•,•) (R n ) + τ j u -ρ ε W s,q(•),p(•,•) (R n ) ≤ δ 2 + δ 2 ≤ δ.
Since δ can be taken arbitrarily small, this concludes the proof. Theorem 3.3. Let q ∈ P log (Ω), p ∈ P log (Ω × Ω) and suppose that Ω is a W s,q(•),p(•,•) -extension domain. Then C ∞ (Ω) is dense in W s,q(•),p(•,•) (Ω).

Proof. Let u ∈ W s,q(•),p(•,•) (Ω). Since Ω is a W s,q(•),p(•,•) -extension domain, we find that ũ ∈ W s,q(•),p(•,•) (R n ) with ũ(x) = u(x) for all x ∈ Ω and ũ W s,q(•),p(•,•) (R n ) ≤ C u W s,q(•),p(•,•) (Ω) .

Due to Theorem 3.2, we can choose ũε ∈ C ∞ 0 (R n ) with ũε -→ u in W s,q(•),p(•,•) (R n ).

We set u ε := ũε | Ω . Then u -u ε W s,q(•),p(•,•) (Ω) ≤ ũ -ũε W s,q(•),p(•,•) (R n ) -→ 0.

Hence u ε ∈ C ∞ (Ω) are the required approximating functions.

  [6], [8]-[10], [15], [21], [27] and the references therein).

  by putting v(x, y) := |u(x) -u(y)| p(x,y) |x -y| n p(x,y) +sin Lemma 2.4, we can assume that|u(x) -u(y)| p(x,y) |x -y| n+sp(x,y) ∈ L ∞ (R n × R n ). Therefore R n R n |(u -u δ )(x) -(u -u δ )(y)| p(x,y)

Department of Mathematics and Computer Science, Hassan II University of Casablanca, Faculty of Sciences Aïn Chock, B.P. 5366 Maarif, Casablanca, Morocco.E-mail address: moh.berghout@gmail.com

Acknowledgment. The authors thank Lars Diening for insightful comments and for giving some very useful suggestions for the improvement of our paper.