Analytical properties of Graetz modes in parallel and concentric configurations

TL;DR

This paper develops an analytical framework for solving the generalized Graetz problem in layered cylindrical or parallel domains, enabling explicit, mesh-less computation of heat transfer modes for heat exchanger applications.

Contribution

It introduces a new analytical approach to compute Graetz modes explicitly in layered configurations, simplifying the spectral problem and facilitating mesh-less solutions.

Findings

Explicit formulas for eigenvalues via zeros of an analytical series

Recursive computation of eigenfunctions (Graetz modes)

Mesh-less, software-compatible solution method

Abstract

The generalized Graetz problem refers to stationary convection-diffusion in uni-directional flows. In this contribution we demonstrate the ana-lyticity of generalized Graetz solutions associated with layered domains: either cylindrical (possibly concentric) or parallel. Such configurations are considered as prototypes for heat exchanger devices and appear in numerous applications involving heat or mass transfer. The established framework of Graetz modes allows to recast the 3D resolution of the heat transfer into a 2D or even 1D spectral problem. The associated eigefunctions (called Graetz modes) are obtained with the help of a sequence of closure functions that are recursively computed. The spectrum is given by the zeros of an explicit analytical serie, the truncation of which allows to approximate the eigenvalues by solving a polynomial equation. Graetz mode computation is henceforth made explicit and can be performed using standard softwares of formal calculus. It permits a direct and mesh-less computation of the resulting solutions for a broad range of configurations. Some solutions are illustrated to showcase the interest of mesh-less analytical derivation of the Graetz solutions, useful to validate other numerical approaches.

Equations

$${\rm div}_{(r,\varphi,z)}\big{(}\mathbf{v}T-k\nabla_{(r,\varphi,z)}T\big{)}=0.$$

$$\mathrm{div}(k\nabla T)+k{\partial }_z^{2}T=v{\partial }_{z}T\text{in}\Omega \times (a,b),$$

$$T=g(z)\text{or}k\nabla T=g(z)\text{on}\partial \Omega \times (a,b).$$

$$\partial_{z}\Psi=A\Psi\quad\text{with}\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \Psi=(T,\mathbf{p})\leavevmode\nobreak\ ,\quad A=\left(\begin{array}[]{cc}vk^{-1}&-k^{-1}{\rm div}(\cdot)\\ k\nabla\cdot&0\end{array}\right),$$

$$-\infty \leftarrow \le {\lambda }_i\le \cdots \le {\lambda }_1<0<{\lambda }_{-1}\le \cdots \le {\lambda }_{-i}\to +\infty .$$

$$\mathrm{div}(k(r)\nabla \Theta )+{\lambda }^{2}k(r)\Theta =\lambda v(r)\Theta \text{on}\text{leavevmode}\text{leavevmode}\text{leavevmode}\Omega ,$$

$$T(r,\phi ,z)=\sum _{i\in {\mathbb{Z}}^{\star }}{c}_i(z){\Theta }_i(r,\phi )e^{{\lambda }_{i}z}.$$

$$T(r,\phi ,z)=\sum _{i\in {\mathbb{Z}}^{+}}{c}_i{\Theta }_i(r,\phi )e^{{\lambda }_{i}z}.$$

$$T^i=\sum _{i\in {\mathbb{Z}}^{+}}{c}_i{\Theta }_i.$$

$$T^i=\sum _{i\in {\mathbb{Z}}^{+}}{c}_i{\Theta }_i+\sum _{i\in {\mathbb{Z}}^{-}}{c}_i{\Theta }_i{e}^{-{\lambda }_{i}L},T^o=\sum _{i\in {\mathbb{Z}}^{+}}{c}_i{\Theta }_i{e}^{{\lambda }_{i}L}+\sum _{i\in {\mathbb{Z}}^{-}}{c}_i{\Theta }_i,$$

$${\Delta }_{n}f=\frac{1}{r}\text{diff}r(r\text{diff}rf)-\frac{n^2}{r^2}f.$$

$${\lambda }^2{k}_{j}T+k_j{\Delta }_{n}T=\lambda v_{j}T,\text{on}(r_{j-1},r_j),j=1\dots m,$$

$$T(r_j^{-})=T(r_j^{+})\text{leavevmode},k_j\text{diff}rT(r_j^{-})=k_{j+1}\text{diff}rT(r_j^{+}),j=1\dots m-1.$$

$$T_{n,\lambda }(r)\sim r^n\text{as}r\to 0^{+}.$$

$$\Lambda =\bigcup _{n\in \mathbb{N}}{\Lambda }_n,{\Lambda }_n=\left\{\lambda \in \mathbb{C},T_{n,\lambda }(R)=0\right\},$$

$$\Lambda =\bigcup _{n\in \mathbb{N}}{\Lambda }_n,{\Lambda }_n=\left\{\lambda \in \mathbb{C},\text{diff}r{T}_{n,\lambda }(R)=0\right\}$$

$$\Delta \Theta =\frac{1}{k_j}\left(\lambda v_j(r)\Theta -{\lambda }^2\Theta \right).$$

$$\Theta =\sum _{n\in \mathbb{Z}}{\theta }_n(r)e^{-in\phi }.$$

$$\Delta \Theta =\sum _{n\in \mathbb{Z}}{\Delta }_n\left({\theta }_n(r)\right)e^{-in\phi }.$$

$$\Theta =\sum _{n\ge 0}{\beta }_{|n|}{T}_{n,\lambda }(r)\cos (n\phi )+\sum _{n>0}{\gamma }_{|n|}{T}_{n,\lambda }(r)\sin (n\phi ).$$

$$T_{n,\lambda }(r)=\sum _{p\in \mathbb{N}}{t}_{n,p}(r){\lambda }^p,$$

$$k_j{\Delta }_n{t}_{n,p}+k_j{t}_{n,p-2}=v(r)t_{n,p-1}.$$

