Julien Bouyssier

Charles Pierre

Franck Plourabou

MATHEMATICAL ANALYSIS OF PARALLEL CONVECTIVE EXCHANGERS WITH GENERAL LATERAL BOUNDARY CONDITIONS USING GENERALISED GRAETZ MODES

Keywords: February, 2012 Heat and mass transfer, convection diffusion, variational formulation, Graetz mode decomposition, mixed formulation, convective exchangers

We propose a mathematical analysis of parallel convective exchangers for any general but longitudinally invariant domains. We analyse general Dirichlet or Neumann prescribed boundary conditions at the outer solid domain. Our study provides general mathematical expressions for the solution of convection/diffusion problems. Explicit form of generalised solutions along longitudinal coordinate are found from convoluting elementary base Graetz mode with the applied sources at the boundary. In the case of adiabatic zero flux counter-current configuration we recover the longitudinally linearly varying solution associated with the zeroth eigenmode which can be considered as the fully developed behaviour for heat-exchangers. We also provide general expression for the infinite asymptotic behaviour of the solutions which depends on simple parameters such as total convective flux, outer domain perimeter and the applied boundary conditions. Practical considerations associated with the numerical precision of truncated mode decomposition is also analysed in various configurations for illustrating the versatility of the formalism. Numerical quantities of interest are investigated, such as fluid/solid internal and external fluxes.

1. Introduction 1.1. Applicative context. Heat exchangers are ubiquitous in many industrial processes where heat is to be to recovered or, on the contrary, disposed, from one fluid onto another. Applications might be associated with heating or cooling systems, but can also involve other processes such as pasteurisation, crystallisation, distillation, concentration or separation of some substances [START_REF] Dušan | Fundamentals of heat exchanger design[END_REF][START_REF] Ho | Improvement in performance on laminar counterflow concentric circular heat exchangers with external refluxes[END_REF][START_REF] Ho | Double-pass Flow Heat Transfer In A Circular Conduit By Inserting A Concentric Tube For Improved Performance[END_REF]. Similarly mass exchangers are also important either in natural biological organs such as kidney and/or biotechnological applications such as devices for continuous extra corporeal blood purification associated with hemo-dialysis [START_REF] Gostoli | Mass transfer in a hollow fiber dialyzer[END_REF], hemo-filtration or extracorporeal oxygenation.

For both mass or hear exchangers, the exchange takes place from coupled convection/diffusion processes without any direct contact between the input and the output fluids, for obvious contamination purposes. Many industrial examples of such devices are possible to find such as radiators, condensers, evaporators, air pre-heaters, cooling towers as well as extra corporeal membrane oxygenation and blood micro-filters [START_REF] Chen | Microfluidic chip for blood cell separation and collection based on crossflow filtration[END_REF]. Also found in exchangers as a generic, although not systematic, common feature, is parallel flow design configuration. This is the class of exchangers that we are going to consider in this paper, with the hypothesis that there is no longitudinal variation of the fluid-velocity along the exchanger axis.

In previous contributions a similar Graetz decompositions for solving exchange configurations has already been used [START_REF] Li'n | Laminar counteflow parallel-plate heat exchangers: Exact and approximate solutions[END_REF][START_REF] Ho | An analytical study of laminar counterflow double-pass heat exchangers with external refluxes[END_REF][START_REF] Weigand | The extended graetz problem with piecewise constant wall heat flux for pipe and channel flows[END_REF]25]. Graetz decomposition applied to stationary convective exchanges problems provides elegant and compact solutions. Furthemore, the obtained familly of exponentially decaying modes also permits to set up a hierarchy of modes in fully developped configurations [START_REF] Ho | An analytical study of laminar counterflow double-pass heat exchangers with external refluxes[END_REF]. Nevertheless, there is a number of limitations that have preclude the systematic and intensive use of such decomposition in more realistic configurations [START_REF] Bouyssier | Analysis of parallel convective exchangers with versatile inlet/outlet conditions[END_REF] It is much straighforward to used them to convective dominated situations (where the Péclet number is large) [START_REF] Li'n | Laminar counteflow parallel-plate heat exchangers: Exact and approximate solutions[END_REF][START_REF] Ho | An analytical study of laminar counterflow double-pass heat exchangers with external refluxes[END_REF] (2) Their used has been restricted to two dimensional [START_REF] Nunge | Analysis of heat or mass tranfer in some contercurrent flows[END_REF][START_REF] Nunge | An analytical study of laminar counter flow double-pipe heat exchangers[END_REF][START_REF] Li'n | Laminar counteflow parallel-plate heat exchangers: Exact and approximate solutions[END_REF][START_REF] Telles | Solutions of the extended Graetz problem[END_REF] or possibly concentric [START_REF] Telles | Solutions of the extended Graetz problem[END_REF][START_REF] Ho | An analytical study of laminar counterflow double-pass heat exchangers with external refluxes[END_REF]25] configurations.

(3) They have been used only for constant or piecewise constant prescribed Dirichlet lateral temperature profiles [START_REF] Weigand | The extended graetz problem with piecewise constant wall heat flux for pipe and channel flows[END_REF]25] or homeneous Neumann adiabatic lateral boundary conditions [START_REF] Ho | An analytical study of laminar counterflow double-pass heat exchangers with external refluxes[END_REF]. (4) Input/output conditions where generally considered as prescribed uniform temperature [START_REF] Li'n | Laminar counteflow parallel-plate heat exchangers: Exact and approximate solutions[END_REF][START_REF] Ho | An analytical study of laminar counterflow double-pass heat exchangers with external refluxes[END_REF][START_REF] Weigand | The extended graetz problem with piecewise constant wall heat flux for pipe and channel flows[END_REF]25] whitout considering the possible coupling with Inlet/Outlet conditions in realistic configurations.

Limitation [START_REF] Bouyssier | Analysis of parallel convective exchangers with versatile inlet/outlet conditions[END_REF] can be overcomed from considering the proper set of orthogonal modes as first noticed by [START_REF] Papoutsakis | The extended Graetz problem with Dirichlet wall boundary conditions[END_REF], so that axially diffusion can also be included in coupled co-current or counter-current problems. Nevertheless, limitation (2) has only been overcomed recently in [START_REF] Pierre | Numerical analysis of a new mixed formulation for eigenvalue convection-diffusion problems[END_REF] for longitudinally infinite exchangers with homogeneous lateral Dirichlet boundary conditions or for finite exchangers [START_REF] Fehrenbach | The Generalized Graetz problem in finite domains[END_REF], again for the homogeneous Dirichlet boundary condition at the outer solid surface.

In a effort to obtain general formulation to reach realistic applicative configurations, limitations (3) and ( 4) are still pending. The contribution of this work is to remove restriction [START_REF] Chen | Microfluidic chip for blood cell separation and collection based on crossflow filtration[END_REF]. In the following we provide the necessary mathematical theory and numerical implementation to permit the use of generalysed Graetz decomposition for any general lateral boundary conditions. The goal of a next paper will be to remove restriction [START_REF] Fehrenbach | The Generalized Graetz problem in finite domains[END_REF].

In order to realise that restriction (3) is important in application, it is interesting to mention that the heat pipes literature has considered a number of different lateral boundary conditions [START_REF] Shah | Laminar Flow Forced Convection in Ducts[END_REF][START_REF] Papoutsakis | The extended Graetz problem with Dirichlet wall boundary conditions[END_REF][START_REF] Telles | Solutions of the extended Graetz problem[END_REF][START_REF] Mahulikar | A new classification for thermal development of fluid flow in a circular tube under laminar forced convection[END_REF][START_REF] Weigand | The effect of wall conduction for the extended Graetz problem for laminar and turbulent channel flows[END_REF]25].

As described in [START_REF] Shah | Laminar Flow Forced Convection in Ducts[END_REF] one founds, in heat pipes, lateral boundary conditions with uniform profile (uniform Dirichlet) in transverse and longitudinal directions, uniform profile along longitudinal direction only, radiative boundary conditions, prescribed uniform flux (uniform Neumann), or exponentially varying profile along longitudinal direction. It is interesting to mention that exponentially varying lateral boundary conditions in the longitudinal direction permits to take into account some convection/diffusion coupling between the fluid and the solid, as described by fully developed Graetz modes which are indeed exponentially decaying solutions [START_REF] Graetz | Uber die Wärmeleitungsfähigkeit von Flüssigkeiten[END_REF][START_REF] Nusselt | Die abhängigkeit der wärmeübergangszahl von der rohrlänge[END_REF].

Hence, each operating conditions thus necessitates a case-specific theoretical treatment without any generally theoretical framework which could describe the complete coupling between convection arising inside the fluid coupled with the diffusion inside the solid.

The purpose of this contribution is to provide such theoretical framework for any general set of prescribed temperature profile or applied flux around the exterior solid boundary of the exchanger. This work is an extension of two previous contributions which have permitted to generalise standard Graetz eigenmodes decomposition to any, possibly complicated, configuration in the transverse direction, whilst longitudinally invariant. One considerable advantage of the developed formalism is to provide a two-dimensional formulation of a fully tri-dimensional problem.

In [START_REF] Pierre | Numerical analysis of a new mixed formulation for eigenvalue convection-diffusion problems[END_REF] longitudinally infinite exchangers are considered with homogeneous Dirichlet boundary conditions. In a second contribution, the extension of two-dimensional formulation to finite configurations [START_REF] Fehrenbach | The Generalized Graetz problem in finite domains[END_REF] has also been considered, again for the homogeneous Dirichlet boundary condition at the outer solid surface.

In this paper, we generalise previous approaches for the very general case of any applied Dirichlet or Neumann boundary condition. Of particular interest is the Neumann case for which, two distinct class of problems emerges from its mathematical properties, as discussed in [START_REF] Li'n | Laminar counteflow parallel-plate heat exchangers: Exact and approximate solutions[END_REF]. Given the fluid density ρ i , the heat capacity c i and the flux Qi of a fluid i, one should distinguish the case where the total heat capacity flow rate Σ i ρ i Qi c i over all fluid i involved in the exchanger is equal or distinct from zero. It is, for exemple, distinct from zero for convective heat pipes for which convection is driving the heat from inlet to outlet, and in that case, any generalised Graetz mode is exponentially varying in the longitudinal direction. In the case where the total heat capacity flow rate Σ i ρ i Qi c i = 0 as encountered in balanced counterflow heat exchangers, we show that an additional linearly longitudinally varying mode has to be considered as already discussed in two-dimensional convection-dominated configurations [START_REF] Li'n | Laminar counteflow parallel-plate heat exchangers: Exact and approximate solutions[END_REF].

1.2. Physical problem and state of the art. This paper considers stationnary convection diffusion in a tube, i.e. in a domain Ω × I with Ω ⊂ R 2 a smooth bounded domain and I ⊂ R an interval (possibly unbounded). A point M ∈ Ω × I has for coordinates M = (x, z) with x = (x 1 , x2 ) ∈ Ω and z ∈ I.

Inside the tube a moving fluid convects a passive tracer, whilst it diffuses inside the immobile solid part. Two physical assumptions are considered. Firstly the fluid velocity ṽ in the tube is independent of z and directed along the z-direction: ṽ(x, z) = ṽ(x)e z . Moreover we adopt the natural convention that ṽ = 0 in the solid part of the domain Ω. Secondly the thermal conductivity k is assumed to be isotropic and independent of z: k = k(x) ∈ R (anisotropic conductivity however could be considered with the condition that e z is one principal direction of the conductivity tensor).

In this setting, the stationary convection-diffusion equation for the temperature T on Ω × I reads [START_REF] Bouyssier | Analysis of parallel convective exchangers with versatile inlet/outlet conditions[END_REF] div( k ∇ T ) + k∂ 2 z T = ρcṽ∂ z T , where div = div x and ∇ = ∇ x, ρ the fluid density, c the heat capacity, ṽ the fluid velocity and k the conductivity.

In the following, we discuss a dimensionless form of this convection/diffusion equation. Following previous conventions in the heat exchanger litterature, considering fluid pipes of radius R, the dimensionless coordinates (x, z) are defined as x = (x 1 /R, x2 /R) and z = z/R so that the pipe radius is unity. Furthermore, the conductivity k is non-dimensionalised by the fluid conductivity kf so that k = k/ kf is unity in the fluid. The dimensionless temperature T is obtained from considering a reference temperature T 0 , T = T /T 0 . The fluid non-dimensional velocity is defined such that, [START_REF] Gostoli | Mass transfer in a hollow fiber dialyzer[END_REF] v = P e 2 ṽ V where V is the average fluid velocity in the fluid pipe, and P e is the Péclet number usually considered as ( 3)

