Bernard Helffer
email: bernard.helffer@math.u-psud.fr

Rola Kiwan
email: rkiwan@aud.edu

Dirichlet eigenfunctions in the cube, sharpening the Courant nodal inequality

Keywords: 2010 Mathematics Subject Classification. 35B05, 35P20, 58J50. Keywords. Courant Theorem, Nodal lines, Nodal domains, Dirichlet, Cube

This paper is devoted to the refined analysis of Courant's theorem for the Dirichlet Laplacian in a bounded open set. Starting from the work by Å. Pleijel in 1956, many papers have investigated in which cases the inequality in Courant's theorem is an equality. All these results were established for open sets in R 2 or for surfaces like S 2 or T 2 . The aim of the current paper is to look for the case of the cube in R 3 . We will prove that the only eigenvalues of the Dirichlet Laplacian which are Courant sharp are the two first eigenvalues.

Introduction and main result

Consider the Dirichlet eigenvalues of the Laplacian in a bounded domain Ω with piecewise C 1 -boundary.

-∆u = λu in Ω u = 0 on ∂Ω.

We denote by {λ k } k≥1 = {λ k (Ω)} k≥1 the sequence of eigenvalues (counted with their multiplicity):

λ 1 < λ 2 ≤ λ 3 ≤ ... ≤ λ k ≤ ...
It is well known that the first eigenvalue is simple and that the eigenfunction u 1 has a constant sign in Ω. All the higher order eigenfunctions must change sign inside Ω and, consequently, must vanish inside Ω. We call nodal set of an eigenfunction u k associated with λ k the closure of the zero set of u k ,

N (u k ) = {x ∈ Ω; u k (x) = 0}.
This nodal set cuts the domain Ω \ N (u k ) into µ k = µ(u k ) connected components called "nodal domains". The famous Courant nodal theorem [START_REF] Courant | Ein allgemeiner Satz zur Theorie der Eigenfunktionen selbstadjungierter Differentialausdrücke[END_REF] of 1923 states that µ(u k ) ≤ k.

We will say that an eigenvalue λ is Courant sharp if λ = λ k and if there exists an associate eigenfunction with k nodal domains. While it is always true in the case of dimension 1 by the Sturm-Liouvillle theory, Pleijel's theorem [START_REF] Pleijel | Remarks on Courant's nodal theorem[END_REF] asserts in 1956 that equality can only occur for a finite set of k's, when the dimension is at least two.

Since we know that the first eigenfunction does not vanish and that the second eigenfunction has exactly two nodal domains, λ 1 and λ 2 are Courant sharp (µ 1 = 1 and µ 2 = 2). We are now interested in checking if other eigenvalues are Courant sharp.

Starting from the founding paper by Å. Pleijel [START_REF] Pleijel | Remarks on Courant's nodal theorem[END_REF], many papers (and some of them quite recent) have investigated in which cases this inequality is sharp: Helffer-Hoffmann-Ostenhof-Terracini [START_REF] Helffer | Nodal domains and spectral minimal partitions[END_REF][START_REF] Helffer | On spectral minimal partitions : the case of the sphere. Around the Research of Vladimir Maz'ya III[END_REF], Helffer-Hoffmann-Ostenhof [START_REF] Helffer | Minimal partitions for anisotropic tori[END_REF][START_REF] Helffer | A review on large k minimal spectral k-partitions and Pleijel's Theorem[END_REF], Bérard-Helffer [START_REF] Bérard | Dirichlet eigenfunctions of the square membrane: Courant's property, and A. Stern's and Å. Pleijel's analyses[END_REF][START_REF] Bérard | Stern's analysis of the nodal sets of some families of spherical harmonics revisited[END_REF][START_REF] Bérard | Courant sharp eigenvalues for the equilateral torus, and for the equilateral triangle[END_REF], Helffer-Persson-Sundqvist [START_REF] Helffer | Nodal domains in the square-the Neumann case[END_REF], Léna [START_REF] Léna | Courant-sharp eigenvalues of a two-dimensional torus[END_REF], Leydold [START_REF] Leydold | Knotenlinien und Knotengebiete von Eigenfunktionen[END_REF][START_REF] Leydold | On the number of nodal domains of spherical harmonics[END_REF][START_REF] Leydold | On the number of nodal domains of spherical harmonics[END_REF]. All these results were devoted to (2D)-cases in open sets in R 2 : squares, rectangles, equilateral triangles, disks,.. or in surfaces like S 2 or T 2 The aim of the current paper is to start the analysis of analogous results for domains in R 3 . Looking at explicitly solvable models and avoiding easy cases where the multiplicity of each eigenvalue is one (making the analysis of the nodal set easy), it is natural to first consider the case of the cube. More precisely, we will prove: Theorem 1.1. In the case of the cube C := (0, π) 3 the only eigenvalues of the Dirichlet Laplacian which are Courant sharp are the two first eigenvalues: λ 1 = 3 and λ 2 = 6.

Coming back to Pleijel's paper

