Pierre Berthonnaud

Komla Domelevo

Existence of travelling waves and large activation energy limits for a onedimensional thermo-diffusive lean spray flame model

Keywords: reaction-diffusion equations, spray flames, travelling waves, large energy activation, singular limits

. The transition parameter is proved to be the ratio between the speed of vaporisation and the speed of combustion, and we are able to prove and estimate the thickening of the combustion zone in the case of slow vaporisation, which is again a novelty compared to purely gaseous combustion. The model at task is therefore on the one hand simple enough to allow for explicit asymptotic limits, and on the other hand rich enough to capture some particular aspects of spray combustion. Finally, we briefly consider situations where the vaporisation is infinitely fast, where the spray is polydisperse, or where the geometry of the problem is different from that of the onedimensional travelling wave.

Introduction

This paper provides a rigourous mathematical analysis of some aspects of spray combustion, including the analysis of the so-called large activation energy limit for a spray flame model. This notion of large activation energy limit was first introduced in the pioneering work of Zeldovich and Frank-Kamenetskii [START_REF] Ya | A theory of thermal propagation of flame[END_REF], and refers to the limit where the combustion rate is much faster than any other physical phenomenon, in particular diffusion. Gas-vapor-droplets systems have many applications in industry or everyday's life, such as diesel or propulsion engines. When trying to understand some specific features, it appears that the structure of these two-phase flames as well as their speed or stability are greatly affected by the presence of vaporising liquid droplets possibly interacting with the combustion zone.

The behaviour of spray flames has been investigated a lot in the physics literature and a wide variety of regimes were considered. Dating from the 70s and early 80s, we can quote the works of Polymeropoulos et al. [START_REF] Polymeropoulos | Flame propagation in a one-dimensional liquid fuel spray[END_REF] [START_REF] Polymeropoulos | The effect of droplet size on the burning velocity of a kerosenair spray[END_REF], Mizutani et al. [START_REF] Mizutani | Combustion of fuel vapor-drop-air systems: Part i -open burner flames[END_REF] [START_REF] Mizutani | Combustion of fuel vapor-drop-air systems: Part ii -spherical flames in a vessel[END_REF], Hayashi et al. [START_REF] Hayashi | Flame propagation in droplet-vapor-air mixtures[END_REF] [START_REF] Hayashi | Propagation velocity and structure of flames in droplet-vapor-air mixtures[END_REF]. Ballal et al. [START_REF] Ballal | Flame propagation in heterogeneous mixtures of fuel droplets, fuel vapour and air[END_REF]. They present studies of the propagation of liquid fuel sprays, study the influence of the size of the droplets, the type and geometry of the spray flame as well as its structure. The investigation of those vapor-drop-air systems was continued in the 90s. To quote but a few works, we refer to the work of Aggarwal and Sirignano [START_REF] Aggarwal | Unsteady spray flame propagation in a closed volume[END_REF] as well as the papers of Greenberg, Tambour and Silverman [START_REF] Silverman | Asymptotic analysis of a premixed polydisperse spray flame[END_REF] [START_REF] Silverman | Solutions of polydisperse spray sectional equations via a multiple-scale approach: an analysis of a premixed polydisperse spray flame[END_REF], where the structure of spray flames as well as the influence of parameters such as droplet size, fuel volatility, or equivalence ratio are also investigated analytically. More elaborate situations appearing in propulsion engines are for example pulsating or acoustic instabilities. We refer to [START_REF] Greenberg | Instability of a flame front propagating through a fuelrich droplet-vapour-air cloud[END_REF] [START_REF] Greenberg | Linear stability analysis of laminar premixed spray flames[END_REF] [START_REF] Greenberg | Stability boundaries of laminar premixed polydisperse spray flames[END_REF][10] [START_REF] Greenberg | Opposed flow polydisperse spray diffusion flames: steady state and extinction analysis[END_REF][7] and more recently [START_REF] Constantin | Existence of pulsating waves in a model of flames in sprays[END_REF] [START_REF] Greenberg | Stability of rich premixed spray flames[END_REF][17] [START_REF] Kats | Application of a non-asymptotic approach to prediction of the propagation of a flame through a fuel and/or oxidant droplet[END_REF] for studies in that direction.

Only few of these studies of spray flames (e.g. [START_REF] Constantin | Existence of pulsating waves in a model of flames in sprays[END_REF]) involve rigorous-in-the-mathematical-sense analysis of the model at task where existence, uniqueness, or asymptotic limits are derived. This is in contrast with purely gaseous combustion, where a lot of results exist in the mathematical literature for various regimes and asymptotics. In particular, a complete study of thermo-diffusive lean gaseous flame fronts and the large activation energy limit is in the paper of Berestycki, Nicolaenko and Scheurer [START_REF] Berestycki | Traveling wave solutions to combustion models and their singular limits[END_REF], and in the paper of Berestycki and Larrouturou [START_REF] Berestycki | Quelques aspects mathématiques de la propagation des flammes prémélangées[END_REF]. The originality of the present work is to provide a complete mathematical analysis of the counterpart of those systems, namely a lean spray flame model that on the hand is simple enough to allow for explicit asymptotic limits, and on the other hand rich enough to capture some particular aspects of spray combustion. We prove in particular the existence of two distinct combustion regimes in the large activation energy limit. We prove also that the speed of the spray flame decreases when the time for complete vaporisation of the droplets exceeds a critical value. Those results should be compared with the very interesting work of Suard, Nicoli and Haldenwang [START_REF] Suard | Vaporisation controlled regime of flames propagating in fuel-lean sprays[END_REF]. Using numerical simulations, these authors investigate the scaling laws of the spray flame velocity versus the Damkhoeler number (i.e. the ratio of the typical time for vaporisation versus the typical time for combustion). They also studied the structure of the reaction zone, of which they prove that the complexity increases in the so-called "vaporisation-controlled regime". Our mathematical analysis shows similar results, with however a sharp transition in the large activation energy limit. Also, an important result of the present work is that we can characterize completely the limiting profiles, with explicit expressions.

The model we consider is built as an elaboration of the purely gaseous thermo-diffusive lean flame (see [START_REF] Berestycki | Traveling wave solutions to combustion models and their singular limits[END_REF], [START_REF] Berestycki | Quelques aspects mathématiques de la propagation des flammes prémélangées[END_REF]) that originally involves only the temperature T of the mixture and the fraction Y of the gaseous fuel. We couple this thermo-diffusive model to vaporising droplets that act as a reservoir of liquid fuel. That is, the vaporisation of the liquid droplets produces (part of) the fuel vapor needed for the combustion. There is no combustion of the liquid phase. This leads to the following minimal model for a lean spray flame:

     -T ′′ + cT ′ = f (T )Y , -dY ′′ + cY ′ = -f (T )Y + 2πng(T ) R √ , cR ′ = -g(T )H(R), (1)
supplemented with conditions at infinity

   T (-∞) = 0, T (+∞) = T b , Y (-∞) = Y u , Y (+∞) = 0, R(-∞) = R u .
where T (x, t) is the nondimensional temperature, Y (x, t) the fraction of reactant, R(x, t) the square radius (i.e. 4 π R is the surface of the droplet), and H is the Heaviside function. Also, we have

T b := Y u + 4/3 π n u R u 3/2
. Finally, f (T ) is the reaction rate and g(T ) is the vaporisation rate. The reaction term is of ignition type, that is f is nonnegative, piecewise Lipschitz continuous, and there exists an ignition temperature θ i > 0 such that f (T ) ≡ 0 for T ≤ θ i and f (T ) > 0 for T > θ i . The function g satisfies exactly the same hypotheses but with a boiling temperature θ v < θ i , which means that the droplets start vaporising before the combustion starts. Finally, the constant d is the inverse of the Lewis number. At constrast with the classical thermo-diffusive system for gas flames, the system above in the case d = 1 does not reduce to a single scalar equation. Finally, notice that this model does not include the effect of the latent heat. The latent heat would appear as a loss term proportional to the vaporisation term in the equation for the temperature. Neglecting this phenomena is certainly an important simplification in the case where the latent heat is large. We therefore assume implicitely in thius study that this effect is negligible.

The plan of the paper is as follows. Section 2 is devoted to the proof of the existence of travelling waves for the system (1) above as stated in Theorem 1 for the monodisperse case, that is in the case where the droplets at a given position all have the same size. The extension to the case of polydisperse sprays is briefly discussed in section 4.2. Section 3 is devoted to the study of the high activation energy limit, where one assumes that combustion is an infinitely fast process. Suitable scalings for the reaction rate give rise to a limiting problem where combustion occurs in a possibly infinitely small region:

     -T ′′ + c T ′ =cT b δ x=0 , -d Y ′′ + c Y ′ =-cT b δ x=0 + 2 π n u g(T ) R √ , c R ′ =-g(T ) H(R),    T (-∞) = 0, T (∞) = T b Y (-∞) = Y u , Y (∞) = 0 R(-∞) = R u , (2)
As announced before, we will prove that this system obeys two different regimes depending on the size of the droplets (i.e. depending on the Damkhoeler number). With the notations above, when the size of the droplets is below a critical size R c , the droplets finish vaporising before the reaction front, and the velocity of the flame is a constant equal to the velocity of the equivalent gas flame. On the contrary, when the size of the droplets is larger than the critical value R c , the velocity of the flame starts decreasing. This is the so-called "vaporisation-limited" regime (see [START_REF] Suard | Vaporisation controlled regime of flames propagating in fuel-lean sprays[END_REF]). Moreover, by zooming on the reaction zone in the large activation energy limit, we are able to estimate the size of the vaporisation region that intersects the combustion region. These results are stated in Theorem 11.

Existence of Travelling Waves

We want to prove here Theorem 1. There exists at least one travelling wave

