Patrick Dutournie

Jacques Mercadier
email: jacques.mercadier@univ-pau.fr

Cyril Aymonier

Aure

Lia Gratias

Franc ¸ Ois Cansell

Determination of Hydrothermal Oxidation Reaction Heats by Experimental and Simulation Investigations

The purpose of this study is to develop a procedure for the prediction of oxidation reaction parameters. The primary aim will be the determination of reaction heats and kinetics for any oxided wastewater substance. This even includes unknown waste compositions. The strategy of the research is, first of all, to build an experimental nearly adiabatic reactor for the measurement of temperature profiles during the reaction. The second part is the numerical simulation of the process, where two procedures are proposed for the determination reaction heats and kinetic parameters. To validate the overall project, two model oxidations were used, namely, the oxidations of acetic acid and methanol by hydrogen peroxide. Results obtained for the modeled molecules are close to those of other researchers, which leads to the conclusion that the developed model works correctly and will be applicable as well for different waste compositions.

Introduction

The evolution of environmental regulations and the increasing wastewater disposal costs lead to the development of new concepts for the complete destruction of toxic substances and sludges. The hydrothermal oxidation of wastes is developed as an alternative technique in order to limit the risks of secondary pollution and the energy supply. This process can be considered as a clean technology.

Hydrothermal oxidation concerns the oxidation reaction of organic matter in water under high pressure and high temperature. Supercritical water oxidation (SCWO, P g 22.1 MPa and 374 °C e T e 700 °C), introduced by M. Modell, 1 is based on high reactivity under these conditions because of the unique properties of supercritical water. Some of these properties include a low dielectric constant, low viscosity, high diffusivity, and complete miscibility of gases such as O 2 and organic compounds, which promotes high mass transport rates with no interphase transport limitations. This leads to high destruction rates (>99.99%) in relatively short reaction times (few seconds). The salt solubility decrease is important because of decreasing density and dielectric constant. When the operating conditions are optimized, the end products are water, carbon dioxide, nitrogen, and eventually mineral acids.

The development of simulation tools is necessary for scaling up SCWO processes. Different approaches are known. Some authors [START_REF] Oh | Numerical analysis and data comparison of a supercritical water oxidation reactor[END_REF][START_REF] Oh | Modelling of thermal characteristics in supercritical water oxidation reactors[END_REF][START_REF] Dutournie | Simulation of a tubular reactor for supercritical water oxidation[END_REF] use trading simulation computer softwares (such as Modar or Fluent) to describe supercritical oxidation in specific reactors. Other works [START_REF] Chen | Simulation of a concentrictube reactor for supercritical water oxidation[END_REF][START_REF] Petrich | Model and simulation of supercritical water oxidation[END_REF][START_REF] Dutournie | One-dimensional simulation of a tubular reactor for waste oxidation[END_REF] devoted to supercritical water oxidation simulations based on one-dimensional models concern the modeling of stationary or nonstationary tubular reactors in turbulent flow. Knowledge of kinetic parameters and reaction heats is very important for the development of simulation tools. Currently, there is a great deal of data in the literature on kinetic parameters, for example, the works of Eckert et al., [START_REF] Eckert | Tuning fluid solvents for chemical reactions[END_REF] Gloyna et al., [START_REF] Li | Kinetic model for wet oxidation of organic compounds in subcritical and supercritical water[END_REF] Savage et al. [START_REF] Martino | Supercritical Water Oxidation Kinetics, Products, and Pathways for CH3-and CHO-Substituted Phenols[END_REF] and Tester et al. 11. Nevertheless, works on the determination of reaction heats are very limited. P. Chen et al. [START_REF] Chen | Simulation of a concentrictube reactor for supercritical water oxidation[END_REF] propose a hydrothermal oxidation reaction heat of -435 kJ mol -1 for most organic compounds and -870 kJ mol -1 for acetic acid.

Some authors calculate formation enthalpies with thermodynamic tables for simple compounds. However, few compounds are represented in this pressure and temperature domain in these tables.

Enthalpies can be estimated with computer software packages such as aspen plus [START_REF] Oh | Numerical analysis and data comparison of a supercritical water oxidation reactor[END_REF][START_REF] Oh | Modelling of thermal characteristics in supercritical water oxidation reactors[END_REF] or Prophy Plus. [START_REF] Prosim | Prosim, 132 route d′Espagne[END_REF] These software programs calculate physical properties of the chemical compounds and determine enthalpy variations between inlet reagents and outlet products of reactions. These calculations are realized with an equation of state valid in this temperature and pressure region. For the hydrothermal oxidations of methanol and acetic acid with oxygen, Prophy Plus gives, respectively, -686 and -847 kJ mol -1 .

However, these software programs do not permit the calculation of some reaction heats of real wastes because their compositions are unknown or very complex. Because of these problems, we have developed a new approach to the determination of reaction heats and kinetic parameters of hydrothermal oxidation reactions.

In this paper, we present the ICMCB pilot plant facility equipped with a quasi-adiabatic reactor, which permits temperature profiles along the tubular reactor to be obtained. Starting from these experimental temperature profiles, reaction heats and kinetic parameters are calculated with two numerical programs developed in this study. Then, we discuss the results obtained for two model molecules: acetic acid and methanol.

Experimental Section

ICMCB Pilot Plant Facility. The ICMCB pilot plant facility is able to treat 3 kg/h of aqueous wastes in a pressure range from 0.1 to 50 MPa and a temperature range from 25 to 600 °C. Figure 1 is a schematic view of this pilot plant.

This pilot plant facility is composed of a waste injection line with an high-pressure pump and an oxidant injection line (hydrogen peroxide diluted in water) with another high-pressure pump (LEWA). These two lines go through preheaters and join at the beginning of the reactor. At the exit of the reactor, the effluent is brought to room temperature and atmospheric pressure. The resulting liquid and vapor phases are separated in a gas/liquid separator.

The composition of the gas phase is analyzed by online gas-phase chromatography. The organic matter concentration in both the liquid feed and the liquid effluent is determined by measurement of the chemical oxygen demand (COD) with a colorimetric method. [START_REF] Aymonier | Hydrothermal oxidation of a nitrogen-containing compound: the fenuron[END_REF] The experimental procedures for the pilot plant, reproducibility data, and experimental error are detailed in the work of P. Beslin. 14 Quasi-Adiabatic Reactor. The quasi-adiabatic reactor has been developed to determine heat quantities and kinetic parameters (activation energies and preexponential factors of the Arrhenius law) [START_REF] Aymonier | Hydrothermal treatment of waste and environmental and energetic acceptable process[END_REF][START_REF] Aymonier | Hydrothermal oxidation of Waste: an environmental acceptable process[END_REF][START_REF] Dutournie | Estimation of reaction heat by supercritical reactor simulation[END_REF] of hydrothermal oxidation reactions of model molecules and of industrial wastes.