Outside the proof of Pleijel's theorem in 2D, Pleijel [START_REF] Pleijel | Remarks on Courant's nodal theorem[END_REF] (see also [START_REF] Bérard | Dirichlet eigenfunctions of the square membrane: Courant's property, and A. Stern's and Å. Pleijel's analyses[END_REF] for a more detailed analysis) considers as an example the case of the square which reads Theorem 2.1. In the case of the square the only eigenvalues which are Courant sharp for the Dirichlet Laplacian are the two first eigenvalues and the fourth one.

The proof was based on a first reduction to the analysis of the eigenvalues less than 68 (the argument will be extended to the (3D)-case below and this is a quantitative version of the proof of Pleijel's theorem), then all the other eigenvalues were eliminated using this time a more direct consequence of Faber-Krahn's inequality, except three remaining cases for which Pleijel was rather sketchy which have to be treated by hand.

At the end of his celebrated paper Å. Pleijel wrote: " In order to treat, for instance the case of the free three-dimensional membrane [0, π] 3 , it would be necessary to use, in a special case, the theorem quoted in [8], p. 394 1 . This theorem which generalizes part of the Liouville-Rayleigh theorem for the string asserts that a linear combination, with constant coefficients, of the n first eigenfunctions can have at most n nodal domains. However, as far as I have been able to find there is no proof of this assertion in the literature." Pleijel was indeed speaking of a result presented in [8] as being proved in the thesis defended in 1932 at the University of Göttingen by Horst Herrmann (with Richard Courant as advisor). This so called Courant-Herrmann conjecture was asserting that, for a given k ∈ N, Courant's theorem holds also for linear combinations of eigenfunctions associated with eigenvalues λ j with j ≤ k. It is said in [START_REF] Gladwell | The Courant-Herrmann conjecture[END_REF] that the authors can not find any mention of the result in the thesis itself. We learned in a recent article by N. Kuznetsov [20] that V. Arnold [START_REF] Arnold | Topological properties of eigenoscillations in mathematical physics[END_REF] has discussed this question with R. Courant in the seventies and that he can prove that this theorem cannot be true in general.

Pleijel is not explicitly saying why he needs this result but one could think that he is interested, because he speaks about the "free problem" (i.e. the Neumann problem), in counting the number of components of the restriction of an eigenfunction to a face of the cube (0, π) 3 . Looking for example to the zeroset of (x, y, z) → a cos x cos y cos nz + b cos y cos z cos nx + c cos z cos x cos ny , one gets for fixed z = 0, a linear combination of the eigenfunctions of the square cos x cos y, cos y cos nx and cos x cos ny corresponding to two different eigenspaces for the Neumann Laplacian in the square (0, π) × (0, π). We will not go further in this paper on the Neumann problem but similar questions could also occur in the Dirichlet problem and we typically meet below the eigenfunction (x, y, z) → a sin x sin y sin nz + b sin y sin z sin nx + c sin z sin x sin ny , and will be interested for example in the intersection of its zero set with the hyperplace {z = π 2 } inside the cube (in the case n = 3).

Reminder on Pleijel's theorem in 3D

Let us first prove that there are only a finite number of eigenvalues that satisfy µ k := µ(u k ) = k. This proof was given in dimension n by Bérard-Meyer [START_REF] Bérard | Inégalités isopérimétriques et applications[END_REF].

Proposition 3.1. If λ k is an eigenvalue of (1) such that λ k-1 < λ k , and u k is an associated eigenfunction then:

λ 3 2 k |Ω| ≥ µ(u k ) 4 3 π 4 , (2)
where |Ω| denotes the volume of Ω.

Proof. Assume that the nodal set cuts the domain Ω in µ k connected components and let us denote them by Ω i ,

1 ≤ i ≤ µ k . Since u k does not vanish inside Ω i , it
is equal to its first eigenvalue and now using the (3D)-Faber-Krahn inequality on each component (see for example Bérard-Meyer [START_REF] Bérard | Inégalités isopérimétriques et applications[END_REF]):

λ 3 2 k |Ω i | ≥ 4 3 π 4 for 1 ≤ i ≤ µ k .
Adding together all the equations we get (2).

Theorem 3.2.

lim sup k→+∞ µ k k ≤ 9 2π 2 < 1 . (3)
In particular, there exists only a finite number of eigenvalues satisfying µ k = k.

Proof. We start from the Weyl asymptotics for the counting function

N (λ) := #{k , λ k < λ}, (4)
which reads

N (λ) ∼ 1 6π 2 |Ω|λ 3 2 . ( 5
)
For an eigenvalue

λ k such that λ k-1 < λ k , we have N (λ k ) = k -1 . Then from λ 3 2 k ∼ 6π 2 |Ω| k
together with (2), we get (3).

Remark 3.3. It is clear from (3) that we cannot have an infinite number of eigenvalues satisfying µ k = k.

The case of the cube

Let us consider the cube (0, π) 3 for which an orthogonal basis of eigenfunctions for the Dirichlet problem is given by:

u ,m,n (x, y, z) = sin( x) • sin(my) • sin(nz) , λ ,m,n = 2 + m 2 + n 2 ,
for , m, n ≥ 1. Applying Proposition 3.1 for this domain, we get Proposition 4.1. If u k is an eigenfunction associated with λ k such that u k has k nodal domains and if λ k-1 < λ k we have:

λ k 3 2 k ≥ 4 3 π . (6)
Here we will try to find a lower bound for the number N (λ), since we know the λ's are equal to 2 + m 2 + n 2 where , m, n are integers, so we need to count the number of the lattice points of R 3 inside the sphere of radius √ λ.

Lemma 4.2. If λ ≥ 3, then N (λ) > π 6 λ 3 2 - 3π 4 λ + 3 √ λ -2 -1 . ( 7
)
The proof is given in the appendix.

Lemma 4.3. If u k is an eigenfunction associated with λ k such that u k has k nodal domains we have:

3 4π - π 6 λ 3 2 + 3π 4 λ -3 √ λ + 3 > 0 . (8)
Proof. First by Courant theorem, we have necessarily λ k-1 < λ k . Applying [START_REF] Courant | Ein allgemeiner Satz zur Theorie der Eigenfunktionen selbstadjungierter Differentialausdrücke[END_REF], we have

k -1 = N (λ) > π 6 λ 3 2 - 3π 4 λ + 3 √ λ -2 -1 . i.e. k > π 6 λ 3 2 - 3π 4 λ + 3 √ λ -2 . (9)
Together with [START_REF] Bérard | Inégalités isopérimétriques et applications[END_REF], this implies:

3 4π - π 6 λ 3 2 + 3π 4 λ > 3 √ λ -2 .
One immediately sees that, for λ ≥ 3 ,

√ λ -2 - √ λ = - 2 √ λ + √ λ -2 ≥ - 2 1 + √ 3 > -1 .
Now setting µ = √ λ we get the third order inequality:

3 4π - π 6 µ 3 + 3π 4 µ 2 -3µ + 3 > 0 . (10)
Now, consider the cubic function h(x) = 3 4π -π 6 x 3 + 3π 4 x 2 -3x + 3 on (0, +∞). The coefficient of x 3 is strictly negative. One can check easily that h(7) < 0 , h (7) < 0 and that h (x) < 0 for all x > 7. This implies that h(x) is negative for x > 7. Hence Inequality [START_REF] Gladwell | The Courant-Herrmann conjecture[END_REF] for a positive µ implies µ ∈ (0, 7) and coming back to λ: λ < 49 . So we have finally proved: Proposition 4.4. If u k is an eigenfunction associated with λ k such that u k has k nodal domains and if λ k-1 < λ k , we have:

λ k ≤ 48 . ( 11
)

The list

In this section, we establish the list of the eigenvalues which are less than 48 and determine which of these eigenvalues satisfy the necessary condition (6) for being Courant sharp.

k ( , m, n) λ k λ 1 (1,1,1) λ 2 = λ 3 = λ 4 (1,1,2) λ 5 = λ 6 = λ 7 (1, 2, 2) λ 8 = λ 9 = λ 10 (1,1,3) λ 11
(2,2,2)

λ 12 = λ 13 = λ 14 = λ 15 = λ 16 = λ 17 (1,2,3) λ 18 = λ 19 = λ 20
(2,2,3) Table 1: List of all eigenvalues less than 49

λ 21 = λ 22 = λ 23 (1,1,4) λ 24 = λ 25 = λ 26 (1,3,3) λ 27 = λ 28 = λ 29 = λ 30 = λ 31 = λ 32 (1,2,4) λ 33 = λ 34 = λ 35 (2,3,3) λ 36 = λ 37 = λ 38 (2, 2, 4) λ 39 = λ 40 = λ 41 = λ 42 = λ 43 = λ 44 (1, 3, 4) λ 45 = λ 46 = λ 47 = λ 48 (3, 3, 3) & (1, 1, 5) λ 49 = λ 50 = λ 51 = λ 52 = λ 53 = λ 54 (2, 3, 4) λ 55 = λ 56 = λ 57 = λ 58 = λ 59 = λ 60 (1, 2, 5)
Coming back to the consequences of Faber-Krahn's inequality, one can check that among all the values on Table 1, the only eigenvalues that satisfy inequality (6) and λ k-1 < λ k are λ 1 , λ 2 , λ 5 , λ 8 and λ 12 .

Proposition 5.1. The only eigenvalues which can be "Courant sharp" are the eigenvalues λ k with k = 1, 2, 5, 8 and 12.

As λ 1 and λ 2 are Courant sharp, the only remaining cases to analyze correspond to k = 5, 8, 12.

In the next section we will by a finer analysis involving symmetries eliminate other cases.

Courant theorem with symmetry

We first recall some generalities which come back to Leydold [START_REF] Leydold | Knotenlinien und Knotengebiete von Eigenfunktionen[END_REF], and were used in various contexts [START_REF] Leydold | On the number of nodal domains of spherical harmonics[END_REF][START_REF] Leydold | On the number of nodal domains of spherical harmonics[END_REF][START_REF] Helffer | Nodal domains in the square-the Neumann case[END_REF][START_REF] Helffer | On spectral minimal partitions : the case of the sphere. Around the Research of Vladimir Maz'ya III[END_REF]. Suppose that there exists an isometry g such that g(Ω) = Ω and g 2 = Id . Then g acts naturally on L 2 (Ω) by gu(x) = u(g -1 x) , ∀x ∈ Ω , and one can naturally define an orthogonal decomposition of L 2 (Ω)

L 2 (Ω) = L 2 odd ⊕ L 2 even ,
where by definition L 2 odd = {u ∈ L 2 , gu = -u}, resp. L 2 even = {u ∈ L 2 , gu = u}. These two spaces are left invariant by the Laplacian and one can consider separately the spectrum of the two restrictions. Let us explain for the "odd case" what could be a Courant theorem with symmetry. If u is an eigenfunction in L 2 odd associated with λ, we see immediately that the nodal domains appear by pairs (exchanged by g) and following the proof of the standard Courant theorem we see that if λ = λ odd j for some j (that is the j-th eigenvalue in the odd space), then the number µ(u) of nodal domains of u satisfies µ(u) ≤ j. We get a similar result for the "even" case (but in this case a nodal domain D is either g-invariant or g(D) is a distinct nodal domain). These remarks may lead to improvements when each eigenspace has a specific symmetry. As we shall see, this will be the case for the cube with the map (x, y, z) → (π -x, π -y, π -z). We observe indeed that u ,m,n (π -x, π -y, π -z) = (-1) +m+n+1 u ,m,n (x, y) ,

and that 2 + m 2 + n 2 ≡ + m + n (mod.2) .
Hence, for a given eigenvalue the whole eigenspace is even if + m + n is odd and odd if +m+n is even. Equivalently, the whole eigenspace is even if the eigenvalue is odd and even if the eigenvalue is odd.

Application.

λ 5 is not Courant sharp. The eigenspace associated with λ 5 = 9 is even. This is the second one (in this even space). Hence it should have less than four nodal domains by Courant's theorem with symmetry and has labelling 5. λ 12 is not Courant sharp. λ 12 = 14 is the fifth eigenvalue in the odd space with respect to σ. It should has less than 10 nodal domains and has labelling 12.

7. The remaining value: k = 8 7.1. Main result. The proof of our main theorem relies now on the analysis of the last case which is the object of the next proposition. Proposition 7.1. In the eigenspace associated with λ 8 the eigenfunctions have either 2, 3 or 4 nodal domains. In particular λ 8 cannot be Courant sharp.

Preliminaries.

For the value λ 8 = 11 we have to analyze the zeroset of Φ a,b,c (x, y, z) := a sin x sin y sin 3z + b sin y sin z sin 3x + c sin z sin x sin 3y , for (a, b, c) = (0, 0, 0) . This looks nice because we can divide by sin x sin y sin z and by making the change of coordinates u = cos x, v = cos y, w = cos z, we get for the zero set of Φ a,b,c in the new coordinates a quadric surface Q a,b,c to analyze in the cube C = (-1, 1) 3 , whose equation is

(Q a,b,c ) 4 (au 2 + bv 2 + cw 2 ) -(a + b + c) = 0 , for (a, b, c) = (0, 0, 0) .
When a + b + c = 0, we immediately see that there are no critical points inside the cube, so the nodal set is simply an hypersurface (cylinder, ellipsoid or hyperboloid with one or two sheets). In this case, this is the analysis at the six faces of the cube which will be decisive for analyzing possible changes in the number of connected components. In the case when a + b + c = 0, we have a double cone with a unique critical point at (0, 0, 0). In the next subsections, we discuss the different cases.

Cylinder.

This corresponds to the case abc = 0 . We can use the (2D)-analysis as done in [START_REF] Bérard | Dirichlet eigenfunctions of the square membrane: Courant's property, and A. Stern's and Å. Pleijel's analyses[END_REF]. It is known that the number of nodal domains can only be 2,3 or 4 (See Section 3.1 and figure 2.1 there). See figure 1.

Double cone.

This corresponds to abc = 0 and a + b + c = 0 . The equation of Q a,b,c is:

au 2 + bv 2 = -cw 2 .
One can verify that the intersection of this cone with each horizontal side w = ±1 is exactly at the vertices of the cube u 2 = v 2 = 1, and that the intersection with each vertical face is a hyperbola. Therefore there are three connected components of C \ Q a,b,c . See figure 2 2 .

2 All the pictures in this paper were produced using MATLAB.

au 2 + bv 2 + (1 -a -b)w 2 = 1 4 . (12)
We denote by Ω a,b,c the open full ellipsoid delimited by Q a,b,c .

Let us look at the intersection of Q a,b,c with the horizontal faces. We have

au 2 + bv 2 ≤ - 3 4 + 2 3 < 0 .
We deduce that in this case there are no possible intersections with the horizontal faces, and therefore two subcases can occur depending on the intersection of Q a,b,c with the vertical edges. This set is determined by We note that this implies that Q a,b,c ∩ {w = 0} is an ellipse contained in the cube.

(1 -a -b)w 2 = 1 4 -(a + b) , w ∈ (-1, +1) . (13)
The equation of Q a,b,c can be written as:

au 2 + bv 2 = 1 4 -cw 2 .
Q a,b,c cuts R 3 into two components Ω + a,b,c and Ω - a,b,c where Ω + a,b,c contains (0, 0, 0). But we have to look inside the cube.

We first observe that Q a,b,c has empty intersection the vertical edges. We have indeed

a + b = 1 -c > 1 4 -c ≥ 1 4 -cw 2 .
We now look at the intersection with w = 0. We get an ellipse E 0 a,b,c := Q a,b,c ∩ {w = 0}, whose equation is

au 2 + bv 2 = 1 4 .
We observe that this ellipse could be included in the cube, if a > 1 4 or not if a ≤ 1 4 . We also look at the intersection with the upper horizontal face. We note that the ellipse E cut the cube in three components.

The three last subcases are under the condition that a > 1 4 . We note that this condition implies that E 0 a,b,c is strictly included in the square (-1, +1) × (-1, +1) and the discussion continues according to the position of E 1 a,b,c in the horizontal face.

Conclusion

In this paper we have analyzed the problem in the simplest example proposed by Pleijel. One can of course ask for similar questions for other geometries starting with the parallelepipeds, the ball (this case has been solved in [START_REF] Helffer | On nodal domains in Euclidean balls[END_REF]), the flat torus (C. Léna has announced the characterization in this case). It is probably difficult to find other solvable cases. The situation for (0, απ) × (0, βπ) × (0, γπ) is in principle easier in the "irrational" case when

α 2 + βm 2 + γn 2 = α 2 1 + βm 2 1 + γ 1 n 2 1 implies ( , m, n) = ( 1 , m 1 , n 1 )
. Each eigenvalue α 2 + βm 2 + γn 2 is indeed of multiplicity 1 and the corresponding eigenfunction has mn nodal domains. One can also think of analyzing "thin structures" (for example γ small or β + γ small, where previous results in lower dimension can probably be used) in the spirit of [START_REF] Helffer | Minimal partitions for anisotropic tori[END_REF] and get partial results. Another interesting question would be to analyze the Neumann problem for the cube in the spirit of [START_REF] Helffer | Nodal domains in the square-the Neumann case[END_REF].

A. The proof of Lemma 4.2

We follow an idea appearing in the (2D) case in a course of R. Laugesen [START_REF] Laugesen | Spectral theory of Partial Differential Equations[END_REF]. We start by assuming that λ is not an eigenvalue. With each triple ( , m, n) with ≥ 1, m ≥ 1, n ≥ 1, we associate the cube

Q ,m,n = ( -1, ) × (m -1, m) × (n -1, n) .
We observe that where we recall that |Ω| denotes the volume of Ω. We are interested in the lower bound. The claim of Laugesen is that

N (λ) =
N (λ) > |B λ | , (14)
where B λ := {(x + 1) 2 + (y + 1) 2 + (z + 1) 2 < λ, x > 0, y > 0, z > 0}.

The observation is that B λ ⊂ ∪ 2 +m 2 +n 2 <λ, ≥1,m≥1,n≥1 Q ,m,n .

For t > 0, [t] + denotes the smallest integer ≥ t. Let (x, y, z) ∈ B λ , then it is immediate to see that (x, y) ∈ Q Coming back to [START_REF] Helffer | On spectral minimal partitions : the case of the sphere. Around the Research of Vladimir Maz'ya III[END_REF], we have to find a lower bound for the area of B λ . We note that by translation by the vector (1, 1, 1):

|B λ | = |C λ | , (15)
where

C λ := D(0, √ λ) ∩ {x > 1} ∩ {y > 1} ∩ {z > 1} .
Let χ the characteristic function of the interval (0, 1). We have to compute the integral

|C λ | = D(0, √ λ)
(1 -χ(x))(1 -χ(y))(1 -χ(z)) dxdydz .

Developing the formula and using the symmetry by permutation of the variables, we get, if λ ≥ 3 , It is then immediate to get the lemma by observing that

|C λ | = D(0, √ λ) dxdydz (16)
D(0, √ λ) χ(x)χ(y)dxdydz > √ λ -2 . ( 17
)
We have assumed till now that λ was not an eigenvalue. But if λ is an eigenvalue > 3, we can apply the previous result for an increasing sequence λj such that λj → λ (where λj > 3 is not an eigenvalue). According to our definition of N (λ) in (4), we can pass to the limit and observing that in [START_REF] Helffer | On nodal domains in Euclidean balls[END_REF] the inequality is uniformly strict when applied to the sequence λ j , we keep the strict inequality when passing to the limit. The case λ = 3 can be verified directly.

  λ 61 = λ 62 = λ 63 = λ 64 = λ 65 = λ 66 (1, 4, 4) & (2, 2, 5) λ 67 = λ 68 = λ 69 (3, 3, 4) λ 70 = λ 71 = λ 72 = λ 73 = λ 74 = λ 75 (1, 3, 5) λ 76 = λ 77 = λ 78 (2, 4, 4) λ 79 = λ 80 = λ 81 = λ 82 = λ 83 = λ 84 = λ 85 = λ 86 = λ 87 (1, 1, 6) & (2, 3, 5) λ 88 = λ 89 = λ 90 = λ 91 = λ 92 = λ 93 = λ 94 = λ 95 = λ 96 (1, 2, 6) & (3, 4, 4) λ 97 = λ 98 = λ 99 = λ 100 = λ 101 = λ 102 (1, 4, 5) λ 103 = λ 104 = λ 105 (3, 3, 5) λ 106 = λ 107 = λ 108 (2, 2, 6) λ 109 = λ 110 = λ 111 = λ 112 = λ 113 = λ 114 (2, 4, 5) λ 115 = λ 116 = λ 117 = λ 118 = λ 119 = λ 120(1, 3, 6) λ 121[START_REF] Bérard | Stern's analysis of the nodal sets of some families of spherical harmonics revisited[END_REF][START_REF] Bérard | Stern's analysis of the nodal sets of some families of spherical harmonics revisited[END_REF][START_REF] Bérard | Stern's analysis of the nodal sets of some families of spherical harmonics revisited[END_REF]

  b, c) = (1, -1, 0) (c) (a, b, c) = (1, 0, 0)

Figure 1 :Figure 2 :

 12 Figure 1: Cylinders

See figure 3 . 4 -

 34 Subcase (a + b) > 1 4 . The ellipsoid Q a,b,c does not touch the vertical edges and in this case C ∩ Ω a,b,c c is connected and C \ Q a,b,c has exactly two connected components. Subcase (a + b) ≤ 1 4 . Q a,b,c cuts each vertical edge along a segment [-w 0 , +w 0 ] with w 0 = 1 (a + b) /(1 -a -b).The intersection of Q a,b,c with each vertical face of the cube is the union of two arcs of an ellipse. In this case it is clear that C \ Q a,b,c has three connected components.

  Figure 3: Ellipsoids

  1 a,b,c := Q a,b,c ∩ {w = 1} has always a non empty intersection with this face.

Figure 4 :Subcase a ≤ 1 4± 1 4-

 411 Figure 4: One Sheet Hyperboloid

7 . 7 . 4 -cw 2 .

 7742 ,c is contained in the horizontal face and Q a,b,c cuts the cube in two connected domains. ,c ∩ ∂C consists of two curves but Q a,b,c continue to cut the cube in two domains. For joining two points of Ω - a,b,c ∩∂C one can always go to a point in {w = 0} outside of E 0 a,b,c and use the connexity (inside the square C ∩ {w = 0}) of the complementary of the full ellipse. Subcase 3 4 ≤ a E 1 a,b,c ∩ C consists of four curves. Q a,b,c continue to cut the cube in two domains. Two sheets hyperboloid. This corresponds to abc = 0, a, b, c not of the same sign and (abc)(a + b + c) > 0 . We can assume b ≥ a > 0, c < 0 and a + b + c = -1 . The equation of Q a,b,c can be written as: au 2 + bv 2 = -1 The hyperplane {w = 0} is contained in one connected component. Hence looking at the symmetry w → -w, we get that necessarily an odd number (≥ 3) of nodal domains and ≤ 8 by Courant's theorem. Hence we know that it cannot be Courant sharp. More precisely, Q a,b,c meets the hyperplane {w = 1} along the ellipse E a,b,c which this time contains the horizontal upper face of the cube. The analysis of the intersection along each of the vertical faces (two symmetric curves by w → -w) shows that we always have exactly three connected components. See figure 5.

Figure 5 :

 5 Figure 5: Two Sheets Hyperboloid. (a, b, c) = (0.8, 0.8, -2.6)

2 +m 2

 2 +n 2 <λ, ≥1,m≥1,n≥1 |Q ,m,n | ≤

  [x]+,[y]+,[z]+ . It remains to verify that Q [x]+,[y]+,[z]+ ⊂ D(0, √ λ) , where D(0, ρ) denotes the open ball in R 3 of radius ρ. But we have, for (x, y, z) ∈ B λ , [x] 2 + + [y] 2 + + [z] 2 + ≤ (x + 1) 2 + (y + 1) 2 + (z + 1) 2 < λ .

  )χ(y)χ(z)dxdydz .