(T , Y , R, c) in X = C 1 (R) × C 1 (R) × C(R) × R solution to Problem ( 1) above .
The study of the existence of travelling waves for system (1) is done following the strategy of H. Berestycki, B. Nicolaenko and B. Scheurer [START_REF] Berestycki | Traveling wave solutions to combustion models and their singular limits[END_REF]. Namely : in Step 1, we state the problem in a bounded domain. In Step 2 we rewrite the problem as a fixed point problem involving an homotopy parameter. In Step 3, we prove a priori estimates independent of the size of the domain. In Step 4, we prove the existence of a travelling wave in a bounded domain. Finally, in Step 5, we let the size of the domain go to infinity. The major step is Step 3 whose proof is postponed to the next section. The proof of Step 5 being very similar to the classical gaseous flame case, we will not detail it here (see [START_REF] Berestycki | Traveling wave solutions to combustion models and their singular limits[END_REF]).

Steps 1 & 2: Problem in a bounded domain & homotopy technique

Let a be a fixed positive number. We consider the following boundary value problem on [-a, a]:

Problem 2. Find (T , Y , R, c) satisfying in ] -a, a[:          -T ′′ + c T ′ =τ [f (T ) Y ], -d Y ′′ + c Y ′ =τ -f (T ) Y + 2 π n g(T ) R √ , c R ′ =τ [-g(T ) H(R)], c =T (0) -θ i + c    -T ′ (-a) + c T (-a) = 0, T (a) = T b , -d Y ′ (-a) + c T (-a) = Y u , Y (a) = 0, R(-a) = R u , (3)
We will note R,c) in the Right-Hand-Side of the system, integrating and solving for (T , Y , R, c) in the Left-Hand-Side defines a fixed point problem that we denote by F τ ≡ I -K τ , so that proving the existence of a solution (T , Y , R, c) to system (3) is equivalent to finding (T , Y , R, c) such that F τ (T , Y , R, c) cancels. Thanks to the boundedness of the domain and the regularizing properties of the elliptic operator, one verifies easily that we have: Lemma 3. Assume that the functions f and g are continous and Lipshitz. The mapping F : R,c) is well-defined, compact and uniformly continous with respect to τ.

