Caroline Kalla
email: kalla@crm.umontreal.ca

Christian Klein
email: christian.klein@u-bourgogne.fr

COMPUTATION OF THE TOPOLOGICAL TYPE OF A REAL RIEMANN SURFACE

published or not. The documents may come

Introduction

Riemann surfaces have many applications in physics and mathematics as in topological field theories and in the theory of integrable partial differential equations (PDEs). In concrete applications such as solutions of PDEs, e.g. Korteweg-de Vries and nonlinear Schrödinger (NLS) equations, see e.g. [START_REF] Belokolos | Algebro-geometric approach to nonlinear integrable equations[END_REF] and references therein, physical quantities as for instance the amplitude of a water wave are real. Thus reality conditions on the solutions are important in practice. The corresponding solutions have to be constructed on real Riemann surfaces, i.e., surfaces with an anti-holomorphic involution τ acting as the complex conjugation on a local parameter on the surface. Regularity conditions for these solutions depend on the topological type of the surface, i.e., whether there are connected sets of fixed points of the involution τ , the real ovals, and whether these ovals separate the surface into two connected components.

It is well known that all compact Riemann surfaces can be realized via nonsingular algebraic curves in P 2 (C). F. Klein [START_REF] Klein | On Riemann's theory of algebraic functions and their integrals[END_REF] observed that a real Riemann surface can be obtained in an analogous way from a nonsingular real plane algebraic curve R with an affine part of the form

(1.1) f (x, y) = N n=1 M m=1
a mn x m y n = 0, x, y ∈ C, a mn ∈ R.

The focus of this paper is on real compact Riemann surfaces. For curves of the form (1.1) the action of the complex conjugation gives rise to an antiholomorphic involution τ defined on R by τ (x, y) = (x, ȳ). The set of fixed points of τ is denoted by R(R) and is called the real part of R. The connected components of R(R) are called real ovals. Historically, the first result in the topology of real algebraic curves was obtained by Harnack [START_REF] Harnack | Ueber die Vieltheiligkeit der ebenen algebraischen Curven[END_REF]: the number k of real ovals of a curve R of genus g satisfies 0 ≤ k ≤ g + 1. In other words, the Date: April 21, 2012. We thank V. Shramchenko for useful discussions and hints. This work has been supported in part by the project FroM-PDE funded by the European Research Council through the Advanced Investigator Grant Scheme, and the ANR via the program ANR-09-BLAN-0117-01. number of connected components of the real part of a real nonsingular plane algebraic curve cannot exceed g + 1. Curves with the maximal number of real ovals are called M-curves.

The complement R \ R(R) has either one or two connected components: if R \ R(R) has two components, the curve R is called a dividing curve, otherwise it is called non-dividing (notice that an M-curve is always a dividing curve). The topological type of R is usually denoted by (g, k, a) where g denotes the genus, k the number of real ovals, and a = 0 if the curve is dividing, a = 1 if it is non-dividing. This implies that the topological type of a curve without real oval is (g, 0, 1). Notice that the first part of Hilbert's 16th problem is concerned with the relative configuration of real ovals of a plane algebraic curve of given degree in P 2 (R), i.e., how many ovals can lie in interior of another oval. This question has been studied by many authors, see, for instance, [START_REF] Hilbert | Über die reellen Züge algebraischen Kurven[END_REF][START_REF] Rohn | Die Maximalzahl und Anordnung der Ovale bei der ebenen Kurve 6. Ordnung und bei der Fläche 4. Ordnung[END_REF][START_REF] Petrovsky | On the topology of real plane algebraic curves[END_REF][START_REF] Gudkov | Complete topological classification of the disposition of ovals of a sixth order curve in the projective plane[END_REF][START_REF] Ya | Gluing of plane real algebraic curves and constructions of curves of degrees 6 and 7[END_REF], 1] and references therein. However, until now the complete answer is known only for curves of degree 7 and less. We will not discuss this topic here, but we would like to mention that in general, a solution to this problem, namely, the knowledge of the embedding R(R) ⊂ P 2 (R), does not provide any information on the embedding R(R) ⊂ R which is the subject of the present paper. For instance, it is possible to construct two real plane algebraic curves having the same degree and the same configuration of ovals, one of them being dividing and the other non-dividing (see [START_REF] Gabard | Sur la topologie et la géométrie des courbes algébriques réelles[END_REF] p. 8).

The aim of this paper is twofold: to determine the topological type (g, k, a) of a real algebraic curve of the form (1.1) with a numerical approach, and to transform periods of the holomorphic differentials of the curve to a form where the A-periods are real. There exist various algorithms which give the oval arrangements of a given real algebraic curve, see, for instance, [START_REF] Arnon | A polynomial time algorithm for the topological type of a real algebraic curve[END_REF][START_REF] Coste | Thom's lemma, the coding of real algebraic numbers and the computation of the topology of semi-algebraic sets[END_REF][START_REF] Sakkalis | The topological configuration of a real algebraic curve[END_REF][START_REF] Feng | Decomposition and Computation of the Topology of Plane Real Algebraic Curves[END_REF][START_REF] Gonzalez-Vega | Efficient topology determination of implicitly defined algebraic plane curves[END_REF][START_REF] Hong | An efficient method for analyzing the topology of plane real algebraic curves[END_REF][START_REF] Seidel | On the exact computation of the topology of real algebraic curves[END_REF], all of them following the same scheme. But to the best of our knowledge, there exists no algorithm that computes the parameter a in the topological type (g, k, a) of R, which encodes the property of the curve to be dividing or not. The starting point of our algorithm is the work by Deconinck and van Hoeij who developed an approach to the symbolic-numerical treatment of algebraic curves. This approach is distributed as the algcurves package with Maple, see [START_REF] Deconinck | Computing Riemann matrices of algebraic curves[END_REF][START_REF] Deconinck | Computing Riemann Theta Functions[END_REF][START_REF] Deconinck | Computing with plane algebraic curves and Riemann surfaces: the algorithms of the Maple package "algcurves[END_REF]. A purely numerical approach to real hyperelliptic Riemann surfaces was given in [START_REF] Frauendiener | Hyperelliptic theta functions and spectral methods[END_REF][START_REF] Frauendiener | Hyperelliptic theta functions and spectral methods: KdV and KP solutions[END_REF], and for general Riemann surfaces in [START_REF] Frauendiener | Algebraic curves and Riemann surfaces in Matlab[END_REF]. For a review on computational approaches to Riemann surfaces the reader is referred to [START_REF]Computational Approach to Riemann Surfaces[END_REF].

The codes [START_REF] Deconinck | Computing with plane algebraic curves and Riemann surfaces: the algorithms of the Maple package "algcurves[END_REF][START_REF] Frauendiener | Algebraic curves and Riemann surfaces in Matlab[END_REF] compute the periods of a Riemann surface in a homology basis which is determined by an algorithm due to Tretkoff and Tretkoff [36]. This homology basis is in general not adapted to possible symmetries of the curve (as the involution τ of real curves). It means that the action on the computed homology basis of any automorphism of the curve cannot be expressed in a simple way in terms of this basis. However, the choice of a basis, where certain cycles are invariant under the automorphisms, is often convenient in applications. In the context of solutions to integrable PDEs on general compact Riemann surfaces as for the Kadomtsev-Petviashvili (KP) (see [START_REF] Dubrovin | Real theta-function solutions of the Kadomtsev-Petviashvili equation[END_REF]) and the Davey-Stewartson (DS) equations (see [START_REF] Malanyuk | Finite-gap solutions of the Davey-Stewartson equations[END_REF][START_REF] Kalla | New degeneration of Fay's identity and its application to integrable systems[END_REF]), smoothness conditions are formulated conveniently in a homology basis adapted to the anti-holomorphic involution τ defined on the surface. For instance, on a real surface there exists a canonical homology basis (A, B) (that we call for simplicity symmetric homology basis in the following) satisfying the conditions (see [START_REF] Seppäłä | Moduli Spaces for Real Algebraic Curves and Real Abelian Varieties[END_REF][START_REF] Vinnikov | Self-adjoint determinantal representations of real plane curves[END_REF])

τ A τ B = I g 0 H -I g A B ,
where H ∈ M g (Z/2Z) is a g × g matrix which depends on the topological type (g, k, a) (see section 2); here I g denotes the g × g unit matrix.

In [START_REF] Kalla | On the numerical evaluation of algebro-geometric solutions to integrable equations[END_REF] we studied for the first time numerically solutions to integrable equations from the family of NLS equations, namely, the multi-component nonlinear Schrödinger equation and the (2 + 1)-dimensional DS equations. A symplectic transformation of the computed homology basis to the symmetric homology basis was introduced in [START_REF] Kalla | On the numerical evaluation of algebro-geometric solutions to integrable equations[END_REF] and was constructed explicitly for concrete examples. In the present work we give in Theorem 3.3 a complete description of such a symplectic transformation depending on the topological type of the underlying real algebraic curve. The formulas are expressed in terms of period matrices of holomorphic differentials on the curve; these holomorphic differentials must satisfy the condition τ * ν j = ν j , j = 1, . . . , g, where τ * is the action of τ lifted to the space of holomorphic differentials. This allows to give an algorithm to construct explicitly a symplectic transformation of an arbitrary canonical homology basis for a real Riemann surface to the symmetric form. The algorithm permits to systematically study smooth real solutions on general real Riemann surfaces for equations like KP and DS starting from a representation of the surface via an algebraic curve which so far was only possible for hyperelliptic surfaces. In addition, this provides a numerical way to compute the topological type of a real Riemann surface for given periods of the holomorphic differentials.

The paper is organized as follows. In section 2 we introduce a symmetric homology basis which depends on the topological type of R. In section 3 we give explicitly the symplectic transformation between the computed homology basis and the symmetric homology basis. This result will be used in section 4 to construct an algorithm which gives the topological type (g, k, a) of a real algebraic curve for given periods of the holomorphic differentials. In section 5 we discuss examples of real curves for higher genus. Some concluding remarks are added in section 6.

Symmetric homology basis

In what follows R denotes a compact Riemann surface of genus g > 0. A homology basis (A, B) := (A 1 , . . . , A g , B 1 , . . . , B g ) with the following intersection indices