P e = ρcV (2R) k f = V (2R) α ,
where α = k f /ρc is the fluid thermal diffusivity. For notation simplification we will only consider in the following one fluid type and only one pipe radius R. This is nevertheless not a restriction of the presented results, which can easily be generalised to more complex configurations involving tubes of different diameters and different fluids. Note also that we have included the P e number in the definition of dimensionless velocity v for notation simplification. Using dimensionless formulation Equation (1) then reads (4) div(k∇T ) + k∂ 2 z T = v∂ z T, where similarly div = div x and ∇ = ∇ x . The dimensionless velocity v and conductivity k satisfy mandatorily the two following properties:

(5) v ∈ L ∞ (Ω) and 0 < k m ≤ k(x) ≤ k M , x ∈ Ω,
no additional regularity assumptions on v and k are needed. Problem (4) has been reformulated in [START_REF] Pierre | Numerical analysis of a new mixed formulation for eigenvalue convection-diffusion problems[END_REF], as briefly stated below. On the Hilbert space

H = L 2 (Ω) × [L 2 (Ω)] 2
we consider the unbounded operator

A A : D(A) ⊂ H → H , D(A) = H 1 (Ω) × H div (Ω) (6) ∀ (s, q) ∈ D(A) , A(s, q) = (k -1 vs -k -1 div q, k∇s).
In matricial notation, the operator A display the form:

A = k -1 v -k -1 div k∇ 0 Lemma 1.1. Let z ∈ I → ψ(z) = (T (z), q(z)) ∈ H be a differentiable function so that ∀ z ∈ I, ψ(z) ∈ D(A) and d dz ψ(z) = Aψ(z),
then T is a solution of (4) and k∇T = ∂ z q.

Here z → T (z) is once differentiable in L 2 -norm, so ∂ 2 z T only has a weak (distribution) sense and T is a solution to (4) in the same weak sense. A strong solution to (4) is recovered if additionally z → q(z) is differentiable in H div (Ω). Lemma 1.1 is the starting point in [START_REF] Pierre | Numerical analysis of a new mixed formulation for eigenvalue convection-diffusion problems[END_REF][START_REF] Fehrenbach | The Generalized Graetz problem in finite domains[END_REF] to derive solutions to (4) using the spectral properties of A.

In all the sequel Q will denote the total flow flux across Ω

Q = Ω v(x)dx.
The case where no net flux arises, i.e. Q = 0, is singular for the Neumann case. Note that, given the definition of the dimensionless velocity v in (2) this condition also reads in dimensional form as a zero total heat capacity flow rate Σ i ρ i Qi c i = 0 for different fluid i flowing in different pipes with flux Qi . Hence, in the following we will considerer the "zero flux" dimensionless condition Q = 0 whithout mentionning any more that it corresponds, in fact, to an adiabatic zero flux countercurrent configuration, as discussed in [START_REF] Li'n | Laminar counteflow parallel-plate heat exchangers: Exact and approximate solutions[END_REF] for convection dominated regimes. This case is singular because of the existence of a non trivial element in the kernel of A denoted Ψ 0 : [START_REF] Ho | Improvement in performance on laminar counterflow concentric circular heat exchangers with external refluxes[END_REF] Ψ 0 = (1, k∇u 0 ), div(k∇u 0 ) = v and k∇u 0 • n = 0 on ∂Ω.

When Q = 0, u 0 ∈ H 1 (Ω) is well defined (up to a constant) since the compatibility condition Ω div(k∇u 0 )dx = ∂Ω k∇u 0 • ndl = 0 = Ω vdx.
This induces the existence of a special solution T 0 satisfying zero flux condition on ∂Ω,

T 0 (x, z) = C 1 (u 0 (x) + z) + C 2 , with C 1 , C 2 ∈ R.
1.3. Summary of the paper. This paper is organised in two main parts. In the first part, we present theoretical results and derive analytical solutions to problem (4). The second part provides numerical illustration of the obtained results of the first part whilst the efficiency of the analytical solutions here derived to describe heat exchanges in tubes.

In the first theoretical part we start in section 2 with a spectral analysis of the operator A in (6) either considering a Dirichlet or Neumann type boundary condition on ∂Ω. We provide an extension of the results in [START_REF] Pierre | Numerical analysis of a new mixed formulation for eigenvalue convection-diffusion problems[END_REF] dedicated to the Dirichlet case. This extension in particular shows that in the Neumann case the physics of the problem depends on the value of the total flux Q = Ω vdx. The case Q = 0 is quite singular and moreover of great interest in our applicative context: it corresponds to counter-current configuration heat exchanger devices. In all cases our main result in theorem 2.3 states that A is diagonal over a (complete) orthogonal basis. The spectrum moreover is made of a double infinite sequence of eigenvalues going both to +∞ and -∞, each sequence corresponding either to the upstream (z < 0) or downstream (z > 0) region descriptions.

Using these results we derive the solutions to problem (4) for non-homogeneous boundary conditions of Dirichlet and Neumann type. These solutions are studied both for an infinite (Ω × R) or a semi infinite (Ω × R + ) domain in sections 3 and 4 respectively. These solutions are obtained as separate variable series: the variation in the transverse (i.e. Ω) direction is given by the operator A eigenfunctions and the longitudinal variation is explicitly given by a simple integral transformation involving both the boundary data (treated as a source term) and the eigenvalues of A.

Numerical results are given in the last section 5. The analytical solutions in the two previous sections can be approximated by truncating their series expansion and by approximating the eigenvalues and eigenfunctions. This approximation is performed with two-dimensional finite element setting as presented in section 5.1. A first axy-symmetric test case is presented in section 5.2 whose purpose is to validate the method. Finally the method is developed to describe the fluid/solid heat exchange for two more complex configurations: a periodic set of parallel pipes and a counter current heat exchanger.

Spectral analysis

The Hilbert space H is equipped with the scalar product:

∀ Ψ i = (f i , p i ) ∈ H, i = 1, 2, (Ψ 1 |Ψ 2 ) H = Ω f 1 f 2 k(x)dx + Ω p 1 • p 2 k -1 (x)dx,
that is equivalent with the canonical scalar product on H (i.e. taking k=1) thanks to property [START_REF] Graetz | Uber die Wärmeleitungsfähigkeit von Flüssigkeiten[END_REF] on k.

With this definition the operator A satisfies:

∀ Ψ i = (s i , q i ) ∈ D(A), i = 1, 2, ( 8
) AΨ 1 |Ψ 2 H = Ψ 1 |AΨ 2 H + ∂Ω s 1 q 2 • ndl - ∂Ω s 2 q 1 • ndl,
with n the unit normal on ∂Ω pointing outwards Ω.

D(A D ) = H 1 0 (Ω) × H div (Ω) , D(A N ) = H 1 (Ω) × H 0 div (Ω), with H 0 div (Ω) = {q ∈ H div (Ω), q
• n = 0 on ∂Ω} . The operators A D and A N clearly have dense domains in H. Using the property (8) they also are symmetric:

∀ Ψ 1 , Ψ 2 ∈ D(A) : (AΨ 1 |Ψ 2 ) H = (Ψ 1 |AΨ 2 ) H , either with A = A D or A = A N .
Theorem 2.1. The two operators A D and A N in definition 2.1 are self adjoint.

With theorem 2.1 We have:

H = Ker (A D ) ⊥ Ran (A D ) = Ker (A N ) ⊥ Ran (A N ) ,
we now characterise these spaces.

Corollary 2.2. In the Dirichlet case

Ker (A D ) = {(0, q), q ∈ H div (Ω) and div q = 0} , [START_REF] Ho | Improvement in performance on laminar counterflow concentric circular heat exchangers with external refluxes[END_REF],

Ran (A D ) = (f, k∇s), f ∈ L 2 (Ω), s ∈ H 1 0 (Ω) .

In the Neumann case let us consider

K N = (0, q), q ∈ H 0 div (Ω) and div q = 0 , R N = (f, k∇s), f ∈ L 2 (Ω), s ∈ H 1 (Ω) . if Q = 0: Ker (A N ) = K N and Ran (A N ) = R N , if Q = 0: let us consider Ψ 0 defined in
Ker (A N ) = K N ⊥ Span (Ψ 0 ) , R N = Ran (A N ) ⊥ Span (Ψ 0 ) . We always have H = K N ⊥ R N .
In addition to these symmetry properties, the spectrum of the operators A D and A N can be fully characterised. Let us denote Sp ⋆ (A) := Sp(A) -{0}, the spectrum of the operator A without the singleton {0}. Either for A = A D or A = A N , Sp ⋆ (A) displays the same double sequence structure Sp ⋆ (A) = (λ n ) n∈Z ⋆ with:

(9) -∞ ←- n→+∞ λ n ≤ • • • ≤ λ 1 < 0 < λ -1 ≤ • • • ≤ λ -n -→ n→+∞ +∞.
Negative values of λ ∈ Sp ⋆ (A) will be referred to as downstream modes whereas positive values will be referred to as upstream modes.

Theorem 2.3. The two operators A -1 D : Ran (A D ) → H and A -1 N : Ran (A N ) → H are compact. We have Sp ⋆ (A D ) = (λ D n ) n∈Z ⋆ and Sp ⋆ (A N ) = (λ N n ) n∈Z ⋆ , these two sequences satisfy [START_REF] Ho | Double-pass Flow Heat Transfer In A Circular Conduit By Inserting A Concentric Tube For Improved Performance[END_REF]. Elements of Sp ⋆ (A D ) and of Sp ⋆ (A N ) are eigenvalues of finite order, the associated eigenvectors (Ψ D n ) n∈Z ⋆ and (Ψ N n ) n∈Z ⋆ form a Hilbert basis (with orthogonal vectors of norm 1) of Ran (A D ) and Ran (A N ) respectively. We will denote

Ψ D n = (T D n , q D n ) and Ψ N n = (T N n , q N n ) the
ψ ∈ Ran (A) iff ψ 2 H = n∈Z ⋆ | (ψ|Ψ n ) H | 2 < +∞, ( 10
)
ψ ∈ D(A) ∩ Ran (A) iff ψ 2 D(A) := n∈Z ⋆ |λ n (ψ|Ψ n ) H | 2 < +∞, (11)
and

ψ D(A) is a norm equivalent to the H 1 (Ω) × H div (Ω) norm.
Proof of theorem 2.1. The Dirichlet case has already been proved in [START_REF] Pierre | Numerical analysis of a new mixed formulation for eigenvalue convection-diffusion problems[END_REF]. The proof in the Neumann case follows the same arguments and is detailed here.

We have A N = A 0 + V with A 0 : (s, q) ∈ D(A N ) → (-k -1 div q, k∇s) and V : (f, p) → (k -1 vf, 0). Both A N and A 0 are symmetric with dense domains and V is bounded on H. Using the Kato-Rellich theorem (see e.g. [START_REF] Reed | Methods of modern mathematical physics. II. Fourier analysis, self-adjointness[END_REF] p. 163), the self-adjointness of A 0 implies the self-adjointness of A N . To prove the self-adjointness of A 0 , let us show that A 0 + i has range H (see e.g. [START_REF] Reed | Methods of modern mathematical physics. I. Functional Analysis[END_REF]).

Let (f, p) ∈ H. Using the Lax Milgram theorem, there exists a unique q ∈ H 0 div (Ω) so that for all χ ∈ H 0 div (Ω): ( 12)

Ω (q • χ + div q div χ) k -1 dx = Ω -ip • χk -1 -f div χ dx.
More precisely the bilinear form on the left is clearly coercive on H 0 div (Ω) thanks to property (5) and the linear form on the right is continuous on H 0 div (Ω). We introduce the function s so that is -k -1 div q = f . Let us prove that Ψ = (s, q) ∈ D(A N ). With div q = k(is -f ) we get with [START_REF] Nunge | An analytical study of laminar counter flow double-pipe heat exchangers[END_REF]:

∀ χ ∈ H 0 div (Ω), - Ω s div χdx = Ω k -1 (p -iq) • χdx, so that in distribution sense ∇s = k -1 (p -iq) ∈ [L 2 (Ω)] 2
. It follows that s ∈ H 1 (Ω) and therefore Ψ ∈ D(A N ). Finally we have is -k -1 div q = f and k∇s + iq = p which means that (f, p) = A 0 Ψ + iΨ and so Ran (A 0 + i) = H.

Proof of corollary 2.2. We prove the Neumann case only. If Ψ = (s, q) ∈ D(A N ) satisfies A N Ψ = 0 then k∇s = 0 and therefore s is a constant. Now we have sv -div(q) = 0 with s ∈ R.

