Fast heat transfer calculations in supercritical fluids versus hydrodynamic approach

0.4 MB·

by  Vadim Nikolayev

Fast heat transfer calculations in supercritical fluids versus hydrodynamic approach
V. S. Nikolayev,1,∗ A. Dejoan,2,† Y. Garrabos,2 and D. Beysens1
1DSM-DRFMC-Service des Basses Temp´eratures/ESEME, CEA-Grenoble, France‡
2CNRS-ESEME, Institut de Chimie de la Mati`ere Condens´ee de Bordeaux,
6
87, Avenue du Dr. A. Schweitzer, 33608 Pessac Cedex, France
1
(Dated: January 26, 2016)
0
2
Thisstudyinvestigatestheheattransferinasimplepurefluidwhosetemperatureisslightlyabove

n its critical temperature. We propose a efficient numerical method to predict the heat transfer in
a such fluids when the gravity can be neglected. The method, based on a simplified thermodynamic
J approach,iscomparedwithdirectnumericalsimulationsoftheNavier-Stokesandenergyequations

5 performed for CO2 and SF6. A realistic equation of state is used to describe both fluids. The pro-
2 posedmethodagrees with thefullhydrodynamicsolution andprovidesahugegain incomputation
time. The connection between the purely thermodynamic and hydrodynamic descriptions is also

] discussed.
n
y
d I. INTRODUCTION derstood, the calculations required to represent realistic
-
experimentalconditionsaredifficultbecausetheinherent
u
l In fluids near their liquid-gas critical point, the char- nonlinear dynamical behavior of such fluids is compli-
f
. acteristic size of the density fluctuations becomes larger cated by the highly non-linear equations of state (EOS)
s
c than the characteristic size of the molecular structure. used to describe real near-critical fluids. Two computa-
i Consequently, the fluid behavior is ruled by fluctuations tionalapproacheshavebeensuggestedforconfinedfluids
s
y and not by its particular molecular structure. This im- in absence of convection. In one of them, which we will
h plies that most fluids behave similarly near the critical refer as the “thermodynamic approach” [5], the Piston
p point (like the 3D Ising model). This universality makes effect is taken into account by a supplementary term g,
[ the study of near critical fluids very appealing. Due to introduced in the heat conduction equation as follows,

1 their very low thermal diffusivity and to their very large ∂T 1
v thermalexpansionandcompressibility,the studyofheat = (k T)+g(T), (1)
∂t ρc ∇· ∇
6 transfer in such fluids is particularly challenging. When p
4 such a fluid is confined in a heated cavity, a very thin with
6
hot boundary layer develops and induces a fast expan-
6 cv ∂T ∂p
sion that compresses the rest of the fluid. The result- g(T)= 1 , (2)
0 − c ∂p ∂t
. ing pressure wavesspread at the sound velocity (i.e very (cid:18) p(cid:19)(cid:18) (cid:19)ρ
1
rapidly) and adiabaticallycompress the bulk of the fluid whereT isthelocalfluidtemperature,c (c ),thespecific
0 v p
6 which is therefore homogeneously heated. After several heat at constant volume (pressure) per unit mass, ρ the
1 sound wave periods the pressure is already equilibrated localdensityandk the thermalconductivityofthe fluid.
: andcanbe assumedto be nearlyhomogeneousalongthe One can note that the term g(T) is only relevant near
v
cavity. During the initial stage of heating, this process the critical point where c c .
i p v
X of energy transfer, called “Piston effect” [1, 2], is much Thefluidmotionsareneg≫lectedandthepressurep,as-
r moreefficientthantheusualdiffusionscenario. Indeed,if sumed to be homogeneous, is only a function of time t.
a the Pistoneffectwereabsent,the bulk ofthe fluid would The pressurep isdeterminedfromthe fluidmassconser-
remain at the initial temperature. vation and computed via the nonlinear expression [5]
Inthe industrialdomain,thePistoneffectcanbeused
to transfer heat much faster than by conduction. This ∂p = v(∂ρ/∂T)p ∂T/∂tdv (3)
featurecanbereadilyappliedtothedevelopmentofheat ∂t − ρχ dv
R v T
exchangers under microgravity conditions [3, 4] where
whereχ =ρ−1(∂ρ/∂p) isRtheisothermalcompressibil-
heat transfer by naturalconvectionis obviously not pos- T T
ity and v the volume of the fluid sample. The resolution
sible.
of (3) requires an iterative procedure for each time step.
WhilethephysicaloriginofthePistoneffectiswellun-
Thisconsistsincalculatingthe temperaturein the whole
fluid volume using Eqs. (1-2) for some trial value of p
(and thus ∂p/∂t) and determining the other thermody-
∗email:[email protected] namic parameters (ρ, χT,...) by an EOS
†Present address: Aeronautics Department, Imperial College,
Λ(p,ρ,T)=0. (4)
Prince Consort Road, South Kensington, London, England, SW7
2BY,UnitedKingdom
Subsequent computation of the volume integrals in
‡Mailingaddress: CEA-ESEME,Institut deChimiedelaMati`ere
Eq.(3) gives a new value of ∂p/∂t from which the pres-
Condens´ee de Bordeaux, 87, Avenue du Dr. A. Schweitzer, 33608
PessacCedex, France sure p is corrected. The correctionstep is repeated until
2
convergence. This approach has been extensively used variables at the boundaries. Such a numerical method is
by several groups [6, 7] in one-dimensional (1D) calcu- calledBoundaryElementMethod (BEM)andis broadly
lation in conjunction with the restricted cubic EOS [8] used in heat transfer problems. The BEM is expounded
and the finite difference numericalmethod. However,its in Appendix A.
extensionto higher dimensions induces a largecomputa- Several simplifications have to be introduced before
tional effort. First, it requires a sophisticated program- presenting the formulation of the energy and pressure
ming and, second, the computer resources rise steeply equations used in the present method. As stated in the
because the thermodynamic variables have to be evalu- introduction, we are mainly interested in the description
ated at each grid point of the computational domain by of the early stages of the heating, i.e. on times of the
means of the iterative procedure to solve Eq. (3). order of the Piston effect time scale [12]. In this regime,
Teams with expertise in theoretical hydrodynamics the thermodynamic quantities c ,c ,k,... are constant
p v
have developed a rigoroushydrodynamic approach[1] in inthebulkofthefluidandonlyvaryintheverythinhot
which the Navier-Stokes and energy equations are cou- and cold boundary layers. Thus, the following assump-
pledwiththeEOS.Thissetofequationshasbeensolved tions can be made:
analyticallyin1Dbyasymptoticmatchingtechniques[9] (i) The spatially varying parameters (c ,c ,k, etc.) in
p v
and also by direct numerical simulation (DNS) [10] via Eq. (1) can be replaced by the spatially homogeneous
thefinitevolumemethod[11]bothin1Dand2D.Forthe timedependentvalueswhichwillbedenotedhereafterby
sake of simplicity, the van der Waals EOShas been used an over-bar (e.g. c¯ ). These values are calculated with
p
in all these works. Although this simple EOS provides the EOS using density ρ (where the brackets indicate
h i
satisfactory qualitative results, accuracy can only be en- the spatial average) and the pressure p = p(t). We note
suredbyusingarealisticEOS.However,asshowninthe that the density ρ is constant as the system is closed.
h i
presentstudy,insertingarealisticEOSintothehydrody- This assumption is equivalent to the statement that
namicequationsleadstoamuchmoredifficultcomputa- all these quantities should be calculated using the tem-
tionaltask,whichinvolvesprohibitivelylargecalculation perature value T¯ = T¯(t) obtained from the EOS for the
times. For instance, to reach the final steady state in pressurepandthe density ρ . Thequantitiescalculated
h i
one of the 1D runs carried out for the cubic EOS in the in such a way correspond to the over-bar variables men-
present study, one requires about two months on a 800 tionedabove. Ingeneral,for a thermodynamic quantity
MHz PC. X, X¯ = X .
6 h i
The purpose of this article is twofold. We first formu- (ii) During the system evolution, the thermodynamic
late an approximate method which is both simple and quantities are supposed not to vary sharply in time.
rapid(e.g.,thesamecalculationcitedaboverequiredonly (iii) The initial temperature T0 is uniform.
20s!) andthencompareitwiththeDNSformalism. Sec- The assumption(iii) is made for simplicity andcanbe
ond,sincethisworkisthefirstonetousearealisticEOS relaxed if necessary. However, the assumptions (i) and
forthe DNS, wedescribe indetailthe hydrodynamicap- (ii) are essential for this approach.
proachforthenear-criticalfluidswithageneralEOS.We In the following section we formulate the energy and
expectthatthiscompleteandunifieddescriptionmaybe pressure equations that governthe kinetics of the super-
useful for the scientific community, as the hydrodynamic critical fluid in a reduced gravity environment.
method is dispersed over many (some of them not easily
accessible) publications.
The article is organized as follows. In Sec. II we de- A. Energy equation
scribe the fast calculation method. Sec. III presents the
hydrodynamicapproach. Sec.IVdeals withthe compar- Usingtheassumption(i)andtheconstancyof ρ ,one
h i
isonbetweenbothapproaches. Theconclusionsaregiven can write
in Sec. V.
dp=(∂p/∂T) dT¯ (5)
ρ
so that the term (2) reduces to
II. FAST CALCULATION METHOD
c¯ dT¯
g¯(T¯)= 1 v , (6)
Themethodisbasedonthethermodynamicapproach, − c¯ dt
(cid:18) p(cid:19)
i.e. on the energy equation (1). However, a differ-
ent pressure equation will be used instead of Eq. (3). and Eq. (1) can be reduced to the equation
While the latter equation integrates over the fluid vol-
∂ψ
ume, we are looking for a pressure equation that in- =D¯ 2ψ. (7)
tegrates only over the boundaries of the fluid domain. ∂t ∇
Such an equation would accelerate the iterative proce- The thermal diffusion coefficient D¯ = k¯/ ρ c¯ depends
p
dure. However,itsadvantagewouldnotbedecisivewith- on T¯ i.e. on time t only, and h i
out an appropriate numerical method for Eq. (1) that
should only require computation of the thermodynamic ψ(~x,t)=T(~x,t) T E(T¯). (8)
0
− −
3
with ~x the position vector and The use of appropriate thermodynamic relationships
leads to
T¯
c¯
E(T¯)= 1 v dT¯. (9) χT T ∂p ρδs
TZ0 (cid:18) − c¯p(cid:19) δp= (cid:28)cp (cid:16)∂T(cid:17)ρ (cid:29). (15)
χT ρc
The initial condition is ψ = 0 because according to cp v
t=0
| D E
the assumption (iii), In order to use the second law of thermodynamics
T¯ =T . (10) δQ
t=0 0
| ρδs = , (16)
h i vT¯
Finally, a known dependence of D¯ = D f(T¯), where
d
D is a dimensionalconstantandf is a non-dimensional where δQ is the total change of the amount of heat of
d
function, allowsthe time t to be replacedby a new inde- the fluid, one needs to separate out the averages of the
pendent variable τ defined by the equation form YZ in Eq. (15). Under the assumption that Y
h i
(or Z) does not vary sharply over the fluid volume, the
dτ =f(T¯), (11) following approximation holds (see Appendix C):
dt
YZ Y Z . (17)
whose initial condition can be imposed as τ = 0. h i≈h ih i
t=0
Since T¯ is a function of t only, this initial value| problem Among the quantities that appear in Eq. (15), only χ
T
is fully defined. The substitutionofEq.(11)intoEq.(7) and c vary sharply near the critical point and could
p
results in the linear diffusion problem with the constant thus vary strongly across the fluid volume. However,
diffusion coefficient D only their ratio, which remains constant near the crit-
d
ical point, enters Eq. (15). Hence, the average of this
∂ψ
=D 2ψ, ratio as well as the remainder averagesof slowly varying
∂τ d∇ (12) quantities can be separated. By using the expression for
ψ =0.
|τ=0 the total heat change rate
It can be solved with BEM as shown in Appendix A.
δQ ∂T
Usually, the ”temperature step” boundary condition = k dA, (18)
δt ∂~n
has been applied for 1D problems. This heating process ZA
correspondsto a fluid cell, initially ata uniformtemper- where the r.h.s. is simply the integrated heat flux sup-
ature,whichissubmittedtoasuddenincreaseoftemper- plied to the fluid through its boundary A (with ~n the
ature at one of its boundaries, while the other is kept at external normal vector to it), one gets to the final ex-
the initial temperature (Dirichlet boundary conditions). pression
This heating condition is physically unrealistic because
the initial value for the heat flux at the heated bound- dT¯ 1 ∂T
= k dA. (19)
ary is infinite. Instead, in this work we use Neumann- dt ρ vc ∂~n
h i v ZA
Dirichletboundaryconditions: aheatfluxq isimposed
in
In the fastcalculationmethod, Eq.(19) plays the roleof
atoneoftheboundaries,whiletheinitialtemperatureT
0
the pressure equation (3). Equation (19) is both substi-
is maintained at the other boundary.
tuted directly into Eq. (6) and solvedto get the temper-
ature T¯ using the initial condition (10). The obtained
value for T¯ is used to solve Eq. (11) and to calculate all
B. Boundary form of the pressure equation
the fluid properties. Note that T¯ should notbe confused
with T from Eqs. (18, 19). The spatially varying fluid
Let us begin by writing the linearized relationship,
temperature T has to be calculated with Eqs. (8,12).
valid under the assumption (ii):
Substituting T¯ by T , Eq. (19) coincides with the re-
h i
∂ρ ∂ρ sultofOnukiandFerrell[2]whichwasderivedbyadiffer-
δρ= δs+ δp, (13) entway. Eq.(19),writtenintermsof T ,wasemployed
∂s ∂p
(cid:18) (cid:19)p (cid:18) (cid:19)s recently [13, 14] to simulate the gravithatiionalconvection
where s is the fluid entropy per unit mass and δ stands in2Dbythefinitedifferencemethod. However,thefinite
for the variation of the thermodynamic quantity during difference numerical method is not the most efficient for
the time interval δt. From mass conservation it follows the computation of heat transfer problems.
δρ =0,and,fromthepressurehomogeneitythat δp =
h i h i
δp.
By averaging Eq. (13) one obtains III. HYDRODYNAMIC APPROACH
∂ρ ∂ρ Analytical analysis as well as direct simulations were
δs + δp=0. (14)
h ∂s i h ∂p i carried out in previous works. Bailly and Zappoli [15]
(cid:18) (cid:19)p (cid:18) (cid:19)s
4
have developed a complete hydrodynamic theory of den- B. cv-formulation
sity relaxation after a temperature step at the bound-
ary of a cell filled with a nearly supercritical fluid in
In the DNS, the energy equation (22) is re-written in
microgravity conditions. In [15] they describe the dif-
terms of temperature T. This is achieved by expressing
ferent stages of the fluid relaxation towards its complete
theinternalenergyasafunctionofdensityandtempera-
thermodynamic equilibrium, covering the acoustic, Pis-
ture so thatone canmakeuse ofthe wellknownrelation
ton effect and heat diffusion time scale. The analytical
approachleansonthematchedasymptoticexpansionsto
de 1 ∂p dρ dT
shoelavte-dtihffeus1iDngN, anveiaerr--cSrtiotikceaslevqaunatdieornsWfoaralas vfliusicdou(ss,eelo[w9-] dt = ρ2 "p−T (cid:18)∂T(cid:19)ρ# dt +cv dt . (24)
and [15]). The DNS of the Navier-Stokesequations were
performed in 1D and 2D geometries. Some of them take Then,bysubstitutingtheEqs.(20)and(22)intoEq.(24)
into account gravity effects, as for example, the interac- one obtains:
tionofanear-criticalthermalplumewithathermostated
dT ∂p
boundary [16]. Numerical results are also available on ρc = (k T) T ~u+Φ. (25)
v
dt ∇· ∇ − ∂T ∇·
thermo-vibrationalmechanisms [17]. (cid:18) (cid:19)ρ
To date the hydrodynamic approach has been solved
Note that Eq. (25) involvesc and not c as in the ther-
forthe classical,vander Waals,EOS.ThisEOSallowsa v p
modynamic Eq. (1). The “c -formulation” is preferred
considerablereductionincomputationaltimewhencom- v
to the “c -formulation” because the much weaker near-
pared to the restricted cubic EOS. However, it does not p
critical divergence of c (in comparison to c ) allows c
provide a correct description of the real fluids. In par- v p v
to be assumed constant.
ticular, it fails to predict the critical exponents for the
divergencelawsofthe thermodynamicproperties. Inthe The boundary conditions for the Navier-Stokes equa-
present work we use a more realistic cubic EOS to de- tions are ~u = 0 at the walls. The initial conditions are
scribethefluidbehaviorinthenear-criticalregion. Here- given by ~u(t = 0)= 0. For the energy equation (25) the
after,we describethe methodologysuitable fora general boundaryandinitialconditionsareidenticaltothoseap-
EOS. plied in the energy equation (1) (cf. Sec. II). The values
of the physical parameters used in the simulations are
discussed in Appendix B.
A. Problem statement
C. Acoustic filtering
The hydrodynamic description leads to the following
set of equations Heattransfer insupercriticalfluid involvesthree char-
acteristic time scales [12, 18]: the acoustic time scale
dρ defined by t =L/c (where c is the soundvelocityand
+ρ ~u=0, (20) a 0 0
dt ∇· L is the cell size), the diffusion time scale t = L2/D
D
d~u (D being the thermal diffusivity) and the Piston Effect
ρ = p+µ 2~u, (21)
dt −∇ ∇ time scale defined by tPE = L2/[D(cp/cv 1)2], with
−
de ta tPE < tD. The present study is mainly concerned
ρ = (k T) p ~u+Φ, (22) ≪
with time of the same order as the Piston effect time
dt ∇· ∇ − ∇·
scalesothatafinedescriptionoftheacousticphenomena
where e is internalenergyper unitmass,~u=(u ,u ,u ) is not needed. This suggests that one can filter out the
1 2 3
is the fluid velocity at the point ~x=(x ,x ,x ), acousticmotionsofthesetofEqs.(4,20,21,25)andretain
1 2 3
only their integrated effects without altering the physics
of our problem. The removal of the acoustic motions
∂u ∂u ∂u ∂u 2∂u ∂u
i j i i i j
Φ=µ + isachievedbyapplyingtheacousticfilteringmethod[19]
∂x ∂x ∂x ∂x − 3∂x ∂x
i,j (cid:18) j i j j i j(cid:19) whichisbroadlyusedinthecomputationofthelowMach
X
number compressible Navier-Stokesequations because it
is the dissipation function due to the shear viscosity µ avoidsnumericalinstabilitieswhentimesteps,∆t>>ta,
(the bulk viscosity is neglected). The operator d/dt is are used in the simulations. The following presents the
defined as main points of the acoustic filtering method.
Theequationofmomentumisfirstrewrittenbychoos-
d ∂ ing the sound velocity c as the reference velocity scale
= +~u . (23) 0
dt ∂t ·∇ and L/u0 as the reference time scale (here u0 is the
characteristic velocity of large scale fluid motions, in
The set of equations (20,21,22) is closed by adding the our case u = L/t ). Using this time and velocity
0 PE
EOS (4). scale the Mach number Ma=u /c appears in the non-
0 0
5
dimensional momentum equation as follows, SIMPLERalgorithmand by applyingthe Finite Volume
Method (FVM, see Appendix D) on each grid cell of the
∂~u Ma−1p 1
ρ +Ma−1(~u )~u = c p+ 2~u, 1D cell. Near the walls the mesh is refined to properly
∂t ·∇ − c2ρ ∇ Re∇
(cid:20) (cid:21) 0 c resolve the very thin thermal boundary layers.
(26)
In the present work the thermodynamics variables are
where Re = ρ u L/µ is the Reynolds number and the
c 0 determined using the parametric EOS [8]. This uses two
densityandpressurearenon-dimensionalizedbythecrit-
parameters r and θ which both depend on temperature
ical density ρ and critical pressure p taken as the ref-
c c T anddensityρ. Therefore,oneneedstosolvetwoequa-
erence values. For small Mach numbers, one can express
tions
the fluid variables as series of Ma,
Λ (r ,θ ,T )=0
~u=Ma[~u(0)+Ma2~u(1)+o(Ma2)], (27) 1 i i i (34)
Λ (r ,θ ,p)=0
2 i i
p=p(0)+Ma2p(1)+o(Ma2), (28)
insteadofoneEq.(4)foreachvolumeelementiandtime
While ~u in the l.h.s. of Eq. (27) is non-dimensionalized
step.
with c , the term in the square brackets defines the ve-
0 The whole numerical procedure consists in solving by
locity non-dimensionalized with u . This explains the
0 the Newton-Raphson method, at each time step, a set
factor Ma in Eq. (27). The density and temperature are
of equations that includes Eqs. (34), written for each
expanded like p in Eq. (28). By substituting the series
volume element and Eq. (33). This makes a system of
(27, 28) into Eq. (26) and neglecting the terms of or-
2N+1equationstoresolve,N beingthetotalnumberof
der O(Ma), one obtains p(0) = 0, which means that
the volume elements. The local temperature T is given
p(0) depends on time onl∇y. By retaining O(Ma) terms i
by the resolution of Eqs. (29-31) at each iteration of the
in Eqs. (20,25,26) and O(1) terms in the EOS (4), one
SIMPLER algorithm for each time step, as described in
obtains the final (dimensional) form for the governing
Appendix D.ForeachvalueT the 2N+1(r ,θ ,p)vari-
i i i
equations:
ables are computed via the system (34).
dρ(0)
= ρ(0) ~u(0), (29)
dt − ∇·
IV. RESULTS AND DISCUSSION
d~u(0)
ρ(0) = p(1)+µ 2~u(0), (30)
dt −∇ ∇
A brief analysis comparing the c and c formulations
dT(0) ∂p p v
ρ(0)c(0) = T(0) ~u(0)+ (1)and(31)oftheenergyequationallowsustogainmore
v dt − ∂T ∇· insight into the relation between the two approaches.
(cid:18) (cid:19)ρ
(k T(0)), (31) Formally, Eq.(1) and Eq.(31) become equivalent if the
∇· ∇ advection term
Λ(p(0),ρ(0),T(0))=0, (32)
(~u )T (35)
where (pc/c20ρc)p(1) is replaced by p(1) for the sake of ·∇
compactness. The pressure term p(1) has to be inter- is added to l.h.s. of Eq. (1). However, the equivalence
preted as the dynamic pressure that makes the veloc- of the two forms under which the pressure work appears
ity field satisfy the continuity equation (29). This term (seesecondtermofthel.h.sofEq.(1)andEq.(31))isnot
reflects the contribution of the acoustic waves averaged trivial and deserves to be detailed. At the early stage
overseveralwaveperiodstothetotalpressurefield. One of the heating (t < t ) the velocity at the front of
PE
notesthatthevelocityscalec0 isnotpresentanymorein the cold boundary layer being very small, the velocity
Eqs.(29-32) which was the main purpose of the acoustic can be assumed to decrease linearly in the bulk cell as
filtering. ∂u/∂x u /L, where u is the maximum veloc-
max max
The assessment of p(0) requires one more equation to ity loca≃ted−at the front of the hot boundary layer and x
close the set (29-32). This additional equation expresses the distance fromthe hotwallto the coldwall. The rate
the mass conservation: of temperature increasedue to the pressure contribution
1 in Eq. (31) can thus be written as follows,
ρ(0)dv = ρ , (33)
v h i
Zv T ∂p T ∂p u
max
where ρ is a known constant. ~u= (36)
Inthheifollowing,the superscript(0)is droppedtocon- − ρcv (cid:18)∂T(cid:19)ρ∇· ρcv (cid:18)∂T(cid:19)ρ L
form to the notation of Sec. II. By using the expression [12],
1 ∂T δQ
u = , (37)
D. Numerical procedure max
T ∂p Aδt
(cid:18) (cid:19)ρ
For the time integration, the first order Euler scheme and Eqs. (18,19) one concludes that the term (36) is
is used. Equations (29-31) are solved by the iterative equivalent to (6) near the critical point where c c .
p v
≫
6
14 15
6)
&2
(cid:25)
12 (cid:21)
7 7 (cid:14)(cid:20). 7 7 (cid:14)(cid:20).
(cid:19) F
10 (cid:19) F
7(cid:16)7(cid:3) (mK)(cid:19) 68 T2LQ1 .7(cid:21)9s(cid:3):(cid:18)P(cid:21) 17.43s 7(cid:16)7 (mK)0 10 83.05s TLQ (cid:21)(cid:3):(cid:18)P(cid:21)
5
4
57.5s
2 31.94s
8.72s
13.08s
6.39s
0
0 1 2 3 4 5 0
0 1 2 3 4 5
[ (mm) [ (mm)
(a) (a)
5 1.5 40 2
T
35 RXW
4
DNS
30 1.5
7(cid:16)7 (mK)FHQWHU(cid:19) 23 TRXW DNS7FHQWHU New method 01.5 RXW (W/m)T2 7(cid:16)7(cid:3)(mK)FHQWHU(cid:19)122505 7FHQWH45U New metho1.d5 1RXW (W/m)T2
1
3 W
W 10 0.5
1 3( 2 3( 0.5
1
5
0 0
0 0 0 5 10 15 20
0 5 10 15 20 0 0
0 500 1000 1500
W (s) W (s)
(b) (b)
FIG.1: ComparisonoftwoapproachesforSF6at1KaboveTc FIG. 2: Comparison of the two approaches for CO2 at 1K
(reducedtemperature3.1·10−3),qin =2W/m2andhρi=ρc. above Tc (reduced temperature 3.3·10−3), qin = 2 W/m2
SolidcurvesaretheDNSresultsandthedottedcurvesarethe andhρi=ρc. SolidcurvesaretheDNSresultsandthedotted
newmethodresults. (a) Spatialvariation of thetemperature curves are the new method results. (a) Spatial variation of
at different times. (b) Time evolution of the temperature at thetemperatureatdifferenttimes. (b)Timeevolutionofthe
thecellcenterandofthefluxattheexitofthecell. Thevalue temperature at the cell center and of the flux at the exit of
oftPE =7.73sobtainedwithourEOSisshownbyanarrow. thecell. Thevalueof tPE =3.45 sobtained with ourEOSis
shown by an arrow in the insert that presents the short time
evolution.
One can note that in the hydrodynamic approach, the
pressureworkisdirectlyrelatedtothemasstransferfrom
perature gradient is small. At very large times, t > t ,
the hot boundary layer to the bulk fluid via the gradi- D
the Piston effect is not efficient and the velocity tends
ent velocity. It is then very important to asses properly
to zero. We thus do not expect a strong influence of the
the effect of the velocity field in order to compare the
advection effects on the temperature field. This will be
fast calculation and hydrodynamic methods. The above
confirmed by the results presented below. Note that the
analysis has shown that the expressions of the pressure
advectiontermcannotbe neglectedwhenthe flux distri-
workisequivalentfor bothmethods. Hence, the remain-
bution over the heater surface is highly inhomogeneous.
ingpotentialinteractionbetweenthevelocityandenergy
Hot jets [26] can be generated in this case.
fieldscanmanifestitselfonlythroughtheadvectionterm
(35). This term is only relevant when, at the same spot The calculations have been performed for two fluids,
of the fluid, both the fluid velocity and the temperature CO andSF ,confinedinacelloflengthL=5mm. The
2 6
gradientarelarge. Atthe smalltimes, t<t , the tem- initial temperatures 1K and 5K above the critical point
PE
perature gradients are confined very near the wallwhere havebeenconsideredforCO . Thecomputationsrelated
2
the velocity remains small [18]. Later on, the velocity toSF concernonlytheinitialtemperature1Kabovethe
6
maximum shifts to the center of the cell where the tem- critical point. The cell boundary situated at x = 0 has
7
been submitted to the constant heat flux q = 2W/m2 outtwobehaviors. Ononehand,thefluxq appearsto
in out
(forT =T +1K)andq =9.5W/m2(forT =T +5K) fit very well with the DNS over the full time evolution,
0 c in 0 c
and the opposite boundary has been maintained at the at the time scale t as well at the time scale t , see
PE D
constant temperature T(x=L)=T . Figs. 2b, 3b. On the other hand, the temperature at the
0
Thetimeevolutionofthetemperatureprofilesandthe cell center T tends to be lower than the DNS data.
center
temperature at the cell center T = T(x = L/2) as This discrepancy increases with time and is larger when
center
wellastheheatfluxq = k∂T(x=L) arecomparedand the temperatureis closertothe criticaltemperature(see
analyzed. In the caseouotf CO− , th∂extime evolution covers Fig. 2b and Fig. 3b). Both behaviors can be explained
2
not only the Piston effect time scale t but also the by considering how the thermal conductivity k is esti-
PE
large diffusion time scale. mated in each method. In the hydrodynamic approach
k is determined locally,whereas the fast calculationuses
600 the spatial average value of k. Thus, keeping in mind
&2
that the thermal conductivity diverges when approach-
(cid:21)
500 7 7 (cid:14)(cid:24). ing the critical point, the increment in temperature near
(cid:19) F theheatingsurfacetendstobesmallerinthenewmethod
K) 400 642.85s T (cid:28)(cid:17)(cid:24)(cid:3):(cid:18)P(cid:21) tthhaenteimnptehreatDuNreSis(sfiexeeFdiagt. t2hae).inAititatlhteemoppeproastiuteresuarnfdacies
m LQ
7 (0 300 closesttothebulktemperaturesothattheeffectofaver-
7(cid:16) 2892.82s aging k is less influent in this region. One can note that
the thermal diffusivity can be computed locally in the
200
fast calculation method by applying the Kirchhoff sub-
321.42s
stitutionofthedependentvariableψ(definedbyEq.(8))
100
ψ
by φ= k(ψ)dψ.
0
0 A physical interpretationof the temperature T¯, which
R
0 1 2 3 4 5 was formally introduced in Sec. II, can now be given.
[ (mm) Indeed, since in the fluid bulk (i.e., arround the center
(a) of the cell) ∂T = 0 at t < t , according to Eq. (1)
∂x PE
300 10 we have ∂T /∂t = g¯(T¯). As near the critical point
center
T c c , Eq. (6) provides g¯(T¯) ∂T¯/∂t. Finally, one
250 RXW cpan≫convclude that T¯ T . In≈other words, T¯ can be
8 ≈ center
7 considered as the bulk temperature.
K) 200 FHQWHU T Aaafurtherremark,wenotethatfor1DtheEq.(A1)
7(cid:16)7(cid:3)(mFHQWHU(cid:19)150 46RXW (W/m)2 cNeWoreaeuvlliendtryoththeaaenvlteedhsabsi,ttesewinnpeo2spsoDrsleivbfaeelnderdatehnx3eaDtelyunttsshieiceoaonBlflytEtohMbehyirBgseehEmreMiraeisndfeosimxrapediatnvsnsaisgnoioentnans--..
100
geousinresolvinglinearizedproblemswhencomparedto
2 other numerical methods. Its success is based on several
50
factors. One of them is its numerical stability: the nu-
merical solution of the integral equations is much more
0 0
stable than that of the differential equations and allows
0 500 1000 1500
W (s) the use of larger time steps. Another advantage consists
(b) in the possibility of determining analytically the BEM
coefficients,Eqs.(A6,A7). For2Dconfigurations,thedi-
agonalcoefficientsG andH (whichhavethelargest
FF FF
FIG. 3: Comparison of two approaches for CO2 at 5K above absolutevalueandthusarethemostrelevant)canbecal-
the critical temperature (reduced temperature 1.6 · 10−2),
culatedanalytically. The semi-analyticalintegrationcan
qin = 9.5 W/m2 and hρi = ρc. Solid curves are the DNS
be used for the remaining coefficients [22, 23].
results and the dotted curves are the new method results.
According to our EOS, tPE = 25.67 s. (a) Spatial variation
of the temperature at different times. (b) Time evolution of
the temperature at the cell center and of the flux at the exit V. CONCLUSIONS
of the cell.
In this work we propose a thermodynamic method for
ThesetofFigs.1,2,and3exhibitaverygoodqualita- describing the heat transfer in supercritical fluids in ab-
tiveagreementbetweentheDNSandthefastcalculation sence of gravity effects. The method has been compared
results. The thin boundary layers and the homogeneous with the solution of the full hydrodynamic equations,
enhancement of the temperature in the center cell are showing an excellent agreement. In general, a thermo-
very well predicted. The quantitative comparison sets dynamic approach leans on the possibility of expressing
8
tsho,etphreetsrsaunresfewroorfkminodmeepnetnudmendtloyesonfottheneveedlotcoitbyeficeolnds.idI-f ψ(x~′,t′)∂x′G(~x−∂~nx~′,t−t′) dx′A= 21ψ(~x,t). (A2)
#
ered, allowing a large reduction in computational time.
As an example, in calculations carried out for CO and The integration is performed over the surface A of the
2
SF6,the presentthermodynamicmethodwithinminutes fluid volume v, ~x A. The outward unit normal to A
∈
providedthecompleteevolutionoftheheattransferpro- is ~n. The Green function G for the infinite space for the
cess,while the directnumericalsimulationofthe fullhy- equation adjoint to Eq. (A1) reads
drodynamic equations required weeks of CPU time.
~x2
Compared with previous thermodynamic methods [5], G(~x,t)=(4πDt)−d/2exp | | , (A3)
the fast calculation method presented here does not re- −4Dt
(cid:18) (cid:19)
quire the evaluation of the variables at each cell of the
wheredisthespatialdimensionalityoftheproblem(A1).
computation domain. This fact ensures a much better
Only the 1D case with the space variable x (0,L)
performance. Moreover, the proposed method offers the ∈
is considered below, the 2D counterpart being described
possibility to explicitly include the thermal behaviour of
elsewhere [22, 23]. For d = 1, A degenerates into two
thematerialvesselcontainingthefluidbytakingintoac-
points, and Eq. (A2) reduces to two equations
count the heat conduction along the solid walls, see Ref.
[4]. t
∂ψ
The direct numerical simulation of the flow has been 2D dt′ G(x x′,t t′)
usedtoanalyzethevalidityofthemethodproposedhere. " − − ∂x′ −
Z
0
The accuracy of the latter approach is explained by the
fact that the advection of energy remains negligible. ∂G(x x′,t t′) x′=L
ψ(x′,t′) − − =ψ(x,t) (A4)
For the sake of completeness, we have also presented ∂x′
(cid:21)(cid:12)x′=0
a detailed description of the hydrodynamic approach. (cid:12)
writtenforx=0,L. Avarietyo(cid:12)fnumericalmethodscan
While ithasbeenusedforaboutadecade,somepartsof (cid:12)
be applied to solve Eqs. (A4). The simplest way is to
its description for a general equation of state are either
present the integral over (0,t) as a sum of the integrals
dispersed over many literature sources or not published
over (t ,t ), f = 1..F with t = 0 and t = t and as-
at all in the accessible literature. f−1 f 0 F
sume a constant values ψ =ψ(t ) and ψ′ =∂ψ(t )/∂x
Concerning the future development of the present re- f f f f
search, we plan to extend the fast calculation method overeachoftheseintervals. Eqs.(A4)willreducethento
to two- and three-dimensional problems. Finally, we in- the set of 2Fm (tFm is the maximum desired calculation
tend to use this method to investigate the heat transfer time) linear equations
in two-phase fluids.
F
[ψ (0)H (x) ψ (L)H (x L)
f Ff f Ff
− − −
f=1
Acknowledgments X
ψ′(0)G (x)+ψ′(L)G (x L)]=
f Ff f Ff −
This work was partially supported by the CNES. We ψFi/2, x=0,L; F =1,Fm (A5)
thankCaroleLecoutreforherhelpwiththecodelaunch- for 4F variables ψ , ψ′(0,L), 2F of them being de-
ing. A. D. acknowledges Jalil Ouazzani for the helpful m f f m
fined by the boundary conditions. The coefficients in
discussions on the numerical method used in the hydro-
these equations can be calculated analytically:
dynamic approach.
tf
G (x)=D G(x,t t′)dt′ =
Ff F
−
Appendix A: BEM for the diffusion equation Ztf−1
x2
x exp( u) 4D(tF−tf)
| | erfc(√u) − (A6)
thaInt DthiasnAdptpceonrdreixspwonedutsoeDthdeatnrdadτitoiofnEaql.n(o1t2a)t.ioItn,casno 2 (cid:20) − √uπ (cid:21)(cid:12)(cid:12)4D(tFx−2tf−1)
(cid:12)
be shown [21] that the linear diffusion problem and (cid:12)
∂∂ψt =D∇2ψ (A1) HFf(x)=−D tf ∂G(x,∂txF −t′)dt′ =
ψ =0 Ztf−1
|t=0 x2
sign(x) 4D(tF−tf)
withtheconstantthermaldiffusioncoefficientDisequiv- erfc(√u) , (A7)
alent to the boundary integral equation − 2 (cid:12)(cid:12)4D(tFx−2tf−1)
(cid:12)
t where (cid:12)
D dt′ G(~x x~′,t t′)∂x′ψ(x~′,t′) ∞
" − − ∂~n − erfc(x)= exp( u2)du,
Z0 AZ Zx −
9
is the complementary error function and Appendix D: Application of the Finite Volume
Method (FVM) and SIMPLER algorithm
1, x>0,
sign(x)=
1, x<0. According to the FVM, the calculation domain is di-
(cid:26) −
vided into a number of non-overlapping control volumes
One needs to mention that the case x = 0 is special: sothatthereisonecontrolvolumesurroundingeachgrid
H (0)=0 for all f and F and point. Thedifferentialequationsareintegratedovereach
Ff
control volume. The attractive feature of this method is
that the integral balance of mass, momentum, and en-
G (0)= D /π( t t t t ).
Ff 0 F f−1 F f
− − − ergy is exactly satisfied over any control volume (called
p p p below the cell for the sake of brevity), and thus over the
The set of linear equations (A5) can be solved by any
whole calculation domain. The integral formulation is
appropriate method, e.g. by the Gauss elimination.
also more robust than the finite difference method for
problems which present strong variations of properties
observed in a near-critical fluid [10, 20]. The equations
Appendix B: The fluid properties are resolved on a staggered grid. This means that the
velocity is computed at the points that lie on the faces
of the cell while the scalar variables (pressure, density
The thermal conductivity k is deduced from the ther-
andtemperature) arecomputed at the center ofthe cell.
mal diffusivity
This choice is made to avoid pressure oscillations in the
T T ϕ1 T T ϕ2 computations [11]. For the time discretization, the first
c c
D =D1 − +D2 − (B1) orderEulerschemeisused. Forthesakeofsimplicityand
T T
(cid:18) c (cid:19) (cid:18) c (cid:19) clarity we present the finite volume method for the 1D
generalizedtransportequationforavariableY (whereY
andthe constantpressurespecificheatatcriticaldensity
can be substituted by either u or T)
ρ , k =Dρ c . Coefficients values for CO are:
c c p|ρc 2
D =5.89184 10−8m2/s,D =7.98068 10−7m2/s, ∂ρY ∂ρuY ∂ ∂Y
1 × 2 × + = Γ +S (D1)
ϕ1 =0.67 and ϕ2 =1.24. ∂t ∂x ∂x ∂x
(cid:18) (cid:19)
Coefficients values for SF are:
6 where Γ denotes the generalized diffusion coefficient and
D = 6.457 10−7 m2/s, D = 0, ϕ = 0.877 and
1 2 1 Sthegeneralizedsourceterm(volumeforces). Integrated
×
ϕ =0.
2 overtheithcellofthelengthδx,Eq.(D1)takestheform
The specific heat at constant pressure is calculated by
using the thermodynamic relationship ρPYP −ρpPYPp∆x+J J =S ∆x (D2)
e w P
∆t −
2
∂p wherethesuperscriptpdenotesthevalueontheprevious
c =c +T χ (B2)
p v ∂T T time step, the subscript P represents the center of the
(cid:18) (cid:19)ρ
cell, the subscripts e and w its “east” and “west” face
respectively. The calculation of the flux
The isothermal compressibility coefficient χ and the
T
specific heat at constant volume are given by the re- ∂Y
J =ρuY Γ (D3)
strictedcubicmodel[8]. Forthereferencehydrodynamic
− ∂x
DNSweusedaconstantc valuecalculatedfortheinitial
v on the faces requires the knowledge of Y and ρ at the
value of temperature and density. We used a constant
centers of two neighboring “East” and “West” cells de-
value for the viscosity µ: 3.74 10−5 Pas for SF and
3.45 10−5 Pas for CO . · · 6 notedbythecapital lettersEandW.Theirvaluesatthe
· · 2 faces can be found by linear interpolation between their
valuesatthe centers,e.g. Y =0.5(Y +Y )ifthenodes
e P E
are equidistant.
Appendix C The continuity equation integratedon the controlvol-
ume is given by:
According to the integral theorem about the mean ρ ρp
value [24], there is always a point ~xm ∈v so that P∆−t P∆x+Fe−Fw =0 (D4)
with F = ρu. When multiplying Eq. (D4) by Y and
P
Y(~x)Z(~x)d~x=Y(~xm) Z(~x)d~x (C1) subtracting the result from Eq. (D2), one obtains the
Zv Zv equation
if the functions Y,Z are continuous. When the spatial ρp∆x
P (Y Yp)+(J Y F ) (J Y F )=S ∆x
variation of Y in v is small, Y Y(~xm) and Eq. 17 ∆t P − P e− P e − w− P w P
h i ≈
stems from Eq. (C1). (D5)
10
that can be rewritten in the following form whereu representstheneighborvelocities. usat-
nb
isfies
a Y =a Y +a Y +b (D6)
P P W W E E b
p(1)∗ p(1)∗
The tridiagonalsetof linearequations (D6) with respect ue =ue+ P − E . (D11)
a
e
toY issolvedbytheThomasalgorithm[25]. Thestencil
P 3. Compute the pressure p(1)∗ whose equation is
coefficients aP, aW et aE depend on the discretization deduced by appblying the divergence operator to
scheme. Their general expression is
Eq. (D11) and using the continuity equation (D4):
a =B A +max( F ,0),
W w w w
− ρ ρ ρ ρ
aaPE ==aBWeA+e+aEm+axρ(pPF∆e,x0/)∆, t, (D7) (cid:18)aww + aee(cid:19)p(P1)∗ = awwp(W1)∗+ aeep(E1)∗+
b=SP∆x+ρpPYPp∆x/∆t, ρpP −ρP∆x ρ u +ρ u . (D12)
e e w w
∆t −
where B = Γ/∆x. We use the “power law scheme” [11]
that requires b b
4. Solve Eq. (D8) with p(1)∗ used for p(1) and thus
5 obtaining u∗.
0.1F
i
A =max 0, 1 | | , i=e,w.
i
" (cid:18) − Bi (cid:19) # 5. Compute p(1)′ whose equation is obtained analo-
gously to Eq. (D12) from
Theset(D6)shouldbewrittenandsolvedbothforthe
velocity and the temperature. While the above scheme (p(1)′ p(1)′)
canbedirectlyappliedforthetemperaturecase,thecou- u =u∗+ P − E . (D13)
e e a
pling of the velocity and the pressure p(1) (which is de- e
fined implicitly by the continuity equation) requires a
It takes the form
special treatment for the velocity equation as described
below. ρ ρ ρ ρ
w + e p(1)′ = wp(1)′+ ep(1)′+
The non-dimensionalized and discretized Navier- a a P a W a E
Stokes equation (30) (cid:18) w e(cid:19) w e
ρp ρ
P − P∆x ρ u∗+ρ u∗. (D14)
a u = a u +(p(1) p(1))+b, (D8) ∆t − e e w w
e e nb nb P − E
X
wherethesubscriptnbdenotestheneighborsofthepoint 6. Calculate the velocity u using Eq. (D13). Do not
e,canbesolvedonlywhenthepressurefieldisgiven. Un- correctthepressurep(1),p(1)′isusedtocorrectonly
less the correct pressure field is employed, the resulting the velocity field, the pressure being computed by
velocityfieldwillnotsatisfythecontinuityequation. We Eq. (D14).
use the iterative SIMPLER algorithm [11] to couple the
7. Solve the energyequationforT using the obtained
velocity and the pressure fields. This algorithm is based
u values.
on successive corrections of the velocity field and pres-
surefieldatagiventime step. Thevelocityandpressure
8. Calculate the density distribution and p(0) via
variables are decomposed as follows:
Eqs. (33,34).
u=u∗+u′,
p(1) =p(1)∗+p(1)′, (D9) 9. Return to step 2 and repeat until the converged
solution is obtained.
wheretheasteriskdenotestheguessesandprimethecor-
Ithas to be noted thatwhereasthe fractionalstep PISO
rections. The steps of the SIMPLER algorithm are the
algorithm[10]issuccessfulinresolvingtheequations(29-
following:
32)ontheacoustictimescale,itisnotthecasewhenthe
1. Start with a guessed velocity field. acoustic filtering method is used. Due to the different
meanings of pressure (see the subsection IIIC) in the
2. A pseudo-velocity u (without taking into account momentum equation (involving p(1)) and in the energy
the pressure gradient) is first computed and is de- equation (involving p(0)), it appears that only an itera-
fined as b tive algorithm can correctly couple the thermodynamic
field and the velocity field, the PISO algorithm leading
a u +b
nb nb
u = (D10) to unstable solutions.
e
a
P e
b