Keywords: Navier-stokes surface libre, condition de glissement de Navier généralisée, dynamique de ligne triple, impact de gouttelette, Comsol multiphysics 3.2, Micralef

Navier-Stokes ALE surface libre avec conditions de glissement de Navier generalisée. Impact de gouttelettes et tentative de simulations avec Comsol Multiphysics 3.2. Résumé : Nous nous intéressons à la modélisation numérique d'une gouttelette impactant un solide. Les équations sont celles de Navier-Stokes surface libre en formulation ALE. La dynamique de la ligne triple est implicitement modélisée par les équations de Shikhmurzaev. Ces équations introduisent une condition de glissement de Navier généralisée, avec un gradient de tension de surface (effet Marangoni localisé). Les potentialités du logiciel Comsol Multiphysics 3.2 pour résoudre ce problème sont étudiées. Finalement, nous obtenons des résultats basés sur les algorithmes établis mais implémentés au sein de Micralef, un code de calcul spécialement développé pour cette étude.

Introduction

The aim of this work is to elaborate a mathematical and numerical model dealing with the dynamic contact line problem. Involved in the presence of three phases, where at least two of them are fluids (liquid-gas or non miscible liquid-liquid flow), the dynamics of this line significantly influences the interface as well as the bulk motion. The physical phenomenon of dynamical contact line appears in several industrial contexts such as coating and cleaning of solids by liquids. Two features of such slow, but not creeping, flows are revealed by experiments: i) the liquid front advances following a rolling motion: the particles of the liquid-gas interface are progressively projected onto the solide-liquid interface, in the same manner as does a rolling non-sliding sphere on a solide surface; ii) the contact angle depends on both its static value, determined thanks to the classical Young's equation, and the fluid velocity in the bulk. Even if succesfull attempts have been done to prescribe an explicit law of that contact angle, none of them is general because of the need to adjust their coefficients empirically. The contact angle and the triple line velocity should be treated by a more deep model, which describes and formalizes correctly the flow that occurs in the immediate vicinity of the triple line. The mathematical modeling of the moving contact line is hard to deal with. A classical no-slip boundary condition at the solid-liquid interface implies a non-physical singularity: the fluid exerts an infinite force on the solid surface (glass filling paradox). Most of models have then been based on a fully slip description of the solid-liquid interface. Even slip condition removes the singularity, it prevents the fluid from rolling, which is actually not the case as mentioned above. Furthermore, for a normal liquid flowing over an even smooth solid, slippage is usually negligible. So, the solution is expected to be between slipping and no-slipping conditions. In fact, according to Shikhmurzaev's theory [START_REF] Shikhmurzaev | The moving contact line problem on a smooth solid surface[END_REF], [START_REF] Blake | Experimental evidence of nonlocal hydrodynamic influence on the dynamic contact angle[END_REF], , for instance, generalized Navier condition has to be applied in the vicinity the triple line. Shikhmurzaev's theory is however a much deeper translation of the phenomena that happen near the triple line than a rude technique to overcome mathematical inconsistencies or to adjust an arbitrary approximation law. This paper is organized as follows. In section 2, we present the full Navier-Stokes equations (the 2D axisymmetric), the Shikhmurzaev surface equations (written in a planar surface co-ordinate). A link with the surfactant equations is done. In section 3, we recall the basic principles of the ALE method. In section 4, we present how the ALE method is implemented into Comsol (TM) Multiphysics 3.2. Since few important quantities (such as the curvature) are not easily computable in this software, we present some tricks to overcome these difficulties. Since the ALE implementation in Comsol Multiphysics 3.2 is not robust enough to handle large mesh distorsion, we present some preliminary numerical results obtained using our home-developed software Micralef. We compute the spreading phase of a 2D axisymmetric droplet impacting a solid surface. The algorithms implemented (in particular, the spreading algorithm) are those presented previously.

INRIA 2 Mathematical model 2.1 Navier-Stokes equations

The droplet dynamics is modeled by the unsteady incompressible 2D axisymmetric Navier-Stokes equation. We denote by u = (u r , u z ) T the fluid velocity, p its pressure, Σ the stress tensor, D the deformation tensor and Re the Reynolds number. We denote by ( τ , n) the unit tangential and external normal vectors, respectively, such that it is direct. We set Σ n = n • Σ; Σ n = Σ n • n; Σ τ = Σ n • τ and f the body force density. This equation writes down like

ρ D u Dt = ∇ • Σ + ρ f (1)
For a Newtonian fluid, the stress tensor is given by

Σ = -pI + µ ∇ u + (∇ u) T (2)
which gives

ρ ∂ u ∂t + u • ∇ u = -∇p + µ∇ 2 u + ρ f in [0, T ] × Ω f (3)
This equation is completed by the mass conservation equation

Dρ Dt + ρ∇ • u = 0 in [0, T ] × Ω f
which reduces for an incompressible fluid (as Dρ/Dt = 0) to

∇ • u = 0 in [0, T ] × Ω f (4)
We present the modelling of the free surface and its dynamics in the ALE section below. Equations ( 3) and ( 4) are completed by suitable initial and boundary conditions.

Surface equations of the Shikhmurzaev theory

According to the Shikhmurzaev's theory, [START_REF] Shikhmurzaev | The moving contact line problem on a smooth solid surface[END_REF], [START_REF] Blake | Experimental evidence of nonlocal hydrodynamic influence on the dynamic contact angle[END_REF] (see also [6] [8] for a reformulation of the equations) the surface equation to be solved writes down like

∂ρ s i ∂t -ξ i ∇ s • (ρ s i ∇ s ρ s i ) + ∇ s • (ρ s i u i ) + 1 τ * ρ s i = 1 τ * ρ eq i in [0, T ] × Γ m f ( 5
)
where ξ i = ξ for i = 1 and ξ i = λ for i = 2. ρ s i is a surface density at the interface i that is related to the local surface tension by a given law. Each of the surface differential operators has to be written in the surface co-ordinate system (s, φ), where s is an arbitrary azimutal parametrization of the (free) surface and φ is the meridional angle as depicted in figure 1.

Figure 1. Surface co-ordinate system, (s, φ).

Injecting the expressions of the different surface differential operators given by appendix 1, the Shikhmurzaev equation (5) becomes

∂ρ s i ∂t - 1 σ ξ i ∂ ∂s σρ s i ∂ρ s i ∂s + 1 σ ∂ ∂s (σρ s i u is ) + 1 τ * ρ s i = 1 τ * eq ρ eq i , (t, s) ∈ [0, T ] × [0, 1] (6)
where σ = r ∂r ∂s

2 + ∂z ∂s

2 1 2 (7)
Multiplying equation ( 6) by the test function ρs i associated to ρ s i and integrating over the relevant boundary ∂Ω, see figure 2, we get

∂Ω r ρs i ∂ρ s i ∂t ds = INRIA ξ i ρs i (s 1 ) ρ s i (s 1 ) ∂ ρs i ∂s (s 1 ) r (s 1 ) -ξ i ρs i (s 0 ) ρ s i (s 0 ) ∂ ρs i ∂s (s 0 ) r (s 0 ) - ∂Ω ξ i    ∂ ρs i ∂s 1 ∂r ∂s
2 + ∂z ∂s 2

ρ s i ∂ρ s i ∂s ds -[ρ s i (s 1 ) ρ s i (s 1 ) u is (s 1 ) r (s 1 ) -ρs i (s 0 ) ρ s i (s 0 ) u is (s 0 ) r (s 0 )] + ∂Ω    ∂ ρs i ∂s 1 ∂r ∂s
2 + ∂z ∂s 2 Boundary conditions associated to this equation are

ρ s i u is ds - ∂Ω r ρs i 1 τ * ρ s i ds + ∂Ω r ρs i 1 τ * eq ρ eq i ds, t ∈ [0, T ] (8)
ρ s i (s 0 ) = ρ eq 0 ( 9
)
ρ s i (s 1 ) = ρ eq 1 ( 10
)
RR n 6175

These surface equations are fully coupled with the Navier-Stokes equations (macroscopic model) through the generalized Navier slip boundary conditioins presented in next section, see [START_REF] Shikhmurzaev | The moving contact line problem on a smooth solid surface[END_REF], [START_REF] Blake | Experimental evidence of nonlocal hydrodynamic influence on the dynamic contact angle[END_REF] , [START_REF] Monnier | Modélisation Numérique de la Ligne Triple[END_REF], [START_REF] Monnier | Analysis of a Local Hydrodynamic Model with Marangoni Effect[END_REF] .

An analogy with the surfactant equations

The surface equatiions arising in the Shikhmurzaev's theory are similar to those modelling the dynamics of surfactants. We present in this section the surfacatants equations and their weak form. Let us denote Γ the surface concentration of the surfactant, mesured in moles per surface area. According to [START_REF] Slatery | Interfacial transport phenomena[END_REF], the surfactant mass balance over an infinitesimal section of the interface is governed by a convection-diffusion equation, in the presence of a distributed source,

d s Γ dt + ∇ s • (Γ u s ) -w s • ∇ s Γ + Γκ ( u • n) = D s ∇ 2 s Γ, (t, s) ∈ [0, T ] × [0, 1] (11)
u s is the tangent to the interface component of the fluid velocity, i.e.