$$\Lambda =\bigcup _{n\in \mathbb{N}}{\Lambda }_n,{\Lambda }_n=\left\{\lambda \in \mathbb{C},\sum _{p\in \mathbb{N}}{c}_{n,p}{\lambda }^p=0\right\},$$

$$k_j{\Delta }_n{t}_p+k_j{t}_{p-2}=v(r)t_{p-1}\text{on}\text{leavevmode}\text{leavevmode}(r_{j-1},r_j),$$

$$t_p(r_j^{+})=t_p(r_j^{-}),k_j\text{diff}r{t}_p(r_j^{+})=k_{j+1}\text{diff}r{t}_p(r_j^{-}),$$

$$\underset{r\to 0}{\lim }\frac{t_p(r)}{r^n}=0.$$

$$F_j[f](r):=r^n{\int }_{r_{j-1}}^r\frac{1}{x^{2n+1}}{\int }_{r_{j-1}}^x{y}^{n+1}f(y)\mathrm{d}y\mathrm{d}x,$$

$$f_{p-1}:=\frac{v}{k}{t}_{p-1}-t_{p-2}.$$

$$t_p(r)={\alpha }_j{\psi }_1(r)+{\beta }_j{\psi }_2(r)+F_j[f_{p-1}](r).$$

$$t_p(r)=F_1[f_{p-1}](r)\mathrm{on}[0,r_1].$$

Full text

Keywords: Heat and mass transfer, convection-diffusion, reduced problem, separation of variables, analytical solutions.

Introduction

Parallel convective heat exchangers are relevant in many applicative contexts such as heating/cooling systems [23], haemodialysis [9], as well as convective heat exchangers [14]. A number of works devoted to parallel convective heat exchangers in simple two dimensional configurations [15, 16, 11, 12, 27, 24, 28, 25] can be found to cite only a few, whilst many others can be found in a recent review [6].

As quoted in [6] conjugate heat transfer are mixed parabolic/hyperbolic problems which makes them numerically challenging. In many applications the ratio between the solid and fluid thermal conductances is high (larger than one thousand in many cases). The convection is dominating, so that the ratio of convection to diffusion effects provided by the so-called Péclet number is very high (e.g. larger than $10^{5}$ in [21, 20]). When dealing with such highly hyperbolic situations, numerical convergence might be an issue. The increase in computer power has permitted and popularized the use of direct numerical simulations to predict heat exchanger performances [21, 20, 26, 7, 13]. The derivation of analytical mesh-less reference solutions allows to evaluate the accuracy and the quality of the discrete solutions, as done in [18, 8, 2, 17, 5] in a finite-element framework. In most cases, it is interesting to validate the numerical solution in simple configurations as well as being able to test the solution quality for extreme values of the parameters, when rapid variations of the temperature might occur in localized regions. However, few analytic solutions are known, apart from very simplified cases. Namely, such analytic solutions can be obtained for axi-symmetric configurations, when the longitudinal diffusion has been neglected whilst assuming a parabolic velocity profile, as originally studied by Graetz [10]. In this very special case, the Graetz problem maps to a Sturm-Liouville ODE class, and the resulting analytic solutions can be formulated from hypergeometric functions, see [4] or for example [25].

In this contribution we introduce analytical generalized Graetz modes: including longitudinal diffusion, for any regular velocity profile, and for general boundary conditions. The derivation of the generalized Graetz modes follows an iterative process that can be performed using a standard formal calculus software. Then, section 1.1 sets notations (mainly for the cylindrical case) and provides the physical context as well as the constitutive equations under study. Section 1.2 gives the necessary mathematical background for the subject, with an emphasis on most recent results useful for the presented analysis. Section 2.1 shows that discrete mode decomposition also holds for non-axi-symmetric configuration. Section 2.3 gives the central result of this contribution regarding the analyticity of the generalized Graetz modes. Finally section 3 illustrate specific applications obtained with the method with explicit analytical computations.

1 Setting the problem

1.1 Physical problem

We study stationary convection-diffusion in a circular duct made of several concentric layers (fluid or solid). The domain is set to $\Omega\times(a,b)$ with $(a,b)\subset\mathbb{R}$ an interval and $\Omega$ the disk with centre the origin and radius $R$. The longitudinal coordinate is denoted by $z$ and cylindrical coordinates $(r,\varphi)$ are used in the transverse plane. Then $\Omega$ is split into $m$ different compartments $\Omega_{j}$, $j=1\dots m$, either fluid or solid and centered on the origin: $\Omega_{1}$ is the disk of radius $r_{1}$ whereas $\Omega_{j}$, $j\geq 2$, is the annulus with inner and outer radius $r_{j-1}$ and $r_{j}$ for a given sequence $0<r_{1}<\dots<r_{m}=R$. Two such configurations are depicted on Figure 1.

The physical framework is set as follows:

Velocity: $\mathbf{v}(r,\varphi,z)=v(r)\leavevmode\nobreak\ \mathbf{e}_{z}$ with $\mathbf{e}_{z}$ the unit vector along the $z$ direction. We denote $\mathbf{v}_{j}=\mathbf{v}_{|\Omega_{j}}=v_{j}(r)\leavevmode\nobreak\ \mathbf{e}_{z}$ the restriction of the velocity to compartment $\Omega_{j}$. In case this compartment is solid we have $v_{j}(r)=0$.

We make the mathematical assumption that each $v_{j}(r)$ is analytical, though $v(r)$ is allowed to be discontinuous at each interface. 2. 2.

Conductivity: $k(r,\varphi,z)=k(r)$ and moreover $k_{|\Omega_{j}}=k_{j}>0$ is a constant.

The general equation for stationary heat convection-diffusion reads

[TABLE]

With the assumptions we have made, it simplifies to

[TABLE]