In case Q = 0, by integrating sv -div(q) = 0 over Ω and using the divergence formula we get s Ω vdx = sQ = 0 and so s = 0. It follows that div q = 0, as a result Ker (A N ) = K N in this case.

In case Q = 0, s can be a non zero constant and in this case Ψ ∈ Span (Ψ 0 ) with Ψ 0 defined in [START_REF] Ho | Improvement in performance on laminar counterflow concentric circular heat exchangers with external refluxes[END_REF].

Thus Ker (A N ) = K N ⊥ Span (Ψ 0 ) in this case. we clearly have R N ⊂ K ⊥ N . Since H = Ker (A N ) ⊥ Ran (A N ), to end the proof of corollary 2.2 it suffices to show that H = K N R N . Let (f, p) ∈ H. Let s ∈ H 1 (Ω) so that, ∀ ϕ ∈ H 1 (Ω) , Ω k∇s • ∇ϕdx = Ω p • ∇ϕdx
The function s ∈ H 1 (Ω) is defined up to an additive constant by the Lax Milgram theorem. Let q = p -k∇s: for all ϕ ∈ H 1 (Ω) we then have Ω q • ∇ϕdx = 0 and so div q = 0 in distribution sense. Thus q ∈ H div (Ω) and with the Green formula we get for all ϕ ∈ H 1 (Ω), Ω q • ∇ϕdx = ∂Ω ϕq • ndl = 0 and finally q ∈ H 0 div (Ω). We eventually can decompose (f, p) = (f, k∇s) + (0, q) so that H = K N R N which ends the proof.

Proof of theorem 2.3. The compactness has already been proved in the Dirichlet case with additional regularity assumptions on k in [START_REF] Pierre | Numerical analysis of a new mixed formulation for eigenvalue convection-diffusion problems[END_REF]. We here give the proof for k ∈ L ∞ (Ω) and for the Neumann case only, the proof in the Dirichlet case is similar and a bit simpler. We simply denote A N = A along this proof.

We consider (f n , p n ) ∈ Ran (A) a bounded sequence in H-norm. We consider the unique (s n , q n ) ∈ D(A) ∩ Ran (A) so that A(s n , q n ) = (f n , p n ), one has to prove that the sequence (s n , q n ) is compact. Compacity is in H-norm, equivalent with the L 2 -norm.

Let us first prove that the sequence (s n ) is compact in L 2 (Ω) We firstly have that k∇s n = p n and thus ∇s n L 2 is bounded. The Poincaré Wirtinger inequality then ensures that s n -c n L 2 is bounded with c n = Ω s n dx the mean value of s n . If we prove that (c n ) is bounded, we obtain that s n L 2 is also bounded and then that (s n ) is bounded in H 1 (Ω): this will prove that (s n ) is compact in L 2 (Ω) using the Rellich-Kondrachov compactness theorem. If Q = 0, we have vs n -div(q n ) = kf n and integrating over Ω we get that

Ω v(s n -c n )dx + c n Ω vdx - Ω div(q n )dx = Ω v(s n -c n )dx + c n Q = Ω kf n dx,
because q n satisfies a zero flux condition on ∂Ω. Since f n and s n -c n are bounded in L 2 -norm, we get that c n is bounded. If Q = 0, with corollary 2.2 we have the additional constraint here that ((

s n , q n )|Ψ 0 ) H = 0 = Ω ks n dx + Ω q n • ∇u 0 dx. With vs n -div(q n ) = kf n we have Ω q n • ∇u 0 dx = Ω -div(q n )u 0 dx = Ω (kf n -vs n )u 0 dx. The orthogonality constraint then gives Ω ((k -vu 0 )s n + kf n ) dx = 0 . From this last equality we get, c n Ω (k -vu 0 )dx = Ω (s n -c n )(vu 0 -k)dx - Ω ku 0 f n dx.
The right hand side is bounded. Using that div(k∇u 0 ) = v, the pre-factor satisfies Ω (k -vu 0 )dx = Ω k(1 + ∇u 0 • ∇u 0 )dx and is non zero. It follows that c n is bounded.

Let us now prove that (q n ) is compact in [L 2 (Ω)] 2 . With corollary 2.2 we firstly have that q n = k∇u n with u n ∈ H 1 (Ω). We can moreover impose that Ω u n dx = 0. We have that k∇u n L 2 is bounded and with the Poincaré-Wirtinger inequality it follows that u n is bounded in H 1 (Ω) thus compact in L 2 (Ω). We also have that

g n := div(k∇u n ) = vs n -kf n is bounded in L 2 (Ω).
Up to sub-sequence extractions, we then can assume that:

s n -→ s strongly in L 2 (Ω), u n -→ u strongly in L 2 (Ω), q n -⇀ q weakly in [L 2 (Ω)] 2 , g n -⇀ g weakly in L 2 (Ω).
We have to prove that indeed q n -→ q strongly in [L 2 (Ω)] 2 .

Let us first prove that u ∈ H 1 (Ω) with k∇u = q and that q ∈ H 0 div (Ω) with div q = g.

For all test function χ ∈ [C ∞ c (Ω)] 2 , Ω k -1 q • χdx = lim n Ω k -1 q n • χdx = lim n Ω ∇u n • χdx = -lim n Ω u n div χdx = - Ω u div χdx,
in the distribution sense, this means that ∇u = k -1 q ∈ L 2 (Ω) i.e. u ∈ H 1 (Ω) and k∇u = q. Now for all ϕ ∈ H 1 (Ω), since q n ∈ H 0 div (Ω), we have

Ω ϕgdx = lim n Ω ϕg n dx = lim n Ω ϕ div q n dx = -lim n Ω ∇ϕ • q n dx = - Ω ∇ϕ • qdx.
This first ensures that div q = g ∈ L 2 (Ω) in the sense of distributions, so that q ∈ H div (Ω). Now that we nave q ∈ H div (Ω), using the Green formula we moreover have for all ϕ ∈ H 1 (Ω):

Ω ϕgdx = - Ω ∇ϕ • qdx = Ω ϕ div qdx - ∂Ω ϕq • ndl = Ω ϕgdx - ∂Ω ϕq • ndl..
The boundary integral is always equal to zero and thus q ∈ H 0 div (Ω). Finally, we now can conclude than q nq L 2 → 0:

q n -q 2 L 2 = Ω kk -1 (q n -q) • (q n -q)dx ≤ k M Ω k -1 (k∇u n -k∇u) • (q n -q)dx,
using inequality [START_REF] Graetz | Uber die Wärmeleitungsfähigkeit von Flüssigkeiten[END_REF]. With the Green formula it follows that,

q n -q 2 L 2 ≤ k M Ω ∇(u n -u) • (q n -q)dx = -k M Ω (u n -u)(g n -g)dx ≤ k M u n -u L 2 g n -g L 2 , ,
with the Cauchy-Schwartz inequality. We conclude to q nq

L 2 → 0 because u n -u L 2 → 0 and g n -g L 2 is bounded.
We proved that A -1 : Ran (A) → Ran (A) is compact. It is moreover self adjoint by theorem 2.1 and injective by construction. By Hilbert Schmidt theorem there exists an orthogonal Hilbert basis of Ran (A) made of eigenvectors of A -1 , moreover 0 is the only limit points of the associated sequence of (non-zero) eigenvalues.

This proves the remaining part of theorem 2.3 excepted the particular structure (9) displayed by the eigenvalues. To prove this we have to show that the Rayleigh coefficients (ψ, Aψ) are bounded neither above nor below for ψ ∈ D(A) and ψ H = 1, we already know that they are unbounded. We have with ψ = (s, q),

(ψ, Aψ) = Ω vs 2 dx + 2 Ω q • ∇sdx.
The first term on the right is clearly bounded, then the second one is unbounded and moreover changes of sigh when performing the transformation (s, q) → (-s, q). This second term then is unbounded above and below.

Solutions on infinite domains

We here consider the case I = R. Being given a function f : R → R, we look for a solution T to (4) on Ω × R for the two following problems.

Dirichlet problem:

T

(x, z) = f (z) for x ∈ ∂Ω, ( 13
)
Neumann problem: k∇T (x, z) • n = f (z) for x ∈ ∂Ω. ( 14
)
We firstly draw the basic ideas to derive solutions to these two problems. Rigorous statements of the solutions are given in the following sub-sections, they follow from these preliminary formal results given below.

With lemma 1.1, we search for a solution ψ(z) = (T (z), q(z)) to dψ/dz = Aψ under the form,

ψ(z) = n∈Z ⋆ d n (z)Ψ D n or ψ(z) = n∈Z ⋆ d n (z)Ψ N n ,
for the Dirichlet or Neumann cases respectively.

Formally differentiating the sums we get,

d dz ψ(z) = n∈Z ⋆ d ′ n (z)Ψ D n or d dz ψ(z) = n∈Z ⋆ d ′ n (z)Ψ N n , so that d ′ n = Aψ|Ψ D n H or d ′ n = Aψ|Ψ N n H respectively. Using (8) we obtain, d ′ n = λ D n d n + ∂Ω T (z)q D n • ndl or d ′ n = λ N n d n - ∂Ω q(z) • nT N n dl.
For the Dirichlet problem T (z) = f (z) on ∂Ω. For the Neumann problem we have k∇T = ∂ z q so that q(z)

• n = F (z) with F a primitive of f on ∂Ω.
Let us introduce the coefficients α n given either for the Dirichlet or the Neumann problems by:

(15) α D n = - 1 λ D n ∂Ω q D n • ndl, α N n = 1 λ N n ∂Ω T N n dl.
Finally the functions d n satisfy,

d ′ n = λ D n d n -λ D n α D n f (z) or d ′ n = λ N n d n -λ N n α N n F (z)
, respectively for the Dirichlet or Neumann cases.

In the sequel the function d n are sought under the form,

d n (z) = α D n f (z) + α D n c n (z) or d n (z) = α N n F (z) + α N n c n (z), thus with the function c n solution of, c ′ n = λ D n c n -f ′ or c ′ n = λ N n c n -f
, respectively for the Dirichlet or Neumann cases. In the Dirichlet case, the solution Ψ is then sought under the form,

Ψ(z) = f (z) n∈Z ⋆ α D n Ψ D n + n∈Z ⋆ α D n c n (z)Ψ D n ,
whereas in the Neumann case it reads,

Ψ(z) = F (z) n∈Z ⋆ α N n Ψ N n + n∈Z ⋆ α N n c n (z)Ψ N n .

Let us first characterise the functions n∈Z

⋆ α D n Ψ D n and n∈Z ⋆ α N n Ψ N n . 3.1. The coefficients α n .
Lemma 3.1. The coefficients α n defined in (15) satisfy:

(16) n∈Z ⋆ α D n Ψ D n = ϕ D , n∈Z ⋆ α N n Ψ N n = ϕ N , with ϕ D ∈ Ran (A D ) ∩ D(A) uniquely determined by, (17) ϕ
D = (1, k∇u D ), u D ∈ H 1 0 (Ω)
and Aϕ D = 0, and with ϕ N ∈ Ran (A N ) ∩ D(A) uniquely defined by,

ϕ N = (s N , k∇u N ), k∇u N • n = 1 on ∂Ω (18)
and

Aϕ N = 0 if Q = 0 aΨ 0 if Q = 0, a ∈ R ,
for Ψ 0 defined in [START_REF] Ho | Improvement in performance on laminar counterflow concentric circular heat exchangers with external refluxes[END_REF]. In the particular case Q = 0, the constraint ϕ N ∈ Ran (A N ) implies that ϕ N |Ψ 0 H = 0. The Bessel inequality [START_REF] Mahulikar | A new classification for thermal development of fluid flow in a circular tube under laminar forced convection[END_REF]

ensures that n∈Z ⋆ |α D n | 2 < +∞ and n∈Z ⋆ |α N n | 2 < +∞. Mean- while, ϕ D / ∈ D(A D ) and ϕ N / ∈ D(A N )
, we also have with relation [START_REF] Nunge | Analysis of heat or mass tranfer in some contercurrent flows[END_REF]

n∈Z ⋆ |λ D n α D n | 2 = +∞ and n∈Z ⋆ |λ N n α N n | 2 = +∞. Remark 3.1.
Let us precise the definition of the particular functions ϕ D and ϕ N .

• For the Dirichlet case ϕ D = (1, k∇u D ), the function u D is determined by the equation div(k∇u D ) = v and u D = 0 on the boundary ∂Ω. • For the Neumann case when Q = 0, ϕ N = (s N , k∇u N ) and s N is a constant, equal to P/Q with P the perimeter of Ω. The second component is defined by Q div(k∇u N ) = P v and k∇u N • n = 1 on ∂Ω. This equation is well posed as long as Q = 0 and u N is defined up to an additive constant. In the sequel we will fix this constant by imposing that:

(19) Q Ω vu N dx = P Ω kdx.
• For the Neumann case when Q = 0, let us consider the two constants a, b ∈ R:

a Ω (vu 0 -k)dx = P, (20) b
Ω (vu 0 -k)dx = a Ω u 0 (2k -vu 0 )dx + ∂Ω u 0 dl, (21)
with P the perimeter of the domain Ω and u 0 defined in [START_REF] Ho | Improvement in performance on laminar counterflow concentric circular heat exchangers with external refluxes[END_REF]. In this case the function Proof. We first prove that there exists a unique ϕ D ∈ Ran (A D ) ∩ D(A) satisfyng [START_REF] Dušan | Fundamentals of heat exchanger design[END_REF]. The condition Aϕ D = 0 imposes div(k∇u D ) = v which equation has a unique solution u D ∈ H 1 0 (Ω) . We now prove that there exists a unique ϕ N ∈ Ran (A N ) ∩ D(A) satisfyng [START_REF] Reed | Methods of modern mathematical physics. I. Functional Analysis[END_REF]. Let first Q = 0. The condition Aϕ N = 0 imposes s N ∈ R and div(k∇u N ) = s N v. With k∇u N • n = 1 on ∂Ω the compatibility condition has the form Ω dl = P = s N Q and sets the constant s N = P/Q. The equation Q div(k∇u N ) = P v with boundary condition k∇u N • n = 1 on ∂Ω is well posed and defines u N up to an additive constant. Let now Q = 0. The condition Aϕ N = aΨ 0 first imposes that k∇s N = ak∇u 0 and so s N = au 0 + b with b ∈ R. It secondly imposes the equation,

ϕ N = (s N , k∇u N ) satisfies s N = au 0 + b.
v au o + b -div(k∇u N ) = ak.
With k∇u N • n = 1 on ∂Ω the compatibility condition for this equation is independent on b and reads:

Ω (avu o + bv -div(k∇u N )dx = a Ω vu 0 dx -P = a Ω kdx.
This relation uniquely determines the value of a as stated in [START_REF] Shah | Laminar Flow Forced Convection in Ducts[END_REF]: note that a is well defined since

Ω (vu 0 -k)dx = Ω (div(k∇u 0 )u 0 -k)dx = -Ω (k∇u 0 • ∇u 0 + k)dx = 0.
Setting a to this value characterises (for any given value of b) the function u N up to an additive constant. We got ϕ N = (au 0 + b, k∇u N ) ∈ D(A) ∩ R N satisfying Aϕ N = aΨ 0 . We still have to fulfill the constraint ϕ N ∈ Ran (A N ), in the present case Q = 0 (see corollary 2.2) we only have to impose ϕ N |Ψ 0 H = 0:

ϕ N |Ψ 0 H = Ω k(au 0 + b)dx + Ω k∇u N • ∇u 0 dx = Ω k(au 0 + b)dx - Ω div(k∇u N )u 0 dx + ∂Ω u 0 dl = Ω k(au 0 + b)dx + Ω ka -(au 0 + b)v u 0 dx + ∂Ω u 0 dl = b Ω (k -vu 0 )dx + a Ω u 0 (2k -vu 0 )dx + ∂Ω u 0 dl = 0
This sets the value of b to [START_REF] Telles | Solutions of the extended Graetz problem[END_REF].

With ϕ D and ϕ N so defined, let us prove [START_REF] Pierre | Generalised graetz problem: analytical solutions for concentric or parallel configurations[END_REF]. For this, since ϕ D ∈ Ran (A D ) and ϕ N ∈ Ran (A N ), it suffices to prove that

ϕ D |Ψ D n H = α N n and ϕ N |Ψ D n H = α N n .
We do it in the Neumann case only. Using property (8) we get,

λ N n ϕ N |Ψ N n H = ϕ N |AΨ N n H = Aϕ N |Ψ N n H + ∂Ω T N n k∇u N • ndl.
We always have

Aϕ N |Ψ N n H = 0: either Aϕ N = 0 if Q = 0 or Aϕ N = aΨ 0 ∈ Ran (A N ) ⊥ if Q = 0.
With k∇u N • n = 1 on ∂Ω we obtain the result:

ϕ N |Ψ N n H = ∂Ω T N n dl/λ N n = α N n . 3.2. The Dirichlet problem. We simply denote here λ n = λ D n , Ψ n = Ψ D n and α n = α D n . We introduce the functions c n (z) for n ∈ Z ⋆ , (22) c n (z) = +∞ z f ′ (ξ)e λn(z-ξ) dξ if n < 0, c n (z) = - z -∞ f ′ (ξ)e λn(z-ξ) dξ if n > 0.
If we assume that f ′ is bounded, these functions are well defined (because λ n < 0 for n > 0 and vice versa), they also are bounded, differentiable and verify c ′ n = λ n c n -f ′ . Proposition 3.2. We assume that f ∈ C 1 (R) with f ′ bounded and we consider the mapping z ∈ R → ψ(z) = (T (z), q(z)) ∈ Ran (A D ), ( 23)

ψ(z) = f (z)ϕ D + n∈Z ⋆ α n c n (z)Ψ n ,
We have,

(24) ψ ∈ C 1 (R, H) ∩ C 0 (R, D(A)), d dz ψ = Aψ on R,
and T is the solution to the Dirichlet problem (4) [START_REF] Nusselt | Die abhängigkeit der wärmeübergangszahl von der rohrlänge[END_REF]. The regularity estimate above simply means that z → T (z) is C 1 in L 2 (Ω) and continuous in H 1 (Ω): T is a weak solution to equation ( 4), as stated in lemma, 1.1.

If we moreover assume that f ∈ C 2 (R) with f ′′ bounded we get the additional regularity,

(25) ψ ∈ C 2 (R, H) ∩ C 1 (R, D(A)). This means that z → T (z) is C 2 in L 2 (Ω), C 1 in H 1 (Ω) and that z → k∇T (z) is continuous in H div (Ω).
With this assumption T is a strong solution to equation ( 4), as stated in lemma, 1.1.

The definition of the temperature T (z) associated to ψ in ( 23) can be precised thanks to remark 3.1. We have,

T (z, x) = f (z) + n∈Z ⋆ α n c n (z)T n (x).
Moreover far-field estimates on the temperature can be derived from this expression under suitable assumptions on f . Roughly speaking, if f (z) → f (+∞) as z → +∞ then we also have T (z) → f (+∞). If we instead have a linear growth of f at +∞, then T (z) verifies a similar asymptote. This is precised in the following corollary.

Corollary 3.3. We assume f ∈ C 1 (R) with f ′ bounded. If +∞ 0
|f ′ |dz < +∞, then f has a limit in +∞ and:

T (z) -→ z→+∞ f (+∞) in L 2 (Ω).
We moreover assume f ∈ C 2 (R) with f ′′ bounded. If +∞ 0 |f ′′ |dz < +∞, then f ′ has a limit in +∞ and:

∂ z T (z) -→ z→+∞ f ′ (+∞), T (z) z→+∞ f (z) + f ′ (+∞)u D in L 2 (Ω), with u D defined in remark 3.1: div(k∇u D ) = v and u D ∈ H 1 0 (Ω).
The regularity assumption on f can be weakened. In particular jumps of f can be taken into account. We still can derive solutions when f ′ = g + δ with g continuous and bounded and δ a Dirac-type distribution. The functions c n (z) in [START_REF] Li'n | Laminar counteflow parallel-plate heat exchangers: Exact and approximate solutions[END_REF] can be defined in this framework as well as the mapping ψ in [START_REF] Weigand | The effect of wall conduction for the extended Graetz problem for laminar and turbulent channel flows[END_REF]. With such a boundary data ψ remains continuous in H but is only differentiable outside the support of δ. These properties are detailed in the corollary 3.4 Similarly jumps on f ′ can be taken into account. Taking now f ′′ = g + δ as previously, then

T (z) remains C 2 in L 2 (Ω), C 1 in H 1 (Ω) and k∇T (z) remains C 0 in H div (Ω) outside the support of δ.
These properties are detailed in the corollary 3.7.

Corollary 3.4. Assume that f (z) = ω(z) with ω(z) = 1 if z < 0 and ω(z) = 0 otherwise, so that f ′ = -δ 0 . The computation of the functions c n (z) in ( 22) leads to the following definition of ψ(z) = (T (z), q(z)):

(26) ψ(z) = ω(z)ϕ D -ω(z) n<0 α n e λnz Ψ n + (1 -ω(z)) n>0
α n e λnz Ψ n .

We have,

ψ ∈ C 0 (R, H) ∩ C ∞ (R ⋆ , D(A)) and d dz ψ = Aψ on R ⋆ ,
and T is the solution to the Dirichlet problem (4) ( 13) with f = ω. It has the following regularity:

z → T (z) is C ∞ on R ⋆ in H 1 (Ω) and z → k∇T (z) is C ∞ on R ⋆ in H div (Ω). It then is a strong solution on R ⋆ . At the origin z = 0, T is continuous in L 2 (Ω).
Proof of proposition 3.2. Let us first prove the regularity estimates for ψ in equation [START_REF] Weigand | The effect of wall conduction for the extended Graetz problem for laminar and turbulent channel flows[END_REF]. The regularity for z → f (z)ϕ D is clear. From [START_REF] Mahulikar | A new classification for thermal development of fluid flow in a circular tube under laminar forced convection[END_REF] [START_REF] Nunge | Analysis of heat or mass tranfer in some contercurrent flows[END_REF], it suffices to prove that the two series [START_REF] Dušan | Fundamentals of heat exchanger design[END_REF] of ϕ D and therefore

|λ n c n (z)| ≤ f ′ L ∞ and |c ′ n | = |λ n c n -f ′ | ≤ 2 f ′ L ∞ . The boundary condition (13) follows from ψ(z) -f (z)ϕ D = n∈Z ⋆ α n c n (z)Ψ n ∈ D(A D ) = H 1 0 (Ω) × H div (Ω). This implies T (z) -f (z) ∈ H 1 0 (Ω) by definition
T |∂Ω = f (z).
Let us now prove that dψ/dz = Aψ. On one hand, since Aϕ D = 0 we have

Aψ = A n∈Z ⋆ α n c n (z)Ψ n = n∈Z ⋆ λ n α n c n (z)Ψ n .
On the other hand,

d dz ψ = f ′ (z)ϕ D + n∈Z ⋆ α n c ′ n (z)Ψ n = f ′ (z)ϕ D + n∈Z ⋆ α n (-f ′ (z) + λ n c n (z))Ψ n = f ′ (z) ϕ D - n∈Z ⋆ α n Ψ n + Aψ.
This gives dψ/dz = Aψ with ( 16). Now assume that f ∈ C 2 (R) with f ′′ bounded. By integrating by part we get for n < 0:

c n (z) = 1 λ n f ′ (z) + 1 λ n +∞ z f ′′ (ξ)e λn(z-ξ) dξ. With c ′ n = λ n c n -f ′ we get, c ′ n (z) = +∞ z f ′′ (ξ)e λn(z-ξ) dξ. It follows that |λ n c ′ n (z)| ≤ f ′′ L ∞ . Since c ′′ n = λ n c ′ n -f ′′ we also get the second upper bound |c ′′ n (z)| ≤ 2 f ′′ L ∞ .
We have the same upper bounds for n > 0, they imply that the two series n∈Z ⋆ |λ n α n c ′ n (z)| 2 and n∈Z ⋆ |α n c ′′ n (z)| 2 are uniformly converging. With (10) [START_REF] Nunge | Analysis of heat or mass tranfer in some contercurrent flows[END_REF], this respectively ensures that

z → n∈Z ⋆ α n c n (z)Ψ n is C 1 in D(A) and C 2 in H.
Proof of corollary 3.3. We assume that

+∞ 0 |f ′ |dz < +∞. It implies that ε(z) = sup [z,+∞) |f ′ | → 0 as z → +∞ .
Let us prove that the c n (z) uniformly converges towards 0 as z → +∞. With n∈Z ⋆ |α n | 2 < +∞ this will ensure that,

n∈Z ⋆ α n c n (z)Ψ n -→ z→+∞ 0 in H, which implies that T (z) → f (+∞) in L 2 (Ω).
First for n < 0 (and so λ n > 0), we have:

|c n (z)| ≤ e λnz +∞ z |f ′ (ξ)|e -λnξ dξ ≤ ε(z)e λnz e -λnξ -λ n +∞ z = ε(z)/λ n ,
that implies uniform convergence to 0 for n < 0. Now for n > 0 (and so λ n < 0), we have for any z 0 < z:

|c n (z)| ≤ e λnz z0 -∞ |f ′ (ξ)|e -λnξ dξ + e λnz z z0 |f ′ (ξ)|e -λnξ dξ.
The second integral is easy to bound: since -λ n > 0, we have 0 < e -λnξ < e -λnz on [z 0 , z] and so:

e λnz z z0 |f ′ (ξ)|e -λnξ dξ ≤ z z0 |f ′ (ξ)|dξ ≤ +∞ z0 |f ′ (ξ)|dξ.
We now bound the first integral,

e λnz z0 -∞ |f ′ (ξ)|e -λnξ dξ ≤ f ′ ∞ e λnz e -λnξ -λ n z0 -∞ = f ′ ∞ e λn(z-z0) |λ n | ≤ f ′ ∞ e λ1(z-z0) |λ 1 | ,
using λ n < λ 1 < 0 for n > 0 and z -z 0 > 0. For a given ε > 0, we can find z 0 so that +∞ z0

|f ′ (ξ)|dξ < ε and we can find z 1 > z 0 so that for all z > z 1 we have f ′ ∞ e λ1(z-z0) /|λ 1 | < ε. It follows that |c n (z)| < 2ε for all z > z 1 and all n > 0.

In case f ∈ C 2 (R) with f ′′ bounded, it has been proved in the proof for proposition 3.2 that,

c ′ n (z) = +∞ z f ′′ (ξ)e λn(z-ξ) dξ, c ′ n (z) = - z -∞
f ′′ (ξ)e λn(z-ξ) dξ, respectively for n < 0 and n > 0. Thus the same arguments as in the previous case prove that,

n∈Z ⋆ α n c ′ n (z)Ψ n -→ z→+∞ 0 in H, so that ∂ z T → f ′ (+∞) as z → +∞. Now we have c n = c ′ n /λ n + f ′ /λ n : n∈Z ⋆ α n c n (z)Ψ n = n∈Z ⋆ α n c ′ n (z) λ n Ψ n + f ′ (z) n∈Z ⋆ α n λ n Ψ n .
The first sum converges to zero as z → +∞. The second one towards

f ′ (+∞)A -1 D ϕ D , with A -1 D ϕ D = n∈Z ⋆ α n Ψ n /λ n by definition. Let us characterise A -1 D ϕ D .
We search for (s, q) ∈ D(A D ) so that A(s, q) = ϕ D , it satisfies s ∈ H 1 0 (Ω) and k∇s = k∇u D so that s = u D . Finally we showed that n∈Z ⋆ α n c n (z)T n → f ′ (+∞)u D as z → +∞ which ends the proof.

Proof of corollary 3.4. The regularity of ψ in equation ( 26) follows from the observation that, since λ n → +∞ as n → -∞ and since n∈Z ⋆ |α n | 2 < +∞, then the series n<0 |λ k n α n e λnz | 2 is uniformly converging for z ∈ (-∞, -ε) for all ε > 0 and all k ∈ N. As a result with [START_REF] Nunge | Analysis of heat or mass tranfer in some contercurrent flows[END_REF], z ∈ (-∞, 0) → n<0 α n e λnz Ψ n is infinitely differentiable in D(A D ). The same result holds for z ∈ (0, +∞) → n>0 α n e λnz Ψ n . The continuity of ψ at the origin in H-norm follows from [START_REF] Pierre | Generalised graetz problem: analytical solutions for concentric or parallel configurations[END_REF]. The proof that dψ/dz = Aψ on R ⋆ and that T = w on ∂Ω is identical with the proof of proposition 3.2.

) c n (z) = +∞ z f (ξ)e λn(z-ξ) dξ if n < 0, c n (z) = - z -∞ f (ξ)e λn(z-ξ) dξ if n > 0, 27
For f bounded they are well defined, bounded. differentiable and verify c ′ n = λ n c n -f . Proposition 3.5. We assume that f ∈ C 1 (R) and that both f and f ′ are bounded. We introduce F a primitive of f .