In the german version, this is p. 454 in the english version.


Romain Casati

Gilles Daviet

Alexandre Derouet-Jourdan

Geometry and Mechanics of fibers: Some numerical models (Invited Talk) Florence Bertails-Descoubes

Keywords: Physics-based simulation, thin elastic rod, frictional contact, hair simulation, inverse physics-based design

In this talk I will give an overview of our work on the simulation of fibers and entangled materials, such as hair, with a specific interest for virtual prototyping and computer graphics applications. I will first introduce a family of high-order, reduced models for discretizing Kirchhoff's equations for thin elastic rods in a both faithful and robust way. Such models are particularly well-suited for simulating inextensible fibers subject to bending and twisting, and featuring an arbitrary curly resting geometry. Then I will show how such models can be coupled to frictional contact using the nonsmooth contact dynamics framework, and I will present a hybrid iterative solver suitable for robustly handling thousands packed fibers at reasonable frame rates. Finally, I will give some insights into the inverse modeling of fibers, consisting in taking an arbitrary curve geometry as input and inferring corresponding geometric and physical parameters of the simulator such that the input geometry corresponds to a stable configuration at equilibrium.

Introduction

Deformable slender structures such as hair fibers, cloth, ribbons, tree branches or leaves, are ubiquitous around us. They often feature an intricate natural shape, ranging from straight to curly, and are characterized by a complex motion involving strongly nonlinear deformations, such as buckling. These complex shapes and motions greatly contribute to the visual richness of the real world. When multiple such structures are coupled together with contact and friction, the range of emerging phenomena is even more exacerbated, giving rise to stick-slip dynamical instabilities, entangling, or spontaneous collective behavior. Human hair, which is typically composed of 150,000 thin fibers, beautifully depicts such complex mechanical behavior when fluttering in the wind.