A i • B j = δ ij A i • A j = B i • B j = 0,
is called a canonical basis of cycles. With A (resp. B) we denote the vector (A 1 , . . . , A g ) t (resp. (B 1 , . . . , B g ) t ). Canonical homology bases are related via a symplectic transformation. Let (A, B) and ( Ã, B) be arbitrary canonical homology bases on R. Then there exists a symplectic matrix

A B C D ∈ Sp(2g, Z) such that (2.1) A B C D Ã B = A B .
Recall that a symplectic matrix M ∈ Sp(2g, Z) satisfies M t J g M = J g , with the matrix J g given by J g = 0 I g -I g 0 , where I g denotes the g × g unit matrix. Symplectic matrices

M = A B C D ∈ Sp(2g, Z
) are characterized by the following system:

A t D-C t B = I g , (2.2) A t C = C t A, (2.3) D t B = B t D. (2.4)
Moreover, the inverse matrix M -1 is given by (2.5)

M -1 = D t -B t -C t A t .
Now let τ be an anti-holomorphic involution defined on R. Recall that (g, k, a) denotes the topological type of R, where k is the number of connected components of R(R) (the set of fixed points of τ ), and a = 0 if the curve is dividing (i.e., if R \ R(R) has two components), a = 1 if it is non-dividing (i.e., if R \ R(R) has just one component). A curve with the topological type (g, g + 1, 0) is called an M-curve. According to Proposition 2.2 in Vinnikov's paper [START_REF] Vinnikov | Self-adjoint determinantal representations of real plane curves[END_REF], there exists a canonical homology basis (A, B) (called for simplicity symmetric homology basis) such that

(2.6) τ A τ B = I g 0 H -I g A B ,
where H is a block diagonal g × g matrix which depends on the topological type (g, k, a) of R and is defined as follows:

-if k > 0 and a = 0,

H =              0 1 1 0 . . . 0 1 1 0 0 . . . 0              , -if k > 0 and a = 1, H =          1 . . . 1 0 . . . 0          ; rank(H) = g + 1 -k in both cases; -if k = 0, H =        0 1 1 0 . . . 0 1 1 0        for even g or H =          0 1 1 0 . . . 0 1 1 0 0          for odd g; rank(H) = g if g is even, rank(H) = g -1 if g is odd.
Remark 2.1. For given information whether there are real ovals (k > 0) or not (k = 0), the matrix H completely encodes the topological type of the real Riemann surface. In [START_REF] Seppäłä | Moduli Spaces for Real Algebraic Curves and Real Abelian Varieties[END_REF] a different, but equivalent form of H was used, which has ones only in the antidiagonal, if the rank is equal to the genus, or in a parallel to the antidiagonal for smaller rank.

Example 2.1. Consider the hyperelliptic curve of genus g defined by the equation

(2.7) y 2 = 2g+2 i=1 (x -x i ),
where the branch points x i ∈ R are ordered such that x 1 < . . . < x 2g+2 . On such a curve, we can define two anti-holomorphic involutions τ 1 and τ 2 , given respectively by τ 1 (x, y) = (x, y) and τ 2 (x, y) = (x, -y). Projections of real ovals of τ 1 on the x-plane coincide with the intervals [x 2g+2 , x 1 ], . . . , [x 2g , x 2g+1 ], whereas projections of real ovals of τ 2 on the x-plane coincide with the intervals

[x 1 , x 2 ], . . . , [x 2g+1 , x 2g+2 ]. Hence the curve (2.7
) is an M-curve with respect to both anti-involutions τ 1 and τ 2 . If all x i are non-real and pairwise conjugate, the curve has no real ovals with respect to the involution τ 2 ; it is dividing for the involution τ 1 .

Symplectic transformation between homology bases

In this section we construct a symplectic transformation between an arbitrary homology basis ( Ã, B) and a symmetric homology basis on a real Riemann surface (R, τ ), where τ denotes the anti-holomorphic involution defined on R. The main result of this paper is contained in Theorem 3.3 which gives the underlying symplectic matrix A B C D in terms of the period matrices of holomorphic differentials satisfying the condition (3.2) in the homology basis ( Ã, B). The key ingredient in this context is the description of the action of τ on the cycles ( Ã, B), given in Proposition 3.1 by the integer matrix R. Then Theorem 3.3 states that the column vectors of the matrix A B t form in fact a Z-basis of the integer kernel of the matrix R t -I 2g . The matrix Q in Theorem 3.3 encodes the degree of freedom in the choice of such a Z-basis. For the ease of the reader, we start recalling some basic facts from the theory of Z-modules used to prove Theorem 3.3, which differs from the usual linear algebra over vector spaces.

3.1. Basic facts from the theory of Z-modules. A Z-module or more generally a Amodule where A denotes a commutative ring, is a natural generalization of vector spaces where the usual scalar field is replaced by the ring A. For a review on the subject we refer to [START_REF] Lang | Algebra, Revised Third Edition[END_REF]. In what follows we assume that A is the principal ring Z and we denote by M a Z-module. Definition 3.1. i. M is of finite type if it admits a finite set of generators. ii. M is said to be free if there exists a Z-basis, namely, a set {x i } i∈I ⊂ M with I ⊂ N such that any element x ∈ M can be written uniquely as x = i∈I α i x i where the scalars α i ∈ Z are non-zero only for a finite number of them.

The following theorems provide important results in the general theory of modules over a principal ring. Theorem 3.1.

i. If M is free and of finite type then all Z-bases of M are finite with the same cardinality called the rank of M. ii. Since Z is integral, the rank of M equals the dimension of the Q-vector space S -1 M where S = Z \ {0}. iii. A submodule of a free Z-module of finite type of rank n is a free Z-module of finite type of rank r ≤ n.

The example we will need in this paper are submodules of Z n which are because of the Theorem 3.1 free modules of finite type with rank r ≤ n.

Remark 3.1. Contrary to the case of vector spaces, if N is a free submodule of the free module M with the same rank, this does not imply that N equals M. Consider the submodule 2Z of Z as an example.

The following theorem, also called the theorem of the adapted basis, gives the classification of modules over the principal ring Z: Theorem 3.2. Let M be a free Z-module of rank n and let N be a submodule of M. Then there exists a basis (e 1 , . . . , e n ) of M and unique non-zero integers (p 1 , . . . , p r ) (with

r ≤ n) such that (1) (p 1 e 1 , . . . , p r e r ) is a Z-basis of N (2) p 1 |p 2 | . . . |p r (this means that ∀i, p i divides all p j with j > i).
In what follows GL n (Z) denotes the set of n × n invertible matrices over Z; it is well known that the determinant of these matrices equals ±1. Notice that if M denotes the matrix formed by the column vectors of a Z-basis of a Z-module of rank n, then other bases are given by the column vectors of matrices of the form M Q with Q ∈ GL n (Z). Therefore, Theorem 3.2 has the following matrix interpretation that we will use in section 3.3: for any non-zero m × n matrix M ∈ M m,n (Z) of rank r, there exist matrices

U ∈ GL m (Z), V ∈ GL n (Z) such that (3.1) U M V = Diag(p 1 , . . . , p r ) 0 0 0 ,
where p i ∈ Z\{0} satisfy p 1 |p 2 | . . . |p r and where Diag(.) denotes the diagonal matrix. This is called the Smith normal form of M . In section 4 we will use a well known algorithm to compute the Smith normal form. In particular, this algorithm provides a Z-basis of the integer kernel of M given by the last nr column vectors of the matrix V in (3.1).

3.2.

Action of τ on an arbitrary homology basis. We denote by (A, B) a symmetric homology basis (i.e. which satisfies (2.6)). Let (ν 1 , . . . , ν g ) be a basis of holomorphic differentials such that

(3.2) τ * ν j = ν j , j = 1, . . . , g,
where τ * is the action of τ lifted to the space of holomorphic differentials: τ * ω(p) = ω(τ p) for any p ∈ R. The matrices P A and P B defined by

(3.3) (P A ) ij = A i ν j , (P B ) ij = B i ν j , i, j = 1, . . . , g
are called the matrices of Aand B-periods of the differentials ν j . From (2.6) and (3.2) we deduce the action of the complex conjugation on the matrices P A and P B :

(3.4)

(P A ) ij ∈ R, (3.5
)

P B = -P B + HP A .
Denote by ( Ã, B) an arbitrary homology basis. From the symplectic transformation (2.1) we obtain the following transformation law between the matrices P Ã, P B and P A , P B defined in (3.3):

(3.6) A B C D P Ã P B = P A P B .
Therefore, by (3.4) one gets

A Re P Ã + B Re P B = P A (3.7) A Im P Ã + B Im P B = 0, (3.8)
and by (3.5)

C Re P Ã + D Re P B = 1 2 HP A (3.9) C Im P Ã + D Im P B = Im (P B ) . (3.10)
From (3.7) and (3.10) it can be checked that the matrices A Re P Ã + B Re P B and C Im P Ã + D Im P B are invertible (the first because P A is, for the second see [START_REF] Kalla | On the numerical evaluation of algebro-geometric solutions to integrable equations[END_REF] for more details). The following Lemma proved in [START_REF] Kalla | On the numerical evaluation of algebro-geometric solutions to integrable equations[END_REF] shows that it is sufficient to know the pairs of matrices A, B or C, D to get the full symplectic transformation (3.6): (3.12)

C t = 1 2 A t H -Re P B A Re P Ã + B Re P B -1
(3.13)

D t = 1 2 B t H + Re P Ã A Re P Ã + B Re P B -1 . (3.14)
The action of τ on an arbitrary homology basis ( Ã, B) can be written as

(3.15) τ Ã τ B = R Ã B ,
where R ∈ M 2g (Z). In the following proposition, we give an explicit expression for the matrix R in terms of the period matrices P Ã and P B only.

Proposition 3.1. The matrix R defined in (3.15) is given by

(3.16) R =   2 Re P B M-1 Im P t à + I g t -2 Re P à M-1 Im P t à 2 Re P B M-1 Im P t B -2 Re P B M-1 Im P t à + I g   ,
where

(3.17) M = Im P t B Re P à -Im P t à Re P B .
Proof. Using (2.6) we deduce the action of τ on (2.1):