X a = C 1 ([ -a, a]) × C 1 ([ -a, a]) × C([ -a, a]) × R. Now, injecting functions (T , Y ,
X a × [0, 1] -→X a defined by (T , Y , R, c) × τ -→ (I -K τ )(T , Y ,
Step 3: Uniform a priori estimates The following result is proved in Section 2.1 below: Lemma 4. Let F τ be defined as above. There exists finite positive constants M , c, c ¯independent of a, such that, for Ω defined as:

Ω := {(T , Y , R, c) ∈ X ; T C 1 (Ia) M , Y C 1 (Ia) M , R C(Ia) M , c < c < c}
Step 4: Application of the Leray-Schauder Topologocal Degree. Assuming the uniform a priori estimates above, we can find positive constants M , c, c ¯such that the solution of the fixed point cannot be found on the frontier ∂ Ω. Also, using Lemma 3, we have:

F τ (∂Ω) = (I -K τ )(∂Ω) = / 0, ∀τ ∈ [0, 1] deg (F τ , Ω, 0) = deg (F 0 , Ω, 0) = / 0, ∀τ ∈ [0, 1]
Finally, one checks easily that the linear problem F 0 has a unique solution in X, so that deg (F 0 , Ω, 0) = / 0. Therefore we have also deg (F 1 , Ω, 0) = / 0 which implies the existence of a solution to the problem in the bounded domain [-a, a].

Remark 5. If the functions f and/or g have a discontinuity at θ v and/or at θ i , then we can prove the same result by using regularisation techniques in order to overcome the fact that K τ is no longer continuous (see [START_REF] Berestycki | Traveling wave solutions to combustion models and their singular limits[END_REF]).

Uniform a priori estimates. Proof of Lemma 4.

The aim of this section is the proof of Lemma 4. Without loss of generality, we can assume that τ = 1. This lemma can be divided into two propositions: Proposition 6. [Estimates on the solutions]For every solution (T , Y , R, c) of problem ( 3) with c ≥ 0, the following holds:

c > 0, 0 ≤ R ≤ R u , 0 < T ≤ T b , 0 ≤ Y < T b , R ′ ≤ 0, 0 < T ′ ≤ cT b , - c d T b ≤ Y ′ ≤ c d T b + 4 3 π n u R u 3/2 .
Proof. The estimates on the size R of the droplets are obvious. The estimates on the other unknowns follow the same lines as in [START_REF] Berestycki | Traveling wave solutions to combustion models and their singular limits[END_REF].

Proposition 7. [Estimates on the velocity]Let

h(s) = (T b -s) f (s), and assume that H(T b ) = θi Tb h(s) d s < ∞. Then, if (T , Y , R, c
) is any solution of ( 3), c satisfies:

min 2 dT b 2 H(T b ), c s 2 ≤ c 2 ≤ max   ln ( 2Tb θv ) a 0 , max 2 M , 2 MT b θ v   (4)
We prove successively lower and upper-bounds independent of the size of the domain for the velocity c:

Proof of the Lower-Bound of Proposition 7:

The simple but important alternative stated in Lemma 8 below allows us to take advantage of the estimates of Proposition 9 below:

Lemma 8. [An important alternative] If c ≤ c s := 1 R ui θv θi g(s) s d s , then R(x i ) ≤ 0 with x i satisfying T (x i ) = θ i . Proposition 9. Every solution (T , Y , R, c) of problem ( 3), with c ≥ 0, satisfy on [-a, a]: |T + Y -T b | ≤ | 1 -d d | (T b -T ) + 4 3 π n u R 3/2 (x) |T + d Y -T b | ≤ |1 -d| (T b -T ) + 4 3 π n u R 3/2 (x) (5)
and consequently

d > 1: 1 d (T b -T ) - 4 3 π n u R 3/2 (x) ≤ Y ≤ (T b -T ) + 4 3 π n u R 3/2 (x) 0 < d < 1: (T b -T ) - 4 3 π n u R 3/2 (x) ≤ Y ≤ 1 d (T b -T ) + 4 3 π n u R 3/2 (x) (6)
Proof of Lemma 8:

Proof. (of Lemma 8) First, we recall that T (x) = θ v e cx on [-a, x i ], so that T ′ (x) = c T (x). The equation for R is:

dR dx = dR dT dT dx = - 1 c g(T (x)). (7)
If we integrate dR dT between the temperatures θ v and θ i , we obtain:

R(x i ) -R u = - 1 c 2 θv θi g(s) s d s,
and, for c ≤ c s , the droplets vaporise before they reach the reaction zone, so that the combustion is monophasic.

Proof. (of Proposition 9 )

The proof follows the same lines as in [START_REF] Berestycki | Traveling wave solutions to combustion models and their singular limits[END_REF]. It relies on suitable integrations of the auxiliary functions W = T + Y and Z = T + d Y . We refer to [START_REF] Berestycki | Traveling wave solutions to combustion models and their singular limits[END_REF] as well as the next section for an exposition of this classical manipulations. Now, we can use the same procedure as in Berestycki, Nicolaenko and Scheurer ( [START_REF] Berestycki | Traveling wave solutions to combustion models and their singular limits[END_REF]) in order to prove lower-bounds for the velocity c.

Proof. (of the lower-bound of Proposition 7)

We give the details for the case d > 1 only, the proof for the case d < 1 being identical. First, if c ≥ c s , then the first inequality in ( 4) is proved. Consider now the case where c ≤ c s , which implies that the combustion . We start now with the identity obtained by integration of -T ′′ + c T ′ = f (T ) Y between 0 and a, after multiplication by T ′ :

- 1 2 |T ′ (a)| 2 + 1 2 c 2 θ v 2 + c 0 a |T ′ (s)| 2 d s = xi a f (T (s)) T ′ (s) Y (s) d s. (8)
But, c ≤ cs , so that R(x) = 0, ∀x ≥ x i , and we deduce:

Y (x) ≥ 1 d (T b -T (x)), ∀x ≥ x i .
We use this lower-bound for Y in the previous integral to obtain:

- 1 2 |T ′ (a)| 2 + 1 2 c 2 θ v 2 + c 0 a |T ′ (s)| 2 d s ≥ 1 d x i a f (T (s)) (T b -T (s)) T ′ (s) d s.
We make the change of variable t := T (s) in the integral and obtain:

- 1 2 |T ′ (a)| 2 + 1 2 c 2 θ v 2 + c 0 a |T ′ (s)| 2 d s ≥ 1 d H(T b ).
Then, similar L 2 estimates can be obtained by integrating against T the equation for the temperature. This allows to control the integral above and therefore to prove that f (s) and a 0 > 0 fixed: then, for each a ≥ a 0 , any solution (T , Y , R, c) of ( 3) satisfies:

2 d T b 2 H(T b ) ≤ c 2 ,
c ≤ max   ln ( 2Tb θv ) a 0 , max 2 M , 2 MT b θ v   . (9)
Proof. We compare the solution T of (3) to the solution T ¯of the following system:

-T ¯′′ + cT ¯′ = T b M H(x) o n] -a, a[ -T ¯′(-a) + cT ¯(-a) = 0, T ¯(a) = T b . ( 10
)
where H is the Heaviside function. By a direct calculation, we have:

     T ¯(x) = T (x), o n [-a, 0], T ¯(x) = T b e cx M c 2 (1 -e -ca ) + (1 - M c a) e -ca , o n[0, a]. (11)
By using maximum principle, we have T ≤ T ¯on [-a, a]. In particular, for x = 0, we have:

θ v ≤ T b e -ca (1 - M c ) + T b M c 2 (1 -e -ca ). (12)
We set c 0 = max (2 M , This completes the proof of the proposition.

Large Activation Energy Limit

We deal now with the most important part of the paper. We want to derive an asymptotic model for the spray flame corresponding to the so-called large activation energy limit, where one assumes that the typical time associated to the chemical reaction is infinitely small compared to any other time-scales, namely in our case time-scales associated to the diffusion process or that associated to the vaporisation process. The most important result is that we are able to prove and quantify in the limit the existence of two different regimes for the propagation of the spray flame, which is a striking and novel feature compared to classical gas flames. On the one hand, when the radius of the droplets in the fresh mixture is small enough, the velocity of the flame is that of the equivalent gas flame, that is the correponding gas flame where all the fuel is injected in gaseous form. On the converse, above a critical radius of the incoming droplets, the velocity of the flame starts decreasing and, following the terms of Suard, Nicoli and Haldenwang ( [START_REF] Suard | Vaporisation controlled regime of flames propagating in fuel-lean sprays[END_REF]), we are in the presence of a "Vaporisation Controlled" combustion regime. Moreover, the asymptotic analysis performed below provides details on the structure of the internal combustion layer in the presence of vaporising droplets. This are new mathematical results, consistent with the results and numerical observations made by the previous authors.

After stating the results, Section 3.0.1 is devoted to a heuristic proof of the limiting profile. It allows us to emphasize the existence of the two combustion regimes and to give more insight in the system. The rigourous proof together with the study of the internal layer are done in Section 3.1.

We start again with system (1) in which the reaction term is now set as:

f ε (T ) :=      0, T < θ i , 1 ε γ e x p T -T b ε , T > θ i . ( 13
)
and we are interested in the existence of travelling waves in the limit where ε goes to 0. Formally, one expects at the limit that the reaction term occurs in an infinitely small region, that we set at the origin

x = 0.      -T ′′ + c T ′ =cT b δ x=0 , -d Y ′′ + c Y ′ =-cT b δ x=0 + 2 π n u g(T ) R √ , c R ′ =-g(T ) H(R). (14)
We want to prove: [START_REF] Aggarwal | Unsteady spray flame propagation in a closed volume[END_REF] tend to a degenerate travelling wave in the sense that c ε tends to 0. * . If γ = 2 then the solutions (T ε , Y ε , R ε , c ε ) to system [START_REF] Aggarwal | Unsteady spray flame propagation in a closed volume[END_REF] tend to a travelling wave solution of the limiting system [START_REF] Greenberg | Instability of a flame front propagating through a fuelrich droplet-vapour-air cloud[END_REF] above, were the limiting velocity c has the expression

Theorem 11. Let M u := 4/3 π n u R u 3/2 fixed be given, such that T b = Y u + M u is kept constant when R u varies. Let γ be given in the intervall [0, 2]. We have that * . If γ < 2 then the solutions (T ε , Y ε , R ε , c ε ) to system
c = min   1 T b 2 d , 1 R u θv Tb g(s) s d s   ( 15
)
As a consequence, the behaviour of the velocity c of the limiting travelling waves shows a bifurcation when the radius R u of the incoming droplets is greater than the critical radius R c :

R c := d T b 2 2 θv Tb g(s) s d s (16)
* . Moreover, when γ = 2 and R u ≥ R c , then the size of the vaporising zone intersecting the reaction zone is at most of order O(ε 2/5 |ln ε|) for small values of ε > 0.

Heuristic asymptotic analysis of spray flames.

Starting from Hypothesis 1&2 below, we show in this section how to solve completely the spray flame problem using simple arguments. All these hypotheses and statement become rigourous in the next section. In this study, the main parameter is the size of the incoming droplet. When this parameter varies, it is crucial to keep constant the temperature T b of the burnt gases. Equivalently, this means that we want to compare flames whose total amount of (liquid and/or gaseous) incoming fuel is the same, which is a natural constraint. The velocity, directly proportional to the burning rate, therefore measures the efficiency of the combustion. For example, the total mass M u := 4/3 π n u R u 3/2 of incoming liquid fuel is kept constant by adjusting the number density n u of the incoming droplets when R u varies.

Hypothesis 1: Combustion occurs instantaneously at the flame front x f := 0 at the temperature T b of the burnt gases.

Hypothesis 2: For small enough sizes of the incoming droplets, the vaporisation finishes strictly before the reaction zone.

Here x f is the position of the flame front, taken as the origin, x v is the position of the beginning of the vaporisation zone, i.e. T (x v ) = θ v , and x vf is the position of the vaporisation front, that is the first point where the size of the droplets vanishes. For small droplets, Hypothesis 2 implies that the velocity of the wave is equal to the Zeldovich velocity c m := 2/d /T b of the classical large activation energy monophasic gas flame. Indeed, the internal structure of the flame involves only the unknowns T and Y as well as the corresponding profiles in the hot gases, all of which are not affected by the presence of the droplets. Increase now the size of the droplets, keeping M u constant. Assuming that the corresponding profiles vary continuously with respect to the size of the droplets, define R c > 0 as the critical smallest size at which the droplets finish vaporising exactely at the flame front. This means in particular that for all R u ∈ ]0, R c [, the velocity of the flame is equal to c m . We have:

Statement. Assume for simplicity that the vaporisation law is the step-function g(T ) = g 0 H (Tθ v ). Let R c > 0 be the critical size of droplets defined above. Then we have:

R c = g 0 c m 2 ln T b θ v (17)
Indeed, according to the vaporisation law, for a given size R u of incoming droplets, the quantities c(R u ), x v (R u ) and x vf (R u ) obey the relation:

x vf = x v + c g 0 R u = 1 c ln (T b θ v ) + c g 0 R u (18)
Then, by continuity of the solution with respect to the parameter R u , on the one hand we must have c(R c ) = c m and on the other hand, the definition of R c implies x vf (R c ) = 0. These two relations give the relation for R c in the Statement above.

Now, what happens now if R u > R c ? We have:

Statement. In the case where R u > R c , we have necessarily x vf (R u ) = 0 and therefore

c(R u ) = c d (R u ) := g 0 R u ln T b θ v ( 19
)
where c d designates the velocity of the flame in the diphasic case as opposed to the monophasic case abovementionned. The second part of the statement is a direct consequence of the first part, after solving the equation x vf (R u ) = 0 for c. Then, it is impossible that x vf > 0. Otherwise, vaporisation, combustion and finally an elevation of temperature would occur after the flame front, which is prohibited by the fact that the temperature is already equal to its maximum value at x = 0. On the converse, x vf (R u ) < 0 is impossible. Indeed, one checks easily that we would have x vf (R) < 0 for all R < R u , a contradiction with x vf R c = 0.

As a conclusion, for any R u , we have the alternative:

Statement. Let R u be given. Then, the limiting system is completely determined in terms of the following alternative involving the data of the problem:

• R u ≤ R c is equivalent to {x vf (R u ) < 0 and c(R u ) = c m } • R u > R c is equivalent to {x vf (R u ) = 0 and c(R u ) = c d (R u ) ≤ c m }
As announced, we observe that the velocity of the spray-flames diminishes for size of the incoming droplets larger than a critical value. This change of regime occurs abruptly at that critical value. We recover here the existence of the so-called "vaporisation-limited" regime emphasized in [HNS]. In particular, it follows from these remarks that substituting some gas fuel for some liquid fuel in the fresh gas can only lead to a decrease of the burning rate, i.e. the velocity, of the flame (at least in the scope of the hypotheses of this Section). If one sets the origin x = 0 in the fresh gas where T (x = 0) = T 0 < T b , then the decrease of the velocity of the flame means an increase of the preheating zone. In other words, large droplets push the reaction front towards the hot gases.

We invite the reader to verify that we would obtain exactely the same conclusions by replacing Hypothesis 2 above by the following:

Hypothesis 3: The velocity of the spray flame is always lower or equal to the velocity c m of the classical gas-flame having the same temperature in the burnt gas.

The rigourous analysis of the structure of the spray flame in the large activation energy limit precisely aims at proving Hypotheses 1 and 3 plus some more information about the structure of the internal combustion layer.

Large activation energy values : existence of flames and their structure.

In the rigourous proof of the convergence of the system to a limiting travelling wave, the main difficulty is to obtain uniform estimates from below as well as from above for the velocity. Thanks to the results of the previous section, this directly yields a solution to the limiting system in the form [START_REF] Greenberg | Instability of a flame front propagating through a fuelrich droplet-vapour-air cloud[END_REF] above. In order to perform these estimates, one has to analyse the internal structure of the reaction zone. Following Zeldovich, this is done by introducing scaled variables and coordinates as we shall see below. First notice that the choice of the value T (0) := T b -Aεlnε provides an information about the velocity of the flame in terms of the temperature profile at the origin:

First expression of the velocity of the flame

Thanks to the remark that, for A large enough, one has

f ε (T ε )Y ε = O(ε A/2 ) whenever T ε ≤ T b + Aεlnε, it follows that T ε (x) = (T b + Aεlnε) e cεx + O(ε A/2 ) on ] -∞, 0].
Integrating the equation for T ε between -∞ and 0:

-T ε ′ (0) + c ε (T b + Aεlnε) = -∞ 0 f ε (T ε )Y ε , (20)
where the Right-Hand-Side is O(ε A/2 ), we have therefore

c ε = T ε ′ (0) T b + O(εlnε). ( 21
)
It remains to estimate the derivative of the temperature profile at the origin, which we do by studying the system on the intervall[0, +∞[.

Notations. Streched system.

The boiling point is denoted by x v,ε and satisfies T ε (x v,ε ) = θ v , the ignition point is denoted by x i,ε and satisfies T ε (x i,ε ) = θ i , the end of vaporisation occurs at x vf,ε satisfies R ε (x vf,ε ) = 0, and we recall that T ε (0) = T b + Aεlnε. Also, we use the auxiliary functions W ε := T ε + Y ε -T b and Z ε := T ε + dY ε -T b . Now, the exponent γ ≤ 2 being given, we define the streched variables by setting:

ξ = x ε γ/2 , p ε (ξ) = T ε (x) -T b ε , q ε (ξ) = Y ε (x) ε , w ε (ξ) = W ε (x) ε , z ε (ξ) = Z ε (x) ε . ( 22
)
p ε (0) = A lnε, p ε (∞) = 0, q ε (∞) = 0, w ε (∞) = 0, z ε (∞) = 0
Moreover, we naturally denoted here

(ξ v,ε , ξ i,ε , ξ vf,ε ) = (x v,ε , x i,ε , x vf,ε )/ε γ/2
. Finally, it is important to recast the expression of the velocity in terms of the streched system, namely, for A fixed large enough:

c ε := ε 1-γ/2 p ε ′ (0) T b + O(εlnε). ( 23
)
and it remains to estimate the derivative of the temperature profile in the streched variables.

First Estimate of p ε ′ (0)
The estimate of p ε ′ (0) goes in two steps. First, we prove here that this quantity is bounded uniformly with respect to ε. We provide an optimal value of the upper-bound later.

It is convenient to write the system in terms of the unknowns (p ε , w ε ) or (p ε , z ε ) defined above. Simple algebra yields:

                 -p ε ′′ + c ε e γ/2 p ε ′ = (w ε -p ε ) e pε = - 1 d (z ε -p ε ), -w ε ′′ + c ε ε γ/2 w ε ′ = 1 -d d c ε ε γ/2 p ε ′′ + 2 π n u c ε √ ε 5γ 4 -1 (ξ vf,ε -ξ) + , -z ε ′′ + c ε ε γ/2 d z ε ′ = 1 -d d c ε ε γ/2 p ε ′ + 2 π n u c ε √ ε 5γ 4 -1 (ξ vf,ε -ξ) + . ( 24
)
First, integrating the equation for the temperature against p ε ′ gives:

c ε ε γ/2 0 ∞ (p ε ′ ) 2 d ξ + (p ε ′ (0)) 2 2 = 0 ∞ (w ε -p ε ) e pε p ε ′ d ξ = 1 d 0 ∞ (z ε -p ε ) e pε p ε ′ d ξ ( 25
)
where the contribution of p ε in the last two integrals can be integrated explicitely. It remains to estimate w ε or z ε in terms of p ε . For that, we distinguish the two cases 0 < d ≤ 1 and d ≥ 1.

Case 0 < d ≤ 1: An integration of the equation for w ε together with the condition at infinity w ε ′ (+∞) = 0 yields, since the contribution of the vaporisation term is nonnegative:

w ε ′ (ξ) ≥ 1 -d d η= ξ η=+∞ exp -c ε ε γ/2 (ξ -η) p ε ′′ (η) dη = 1 -d d exp (c ε ε γ/2 (ξ -η)) p ε ′ (η) η=ξ η=+∞ + η= ξ η=+∞ exp (c ε ε γ/2 (ξ -η)) c ε ε γ/2 p ε ′ (η) d η ≥ε 1 -d d (-p ε ′ (ξ)) (26)
because 0 < d ≤ 1 and the last integral is positive. Finally, since w ε (+∞) = 0, and p ε ′ ≥ 0, we have:

∀ξ ≥ 0, w ε (ξ) ≤ |p ε (ξ)| = -p ε (ξ) (27)
Hence, Eq. ( 25) gives

p ε ′ (0) 2 2 ≤ 1 -d d + 1 = 1 d (28)
Case d ≥ 1: Now, a direct integration of the equation for z ε yields, since the contribution of the vaporisation term is nonnegative:

z ε ′ (ξ) - c ε ε γ/2 d z ε (ξ) ≥ d -1 d c ε ε γ/2 p ε (ξ) (29) Hence z ε (ξ) ≤ d -1 d η= ξ η=+∞ exp c ε ε γ/2 d (ξ -η) c ε ε γ/2 (-p ε (η))dη ≤ d -1 d (-p ε (ξ)) η=ξ η=+∞ c ε ε γ/2 exp ε c ε ε γ/2 d (ξ -η) d η =(d -1)(-p ε (ξ)) (30)
where we have used the monotonicity of p ε . Hence, Eq. ( 25) gives

p ε ′ (0) 2 2 ≤ (d -1) + 1 d = 1 (31)
As a conclusion, p ε ′ (0) is uniformly bounded with respect to ε and the boundedness of c ε follows from [START_REF] Silverman | Asymptotic analysis of a premixed polydisperse spray flame[END_REF] as well as: Lemma 12. Let γ < 2 be given. Then, the velocity c ε tends to zero as ε goes to zero.

Hence, in the rest of this section, we shall only consider the case where γ = 2.

Approximate System: Come back now to the unknowns (p ε , z ε ). It is convenient to introduce the approximate problem:

       -p ˜ε′′ = 1 d (z ˜ε -p ˜ε) e p ˜ε, -z ˜ε′′ = 2 π n u c ε √ ε 3/2 (ξ vf,ε -ξ) + , (32)
with the same boundary conditions as for (p ε , z ε ), but where we omitted the convection terms involving the factor c ε ε. Classical results on elliptic regularity ensure that (p ε , z ε ) -(p ˜ε, z ˜ε) as well as their derivatives (p ε ′ , z ε ′ ) -(p ˜ε′ , z ˜ε′ ) are of order ε, uniformely on [0, +∞[, in the space of continuous functions.

Estimate of the Vaporisation Front ξ vf,ε

It is now possible to estimate the size ξ vf,ε of the intervall of the vaporisation zone intersecting that of the reaction zone (i.e. ξ > 0): Lemma 13. There exists a constant C > 0 such that for ε small enough:

ξ vf,ε ≤ Cε -3/2 |lnε| 2/5 , or equivalently x vf,ε ≤ Cε 2/5 |lnε| 2/5
(33)

Proof. Integrate the equation for z ˜ε taking into account the vaporisation term. For that, recall that we have at the vaporisation front z ˜ε(ξ vf,ε ) = 0 and z ˜ε′ (ξ vf,ε ) = 0. It comes:

z ˜ε′ (ξ) = C ε 3/2 c ε 1/2 (ξ vf,ε -ξ) 3/2 ≥ 0 (34)
where C denotes a generic positive constant independent of ε. Integrating again, we have immediately, since z ε (+∞) = 0,

z ˜ε(ξ) = -C ε 3/2 c ε 1/2 (ξ vf,ε -ξ) 5/2 ≤ 0 (35)
Recalling that z ε = z ˜ε + O(ε), we have

-C ε 3/2 c ε 1/2 ξ vf,ε 5/2 + O(ε) = z ε (0) = p ε (0) + dq ε (0) ≥ A lnε (36)
because q ε (0) ≥ 0, and the result follows thanks to the boundedness of c ε .

Second Estimate of p e ′ (0)

We precise here the upper-bound and give a lower-bound for the velocity and prove: Lemma 14. Let d be given. We have for a certain c > 0, and for ε small enough:

c 2 T b 2 2 ≤ (p ε ′ (0)) 2 2 ≤ 1 d + 0(ε 2/5 |lnε|) (37)
Proof. Since z ˜ε above is negative, it follows that z ε ≤ O(ε), so that we have from [START_REF] Suard | Vaporisation controlled regime of flames propagating in fuel-lean sprays[END_REF] the uniform upper-bound:

∀d, p ε ′ (0) 2 2 ≤ 1 d + O(ε) (38)
This inequality becomes an equality in the case where limsup ε→0 ξ vf ,ε < 0 and we have therefore in that case a lower-bound for the velocity of the limiting flame. This happens in particular in the case where limsup ε→0 x vf ,ε < 0 (which implies ξ vf ,ε → -∞) that is when the droplets finish vaporising strictly before the reaction front x = 0. This remark allows us to use the same alternative (see Lemma 8) as in the previous section. Namely, using the notations of the heuristic analysis, choose 0 < c < c d (R u ). Then, for ε small enough, c ε ≤ c implies that the droplets finish vaporising strictly before the reaction front. Hence, c ε ≃ c m . This completes the proof.

Large activation energy values :

Limiting system as ε → 0.

It remains to prove that the limiting system in the High Energy Activation limit is indeed system [START_REF] Greenberg | Instability of a flame front propagating through a fuelrich droplet-vapour-air cloud[END_REF] with the velocity c as in Theorem 11.

Proof. (of Theorem 11 -Continued) Following [START_REF] Berestycki | Traveling wave solutions to combustion models and their singular limits[END_REF], the uniform smoothness in C 1 (R) of the sequence T ε together with the uniform lower-bound for c ε allow us to extract converging subsequences T ε ′, x vf ,ε ′, x v,ε ′, c ε ′ on any large enough compact intervall around the origin. Then, let x vf ,1 ≤ 0 and x vf ,2 ≤ 0 be two limits of subsequences of the sequence of vaporisation fronts x vf ,ε . If x vf ,1 < x vf ,2 < 0 then necessarily we have for the corresponding velocities c 1 = c m = c 2 , which implies x vf ,1 = x vf ,2 , a contradiction. Now, recalling the definition of R c in the heuristic analysis, if x vf ,1 < x vf ,2 = 0, we must have both R u < R c and R u ≥ 0, a contradiction again. As a conclusion, we have necessarily x vf ,1 = x vf ,2 , hence also c 1 = c 2 , so that the limiting profile is unique, and all the sequence of solutions tend towards their respective limits.

Extensions

Singular limits for vaporisation and/or combustion.

The case of rapid vaporisation: There are situations where the vaporisation of the droplets occurs in a very small region compared to the size of the heating zone, that is situations where the time for complete vaporisation is very small. This may happen in two situations. First, when the droplets are very small, we have R u = η 1 R u,0 where η 1 is a small parameter and R u,0 is a reference size of order 1. In that case, the density number of droplets in the fresh gases is normalised to

n u = n u,0 /η 1 3/2
. This allows us as usual to maintain a constant temperature in the burnt gases

T b = Y u + 4/3 π n u,0 R u,0
3 /2 while letting η 1 go to 0. Second, when the vaporisation rate is very large, me may normalize it as g(T ) = g 0 (T )/η 2 where η 2 is a small parameter. Recall here that g 0 (T ) is an increasing function, positive for temperatures larger than θ v . We have:

Set the origin as x v = 0, such that T (0) = θ v . Then, when η 1 and/or η 2 tend to 0, the solution (T , Y , R, c) η 1 ,η 2 of problem (1) converges to a solution of the following system, where complete vaporisation occurs instantaneously at x v = 0:

   -T ′′ + c T ′ = f (T )Y -d Y ′′ + c Y ′ = -f (T )Y + c 4 3 πn 1 R u,1 3/2 δ(x = 0)
with the usual boundary conditions

Proof. The details of the proof are left to the reader. We refer to [START_REF] Berthonnaud | Contribution à la Modélisation et à l'Etude Mathématique des Ecoulements Diphasiques Turbulents et Réactifs[END_REF] for more details.

The case of rapid vaporisation and combustion: The previous limit is easily combined whith the high activation energy asymptotics. In that case, when rapid vaporisation and reaction occur at the same time, the limiting system has two singular terms. Both the vaporisation and the reaction term are concentrated in infinitely small regions. This is in some sense a simplified situation where vaporisation necessarily ends before the reaction occurs, since θ v < T b . Hence, the velocity of the limiting wave is obviously c = c m .

Results for a polydisperse spray.

All the previous results can be extended in the case of polydisperse sprays, that is in the case where droplets of possibly different sizes are present at a given position. In that case, we need a statistical description of the distribution of size of the droplets as was first pointed out by Williams [?, 26]. For our problem, we assume more precisely that the distribution of size of droplets is known in the fresh gases. Introduce for that the probability density number h u (R) of droplets in the fresh gases, such that h u (R) d R is the probable number of droplets having size in the range [R, R + d R]. We define in the same manner the density number h(x, R) at any position x and we have the limiting condition in the fresh gas h(-∞, R) = h u (R). Also, we assume that h u is continuous on R + and has compact support. Then, the spray flame problem writes:

       -T ′′ + c T ′ =Y f (T ), -d Y ′′ + c Y ′ =-Y f (T ) + 2 π g(T ) R R 1/2 h(x, R) d R, c∂ x h + ∂ R (h( -g(T ) H(R)) =0, (39)
In this setting, the total mass fraction M u of liquid fuel in the fresh gas, and therefore the corresponding temperature T b of the burnt gases are given by:

M u = 0 +∞ 4 3 π R 3/2 h u (R) d R, T b = Y u + M u (40)
Notice finally that a monodisperse spray corresponds to the situation where h u writes as:

h u (R) = n u δ (R -R u ) (41)
Existence of Travelling Waves: The proof of the existence of a travelling wave for the polydisperse spray follows exactely the same lines as in the monodisperse case. In the alternative stated in Lemma 8, page 4, the threshold velocity may now be written in terms of R m := sup {R; h u (R) > 0} instead of R u .

Asymptotic Limits: The limiting profiles are again of two types, but now the critical radius R c must be compared to R m defined above in order to know if the regime of combustion is controlled by vaporisation or not. Also, in the case of infinitely rapid vaporisation, the vaporisation eventually concentrates in an infinitely small region. Hence, there is no difference between a monodisperse and polydisperse spray in this asymptotics. We refer the reader to [START_REF] Berthonnaud | Contribution à la Modélisation et à l'Etude Mathématique des Ecoulements Diphasiques Turbulents et Réactifs[END_REF] for more details.

Other geometries.

The large activation energy limit analysis can be carried out in other onedimensionnal geometries, such as the anchored flame on the half-line, or counter-flow-like configurations, or finally radial in dimension 2 or 3. In each case, we are able to solve explicitely the problem (provided the expression of the vaporisation law is relatively simple), and it is possible to study the influence of different parameters on the existence of a profile or on the value of the burning rate for example. Relevant parameters are the typical velocities of injection of the gas and/or droplets, the value of the vaporisation rate, the dimension. We refer the reader to [START_REF] Berthonnaud | Contribution à la Modélisation et à l'Etude Mathématique des Ecoulements Diphasiques Turbulents et Réactifs[END_REF] for more details.

Perspectives and Concluding Remarks

This paper considered the existence of travelling fronts for a simple onedimensionnal thermodiffusive lean spray flame model. We proved the existence of travelling waves for general combustion and vaporisation laws. As far as qualitative results are concerned, the most important part of the paper is the study of the high activation energy limit for the system. Indeed, the limiting problem involve simple explicit profiles but still retains the important features of the dynamics. Extensions of these results to the cases of fast vaporisation, polydisperse sprays or other onedimensionnal geometries was briefly mentionned. The corresponding details can be found in [START_REF] Berthonnaud | Contribution à la Modélisation et à l'Etude Mathématique des Ecoulements Diphasiques Turbulents et Réactifs[END_REF].

The present work is therefore a first step towards the rigorous derivation of asymptotic models for spray flames. However, much remains to be done in order to understand the effect of droplets on the dynamics of the flame structure, as observed by many physicists and experimentalists (see [START_REF] Greenberg | Instability of a flame front propagating through a fuelrich droplet-vapour-air cloud[END_REF] and the references therein). First, as already noticed in the Introduction, the model under consideration in this paper omits the effects of the latent heat. It would be very interesting to derive some high activation energy models incorporating these effects, in the hope of deriving explicit expressions of the limiting profiles and combustion rates, as well as extinction limits. Second, a striking consequence of the high activation energy asymptotics is that droplets cannot cross the flame front in that limit, but may enter the reaction zone only for large but not infinite values of the activation energy. It would be interesting to derive intermediate asymptotic models making explicit the structure and thickening of the combustion region for large but not infinite .

As far as dynamical phenomena are concerned, the problem of the stability of spray flame system are of particular importance. As a first step, one has to study the nonlinear stability for small perturbations of the thermodiffusive lean spray flame model (see [START_REF] Berthonnaud | Contribution à la Modélisation et à l'Etude Mathématique des Ecoulements Diphasiques Turbulents et Réactifs[END_REF][START_REF] Berthonnaud | Stability of travelling waves solutions in a one-dimensional two-phase flow combustion model[END_REF]). As a second step, we are interested in the problem of acoustic instabilities in spray flame systems. In this direction, the work Clavin an Sun [START_REF] Clavin | Acoustic instabilities in two-phase combustion[END_REF] is of particular interest to us because these authors propose a framework for the mathematical analysis of gas or spray flames coupled to the acoustics. This mathematical analysis relies on the possibility of deriving explicit expressions of the solutions to the problem. This is only possible in some asymptotic limits. For that, the authors consider the large activation energy limit for combustion phenomena and assume an instantaneous vaporisation of the droplets. As a consequence, both the combustion and the vaporisation zones reduce to infinitely small regions, whose internal layer is analysed. However, assuming that the vaporisation zone is very small compared to the preheating zone is a very restrictive assumption in many applications where the droplets can spread into the preheating zone, approach the combustion zone, or even cross the combustion zone. We have shown in the present paper that it is possible to handle analytically such situations where the droplets approach or reach the reaction front, and where the slowly vaporising droplets induce a dramatic change of the combustion rate. It would therefore be intersting to study mathematically and generalise the work of Clavin and Sun [START_REF] Clavin | Acoustic instabilities in two-phase combustion[END_REF] to more general situations regarding the size of the vaporisation zone. It is reasonable to expect that different size of the vaporising region would give rise to different behaviours or stability properties of the system.

  which is the desired result. Proof of the Upper-Bound of Proposition 7 We have the Proposition 10. Let M = sup [0,T b ]

2

 2 MTb θv ), and we compute an upper bound of c independent of a. If c ≤ c 0 , the result is proved. Conversely, if c > c 0 , then since θ v ≤ T b e -ca + θv 2 , we have that e ca ≤ 2Tb θv , and we conclude that: c ≤ ln 2Tb θv a 0 .

Existence of travelling waves and large activation energy limits for a onedimensional thermodiffusive lean spray flame model

Pierre BERTHONNAUD, Komla DOMELEVO

Acknowledgments. It is our pleasure to thank J.-M. Roquejoffre for suggesting the study and many advices. Also, we would like to thank P. Haldenwang for very interesting and stimulating discussions.


Elisabeth Mr Robert
email: erobert@vub.ac.be

Nele Schmitz

Paul Copini

Edo Gerkema

Frank J Vergeldt

Carel W Windt

Hans Beeckman

Nico Koedam

Henk Van As

Elisabeth M R Robert

Elisabeth Robert

Visualization of the stem water content of two genera with secondary phloem produced by successive cambia through Magnetic Resonance Imaging (MRI)

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The secondary phloem tissue produced by successive cambia in stems of Avicennia trees is more abundant with increasing soil water salinity and decreasing inundation frequency within the mangrove forest and woody shrub and tree species with this type of secondary phloem mainly occurs in either dry or saline environments [START_REF] Robert | Successive cambia: a developmental oddity or an adaptive structure[END_REF]. Hence, we hypothesize that secondary phloem tissue produced by successive cambia allows for increased water storage. Phloem water storage may be advantageous in environmental conditions with variable water availability and frequent drought. Through estimation of the effective water content of woody species with secondary phloem tissue produced by successive cambia, and the variation therein across the stem, we may gain insight into the role of this secondary phloem in the water balance of drought and salt tolerant woody plants, such as Avicennia and Bougainvillea [START_REF] Knox | Salt Tolerance of Landscape Plants for South Florida[END_REF][START_REF] Knox | Drought Tolerant Plants for North and Central Florida[END_REF][START_REF] Kobayashi | Bougainvillea. Ornamentals and Flowers: Oct. 2007 -OF-38[END_REF] individuals. Imaging methods based on nuclear magnetic resonance (NMR), called magnetic resonance imaging (MRI), are among the techniques that allow for direct in vivo observation of water content and water movement in

The$Journal$of$Plant$Hydraulics$1:$e70006$$ ! 2! plants ( [START_REF] Van As | MRI of intact plants[END_REF]. Other methods are heat tracing methods (heat balance, heat pulse or heat dissipation) and radioisotope or stable isotope tracing methods [START_REF] Calder | Investigation into the use of deuterium as a tracer for measuring transpiration from Eucalypts[END_REF][START_REF] Lambs | Sap flow measurements of Ceriops tagal and Rhizophora mucronata mangrove trees by deuterium tracing and lysimetry[END_REF]. However, the exceptional aspect of MRI is the possibility for spatially resolved results, allowing for sap flow and water content data on every location of the plant's stem [START_REF] Van As | MRI of intact plants[END_REF]. In the last decades, NMR flow imaging has become a growing discipline in the study of xylem and phloem sap transport over long and short distances in plants (e.g. [START_REF] Kuchenbrod | Measurement of water flow in the xylem vessels of intact maize plants using flowsensitive NMR imaging[END_REF][START_REF] Rokitta | Dynamic studies of phloem and xylem flow in fully differentiated plants by fast nuclear-magnetic-resonance microimaging[END_REF][START_REF] Olt | Na-23 and H-1 NMR micro-imaging of intact plants[END_REF][START_REF] Holbrook | In vivo observation of cavitation and embolism repair using magnetic resonance imaging[END_REF][START_REF] Scheenen | Functional imaging of plants: a nuclear magnetic resonance study of a cucumber plant[END_REF][START_REF] Scheenen | Intact plant magnetic resonance imaging to study dynamics in long-distance sap flow and flow-conducting surface area[END_REF][START_REF] Van As | Intact plant MRI for the study of cell water relations, membrane permeability, cell-to-cell and longdistance water transport[END_REF][START_REF] Windt | Nuclear magnetic resonance imaging of sap flow in plants[END_REF][START_REF] Kaufmann | Functional repair of embolized vessels in maize roots after temporal drought stress, as demonstrated by magnetic resonance imaging[END_REF][START_REF] Mullendore | Sieve tube geometry in relation to phloem flow[END_REF] with few studies on trees [START_REF] Kuroda | Magnetic resonance micro-imaging of xylem sap distribution and necrotic lesions in tree stems[END_REF][START_REF] Windt | MRI of long-distance water transport: a comparison of the phloem and xylem flow characteristics and dynamics in poplar, castor bean, tomato and tobacco[END_REF][START_REF] Homan | 0.7 and 3 T MRI and sap flow in intact trees: xylem and phloem in action[END_REF][START_REF] Merela | A single point NMR method for an instantaneous determination of the moisture content of wood[END_REF][START_REF] Umebayashi | The developmental process of xylem embolisms in pine wilt disease monitored by Multipoint Imaging using Compact Magnetic Resonance Imaging[END_REF][START_REF] Schepper | MRI links stem water content to stem diameter variations in transpiring trees[END_REF]. All NMR experiments are based on the interaction between the magnetic moment of certain nuclei, especially 1 H, and an external magnetic field [START_REF] Köckenberger | Functional imaging of plants by magnetic resonance experiments[END_REF][START_REF] Van As | Intact plant MRI for the study of cell water relations, membrane permeability, cell-to-cell and longdistance water transport[END_REF]. A strong static magnetic field is used to create sample magnetization that is then manipulated by radio frequency (RF) pulses and detected through induction of a weak voltage in a coil placed around the sample [START_REF] Köckenberger | Functional imaging of plants by magnetic resonance experiments[END_REF]. The magnetization is measured in the presence of a magnetic gradient, allowing for one-dimensional projection in the direction of the gradient after Fourier transformation of the detected signal as a function of time [START_REF] Köckenberger | Functional imaging of plants by magnetic resonance experiments[END_REF]. It has already been proven that the bark tissues of Quercus robur [START_REF] Schepper | MRI links stem water content to stem diameter variations in transpiring trees[END_REF] and Populus [START_REF] Windt | MRI of long-distance water transport: a comparison of the phloem and xylem flow characteristics and dynamics in poplar, castor bean, tomato and tobacco[END_REF] trees, 1-2 cm in stem diameter and 1.1-1.5 m in height, have the highest water content when compared to the other tissues within the stems of these trees. Also, the non-conductive bark cells in the same Quercus trees were put forward as contributing most to the water storage capacity of the tree stem [START_REF] Schepper | MRI links stem water content to stem diameter variations in transpiring trees[END_REF]. If the phloem tissue of the bark has an important role in the water household of plants [START_REF] Van As | MRI of intact plants[END_REF] and trees, it is probable that also other living tissues within a tree's stem, especially secondary phloem if present, significantly contribute to the tree's water reserves. It is generally accepted that phloem in plants with a single vascular cambium function for one year only [START_REF] Raven | Biology of plants[END_REF][START_REF] Rosner | Osmotic potential of Norway spruce (Picea abies (L.) Karst.) secondary phloem in relation to anatomy[END_REF][START_REF] Carlquist | Successive cambia revisited: ontogeny, histology, diversity, and functional significance[END_REF]. For plants with successive cambia, a greater area of conducting phloem area has been proposed, based on the continued production of secondary phloem by the internal vascular cambia [START_REF] Carlquist | Successive cambia revisited: ontogeny, histology, diversity, and functional significance[END_REF]. To our knowledge, physiological experiments providing evidence that secondary phloem tissue in plants with successive cambia is still active are lacking, with the exception of one study by [START_REF] Fahn | The organization of secondary conducting tissues in some species of the Chenopodiaceae[END_REF]. Water conduction can take place in more than the outermost stem part of Avicennia trees (Robert & Schmitz, unpublished data). The question then remains: what about the actual water content of secondary phloem tissue produced by successive cambia, if present, in woody species? Are the secondary phloem cells water-filled over the entire stem? Do they play a role in the axial and radial water transport in tree stems? A positive answer would confirm the ecological advantage that we suggested for species with secondary phloem tissue produced by successive cambia in conditions of physiological drought due to high soil water salinity or water shortage [START_REF] Robert | Successive cambia: a developmental oddity or an adaptive structure[END_REF]. In this study, we used MRI to visualize the stem water content of three species with successive cambia, the mangroves Avicennia marina and A. officinalis and the non-mangrove Bougainvillea spectabilis. We tested the hypothesis that not only the outermost phloem tissue has high water content but also the secondary phloem tissues over the entire stem from the bark inward to the pith, herewith serving as water storage sites. One tree of each species has been selected for the study: a Bougainvillea tree, cultivated in a temperate climate and directly subjected to MRI and two Avicennia trees taken from the mangrove forest, transported and measured, after an acclimatisation period in conditions different from its natural environment. The general approach is viewed as a preliminary and pioneering step to assess an extra line of evidence for the ecological significance of secondary phloem tissue produced by successive cambia.

Materials and Methods

Plant material

A Bougainvillea spectabilis WILLD. plant with a stem of 1.20 cm in diameter and ca. 80 cm in height, trimmed when cultivated as a single stemmed small tree, has been obtained from a local plant dealer in Belgium. Two small Avicennia trees have been transported from the mangrove forest of Gazi Bay (Kenya) (Avicennia marina (Forssk.) Vierh. -tree 1) and Leyte island (The Philippines) (Avicennia cf. officinalis Linn. -tree 2) to the Wageningen NMR Centre (Wageningen University, The Netherlands) where MRI measurements were conducted. Tree 1 was ca. 150 cm high and had a stem diameter of 1.25 cm while tree 2 was ca. 65 cm high with a stem diameter of 0.8 cm. During transport, the Avicennia trees were covered with insulation material and between arrival and actual MRI measurements there was an acclimatization period of 2.5 months in the mangrove greenhouse of Burger's Zoo in Arnhem (The Netherlands) during which the tree was watered with fresh water (tree 1), and of one week in a climate chamber (air temperature: 22ºC; relative air humidity: 70%) during which the tree was watered with saline water (20-30 g/l) (tree 2). While the xylem and phloem

The$Journal$of$Plant$Hydraulics$1:$e70006$$ ! 3! tissue in Avicennia trees has a reticulate structure [START_REF] Robert | Successive cambia: a developmental oddity or an adaptive structure[END_REF], numerous strands of vascular increments, composed of both xylem vessels and phloem tissue, are scattered over the stem of Bougainvillea (Figure 1).

Magnetic Resonance Imaging (MRI)

Measurements were done in a 3 T (Bougainvillea tree and Avicennia tree 1) and a 0.7 T MRI system (Avicennia tree 2) at the Wageningen NMR Centre (Wageningen University, The Netherlands). The 3T MRI system (Figure 2) consists of an Avance console (Bruker, Karlsruhe, Germany) and a superconducting magnet with a 50 cm vertical bore (Magnex, Oxford, UK), generating a magnetic field of 3 T (128 MHz proton frequency). Measurements were performed in the centre of a gradient coil that was mounted around the stem (Figure 2b-d) and that generates maximum gradient strengths of 1 T/m. For induction and detection of the signal, a radio frequency (RF) coil, consisting of two half cylinders, was placed around the tree stem at ca. 40 cm height above the pot surface. The 0.7 T MRI system consists of an Avance console and electromagnet (Bruker, Karlsruhe, Germany) with a 10 cm air gap, generating a magnetic field of 0.7 T (30 MHz proton frequency). Maximum gradient strength that can be generated is 1 T/m. For induction and detection of the NMR signal a solenoid RF coil was constructed by wrapping copper wire around a mould slightly wider than the tree stem. The coil was connected to a tuning circuit and electromagnetically shielded with aluminium foil and copper tape. For a more detailed description of the MRI systems: see [START_REF] Homan | 0.7 and 3 T MRI and sap flow in intact trees: xylem and phloem in action[END_REF]. During the measurement periods air temperature ranged from 22ºC to 26ºC. Relative air humidity was 28-40% (3 T MRI) and 50-65% (0.7 T MRI). Trees were subjected to a day-night cycle of 12:12 hours and regularly watered with fresh water (Bougainvillea and Avicennia tree 1) and saline water (20-30g/l) (Avicennia tree 2). These conditions were chosen to maximally allow normal functioning of all three plants. The conditions imposed on the Avicennia trees were intended to represent natural conditions with high demand on the water transport system either imposed by very low relative humidity or by the combination of high relative humidity and saline water. For all quantitative T2 measurements, a multiple spin echo imaging sequence has been used [START_REF] Edzes | Quantitative T2 imaging of plant tissues by means of multiecho MRI microscopy[END_REF]. The following settings were applied: field of view: 20 x 20 mm (Bougainvillea), 18.5 x 18.5 mm (Avicennia 1) and 11.4 x 11.4 (Avicennia 1); matrix: 256 x 256; slice thickness: 3 mm; number of averages: 8 (Bougainvillea) and 16 (Avicennia); echo time: 7.19 ms (Bougainvillea and Avicennia 1) and 6.88 ms (Avicennia 2); number of echos: 96 (Bougainvillea and Avicennia 2) and 64 (Avicennia 1); repetition time: 2500 ms. Differences of values between the trees were to optimise the spatial (field of view) and the signal-tonoise and time resolution (number of averages). They furthermore depend on the longest relevant T2 (number of echoes). Several quantitative T2 measurements were made, both during day and night after which, measured proton density and T2 values were visualized over the stem disc through a colour scale from the smallest (black) to the largest (white) value. Proton density is a measure of the water content, while T2 is a time constant related to transverse relaxation after magnetization and radio frequency pulsing of a sample [START_REF] Köckenberger | Functional imaging of plants by magnetic resonance experiments[END_REF][START_REF] Van As | Intact plant MRI for the study of cell water relations, membrane permeability, cell-to-cell and longdistance water transport[END_REF]. T2 values change with the size of the measured compartment (van der Weerdt et al., 2001[START_REF] Van As | Intact plant MRI for the study of cell water relations, membrane permeability, cell-to-cell and longdistance water transport[END_REF] and are as such an indication of the tissue type in plants. No differences could be observed in water content and T2 patterns between the different measurements, neither in Avicennia, nor in Bougainvillea.

The$Journal$of$Plant$Hydraulics$1:$e70006$$ !

5!

Figure 3: Proton density (water content) on a gradient from 0 (black -low water content) to 100 (white -high water content) and T2 images of intact Bougainvillea spectabilis and Avicennia marina stems obtained through magnetic resonance imaging (MRI). The stem diameters are 1.20 cm and 1.25 cm for the Bougainvillea and Avicennia plant, respectively. The inset in the upper left panel is a zoom of two vascular strands, composed of phloem (Ph) and xylem (X) cells (Figure 1). Proton density is in arbitrary units (a.u.). Only proton density directly links to water content. T2 values represent the relaxation time after magnetisation and radio frequency pulsing and are strongly influenced by the degree of mobility of the protons and their physico-chemical environment with higher values for more mobile and lower values for less mobile protons. T2 values are as such also a measure of cell size/type.

The$Journal$of$Plant$Hydraulics$1:$e70006$$ ! 6!

Results and discussion

Amplitude images showed the water content of the different tissues over the entire stem of both Bougainvillea and Avicennia trees (Figure 3). The bark and the vascular bundles of Bougainvillea had high water content but the highest amount of water could be observed in the zone of the outermost vascular cambium (up to the first lignified xylem cells). This is in agreement with observations in Quercus [START_REF] Kuroda | Magnetic resonance micro-imaging of xylem sap distribution and necrotic lesions in tree stems[END_REF][START_REF] Schepper | MRI links stem water content to stem diameter variations in transpiring trees[END_REF], Populus [START_REF] Windt | MRI of long-distance water transport: a comparison of the phloem and xylem flow characteristics and dynamics in poplar, castor bean, tomato and tobacco[END_REF] and Pinus [START_REF] Umebayashi | The developmental process of xylem embolisms in pine wilt disease monitored by Multipoint Imaging using Compact Magnetic Resonance Imaging[END_REF] tree stems. High water content is expected in the cambial zone since water is crucial in cell formation and especially cell enlargement [START_REF] Woodward | Plant science: tall storeys[END_REF][START_REF] Arend | Seasonal change in the drought response of wood cell development in poplar[END_REF]. When looking in more detail to the water content of the vascular bundles of Bougainvillea, we see that in the outermost zone of many bundles, the water content was higher than in the innermost zone, visible by the lighter yellow-green and the slightly darker green colour respectively (Figure 3 -inset in above left panel). Since each vascular bundle in this species is composed of phloem on the outside and xylem on the inside (Figure 1), the higher water content was situated in the phloem tissue. The phloem cells in Bougainvillea stems are thus full of water, supporting their role in water storage and supply.

In Avicennia, amplitude images showed higher water content in the pith and in the outermost zone of the studied stems (Figure 3; tree 2: data not shown), although this was less pronounced than in the studied Bougainvillea tree. The phloem bands in both studied Avicennia trees had relatively high water content throughout the entire stem, despite the different treatments they got (tree 1: fresh water and low relative humidity versus tree 2: saline water and higher relative humidity). The considerable amount of water in the secondary phloem tissue of both Avicennia and Bougainvillea throughout the entire stem strengthens the importance attributed to the secondary phloem tissue in the water household of trees with successive cambia and other cambial variants. The presence of water storage tissues evenly distributed over the tree stem could indeed be responsible for their ecological advantage in conditions of water shortage as has been proposed earlier [START_REF] Robert | Successive cambia: a developmental oddity or an adaptive structure[END_REF]. T2 values are not radially homogeneous in the studied B. spectabilis and A. marina stems, which is indicated by the presence or absence of light green dots on the T2 images (Figure 3 -right panels) often in concentric zones. Since T2 values scale with compartment size (van der Weerdt et al., 2001;[START_REF] Van As | Intact plant MRI for the study of cell water relations, membrane permeability, cell-to-cell and longdistance water transport[END_REF], narrower vessels in some zones of the stems could possibly be at the basis of lower T2 values in these zones. However, to obtain further understanding of the physiology in terms of hydraulic function behind the observed patterns, sap flow measurements should be conducted.

Conclusion

We can conclude that the secondary phloem tissue produced by successive cambia of both Bougainvillea and Avicennia has high water content. This strengthens the contribution of this secondary phloem tissue to ecological success in conditions of physiological drought caused by water shortage or high soil water salinity that we proposed earlier [START_REF] Robert | Successive cambia: a developmental oddity or an adaptive structure[END_REF]. In this study we could only address young trees with narrow diameter. Further research should address the water content of the secondary phloem at the interior part of large trees with successive cambia in relation to the functionality of the xylem vessels. Our MRI results on Bougainvillea and Avicennia are preliminary and only a very few trees have been measured, yet we consider them significant. The nondestructive and real-time study of phloem transport and water exchange between the xylem and the phloem tissue in mature trees in the field is difficult [START_REF] Sevanto | Effects of the hydraulic coupling between xylem and phloem on diurnal phloem diameter variation[END_REF]. The study of secondary (not bark) phloem tissue in trees presents additional complications. However, further studies of the functioning of secondary phloem tissue produced by successive cambia and its role in the water household of plants could offer important insight in understanding the survival and thriving of plants under physiological drought and the functioning of trees in general. Combining techniques and new approaches may overcome the current difficulties and should therefore be continued. Studies on drought and/or salt tolerant woody plants combining MRI and dendrometer measurements (De Schepper et al., 2012), using portable NMR devices [START_REF] Windt | A portable Halbach magnet that can be opened and closed without force: the NMR-CUFF[END_REF] or specifically addressing the coupling between xylem and phloem through dendrometer analysis [START_REF] Sevanto | Effects of the hydraulic coupling between xylem and phloem on diurnal phloem diameter variation[END_REF] and NMR measurements of radial transport and (radial) tissue permeability [START_REF] Sibgatullin | Quantitative permeability imaging of plant tissues[END_REF] could be a first possible continuation. The specific effect of different treatments, e.g. watering with fresh water versus salty water, should also be addressed in long-term experimental set-ups or by comparing several trees from contrasting locations.

The$Journal$of$Plant$Hydraulics$1:$e70006$$ !

Figure 1 :

 1 Figure 1: Internal stem anatomy of Bougainvillea spectabilis (A, B) and Avicennia marina (C, D), showing the bark tissue (grey bars) (A, C), the secondary phloem tissue produced by successive cambia (white bars) and the xylem tissue (black bars) with water conducting vessels (v) (B, D). Stems were sectioned using a sliding microtome after which a picture was taken via a camera connected to a microscope. The stem anatomy of A. officinalis is highly similar to the one of A. marina.

Figure 2 :

 2 Figure 2: 3T MRI system at Wageningen University (The Netherlands) showing (a) the lifting of a tree into the superconducting magnet, (b, c, d, e) the gradient coil, (b) open and (c, d, e) placed around a tree, (d) Edo Gerkema preparing an oak tree for MRI analysis in the space under the magnet and (f) the 3T MRI building from outside. The grey arrow indicates the length of a 1.80 m high tree, from roots to crown top. Figures a, b and c are adapted from Homan et al., 2007. Picture d is taken by Paul Copini.

Acknowledgements

We are grateful to Hamisi Ali Kirauni and James Gitundu Kairo of the Kenya Marine and Fisheries Research Institute (KMFRI, Mombasa, Kenya), to Mohamed Omar Said Mohamed of Kenya Wildlife Service and to Hilde Robert and Mia Van Mellaert for their help and support in the transport of the Avicennia marina tree; to Joep Wensing and Ernst Kamphuis of Burger's Zoo (Arnhem, The Netherlands) for housing the Avicennia plants in the mangrove greenhouse of Burger's Zoo; to Arie van Hoek and John Philippi of the Wageningen NMR Centre (Wageningen University, The Netherlands) for their assistance during scanning and scanning preparation and to Hilde

The$Journal$of$Plant$Hydraulics$1:$e70006$$ ! 7! Robert for her assistance in making thin-sections. The research was financial supported by the Agency for Innovation by Science and Technology (IWT, Flanders, Belgium), the Research Foundation -Flanders (FWO, Flanders, Belgium), the Company of Biologists (COB) through the Society for Experimental Biology (SEB), the European Cooperation in Science and Technology (COST) on the Terrestrial Biosphere in the Earth System (TERRABITES) via a Short Term Scientific Mission and the "Stichting ter bevordering van wetenschappelijk onderzoek in Afrika" (SBWOA, Belgium). This work has been carried out under the framework of the project 'CREC' (EU IRSES #247514) and the COST FP1106 network STReESS.


Bertrand Neveu
email: bertrand.neveu@enpc.fr

Gilles Trombettoni
email: gilles.trombettoni@lirmm.fr

Ignacio Araya

Node Selection Heuristics Using the Upper Bound in Interval Branch and Bound

Keywords: Global Optimization, Interval Branch and Bound, Node Selection 1

We present in this article a new strategy for selecting the current node in an interval Branch and Bound algorithm for constrained global optimization. The standard best-first strategy selects the node with the lowest lower bound of the objective estimate. We propose in this article new node selection policies where an upper bound of each node/box is also taken into account. The good accuracy of this upper bound achieved by several operators leads to a good performance of the criterion. These new strategies obtain better experimental results than classical best-first search on difficult instances.

Introduction

The paper deals with continuous global optimization (nonlinear programming) deterministically handled by interval branch and bound (B&B). Several works have been performed for finding good branching heuristics [START_REF] Csendes | Subdivision Direction Selection in Interval Methods for Global Optimization[END_REF], but very little work for the node selection itself. The solvers generally follow a best first search strategy (BFS), with some studies for limiting its exponential memory growth [START_REF] Casado | Experiments with a New Selection Criterion in a Fast Interval Optimization Algorithm[END_REF][START_REF] Ninin | A Metaheuristic Methodology Based on the Limitation of the Memory of Interval Branch and Bound Algorithms[END_REF]. Different variants of BFS have been proposed for discrete problems, such as K-Best-First Search [START_REF] Felner | KBFS: K-Best-First Search[END_REF].

To our knowledge, only Casado et al. in [START_REF] Casado | Experiments with a New Selection Criterion in a Fast Interval Optimization Algorithm[END_REF] and Csendes in [START_REF] Csendes | New Subinterval Selection Criteria for Interval Global Optimization[END_REF] proposed node selection heuristics for interval B&Bs. One criterion to maximize, called C 3 in [START_REF] Markot | New Interval Methods for Constrained Global Optimization[END_REF], and suitable for unconstrained global optimization is:

f * -lb ub -lb where [lb, ub] = [f ] N ([x]
) is the interval obtained by the natural interval evaluation of the realvalued objective function f in the current box [x]: [lb, ub] is a range interval including all real images of any point in [x] by f . f * is the optimum. When f * is not known, f , the cost of the best feasible point found so far can be used as an approximation of f * . This criterion favors small boxes (i.e., a small interval width (ublb)) and nodes with good lb. For constrained optimization, another criterion (to maximize) called C 5 is equal to C 3 ×f r. It takes into account a feasibility ratio f r computed from all the inequality constraints. The criterion C 7 proposes to minimize lb C 5 . This paper proposes two other policies taking into account an accurate upper bound of the optimum cost.

Standard Interval B&B

An interval

[x i ] = [x i , x i ] defines the set of reals x i s.t. x i ≤ x i ≤ x i . A box [x] is a Cartesian product of intervals [x 1 ] × ... × [x i ] × ... × [x n ].
The paper deals with continuous global optimization under inequality constraints defined by: min x∈[x] f (x) subject to g(x) ≤ 0, where f : R n → R is the real-valued objective function and g : R n → R m is a vector-valued function. x = {x 1 , ..., x i , ...x n } is a vector of variables varying in a domain/box [x]. x is said to be feasible if it satisfies the constraints.

A standard interval (or spatial) B&B scheme for continuous constrained global optimization (also known as nonlinear programming) is described in algorithms below. The B&B maintains two main types of information during search: f , which is the cost of the best feasible point found so far, and f min the minimum value of the lower bounds lb([x]) of the boxes/nodes [x] to explore. In other terms, in every box [x], there is a guarantee that no feasible point exists with a cost lower than lb([x]).

Let us first ignore the bold-faced instructions corresponding to the new strategy and detailed in Section 3. The algorithm is launched with the set g of constraints, the objective function f and with the initial box put into a list B of boxes to be handled. ǫ obj is the required precision on the objective cost and is used as stopping criterion. We add a variable x obj in the system (the vector x of variables) corresponding to the image (cost) of the objective function, and a constraint f (x) = x obj . The criterion generally used in existing interval B&Bs is denoted x obj ≤ f . There exist two phases in a branch and bound: a phase where we try to find the optimal solution, and a second phase where we prove that this solution is optimal, which requires us expand all the remaining nodes. Therefore the node selection matters only in the first phase.

We define new strategies aggregating two criteria for selecting the current box:

1. LBBox: The well known criterion used by BFS and minimizing lb([x]) (for all boxes [x] in the set B). This criterion is optimistic since we hope to find a solution with cost f min , in which case the search would end. For each box, lb([x]) is computed by the Contract&Bound procedure and the computed value labels the node stored into the set B of boxes.

2. UBBox: This criterion selects the node having the smallest goal upperbound. Thus, if a feasible point was found inside this box, it would more likely improve the best cost found so far.

The UBBox criterion is symmetric to the LBBox one: for every box, ub([x]) is computed by the Contract&Bound procedure and labels the node before storing it in the set B of boxes. In particular, constraint programming techniques like 3BCID [START_REF] Trombettoni | Constructive Interval Disjunction[END_REF] can improve ub([x]) by discarding small slices at the upper bound of [x obj ] (shaving process).

We think that a key of success of the UBBox criterion is that it evaluates more accurately the objective upperbound than the natural interval evaluation would do (i.e., ub computed by

[lb, ub] := [f ] N ([x])).
We propose two main ways to aggregate these two criteria.

LB+UBBox.

This strategy selects the node [x] with the lowest value of the sum lb([x]) + ub ([x]). This corresponds to minimizing both criteria with the same weight, i.e. minimizing the middle of the interval of the objective estimate in the box.

Alternating both criteria.

In this second strategy, the next box to handle is chosen using one of the two criteria. A random choice is made by the criterionChoice function at each node selection, with a probability UBProb of choosing UBBox. If UBBox (resp. LBBox) is chosen and several nodes have the same cost ub([x]) (resp. lb([x])), then we use the other criterion LBBox (resp. UBBox) to tie breaks.

Experiments showed that the performance is not sensitive to a fine tuning of the UBProb parameter provided it remains between 0.2 and 0.8, so that the parameter has been fixed to 0.5. The experiments in Section 5 highlight the positive impact of this criterion on performance.

These results suggest that it is important to invest both in intensification (UBBox) and diversification (LBBox). In other words, the use of a second criterion allows the search to avoid the drawback of using one criterion alone, i.e. (for LBBox) choosing promising boxes with no feasible point and (for UBBox) going deep in the search tree where only slightly better solutions will be found trapped inside a local minimum.

Details on the criterion UBBox

All candidate boxes in the set B, have a UB label depending on the best cost found so far ( f ) when the boxes were handled by Contract&Bound. All labels fall in four main cost ranges categories given by f and explaining with which priority the boxes are chosen using the UBBox criterion. The label is: 1. lower than fǫ obj if the contraction procedure reduced the maximum estimate of the objective in the box, 2. equal to fǫ obj , if the box is a descendant of the box containing the current best feasible point f , 3. equal to f -0.9 ǫ obj if the box was handled after the last update of best cost, 4. greater than f -0.9 ǫ obj in the remaining case.

As shown in the Contract&Bound pseudocode, the additional term 0.1 ǫ obj allows penalizing the boxes that are not issued by a bisection from boxes where the current best feasible point was found.

4.

Implementation of the Set B of Boxes

The set B was initially implemented by a heap structure ordered on the LBBox criterion. The implementation s1-05 keeps this unique data structure for taking into account the two criteria in the randomized strategy, but the node selection using UBBox comes at a linear cost in the number N of nodes. In practice, this takes about 10% of the total cost when N exceeds 50,000.

We have then built a variant s1-05-01 changing dynamically the probability UBProb in the following way: UBProb = 0.5 if N ≤ 50, 000 and UBProb = 0.1 if N > 50, 000.

Finally, we tried a cleaner implementation, s2-nodiv, with two heaps, one for each criterion. All the operations are then in log 2 (N ), except for the heap filtering process launched occasionnally during search. (This "garbage collector" is performed in all our implementations and removes from B all the nodes with lb([x]) > f .)

The s2-nodiv implementation brings less diversification than the first one which reconstructs the heap at each criterion change. Indeed, there often exist several nodes with the same UBBox and LBBox values, so it is useful to periodically rebuild the data structures randomly for breaking the ties. So we propose variants that use a second parameter corresponding to a diversification period. We obtained good results by fixing it to 50 in the s2-50 variant or 100 in the s2-100 variant, the first parameter UBProb being still fixed to 0.5. To be fair, we also applied the first strategy with the minUB+LB aggregative criterion running heap filtering every 100 nodes (s0-100).

Experiments

We have run experiments on 82 problems, issued from the series 1 and 2 of the Coconut benchmark. The best strategies s2-50 and s2-100 obtain a gain of about 40% w.r.t. the total time and 23% on average, meaning that greater gains are obtained on difficult problems.

Summary

The node selection policy is a promising line of research to improve performance of interval B&B. We have obtained good results by taking into account for each node a lower bound but also an upper bound, provided that this upper bound is made accurate by box contraction operations, by a random selection between both criteria and by a work on heap data structures. In a short term, we are going to investigate how relevant components of criteria proposed by Markot and al. in [5] can improve the current policies, especially the feasibility ratio.

IntervalBranch&Bound(B, g, f ) {

While (B = ∅ and ffmin > ǫ obj ) { criterion := criterionChoice(LBBox, UBBox, UBProb)

[x] := bestBox (B, criterion);

in this paper by LBBox. It consists in selecting a node/box [x] with a minimum lower bound estimate of the objective function (lb([x])). The selected box [x] is then split into two sub-boxes [x] 1 and [x] 2 along one dimension (selected by another and more studied branching heuristic). Both sub-boxes are then handled by the Contract&Bound procedure.

A constraint x obj ≤ UBBox is first added to the system for decreasing the upperbound of the objective function in the box. This aims at finding a solution better than the current best feasible point. The procedure then contracts the handled box without loss of feasible part. In other words, some infeasible parts at the bounds of the domain are discarded by constraint programming (CP) and convexification algorithms. Since this contraction works on the extended box including the objective variable x obj , it may improve both bounds of the objective on the current box.

The last part of the procedure carries out upperbounding. FeasibleSearch calls one or several heuristics searching for a feasible point that improves the best cost found so far.

Note that the search tree is traversed in best-first order so that an exponential memory may be required to store the nodes to handle.

New Strategies Using Upper Bounds

In optimization, the selection of the next node to expand is crucial for obtaining a good performance. The best node we can choose is such that it will improve the most the upperbound. Indeed, the upperbound improvement reduces globally the feasible space due to the constraint: