Alexey Ph Polikarpov

Irina A Graur

Heat and mass transfer in a rarefied gas confined between its two parallel condensed phases

The heat and mass transfer in a rarefied gas between its two parallel condensed phases is considered on the basis of linearized and non-linear S-model kinetic equations. The profiles of the macroscopic parameters in the gap between gas-liquid interfaces are obtained for several Knudsen numbers and for the cases of complete and non-complete evaporation and condensation. The linearized Navier-Stokes equations and energy equation, subjected to the temperature and pressure jump boundary conditions, are solved analytically and the expressions for the macroscopic parameters are obtained. The comparison of three approaches allowed us to establish the limit of the application of the linearized approaches in term of the saturation temperature ratio between the condensed phases.

Introduction

Evaporation and condensation phenomena are interesting not only from the fundamental point of view, but they are also largely present in various industrial applications. During the evaporation and condensation the liquid-vapor interface is in non-equilibrium state and its detailed modeling requires the use of the molecular based approaches, like the kinetic theory of gases [START_REF] Sone | Kinetic Theory and Fluid Mechanics[END_REF] or the molecular dynamics [START_REF] Meland | Molecular dynamics simulation of the inverted temperature gradient phenomenon[END_REF]. The adequate description of the evaporation-condensation was in focus of interest from a long time [START_REF] Hertz | Uber die verdungstung der flussigkeiten, inbesondere des quecksilbers, im luftleeren raume[END_REF], [START_REF] Knudsen | Die maximale verdampfungsgeschwindigkeit des que-cksilbers[END_REF], [START_REF] Schrage | A Theoretical Study of Interphase Mass Transfer[END_REF], [START_REF] Pao | Application of kinetic theory to the problem of evaporation and condensation[END_REF], [START_REF] Pao | Temperature and density jumps in the kinetic theory of gases and vapors[END_REF], [START_REF] Aoki | Gas flows caused by evaporation and condensation on two parallel condensed phases and the negative temperature gradient: numerical analysis by using a nonlinear kinetic equation[END_REF], [START_REF] Sone | Kinetic Theory and Fluid Mechanics[END_REF]. However, the complete understanding of these phenomena is just far to be completely achieved. In addition, the implementation of the kinetic models or the molecular dynamics based approaches is not so easy for the practical problems. This is why it is also interesting to use the continuum approaches, like the classical Navier-Stokes equations, subjected with the temperature and pressure jumps boundary conditions [START_REF] Cipolla | Kinetic theory of condensation and evaporation. II[END_REF], [START_REF] Onishi | Kinetic theory of slightly strong evaporation and condensation -hydrodynamic equation and slip boundary condition for finite reynolds numberevaporation and condensation -hydrodynamic equation and slip boundary condition[END_REF] to simulate the evaporation-condensation problems. Recently the alternative macroscopic approach based on the regularized 13 moments equations [START_REF] Struchtrup | Evaporation boundary conditions for the r13 equations of rarefied gas dynamics[END_REF] was also proposed and applied for the simulation of the evaporation and condensation phenomena.

The main objective of this study is the development of the numerical approach to simulate the behavior of the vapor phase between two parallel liquidvapor interfaces. Contrarily to the authors of Ref. [START_REF] Aoki | Gas flows caused by evaporation and condensation on two parallel condensed phases and the negative temperature gradient: numerical analysis by using a nonlinear kinetic equation[END_REF], where the BGK equation is implemented, the linearized and non-linear S-model [START_REF] Shakhov | Generalization of the Krook kinetic relaxation equation[END_REF] kinetic equations are solved numerically to simulate the evaporation-condensation phenomena appearing between these two plane condensed phases. The main advantage of the S-model is its capacity to provide the correct Prandtl number, equals to 2/3, for the monatomic gases, which can be important, when the problems involving the vapor flows and the heat flux are considered. In addition, the linearized Navier-Stokes equations and energy equation subjected to the temperature and pressure jumps boundary conditions [START_REF] Cipolla | Kinetic theory of condensation and evaporation. II[END_REF], [START_REF] Onishi | Kinetic theory of slightly strong evaporation and condensation -hydrodynamic equation and slip boundary condition for finite reynolds numberevaporation and condensation -hydrodynamic equation and slip boundary condition[END_REF] are also solved analytically and the simple relations are obtained for the profiles of the macroscopic parameters of vapor in the gap between two condensed phases. The influence of evaporation-condensation coefficient on the macroscopic parameters is also analyzed. Comparing the results obtained by three approaches we establish the limit of the applicability of the linearized kinetic equation as well as of the linearized Navier-Stokes equations for the simulation of the evaporation and condensation phenomena.

Problem statement

We consider two parallel plane condensed phases at rest, maintained at temperature T 1 and T 2 , on the bottom (y = 0) and top (y = H) interfaces, respectively, y is the coordinate normal to both condensed phases. Let p 1 and p 2 be the saturation gas pressures at temperatures T 1 and T 2 , (T 1 > T 2 ), respectively. We investigate here the behavior of the monatomic gas motion caused by evaporation and condensation on the condensed phases, first on the basis of the kinetic theory.

If we consider that the ratios p 1 /p 2 and T 1 /T 2 are independent, then the evaporation-condensation problem is completely characterized by three parameters: p 1 /p 2 and T 1 /T 2 ratios and by the rarefaction parameter δ. However, in practice, the saturation pressure and temperature are related by the Clausius-Clapeyron equation [START_REF] Landau | Statistical physics[END_REF] and therefore cannot take the arbitrary values. In this case the evaporation-condensation problem depends on two parameters only: the slop of p (T ) function and the rarefaction parameter δ. In the following we consider both cases: first the situation when two parameters are independent, and then, the cases of Argon and several other monatomic gases, where the saturation pressure and temperature are related through the Clausius-Clapeyron equation.

The distance H between two interfaces is taken as the characteristic dimension of the problem, so the rarefaction parameter δ is defined as following: where is the equivalent mean free path:

δ = H , (1)
= µ 2 v 2 p 2 , (2)
which is defined using the reference parameters with subscript 2, corresponding to the upper plate. The choice of the parameters related to the upper plate as the reference parameters was done without lost of generality. In Eq. ( 2)

µ 2 = µ(T 2 )
is the dynamic viscosity of the vapor phase

µ(T ) = µ 2 T T 2 , ( 3
)
v 2 is the most probable molecular speed,

v(T ) = √ 2RT , (4)
calculated at the temperature

T 2 : v 2 = v(T 2 ), R = k B /m
is the specific gas constant, k B is the Boltzmann constant, m is the molecular mass.

Modeling based on the kinetic equation

To model the evaporation and condensation phenomena between two parallel condensed phases maintained at different temperatures first the S-model kinetic equation [START_REF] Shakhov | Generalization of the Krook kinetic relaxation equation[END_REF] is used. The evaporation-condensation problem between two parallel condensed phases is considered here as steady state and one dimensional in physical space, so the S-model kinetic equation is written as

v y ∂f ∂y = υ f S -f , (5)
where f (y , v) is the one particle velocity distribution function, v = (v x , v y , v z ) is the molecular velocity vector, υ is the collision frequency, υ = p /µ. In the frame of S-model the equilibrium distribution function f S in Eq. ( 5) is defined as following

f S (y , v) = f M 1 + 2mVq 15n (y )(k B T (y )) 2 mV 2 2k B T (y ) - 5 2 , (6)
here T (y ) is a gas temperature, n (y ) is a gas number density, u = (0, u , 0) is a bulk velocity vector, V = vu is the peculiar velocity vector, q = (0, q , 0) is a heat flux vector, f M is the Maxwellian distribution function [START_REF] Kogan | Rarefied gas dynamics[END_REF]. The macroscopic parameters are defined as follows:

n (y ) = f (y , v)dv, u (y ) = 1 n f (y , v)v y dv, (7)
T (y ) = m 3k B n f (y , v)V 2 dv, q (y ) = m 2 f (y , v)V 2 (v y -u y )dv . (8)
By integrating Eq. ( 5) multiplied by (1, v y , v 2 i ) we obtain the integrals expressing the number of molecules, y component of momentum and energy transported in y direction per a unit area of a plan, parallel to the condensed phases. The evaporation flow rate, expressed in the number of molecules per time and per surface unit, J n , and the evaporation mass flow rate, expressed in kilogram per time and surface units, J ρ , are defined as followings:

J n = v y f (y , v)dv, J ρ = m v y f (y , v)dv. (9)
The second definition of the evaporation mass flow rate is usually provided from the experiments. The energy flux is defined as

J E = v y v 2 f (y , v)dv. ( 10
)
The constancy of the integrals J n and J E will be further used to estimate the accuracy of the numerical calculations and for the comparison with the results of Ref. [START_REF] Aoki | Gas flows caused by evaporation and condensation on two parallel condensed phases and the negative temperature gradient: numerical analysis by using a nonlinear kinetic equation[END_REF].

Linearized S-model equation

To linearize S-model kinetic equation [START_REF] Schrage | A Theoretical Study of Interphase Mass Transfer[END_REF] we assume that deviations between the temperatures of two condensed surfaces and the corresponding saturation pressures are small:

X P = p 1 -p 2 p 2 1, X T = T 1 -T 2 T 2 1. (11)
In previous expressions X P and X T can be related to the thermodynamic forces [START_REF] Sharipov | Onsager-casimir reciprocity relations for open gaseous systems at arbitrary rarefaction. II Application of the theory for single gas[END_REF]. As it was mentioned in Section 2, for a given gas the pressure and temperature differences are coupled by the relation

p 1 -p 2 = β(T 1 -T 2 ), ( 12
)
where β is a positive constant corresponding to the slop of the Clausius-Clapeyron curve at T 2 , so X P and X T are not independent quantities. However here we consider two forces separately, to see clearly the impact of each force to the evaporation-condensation process.

For further derivation we introduce the following dimensionless quantities:

y = y H , c = v v 2 , u = u v 2 , n = n n 2 , T = T T 2 , p = p p 2 , q = q p 2 v 2 .
(13) When X P and X T are small enough compared to 1 we can linearize the distribution function as following:

f = f M 0 (1 + h P X P + h T X T ), (14)
where h P and h T are the perturbation functions, related to the pressure and temperature difference, respectively, f M 2 is the absolute Maxwellian distribution function with the reference number density, n 2 and temperature, T 2 :

f M 2 = n 2 m 2πk B T 2 3/2 exp(-c 2 ). (15)
For the real gases, for any value of β we can obtain the solution as a superposition of two functions as h = (βh

P + h T )∆T, (16)
where ∆T = T 1 -T 2 . However, as it was mentioned above, we consider here X P and X T as independent quantities. The S-model kinetic equation [START_REF] Schrage | A Theoretical Study of Interphase Mass Transfer[END_REF] in its linearized form can be separated into two equations with reference to the thermodynamic forces X P and X T :

c y ∂h i ∂y = δ ν i + 2c y u i + c 2 - 5 2 τ i + 4 15 q i c y -h i , i = P, T. (17)
The macroscopic parameters in Eq. ( 17) can be calculated through the following relations:

n = 1 + ν P X P + ν T X T , u = u P X P + u T X T , (18)
T = 1 + τ P X P + τ T X T , q = q P X P + q T X T . (19)
The dimensionless moments of the distribution function are given as:

ν i = 1 π 3/2 exp(-c 2 )h i dc , (20)
u i = 1 π 3/2 exp(-c 2 )h i c y dc , (21)
τ i = 1 π 3/2 exp(-c 2 )h i 2 3 c 2 -1 dc , (22)
q i = 1 π 3/2 exp(-c 2 )h i c y c 2 - 5 2 dc , i = P, T. (23)
The distribution function of the molecules outgoing from the condensed phases is supposed to be decomposed into two parts [START_REF] Aoki | Gas flows caused by evaporation and condensation on two parallel condensed phases and the negative temperature gradient: numerical analysis by using a nonlinear kinetic equation[END_REF]: σ e part of molecules evaporates from the surface with saturation parameters while another part, (1-σ c ), reflects diffusively from the surface. We assume here that the evaporation σ e and the condensation σ c coefficients are equal each to other σ = σ e = σ c . The linearized form of evaporation-condensation boundary conditions reads: for the bottom surface, y = 0:

h P = σ + (1 -σ)ν Pw , h T = σ c 2 - 5 2 + (1 -σ) c 2 - 3 2 + ν Tw , (24)
ν Pw = - 2 π cy<0 h P c y dc, ν Tw = - 1 2 - 2 π cy<0 h T c y dc (25)
and upper surface, y = 1:

h P = (1 -σ)ν Pw , h T = (1 -σ)ν Tw , (26)
ν Pw = 2 π cy>0 h P c y dc, ν Tw = 2 π cy>0 h T c y dc. (27)
To reduce the number of variables in a velocity space we use the standard procedure [START_REF] Chu | Kinetic-theoretic description of the formation of a shock wave[END_REF]: we multiplied first Eqs. [START_REF] Broadwell | Shock structure in a simple discrete velocity gas[END_REF] successively by 1 and (c 2 z -1 2 ), and then we integrated them over c z velocity. Doing so we obtain the system of two equations for reduced distribution functions:

c y ∂φ i ∂y = δ ν i + 2c y u i + c 2 x + c 2 y -1 τ i + 4 15 q i c y c 2 x + c 2 y -2 -h i , (28)
c y ∂ψ i ∂y = δ 1 2 τ i + 2 15 q i c y , i = T, P, (29)
where these reduced distribution functions read:

φ i = 1 √ π ∞ -∞ exp(-c 2 z )h i dc z , ψ i = 1 √ π ∞ -∞ exp(-c 2 z ) c 2 z - 1 2 h i dc z , i = T, P . (30)
The system of equations ( 28), (29) with the boundary conditions, Eqs. ( 24) -(27), is solved by the Discrete Velocity Method (DVM) [START_REF] Broadwell | Shock structure in a simple discrete velocity gas[END_REF], see Section 4.

Thermodynamic analysis

The author of Refs. [START_REF] Sharipov | Onsager-casimir reciprocity relations for open gaseous systems at arbitrary rarefaction. I General theory for single gas[END_REF], [START_REF] Sharipov | Onsager-casimir reciprocity relations for open gaseous systems at arbitrary rarefaction. II Application of the theory for single gas[END_REF] have obtained the Onsager-Casimir reciprocity relation in general form for rarefied gas flows. It was shown [START_REF] Sharipov | Onsager-casimir reciprocity relations for open gaseous systems at arbitrary rarefaction. II Application of the theory for single gas[END_REF] that for the case of evaporation-condensation of a gas between two parallel condensed phases the "fluxes" corresponding to the driving "forces" are defined as follows:

J P = n 2 u , J T = 1 k B T 2 q , (31)
where u y and q y do not depend on y . The thermodynamic fluxes are related to the thermodynamic forces in the matrix form:

J P J T = Λ P P Λ P T Λ T P Λ T T × X P X T . ( 32
)
In this case, the Onsager-Casimir relation, Λ P T = Λ T P , yields the coupling between the mass flux caused by temperature drop and the thermal flux caused by the pressure drop [START_REF] Sharipov | Onsager-casimir reciprocity relations for open gaseous systems at arbitrary rarefaction. II Application of the theory for single gas[END_REF].

Non-linear S-model equation

To establish the limit of applicability of linearized S-model kinetic equation to model two condensed surfaces problem we solve numerically the non-linear S-model kinetic equation. By using Eqs. ( 13) the dimensionless form of Eq. ( 5) reads:

c y ∂f ∂y = δn √ T f S -f . (33)
As for the linearized case, see Section 3.1, we assume that only one part σ of the incident molecules evaporates immediately and (1 -σ) part reflects from the surface diffusively. Therefore the dimensionless boundary conditions for the distribution function of the reflected molecules at the liquid-gas interfaces can be written as

y = 0 c y > 0, f = (σn s1 + (1 -σ)n r1 )f M 1 , f M 1 = 1 π 3/2 p 1 p 2 T 2 T 1 5/2 exp - T 2 T 1 c 2 y = 1 c y < 0, f = (σn s2 + (1 -σ)n r2 )f M 2 , f M 2 = 1 π 3/2 exp -c 2 . ( 34
)
Here n si , i = 1, 2 is the number density of the saturation vapor near each plane interface, calculated from the equation of state using the values of the saturation temperature and corresponding saturation pressure. The number density n ri , i = 1, 2, can be calculated from the impermeability condition on the corresponding interface:

n r1 = -2 π T 1 T 2 cy<0 c y f dc, n r2 = 2 √ π cy>0 c y f dc. ( 35
)
To minimize computational efforts, the c z variable is eliminated by introducing the reduced distribution functions as in [START_REF] Graur | Comparison of different kinetic models for the heat transfer problem[END_REF].

Numerical technique

To solve numerically both linearized and non-linear S-model kinetic equation the same Discrete Velocity Method (DVM) [START_REF] Graur | Comparison of different kinetic models for the heat transfer problem[END_REF] is used. First, the Discrete Velocity Method is applied to split the continuum (-∞, ∞) molecular velocity space c x , c y in the governing equation ( 5) into discrete velocity sets c xm , c y k where m, k = -N c , .., -1, 1, 2, .., N c . These velocities, c xm and c y k , are taken to be the roots of the Hermite polynomial of order N c . Then, the set of 2N c kinetic equations, corresponding to 2N c values of discrete velocities, is discretized in space by Finite Difference Method (FDM). In presented here simulations N c is taken to be equal to 20, so the number of points in the molecular velocity space is equal to 40 in each direction. The grid-independence is checked with finer grid with N c = 25 points: the difference in all macroscopic quantities was less than 1%.

The spacial derivatives are approximated by the first-order accurate upwind type scheme, see Ref. [START_REF] Graur | Comparison of different kinetic models for the heat transfer problem[END_REF] for more details. The number of uniformly distributed numerical grid points in physical space, N y , is equal to 1000. This number of grid points insures the calculation accuracy of 1%. The convergence criterion for iteration process is employed in following form

min i q (l+1) i -q (l) i / q (l) i ≤ ε, ( 36
)
where l is the iteration number, i = 1, N y ; ε = 10 -7 is used here. The described above numerical approach was used to solve the linearized and non-linear S-model kinetic equation for the large range of the rarefaction parameter and for several values of the evaporation coefficients. The obtained results are discussed in Section 6.