(3.18) A B C D -1 I g 0 H -I g A B = R A B C D -1 A B , which yields (3.19) R = A B C D -1 I g 0 H -I g A B C D .
In other words, a symplectic matrix which transforms a basis ( Ã, B) to a symmetric form also satisfies (3.19). From (3.19) and (2.5) one gets

(3.20) R = (2 C t B -A t HB + I g ) t 2 D t B -B t HB -2 C t A + A t HA -(2 C t B -A t HB + I g ) .
Replacing

2 C t A = A t HA -2 Re P B M-1 Im P t B , 2 C t B = A t HB + 2 Re P B M-1 Im P t à , 2 D t B = B t HB -2 Re P à M-1 Im P t à ,
with M given by (3.17). Substituting these relations in (3.20) one gets (3.16).

Symplectic transformation.

In this part we present in Theorem 3.3 the main result of the present paper: the symplectic transformation between an arbitrary basis ( Ã, B) on a real Riemann surface and a homology basis adapted to the symmetry is given in terms of the period matrices P Ã, P B defined in the previous section. This result will allow to construct in section 4 an algorithm which computes the topological type of a given real Riemann surface. We start with the following lemma which describes the spectral properties of the matrix R (3.16): Lemma 3.2. The matrix R in (3.16) is diagonalizable over Q with eigenvalues 1 and -1. The dimension of the corresponding eigenspaces equals g.

Proof.

It is straightforward to see that the matrix

I g 0 H -I g is diagonalizable over Q:
the eigenvalues are 1 and -1, and the dimension of the corresponding eigenspaces equals g. Therefore, by (3.19) the same holds for the matrix R.

In what follows we denote by

(3.21) K Z := w ∈ Z 2g ; (R t -I 2g )w = 0
the integer kernel of the matrix R t -I 2g . According to the theory of modules over the principal ring Z, the Z-module K Z admits a Z-basis which can, for instance, be computed from the Smith normal form (see section 3.1).

Theorem 3.3. Let R be a real compact Riemann surface of genus g and ( Ã, B) a canonical homology basis on R. For the given matrices of periods P Ã, P B of holomorphic differentials satisfying (3.2), let S 1 S 2 be the 2g × g matrix formed by a Z-basis of K Z (3.21). Then a symplectic matrix in (2.1) which relates the homology basis ( Ã, B) to a symmetric homology basis (A, B) is given by:

A t B t = S 1 S 2 Q (3.22) C t D t = 1 2 S 1 S 2 QH + -Re P B Re P Ã S t 1 Re P Ã + S t 2 Re P B -1 (Q t ) -1 (3.23)
where the matrix H is defined in section 2, and where Q ∈ GL g (Z) is such that

(3.24) S 1 S 2 QHQ t ≡ 2 -Re P B Re P Ã S t 1 Re P Ã + S t 2 Re P B -1
(mod 2).

Proof. Notice that (3.19) can be rewritten as

(3.25) A B C D R = I g 0 H -I g A B C D ,
which, in particular, gives the following condition for the matrices A and B:

(3.26) A B (R -I 2g ) = 0.
Denote by u 1 , . . . , u g , v 1 , . . . , v g the column vectors of the matrix A B C D t

. By (3.26) the vectors u i for i = 1, . . . , g lie in the integer kernel K Z of the matrix R t -I 2g . Let us prove that (u 1 , . . . , u g ) form in fact a Z-basis of the module K Z . By Lemma 3.2 one has dim(K Q ) = g, where K Q denotes the kernel of R t -I 2g over the field Q, which by Theorem 3.1 yields rank(K Z ) = g. Therefore one has to check that (u 1 , . . . , u g ) are free vectors over Z and generate the Z-module K Z . Notice that here it is important to check that these vectors form a set of generators for the Z-module K Z , as we saw in Remark 3.1.

Since

A B C D

t ∈ GL 2g (Z), the 2g vectors u 1 , . . . , u g , v 1 , . . . , v g form a Z-basis of the module Z 2g , which in particular implies that these vectors are free over Z. Then it remains to prove that the vectors u i , i = 1, . . . , g generate the Z-module K Z which is done by contradiction as follows. Let w ∈ K Z which we write in Z 2g as w = g i=1 α i u i + g j=1 β j v j with α i , β j ∈ Z such that at least one of the β j is non-zero. Since w, u 1 , . . . , u g ∈ K Z , one has v := g j=1 β j v j ∈ K Z and v = 0. We deduce that u 1 , . . . , u g , v are g + 1 free vectors in K Z . This is impossible since rank(K Z ) = g. Thus one has β j = 0 for j = 1, . . . , g which implies that the vectors u i , i = 1, . . . , g generate the Z-module K Z .

Hence we can write (3.27)

A t B t = S 1 S 2 Q, for some Q ∈ GL g (Z)
, where the g column vectors of the matrix S 1 S 2 form a Z-basis of the module K Z . Here the matrix Q encodes the freedom in the choice of such a basis. The matrices C and D are then given by (3.13) and (3.14). It follows that these two matrices are integer matrices if and only if the matrix Q satisfies (3.24), which completes the proof.

Remark 3.2. If the curve is an M-curve, then H = 0 and the matrices A, B, C, D in Theorem 3.3 are given by

A t B t = S 1 S 2 Q (3.28) C t D t = -Re P B Re P Ã S t 1 Re P Ã + S t 2 Re P B -1 (Q t ) -1 (3.29)
where Q ∈ GL g (Z) is arbitrary.

Remark 3.3. As explained below, if the curve is not an M-curve, namely, H = 0, one can construct explicitly a matrix Q ∈ GL g (Z/2Z) such that (3.24) holds. This construction of Q is based on the Smith normal form of the matrix S 1 S 2 (see section 3.1) which allows the simplification of the system (3.24).

Lemma 3.3. There exist

U ∈ GL 2g (Z), V ∈ GL g (Z) such that (3.30) U S 1 S 2 V = E 0 ,
where E is a g × g diagonal matrix with elements E ii = ±1 for i = 1, . . . , g.

Proof. By (3.1) there exist U ∈ GL 2g (Z), V ∈ GL g (Z) and p 1 , . . . , p g ∈ N \ {0} satisfying

p 1 |p 2 | . . . |p g such that (3.31) U S 1 S 2 V = D 0 ,
where D := Diag(p 1 , . . . , p g ). The fact that D = E can be deduced from the following equalities:

1 = det A B C D = det S 1 C t S 2 D t det(Q) = ± det D F 1 0 F 2 = ± det(D) det(F 2 ) (3.32)
where we multiplied the matrix in the determinant from the left by U and from the right by V 0 0 I g , and used

det(Q) det(U ) det(V ) = ±1 since Q, V ∈ GL g (Z) and U ∈ GL 2g (Z); here F 1 F 2 = U C t D t .
We deduce that det(D) = ±1 since the right-hand side of (3.32) is a product of determinants of matrices with integer coefficients. This completes the proof. Now let us define matrices N 1 , N 2 ∈ M g (Z) as follows:

(3.33) N 1 N 2 := 2 V U -Re P B Re P Ã S t 1 Re P Ã + S t 2 Re P B -1 .
Then one has:

Proposition 3.2. Q ∈ GL g (Z)
(mod 2),
which is equivalent to

I g 0 QHQ t ≡ 2 V U -Re P B Re P Ã S t 1 Re P Ã + S t 2 Re P B -1 (mod 2).
Using the definition (3.33) one gets (3.34). Now to check that a transformation of the form Q → QQ 0 where Q 0 solves (3.35) is the only one which preserves (3.34), notice that such a transformation corresponds to a symplectic transformation between the symmetric homology basis obtained from Theorem 3.3 and another symmetric homology basis, since this coincides with a change of the Zbasis in (3.22). Hence from (2.6) it is straightforward to see that the symplectic matrix which relates two symmetric homology bases is given by

(3.36) Q t 0 0 1 2 (HQ t 0 -Q -1 0 H) Q -1 0 , with Q 0 ∈ GL g (Z) satisfying Q 0 HQ t 0 ≡ H (mod 2)
. This completes the proof.

4.

Algorithm for the computation of the topological type (g, k, a)

The results of the previous section allow us to formulate an algorithm to transform an arbitrary canonical homology basis ( Ã, B), for instance obtained via the algorithm [START_REF] Tretkoff | Combinatorial group theory, Riemann surfaces and differential equations[END_REF], to a symmetric basis (A, B) satisfying (2.6). The key task in this context is the computation of the matrix Q in (3.34). In the process of computing Q, the matrix H giving the topology of the real Riemann surface can be determined.

The starting point of the algorithm are the periods P Ã and P B of a basis of differentials (ν 1 , . . . , ν g ) satisfying (3.2). Notice that condition (3.2) is important. The Maple algcurves package generates such differentials for real curves by default. For the Matlab code used in this paper this is in general not the case if a numerically optimal approach is used to determine the holomorphic differentials, see [START_REF] Frauendiener | Algebraic curves and Riemann surfaces in Matlab[END_REF] for details. However, it is possible to determine a basis of the holomorphic differentials with rational coefficients which will satisfy condition (3.2). For the examples in the following section, we always choose this option.

With these periods the code computes the matrix R via (3.16). This matrix will have integer entries up to the used precision (by default 10 -6 in Maple and 10 -12 in Matlab). Rounding has to be used to obtain an integer matrix. Computing the Smith normal form (3.1) of the matrix R t -I 2g , we obtain from the last g vectors of the resulting matrix V a Z-basis of the integer kernel K Z (3.21), i.e., the column vectors of the matrix S 1 S 2 in (3.27). An algorithm to compute the Smith normal form U M V = S of an integer matrix M is implemented in Maple. This algorithm can be called from Matlab via the symbolic toolbox. We use here an own implementation of the standard algorithm to compute the Smith normal form which we briefly summarize: we always work on column j starting with j = 1. If there is no nonzero element in this column, it is swapped by multiplication with an appropriate matrix U from the left with the last column with a nonzero element. If the element M jj = 0, a row with a nonzero element in position j is added (to avoid clumsy notation, the transformed matrix is still called M ). Then all nonzero elements M jk for k = j are eliminated by adding row j with appropriate multipliers obtained via the Euclidean algorithm. In the same way the row with index j is cleaned by acting on M via multiplication by a matrix V from the right. If the resulting element M jj does not divide all other elements of M , one of these elements is added by using the Euclidean algorithm to M jj in a way that the latter becomes smaller. This destroys possibly the nullity of the remaining elements in column j and row j which thus have to be cleaned as before. This process is repeated until M jj divides all other elements of M . Then the index j is incremented by 1. The procedure is repeated until the Smith normal form is obtained.