As the essence of Computer Graphics is to represent the visual appearance of the real world with the highest fidelity, it is important for practitioners to be able to capture all the relevant details of a given scene. In the case of a dynamical scene involving passive objects such as hair, cloth, or natural phenomena, physics-based simulation has proven over the years to be a method of choice for capturing resulting visual effects. Unlike phenomenological methods which develop descriptive models for reproducing a given emerging phenomenon, physics-based methods provide generative models whose goal is to explain the physical causes of the phenomenon. From a set of initial conditions as well as a few physical parameters (e.g., the mass, the stiffness, the natural shape), a physics-based simulator may generate not just a single effect, but a wide range of emerging phenomena revealing the whole complexity of the underlying physics.

When designing a physics-based simulator for computer graphics, one has inevitably in mind the four following criteria: 1. Realism: Ingredients which are necessary to capture relevant visual effects should be identified, and translated numerically with as few quality loss as possible. 2. Robustness: The simulator should converge properly for a subsequent range of parameters. 3. Efficiency: The simulator should be fast enough for allowing complex scenes to be simulated in reasonable timings (in our case, a few days of computation for a given scene is considered as an upper-bound). 4. Control: The simulator should provide the user with some handles to control the shape and motion of the object in an intuitive way.

Over the ten past years, we have been striving to develop some numerical models satisfying all four criteria at the same time. Our work has focused on the simulation of slender structures prone to contact and dry friction, and especially on the dynamics of thin fiber assemblies, with some direct applications in hair simulation and inverse physics-based design.

High-order reduced models for simulating dynamic fibers

The first part of the talk will be devoted to the presentation of the numerical models that we have been developing for simulating the dynamics of flexible fibers, the so-called super-helix and super-clothoid models [START_REF] Bertails | Super-helices for predicting the dynamics of natural hair[END_REF][START_REF] Bertails-Descoubes | Super-clothoids[END_REF][START_REF] Casati | Super space clothoids[END_REF].

I will first introduce the mechanical equations for inextensible thin elastic rods, namely the Kirchoff equations, which take the form of second order partial differential equations subject to boundary conditions. Noting that curvatures and twist play a major role both in the geometric and dynamic description of this model, we have come up with a spatial rod discretization based on elements that are polynomial in such quantities. Our first scheme relied on piecewise constant curvatures [START_REF] Bertails | Super-helices for predicting the dynamics of natural hair[END_REF], and was then extended to piecewise linear curvatures [START_REF] Bertails-Descoubes | Super-clothoids[END_REF][START_REF] Casati | Super space clothoids[END_REF]. One major advantage of such curvaturebased formulations is that the kinematics of the discretized rod remains, by construction, perfectly inextensible. Such intrisic inextensibility thus removes the burden of adding subsequent (stiff) inextensibility contraints when solving the dynamics. The price to pay, however, is that the geometry of the rod is not readily available but has to be computed iteratively from the curvatures.

In the piecewise constant case [START_REF] Bertails | Super-helices for predicting the dynamics of natural hair[END_REF], each element turns out to be a perfect circular helix (hence the super-helix name for the model), leading to a cheap and exact evaluation of all the terms of the discrete dynamic equations. For higher orders however, one unfortunately loses such a closed-form formula and a both accurate and efficient spatial integration scheme has to be designed.