Let us consider the mapping z ∈ R → ψ(z) ∈ Ran (A N ), ( 28)

ψ(z) = F (z)ϕ N + n∈Z ⋆ α n c n (z)Ψ n ,
we have

ψ ∈ C 2 (R, H) ∩ C 1 (R, D(A)) and d dz ψ = Aψ on R,
and T is a strong solution to the Neumann problem (4) [START_REF] Papoutsakis | The extended Graetz problem with Dirichlet wall boundary conditions[END_REF].

The regularity of the solution to the Neumann problem is increased of one degree comparatively to the Dirichlet case. This comes from the definition of the function c n that are defined with respect to f (equation ( 27)) in the Neumann case whereas they are defined with the help of f ′ (equation ( 22)) for the Dirichlet problem.

Another interesting difference with the Dirichlet case is that the temperature now is defined up to an additive constant and we have an infinite set of solutions T . Precisely with remark 3.1, the temperature T can be written as:

(29) T (z, x) = P Q F (z) + n∈Z ⋆ α n c n (z)T n (x),
where F is defined up to an additive constant. This was expected since any constant C is solution to (4) with homogeneous Neumann boundary condition on R × ∂Ω. We however have uniqueness for the gradient of T that describes the heat exchanges. To have a unique determination of the temperature, the constant in F has to be set. This means adding some normalisation condition on the temperature (indeed in the Dirichlet case this normalisation also is present but implicitly). Such a normalisation can be done considering the far field temperature with suitable conditions on f (roughly f → 0 at one end of the duct at least). This is precised in the following corollary Corollary 3.6. We assume as in proposition 3.5 that f ∈ C 1 (R) and that both f and f ′ are bounded. If +∞ 0

|f |dz < +∞, then F has a limit in +∞ and:

T (z) -→ z→+∞ F (+∞) P Q in H.
The constant F (+∞) then can be fixed by a normalisation condition on T at z = +∞.

If +∞ 0
|f ′ |dz < +∞, then f has a limit in +∞ and:

∂ z T (z) -→ z→+∞ f (+∞) P Q , T (z) z→+∞ F (z) + f (+∞)u N in L 2 (Ω),
with u N defined in remark 3.1: Q div(k∇u N ) = P v and k∇u N •n = 1 on ∂Ω with the normalisation condition [START_REF] Reed | Methods of modern mathematical physics. II. Fourier analysis, self-adjointness[END_REF].

Similar results could of course also be obtained with asymptotic assumptions on f at z = -∞.

Solution can also be obtained with weaker regularity on the boundary data f : precisely with f ′ = g + δ with g continuous and bounded and with δ a Dirac-type distribution. The functions c n (z) in ( 27) can be defined in this framework as well as the function ψ in (28). With such a boundary data ψ remains C 1 in H on R but is only C 2 in H and C 1 in D(A) outside the support of δ. These properties are detailed in the corollary 3.7. Transposed to the Dirichlet context this corresponds to the case f ′′ = g + δ.

Corollary 3.7. Assume that f (z) = ω(z) with ω(z) = 1 if z < 0 and ω(z) = 0 otherwise, so that f ′ = -δ 0 . The computation of the functions c n (z) in (27) in this case leads to the following definition of ψ(z) = (T (z), q(z)):

(30) ψ(z) = A -1 N ϕ N + zϕ N - n<0 α n e λnz Ψ n /λ n if z < 0 n>0 α n e λnz Ψ n /λ n if z > 0 , where A -1 N ϕ N = n∈Z ⋆ α n Ψ n /λ n is well defined since ϕ N ∈ Ran (A N ). We have, ψ ∈ C 1 (R, H) ∩ C 0 (R, D(A)) ∩ C ∞ (R ⋆ , D(A)),
and T is a solution to the Neumann problem (4) ( 14) with f = ω. The regularity result above neans:

z → T (z) is C ∞ on R ⋆ in H 1 (Ω) and z → k∇T (z) is C ∞ on R ⋆ in H div -norm. At the origin z = 0, T is C 1 in L 2 (Ω) and k∇T is continuous in L 2 (Ω).
Proof of proposition 3.5. Let ψ be given by equation ( 28). The function F in this section plays a symmetric role with f in the previous section (Dirichlet case). Here F ∈ C 2 (R) with F ′ and F ′′ bounded: with proposition 3.2, relations ( 24) and (25) hold for ψ. It only remains to prove that T satisfies [START_REF] Papoutsakis | The extended Graetz problem with Dirichlet wall boundary conditions[END_REF] We have

∂ z ψ = f (z)ϕ N + n∈Z ⋆ α n c ′ n (z)Ψ n .
In the proof of proposition 3.2 we showed that [START_REF] Reed | Methods of modern mathematical physics. I. Functional Analysis[END_REF]. With ∂ z q = k∇T we obtain the desired boundary condition [START_REF] Papoutsakis | The extended Graetz problem with Dirichlet wall boundary conditions[END_REF].

|λ n c ′ n (z)| ≤ F ′′ L ∞ , consequently n∈Z ⋆ α n c ′ n (z)Ψ n ∈ D(A N ). It follows that on ∂Ω we have ∂ z q • n = f (z)k∇u N • n = f (z) using
Proof of corollary 3.6. Replacing f by F in the proof of corollary 3.3 gives that:

• if +∞ 0 |f |dz < +∞ then n∈Z ⋆ α n c n (z)Ψ n → 0 as z → +∞, • if +∞ 0 |f ′ |dz < +∞ then n∈Z ⋆ α n c ′ n (z)Ψ n → 0 and n∈Z ⋆ α n c n (z)Ψ n → f (+∞)A -1 N ϕ N as z → +∞.
It remains to characterise (s, q) = A -1 N ϕ N . We use the determination of ϕ N in remark 3.1: k∇s = k∇u N so that s = u N + C with C a constant. We now have v(u N + C) -div(q) = kP/Q. Integrating over Ω and since q ∈ H 0 div (Ω) we get:

Q Ω v(u N + C)dx = P Ω kdx,
so that with the chosen normalisation (19) C = 0.

Proof of corollary 3.7. Since F ′ = ω ∈ L ∞ (R) we can apply the first part of proposition 3.2 so that (24) holds for ψ in equation ( 30). The regularity estimates on R ⋆ are identical to the ones n the Dirichlet case.

Let us examine the boundary condition.

d dz ψ = ϕ N - n<0 α n e λnz Ψ n if z < 0 n>0 α n e λnz Ψ n if z > 0 .
We clearly have dψ/dz -ω(z)ϕ N ∈ D(A N ) for z = 0. As a result on the boundary ∂Ω, k∇T We recall that the definition of ϕ N = (s N , k∇u N ) = n∈Z ⋆ α n Ψ n is singular in this case. As stated in lemma 3.1 and in remark 3.1, s N = au 0 + b with a and b two constants given in (20) ( 21) and u 0 defined in [START_REF] Ho | Improvement in performance on laminar counterflow concentric circular heat exchangers with external refluxes[END_REF]. We have ϕ N ∈ Ran (A N ), so that ϕ N |Ψ 0 H = 0, and Aϕ N = aΨ 0 . We also recall that Ψ 0 = (1, k∇u 0 ) and the definition of Ran (A N ) in corollary 2.2: R N = Ran (A N ) ⊥ Span (Ψ 0 ).

• n = ∂ z q • n = ω(z)k∇u N • n = ω(z) with (18).
Proposition 3.8. The results of proposition 3.5 extends to the case Q = 0 with the alternative definition of ψ : z ∈ R → R N :

ψ(z) = aG(z)Ψ 0 + F (z)ϕ N + n∈Z ⋆ α n c n (z)Ψ n , with G : R → R satisfying G ′ = F .
The case Q = 0 present singular characteristics that deserves our attention. The temperature reads:

T (z) = aG(z) + F (z)(au 0 + b) + n∈Z ⋆ α n c n (z)T n .
When compared to (29) we see that the leading term as z → ±∞ in the temperature are different: if Q = 0 it is F (z)P/Q whereas when Q = 0 it is aG(z). Assume for example that f = 0 in a neighbourhood of +∞: in this case T (z) will converge to some limit if Q = 0 whereas for Q = 0 it will present a linear growth. Similarly if f = L = 0 in a neighbourhood of +∞: in this case T (z) will present a linear growth if Q = 0 whereas for Q = 0 this growth will instead be parabolic. These two properties being consequences of corollary 3.6.

A second important difference is that the solution now is defined up to two constants, which was expected since any function of the form C 1 (z + u 0 ) + C 2 is solution of the homogeneous Neumann problem. Thus two solutions of the problem may have different gradient and then correspond to different heat exchanges. To clarify this we rewrite the temperature as,

