International Journal of Thermal Sciences xxx (2014) 1e12
Contents lists available at ScienceDirect
International Journal of Thermal Sciences journal homepage: www.elsevier.com/locate/ijts
Modeling supercritical heat transfer in compressible fluids Patrícia C. Teixeira a, Leonardo S. de B. Alves b, * a ^mica e Hiperso rio de Aerotermodina ^nica e LAH, Instituto de Estudos Avançados e IEAv, Centro T Laborato ecnico Aeroespacial e CTA, ~o Jos Trevo Coronel Aviador Jos e Albano do Amarante 1, Putim, Sa e dos Campos, SP 12228-001, Brazil b ^nica Teo ^nica e TEM, Universidade Federal Fluminense e UFF, rio de Meca rica e Aplicada e LMTA, Departamento de Engenharia Meca Laborato tria 156, bloco E, sala 216, Nitero i, RJ 24210-240, Brazil Rua Passo da Pa
a r t i c l e i n f o
a b s t r a c t
Article history: Received 16 November 2013 Received in revised form 9 August 2014 Accepted 9 August 2014 Available online xxx
Novel solutions have been derived for both thermodynamic and hydrodynamic models of the heat transfer inside a cavity containing supercritical fluid in zero gravity. A fully analytical solution of the thermodynamic model was obtained through a combination of the Generalized Integral Transform Solution and the Matrix Exponential Method. Its accuracy is entirely controlled by a single user prescribed parameter. Furthermore, a low Mach Preconditioned Density-Based Method was employed to generate a numerical solution of the hydrodynamic model, avoiding acoustic filtering and the need to resolve acoustic time scales without it. A proper model for the piston effect evolution in pseudo-time must be included to generate a physically correct description of its physical-time evolution. Furthermore, both models generate graphically identical results, but only upon a thermodynamically consistent selection of fluid properties and equation of state. Finally, the theoretical expression for the piston effect relaxation time underestimates the actual value of this characteristic time estimated from a simulation of the same model used to derive this expression. This feature is not an artifact of boundary condition choice. © 2014 Elsevier Masson SAS. All rights reserved.
Keywords: Thermodynamic critical point Differential equation of state Piston effect Low Mach preconditioning Integral transforms Analytic solution
1. Introduction All fluids are subject to a universal divergence of their thermodynamic properties as their critical points are approached. Even though the thermal diffusion relaxation time diverges near the critical point, a very fast temperature equilibration is still observed in enclosed samples on both ground-based [1] and microgravity [2] experiments. Such a critical speeding-up has been explained using simulations of a thermodynamic model [3,4] as well as a hydrodynamic model [5]. It is caused by the ability of small temperature perturbations to create severe compression in near-critical fluids, which generates thermo-acoustic waves. When entrapped within cavity walls, their propagation and reflection induces a rapid heating of the entire fluid, resulting in a homogeneous increase of its bulk temperature. This fast temperature relaxation phenomenon is known today as piston effect and it has been the subject of several reviews in the literature [6e9]. The original thermodynamic model developed for supercritical heat transfer under microgravity conditions [4,3] is essentially the heat conduction equation with a source term proportional to the bulk temperature time derivative. This source term models the
* Corresponding author. Tel.: þ55 21 2629 5576; fax: þ55 21 2629 5588. E-mail address:
[email protected] (L.S.B. Alves).
adiabatic compression mechanism responsible for the piston effect. It becomes dominant whenever g [ 1, i.e. closer to the critical point. On the other hand, it vanishes in the incompressible limit, where g ¼ 1. A solution to this equation, including variable properties obtained from available data in the literature, was first generated with numerical methods [4]. However, bulk temperatures are obtained by numerical integration of the temperature field over the entire simulated volume at each time step. Such an approach can become costly for two and three-dimensional domains. In order to circumvent this problem, an alternative approach was proposed that replaces the volume integral by a surface integral using the Boundary-Element Method [10]. An additional solution to this equation, now considering constant properties, was first obtained with an approximate Fourier transformation procedure [3], where a detailed derivation is provided elsewhere [11]. Although both studies considered steady heating at the boundaries, the approximate Fourier procedure has been extended towards pulsed [12] and unsteady [13] heating as well. The original work [11] also considered a conjugate problem, coupling the thermodynamic model for the fluid with a heat conduction model for the solid walls containing a composite material. It was extended later to include curvature effects due to cylindrical container walls [14], although a different approximate analytical approach was utilized. Separation of variables was applied to the solid wall heat conduction problem whereas a Laplace transform with numerical
http://dx.doi.org/10.1016/j.ijthermalsci.2014.08.011 1290-0729/© 2014 Elsevier Masson SAS. All rights reserved.
Please cite this article in press as: P.C. Teixeira, L.S.B. Alves, Modeling supercritical heat transfer in compressible fluids, International Journal of Thermal Sciences (2014), http://dx.doi.org/10.1016/j.ijthermalsci.2014.08.011
2
P.C. Teixeira, L.S.B. Alves / International Journal of Thermal Sciences xxx (2014) 1e12
Nomenclature b A
A Ai,j B Bi,j BDF c b c CP CV e E Ei Ev F h hP hT H H k L c M Ma N Ni NT Nx P PH PT Q b Q t tD tPE T T u x
inviscid flux Jacobian based on primitive variables matrix integral transformed matrix integral transformed matrix coefficient integral transformed matrix integral transformed matrix coefficient backwards difference formula sound speed preconditioned sound speed specific heat at constant pressure specific heat at constant volume internal energy per unit mass total internal energy per unit mass inviscid flux vector viscous flux vector dimensionless temperature filter enthalpy per unit mass preconditioned enthalpy Jacobian with respect to temperature preconditioned enthalpy Jacobian with respect to pressure total enthalpy per unit mass source term vector thermal conductivity cavity length preconditioned inviscid flux Jacobian matrix Mach number number of terms in summation series solution norm number of terms in temporal grid number of terms in spatial grid pressure hydrodynamic pressure thermodynamic pressure conservative dependent variable vector primitive dependent variable vector physical-time coordinate thermal diffusion relaxation time piston effect relaxation time temperature conservative to primitive dependent variable Jacobian matrix velocity spatial coordinate
inversion was applied to the supercritical fluid problem. Recently, the Generalized Integral Transform Technique was applied to a variation of this thermodynamic model where all fluid properties have a linear dependence on both temperature and pressure [15], confirming that variable properties have little impact on the bulk temperature evolution [16]. This hybrid technique analytically transforms the thermodynamic model into an unsteady system of ordinary differential equations, which is then numerically marched forward in time. Such simplified models are able to predict bulk temperature and density evolutions on piston effect and thermal diffusion time scales, respectively [10]. However, pressure is space averaged since the speed of sound is assumed much faster than any other local fluid speeds. If an accurate description of supercritical heat transfer is desired on acoustic time scales, the propagation
Greek symbols thermal diffusivity isobaric thermal expansion coefficient eigenvalue increment Kronecker delta absolute error ε preconditioned sound speed control parameter hi integral transform coefficient g ratio between specific heats G preconditioning matrix kT isothermal compressibility m dynamic viscosity ji eigenfunction ~i j normalized eigenfunction q homogeneous dimensionless temperature q integral transformed homogeneous dimensionless temperature vector qi integral transformed homogeneous dimensionless temperature Q dimensionless temperature Qb dimensionless bulk temperature r density rP preconditioned density Jacobian with respect to temperature rT preconditioned density Jacobian with respect to pressure t dimensionless pseudo-time coordinate tPE dimensionless piston effect relaxation time x dimensionless spatial coordinate
a aP bi d di,j D
Subscripts 0 initial state, right wall 1 left wall b bulk state c critical point D thermal diffusion H hydrodynamic i series summation index j series summation index P isobaric PE piston effect T isothermal, thermodynamic V isovolumetric
of pressure waves must be resolved [9], as demonstrated in early studies using a hydrodynamic model [5]. This particular model solved the NaviereStokes equations coupled with the van der Waals equation of state using the PISO method [17]. Improved numerical results at both acoustic and thermal diffusion time scales have been obtained from this model for single [18,19] and binary [20] supercritical fluids using the explicit Mac-Cormack method with the FCT algorithm [21] to minimize spurious numerical oscillations. The same method has also been applied to an alternative version of the above hydrodynamic model for supercritical fluids [22e24] that uses instead a real gas equation of state as well as characteristics-based boundary conditions [25]. This model was solved recently using a CrankeNicolson scheme for temporal discretization and a central differences scheme for spatial discretization [26], coupling velocity and
Please cite this article in press as: P.C. Teixeira, L.S.B. Alves, Modeling supercritical heat transfer in compressible fluids, International Journal of Thermal Sciences (2014), http://dx.doi.org/10.1016/j.ijthermalsci.2014.08.011
3
P.C. Teixeira, L.S.B. Alves / International Journal of Thermal Sciences xxx (2014) 1e12
pressure fields through the SIMPLEC algorithm [27]. These latter studies using real gas equations of state are able to capture the propagation of pressure waves with excellent numerical resolution and very good agreement with experimental data. An important limitation of the numerical methods utilized in these aforementioned studies is their practical restriction to short simulation times, on the order of the acoustic time scale [9]. The time span of this scale is a few micro seconds whereas the thermal diffusion time span is a few minutes in the same respective studies, effectively rendering simulations towards steady-state too cpu-time consuming even when restricted to one-dimensional domains. As far as the authors are aware, all numerical simulations using hydrodynamic models for heat transfer inside cavities containing supercritical fluids that have been reported so far in the literature overcome this difficulty with the same approach [9], known as either acoustic filtering [28] or low Mach number asymptotic equations [29]. It works by asymptotically expanding all dependent variables with respect to the Mach number, considered a perturbation parameter, and neglecting high-order terms in the flow governing equations to filter high frequency acoustic behavior, leaving only low frequency behavior such as the adiabatic compression of the bulk fluid. This is the reason why this approach cannot capture pulse like sound waves due to strong nonlinearities that might be generated by sudden boundary temperature increases [30]. Many nard convection of fluids near their simulations of RayleigheBe critical point have been performed with this approach, where the pressureevelocityedensity coupling was enforced using either SIMPLE [31,30], PISO [32e34] or SIMPLER [35e39] algorithms. Despite its quite common use in the literature, this approach has another significant drawback. The volumetric rate of change of a fluid particle depends on both temperature and pressure material derivatives. This latter dependence is proportional to the ratio between specific heats and the Mach number squared for problems with prescribed wall temperatures [40]. In fact, these two parameters have been used in the asymptotic expansion employed by earlier studies [32e34]. However, the ratio between specific heats diverges as the critical point is approached, requiring the Mach number to be even smaller for the overall dependence to be small. Furthermore, the speed of sound goes to zero in the same limit, which effectively increases the Mach number. Hence, more terms are required in the asymptotic expansion as the distance to the critical point decreases. Others have noted similar limitations of the asymptotic expansion approach as well [9]. This explains why a higher-order buoyancy term must be added to the momentum equations in order for this approach to satisfy the Schwarzschild criterion [32e34], which becomes dominant close to the critical point. At this stage, it is clear that numerical simulations of hydrodynamic models face an interesting dilemma. On one hand, time marching is restricted by acoustic time scales without acoustic filtering, rendering these simulations computationally quite expensive. On the other hand, acoustic filtering requires the inclusion of more terms in the low Mach asymptotic expansion as the critical point is approached because of the fluid property divergence phenomenon. Hence, an important question naturally arises: is it possible to bypass acoustic time scales in order to simulate both piston effect and thermal diffusion time scales without acoustic filtering? The present paper provides an answer to this issue using a low Mach number Preconditioned DensityBased Method [41e43]. It has an additional advantage over acoustic filtering: it uses no asymptotic expansions, allowing it to simulate arbitrary Mach number flows. In doing so, it should have no difficulties trying to predict the Schwarzschild criterion, since it has been quite capable of predicting density stratification in the
nard problem [44e46]. The authors believe classical RayleigheBe this is the first time such a method is applied to the hydrodynamic models discussed here. Numerical results are validated against a novel solution of the thermodynamic model obtained with the Generalized Integral Transform Technique. Once again, the authors are not aware of any other machine-precision accurate fully analytical solution of this model in the literature. Such a comparative analysis between hydrodynamic and thermodynamic models has been successfully employed in the past for onedimensional cavities in the absence of buoyancy forces [47,10], where filtering has little impact on solution accuracy at large time scales [9]. A similar comparison for the heat transfer within supercritical fluids is pursued in the present paper and discussed in detail next. 2. Hydrodynamic model 2.1. Compressible governing equations with real gas equation of state The one-dimensional governing equations are written in conservative form as
vQ vEv ¼ vx vt
vEi ; vx
(1)
where independent variables t and x represent the physical-time and spatial coordinate, respectively. Furthermore, the conservative dependent variable vector is defined as
(2)
Q ¼ fr; ru; rEg; and the inviscid and viscous flux vectors are defined as
n o Ei ¼ ru; ru2 þ P; ðrE þ PÞu and 4m vu 4m vu vT ; u þk ; ¼ 0; 3 vx 3 vx vx
Ev (3)
respectively. In the above definitions, r stands for density, u for velocity, E ¼ e þ u2/2 for total internal energy per unit mass, e for internal energy per unit mass, P for pressure, m for dynamic viscosity, k for thermal conductivity and T for temperature. The Newtonian fluid hypothesis is assumed to hold, i.e. bulk viscosity is neglected [10]. Relations (1) to (3) form a system of three equations and four unknowns. Hence, closure is necessary. This is achieved with the introduction of an equation of state. In the present study, the thermodynamic identity
dr ¼ rkT dP
(4)
raP dT;
is utilized for this purpose, where the isothermal compressibility and isobaric thermal expansion coefficient are defined as
kT ¼
1 vr r vP T
and
aP ¼
1 vr ; r vT P
(5)
respectively. These thermodynamic properties can be extracted directly from tables or correlations for real gases. Otherwise, they can be calculated through their definitions in (5) when algebraic equations of state are accurate enough. In other words, Eq. (4) represents a general relation for arbitrary gases. It has already been introduced in the context of hydrodynamic models for supercritical heat transfer [18e20] and its use is more common nowadays [22e24].
Please cite this article in press as: P.C. Teixeira, L.S.B. Alves, Modeling supercritical heat transfer in compressible fluids, International Journal of Thermal Sciences (2014), http://dx.doi.org/10.1016/j.ijthermalsci.2014.08.011
4
P.C. Teixeira, L.S.B. Alves / International Journal of Thermal Sciences xxx (2014) 1e12
2.2. Low Mach preconditioned Density-Based Method Any attempt to generate numerical solutions for the compressible flow Eqs. (1)e(4) is going to be quite difficult at small Mach numbers Ma, since both accuracy and efficiency deteriorate. A well known criterion for these problems to appear is Ma < 0.3, but it is strictly valid for isothermal flows of perfect gases only. Stiffness has long been identified as the main agent responsible for these losses. It manifests in two different ways. The first one is related to time step restrictions. They are caused by discrepancies between acoustic and convective characteristic time scales. Preconditioning is a popular way of minimizing this problem [41e43], since it forces all eigenvalues to have the same order of magnitude. On the other hand, the second one is attributed to cancelation errors that increase as the Mach number decreases. They amplify round-off errors due to discrepancies between thermodynamic and hydrodynamic contributions to primitive and conservative flow variables. Hence, this problem is alleviated by introducing known reference states and solving for the remaining fluctuations [48e50], leading to reduced cancelation errors for pressure and temperature residuals as well [51,52]. As far as the authors are aware, this is the first time such a methodology is applied to the heat transfer problem in a cavity containing a supercritical fluid. In the present work, it is sufficient to control only pressure cancelation errors. This is achieved by re-writing Eq. (1) as
b vQ vEv ¼ T vt vx
vEi vx
H
dPT ; dt
vr B vP B B vQ vr B ¼B T¼ u b B vP vQ B @ vr vh þr H vP vP
0 r 1
ru
for real gases, where CP stands for the specific heat at constant pressure and can be also extracted from tables or correlations. Finally, source term H is simply the first column of matrix T, i.e. its coefficients are obtained from relation Hi ¼ Ti,1. Nevertheless, neither Eq. (1) nor Eq. (6) can be solved in the limit of very small Mach numbers without further modifications. It is necessary to precondition these equations to force their eigenvalues to have the same order of magnitude at low Mach numbers. This leads to re-writing Eq. (6) as
G
b vQ vEv ¼ vt vx
vEi vx
H
dPT ; dt
(11)
where the low Mach preconditioning matrix is based on T and defined as
0
rP urP HrP þ rhP
0 r ru
1 rT A; urT HrT þ rhT
(7)
(8)
1 vr C vT C C vr C u C; C vT C vr vh A H þr vT vT
(10)
by modifying both density and enthalpy dependencies on both pressure (rP and hP) and temperature (rT and hT). These new parameters are selected to modify the eigenvalues of the precondib where tioned inviscid flux Jacobian c M ¼ G 1 A,
which is possible because PT was assumed spatially uniform. It is important to emphasize that such a pressure splitting does not introduce additional approximations because not a single term in the governing equations is neglected on the basis that the Mach number is small, which is an approximation at the core of the acoustic filtering procedure. Hence, the present method is applicable to flows with arbitrary Mach numbers whereas acoustic filtering can only be used for small Mach number flows. Although b , temperature both pressure and velocities are often present in Q can be replaced by enthalpy [53] or entropy [54], among others. The conservative to primitive variable Jacobian is defined as
0
vh 1 TaP ; ¼ r vP T
and
(6)
where pressure is decomposed into P ¼ PT þ PH, which is the sum of its hydrodynamic (PH ~ ru2) and thermodynamic (PT ~ rc2) contributions, with c being defined as the sound speed. The latter is a zero-dimensional, but possibly time-dependent, variable and orders of magnitude higher than the former, when the Mach number is very small. This step is essential for controlling pressure cancelation errors in the momentum conservation equations, where P must be replaced by PH. In other words,
n o Ei ¼ ru; ru2 þ PH ; ðrE þ PÞu ;
vh ¼ CP vT P
G¼@
and solving for the primitive dependent variable vector
b ¼ fP ; u; Tg; Q H
where H ¼ h þ u2/2 stands for total enthalpy per unit mass and h ¼ e þ P/r for enthalpy per unit mass. Density dependencies on pressure and temperature are extracted from the isothermal compressibility and isobaric thermal expansion coefficient according to relations in (4), respectively, whereas enthalpy dependencies are obtained from
(9)
1
0
vr u B vP B B vr b ¼ vEi ¼ B 1 þ u2 A B vP b B vQ B @ vr vh u H þr vP vP
r 2ru r u2 þ H
(12)
1 vr C vT C C vr C 2 u C; vT C C vr vh A þr u H vT vT u
(13)
is the inviscid flux Jacobian based on primitive variables. They are given by
rP ¼
1
2 b c vh ; ¼ vT
rT ð1
rhP Þ=ðrhT Þ; rT ¼
vr vh ; hP ¼ vT vP
and
hT
(14) where the preconditioned sound speed is defined as
h hpffiffiffiffiffiffi i i b c ¼ min max u2 ; ε ; c ;
(15)
in order to maintain the maximum eigenvalue ratio near unity. The parameter ε, defined elsewhere [43], is employed to minimize eigenvector orthogonality loss [55]. As a direct consequence of preconditioning, source term H is now given by the first column of matrix G, i.e. its coefficients are obtained from relation Hi ¼ Gi,1. Finally, the eigenvalue modification introduced by low Mach preconditioning changes the natural unsteady evolution of Eq. (11),
Please cite this article in press as: P.C. Teixeira, L.S.B. Alves, Modeling supercritical heat transfer in compressible fluids, International Journal of Thermal Sciences (2014), http://dx.doi.org/10.1016/j.ijthermalsci.2014.08.011
5
P.C. Teixeira, L.S.B. Alves / International Journal of Thermal Sciences xxx (2014) 1e12
which occurs now under pseudo-time t. An unsteady simulation can still be performed with dual-time-stepping, i.e.
G
b vQ vQ vEv þ ¼ vt vt vx
vEi vx
H
dPT ; dt
(16)
where pseudo-time iterations must be performed until pseudotime steady-state, defined through an user prescribed tolerance, at each physical-time step to recover Eq. (1). 2.3. Adiabatic compression mechanism Cancelation errors required a pressure decomposition to be minimized, which led to the introduction of PT as an additional dependent variable. There are different ways to obtain this thermodynamic pressure contribution. Since PT is a spatially homogeneous variable, a quite common procedure enforces one of the conservation laws in integral form. For instance, flows completely bounded by impermeable walls require total mass conservation through time. In other words,
ZL
rðx; tÞdx ¼
0
ZL
rðx; t ¼ 0Þdx ¼ const;
(17)
0
where L stands for the one-dimensional cavity length. Relation (17) can be re-written in a more appropriate form in the present paper, since the total differential expressed in (4) can be divided by dt to become
vr vP ¼ rkT vt vt
raP
vT ; vt
(18)
and integrated over the entire cavity length to yield
0
dPT @ ¼ dt
ZL 0
raP
vT dx vt
ZL 0
rkT
1,
vPH A dx vt
ZL
order Simpson's rule. Boundary conditions were defined as follows: i) cavity walls are assumed impermeable, ii) the right wall temperature is maintained at its initial value whereas its left counterpart is prescribed at a higher value at t > 0þ in order to induce an impulsive heating, iii) the pressure spatial derivative is assumed zero at both walls, iv) with their respective densities extracted from equation of state (4). Such a boundary condition for pressure, although an approximate one, is consistent with the spatial homogeneity of the thermodynamic pressure contribution, which is dominant compared to the hydrodynamic contribution. This code was originally developed for unsteady threedimensional simulations of transverse jets, and it was validated for two-dimensional planar mixing-layers and three-dimensional free jets as well [56]. Furthermore, it was recently utilized to generate steady-states for convectively unstable planar-mixing layers [57].
3. Thermodynamic model 3.1. Heat equation with adiabatic compression source term In order to validate the hydrodynamic model proposed in the previous section for the heat transfer inside a one-dimensional cavity filled with supercritical fluid under zero gravity, a comparison with the thermodynamic model [4,3] is performed. This model neglects the fluid flow and considers an energy balance in the form of
rCP
vT vt
CV Þ
rðCP
(20)
where pressure, assumed spatially uniform, is obtained by simplifying Eq. (19) to
dP ¼ dt
ZL
raP
vT dx vt
0
rkT dx;
kT dP v vT ¼ k ; vx vx aP dt
, ZL
rkT dx;
(21)
0
(19)
0
which is employed as the additional closure equation for this hydrodynamic model. It is worth noting that the second term within parenthesis in Eq. (19) can be neglected based on the assumption that PT [ PH, valid when PH/PT ~ O(gMa2), where g ¼ CP/CV stands for the ratio between specific heats and CV for the heat capacity at constant volume. 2.4. Temporal and spatial resolution Eq. (16) is marched in physical-time with the second-order Lstable implicit backwards difference formula (also known as BDF) scheme. On the other hand, pseudo-time integration is performed with the first-order L-stable implicit Euler scheme [43]. Accurate unsteady simulations are achieved at each physical-time step by requiring a six orders of magnitude reduction in the maximum b ¼Q b pþ1 Q b p , where p is the pseudo-time solution increment d Q pseudo-time iteration counter. As discussed and justified in the next section, fluid properties are assumed constant. Hence, Eq. (19) can be integrated analytically in time. Initial conditions assume a motionless fluid at a fixed and uniform temperature and pressure. Respective initial densities are extracted from equation of state (4). Viscous fluxes in Eq. (16) are spatially resolved with a secondorder conservative central-difference scheme. On the other hand, inviscid fluxes in this same equation are spatially resolved with a third-order preconditioned flux-difference scheme [43]. All the numerical integrations in Eq. (19) are performed using the fourth-
which represents a mass balance inside the cavity. Hence, an integraledifferential equation can be constructed for this model by combining Eqs. (20) and (21). If all fluid properties are assumed constant and evaluated at the initial state, this model takes the form
vT vt
1
1 g0
dTb v2 T ¼ a0 2 ; dt vx
(22)
where the thermal diffusivity is defined as a0 ¼ k0/(r0CP0 ) and the bulk temperature of the fluid is calculated using definition
Tb ðt Þ ¼
1 L
ZL
(23)
T ðx; t Þdx;
0
and the initial boundary conditions for temperature are
Tðx; t ¼ 0Þ ¼ T0 ; Tðx ¼ 0; tÞ ¼ T1
and
Tðx ¼ L; tÞ ¼ T0 ;
(24)
respectively. Density can now be obtained from an equation of state. In this particular case where all properties are assumed constant, integrating Eq. (4) yields
r ¼ 1 þ kT0 ðP r0
P0 Þ
aP0 ðT
T0 Þ;
(25)
where P0 and r0 are the initial pressure and density.
Please cite this article in press as: P.C. Teixeira, L.S.B. Alves, Modeling supercritical heat transfer in compressible fluids, International Journal of Thermal Sciences (2014), http://dx.doi.org/10.1016/j.ijthermalsci.2014.08.011
6
P.C. Teixeira, L.S.B. Alves / International Journal of Thermal Sciences xxx (2014) 1e12
3.2. Integral transformation A solution for the thermodynamic model described in the previous subsection can be obtained with the Generalized Integral Transform Technique [58,59]. This technique is a natural extension of the Classical Integral Transform Technique, born out of Separation of Variables [60], required in the presence of otherwise nontransformable terms. It has been successfully used in the past by the author to simulate natural convection in porous media [61], nard instanonlinear convectionediffusion flows [62], DarcyeBe bility [63] and Prats flow of a power-law fluid [64]. Recently, Eqs. (20) and (21), subjected to initial and boundary conditions (24), were solved with this technique for the particular case where all fluid properties were allowed to vary linearly with both pressure and temperature [15]. Since this nonlinearity was shown to have very little impact on the piston effect, fluid properties are going to be considered constant in this work. Furthermore, it is shown for the first time here that such a simplification allows the derivation of a novel fully analytic and machine-precision accurate solution for this linear thermodynamic model. Before deriving such a solution, it is useful to re-write the problem in dimensionless form using the following transformation
t x t¼ ; x¼ tD L
and
T Q¼ T1
T0 ; T0
(26)
where tD ¼ L2/a0 stands for the thermal diffusion relation time. Eq. (22) can be re-written for the dimensionless temperature as
vQ vt
1 g0
1
Z1
(27)
(28)
Qðx; tÞdx;
0
Qðx; t ¼ 0Þ ¼ 0; Qðx ¼ 0; tÞ ¼ 1
and
Qðx ¼ 1; tÞ ¼ 0:
Qðx; tÞ ¼ F ðxÞ þ qðx; tÞ;
(31)
x;
and q stands for the homogeneous dimensionless temperature, which is the new dependent variable of the model. Substituting Eq. (30) into Eq. (27) yields
1
(33)
respectively. Doing the same to initial and boundary conditions (29) leads to
qðx; t ¼ 0Þ ¼ x
1
and
qðx ¼ 0; tÞ ¼ qðx ¼ 1; tÞ ¼ 0:
(34)
Now, spatial and temporal dependences of the model problem unsteady representation are separated by proposing a solution to Eq. (32) in the form
qðx; tÞ ¼
∞ X i¼1
~ i ðxÞqi ðtÞ; j
(35)
where the eigenfunction is provided by eigensystem
d2 ji dx2
þ b2i ji ðxÞ ¼ 0;
(36)
with boundary conditions
ji ð0Þ ¼ 0
and
ji ð1Þ ¼ 0;
(37)
which yield the eigenfunctions and eigenvalues
and
(38)
bi ¼ ip;
respectively, where i ¼ 1, 2,…,∞. Since these eigenfunctions are orthogonal, the integral transformed homogeneous temperature can be defined according to relation
qi ðtÞ ¼
Z1
~ i ðxÞqðx; tÞdx; j
(39)
1 g0
dqb v2 q ¼ 2; dt vx
based on the above eigensystem, where the norm
Ni ¼
Z1
ji ðxÞ2 dx ¼
1 ; 2
(40)
0
pffiffiffiffiffi ~ i ¼ ji = Ni . is used by the normalized eigenfunctions j ~ i, integrating the result over the Multiplying Eq. (32) by j dimensionless cavity length and applying transformation (39) to the time derivative term yields
(30)
where the filter is defined as the steady-state temperature distribution
qðx; tÞdx;
0
(29)
An important first step in the solution procedure is to make the boundary conditions homogeneous before initiating the integral transformation procedure. Quite significant improvements in summation series convergence rates are achieved due to this step. It can be accomplished by introducing a filter for the dimensionless temperature, given by
vq vt
qb ðtÞ ¼
Z1
0
with dimensionless initial and boundary conditions written as
F ðxÞ ¼ 1
and
ji ðxÞ ¼ sin½bi x
dQb v2 Q ¼ 2 ; dt vx
where the dimensionless bulk temperature is given by
Qb ðtÞ ¼
1 Fb ¼ 2
dqi dt
hi 1
and into Eq. (28) yields the bulk filter and dimensionless temperatures
dqb ¼ dt
Z1
2
~ i ðxÞ v qdx; j vx2
(41)
0
where integral transform coefficient
hi ¼
Z1 0
(32)
1 g0
~ i ðxÞdx ¼ j
pffiffiffi1 2
cos½bi ; bi
(42)
is defined for simplicity. Integrating the r.h.s. of Eq. (41) by parts, applying both boundary conditions in (34) as well as (37), substituting Eq. (36) and then applying integral transform definition (39) into the result yields
Please cite this article in press as: P.C. Teixeira, L.S.B. Alves, Modeling supercritical heat transfer in compressible fluids, International Journal of Thermal Sciences (2014), http://dx.doi.org/10.1016/j.ijthermalsci.2014.08.011
7
P.C. Teixeira, L.S.B. Alves / International Journal of Thermal Sciences xxx (2014) 1e12
Z1
2
~ i ðxÞ v qdx ¼ j vx2
0
Z1 0
~i d2 j dx2
qðx;tÞdx ¼
b2i
Z1
~ i ðxÞqðx;tÞdx ¼ b2 qi ðtÞ; j i
which can be substituted into Eq. (41), using inverse definition (35) in the non-transformed bulk temperature term, to generate
1
∞ dqj 1 X hi hj þ b2i qi ðtÞ ¼ 0; g0 dt
(44)
j¼1
for i ¼ 1, 2,…,∞, or, in a more compact form, ∞ X j¼1
Ai;j
xþ
0
(43)
dqi dt
Qðx; tÞ ¼ 1
dqj þ b2i qi ðtÞ ¼ 0 dt
where Ai;j ¼ di;j
1
1 hh; g0 i j
N X i¼1
~ i ðxÞqi ðtÞ; j
(50)
with its also analytical bulk value, defined in (28), yielding
Qb ðtÞ ¼ F b þ qb ðtÞ ¼
N 1 X þ h q ðtÞ; 2 i¼1 i i
(51)
where both expressions can be machine-precision accurate for a high enough N. Its value was chosen to impose a temperature absolute error smaller than six orders of magnitude in the present calculations for each test case presented here and discussed next. 4. Results and discussion 4.1. Fluid property selection
(45)
is the integral transform matrix coefficient, with di,j representing the Kronecker delta. This integral transformed equation is subject to the transformed initial condition, obtained by applying a similar transformation procedure to the initial condition in Eq. (34),
qi ð0Þ ¼
Z1
~ i ðxÞdx ¼ 1Þj
ðx
0
pffiffiffi 2 : bi
(46)
Five test cases with g ¼ 2, 5, 10, 15 and 20 were employed to demonstrate the low Mach Preconditioned Density-Based Method capability to simulate zero gravity heat transfer within supercritical fluids. In all cases the fluid is assumed initially quiescent until the left wall temperature is raised to T1 ¼ T0 þ 10 mK. CO2 was chosen as a test fluid and its properties were obtained from the NIST database at P0 ¼ 7.4 MPa with T0 < Tc for all 5 test cases, where its critical point is located at Pc ¼ 7.3773 MPa, Tc ¼ 304.1282 K and rc ¼ 467.6 kg/m3. These properties are shown in Table 1. Thermodynamic consistency was verified by comparing NIST data for CV0 and c0 at these conditions to their respective values obtained from equation
3.3. Temporal and spatial resolution In order to obtain a solution for the system of Eqs. (45) and (46), the summation series must be truncated. Since Eqs. (36)e(38) represent an eigensystem from the SturmeLiouville class, the series solution proposed in Eq. (35) is convergent. Hence, both temporal and spatial accuracies of this proposed solution for the thermodynamic model are controlled entirely by the number of terms N it utilizes, as long as the error introduced by the numerical integration of Eqs. (45) and (46) is kept small enough. This additional error can be avoided altogether by inverting matrix A in Eq. (45), leading to N dqi X þ Bi;k qk ðtÞ ¼ 0 dt k¼1
where Bi;k ¼ Ai;j1 dj;k b2k ;
(47)
is a new integral transform matrix coefficient. Matrix B also has an analytical expression, since the inverse matrix A 1 exists and is given by
Ai;j1 ¼ di;j þ
ðg0 g0
ðg0
1Þhi hj ; P 2 1Þ N k¼1 hk
(48)
noting that both matrixes become diagonal when g0 ¼ 1, as expected for pure thermal diffusion problems. Eq. (47), subjected to initial condition (46), has an analytical solution written here in vector form as
CV0 ¼ CP0 þ
qðtÞ ¼ qð0Þe ;
(49)
based on the matrix exponential method [65]. Combining the above solution with Eqs. (30), (31) and (35) generates the analytical solution for the dimensionless temperature
(52)
which is Mayer's generalized relation, and
vr vr 1 c20 ¼ 1 þ vP 0 vT 0
r0
vh vh r ; 0 vP 0 vT 0
(53)
which is an alternative definition for the sound speed [43]. These two consistency checks are vital for an adequate comparison between hydrodynamic and thermodynamic models because neither CV0 nor c0 are provided as input data for the compressible flow code used to simulate the hydrodynamic model. If the values of CV0 taken from NIST and calculated from Eq. (52) are different, it is not possible to guarantee that simulations with both models use the same value of g. Since Mayer's relation is derived from an entropy balance, such errors represent either generation or destruction of entropy. Furthermore, if the values of c0 taken from NIST and calculated from Eq. (53) are different, the adiabatic compression
Table 1 CO2 supercritical fluid properties at P0 ¼ 7.4 MPa for all 5 test cases. Units are [K], [kg/ m3], [J/(kg K)], [Pa/s], [W/(m K)], [kg/(m3 Pa)], [kg/(m3 K)], [Pa/K] and [J/kg], respectively.
g T0
r0 CP0 10
Bt
T0 vr vP ; r20 vT 0 vT 0
m0 106
3
k0 103 vr/vPj0 106 vr/vTj0 vP/vTj0 10 7 h0 10 3
2
5
10
15
20
227.751653 1153.24501 1.93354539 224.773388 171.870551 2.27999999 3.58772614 15.7357075 103.563563
300.674704 716.592550 5.14963289 57.9506113 79.7014590 56.1890000 19.8831281 3.53858608 277.542995
303.396869 635.674676 11.5249917 48.1881349 80.1302023 195.920000 52.0253582 2.65540461 296.872420
303.849596 605.776395 18.6016283 45.0487798 84.8271279 366.299999 87.6389969 2.39255367 303.351139
304.012817 588.786850 26.3001728 43.3595669 89.5621628 560.300001 126.345769 2.25500021 306.922011
Please cite this article in press as: P.C. Teixeira, L.S.B. Alves, Modeling supercritical heat transfer in compressible fluids, International Journal of Thermal Sciences (2014), http://dx.doi.org/10.1016/j.ijthermalsci.2014.08.011
8
P.C. Teixeira, L.S.B. Alves / International Journal of Thermal Sciences xxx (2014) 1e12
mechanism implicit in each model will not be the same because they are a consequence of the propagation and reflection of acoustic waves. As a result, the bulk temperature increase induced by the piston effect in each model will be different. The use of Eq. (25), which is the constant property version of Eq. (4), in conjunction with the experimental data presented in Table 1 guarantees between four and six digits of accuracy for CV0 and c0 when comparing Eqs. (52) and (53) with their respective NIST data. This comparison is shown in Table 2, which includes relative errors for both properties.
10
2
10
3
10
4
10
5
1 5
2
2
4.2. Verification analysis Before the solutions obtained for the hydrodynamic and thermodynamic models can be compared, their numerical behavior must be first verified against theoretical predictions based on the methodologies used to derive them. This verification step is shown in Figs. (1) and (2) for the spatial and temporal resolution of the numerical schemes implemented with the low Mach Preconditioned Density-Based Method, respectively. They show the dimensionless temperature absolute error versus dimensionless grid size and time step, respectively, for g ¼ 2 and 5. Their slope represents the actual order measured in the numerical results whereas dashed lines in both figures represent the theoretical slopes of orders 1 and 2. As one would expect, low resolution, i.e. large grid sizes and time steps, decreases actual order to 1. Furthermore, actual order tends to the theoretical value of 2 as resolution is increased in both figures. Another feature is worth mentioning in these figures. Error increases as g increases from 2 to 5. This is also expected, since a higher g value leads to stronger temporal gradients as well as a thinner thermal boundary-layer thickness at early simulation times. These phenomena require improved resolution in time and space, respectively, to contain error propagation. In order to avoid it, a non-uniform grid spacing is employed with grid points clustered near both walls and a minimum of 100 time steps are taken before reporting a solution early in the piston effect regime. Once both temporal and spatial resolution requirements have been met, convergence rates can be properly analyzed. Fig. 3 shows the number of pseudo-time steps per physical-time step required to reduce the maximum dimensionless temperature pseudo-time increment in six orders of magnitude for g ¼ 2, 5, 10, 15 and 20. Temperature is chosen instead of pressure and velocity because, among all three dependent variables, it yields the largest dimensionless increments. Three features are clear in this figure. First, the number of pseudo-time steps is higher for n ¼ 1 in all cases. This problem is caused by the discontinuity between initial and left boundary conditions. It models the infinitely fast boundary heating imposed at t ¼ 0þ when the left wall temperature is raised from T0 to T1. Nevertheless, this number decreases drastically for n ¼ 2 in all cases. From this point onwards, the number of pseudo-time
10
3
2 10
3
5 10
3
10
2
2 10
2
Fig. 1. Dimensionless temperature absolute error DQ versus dimensionless grid size Dx for g ¼ 2 and 5, both measured at x ¼ 0.1 and t ¼ 0.05. Spatial grid has a uniform distribution with Nx ¼ 51, 101, 201, 401, 801 and 1601 points whereas a constant time step with NT ¼ 20,000 points per piston effect relaxation time was employed. Dashed lines indicate theoretical order behavior.
iterations slowly decreases towards 1 as physical-time steady-state is approached. This is the second important feature implied in this figure, since it only shows results up to n ¼ 50. Finally, the number of pseudo-time steps required for convergence increases with g. This third feature is better illustrated in Fig. 4, which presents the dimensionless temperature increment convergence in pseudo-time at n ¼ 2 for the same five cases with g ¼ 2, 5, 10, 15 and 20. It was generated using the optimal pseudo-time CFL number for each case. The convergence slow down is caused by the stronger adiabatic heating caused by higher g values, which induces stronger temporal gradients at early simulation times due to a decreasing piston effect relaxation time as well as stronger thermal boundarylayer spatial gradients that are long lasting due to an increasing thermal diffusion relaxation time. Finally, having verified the numerical behavior of the hydrodynamic model numerical solution, the same must be done to the thermodynamic model analytic solution. The convergence characteristics of the integral transform series (51) are evaluated in Fig. 5 for g ¼ 2, 5, 10, 15 and 20. This figure shows the dimensionless bulk temperature absolute error as a function of the number of terms utilized in the summation series, which is the only parameter controlling accuracy. All temperature data is measured after 100
3 10
4
1 10
5
4
2 Table 2 Comparison between NIST data and Eq. (52) for specific heat at constant volume CV0 in [J/(kg K)] as well as Eq. (53) for sound speed c0 in [m/s], respectively, obtained using the conditions specified in Table 1. jR.E.j means absolute value of the relative error.
g
2
CVNIST 0 Eq:ð52Þ C V0
5
10
15
966.772706 1029.92485 1152.50734 1240.10180 1315.01809 966.770313 1029.92152 1152.40540 1239.89589 1314.96426
5
10
5
2
20
jR.E.j 106 2.47488219 3.23624407 88.4520310 166.039205 40.9347205 cNIST 936.587859 298.303820 225.921799 202.362479 188.933769 0 Eq:ð53Þ 936.586623 298.309926 225.945398 202.374573 188.907927 c 0
3 10
jR.E.j 106 1.31878094 20.4679299 104.456944 59.7659034 136.779605
5 10
5
10
4
2 10
4
5 10
4
10
3
Fig. 2. Dimensionless temperature absolute error DQ versus dimensionless time step Dt for g ¼ 2 and 5, both measured at x ¼ 0.1 and t¼0.0015. Constant time steps were
obtained with NT ¼ 625, 1250, 2500, 5000 and 10,000 points per piston effect relaxation time whereas a uniform spatial grid with Nx ¼ 1601 was employed. Dashed lines indicate theoretical order behavior.
Please cite this article in press as: P.C. Teixeira, L.S.B. Alves, Modeling supercritical heat transfer in compressible fluids, International Journal of Thermal Sciences (2014), http://dx.doi.org/10.1016/j.ijthermalsci.2014.08.011
9
P.C. Teixeira, L.S.B. Alves / International Journal of Thermal Sciences xxx (2014) 1e12
p
b
0.1
350 300
0.01 250
0.001
200 150 100
20 15 10 5 2
50 0 10
20
30
40
50
n
Fig. 3. Number of pseudo-time iterations p per physical-time step n for g ¼ 2, 5, 10, 15 and 20. All simulations were performed with Nx ¼ 801 and NT ¼ 10,000 points per piston effect relaxation time.
time steps or, in other words, 1% of each respective piston effect relaxation time. Convergence patterns are qualitatively similar in all cases, but error increases with g as expected. Nevertheless, such an error is smaller than 10 4 when N ¼ 450. Since the steady filter (31) is employed, accuracy increases with time. Hence, bulk temperature error will be smaller than 10 4. 4.3. Validation analysis Now that all verification steps for both numerical and analytical solutions are finalized, a comparative analysis between them can be undertaken to validate the numerical solution obtained with the low Mach preconditioned Density-Based Method. This is performed in Fig. 6, which shows the dimensionless temperature distribution over the dimensionless distance at three different dimensionless times for g ¼ 2, 5, 10, 15 and 20. Lines represent the thermodynamic model analytic solution given by Eq. (50) whereas all points represent the hydrodynamic model numerical solution of Eq. (16), but only a few data points are shown to facilitate visualization. Circular points, and their respective lines, are solutions that include the piston effect model from Eq. (19). On the other hand, square points, and their respective lines, are solutions that include thermal diffusion only, since they were generated by imposing PT ¼ P0 instead of Eq. (19). There is an excellent agreement between
10
5
10
6
10
7
10
8
10
9
10
10
10
11
2 0
20
5 40
10 60
80
15 100
20 120
140
p
Fig. 4. Maximum dimensionless temperature pseudo-time increment dQ per pseudotime step p for g ¼ 2, 5, 10, 15 and 20. Data extracted from physical-time step n ¼ 2 in Fig. 3.
10
4
10
5
20 15 10 5 2 0
100
200
300
400
N 500
Fig. 5. Dimensionless bulk temperature absolute error DQb versus number terms N used in the summation series solution, measured at 1% of the piston effect relaxation time for g ¼ 2, 5, 10, 15 and 20.
both solutions in all cases, demonstrating that low Mach preconditioning is a viable option to model the piston effect phenomenon. In previous comparisons using a cubic equations of state [10], there is good agreement between both model solutions but non-negligible differences can still be observed. The source of this discrepancy seems to be the approximate algebraic equation of state employed instead of Eq. (25), coupled with the fact that Eq. (52) was verified but Eq. (53) was not. When we attempted to do the same in our study, similar discrepancies appeared between both model solutions. Another important aspect of the low Mach preconditioning method employed here to solve the hydrodynamic model is the inclusion of source term H in Eq. (16). This term is multiplied by a pseudo-time derivative, which effectively disappears once pseudotime steady-state is reached within each physical-time step. Hence, one might be tempted to exclude this term altogether from this hydrodynamic model solution methodology. However, its presence is essential for the generation of a correct unsteady solution. Without this term, the march in pseudo-time is not capable of leading the previous physical-time step solution towards its next correct physical-time step. The cumulative effect of this error over physical-time leads to an unphysical bulk temperature evolution. There are a few important differences between the plots presented in Fig. 6 that are worth discussing. Although all plots used Nx ¼ 801 non-uniformly distributed grid points, only the plots for g ¼ 2 and 5 were generated using NT ¼ 10,000. Both model solutions are indistinguishable for all three times shown in this figure when g ¼ 2. However, a small deviation can be observed at the later time between both model solutions including the piston effect when g ¼ 5. The stronger early temporal gradients in this latter case lead to bigger errors, which cumulate over time and become graphically discernible at t ¼ 0.2. This effect increases with g. In order to avoid it, the numerical data shown for g ¼ 10, 15 and 20 was generated using NT ¼ 20,000, 30,000 and 40,000, respectively. An alternative solution to this problem would employ higher-order schemes in physical-time instead, which would decrease the number of necessary physical-time steps. Nevertheless, no discernible differences can be observed between both model solutions for these latter three cases, where g ¼ 10, 15 and 20. Fig. 6 indicates that both strong thermal boundary-layer and fast bulk temperature increase are captured by the hydrodynamic model numerical solution provided by the low Mach Preconditioned Density-Based Method. As a final comparison, bulk temperature evolution in time is shown in Fig. 7 for g ¼ 2, 5, 10, 15 and 20. Points show the hydrodynamic model numerical solution
Please cite this article in press as: P.C. Teixeira, L.S.B. Alves, Modeling supercritical heat transfer in compressible fluids, International Journal of Thermal Sciences (2014), http://dx.doi.org/10.1016/j.ijthermalsci.2014.08.011
10
P.C. Teixeira, L.S.B. Alves / International Journal of Thermal Sciences xxx (2014) 1e12
Fig. 6. Spatial distribution of the dimensionless bulk temperature Q at three different dimensionless times for g ¼ 2, 5, 10, 15 and 20. Points represent the hydrodynamic model numerical solution whereas lines represent the thermodynamic model analytical solution. Circles include the piston effect model whereas squares do not.
whereas solid lines show the thermodynamic model analytic solution. Furthermore, dashed lines indicate the location of the piston effect relaxation time. Since this characteristic time represents the time it takes the bulk temperature to relax towards steady-state, its dimensionless value tPE can be estimated from either solution using a standard boundary-layer theory approximation [66]
Qb ðtPE Þx0:99F b ;
(54)
since F b is the dimensionless bulk temperature steady-state. Using Eq. (51), it is easy to prove that Eq. (54) is equivalent to 2qb(tPE) x 0.01, which is plotted instead in Fig. 7. This figure shows that tPE decreases from its thermal diffusion value of tPE ¼ tD x 0.44532 as g increases from its thermal diffusion value of g ¼ 1, as expected. In fact, its values are tPE x 0.25399, 0.13736, 0.093965, 0.076385 and 0.065860 when g ¼ 2, 5, 10, 15 and 20, respectively. The classical theoretical expression for the piston effect relaxation time, derived using an approximate Fourier transform procedure [3], yields tPE x 1, 0.062500, 0.012346, 0.0051020 and 0.0027701 for the same values of g. Hence, the theoretical expression appears to underestimate the piston effect relaxation time whenever g > 2. Mid point temperature data from a previous study obtained while imposing a prescribed heat flux at the hotter wall [10], instead of the prescribed temperature utilized here, indicate the same trend. This means such an underestimation by
the theoretical expression is not due to the use of alternative boundary conditions. 5. Conclusions and future work This study presented simulations of thermodynamic as well as hydrodynamic models for the heat transfer taking place in a onedimensional cavity filled with supercritical fluid under zero gravity in order to capture the piston effect induced by differentially heated walls. All major conclusions reached here are discussed in detail below: When all fluid properties are considered constant, it is possible to derive an arbitrary precision analytical solution of the thermodynamic model. This is done by combing the Generalized Integral Transform Technique with the Matrix Exponential Method. Accuracy is controlled by a single parameter, which represents the number of terms employed in the convergent summation series solution. The authors believe this is the first fully analytical, essentially machine-precision accurate, solution of the thermodynamic model of this problem reported in the literature. It is possible to generate a numerical solution of the hydrodynamic model without acoustic filtering but also without having to resolve acoustic scales. In order to do so, a low Mach number
Please cite this article in press as: P.C. Teixeira, L.S.B. Alves, Modeling supercritical heat transfer in compressible fluids, International Journal of Thermal Sciences (2014), http://dx.doi.org/10.1016/j.ijthermalsci.2014.08.011
P.C. Teixeira, L.S.B. Alves / International Journal of Thermal Sciences xxx (2014) 1e12
Fig. 7. Dimensionless homogeneous bulk temperature qb normalized by dimensionless bulk filter F b versus dimensionless time t for g ¼ 2, 5, 10, 15 and 20. Points represent the hydrodynamic model numerical solution whereas lines represent the thermodynamic model analytical solution. Dashed lines and respective points locate the values of the dimensionless piston effect characteristic time tPE obtained from Eq. (54).
Preconditioned Density-Based Method is employed instead. This can only be achieved when the thermodynamic pressure evolution in pseudo-time is properly included in the simulation. As far as the authors are aware, this is the first application of a low Mach preconditioning method to the hydrodynamic model of this problem reported in the literature. The selection of an equation of state is dependent on the imposed fluid properties. They must satisfy both Mayer's generalized relation and the sound speed definition. It is only in this scenario that thermodynamic and hydrodynamic models yield the same results. Otherwise, the entropy balance as well as the adiabatic compression mechanism are not equally accounted for by both models, leading to discrepancies between their respective results. The classical theoretical expression for the piston effect relaxation time appears to underestimate the piston effect relaxation time estimated from a simulation of the same thermodynamic model used to derive this expression. Combining data from an earlier study [10] with the present data, it is possible to state that this feature is not an artifact of boundary condition choice. Current efforts aim at investigating the effects of different boundary conditions on the piston effect dynamics using both models. Furthermore, nonlinear variations are being investigated by accounting for property dependences on both pressure and temperature. Acknowledgments The authors would also like to acknowledge the financial support from the United States Air Force, through the AFOSReSOARD grant FA9550-13-1-0178, as well CNPq, through grant MCT/ CNPq14-481072/2012-8. References [1] D. Dahl, M.R. Moldover, Thermal relaxation near the critical point, Phys. Rev. A 6 (1972) 1915e1920. [2] K. Nitsche, J. Straub, The critical hump of cv under microgravity, results from d-spacelab experiment, in: Proceedings of the 6th European Symposium on Material Sciences under Microgravity Conditions, ESA SP-256. [3] A. Onuki, H. Hao, R.A. Ferrell, Fast adiabatic equilibration in a singlecomponent fluid near the liquid-vapor critical point, Phys. Rev. A 41 (1990) 2256e2259.
11
[4] H. Boukari, J.N. Shaumeyer, M.E. Briggs, R.W. Gammon, Critical speeding up in pure fluids, Phys. Rev. A 41 (1990) 2260e2263. [5] B. Zappoli, D. Bailly, Y. Garrabos, B.L. Neindre, P. Guenoun, D. Beysens, Anomalous heat transport by the piston effect in supercritical fluids under zero gravity, Phys. Rev. A 41 (1990) 2264e2267. [6] B. Zappoli, Near-critical fluid hydrodynamics, C. R. Mec. 331 (2003) 713e726. [7] M. Barmatz, I. Hahn, J.A. Lipa, R.V. Duncan, Critical phenomena in microgravity: past, present and future, Rev. Mod. Phys. 79 (2007) 1e52. s, A brief review of the thermophysical properties of supercritical [8] P. Carle fluids, J. Supercrit. Fluids 53 (2010) 2e11. [9] B. Shen, P. Zhang, An overview of heat transfer near the liquidegas critical point under the influence of the piston effect: phenomena and theory, Int. J. Therm. Sci. 71 (2013) 1e19. [10] V.S. Nikolayev, A. Dejoan, Y. Garrabos, D. Beysens, Fast heat transfer calculations in supercritical fluids versus hydrodynamic approach, Phys. Rev. E 67 (2003) 061202-1e061202-11. [11] A. Onuki, R.A. Ferrell, Adiabatic heating effect near the gaseliquid critical point, Phys. A 164 (1990) 245e264. [12] R.A. Ferrell, H. Hao, Adiabatic temperature changes in a one-component fluid near the liquidevapor critical point, Phys. A 197 (1993) 23e46. €hlich, P. Carle s, B. Zappoli, [13] Y. Garrabos, M. Bonetti, D. Beysens, F. Perrot, T. Fro Relaxation of a supercritical fluid after a heat pulse in the absence of gravity effects: theory and experiments, Phys. Rev. E 57 (1998) 5665e5681. s, F. Zhong, M. Weilert, M. Barmatz, Temperature and density relax[14] P. Carle ation close to the liquid-gas critical point: an analytical solution for cylindrical cells, Phys. Rev. E 71 (2005) 041201. [15] L.S. de B. Alves, An integral transform solution for unsteady compressible heat transfer in fluids near their thermodynamics critical point, Therm. Sci. 17 (2013) 673e686. [16] Y. Masuda, T. Aizawa, M. Kanakubo, N. Saito, Y. Ikushima, One dimensional heat transfer on the thermal diffusion and piston effect of supercritical water, Int. J. Heat Mass Transfer 45 (2002) 3673e3677. [17] R.I. Issa, A.D. Gosman, A.P. Watkins, The computation of compressible and incompressible recirculating flows by a non-iterative implicit scheme, J. Comput. Phys. 62 (1986) 66e82. [18] A. Nakano, M. Shiraishi, Numerical simulation for the Piston effect and thermal diffusion observed in supercritical nitrogen, Cryogenics 44 (2004) 867e873. [19] A. Nakano, M. Shiraishi, Piston effect in supercritical nitrogen around the pseudo-critical line, Int. Commun. Heat Mass Transfer 32 (2005) 1152e1164. [20] A. Nakano, Studies on piston and soret effects in a binary mixture supercritical fluid, Int. J. Heat Mass Transfer 50 (2007) 4678e4687. [21] C.A.J. Fletcher, Computational Technique for Fluid Dynamics, in: Scientific Computation, second ed., vol. 2, Springer-Verlag, Berlin, Germany, 1988. [22] P. Zhang, B. Shen, Thermoacoustic wave propagation and reflection near the liquidegas critical point, Phys. Rev. E 79 (2009) 060103. [23] B. Shen, P. Zhang, On the transition from thermoacoustic convection to diffusion in a near-critical fluid, Int. J. Heat Mass Transfer 53 (2010) 4832e4843. [24] B. Shen, P. Zhang, Thermoacoustic waves along the critical isochore, Phys. Rev. E 83 (2011). [25] T.J. Poinsot, S.K. Lele, Boundary conditions for direct simulations of compressible viscous flows, J. Comput. Phys. 101 (1992) 104e129. [26] N. Hasan, B. Farouk, Thermoacoustic transport in supercritical fluids at nearcritical and near-pseudo-critical states, J. Supercrit. Fluids 68 (2012) 13e24. [27] J.P. van Doormaal, G.D. Raithby, Enhancements of the simple method for predicting incompressible fluid flows, Numer. Heat Transfer 7 (1984) 147e163. [28] S. Paolucci, On the Filtering of Sound from the NaviereStokes Equations, Technical Report Technical Report SAND 82-8257, Sandia National Laboratories, 1982. [29] S. Klainerman, A. Majda, Compressible and incompressible fluids, Commun. Pure Appl. Math. 35 (1982) 629e651. [30] B. Shen, P. Zhang, RayleigheBenard convection in a supercritical fluid along its critical isochore in a shallow cavity, Int. J. Heat Mass Transfer 55 (2012) 7151e7165. [31] B. Zappoli, Influence of convection on the piston effect, Int. J. Thermophys. 19 (1998) 803e815. [32] G. Accary, I. Raspo, P. Bontoux, B. Zappoli, An adaptation of the low Mach number approximation for supercritical fluid buoyant flows, C. R. Mec. 333 (2005) 397e404. [33] G. Accary, I. Raspo, A 3d finite volume method for the prediction of a supercritical fluid buoyant flow in a differentially heated cavity, Comput. Fluids 35 (2006) 1316e1331. nard convection in a [34] G. Accary, P. Bontoux, B. Zappoli, Turbulent RayleigheBe near-critical fluid by three-dimensional direct numerical simulation, J. Fluid Mech. 619 (2009) 127e145. s, J. Ouazzani, Thermoacoustic and buoyancy [35] B. Zappoli, S. Amiroudine, P. Carle driven transport in a square side-heated cavity filled with a near-critical fluid, J. Fluid Mech. 316 (1996) 53e72. [36] B. Zappoli, A. Jounet, S. Amiroudine, A. Mojtabi, Thermoacoustic heating and cooling in near-critical fluids in the presence of a thermal plume, J. Fluid Mech. 388 (1999) 389e409. [37] S. Amiroudine, P. Bontoux, P. Larroude, B. Gilly, B. Zappoli, Direct numerical simulation of instabilities in a two-dimensional near-critical fluid layer heated from below, J. Fluid Mech. 442 (2001) 119e140.
Please cite this article in press as: P.C. Teixeira, L.S.B. Alves, Modeling supercritical heat transfer in compressible fluids, International Journal of Thermal Sciences (2014), http://dx.doi.org/10.1016/j.ijthermalsci.2014.08.011
12
P.C. Teixeira, L.S.B. Alves / International Journal of Thermal Sciences xxx (2014) 1e12
[38] S. Amiroudine, B. Zappoli, Piston-effect-induced thermal oscillations at the nard threshold in supercritical he, Phys. Rev. Lett. 90 (2003) 1e4. RayleigheBe [39] I. Raspo, B. Zappoli, P. Bontoux, Unsteady two-dimensional convection in a bottom heated supercritical fluid, C. R. Mec. 332 (2004) 353e360. [40] R.L. Panton, Incompressible Flow, second ed., John Wiley & Sons, United States, 1996. [41] E. Turkel, Review of preconditioning methods for fluid dynamics, Appl. Numer. Math. 12 (1993) 257e284. [42] E. Turkel, Preconditioning techniques in computational fluid dynamics, Annu. Rev. Fluid Mech. 31 (1999) 385e416. [43] S. Venkateswaran, C.L. Merkle, Analysis of Preconditioning Methods for the Euler and NaviereStokes Equations, in: Von Karman Institute Lecture Series, 1999. [44] Y.H. Choi, C.L. Merkle, The application of preconditioning in viscous flows, J. Comput. Phys. 105 (1993) 207e223. [45] J.M. Weiss, J.P. Maruszeski, W.A. Smith, Implicit solution of preconditioned NaviereStokes equations using algebraic multigrid, AIAA J. 37 (1999) 29e36. [46] I. Mary, P. Sagaut, M. Deville, An algorithm for low Mach number unsteady flows, Comput. Fluids 29 (2000) 119e147. s, B. Zappoli, Numerical solutions of 1D [47] S. Amiroudine, J. Quazzani, P. Carle unsteady near-critical fluid flows using finite volume methods, Eur. J. Mech. B Fluids 16 (1997) 665e680. [48] C.L. Merkle, Y.H. Choi, Computation of low speed compressible flows with time-marching methods, Int. J. Numer. Meth. Eng. 25 (1988) 293e311. [49] E. Turkel, A. Fiterman, B. van Leer, Preconditioning and the Limit to the Incompressible Flow Equations, Technical Report 93e42, ICASE, 1993. [50] J. Sesterhenn, B. Müller, H. Thomann, On the cancelation problem in calculating compressible low Mach number flows, J. Comput. Phys. 151 (1999) 597e615. [51] S.-H. Lee, Convergence characteristics of preconditioned Euler equations, J. Comput. Phys. 208 (2005) 266e288. [52] S.-H. Lee, Cancellation problem of preconditioning method at low Mach numbers, J. Comput. Phys. 225 (2007) 1199e1210.
[53] J.S. Shuen, K.H. Chen, Y.H. Choi, A coupled implicit method for chemical nonequilibrium flows at all speeds, J. Comput. Phys. 106 (1993) 306e318. [54] E. Turkel, Preconditioned methods for solving the incompressible and low speed compressible equations, J. Comput. Phys. 72 (1987) 277e298. [55] D.L. Darmofal, P.J. Schmid, The importance of eigenvectors for local preconditioners of the Euler equations, J. Comput. Phys. 127 (1996) 346e362. [56] L.S.B. Alves, Transverse Jet Shear-layer Instabilities: Linear Stability Analysis and Numerical Simulations (Ph.D. thesis), University of California at Los Angeles, Los Angeles, USA, 2006. [57] R. de S. Teixeira, L.S. de B. Alves, Modeling far field entrainment in compressible flows, Int. J. Comput. Fluid Dyn. 26 (2012) 67e78. [58] R.M. Cotta, Integral Transforms in Computational Heat and Fluid Flow, CRC Press, Boca Raton, FL, 1993. [59] R.M. Cotta (Ed.), The Integral Transform Method in Thermal and Fluids Science and Engineering, Begell House, Inc., New York, 1998. € Heat Conduction, second ed., Wiley Interscience, New York, 1993. [60] M.N. Ozisik, [61] L.S. de B. Alves, R.M. Cotta, Transient natural convection inside porous cavities: hybrid numericaleanalytical solution and mixed symbolic-numerical computation, Numer. Heat Transfer Part A Appl. 38 (2000) 89e110. [62] L.S. de B. Alves, R.M. Cotta, M.D. Mikhailov, Covalidation of hybrid integral transforms and method of lines in nonlinear convectionediffusion with mathematica, J. Braz. Soc. Mech. Sci. RBCM 23 (2001) 303e320. [63] L.S. de B. Alves, R.M. Cotta, J. Pontes, Stability analysis of natural convection in porous cavities through integral transforms, Int. J. Heat Mass Transfer 45 (2002) 1185e1195. nard [64] L.S. de B. Alves, A. Barletta, Convective instability of the DarcyeBe problem with through flow in a porous layer saturated by a power-law fluid, Int. J. Heat Mass Transfer 62 (2013) 495e506. [65] C. Moler, C.V. Loan, Nineteen dubious ways to compute the exponential of a matrix, twenty-five years later, SIAM Rev. 45 (2003) 1e46. [66] H. Schlichting, Boundary-layer Theory, seventh ed., MacGraw Hill, Inc., New York, 1986.
Please cite this article in press as: P.C. Teixeira, L.S.B. Alves, Modeling supercritical heat transfer in compressible fluids, International Journal of Thermal Sciences (2014), http://dx.doi.org/10.1016/j.ijthermalsci.2014.08.011