The quasi-adiabatic reactor is made of Inconel 718 tube of 2.4 mm inside diameter and 8.4 m length. It is equipped with 28 external thermocouples (type K) positioned along the reactor.

The tube is isolated with ceramic fibers and surrounded with a thermal shield. The thermal shield is constituted of a copper sheet, surrounded by an electrical heating material in order to reduce thermal loss. The temperature of the thermal shield is uniform and regulated at the temperature of the mixing point between the waste feed solution and the hydrogen peroxide feed solution.

The quasi-adiabatic reactor allows the temperature profile along the tube to be obtained, which characterizes the exothermic oxidation reaction.

Procedure of Reactor Temperature Profile Treatment. A program has been developed to treat reactor temperature profiles. This one-dimensional program where the term of thermal expansion represents the internal energy variations of a compressible fluid subject to a pressure gradient. In eqs 3 and 4, the reaction rate term (A e -Ea/RT C) represents a first-order reaction with respect to the organic molecule. These hypotheses are used by most of the authors. [START_REF] Cansell | Fluides Supercritiques et Mate ´riaux[END_REF] The first method, the so-called "integral method", for the determination of the reaction heat involves the integration of the experimental temperature profile from the association of eqs 3 and 4. The first integration shows that the effects of the thermal viscous dissipation and the pressure gradient were negligible. The molecular diffusion characteristic time was very high with respect to the convective transport characteristic time.

In neglecting these terms, the reaction heat can be expressed as follows: This program solves the oxidation of model compounds in an adiabatic reactor at constant pressure or in a reactor cooling by convective heat transfer. It allows for the simulation, in particular, of the temperature profile along the reactor. The quantity of organic compounds in the reactor being less than 5 wt %, the mixture is considered to be pure water for the thermodynamic behavior of the fluid. [START_REF] Schmidt | Properties of Water and Steam in SI Units[END_REF] The modeling of the flow field in the reactor is described by the mass, momentum, species, and energy conservation equations.

where the values of C p were considered to be those of pure water, [START_REF] Dutournie | Estimation of reaction heat by supercritical reactor simulation[END_REF] h was determined from heat losses along the tubular reactor without chemical reaction, and C was calculated from chemical analysis at the reactor inlet and at the reactor outlet.

∆H values were calculated from eq 5 and from the obtained experimental temperature profiles. From the ∆H value, the concentration profile, C(z), was determined by the integration of eq 5 from 0 to L. The kinetic rate was then calculated with eq 3. Then, the activation energy was obtained from an Arrhenius plot.

We can also determine the reaction heat and the kinetic parameters (activation energy and preexponential factor) by using the simulation method, which solves eqs 1-4 simultaneously. This second method, the socalled "simulation method", calculates simulated profiles of u, C, and T as functions of z by Euler's method. Because of the considerable variations in the water properties, the program is iterative and uses relaxation coefficients. First, the adjustment of the reaction heat is performed with experimental data for total conversion (independent of the reaction rate). Then, the kinetic parameters are adjusted in order to superpose the simulated [START_REF] Aymonier | Hydrothermal oxidation of Waste: an environmental acceptable process[END_REF][START_REF] Dutournie | Estimation of reaction heat by supercritical reactor simulation[END_REF] and experimental profiles.

Acetic acid and methanol were studied as model molecules with these two treatment procedures.

Results and Discussion

Determination of Reaction Heat and Kinetic Parameters of Acetic Acid SCWO. The kinetic rate expression is given as follows:

where [CH 3 COOH] is the molar concentration of acetic acid (mol m -3 ).

The unpreheated hydrogen peroxide was brought into the reactor in stoechiometric quantity according to the following equation:

The reactor is able to treat aqueous wastes in the range from 380 to 600 °C for supercritical oxidation. Indeed, the reactor material could be damaged for temperatures higher than 600 °C. Under 400 °C, the kinetic energy is too low to start the reaction. Thus, the chosen inlet temperature is around 400 °C. Inlet concentrations have then been selected so that different conversion yields can be obtained (20-100%). (See Table 1.)

The experimental temperature profiles obtained for tests 1-3 are presented in Figure 2. The mass flow rate of test 2 is more important than that of test 1. Hence, the average residence time of waste in the reactor is different around 34 s for test 1 and 24 s for test 2. The shield temperature profiles in tests 1 and 2 are different. For these reasons, the temperature profile of test 1 decreases at the end of the reactor, even if the conversion is not complete. Test 3 has an inlet temperature and concentration more important than those of tests 1 and 2.

From these experimental temperature profiles, the reaction heat and kinetic parameters have been determined with the integral and simulation methods. The results are given in Table 2.

The results obtained with the integral method and simulation method are very close and lead to the average results: ∆H ) -909 kJ mol -1 , E a ) 170 kJ mol -1 , and A ) 85 × 10 9 s -1 .

The comparison between experimental and simulated temperature profiles is presented in Figure 3.

The obtained kinetic parameters are similar to those obtained by Lee and co-workers [START_REF] Lee | Efficiency of H 2O2 and O2 in Supercritical Water Oxidation of 2,4-Dichlorophenol and Acetic Acid[END_REF] (for the oxidation of acetic acid by hydrogen peroxide). Indeed, they studied the SCWO of acetic acid in a flow reactor using the isothermal plug flow design equation. They found for the working conditions 22.5 e P e 31 MPa and 415 e T e 525 °C, a preexponential factor of 3.16 × 10 11 s -1 and an activation energy of 180 kJ mol -1 .

Determination of Reaction Heat and Kinetic Parameters of Methanol SCWO. The operating procedure was the same as previously described. The unpreheated hydrogen peroxide was brought into the a Tinlet is the temperature at the mixing point between the acetic acid feed solution and the hydrogen peroxide feed solution. b Dw and Dox are given at 0.1 MPa and 25 °C. c The mass fraction is the mass ratio of acetic acid in the initial feed solution (Dw). d ∆COD is defined by: ∆COD ) reactor in stoechiometric quantity according to the following equation:

∆H ) ∫ 0 L FuC p dT dz dz + ∫ 0 L4h φ (T -T ∞ ) dz ∫ 0 Ld(uC) dz dz (5) k ) A e -E a /RT [CH 3 COOH] (6) CH 3 COOH + 4H 2 O 2 f 2CO 2 + 6H 2 O (7)
(D w × COD i ) -(D w + D ox ) × COD f D w × COD i ( 8
)
The results are presented in Table 3. We performed calculations by the simulation and integral method obtained the following average results: ∆H ) -670 kJ mol -1, E a ) 57.6 kJ mol -1 , and A ) 4250 s -1 .