Solution of the linearized Navier-Stokes equations with the jump boundary conditions

The numerical solution of the model kinetic equation provide very accurate information about the evaporation-condensation process. However, from the practical point of view it is interesting to have the simple and easy to use estimation of the vapor behaviors. Below the analytical solution of the linearized Navier-Stokes equations with the energy equation subjected with the temperature and pressure jump boundary conditions is derived.

For one dimensional in the physical space case, considered here, the linearized Navier-Stokes equations, completed with the energy equations read:

du y dy = 0, dp dy = 0, d 2 T dy 2 = 0. ( 37
)
To take into account the evaporation-condensation on the upper and down interfaces, equations (37) have to be subjected to the special boundary conditions [START_REF] Cipolla | Kinetic theory of condensation and evaporation. II[END_REF], [START_REF] Lang | Evaporation and condensation for general gas-liquid surface scattering[END_REF], [START_REF] Sone | Kinetic theory of evaporation and condensation -hydrodynamic equation and slip boundary condition[END_REF], [START_REF] Onishi | Kinetic theory of slightly strong evaporation and condensation -hydrodynamic equation and slip boundary condition for finite reynolds numberevaporation and condensation -hydrodynamic equation and slip boundary condition[END_REF]. We will use here the most complete form of the boundary conditions suggested by the authors of Ref. [START_REF] Cipolla | Kinetic theory of condensation and evaporation. II[END_REF] and which relate the temperature and pressure jumps on the liquid-vapor interface, the driving "forces", with the "fluxes", evaporation rate and heat flux through interface. For the down plate this relation between "fluxes" and "forces" reads:

(p 1 -p )/p 1 (T 1 -T )/T 1 = a 11 a 12 a 21 a 22 × n 1 u q /(k B T 1 ) . ( 38
)
The heat flux in system (38) reads:

q = -κ dT dy , ( 39
)
where κ is the thermal conductivity coefficient, which can be calculated as:

κ = γ γ -1 1 P r Rµ, ( 40
)
where γ is the heat capacity ratio, P r is the Prandtl number. For the monatomic gas considered here, γ = 5/3 and P r = 2/3. The matrix in Eq. ( 38) is the inverse Onsager matrix, so it is symmetric and a 12 = a 21 . For further derivation it is convenient to introduce the reduced form:

a ij = a ij /(n 1 v 1 ) (41)
The values of these coefficients were calculated in Refs. [START_REF] Cipolla | Kinetic theory of condensation and evaporation. II[END_REF], [START_REF] Lang | Evaporation and condensation for general gas-liquid surface scattering[END_REF], [START_REF] Sone | Kinetic theory of evaporation and condensation -hydrodynamic equation and slip boundary condition[END_REF], [START_REF] Onishi | Kinetic theory of slightly strong evaporation and condensation -hydrodynamic equation and slip boundary condition for finite reynolds numberevaporation and condensation -hydrodynamic equation and slip boundary condition[END_REF] for the monatomic gases and for the case of the diffuse reflection of the molecules which do not condense on the surface:

a 11 = 2 √ π σ -1.4195228, a 12 = a 21 = 0.446658, a 22 = 1.0422028. ( 42
)
It is worth to note that only a 11 coefficient depends on the evaporation coefficient σ.

In our case of the evaporation-condensation between two parallel condensed surfaces maintained at different temperatures the boundary conditions, Eqs. (38), are write as following:

for the bottom condensed surface y = 0:

p = p 1 1 -a 11 u v 1 -a 12 q v 1 p 1 , T = T 1 1 -a 21 u v 1 -a 22 q v 1 p 1 , (43)
for the upper condensed surface y = 1:

p = p 2 1 + a 11 u v 2 + a 12 q v 2 p 2 , T = T 2 1 + a 21 u v 2 + a 22 q v 2 p 2 , (44)
Taking into account the conditions of the smallness of the pressure and temperature gradients, Eqs. [START_REF] Struchtrup | Evaporation boundary conditions for the r13 equations of rarefied gas dynamics[END_REF], linearized system of the Navier-Stokes equations (37) with the boundary conditions (43), (44) was solved analytically. The analytical dimensionless expressions read:

T (y) = Ay + B, (45)
A = a 12 a 11 CX P -CX T , B = 1 - a 12 2a 11 CX P -1 + η 2δ M CX T , (46)
p = 1 + 1 2 X P , (47)
u = 1 2a 11 + a 22 a 11 η 2δ CX P - a 12 a 11 η 2δ CX T , (48)
q = - a 12 a 11 η 2δ CX P + η 2δ CX T , (49)
J n = u, (50)
η = γ γ -1 1 P r , M = a 22 - a 2 12 a 11 , C = 1 + η δ M -1 . (51)
From the analysis of analytical solutions, Eqs. ( 45) -( 49), it is clear that the temperature profile is linear while the evaporation velocity, pressure and heat flux in the gap are constant. The form of the analytical solutions is analogous to the form of Eqs. [START_REF] Sharipov | Onsager-casimir reciprocity relations for open gaseous systems at arbitrary rarefaction. I General theory for single gas[END_REF], where the temperature, velocity and heat flux are split into two parts: one part related to the pressure gradient and the second one related to the temperature gradient. In Section 6 both solutions: the analytical expressions, Eqs. ( 45) -(49), and the numerical solutions of the linearized Smodel equation are compared in their separated (according to the corresponding gradients) form.

It is worth to note that the Onsager-Casimir relation leads here to the equality between the mass flux caused by the temperature drop u T and the thermal flux caused by the pressure drop q P [15]:

u T = q P , (52)
see analytical expressions, Eqs. ( 48), (49).

From the boundary conditions on the down plate, Eqs. (43), we can see that the temperature jump on the liquid-vapor interface depends on both evaporation velocity and heat flux. According to Eqs. (43) the temperature of the gas over an evaporating surface is lower than that of the condensed phase [START_REF] Sone | Kinetic theory of evaporation and condensation -hydrodynamic equation and slip boundary condition[END_REF]. Because of this temperature drop the "negative temperature gradient" phenomenon is observed when the rate of evaporation (or X P /X T ) is sufficiently large [START_REF] Pao | Application of kinetic theory to the problem of evaporation and condensation[END_REF], [START_REF] Sone | Kinetic theory of evaporation and condensation -hydrodynamic equation and slip boundary condition[END_REF]. However, when X P /X T is very small, the mass flow rate is directed from the colder to the hotter condensed phase. This situation was observed in several recent measurements, see [START_REF] Ward | Interfacial conditions during evaporation or condensation of water[END_REF], [START_REF] Badam | Experimental and theoretical invesigations on interfacial temerature jumps during evaporation[END_REF], [START_REF] Kazemi | Experimental and numerical study of the evaporation of water at low pressures[END_REF]. However in all these papers the positive temperature jump the vapor temperature is higher as that of the liquid phase, was observed while the linearized approach predicts the negative one.

Results and discussion

We present the results in following manner. First, the results obtained numerically from the linearized S-model kinetic equation and the analytical solution of the linearized Navier-Stokes equations subjected to the jumps boundary conditions are considered in order to show the limits of the applicability of this analytical solution in term of the rarefaction parameter. The cases of complete and non-complete evaporation-condensation are also presented and compared. Then, the linearized and non-linear numerical solutions of the S-model kinetic equation are compared to show the limits of the linearized approach application. Finally, the numerical solutions obtained from the S-model kinetic equation are compared with those giving by the BGK model in linearized and non-linear formulations.

Linearized approaches, complete evaporation

In this Section we will present the results obtained from the numerical solution of the linearized S-model kinetic equation and the analytical solution of the linearized Navier-Stokes equations completed by the pressure and temperature jump boundary conditions. Both solutions are determined only by two parameters: the gas rarefaction δ and the evaporation-condensation coefficient σ. Figure 2 shows the variations of the macroscopic vapor parameters: dimensionless number density and temperature in the gap between two condensed phases, obtained from the numerical solution of the linearized S-model equation. The case of the complete evaporation-condensation is considered, σ = 1, for four values of the rarefaction parameter: δ = 0.1, 1, 10 and 100. The macroscopic parameters are separated in two parts: one part due to the pressure gradient (subscript P ) and another one due to the temperature gradient (subcript T ), so we consider here these two gradients to be independent. From Figure 2 we can see that both dimensionless number density and temperature vary quasi linearly in the gap. The profiles of the number density, related to the pressure gradient, ν P , are almost coincide for the values of the rarefaction parameter δ = 1, 10 and 100, while the profiles of the number density variation due to the temperature gradient, ν T , are different for different values of the rarefaction parameter. The profiles of the dimensionless temperatures, τ P and τ T , show the "mirror" behavior in the gap: τ P increases while τ T decreases from down to upper condensed phases.