Notice that this standard algorithm has a well known problem: in general for larger matrices the entries of the matrices U and V become very large. Though these are integer matrices, this is problematic once some of the entries are of the order of 10 16 (machine precision in Matlab is 10 -16 which implies that integers of the order of 10 16 , which are internally treated as floating point numbers, can no longer be numerically distinguished). There are more sophisticated algorithms to treat larger matrices as the one given in [START_REF] Hafner | Asymptotically fast triangularization of matrices over rings[END_REF]. In practice the standard algorithm works well for examples of a genus g ≤ 6 which is sufficient for our purposes. Only for higher genus, the algorithm [START_REF] Hafner | Asymptotically fast triangularization of matrices over rings[END_REF] would be needed.

By computing the Smith normal form of the vector S 1 S 2 in (3.30), we get the needed quantities to compute the matrix N 1 in (3.33). The main task is then to determine the matrix Q in (3.24) since for given periods, the whole symplectic matrix (2.1) follows from equations (3.22) and (3.23) for given S 1 S 2 . The matrices Q and H can be determined from relation (3.34) for a given matrix N 1 by standard Gaussian elimination in Z/2Z and by imposing the block diagonal form of section 2 on H as we will outline below.

Remark 4.1. It was shown in section 2 that the matrix H can be chosen to be either diagonal or to consist of blocks of the form

H 0 = 0 1 1 0 .
A unique determination of the matrix H in the computation of Q from (3.34) is only possible, if this block H 0 cannot be related through a similarity transformation in Z/2Z to I 2 . In fact this is the case since H 0 cannot be diagonalized in Z/2Z. The same reasoning applies if the matrix H consists of several blocks H 0 and zeros otherwise. However, a block of the form

H 1 =   1 0 0 0 0 1 0 1 0  
can be diagonalized in Z/2Z by multiplication from the left and the right by a matrix of the form

Q 1 =   1 1 1 1 0 1 1 1 0   .
It follows that if the matrix H determined in the computation of Q from (3.34) has a non-zero diagonal element, then H can be diagonalized in several steps: if the matrix H has the form H below ( Hij = δ ij for i, j < k and a block H 0 for i, j = k, k + 1)

H =             1 . . . 1 0 1 1 0             Q =             1 1 1 1 0 1 1 1 0            
then multiplication from the left and the right with a matrix Q (not shown elements of this matrix are 0) gives the identity matrix for i, j ≤ k + 1. Applying this procedure several times will lead to a diagonal matrix H.

The algorithm for the computation of Q and H via the similarity relation (3.34) for given N 1 by imposing a block diagonal form for the matrix H as in section 2 uses in principle standard Gauss elimination on the rows and columns of N 1 over the field Z/2Z with minor modifications as detailed below. We only describe the action on the columns via a matrix q from the right, since the action on the rows follows by symmetry by multiplication with q t from the left: -if N 1 ≡ 0 (mod 2), put H = 0 and Q = I g and end the algorithm, otherwise, put the index j of the column under consideration equal to 1; -if column j contains only zeros, it is swapped with the last column with non-zero entries; -if there is a 1 in position j of the column, all further non-zero entries in the column are eliminated in standard way; -if there is a 1 in the column, but not in position j, rows are swapped in a way that it appears in the position j + 1 of the column (it cannot be put to position j as explained in Remark 4.1). Further ones in the column are eliminated; -if there was a non-zero diagonal entry, the column index j is incremented by 1, if there was a block 0 1 1 0 , the index j is incremented by 2. Then the algorithm is repeated with column j until j = g or until the columns with index j and higher only contain zeros; -if there are blocks of the form H in Remark 4.1, then H will be diagonalized by multiplication with the corresponding matrix Q as explained in 4.1.

Examples

In this section we study examples of real algebraic curves, provide the computed periods and the application of the algorithm to obtain a symmetric homology basis as well as the matrix H encoding the topological information of the curve. For convenience we use here the Matlab algebraic curves package, but the same examples can be of course studied with the Maple package. We also give graphical representations of the real variety of an algebraic curve if there is any which is generated via contour plots of f (x, y) = 0 for real x and y (this corresponds to the command plot_real_curve in the Maple algcurves package). This is not identical with the set R(R) of real ovals of the Riemann surface since the curves may have singularities, whereas the Riemann surface is defined by desingularized such curves (see [START_REF] Frauendiener | Algebraic curves and Riemann surfaces in Matlab[END_REF] for how this is done in the Matlab package). Thus there may be cusps and selfintersections in the shown plots. Moreover these plots are not conclusive if the curves come very close, and if there are self-intersections as in Fig. 4 below. In addition we only show the curves for finite values of x and y from which it cannot be decided which curves cross at infinity and which lines belong to the same ovals as in Fig. 5. They just serve for illustration purpose, for more sophisticated approaches, see [START_REF] Arnon | A polynomial time algorithm for the topological type of a real algebraic curve[END_REF][START_REF] Coste | Thom's lemma, the coding of real algebraic numbers and the computation of the topology of semi-algebraic sets[END_REF][START_REF] Sakkalis | The topological configuration of a real algebraic curve[END_REF][START_REF] Feng | Decomposition and Computation of the Topology of Plane Real Algebraic Curves[END_REF][START_REF] Gonzalez-Vega | Efficient topology determination of implicitly defined algebraic plane curves[END_REF][START_REF] Hong | An efficient method for analyzing the topology of plane real algebraic curves[END_REF][START_REF] Seidel | On the exact computation of the topology of real algebraic curves[END_REF]. The computed number of real ovals via the algorithm is, however, unique: as outlined in section 2, it follows from the rank of the matrix H (for k = 0 one has k = g + 1rank(H)). We always assume in the following that it is known whether there are any real ovals. This allows the unique identification of the topological type (g, k, a) via the matrix H.

The Trott curve [START_REF] Trott | Applying Groebner Basis to Three Problems in Geometry[END_REF] given by the algebraic equation (5.1) 144 (x 4 + y 4 ) -225 (x 2 + y 2 ) + 350 x 2 y 2 + 81 = 0 is known to be an M-curve of genus 3 (it has the maximal number g + 1 = 4 of real ovals, as can be seen in Fig. 1). Moreover, this curve has real branch points only (and 28 real bitangents, namely, tangent lines to the curve in two places). Our computed matrices of à and B-periods denoted by aper and bper respectively read 1 aper = -0.0000 + 0.0235i -0.0000 + 0.0138i -0.0000 + 0.0138i 0 + 0.0000i 0.0000 + 0.0277i 0 + 0.0000i -0.0315 -0.0000 + 0.0000i 0.0250 -0.0000i bper =

1 For the ease of representation we only give 4 digits here, though at least 12 digits are known for these quantities.

-0.0315 + 0.0235i -0.0000 + 0.0138i -0.0250 + 0.0138i -0.0000 + 0.0000i -0.0250 + 0.0277i 0.0250 -0.0000i -0.0000 -0.0235i 0.0000 + 0.0138i 0 + 0.0138i. For this the algorithm produces as expected H = 0 and Q of (3.34) the identity matrix. The symplectic transformation found by the algorithm via (3.22) and (3.23) has the form

[A,B,C,D] = 0 1 0 0 -1 0 0 1 0 0 0 0 1 0 0 -1 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0
1. Since we will not actually use the matrices A, B, C, D in this article, and since they follow for given periods from H and Q via (3.22) and (3.23) we will only give them for this example. However they are needed e.g. for the study of algebro-geometric solutions to integrable equations as KP and DS as in [START_REF] Kalla | On the numerical evaluation of algebro-geometric solutions to integrable equations[END_REF].

The Klein curve given by the equation (5.2)

y 7 -x(x -1) 2 = 0
has the maximal number of automorphisms (168) of a genus 3 curve. The computed periods read aper = -0.9667 + 0.7709i 0.9667 + 0.2206i 0.9667 -2.0073i -1.2054 -0.2751i -0.4302 + 0.8933i -1.7419 + 1.3891i -0.4302 -0.8933i 1.7419 + 1.3891i -1.2054 + 0.2751i bper = -2.7085 -0.6182i -0.2387 + 0.4958i 1.3969 -1.1140i -2.1721 -1.7322i 0.5365 -0.1224i -0.7752 -1.6097i 0.9667 + 0.2206i -0.9667 + 2.0073i -0.9667 + 0.7709i.

The algorithm finds

H = 1 0 0 0 1 0 0 0 1 Q = 1 1 1 0 0 1 0 1 0.
Therefore, the topological type of the curve is (3, 1, 1), namely, the curve has genus 3, one real oval (as can be also seen in Fig. 2) and is non-dividing.

The Fermat curve (5.3) has for n = 4 the topological type (3, 0, 1). This is confirmed by the algorithm. For the periods aper = 0.9270 + 0.0000i 0.0000 -0.9270i 0.0000 -0.9270i -0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 -1.8541i 0 + 0.9270i -0.9270 + 0.0000i 0.0000 -0.9270i bper = 0.9270 + 0.9270i 0.9270 -0.9270i 0.0000 + 0.0000i 0.0000 -0.9270 + 0.9270i 0.9270 -0.9270i -0.9270 + 0.0000i 0.0000 -0.9270i 0.0000 -0.9270i

y n + x n + 1 = 0 x y -5 -4 -3 -2 -1 0 1 2 3
we find

H = 0 1 0 1 0 0 0 0 0 Q = 1 0 0 0 0 1 0
1 0 in accordance with the expectation. For n = 5 the curve has genus 6. We find with aper = Columns 1 through 4 -0.1623 + 0.4995i -0.1890 + 0.5817i -0.1623 + 0.4995i 0.4948 + 0.3595i -0.3246 -0.0000i -0.3780 + 0.0000i -0.3246 + 0.0000i 0.9896 -0.0000i -0.4249 + 0.3087i -0.0722 -0.2222i 0.5252 -0.0000i 0.6116 -0.0000i -0.4249 -0.3087i