In the piecewise linear case (where each element takes the form of a 3D clothoid), we were able to build an accurate integration scheme which proved to be orders of magnitude faster compared to classical integration methods [START_REF] Casati | Super space clothoids[END_REF]. The key of our approach was to leverage the fact that the solution takes the form of a power series expansion, while avoiding the pitfall of catastrophic cancellation through an adaptive integration strategy. With this tool in hand, we were able to demonstrate that the super-clothoid model could capture intricate shapes both robustly and efficiently, with better spatial accuracy and geometric fairness compared to state-ofthe-art methods (see Figure 1).

Robust frictional contact model for fiber assemblies

The second part of the talk will be focused on the dynamic simulation of fiber assemblies, where individual fibers are coupled to each other through contact and friction. Fig. 2: Simulation of a fast head movement without (top) and with (bottom) Coulomb friction, using our robust solver from [START_REF] Daviet | A hybrid iterative solver for robustly capturing Coulomb friction in hair dynamics[END_REF]. In the latter case, hair remains much more coherent and depicts typical stick-slip effects.

I will first illustrate why capturing threshold effects in friction is key to realism (see Figure 2), before introducing the nonsmooth Signorini-Coulomb friction model and its various formulations. We shall see what the numerical counter-parts are for each formulation, and how each of them performs in terms of efficiency, robustness, and scalability [START_REF] Bertails-Descoubes | A nonsmooth Newton solver for capturing exact Coulomb friction in fiber assemblies[END_REF].

Then I will explain how we managed to design a robust and scalable frictional contact solver by combining an iterative Gauss-Seidel strategy together with an extremely robust one-contact solver [START_REF] Daviet | A hybrid iterative solver for robustly capturing Coulomb friction in hair dynamics[END_REF]. Our global solver proved to converge well in scenarios involving thousands fibers subject to tens thousands frictional contact points, and thus allowed us to enhance considerably the realism of hair simulations. Our method has been adopted by the special effects industry for simulating hair and fur accurately [START_REF] Kaufman | Adaptive nonlinearity for collisions in complex rod assemblies[END_REF].

3 From geometry to physics: Inverse design of fiber assemblies Finally, I will present some new important challenges regarding inverse physicsbased design. While current simulators may succeed in reaching a good level of realism, they remain difficult to control in order to achieve a precise artistic goal. More precisely, to generate some desirable shapes and motions, one should be able to feed a simulator with the "right" parameters. Finding such parameters remains a difficult task, which is often performed through a tedious trial and error process. To make this task fully automatic, we have started looking at inverse solutions in the case where a static shape is provided as input: the inverse model should be able to interpret this shape automatically as a stable equilibrium of the simulator under gravity and other external forces, such as contact and friction.

In the case of an isolated fiber, we have shown that inverting any of our supermodel [START_REF] Bertails | Super-helices for predicting the dynamics of natural hair[END_REF][START_REF] Casati | Super space clothoids[END_REF] boils down to two decoupled problems, both of them being easy to solve [START_REF] Derouet-Jourdan | Stable inverse dynamic curves[END_REF][START_REF] Derouet-Jourdan | 3D Inverse Dynamic Modeling of Strands[END_REF]: First, an equilibrium condition which appears to be linear in the natural shape of the fiber, thanks to the curvature-based parameterization of our fiber models; Second, a sufficient stability condition that can be simply set by fixing a lowerbound for the ratio of stiffness over mass. Actually, the only remaining difficulty is to solve a merely geometric fitting problem -converting a curve as a piecewise helix or clothoid. In the case of helical fitting, we have already brought some efficient and robust solutions to this geometric approximation problem [START_REF] Derouet-Jourdan | Stable inverse dynamic curves[END_REF][START_REF] Derouet-Jourdan | Floating tangents for approximating spatial curves with G 1 piecewise helices[END_REF].

In the presence of contact and friction, Coulomb sticking constraints have to be considered, which makes the overall inverse problem nonsmooth and ill-posed. We have shown that assuming known mass and stiffness, as well as a simplified inverse model, it is possible to recover a plausible natural shape as well as frictional contact forces at play [START_REF] Derouet-Jourdan | Inverse dynamic hair modeling with frictional contact[END_REF]. This work allowed us, for the first time, to animate in a plausible way a few hair geometries stemming from recent hair captures (see Figure 4). [START_REF] Luo | Structure-aware hair capture[END_REF], inverted by our method in [START_REF] Derouet-Jourdan | Inverse dynamic hair modeling with frictional contact[END_REF] and physically animated (middle) and trimmed (right).

Conclusion

Throughout this long-term work on the numerical modeling of fibers and frictional contact, we have learnt that systematically concentrating the efforts on the upstream modeling and formulation of problems often pays off: even for very complex problems, the resulting numerics may be greatly simplified, and thus solved more easily and robustly. Keeping in mind this key lesson, we are starting to investigate the case of 2D slender structures (namely plates and shells), for which many exciting challenges remain open.

Fig. 1 :

 1 Fig. 1: Many physical strands exhibit a smooth curled geometry with linear-like curvature profile, which is captured and deformed accurately thanks to our super-clothoid model [4]. From left to right and top to bottom, three examples of real strands whose shapes are synthesized and virtually deformed in real-time using a very low number of 3D clothoidal elements: a vine tendril (4 elements), a hair ringlet (2 elements), and a curled paper ribbon (1 single element). Left photograph courtesy of Jon Sullivan, pdphoto.org.

Fig. 4 :

 4 Fig.4: Real curly wig (left) captured from[START_REF] Luo | Structure-aware hair capture[END_REF], inverted by our method in[START_REF] Derouet-Jourdan | Inverse dynamic hair modeling with frictional contact[END_REF] and physically animated (middle) and trimmed (right).


Alessandro Ferreira Leite

Azzedine Boukerche

Alba Cristina

Magalhaes Alves De Melo

Christine Eisenbeis

Claude Tadonki

Célia Ghedini Ralha

Power-Aware Server Consolidation for Federated Clouds

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


Julien Colin
email: julien.colin@centralesupelec.fr

Romain Rémond

Patrick Perre

P Perré

Design and optimization of industrial woody

Keywords: biomass pretreatment, wood energy, continuous dryers, multiscale modelling, design and optimization

INTRODUCTION

The bioeconomyi.e. the use of biomass for food, feed, industrial and energy purposes [START_REF]European Commission Innovating for sustainable growth: A bioeconomy for Europe[END_REF] has the potential to contribute to a more sustainable use of resources. In particular, the bioeconomy has a positive effect to reduce global warming, and lower the energy dependence of France. In this context, thermochemical conversion processes are promising for a wide range of applications such as heat, biomolecules or biofuels production [START_REF] Damartzis | Thermochemical conversion of biomass to second generation biofuels through integrated process design -A review[END_REF].

These dry conversion routes require increasingly high quality of the raw materials, particularly regarding moisture content (MC). The use of continuous convective dryers for the pretreatment of wood energy is attractive because of their low operating cost and their ease of integration in the production line.

However, the use of such dryers is not without pitfall. For example, heterogeneity of the climatic conditions inside the dryer, together with biomass variability, make the conception and optimization of these industrial plants complex. Identifying the best technical solutions and estimating their gains through empirical knowledge may be costly in both time and money. Simulation models may therefore be useful tools, either for optimal design or for on-line control of industrial plants. However, the various spatial scales involved in such processes, together with the huge variability of biomass, render such a modeling approach challenging.

At the dryer scale, the air parameters change along the flow. This profile of climatic conditions depends on the heat and mass transfer at the exchange surfaceswood particles and wall of the dryerand vice versa: the relative humidity (RH) and temperature of the air affect the drying rate of each particle. Therefore the two-way coupling between two spatial scales (particle ↔ dryer) must be understood in order to simulate the continuous drying of a bed of wood particles [START_REF] Perré | Multiscale aspects of heat and mass transfer during drying[END_REF][START_REF] Perré | A Dual scale computational model of kiln wood drying including single board and stack level simulation[END_REF]. Some numerical approaches use the concept of distributed micro-models. It allows the two-way coupling to be comprehensive even when all particles are different. Thus it is possible to assess the dryer performance for reducing the heterogeneity of final MC of particles, in spite of the initial variability.