The evaporation velocity and the heat flux are constant in the gap, so their average over the gap values for different gas rarefaction are presented in Table 1. It is to note that the gas velocity generated by temperature gradient, u T , is oriented in the direction opposite to the evaporation flow generated by the pressure gradient, u P . In addition, the values of u T velocity is much lower than that of u P velocity, especially in the case of large value of rarefaction parameter, δ = 100. In this regime, the evaporation rate is engendered essentially by the pressure gradient. The evaporation velocities, u P and |u T |, decrease with rarefaction parameter increases. It can be seen from Table 1 that the heat flux generated by the temperature gradient, q T , is positive and directed from the hot bottom plate to the cold upper plate, while the heat flux generated by the pressure gradient, q P , has the opposite sign and is much smaller in absolute value. The values of the heat flux related to the pressure gradient, q P , are not provided in Table 1, because their values can be found from the Onsager-Casimir reciprocity relations q P = u T , see Eq. ( 52). Both heat fluxes, |q P | and q T , decrease with δ increasing, see Table 1.

S-model

By analyzing the results provided in Table 1 one can conclude that the difference between the analytical solution of the linearized Navier-Stokes equations subjected by the jump boundary conditions and the numerical solution of the S-model kinetic equation is very small. The difference between both solutions for δ = 10 is of the order of 0.1%. For larger value of the rarefaction parameter, δ = 1, this difference keeps surprisingly small, around 5%. Therefore, analytical expressions, Eqs. (45) -(49), can be used for estimation of the evaporation rate with a good accuracy up to δ = 10.

Linearized approaches, non-complete evaporation

Two approaches, the numerical solution of the S-model kinetic equation and the analytical solution of the Navier-Stokes equations subjected to the jump boundary conditions, are analyzed now for the case of non-complete evaporation. The cases of σ = 0.2, 0.5, 0.8 and 1 are considered. Figure 3 shows the influence of the non-complete evaporation-condensation on the macroscopic profiles of the number density and temperature, obtained from the numerical solution of the linearized S-model kinetic equation. The results for four values of the evaporation coefficient, σ = 0.2, 0.5, 0.8 and 1, are presented on Fig. 3 for rarefaction parameter δ = 10. It is worth to note that the profiles of the number density and temperature, caused by the temperature gradient, ν T and τ T , are independent from the evaporation coefficient, so they depend only on the temperature gradient between the condensed surfaces.

As it was mentioned above, the evaporation velocity and the heat flux are constant in the gap, so their averaged over the gap values, obtained from the numerical solution of the linearized S-model kinetic equation and the analytical solution of the linearized Navier-Stokes equations, Eqs. (45) -(49), are presented in Table 2 for the case of non-complete evaporation, σ = 0.5, and different values of the rarefaction parameter δ. Comparing the values of the evaporation velocity, provided in Tables 1 and2 one can conclude that the decrease in the evaporation coefficient decreases the absolute value of both evaporation velocities. The heat flux caused by temperature gradient, q T , remains practically the same with σ decreases, while the values of the heat flux caused by pressure gradient, q P = u T , see Table 2, decrease in approximately two times when the evaporation coefficient reduces from 1 to 0.5, see also Table 1. This behavior of the heat flux is in agreement with the temperature behaviors, presented on Fig. 3: temperature τ T , due to the temperature gradient, remains the same when the evaporation coefficient changes. It is worth to note that as for the case of the complete evaporation, the agreement between numerical solution of the linearized S-model kinetic equation and the analytical solution of the Navier-Stokes equation is very good: the difference between both solutions for δ = 10 is of the order of 1%.

S-model

Navier

Comparison of the BGK and S-model numerical solutions

As it was underlined in the Introduction, the S-model is chosen here because it provides the correct Prandtl number, equal to 2/3 for monatomic gases, contrary to the BGK model, for which the Prandtl number is equal to 1. In this Section both linear and non-linear BGK and S-model solutions are compared to see the influence of the applied models on the macroscopic profiles.

First, the numerical solutions of the linearized BGK and S-model equations are shown on Fig. 4 for the case of weak evaporation and condensation, T 1 /T 2 = 1.001, p 1 /p 2 = 1.01 and for the large value of the rarefaction parameter, δ = 250. It can be seen from Fig. 4 that both linearized solutions are in good agreement, except for the heat flux profiles, as expected. However, the heat flux profiles obtained by the linearized S-model equation and linearized Navier-Stokes equations coincide, because both models have the same Prandtl number.

Let us compare now the solutions of the non-linear BGK and S-model equations. Tables 3 and4 provide the numerical values of the dimensionless evaporation rate, Eq. ( 9), J n /(n 2 v 2 ) and energy flux, Eq. ( 10), J E /(p 2 v 2 ). Both quantities are obtained by integration of the kinetic equation, Eq. ( 5), and show the conservation of momentum and energy. Both quantities are independent from the type of the kinetic model and show the reliability of applied numerical methods to obtain the numerical solutions. The temperature and pressure ratios were chosen according to Ref. [START_REF] Aoki | Gas flows caused by evaporation and condensation on two parallel condensed phases and the negative temperature gradient: numerical analysis by using a nonlinear kinetic equation[END_REF]. As it is clear from Tables 3 and4 very good agreement was found between the presented here results (BGK and S-model) and the BGK results from Ref. [START_REF] Aoki | Gas flows caused by evaporation and condensation on two parallel condensed phases and the negative temperature gradient: numerical analysis by using a nonlinear kinetic equation[END_REF] for the dimensionless evaporation rate, J n and energy flux J E . The difference between three solutions is of the order of 1%.

The profiles of the macroscopic parameters, obtained from the numerical solution of the BGK and S-model equations are however slightly different, especially, as it was expected, the profiles of heat flux, which was also observed when comparing the linearized versions of both models. Figure 5 illustrates this difference in the macroscopic profiles of the evaporation velocity and heat flux in the gap between condensed phases. The comparison with the measurements are needed to conclude definitively about the advantages of one or other models.

Evaporation-condensation of Argon

In this Section we establish the applicability domain of the linearized S-model kinetic equation and of the analytical expressions (45) -(49). Contrarily to the linearized equation, to solve the non-linear S-model equation we need to fix the pressure and temperature ratios between two condensed surfaces. Usually the authors [START_REF] Aoki | Gas flows caused by evaporation and condensation on two parallel condensed phases and the negative temperature gradient: numerical analysis by using a nonlinear kinetic equation[END_REF] consider the variations of the saturation temperature and pressure independently, but here we consider more realistic case.

The S-model kinetic equation describes the behavior of the monatomic gases. Therefore, we fixed three temperature ratios, T 1 /T 2 = 1.001, 1.01 and 1.05 and estimated the corresponding saturation pressures for some monatomic gases, taking into account the Clayperon-Clausius equation [START_REF] Landau | Fluid Mechanics[END_REF]. The fixed temperature ratios of two condensed phases, T 1 /T 2 = 1.001, 1.01 and 1.05, lead to corresponding pressure ratios p 1 /p 2 = 1.01, 1.1, 1.55 for Krypton, starting with temperature 120K, to pressure ratios p 1 /p 2 = 1.01, 1.1, 1.57 for Argon, if we start with the temperature T 2 = 84K and to pressure ratios p 1 /p 2 = 1.004, 1.04, 1.2 for Helium. Some approximative tendency is so deduced: the increase of saturation temperature in 1% causes the increase of saturation pressure in approximately 10%. Therefore, the three sets of saturation parameters are chosen to solve numerically the non-linear S-model kinetic equation: T 1 /T 2 = 1.001, 1.01, 1.05 and p 1 /p 2 = 1.01, 1.1, 1.5 (see also two first lines in Table 5) to analyze the eventual limits of application for the linearized models.

δ T 1 /T 2 = 1.001 1.01 1.05 p 1 /p 2 = 1.01 1.1 1.5 Lin S
Lin NS Non-l S Lin S Lin NS Non-l S Lin S Lin NS Non-l S 1 0.00250 0.00245 0.00248 0.02599 0.02545 0.02469 0.15235 0.14842 0.12092 10 0.00239 0.00238 0.00237 0.02486 0.02480 0.02367 0.14573 0.14473 0.11684 100 0.00236 0.00236 0.00235 0.02458 0.02452 0.02344 0.14410 0.14310 0.11654 Table 5: Dimensionless evaporation rate Jn obtained by: the numerical solution of the linearized S-model kinetic equation (Lin S), the analytical expression Eq.( 50), (Lin NS), and the numerical solution of the non-linear S-model kinetic equation (Non-l S).

Table 5 provides the values of dimensionless evaporation rate, J n , obtained by: numerical solution of the linearized S-model kinetic equation, analytical expression, Eq.( 50), and numerical solution of the non-linear S-model kinetic equation. By comparing the results obtained from these three approaches it is clear that for two first pairs of saturation temperature/pressure values all the approaches gives very similar results. The difference between two linearized approaches does not exceed 0.2%. It is to note that the values of analytical solution of the linearized Navier-Stokes equations are provided in Table 5 also for δ = 1, even if theoretically this approach is not valid for δ < 10. But for given conditions the good agreement was found between all three approaches also for δ = 1. However, the difference between numerical solutions of the linearized and non-linear S-model kinetic equations becomes of the order of 4.6% for the temperature ratio 1.01 and this difference increases up to 20%, when the saturation temperature ratio increases up to 1.05. In the same time the agreement between two linearized approaches remains very good, of the order of 1% even for relatively large temperature ratio of 1.05. From the presented analysis, we can conclude that the linearized approaches work well for the real monatomic gases, i.e. taking into account their respective Clausius-Clayperon relations, only up to very small temperature ratio of the order of 1%, and therefore for the larger temperature ratios the non-linear approaches have to be applied.