= 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 Q = 1 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 1 0 0 0 1 0 1 0 0 1 1 0 1 0 0 0 1 0 1
1. This implies that there is one real oval as can be also seen in Fig. 3, and the curve is non-dividing. This corresponds to (g, k, a) = (6, 1, 1). This gives the topological type (3, 1, 1). In particular, the number of real ovals equals one. The real variety of the curve, which has a self intersection and a cusp, can be seen in Fig. 4. The curve (5.5) -180x 5 +396yx 4 -307x 3 y 2 +107x 2 y 3 +273x 3 -318x 2 y-17xy 4 +117xy 2 -68x+y 5 -12y 3 +19y = 0 is known (see [START_REF] Dubrovin | Matrix finite-zone operators[END_REF]) to be a dividing curve of genus 6. In fact we get for aper = Columns 1 through 4 0.0414 + 0.0278i -0.0345 + 0.0272i -0.0979 + 0.0264i -0.3041 + 0.0603i -0.1149 -0.0446i 0.0000 -0.0000i 0.0532 -0.0226i 0.0000 -0.0000i 0.1149 + 0.0169i -0.0000 + 0.0544i -0.0532 -0.0805i -0.0000 + 0.1206i 0.0000 -0.0111i 0.0000 + 0.0183i 0.0000 -0.0303i 0.0000 + 0.1114i 0.0820 -0.0278i 0.0000 + 0.0000i -0.0527 -0.1031i 0.0000 + 0.0000i -0.0820 -0.0000i -0.0000 -0.0544i 0.0527 -0.0000i -0.0000 -0.1206i

Columns 5 through 6 -0.5369 + 0.0872i -2.8149 + 0.5433i 0.2427 + 0.0099i 0.7641 + 0.0357i -0.2427 -0.2372i -0.7641 -0.5381i 0.0000 -0.1843i 0.0000 -1.1224i -0.0881 -0.2274i -0.1275 -0.5024i 0.0881 -0.0000i 0.1275 -0.0000i bper = Columns 1 through 4 0.0414 + 0.0278i -0.0345 -0.0094i -0.0979 + 0.0264i -0.3041 -0.1626i -0.0089 -0.1009i -0.0666 + 0.0180i 0.1091 -0.0733i -0.1162 + 0.0046i -0.0089 -0.0563i -0.0666 -0.0180i 0.1091 -0.0507i -0.1162 -0.0046i 0.0320 -0.0000i 0.0000 -0.0183i 0.1425 -0.0000i 0.0000 -0.1114i 0.1060 + 0.0286i 0.0666 + 0.0724i 0.0559 -0.0525i 0.1162 + 0.1252i -0.0580 + 0.0160i 0.0666 -0.0363i 0.1614 + 0.0751i 0.1162 -0.1160i

Columns 5 through 6 -0.5369 + 0.0872i -2.8149 + 0.5433i 0.3380 -0.1055i 0.9555 -0.2619i 0.3380 -0.1154i 0.9555 -0.2976i 0.8311 -0. This means that it is a dividing curve with 3 real ovals, which is, however, not obvious from Fig. 5. The topological type is (6, 3, 0).

Outlook

The algorithm presented in this paper allows the transformation of an arbitrary canonical homology basis to a form adapted to the underlying symmetry, here the anti-holomorphic involution. This permits to find a basis satisfying relation (2.6). A similar condition can be imposed for any involution, and the algorithm presented here can be easily adapted to that case.

In general, the presence of symmetries allows to significantly simplify the Riemann matrix of a surface, but only in a homology basis adapted to the symmetries, for instance in a basis such that the A-cycles are invariant under symmetry operations, see [START_REF] Belokolos | Algebro-geometric approach to nonlinear integrable equations[END_REF] and [5] for the Klein curve. The latter reference uses an approach to this problem based on Comessati's theorem [START_REF] Comessati | Sulla connessione delle superficie algebriche reali[END_REF] via two pieces of software, extcurves and CyclePainter2 . It will be the subject of further work to generalize the algorithm to symmetry groups beyond involutions. In a first step it would be interesting to extend the approach presented in this paper to automorphisms τ satisfying τ n = id with n > 2.

Lemma 3 . 1 . 1 ( 3

 3113 The matrices A, B, C, D ∈ M g (Z) solving (3.7)-(3.10) satisfy: A t = Im P B C Im P Ã + D Im P B -.11) B t = -Im P Ã C Im P Ã + D Im P B -1

Figure 1 .

 1 Figure 1. Real ovals of the Trott curve (5.1), an M-curve of genus 3.

Figure 2 .

 2 Figure 2. Real variety of the Klein curve (5.2), the curve with the maximal number of automorphisms in genus 3.

Figure 3 .

 3 Figure 3. Real variety of the Fermat curve (5.3) for n = 5. For the curve (5.4) y 3 -2x 3 yx 9 = 0 of genus 3 we have aper =

Figure 4 .

 4 Figure 4. Real variety of the curve (5.4).

Figure 5 .

 5 Figure 5. Real variety of the curve (5.5).

  A via(3.11) and B via (3.12) in (3.13) and then using(3.11) again to eliminate the factor C Im P Ã + D Im P B leads to a relation for C t A only (using in the last step (3.12) instead of (3.11) gives a relation for C t B). Similarly one gets a relation for D t B,

  where N 1 ∈ M g (Z) is defined in(3.33). Moreover, if Q is a particular solution of (3.34) then the general solution can be written as QQ 0 where Q 0 solves

	(3.35)		Q 0 HQ t 0 ≡ H (mod 2).
	Proof. Multiplying equality (3.24) from the left by the matrix U of Lemma 3.3 and using
	(3.30) one gets
	I g 0	V -1 QHQ t ≡ 2 U	-Re P B Re P Ã	S t 1 Re P Ã + S t 2 Re P B	-1

satisfies (3.24) if and only if it solves (3.34) QHQ t ≡ N 1 (mod 2),

Located at http://gitorious.org/riemanncycles.


Hua Sun

Ru Li

Eric Chénier

Guy Lauriat
email: lauriat@univ-mlv.fr

On the modeling of aiding mixed convection in vertical channels

Keywords: Aiding mixed convection, Vertical open channel flows, Computational fluid dynamics, Boundary conditions Nomenclature

This paper aims at showing that to prescribe a flow rate at the inlet section of a vertical channel with heated walls leads to surprising and counterintuitive physical solutions, especially when the problem is modeled as elliptical. Such an approach can give rise to the onset of recirculation cells in the entry region while the heat transfer is slightly increased under the influence of buoyancy force. We suggest an alternative model based on more realistic boundary conditions based on a prescribed total pressure at the inlet and a fixed pressure at the outlet sections. In this case, the pressure and buoyancy forces act effectively in the same direction and, the concept of buoyancy aiding convection makes sense. The numerical solutions emphasize the large differences between prescribed inlet velocity and solutions provided with our pressure based boundary conditions.

a thermal diffusivity [m 2 s -1 ] A aspect ratio, A = H/D c p specific heat [J K -1 kg -1 ] c x ,
Mixed convection may be defined as heat transfer situations where both pressure and buoyancy forces interact. In vertical channels, the bulk fluid flow can be either upward or downward, and the solutal and thermal buoyancy forces may be either assisting or opposing the forced flow. As a result, several hundred of papers were published just for mixed convection flows in vertical ducts of various shapes and thermal conditions, including two phase flows with surface condensation or evaporation, and flows in porous filled ducts. Most studies considered buoyancy assisted flow, i.e. upward flows with heating or downward flows with cooling. Under these conditions, the axial velocities may increase near the channel walls and decrease within the core region, with possible occurrence of a flow recirculation in the entry region for large effects of the buoyancy force [START_REF] Ingham | Recirculating pipe flows[END_REF][START_REF] Hammou | Laminar mixed convection of humid air in a vertical channel with evaporation or condensation[END_REF][START_REF] Desrayaud | Flow reversal of laminar mixed convection in the entry region of symmetrically heated, vertical plate channels[END_REF][START_REF] Laaroussi | Effects of variable density for film evaporation on laminar mixed convection in a vertical channel[END_REF] (not to be confused with flow reversal occurring close to the outlet section in asymmetrically heated channels [START_REF] Aung | Theory of fully developed, combined convection including flow reversal[END_REF][START_REF] Aung | Developing flow and flow reversal in a vertical channel with asymmetric wall temperatures[END_REF][START_REF] Gau | Measurements of heat transfer and flow structure in heated vertical channels with buoyancy assisted and opposed flows[END_REF][START_REF] Gau | Reversed flow structure and heat transfer measurements for buoyancy-assisted convection in a heated vertical duct[END_REF][START_REF] Cheng | Enhancement of flare method for predicting buoyancy-induced flow reversal in vertical ducts via parabolic model[END_REF][START_REF] Yang | Numerical and analytical study of reversed flow and heat transfer in a heated vertical duct[END_REF]). In these cases, it is generally expected that heat transfer for mixed convection is larger than for forced convection. This intuition is obviously well physically based. However, numerical solutions may lead to opposite effects, not owing to the numerical accuracy but because the boundary conditions are not well specified. This paper is focused on this issue.

The extensive interest borne upon mixed convection in ducts stems from its wide range of practical applications including the design of compact heat exchangers, cooling of electronic equipment and, solar energy collectors. Numerous theoretical and experimental investigations were reviewed in many textbooks, such as in Aung [START_REF] Aung | Mixed convection in internal flow, Handbook of Single-Phase Convective Heat Transfer[END_REF], Gebhart et al. [START_REF] Gebhart | Buoyancy Induced Flows and Transport[END_REF] or Bejan [START_REF] Bejan | Convection Heat Transfer[END_REF]. The present study is focused on mixed convection of air as the working fluid for typical geometries used for active cooling from the back of vertical photovoltaic cell panels operating at high fluxes.

This work considers the effect of thermal buoyancy force only (the effect of solutal force being similar, provided it acts in the same direction as the thermal force), and in the case of upward forced or natural flow direction (symmetrically heated vertical channel), commonly termed as buoyancy aiding or assisting mixed convective flows.