where we denoted by ${\rm div}={\rm div}_{(r,\varphi)}$ and $\nabla=\nabla_{(r,\varphi)}$ the gradient and divergence operators restricted to the transverse plane. The following boundary conditions, either of Dirichlet or Neumann type, are considered

[TABLE]

1.2 Mathematical background

Problem reformulation

Adding a supplementary vector unknown $\mathbf{p}:\leavevmode\nobreak\ \Omega\to\mathbb{R}^{2}$, problem (1)-(2) has been reformulated in [18, 8, 2, 17] into a system of rwo coupled PDEs of first order:

[TABLE]

on the space $\mathcal{H}={\rm L}^{2}(\Omega)\times[{\rm L}^{2}(\Omega)]^{2}$ and involving the differential operator $A:\leavevmode\nobreak\ D(A)\subset\mathcal{H}\rightarrow\mathcal{H}$. The definition of the domain $D(A)$ of the operator $A$ depends on the chosen Dirichlet or Neumann boundary condition.

For simplicity we briefly recall the properties of operator $A$ in the Dirichlet case, as presented in [18, 8]. These properties have been extended to the Neumann case in [17] and to the Robin case in [3]. For a Dirichlet boundary condition, we set $D(A)={\rm H}^{1}_{0}(\Omega)\times{\rm H}_{\operatorname{div}}(\Omega)$. Then $A$ is self-adjoint with compact resolvent. Apart from the kernel space $K:={\rm ker}A=\{(0,\mathbf{p}),\leavevmode\nobreak\ \mathbf{p}\in{\rm H}_{\operatorname{div}}(\Omega),\leavevmode\nobreak\ \operatorname{div}\mathbf{p=0}\}$ the spectrum of $A$ is composed of a set $\Lambda$ of eigenvalues of finite multiplicity. It has been shown in [18] that $\Lambda$ decomposes into a double sequence of eigenvalues $\lambda_{i}$,

[TABLE]

We call upstream eigenvalues the positive eigenvalues $\{\lambda_{i},\leavevmode\nobreak\ i<0\}$ and downstream eigenvalues the negative ones $\{\lambda_{i},\leavevmode\nobreak\ i>0\}$. The associated eigenfunctions $(\Psi_{i})_{i\in\mathbb{Z}^{\star}}$ form an orthogonal (Hilbert) basis of $K^{\perp}$.

Eigenmodes

Let us write $\Psi_{i}=(\Theta_{i},\mathbf{p}_{i})$ the eigenfunctions. Their vector component satisfies $\mathbf{p}_{i}=k\nabla\Theta_{i}/\lambda_{i}$. It is important to understand that $\Theta_{i}:\leavevmode\nobreak\ \Omega\mapsto\mathbb{R}$ only is the scalar component of the associated eigenfunction $\Psi_{i}$. As a result the $(\Theta_{i})_{i\in\mathbb{Z}^{\star}}$ are not eigenfunctions themselves, they are neither orthogonal nor form a basis of ${\rm L}^{2}(\Omega)$. To clarify this distinction we refer to $\Theta_{i}$ as an eigenmode associated with $\lambda_{i}$.

Eigenmodes can be directly defined through a generalized eigenvalue problem. A function $\Theta:\Omega\rightarrow\mathbb{R}$ is an eigenmode if $\Theta\in{\rm H}^{1}(\Omega)$, $k\nabla\Theta\in{\rm H}_{\operatorname{div}}(\Omega)$ and their exists a scalar $\lambda$ so that

[TABLE]

with $\Theta=0$ or $\nabla\Theta\cdot\mathbf{n}=0$ on $\partial\Omega$ depending on the considered Dirichlet or Neumann boundary condition. In that situation, $\lambda$ is an eigenvalue of $A$ associated with the eigenfunction $\Psi_{i}=(\Theta,k\nabla\Theta/\lambda)$. As a consequence, the eigenmodes always are real functions since the operator $A$ is symmetric.

The upstream and downstream eigenmodes have the following important property (proved in [8]):

• •

The upstream eigenmodes $\{\Theta_{i},\leavevmode\nobreak\ i<0\}$ form a (Hilbert) basis of ${\rm L}^{2}(\Omega)$.

• •

The downstream eigenmodes $\{\Theta_{i},\leavevmode\nobreak\ i>0\}$ also form a basis of ${\rm L}^{2}(\Omega)$.

Problem resolution

The problem (1)-(2) can be solved by separation of variables. General solutions for non-homogeneous boundary conditions of Dirichlet, Neumann or Robin type have been derived in [18, 8, 17, 2, 3]. Such solutions are detailed in section 3. We simply recall their formulation for a homogeneous Dirichlet boundary condition:

[TABLE]

The functions $c_{i}(z)$ are determined with the help of the eigenmodes, of the boundary condition $g(z)$ and of the inlet/outlet conditions. As an illustration, we precise that derivation in two cases. In the case of a homogeneous boundary condition $g(z)$=0 in (2), then $c_{i}(z)=c_{i}\in\mathbb{R}$ are constant scalars. On a semi-infinite domain $\Omega\times(0,+\infty)$, the upstream coefficients are zero, $c_{i}=0$ for $i<0$, and

[TABLE]

The coefficients $c_{i}$ for $i>0$ are given by the inlet condition $T^{i}=T_{|z=0}$

[TABLE]

If the domain is finite, equal to $\Omega\times(0,L)$, then the upstream coefficients are no longer equal to zero, the upstream and downstream coefficients $c_{i}$ satisfy

[TABLE]

where $T^{o}=T_{|z=L}$ is the outlet condition.

2 Analyticity of the generalized Graetz modes

2.1 Series decomposition

To take advantage of the azimuthal symmetry of the physical problem we perform the Fourier decomposition of the eigenmodes. Their Fourier series expansion is composed by terms of the form $T(r)\cos(n\varphi)$ or $T(r)\sin(n\varphi)$. We prove here that we have a finite number of such terms and characterize $T(r)$. Let us introduce the operator $\Delta_{n}$