u s = n × ( u × n) (12)
The derivative d s /dt expresses the rate of change of a variable following an interfacial marker point moving with the velocity w s . In other words

d s dt = ∂ ∂t + w s • ∇ s ( 13
)
w s is definetely arbitrary and it is related to the markers velocity to the interface by the relation

v = u • n + w s (14)
When w s is taken to be equal to u s , the markers will represent material fluid particles.

As done for Shikhmurzaev equation, equation [START_REF] Shikhmurzaev | The moving contact line problem on a smooth solid surface[END_REF] has to be written in the surface coordinate system (s, φ). According to (77), (78), (79) and [START_REF] Qian | Generalized Navier boundary condition for the moving contact line[END_REF], this is straightforward and gives

∂Γ ∂t + w s • ∇ s Γ + 1 σ ∂ ∂s (σΓu s ) -w s ∂Γ ∂s + Γκ( u • n) = D s 1 σ ∂ ∂s σ 3 r 2 ∂Γ ∂s , (t, s) ∈ [0, T ] × [0, 1] (15)
INRIA σ being given by expression [START_REF] Monnier | Flow patterns in the vicinity of triple line dynamics arising from a local surface tension model[END_REF]. Equation (15) has to be multiplied by the test function Γ associated to Γ and integrated over the whole free interface so to obtain

Thus, if we manage to implement and solve correctly the Shikhmurzaev's equations, next we can consider the coupling with the present surfactant model.

Droplet geometry and mesoscopic boundary conditions

In the Shikhmurzaev's theory, [START_REF] Blake | Experimental evidence of nonlocal hydrodynamic influence on the dynamic contact angle[END_REF] , the previous "local surface model" of the Shikhmurzaev theory are fully coupled with the present boundary conditions for the Navier-Stokes equations . We refer to [START_REF] Blake | Experimental evidence of nonlocal hydrodynamic influence on the dynamic contact angle[END_REF] , [START_REF] Monnier | Modélisation Numérique de la Ligne Triple[END_REF], [START_REF] Monnier | Analysis of a Local Hydrodynamic Model with Marangoni Effect[END_REF] for more details.

Let us consider the following droplet geometry. The generalized boundary conditions on the free surface (liquid-gas) are:

Σ n = (-p ext + κ Ca ) n + h τ in [0, T ] × Γ f ( 19
)
where Ca is the Capillary number, κ is the mean curvature and p ext is the external pressure.

The extra term h is given (it is a surface tension gradient). We have classical boundary conditions on the symmetry axis. We decompose the liquid-solid interface into two parts Γ ad and Γ sl . Γ sl denotes a "small" part of the liquid-solid interface near the triple point. We consider a generalized Navier slip type boundary condition (local slipping)

u. n = 0 in [0, T ] × Γ sl Σ τ = -[β u + g]. τ in [0, T ] × Γ sl ( 20
)
where β is a given sliding coefficient. The extra term g is given, it models either a surface tension gradient in the Shikhmurzaev's theory or a uncompensated Young stress in the Qian-Wang-Sheng theory, see below. On Γ ad , we impose adherence boundary conditions u = 0.

Contact line dynamics modeling

We consider two different types of model for the contact line dynamics an explicit model (Tanner type law) and an "implicit" one deriving from the Shikhmurzaev theory.

Mobility relation (Tanner type Law). The contact line velocity and the wetting angle are related by