Most of the analytical and numerical studies on mixed convection were based on approximate problem formulations in order to easily solve the governing equations thanks to parabolized forms (or boundary-layer type equations for external flows with the streamwise pressure gradient dropped out [START_REF] Chen | Natural Convection: Fundamentals and Applications[END_REF]) which allow not to prescribe boundary conditions at the outlet duct section. On the other hand, the increasingly use of commercial computational fluid dynamics packages during the last decade has led to solve the problem when it is based on fully elliptic formulation, mostly by considering again prescribed inlet flow rate and outlet free boundary conditions. However, despite the large number of papers available in the current literature, the fundamental question of inlet/outlet flow conditions seems not to have been properly addressed because it appears that mixed convection with fixed inlet velocity (and thus Reynolds number) was very often considered as the relevant formulation for solving mixed convection.

The aim of this paper is at showing that the use of prescribed inlet velocity condition (constant flow rate whatever the importance of the buoyancy force) leads to completely different solutions in comparison with those obtained using prescribed total pressure inlet and pressure outlet conditions. For our viewpoint, such pressure boundary conditions appear to provide more realistic physical behaviors in practical applications, as it will be shown in what follows.

We are considering the practical case of air flowing in vertical flat-plate channels submitted to uniform wall temperatures, and maximum temperature differences compatible with the Boussinesq approximation. Although most of the results are presented in dimensional form, we also analyze the relevance of scaling predictions based on dimensionless parameters. It should be added that the extension of the present work to weakly compressible formulation (i.e. large temperature differences) is straightforward when using our in-house computational code, but the results obtained have shown that the main conclusions drawn in this paper remain unchanged. We consider two-dimensional, incompressible and laminar buoyancy-assisting flows inside a vertical parallel-plate channel, as shown in Fig. 1. The fluid enters the channel of height H at ambient temperature and traverses upward, being heated by the hot walls at uniform temperature T h . The fully developed region may be eventually reached at the outlet (z = H), after a development length depending on the value of the plate spacing as well as on importances of the pressure and buoyancy forces. On account of the maximum temperature difference invoked, it is assumed that all physical properties are constant except for the density in the buoyancy force in the vertical direction (Boussinesq approximation). The reference temperature has been taken to be the inlet temperature, T 0 , like in most of the solutions reported in the current literature. However, this choice remains an open question, as discussed by Barletta and Zanchini [START_REF] Barletta | On the choice of the reference temperature for fullydeveloped mixed convection in a vertical channel[END_REF], especially for fully developed mixed-convection. Owing to the stability results by Chen and Chung [START_REF] Chen | The linear stability of mixed convection in a vertical channel flow[END_REF], the governing equations are written in transient form in order to capture possible transitions to unsteady flows. With the z-axis pointing upwards and the origin of coordinates placed at the center of the inlet section, the conservation equations based on a fully elliptic model are

Governing equations

x z H u = 0 T = T 0              u = 0 ∂w ∂z = 0 ∂T ∂z = 0 D T h , v = 0
∇.v = 0 (1) ρ 0 ( ∂v ∂t + ∇.(v ⊗ v)) = -∇(p + ρ 0 gz) + µ 0 ∇ 2 v + ρ 0 gβ(T -T 0 )e z ( 2
)
∂T ∂t + ∇.(vT ) = a 0 ∇ 2 T ( 3
)
where ρ 0 , µ 0 and a 0 are the fluid density, dynamic viscosity and thermal diffusivity evaluated at T 0 .

Boundary and initial conditions

The initial condition is a fluid at rest at uniform temperature T 0 . The boundary conditions are written as follows:

• Along the vertical isothermal walls:

v = 0, T = T h at x = ±D/2 and for 0 < z < H (4)
where D and H are the channel width and height, respectively. • At the inlet (z = 0) and outlet (z = H) sections for -D/2 < x < D/2:

(a) natural convection

u = 0, ∂w ∂z = 0, p = -ρ 0 w 2 2 , T = T 0 at z = 0 u = 0, ∂w ∂z = 0, p + ρ 0 gH = 0, ∂T ∂z = 0 at z = H (5)
(b) forced or mixed convection b1 : fixed flow rate (w 0 prescribed)

u = 0, w = w 0 , T = T 0 at z = 0 u = 0, ∂w ∂z = 0, p + ρ 0 gH = 0, ∂T ∂z = 0 at z = H (6)
b2 : fixed inlet total pressure and outlet pressure (p s )

u = 0, ∂w ∂z = 0, p = -ρ 0 w 2 2 , T = T 0 at z = 0 u = 0, ∂w ∂z = 0, p + ρ 0 gH = p s , ∂T ∂z = 0 at z = H (7)
With boundary conditions Eq. 7 applied to forced or mixed convection, a pressure drop (p s < 0) is prescribed at the outlet section (it could be as well an increase in pressure at the inlet section). Therefore, the boundary conditions are inlet temperature and outlet pressure, zero z-derivative for the vertical velocity component and zero horizontal velocity-component. It should be emphasized that these flow B.C. differ somehow from those used in most of previous works.

Let us emphasize that our results achieved using Eq. 6 or Eq. 7 are almost identical, provided the average velocities w(0) are the same. The choice of the particular boundary condition Eq. 6 may lead to wrong physical interpretations, because the imposed velocity w 0 is completely uncorrelated to the velocity induced by the natural convection flow. In contrast, the use of Eq. 7 preserves the coherence of the fluid flow with respect to free convection (Eq. 5) as it will be shown in what follows.

The governing equations are generally cast in dimensionless form by using either the set of dimensionless variables 3 and P r = ν 0 /a 0 , the channel aspect ratio, the Reynolds, Grashof and Prandtl numbers, respectively. The relative strength of the buoyancy force is then characterized either by the Richardson number, Ri = Gr/Re 2 , for the first set of dimensionless variables or by the product Ri × Re = Gr/Re for the second set. Some controversy still exist about the best pertinent parameter, e.g. Ri or Ri × Re.

X = x/D h , Z = z/D h , θ = (T -T 0 )/(T h - T 0 ), U = u/w 0 , W = w/w 0 and P = p/ρ 0 w 2 0 or the set X, Z = z/(D h Re), θ, U = D h u/(ν 0 ), W and P where D h = 2D is the hydraulic diameter.

It follows that the problem formulation depends on four parameters:

A = H/D, Re = w 0 D h /ν 0 , Gr = gβ(T h -T 0 )D 3 h /ν 2 0 = 8Gr H /A

Heat transfer

The total heat transfer rate released by two walls of depth L and height H, denoted as Q 2w , is calculated as follows:

Q 2w = 2Lk 0 ∫ H 0 ∂T ∂x x=D/2 dz = 2h(LH)(T h -T 0 ) (W ) ( 8
)
where L is the length of the channel in the y-direction, h is the mean heat transfer coefficient defined by the right equality in Eq. 8. For natural or mixed convection, a mean Nusselt number expressed by

N u 2w = Q 2w /[2(LH)(k 0 ∆T /D)
] may be introduced [START_REF] Bar-Cohen | Thermally optimum spacing of vertical, natural convection cooled, parallel plates[END_REF].

The enthalpy flux removed from the two heated surfaces by the stream is

Q en = ρ 0 c p S c [(wT ) H -(wT ) 0 ] (W ) (9)
where (wT ) z is the average of (wT ) over the channel cross-section of area S c = DL. By using again 2(LH)(k 0 ∆T /D) as a reference heat flux, we can define a Nusselt number

N u en = Q en /[2(LH)(k 0 ∆T /D)].
It should be emphasized that Q 2w is always greater than Q en because the heat losses by axial diffusion through the channel inlet are not included into the definition of Q en . When the wall temperatures are kept fixed, Q 2w and Q en merge provided that the Reynolds number is large enough in order that axial diffusion may be assumed negligibly small. The difference Q cond = Q 2w -Q en yields the effect of axial diffusion which can be written as

Q cond = k 0 ∫ D/2 -D/2 ∂T ∂z z=0 dx = 2k 0 ∫ H 0 ∂T ∂x x=D/2 dz -ρ 0 c p ∫ D/2 -D/2 [w(x, H)T (x, H) -w(x, 0)T 0 ]dx (10
) For negligible axial heat diffusion, equation (8) reads:

2k 0 L ∫ H 0 ∂T ∂x x=D/2 dz = 2h(HL)(T h -T 0 ) = ṁc p (T b (H) -T 0 ) (11)
where ṁ = ρ 0 w 0 (DL) and T b (H) = (wT ) H /w 0 is the bulk temperature at the outlet section. Such an approximation is one of the foundations introduced for establishing the analytical solutions reported in Ref. [START_REF] Bejan | Convection Heat Transfer[END_REF][START_REF] Bar-Cohen | Thermally optimum spacing of vertical, natural convection cooled, parallel plates[END_REF][START_REF] Bejan | The optimal spacing of parallel plates cooled by forced convection[END_REF].

Numerical method and validation

The mass, momentum and energy conservation equations were spatially discretized by using the collocated finite volume method presented in [START_REF] Chénier | Numerical results using a colocated finite volume scheme on unstructured grids for incompressible fluid flows[END_REF][START_REF] Touazi | Simulation of natural convection with the collocated clustered finite volume scheme[END_REF], applied here for structured meshes. The key features of this method are:

the normal fluxes to cell faces are defined by a consistent two-points approximation,

the variables transported by the fluid flow are expressed on cell faces in a centered manner, not by using a linear interpolation at the faces, the "pressure gradient" is defined by the dual expression of the velocity divergence and therefore, it can be non-consistent with a gradient.

An implicit second-order Euler scheme was adopted for time derivatives at time t = (n + 1)∆t, with an implicit treatment of the diffusion terms and an Adams-Bashforth extrapolation procedure for the transport terms. A stabilization method is necessary to prevent from the onset of checkerboard oscillations associated with collocated schemes. This is performed through the velocity-pressure decoupling which is handled by a projection method [START_REF] Goda | A multistep technique with implicit difference schemes for calculating two-or three-dimensional cavity flows[END_REF].