[TABLE]

Consider $\Theta$ an eigenmode associated with $\lambda\in\Lambda$ and assume that $\Theta(r,\varphi)=T(r)\cos(n\varphi)$ or $\Theta(r,\varphi)=T(r)\sin(n\varphi)$. Then $T$ is a solution of the following ODEs

[TABLE]

that are coupled with the transmission conditions

[TABLE]

Lemma 1**.**

For all $\lambda\in\mathbb{C}$ and all $n\in\mathbb{N}$ there exists a unique function $T_{{n},\lambda}(r):\leavevmode\nobreak\ (0,R)\rightarrow\mathbb{R}$ that satisfies (5)-(6) together with the normalisation condition

[TABLE]

An eigenmode $\Theta$ associated with the eigenvalue $\lambda$ decomposes as a finite sum of terms of the form $T_{{n},\lambda}(r)\cos(n\varphi)$ or $T_{{n},\lambda}(r)\sin(n\varphi)$.

The eigenvalue set $\Lambda$ decomposes in the Dirichlet case as

[TABLE]

and in the Neumann case as

[TABLE]

Finally, if $\lambda\in\Lambda_{n}$, then the associated eigenmodes are $T_{{n},\lambda}(r)\cos(n\varphi)$ and $T_{{n},\lambda}(r)\sin(n\varphi)$.

Proof of lemma 1.

The well poseness of the function $T_{{n},\lambda}$ definition is obtained by induction on the intervals $[r_{j-1},r_{j}]$. Assume that $T_{{n},\lambda}$ is given on $[r_{j-1},r_{j}]$ for some $j\geq 1$. On $[r_{j},r_{j+1}]$ the ODE (5) is regular and has a space of solution of dimension two, therefore $T_{{n},\lambda}$ is uniquely determined by the two initial conditions (6).

Now on $[0,r_{1}]$: the ODE (5) is singular at $r=0$. The Frobenius method (see e.g. [22]), with the assumption that $v(r)$ is analytical on $[0,r_{1}]$, states that the space of solutions is generated by two functions whose behavior near $r=0$ can be characterized:

• •

for $n>0$, one solution is $O(r^{n})$ at the origin and the second is $O(r^{-n})$,

• •

for $n=0$, one solution is $O(1)$ at the origin and the second is $O(\log(r))$,

Therefore condition (7) ensures existence and uniqueness for $T_{{n},\lambda}$.

Let $\Theta$ be an eigenmode for $\lambda\in\Lambda$. On each sub-domain $\Omega_{j}$, equation (4) can be rewritten as

[TABLE]

Using the assumption that $v_{j}(r)$ is analytical on $[r_{j-1},r_{j}]$, elliptic regularity properties imply that $\Theta\in C^{\infty}(\overline{\Omega_{j}})$. Moreover, since $\Theta\in{\rm H}^{1}(\Omega)$ and $k\nabla\Theta\in{\rm H}_{\operatorname{div}}(\Omega)$, it follows that $\Theta$ and $k\nabla\Theta\cdot\mathbf{n}$ are continuous on each interface between $\Omega_{j}$ and $\Omega_{j+1}$. We consider the Fourier series expansion for $\Theta$

[TABLE]

Since $\Theta\in C^{\infty}(\overline{\Omega_{j}})$ we can differentiate under the sum to obtain

[TABLE]

and so equation (4) ensures that each Fourier mode $\theta_{n}(r)$ satisfies the ODEs (5). It also satisfies the transmission conditions (6) because of the continuity of $\Theta$ and of $k\nabla\Theta\cdot\mathbf{n}$ at each interface. We already studied the behavior of the solution of (4) at the origin. Among the two possible behaviors characterized by the Frobenius method, $\Theta\in{\rm H}^{1}(\Omega)$ and $\nabla\Theta\in{\rm L}^{2}(\Omega)$ ensure that $\theta_{n}(r)=O(r^{|n|})$. As a result we have $\theta_{n}(r)=\alpha_{|n|}T_{|n|,\lambda}(r)$.

Finally, $\Theta$ being a real function, we can recombine the Fourier modes to get,

[TABLE]

We also proved that each term $T_{{n},\lambda}(r)\cos(n\varphi)$ or $T_{{n},\lambda}(r)\sin(n\varphi)$ itself is an eigenmode for $\lambda$ which obviously are linearly independent. But each eigenvalue $\lambda\in\Lambda$ being of finite multiplicity, the sums above are finite. ∎

2.2 Closure functions

Assuming the following decomposition:

[TABLE]

and formally injecting this expansion into problem (5) provides recursive relations on $t_{{n},p}(r)$,

[TABLE]

which allows an explicit analytical computation of the functions $t_{n,p}(r)$. We prove in section 2.3 that such a decomposition exists. The functions $t_{n,p}$ are called the closure functions. They are precisely defined in section 2.2.1 and their construction with the help of closure problems is given in section 2.2.2.

A consequence is that the spectrum in (8) and (9) are given by the zeros of the following analytical series

[TABLE]

where the coefficients $c_{n,p}$ are given by $c_{n,p}=t_{{n},p}(R)$ in the Dirichlet case or by $c_{n,p}=\diff{}{r}t_{{n},p}(R)$ in the Neumann case. In practice:

By truncating the series in equation (11) at order $M$, we can compute approximate eigenvalues by searching the zeros of the polynomial in $\lambda$ $\sum_{p=0}^{M}c_{n,p}\lambda^{p}=0$. 2. 2.

If $\bar{\lambda}$ is an approximate eigenvalue, the corresponding approximate eigenmode is $\sum_{p=0}^{M}t_{{n},p}(r)\bar{\lambda}^{p}$.

For more simplicity we fix in the sequel the value of $n\in\mathbb{N}$ and denote $t_{{n},p}=t_{p}$ and $T_{{n},\lambda}=T_{\lambda}$.