The general tendency of the evaporation rate behavior can be derived by analyzing the results from Table 5: the evaporation rate decreases slightly when the rarefaction parameter increases for the fixed saturation temperature/pressure ratios. The evaporation rate increases more than in 10 times when the saturation temperature ratio increases in 10 times, this tendency of the proportionality remains up to temperature ratio equal to 1.05.

The numerical solutions of the linearized and non-linear S-model equations and the analytical solution of the linearized Navier-Stokes equations subjected to the pressure and temperature jumps boundary conditions, Eqs. ( 45) -( 49), are presented on Figs. 6789for two values of the rarefaction parameter δ = 10 and 100, for the case of the complete evaporation, σ = 1, and for three values of the saturation temperature/pressure ratios, presented in Table 5. The macroscopic profiles of the temperature, pressure, evaporation velocity and heat flux in the gap are provided.

It is seen from these Figures that all the temperature and pressure profiles for two temperature ratios T 1 /T 2 = 1.001 and 1.01 practically coincide for two values of the rarefaction parameter, δ = 10 and 100, except, the case, shown on Fig. 8(c), where the temperature profile becomes non-linear, in agreement with the findings of Ref. [START_REF] Aoki | Gas flows caused by evaporation and condensation on two parallel condensed phases and the negative temperature gradient: numerical analysis by using a nonlinear kinetic equation[END_REF]. In this case, the maximum of the temperature/pressure ratio is smaller than the rarefaction parameter and the linearized solutions do not coincide with the non-linear one, even if the deviations of the macroscopic parameters from their equilibrium state is small enough and normally the linearized theory should work.

With further saturation temperature increasing, up to T 1 /T 2 = 1.05, the temperature profiles, provided by the non-linear S-model equation, start to deviate from the linear ones for both value of δ, while the pressure profiles keep their shape and only absolute values of pressure change. The profiles of the evaporation velocity and of the heat flux differ from those provided by the linearized approaches even for small temperature ratio. For larger temperature ratios, T 1 /T 2 = 1.01 and 1.05, the behaviors of the linearized and nonlinear models differ not only quantitatively, but also qualitatively. Especially it is clearly seen near the liquid boundaries for heat flux distributions, Figs. 7 and 9 and for temperature distributions, Figs. 6 and8. The pressure profiles, obtained from all models and for all considered rarefaction levels, have the same qualitative behavior: fully coincide for T 1 /T 2 = 1.001 and have a difference in 2% for T 1 /T 2 = 1.05.

It is to note that the analytical solution of the linearized Navier-Stokes equations, even subjected with the temperature and pressure jump boundary conditions, cannot provide the correct profiles of the macroscopic parameters inside thin layer near condensed phases, the Knudsen layer, see the temperature, evaporation velocity and heat flux profiles. Even for very small pressure and temperature ratios, 1% and 0.1%, respectively, the profiles given by the analytical solutions are different, see Figs. 6789(a) and (b). However, the analytical results provide acceptable accuracy of the macroscopic parameters outside of the Knudsen layer up to temperature ratio T 1 /T 2 = 1.01 and can be used for rough estimations.

An interesting remark can be done concerning the limits of validity of the linearized approach: the consideration of the evaporation and condensation phenomena reduces considerably this limit. It was found in [START_REF] Chernyak | Non-linear effects in gases in couette problem[END_REF] that the heat flux between two parallel solid plates, obtained from the numerical solution of the non-linear S-model equation, still keeps its linear behavior up to temperature ratio equal to T 1 /T 2 ≈ 3, which is much larger compared to the limit found in present analysis.

It is worth to note that the "inversion of the temperature gradient", mentioned for the first time in Refs. [START_REF] Pao | Application of kinetic theory to the problem of evaporation and condensation[END_REF], [START_REF] Pao | Temperature and density jumps in the kinetic theory of gases and vapors[END_REF], is observed in the gap for all considered here cases. With the temperature gradient inversion the temperature near the hotter condensed surface becomes lower than that near the colder one. This phenomenon is obtained in numerical simulations, when the ratio between the pressure and temperature differences, X P /X T becomes larger than 4.772337 [START_REF] Sone | Kinetic theory of evaporation and condensation -hydrodynamic equation and slip boundary condition[END_REF]. Here this criterion is satisfied for all calculated cases, so the vapor temperature is larger near the colder condensed surface, see Figs. 6 and8.

Conclusion

The linearized and non-linear S-model equations are solved for different values of the rarefaction parameters and several saturation pressures and temperatures ratios. The case of the non-complete evaporation-condensation is also considered. The comparison between the BGK and S-model numerical solutions shows the difference only in the heat flux profiles while the other macroscopic parameters have not beed impacted by the model changing. The analytical solution of the linearized Navier-Stokes equations completed by the pressure and temperature jumps boundary conditions is obtained. The numerical solutions of the linearized and non-linear kinetic equations are compared together and with the analytical solutions of the continuum equations. It was found that the analytical solutions of the Navier-Stokes equations provide the excellent agreement with the numerical solutions of the linearized BGK and S-model kinetic equations up to rarefaction parameter equal to 10, in the case of complete and non-complete evaporation-condensation. The case of the real relations between saturation temperature/pressure for some monatomic gases are considered and it was found that the linearized approaches give satisfactory results only for the temperature deviation between condensed phases of the order of 1%, but for the larger saturation temperature ratio the non-linear approach has to be applied.

Acknowledgements

The research work of one of the authors (A. Polikarpov) was funded by RFBR according to the research project no. 18-31-00194. A.Polikarpov also thanks the Ministry of Education and Science of the Russian Federation Agreement no. 02.A03.21.0006 for travel grant. Irina Graur thanks Prof. Aldo Frezzotti for the helpful discussions.

Figure 1 :

 1 Figure 1: Sketch of 1D geometry

δT

  H /T C = 1.01 1.2 P H /P C = 1.02 2 BGK S-mod BGK [8] BGK S-mod BGK [8] 0.090 0.00421 0.00420 0.00421 0.22936 0.22925 0.22891 0.903 0.00427 0.00422 0.00429 0.22151 0.22103 0.22109 9.027 0.00454 0.00450 0.00457 0.22025 0.22006 0.22030

Figure 2 :Figure 3 :

 23 Figure 2: Numerical solution of the linearized S-model equation: the variations of macroscopic quantities, dimensionless number density and temperature, both separated in two parts: due to the pressure gradient (with superscript P , (a) and (c)) and due to the temperature gradient (with superscript T , (b) and (c)) in the gap between two condensed phases for different values of the rarefaction parameter δ and for the case of the complete evaporation, σ = 1.

Figure 4 :Figure 5 :

 45 Figure 4: Temperature (a) and pressure (b), obtained from the numerical solution of linearized BGK and S-model equations, and non-linear S-model equation, T 1 /T 2 = 1.001, p 1 /p 2 = 1.01, δ = 250.

T 1 /T 2 = 1 .05, p 1 /p 2 = 1 5 Figure 6 :

 12112156 Figure 6: Comparison of three results for δ = 10 and complete evaporation, σ = 1: the numerical solution of the linearized S-model kinetic equation, the analytical solution of the linearized Navier-Stokes equations with the temperature and pressure jumps boundary conditions and the numerical solution of the non-linear S-model kinetic equation. The temperature and pressure profiles in the gap are compared for three sets of saturation temperature and pressure ratios: (a), (b) T 1 /T 2 = 1.001, p 1 /p 2 = 1.01, (c), (d) T 1 /T 2 = 1.01, p 1 /p 2 = 1.1; (e), (f) T 1 /T 2 = 1.05, p 1 /p 2 = 1.5.

Table 1 :

 1 Macroscopic parameters obtained from the numerical solution of the linearized Smodel kinetic equation and the analytical solution of the Navier-Stokes equations with the jump boundary conditions for different values of the rarefaction parameter δ and for complete evaporation, σ = 1.

	Navier-Stokes

Table 2 :

 2 Macroscopic parameters obtained from the numerical solution of the linearized Smodel kinetic equation and the analytical solution of the Navier-Stokes equations with the jump boundary conditions for different values of the rarefaction parameter δ and for noncomplete evaporation, σ = 0.5.

	-Stokes

Table 3 :

 3 Dimensionless evaporation rate Jn

		T H /T C = 1.01		1.2
	δ	P H /P C = 1.02		2
	BGK	S-mod BGK [8]	BGK	S-mod BGK [8]
	0.090 0.01395 0.01397 0.01391 0.66267 0.66293 0.66141
	0.903 0.01308 0.01320 0.01296 0.63507 0.63633 0.63287
	9.027 0.01206 0.01217 0.01201 0.62106 0.62161 0.62105

Table 4 :

 4 Dimensionless energy flux J E


Mustafa Hadj-Nacer
email: mhadjnacer@unr.edu

Dilesh Maharjan
email: dileshz@gmail.com

Minh-Tuan Ho
email: minh-tuan.ho@strath.ac.uk

Professor Stefan K Stefanov
email: stefanov@imbm.bas.bg

Irina Martin Graur

Professor, Miles Greiner
email: greiner@unr.edu

Continuum and kinetic simulations of Heat Transfer trough Rarefied Gas in Annular and Planar Geometries in the Slip Regime

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

Introduction