The simulation results are very close to the experimental data (Figure 4). The reaction heats obtained with methanol at 6 and 3 wt % in the initial solution are approximately equal. The activation energy of 57.6 kJ mol -1 and the preexponential factor of 4250 s -1 were similar to those obtained by P. Beslin and co-workers [START_REF] Cansell | Hydrothermal oxidation of model molecules and industrial wastes[END_REF] (for methanol oxidation by hydrogen peroxide).

Conclusions

As the overall purpose of our study was to develop a procedure for predicting oxidation reaction parameters, we built a quasi-adiabatic reactor in order to obtain its temperature profile during the hydrothermal oxidation reactions of organic compounds.

Furthermore, we developed two numerical procedures that use the experimentally obtained temperature profiles to determine reaction heats and kinetic parameters. Reaction heats of -910 and -670 kJ mol -1 were calculated for the hydrothermal oxidations of acetic acid and methanol, respectively. Preexponential factors were also determined for acetic acid (8 5 × 10 10 s -1 ) and methanol (4250 s -1 ). For the activation energies, we found 170 kJ mol -1 (for acetic acid) and 57.6 kJ mol -1 (for methanol). These data allow for improvements in the reactor design. Furthermore, the obtained reaction heats permit the calculation of the energetic balance of the processes under additional operating conditions. Acknowledgment L'Electrolyse is gratefully acknowledged by the authors for financial support for C.A. P.D. gratefully acknowledges the financial support provided by Anjou Recherche, CNIM, and ECODEV.

Notation

A ) preexponential constant (s -1 ) C ) molar concentration of waste (mol m -3 ) COD ) chemical oxygen demand (mg/L) C p ) specific heat (J kg -1 K -1 ) D ) diffusion coefficient of compound in the mixture (m 2 s -1 ) D w ) mass flow rate (kg s -1 ) E a ) activation energy (kJ mol -1 ) a The mass fraction is the mass ratio of acetic acid in the initial feed solution (Dw). b Tinlet is the temperature at the mixing point between the acetic acid feed solution and the hydrogen peroxide feed solution.

Figure 1 .

 1 Figure 1. ICMCB pilot facility for hydrothermal waste treatment.

Figure 3 .

 3 Figure 3. Experimental and simulated temperature profiles of hydrothermal oxidation of acetic acid.

Figure 4 .

 4 Figure 4. Experimental and simulated temperature profiles (°C) for test 1 and test 2 vs the reactor length (m)Table 3. Methanol Hydrothermal Oxidation at 25 MPa

Table 1 . Operating Conditions and Results of Hydrothermal Oxidation of Acetic Acid at 25 MPa

 1

				mass				shield
		Tinlet	Dw/Dox b	fraction c	CODi	CODf	∆COD d	temp
	test	(°C)	(kg/h)	(%)	(mg/L)	(mg/L)	(%)	(°C)
	1 396 0.31/0.23	6	63679 27440 24.94 400
	2 396 0.41/0.32	8	83389 34650 26.01 400-440
	3 403 0.25/0.26	10	105657 1586 96.94 400-440

Table 2 . Reaction Heats and Kinetic Parameters for Acetic Acid Hydrothermal Oxidation

 2

		simulation method	integral method
		∆H	Ea	A	∆H	Ea
	test	(kJ mol -1 )	(kJ mol -1 )	(s -1 )	(kJ mol -1 )	(kJ mol -1 )
	1	-909	166.8	87 × 10 9	-913	171.9
	2	-910	170.7	82 × 10 9	-902	170
	3	-911	171.2	85 × 10 9	-912	169.3
	Figure 2. Experimental temperature profiles in terms of reactor
	length.


Keywords: Amokrane, Ayeley Tchangani, François

published or not. The documents may come

Evaluation process in end-of-life systems management using BOCR analysis.

This is an author-deposited version published in: http://oatao.univ-toulouse.fr/ Eprints ID: 10694

To cite this version:

Bouzarour-Amokrane, Yasmina and Tchangani, Ayeley and Peres, François Evaluation process in end-of-life systems management using BOCR analysis.

Open Archive Toulouse Archive Ouverte (OATAO)

OATAO is an open access repository that collects the work of Toulouse researchers and makes it freely available over the web where possible.

Any correspondence concerning this service should be sent to the repository administrator: staff-oatao@listes-diff.inp-toulouse.fr

  (2013) In: IFAC Conference on Manufacturing, Modeling and Control, 19 June 2013 -21 June 2013 (Saint Petersburg, Russian Federation).


Martin Durand
email: martin.durand@lip6.fr

Fanny Pascual
email: fanny.pascual@lip6.fr

Efficiency and Equity in the Multiple Organization Scheduling Problem

Keywords: multi agents scheduling, fairness, makespan minimization

Introduction

The Multi Organization Scheduling Problem (MOSP) (Pascual et al. 2007) deals a set of n organizations {O 1 , . . . , O n } which each owns both a set of identical parallel machines, and a set of sequential tasks to execute. The objective is to minimize the completion time of the last task completed on the machines shared by the organizations (the makespan), under an additional constraint: no organization should increase the last completion time of its tasks in the shared system, compared to the case where it executes its own tasks on its own machines. This last constraint is called the rationality constraint, and ensures that all the organizations have incentive to share their machines. More formally, let us denote by C i loc is the makespan of Organization O i if it schedules its own tasks on its owns machines -this scheduled is assumed to be given by the organization (it can minimizes the makespan of Organization O i ,or not) -and is called the local makespan of O i . Given any schedule S of all the tasks on all the machines, we denote by C i (S) the makespan of O i , i.e. the maximum completion time of a task of O i in S. Our problem is the following one: minimize C max (S) such that, for each i ∈ {1, . . . , n}, C i (S) ≤ C i loc .

Interest of cooperation. Let us first show that sharing machines does not only allow organizations which have many tasks and few machines to decrease theirs makespans by given tasks to the other organizations, but that each organization may decrease its makespan. Let us consider the following instance, in which all organizations can benefit from cooperating. There are n = 3 organizations, each one having only one machine. All the tasks are of length 1. Organization O 1 owns 3 tasks, O 2 owns 6 tasks, and O 3 owns 12 tasks. On figure 1 we can see on the left the local schedules (schedules in which each organization schedules its own tasks on its owns machines), and on the right a schedule in which the organizations share their machines. In this example, sharing the machines allow each organization to decrease its makespan. This example can be extended for higher values of n. We have shown that the best improvement which can be obtained for each organization simultaneously is a factor n, and that this ratio can be obtained on some instances.

Focus of this paper and map of the paper. Besides analyzing the best possible benefit that organizations can mutually have by sharing their machines, our aim is to focus on the efficiency of algorithms (where the efficiency is thought in term of makespan -the date at which all the tasks have been computed), and on the equity of algorithms for MOSP (it is not suitable that, even if the returned schedule fulfills the rationality constraint, the machines which are free are used only for the tasks of a single organization while some tasks of the other organizations are waiting). These two aspects may be antagonist, and our aim is to see to which extent, since what we want would be a schedule with a small makespan and in which machines are shared with equity.