2.2.1 Definition

We consider the ODEs, for $j=1\dots m$,

[TABLE]

together with the transmission conditions for $j=1\dots m-1$,

[TABLE]

and the normalisation condition at the origin,

[TABLE]

Lemma 2**.**

Setting $t_{-1}=0$ and $t_{0}=r^{n}$, then the closure functions $(t_{p}(r))_{p\geq 1}$ satisfying (12), (13) and (14) for $p\geq 1$ are uniquely defined.

The proof is set-up by construction in section 2.2.2.

2.2.2 Construction

We assume that for some $p\geq 1$, $t_{p-2}(r)$ and $t_{p-1}(r)$ are known. We hereby derive $t_{p}(r)$. Let us first introduce the operators $F_{j}$ for $j=1,\dots m$, defined for a function $f$

[TABLE]

which is the inverse of operator $\Delta_{n}$. We denote $\psi_{1}(r)=r^{n}$ and $\psi_{2}(r)=r^{-n}$ if $n>0$ or $\psi_{2}(r)=\ln(r)$ if $n=0$, that are the basis solution of $\Delta_{n}f=0$. We consider the right hand side $f_{p-1}$

[TABLE]

Then on each compartment $(r_{j-1},r_{j})$, $t_{p}(r)$ is solution of (12) and therefore reads,

[TABLE]

We finally show how to compute the constants $\alpha_{j}$ and $\beta_{j}$

First compartment $[0,r_{1}]$

Assume that $t_{p-1}=O(r^{n})$ and $f_{p-1}=O(r^{n})$ at r=0, which is true for $p=1$.

We get that $F_{1}[f_{p-1}]=O(r^{n+2})$ and the normalisation condition (13) sets $\alpha_{1}=\beta_{1}=0$. We then have,

[TABLE]

It follows that $t_{p}=O(r^{n+2})=O(r^{n})$ and $f_{p}=O(r^{n})$.

Further compartments $[r_{j},r_{j+1}]$, $j\geq 1$

We assume that $t_{p}(r)$ has been computed on the compartment $[r_{j-1},r_{j}]$ and determine $t_{p}(r)$ on $[r_{j},r_{j+1}]$, $j\geq 1$.

We clearly have $F_{j+1}[f](r_{j})=0$ and $\diff{}{r}F_{j}[f](r_{j})=0$. Then equation (13) at $r_{j}$ reformulates as

[TABLE]

which equation has a unique solution since $\psi_{1}$ and $\psi_{2}$ form a basis for the solutions of the homogeneous equation $\Delta_{n}f=0$.

2.3 Series expansion of the eigenmodes

Our main result is the following.

Theorem 1**.**

The functions $T_{{n},\lambda}$ satisfy,

[TABLE]

where the $(t_{{n},p}(r))_{p\in\mathbb{N}}$ are the closure functions introduced in the previous section.

Proof of theorem 1

We fix the value of $n\in\mathbb{N}$ and simply denote $t_{{n},p}=t_{p}$ and $T_{{n},\lambda}=T_{\lambda}$.

Assume that the three functions $T_{\lambda}(r)$, $\diff{}{r}T_{\lambda}(r)$ and $\Delta_{n}T_{\lambda}(r)$ are analytical for $r\in[r_{j-1},r_{j}]$ and $\lambda\in\mathbb{C}$. We can write $T_{\lambda}(r)=\sum_{p\geq 0}s_{p}(r)\lambda^{p}$. The derivation theorem imply that $\diff{}{r}T_{\lambda}(r)=\sum_{p\geq 0}\diff{}{r}s_{p}(r)\lambda^{p}$ and that $\Delta_{n}T_{\lambda}(r)=\sum_{p\geq 0}\Delta_{n}s_{p}(r)\lambda^{p}$. Injecting this in (5) shows the $s_{p}(r)$ satisfy (5). Similarly the transmission and renormalisation conditions (6) (7) imply that the $s_{p}(r)$ satisfy (13) (14). Uniqueness in lemma 2 then imply that $s_{p}(r)=t_{p}(r)$.

Let us then prove that $T_{\lambda}(r)$, $\diff{}{r}T_{\lambda}(r)$ and $\Delta_{n}T_{\lambda}(r)$ are analytical for $r\in[r_{j-1},r_{j}]$ and $\lambda\in\mathbb{C}$ for all $j=1\dots m$. We proceed by induction.

Assume that this is true on $[r_{j-1},r_{j}]$. Then the initial data $\lambda\rightarrow T_{\lambda}(r_{j})$ and $\lambda\rightarrow\partial_{r}T_{\lambda}(r_{j})$ are analytical. On $[r_{j},r_{j+1}]$, $T_{\lambda}$ is the solution of the regular ODE (5) that analytically depends on $\lambda$, $r$ and whose initial conditions (6) at $r_{j}$ also analytically depend on $\lambda$. Classical results on ODEs (see e.g. [1, section 32.5]) state that $T_{\lambda}(r)$ analytically depends on $\lambda$ and $r$ on $[r_{j},r_{j+1}]$. This is also true for $\Delta_{n}T_{\lambda}$ since $\Delta_{n}T_{\lambda}=-\lambda^{2}T_{\lambda}+\lambda v/kT_{\lambda}$. Finally this is also true for $\diff{}{r}T_{\lambda}$ by integration.

It remains to prove the result for $r\in[0,r_{1}]$. This is harder because of the singularity at $r=0$. The problem being local at $r=0$, we can assume $r_{1}\leq 1.$ We formally introduce the series,

[TABLE]

Let us denote $F[f]=F_{1}[f]$ for $F_{1}[f]$ defined in (15). We introduce $F^{(i)}=F\circ\dots\circ F$ the $i^{th}$ iterate of $F$. Let us define $\tau_{i}=F^{(i)}[t_{0}]$ for $t_{0}(r)=r^{n}$ the 0th closure function. It is easy to compute $\tau_{i}$,