Attention was also paid to the temperature discontinuities at the corner of the inlet section (i.e. at x = ±D/2, z = 0). In the case of pure conduction, the analytical heat fluxes along the vertical walls are indeed infinite. Therefore, the wall heat fluxes calculated numerically, Q 2w , cannot converge when refining the meshes (see [START_REF] Sun | Optimal plate spacing for mixed convection from an array of vertical isothermal plates[END_REF] for details). To avoid this problem, it is then necessary to smooth the thermal boundary conditions in order to recover some regularity. To this end, we have substituted the uniform temperature imposed at the walls by a temperature distribution which decreases abruptly towards T 0 close to z = 0. The steep function retained in this work writes:

T p (z) = T 0 + (T h -T 0 ) tanh(Γ z H ) (12)
For example, by choosing Γ = 1000, the wall temperature T p (z) for a channel height H = 1 m varies from 0.5T h at z ≃ 0.55 mm to 0.9T h at z ≃ 1.5 mm and then becomes larger than 0.99T h for z > 2.65 mm. We checked that such a steep variation in T p (z) does change significantly the temperature field while it allow second order convergence of the wall heat fluxes, Q 2w .

The velocity components, the temperature field and the pressure correction were calculated by using the Bi-Conjugate Gradient Stabilized (BCGS) method, preconditioned by an incomplete LU decomposition. Owing to the axial symmetry of the flow, the problem was solved on the half-width of the channel. The faces of the mesh x f i and z f k are defined as follows:

x f i D/2 = i n x - c x 2π sin ( i n x π ) , 0 ≤ i ≤ n x z f k H = tanh ( c z ( k n z -1
))

tanh(c z ) + 1, 0 ≤ k ≤ n z (13)
The coefficients c x and c z define the mesh refinements in the entrance region and along the isothermal wall. By choosing c x = 1 and c z = 1.82, the size ratios between the largest cell and the smallest cell are equal to 3 and 10 in horizontal and vertical directions, respectively. The time integration was performed with the Courant-Friedrich-Levy number kept fixed to CF L = 0.5. The asymptotic behavior of the numerical solution was investigated using a Richardson extrapolation on one set of the three meshes defined by n x ∈ {10, 20, 40} with n z = 20n x . This method leads to the extrapolated value of f which is written

f Extrap. = f n x + c n x (n x ) s ( 14
)
where c n x is a coefficient depending on n x and f n x is the numerical result on grid n x × n z . For sufficiently fine meshes, the value of s must tend towards the order of consistency of the numerical scheme (i.e. s = 2), f Extrap. must then become independent of the mesh. For natural convection flow, the mesh 20 × 400 defined for a half width of the channel provides accurate results within 3 digits (Tab. 1). This grid size is then retained for all computations.

Results

The influence of the dynamical boundary conditions at the inlet and outlet sections is discussed by considering air entering at T 0 = 300 K into a vertical channel with height walls H = 1 m at uniform hot temperature T h = 320 K (or H = 1.5 m, just for allowing comparisons with previously published works [START_REF] Hammou | Laminar mixed convection of humid air in a vertical channel with evaporation or condensation[END_REF][START_REF] Desrayaud | Flow reversal of laminar mixed convection in the entry region of symmetrically heated, vertical plate channels[END_REF]). The channel width was D = 2 cm or D = 3 cm. With the air properties taken at T 0 = 300 K (ρ 0 = 1.176 kg.m -3 , µ 0 = 1.85 10 -5 P a.s, c p,0 = 1006

J.kg -1 .K -1 , k 0 = 0.0261 W.m -1 .K -1
) the Grashof number is then Gr H = 2.64 10 9 H 3 . Channel heights less than about two meters are thus considered as maximum in order to keep a laminar flow.

For pure natural convection (Tab. 1), the computations yield an average velocity w cn = 0.38 m/s and a pressure drop at the inlet section (Eq. 5) p cn,in = -0.087 P a. The heat flux transferred along one of the isothermal walls is Q 1w = Q 2w /2 = 65.2 W by unit of depth length (direction perpendicular to the plane of Fig. 1). These values are the basis for the discussions of mixed convection results presented in what follows.

Before getting further into the result section, let us recall that flow solutions are similar provided that the average velocities, imposed or resulting from computation, are identical. However, when a constant flow rate is prescribed, counterintuitive results or surprising phenomena may arise. Such situations, developed in details in Sec.4.1, are now briefly exemplified. By using Eq. 6, astonishing flow recirculations are highlighted close to inlet section. Their occurrences are simply due to the incompatibility between the imposed flow rate and the necessity of feeding the dynamic boundary layers produced by the buoyancy along the vertical heated walls. Similarly from a purely thermal viewpoint, it is also expected that the buoyant force increases the velocities in the boundary layers and then improves the heat transfer at the solid walls. However this behavior is not observed when the inlet velocity is prescribed.

The next two subsections are then devoted to the detailed analysis of the fluid flows and heat transfers when using boundary conditions 6 and 7. For both sets of boundary conditions, comparisons between natural, forced and mixed convection flows are performed. These comparisons emphasize clearly that only the pressure boundary conditions make sense.

Constant flow rate

We are first considering prescribed inlet velocity in the range 0.1 m/s ≤ w 0 ≤ 1 m/s in order to satisfy the assumption of laminar flow (Re = 2500 for w 0 = 1 m/s). For pure forced convection, it results that the difference in pressure between the inlet and outlet sections is within ∆p = 0.061 P a and ∆p = 0.96 P a.

The ∆p-variations versus the inlet velocity are shown in Fig. 2 for forced convection and for mixed convection. Since the outlet pressure has been fixed to p s = 0 (see Eq. 6), that implies ∆p = p(0). For forced convection, the inlet pressure is always positive while it is negative for mixed convection if the inlet velocity is smaller than the inlet velocity corresponding to pure natural convection (w nc = 0.38 m/s). Therefore, the pressure and buoyancy forces are opposite if w 0 < w nc .

Let us now consider computations based on the dimensionless form of the governing equations. In most of the papers dealing with aiding mixed convection, the Reynolds and Grashof numbers were arbitrarily fixed: the relative importance of one convective mode was discussed by considering only the value of Ri or Ri×Re with the same length scale for Re and Gr. For the present computations (D = 2 cm), it is found that Ri = 2.62 10 -2 /w 2 0 (w 0 in m/s) and Ri × Re = 66.4/w 0 . Therefore, either the first or the second criterion (e.g Ri > O(1) [START_REF] Bejan | Convection Heat Transfer[END_REF] or Ri × Re > O( 103 ) [START_REF] Padet | Minimization of entropy production in fully developed mixed thermal convection[END_REF]) leads to assume that natural convection dominates for w 0 less than about 0.15 m/s, a value smaller than w cn = 0.38 m/s. For the practical case considered, the average velocity 0.15 m/s can only be reached for a negative pressure gradient (Fig. 2), namely for a "forced convection" opposing to the natural convection. Thus, concluding that natural convection is dominating has no physical meaning since the pressure gradient plays a dominant part into the flow rate.

The streamlines plotted in Fig. 3 shows another view of the results linked to a fixed flow rate smaller than the natural inlet flow rate. For this case, we consider a channel height H = 1.5 m and a channel width D = 3 cm in order to check that the present results are in full agreement with those published in [START_REF] Hammou | Laminar mixed convection of humid air in a vertical channel with evaporation or condensation[END_REF][START_REF] Desrayaud | Flow reversal of laminar mixed convection in the entry region of symmetrically heated, vertical plate channels[END_REF][START_REF] Laaroussi | Effects of variable density for film evaporation on laminar mixed convection in a vertical channel[END_REF].

Figure 3 raises obviously the question of the origin of existence of flow recirculations within the entry region for aiding mixed convection in vertical channels. As it is well established nowadays, the separation bubble (the axis of which being the channel axis) cannot be predicted when using parabolized formulations. On the other hand, elliptical models may lead to such flow topologies provided that the inlet velocity is prescribed without any reference to the natural convection velocity. The case shown in Fig. 3 corresponds to a clear recirculating flow, easily reproducible by using the simulation conditions reported in the legend of this figure. The conclusion is that the flow field referenced as "buoyancy-assisted mixed flow" is similar to the flow predicted for "opposing mixed convection" with p s = 1.4 P a (∆p < 0).

Figure 4 shows the variations of the vertical velocity component along the centerline (e.g. w(0, z)) for the three convective modes. is the maximum flow velocity at any z-section and it increases upstream until the flow regime is fully established (before z = H if w 0 < 0.2 m/s in the case considered). For mixed convection, the buoyancy force produces maximum in velocities near the walls. Therefore, w(0, z) for mixed convection is lower than for forced convection when the flow rate is prescribed, the difference increasing upstream with w 0 . If w 0 < w nc (w nc being the average velocity for natural convection) the decrease in w(0, z) at the bottom of the channel indicates a possible onset of flow recirculation (see Fig. 3). For w 0 = 0.1 m/s and 0.2 m/s, a stagnant zone characterized by a decrease in the axial velocity spreads within the lower half-region. This stagnation zone becomes more and more confined close to the channel inlet as the difference between w 0 and the natural convection velocity increases.

Figure 5 displays the variation in the total heat flux transferred at one of the channel walls (Q 1w = Q 2w /2) as function of the inlet velocity, the horizontal dashed line being for pure natural convection. As can be seen, mixed convection with a prescribed inlet velocity lower than that for pure natural convection leads to a lower heat transfer rate. In addition, the very small differences between the results for forced and mixed convection clearly show that the buoyancy force does not assist significantly the pressure force. For the lowest inlet velocity considered here (w 0 = 0.1 m/s), the Péclet number is P e ≈ 259. Therefore, axial conduction is negligibly small (Q cond ≪ Q 2w , Eq. 8 and 10).

Pressure driven mixed convection

When a pressure drop is prescribed at the outlet section (p s < 0), with an inlet total pressure fixed, the pressure and the velocity profile at the inlet as well as the flow rate are not a priori known. Therefore, it is not possible to calculate the Reynolds and Richardson numbers until the end of the computations. However, the flow rate is necessary greater than that for pure natural convection. For the practical case investigated (H = 1 m, D = 2 cm, ∆T = 20 K), the maximum pressure (at the center of the inlet section) for natural convection is p(0, 0) = -0.143 P a. The inlet pressure for mixed convection should be lower than this value since the flow rate is assumed to increase (p(0) ≈ -ρ 0 w 2 /2). Figure 6 shows the axial variation of pressure along the channel axis. For natural convection the pressure minimum is at z = 0.24 m while it is located at the outlet section for forced convection. For mixed convection, the pressure curves exhibit a minimum, whose location is indicated by an arrow in Fig. 6. Therefore, ∂p/∂z is zero at this z-coordinate where the buoyancy force compensates viscosity and advection. When the pressure drop at the outlet section exceeds 0.6 P a, the pressure decreases monotonously from the inlet to the outlet section, indicating the predominance of forced convection.

