Imed Kacem
email: kacem@univ-metz.fr

Vangelis Th

Paschos
email: paschos@lamsade.dauphine.fr

Weighted completion time minimization on a single-machine with a fixed non-availability interval: differential approximability

This paper is the first successful attempt on differential approximability study for a scheduling problem. Such a study considers the weighted completion time minimization on a single machine with a fixed non-availability interval. The analysis shows that the Weighted Shortest Processing Time (WSPT) rule cannot yield a differential approximation for the studied problem in the general case. Nevertheless, a slight modification of this rule provides an approximation with a differential ratio of 3- √ 5 2 ≈ 0.38.

Introduction

In this paper, we study the differential approximability of a well-known scheduling problem. Our work is motivated by the fact that the differential approximability concept has not yet been investigated for scheduling problems. Contrary to the standard approximability based on the comparison in the worst case of a heuristic solution to the optimal one, the differential approximability principle consists in comparing the heuristic solution to both the optimal and the worst solutions. More precisely, we say that heuristic A is an α-differential approximation for problem Π if for every instance I of Π the following relation holds f (A(I)) ≤ αf (opt(I)) + (1α)f (wst(I)), where f is the objective function to be minimized in problem Π and the values opt(I), A(I) and wst(I) denote respectively the values of an optimal solution, of an approximate solution and of a worst solution. This last solution is defined as the optimal solution of a problem having the same instances and set of constraints with the initial problem but the opposite goal (i.e., max, if the initial problem is a minimization one and min if the initial problem is a maximization problem. Let us also note that worst solutions are not always easy to compute. For instance, for the minimization version of travelling salesman problem, the worst solution is a Hamiltonian cycle of maximum total distance, i.e., the optimum solution of maximum travelling salesman problem. The computation of such a solution is not trivial since the latter problem is as hard as the former one. On the contrary, examples of problems for which a worst solution is easily computed are maximum independent set where the worst solution is the empty set, minimum vertex cover, where this solution is the whole vertex-set of the input graph, or, even, minimum graph-coloring, where the worst solution consists of taking a color per vertex of the input graph. The value α is called the differential ratio and it belongs to (0, 1). For more details on these approaches, the reader is invited to consult Ausiello and Paschos [START_REF] Ausiello | Reductions, completeness and the hardness of approximability[END_REF] and Demange and Paschos [START_REF] Demange | On an approximation measure founded on the links between optimization and polynomial approximation theory[END_REF].

In the studied problem we have a set of independent jobs to be performed on a single machine under the constraint of a fixed non-availability interval. The objective is to minimize the total weighted completion time under the non-resumable scenario. This problem has been proved to be NP-Hard by Adiri et al. [START_REF] Adiri | Single machine flow-time scheduling with a single breakdown[END_REF] and Lee [START_REF] Lee | Machine scheduling with an availability constraints[END_REF] and it has been studied in the literature under various criteria. Several standard approximations have been proposed. A sample of them include the worst-case analysis of heuristic methods (see for example Adiri et al. [START_REF] Adiri | Single machine flow-time scheduling with a single breakdown[END_REF]; Lee and Liman [START_REF] Lee | Single machine flow-time scheduling with scheduled maitenance[END_REF]; Sadfi et al. [START_REF] Sadfi | An improved approximation algorithm for the single machine total completion time scheduling problem with availability constraints[END_REF]; He et al. [START_REF] He | Improved algorithms for two single machine scheduling problems[END_REF]; Wang et al. [START_REF] Wang | Preemptive scheduling with availability constraints to minimize total weighted completion times[END_REF] and Breit [START_REF] Breit | Improved approximation for non-preemptive single machine flow-time scheduling with an availability constraint[END_REF]; Kacem and Chu [START_REF] Kacem | Worst-case analysis of the WSPT and MWSPT rules for single machine scheduling with one planned setup period[END_REF]; Kacem [START_REF] Kacem | Approximation algorithm for the weighted flowtime minimization on a single machine with a fixed non-availability interval[END_REF]; Kellerer and Strusevich [START_REF] Kellerer | Two simple constant ratio approximation algorithms for minimizing the total weighted completion time on a single machine with a fixed non-availability interval[END_REF]). Efficient standard approximation schemes were also published in Kellerer and Strusevich [START_REF] Kellerer | Fully Polynomial Approximation Schemes for a Symmetric Quadratic Knapsack Problem and its Scheduling Applications[END_REF]; Kacem and Mahjoub [START_REF] Kacem | Fully Polynomial Time Approximation Scheme for the Weighted Flow-time Minimization on a Single Machine with a Fixed Non-Availability Interval[END_REF] and He et al. [START_REF] He | Improved algorithms for two single machine scheduling problems[END_REF]. Other exact methods to solve this problem have been proposed in Kacem et al. [START_REF] Kacem | Single-machine scheduling with an availability constraint to minimize the weighted sum of the completion times[END_REF]- [START_REF] Kacem | Efficient branch-and-bound algorithm for minimizing the weighted sum of completion times on a single machine with one availability constraint[END_REF]. For more detail on scheduling problems under nonavailability constraints, we refer to the state-of-the-art papers by Lee [START_REF] Lee | Machine scheduling with an availability constraint[END_REF] and Schmidt [START_REF] Schmidt | Scheduling with limited machine availability[END_REF].

The review of the related literature shows that no differential approximation has been proposed to this problem according to the best of our knowledge. In a more general way, we did not find any work dedicated to the differential approximation to scheduling problems. For these reasons, this paper is a first successful attempt to develop a polynomial 3- √ 5

2 -differential approximation algorithm for the studied problem.

The paper is organized as follows. In Section 2, we present a description of the studied problem. Section 3 provides the differential analysis of the Weighted Shortest Processing Time heuristic (W SP T ). In Section 4, we show that the modification of the above heuristic yields a differential ratio of 3- √ 5

2 . Finally, Section 5 concludes the paper.

Problem Formulation

We have to schedule a set of n jobs J = {1, 2, . . . , n} on a single machine.

Every job i has a processing time p i and a weight w i . The machine is unavailable between t 1 and t 2 and it can process at most one job at a time. The fixed non-availability interval length is denoted by ∆t where ∆t = t 2t 1 .

Let C i (σ) denote the completion time of job i in a feasible schedule σ.

The aim is to find a schedule σ * that minimizes the total weighted completion time n i=1 w i C i (σ * ). With no loss of generality, we consider that all data are integers and that jobs are indexed according to the W SP T rule (i.e., p 1 w 1 ≤ p 2 w 2 ≤ . . . ≤ pn wn ). Due to the dominance of the W SP T order (see Smith [START_REF] Smith | Various optimizers for single stage production[END_REF]), an optimal schedule is composed of two sequences of jobs scheduled in nondecreasing order of their indexes (one sequence will be performed before t 1 and another after t 2 ).

If all the jobs can be inserted before t 1 , the problem studied (P) has obviously a trivial optimal solution obtained by the W SP T rule (Smith [START_REF] Smith | Various optimizers for single stage production[END_REF]). We therefore consider only the problems in which all the jobs cannot be scheduled before t 1 . The worst solution can be naturally defined as the solution W ST consisting of scheduling all the jobs after t 2 in the W SP T order. Figure 1 illustrates the two sequences W SP T and W ST and the related notations. p i for every k ∈ {1, 2, . . . , n}. Job g is identified by Z g ≤ t 1 and Z g+1 > t 2 . Variable δ denotes the idle-time between t 1 and the completion time of g in the W SP T sequence (i.e., δ = t 1 -Z g ). Moreover, ϕ * (P) denotes the minimum weighted sum of the completion times for problem P and ϕ σ (P) is the weighted sum of the completion times of schedule σ for problem P.

WSPT Analysis

In this section, we are interested in the differential approximability of the W SP T rule. We recall that the absolute worst-case performance ratio of this rule can be arbitrarily large [START_REF] Lee | Machine scheduling with an availability constraints[END_REF], but not smaller than 3 under some conditions [START_REF] Kacem | Worst-case analysis of the WSPT and MWSPT rules for single machine scheduling with one planned setup period[END_REF]. Our analysis is based upon the comparison of W SP T to W ST and it uses a lower bound introduced in Kacem and Chu [START_REF] Kacem | Worst-case analysis of the WSPT and MWSPT rules for single machine scheduling with one planned setup period[END_REF] (page 1083, Equation ( 6)).

Lemma 1 [START_REF] Kacem | Worst-case analysis of the WSPT and MWSPT rules for single machine scheduling with one planned setup period[END_REF] The following relation holds:

ϕ * (P) ≥ g i=1 w i Z i + w g+1 ∆t p g+1 (Z g + ∆t) + n i=g+1 w i (Z i + ∆t) -w g+1 ∆t p g+1 t 2
From the definition of W SP T and W ST solutions (Figure 1), the following proposition can directly be established.

Proposition 2

ϕ W ST (P) -ϕ W SP T (P) = t 2 g i=1 w i + Z g n i=g+1 w i Proposition 3 ϕ W ST (P) -ϕ * (P) ≤ t 2 g i=1 w i + t 1 n i=g+1 w i + w g+1 ∆t p g+1 δ Proof.
From Figure 1, we can establish that:

ϕ W ST (P) = n i=1 w i (t 2 + Z i ) (1)
By combining Lemma (1) and Equation (1), we obtain:

ϕ W ST (P) -ϕ * (P) ≤ n i=1 w i (t 2 + Z i ) - g i=1 w i Z i -w g+1 ∆t p g+1 (Z g + ∆t) - n i=g+1 w i (Z i + ∆t) + w g+1 ∆t p g+1 t 2 = t 2 g i=1 w i + t 1 n i=g+1 w i + w g+1 ∆t p g+1 (t 2 -(Z g + ∆t)) = t 2 g i=1 w i + t 1 n i=g+1 w i + w g+1 ∆t p g+1 δ as claimed.
The following well-known lemma will be used in what follows (see also Kacem [9]).

Lemma 4 Let a i and b i be positive numbers (i = 1, 2, . . . , k) such that b i > 0 for every i and

a 1 b 1 ≥ a 2 b 2 ≥ . . . ≥ a k b k .
The following relation holds:

k-1 i=1 a i k-1 i=1 b i ≥ a k b k (2)
Theorem 5 Let ρ be a positive number such that ρ ∈ (0, 1).

If δ ≤ ρt 1 , then W SP T is a (1 -ρ)-differential approximation for P, i.e., ϕ W SP T (P) ≤ (1 -ρ) ϕ * (P) + ρϕ W ST (P) (3)
Proof. Let X and Y be defined as follows:

X = t 2 g i=1 w i + Z g n i=g+1 w i Y = t 2 g i=1 w i + t 1 n i=g+1 w i + w g+1 ∆t p g+1 δ Hence, Y X = t 2 g i=1 w i + t 1 n i=g+1 w i + w g+1 ∆t p g+1 δ t 2 g i=1 w i + Z g n i=g+1 w i = 1 + δ n i=g+1 w i + w g+1 ∆t p g+1 δ t 2 g i=1 w i + Z g n i=g+1 w i ≤ 1 + ρt 1 n i=g+1 w i + w g+1 ∆t p g+1 δ (1 -ρ) t 1 n i=g+1 w i + t 2 g i=1 w i = 1 + ρt 1 n i=g+1 w i + w g+1 ∆t p g+1 δ (1 -ρ) t 1 n i=g+1 w i + t 2 g i=1 p i       g i=1 w i g i=1 p i      
By Lemma 4, it can be established that

w g+1 p g+1 ≤ w g p g ≤       g i=1 w i g i=1 p i      
Moreover, we know that ∆t ≤ t 2 and g i=1 p i = Z g . Hence, we deduce that:

Y X ≤ 1 + ρt 1 n i=g+1 w i + t 2 w g+1 p g+1 δ (1 -ρ) t 1 n i=g+1 w i + t 2 w g+1 p g+1 Z g ≤ 1 + ρt 1 n i=g+1 w i + t 2 w g+1 p g+1 ρt 1 (1 -ρ) t 1 n i=g+1 w i + t 2 w g+1 p g+1 (1 -ρ) t 1 ≤ 1 + ρt 1   n i=g+1 w i + t 2 w g+1 p g+1   (1 -ρ) t 1   n i=g+1 w i + t 2 w g+1 p g+1   = 1 (1 -ρ)
From the last result:

ϕ W ST (P) -ϕ W SP T (P) ϕ W ST (P) -ϕ * (P) ≥ X Y ≥ 1 -ρ ⇒ϕ W SP T (P) ≤ (1 -ρ) ϕ * (P) + ρϕ W ST (P) as claimed.
Let us note that, in the general case, W SP T can be arbitrarily bad when δ is large (compared to t 1 ). As an example, let us consider the towjob instance with p 1 = ε, w 1 = ε, p 2 = t, w 2 = tε, t 1 = t and t 2 = t + t 2 (with t >> ε). Figure 2 illustrates the WSPT, WST and OPT solutions. For this instance we have ϕ W SP T (P) ≈ t 3 , ϕ W ST (P) ≈ t 3 whereas ϕ * (P) ≈ t 2 . In other words, the differential approximation ratio in this case tends to 0.

-Differential Approximation

Based upon the previous analysis of W SP T rule, it appears that the wrong scheduling of job g + 1 can be the origin of its possible weakness. Hence, we investigate the modification of the W SP T sequence based upon the following algorithm H, which tests the two possibilities of scheduling job g + 1 before and after the non-availability interval. This algorithm generates two sequences. The first one is the W SP T sequence. In the second one (denoted as W SP T 2), job g + 1 is scheduled before t 1 and the other jobs are scheduled in the W SP T order. The output of this algorithm is the best generated sequence.

HEURISTIC H

(i) Construct the sequence W SP T and compute ϕ W SP T (P).

(ii) Let y be the y th job in J -{g + 1} according to the WSPT order (y < g) such that y i=1 p i + p g+1 ≤ t 1 and y+1 i=1 p i + p g+1 > t 1 . Construct the sequence: W SP T 2 = 1, 2, . . . , y, g + 1, y + 1, y + 2, . . . , g, g + 2, g + 3, . . . , n and, if feasible, compute ϕ W SP T 2 (P).

(iii) Output the best among the solutions obtained in steps (i) and (ii) of value ϕ H (P) = min {ϕ W SP T (P) , ϕ W SP T 2 (P)}.

It can be easily seen that Heuristic H can be implemented in O (n log(n)) time.

To illustrate Heuristic H, let us consider the following four-job instance: p 1 = 1; w 1 = 2; p 2 = 2; w 2 = 3; p 3 = 3; w 3 = 4; p 4 = 1; w 4 = 1; t 1 = 4; t 2 = 7. Given this instance, we have: ϕ W SP T (P) = 62; ϕ W SP T 2 (P) = 55 and ϕ H (P) = 55. Figure 3 illustrates the obtained schedules. In this case, we have g = 2 and y = 1.

Proof. From Figure 4, the value ϕ W ST (P)ϕ W SP T 2 (P) can be computed:

ϕ W ST (P) -ϕ W SP T 2 (P) =
(P) ≤ ρ (1 + ρ) ϕ * (P) + 1 (1 + ρ) ϕ W ST (P) (5)
Proof. By definition:

ϕ H (P) = min {ϕ W SP T (P), ϕ W SP T 2 (P)} ≤ ρ (1 + ρ) ϕ W SP T (P) + 1 (1 + ρ) ϕ W SP T 2 (P)
Hence,

ϕ W ST (P) -ϕ H (P) ≥ ρ (1 + ρ) (ϕ W ST (P) -ϕ W SP T (P)) + 1 (1 + ρ) (ϕ W ST (P) -ϕ W SP T 2 (P))
From Propositions 2 and 6, we can deduce the following inequality:

ϕ W ST (P) -ϕ H (P) ≥ ρ (1 + ρ)   t 2 g i=1 w i + Z g n i=g+1 w i   + 1 (1 + ρ)   w g+1 δ p g+1 ∆t + δ n i=g+1 w i  
By assumption, δ > ρt 1 . Therefore, we deduce that:

ϕ W ST (P) -ϕ H (P) > ρ (1 + ρ)   t 2 g i=1 w i + Z g n i=g+1 w i   + 1 (1 + ρ)   w g+1 δ p g+1 ∆t + ρt 1 n i=g+1 w i   > ρ (1 + ρ) t 2 g i=1 w i + ρ (1 + ρ)   w g+1 δ p g+1 ∆t + t 1 n i=g+1 w i   > ρ (1 + ρ)   t 2 g i=1 w i + t 1 n i=g+1 w i + w g+1 δ p g+1 ∆t  
Finally, from Proposition 3 we obtain:

ϕ W ST (P) -ϕ H (P) ≥ ρ (1 + ρ)
(ϕ W ST (P)ϕ * (P)) , and then, Equation ( 5) is proved.

Theorem 8 Heuristic H is a 3- √ 5

2 -differential approximation for problem P.

Proof. Let ρ be a positive number such that ρ ∈ (0, 1). By combining Theorems 5 and 7, Heuristic H is a ρ (1+ρ) -differential approximation for P (if δ > ρt 1 ) and a (1ρ)-differential approximation for P (if δ ≤ ρt 1 ). Hence, by taking ρ ∈ (0, 1) such that ρ (1+ρ) = 1ρ, we obtain ρ = that completes the proof.

Conclusion

Motivated by the absence in the literature of differential approximability analysis for scheduling problems, this paper aims to investigate this new direction. The considered study is related to the weighted completion time minimization on a single machine with a fixed non-availability interval. The analysis of the Weighted Shortest Processing Time (W SP T ) rule shows that this rule cannot yield a differential approximation for the studied problem in the general case. Nevertheless, a slight modification of this rule provides a 3- √ 5

Figure 1 :

 1 Figure 1: Illustration of W SP T and W ST solutions

Figure 2 :

 2 Figure 2: Worst-case example for W SP T

Figure 3 : 6 ϕ

 36 Figure 3: Illustration of H

-differential approximation. Ongoing research will aim at designing more efficient differential approximations for scheduling problems (in particular, differential approximation schemes).

* Research partially supported by the French Agency for Research under the DEFIS program TODO, ANR-09-EMER-010


F Doisneau

J Dupays

A Murrone

F Laurent

M Massot

Eulerian VS Lagrangian simulation of unsteady two-way coupled coalescing two-phase flows in solid propellant combustion

Keywords: Solid Rocket Motor, Eulerian Multi-Fluid, Coalescence, Moderately dense spray Moteur à Propergol Solide, Multi-Fluide Eulérien, Coalescence, Spray modérément dense

Solid rocket motor flows depend strongly on alumina granulometry, which spreads up with coalescence. Solving droplet sizes in unsteady flows is a modeling and numerical challenge, specially when high loadings generate retrocoupling. As an alternative to Lagrangian approaches, the Eulerian Multi-Fluid model describes polydispersity with size-sorted droplet fluids called sections. A two size moment method, that spares sections, is implemented in the CEDRE code. A new transport and coupling method allows to capture moderately dense spray structures with low diffusion. The aim is threefold : compare CEDRE's Lagrangian and new Eulerian two-phase approaches, prove the feasibility of unsteady and complex motor flow simulations, and show the effect of coalescence on instabilities.

Résumé

Simulation eulérienne et lagrangienne d'écoulements diphasiques instationnaires avec coalescence et retrocouplage en propulsion solide. L'écoulement d'un moteur à propergol solide dépend fortement de la granulométrie d'alumine, qui s'étale avec la coalescence. Résoudre la taille des gouttes dans les écoulements instationnaires est un défi numérique et de modélisation, surtout quand les forts chargements génèrent du retrocouplage. Une alternative à l'approche lagrangienne est le modèle Multi-Fluide eulérien qui décrit la polydispersion en résolvant des fluides de gouttes triés par tailles, appelés sections. Une méthode à deux moments en taille, qui économise des sections, est implémentée dans le code CEDRE. Une nouvelle méthode pour le transport et le couplage permet de capturer les structures des sprays modérément denses avec une diffusion faible. L'objectif est triple : comparer les approches diphasiques lagrangienne et eulérienne améliorée de CEDRE, prouver la faisabilité de simulations d'écoulements complexes dans les moteurs et exhiber l'effet de la coalescence sur les instabilités.

Introduction

In aluminized propellant combustion the droplet size distribution, which depends on propellant formulation and combustion and on its own distributed combustion, is strongly polydisperse [START_REF] Dupays | Two-phase unsteady flow in solid rocket motors[END_REF]. Moreover, it can dramatically change because of coalescence, that is the amalgamation of two droplets after a collision. But coalescence intensity depends on droplet velocities, which are conditioned by size. So droplet size and velocity are intimately intricate, which will be referred to as the "size/velocity coupling". The crucial importance of polydispersity and size/velocity coupling has been demonstrated in general cases in [START_REF] Kah | Eulerian quadrature-based moment models for dilute polydisperse evaporating sprays[END_REF] and in solid rocket motor (SRM) cases in [START_REF] Doisneau | Eulerian Multi-Fluid models for the simulation of dynamics and coalescence of particles in solid propellant combustion[END_REF][START_REF] Doisneau | Two-way coupled simulation of acoustic waves in polydispersed coalescing two-phase flows : application to Solid Rocket Motor instabilities[END_REF]. Accurate treatment of polydispersity, and of what controls it, is required since the size parameter directly drives the droplet-gas interactions and therefore instabilities. The Williams-Boltzmann equation allows to describe polydisperse, dilute and moderately dense liquid sprays by solving the evolution of the dispersed phase number density function (NDF) due to advection, drag force and heat transfer that couple the gaseous phase, and droplet-droplet interactions such as coalescence [START_REF] Hylkema | A random particle method to simulate coalescence phenomena in dense liquid sprays[END_REF][START_REF] Laurent | Eulerian Multi-Fluid modeling for the numerical simulation of coalescence in polydisperse dense liquid sprays[END_REF][START_REF] Doisneau | Eulerian Multi-Fluid models for the simulation of dynamics and coalescence of particles in solid propellant combustion[END_REF]. We focus here on moderately dense sprays, which have a significant enough mass fraction to strongly perturbate the gas but a small enough volume fraction to avoid particle volume effects. Moreover, coalescence occurs in many moderately dense sprays and is an issue in SRM.

A first strategy is to approximate the NDF by discrete numerical parcels of particles of various sizes through a Lagrangian Monte-Carlo approach [START_REF] Bird | Molecular gas dynamics and the direct simulation of gas flows[END_REF], which is generally considered to be the most accurate for solving the Williams-Boltzmann equation. But it introduces diffusion when coupled to a Eulerian fluid and is difficult to parallelize. Collisions and coalescence can be treated in this context, usually with a stochastic approach [START_REF] Hylkema | A random particle method to simulate coalescence phenomena in dense liquid sprays[END_REF][START_REF] Laurent | Eulerian Multi-Fluid modeling for the numerical simulation of coalescence in polydisperse dense liquid sprays[END_REF] but here again, accuracy and performance are difficult to reconcile. As an alternative the Eulerian Multi-Fluid model discretizes the NDF using a finite volume approach of the size variable but assumes size/velocity correlations to reduce the phase space dimension [START_REF] Laurent | Multi-fluid modeling of laminar poly-dispersed spray flames: origin, assumptions and comparison of the sectional and sampling methods[END_REF]. After integration on size intervals, the so-called "sections", the resulting conservation equations are similar to those of the pressureless gas dynamics (PGD), with source terms accounting for gas-particle and particle-particle interactions. The resolution of spray dynamics and coalescence coupled to the gaseous flow field can be performed with one size moment [START_REF] Laurent | Eulerian Multi-Fluid modeling for the numerical simulation of coalescence in polydisperse dense liquid sprays[END_REF] or two size moments to approximate the size distribution in each section. The latter, introduced in [START_REF] Dufour | A second-order Multi-Fluid model for evaporating sprays[END_REF] in order to treat evaporation, is well-suited in SRM cases as demonstrated in [START_REF] Doisneau | Eulerian Multi-Fluid models for the simulation of dynamics and coalescence of particles in solid propellant combustion[END_REF]. It is indeed a good compromise between high-order moment methods yielding accuracy but algebraic complexity -for instance on size-velocity coupling [START_REF] Vié | Size-velocity correlations in high order moment methods for polydisperse evaporating sprays : modelling and numerical issues[END_REF]-and modeling flexibility allowing to solve complex configurations. To solve the complex interactions of acoustic waves in a polydisperse spray without being limited by the smallest droplet dynamic time scales, a new splitting algorithm for which accuracy and stability were proven, was implemented for the Eulerian approach in the CEDRE code [START_REF] Doisneau | Numerical strategy for two-way coupling in unsteady polydisperse moderately dense sprays[END_REF][START_REF] Doisneau | Numerical strategy for two-way coupling in unsteady polydisperse moderately dense sprays[END_REF], a multi-physics platform on general unstructured meshes, for both research and industrial applications, in the fields of energetics and propulsion [START_REF] Refloch | CEDRE software[END_REF]. A low diffusion transport scheme for the dispersed phase has recently been implemented in CEDRE, that allows to capture moderately dense spray structures. But the impact of the entire strategy in SRM has never been assessed.

The long term purpose of our research is to reach a predictive Euler-Euler solver for the large eddy simulation (LES) of SRM in order to resolve the complex dynamics of polydisperse sprays, coalescence, combustion, and acoustic waves interaction, which are at the origin of the pressure oscillations. The objective of the present paper is to qualify the methods implemented in CEDRE in a realistic configuration of a P230 SRM featuring instabilities. We will focus on the interaction of hydrodynamic gaseous flow instabilities coupled to a polydisperse spray segregation dynamics and coalescence, where the resolution of the spray size spectrum is crucial. Such a study is conducted on an under-resolved mesh, where the main flow structures and the main numerical peculiarities are still featured, with a special attention devoted to the fact that the conclusions to be drawn for this study are valid for a larger simulation, be it in the framework of DNS or LES. Besides, since a Lagrangian as well as the Two Size Moment Eulerian Multi-Fluid method have been implemented in the CEDRE code, we perform a comparison of the two approaches regarding cost, ease of use and accuracy. Let us insist on the fact that the proposed Euler-Euler formalism and strategy provide for the first time a reliable resolution of all the physical phenomena involved, with algorithms that have the potential to be extended to realistic cases.

The paper is organized as follows : section 1 describes Eulerian and Lagrangian modeling of a coalescing spray from their common basis at the kinetic level. Section 2 describes recent improvements in Eulerian simulations to capture accurately the two-way coupling and the spray structures despite the high inhomogeneities that are vacuum and dense zones. Section 3 provides quantitative accuracy and cost comparisons between the Euler-Euler and Euler-Lagrange frameworks on a SRM without and with coalescence. Finally section 4 compares the advantages and drawbacks of both approaches for SRM simulation.

Eulerian and Lagrangian modeling for coalescing sprays

A common kinetic basis : the Williams equation

Let us define the number density function (NDF) f of the spray, where f (t, x, u, T, S)dxdSdT du denotes the average number of droplets (in a statistical sense) at time t, in a volume of size dx around a space location x, with a velocity u, internal temperature T , and size S in a dudT dS-neighborhood of (u, T, S). The evolution of the NDF is described by the classical Williams-Boltzmann equation. Considering local momentum and heat transfer with the gas due to drag forces and conducto-convective fluxes and considering collisions between droplets but no evaporation nor break-up, it reads :

∂ t f + u • ∂ x f + ∂ u • (Ff ) + ∂ T H c p (T ) f = Γ ( 1
)
where F is the drag force per unit mass, H is the heat transfer per unit mass and Γ is the collision source term. These terms require models which should take into account physical phenomena at the droplet scale, which specific heat capacity is c p (T ). They are typically modeled by Stokes' laws with eventually Schiller-Nauman or Ranz-Marshall corrections, yielding characteristic dynamic and thermal times τ u (S) and τ T (S) strongly depending on droplet size. Regarding collisions, we assume [START_REF] Laurent | Eulerian Multi-Fluid modeling for the numerical simulation of coalescence in polydisperse dense liquid sprays[END_REF][START_REF] Doisneau | Eulerian Multi-Fluid models for the simulation of dynamics and coalescence of particles in solid propellant combustion[END_REF] :

[C1] We only take binary collisions into account.

[C2] The mean collision time is very small compared to the intercollision time.

[C3] Mass, momentum and droplet enthalpy are preserved during collisions.

[C4] Every collision leads to coalescence. The collision operator resulting from these hypotheses is developed as Γ = Q + -Q -, respectively corresponding to creation and destruction of droplets due to coalescence. Considering two precursor droplets of volumes v and v colliding to form a new droplet of volume v, the operators read [START_REF] Hylkema | A random particle method to simulate coalescence phenomena in dense liquid sprays[END_REF]:

Q + = 1 2 T u ,v ∈[0,v] f (t, x, u , T , v )f (t, x, u , T , v )B(|u -u |, v , v )Jdv du dT Q -= T u ,v f (t, x, u, T, v)f (t, x, u , T , v )B(|u -u |, v, v )dv du dT (2)
where

v = v -v , v u = vu -v u , v h(T ) = vh(T ) -v h(T ) are pre-collisional parameters, J = (v/v ) n d +1 c p /c p
is the Jacobian of the mapping (v, T, u) → (v , T , u ) with n d the dimension of the velocity phase space, c p = c p (T ) and c p = c p (T ) and B is the collision/coalescence probability kernel :

B(|u -u |, v, v ) =|u -u | β(v, v )E(|u -u |, v, v ).
In this kernel, β(v, v ) = π (r + r ) 2 is the impact parameter or geometric cross-section and E is the collision efficiency. This efficiency accounts for the effects of the surrounding gas flow. Experimental results in SRM conditions suggest to use such models [START_REF] D'herbigny | Etude expérimentale et numérique pour la validation d'un modèle de coalescence[END_REF] to avoid significant over-estimations on coalescence intensity. In the following, we use a Langmuir-Blodgett law to get a reasonable coalescence intensity and to prove that these models are compatible with Eulerian approaches. The coalescence terms (2) are referred to as ballistic since their collision probability kernel is based on the resolved particle velocities. They are typical of DNS. In a LES context, sub-grid scale models (SGS) models should include coalescence terms linked to SGS velocity differences, often described as a statistical dispersion.

Lagrangian approach

A Lagrangian or stochastic particle method is implemented in the CEDRE code and consists in approximating the NDF with a finite number I of parcels i representing several droplets with the same position x pi , size S pi , velocity u pi and temperature T pi , as follows :

f (t, x, u, T, S) = i∈[1,I] w pi δ(x -x pi (t))δ(S -S pi (t))δ(u -u pi (t))δ(T -T pi (t)) (3)
where w pi is a weight. Practically in CEDRE, Eq. ( 1) is not solved but splitted in three steps [START_REF] Murrone | Numerical modeling of dispersed two-phase flows[END_REF] that are transport in the physical and phase space, collision and fragmentation, which is not considered here.

In the transport step, the r.h.s. of (1) decouples for each parcel i in the following system of ODE :

∂ t x pi = u pi ; ∂ t S pi = 0; ∂ t u pi = F pi ; ∂ t T pi = H pi . (4)
which has to be solved with a method that respects the coupling to the gaseous phase. An unconditionally stable first order scheme is chosen to do so. In the collision step, pairs of particles are randomly chosen among each cell and their state is modified according to an algorithm described in [START_REF] Hylkema | A random particle method to simulate coalescence phenomena in dense liquid sprays[END_REF][START_REF] Laurent | Eulerian Multi-Fluid modeling for the numerical simulation of coalescence in polydisperse dense liquid sprays[END_REF].

On the one hand, particle trajectories can be accurately integrated because they depend only on time so numerical diffusion is very limited. Moreover the method allows any type of particle trajectory crossings (PTC), especially those occurring on symmetry axes, close to a wall or at the rim of vortices. On the other hand, more particles are needed when degrees of freedom such as polydispersity or unsteadiness are solved which is costly in terms of CPU and memory. Moreover numerical diffusion is still introduced by the gas two-way coupling and the collision algorithm since source terms have to be evaluated on the Eulerian mesh. Finally, in a high performance computing context, possibilities of vectorization/parallelization or implicitation are limited. The particle load balance is particularly difficult to achieve [START_REF] Garcia | Développement et validation du formalisme Euler-Lagrange dans un solveur parallèle non-structuré pour la simulation aux grandes échelles[END_REF].

Two size Moment Eulerian Multi-Fluid approach

To solve Eq. ( 1) with an affordable Eulerian approach, we need to reduce the size of the phase space so we only consider the moments of f of order zero and one in the velocity and temperature variables. To close the system the following assumptions are introduced for a given droplet size S, at (t, x) :

[H1] There is only one velocity ū(t, x, S) and no dispersion.

[H2] There is only one temperature T (t, x, S) and no dispersion. Hypothesis [H1] is known as the monokinetic hypothesis [START_REF] Laurent | Multi-fluid modeling of laminar poly-dispersed spray flames: origin, assumptions and comparison of the sectional and sampling methods[END_REF]. It is equivalent to reducing the velocity distribution support to a one dimensional sub-manifold parameterized by droplet size. Hypothesis [H2] is similarly introduced so that temperatures are locally correlated for droplets of a given size, which is, as soon as the monokinetic hypothesis holds, correct when τ T (S) is close or smaller than τ u (S). We then choose a discretization 0 = S 0 < S 1 < • • • < S Nsec = ∞ for the droplet size phase space and we average the conservation laws for number, velocity and enthalpy over each fixed size interval [S k-1 , S k [, called section. The set of droplets in one section can be seen as a "fluid" for which conservation equations are written, the sections exchanging mass, momentum and enthalpy. We also assume :

[H3] In each section, the velocity does not depend on the size of the droplets.

[H4] In each section, the temperature does not depend on the size of the droplets.

[H5] In each section, the form of n as a function of S is presumed. We choose for assumption [H3] the notation ū(t, x, S) = ūk (t, x) to designate the constant velocity distribution in section k. Hypothesis [H4] is similarly introduced here and is equivalent to presuming h(t, x, S) = hk (t, x) as an enthalpy distribution in section k, and the corresponding effective temperature Tk (t, x) which allows to define the sectional specific heat capacity ck = c p ( Tk ). Finally, all these hypotheses are equivalent to presume the following form of NDF :

f (t, x, u, T, S) = Nsec k=1 κ k (t, x, S)δ(u -ūk (t, x))δ(T -Tk (t, x)).
The Two Size Moment Multi-Fluid model, described in [START_REF] Dufour | A second-order Multi-Fluid model for evaporating sprays[END_REF] and extended to coalescence in [START_REF] Doisneau | Eulerian Multi-Fluid models for the simulation of dynamics and coalescence of particles in solid propellant combustion[END_REF], is based on a two-coefficient exponential reconstruction :

κ k (t, x, S) = a k (t, x) exp(-b k (t, x)S) for S ∈ [S k-1 , S k [, where (a k (t, x), b k (t, x)) k are chosen to preserve the following set of moments (n k (t, x), m k (t, x)) k .
Finally the two-phase flow is described by a system including gaseous Navier-Stokes equations coupled to the conservation equations for the chosen moments. It is written here with Euler equations for simplicity :

                                         ∂ t ρ g + ∂ x •(ρ g u g ) = 0 ∂ t (ρ g u g ) + ∂ x •(ρ g u g ⊗u g + pI) = - Nsec k=1 m k Fk ∂ t (ρ g e g ) + ∂ x •(ρ g e g u g ) + p∂ x u g = - Nsec k=1 m k Hk - Nsec k m k Fk (u g -ūk ) ∂ t n k + ∂ x •(n k ūk ) = C n+ k -C n- k ∂ t m k + ∂ x •(m k ūk ) = C m+ k -C m- k ∂ t (m k ūk ) + ∂ x •(m k ūk ⊗ ūk ) = m k Fk + C mu+ k -C mu- k ∂ t (m k hk ) + ∂ x •(m k hk ūk ) = m k Hk + C mh+ k -C mh- k                k = 1, N sec (5)
where Fk and Hk are averaged on the section from their kinetic expressions. The coalescence source terms come from the expressions (2) that have been integrated on velocity thanks to the monokinetic hypotheses and on the two colliding partner section size intervals thanks to the Multi-Fluid reconstruction hypotheses. They couple the sections together but they are local in the physical space, analogous to chemical terms for instance. The size-integration method is crucial and has been carefully developed in the CEDRE code regarding computational performance : it is based on an adaptive quadrature technique that can handle smooth as well as steep and even singular size distributions [START_REF] Doisneau | Eulerian Multi-Fluid models for the simulation of dynamics and coalescence of particles in solid propellant combustion[END_REF].

Numerical methods for Euler-Euler simulations

Achieving two-way coupling in moderately dense sprays

A numerical method based on operator splitting is suggested in [START_REF] Doisneau | Numerical strategy for two-way coupling in unsteady polydisperse moderately dense sprays[END_REF][START_REF] Doisneau | Numerical strategy for two-way coupling in unsteady polydisperse moderately dense sprays[END_REF] to ensure a robust and accurate two-way coupling for moderately dense sprays. It is referred to as Acoustic-Convective Splitting (ACS) since it has been designed to treat efficiently the fast scales of acoustics in a two-phase context. Its benefits on cost and accuracy have been quantified on two-phase acoustics in [START_REF] Doisneau | Numerical strategy for two-way coupling in unsteady polydisperse moderately dense sprays[END_REF][START_REF] Doisneau | Numerical strategy for two-way coupling in unsteady polydisperse moderately dense sprays[END_REF]. In complex flows however, twoway coupling unsteadiness has a similar sensitivity to polydispersity so the method is a priori well adapted for SRM cases. The method is briefly described here, as well as the treatment of coalescence in this context. The ACS is a two stage splitting based on the following timescale hierarchy : liquid convection generally yields the slowest scales while the fastest scales are gas transport scales, which include acoustic waves. The relaxation scales, that are drag and heat transfers between the liquid and the gas, are to be linked to the acoustic scales in order to account for two-phase acoustics. System ( 5) is considered as composed of transport operators T g and (T k ) k for the fluids, a relaxation operator R and a coalescence operator C. These operators apply on

U(t) = ρ g , ρ g u g , ρ g h g , (n k , m k , m k ūk , m k hk ) k=1,Nsec t (t) ∈ R 4Nsec+3 .
The transport operators read :

T g          ∂ t ρ g + ∂ x •(ρ g u g ) = 0 ∂ t (ρ g u g ) + ∂ x •(ρ g u g ⊗u g + pI) = 0 ∂ t (ρ g e g ) + ∂ x •(ρ g e g u g ) + p∂ x u g = 0 T k                ∂ t n k + ∂ x •(n k ūk ) = 0 ∂ t m k + ∂ x •(m k ūk ) = 0 ∂ t (m k ūk ) + ∂ x •(m k ūk ⊗ ūk ) = 0 ∂ t (m k hk ) + ∂ x •(m k hk ūk ) = 0
The local operators that are relaxation and coalescence ensure the coupling between all the "fluids" and all the sections respectively :

R                                          ∂ t ρ g = 0 ∂ t (ρ g u g ) = - Nsec k=1 m k Fk ∂ t (ρ g e g ) = - Nsec k=1 m k Hk - Nsec k m k Fk (u g -ūk ) ∂ t n k = 0 ∂ t m k = 0 ∂ t (m k ūk ) = m k Fk ∂ t (m k hk ) = m k Hk                k = 1, N sec C                ∂ t n k = C n+ k -C n- k ∂ t m k = C m+ k -C m- k ∂ t (m k ūk ) = C mu+ k -C mu- k ∂ t (m k hk ) = C mh+ k -C mh- k k = 1, N sec
We can rewrite system (5) as an operator equation

U exact (t + ∆t c ) = C + R + T g + Nsec k=1
T k (∆t c ) U(t) between t and t+∆t c . We now introduce the structure of the first stage of ACS, performed here according to a Strang splitting :

U Strang (t + ∆t c ) = C + R + T g ∆t c 2 Nsec k=1 T k (∆t c ) C + R + T g ∆t c 2 U(t).
The convective splitting time step is chosen to allow liquid convection to occur at its physical timescale which matches numerical criterion of the transport scheme ∆t c ∼ τ c k . We then suggest to consider the relaxation and the gas transport as a two-phase acoustics operator A. Compared to the strategy presented in [START_REF] Doisneau | Numerical strategy for two-way coupling in unsteady polydisperse moderately dense sprays[END_REF][START_REF] Doisneau | Numerical strategy for two-way coupling in unsteady polydisperse moderately dense sprays[END_REF], the coalescence operator has been included at the acoustic level since it can be stiff in dense zones and it modifies the droplet sizes, which condition the other time steps. Since C has a high numerical cost, it has been included through a first order Lie splitting to spare resources, while the acoustic operator is built on a Strang splitting. On a time step ∆t a respectful of acoustics, we finally have :

C + R + T g (∆t a ) ≈ C(∆t a ) A(∆t a ) with A(∆t a ) = R ∆t a 2 T g (∆t a ) R ∆t a 2 .

Capturing spray structures

Since the spray follows the PGD in DNS or has a weak pressure in LES, strong gradients can appear, depending on the flow topology. For instance vortices yield, as soon as droplets are inertial enough, a remarkable segregation corresponding to strong density gradients. It is crucial to preserve these spray structures to remain accurate and to achieve predictive simulations. Therefore, a high quality solution is required, either with refined meshes or with an accurate method [START_REF] Fréret | On the role of preferential segregation in flame dynamics in polydisperse evaporating sprays[END_REF][START_REF] Martinez | Simulation aux grandes échelles de l'injection de carburant liquide dans les moteurs à combustion interne[END_REF]. To spare degrees of freedom and cost, a high order transport scheme, yielding low diffusion, can be used. This is rarely achieved in Euler-Euler simulations : for instance most high order numerical methods developed for Euler or Navier-Stokes equations are ill-suited because of the peculiarities of PGD that are vacuum zones and δ-shocks [START_REF] De Chaisemartin | Polydisperse evaporating spray turbulent dispersion : Eulerian model and numerical simulation[END_REF]. A specific method has been developed for spray transport on generic unstructured meshes in CEDRE. It is based on a second order MUSCL reconstruction with unstructured gradient limitation [START_REF] Courbet | Space discretization methods[END_REF] and with a Godunov-type flux scheme, but adapted to PGD. A first order method has been used for most simulations in this paper but the second order method is used and specifically discussed in section 4.2.

Towards Large Eddy Simulation

It is now acknowledged that, in many problems, LES provides more predictive solutions compared to RANS, whereas DNS is most of the time out of reach. For SRM studies, LES requires refined 3D meshes and a SGS modeling [START_REF] Simoes | Modélisation eulérienne de la phase dispersée dans les moteurs à propergol solide, avec prise en compte de la pression particulaire[END_REF][START_REF] Moreau | Development of gas-particle Euler-Euler LES approach: A priori analysis of particle sub-grid models in Homogeneous Isotropic Turbulence[END_REF] that is specific to two-phase flows. Some models are available in the community, even if they cover only a partial range of Stokes numbers [START_REF] Masi | Étude théorique et numérique de la modélisation instationnaire des écoulements turbulents anisothermes gaz-particules par une approche Euler-Euler[END_REF]. However, we claim that their implementation is straightforward once the PGD equations can be resolved with enough numerical efficiency, accuracy and robustness, which is the purpose of the present contribution. Remarkably a hybrid strategy between PGD and velocity-dispersed zones has already been suggested [START_REF] Boileau | Robust numerical schemes for Eulerian spray DNS and LES in two-phase turbulent flows[END_REF]. In this paper, we do not achieve predictive simulations regarding models and meshes but we perform DNS-like simulations : this to validate methods that are a prerequisite to any LES. The question of a robust and high order transport scheme is crucial for resolved turbulent scales to be effectively captured. The coupling method as introduced above aims at solving robustly the transfers between the two phases' resolved scales and is designed to be directly applicable in the framework of LES. And ballistic coalescence -based on the resolved particle velocity differences as above-remains significant even when SGS collisions appear.

Numerical simulation of instabilities in the P230 SRM

In this section, we give numerical results from the CEDRE code on a 2D configuration which reproduces the P230 motor that equips Ariane 5. This motor features obstacle vortex shedding (VSO) due to the thermal protection and parietal vortex shedding (VSP), which occurs in long segmented motors [START_REF] Lupoglazoff | Parietal vortex shedding as a cause of intstability for long solid propellant motors. numerical simulations and comparisons with firing tests[END_REF]. The chamber can host acoustic modes and flows through a supersonic nozzle. The propellant is aluminized so the flow contains liquid droplets of aluminum oxide from a zone that is close to the wall down to the nozzle. The droplet size distribution is often considered as bimodal after combustion [START_REF] Price | Combustion of Metallized Propellants[END_REF] so we choose to model the spray injection with non-burning aluminum oxide droplets of size d 1 and d 2 .

Numerical configuration

We consider a 2D axisymmetric mesh of the P230 at a specific moment of the grain combustion which is shown in Fig. 1. It is an unstructured 45, 500 generic cell mesh and the computation is paralellized on 32 domains. We use a coarse mesh, compared to typical SRM meshes with several million points, only to assess the feasibility of complex simulations with an Euler-Euler approach. For this, we now consider either the possibility of coalescence or not. The simulation starts from a steady state of the single phase gaseous flow field. When the particles are injected, the pressure builds up to a new value of average pressure. The simulated time is long enough to perform a harmonic study of the instabilities. Gas and particles are injected at a fixed surface flow rate from the two walls representing the two propellant segments. The injected liquid flow rate is fairly the same for the two segments. The loading is high since the liquid mass fraction reads 30% at the injection and can reach much more in the accumulation zones so a two-way coupling approach is required. The small droplets (diameter d 1 ) represent 90% of the injected mass and the big droplets (diameter d 2 ) the remaining 10%. The number concentration is around 10 15 at the wall so coalescence is to be considered. Because of the high density of alumina, the liquid volume fraction is much smaller than one so the liquid volume occupation effect can be neglected, the spray is moderately dense and system ( 5) is relevant. ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? The case is numerically costly because of the extreme stiffness yielded by the various space and time scales. Indeed the motor is more than 25m long while the smallest hydrodynamic scales are a few centimeters. Its main instability mode is low frequency around 20Hz while the small particle response time is close to the microsecond. The injection and nozzle output velocities, respectively a few m/s and Mach 3, introduce time scales with a ratio of about 3, 000. The case is also costly because of the polydisperse coalescing liquid phase. First, it features a strongly bimodal disperse phase : because droplet characteristic times depend on the square of their diameters, this introduces two very different time scales, with a ratio of 10 4 . Second, the flow is injected at a very low velocity, encounters hydrodynamic and acoustic instabilities and finally crosses a supersonic nozzle with a high load ratio. We recapitulate all the transport characteristic times in Table 3.1. Since coalescence is a quadratic operator that is intrinsically costly, treating the transport stiffness by solving the smallest scales is definitely out of reach while reducing the cost with an implicit integration for instance would deteriorate accuracy on the unsteady phenomena. Thanks to the ACS splitting approach, the main time step can be chosen above the small particle response time : this allows to spare an order of magnitude on computational cost without compromising accuracy on the physics of small droplets while most classical coupling methods require a time step below the smallest particle scale for stability reasons [START_REF] Doisneau | Numerical strategy for two-way coupling in unsteady polydisperse moderately dense sprays[END_REF][START_REF] Doisneau | Numerical strategy for two-way coupling in unsteady polydisperse moderately dense sprays[END_REF]. The splitting strategy moreover decouples the coalescence intrinsic stiffness from the other time scales, allowing to treat it by a separate integration method if needed. We now detail the numerical configuration for the size variable in each framework. In the Eulerian Two Size Moment method, the section discretization conditions the simulation's accuracy. For the inert case two sections are required to solve perfectly the bimodal distribution, no matter what bounds they have, thanks to the Two Size Moment method which can accurately approach the dynamics of any Dirac size distribution. And the two kinds of droplets will not change size during their transport. For the coalescing simulation, the monokinetic hypotheses prevents coalescence from occuring among droplets of the same section. This is relevant if the flow structure prevents PTC for droplets with fairly the same sizes. But the coalescence phenomenon does happen between different sections, which particles have different dynamics.

Since the injected droplets are strongly bimodal and the flow yields particle acceleration, coalescence will occur between the two sections. The phase space is then discretized as follows. The first section is defined like in the inert case since no growth can occur. The second section now ranges from d 2 to d 3 = 1.3d 2 to account for the big droplet size increase. A third section is also created to host even bigger droplets that can be created above d 3 . A discretization as coarse as 3 sections has proven to be sufficient to capture basic coalescence size/dynamics coupling with the Two Size Moment Multi-Fluid method [START_REF] Doisneau | Eulerian Multi-Fluid models for the simulation of dynamics and coalescence of particles in solid propellant combustion[END_REF]. In the Lagrangian computations, the parametrization is the same as in [START_REF] Estivalezes | Simulations lagrangiennes polydisperses avec collisions et coalescence sur une configuration Ariane V 2D à l'instant 107s[END_REF]. Numerical particles of the two initial sizes are injected at each wall cell with a period of 4.10 -4 s. The number of particles per computational domain is however limited to 300, 000. Since we use 32 domains, this allows an average of 200 particles per cell which is said to be enough for the coalescence algorithm to converge [START_REF] Estivalezes | Simulations lagrangiennes polydisperses avec collisions et coalescence sur une configuration Ariane V 2D à l'instant 107s[END_REF].

Flow topology

Unsteadiness is visible on the gaseous vorticity field of the inert case (Fig. 2), where one retrieves a strong VSO and VSP along the wall, here tuned on the 1L chamber mode. Both hydrodynamic instabilities are strongly coupled to the dispersed phase. The structure is the same with coalescence. In the inert case we give Eulerian volume fractions in Fig. 3 and Lagrangian particle repartition in Fig. 4. The big droplets are segregated out of the vortices, while the small droplets fill most of the chamber. For the same reasons, small and big droplets flow differently throughout the nozzle. This dependence on size motivates a polydisperse approach. In the coalescing case and for the Eulerian simulation, the instantaneous volume fractions per section are in Fig. 5. This time, mass is transferred towards the third section, the big droplet mass being completely hosted in this section at the nozzle entrance. The Eulerian third section d 30 and the Lagrangian overall d 10 are presented as an illustration in Fig. 7. In both cases, the most intense coalescence zones are the edges of the vortices, the impact zone of the nozzle wall and the nozzle itself. Moreover coalescence can become so intense that droplet average sizes strongly vary through the domain. The droplet growth can however not be precisely compared since we have different average diameters on different size intervals. Finally Euler-Euler and Euler-Lagrange simulations give comparable droplet repartition and size. However quantitative comparison is made difficult by the radically different structure of output data.

Harmonic study

The pressure oscillation depends strongly on the flow so we can use it to quantify the effect of coalescence and the differences between the two approaches. The time evolution of the front and rear pressure signals is now described : when the particle injection starts, the gaseous flow encounters an additional inertia so the chamber pressure drops whatever the evolution of droplet sizes. But when the droplets reach the nozzle, the flow slows down so that the pressure builds up. It reaches its initial level in the non-coalescing but stabilizes about 2bar lower for coalescing particles. We now focus on this permanent regime where the average pressure does not vary anymore and only the instability remains. The main values of the pressure signal at permanent regime are given in Tab. 2. Both Eulerian and Lagrangian approaches well agree on average pressure in the inert case and they predict a lower mean pressure at permanent state when coalescence has occurred. This pressure drop can be related to the particle growth. Small particles tend to be in equilibrium with the gas and to increase the chamber pressure while, for the same loading, larger particles tend to stress the deviation from equilibrium and have reduced effects on the pressure. In the asymptotic limit, i.e. for very large particles, the two-phase pressure is similar to the single-phase one. Having said that, a pressure drop of two bars seems to be quite important and may be the sign that coalescence is over-estimated. But regarding oscillation levels, they strongly disagree in the inert case, though they well agree in the coalescing case. The difference is such that the ability to be predictive is still at stake. Further studies are required to understand the differences on instability levels between Eulerian and Lagrangian simulations and to complete the cross validation. The ease to set and exploit a computation is assessed here. Lagrangian and Eulerian approaches have common points but they differ in their treatment of the size variable and on the form of the output data. Regarding size parametrization, an equivalent diameter approach is unsuitable in SRM as shown in [START_REF] Doisneau | Two-way coupled simulation of acoustic waves in polydispersed coalescing two-phase flows : application to Solid Rocket Motor instabilities[END_REF] so it is compulsory to account for polydispersity. In this case of strongly bimodal distribution though, the size distribution is easy to model with two classes or two sections. Yet coalescence induces a size shift which has to be captured. The Multi-Fluid approach requires section bounds, which are determined from an a priori estimation of the droplet size evolution in the flow. This choice is generally difficult but any coarse but reasonable discretization is often satisfactory to get the dynamic trend as shown in [START_REF] Doisneau | Eulerian Multi-Fluid models for the simulation of dynamics and coalescence of particles in solid propellant combustion[END_REF]. The Lagrangian approach only requires enough particles to converge on the size distribution. Here again, this is difficult to determine a priori and it can only be evaluated thanks to numerical noise. Regarding coalescence computation now, Eulerian approaches compute deterministic integrals so they do not require particular care. On the contrary, Lagrangian collision algorithms are stochastic so they require sufficient particles per cell to converge. Finally most collision efficiency models can be used in both cases, as shown in [START_REF] Doisneau | Eulerian Multi-Fluid models for the simulation of dynamics and coalescence of particles in solid propellant combustion[END_REF] for the Eulerian method. Regarding post-processing, we highlight that Eulerian fields are convenient and mostly used so that Lagrangian output data require averaging operations.

Numerical diffusion assessment

Nonetheless is the flow better solved when the liquid structures are captured, but coalescence, which strongly depends on concentration, absolutely requires spatial accuracy. A transport scheme that is second order in space is currently developed in CEDRE for the Eulerian approach of the liquid phase. Since it has not been fully validated, it has not been used for the previous results but we have performed additional computations to make a qualitative comparison: the diffusion level improvement for the Eulerian approach compared to the Lagrangian approach is visible in Fig. 8, which shows a very good agreement on the particle structuration by the vortices for the high order Eulerian scheme. We insist on the fact that the finer the structures and the steeper the gradients are resolved, the more intense the coalescence and retrocoupling source terms are. The overall robustness of the two-phase numerical method is therefore challenged when higher accuracy is achieved for moderately dense sprays.

Particle Trajectory Crossing treatment with Eulerian methods

PTC is here considered for particles of the same size. It occurs in high strain zones, where the flow accelerates or curves faster than the particles can react. Critical Stokes numbers can be defined, depending on the configuration, above which PTC occurs. This has been tackled for Taylor Green vortices and Homogeneous Isotropic Turbulence [START_REF] De Chaisemartin | Polydisperse evaporating spray turbulent dispersion : Eulerian model and numerical simulation[END_REF][START_REF] Fréret | On the role of preferential segregation in flame dynamics in polydisperse evaporating sprays[END_REF][START_REF] Kah | Eulerian quadrature-based moment models for dilute polydisperse evaporating sprays[END_REF][START_REF] Vié | On the anisotropic gaussian closure for the prediction of inertial-particle laden flows[END_REF]. The Eulerian Multi-Fluid approach fails to treat PTC because of the monokinetic hypothesis [H1] [START_REF] De Chaisemartin | Polydisperse evaporating spray turbulent dispersion : Eulerian model and numerical simulation[END_REF]. In the P230 chamber, PTC does not occur very much on the axis and the unresolved PTC effect is moderate at the edge of the vortices so the Eulerian approach is satisfactory. In the nozzle however an unphysical δ-shock appears on the axis because of momentum averaging of the two symmetric particle jets. Moreover δ-shocks densify as the mesh refines or as the scheme's diffusion decreases, which is well illustrated in Fig. 8 and is a major issue of both numerical convergence and physical meaning. Specific impulse loss can therefore be considered as mis-estimated. The monokinetic hypothesis generates the same kind of limitations at the wall than for PTC : the disperse phase can only slip along the wall with no normal velocity component because of momentum averaging. So the PTC failure of Eulerian mono-kinetic approaches is a major drawback for certain SRM applications at the wall and in the nozzle. It can only be treated by implementing a polykinetic approach, either with multiple velocities [START_REF] Kah | Eulerian quadrature-based moment models for dilute polydisperse evaporating sprays[END_REF][START_REF] Yuan | Conditional quadrature method of moments for kinetic equations[END_REF] or with pressure-like terms that account for an average effect of PTC [START_REF] Simoes | Modélisation eulérienne de la phase dispersée dans les moteurs à propergol solide, avec prise en compte de la pression particulaire[END_REF][START_REF] Vié | On the anisotropic gaussian closure for the prediction of inertial-particle laden flows[END_REF], depending on the range of Stokes numbers.

Computational cost and parallelizability

Finally, let us consider the cost of the two approaches as implemented in CEDRE on a complex SRM case. They are given in Table 3 and appear to be comparable. The major difference, highlighted here, is the natural scalability of the Eulerian method compared to the difficulty to increase the number of CPUs for Lagrangian computations as said previously. This is a major advantage of Euler-Euler approaches for two-phase flows. The cost of such simulations becomes indeed prohibitive on 3D meshes [START_REF] Estivalezes | Simulations lagrangiennes polydisperses avec et sans collision-coalescence sur une configuration Ariane V tridimensionnelle à l'instant 107s[END_REF]. Moreover, convergence on particle repartition and two-way coupling is naturally achieved with Eulerian approaches while no criterion on particle number or integration time step is known for Lagrangian approaches.

Conclusion

This paper brings three salient facts to light. First, considering coalescence in aluminized SRM simulations is incontrovertible because the droplet sizes impact on the liquid repartition, on instability levels, on specific impulse loss, and probably on slag weight. Second, the Euler-Euler approach is valid when compared to Lagrangian computations even though it has modeling weaknesses such as PTC treatment, as long as PTC occurs moderately. Third, the Euler-Euler approach is accessible and recommendable for two-way coupled simulations, for parallelization purpose. Finally, further validations are required to completely match with Lagrangian references in such complex configurations. Then the model and method can bear a transposition to LES in order to achieve accessible and predictive simulations.

  wall injection (d 1 and d 2 )

  wall injection (d 1 and d 2 )

Figure 1 .

 1 Figure 1. 2D axisymmetric geometry and domains of the P230 SRM (• sensor location).

Figure 2 .

 2 Figure 2. Vorticity (rad/s) in the P230 at t = 1.46s without coalescence.

Figure 3 .

 3 Figure 3. Volume fraction per section (Top : d 1 ; Bottom : d 2 ) at t = 1.09s without coalescence.

Figure 4 .

 4 Figure 4. Lagrangian parcel repartition colored by diameter (m) at t = 1.09s without coalescence.

Figure 5 .

 5 Figure 5. Volume fraction per section (Top : d 1 ; Middle : d 2 to d 3 ; Bottom : above d 3 ) at t = 0.81s with coalescence.

Figure 6 .

 6 Figure 6. Lagrangian parcel repartition colored by diameter (m) at t = 0.81s with coalescence.

Figure 7 .

 7 Figure 7. Droplet sizes with coalescence at t = 0.46s -Top : Third section d 30 diameter (m) for the Eulerian simulation; Bottom : Overall d 10 diameter (m) for the Lagrangian simulation.

Figure 8 .

 8 Figure 8. Volume fractions at t = 0.278s without coalescence -Top : d 1 and d 2 with first order transport scheme; Middle : d 1 and d 2 with second order transport scheme; Bottom : both sizes with Lagrangian.

Table 1 .

 1 Timescale hierarchy (µs), exhibiting stiffness.

	d 1 droplet Stokes time	1
	Acoustic CFL time	10
	Convective CFL time (nozzle)	30
	d 2 droplet Stokes time	10, 000
	First acoustic mode	50, 000
	Eddy revolution time	50, 000
	Convective CFL time (injection)	90, 000
	Typical computation time	1, 000, 000

Table 2 .

 2 Non-dimensionned pressure data (Lagrangian values from ONERA simulations[START_REF] Estivalezes | Simulations lagrangiennes polydisperses avec collisions et coalescence sur une configuration Ariane V 2D à l'instant 107s[END_REF]).

			Without coalescence	With coalescence
			Lagrangian	Eulerian	Lagrangian	Eulerian
	Frequency		1.00 ± 1%	1.02 ± 4%	0.99 ± 2%	1.10 ± 4%
	Front sensor Average	1.00	0.99	0.98	0.95
	pressure	RMS	1.0	4.7	1.8	1.6
	4. Eulerian VS Lagrangian computations
	4.1. Ergonomy of both approaches

Table 3 .

 3 Computational cost per time step (32 Nehalem CPU).

		Eulerian	Lagrangian
	Inert	84ms	124ms
	Coalescence	144ms	236ms

Acknowledgements

This research was done thanks to a Ph. D. grant from DGA, Ministry of Defence (M. S. Amiet, Technical Monitor) and a France-Stanford grant from the French Ministry of Foreign Affairs (M. Massot, Technical Monitor). The authors also thank J.-L. Estivalezes (ONERA) for preliminary Lagrangian computations and SME for the P230 mesh.


Nicolas Martin

Jérémie Sauget

Tomas Nyberg

Anisotropic electrical resistivity during annealing of oriented columnar titanium films

Keywords: GLAD, titanium, resistivity, anisotropy, annealing

We report on the evolution of anisotropic electrical resistivity versus temperature of titanium thin films. An oriented titanium film (1 µm thick) is sputter deposited by GLancing Angle Deposition (GLAD) using an incident angle α = 80° of the particle flux. Two parallel Ti electrodes cover this GLAD film. We measure the components of the conductivity tensor by the van der Pauw method during an annealing treatment in air ranging from 298 to 873 K.

The average DC electrical resistivity ρ changes from 6.03×10 -6 to more than 2.40×10 -1 Ωm with the increasing temperature. The anisotropy ratio is A = 1.39 before annealing and reaches 12.4 for the highest temperatures. This enhanced anisotropy is interpreted assuming the oxidation of the porous GLAD titanium film.

Introduction

DC electrical resistivity of metal and semiconductor thin films can be determined by the van der Pauw (vdP) method [1]. It is a convenient method, which requires only four contacts sufficiently small and located at the circumference of the sample. This powerful method is based on averaging measurements of the potential difference between two electrodes when a known current flows between two other electrodes. One of the strong hypotheses of this method assumes an isotropic conductive medium. However, most thin films produced by vacuum processes, often exhibit anisotropic behaviors. In order to get electronic transport properties in such anisotropic thin films, components of the resistivity tensor have to be determined. For this purpose, the theory of the vdP method has been extended [2][3][4][5]. Price [6] showed that the method may be applied to anisotropic materials provided that one of the principal axes of the resistivity tensor is perpendicular to the plane of the sample. Recently, Bierwagen et al. [7] proposed an original procedure to obtain anisotropic transport properties, especially from the results of the vdP measurements. This letter demonstrates that the anisotropic DC electrical resistivity of titanium thin films can be measured from the Bierwagen et al. method during an annealing treatment. a) Corresponding author: nicolas.martin@femto-st.fr; Tel.: +33 (0)3 81 85 39 69; Fax: +33 (0)3 81 85 39 98

Material and methods

The titanium films were deposited on glass by DC reactive magnetron sputtering from a titanium target (purity 99.6 %) in a pure argon atmosphere using a sputtering pressure P sput = 0.10 Pa. The target was sputtered with a current density J Ti = 100 Am -2 . Substrates were grounded and unheated during deposition. A titanium film 1 µm thick was first deposited on the glass substrate by GLancing Angle Deposition (GLAD) [8] using an incident angle α of the particle flux of 80°. An oriented columnar structure was produced with a column angle β close to 40°. The latter is systematically lower than the incident angle because of the shadowing effect produced by the growing columns [8]. Afterwards, two parallel strips extending along the x axis (by means of a mask) of titanium 500 nm thick and separated by 9 mm were deposited by conventional sputtering (α = 0°). Two probes were placed on each titanium strip leading to a square geometry with L x = L y = 10 mm (Fig. 1). Then, the vdP configuration was used to measure the average and anisotropic electrical resistivity. The annealing treatment in air was performed following 9 cycles of temperature. The first cycle started from 298 up to 473 K at 5 K min -1 followed by a return to 298 K with the same ramp of 5 Kmin -1 . The second cycle was similar to the previous one up to 523 K (increment of 50 K). The third cycle was up to 573 K and so on until the last cycle (9 th ) reached 873 K to finally come back to 298 K.

Results and discussion

Assuming that the sides of the sample are parallel to the principal axes of the tensor, the conductivity tensor can be expressed in the following form:

      = yy xx σ σ σ 0 0 ( 1 )
Where σ xx (Sm -1 ) and σ yy (Sm -1 ) are the principal conductivities in the mutually perpendicular directions x,y of the two-dimensional plane, respectively. In order to get these σ xx and σ yy components, the resistances R xx (Ω) and R yy (Ω) are first measured from:

I V I V R R R xx + = + =
I V I V R R R yy + = + = (2)
Where V 23 (V) is the voltage measured across contacts 2 and 3 (inverse for V 32 ) and I 14 (A) is the current flowing through the contacts 1 and 4 (inverse for I 41 ), and analogously for R yy . The average DC electrical resistivity ρ (Ωm) can be calculated using the classical vdP approach with [1]:

1 exp exp =         - +         - ρ π ρ π yy xx R d R d ( 3 )
Where d (m) is the thickness of the sample. The anisotropy in our vdP configuration A vdP is defined as [7,9]:

( ) [ ] ( ) [ ] ∑ ∑ + + = - = - - × × = = odd n eff odd n eff xx yy vdP n A n n A n R R A 1 1 1 sinh sinh π π ( 4 )
This equation ( 4) is used to calculate the effective anisotropy A eff . The anisotropic coefficient A is finally deduced from:

2 yy y xx eff xx x y y L AA L ρ σ ρ σ  =× = =   ( 5 )
And the σ xx and σ yy components are obtained with:

1 1 - × = = A xx xx ρ ρ σ and A yy yy × = = ρ ρ σ 1 ( 6 )
Where ρ xx (Ωm) and ρ yy (Ωm) are the resistivities in the mutually perpendicular directions x,y of the twodimensional plane, respectively.

The average DC electrical resistivity ρ was first measured as a function of the annealing treatment (Fig. 2).

Before heating, the resistivity at room temperature is ρ 298 = 6.03×10 -6 Ωm, which is more than one order of magnitude of the bulk titanium (ρ 298 bulk = 3.91×10 -7 Ωm). It is mainly attributed to the scattering of electrons at the grain boundaries [10]. Due to the high porous structure of the GLAD films [11], scattering is particularly favored at the column interfaces reducing the electronic transport properties. A typical metallic-like behavior is measured as the temperature rises up to T = 473 K since ρ 473 linearly reaches 8.15×10 -6 Ωm. A linear interpolation of ρ versus T is calculated up to 370 K, leading to a temperature coefficient of resistance (TCR) of 1.19×10 -3 K -1 , which is again lower than that of the bulk value (TCR bulk = 5.50×10 -3 K -1 ).

A further increase of the temperature up to 473 K leads to a non linear evolution of the resistivity versus temperature since ρ 473K reaches 8.11×10 -6 Ωm before the temperature drops. For the first cycle, the return to the room temperature gives rise to a linear evolution of ρ versus T and the resistivity is down to ρ 298 = 7.24×10 -6 Ωm before starting the second cycle. This increase of resistivity is mainly attributed to the initial oxidation of the GLAD porous titanium film. Further annealing treatments (2 nd to 5 th cycles) reduce the metallic-like character of the film since the resistivity reaches ρ 298 = 1.83×10 -5 Ωm at the end of the fifth cycle. In addition, the TCR decreases and finally becomes negative (TCR = -1.07×10 -4 K -1 ) due to the emphasized oxidation of the GLAD film. We can suggest that metallic titanium monoxide is produced since negative TCRs are obtained in TiO x films (0.80 < x < 1.23) containing a high concentration of defects [12]. Assuming the Mooij criterion [13], the increasing resistivity vs. annealing cycle is mainly due to a disordering effect (coming from the oxidation), inducing a reduction of the electron mean free path from several to a few interatomic spacings. This trend is even more significant for the 8 th and 9 th cycles, for which resistivity values and their evolution as a function of the temperature correspond to a semi-conducting-like behavior. At the end of the final cycle (9 th ) and as the temperature decreases, ρ 500 reaches 2.40×10 -1 Ωm. Below 500 K the film is too resistive to be accurately measured with our vdP system. This last annealing cycle leads to the most important jump of the electrical resistivity since 3 orders of magnitude separate resistivity values measured between the 8 th and 9 th cycle (axis break in Fig. 2). This smooth metal-to-semiconductor transition has previously been reported for other compounds [14], but rarely for GLAD thin films [15]. It correlates with a change from a low resistivity behavior (metallic-like) of oxygen-rich TiO 1+ε exhibiting a negative TCR [16], to a semiconducting-like behavior observed for oxygen-deficient TiO 2-δ with a low activation energy [17]. In this study, this transition is mainly ascribed to the oxidation of the GLAD porous titanium film, which is emphasized for the 8 th cycle, and becomes even more significant at the final annealing treatment. As a result, it can be claimed that an oxygendeficient titanium dioxide compound is formed. The generalized model of conduction proposed by Goldfarb et al. [18] for amorphous transition-metal oxides is then appropriate to interpret this increase of resistivity vs. annealing treatment. A continuous oxidation causes depletion of the Ti d-band charge carriers in favor of the O 2p valence band, inducing a transition from metallic to hopping conduction, and finally down to the percolation threshold. This oxidation occurs up to a point where the hopping sites tend to be few and so distant from each other that only percolation conduction can take place. In the end, the concentration of hopping sites falls below the threshold for percolation and the resistivity becomes too high to be measured by the vdP technique.

The conduction anisotropy during the annealing treatment has been investigated from calculations of the ρ xx and ρ yy components using equations (1) to (6). Their evolution and the anisotropic coefficient vs. annealing cycle number are measured at 473K (Fig. 3). It is worth noting that the resistivity ρ xx along the x-direction is systematically lower than that of the y-direction ρ yy for any annealing cycle. This difference is mainly related to the GLAD Ti film microstructure. Since a very high porous architecture is especially favored for grazing incident angles of the particle flux (α > 70°), the number of electrical pathways available for the free carriers are reduced leading to a decrease of the electron mean free path. In addition, GLAD deposition of columnar films with incident angles close or higher than 60° gives rise to a significant structural and uniaxial anisotropy in the substrate plane induced by the shadowing effect [19]. This effect mainly occurs in the direction of the incident vapor flux. Islands grow with connections to each other by chains perpendicular to the plane of incidence or to the direction of shadowing [20].

Films produce a columnar structure with columns exhibiting an elliptical section versus the incident angle, which is particularly marked for α ≥ 60°. As a result, ρ xx = 6.90×10 -6 Ωm and ρ yy = 9.62×10 -6 Ωm are measured with a corresponding anisotropic coefficient A = 1.39 for the first cycle. These values seem to be contradictory with results previously obtained by others [21,22], who reported the lowest resistivity in the particle flux direction (y-direction) for GLAD metallic films. Besides some discrepancies are related to the measurement method and the definition of the anisotropic coefficient, the authors also reported that the anisotropic transport properties depend on the local microstructure of the films and the nature of deposited materials. However, it should be reminded that our system is composed of Ti stacks where two parallel strips of Ti deposited by conventional sputtering (α = 0°) cover the GLAD Ti film (α = 80°). Thus, resistivity of our system according to the x-direction is systematically lower than that of the y-direction.

Until the 7 th cycle, ρ xx and ρ yy gradually increase up to 1.66×10 -5 and 2.64×10 -5 Ωm, respectively and inversely the anisotropic coefficient decreases down to A = 1.59. This range of resistivities corresponds to a metallic-like behavior in both x and y directions in spite of a negative TCR value calculated from the average resistivity. The oxidation only occurs on the surface of the conventional Ti film (α = 0°) whereas it becomes significant through the thickness of the GLAD Ti film due to its high porous structure [15]. It supports the formation of a TiO 1+ε -like compound since the resistivity ρ yy remains in the order of magnitude of metallic compounds.

The sudden increase of the average DC electrical resistivity ρ previously measured from the 8 th cycle (Fig. 2) correlates with the beginning of the anisotropy change (Fig. 3). The anisotropic coefficient rises to A = 12.4 as the 9 th cycle ends and resistivities reach ρ xx = 6.54×10 -2 and ρ yy = 8.13×10 -1 Ωm. As a result, these last annealing cycles do not solely lead to a smooth metal-to-semiconductor transition. They also emphasize the anisotropy of the electronic transport properties. Oxidation phenomenon proceeds in the GLAD Ti thin film and becomes substantial in the conventional Ti thin film strips. Taking into account the fact that the oxidation kinetics of pure titanium becomes especially significant from about 830 K [23], the conventional Ti film is certainly transformed into a mixture of TiO and TiO 2 compounds since the temperature of the final cycle reaches 873 K. In addition, oxidation strongly depends on the structure of the film, which is especially favored by voids in the GLAD film. For our titanium stack films, this leads to a reinforcement of the anisotropic conductivity as the temperature rises. This trend agrees with the few reports devoted to the role of temperature on the anisotropic electrical transport in thin films [24,25]. This dependence has been attributed to the anisotropic distribution of the grain boundary (GB) potential barrier heights. A similar approach can be used for our films based on the GB trapping model [26]. Anisotropic growth due to the GLAD process intrinsically leads to the highest GB potential barriers in the y-direction in comparison with the x-ones. The annealing in air enhances this effect, especially in the porous GLAD Ti film which is very sensitive to the oxidation phenomenon.

Conclusions

Anisotropic electrical resistivity vs. temperature of titanium thin films has been measured. An extension of the van der Pauw method has been applied to get resistivities ρ xx and ρ yy in the mutually perpendicular directions x,y of the two-dimensional plane and the anisotropic coefficient during incremental annealing cycles in air (from 298 to 873K). As-deposited films exhibit a metallic-like behavior with an average DC resistivity ρ = 6.03×10 -6 Ωm and an intrinsic anisotropic coefficient A = 1.39. The metallic character remains as the temperature rises up to 773 K and the anisotropy is slightly enhanced to A = 1.59 assigned to the formation of a TiO compound in the GLAD Ti film.

Similarly, the average DC resistivity is close to 2×10 -5 Ωm and the TCR reduces to become negative (-1.07×10 -4 K - 1 ). The metallic character remains for this range of temperatures in spite of the oxidation of the GLAD Ti film, which initially occurs because of its highly porous structure. An increase of the temperature up to 800 K and finally 873 K induces a further oxidation of the titanium stack films. Resistivities ρ xx and ρ yy are especially increased from 800 K to show semiconducting-like behaviors with an emphasized anisotropy coefficient A = 12.4. This anisotropic electrical resistivity is attributed to difference of the grain boundary potential barriers in the x and y directions produced by the anisotropic growth of the GLAD process. This effect is promoted by the annealing treatment in air, especially in the porous GLAD Ti film. Average DC electrical resistivity ρ as a function of the temperature T of the titanium stack films. The annealing treatment in air was performed following 9 cycles of temperature. The first cycle started from 298 up to 473 K followed by a return to 298 K with the same ramp of 5 Kmin -1 . The second cycle was similar to the previous one up to 523 K (increment of 50 K) and so on until the final cycle (9 th ) reached 873 K to finally come back to 298 K.

Figure captions

The annealing dependence of orthogonal DC electrical resistivities ρ xx and ρ yy at 473 K and anisotropic coefficient A calculated from equation (5) of titanium stack films. The conduction anisotropy is especially emphasized from the 8 th cycle, which corresponds to an annealing temperature close to 800 K.

Figure 1 A

 1 Figure 1A schematic representation of the titanium stack films and the setup for the DC electrical resistivity measurements.The GLAD Ti film 1 µm thick is deposited by sputtering with an incident angle α = 80° of the particle flux. A crosssection view by scanning electron microscopy exhibits the inclined columnar structure. The two parallel strips of Ti (500 nm thick and separated by 9 mm) are deposited by conventional sputtering (α = 0°). The four contacts are set according to a square geometry of size L x = L y = 10 mm.

Figure 2

 2 Figure 2

Figure 1 Figure 2

 12 Figure 1

Figure 3

 3 Figure 3

Acknowledgements

The authors thank the Region of Franche-Comté for financial support. This work has been supported by the Labex ACTION program (contract ANR-11-LABX-01-01).