Inside the biomass particle, Finite Element or Control Volume based computational models are very efficient in predicting the evolution of the variable fields (temperature, MC, pressure, strain, stress). They are also well adapted to deal with non-linearities and couplings, and for testing various formulations or parameter values. In this sense, they are able to provide predictive simulations, and can also assist the improvement of the fundamental understanding [START_REF] Patankar | Numerical heat transfer and fluid flow[END_REF][START_REF] Perré | Simulation par volumes finis des transferts couplés en milieux poreux anisotropes: séchage du bois à basse et à haute température[END_REF][START_REF] Prat | 2D modelling of drying of porous media: influence of edge effects at the interface[END_REF][START_REF] Turner | A two-dimensional orthotropic model for simulating wood drying processes[END_REF]. As drawbacks, this approach usually leads to significant computational times, which limits the total number of particles at the drier scale. In addition, particular phenomena such as the effect of checks on transfer are difficult to consider. That is why analytical models [START_REF] Chen | A discussion on a generalized correlation for drying rate modeling[END_REF][START_REF] Sander | Thin layer drying of porous mate-rials: selection of the appropriate model and relationships between thin-layer models parameters[END_REF][START_REF] Van Meel | Adiabatic convection batch drying with recirculation of air[END_REF] are often employed to simulate the drying of a large number of particles, as it is the case for the pretreatment of biomass. The concept of characteristic drying curve defines the drying rate as a proportion of the maximum drying rate which is obtained during the so-called constant drying rate (CDR) period, when the entire exchange surface is supplied with free water. This proportion is a dimensionless function of MC. However, a usual assumption in the approach proposed by van Meel [START_REF] Van Meel | Adiabatic convection batch drying with recirculation of air[END_REF] is that all of the energy supplied to the sample is used to evaporate water. This strong assumption leads to a poor agreement between the simulations and the experimental results during certain drying periods, namely during transient periods such as heating-up or during variable drying conditions [START_REF] Chen | A discussion on a generalized correlation for drying rate modeling[END_REF][START_REF] Saastamoinen | Drying of biomass particles in fixed and moving beds[END_REF]. Without considering the sensible heat, such models can neither faithfully simulate the evolution of wood MC if the airflow temperature varies, as it does along the tunnel of drying, nor predict the phenomena of condensation or heating-up. This paper proposes DryKiln_CRP, a multiscale computational model to design and optimize industrial convective dryers used for wood energy pretreatment. The physical formulation of this model at the particle and the dryer scale is presented in the following part. Then the relevance of the multiscale model as an operational tool will be illustrated via two case studies: the optimization of a fixed-bed dryer and, the design of a conveyor dryer. These dryers may involve simpleparallel-and counter-flowor combined airflow configurations. At this scale, we propose an extended van Meel model, able to account for the coupling between heat and mass transfer in an explicit way [START_REF] Colin | Séchage en continu du bois énergie comme moyen de préconditionnement en vue de sa conversion thermochimique : approches expérimentale et numérique[END_REF]:

AN INTERACTIVE, FLEXIBLE AND REALISTIC MULTISCALE MODEL

,0 ln v wood wood atm v m v wood surf atm v j A X m PP k C M A PP              (1)           , 0 , , 0 , h wood wood P wood P l surf Pl su v db wood rf j A m c X c T m X L T h T A c T T T                   (2)
In these equations, m 0 is the dry mass of the particle and C is the molar concentration in the boundary layer.

In equation ( 1), the vapor flux j v, wood at the exchange surface of the particle A wood is directly tied to the mass variations as a result of the mass balance over the entire particle. This flux is determined by the product of a mass transfer coefficient, the surface area of the particle and the driving force (film theory).

In equation ( 2), the heat flux j h, wood through the exchange surface of the particle A wood is directly tied to the sensible and latent heat variations of the moist wood due to the energy balance of the particle; the driving force is proportional to the difference of temperature between the air and the particle.

The vapor pressure at the surface of the particle is evaluated through the surface activity a*:

    *, surf v vs surf P a R P T H   (3)
Contrary to the water activity, which assumes a local thermodynamic equilibrium, the surface activity a* is a function allowing for the vapor pressure at the surface to be evaluated during drying through values averaged over the entire particle. Thus, the surface activity is no longer an equilibrium value but a kinetic parameter of the model.

During drying, the value of a* decreases from 1, when free water exists at the surface, to RH, when the particle is in equilibrium with the surrounding air. Following the concept of characteristic drying curve, the surface activity is a function of RH and ϕ, the dimensionless MC of wood:

1 cr eq cr cr eq if X X XX if X X XX            (4)
The main assumption of this model is that the temperature of the studied particles is uniform inside the particle:

surf TT  (5)
The assumption of temperature homogeneity is obviously relevant for small particles (small Biot number). It remains valid for large particle (characteristic dimension in the order of 10 cm) because the characteristic time of heat transfer is much lower than the one of mass transfer. This was reported several times in the literature for convective drying, either for measurements or simulation results. Indeed, during convective drying, the particle temperature evolves slowly as consequence of the evolution of surface moisture content. Due to this slow variation, the temperature field has enough time to diffuse inside the product, giving rise to a quasi-uniform temperature.

Equations ( 1), ( 2) and (3) become:

    0 ln *, atm v m v wood atm vs PP X m k C M A P a RH P T             (6)       0 , , 0 db wood P wood P l v h T T A m c X c T m X L T                (7)
This extended van Meel model being a characteristic drying curve based analytical model, all the internal mechanismsheat and mass transferare implicitly considered through the surface activity, a*.

In order to faithfully simulate the evolution of the particle parameters, a* must be previously identified by inverse analysis of drying kinetics of single particles under different drying conditions (temperature and RH).

These drying kinetics can be obtained by performing numerous experiments for any geometry or any wood species. However, this tedious task was reduced in our case by using the 3D version of TransPore [START_REF] Perré | A 3-D version of TransPore: a comprehensive heat and mass transfer computational model for simulating the drying of porous media[END_REF], a computational model simulating heat and mass transfer in porous media, to generate drying kinetics under a wide range of climatic conditions and for any particle size. This is however limited to the species for which TransPore is currently validated (beech and spruce), and for quite simple geometries.

Equations ( 6) and ( 7) are discretized in time using a Taylor series formulation, truncated at the first order, via a fully explicit scheme. The state of the system and the flux are assumed constant between times t and t+Δt.

From the particle to the dryer: a multiscale approach

The global model calculates the evolution of the air and wood parameters along the drying chamber using all fluxes exchanged between each particle. The system of equations describing the stack configuration [START_REF] Perré | A Dual scale computational model of kiln wood drying including single board and stack level simulation[END_REF] takes into account the two-way coupling between particles and dryer [START_REF] Colin | Séchage en continu du bois énergie comme moyen de préconditionnement en vue de sa conversion thermochimique : approches expérimentale et numérique[END_REF]

q j A j A x       (8)           , , , , , , . ., , , , , . .,
T q c q c A j j c T T x A j j c T T                      (9)
A lin.spe. refers to a linear specific areain square meter of exchange surface of wood/air or wall/air per linear meter of kiln along its main axis i.e. the direction of biomass displacement.

In this model, it is assumed that the climatic conditions are homogeneous in any cross-section of the dryer. 8) and ( 9) are discretized in space using a Taylor series formulation, truncated at the first order, via a fully explicit scheme. The dryer is then segmented along its main axis into N Control Volumes (CV), each containing M particles (Fig. 2). For the k th CV, the state of the system and the flux are assumed constant between the positions (k-1)•Δx and k•Δx.

position along the dryer k•Δx (k-1)•Δx CV k , in h wall j , in v wall j , out h wall j , h wood j         , , , 1
1 in db m v m a T k x q k x q            , , , in db m v m a T k x q k x q      , v wood j , out db T  wall conveyor belt
The climatic conditions along the dryer are evaluated from the local heat and mass fluxes simulated for each particle. Thus, the global model is connected to thousands of local models (extended van Meel models), which emphasizes the need for a low computational time for the simulation of each particle behavior during the drying.

In order to improve its predictive potential, the global model takes into account the heat losses as well as the possible condensation at the kiln wall.

The computational solution of the multiscale model is encoded in a custom software -DryKiln_CRPwritten in Fortran 90.

Particular attention is paid to the structure of the model so that it is suitable for a wide range of continuous dryers: it can handle parallel-and counter-flow configurations, several air injection or exhaust points and air recirculation.

This new multiscale model is then validated by comparing simulations and experimental kinetics collected using an original device [START_REF] Colin | Séchage en continu du bois énergie comme moyen de préconditionnement en vue de sa conversion thermochimique : approches expérimentale et numérique[END_REF].

With the help of twofictitious but relevantcase studies, the following sections aim at demonstrating the potential of DryKiln_CRP as an engineering tool.