U CL (t) = k (θ(t) -θ eq ) ι (θ eq ) ι for t ∈]0, T [ (21)
where U CL is the contact line velocity, θ eq is the wetting angle at equilibrium (Young's law), θ is the (dynamic) wetting angle. k and ι are parameters usually determined using

INRIA experimental results.
Local flow modelling and Shikhmurzaev's theory. The "implicit" model does not impose the contact line velocity and the wetting angle but consider them as a response of the model. This model is based on the Shikhmurzaev's theory, [START_REF] Shikhmurzaev | The moving contact line problem on a smooth solid surface[END_REF], [START_REF] Blake | Experimental evidence of nonlocal hydrodynamic influence on the dynamic contact angle[END_REF], which makes introduce the generalized Navier slip condition (20) and the condition (19). These conditions are local since the extra terms g and h vanish except in a vicinity of the triple point. The basic idea of this theory is to consider that the rolling motion observed in experiments, [START_REF] Dussan | On the motion of a fluid-fluid interface along a solid surface[END_REF], implies that particles of the liquid-gas interface become an element of the solid-liquid in a finite time. Then, the surface tension value associated to this particle must change to its new equilibrium value relative to the solid-liquid interface. This process would gives rise to a surface tension gradient in a small vicinity of the advancing contact line (hence a local Marangoni effect). In other respect, the Young equation would remain valid at any time. In this theory,

g = - 1 2Ca ∇γ and h = 1 Ca ∇γ LG . τ (22)
where γ and γ LG are the liquid-solid and the liquid-gas surface tension coefficient respectively. In [START_REF] Monnier | Analysis of a Local Hydrodynamic Model with Marangoni Effect[END_REF] and [START_REF] Monnier | Flow patterns in the vicinity of triple line dynamics arising from a local surface tension model[END_REF], a mathematical and numerical study presents some qualitative behaviors of g and h arising from Shikhmurzaev's theory.

A connection with Qian-Wang-Sheng theory. From molecular dynamics simulations on immiscible fluids, [START_REF] Qian | Generalized Navier boundary condition for the moving contact line[END_REF] shows that the relative slipping between the fluid and the solid wall follows a generalized Navier slip b.c. similar to (20) (if one phase fluid only is considered). In this theory, the extra term g in (20) would model an interfacial uncompensated Young stress. The extra tangential stress g is defined as follows Γint g dy = γ(cosθ-cosθ eq ), where Γint dy denotes the integral across the interface Γ f .

Basic principles of the ALE method

In this section is presented the basic principles of the the so-called Arbitrary Lagrangian Eulerian (ALE) method for tracking moving boundaries. The computational grid therefore follows the physical interface. Although other methods have proven their high efficiency to simulate fluid flows with moving interfaces (level set method, VOF method, front-tracking method), the choice of the ALE approach is motivated here by two arguments. First, no topological changes are expected to happen in the problem we are dealing with. Second, a high precision to achieve on the position of the interface is of major interest in triple line simulations, and one must impose some extra physical equations on this surface (surface equations arising from the Shikhmurzaev's theoery or from surfactant modelling). Thus, we believe that this can be achieved more easily and more accurately with ALE methods rather than with a level set approach for instance.

The ALE formulation for one fluid with a free interface has been widely used by many authors, see e.g. [START_REF] Hughes | Lagrangian-Eulerian Finite Element Formulation for incompressible Viscous Flows[END_REF], [START_REF] Maury | Characteristics ALE method for unsteady free surface flows with surface tension Z[END_REF].

We now introduce some notations. For t > 0, we define a reference configuration Ω. The liquid-gas interface is denoted Γ LG and the solide-liquid interface is denoted Γ SL . The normal n to each interface is directed outward the liquid domain. We define a reference configuration Ω = Ω t=0 and we consider a mapping Ãt which associates to a point x ∈ Ω a point x ∈ Ω. Throughout, the mapping Ãt will be supposed to be smooth egough in x, invertibble with a smooth inverse, and differentiable with respect to time, t. For a function ψ(•, t) defined on Ω, we denote by ψ(•, t) the function defined on ψ satisfying ψ(x, t) = ψ(x, t), with x = Ãt (x). By a classical abuse of notation, we will denote ∂ψ ∂t | x the time derivative on the ALE frame

∂ψ ∂t x (x, t) = ∂ ψ ∂t (x, t), with x = Ã-1 t (x) (23)
The domain velocity v is defined by

v(x, t) = ˜ v(x, t) = ∂A t ∂t (x), with x = Ã-1 t (x) (24)
Since, according to (14), we have

u • n = v • n on ∂Ω, (25)
which means that the domains Ω t follows the free interface. Thus,

for x ∈ Ω, x = Ãt (x) ∈ Ω t (26)
Although the domains Ω and Ω are in fact the same domain of R d (d is the space dimension), we will keep both notations for the sake of clarity. We now recall some standard formulae whith will be useful in the sequel. In other words, all the differential operators (∇, ∇•, ∇ 2 ) will be take with respect to the Eulerian variable x. We have

∂ψ ∂t x (x, t) = ∂ ψ ∂t (x, t) = ∂ψ ∂t (x, t) + v • ∇ψ(x, t) (27)
We denote by Jt the Jacobian matrix of Ãt ,

Jt = ∂ Ãt ∂ xj
and by Jt its determinant. Then we have the Euler formula,

∂ Jt ∂t Ã-1 t (x), t = Jt Ã-1 t (x), t ∇ • v(x, t) (28)

INRIA

Using this relation, we get in particular

d dt Ωt ψ(x, t)dx = d dt Ωt ψ(x, t) Jt dx = Ωt ∂ ∂t ψ(x, t) Jt dx = ωt ∂ ψ ∂t (x, t) Jt + ψ(x, t) Jt ∇ • v dx (29)
and thus

d dt Ωt ψ(x, t) dx = Ωt ∂ψ ∂t x dx + Ωt ψ(x, t)∇ • v dx (30)
This equation will be used when replacing the mapping ψ by any defined given physical mapping (velocity, pressure, etc.). We now propose a weak formulation of equations ( 3) and ( 4). The following functional spaces will be needed

H 1 0 (Ω) = w ∈ H 1 (Ω) d , P 1 0 (Ω) = q ∈ H 1 (Ω) . We define: W = L 2 (0, T ; H 1 0 (Ω)), P = L 2 (0, T ; P 1 0 (Ω)).
In order to express the weak form associted to the mass and momentum equations in the current frame, we introduce the test function spaces defined on the reference frame: W = H 1 0 ( Ω), P = H 1 0 ( Ω). On the current domain, the test function spaces are defined by

W 0 = { w : Ω t × [0, T ] → R d , w(x, t) = ˜ w Ã-1 t (x) , ˜ w ∈ W }, P 0 = {q : Ω t × [0, T ] → R, q(x, t) = q Ã-1 t (x) , q ∈ P }.
It is worthwhile to notice that the test functions do not depend on time on the reference frame whereas they do on the current one. More precisely, denoting by (w i ) i=1,...,d the components of w ∈ W 0 , we have

∂w i ∂t x = ∂w i ∂t + v • ∇w i = 0 (31)
Similar relation holds for the functions q ∈ P 0 . We now give the formulation that will be used in the numerical simulations.

Weak ALE formulation. Suppose that the domain is moving such relation (14), beyong others, is satisfied. We look for ( u, p) in W × P such that, ∀( w, q) ∈ W 0 × P 0 :

d dt Ω ρ u • w dx + Ω ρ( u -v) • ∇ u • w dx - Ω ρ u • w∇ • v dx = Ω -∇p + µ∇ 2 u + ρ f • w dx (32) Ω q∇ • u dx = 0 (33)
Proposition System (3)-( 4) is formally equivalent to the weak ALE formulation (32)-(33).

Proof. Let ( u, p) be a solution to (3)-( 4). For the momentum balance equation (32), using formula (27), we have from ( 3)

ρ ∂ u ∂t x + ( u -v) • ∇ u = h
where h = -∇p + µ∇ 2 u + ρ f . By multiplying this equation by w ∈ W 0 and integrating over Ω, we get

Ω ρ ∂ u ∂t x • w dx + Ω ρ( u -v) • ∇ u • w dx = Ω h • w dx (34)
using ( 30) and (31), we can write the first term as follows

Ω ρ ∂ u ∂t x • w dx = Ω ρ ∂ ∂t x ( u • w) dx = d dt Ω ρ u • w dx - Ω ρ u • w∇ • v dx (35)
Thus, the eventual momentum equation in the weak form reads

d dt Ω ρ u • w dx + Ω ρ( u -v) • ∇ u • w dx - Ω ρ u • w∇ • v dx = Ω h • w dx (36)
In addition, multiplying by q ∈ P 0 and integrating over Ω, the weak equation associated to the mass conservation equation ( 4) reads as given by (33). Conversely, and independently of boundary conditions, if ( u, p) is a solution to the weak ALE formulation (36)-(33) for all ( w, q) ∈ W 0 × P 0 then ( u, p) is a solution to the system (3)-( 4). The demonstration is however out of the scope of this paper.

Mesh velocity law. All we did so far is to specify the mesh velocity v by the condition (14), where w s is an arbitrary tangential to the interface velocity. Perhaps the more intuitive reaction is to take this later equal to zero, that is, no mesh sliding is allowed at the boundaries. Even if it is shown it has a good mass conservation property, using normal velocity component for free interfaces is however not always of benefit on surfaces that are not supposed to deform. Let us suppose, for any reason, we have made our mind up about this surface tangential velocity component. The next question is how will the mesh deform INRIA inside the domain? The answer is given by the name of the method: arbitrarily (once more). One can take virtually any equation to be satisfied on Ω. Virtually any one and not any one, since any point from the domain Ω has to remain inside its (moving) boundary ∂Ω as time goes on. A simple way to satisfy this condition is to choose, among others, an elliptical equation to be satisfied by the mesh velocity on Ω. More precisely, a better choice is to have the mesh velocity in such a way to preserve a quality equivalent to this of the initial mesh during the whole process of deformation. Elasticity equations are a common and efficient way to achieve such a goal.

We refer to [START_REF] Monnier | Flow patterns in the vicinity of triple line dynamics arising from a local surface tension model[END_REF] for a description of the ALE algorithm implemented into Micralef, computational code used in the sequel.

Comsol Multiphysics 3.2 and its incompatibilities with droplet modelling

A general presentation of COMSOL Multiphysics is made in Appendix B. In what follows we focus our attention on the possibilities given by the ALE new mode that we chiefly are concerned. We present below briefly the ALE implementation in Comsol Multiphysics 3.2.

Then some tricks useful to overcome some difficulties when using Comsol Multiphysics 3.2 (difficulties inherent in this software as being not a fully open package).

In summary, unfortunately the present implementation of ALE method into Comsol Multiphysics (version 3.2) does not allow to consider large mesh distortions which necessarily appear when modeling free surface flows (such as a droplet impact). Furthermore, in order to simulate a free surface flow driven by surface tension forces one must be able to compute the curvature forces. The present implementation in Comsol Multiphysics prevent to do it straightforwardly. Also, we encoutered serious difficulties to implement an algorithm to treat the triple line dynamics. Nevertheless, we present below some tricks useful to overcome some of these difficulties.

For the description of an efficient implementation of the ALE method and the description of a triple line dynamics algorithm as required for a free surface droplet flow, we refer to [START_REF] Monnier | Flow patterns in the vicinity of triple line dynamics arising from a local surface tension model[END_REF]. The latter approach has been implemented into our home-developed software Micralef, which led to the preliminary numerical results presented in next section.

The ALE implementation in Comsol Multiphysics 3.2

The ALE algorithm is implemented as a new physics mode strongly coupled with the original mechanical problem since Navier-Stokes equations are expressed in a domain which is determined by the ALE mode and the ALE physics is completed by kinematic boundary conditions supplied by Navier-Stokes equations mode. By default the mesh velocity is governed by Laplace equation:

∇ 2 v = 0 ( 37
)
which is completed with suitable boundary conditions

v| interf ace = ( u • n) n| interf ace (38)
COMSOL Multiphysics also offers the possibility to implement one's own mesh governing equations but many difficulties arose to do so in practice.

As mentioned previously, one must define the mesh velocity in such a way to preserve a good quality of the mesh during the whole process of deformation. An efficient way to acheive this goal is to define the mesh velocity by solving a simple linear elasticity system. Suprisingly, this is has not been foreseen into the ALE new mode of Comsol Mulitphysics. Comsol uses instead the Laplace equation or Winslow equations which is not suitable for large distorsions as this appears in the present application (and in free surface flows in general). Such an approach is probably robust enough for fluid-structure interactions since in that case, mesh distortions are much smaller, but not in our case. Furthermore, in order to implement the present mathematical model (eg the curvature term , and the dynamics of the triple point), one needs to define and / or obtain some essential informations. We present below how we circumvented the standard use of Comsol in order to reach (partially only) our goal.

Some tricks to overcome difficulties encountered

Comsol Multiphysics is reputed to be significantly an open scientific software, although some information are not straightforward or even impossible to obtain. Among these are the curvature of the boundaries, the maximum of a given space-dependent function and the current time step. In the following we explain the different stategies we adopted to get most of those quantities.

Curvature computation

If for some reason we do need to compute the curvature, or other boundary differential quantities, and not its weak integral over a given boundary, the simplest (not the most elegant, though) way is to solve simultanously the following trivial physics problems in the weak form (for simplicity, we suppose we are dealing with a 2D configuration using the co-ordinates r,z):

coordr (coordr -r) = 0 (39) coordz (coordz -z) = 0 (40)
coordrT zz = ∂ 2 coordr ∂s 2 tz 2 (48)
The quantities corresponding to coordz are given in a similar manner. Hence the planer (first principal) curvature is obtained by

κ 1 = ∂ 2 r ∂s 2 ∂z ∂s -∂ 2 z ∂s 2
∂r ∂s

∂r ∂s

2 + ∂z ∂s

2 3/2 (49)
Bearing in mind that due to axisymmetry, the second principal curvature is given by

κ 2 = n z r (50)
the total curvature, given by

κ = κ 1 + κ 2 , (51)
is then definetely determined.

Capillary effects

To introduce forces due to capillary effects, which are surface-located, into the volumic kinetic balance, two ways are conceivable. The first is to compute directly the curvature as given by the previous subsection 4.2.1 and take into account the term κγ n in the surface balance applied to the interface boundary condition. Let us show how the surface balance is injected into the COMSOL interface boundary condition. For accommodation, we rewrite the normal projection of equation ( 19), when the external pressure is equal to zero and no tangential trailing is imposed (we delay this treatment for a next step), in the form

Σ n -κγ = 0 (52)
By default, COMSOL does not propose such a boundary condition. To make it support this condition, we have to relax the boundary condition, by imposing "initially" neutral condition, and then use the Lagrange multipliers technic in order to constrain the desired condition, namely equation ( 52), to be satisfied. Let λ be the Lagrange multiplier we are using. The equation this new unkown has to satisfy (and we have to introduce) is

λ (Σ n -κγ) -λ u • n = 0 (53)
We remind that ˆ(hat) symbol denotes the associated test function of the underlying quantity. This equation expresses only the fact that the condition Σ n -κγ is applied to the normal projection of the momentum equation by means of the associated quantity u • n.

The resolution of fluid dynamic problem is hence done with boundary conditions (at least some of them) handled as weak constraints, since equation ( 53) is by definition a weak form.

Maximum of a space-dependent function

Perhaps one of the most tedious tasks was to get the maximum of a space-dependent function. Let us first explain why have we to compute such a value and then how will we do it.

In the case of a spreading droplet over a solid surface, we have to apply different boundary conditions whether solid-liquid surface or liquid-gas interface are concerned. For instance on the liquid-gas interface:

Σ n -κγ = 0 (54)
and on the solid-liquid surface: The point is during the droplet spreading a given part of the liquid-gas interface will rolls up and sticks up on the solid-liquid interface (see figure 5). This what we call projection of the liquid-gas interface. So, at time t, the boundary conditions applied on the red element are given by (56) while at time t+ ∆t they will be given by ( 54) and ( 55). Two questions have to be answered. The first is how to switch between one type of boundary conditions to another.

u • n = 0 (55) Σ t = 0 (56) INRIA
The second is what is the criterion that permits to choose one boundary condition rather than the other. We will answer to the second question first. The criterion is merely belonging to the solid-liquid interface or not. In other words, a point belong to the solid-liquid surface if it satisfies:

r -r p < 0 (57) z -z p < ǫ (58)
where (r p , z p ) are the coordinates of the triple-point and ǫ is a tolerence much smaller than the smallest mesh size. Now the answer to the first question. The easiest way in COMSOL Multiphysics to switch from a boundary condition to another is to use weak-constraint handling through Lagrange multipliers in the following way:

λ [(Σ n -κγ) H(z -z p -ǫ) + ( u • n) H(-z + z p + ǫ)] -λ u • n = 0 ( 59
)
where H is the Heaviside function (or a smoothed function that behaves alike). The triplepoint coordinate, given by (r p , z p ), have to be evaluated. The triple line z-co-ordinate is obtained in straightforward manner, while the r-co-ordinate has to be conputed in the following way:

r p = max {r/P (r, z) ∈ Γ ∧ z -z p < ǫ} ( 60
)
Unfortunately there is no way to directly obtain a maximun of a space dependent-function using COMSOL Multiphysics. However an indirect tricky (and time consuming) way to do so exists using the following property of positive functions (a < b):

lim n→+∞ 1 b -a b a |f (x)| n dx 1/n = max x∈[a,b] (|f (x)|) ( 61
)
In our case the function to be maximized along the concerned boundaries is rH(z p -z + ǫ). From a numerical point of view, and according to (61), the larger n is the better the maximum is computed. A bad news is that this argument is no longer correct because round off error accumulated during computing increases as n gets larger. The other bad news is that error increases with respect to n before getting a good approximation of the desired maximum, as shown by figure 6.

Figure 6. Variations of the relative error with respect to the power n used in expression (61). Now the good news (even here there are some). Another way to get the maximum is to use the expression (61) with a relative small value of n not once only, but many times in a iterative way. Figure 7 shows the variations of the relative error with the iteration number when n = 64. Unfortunately, we did not manage to obtain neither large mesh distortions nor valid result of the spreading of the droplet when using Comsol Multiphysics 3.2.. Therefore, we do not present any figures obtained.

INRIA

5 Numerical results obtained using Micralef software: spreading phase

In the present section, we present some preliminary numerical results of the spreading of a droplet impacting a solid surface. These results are obtained when solving the model presented previously (see also [START_REF] Monnier | Flow patterns in the vicinity of triple line dynamics arising from a local surface tension model[END_REF] for more details). The numerical scheme and algorithms have been implemented into Micralef, a home-developed software.

We consider a water droplet. The reference length and velocity are:

L ref = 2.3 mm, U ref = 0.98 m/s. Then, Re = ρU ref L ref µ ≈ 46, Weber number W e = ρU 2 ref L ref σ ≈ 68 and C a = σ µU ref ≈ 1.5.

Spreading phase using a Tanner type law

We start by testing the efficiency of the algorithm implemented into Micralef. To this end, we consider the Tanner type law as triple line dynamic modelling, and we make fit this law with available experimental data, see e.g. [START_REF] Allaman | Impact and spreading of microdrops onto solid substrates[END_REF], related to a spreading phase. The triple point position is imposed at each time. We focus on the volume conservation, the height and diameter of the spreading (splat radius) and the deformation of the mesh (number of re-meshing necessary), see Fig. 5.1. With a time step dt = 5.10 -5 s and a coarse mesh (≈ 600 elements), volume lost after 650 iterations is roughly 3%. The loss occurs mainly at the very beginning of spreading. Mesh transport is efficient hence re-meshing occurs only when a projection occurs (see previous section). From a qualitative point of view, numerical results are similar to experimental data presented in [START_REF] Allaman | Impact and spreading of microdrops onto solid substrates[END_REF].

Spreading phase using an implicit triple line dynamics

We consider the spreading phase using the Shikhmurzaev's model and the algorithm presented previously. The source term g in the generalized Navier slip condition acts like a control on the position of the triple point. A decrease for g < 0 translates in a faster spreading. If g is small, this influence is likely not to be seen in the first part of the spreading, when the inertial forces dominate all other forces in presence. By increasing |g| over a threshold, which depends on β, it accelerates the spreading and modifies the flow nature. For |g| large, large interface distortions occurs and mesh generator fails to re-mesh.

In order to observe the influence of the g-term, we fix β = 1000, and we perform 1000 iterations with a step size of dt = 5.0 × 10 -6 for g ∈ {-10, -100}, Fig. 5.2. One can observe that after the inertial phase, the triple point position is farther to the right with an increasing |g|. Fig. 5.3 shows a droplet subjected to large distortion (similar to breakup) when |g| is getting larger. Simulation can go no further because of the change of topology it implies is not supported by the ALE method we are using.

A first attempt to simulate the recoiling phase. We try now to see the influence of the sign of parameter g on the droplet behavior. Since we have seen that parameter g monitor the spreading, it is also supposed to do so with recoiling. In Fig. 5.4, simulation is started with initially spread droplet. Also here, the ratio -g β significantly influences the behavior of the droplet while parameter h influences more particulary on the droplet curvature than on the triple-point position.

Conclusion

We presented the Navier-Stokes free surface equa-tions with a generalized Navier slip condition as it appears both in the Shikhmurzaev theory and in the Qian-Wang-Sheng theory. Furthermore, we wrote the surface equations arising in the Shikhmurzaev theory in the surface co-ordinate system and in weak form. An anlmogy with the surfactant equations si done; and these equations are detailled. Then, we consider as test case the droplet impact onto a solid surface. An attempt to model and simulate using Comsol Multiphysics 3.2 is done. Unfortunately, we show that the ALE implementation in Comsol Multiphysics 3.2 is not robust enough to handle large mesh distorsion (as it appears in many surface flows). In addition, some important quantities such as curvature, are not easily computable. We present some tricks to circumvent this latter problem. Finally, we present some preliminary numerical results using our freeware Micralef. Miocralef has been developed especially to solve the present equations and the algorithms we elaborated (ALE, spreading and recoiling algorithms). We obtain qualitatively good results when siumulating the spreading phase. An efficient numerical model for the recoiling phase is still under progress. Droplet profiles and triple point position (right) for (β, g, h) = {(500, 100, 0), (100, 100, -10)} .

INRIA

A Surface differential operators expressions

Let the surface S be affected by an arbitrary co-ordinate system (u 1 , u 2 ) and let ϕ be a scalar function and A α a contravariant vector field both defined in S. The metric tensor of S is defined by the expression, Aris,

a αβ = 3 i=1 ∂y i ∂u α ∂y i ∂u β , y ∈ S ( 62
)
y i is the i-th co-ordinate of the current point of S, expressd in the Reimannian space E 3 . It is shown that the metric tensor a αβ is a symmetric covariant second order tensor. We define the short-hand determinant of a αβ by

a = det|a αβ | = a 11 a 22 -(a 12 ) 2 (63)
The Kronecker delta δ α β which is 1 if α = β and zero otherwize, permits to define the conjugate metric tensor by

a αγ a γβ = a γα a γβ = a αγ a βγ = δ α β (64)
1. The surface divergence of A α is given by

∇ s • A = A α ,α = 1 a 1 2 ∂ ∂u α a 1 2 A α (65)
2. The surface gradient of ϕ is given by

∇ s ϕ = ϕ ,α = ∂ϕ ∂u α e α ( 66
)
where e α is the α-th tangent vector to S.

3. The surface Laplacian of ϕ is given by

∇ 2 s ϕ = a αβ ϕ ,αβ = 1 a 1 2 ∂ ∂u β a 1 2 a αβ ∂ϕ ∂u α (67)
4. The surface curl of A i is given by

∇ s × A = ǫ αβ A α,β n = a -1 2 ∂A 2 ∂u 1 - ∂A 1 ∂u 2 n ( 68
)
ǫ αβ is the (two dimensional) permutation tensor defined by

ǫ 11 = ǫ 22 = 0, ǫ 12 = 1 and ǫ 21 = -1 (69)
and n is the normal vector to S. 5. The elementary surface is given by

ds 2 = a αβ du α du β (70)

INRIA

As surface equations are involved, the whole issue is first to choose the surface co-ordinate system that maps the surface S and second to get the associated metric tensor a αβ in order to express correctly each differential operator according to the mapping co-ordinates as discribed formerly.

Application: axisymmetric surface S In the case where u 1 , u 2 = (s, φ), as shown in figure 1, and considering the axisymmetry of the geometry, then we can write

x = r cos φ (

B Short presentation of Comsol Multiphysics 3.2 software

This section is dedicated to describe the numerical tool used to implement and simulate the triple line dynamics, to name COMSOL Multiphysics 3.2b. COMSOL Multiphysics is an interactive environment for modeling and solving many kinds of scientific and engineering problems based on partial differential equations (PDEs). One accesses the power of COMSOL Multiphysics as a standalone product through a flexible graphical user interface, or by script programming in the COMSOL Script language or in the MATLAB language.

As noted, the underlying mathematical structure in COMSOL Multiphysics is a system of partial differential equations. Three ways of describing PDEs are provided through the following mathematical application modes:

1. Coefficient form, suitable for linear or nearly linear models;

2. General form, suitable for nonlinear models;

3. Weak form, for models with PDEs on boundaries, edges, or points, or for models using terms with mixed space and time derivatives.

Using these application modes, one can perform various types of analysis including:

1. Stationary and time-dependent analysis 2. Linear and nonlinear analysis

Eigenfrequency and modal analysis

INRIA

When solving the PDEs, COMSOL Multiphysics uses the finite element method (FEM).

The software runs the finite element analysis together with adaptive meshing and error control using a variety of numerical solvers. PDEs form the basis for the laws of science and provide the foundation for modeling a wide range of scientific and engineering phenomena. Therefore one can use COMSOL Multiphysics in many application areas.

The time-dependent solver As mentioned before, COMSOL Multiphysics proposes many solver modes, such stationary linear, stationary nonlinear, eigenvalue and time-dependent solvers. Each solver uses a specific algorithm. The time-dependent solver algorithm, which we are burdened with, and the setting parameters are presented in this section. To use the time-dependent solver, one has to go to the Solver Parameters dialog box, find the Solver list, and select Time dependent. Specifying output times. In the Time stepping area, one has to enter the output in the Times edit field, which contains a vector of times at which one wants the solution to the time-dependent model. The matlab syntax 0:0.1:1 represents a vector of times starting at 0 steps of 0.1 up to 1. The relevant time span depends on the model's dynamics.

Setting tolerances. To specify the accuracy and sometimes influence the performance, one has to set absolute and relative tolerance parameters for the time-dependent solver. In the Relative tolerance edit field, one can enter a positive number (default = 0.01). In the Absolute tolerance edit field, one can enter a single positive number (default = 0.001) or a space-separated list whose entries alternate between degree-of-freedom names and positive scalars. More precisely, let U be the solution vector corresponding to the solution at a certain time step, and let E be the solver's estimate of the local error in U commited during the current time step. The step is accepted if

   1 N i   E i A i + R |U i |   2    1/2 < 1
where A i is the absolute tolerance, R is the relative tolerance and N is the number of degrees of freedom.

The time-dependent solver algorithm. The finite element discretization of the time-dependent system of partial differential equations problem is 0 = L U, U , Ü , t -N (U, t)

t Λ 0 = M (U, t)
which is reffered to as the method of lines. Before solving this system, the algorithm eliminates the Lagrange multipliers Λ. If the constraints 0 = M are linear and time independent, the algorithm also eliminates them from the system. otherwise it keeps the constraints, leading to a algebraic-differential system. To solve the above ordinary/algebraic -differential equations system, COMSOL Multiphysics uses a version of the Differential Algebraic equation Solver Package, DASPK, originally proposed by L. Petzolds. This algorithm uses variable-order variable-stepsize backward differential formulae. Thus the solver is an implicit time-stepping scheme, which implies that it must solve a possibly nonlinear system of equations at each time step. It solves the nonlinear system using Newton iteration, and it then solves the resulting sytems with an arbitrary COMSOL Multiphysics linear system solver. The linearization of the above system used in the Newton iteration is

E V + D V + KV = L -N t Λ N V = M
where E = -∂L/∂ Ü is the mass matrix, D = -∂L/∂ U is the damping matrix, K = -∂L/∂U is the stiffness matrix. When E = 0, D is called the mass matrix.

C Short presentation of our home-developed software Micralef

Micralefe is a C++ finite element software (home-developed) solving 2D axisymmetric Navier-Stokes flow with free surface dynamics and curvature forces. It treats of a droplet impact on a solid substrate. The dynamics of the free surface is described using an ALE formulation (Arbitrary Lagrangian Eulerian), while the dynamic triple line modelling is based on the Shikhmurzaev's theory.The code is based on the Rheolef C++ finite element library. Micralefe capabilities include the dynamics of the contact angle, the triple point dynamics (liquid-gas-solid), variable surface tension coefficients, slip type boundary conditions and automatic mesh refinement. It has been developed in order to obtain numerical results presented in the present article.

Figure 2 .

 2 Figure 2. Mesoscopic domain ∂Ω, near the triple point, over which Shikhmurzaev equation is applied.

Figure 3 .

 3 Figure 3. A droplet put down on a solid surface.

Figure 1 .

 1 Figure 1. Switch from free boundary condition to sliding boundary condition, as the droplet spreads and its free interface rolls up.

Figure 7 .

 7 Figure 7. Variations of the relative error with the number of iteration of expression (61)with n = 64.

Figure 5 . 1 :

 51 Figure 5.1: Spreading: droplet profiles and volume in function of time steps

Figure 5 . 2 :

 52 Figure 5.2: Droplet profiles and triple point position (right) for (β, g) = {(10 3 , -10), (10 3 , -10 2 )} .

Figure 5 . 3 :

 53 Figure 5.3: Large interface distortion, similar to a droplet breakup g -10 4 } .

  Thus, the variable coordr (resp. coordz) should give an approximation of the r (resp. z) coordinate. Once evaluated, the first and second order space derivatives are obtained according

	to
	∂r ∂s	=	∂coordr ∂s	= coordrT r × tr + coordrT z × tz	(41)
	∂z ∂s	=	∂coordz ∂s	= coordzT r × tr + coordzT z × tz	(42)
			∂ 2 r ∂s 2 =	∂ 2 coordr ∂s 2	= coordrT rr + coordrT zz	(43)
		∂ 2 z ∂s 2 =	∂ 2 coordz ∂s 2	= coordzT rr + coordzT zz	(44)
	where tr and tz are the components of the tangential vector to the interface, and coordrT r,
	coordrT z, coordrT rr, coordrT zz are COMSOL Multiphysics syntax for the (output) quan-
	tities
					coordrT r =	∂coordr ∂s	tr	(45)
					coordrT z =	∂coordr ∂s	tz	(46)
			coordrT rr =	∂ 2 coordr ∂s 2 tr 2	(47)
						INRIA

  The code modularity allows us to plan further developments such as coupling with a surfactant model or thermal model. Micralef is avalaible under request , see URL: http://ljk.imag.fr/membres/Jerome.Monnier. The C++ source codes are commented while a short user's documentation is available. Potential users must be expert in computer sciences (C++ programming, compilers, Linux and open-source libraries).

	INRIA

RR n 6175

Benselama & Monnier

Acknowledgements. The authors would like to thank P. Pham, A. Glière from LETI -CEA-Grenoble, France, for their numerous and fruitful discussions and remarks. This study has been co-funded by a research program led by Prof. A. Soucemarianadin (UJF); the authors thank him very much.

INRIA


Minh Châu

Phan Huy

Philippe Delaye

Gilles Pauliat

Sylvie Lebrun

Efficient Raman converter in the yellow range with high spatial and spectral brightness

Keywords: Stimulated Raman scattering, hollow core photonic crystal fibers, wavelength converter

We present a Raman converter emitting at 583 nm on the second Stokes order of a line of propan-2-ol pumped by a microlaser at 532 nm in the sub-nanosecond regime. We used a mixture of liquids to adapt the transmission band of a photonic bandgap fiber. The internal conversion efficiency is 67% in photon numbers and the output power is 1.06 mW, corresponding to a maximum peak power of 338 W. The beam delivered by the converter presents a Gaussian spatial structure and a high spectral brightness, typically more than 5 times higher than supercontinuum sources in this spectral range.

Introduction

Today there is a growing demand for pulsed visible and UV laser sources in biological applications using fluorescence such as bioimaging or DNA sequencing or environmental applications such as atmospheric monitoring. In these areas compact laser systems emitting in the microjoule range one or several "exotic" wavelengths with a good beam quality and a high spectral purity would offer new opportunities for instance in order to excite the increasing number of existing fluorophores [1] or for the realization of new microlidars [START_REF] Spinhirn | Micro pulse lidars[END_REF]. Solutions such as supercontinuum sources [3] or OPOs [4] are expensive, cumbersome and oversized when only a few fixed discrete wavelengths are useful. Laser diodes can emit over a large number of wavelengths in the visible range but have many drawbacks such as the poor beam quality and the low peak power in the pulsed regime, limited to a few tens of watts. Moreover to sample the visible spectrum one needs to provide several laser diodes at the expense of extra costs. Microlasers are diode pumped solid-state lasers emitting pulses in the microjoule range with a repetition rate that can vary from Hz to a few hundred of kHz. These lasers are very popular for their compactness and their user-friendly aspect but commercial ones are mostly working at 1064 nm. They can be followed by a non linear stage to emit at 532 nm and possibly at 355 nm. Extending these achievable wavelengths to "on demand" discrete visible wavelengths would be of great interest for non-laser-specialized researchers. Raman effect has been widely used to generate new wavelengths. Starting from a pump source at a fixed wavelength reaching higher wavelengths is quite easy since the phase-matching condition is automatically satisfied. The Raman medium can be a crystal, a gas or a liquid. With a crystal, the systems are quite complex and not versatile. For example, in [START_REF] Omatsu | Passively Q-switched yellow laser formed by a self-Raman composite Nd:YVO4/YVO4 crystal[END_REF] in order to generate a given wavelength in the yellow range the authors need to develop a special pump source in the infrared based on Raman effect at first and then insert a nonlinear crystal in the cavity for frequency doubling. Changing the emitted wavelength needs to change several components and a realignment of the whole system. The conversion efficiency remains quite low, typically a few percents. In order to consequently increase the density of pump power over large interaction lengths capillaries or special fibers filled with nonlinear liquids have been studied [START_REF] Fanjoux | Stimulated Raman-Kerr scattering in an integrated nonlinear optofluidic fiber arrangement[END_REF][START_REF] Yiou | Stimulated Raman scattering in an ethanol core microstructured optical fiber[END_REF]. Even if the total conversion efficiency of the pump towards the Stokes lines can be high the energy of the pump is spread over all these lines generated through Raman cascade effect.

In order to stop this cascade effect and to optimize the transfer of energy of the pump to a

given Stokes order we propose to use hollow-core photonic crystal fibers (HC-PCF) filled with liquids. In such fibers, light can propagate in the core over long distances on small areas with low losses even if the refractive index of the core is smaller than the one of silica [START_REF] Cregan | Single-Mode Photonic Band Gap Guidance of Light in Air[END_REF].

Multiline UV and visible Raman converters in hydrogen-filled HC-PCF have been fabricated [START_REF] Wang | Compact and portable multiline UV and visible Raman lasers in hydrogen-filled HC-PCF[END_REF]. In this reference, many Raman lines and Stokes orders are excited at the same time giving birth to 8 to 23 new wavelengths according to the fiber that is used. The use of liquids with refractive index different from unity offers another degree of freedom. In HC-PCF the photonic bandgap guidance opens spectral transmission bands that can be adapted through the control of the refractive index of the material that fills the holes [START_REF] Antonopoulos | Experimental demonstration of the frequency shift of bandgaps in photonic crystal fibers due to refractive index scaling[END_REF]. In previous works, we demonstrated the efficient conversion of the pump to a single Stokes order in a HC-PCF filled with ethanol by stopping the Raman cascade at the first Stokes order [START_REF] Lebrun | High-efficiency single-mode Raman generation in a liquid-filled photonic bandgap fiber[END_REF][START_REF] Lebrun | Optical characterizations of a Raman generator based on a hollow core photonic crystal fibre filled with a liquid[END_REF][START_REF] Lebrun | Efficient stimulated Raman scattering in hybrid liquid-silica fibers for wavelength conversion[END_REF]. The pump source emitted at 532 nm sub-nanosecond µJ pulses and the conversion efficiency at high pump power to the first Stokes order of ethanol at 630 nm was 70% in photon numbers [START_REF] Lebrun | Optical characterizations of a Raman generator based on a hollow core photonic crystal fibre filled with a liquid[END_REF].

In particular we showed that despite the liquid filled fiber was slightly multimode (the used fiber structure is optimized for good transmission properties when empty) an emission on a single quasi Gaussian mode was possible [START_REF] Lebrun | High-efficiency single-mode Raman generation in a liquid-filled photonic bandgap fiber[END_REF]. The real mode shape revealed the hexagonal symmetry of the cladding structure. We also showed that the temporal shape of the pulses can be perfectly explained by a quasi-steady state regime modelling of the Raman mechanisms i.e. instantaneous response of the Raman effect at the pulse duration time scale [START_REF] Lebrun | Optical characterizations of a Raman generator based on a hollow core photonic crystal fibre filled with a liquid[END_REF][START_REF] Lebrun | Efficient stimulated Raman scattering in hybrid liquid-silica fibers for wavelength conversion[END_REF]. We experimentally measured a long-term stability only limited by the pump beam stability [START_REF] Lebrun | Optical characterizations of a Raman generator based on a hollow core photonic crystal fibre filled with a liquid[END_REF]. All these characteristics have been observed for different liquids at different Stokes wavelengths [START_REF] Lebrun | Efficient stimulated Raman scattering in hybrid liquid-silica fibers for wavelength conversion[END_REF]. After these proofs of principle we propose to build in the present work an efficient Raman converter emitting in the yellow range at the target wavelength of 583 nm. Indeed this wavelength is interesting for biological applications since it matches an absorption line of several fluorophores, such as Bodipy 581/591 or X-Rhodamine for example [1]. We describe firstly the design of this new converter. The choice of the active Raman liquid and the mixture of liquids to be used to adjust the refractive index and the band position is discussed. Then we present the experimental setup and the results. The internal conversion efficiency from the pump at 532 nm to the target wavelength at 583 nm is high (67% in photon numbers) despite working on the second Stokes order of the Raman cascade. The new wavelength is emitted in a quasi-Gaussian spatial beam with a high spectral brightness.

Design of the Raman converter emitting at 583 nm

The Raman converter is based on the Amplified Spontaneous Emission principle, photons emitted through spontaneous Raman scattering being further amplified by stimulated Raman scattering. It is thus very simply implemented as requiring just a pump operating in a single pass scheme in the nonlinear medium. Liquids as any material present complex Raman spectra with many lines among which the one with the highest Raman gain coefficient gR1 is called the main Raman line. Other lines with lower Raman gain coefficients gRi are called secondary Raman lines. For a given line, the Raman threshold can be defined when the

parameter 𝛾 = 𝑔 $ % & ×( ) *++
is approximately equal to 25 [START_REF] Phan Huy | Lowering backward Raman and Brillouin scattering in waveguide Raman wavelength converters[END_REF]. gR is the Raman gain coefficient, Pp is the internal peak pump power, Aeff the effective area and L the interaction length. When increasing Pp and assuming no losses, before its saturation the Stokes power varies as 𝑒𝑥𝑝(𝛾).

So, around threshold, if

? @A ? @B ≾ 0.7 , the Stokes power of the main line will be more than 1800 times higher than the Stokes power of the secondary Raman line. So we propose this value of 0.7 as a good criterion to say that, as we are working in the steady-state regime, Stimulated Raman Scattering (SRS) occurs on the main Raman line only, secondary lines emission being neglected. Above threshold the pump begins to be depleted by energy transfer toward the main line and no more power is available to amplify the secondary lines, increasing again the effect. When ? @A ? @B becomes close to unity, it is difficult to identify the main Raman line from a secondary line. In this case competition occurs between the two lines and a solution to favor one line is to put extra-losses on the other one. Indeed the amplification gain of the unwanted line varies as 𝑒𝑥𝑝(𝛾 -𝛼𝐿) where 𝛼 represents the linear losses and hence vanishes when 𝛾 < 𝛼𝐿.

In this work we want to shift the pump wavelength at 532 nm to the target wavelength of 583 nm by using SRS in a liquid. Such a wavelength shift corresponds to a Raman shift around 1640 cm -1 which is not directly achievable with common and not harmful Raman liquids we want to use [START_REF] Martin | Infrared difference frequency generation[END_REF][START_REF]Spectral Database for Organic Compounds SDBS[END_REF]. It also corresponds to the second Stokes order of a line with a Raman shift around 820 cm -1 . This shift matches with a Raman line of alcohols. To choose the most convenient alcohol, we have to calculate the Raman gain coefficients. For that, we have used the formula demonstrated in [START_REF] Maier | Backward Stimulated Raman scattering[END_REF] and considered the depolarization ratio 𝜌 L which adds a correction factor M MNO P to it [START_REF] Porto | Angular dependence and Depolarization Ratio of the Raman effect[END_REF]. This formula becomes :

𝑔 $ = 1 1 + 𝜌 L 8𝜋 Δ𝜐 WXYY 𝑁𝑐 \ 𝜔 L ^ℏ𝑛 L \ a 𝜕𝜎 𝜕Ω e (1)
For the liquids and lines we use in this work, 𝜌 L is much smaller than one and brings only a minor correction, but has to be considered in some other liquids or lines. Δ𝜐 WXYY is the full width at half height of the spontaneous Raman line, N is the molecular density, c is the velocity of light, ws and ns are the Stokes frequency and refractive index.

fg fh is the total differential scattering cross-section [START_REF] Shan | Design of nanofibres for efficient stimulated Raman scattering in the evanescent field[END_REF]. In [START_REF] Griffiths | Raman scattering cross sections in strongly interacting liquid systems : CH3OH, C2H5OH, i-C3H7OH, (CH3)2CO, H2O, and D2O[END_REF] this cross-section is given for a pump wavelength of 488 nm. As this parameter is wavelength dependent, we have to calculate it for the wavelengths we are interested in according to the following formula [START_REF] Bischel | Wavelength dependence of Raman scattering cross sections from 200-600 nm[END_REF] :

𝜕𝜎 𝜕Ω = 𝐴𝜔 L j k𝜔 l \ -𝜔 m \ n \ (2)
A is a constant of proportionality, wp is the pump frequency, ws is the Stokes frequency and wi a UV material resonance frequency. In the liquids that we study, 𝜔 m ≪ 𝜔 l , so we can consider that fg fh is proportional to the fourth power of ws.

In alcohols there are usually two strong Raman bands, one around 2800-3000 cm -1 (denoted band I in the following) and one around 800-1040 cm -1 (denoted band II). Each band contains several Raman lines. In tables 1 and 2, we have computed the Raman gain coefficients gR in several alcohols for the strongest Raman line in band I (denoted gR,I) and for the strongest Raman line in band II (denoted gR,II) from the figures published in [START_REF] Griffiths | Raman scattering cross sections in strongly interacting liquid systems : CH3OH, C2H5OH, i-C3H7OH, (CH3)2CO, H2O, and D2O[END_REF][START_REF] Colles | Relative and absolute Raman scattering cross sections in liquids[END_REF] at the pump wavelength lp of 532 nm. We calculated gR,II for the first Stokes order emitting at l1s (gR,II at lp) and for the second Stokes order emitting at l2s (gR,II at l1s, l1s being the pumping wavelength) and finally the ratio ? @,qq rs t Bu ? @,q rs t & .

In methanol and ethanol, the ratio ? @,qq rs t Bu ? @,q rs t & is largely below 0.7, so we expect that SRS will occur on the Raman line of band I only. Propanol has two isomers, propan-1-ol and propan-2ol. In the literature we did not find values of Raman cross sections in propan-1-ol, so we are not able to calculate the Raman gain. However, in [START_REF] Krishnan | The Raman spectra of organic compounds, Part. I Methyl, Ethyl, n-Propyl and n-Butyl Alcohols[END_REF] and [START_REF] Wood | Raman spectra of a Series of Normal Alcohols and Other Compounds[END_REF] the authors both evaluated that the intensity of the Raman line at 2876 cm -1 is higher than the intensity of the Raman line at 860 cm -1 . This assumption is also stressed by a visual comparison of the Raman spectrum of propan-1-ol with the Raman spectrum of ethanol, showing the same global trends [START_REF]Spectral Database for Organic Compounds SDBS[END_REF].

This makes us believe that the strongest line in band I is the main Raman line in propan-1-ol, as for methanol and ethanol, and will be strongly favored at the expense of band II.

For propan-2-ol, more data are available in the literature and enable to calculate the Raman gains [START_REF] Griffiths | Raman scattering cross sections in strongly interacting liquid systems : CH3OH, C2H5OH, i-C3H7OH, (CH3)2CO, H2O, and D2O[END_REF][START_REF] Colles | Relative and absolute Raman scattering cross sections in liquids[END_REF]. As can be seen in Raman converter on the 2928 cm -1 line [START_REF] Lebrun | High-efficiency single-mode Raman generation in a liquid-filled photonic bandgap fiber[END_REF][START_REF] Lebrun | Optical characterizations of a Raman generator based on a hollow core photonic crystal fibre filled with a liquid[END_REF], we are confident that it will be high enough to generate SRS in this new converter. So we decided to use propan-2-ol to generate SRS on the second Stokes order of the line at 819 cm -1 . Our objective is then twofold : 1) to stop the Raman cascade at the second Stokes order of the line at 819 cm -1 in order to optimize the pump transfer energy on the target wavelength and 2) to prevent SRS on the line at 2882 cm -1 by adding extra-losses on it.

We have used a 1 m long HC-PCF with a transmission band centered around 1 µm (ref.

NKT Photonics, HC 1060-02). In order to shift the transmission band around 580 nm, the fiber has to be filled with a liquid having a refractive index around 1.39 [START_REF] Antonopoulos | Experimental demonstration of the frequency shift of bandgaps in photonic crystal fibers due to refractive index scaling[END_REF][START_REF] Lebrun | Efficient stimulated Raman scattering in hybrid liquid-silica fibers for wavelength conversion[END_REF]. As the refractive index of propan-2-ol is 1.37, one needs to mix it with a high index liquid. We have chosen to use DMSO (Dimethyl Sulfoxyde) as the diluent, as it mixes well with alcohols and as its refractive index being 1.48, allowing to keep a high proportion of the alcohol in the mixture.

The proportions of the liquids were 11% of DMSO and 89% of propan-2-ol in volume. This weak dilution ensures that the Raman gain of DMSO in the mixture remains much smaller than the Raman gain of propan-2-ol in the mixture, so that the converter will not emit on a

Raman line of DMSO. The refractive index of the mixture is firstly measured with an Abbe refractometer at 589 nm. Then a small piece of fiber (5 cm) is filled with the mixture to measure the transmission band by injecting the light of a supercontinuum laser in it [3]. At the output of this small piece of fiber, we collect and filter the light that comes from the core before sending it to the optical spectrum analyzer (OSA). Some adjustments of the liquid index are sometimes necessary and once the transmission band is satisfactory we fill the fiber of the converter with the liquid thanks to a homemade microfluidic system. Both ends are hermetically closed with micro-tanks [START_REF] Lebrun | Optical device having a liquidcore optical fibre and method for producing such a device[END_REF]. The whole device "liquid HC-PCF + tanks" is called the "Raman converter" (see insert of fig. 2). The transmission bandwidth of the converter may be slightly narrower than the one of the small piece of fiber because of the filtering of the modes along the length of propagation (typically 1 m). Another fine adjustment of the liquid index may be necessary. In fig. 1, we show three transmission bands of the same fiber, being filled with a liquid mixture having successively a refractive index of 1.384, 1.388 and 1.394. Between two fillings our microfluidic system enables to rinse the fiber to ensure the mixture of liquids is the expected one. The transmission bands were registered by optimizing the level of the output signal of the converter at 532 nm. We can estimate from fig. 1 that the right band edge of the transmission band shifts with a rate of roughly 5 nm for a change in the refractive index of 10 -3 in the range of indices we use. This estimation is helpful to optimize the position of the transmission band with the refractive index of the liquid mixture. The band is optimized to contain the pump line at 532 nm, the first Stokes order line at 556 nm and the second Stokes order line at 583 nm. The third Stokes order at 612 nm and the first Stokes order at 628 nm of the line at 2882 cm -1 have to be out of the transmission band, so that these lines will not be generated efficiently. From fig. 1, we can see that converter with n=1.394 is not optimized since the line at 583 nm suffers from too much losses. We have tested the two other converters (n=1.384 and n=1.388) and will present the best results obtained with the converter having n=1.388.

Experimental set-up

The experimental set-up is depicted in figure 2. We have used a frequency doubled microlaser emitting at 532 nm as the pump source, with a pulse duration of 560 ps (FWHM). The frequency rate repetition is 5.6 kHz and the maximum mean pump power incident on the converter is 5.75 mW (corresponding to a peak power of 1.8 kW and an energy per pulse of 1µJ). The pump beam is focused inside the Raman converter through a microscope objective (x6.3, NA 0.2). At the output of the Raman converter, we collect the light through a second microscope objective (x6.3, NA 0.2). The light is sent in an OSA through a multimode fiber.

We measured the power of the different lines emitted by the Raman converter by sending the output beam through a multimode fiber in an OSA used with a low resolution of 10 nm, much larger than the Raman linewidths. The power calibration of the system "multimode fiber + OSA" is firstly made by measuring the power in mW of the different wavelengths taken separately thanks to coloured filters, put just before the multimode fiber. Knowing the transmission of the coloured filters, we can deduce the power in mW of the different wavelengths at the output of the converter. Then for each wavelength a mean power in arbitrary unit is calculated from the OSA and can be converted into the power in mW at the output of the converter. Within this method we can have simultaneous measurements of the powers of each wavelength delivered by the Raman converter.

Results

In fig. 3 we show an example of the transmitted powers of the different lines delivered by the converter versus the incident pump power Pin at 532 nm. The injection of the pump beam in the converter was firstly optimized for low pump powers, before the threshold of the first Stokes order at 556 nm. The pump injection efficiency is then estimated to be 30%. Then we increased Pin until the appearance of the first Stokes order at 556 nm for Pin =1.5 mW. Moving slightly the injection microscope enables to optimize the output powers but sometimes led to irreversible damage of the converter (probable burning of impurities or dusts at the input of the converter) so we chose not to optimize the output power at 556 nm as we were not interested in. This explains why the power at 556 nm does not grow higher. By continuing to increase Pin we observed the appearance of the second Stokes order at 583 nm for Pin =3 mW.

The first Stokes order saturated, and the second Stokes order then grew without any saturation. Just before the maximum incident pump power we carefully optimized the output power at 583 nm and obtained a power of 1.06 mW for the maximum pump power of 5.75 mW. It corresponds to a maximum peak power of 338 W at 583 nm and an energy of 0.19 µJ per pulse. The internal conversion efficiency from the internal pump power to the second Stokes order at 583 nm is 61 % (or 67% in photon numbers). For high pump powers, we observed that the pump and the first Stokes order at 556 nm are not completely depleted which has been previously explained by the Gaussian temporal shape of the pump pulse [START_REF] Lebrun | Optical characterizations of a Raman generator based on a hollow core photonic crystal fibre filled with a liquid[END_REF][START_REF] Lebrun | Efficient stimulated Raman scattering in hybrid liquid-silica fibers for wavelength conversion[END_REF]. At maximum pump power we observed the onset of the third Stokes order at 612 nm in a core mode and the onset of the first Stokes order of the line in band I at 628 nm in a cladding mode but the powers were very weak (about a few µW). This illustrates the difficulty to adjust the transmission band that may still be improved but also shows the efficiency of the SRS in the converter.

The highest achievable resolution of the OSA we used in the present experiment was 50 pm which is in the range of the Raman spectral linewidths of liquids. So even if it was not possible to precisely measure the linewidths of the emitted lines we had an upper limit of these values. Insert in fig. 3 shows an example of the spectrum of the line at 583 nm obtained at the highest resolution of the OSA. To get a more precise value we used the fact that the spectral width of the Stokes output Δ𝜔 M/\ (FWHH) narrows with the gain according to the following formula [START_REF] Wang | Theory of stimulated Raman scattering[END_REF] :

Δ𝜔 M/\ ≅ a ln 2 2𝛾 × 𝐿 e M/\ Γ (3)
L is the length of the converter, G is the spontaneous spectral linewidth (FWHH). We took L=1 m and G=6.6 cm -1 [START_REF] Colles | Relative and absolute Raman scattering cross sections in liquids[END_REF].

To calculate we took Pp=59 W (see fig. 3). We approximate the Raman gain of propan-2-ol in the mixture by simply multiplying its value obtained in table 2 for its pure form by its proportion in the mixture, i.e. 89%. We obtained gR=2.11x10 a spectral power density of 4 mW/nm has been obtained but in the ps regime).