[TABLE]

We consider the constant $M=\max(\|v/k\|_{\infty},1)\geq 1$.

Lemma 3**.**

If $r_{1}\leq 1$, then on $[0,r_{1}]$ we have,

[TABLE]

with $\alpha_{p}=(2M)^{p}K_{i-1}$ for $p=2i$ or $p=2i+1$.

With definition (18) of the coefficients $K_{i}$, it is clear that the series $\sum_{p\geq 0}\alpha_{p}\lambda^{p}$ converges over $\mathbb{C}$. The three series $A_{\lambda}(r)$, $B_{\lambda}(r)$ and $C_{\lambda}(r)$ therefore are normally converging for $r\in[0,r_{1}]$ and for $\lambda$ in any compact in $\mathbb{C}$. As a result the integration theorem implies that $B_{\lambda}=\partial_{r}A_{\lambda}$ and $C_{\lambda}=\Delta_{n}A_{\lambda}$. Relation (12) ensures that $A_{\lambda}$ satisfies (5) whereas relation (14) together with $t_{0}=r^{n}$ ensures that $A_{\lambda}$ satisfies (7). Uniqueness in lemma 1 then implies that $A_{\lambda}=T_{\lambda}$. This proves theorem 1 on $[0,R]$ and ends this proof.

Proof of lemma 3.

We will systematically use that $r\leq 1$, that $K_{i}$ and $\tau_{i}(r)$ in (18) are decreasing and that the operator $F$ satisfies,

[TABLE]

With definitions (16)-(17) we have the upper bound,

[TABLE]

By recursion, we obtain an upper bound involving the $\tau_{i}=F^{(i)}(t_{0})$ of the form

[TABLE]

The number of terms in the sum is less than $2^{p}$. Index $n_{k}$ is smaller than $p$ and $M^{n_{k}}\leq M^{p}$. The minimal value for $m_{k}$ is $i$ if $p=2i$ or $i+1$ if $p=2i+1$, so that $\tau_{m_{k}}\leq\tau_{i}$ or $\tau_{m_{k}}\leq\tau_{i+1}$ respectively. Therefore,

[TABLE]

which upper bound ensures the first inequality in lemma 3.

From that last inequality it is easy to check that $|t_{p}|+|t_{p-1}|\leq 2(2M)^{p}K_{i}r^{n+2i}$ if $p=2i$ or $p=2i+1$. For $p=2i$ or $p=2i+1$ it follows that

[TABLE]

This gives the third inequality in lemma 3.

By differentiating (17) we get,

[TABLE]

It follows that

[TABLE]

with,

[TABLE]

implying the second inequality in lemma 3. ∎

2.4 Extension to planar configurations

We consider layered planar configurations as depicted on Figure 2. The transverse coordinate perpendicular to the layers is denoted by $x$. The coordinate $x$ is homologue to the radial coordinate $r$ in the cylindrical case. The origin is set at the center so that $-R\leq x\leq R$ with $2R$ the total thickness of the geometry.

Actually, the results that we obtained for concentric cylindrical configurations are easier to establish in the case of layered planar configurations. This is because the operator \Delta_{n}:=\big{(}\diff{{}^{2}}{x^{2}}+n^{2}\big{)} associated with the $y$-periodic decomposition

[TABLE]

is no more singular in Cartesian coordinates. Hence, the technical issues associated with the proof of analyticity in the variable $\lambda$ for the functions $T_{\lambda}$, $\,\mathrm{d}T_{\lambda}/\,\mathrm{d}r$ and $\Delta_{n}T_{\lambda}$ are no longer present in this case. Furthermore, each step of the proofs provided in sections 2.1 and 2.2 directly apply to the planar case, so that theorem 1 also holds.

3 Exemples of applications

In this section we develop various exemples of solutions so as to illustrate the versatility and usefulness of the previously presented theoretical results. In section 3.1 we first give explicit general solutions adapted for two families of geometries, i.e planar or cylindrical, for general boundary conditions. We pursue towards illustrating interesting and relevant solutions considering two idealyzed but non trivial configurations in the subsequent sections. In section 3.2 we showcase how a localized heat source can lead to a ’hot spot’ of temperature in its neibourghood, and illustrate how our mesh-less analytical method can effectively capture the temperature peak. A second exemple is provided in section 3.3 where we examine a double-pass configuration in the planar framework for which, again, a localized heat source is imposed nearby the origin.

3.1 Explicit families of solutions

As in equation (2), we will consider symmetric boundary conditions (only depending on $z$). Thus we will consider the spectrum $\Lambda_{0}$ in definitions (8)-(9) for $n$=0. In the Dirichlet case $\Lambda_{0}=\left\{\lambda\in\mathbb{C},\quad T_{0,\lambda}(R)=0\right\}$ and in the Neuman case $\Lambda_{0}=\left\{\lambda\in\mathbb{C},\quad\,\mathrm{d}T_{0,\lambda}/\,\mathrm{d}r(R)=0\right\}$. The spectrum is computed with the closure functions as in equation (11). It decomposes as in equation (3): $\Lambda_{0}=\{\lambda_{+i},\leavevmode\nobreak\ \lambda_{-i},\leavevmode\nobreak\ i\in\mathbb{N}^{\star}\}$ with $\lambda_{+i}<0$ the upstream modes and $\lambda_{-i}>0$ the downstream modes. We will simply denote $T_{\pm i}=T_{\lambda_{\pm i},0}$. Remember that $\lambda_{\pm i}$ and $T_{\pm i}$ depend on the nature of the boundary condition (Dirichlet or Neumann).

Dirichlet boundary condition

For the lateral Dirichlet boundary condition in equation (2), the temperature solution are given in [2] for the cylindrical case

[TABLE]

with, denoting $k$ the conductivity in the boundary annulus:

[TABLE]

This adapts to the parallel planar configuration with

[TABLE]

In both cylindrical and planar cases, the functions $c_{i}(z)e^{\lambda_{i}z}$ are given by the convolution between $\,\mathrm{d}g/\,\mathrm{d}z$ and the exponentially decaying modes

[TABLE]

for the upstream modes and downstream modes respectively.

Neumann boundary condition &

non-equilibrated case

Consider now a Neumann boundary conditions (2) in the case where $Q:=\int_{\Omega}v\,\mathrm{d}x\neq 0$, i.e. the total convective flux is not zero. Then from [2] the solution reads

[TABLE]

with $G(z)=\int_{-\infty}^{z}g(\xi)d\xi$ the primitive of the heat source $g(z)$ and $P$ the perimeter of the external cylinder. Note that the temperature indeed is defined up to an additive constant that has been fixed by setting $T_{-\infty}=0$ here.

For the cylindrical configuration, we choose a Poiseuille velocity profile $v(r)={\rm Pe}(1-(r/r_{0})^{2})$, where ${\rm Pe}$ is the Péclet number that quantifies the ratio between convection and diffusion (here based on the maximal velocity in the tube). We have $P/Q=4R/({\rm Pe}r_{0}^{2})$ and

[TABLE]

Whereas, for parallel planar configurations:

[TABLE]

In both cases, the functions $c_{i}(z)e^{\lambda_{i}z}$ are given by the convolution between the imposed flux at the boundary and the exponentially decaying modes

[TABLE]

for the upstream modes and downstream modes respectively.

Neumann boundary condition & equilibrated case

Consider now a Neumann boundary conditions (2) in the case where the total convective flux cancels out: $Q:=\int_{\Omega}v\,\mathrm{d}x=0$ . This is the case of an equilibrated exchanger. In this case, the solution displays a distinct form (see [2]) involving the (adiabatic) kernel $T_{0}$ solution of

[TABLE]

In the section 3.3 we will consider an exemple of such a configuration for which we will give an explicit solution of the kernel $T_{0}$. In general form, the complete solution associated with equilibrated case $Q=0$, reads

[TABLE]

with $\mathcal{G}(z)=\int_{-\infty}^{z}G(\xi)d\xi$, the second primitive of the heat source $g(z)$, $\alpha_{i}$ and $c_{i}(z)$ again given by (20) and (23) and where $a$ and $b$ are two constants characterizing the heat exchange with values detailed below. Note that for this configuration the temperature field is defined up to $C_{1}(z+T_{0})+C_{2}$, see details in [2].

In cylindrical configuration the parameters $a$ and $b$ are given by

[TABLE]

whereas for parallel planar configuration the parameters $a$ and $b$ read

[TABLE]

3.2 Locally heated pipe & non-equilibrated case $Q\neq 0$

We illustrate the use of explicit computation of the eigenmode decomposition, through the recursive relations (12) and (13), in a simple and classical configuration: sometimes referred to as ’generalized Graetz’ configuration. Two concentric cylinders are thus considered. A central one, for which $r\in[0,r_{0}]$ and whereby the fluid convects the temperature, and an external one, $r\in[r_{0},R]$ where temperature conduction occurs. The dimensionless axi-symmetric longitudinal velocity $v(r)$ inside the inner cylinder is chosen such as $v(r)={\rm Pe}(1-(r/r_{0})^{2})$, where ${\rm Pe}$ is the Péclet number which quantifies the ratio between convection and diffusion. The domain dimensions are $r_{0}=1$ and $R=2$. The conductivity is set to $k=1$. The solution is defined up to an additive constant that is fixed by setting $T_{-\infty}=0$. A Neumann boundary condition $k\nabla T=g(z)$ is set. The applied boundary condition is chosen so as to present a localyzed (and regular) heat flux nearby the origin, with $z_{0}=1/2$ here:

[TABLE]

and $g(z)=0$ otherwise. With these conditions, a simple balance on the domain allows to compute $T_{+\infty}$:

[TABLE]

so that $T_{+\infty}=8/{\rm Pe}$ here. Using Neumann boundary condition (28) and equation (20) one is able to provide a mesh-less explicit analytical solution for the temperature, illustrated in Figure 3 for various values of ${\rm Pe}$ varying between $100$ to $0.1$ so as to show-case the drastic effect of convection on the temperature profiles. Figure 3a exemplifies that, when convection dominates in the centerline $r=0$, the effect of the heat source nearby the origin is weak. The local temperature is almost zero at $r=0$ for $z\in[-1,0]$, since the prescribed temperature at $z\rightarrow-\infty$ is zero. Nevertheless, a slight tilt of the centerline temperature profile is noticeable as $z>0$ so that it barely reaches the non-zero asymptotic downstream constant temperature $T_{+\infty}$ at $z=10$. On the contrary to the centerline profile, the wall profile at $r=R$ displays a strong deflection with a maximum located at the heat source maximum $z=1/2$, and both upstream and downstream decay from this maximum. The typical downstream decay length is related to the convection ability to transport the heat flux downstream. Hence the larger the Péclet, the longer the downstream decay length. The upstream decay length, on the contrary both depends on the solid conduction and the wall radius. In the case of small solid walls thickness, some asymptotic behavior have been documented [19]. The other radially intermediate temperature profiles shown in Figure 3a display a medium behavior between the centerline and the wall profile. The closer to the outer cylinder wall, the closer the temperature peak to the wall profile. Figures 3b, 3c and 3d display the effect of decreasing the convection on the temperature profile.

From one hand, these profiles display smoother and smaller peaks at the heat source as convective effects are weakened. On the other hand, the profiles are increasingly non-symetric at smaller Péclet numbers, with an increasing downstream temperature $T_{+\infty}=8/{\rm Pe}$.

3.3 Parallel configuration & equilibrated case $Q=0$