The average velocity is displayed in Fig 7 versus the outlet pressure drop. In contrast to what has been found for prescribed inlet velocity, w for mixed convection is systematically greater than for forced convection: buoyancy force assists thus the pressure force whatever p s < 0 and, we can effectively speak ) clearly exhibits the large differences in the wall heat flux according to the type of boundary conditions. It should be noted that the forced convection curve does not tend to zero as p s → 0 because the problem turns then into a pure conduction problem. From Eq. 10 it is found that Q cond = 4.8 W for p s = 0.

Let us now discuss the above results in terms of Ri or Ri×Re (Re and Gr being both based on the same length scale, i.e D h ). When increasing the pressure drop, Ri and Ri × Re decrease since the flow rate is increased. From p s = -0.1 P a to p s = -1.5 P a, Ri decreases from 0.134 to 0.022 while Ri×Re decreases from about 151 to about 57. The criterion for the transition from dominant forced convection to dominant natural convection being Ri > O(1) [START_REF] Bejan | Convection Heat Transfer[END_REF] or Ri × Re > O( 103 ) [START_REF] Padet | Minimization of entropy production in fully developed mixed thermal convection[END_REF], it could be assumed that natural convection plays a negligible role for the results discussed herein. In fact, the heat flux at one of the vertical walls is Q 1w ≈ 65 W for pure natural convection. That substantially differs from Q 1w ≈ 75 W for mixed convection when p s = -0.2 P a (Q 1w ≈ 50 W for forced convection). Therefore, any of these two criteria based on approximate analyses appears to be relevant. The temperature distribution along the channel axis is plotted in Fig. 9 according to the pressure at the channel outlet section. Small increases in the axial temperature indicate that the channel flow regime turns into a boundary layer type regime. Asymptotic solutions for assisting mixed convection along a vertical isothermal plate may be thus useful. The present results show that it is not the case, at least for laminar flows (i.e. |p s | < 1.5 P a). In addition, temperature as well as velocity profiles do not exhibit boundary layer-type behavior. When p s is decreased from p s = 0 (pure natural convection) to p s = -1.5 P a (dominant forced convection), Figure 9 indicates a significant decrease in the axial temperature which is closely linked to an increase in heat transfer at the walls. Since the temperature increases up to the outlet section, an established velocity field cannot be predicted, as shown in Fig. 10. This figure exhibits also that the dynamical entry length augments with the Reynolds number. However, the usual criterion for defining the dynamical and thermal entry lengths for forced convection in a flat plate channel (z d ≈ 4 10 -4 Re D h and z t /D h ≈ 0.01Re D h P r, respectively) are not met for the largest pressure drop considered.

Conclusion

Mixed-convection for assisting flows occurring in isothermal flat-plate channels was studied by solving numerically the system of conservation equations written in elliptical form. The computations were carried out by using a finite-volume method implemented in an in-house code. The spatial order of convergence of the method and the accuracy of the results have been checked thanks to the Richardson extrapolation method and regularized inlet thermal conditions in order to properly solve cases of low flow rates linked to small plate spacings and/or small temperature differences.

As an example, we have considered air at ambient temperature flowing in vertical channels (height H = 1 m or H = 1.5 m, width D = 2 cm or D = 3 cm). For a maximum temperature difference compatible with the Boussinesq approximation (∆T = 20 K), leading to Gr H = 2.65 10 9 H 3 , the entire laminar mixed-convection regime has been covered by changing either the inlet velocity or the outlet pressure drop. Re D h was varied from Re D h ≈ 200 to Re D h ≈ 2800. Although most of the results have been discussed in their dimensional forms, we have also analyzed the reliability of scaling predictions based on dimensionless numbers.

The results for a prescribed inlet velocity w 0 show that flow recirculations may appear in the entry region when the inlet velocity is lower than that for pure natural convection. In that case, the pressure force opposes to the buoyancy force and the problem is, in fact, similar to that of opposing natural convection with a downward bulk flow. In addition, we have shown that to prescribe w 0 yields no significant difference between the wall heat fluxes calculated either for forced convection or for mixed convection.

The modeling of mixed convection based on a pressure drop at the outlet section (p s < 0) and a total pressure at the inlet section leads to a completely different analysis of the results, with a better accordance with the common physical sense. When imposing an outlet pressure drop, natural convection assists effectively forced convection and the flow rate as well as the wall heat fluxes are systematically increased over those found for pure forced convection.

Fig. 1

 1 Fig.1The physical model of mixed convection in a vertical channel. The full boundary conditions are written in Eqs. (5-7) for the three cases investigated.

Fig. 2

 2 Fig.2Variations of the pressure difference between inlet and outlet sections according to the prescribed inlet velocity for forced, mixed and natural convection.

Fig. 3

 3 Fig. 3 Flow recirculation for mixed convection (a) prescribed inlet velocity , w 0 = 0.05 m/s (b) overpressure at the outlet section ps = 1.4 P a (D = 3 cm, H = 1.5m, ∆T = 20 K).

Fig. 4

 4 Fig. 4 Vertical velocity variations along the channel axis for forced and mixed convection according to the prescribed inlet velocity. The natural convection curve is also shown (D = 2 cm, H = 1 m, ∆T = 20 K).

  Figure 8 shows that natural convection dominates for the smallest values of -p s and that forced convection dominates for the largest values of -p s . The comparison between Fig. 5 (B.C. 6) and Fig. 8 (B.C. 7

Fig. 5

 5 Fig. 5 Total heat flux transferred at one of the hot walls for forced and mixed convection according to the prescribed inlet velocity. The value for natural convection is also shown (D = 2 cm, H = 1 m, ∆T = 20 K).

Fig. 6

 6 Fig. 6 Axial variation of pressure along the channel axis according to the outlet pressure drop (D = 2 cm, H = 1 m, ∆T = 20 K). The arrows indicate the positions (in m) of the pressure minimum.

Fig. 7 Fig. 8 Fig. 9

 789 Fig. 7 Variations of the average velocity according to the pressure drop at the outlet section for forced and mixed convection (D = 2 cm, H = 1 m, ∆T = 20 K).

Fig. 10

 10 Fig. 10 Variations of the vertical velocity component along the channel axis according to the pressure drop at the outlet section for natural and mixed convection (D = 2 cm, H = 1 m, ∆T = 20 K).

5

  Total heat flux transferred at one of the hot walls for forced and mixed convection according to the prescribed inlet velocity. The value for natural convection is also shown (D = 2 cm, H = 1 m, ∆T = 20 K). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 6 Axial variation of pressure along the channel axis according to the outlet pressure drop (D = 2 cm, H = 1 m, ∆T = 20 K). The arrows indicate the positions (in m) of the pressure minimum. . . . 14 7 Variations of the average velocity according to the pressure drop at the outlet section for forced and mixed convection (D = 2 cm, H = 1 m, ∆T = 20 K). . . . . . . . . . . . . . . . . . . . . . . . 15 8 Variations of the heat flux transferred at one of the walls according to the pressure drop at the outlet section for forced and mixed convection (D = 2 cm, H = 1 m, ∆T = 20 K). . . . . . . . . . . . 15 9 Variations of the temperature along the channel axis according to the pressure drop at the outlet section for natural and mixed convection (D = 2 cm, H = 1 m, ∆T = 20 K). . . . . . . . . . . . . . 16 10 Variations of the vertical velocity component along the channel axis according to the pressure drop at the outlet section for natural and mixed convection (D = 2 cm, H = 1 m, ∆T = 20 K). . . . . . . . 17

  Re Reynolds number based on D h , Re = w 0 D h /ν 0 Ri Richardson number, Ri = Gr/Re 2 S c area of the channel cross section, S c = DL [m 2 ]

	k	thermal conductivity [W m -1 K -1 ]
	L	channel length in the spanwise direction [m]
	ṁ	mass flow rate [kg s -1 ]
	n x , n z	numbers of grid points in x-and z-directions
	p	pressure [P a]
	p s	pressure at the outlet section [P a]
	P r	Prandtl number, P r = ν 0 /a 0
	Q	heat flux, [W ]
	Q en	enthalpy heat flux, [W ]
	Q 2w	convective heat flux along the two channel walls, [W ]
	t	time [s]
	T	temperature [K]
	u, w	velocity components [m s -1 ]
	x, z	coordinates [m]
	Greeks
	β coefficient of thermal expansion, β = 1/T τ dimensionless time
	Subscripts
	a, b	analytical solutions
	h	hot wall
	H	quantity based on channel height
	nc	natural convection
	0	inlet section
	Superscripts
	-	averaged quantity
	*	dimensionless quantity
	c z 1 Introduction stretching parameters, Eq. 13 and 13
	D	plate spacing [m]
	D h	hydraulic diameter, D h = 2D [m]
	g	gravitational acceleration [m s -2 ]
	Gr H	Grashof number based on H,
		Gr H = gβ 0 ∆T H 3 /ν 2 0
	h	heat transfer coefficient [W m -2 K -1 ]
	H	channel height [m]

0 [K -1 ] ∆T temperature difference, ∆T = (T h -T 0 ) [K] µ dynamic viscosity [P a.s] ν kinematic viscosity [m 2 s -1 ] ρ density [kg m -3 ] θ dimensionless temperature ratio, θ = (T -T 0 )/∆T

Table 1

 1 Convective heat flux, enthalpy heat flux and mean velocity according to the mesh refinement (nz = 20nx), for D/2 = 1 cm and H = 1 m. Order s of the numerical scheme and Richardson-extrapolated values (see Eq. 14).

	nx	10	20	40	s fExtrap.
	Q 2w (W ) 130.39 130.31 130.30 2.66 130.30
	Qen (W ) 130.32 130.22 130.20 2.01 130.19
	wnc(m/s) 0.3840 0.3833 0.3831 2.01 0.3830

List of Tables