The output beam of the converter was sent on a diffraction grating after which we observe the spatial patterns of the modes projected on a screen (Fig. 4 (a)). On this saturated image we can see the typical shape of the modes of HC-PCF fibers with an hexagonal symmetry due to the structure of the cladding. The periodic structure of the holes in the cladding is also responsible for the presence of petal lobes around the core. Theses petal lobes can be seen as part of the fundamental mode and carry only a very small percentage of the energy. On the far-field patterns of Fig. 4(b) and (c) we can see that the Stokes lines at 556 nm and 583 nm are emitted in quasi-Gaussian modes.

From simulations of the Raman cascade [START_REF] Lebrun | Optical characterizations of a Raman generator based on a hollow core photonic crystal fibre filled with a liquid[END_REF] and from previous experimental studies [START_REF] Shan | Stimulated Raman scattering in the evanescent field of liquid immersed tapered nanofibers[END_REF] we can estimate that the temporal pulse duration at 583 nm is shorter than the pump pulse duration (i.e. 560 ps FWHM), an order of magnitude of this pulse duration at maximum pump power is expected to be around 100-200 ps FWHM by the numerical simulations.

In this work we have shown that the exploitation of the various Raman lines of common liquids and the adjustment of the refractive index of the liquid mixtures to finely control the position of the transmission bands broaden the possibilities of our Raman converters. A simple microlaser at 532 nm associated with different Raman converters will enable to access "on demand" discrete wavelengths with high conversion efficiencies in beams with high spectral brightness and good spatial quality, opening new opportunities in biological applications for instance.