T (z) =C 1 (z + u 0 ) + C 2 + aG(z) + F (z)(au 0 + b) + n∈Z ⋆ α n c n (z)T n , (31)
∂ z T (z) =C 1 + aF (z) + f (z)(au 0 + b) + n∈Z ⋆ α n c ′ n (z)T n .
imposing F (0) = G(0) = 0 and with C 1 , C 2 two constants. Assume that f (+∞) = 0 and +∞ 0 |f ′ |dz < +∞: physically we could say f = 0 at one end of the duct. As showed in the proof of corollary 3.6, it implies that n∈Z ⋆ α n c ′ n (z)Ψ n → 0 as z → +∞. Here F (+∞) = +∞ 0 f dz and we get the following limit,

∂ z T (z) -→ z→+∞ C 1 + a +∞ 0 f dz.
This is another important difference with the case Q = 0 where this limit would be fixed here, equal to f (+∞)P/Q. In the case Q = 0 this limit is free. We can impose the heat flux at +∞ : such an additional condition determines C 1 . With this supplementary condition, we conserve an infinite set of solutions depending on C 2 : but two such solutions have equal gradients and now correspond to equal heat exchanges.

Proof of proposition 3.8. All the regularity estimates will still hold in this case: they only depend on the c n (z) whose definition remained unchanged. The boundary condition [START_REF] Papoutsakis | The extended Graetz problem with Dirichlet wall boundary conditions[END_REF] will also be satisfied since Ψ 0 = (1, k∇u 0 ) and u 0 satisfies a zero flux condition on ∂Ω.

We then only have to prove that ∂ z ψ = Aψ. On one hand, since Aϕ N = aΨ 0 and AΨ 0 = 0 we have

Aψ = aF (z)Ψ 0 + n∈Z ⋆ λ n α n c n (z)Ψ n .
On the other hand,

d dz ψ = aF (z)Ψ 0 + f (z)ϕ N + n∈Z ⋆ α n c ′ n (z)Ψ n = aF (z)Ψ 0 + f (z)ϕ N + n∈Z ⋆ α n (-f (z) + λ n c n (z))Ψ n = f (z) ϕ N - n∈Z ⋆ α n Ψ n + Aψ.
This gives dψ/dz = Aψ with ( 16).

Solutions on semi-infinite domains

We here consider the case I = R + = (0, +∞).

Being given a function f : R + → R, we look for a solution T either to the Dirichlet problem (4) [START_REF] Nusselt | Die abhängigkeit der wärmeübergangszahl von der rohrlänge[END_REF] or to the Neumann problem (4) ( 14) on ∂Ω × R + . These two problems need an additional entry condition on Ω × {0} to be closed. The theoretical and numerical bacground to deal with Dirichlet entry condition has been developped in [START_REF] Fehrenbach | The Generalized Graetz problem in finite domains[END_REF], they will be used and briefly discussed in the next numerical section 5. A second methodology has been developped in [START_REF] Bouyssier | Analysis of parallel convective exchangers with versatile inlet/outlet conditions[END_REF] in order to deal with quite general entry condition: a mixed combination of Dirichlet, Neumann and Robin type, moreover allowing the coupling with additional tubes. Since the focus in this paper is on the lateral boundary conditions, entry condition will not be precised in this section.

4.1. The Dirichlet problem. We again simply denote λ n = λ D n , Ψ n = Ψ D n and α n = α D n . We assume that f is differentiable and that f ′ is bounded. The functions c n (z) now are defined for n ∈ Z ⋆ by,

c n (z) = +∞ z f ′ (ξ)e λn(z-ξ) dξ if n < 0, c n (z) = - z 0 f ′ (ξ)e λn(z-ξ) dξ if n > 0.
They are well defined, bounded, differentiable and satisfy c ′ n = λ n c n -f ′ . These functions only differs from the ones in equation ( 22) when n > 0, in particular we have c n (0) = 0 for n > 0.

Consider the mapping

z ∈ R + → ψ(z) = (T (z), q(z)) ∈ Ran (A N ), (32) ψ(z) = f (z)ϕ D + n∈Z ⋆ α n c n (z)Ψ n + n>0 β n e λnz Ψ n , for a given sequence (β n ) n>0 satisfying n∈Z ⋆ |β n | 2 < +∞.
When comparing (32) with the solution on infinite domain [START_REF] Weigand | The effect of wall conduction for the extended Graetz problem for laminar and turbulent channel flows[END_REF] we point out that here the functions c n (z) have an alternative definition for n > 0. We also underline that the additional term n>0 β n e λnz Ψ n ∈ Ran (A D ) is well defined since λ n < 0 for n > 0.

The results of proposition 3.2 extends in this context.

Proposition 4.1. If we assume that f ∈ C 1 (R + ) with f ′ bounded then, (33) ψ ∈ C 1 (R + , H) ∩ C 0 (R + , D(A)), d dz ψ = Aψ on R + ,
and T is a solution to the Dirichlet problem (4) [START_REF] Nusselt | Die abhängigkeit der wärmeübergangszahl von der rohrlänge[END_REF] in the weak sense stated in lemma 1.1. With remark 3.1 the temperature T reads:

T = f (z) + n∈Z ⋆ α n c n (z)T n + n>0
β n e λnz T n .

If we moreover assume that f ∈ C 2 (R + ) with f ′′ bounded we get the additional regularity,

ψ ∈ C 2 (R + , H) ∩ C 1 (R + , D(A)).
and T is a strong solution of the Dirichlet problem (4) [START_REF] Nusselt | Die abhängigkeit der wärmeübergangszahl von der rohrlänge[END_REF]. We refer to proposition 3.2 for the translation of these regularity results in term of temperature T and of heat flux k∇T .

Remark 4.1. The meaning of the supplementary term ψ 2 (z) = n>0 β n e λnz Ψ n is the following. The function ψ 2 is a solution to the homogeneous Dirichlet problem (4) (13) (ie taking f = 0, see the proof below). This is true for any sequence (β n ) n>0 satisfying n∈Z ⋆ |β n | 2 < +∞. Its role is to satisfy a prescribed entry condition at z = 0. Let us for instance consider the Dirichlet entry condition T |z=0 = E ∈ L 2 (Ω). Then the sequence (β n ) n>0 must verify, since c n (0) = 0 for n > 0:

f (0) + n<0 α n c n (0)T n + n>0 β n T n = E.
It has been proven in [START_REF] Fehrenbach | The Generalized Graetz problem in finite domains[END_REF] that such a sequence exists and is unique. A numerical method to approximate (β n ) n>0 is also developped in this work: more details will be given in the numerical section 5.1.

As in the infinite domain case, we could weaken the regularity assumptions on f , allowing in particular jumps of f or of f ′ . Results in corollary 3.4 can easily be adapted here and we will not detail this matter.

Far field estimates can be derived exactly as in corollary 3.3 (because z → n>0 β n e λnz Ψ n goes to 0 as well as all its derivatives as z → +∞). We only recall these estimates:

• if +∞ 0 |f ′ |dz < +∞, then T (z) → f (+∞) as z → +∞ in L 2 (Ω), • if +∞ 0 |f ′′ |dz < +∞, then ∂ z T (z) → f ′ (+∞) and T (z) z→+∞ f (z)+f ′ (+∞)u D as z → +∞ in L 2 (Ω) (with u D defined in remark 3.1).
Proof. Let us decompose ψ(z) = ψ 1 (z) + ψ 2 (z) with,

ψ 1 (z) = f (z)ϕ D + n∈Z ⋆ α n c n (z)Ψ n , ψ 2 (z) = n>0 β n e λnz Ψ n .
The first term ψ 1 can be analysed exactly as in the proof of proposition 3.2: in particular

ψ 1 ∈ C 1 ([0, +∞), H) ∩ C 0 ([0, +∞), D(A)) and if f ∈ C 2 ([0, +∞)) with f ′′ bounded then ψ 1 ∈ C 2 ([0, +∞), H) ∩ C 1 ([0, +∞), D(A)). Moreover ∂ z ψ 1 = Aψ 1 , and denoting ψ 1 = (T 1 , q 1 ) then T 1 = f (z) on ∂Ω.
The second term ψ 2 in turn also can be analysed exactly as in the proof of corollary 3.4 (since

n>0 |β n | 2 < +∞): ψ 2 ∈ C ∞ (R + , D(A)) and ∂ z ψ 1 = Aψ 1 . We also have ψ 2 = (T 2 , q 2 ) ∈ D(A D ) so that T 2 = 0 on ∂Ω.

The Neumann problem for