Here we consider a parallel planar geometry in a double-pass configuration for which the upper fluid is re-injected into the lower one at one end as in [12]. An exchanger with total thickness $2R$ is considered. A fluid is flowing for $|x|\leq x_{0}$ surrounded by solid walls for $x_{0}\leq|x|\leq R$. We consider the zero total flux for which the upper fluid is convected along $+z$ direction for $x\in[0,x_{0}]$, and on the opposite one for $x\in[-x_{0},0]$. Within $[-x_{0},x_{0}]$, the velocity profile reads

[TABLE]

$\sigma{(x)}$ being the sign of $x$, with Péclet number ${\rm Pe}=\bar{v}x_{0}/D$ (built from the average velocity $\bar{v}=\int_{0}^{x_{0}}v\,\,\mathrm{d}x/x_{0}$, $x_{0}$ the fluid channel half-gap and the diffusivity $D$). At $x=\pm R$, adiabatic conditions are prescribed, (i.e $\nabla T\cdot{\bf n}|_{R}=0$) for $|z|>1/2$ whereas the flux (28) (with $z_{o}=0$ here) is imposed for $z\in[-1/2,1/2]$. In this case the adiabatic kernel $T_{0}$ solution of (24) is given by: for $|x|\leq x_{0}$

[TABLE]

whereas for $|x|\geq x_{0}$

[TABLE]

The two constants $a$ and $b$ defined in (27) read

[TABLE]

Figure 4 illustrates the temperature profiles along the longitudinal direction $z$ at various transverse heights $x$, either in the center of the channel ($x=0$), at the interface between the liquid and the solid ($x=x_{0}$) or at the solid exterior edge ($x=R$). One can observe that the “hot-spot” temperature located very close at $z=0$ at the surface $x=R$, is weakly affected by the increase of the Péclet number expect for small Péclet (where one has to translate back the reference temperature chosen at $\pm\infty$, so as to obtain true physical “hot-spot” temperature). Nevertheless, further-down inside the solid the temperature rise is weakened by increasing fluid convection, as expected. Also, convection drops down the outlet temperature, as expected from heat-flux balance argument.

4 Conclusion

This contribution has provided the mathematical proof, as well as the effective algorithmic framework for the computation of generalized Graetz mode decomposition in cylindrical or parallel configurations. We have shown that, in these special configurations, the Graetz functions analyticity enables mesh-less explicit computation of the steady-state temperature even when boundary condition with source terms are considered. The method has been illustrated in two complementary cases (cylindrical/non equilibrate and parallel/equilibrated) in order to showcase its various aspects.

5 Appendix

5.1 Cylindrical heated pipe case $Q\neq 0$

The solution provided by (20) is

[TABLE]

Rewritting (28) as

[TABLE]

with $H(z)$ the Heaviside function, and using integration by parts leads to the primitive $G(z)=\int_{-\infty}^{z}g(z^{\prime})dz^{\prime}$ equals to

[TABLE]

The function $c_{i}(z)$ in (33) are given by (23), the integration by part of which gives

[TABLE]

The eigenfunctions $T_{i}$ are provided by the $\lambda$-analytical decomposition (10) upon functions $t_{p}(r)$ such that

[TABLE]

where each eigenvalue $\lambda_{i}$ of the discrete spectrum sets its eigenfunctions $T_{i}$ from (34). We hereby provide the first three elements of both dowstream and upstream spectrum computed with a finite truncation of $N_{p}=20$ in (34) and parameter $Pe=1$, with a formal calculus Maple software: $\lambda_{1}=0.674240$, $\lambda_{2}=3.306258$, $\lambda_{3}=4.936416$, $\lambda_{-1}=0$, $\lambda_{-2}=-1.027741$, $\lambda_{-3}=-2.35726$. Function $t_{p}(r)$, $p\in\{0,5\}$ are also hereby given by the following piecewise continuous analytical functions of $r$ along the fluid-solid domains $r\in[0,1]\cup[1,2]$

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Finally, each parameter $\alpha_{i}$ of (33) is given by (22) using the closure function $T_{i}(R=2)$ and its corresponding eigenvalue $\lambda_{i}$.

5.2 Parallel configuration & equilibrated case $Q=0$

The theoretical solution detailed in section 3.1 is hereby detailled. From (25) we recall the temperature solution

[TABLE]

involving the constant $a$ given in (32) and the function $g(z)$ given in (28). Rewritting (28) as

[TABLE]

with $H(z)$ the Heaviside function, and using integration by parts leads to a primitive $G(z)=\int_{-\infty}^{z}g(z^{\prime})dz^{\prime}$ equals to

[TABLE]

Again, the functions $c_{i}(z)$ are given by (23), the integration by part of which gives

[TABLE]

The eigenfunctions $T_{i}$ are provided by the $\lambda$-analytical decomposition (10) upon functions $t_{p}(r)$ such that

[TABLE]

where each eigenvalue $\lambda_{i}$ of the discrete spectrum sets its eigenfunctions $T_{i}$ from (34). We hereby provide the five first elements of these spectrum computed with a finite truncation of $N_{p}=20$ in (34) and parameter $Pe=50$, computed with a formal calculus Maple software. $\lambda_{1}=-1.738793$, $\lambda_{2}=-1.738793$, $\lambda_{3}=-1.585275$, $\lambda_{4}=-1.3093020$, $\lambda_{5}=-1.011529$. Function $t_{p}(r)$, $p\in\{0,5\}$ are also hereby given by the following piecewise continuous polynomial functions of $r$ along the various solid-fluid domains $[-2,-1]\cup[-1,0]\cup[0,1]\cup[1,2]$. Starting with $t_{0}=1$ identically equal to 1, we recursively compute the following functions $t_{i}$ and obtained

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Parameter $\alpha_{i}$ of (36) is given by (22) using closure function $T_{i}(R=2)$ and its corresponding eigenvalue $\lambda_{i}$.