Conclusion

We have designed, built and characterized a Raman converter emitting at 583 nm on the second Stokes order of a line of propan-2-ol pumped by a microlaser at 532 nm in the subnanosecond regime. The internal conversion efficiency is 67% in photon numbers and the output power is 1.06 mW, corresponding to a maximum peak power of 338 W. The beam delivered by the converter presents a Gaussian spatial structure and a spectral brightness more than 5 times higher than supercontinuum sources in this spectral range. The principle of this Raman converter can be adapted to other Raman liquids to deliver other wavelengths in the visible range. These Raman converters have applications in areas such as fluorophore excitation where pulsed low energies at specific discrete wavelengths are needed. Moreover the liquid mixtures method offers a new degree of freedom for the study of other nonlinear effects. For example it may enable to adjust the position of the zero-dispersion wavelength to realize a phase-matching condition in a four-wave mixing process for the generation of correlated photon pairs [START_REF] Barbier | Spontaneous four-wave mixing in liquid-core fibers : towards fibered Raman-free correlated photon sources[END_REF].

Figure captions

Figure 1 .

 1 Figure 1. Optimization of the transmission band of a 1 m long HCPCF 1060-02 filled with a mixture of propan-2-ol and DMSO. The refractive index of the liquid mixture is respectively 1.384 (dashed line), 1.388 (full line) and 1.394 (dotted line).