Q = 0. We simply denote λ n = λ N n , Ψ n = Ψ N n and α n = α N n .
We assume that f ∈ C 1 ([0, +∞) with both f and f ′ bounded. The functions c n (z) now are defined for n ∈ Z ⋆ by,

c n (z) = +∞ z f (ξ)e λn(z-ξ) dξ if n < 0, c n (z) = - z 0 f (ξ)e λn(z-ξ) dξ if n > 0.
They are well defined, bounded, differentiable and satisfy c ′ n = λ n c n -f . These functions only differs from the ones in equation ( 27) when n > 0, in particular we have c n (0) = 0 for n > 0.

We consider the mapping z ∈ R + → ψ(z) = (T (z), q(z)) ∈ Ran (A N ), ( 34)

ψ(z) = F (z)ϕ N + n∈Z ⋆ α n c n (z)Ψ n + n>0 β n e λnz Ψ n ,
for F ′ = f and for some sequence

(β n ) n>0 satisfying n∈Z ⋆ |β n | 2 < +∞.
This mapping can be analysed exactly as in the Dirichlet problem: ψ ∈ C 2 (R, H) ∩ C 1 (R, D(A)) and ∂ z ψ = Aψ. Moreover the temperature T is a strong solution to the Neumann problem (4) [START_REF] Papoutsakis | The extended Graetz problem with Dirichlet wall boundary conditions[END_REF]. With remark 3.1 the temperature T has the form, For instance for the Dirichlet entry condition T |z=0 = E, the sequence (β n ) n>0 must satisfy,

T = P Q F (z) + n∈Z ⋆ α n c n (z)T n + n>0 β n e λnz T
P Q F (0) + n<0 α n c n (0)T n + n>0 β n T n = E.
Solutions can be obtained with weaker regularity assumptions on the boundary data f , as in corollary 3.7. Far field estimates can be derived exactly as in corollary 3.6; because z → n>0 β n e λnz Ψ n goes to 0 as well as all its derivatives as z → +∞.

4.3.

The Neumann problem for Q = 0. We keep here the notations and definition of the previous section 4.2.

We introduce a primitive F of f and a primitive G of F . We consider the mapping

z ∈ R + → ψ(z) = (T (z), q(z)) ∈ R N , (35) ψ(z) = aG(z)Ψ 0 + F (z)ϕ N + n∈Z ⋆ α n c n (z)Ψ n + n>0 β n e λnz Ψ n ,
with a defined in equation ( 20) and for some sequence

(β n ) n>0 satisfying n∈Z ⋆ |β n | 2 < +∞. As in section 4.2, if f ∈ C 1 ([0, +∞) with both f and f ′ bounded then ψ ∈ C 2 (R, H) ∩ C 1 (R, D(A)) and ∂ z ψ = Aψ.
Moreover the temperature T is a strong solution to the Neumann problem (4) [START_REF] Papoutsakis | The extended Graetz problem with Dirichlet wall boundary conditions[END_REF]. The temperature T reads,

T = aG(z) + F (z)(au 0 + b) + n∈Z ⋆ α n c n (z)T n + n>0 β n e λnz T n ,

Numerical results

Our original purpose is the description of the heat exchanges in heating pipes and heat exchangers. In the previous two sections we derived analytical solutions for the temperature and the heat flux on such devices. In this section we provide numerical illustrations and we analyse the efficiency for these analytical solutions to describe the heat exchanges between a fluid flowing in a tube and the surrounding solid. We consider a tube-like geometry Ω × I either with I = R or I = R + . The fluid is assumed to flow in a circular duct. Following adimentionalisation process in section 1.1, the external radius of the duct is taken equal to one. The velocity v has the Poiseuille profile:

v(x) = P e(1 -x -x 0 2 ),
with x 0 the centre of the circular duct and with P e the Péclet number. The thermal diffusivity is taken as homogeneous, k(x) = 1. We will consider four test cases. The two first ones have an axisymmetric configuration. The third one is a periodic configuration describing a collection of parallel circular ducts. The last test case is a counter current configuration where Q = 0. 5.1. Implementation. For the hereby developed analytical solutions, their numerical approximation follows the same two steps. Firstly truncate the series for -N ≤ n ≤ N . Secondly approximate the N th first eigenvalues λ n and eigenfunctions Ψ n (we omitted here the indexes D and N relatively to the Dirichlet or Neumann boundary condition). Once obtained these approximations, the coefficients α n in (15) (only depending on the Ψ n ) and the functions c n (z) (that only depend on the λ n and on the boundary data on ∂Ω × I) can in turn be approximated.

The numerical approximation of the λ n and of the Ψ n has been presented in [START_REF] Fehrenbach | The Generalized Graetz problem in finite domains[END_REF] for the Dirichlet case. We present the adaptation of the method to the Neumann case.

d = 2 l = 4 l = 4 2 R = 1 3 l d = 2
Figure 1. Geometrical configurations. On the left, periodic test case: the parallel circular ducts have diameter d = 2 and are embedded in square cells of size l = 4 and of same centre. On the Right, counter-current test case: two circular ducts of same diameter d = 2 are embedded in a solid matrix, also circular, with diameter 2R = 13. The two ducts are symmetrically located on each side of the matrix center, l denotes their distance that will vary. Definition 5.1 (weak formulation). The problem to be solved is, find

(Ψ, λ) ∈ D(A N ) × R * so that A N Ψ = λΨ. It is equivalent with: find (u, s) ∈ H 1 (Ω) × H 1 (Ω) and for λ ∈ R * so that for all (ũ, s) ∈ H 1 (Ω) × H 1 (Ω) we have, (36) a 1 (u, s), (ũ, s) = λ a 2 (u, s), (ũ, s) ,
where the bilinear products a 1 and a 2 are defined by, a

Ω v uũ + k∇u • ∇s + k∇ũ • ∇s dx, a 2 (u, s), (ũ, s) = Ω kuũ + k∇s • ∇s dx. 1 (u, s), (ũ, s) =
Then the eigenfunction Ψ is given by Ψ = (u, k∇s).

Here the unkmown Ψ has been replaced by (u, s) ∈ H 1 (Ω) × H 1 (Ω) so that Ψ = (u, k∇s). This is possible because Ψ ∈ R N (see corollary 2.2): such a change of variable avoids any problem eventually caused by the kernel of A N with the numerical methods.

Proof. Let Ψ ∈ D(A N ) satisfy A N Ψ = λΨ for some λ ∈ R ⋆ . Then Ψ ∈ R N and Ψ = (u, k∇s) for some (u, s) ∈ H 1 (Ω) × H 1 (Ω) (see corollary 2.2). We have for all (ũ, s) ∈ H 1 (Ω) × H 1 (Ω) that A(u, k∇s)|(ũ, k∇s) H = λ ((u, k∇s)|(ũ, k∇s)) H .
Developping the H-scalar product and using the Green formula exactly gives (36).

Conversely consider a solution (u, s) to (36) and form Ψ = (u, k∇s). One has to show that Ψ ∈ D(A N ) and that A N Ψ = λΨ. Writing (36) for s = 0 and for a smooth test function ũ with compact support in Ω gives,

-div(k∇s), ũ = Ω k∇s • ∇ũdx = Ω (λk -v)uũdx,
where •, • denotes the distribution product. This means that k∇s ∈ H div (Ω) and that k -1 (vudiv(k∇s)) = λu. Meanwhile writing (36) for ũ = 0 gives that k∇u = k∇s. As a result we have Ψ ∈ D(A) and AΨ = λΨ. It only remains to prove that k∇s • n = 0 on ∂Ω to ensure that Ψ ∈ D(A N ). This can be seen simply by rewriting (36) for s = 0 and ũ ∈ H 1 (Ω),

Ω k∇s • ∇ũdx = Ω (λk -v)uũdx,
and performing an integration by part it follows that, ∂Ω k∇s • nũds = Ω (λku -vu + div(k∇s)) ũdx = 0, ensuring that k∇s • n = 0 on ∂Ω.

The weak formulation (36) is used in practice for the approximation of the eigenvalues/eigenfunctions (λ n , Ψ n ) by considering a P1 -Lagrange finite element space P 1 (M) over a mesh M of the domain Ω. Considering a classical basis of P 1 (M), problem (36) leads to the resolution of the eigenvalue problem: find λ ∈ R and X ∈ P 1 (M) so that,

S h X = λM h X,
where S h and M h respectively are the stiffness and mass matrices associated with the products a 1 and a 2 written on the considered bases of P 1 (M). The matrix M h is symmetric positive definite and the matrix S h as well excepted when Q = 0, in this case S h is semi definite positive with a one dimensional kernel. The assembling for these two matrices is done using the finite element library FreeFem++ [6], being given a mesh of Ω also built using FreeFem++. In practice this assembling only involves to build sub-block-matrices that simply are classical mass and stiffness matrices, ie matrices for the L 2 products (u, v) → Ω uvdx and (u, v) → Ω k∇u • ∇vdx, which are built by FreeFem++. The resolution of the spectral problem S h X = λM h X uses the arpack++ library 1 .

Eventually, when considering a semi-infinite problem, we also have to determine the coefficients (β n ) n>0 in (32) or (34). This question is the topic of [START_REF] Bouyssier | Analysis of parallel convective exchangers with versatile inlet/outlet conditions[END_REF] for general entry conditions. For Dirichlet entry conditions it has primary been addressed in [START_REF] Fehrenbach | The Generalized Graetz problem in finite domains[END_REF]: we will adopt this strategy here. The sequence (β n ) n>0 is approximated by a vector b ∈ R N solution of a (symmetric positive definite) K h b = E h . The right hand side E h depends on the entry condition. The matrix K h = [k ij ] is computed using the eigen functions T n : simply k i j = Ω T i T j dx.

Axisymmetric configuration.

We first consider an axisymmetric configuration: the domain Ω is the circle with radius R = 2 and centre 0, the fluid part is is the circle of radius 1 and same centre 0. ) with respect to the mesh size h is here depicted on the left (resp. right) for three values of the Péclet number (P e =0.1, 1 and 10). The relative error is represented as a function of the mesh size h using a (decimal) Log/Log scale. Each plot displays the same linear behaviour with slope 2.

5.2.1.

Finite element solver evaluation. For this geometry, both the eigenvalues and the eigenfunctions have an analytical definition following a technique presented in [START_REF] Pierre | Generalised graetz problem: analytical solutions for concentric or parallel configurations[END_REF]. These analytical solutions can be computed rapidly with an arbitrary small accuracy using a Maple code. This code results will be considered as reference solutions. On the other hand we presented in the previous section our strategy to approximates the λ n and the Ψ n on an arbitrary domain Ω using a mesh of Ω and a finite element solver.

The purpose of this sub section is to evaluate the accuracy of this finite element solver by comparing solutions obtained with the finite element solver with the reference ones. We set the boundary condition to the Neumann case. A series of 16 meshes has been considered with a mesh size h varying between 0.157 and 0.0157. We first analyse the convergence of the computed eigenvalues: results are depicted on figure 2. The relative error of the first downstream and upstream eigenvalues with respect to their exact values λ 1 and λ -1 has been computed on each of these meshes. We considered three values of the Péclet number: P e = 0.1 (dominant diffusion), P e =1 and P e =10 (dominant convection). In all cases, the error goes to zero with an order two convergence with h.

- With the same setting we analysed the convergence of the corresponding coefficients α ±1 that involve the boundary integral of T ±1 . The results are depicted on figure 3. The convergence for the α n is of order 2 with the mesh size.

5.2.2.

First axisymmetric test case. As a first test case we consider the Neumann problem ( 14) on the infinite domain Ω × R. We impose the heat flux f (z) = 1 for z < 0 and f (z) = 0 for z > 0 on the boundary ∂Ω × R, whose solution is given in corollary 3.7. The Péclet number is set to P e = 10. We denote by T the solution given by equation (30) and we denote T N its approximation considering the N th -first computed eigenmodes. We are interested with the computation of the fluid/solid heat flux ϕ(N ), (37)

ϕ(N ) = +∞ 0 ∂O -k∇T N • ndl dz.
with O the fluid domain. So defined the fluid/solid heat flux for the exact solution T is given by ϕ = ϕ(+∞). This limit value ϕ has been evaluated using an extrapolation procedure. Using this evaluation of ϕ we computed the relative error e ϕ (N ) on the fluid/solid heat flux computation,

e ϕ (N ) = ϕ(N ) -ϕ ϕ ,
its behaviour is depicted on figure 4 (on the left). It can be seen on this graph that the error goes to zero with N and that the odd and even values of the error follow two different curves. This two curves however display the same asymptotic behaviour: a geometric convergence towards zero of order 2 with N . In addition to the convergence asymptotic, it is quite important to notice that one can get very good approximations on the fluid/solid heat flux using very few modes. We always have a prediction with less than 10 % of accuracy and using only 3 modes this accuracy is less than 1 %. The exact solution T is given by equation (34). The function F (z) is any primitive of f (z) and is defined up to a constant. We set F 0 (z) the primitive of f (z) so that F 0 (0) = 0 and we write F (z) = F 0 (z) + β 0 Q/P with β 0 a constant. By doing so the constant on F (z) is considered as a supplementary unknown β 0 . Introducing the supplementary eigenmode T 0 = 1 and λ 0 = 0, We can rewrite the temperature T as:

T = F 0 (z)P/Q + n∈Z ⋆ α n c n (z)T n + n≥0 β n e λnz T n ,
We approximate T by T N using a truncation at order N of the series and by computing the N th -first eigenmodes using the maple code. The computation of the constants (β n ) 0≤n≤N is done following [START_REF] Fehrenbach | The Generalized Graetz problem in finite domains[END_REF] as briefly presented in section 5.1. We numerically observed that the entry condition T = 0 ensures the existence and unicity of the constants β 0 , . . . , β N (also depending on N ).

With this setting we computed the fluid/solid heat flux ϕ(N ) as previously defined in (37). The exact heat flux ϕ is again computed by extrapolation on the sequence ϕ(N ) allowing to compute the relative error e ϕ (N ) on the computed flux. The relative error is depicted on figure 4 on the right. To obtain the asymptotic regime we again had to distinguish between the odd and even values on the flux. The convergence speed is really fast and is more than algebraic. We observed an exponential convergence e ϕ (N ) = O(exp(-cN )) with 1/5 ≤ c ≤ 1/4, precisely c has been evaluated to c ≃ -0.22. An important remark is that we have an easy evaluation of the limit temperature as z goes to +∞. For a given N , we have:

T ∞ (N ) := lim z→+∞ T N (z) = P Q +∞ 0 f (z)dz + β 0 .
The exact temperature at +∞ T ∞ := lim z→+∞ T (z) is evaluated by extrapolation on the sequence T ∞ (N ) and we define the relative error e T∞ on the temperature at infinity as,

e T∞(N ) = T ∞ (N ) -T ∞ T ∞ .
Its behaviour is depicted on figure 5 on the left. We observed a extremely fast onvergence of T ∞ (N ) towards T ∞ of algebraic type e T∞(N ) = O(N -c ) and with order c ≃ 5.7. Again, beyond the convergence asymptotic, the important fact is that we already have a precision better than 0.1 % using only one Graetz mode ! It is quite interesting that the high precision and the fast convergence for the far field temperature and the fluid/solid heat flux computation are here obtained even considering a non-regular (discontinuous) boundary heat source term f (z). Eventually, we also plotted the temperature profiles on this configuration. The three profiles z → T (r, z) for r = R, r = 1 and r = 0 have been plotted on figure 5. They represent the temperature at the solid wall surace, at the fluid/solid interface and at the duct centre. One can check on these profiles that the entry condition T = 0 at z = 0 is well respected. The temperature secondly increases between the entry z = 0 and the strating point of the heated region z = a. On the right of the heated region, the temperature reaches fastly the limit temperature T ∞ as z → +∞. 5.3. Periodic test case. We consider in this test case a periodic geometry depicted on figure 1. It consists in a series of parallel circular ducts, in which a fluid is flowing, disposed inside a solid. The domain Ω is not bounded but periodic in the horizontal direction. It is composed of a series of squares of size l = 4. Each squre is made of a fluid domain: the circle of radius 1 and with centre the square centre, and of a solid matrix. The distance between two successive ducts then also is l = 4. We consider an infinite configuration Ω×R with an imposed boundary flux given by set to f (z) = 1 if z < 0 and f (z) = 0 for z > 0. Considering one elementary square, a Neumann boundary condition is considered on its top and bottom edges whereas periodic conditions on T and ∂ x T are imposed on its left and right edges. For symmetry reasons, this boundary condition is equivalent with an homogeneous Neumann one on the whole square cell boundary. The analytical solution for this problem is in corollary 3.7. The Péclet number for each duct is constant P e = 10. The tentth first Graetz modes are depicted on figure (6). One can note on this figure that the apparent structural complexity of the Graetz modes increases with N . Our focus concerns the influence of modes truncation on exchanges estimation in order to evaluate exchanger's performance. We computed on this configuration the fluid/solid heat flux ϕ(N ) as previously defined in (37). Figure 6 illustrates the contribution of the first ten contributing modes to the exchange flux. For the chosen convection dominated examined situation, P e = 10, one can see that few Graetz modes are enough to obtain accurate estimates of the fluid/solid flux. This is more precisely stated in table 1 where the fluxes values ϕ(N ) relatively to figure 6 have been reported together with the asoociated relative error e ϕ (N ). Though these results (on the contrary of the previous computation in section 5.2) might be blurred by some numerical error induced by discretisation, they clearly indicate that with very few Graetz modes one obtains accurate estimations on the fluid/solid flux ϕ with less than 1% of relative error. This is a very interesting observation, that convection dominated configuration provide an excellent performance for the proposed mode decomposition, so that it is quite easy to get fast and accurate estimate of the exchange performance using the proposed formulation, and this even for a more complex geometrical configuration.

5.4.

Counter current case. We finally consider the counter current configuration, for which the total debit Q = 0. It is depicted on the right side of figure 1 and consists of two parallel circular ducts, where a fluid is flowing in opposite directions. The two ducts are encapsulated inside a cylindrical solid with diameter 2R = 13. The two ducts are symmetrical with the solid centre and their centres distance l is variable. Various values for the Pechlet number will be considered.

We adopt an infinite configuration where the outer solid wall is heated for z < 0 and verifies a zero flux condition for z > 0 as described in corollary 3.7. The temperature in this case is given by (31): in the z > 0 region, the functions F and G are equal to zero. We concentrate on the heat flux on the two internal duct boundaries for z > 0, as defined in (37). These fluxes only depend on the constant C 1 in (31), this constant represent the heat flux ∂ z T as z → +∞ and is set to C 1 = 0 here. The roles of the two ducts is absolutely non symmetric. The fluid in the left tube flows towards the z > 0 region. This fluid is heated in the z < 0 region by heat diffusion process in the solid. It then brings heat by convection to the z > 0 region: then this tube can be considered as the input duct. On the contrary the fluid in the right tube flows towards the z < 0 region and it evacuates heat by convection from the z > 0 region: it can be considered as the output duct. or P e = 10 and for a distance l = 1.5d = 3 between the two ducts centres above. Below, the Pechlet number is set to P e = 10 and the distance l between two ducts centres varies from l = 1.5d = 3 and l = 3d = 6 to l = 4.5d = 9. It is interesting to mention that the convergence rate of the flux according to the number of considered eigenmodes is sensitive to the chosen geometrical parameters as well as to the Péclet number. Qualitatively the closer the tubes, the faster mode truncation converges to the exchange flux. On the other hand, increasing the Péclet number provides a slower mode convergence as observed on the upper part of figure 7. Nevertheless, an estimate of the exchange flux accurate within a few percent is obtained in every configurations when aggregating the contribution of less than ten modes. Another interesting observation is provided on figure 8 where one can observe the spatial structure of the most contributing modes to the exchange flux. We focus for this figure on the case P e = 10, l = 1.5d and on the left tube. The spectral convergence of the fluid/solid flux on this left tube is displayed together with the visualisation of the controbuting Graetz modes. As for figure 6, it can be observed that the spatial structure of the modes increases in complexity as their contribution to the exchange flux decreases. For example, the first mode is mostly of zeroth azimuthal order, the second and third modes are of first azimuthal order, the fourth to six modes are principally of azimuthal order three, etc... Nevertheless this observation is not a golden rule since the seventh mode has a azimuthal order one, with an horizontal symmetry as opposed to the second mode which also has a first azimuthal order but with a vertical symmetry. This observation indicates that the chosen configurations favours some symmetries.

Finally the modal convergence of the flux ϕ(N ) displayed on figure 8 is detailed in table 2. It can be observed thet considering a few Graetz modes provides an accurate estimation of the fluid/solid flux: 5 to 10 Graetz modes are sufficient to get a 1% accurate evaluation. This confirm the observation made for the periodic test case, even consideringa convection dominated configuration and a complex geometry, the Graetz mode decomposition remains efficient to accurately capture the physically important features of the heat transfer.

Conclusion

. This work has permitted to extent the two-dimensional mapping of longitudinally invariant convection/diffusion problems to very general configurations with either prescribed field or fluxes at the outer boundary. In the case of prescribed fluxes, it is necessary to distinguish the case of zero total convective flux (typically encountered in counter-current configurations) from the case of non-zero convective flux. In both cases, we found general analytical expression for the longitudinal variation of the solution, which depends on the applied boundary condition. Those considerations apply to convective exchangers and have been illustrated in some non-trivial configuration to illustrate the versatility and the numerical efficiency of the method for studying complex configurations. This analysis opens new perspective for a systematic and accurate study of convective exchangers and towards their optimisation. E-mail address: fplourab@imft.fr

Definition 2 . 1 .

 21 We respectively introduce two restrictions A D and A N of the operator A relatively to a homogeneous Dirichlet (D) or homogeneous Neumann (N ) boundary condition with domains D(A D ) and D(A N ):

.

  eigenfunctions for the Dirichlet and Neumann cases respectively. We have the relation q D,Two important consequences are (skipping the indexes D and N ):

  The function u N satisfies the elliptic equation v au o + b -div(k∇u N ) = ak together with the boundary condition k∇u N • n = 1 on ∂Ω. The justification of the well posedness of a, b and u N is detailed in the following proof.

3. 3 .

 3 The Neumann problem for Q = 0. We simply denote here λ n = λ N n , Ψ n = Ψ N n and α n = α N n . The functions c n (z) for n ∈ Z ⋆ are alternatively defined as, (

3. 4 .

 4 The Neumann problem for Q = 0. Let us adapt the previous results to the particular case Q = 0: notations are unchanged as for the Neumann problem with Q = 0.

Figure 2 .

 2 Figure2. Eigenvalue convergence. The convergence of the first downstream (resp. upstream) computed eigenvalue towards its exact value λ 1 (resp. λ -1 ) with respect to the mesh size h is here depicted on the left (resp. right) for three values of the Péclet number (P e =0.1, 1 and 10). The relative error is represented as a function of the mesh size h using a (decimal) Log/Log scale. Each plot displays the same linear behaviour with slope 2.

Figure 3 .

 3 Figure3. Convergence of the coefficients α n . The convergence of the first downstream (resp. upstream) computed coefficients α n towards its exact value α 1 (resp. α -1 ) with respect to the mesh size is here depicted on the left (resp. right) for three values of the Péclet number (P e =0.1, 1 and 10). The relative error is represented as a function of the mesh size h using a (decimal) Log/Log scale. Each plot displays the same linear behaviour with slope 2.

Figure 4 .

 4 Figure 4. Relative error e ϕ (N ) on the computed fluid/solid heat flux ϕ(N ) according to the number of considered eigenmodes N . On the left: infinite domain configuration. The ood and even values of the error have been plotted separately. By doing so we observe the same algebraic convergence of order 2 with N of the error: e ϕ (N ) = O(N -2 ). On the right: semi-infinite configuration. Again we distinguished between the odd and even values of the error to observe an exponential convergence e ϕ (N ) = O(exp(cN )) where c has been evaluated to c ≃ -0.22.

5. 2 . 3 .

 23 Second axisymmetric test case. For the second test case we now consider the Neumann problem (14) on a semi-infinite domain. We impose the heat flux f (z) on the boundary ∂Ω×(0, +∞) with, f (z) = 1 for z ∈ [a, b] and f (z) = 0 otherwise, which means that we consider here a tube insulated outside the heated region [a, b] and with a homogeneous heating f (z) = 1 inside [a, b]. We set the heated region defining a = 2R = 4 and b = 3R = 6. The homogeneous Dirichlet condition T = 0 is considered at the entry z = 0 modelling a cold fluid injection.

Figure 5 .

 5 Figure 5. Semi infinite test case. Left: relaive error e T∞(N ) on the predicted temperature at z = +∞. We observe a geometric convergence towards zero with a high order evaluated to 5.7. Right: representation of the temperature profiles z → T (r, z) for three fixed values of r: r = R (boundary), r = 1 (fluid/solid interface) and r = 0 (duct centre). The heated region (interval [a, b] where f (z) = 0) is between z = 2R = 4 and z = 3R = 6. The temperature T ∞ at z = +∞ is also plotted.

Figure 6 .

 6 Figure 6. Periodic configuration. Fluid-solid flux computed with the ten first eigenmodes and their visualisations

N

Figure 7 .

 7 Figure 7. Counter current configuration. Fluid/solid heat exchange on the left and right duct boundaries. These fluxes are computed with an increasing number N of eigen modes, the dependence of the computed flux with N is here depicted for various values of the Péclet number (above), and for a varying distance l between the two pipes (below). On figure7is displayed the evaluation of the fluid/solid exchange on both left and right pipe boundaries. These fluxes are computed for three values of the Pechley number, P e = 0.1, P e = 1 or P e = 10 and for a distance l = 1.5d = 3 between the two ducts centres above. Below, the Pechlet number is set to P e = 10 and the distance l between two ducts centres varies from l = 1.5d = 3 and l = 3d = 6 to l = 4.5d = 9. It is interesting to mention that the convergence rate of the flux according to the number of considered eigenmodes is sensitive to the chosen geometrical parameters as well as to the Péclet number. Qualitatively the closer the tubes, the faster mode truncation converges to the exchange flux. On the other hand, increasing the Péclet number provides a slower mode convergence as observed on the upper part of figure7. Nevertheless, an estimate of the exchange flux accurate within a few percent is obtained in every configurations when aggregating the contribution of less than ten modes.

Figure 8 .

 8 Figure 8. Counter current configuration. Fluid/solid flux on the left tube boundary for z > 0 computed with the ten first Graetz modes and their visualisations (for P e=10 and l = 1.5d).

N

[

  25] B. Weigand and F. Wrona. The extended Graetz problem with piecewise constant wall heat flux for laminar and turbulent flows inside concentric. Heat and Mass Transfer, 39(4):1432-1181, 2003. (Julien Bouyssier) Université de Toulouse, INPT, UPS, IMFT (Institut de Mécanique des Fluides de Toulouse), Allés Camille Soula, F-31400 Toulouse, France, and CNRS, IMFT, F-31400 Toulouse, France. E-mail address: jbouyssi@imft.fr (Charles Pierre) Laboratoire de Mathématiques et Applications Université de Pau et du Pays de l'Adour av. de l'Université BP 1155 64013 PAU Cedex -FRANCE E-mail address: charles.pierre@univ-pau.fr (Franck. Plouraboué) Université de Toulouse, INPT, UPS, IMFT (Institut de Mécanique des Fluides de Toulouse), Allés Camille Soula, F-31400 Toulouse, France, and CNRS, IMFT, F-31400 Toulouse, France.

  n∈Z ⋆ |λ n α n c n (z)| 2 and n∈Z ⋆ |α n c ′ n (z)| 2 are uniformly converging. We already have from lemma 3.1 that n∈Z ⋆ |α n | 2 < +∞. The uniform convergence then follows from the upper bounds

  Remark 4.2. As stated in remark 4.1, the role of ψ 2 (z) = n>0 β n e λnz Ψ n is to satisfy a prescribed entry condition. The function ψ 2 is a solution to the homogeneous Neumann problem (4) (14) (as long as n∈Z ⋆ |β n | 2 < +∞).

n ,

Table 1 .

 1 Periodic

		1	2	3	5	8	10
	ϕ(N )	-7.10	-7.28	-7.11	-6.98	-7.02	-7.04
	e ϕ (N )	0.011	0.036	0.012	0.007	0.002	0.001

configuration. Fluid/solid flux ϕ(N ) with respect to N and the associated relative error e ϕ (N ), the reported figures are relative to figure 6.

Table 2 .

 2 Counter current configuration. Fluid/solid flux ϕ(N ) for the left tube with respect to N and the associated relative error e ϕ (N ), the reported figures are relative to figure 8.

		1	2	3	5	8	10
	ϕ(N )	-25.84	-39.32	-38.03	-36.56	-37.02	-37.10
	e ϕ (N )	0.30	0.058	0.024	0.016	0.003	0.001

J. BOUYSSIER, C. PIERRE, AND F. PLOURABOU É

Arpack software: http://www.caam.rice.edu/software/ARPACK/

Acknowledgments. The authors thank F. De. Gournay and J. Ferhenbach for very fruitful discussions. They led us to a deeper understanding of the physics and of the mathematical properties of the problem, as well as to an efficient implementation to obtain the numerical results.. The authors are indebted to F. De Gournay for insinghful mathematical remarks concerning the extension of the compactness properties.