In Section 2, we focus on efficiency: we analyze the highest possible increase of a local makespan which may occur if we want a schedule which minimizes the (global) makespan.

We then look at the problem where each organization agrees to increases its makespan (compared to its local makespan) by a factor (1 + ).

In Section 3, we focus on fairness: we introduce a new problem which consists in maximizing the minimal gain (decrease of its makespan) of an organization. Before presenting these results, we start by reviewing existing work on MOSP. (Pascual et al. 2007) has been introduced with parallel rigid tasks (tasks that need to be executed in parallel on several machines) and has mainly been studied from an approximation viewpoint. The best approximate algorithm is a 3-approximation algorithm when the organizations schedule locally the tasks in decreasing order of their heights (the height of a task is the number of machines needed to execute the task), or a 4-approximation algorithm in the general case (Dutot et. al. 2011). For sequential tasks (tasks that need to be executed on one machine only), the best known algorithm is a 2-approximate algorithm [START_REF] Cohen | Analysis of Multi-Organization Scheduling Algorithms[END_REF]) (in the sequel, all the papers -as well as our results -deal with sequential tasks).

State of art. The Multi Organization Scheduling Problem

Some papers also consider a relaxed version of MOSP: it is assumed that the organizations tolerate a bounded degradation on the makespan of their own tasks, and the aim is to minimize the global makespan. This problem is denoted by (1 + α)-MOSP [START_REF] Ooshita | A Generalized Multi-Organization Scheduling on Unrelated Parallel Machines[END_REF] when it is assumed that each organization accepts to increase the maximum completion time of its tasks by a factor at most (1+α). A 3 2 -approximate algorithm for 2-MOSP has been given (Cordeiro et. al. 2011). The closest work in spirit to what we will do in Section 2 is a study of (1 + α)-MOSP on unrelated machines [START_REF] Ooshita | A Generalized Multi-Organization Scheduling on Unrelated Parallel Machines[END_REF][START_REF] Ooshita | The Price of Multi-Organization Constraint in Unrelated Parallel Machine Scheduling[END_REF]. In this setting, Ooshita et al. show that, when there is no cooperation (α = 0), the makespan can be m times higher than in the optimal makespan without the rationality constraint. When α > 0, the authors also give a (2 + 2 α )-approximate algorithm for (1 + α)-MOSP.

Efficiency vs. increase of the local makespans

In this section, we study how the aim of minimizing the makespan is in opposition with the rationality constraint.

2.1 Necessary trade-off between the (global) makespan and the increase of local makespans.

We have shown that in an optimal schedule for the makespan minimization, an organization may increase its makespan up to a factor m, where m is the number of machines (due to lack of space the proof is omitted). This value, that could be called "the price of efficiency", is high. We will now assume that organizations may accept to increase their makespans in order to get an efficient schedule, but only if this does not increase to much their makespans. We will now assume that each organization agrees to increase a little bit its makespan: it will accept a schedule in which its makespan is increased by a factor at most (1 + ) compared to its local makespan.

Let ≥ 0. We assume that each organization O i agrees to have a makespan at most equal to (1 + )C i loc . If = 0, this is the MOSP. Otherwise, each organization agrees to increase a little bit its makespan (the higher α is, the higher an organization agrees to increase its makespan). We call (1 + )-MOSP, the problem where we wish to minimize the makespan with these relaxed constraints:

minimize C max (S) such that, for each i ∈ {1, . . . , N }, C i (S) ≤ (1 + )C i loc .
Thanks to a specific instance, we give a lower bound on the approximation ratio of any algorithm for (1 + )-MOSP with respect to the optimal makespan when there is no rationality constraints: this shows what we loose, in term of makespan, due to the relaxed rationality constraint. When = 0, we obtain the following proposition.

Proposition 1 There is no algorithm which returns schedules which fulfill the rationality constraint, and which is less than 2-approximate with respect to the global makespan.

This bound improves the previous one, 3 2 , which had been given (Pascual et al. 2007) first for two organizations and then [START_REF] Cohen | Analysis of Multi-Organization Scheduling Algorithms[END_REF] for more than two organizations. Furthermore, in [START_REF] Cohen | Multi-organization scheduling approximation algorithms[END_REF] the authors show that no approximation algorithm for MOSP has a ratio asymptotically better than 2 w.r.t. the global makespan (when m tends towards the infinity) when we add the constraint that on the returned schedule, each machine schedules the tasks of its organization (if any) before the tasks of other organizations. This constraint is thus not necessary to obtain the asymptotic ratio of 2.

A PTAS for the makespan minimization with a bounded increase on the local makespan

We adapt the polynomial approximation scheme (PTAS) presented [START_REF] Hall | Approximation Schemes for Constrained Scheduling Problems[END_REF] for a scheduling problem (makespan minimization with delivery times), to get a PTAS with resource augmentation for our problem. More precisely: given a fixed > 0, and a fixed number of organizations n, we will get a polynomial time algorithm which returns a schedule with a makespan at most (1 + ) times the optimal makespan, and in which the makespan of each organization is at most (1 + ) times its local makespan. The rationality constraint may thus be violated, but the increase of the makespans of the organizations is bounded, and may be acceptable if is small.

In the previous sections, we have assumed either that the rationality constraint should be fulfilled (but we then had as only objective function to minimize the makespan, and the gains for the organizations -the decrease of their makespans -in the returned schedule could be very different), or we have even assumed than we can relax (in a bounded way) the rationality constraint to get a schedule with an even smaller makespan. In the following section, we focus on fairness issues: we will keep the rationality constraint, and our focus will not be to decrease the makespan, but to get schedule in which all the organizations decrease their makespans by a factor as large as possible.

Maximizing the Minimal Decrease of the Local Makespans

Given a schedule S , the gain g i (S) of Organization O i represents how much Organization O i has decreased its makespan in S in comparison to its local schedule:

g i (S) = C i loc C i (S) .
The Maximal Minimal Gain problem, denoted as MaxMinGain, takes the same input as MOSP. It builds a schedule of all the tasks of all the organizations on the m machines of the organizations, in order to maximize the minimum gain among the organizations. The returned schedule is thus S = arg max S min i∈{1,...,N } g i (S) Problem MaxMinGain can be solved in polynomial time when all the tasks have the same length. Moreover, in this case, it is possible to find a schedule S which is optimal for MaxMinGain and which is optimal for problem (P ||C max ): the global makespan is minimized while the minimal gain of an organization is maximized. The algorithm is very simple: it consists in scheduling the tasks greedily, by increasing local makespans.