Figure 2 .

 2 Figure 2. Experimental setup. The l/4 plate enables to optimize the conversion efficiency of the IR beam to the 532 nm in the type 2 KTP crystal. L1 : focalizing lens. L2 : collimating lens. The l/2 plate associated with the polarizing beam splitter (PBS) enables to vary the incident pump power. The rejected beam is sent to a powermeter (PM reference) to serve as a pump power reference. M1, M2 : mirrors. O1, O2 : microscope objectives. MMF : multimode fiber. OSA : optical spectrum analyzer. In dashed line, removable coloured filter and powermeter (PM calibration) used for the power calibration. Insert : photo of the Raman converter.

Figure 3 .

 3 Figure 3. Mean output powers of the different lines emitted by the Raman converter at 583 nm versus the mean incident pump power. Insert : example of spectrum of the line at 583 nm at the highest resolution of the OSA.

Figure 4 .

 4 Figure 4. (a) Spatial modes at 532 nm, 556 nm and 583 nm (from left to right) at the output of the Raman converter around maximum incident pump power. (b) Example of far field pattern of the line at 556 nm (c) Example of far field pattern of the line at 583 nm (dotted line : experimental data, full line : Gaussian fit).

Table 1 . Estimated Raman gains (in x10 12 W -1 m) of different alcohols for the strongest Raman line in band I. The pump wavelength is l p=532 nm. l 1s is the wavelength of the first Stokes order.

 1 table 1 the Raman line of propan-2-ol at 2882 cm -1 has a gain of 2.43 x10 12 W -1 m. As can be seen in table2the Raman line at 819 cm -1 has a gain of 2.50 x10 12 W -1 m for the first Stokes order and a gain of 2.37 x10 12 W -1 m for the second Stokes order we want to generate.

	Liquid	Raman shift (cm -1 )		l1s (nm)		gR,I at lp
	Methanol	2837		626		2.16
	Ethanol		2928		630		2.92
	Propan-2-ol	2882		628		2.43
	Liquid	Raman shift (cm -1 ) l1s (nm) l2s (nm) gR,II at lp gR,II at l1s	g w,xx at λ Mz g w,x at λ {
	Methanol	1037	563	598	0.57	0.54	0.25
	Ethanol	881	558	587	1.04	0.99	0.34
	Propan-2-ol	819	556	583	2.50	2.37	0.97

Table 2 . Estimated Raman gains (in x10 12 W -1 m) of different alcohols for the strongest Raman line in band II. The pump wavelength is l p=532 nm. l 1s (l2s) is the wavelength of the first (second) Stokes order.

 2 is very close to unity. So in comparison with ethanol, methanol and propan-1-ol this liquid is the most favorable for SRS emission on the 819 cm -1 line. As gR,II in propan-2-ol is almost the same as gR,I in ethanol for which we already demonstrated SRS in ethanol in a

	Contrary	to	the	other	alcohols	studied,	in	propan-2-ol,
	? @,qq rs t Bu
	? @,q rs t &

  -12 W -1 m. For Aeff we took the value measured in a previous work in a Raman converter built with the same HC-PCF and filled Δ𝜔 M/\ = 1.2 𝑐𝑚 'M (corresponding to 47 pm). The estimated spectral power density at 583 nm is 22 mW/nm, which is well above the one obtained in this spectral range by sub-ns supercontinuum sources, typically 700 µW/nm in the MHz range (among the references of[3]

	with	ethanol	[11],	i.e.	12x10 -12	m².	We	found	then

Acknowledgements

The authors thank the "Fondation Coopération Scientifique Campus Paris Saclay" for supporting the work of Minh Châu Phan Huy.