Used nuclear fuel assemblies consists primarily of square arrays of zircaloy cladding tubes containing highly-radioactive solids and gases [START_REF]0184, december[END_REF]. They are initially stored underwater while their heat generation and radioactivity decrease [START_REF] Saling | Radioactive Waste Management[END_REF]. After sufficient time individual assemblies are placed in square openings within helium-fill canisters. Before a canister is sealed, vacuum drying is used to remove remaining moisture [START_REF] Colmont | Improved thermal modeling of snf shipping cask drying process using analytical and statistical approaches[END_REF]. During periods of time of this process the helium within the canister's narrow spaces is in the rarefied-gas slip regime. This causes the cladding temperature to increase compared to continuum gas conditions. The cladding temperature must be accurately predicted to assure its temperature does not exceed certain important-to-safety limits [START_REF] Commission | Cladding considerations for the transportation and storage of spent fuel[END_REF].

The ultimate goal of this research is to develop and experimentally-benchmark computational methods to predict surface temperatures of complex void spaces containing slip-regime gases. For moderately low surface temperature differences and for surfaces with radii of curvature that are larger than the mean free path (mean distance traveled by molecules between two successive collisions), transport across slip-regime gases may be approximated using continuum Fourier conduction in the gas bulk, and temperature-jumps boundary condition at its interfaces with heated surfaces [START_REF] Kennard | Kitetic theory of gases[END_REF].

In the current paper, conduction across annular gaps filled with slip regime helium calculated from two kinetic gas methods [START_REF] Graur | Comparison of different kinetic models for the heat transfer problem[END_REF][START_REF] Graur | Simulation of the transient heat transfer between two coaxial cylinders[END_REF] and Direct Simulation Monte Carlo [START_REF] Bird | Molecular Gas Dynamics and the Direct Simulation of Gas Flows[END_REF]) are compared to calculations from continuum models with temperature-jump boundary conditions. Annular gaps with a range of radius ratios are simplified models for regions within used fuel canister and other engineering systems. In this paper, comparisons are made between the continuum and kinetic models for a range of radius ratios, surface temperature ratios, gas rarefaction parameters and surface accommodation coefficients. These comparisons are used to determine the conditions where the continuum models with temperature-jump boundary conditions may be used with confidence for the simple annular gap geometry. Future work will compare continuum and kinetic models for calculating heat transfer within an enclosed array of heated rods with experiment data in order to determine when the continuum models may be used in complex configurations.

Problem formulation

The problem of conduction heat transfer between parallel plates and coaxial cylinders separated a gas at rest is considered here (see Fig. 1). The radii and temperatures of the inner and outer cylinders are (R 1 , T w 1 ) and (R 2 , T w 2 ), respectively, with R 1 < R 2 and T w 1 > T w 2 . For simplicity, the same notation is used for parallel plates: the r axis is normal to both plate's surfaces. The bottom and upper plates have the locations and temperatures (R 1 , T w 1 ) and (R 2 , T w 2 ), respectively, with R 2 = -R 1 = H /2. The cylinders and plates are assumed to have infinite length in the direction perpendicular to the figure.

Both problems are similar, however, the coaxial cylinder problem is governed by three physical parameters (temperature and radius ratios, and rarefaction parameter; see definitions below). The parallel plate problem is defined by only two parameters (temperature ratio and rarefaction parameter). It is convenient to introduce the following parameters:

• The temperature ratio between the cylinders or plates:

T = T w 1 /T w 2 ,
• The aspect ratio for the cylinders and plates: R

= |R 1 /R 2 |, • The rarefaction parameter δ = R 0 , where = µ 0 υ 0 p 0 , υ 0 = 2k B T 0 m (1)
In these expressions, is the equivalent mean free path at reference pressure p 0 , m is the molecular mass of the gas, k B is the Boltzmann constant, and υ 0 is the most probable molecular velocity at reference temperature T 0 .

In expression (1), R 0 and T 0 are the reference space dimension and temperature, respectively, and µ 0 is the gas viscosity calculated at the reference temperature; µ 0 = µ(T 0 ).

It is convenient to take the distance between the cylinders (or plates) as the reference space dimension R 0 = R 2 -R 1 . The temperature of the external cylinder (or upper plate) is used as the reference temperature, T 0 = T w 2 . The problem of heat transfer considered here is one dimensional, therefore, the heat flux has only one component in the direction normal to the surfaces, denoted by q .

For practical reasons, the following dimensionless variables are introduced

r = r R 0 , p = p p 0 , T = T T 0 and q = q p 0 υ 0 . ( 2
)
The influence of the gas-surface interaction is taken into account by the thermal accommodation coefficient, denoted as

α = T i -T r T i -T w , (3)
where T i and T r are the temperature of the incident and reflected molecules, respectively.

Continuum approach

The analytical expressions for heat flux and temperature can be obtained from an energy balance. The hypothesis of zero macroscopic velocity and constant pressure between the plates or cylinders are used to derive the analytical expressions. The Fourier law can be applied to calculate the heat flux

q = -κ dT dr . (4)
For monatomic gases, the gas thermal conductivity κ is related to the gas viscosity µ as follows [START_REF] Bird | Molecular Gas Dynamics and the Direct Simulation of Gas Flows[END_REF]

κ = 15 4 k B m µ . (5)
In order to define the dependence of viscosity on temperature, the molecular interaction potential must be specified. In this work, the inverse power law potential [START_REF] Bird | Molecular Gas Dynamics and the Direct Simulation of Gas Flows[END_REF] is employed, which leads to a power law temperature dependence for the viscosity,

µ = µ 0 T T 0 ω , ( 6
)
where ω is the viscosity index, which is equal to 0.5 for the Hard-Sphere (HS) model and 1 for the Maxwell model. Its value varies as a function of the gas nature in the case of Variable Hard-Sphere (VHS) model, see [START_REF] Bird | Molecular Gas Dynamics and the Direct Simulation of Gas Flows[END_REF]. In this paper, the Hard-Sphere model is retained for all calculations, so ω = 0.5 for all gases.

In the continuum regime, the temperature continuity condition may be assumed on the surfaces. However, in the slip regime, the temperature-jump conditions [START_REF] Cercignani | Mathematical methods in kinetic theory[END_REF] must be used as boundary conditions,

T g = T w + ξ T dT dr w , (7)
where T g is the gas temperature near the wall and ξ T is the temperature-jump coefficient [START_REF] Kogan | Rarefied gas dynamics[END_REF]. The dimensionless form of the temperature-jump boundary condition at the cylinders and plates walls can be written as

T =      T w1 + ξ T 1 δ T ω+1/2 w1 dT dr , r = R 1 T w2 -ξ T 2 δ T ω+1/2 w2 dT dr , r = R 2 (8)
Here ξ T 1 , ξ T 2 are the temperature-jump coefficients of the hotter and colder surfaces, respectively. The assumption of the constant pressure between the cylinders and plates is used to obtain the previous expressions. The coefficient ξ T in Eq. ( 7) and [START_REF] Graur | Simulation of the transient heat transfer between two coaxial cylinders[END_REF] depends on the gas nature and surface state through the thermal accommodation coefficient. Lin and Willis, 1972 [START_REF] Lin | Kinetic theory analysis of temperature jump in a polyatomic gas[END_REF] proposed the following expression for polyatomic gases,

ξ T = 2 -α α + 0.17 √ π Pr γ γ + 1 , (9)
obtained by applying a variational method to the Morse equation [START_REF] Morse | Kinetic model for gases with internal degrees of freedom[END_REF] and to the Holway model [START_REF] Holway | New statistical models in kinetic theory: methods of construction[END_REF]. In Eq. 9 γ is the gas specific heat ratio and Pr is the Prandtl number. For monatomic gases, this expression is reduced to the one proposed by Welander [START_REF] Welander | On the temperature jump in a rarefied gas[END_REF]. For the case of complete accommodation (α = 1) and a monatomic gas (γ=5/3 and Pr=2/3), the value of this coefficient is ξ T ∼ 1.95. In this paper, Eq (4) subjected to the boundary condition ( 8) is solved numerically using an iterative method for parallel plates and coaxial cylinders having different temperature and aspect ratios at different values of the rarefaction parameter δ in the slip regime, and various values of the thermal accommodation coefficient α.

Analytical solution

To obtain the completely explicit expression of the temperature distribution between the cylinder and plates, a linearization of the temperature is carried out and the terms of the order of ε 2 are neglected, where ε = T w -T g T w (Eq. 7), so the temperature profile between the cylinders and plates is obtained analytically with accuracy of ε 2 .

Parallel plates

The following dimensionless analytical expressions are obtained for the temperature profile

T (r) = 1 2 T ω+1 g1 + T ω+1 g2 +C (ω + 1) (R 1 + R 2 -2r) 1 ω+1 (10)
and heat flux profile

q (r) = 15 8δ C, (11)
where

C = T ω+1 w1 -T ω+1 w2 / (ω + 1) R 2 -R 1 + ξ T 1 δ T ω+1/2 w1 + ξ T 2 δ T ω+1/2 w2 . ( 12
)
The gas temperatures near the walls are

T ω+1 g1 = T ω+1 w1 [1 -η 1 (ω + 1)] , T ω+1 g2 = T ω+1 w2 [1 -η 2 (ω + 1)] , (13)
where