When tasks can have different lengths, MaxMinGain is strongly NP-hard and hard to approximate. We also show that there is no algorithm which is optimal for MaxMinGain and which has an approximation smaller than 2 for MOSP. Naturally, this implies that no algorithm can be optimal for MaxMinGain and have an approximation ratio smaller than 2 for (P ||C max ).

As seen in the previous section, a list scheduling by increasing local makespan is optimal for MaxMinGain when tasks all have the same length. However, in the general case, such a schedule can break the rationality constraint. We complete our results with a heuristic that aims at returning a schedule as close as possible from a list schedule by increasing local makespan but respecting the rationality constraint. On tested instances, this heuristic returns a schedule with an average makespan below 1.07 times the optimal makespan and an average minimum gain above 0.92 times the optimal one.

Fig. 1 .

 1 Fig. 1. An example in which each organization benefits from cooperating


A Audouin

P Ern

F Charru

Droplet entrainment from a single wave propagating in a stratified air-water pipe flow

We investigate the distribution of sizes and velocities of droplets detached from a single roll wave propagating in a pipe. The wave is induced by a pulse in liquid flowrate in a stratified air-water flow. Due to the shear applied on the crest of the wave by the faster gas stream, droplets are formed and entrained by the gas. Detection and tracking of the droplets is performed by image processing. The distributions of sizes and velocities of the droplets are determined in the center of the pipe for various distances to the roll wave and compared to those measured for higher gas flowrate when droplets are produced from all the roll waves formed naturally at the gasliquid interface without external forcing.

Introduction

In this paper, we investigate the atomization and redeposition of droplets in a horizontal stratified airwater flow in a pipe. Due to the shear applied on the liquid layer by the faster gas flow, waves appear at the interface and propagate on the liquid surface. For sufficiently high velocities of the gas, large-amplitude waves arise and droplets detach from the crest of these waves. It has recently been shown that primary atomization mainly occurs through two mechanisms: bag-and ligament-breakup, as illustrated in Figure 1. Since the characteristic sizes of the ligaments and bags vary from a few millimeters to a few centimeters, the range of sizes of the droplets produced is large [START_REF] Boulesteix | Size and velocity distributions of droplets in an air-water horizontal pipe flow[END_REF]. Moreover, it turns out that the droplets size distributions are more peaked and have heavier tails than the lognormal distribution usually used to account for these distributions. Notably, the distribution of droplets sizes presents a heavy tail at large drop diameters (> 0.5 mm), which gives the predominant contribution to the mass of entrained droplets. Analysis of the joint distributions of sizes and horizontal velocities also revealed that the most numerous droplets (diameters about 0.1 mm) travel at a velocity close to the gas speed, while larger drops travel at much lower velocities. Furthermore, the associated distribution of vertical velocities is non-Gaussian (displaying exponential power tail-off), large-size droplets being responsible for the larger weight of positive vertical velocities. However, most droplets present negative vertical velocities. As discussed in [START_REF] Boulesteix | Size and velocity distributions of droplets in an air-water horizontal pipe flow[END_REF], the following sketch then emerged concerning the contributions to the dispersed phase provided by the various roll waves upstream the volume of measure: drops that are still accelerating and in particular larger drops that are still not redeposited are provided by roll waves closer to the measurement volume, while smaller drops also originate from further away waves, since they rapidly reach their terminal velocity and are carried by the gas stream over longer distances before they redeposit. A better understanding of the effect of acceleration and transport by the gas and the effect of sedimentation by gravity on the droplets detaching from a roll wave is therefore needed. The aim of this work is thus to focus on the contribution provided by a single wave and to determine the size and velocity distributions of the droplets according to the distance between the measurement volume and the isolated wave.

Experimental set-up

The experiments were performed for an air-water flow in a horizontal pipe of diameter 5 cm and length 5 m. The configuration considered here corresponds to the superficial gas velocity U gs = 10 m/s and the superficial liquid velocity U ls = 0.07 m/s. The roll wave is induced by increasing the liquid flowrate by a factor 2.4 during about 1s. The passage of the wave in a tube section is detected with a differential pressure transducer measuring the pressure loss between two sections of the pipe located 3.75 m-apart and is used as trigger for image recording. The characteristics of the droplets produced from the wave were investigated in the center of the pipe for two distances from the volume of observation to the wave. The shortest one d s is 35 cm and the largest one d l is 80 cm. We used a high-speed camera (5000 fps) and a 20ns flash lamp to get motion-blur-free images (especially for the smaller droplets). The field of view was set to a length of 12.47 mm, a height of 6.50 mm, along with a depth of field 1.9 mm and a resolution of 9.28 µm/pix. This enabled us to detect droplets as small as 30 µm and as large as 2 mm in diameter. Image calibration and processing were used to discriminate the droplets present in the volume of observation from those which were out-of-focus thanks to a gradient-based method [START_REF] Marmottant | On spray formation[END_REF] and to retrieve their equivalent diameter from their area in the image. Together with the droplet detection program, a tracking routine employing a multi-hypothesis method [START_REF] Reid | An algorithm for tracking multiple targets[END_REF] was used to follow the droplets along a sequence of images and to determine their vertical and horizontal velocities. All programs were implemented using Matlab and its Image Processing Toolbox.

Celerity of the roll wave

It is generally assumed that roll waves propagate on the liquid surface with a constant celerity unlike the capillary-gravity ripples [START_REF] Belt | On the liquid film in inclined annular flow[END_REF]. To confirm this, we recorded the passage of sixty isolated roll waves at two differents locations of the pipe. We used two high-speed cameras (500 fps) situated 1.1 m apart to record the position of the air-liquid interface using a field of view of length about 58 mm and height 46 mm. Image processing then provided an estimation of the liquid layer thickness as well as the celerity of the wave by crosscorrelating the temporal signals of the liquid-layer thickness obtained at two points (located 5 cm apart) of an image of a given camera and between two points of each camera. The celerity of the roll wave U w is then defined by the ratio between the distance separating the two points of reference used in the cross-correlation and the time delay corresponding to its maximum value. For a given roll wave, three different estimations were thus obtained, yielding the celerity U w =1.5 ± 0.25 m/s. No acceleration of the roll wave was observed along the pipe, but the differences in the measured values indicate that the wave celerity fluctuates as the wave propagates in the pipe, as do the structure and shape of the wave.

Size distributions