OPTIMIZATION OF A FIXED BED DRYER

Description of the studied case

A company needs to produce one ton per hour (anhydrous mass flow-rate) of dry beech with an average MC less than 20 % and a heterogeneity of MC lower than 10 %.

After harvesting, woody biomass is ground into large wood chips (60×40×20 mm). The macroporosity of the resulting granular medium is equal to 0.6. It is then stored in bulk to let it dry in open air for a few weeks, bringing its MC down to an average of 60 % with a standard deviation of 5 %. Finally the drying process is finished in an industrial plant (Fig. 3). The company has a fixed-bed dryer, 4 m high and 3 m in diameter. The wood goes through the vertical drying tower from top to bottom, and the air from bottom to top: this is a counterflow configuration. The dryer is supplied with warm air through a waste heat source from a nearby industrial facility. The inlet temperature of the air is 60 °C, its RH is 4.3 %, and its volume flow-rate is 5 m 3 •s -1 . The value of the heat transfer coefficient (49 W•m -2 •K -1 ) is calculated thanks to the formula proposed by Nganhou [START_REF] Nganhou | Heat and mass transfer through a thick bed of cocoa beans during drying[END_REF] and the value of the mass transfer coefficient (0.046 m•s -1 ) is obtained by transfer analogy [START_REF] Incropera | Fundamentals of heat and mass transfer[END_REF][START_REF] Salin | External heat and mass transfert[END_REF]. For these climatic conditions and this geometry of particles, the critical moisture contentidentified from simulated kinetics (TransPore 3D)is equal to 103 %.

wet wood chips (beech): = 60 % 5 % = 15 C dry wood chips; objectives:

< 20 %  Δ < 10 %  q m0, wood = 1 t•h -1 warm air:  T ∞ = 60 C  RH = 4.3 %  q v = 5 m 3 •s -1 cold air 4 m 3 m X Fig. 3

. Schematic description of the fixed bed dryer

The wall consists of a single steel sheet (Table 1).

ρ (kg•m -3 ) c P (J•kg -1 •K -1 ) λ (W•m -1 •K -1 ) 1 7500 450 45
The climatic conditions outside of the dryer are as follows: air temperature = 15 °C, RH = 50 % and heat transfer

coefficient = 8 W•m -2 •K -1 .
The initial drying process does not therefore meet the specifications required. The multiscale model is then used test different ways which may be employed to overcome this problem, and finally to propose adequate technical solutions.

The initial configuration

The flow-rate is set to meet the company needs: one ton of dry mass per hour.

m•h -1 (311 kg•h -1 )
With a wood flow-rate equal to 311 kg•h -1 , the average MC of the biomass reaches 20 % at the outlet of the dryer, and the dispersion of these values is about 2.7 % (Fig. 5). This configuration achieves the qualitative objectives. However, from a quantitative point of view, the mass flow-rate of dry wood is only one-third of the requirement.

Thermal insulation of the dryer wall

Now we thermally insulate the drying tower by covering its outer face with a 110 mm thick layer made of glass wool in order to decrease heat loss through the wall (Table 2). The saved energy is expected to improve the plant performance.

Table 2. Characteristics of the insulated drying tower wall The energy yield is defined as the share of the heat source power used for evaporating the water from the wood.

thickness (cm) ρ (kg•m -3 ) c P (J•kg -1 •K -1 ) λ (W•m -1 •K -1 )
DryKiln_CRP is a useful tool to test the performance of the wall insulation. With a wood flow-rate equal to one ton per hour, the gain is perceptible but remains low (Fig. 6). Regarding the outlet MC of the wood, the difference between the initial configuration and the current one is lower than 0.5 % (Table 3). With the insulated dryer, we need to decrease the wood flow-rate to 325 kg•h -1 to reach an average MC of 20 %. The increase of productivity obtained thanks to the thermal insulation of the wall is not negligible (14 kg•h -1 ). Because of its low cost and its long-term benefit, thermal insulation of the wall is implemented. This however is clearly far not enough to meet the company requirements, and it must be coupled with other technical solutions.

Reduction of the particle size

The analysis of Fig. 4-b tells us that the air is not fully saturated with moisture at the outlet of the dryer. Thus the evaporation capacity of the fluid is not totally exploited. The resistance to internal transfer within large wood chips certainly limits the drying rate. In order to test this hypothesis, we examine the drying of smaller wood chips: 20×20×5 mm. This reduction of the particle size affects the surface activity and the critical moisture content rises up to 107 %. The quality requirements limit the flow-rate of small wood chips to 639 kg•h -1 (Fig. 7-a). Of course this is lower than the target flow-rate (1 t•h -1 ) but is almost twice as high as the maximum flux previously identified for the large particles (325 kg•h -1 ). We can notice that the MC decreases very slowly for the two first meters of the dryer. The air is indeed almost saturated in this part of the dryer (Fig. 7-b), and the tower could be smaller in height without affecting the quality of drying. With the initial size of wood chips, the drying rate was indeed controlled by the resistance to internal transfer. It is then not possible to exploit the full evaporation potential of the air (Table 4). In the case of the configuration with a wood mass flow-rate optimized for large particles (325 kg•h -1 ), the air at the outlet remains warm (41.4 °C) and far from saturation (24.4 %). Thus it is not surprising that this open system (no air recirculation device) has a very low energy yield (37.4 %). The reduction of the wood chips size is therefore a simple and efficient solution to optimize the drying performance. However, before implementing this solution, the company must ensure that this particle size is compatible with the downstream processes. Furthermore, the additional cost of grinding should be evaluated ex ante.

Additional heat source

Whatever the particle size, an additional heat source is needed to reach the target wood flow-rate

(1 t•h -1 )
through increasing the evaporation potential of air. DryKiln_CRP is then useful to define the required power.

Two parameters of the heat source can be adjusted: flow-rate and temperature. For the drying of small wood chips (20×20×5 mm), we increase the air flow-rate at the same temperature (60 °C). This flow-rate must reach 7.8 m 3 •s -1to be compared to the initial value of 5 m 3 •s -1in order to meet all the specifications. The corresponding additional power required is 136 kW (Fig. 8). Compared to the drying of small particles, the temperature of the outlet air is higher in this caseat respectively 24.0 °C and 45.3 °C -and its relative humidity is farther from saturationat respectively 99.4 % and 32.2 % (Fig. 10). Thus the energy yield of small particles is higher than in the case of large particlesat respectively 72.5 % and 49.9 %. The results of simulation allow us also to estimate the energy efficiency, i.e. the specific energy supplied by ton of evaporated water by differentiating the waste and the additional heat sources (Table 5). In both case, the configurations achieve the qualitative requirements for heterogeneity of MC: the dispersion of values is respectively about 3.3 % and 2.8 %.

On the one hand, Fig. 8 and Fig. 10 confirm that internal transfers are intensified by the grinding of the wood into small particles. In this case, because of the moisture-saturated air, the two first meters from the top of the dryer are useless. The wood flow-rate is driven by evaporation potential of air. On the other hand, the drying performance of large particles remains limited by internal transfers (Fig. 9 and Fig. 10 ). The wood flow-rate is now driven by the drying kinetic of the single particle. There is therefore an intermediate particle size that enables to exploit the full evaporation potential of the air while using the entire height of the dryer. The drying performance is then comparable to the value obtained for small particles, but the energy required for the grinding of the wood is lower. Thus, the particle size is a key parameter for the energy and economic optimization of the whole production line.

Whether by determining the power of the additional heat source needed, or by calculating the energy performance of each configuration, DryKiln_CRP is a powerful tool for the improvement of existing dryers.

DESIGN OF A CONVEYOR DRYER

Description of the studied case

A company needs to produce 450 kg•h -1 (anhydrous flow-rate) of dry wood for the production of charcoal.

Before carbonization, the mean and the heterogeneity of the wood MC must be lower than respectively 15 % and 10 %. The woody biomass is composed of eucalyptus logs (assumed to be cylindrical, Table 6): The company wants to build an industrial continuous dryer to pretreat this wet raw material. In order to reduce the operating costs, the dryer will be supplied with warm air through a waste heat source from a nearby industrial facility: the inlet temperature of the air is 120 °C, its RH is 0.59 %, and its volume flow-rate is 4.5 m 3 •s -1 . The use of a conveyor dryer is appropriate due to the large particle size. The conveyor belt carries a 50 cm thick bed of logswhich is thin compared to the particle sizethrough a tunnel. The cross-section inside this tunnel is 1 m high and 3 m wide. Considering this bed configuration and the heat source characteristics, the values of the heat and mass transfer coefficients are respectively equal to 29 W•m -2 •K -1 and 0.032 m•s -1 [START_REF] Nganhou | Heat and mass transfer through a thick bed of cocoa beans during drying[END_REF][START_REF] Incropera | Fundamentals of heat and mass transfer[END_REF][START_REF] Salin | External heat and mass transfert[END_REF]. For these climatic conditions and this geometry of particles, the critical moisture contentidentified from experimental kineticsequals to 78.0 %. The wall consists of a 200 mm thick inner layer of concrete and a 100 mm thick outer layer of glass wool.