η 1 = ξ T 1 δ C 1 T 1/2 w1 , η 2 = - ξ T 2 δ C 1 T 1/2 w2 . (14)
Here, the dimensionless coordinates R 1 and R 2 for parallel plates are: R 1 = -0.5 and R 2 = 0.5.

Coaxial cylinders

By implementing the same technique, the temperature and heat flux profiles for coaxial cylinders were obtained, as in Ref. [START_REF] Graur | Rarefied gas flow through a long rectangular channel of variable cross section[END_REF]. The dimensionless temperature profile with an accuracy of ε 2 is

T (r) = 1 2 T ω+1 g1 + T ω+1 g2 +C (ω + 1) ln R 1 r + ln R 2 r 1 ω+1 , (15)
and the dimensionless heat flux profile is

q (r) = 15 8δ C r , (16)
where

C = T ω+1 w1 -T ω+1 w2 / (ω + 1) ln R 2 R 1 + ξ T 1 δ T ω+1/2 w1 R 1 + ξ T 2 δ T ω+1/2 w2 R 2 . (17)
In this paper, two aspect ratios are considered, R = 0.5 and R = 0.1. The corresponding dimensionless coordinate [R 1 ,

R 2 ] for each case are: [START_REF]0184, december[END_REF][START_REF] Saling | Radioactive Waste Management[END_REF] and [0.11, 1.11], respectively. The gas temperature in the vicinity of the walls is calculated with the same expressions (13) as for the two parallel plates, with the parameters ε 1 and ε 2 calculated as

η 1 = ξ T 1 δ C R 1 1 T 1/2 w1 , η 2 = - ξ T 2 δ C R 2 1 T 1/2 w2 . ( 18
)
The analytical model is more convenient to calculate the temperature profile between the cylinders and plates than the numerical model as it needs no iteration.

Kinetic approach

Two kinetic models are used in this paper to simulate heat transfer between parallel plates and coaxial cylinders in the slip regime.

S-model kinetic equation

In the steady-state case the S-model kinetic equation [START_REF] Shakhov | Generalization of the Krook kinetic relaxation equation[END_REF] reads

v ∂ f ∂r = ν f S -f . ( 19
)
where, f (r , v) is the one particle molecular velocity distribution function, v is the molecular velocity vector, r is the position vector, and ν is the molecular collision frequency. It should be underlined that Eq.( 19) is valid if a Cartesian frame is used in the velocity space. In the case of the two coaxial cylinders, an axisymmetric problem, the transport part of Eq.( 19) can be written in form used in Refs. [START_REF] Sone | Evaporation of a rarefied gas from a cylindrical condensed phase into a vacuum[END_REF], [START_REF] Graur | Simulation of the transient heat transfer between two coaxial cylinders[END_REF].

The equilibrium distribution function f S is

f S = f M 1 + 2mVq 15n (k B T ) 2 mV 2 2k B T - 5 2 , f M (n , T ) = n m 2πk B T 3/2 exp - m(v -u ) 2 2k B T , (20)
where f M is the local Maxwellian distribution function, u is the bulk velocity vector, V = vu is the peculiar velocity vector, and q is the heat flux vector.

In the frame of this model, the molecular collision frequency is assumed to be independent of the molecular velocities and may be found as follows [START_REF] Shakhov | Generalization of the Krook kinetic relaxation equation[END_REF]

ν = p µ . ( 21
)
The macroscopic parameters are calculated from the solution of the molecular distribution function f obtained by solving eq. ( 19), as

T (r ) = m 3n k B V 2 f (r , v)dv, q (r ) = m 2 VV 2 f (r , v)dv. (22)
Considering the axial symmetry of the problem, the S-model kinetic equation in the completely conservative form is found in Refs. [START_REF] Graur | Simulation of the transient heat transfer between two coaxial cylinders[END_REF], [START_REF] Larina | A numerical method for calculationg axisymmetric rarefied gas flows[END_REF], and [START_REF] Shakhov | Numerical study of the generalized cylindrical couette flow of rarefied gas[END_REF]. The generalization of a kinetic model for the case of the polyatomic gases is presented in [START_REF] Hsu | Kinetic theory of parallel plate heat transfer in a polyatomic gas[END_REF].

The Discrete Velocity Method (DVM) is used to divide the continuum-molecular velocity space c p in the system of kinetic equations (see Ref. [START_REF] Graur | Simulation of the transient heat transfer between two coaxial cylinders[END_REF]) into a discrete velocity set c p k . The system of kinetic equations with discrete velocity set c p k is then discretized in space by the Finite Difference Method (FDM). The spatial derivatives are approximated by the second order-of-accuracy upwind-type numerical scheme. The Gauss-Hermite quadrature formulas are used to evaluate the integrals for calculating the macroscopic parameters (Eq. 22).

In the physical space, the reduced distribution functions depend only on the spatial variable r, which is either the distance from the cylinder's axis or the distance from midway between the plates. For coaxial cylinders, the distance between the walls is split into N r intervals, N r = 800 for δ ≤ 10 and N r = 6400 for δ > 10. For parallel plates, N r = 4000 for all δ. In the velocity space, the distribution functions depend on two variables; the magnitude and orientation of the molecular velocity vector c p k and ϕ, respectively. The velocity vector magnitude c p k is distributed according to the Gaussian quadrature rule, which is characterized by Gaussian abscissas N c p and their corresponding weights. The number of implemented points N c p depends on gas rarefaction; N c p = 25 for δ < 0.5 and N c p = 12 for δ ≤ 0.5. The range of molecular velocity orientations (0 ≤ ϕ ≤ π) is divided into N ϕ equal intervals. For coaxial cylinders, N ϕ = 6400 for δ ≤ 10 and N ϕ = 100 for larger δ. For parallel plates, N ϕ = 100 for all δ. Moreover, this range of ϕ is divided into two sub-domains, according to the sign of the molecular velocity components. This set of the numerical grid parameters guarantees the accuracy for energy conservation law of the order of 0.1%.

Direct Simulation Monte Carlo (DSMC)

The traditional DSMC technique proposed by Bird [START_REF] Bird | Molecular Gas Dynamics and the Direct Simulation of Gas Flows[END_REF] is employed in this paper for the calculation of heat transfer between parallel plates and coaxial cylinders. The DSMC technique enables gas flows to be modeled on a molecular level by simulating the motion of individual particles according to their physical properties. This technique can be viewed as a Monte Carlo method for solving the time-dependent non-linear Boltzmann equation. Within each time step t of the simulation, this method combines deterministic aspects for modeling particle motions with statistical aspects for computing collisions between particles. The collision technique that is used is the NoTime Counter scheme suggested by Bird, with a slight modification in the calculation of maximum number of collision in a cell, as described in Ref. [START_REF] Stefanov | Monte carlo simulation and navier-stokes finite differnce solution of rarefied gas flow problems[END_REF].

The procedure for axis-symmetric flows is employed as given in [START_REF] Bird | Molecular Gas Dynamics and the Direct Simulation of Gas Flows[END_REF]. For each particle, only the radial coordinate r and the three velocity components v x , v r , and v θ are stored. Then, the molecules are moved according to their velocities and acquire a new radial coordinate r + given by

r + = (r + v r ∆t) 2 + (v θ ∆t) 2 . ( 23
)
Following this procedure, the new velocity components, v + x ,v + r , and v + z , are

v + x = v x , v + r = v r (r + v r ∆t) -v 2 θ ∆t r +
, and

v + z = v r v θ ∆t -v θ (r + v r ∆t) r + . ( 24
)
The domain is discretized in the radial direction using N r computational cells (N r ranges from 400 to 4000) and the particles are initially distributed in such a way that the minimum number of particles in each cell of the domain is N p = 50. Since no weighting factors are used, the region where the number of particles per cell is the smallest corresponds to the area near the inner cylinder for geometrical reasons. The time step ∆t is chosen to be less or equal to one third the cell traversal time, ∆r/ 2kT 0 /m. The macroscopic properties are obtained for the steady state solution by time averaging over N T = 10000 time steps. The heat flux and other macroscopic quantities are volume based calculated by averaging the microscopic values of the particles at a given cell, q r = v r v 2 . For these grid and time parameters, the total computational error (systematic and statistical) is less than 1% in all the computed cases. Here, the available analyses of the splitting scheme and domain discretization [START_REF] Gallis | Convergence behavior of a new dsmc algorithm[END_REF][START_REF] Stefanov | On DSMC calculations of rarefied gas flows with small number of particles in cells[END_REF] suggest a first order systematic approximation error of (O(∆t, ∆r)), with respect to the Boltzmann equation. The statistical error [START_REF] Hadjiconstantinou | Statistical error in particle simulations of hydrodynamic phenomena[END_REF] is estimated as O(1/ N T N p ).

It should be noted here that the computer code described above was employed for two coaxial cylinders with R = 0.1 and 0.5, and two parallel plates with R = 0.999.

Results and Discussions

In this section, results of steady-state heat transfer, between parallel plates and coaxial cylinders in the slip regime, obtained from kinetic and continuum approaches, are presented. The results are showed in graphical form for the macroscopic quantities (temperature, heat flux, and pressure) in term of all the quantities involved in the problem. For all three geometries, parallel plates (R = 1) and coaxial cylinders (R = 0.5 and 0.1), two temperature ratios T = 1.1 and 2 are examined.