The probability density functions (pdf) of droplets sizes and velocities were obtained for a total number of 16.000 droplets for each distance to the wave. In order to compare the different results, all the pdf are normalized, so that their integral is equal to 1. For both distances to the waves, the measured distributions are large, with droplets diameters ranging from about 40 µm to 3 mm with a predominance of small droplets. Nevertheless, the contribution to the liquid volume in the dispersed phase provided by larger droplets is considerable. It also appears that larger drops disappear when increasing the distance between the wave and the measurement volume as a consequence of deposition while the amount of small drops is kept almost constant (figure 2). The size and velocity distributions were investigated in a previous work [START_REF] Boulesteix | Size and velocity distributions of droplets in an air-water horizontal pipe flow[END_REF] for a wavy-stratified regime characterized by nearly the same superficial liquid velocity U ls = 0.067 m/s and a higher superficial gas velocity U gs = 16.9 m/s. In this regime, roll waves are formed naturally at the gas-liquid interface without external forcing of the liquid flow rate. Comparison of the distributions measured in the present work to those obtained in the wavy-stratified regime with natural atomization indicates that the relative weight in the distributions of smaller drops and of larger drops is stronger in the latter regime, which supports the idea of an accumulation in the dispersed phase of smaller droplets originating from a large range of roll waves.

For the wavy-stratified flow and for the isolated wave, figure 3 shows the comparison between the experimental distributions of droplets sizes and the Fréchet law (Fr) and the lognormal law (LN), which is commonly used in this context, e.g. [START_REF] Azzopardi | Drops in annular two-phase flow[END_REF]. The expressions of these laws are

€ ν Fr (d) = kd -k -1 d -k e d d $ % & ' ( ) -k and € ν LN (d) = 1 σ d 2π e - 1 2σ

Horizontal and vertical velocities distributions

We now consider the distributions of horizontal and vertical velocities associated to the size distributions of the droplets. Figure 4 presents the distributions of horizontal velocities V x of the droplets for both distances to the isolated wave and for the wavy-stratified regime investigated in [START_REF] Boulesteix | Size and velocity distributions of droplets in an air-water horizontal pipe flow[END_REF]. The distributions display a similar shape, with a long tail for lower horizontal velocities. It also appears that the pdf for the largest distance to the wave is slightly shifted towards the higher velocities. This result is consistent with an acceleration of the droplets by the gas stream. As shown in figure 5, two different contributions can be distinguished in the distributions. The first is a gaussian contribution which includes smaller droplets whose horizontal velocities are distributed around the mean gas velocity U g (determined by considering the mean liquid-layer thickness) because of the turbulent velocity fluctuations of the gas and of the variations in liquid layer thickness. The second component of the distribution which corresponds mainly to the lower velocities includes both large and small droplets that are still accelerating. This result is more clearly seen in Figure 6, which displays the distributions of droplets diameters by classes of horizontal velocities (for V x < 10 m/s, 10 < V x < 14 m/s, 14 < V x < 25 m/s): only small droplets reach high horizontal velocities (14 < V x < 25 m/s), while droplets with low horizontal velocities (V x < 10 m/s) present a broad range of diameters. Note that the lognormal law also suits to these distributions.

The distributions of vertical velocities V y of the droplets display a non-gaussian behavior (figure 7). Owing to the gravity, most of the droplets are settling (around 65%) and in particular the largest ones. Comparison of the distributions measured in this work to those obtained for the wavy-stratified with natural atomization (with U ls = 0.067 m/s and U gs = 16.9 m/s) indicates that the peak of the distributions as well as the mean vertical velocity are negative, with a more negative peak for the wavy-stratified flow with natural atomization. Nevertheless, the distributions are asymmetric with heavier tails for the positive velocities. As can be seen in Figure 8, most droplets with horizontal velocities close to the gas velocity have a negative vertical velocity while droplets with positive vertical velocities have lower horizontal velocities and can be considered to be still under the inertia of their ejection from the liquid layer.

Conclusions

Digital image processing was used to investigate the probability density functions of the sizes and velocities of the droplets detached from an isolated roll wave in a horizontal stratified air-water flow. For the two distances to the wave that were investigated, the size and velocity distributions of the droplets are large. However, we observed that larger drops are present in the vicinity of the isolated wave and disappear when increasing the distance between the measurement volume and the wave, as a consequence of redeposition. On the other hand, the amount of small drops is kept almost constant for both distances to the wave. Furthermore, it turns out that the droplets size distribution can be satisfactorily fitted by a lognormal law in both cases. Comparison of the distributions measured in this work to those obtained for a wavy-stratified regime with natural atomization [START_REF] Boulesteix | Size and velocity distributions of droplets in an air-water horizontal pipe flow[END_REF] indicates that the relative weight in these distributions of smaller drops and of larger drops is stronger in the latter regime, which supports the idea of an accumulation in the dispersed phase of smaller droplets originating from a large range of roll waves. Analysis of the distribution of horizontal velocities also revealed that the droplets have a lower speed for the smallest distance to the wave, confirming the acceleration effect of the gas on the droplets. Our goal is now to obtain a better description of the atomization and deposition distributions that could be used in the closure models for the momentum conservation equations in the liquid and gas phases.

Figure 1 .

 1 Figure 1. Example of bag-breakup at the crest of a roll wave.

Figure 2 .

 2 Figure 2. Distributions of droplets diameters for two distances from an isolated wave (with U ls = 0.07 m/s and U gs = 10 m/s) and for a wavy-stratified regime with natural atomization (with U ls = 0.07 m/s and U gs = 16.9 m/s).

d

  the mean diameter of the droplets, € σ the standard deviation, and k an ajustable parameter. It turns out that the droplets sizes distributions for the wavy-stratified flow are clearly better fitted by the Fréchet law (figure3.b), as shown in[START_REF] Boulesteix | Size and velocity distributions of droplets in an air-water horizontal pipe flow[END_REF], whereas the lognormal law is more appropriate in the case of the isolated wave it for both distances d s and d l (figure3.a).

Figure 3 .

 3 Figure 3. Comparison between the experimental pdf of droplets diameters (dotted blue) and a lognormal fit (dashed red) and a Fréchet fit (solid black), for (a) the isolated wave at the distance d s and (b) the wavy-stratified regime (with U ls = 0.067 m/s and U gs = 16.9 m/s).

Figure 4 .

 4 Figure 4. Horizontal velocity distribution of droplets for both distances to the isolated wave (U ls = 0.07 m/s and U gs =10 m/s) and for the wavy-stratified regime (with U ls = 0.067 m/s and U gs = 16.9 m/s).

Figure 5 .

 5 Figure 5. Horizontal velocity distribution for the short distance to the wave d s compared to a gaussian fit (dashed red line); the dotted black line is the residual component; (U ls = 0.07 m/s and U gs = 10 m/s).

Figure 6 .

 6 Figure 6. Distributions of droplets diameters by classes of horizontal velocities, for the short distance d s and for U ls = 0.07 m/s and U gs = 11 m/s.

Figure 7 .

 7 Figure 7. Vertical velocity distributions of droplets for both distances to the isolated wave (U ls = 0.07 m/s and U gs = 11 m/s) and for the wavy-stratified regime (U ls = 0.067 m/s and U gs = 16.9 m/s).

Acknowledgements