By using the multiscale model, the best airflow configuration will be identified. The objective is to meet the specifications while designing a dryer as short as possible in order to reduce the investment costs. The c-f configuration is the most efficient (Table 7): the tunnel is 10 m long, i.e. less than half the length for the case of the p-f configuration (20.75 m). In light of the previous results, it seems possible to reduce the length of the dryer if we succeed in obtaining a high rate of drying at both the beginning and the end of the tunnel. To this end, air is injected at each end of the dryer with equal flow-rates, and the gas exhaust is located at the mid-length of the tunnel. This constitutes a p-f + c-f configuration. The opposite configuration (c-f + p-f) is also studied (Fig. 12). 13 shows that the expected gain is not obtained. The combined configurations do not significantly reduce the length of the dryer compared to the optimal length determined for each simple configuration. Moreover the p-f + c-f configuration is less efficient than the c-f + p-f. These results are against the intuition that led to the study of the combined configurations; they confirm the usefulness of the model for validating or invalidating preconceived ideas and at a negligible cost compared with the building of a pilot plan.

Combined airflow configurations

DryKiln_CRP is also useful for studying the impact of an unequal distribution of the air flow-rate within the c-f and p-f sections, or the effect of the displacement of the limit between these sections along the dryer. In both cases, the results are similar to the previous observation: sharing the heat source does not appear to be a promising solution, and the c-f configuration remains the most efficient.

Air recirculation system

, , 0 A recirculation system can be implemented in order to exploit the full evaporation potential of the air (Fig. 14).

m a q , ,1 m a q , , 2 m a q , , 2 m a q   , , 2 (1
The air flow-rate inside the tunnel increases with the recirculation rate β whose value can vary between 0 and 1.

The heat source parameters are unchanged. Moreover, we assume that there is neither heat loss through the recirculation pipe nor heat gain because of the recycling fan. The heat and mass balances at the mixing point allow us to calculate the injected air parameters from the heat source and exhausted air parameters, and the recirculation rate: Although the air recirculation cannot be used to make the tunnel shorter, it reduces the risk of fire because the temperature of the air injected into the tunnel decreases as the recirculation rate increases (from 120 °C when β = 0 to 78.9 °C when β = 2/3).

, ,1 , , 2 , , 0 1 1 m a m a m a q q q      (10) , ,1 , , 0 , , 2 m v m v m v q q q     (11)     , , , 0 2 0 , , , 0 0 , , 2 2 1 , , , 0 , , , 0 , , 2 1 1 1
P a m a P v m v m v P a m a P v m v m v c q T T c q T q T T c q c q q                            (

Optimization of the MC heterogeneity

The c-f configuration meets the objective in terms of average MC (15 %) but the heterogeneity remains too high (14 %, Fig. 16). In order to obtain a MC heterogeneity lower than 10 %, we must extend the tunnel by 4.5 m, giving a total length of 14.5 m. The average MC of wood is then 5.3 %. An air recirculation rate equal to 2/3 rds allows the size of the dryer with a p-f configuration to be limited to 14.75 m (Fig. 15). In this case, the MC heterogeneity is about 9 %. Both dryers are similar regarding their length but the c-f configuration -14.5 m long does not need an air recirculation system. Thus the c-f configuration appears as the most suitable design because of its lower investment and operating costs.

CONCLUSION

The model DryKiln_CRP presented in this paper, based on an extended van Meel model at the particle level, is able to faithfully simulate the continuous drying of wood for energy purposes. As main feature, the multiscale structure of the model allows the coupling between the different spatial scales to be accounted for with a large number of particles (up to several thousands) with different properties. Thus, it is able to account for average fields, standard deviations and local values. Beyond its application to wood, this modelwith its multiscale structurecan be applied to any other macroscopic granular medium that presents variability and for which the concept of characteristic drying curves is applicable. Finally, at the drier scale, the calculations of energy efficiency take into account the heat flows through the wall as well as any possible condensation that may occur at its surface.

Both case studies have demonstrated the potential of DryKiln_CRP as an engineering tool to design or optimize industrial facilities. Indeed, after having identified the surface activity of the product to dryfrom drying kinetics obtained by labscale experiments or TransPore 3D simulations -, DryKiln_CRP is suitable for a wide range of configurations thanks to its ability to simulate parallel-and counter-flow configurations, several air injection or exhaust points, air recirculation and the wall thermal insulation of the wall. Such a tool allows us to reduce the investment and operating costs while saving energy while improving the quality of drying.

Fig. 1 .

 1 Fig. 1. Heat and mass fluxes at the surface of the particle

Fig. 2 .

 2 Fig. 2. Heat and mass fluxes inside a Control Volume for a conveyor dryer

Fig. 4 .Fig. 5 .

 45 Fig. 4. Simulated evolution of the wood particles MC (a) and of the climatic conditions (b) along the dryer for the initial configuration

Fig. 6 .

 6 Fig. 6. Simulation of the drying quality improvement thanks to the thermal insulation of the wall

Fig. 7 .

 7 Fig. 7. Simulated evolution of the MC of small wood chips (a) and of the climatic conditions (b) along the dryer

Fig. 8 .

 8 Fig. 8. Simulated evolution of the MC of small wood chips along the dryer supplied with an additional heat source (136 kW)

Fig. 9 .Fig. 10 .

 910 Fig. 9. Simulated evolution of the MC of large wood chips along the dryer supplied with an additional heat source (318 kW)If the downstream processes are not compatible with a fine grinding, an increase in the air temperature has two positive effects; a better evaporation potential and a thermal activation of internal transfer which is very useful for the drying of large particles (60×40×20 mm). With an air flow-rate equal to 7.8 m 3 •s -1 , the temperature should reach 87 °C to ensure the quality of the dry wood. The additional heating power is then 318 kW (Fig.9). 0 20 40 60 80 100 0 1 2 3 4

Fig. 11 .Table 7 .

 117 Fig. 11. Simple airflow configurations: parallel-flow, (p-f) and counter-flow (c-f)

Fig. 12 .

 12 Fig. 12. Combined airflow configurations: parallel-flow followed by counter-flow (p-f + c-f) and counter-flow followed by parallel-flow (c-f + p-f)

Fig. 13 .

 13 Fig. 13. Simulated evolution of the average MC of the bed along the tunnel for simple and combined configurations

Fig. 14 .

 14 Fig. 14. Schematic diagram of the air recirculation system

Fig. 15 .

 15 Fig. 15. Influence of the air recirculation rate on the length of the tunnel with p-f and c-f configurations

Fig. 16 .

 16 Fig. 16. Unsatisfying decrease of the MC heterogeneity with the c-f configuration

Table 1 .

 1 Characteristics of the drying tower wall

	thickness
	(cm)

Table 3 .

 3 Influence of the wall insulation on the energy performance of the dryer

	insulated wall	outlet MC (%)	exhausted air T ∞ RH (°C) (%)	yield (%)	heat loss (%)
	without	40.4	27.0 69.6 55.8 2.66
	with	40.0	27.6 68.0 56.9 0.17

Table 4 .

 4 Summary of qualitative, quantitative and energy aspects of the studied configurations

	wood chips	flow-rate (kg•h -1 )	outlet MC (%)	outlet air T ∞ RH	yield (%)
				(°C)	(%)
	large	1000	40.0	27.6	68.0	56.9
		325	20.0	41.4	24.4	37.4
	small	1000	36.1	23.2	100.0	67.9
		639	20.0	24.0	99.9	72.7

Table 5 .

 5 Energy efficiency of the drying with additional heat sources

	wood	global energy	from waste heat	from additional heat
	chips	efficiency (kWh/kg)	source (kWh/kg)	source (kWh/kg)
	small	0.95	0.61	0.34
	large	1.41	0.80	0.61

Table 6 .

 6 Characteristics of eucalyptus logs to dry

	parameter	mean	standard deviation
	length (mm)	300	20
	diameter (mm)	85	5
	density (kg•m -3 )	600	/
	initial MC (%)	90	10
	temperature (°C)	20	/
	bed macroporosity	2/3	/