The value of the rarefaction parameter δ is varied from 500 to 3, which covers the continuum, slip, and near transitional regimes [START_REF] Schaaf | Flow of rarefied gases[END_REF]. In this work, we wish to determine the conditions where the difference between the continuum and kinetic results of the same order as the difference between the S-model and DSMC results.

Dimensionless temperature profiles

In order to assess the limit of applicability of the continuum models, first, a comparison of the dimensionless temperature profiles obtained from the S-model and DSMC technique is performed. Figure 2 shows the dimensionless temperature distribution along the r-axis, obtained from the kinetic (S-model and DSMC) approaches and the continuum (numerical and analytical) models for the temperature and radius ratios considered in this work, in the case of full accommodation of the molecules at the walls (α 1 = α 2 = 1) and for different values of rarefaction parameter δ = 100, 50, 10 and 3.

Comparison between the profiles of dimensionless temperature obtained from the kinetic model shows that the S-model and DSMC technique are in good agreement for all T , R , and δ employed in this paper. Both models exhibit the same profile even in the Knudsen layers (local non-equilibrium region with a thickness ∼ O( ) from the wall).

As expected, for T = 1.1, the profiles of dimensionless temperature from the continuum numerical and analytical [START_REF] Cercignani | Mathematical methods in kinetic theory[END_REF][START_REF] Welander | On the temperature jump in a rarefied gas[END_REF] models exhibit a nearly linear trend for parallel plates and a nearly logarithmic trend for coaxial cylinders. This is due to the relatively small thermal conductivity variation across the gaps. For the smaller temperature ratio (T = 1.1, left plots), the numerical and analytical temperature profiles are very similar for all values of δ, but the difference increases as R decreases.

For the larger temperature ratio (T = 2, right plots), the difference between the two continuum models is larger and increases as δ decreases. In all cases, the temperatures predicted by the numerical model are systematically larger than that from the analytical model, but within 2%. This difference comes from the linearization of the temperature used in the analytical model, see Sections. 3.1.1, 3.1.2. The comparison between dimensionless temperature profiles obtained from the continuum and kinetic models is divided into two categories according to the temperature ratio T . • For T =1.1 and all R , there is a good agreement between the two types of models for δ =100 and 50. For δ=10, the agreement is good in the core of the gap (away from the walls), except for the case R =0.1. Deviations are observed close to the walls. For smaller δ, there is disagreement even within the core of the gap for all aspect ratios R . • For T =2, the agreement between the continuum and kinetic models for δ =100 and 50 is not as close as it is for T =1.1.

However, for smaller δ ≤ 10, the disagreement increases as δ decreases. Figure 2 shows that the temperature-jump at the walls increases as δ decreases. The temperature-jump at the hotter wall is larger than the corresponding one at the cooler wall and it increases as R decreases. Furthermore, as T increases, the temperature-jump at both walls increases. However, the increase of T has a larger effect on the jump at the inner wall than the outer wall. This is because the gas becomes more rarefied as its temperature increases.

Overall, the temperature-jumps at the walls predicted by the continuum models increase as δ and R decrease, and T increases and they are larger than those predicted by the kinetic models. (q-q S-model )/q δ T =2 Fig. 5. Percent deference of the dimensionless heat flux q between S-model and DSMC (black), and Numerical-continuum (red) models as function of the rarefaction parameter δ obtained for all combinations of T and R and different values of thermal accommodation coefficient α.

Dimensionless heat flux

Figure 3 shows the dimensionless heat flux q, which is the quantity of practical interest in this paper, plotted as a function of rarefaction parameter δ for all combinations of T and R considered in this paper and for α=1, 0.6 and 0.2 with

α = α 1 = α 2 .
Results from kinetic (S-model and DSMC) and continuum (numerical and analytical) approaches are included.

It is clear from this figure that the value of q increases as δ decreases and α and T increase. The values of q obtained from both continuum models (numerical and analytical) are in a good agreement for case T =1.1, regardless R and α. For T =2, the discrepancy between the two models increases as δ, α, and R decrease. The value of q for the continuum numerical model is systematically larger than the continuum-analytical model. At T =1.1, the continuum and kinetic results are in good agreement for all values of δ, except for cases R =0.1 and α=1. The reason of this discrepancy is not clear and may be investigated further. For T =2, the agreement is good for R =0.1 and α=1, but is less good as α decreases and R increases.

The disagreement between the continuum and kinetic models, especially at large T , may be explained from the profiles of dimensionless pressure along the r-axis. Figure 4 shows these profiles between the cylinders and plates obtained from the S-model kinetic equation in the case of full accommodation of the molecules at the walls. The results are plotted for T =1.1 and 2, R =1, 0.5, and 0.1, and δ=100, 50, 10, and 3. The r-axis was scaled to be between 1 and 2 as described in the caption for convenience of comparison. From this figure, it can be seen that, in all cases, the pressure varies along the r-axis between the cylinders and plates and this variation is larger as δ and R decrease, and T increases. As mentioned earlier, the continuum models assume that the pressure is constant along the r-axis. Figure 4 shows that the pressure variation is very small for T =1.1, which explains the good agreement obtained between the continuum and kinetic approaches. However, for T =2, the pressure variation is significant for δ ≤ 10, therefore, the assumption of the constant pressure is not valid. This may explain the disagreement obtained in these cases. Figure 5 gives the percentage difference of dimensionless heat flux q between DSMC and S-model (shown in black symbols) and the continuum-numerical and S-model (shown in red symbols) for both temperature and aspect ratios, and for three values of the thermal accommodation coefficient considered in this work. It can be seen from this figure that the heat fluxes obtained from the DSMC and S-model are within 3% of each other in all configurations. This difference decreases as δ increases and the gas approaches the continuum regime. It is also clear from Fig. 5 that, for the case T =1.1, the percentage difference between the continuum-numerical and S-model is smaller than 4% for all values of R , α, and δ, which is of the order of the difference between the DSMC and S-model approaches. However, in case T =2, the percentage difference between the continuum-numerical and S-models is significantly larger than the corresponding difference between DSMC and S-model and of the order of 50% at small δ. This quantitative comparison between the continuum and kinetic approaches clearly demonstrates that the continuum models are appropriate for simulating heat transfer for small temperature ratio, however, for large ratio their ability is largely reduced even at large value of δ ≥ 50.

In general, it can be concluded that the continuum models are accurate for δ ≥ 10 and T =1.1, and for δ ≥ 50 and T =2, regardless the aspect ratio.

Conclusion

The steady state heat transfer through rarefied gas confined between coaxial cylinders and parallel plates is studied using continuum and kinetic approaches in the slip regime. The calculations were conducted for different temperature and radius ratios. Both walls were assumed to have the same value of thermal accommodation coefficient α.

The dimensionless temperature and heat flux were obtained from S-model and DSMC were compared to continuum numerical and analytical solutions of the Fourier equation subjected to the Lin & Willis temperature-jump boundary condition at the walls.

The two continuum models were in good agreement for all values of δ and α, with maximum difference less than 2%, observed for case T =2 and R =0.1. The comparison between the continuum and kinetic approaches showed that the continuum models are valid in the slip regime for a temperature ratio of T =1.1, however, for T =2, they are only valid for δ ≥ 50.

It should be noticed that here the temperature boundary conditions are employed in the paper. However sometimes in practice it is important to be able to employ constant heat flux boundary conditions so that the heat exchange between the system and the surroundings can be controlled. The examples of the implementation of the heat flux boundary conditions are given in [START_REF] Akhlaghi | A newiterative wall heat flux specifying technique in dsmc for heating/cooling simulations of mems/nems[END_REF], [START_REF] Meng | Numerical simulation of rarefied gas flows with specified heat flux boundary conditions[END_REF] for the DSMC and kinetic approaches, respectively.

Fig. 1 .

 1 Fig. 1. Cross-section of (a) two coaxial cylinders and (b) two parallel plates configurations: dimensions (r, θ) in physical space, dimensions (υ r , υ θ ) (or (υ p , ϕ)) in molecular velocity space.

Fig. 2 .

 2 Fig. 2. Dimensionless temperature profiles between plates and cylinders for all combination of R and T and different values of rarefaction parameters δ = 100, 50, 10 and 3 in case α = 1.

Fig. 3 .

 3 Fig. 3. Dimensionless heat flux as function of the rarefaction parameter δ obtained for all combination of R and T and different values of thermal accommodation coefficient α.

Fig. 4 .

 4 Fig. 4. Dimensionless pressure profiles between plates and cylinders for all combination of R and T and different values of rarefaction parameters δ = 100, 50, 10 and 3 in case α = 1. Black lines are for parallel plates (R = 1), blue and green lines are for coaxial cylinders with R = 0.5 and R = 0.1, respectively. l = r for R = 0.5 and l = 1 -R 1 + r for R = 1 and 0.1.

Acknowledgments

This work was supported by the US Department of Energy Office of Nuclear Energy University Program. This work was granted access to the HPC resources of Aix-Marseille Université, under the project number 13b015, financed by the project Equip@Meso (ANR-10-EQPX-29-01) of the program "Investissements d'Avenir" supervised by the Agence Nationale pour la Recherche.