This work was supported by TOTAL Exploration & Production and Scientific Development Division. We would like to thank in particular D. Larrey, F. Luck and J.M. Munoz for their involvement in the project. We also thank S. Cazin, G. Ehses, H. Ayroles and M. Marchal from IMFT for the technical support. Figure 8. Distributions of droplets horizontal velocities for positive (dashed red) and negative (solid blue) vertical velocities, for the short distance d s , U ls = 0.07 m/s and U gs = 11 m/s.


Gianluca Lavalle

Claire Laurent

François Charru

Jean Paul Vila

Study of the dynamics of thin liquid films sheared by gas flows

Keywords: Saint-Venant equations, thin film, gas flow, stability analysis

We investigate the long-wave instability of a thin liquid film sheared by a concurrent gas flow. The film flows down an inclined plane and is driven by a constant shear stress and a pressure gradient. The height of the film is small compared to the wavelength, which justifies a long-wave approach. We adopt this procedure to approximate the Saint-Venant (or shallow water) equations. A linear stability analysis is thus developed and the limit between stable and unstable conditions is analyzed in terms of a critical Reynolds number. It comes out that the constant shear stress, applied at the free surface of the film, tends to stabilize the liquid layer to long waves. This is true when the mass flow of the film is kept unaltered. On the other side, the role of the pressure gradient is the same as the longitudinal component of gravity, which is known to be a destabilizing factor for the liquid film.

Introduction

Instabilities of thin liquid films have been studied for many years, since the pioneering works of Kapitzas. When a film flows down an inclined plane a wide range of wave evolution can appear at the free surface. The development of long waves manifests itself with the inception of linear structures and the transition to non-linear waves, such as solitons and periodic waves. When a gas flows above the film, its dynamics involves more phenomena, and a strong influence between the two phases rises at the interface. Both natural phenomena and industrial applications are related with the film dynamics. Particularly, in the domains of aeronautics and space, there exists many typical configurations, such as: water ingestion in aircraft engines, primary atomization in fuel combustion, and de-icing of aircraft systems. In order to model the latter, for example, it is necessary considering that the flow around a wing profile manifests strong pressure gradient at the leading edge and relevant shear stresses all along the profile. Therefore, a deep comprehension of the physics of long-wave instabilities of thin liquid layers is still required. When the height of waves at the interface remains smaller than their wavelength, an asymptotic approach based on long wave approximation results to be more suitable to model the film dynamics, due to its low-dimensional construction. This procedure was firstly proposed by Benney [Benney]. Furthermore, the use of Saint-Venant equations, integral form of the full Navier-Stokes, has certain benefits, especially for low Reynolds numbers [Tseluiko and Kalliadasis]. Shkadov firstly proposed a model of Saint-Venant equations where all the variables are enslaved to the film thickness and average velocity [Shkadov]. The limitation of his model is the lack of consistency, because the disagreement of the dispersion relation with the Orr-Sommerfeld theory. Later, Ruyer-Quil and Manneville have improved Shkadov's model, expanding the velocity field by using polynomial functions [Ruyer-Quil and Manneville 1/2]. In this way, they were able to ensure the consistency of their Saint-Venant equations, that is, the linearized model provides the right dispersion relation and the correct critical Reynolds number , above which the liquid film is unstable. Subsequently, Boutounet et al. have proposed a frame to obtain several models at given order accurate [Boutounet et al.]. However, in order to develop such first-order models, it is necessary to calculate the whole velocity field, that could be very tricky. Luchini and Charru have thus presented a way to avoid the calculation of the first-order velocity field [Luchini and Charru 1/2]. Their development is based on the use of the energy equation rather than the momentum equation, in order to get the proper reductions to prevent the hard calculation of the velocity correction. All these already cited Saint-Venant models are consistent and expressed in a non-conservative form. In this paper we present a consistent system of Saint-Venant equations, combined with a Benney-like asymptotic development, which satisfies two more properties, namely it is expressed in a conservative form and the concerning numerical flux does not depend on the possible bottom velocity. The former is relevant in numerical developments, i.e., a conservative numerical flux provides benefits for finite volume method. The latter, instead, is significant when computing industrial calculations involving rotating components, such as turbines and fans. In these configurations, the bottom velocity, represented by the blade speed, reaches very high values; if it is contained into the numerical flux, besides in the source terms, this velocity can be greatly high compared to the other terms of the expression. 1982).

Later, Smith gave a critical condition for the instability of such a film [Smith 1990]. In his configuration the film thickness is supposed to be fixed. Hereby, we provide a critical Reynolds number depending on a prescribed shear stress, when the mass flow of the film is kept unchanged. Our result is different from the cited one, in a way that the effect of the shear stress is to stabilize the film to large wavelength disturbances. Furthermore, the effect of the pressure gradient in the stability condition of the liquid layer has been analyzed. Its role is the same as the longitudinal component of gravity, thus destabilizes the film interface. where and are the streamwise ( ) and cross-stream ( ) velocity components, the pressure, g the gravity, the angle of the inclined plane and the kinematic viscosity. These equations must be completed with boundary conditions at the wall and at the free surface of the film, where and respectively. The no-slip condition reads

00 | 0 , | 0 uv (2)
At the free surface, boundary conditions state the continuity of y-velocity and the balance of tangential and normal stress components. The first of such conditions sets the free surface a material line, and reads

| | . t h x h h u h v (3)

With reference to Figure 1, the two other conditions yield

() • ( , ) n e p x t nT (4) () • ( , ) n xt tT (5)
where is the normal stress vector, with stress tensor, the gas pressure and shear stress of the gas at the free surface. By the development of conditions (4), ( 5) one gets for normal and tangential stresses respectively Hence, the equations (1) with boundary conditions (2), ( 3) and ( 6), ( 7) describe the whole film dynamics. where is the dimensionless average film velocity. Note that the y-momentum equation has been already replaced into the equations ( 13). These are expressed in terms of the unknowns and film thickness and flow rate respectively, except for the momentum transport and wall shear stress terms. In order to enslave all the terms of the equation to and only, two closure models are required. This consists in assuming a proper shape for the velocity profile.

2 2 2 ( | | ) | | 1 |0 () x y h x h x h x y h x he h u v u h v h pp (6) 2 2 1 2 ( | | ) (1 )( | | ) 1 ( , ) 0 () x y h x h x y h x h x h v u h u v h xt (7)
( | | ) | | ( | )(1 ) 0 cos 1 2 ( | | ) 1 11 (1 ) | | 0 () ( ( )) t h x h x y h x h y h x x h h e x x y h x h x x y h x h h u h v h u v v h u Re p p h Fr h v u h h u v (11)

The base state

A trivial solution of system (9-11), in line with the wavyless solution of Nusselt, is based on the assumption that the film thickness is steady and constant, in a way that both the derivatives with respect to and are zero. By developing the calculations, the equilibrium velocity and pressure profiles read

(0) (1) (0) (1) (0) (1) u u u v v v p p p ( 17
)
where zeroth-order development coincides with the base state solution. According to the dimensionless pressure (8), the basic state consists in a constant pressure inside the liquid layer. The classic hydrostatic behavior results just shifted to the next order. The first-order variables can be found by substitution and then development of ( 17) into (9-11). Once the quantities (17) are known until the first order, the closure models for the momentum transport and wall shear stress terms can be obtained. (23) Unlike the previously cited works, this system of consistent shallow water equations is expressed in conservative form. This turns to have certain advantages when using a finite volume method. Furthermore, if we introduce a bottom velocity in the analysis, the numerical flux of equations ( 23) remains unaltered. The bottom velocity will appear as a source term only. When solving industrial calculations with rotating components, typical of an aerospace environment, this properties results to be greatly relevant for numerical stability. .

Linear stability analysis

Let analyze the linear stability associated to the system of shallow water equations defined previously. Let thus calculate the dispersion relation and then the conditions setting a stable or unstable liquid layer. We admit plane wave solutions, such that both film thickness and flow rate can be written in the form 0 0

ˆexp[ ( )] , ˆexp[ ( )]

h h h i kx t q q q i kx t

(24)
where the subscript refers to the base equilibrium solution defined before, and and are the wavenumbers and the angular velocity respectively. We admit real wavenumbers and complex , so that the angular velocity , where is the complex wave velocity, can be written as

r i r i i kc ikc (25)
The sign of the amplification or damping rate defines the conditions for a stable or unstable film. When is negative the solution is stable to large wavelengths, whereas represents the transition and thus defines the neutral-stability curve. Inversely, when is positive the film is unstable to large wavelength disturbances. Let replace the wavy solutions ( 24) into the quasi-linear form of the SW equations ( 23

q qi qik hik h ikh hh h ikh h ikh h ikh Fr q h q h Re Reh Re h (28)
The system of equations ( 26) and ( 28) allows to find the dispersion relation, which governs the relationship between and , and subsequently the neutral-stability curve. Let replace the equation ( 26) into (28). With the help of ( 27), the dispersion relation reads

() i h ik h ik Reh k Re h h h h ik Fr (29) )
Note that the assumptions of and equal to unity permits to semplify the calculations, without loss of generality. By an asymptotic expansion of in powers of the wavenumber (thus, in the limit , which corresponds to large wavelengths), the solution of (29) can be written as (we express here the wave velocity ) 22 2 2 1 () 3 5 5 ()

Re c i k O k Fr (30)
Note that the wave velocity is complex. As stated before, the real part represents the propagation velocity, and the imaginary part the growth rate.

Different cases have been studied. When only gravity acts on the liquid layer, such as a falling film along an inclined plane, the wave velocity reduces to

If

, or , the film is stable to large wavelength disturbances. It is manifest that a horizontal film is always stable, whereas a vertical film is always unstable. When a constant shear stress supports the film, such as the case of a falling film sheared by a concurrent gas flow, the wave velocity changes both in the real and imaginary part. Now we obtain . Hence, the propagation velocity of the waves, related to the film velocity, decreases in respect to the falling film configuration, and this happens because the shear stress increases the film velocity more than the wave speed. The stability condition now reads 5 / 6cot 5 /18 , 1 / 6

(1 / 2)(1 / 6) Therefore, the value of , as for , increases when a constant shearing is applied in the direction of the film flow. Its effect is thus to stabilize the film to large wavelength disturbances.

Let rewrite now the expression of using the equilibrium velocity. It yields For any , the value of is thus confined in the range . This makes the critical Reynolds number lying in the range , as sketched in Figure 2. Therefore, although increases, it is limited up to . The expression of as function of the dimensionless shear stress depends on the definition of the Reynolds number. In our analysis it depends on the average velocity of the liquid layer. Apart from the reference velocity, the results of Smith [Smith 1990] differs from this one in the definition of the Reynolds number. In his work, it does not depend on the shear stress. Hence, his results shows a destabilizing effect of , but when the film thickness is kept unchanged in the experience. Another way to quantify the effect of the shear stress over the free surface of the film is to redefine the dimensionless setting of and introduce a dimensionless film thickness , By expressing both the Reynolds number and its critical value in terms of , a stability condition for the thickness of the film can be found, which reads

In this way, we can easily compare this case to the configuration of falling films. The analysis of provides the relationship 6

Conclusions

Combined with a numerical flux not depending on the possible bottom velocity, the conservative form of Saint-Venant equations expressed above turns to be very useful in many numerical experiences. Among the others, this model results to be the only consistent one to possess these two properties.

The results of the linear stability analysis for a liquid film sheared by a constant shear stress show that the effect of this stress consists in increasing the critical Reynolds number. By its definition, it is manifest that the presence of the shear stress stabilizes the film to large wavelength disturbances. This is not in agreement with the work of Smith. The difference comes out from the definition of the Reynolds number. In his work, it is defined as function of film thickness only. In our work, the Reynolds number depends on the average longitudinal velocity, which contains also the shear stress. In conclusion, his experience is obtained as the film thickness is kept fixed, whereas in our is the mass flow to be unchanged. This result is very important for academic experiments and numerical calculation, as also for industrial applications. Since the evident difficulties to express a unique stability condition for this experience, a critical condition of stability based on dimensionless thickness has been developed. The second main conclusion consists in the critical condition for stability when the film is driven by a pressure gradient too. The modification of the longitudinal component of gravity consists in a fictitious alteration of the inclination angle and the gravity norm . The expression of the critical Reynold number does not change its form, but the angle. Geometrical considerations drive to the conclusion that the decreases compared to the case without pressure gradient. Its role is thus to destabilize the film, in contrast with the shear stress.

  incompressible flow, such as a thin liquid film flowing down an inclined plane and sheared by a concurrent gas flow (Figure1), Navier-Stokes equations reduce to

Figure 1 :

 1 Figure 1: Inclined plane and boundary conditions.

  s development [Benney], let expand the velocity and pressure fields as

Figure 2 :

 2 Figure 2: Critical Reynolds number versus . .

  On the other hand, the momentum equation produces two contributions: one is the non-exponential part, which gives the equilibrium flow rate, and the other is the exponential term.

  One can thus summarize that a negative tends to destabilize the film, because the tangential component of gravity increases.

Figure 3 :

 3 Figure 3: Scheme of fictitious when adding to the longitudinal gravity component.

  Secondly, a linear stability analysis is developed for a thin liquid film falling on an inclined plane, with a constant shearing applied at the free surface. The instability of such a liquid layer has been developed byMiles [Miles 1960] and Benjamin [Benjamin 1959], then improved by Smith and Davis [Smith and Davis

	th International Conference on Multiphase Flow
	ICMF 2013, Jeju, Korea, May 26 -31, 2013
	2