A new formulation of variable turbulent Prandtl number for heat transfer to supercritical fluids

International Journal of Heat and Mass Transfer 92 (2016) 792–806

Contents lists available at ScienceDirect

International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt

A new formulation of variable turbulent Prandtl number for heat transfer to supercritical fluids Yoon Y. Bae Korea Atomic Energy Research Institute, 989-111 Daedeokdaero, Yuseong, Daejeon, 34057, Republic of Korea

a r t i c l e

i n f o

Article history: Received 6 February 2015 Received in revised form 15 September 2015 Accepted 15 September 2015 Available online 30 September 2015 Keywords: Turbulent Prandtl number Reynolds analogy Mixed convection Supercritical pressure Strong property variation

a b s t r a c t When a fluid at supercritical pressure approaches the pseudo-critical temperature it experiences a strong variation in physical properties putting applicability of various turbulent flow modelings in question. Earlier numerical calculations showed, without exception, unrealistic over-predictions, as soon as the fluid temperature approached the pseudo-critical temperature. The over-predictions might have been resulted either from an inapplicability of widely used turbulence models or from an unrealistic treatment of the turbulent Prandtl number (Prt) as a constant. Recent research, both numerical and experimental, indicates that Prt is very likely a function of fluid–thermal variables as well as physical properties, when the gradients of physical properties of a fluid are significant. This paper describes the procedure for a new formulation of Prt which varies with physical properties and fluid–thermal variables. The application of the variable Prt was surprisingly successful in reproducing the fluid temperature in supercritical fluids flowing in small-diameter vertical tubes ranging from 4.57 to 20 mm. Ó 2015 Elsevier Ltd. All rights reserved.

1. Introduction An accurate estimation of the heat transfer rate or temperature of the coolant channel is essential for the development of a supercritical pressure water cooled reactor (SCWR) [1]. Methods for predicting the heat transfer rate to or from supercritical fluids flowing in a very narrow passage are not satisfactory and have yet to be established. The two kinds of fluid, water and carbon dioxide (CO2), are mediums of interest and lot of works for the investigation are being conducted for applications in areas such as SCWR, Brayton cycle and compact printed circuit heat exchangers. A number of correlations for the prediction of the heat transfer rate in fluids at supercritical pressures have been proposed by various researchers, but most of them are applicable fluids in a forced convection regime, as shown in the review papers by Cheng and Schulenberg [2] and Pioro and Duffey [3]. The correlations available in literature predict the heat transfer rate with a reasonable accuracy in a forced convection regime; however, in a mixed convection regime, all of those correlations fail or partially succeed to produce accurate predictions, and the variation is so large that their application to the design needs to be very cautious. Since most of the earlier works have been summarized by Pioro and Duffey [3], several selected recent works are introduced here. Efforts, both experimental and analytical, have been made to for-

E-mail address: [email protected] http://dx.doi.org/10.1016/j.ijheatmasstransfer.2015.09.039 0017-9310/Ó 2015 Elsevier Ltd. All rights reserved.

mulate a reliable correlation for a mixed convection heat transfer by researchers such as Watts and Chou [4], Jackson and Hall [5], Jackson et al. [6], Bae and Kim [7], Bae et al. [8], Bae [9] and Jackson [10]. Zhu et al. [11] investigated the heat transfer characteristics of steam–water flowing upward in tubes at sub- and super-critical pressures in the range of 13–30 MPa. Yang et al. [12] performed an experiment on heat transfer to supercritical water flowing in vertical annular channels, and evaluated four correlations against the data. Li et al. [13] reported recent experimental results from the supercritical water heat transfer test facility SWAMUP at Shanghai Jiatong University. Zhao et al. [14] reported experimental results from the same research group with different conditions only to reveal that the existing heat transfer correlations did not correctly reproduce the heat transfer rate. In addition to the experimental efforts, a large number of numerical works have been performed to simulate the flow and thermal field in a fluid at supercritical pressures, and in doing so, the applicability of various turbulence models was examined. For both forced and mixed convection regimes, experimental and numerical investigations of the thermal and flow field at supercritical pressure was performed by Licht et al. [15]. They confirmed that for the simple case of deterioration, numerical simulations using the commercial CFD code Fluent offered a qualitative insight into changes in fluid temperature and turbulent velocities responsible for the axial evolution of the wall temperature. Cho et al. [16] examined three turbulence models, RNG k–e, SST k–x and one type

Y.Y. Bae / International Journal of Heat and Mass Transfer 92 (2016) 792–806

793

Nomenclature A+ Cl Ce1, Ce2 cp D G Gk h k p Pk Pr Prt Prt,o q r R Re T u, v u+ x y y+

effective viscous sublayer thickness constant in the turbulent viscosity constants in transport equation for e specific heat tube diameter mass flux production of turbulence due to buoyancy enthalpy turbulent kinetic energy pressure production of turbulence due to shear Prandtl number turbulent Prandtl number (variable) Prt before adjustment with additional functions heat flux radial coordinate tube radius Reynolds number temperature velocity in x and r direction non-dimensional u, u=us axial coordinate distance from the wall non-dimensional distance from wall, yus =m

of low-Reynolds number model, against the experimental data obtained for a tube and annulus with an equivalent hydraulic diameter of 4.4 mm, and reported that the performance of the three models was partially successful. He et al. [17] thoroughly investigated low-Reynolds number turbulence models and concluded that both the low Reynolds number k–e models and the V2F models were able to capture the general trends of the interesting wall temperature behavior observed with an upward flow in some experiments with a fluid at a pressure just above the critical value, while the detailed variation of the wall temperature predicted using each model was rather different from that in the experiments. They also found that the effect on the heat transfer was almost entirely due to the shear production effect caused by the distortion of the mean flow as a result of the strong influence of the buoyancy. Zhang et al. [18] successfully reproduced using a modified version of a low-Reynolds turbulence model the data from a DNS calculation and an experiment by employing an algebraic flux model in calculating the turbulence production based on the buoyancy. However, its application to other conditions is still to be proven, and the calculation domain was too small to generalize the results. Zhang et al. [19] compared the experimental heat transfer data in supercritical fluids in a circular tube with the calculation results obtained by employing six different turbulence models and found that the Reynolds stress model (RSM) gave the best agreement with the experimental data, especially with the deteriorated ones. The result of RSM was not much different from that of RNG k–e, and its applicability should be considered in parallel with the fact that it requires solving additional equations. Jaromin and Anglart [20] numerically performed a sensitivity analysis of the heated wall temperature and velocity distribution in the CFD simulation of the upward flow of supercritical water. They claimed that k–x turbulence model successfully simulated the initial temperature peak near the inlet and onset of deterioration, but without the recovery of heat transfer from deterioration at the pipe exit. In summary, the numerical works performed thus far only to prove that all current turbulence models are applicable to the cases with limiting conditions. In this context, a break-through is needed to simulate the fluid thermal behavior in a fluid with severe property

yþ TBL yþ r¼0:8R yþ upeak
a ~ a

y+ at the turbulent boundary layer edge (3 0 0) y+ at r = 0.8R y+ at @u=@y ¼ 0 Reynolds average quantity (a: dummy) Favre average quantity (a: dummy)

Greek symbols turbulent thermal diffusivity b volumetric expansion coefficient e dissipation rate of turbulent kinetic energy j von Karman constant l, lt molecular and turbulent viscosity m, mt molecular and turbulent kinematic viscosity q density rk, re model constants for turbulent diffusion of k, e rt standard turbulent Prandtl number (=0.9) st shear stress

at

Subscripts e effective (molecular + turbulent) o inlet w wall

variation, and focus has been given to Prt rather than trying to improve the turbulence modeling. Without exceptions all numerical works performed thus far struggled in simulating flow and thermal fields with strong property variation, that is, strong buoyancy or acceleration. Most of the turbulence modelings were developed based on an incompressible and constant-property flow. Afterwards, many efforts have been made to extend the models developed for constant property flows to variable property flows, however, they focused on high speed flow or compressible flows. Since the properties of fluids at supercritical pressure change with temperature than with pressure, a direct application of the theory developed for high speed or compressible flow to the cases of supercritical fluids should be very cautious. As expected, the application of the variable property (mainly density) version of turbulence modeling was not used in simulating a highly buoyant flow of fluids at supercritical flow, especially in reproducing fluid temperatures. Prt is a product of a pure perspective of similarity in appearance between the momentum and energy equations, and an assumption that their behavior would be the same or, at least, very similar to each other. Accordingly, it was treated as unity or an experimentallyobtained value slightly smaller than unity of 0.8–0.9. For most of fluid flows with barely varying fluid property, it has successfully worked and produced reasonable results, but never in the cases of strong property variation. In this regard, the author decided to revisit the Reynolds analogy, which connects the momentum and energy equations in the name of the Prt, and tried to find any possibility of extending it to be applied to the case of strong property variation. In all numerical works introduced above, including earlier works not mentioned here, Prt was treated as a constant or function of Prandtl number. A constant Prt, which was not successful in predicting the heat transfer in a deterioration regime, does not seem to properly represent the physics in the supercritical fluid in a vertical tube under a heat flux high enough to cause deterioration. The Prt is highly unlikely to be a constant when the fluid properties experience substantial variations. Quarmby and Quirk [21] measured the eddy diffusivity in air flowing through a plain tube

794

Y.Y. Bae / International Journal of Heat and Mass Transfer 92 (2016) 792–806

under high heat flux, and suggested an equation for Prt. Dai et al. [22] performed experiments of buoyant turbulent plume with carbon dioxide and hexafluoride, in which the density ratio between plume and atmosphere was 1.51 and 5.06, and found that Prt decreased from 0.8 at the centerline to a value smaller than 0.1 near the plume edge. It probably tells us that the density ratio plays an important role in the turbulent heat transfer. Kang et al. [23] experimentally showed that the value of Prt varied between 0.4 and 1.4. Kang and Iaccarino [24]and numerically by using LES revealed that Prt showed a wide variation ranging from 0.1 to 20. The effect of Prt has been studied by Mohseni and Bazargan [25], showing that a value of Prt smaller than unity results in a better agreement with the experimental data, for a deteriorated heat transfer. Liu et al. [26] proposed a new formulation of Prt, which depends on the temperature gradient, for the computation of the film cooling effectiveness. Their calculations clearly showed that the temperature-dependent Prt produced better results than using a constant Prt. Considering the above evidences of Prt being far from unity in the case of steep property gradient, there is sufficient reason to believe that Prt might be a function of flow and thermal variables, as well as the physical properties. In this paper, an attempt was made to develop a new formulation for Prt, which varies depending on the flow and thermal variables as well as the physical properties such as density, specific heat and temperature. After the completion of the development of Prt its applicability was examined against the experimental data for water, carbon dioxide and R22.

2.1. Governing equations and turbulence model In the present study, a vertically upward flowing fluid in a uniformly heated tube was considered. The flow was assumed to be steady and 2-D axi-symmetric. The vertical upward direction was aligned in the positive x direction and the radial coordinate was r. In the circumferential direction, the velocity components as well as the gradient of all flow and thermal properties were assumed to be zero. That is, there is no swirl. The governing equations employed in the present study were continuity equation, ensemble averaged Navier–Stokes equation, energy equation, two-equation turbulence model. The governing equations for a velocity field in a cylindrical coordinate with a two-equation k–e turbulence model are as follows [17]:


 
 ~ 1 @ @ @p 2 @ @u
u
u
gx þ
e ~ 2 Þ þ ðr q ~ v~ Þ ¼  þ q ðr q rl r @x @r @x r @x @x
  ~ @ v~ 1 @ @u
e rl þ þ r @r @r @x

   


 1 @ @ @ l
@ ~e 1 @ l
@ ~e
u
v~ ~eÞ ¼
þ t ~~eÞ þ ðrq ðr q r l l
þ t þ r @x @r @x re @x r @r re @r ~e ~e2
C e1 f e1 ðP k þ Gk Þ  q
C e2 f e2 ð6Þ þq ~ ~ k k The tilde and bar above the variables indicate the ensemble and 0 00 ~ F~ ¼ qF=q). For a Favre average, respectively (f ¼ F  F, f ¼ F  F, vertical upward flow in a tube, g x ¼ g and g r ¼ 0. The eddy viscosity was modeled by the model adopted earlier in Ref. [27]. The Reynolds stress tensor qu0i u0j (considered practically not 00 00
ug much different from q i uj ) was modeled through a Boussi-

nesque approximation.

qu0i u0j ¼ lt

  ~j @ u ~i @u 2 ~  dij qk þ 3 @xi @xj

ð7Þ

The production of turbulence energy by interaction with the mean flow and turbulence Pk was defined as

( "   2  2 #   ) ~ 2 ~ @ v~ 2 @u @ v~ v~ @u þ P k ¼ mt 2 þ þ þ @x @r @x @r r

ð8Þ

~ k @T u0i t0 ¼ ch u0i u0j ~e @xj

ð9Þ

In this work, the low-Reynolds number turbulence model proposed by Myong and Kasagi [29] (will be referred to hereinafter as MK) was used. In the low-Reynolds turbulence model, the eddy viscosity is defined as follows:

~2 k

l
t ¼ q
C l f l ~ e

ð10Þ

The coefficients appearing in the eddy viscosity and dissipation equation are defined as follows:

fl ¼

1þ

ð1Þ

3:45

!


Re1=2 t

 þ  y 1  exp  þ A

ð11Þ

f e1 ¼ 1

ð12Þ

( f e2 ¼ ð2Þ


   ~ 1 @ @ @p 1 @ @ v~ @ u
e ~ v~ Þ þ ðr q
u
v~ 2 Þ ¼  þ q
gr þ ðr q rl þ r @x @r @r r @x @x @r
 
e v~ 2 @ @ v~ l
e rl ð3Þ 2 2 þ r @r @r r


( "   #

1 @ l
t @ h~ ~ þ @ ðr q ~ ¼ þ1 @ r l
u
v~ hÞ ~ hÞ ðr q þ r @x @r r @x Pr Prt @x "   #) @ l
l
t @ h~ þ r þ @r Pr Prt @r

ð5Þ

The gravitational production Gk ¼ q0 u0i g i =q or Gk ¼ bg i u0i t0 was modeled using the generalized gradient diffusion hypothesis (GGDH) [28] as shown below.

2. Numerical method


 1 @ @ ~ Þ þ ðrq
u
v~ Þ ¼ 0 ðr q r @x @r

"
 #   @  1 @  @ l
t @ k~ ~~ ~ ~ ~


lþ r quk þ r qv h ¼ rk @x r @x @r @x "   # 1 @ l
@ k~
þ t
Gk  q
~e
Pk þ q r l þ þq r @r rk @r

1

"   #)
 þ 2 2 2 Ret y 1  exp  exp  9 6 5

ð13Þ

In the Myong–Kasagi model [29], A+, the modified effective sublayer thickness, in Eq. (11) is 70. The constants for the MK turbulence model are summarized in Table 1. Although there is no unique definition of the normalized wall distance y+, most of the computational works performed so far, regardless of compressible or incompressible cases, employed the

Table 1 Constants in the Myong–Kasagi turbulence model.

ð4Þ

Cl

Ce1

Ce2

rt

rk

re

Ch

0.09

1.40

1.80

0.9

1.4

1.3

0.3

795

Y.Y. Bae / International Journal of Heat and Mass Transfer 92 (2016) 792–806

definition yþ ¼ us y=mw , where us ¼ ðsw =qw Þ1=2 using wall properties. Instead, in this work, an altered definition based on the local properties such as yþ ¼ ðsw =qÞ1=2 y=m, which is called the semilocal scaling, was used. The only reason for the choice was that the computational results obtained in this work with the latter definition agreed with the experimental results better than those with the former definition. When only a narrow region near wall is heated to high temperature, the use of us ¼ ðsw =qw Þ1=2 will results

in an unrealistically high value of y+, making fl approach unity more rapidly as y+ increases. When the physical properties considerably vary the definition of us and y+ greatly influences the numerical results through the damping function and coefficients involving y+ as an independent variable. An exact assessment of the effect of the definition of us and y+ is not available at the moment and left for a further research. Huang et al. [30] also indicated that the semi-local scaling was the best in the interpretation of their DNS calculation results of compressible turbulent channel flows. Foysi et al. [31] reported that a use of the semi-local scaling improved the similarity in turbulence, although global similarity was not achieved.

2.2. Turbulent Prandtl number

Fig. 1. Variation of Prt with y+ at three axial locations. G = 166.62 kg/m2 s, q = 30.87 kW/m2, P = 8.0 MPa, d = 2 mm. After [36].

Fig. 2. Configuration of the tube used in this work.

Reynolds [32] reviewed more than thirty ways for predicting Prt and the Schmidt number (Sct). Kays [33] examined the thenavailable experimental data on Prt for the two-dimensional TBL. Prt has been treated as a constant around 0.9 or unity in most of earlier numerical works. However, there are other cases where Prt is far from unity. Prt becomes approximately 0.7 for an axisymmetric case of a heated jet, while planar jet data indicate a value of 0.5; in a thermally developing wall bounded TBL, Prt is around unity only in the core boundary layer, and it decreases to a value less than 0.5 [34]. In a TBL, Prt decreases from around 1.5 in the sub-layer to 0.7 at the outer edge of the boundary layer [35]. According to the experiment data provided by Dai et al. [22], the value of Prt reaches as small as 0.05. The above evidence strongly implies that Prt can hardly be unity or a constant very close to it at least in the case of heating or cooling of fluid experiencing substantial property variation. Prt is purely a product of the Reynolds analogy, which claims that the mechanism of turbulent heat transfer would be very similar to that of the turbulent momentum transfer. An assessment of Prt for the case (C) published in [36] was made and the result is shown in Fig. 1. Evidently, Prt has values much less than unity inside the turbulent boundary layer, especially in the buffer and viscous sublayer. The peaks are due to the strong gradient of velocity and temperature, which is inevitable irregularity originated from the application of gradient diffusion hypothesis in the turbulent momentum and energy diffusion. It should be noted that the values of Prt continuously decrease towards the wall

Table 2 Flow inlet conditions for the cases studies in this paper. Cases

Fluid

Pressure (MPa)

Inlet temp (K)

Mass flux (kg/m2 s)

Heat flux (kW/m2)

D (mm)

Reo

W1-a Li et al. [13] W1-b W1-c W1-d W2, Li et al. [13] W3, Shitsman [42] W4, Shitsman [43] W5 Vikhrev et al. [44] C1 [45] C2 [45] C3, unpublished KAERI data1 C4, unpublished KAERI data1 F1, Mori et al. [46] F2, Mori et al. [46] F3, Mori et al. [46]

Water

25.0

574

65,663

574 598 398.4 334.1 282 282 283

542.5 677 759 918 471 320 360 570 30 50 90

7.6

25.0 24.5 23.3 26.5 8.12 8.12 8.12

789.8 789.8 789.8 789.8 473 430 380 493 400 400 600

7.6 8 16 20.4 4.57 4.57 4.57

39,324 41,732 26,714 21,586 19,303 19,303 29,478

7.75

278.4

500

60

4.4

21,996

5.5 5.5 5.5

365.6 291.3 312.4

400 400 1000

10 25 90

4.4 4.4 4.4

25,694 9297 29,026

CO2

R22

1 These data have been produced with those published in [45] under the same research program executed by the author. The data were deposited in the OECD/NEA data bank after a thorough examination by the participants in the coordinated research program (CRP) supported by IAEA. The results of the CRP was documented in the IAEA TECDOC series, IAEA-TECDOC-1746 in 2014.

796

Y.Y. Bae / International Journal of Heat and Mass Transfer 92 (2016) 792–806

Fig. 3. Grid dependency for the case of carbon dioxide at supercritical pressure flowing upward in a circular tube [45].

contrary to the well-known experimental evidences, were, without exception Prt increased toward the wall. First of all, let us examine the steps taken to derive the Reynolds analogy for the turbulent heat transfer (the discussion of Kays et al. [37] was largely followed.). Imagine an element of fluid of mass dm that moves in the y direction at distance ‘ (which is a ‘‘mixing length”). Let us assume that the effective velocity of fluid in the y pffiffiffiffiffiffiffi (=R  r) direction is C v 02 . According to the momentum theorem, the shear force is equal to the rate of momentum transfer. Then, the effective shear stress is

pffiffiffiffiffiffiffi

pffiffiffiffiffiffiffi

F A

st ¼ ¼ C v 02 dðq
u~Þ ¼ C v 02 ðq
du~ þ u~dq
þ du~dq
Þ

ð14Þ

Similarly, the effective heat flux is

pffiffiffiffiffiffiffi

~ v 02 dðq
cp TÞ pffiffiffiffiffiffiffi  ~ p þ cp Td ~ q ~ q
cp dT~ þ q
þq
Tdc
dcp dT~ þ cp dq
dT~ þ Td
dcp ¼ C v 02 q

q_ 00t ¼ C

ð15Þ If the functional relations of incremental variation of the vari~ , dq
, dcp and dT~ with respect to each other across distance ables, du ‘ are known, the exact value of Prt can be calculated. In a conventional approach, the property variation was generally neglected under the assumption of constant properties or negligible variation of them. However, for a fluid at supercritical pressure, the extra ~ p þ cp Td ~ q ~ dq ~ dq
þ du
þ
and q
Tdc terms in Eqs. (14) and (15), u ~ ~ ~ q
dcp dT þ cp dq
dT þ Tdq
dcp , are no longer negligible, and can possibly be greater than the remaining terms. The second-order terms should not be neglected, if an extremely accurate prediction is required, since they amount to greater than 1% at some locations. When the mixing distance is sufficiently small, the increments ~ dq ~ , dT,
and dcp can safely be expressed as du

~¼‘ du

~ @u ; @y

dT~ ¼ ‘

@ T~ ; @y


¼‘ dq


@q ; @y

and dcp ¼ ‘

@cp @y

ð16Þ

In this work, the radial gradients are much larger than the axial ones, only radial derivatives are retained. When we insert Eq. (16) into Eqs. (14) and (15), we obtain after little rearrangement,









pffiffiffiffiffiffiffi ~ dq ~
jdq
j @u u st ‘ ¼ C v 02 1 þ þ @y q q
du~ q


ð17Þ

Similarly, the effective heat flux is

~ ! pffiffiffiffiffiffiffi
T dcp jdcp j jdqj q_ 00t T~ dq T~ dqdcp @ T~ ‘ ¼ C v 02 1 þ þ þ þ þ cp q
cp q
dT~ cp dT~ q cp q dT~ @y ð18Þ Please note that the incremental values in Eqs. (17) and (18), were taken as absolute values. This somewhat ad hoc arrangement can be justified if we recall that the turbulent mixing has no directional preference. Positive and negative gradient of relevant variables will equally play a mixing or diffusing role. The terms expressed as a ratio between incremental variation of the fluid property or flow variable can be easily replaced with derivatives with respect to the radial coordinate. However, the incremental ~ , dq
, dcp and dT~ need to be modeled. variations such as du Now, we begin to derive an expression for the mixing length, which holds from the wall to centerline of a tube or pipe. An empirical equation for mt =m for the entire region outside the viscous sublayer in a tube [37] is




 r 2 mt jyþ  r 1þ2 ¼ 1þ R R m 6



ð19Þ

At y ¼ 0, mt =m ¼ jyþ c =3. This value and corresponding mixing length will prevail in most of the TBL except viscous sublayer. In a viscous and log-law region, the mixing length is simply determined by the well-known function proposed by Van- Driest, ‘ ¼ jy½1  expðyþ =Aþ Þ, and em =m is defined by simple mixing length theory as





mt ‘2 @ u~ ¼ m m @y

ð20Þ

From Eqs. (19) and (20), ‘max , which is valid in most part of the TBL, is

‘max ¼

"

 mjyþc @ u~ 1 3

@y

#1=2 ð21Þ

Finally, the mixing length in a pipe is determined as follows.

797

Y.Y. Bae / International Journal of Heat and Mass Transfer 92 (2016) 792–806

‘ ¼ min ð‘Van Driest ; ‘max Þ

ð22Þ

By definition of momentum and thermal diffusivity,

~ q_ 00t st @u @ T~ ¼ mt ¼ at ; @y @y q q
cp

ð23Þ

Without the extra terms in the brackets of the right-hand side of pffiffiffiffiffiffiffi Eqs. (17) and (18), Prt ¼ mt =at would be unity since mt ¼ at ¼ C v 02 . However, with the extra terms Prt would not simply become unity, but ends up with a quite different value from unity. ~ , dq
, dcp and dT~ The incremental variations of the variables du

can be expressed as ‘ multiplied by corresponding differentials with respect to y. The insertion of those expressions into Eqs. (17) and (18), and from the definition of the turbulent Prandtl number, results in the basic variable turbulent Prandtl number Prt;o expressed as Prt;o

 . 

~ @u þ q
‘ @@yq 1 þ qu~
@@yq @y ¼  .   .  @c
~ ~ @c ~

@ T~ @ T~ þ cTp @yp þ c‘p @yp þ q
‘ @@yq þ q
‘cTp @@yq 1 þ qT
@@yq @y @y

@cp @y

 1 @ T~ @y ð24Þ

2

In a practical sense, the terms including the mixing length in Eq. (24) can be neglected. The value of Prt,o approaches unity or even a large number when the value of y+ is smaller than 10. In the viscosity-dominated region, turbulent heat flux will also be suppressed as much as the turbulent momentum transfer. Therefore, it was assumed that the Prt dies down toward the wall, following the same manner of the turbulent momentum transfer. The above argument allows us to introduce the following function, which is exactly the same as the damping function originally proposed by Van Driest, as the first factor to be added to the variable Prt,o.

 þ y f 1 ¼ 1  exp  þ A

ð25Þ

Considering that all gradients in flow variables and fluid properties will be reduced to zero around the tube center line, Prt is expected to approach a standard value, say, rt , in this direction as it does toward the wall. Therefore, the following function is introduced to take this assumption into account.


  B  yþ f 2 ¼ 0:5 1 þ tanh 10

ð26Þ

2

Water, Upward, P = 26 MPa, G = 789.9 kg/m s, q = 542.5 kW/m Grid: 100 (L) x 50 (R), GGDH, JPRT=9 , L 0=0.5 m, LQ=3.5 m, MK 750

2

2

Water, Upward, P = 26 MPa, G = 789.9 kg/m s, q = 677 kW/m Grid: 100 (L) x 50 (R), GGDH, JPRT=9 , L 0=0.5 m, LQ=2.5 m, MK

700

750

700

650

Tpc, 661.53 K

T,K

T,K

Tpc, 661.53 K

600

Experiment, T w

Experiment, T b

Calculation, T w

Calculation, Tb

650

Calculation, T axis

600

550 0.0

0.5

1.0

1.5

2.0

2.5

3.0

3.5

4.0

Experiment, T w

Experiment, Tb

Calculation, T w

Calculation, Tb

Calculation, T axis

x, m

550 0.0

(a)

0.5

1.0

1.5

2.0

2.5

3.0

2.5

3.0

x, m

(a) y+ at umax

10000

y+ at tube axis y+ at TBL edge (r=0.8R) Sublayer thickness

y+ at umax

10000

y+ at tube axis y+ at TBL edge (r=0.8R) Sublayer thickness

1000

1000

y+ y+ 100

100

10 0.0

0.5

1.0

1.5

2.0

2.5

3.0

3.5

4.0

x, m

(b) Fig. 4. Case W1-a, (a) temperature distribution, (b) distribution of y+ along the tube wall and axis. p = 26 MPa, G = 789.9 kg/m2 s, q = 542.5 kW/m2. Data from Li et al. [13].

10 0.0

0.5

1.0

1.5

2.0

x, m

(b) Fig. 5. Case W1-b, (a) temperature distribution, (b) distribution of y+ along the tube wall and axis. p = 26 MPa, G = 789.9 kg/m2 s, q = 677 kW/m2. Data from Li et al. [13].

798

Y.Y. Bae / International Journal of Heat and Mass Transfer 92 (2016) 792–806

The constant 10 is an arbitrary value given to make sure that the function f 2 varies smoothly around yþ ¼ B. The value of B should be carefully determined to make sure that beyond the point of yþ ¼ B, the flow is in the wake region. Prt derived here is valid only in the TBL and not in the wake region, where the mixing length theory is no longer valid. In a TBL, the log-law and wake region are separate from each other at around yþ ¼ B, and the momentum diffusivity asymptotically reaches a constant value as it enters the wake region, as can be found in the textbook dealing with TBL [34,37]. There is no reason not to believe that the mechanism of turbulent heat exchange differs from that of the momentum exchange. The experimental data provided by Quarmby and Quirk [21] and Dai et al. [22] show that it may be natural to assume that Prt asymptotically converges to 0.9 as y approaches R, say the tube centerline or axis. In the case of TBL on a flat plate, the boundary condition for e at the TBL outer edge is given as a free stream value. This cannot be the case for a TBL in a pipe, since e is calculated as a part of the solution rather than being defined or given. As can be seen in references [38,39], mt =m asymptotically approaches a constant value, where theoretically every property or variable should have nearly zero gradients. It is strongly indicated that Prt has a constant value

being close to a conventional value of 0.9; this is consistent with the fact that Prt asymptotically approaches a constant in TBL on a flat plate with zero pressure gradient [40]. Finally, with the incorporation of the functions introduced above, the variable Prt takes the following form.

Prt ¼ rt  f 1 f 2 ðrt  Prt;o Þ

ð27Þ

Eq. (27) was tested in the present numerical simulation of various experimental data with water, carbon dioxide and R22, especially those data plagued with heat transfer deterioration accompanying a sudden increase in wall temperature. 2.3. Numerical procedure The present numerical study was performed using an in-house code, a version modified from the one provided by Ferziger and Peric [41]. Basically, the SIMPLE algorithm with a single pressure correction step was applied. All variables were assigned to the collocated grids. Diffusive fluxes were discretized by busing the power law scheme; while, for convective fluxes, linear interpolation was allowed to be blended with an upwind approximation. The resulting matrix of the variables was iteratively solved by 2

2

Experiment, T w

800

Experiment, Tw

+

900

Calculation, Tb

Calculation, Tb

Calculation, T axis

Calculation, Tw (TBL_edge: r=0.8R)

Calculation, T axis

750

Experiment, T b

Calculation, Tw (TBL_edge: y =300)

Experiment, Tb

Calculation, T w

2

Water, Upward, P = 26 MPa, G = 789.9 kg/m s, q = 918.1 kW/m Grid: 100 (L) x 50 (R), GGDH, JPRT=9, L 0=0.5 m, LQ=2.5 m, MK

2

Water, Upward, P = 26 MPa, G = 789.9 kg/m s, q = 759 kW/m Grid: 100 (L) x 50 (R), GGDH, JPRT=9 , L 0=0.5 m, LQ=2.5 m, MK

T,K

T,K

800

700 Tpc, 661.53 K

700

650

Tpc, 661.53 K

600 600

550 0.0

0.5

1.0

1.5

2.0

2.5

0.0

3.0

1.0

1.5

x, m

(a)

(a)

y+ at umax

10000

0.5

x, m

1000

2.5

3.0

y+ at umax

10000

y+ at tube axis y+ at TBL edge (r=0.8R) Sublayer thickness

2.0

y+ at tube axis y+ at TBL edge (r=0.8R) Sublayer thickness

1000

y+

y 100

+

100

10

10 0.0

0.5

1.0

1.5

2.0

2.5

3.0

0.0

0.5

1.0

1.5

x, m

x, m

(b)

(b)

Fig. 6. Case W1-c, (a) temperature distribution, (b) distribution of y+ along the tube wall and axis. p = 26 MPa, G = 789.9 kg/m2 s, q = 759 kW/m2. Data from Li et al. [13].

2.0

2.5

3.0

Fig. 7. Case W1-d, (a) temperature distribution, (b) distribution of y+ along the tube wall and axis. p = 26 MPa, G = 789.9 kg/m2 s, q = 918 kW/m2. Data from Li et al. [13].

799

Y.Y. Bae / International Journal of Heat and Mass Transfer 92 (2016) 792–806

Stone’s SIP (strongly implicit procedure) method. The computational object was a vertical upward flow of water or carbon dioxide or R22 in uniformly heated circular tubes with several different diameters. The configuration of the tube is shown in Fig. 2. To obtain a thermally and dynamically fully developed turbulent velocity distribution at the inlet of the heated section, a velocity distribution proportional to a 1/7 power of reduced radius was assumed at the tube entrance; and a sufficiently long flow developing region with a length of 0.5 m was provided in front of the heated section. The length of the heated section was controlled to insure that the fluid temperature at the end is sufficiently higher than the critical temperature. The calculation domain was discretized into rectangular grids. The grid was refined into the wall in the radial direction. After try+ ing several values of yþ P (the value of y at the first node from the wall), yþ < 0:5 was found to be optimum with a result of reasonP ably accurate converged solutions. When the heat flux is very high, the situation gets more complex. As the wall temperature approaches the pseudo-critical temperature, an oscillation appears. This is considered purely numerical noise and has nothing to do with physics. A finer mesh as small as yþ P ¼ 0:002 reduced the amplitude of the oscillation, but did not remove it completely,

resulting in a convergence problem. The oscillation was suppressed by intentionally limiting the unrealistic sudden temperature jump in the immediate neighbor cell, such that, ½T  0:5ðT front þ T rear Þ < 5 K. The adjustment like this is purely arbitrary without any physical background, and should be determined case by case. However, this did not affect the overall calculations, while considerably alleviating the oscillation. The number of radial nodes greater than 50 did not conspicuously improve the predictions, while requiring substantially a longer time for convergence. The axial grid dependency was investigated by comparing 50, 100 and 150 nodes for the case of carbon dioxide with mass flux of 400 kW/m2 and heat flux of 50 kW/m2 (case C2 in Table 2), and the results are shown in Fig. 3. The results with 100 and 150 axial nodes did not show any noticeable improvement as compared with the result with 50 nodes. Since the purpose of this work was not to obtain accurate numerical results but to examine the applicability of newly developed Prt, a grid, 100(A)  50(R), which allows reasonable accuracy while requiring minimum computational works, was chosen. Grid dependency checked for the case of water was nearly similar to the case of carbon dioxide and R22. Generally speaking, the

2

2

1000

Experiment, Tw

700

Experiment, T b

Calculation, T w , modified TBL outer edge Calculation, T w, TBL outer edge: y

900

Experiment, T w

Experiment, T b

Calculation, T w

Calculation, Tb

Tpc, 663.4 K

Calculation, T axis

+ r=0.8R

600

Calculation, Taxis

T,K

Calculation, T b

2

Water, Upward, P = 24.5 MPa, G = 380 kg/m s, q = 360 kW/m Grid: 100 (L) x 50 (R), GGDH, JPRT=9, L 0=0.5 m, LQ=1.5 m, MK

2

Water, Upward, P = 25 MPa, G = 473 kg/m s, q = 471 kW/m Grid: 100 (L) x 50 (R), GGDH, JPRT= 9, L 0=0.5 m, LQ=3.0 m

T,K

800

500

700 Tpc, 661.53 K

600

400

0.0

500 0.0

0.5

1.0

1.5

2.0

2.5

3.0

0.5

1.0

1.5

2.0

x, m

3.5

x, m

(a)

(a) y+ at umax

y+ at umax

10000

y+ at tube axis y+ at TBL edge (r=0.8R) Sublayer thickness

y+ at tube axis y+ at TBL edge (r=0.8R) Sublayer thickness

1000

1000

y+

y+

100 100

10 0.0

10 0.0

0.5

1.0

1.5

2.0

2.5

3.0

3.5

4.0

x, m

(b) Fig. 8. Case W2, (a) temperature distribution, (b) distribution of y+ along the tube wall and axis. p = 25 MPa, G = 473 kg/m2 s, q = 471 kW/m2. Data from Li et al. [13].

0.5

1.0

1.5

2.0

x, m

(b) Fig. 9. Case W3, (a) temperature distribution, (b) distribution of y+ along the tube wall and axis. p = 24.5 MPa, G = 380 kg/m2 s, q = 360 kW/m2. Data from Shitsman [43].

800

Y.Y. Bae / International Journal of Heat and Mass Transfer 92 (2016) 792–806

calculations strongly depend on the radial number of nodes and weakly on the axial one. The flow conditions examined in this work are summarized in Table 2. All cases selected for the numerical experiment in this paper are for heat transfer deterioration, since it has been proved that the experimental data without heat transfer deterioration can be simulated quite accurately with the well-known turbulence models employing a constant Prt [16,19]. It should be noted that the Reynolds numbers at the inlet for all cases are considered in the range of fully turbulent flow. The following convergence criterion was set as the sum of the net flux of each flow variable (except pressure) normalized by the sum of local flux being less than 106. However, for some cases, the convergence criterion is relieved to a larger value of 103 since it could not be satisfied owing to the oscillation. Nevertheless, at the outlet, the velocity and enthalpy were continually renewed so that the mass and energy conservations were satisfied within 1% error. The fluid properties were calculated from a table constructed from the NIST standard reference database [47] with a suitable interval of enthalpy so that the temperature at the vicinity of the

pseudo-critical temperature (Tpc) is resolved as finer than 0.1 K. The fluid properties were linearly interpolated from the table.

3. Results and discussions 3.1. Water In this section, the results of a numerical simulation on the experiments performed using water as a working fluid are presented. The data were digitized from published papers. Fig. 4 shows the temperature distribution and the value of y+ along the tube wall and axis under the conditions of p = 26.5 MPa, G = 789.9 kg/m2 s, q = 542.5 kW/m2 s. In the title of plate (a) the acronym GGDH implies the generalized gradient diffusion hypothesis, and JPRT = 9 the currently used formula for Prt. In all graphs similar to Fig. 4, the solid lines from the top represent the wall, bulk and centerline temperature, respectively. Considering that the value of A+ is about 26 for a fully developed TBL on a flat plat without a pressure gradient and 70 in a lowReynolds number turbulence model, the value of A+ is given as

2

2

Experiment, T w

Experiment, Tw

Experiment, T b

Calculation, Tw

Calculation, Tb

+

Calculation, T w, A =70

900

Calculation, Taxis

Calculation, T b

Calculation, T axis

800

T,K

T,K

Experiment, Tb

Calculation, T w, modified viscous sublayer

1000

800

2

Water, Upward, P = 26.5 MPa, G = 493 kg/m s, q = 570 kW/m Grid: 100 (L) x 50 (R), GGDH, JPRT=9, L 0=0.5 m, LQ=6.0 m, MK

2

Water, Upward, P = 23.3 MPa, G = 430 kg/m s, q = 320 kW/m Grid: 100 (L) x 50 (R), GGDH, JPRT=9, L 0=0.5 m, LQ=2.0 m, MK

700

Tpc, 663.4 K

700 600

Tpc, 651.7 K

500

600

400 300

0.0

0.5

1.0

1.5

2.0

2.5

0

1

2

3

x, m

x, m

(a)

(a)

y+ at umax

y+ at umax

y+ at tube axis y+ at TBL edge (r=0.8R) Sublayer thickness

y+ at tube axis y+ at TBL edge (r=0.8R) Sublayer thickness

4

5

4

5

6

1000

1000

y+

y+

100

100

10

10 0.0

0.5

1.0

1.5

2.0

2.5

0

1

2

3

x, m

x, m

(b)

(b)

Fig. 10. Case W4, (a) temperature distribution, (b) distribution of y+ along the tube wall and axis. p = 23.3 MPa, G = 430 kg/m2s, q = 320 kW/m2. Data from Shitsman [43].

6

Fig. 11. Case W5, (a) temperature distribution, (b) distribution of y+ along the tube wall and axis. p = 26.5 MPa, G = 493 kg/m2 s, q = 570 kW/m2. Data from Vikhrev et al. [44].

801

Y.Y. Bae / International Journal of Heat and Mass Transfer 92 (2016) 792–806 30

x, m 0.3 (unheated) 0.6 u+ = 5.6 log(y+)+4.9 0.9 (u+ = ln(y+)/0.41+5.0) 1.2 1.5 1.8 2.1 2.4 u+ = y+ 2.7 3.0

25

20

4

u+

u, m/s

6

15

10

2

1.0 0.8 0.6 Prt

8

0.4 0.2

5 0

0.0 0

0.000

0.004

0.1

1

10

100

r, m

y+

(a)

(b)

0.1

1000

1

10

100

1000

y+

(c)

Fig. 12. Distribution of u, u+ and Prt along the tube. p = 26 MPa, G = 789.9 kg/m2 s, q = 918.1 kW/m2. Data from Li et al. [13].

26 and 70 in the regions of the unheated and heated parts of the test section for all calculations. The agreement between the calculation and measurement is extremely good. The agreement was relatively poor in the region near the point of the heating front 2

CO2, Upward, P = 8.12 MPa, G = 400 kg/m s, g = 30 kW/m 320

end. This is probably due to the fact that the thermal boundary begins developing from the point of heat input in the calculation, while it is already developed to some extent due to the thermal conduction in the direction of the tube inlet. As the thermal bound-

2 2

CO2, Upward, P = 8.12 MPa, G = 400 kg/m s, g = 50 kW/m

Grid: 100 (L) x 50 (R), GGDH, JPRT=9 , L 0=0.5 m, LQ=2.0 m, MK

2

Grid: 100 (L) x 50 (R), GGDH, JPRT=9 , L 0=0.5 m, LQ=2.0 m, MK 400

310

Tpc, 308.4 K 380

Experiment, Tw

Experiment, Tb

Calculation, T w

Calculation, T b

300

T,K

T,K

Calculation, T axis 360

290

Experiment, T w

Experiment, T b

Calculation, T w

Calculation, T b

320

Calculation, T axis

280 0.0

0.5

1.0

1.5

2.0

340

Tpc, 661.53 K

300

2.5

3.0

280

x, m

0.0

0.5

1.0

(a)

2.0

2.5

2.0

2.5

(a)

y+ at umax

y+ at umax

y+ at tube axis y+ at TBL edge (r=0.8R) Sublayer thickness

1000

1.5

x, m

y+ at tube axis y+ at TBL edge (r=0.8R) Sublayer thickness

1000

y+

y+

100

100

10

10

0.0

0.5

1.0

1.5

2.0

2.5

3.0

0.0

0.5

1.0

1.5

x, m

x, m

(b)

(b)

Fig. 13. Case C1, (a) temperature distribution, (b) distribution of y+ along the tube wall and axis. p = 8.12 MPa, G = 400 kg/m2 s, q = 30 kW/m2. Data from Bae et al. [45].

Fig. 14. Case C2, (a) temperature distribution, (b) distribution of y+ along the tube wall and axis. p = 8.12 MPa, G = 400 kg/m2 s, q = 50 kW/m2. Data from Bae et al. [45].

802

Y.Y. Bae / International Journal of Heat and Mass Transfer 92 (2016) 792–806

ary layer develops, the wall temperature increases. The outer edge of the TBL was defined as the line of y+ corresponding to r = 0.8R, hereinafter referred to as yþ r¼0:8R [38]. Beyond this line, the region becomes a wake region where the theory of the TBL no longer holds, and accordingly, Prt developed here cannot be applied. This definition of the TBL outer edge works very well when the flow remains in a state of normal or mild heat transfer deterioration. When heat transfer deterioration becomes severe, yþ r¼0:8R no longer represents the TBL outer edge, and the edge moves into the TBL. It can clearly be seen in Fig. 5(b) when compared with Fig. 4(b). An open square represents the location of @u=@y ¼ 0, hereinafter referred to as yþ upeak . As the heat flux increases from 542.5 to 677 kW/m2 s, as shown in Fig. 5, the heat transfer deterioration þ is seen apparently; and yþ upeak gets closer to yr¼0:8R . The result with a stronger heat flux of q = 759 kW/m2 s with the same mass flux is shown in Fig. 6. As expected, the heat transfer deterioration is more apparent. It should be note that yþ r¼0:8R still plays a role of the TBL outer edge. However, when the heat flux increases to q = 918 kW/m2 s (see Fig. 7), yþ upeak closely approaches yþ r¼0:8R . In this case the wake region penetrates well beyond yþ r¼0:8R ; and this line no longer represents the TBL outer edge. Once 2

CO2, Upward, P = 8.12 MPa, G = 600 kg/m s, g = 90 kW/m

a velocity peak penetrates the TBL, the downstream will experience the disturbances resulting from the penetration. For a remedy of this unacceptable situation, yþ TBL ¼ 300 was chosen as the line of the TBL outer edge once the TBL was penetrated þ by yþ upeak . The case with yTBL ¼ 300 apparently agrees with the experimental data better than the case with the line of r = 0.8R. þ The lines of yþ TBL ¼ 300 and yr¼0:8R will merge sufficiently far downstream where the flow field returns to the one corresponding to a typical constant-property turbulent flow in a circular tube. If the boundary between TBL and wake region is precisely defined, a more accurate prediction of the wall temperature can be achieved. However, at the moment, no experimental information is available for an accurate estimation of the boundary between the TBL and wake region. In cases W1 b–d, the reproduction of the wall temperature recovery after a heat transfer deterioration was in surprisingly good agreement with the experimental data. Thus far, we have discussed the cases of normal or mild heat transfer deterioration. We now examine the cases with relatively severe heat transfer deterioration. Fig. 8 shows the results for the case of p = 25 MPa, G = 473 kg/ m2 s and q = 471 kW/m2. Severe heat transfer deterioration is

2

2

CO2, Upward, P = 7.75 MPa, G = 500 kg/m s, g = 60 kW/m

Grid: 100 (L) x 50 (R), GGDH, JPRT=9 , L 0=0.5 m, LQ=2.0 m, MK

400

Experiment, T w

Experiment, T b

Calculation, T w

Calculation, T b

2

Grid: 100 (L) x 50 (R), GGDH, JPRT=9 , L 0=0.5 m, LQ=2.0 m, MK 400

Calculation, T axis

Experiment, Tw

Experiment, Tb

Calculation, T w

Calculation, T b

Calculation, T axis

T,K

T,K

350

350

300

Tpc, 661.53 K

Tpc, 306.4 K

300

0.0

0.5

1.0

1.5

2.0

0.0

2.5

0.5

1.0

1.5

x, m

x, m

(a)

(a)

y+ at umax

y+ at umax

y+ at tube axis y+ at TBL edge (r=0.8R) Sublayer thickness

y+ at tube axis y+ at TBL edge (r=0.8R) Sublayer thickness

1000

2.0

2.5

2.0

2.5

1000

y+

y+

100

100

10

10 0.0

0.5

1.0

1.5

2.0

2.5

0.0

0.5

1.0

1.5

x, m

x, m

(b)

(b)

Fig. 15. Case C3, (a) temperature distribution, (b) distribution of y+ along the tube wall and axis. p = 8.12 MPa, G = 600 kg/m2 s, q = 90 kW/m2. Data from KAERI’s unpublished data.

Fig. 16. Case C4, (a) temperature distribution, (b) distribution of y+ along the tube wall and axis. p = 7.75 MPa, G = 500 kg/m2 s, q = 60 kW/m2. Data from KAERI’s unpublished data.

803

Y.Y. Bae / International Journal of Heat and Mass Transfer 92 (2016) 792–806

apparent. As can be seen in the lower plate (b), the line yþ upeak penetrated deep into the TBL. This penetration will obviously alter the velocity profile in the TBL to result in the modification of corresponding thermal behavior. The wall temperature distribution þ based on yþ TBL ¼ yr¼0:8R did not agree well with the experimental þ data. This discrepancy might stem from the fact that yþ TBL ¼ yr¼0:8R corresponds to the location well beyond the TBL edge and deep in the wake region. The effect of the TBL location will be discussed further in Sections 3 and 4. In Fig. 9, case W3 of p = 24.5 MPa, G = 380 kg/m2 s and q = 360 kW/m2 is shown. The wall temperatures remained under Tpc. The calculated result qualitatively agrees with the experimental data. It should be noted that the calculation did not pick up the first peak of the experiment. It is presumed that the peak is weakly related to the strong property variation, since it appeared well below Tpc. In Fig. 10, case W4 of p = 23.3 MPa, G = 430 kg/m2 s, q = 320 kW/m2 is shown. The agreement is qualitatively good, but it was not able to capture the extremely high value. The case of the most severe heat transfer deterioration, case W5 of p = 26.5 MPa, G = 493 kg/m2 s, q = 570 kW/m2, is shown in Fig. 11. Again, as in case W3 shown in Fig. 9, the first peak was not captured, while the second peak was captured qualitatively well, although the location of the peak did not coincide. An extra simulation was conducted for case W5 with a domain including only the 2

first peak, reducing the axial heated length of the domain from 6 m to 1 m. It is commensurate to a grid refinement since we maintained the number of axial grids. This grid refinement did not capture the first either, and it may indicate that the peak appearing in this case did not originate from a severe property variation. It should be noted that the first peak appeared in cases W3 and W5 did not appear either in other water experiments or in the CO2 experiments. The typical distributions of u, u+ and Prt for case W1-d are shown in Fig. 12. As can be seen in plate (a), the velocity u at the exit still shows an overshoot. The low-density region with an overshoot initially appeared in the region very close to the viscous sublayer expands in the radial direction as the fluid moves downstream. The velocity distribution will eventually resume a typical one for a turbulent boundary layer in the region sufficiently far downstream. In plate (b), the distribution of u+ as a function of y+ is shown. As fluid enters the heated section, the u+ profile suddenly starts to deviate from the standard log-law profile. As the fluid moves further downstream to the exit, the u+ profile returns to the standard one, while in this case, the u+ profile at the exit still showed the influence of the velocity overshoot. The decrease in the value of u+ near the centerline corresponds to the velocity overshoot and increase of friction velocity due to density decrease. The distribution of Prt is shown in plate (c). The Prt showed a very

2

R22, Upward, P = 5.5 MPa, G = 400 kg/m s, q = 10 kW/m Grid: 100 (L) x 50 (R), GGDH, JPRT= 9, JTBL=1 , L 0=0.5 m, LQ=3.0 m

2

385

450

375

Experiment, T b

Calculation, T w

Calculation, Tb

Experiment, Tw

Experiment, Tb

Calculation, Tw

Calculation, T b

Calculation, Taxis

Calculation, T axis

400 Tpc 374.5K

Tpc 374.5 K

T,K

T,K

380

Experiment, T w

2

R22, Upward, P = 5.5 MPa, G = 400 kg/m s, q = 25 kW/m Grid: 100 (L) x 50 (R), GGDH, JPRT= 9, JTBL=1 , L 0=0.5 m, LQ=4.0 m

350 370

365

300 0.0

0.5

1.0

1.5

2.0

2.5

3.0

3.5

0

x, m

1

2

3

4

3

4

x, m

(a)

(a)

umax did not appear

y+ at umax

y+ at tube axis y+ at TBL edge (r=0.8R) Sublayer thickness

y+ at tube axis y+ at TBL edge (r=0.8R) Sublayer thickness

1000

y+

y+

1000

100

100

10

10

0.0

0.5

1.0

1.5

2.0

2.5

3.0

3.5

0

1

2

x, m

x, m

(b)

(b)

Fig. 17. Case F1, (a) temperature distribution, (b) distribution of y+ along the tube wall and axis. p = 5.5 MPa, G = 400 kg/m2 s, q = 10 kW/m2. Data from Mori et al. [46].

Fig. 18. Case F2, (a) temperature distribution, (b) distribution of y+ along the tube wall and axis. p = 5.5 MPa, G = 400 kg/m2 s, q = 25 kW/m2. Data from Mori et al. [46].

804

Y.Y. Bae / International Journal of Heat and Mass Transfer 92 (2016) 792–806

small value in the log-law region, and it gradually returned to a standard value of 0.9 toward the wall and centerline. The inner and outer boundaries of the low Prt region strongly influenced the numerical simulation, and a further study is required for a more robust simulation, although currently used boundaries, A+ = 70 and yþ r¼0:8R , served very well for the purpose of this paper.

3.2. Carbon dioxide In this section, the results of a numerical simulation for the experiments performed using carbon dioxide as a working fluid are presented. The calculation results for carbon dioxide are little different from those for water. A case of nearly normal heat transfer with the conditions of p = 8.12 MPa, G = 400 kg/m2 s, q = 30 kW/m2 is shown in Fig. 13. The calculation closely followed the experimental data with a slight overestimation in the region where Tbulk is close Tpc. As can be seen in the lower plate (b), a velocity peak appeared. However, the values of yþ upeak were too larger than that for r = 0.8R to exert any influence on the thermal behavior. Fig. 14 shows a case with a more severe heat transfer deterioration. As the heat flux increased to 50 kW/m2, with the same mass flux, a severe wall 2

temperature rise appeared. In the region of the front end heating, þ yþ upeak reached well below 70. Such a small value of yupeak will influence the fluid–thermal behavior in the TBL anyway by modifying it, þ as discussed earlier for case W2. yþ TBL ¼ yr¼0:8R worked very well, and no modification was needed, different from case W2. To see the effect of different mass and heat flux, case C3 of G = 600 kg/ m2 s, q = 90 kW/m2 was simulated, and the results are shown in Fig. 15. The calculated wall temperatures agree very well with the experimental data. As can be seen in cases C2 and C3, the accuracy of reproducing the recoveries from the heat transfer deterioration was surprising. Case C4, whose results are shown in Fig. 16, was chosen to investigate the effect of pressure. The pressure in case C4 is 7.75 MPa, which is closer to the Tpc than that of cases C1–C3. The results agreed well with the experimental data, as in cases C1–C3.

3.3. Freon (HCFC22) In this section, the results of a numerical simulation for the experiments performed using R22, a type of Freon, as the working fluids are presented. The data were digitized from the paper pre2

440

Experiment, T w

Experiment, T b

Calculation, T w

Calculation, Tb

+

900

Calculation, Tw, Original Calculation, Tb

850

Calculation, T axis

800

T,K

T,K

Experiment, Tb

Calculation, Tw, TBL and A modified

400 380

Experiment, T w

950

Calculation, T axis

420

2

Water, Upward, P = 25 MPa, G = 473 kg/m s, q = 471 kW/m Grid: 100 (L) x 50 (R), GGDH, JPRT= 9, L 0=0.5 m, LQ=3.0 m

2

R22, Upward, P = 5.5 MPa, G = 1000 kg/m s, q = 90 kW/m Grid: 100 (L) x 50 (R), GGDH, JPRT= 9, JTBL=1 , L 0=0.5 m, LQ=2.5 m

Tpc, 374.5 K

750 700

360

650

340

Tpc, 661.53 K

600 320

550 300

0.0 0.0

0.5

1.0

1.5

2.0

2.5

0.5

1.0

1.5

3.0

2.0

2.5

3.0

3.5

4.0

2.5

3.0

3.5

4.0

x, m

x, m

(a)

(a) y+ at umax

y+ at umax

y+ at tube axis y+ at adjusted TBL edge y+ at 0.8R Sublayer thickness

y+ at tube axis y+ at TBL edge (r=0.8R) Sublayer thickness 1000

y+

1000

y+

100

100

10 0.0

0.5

1.0

1.5

2.0

2.5

3.0

x, m

(b) Fig. 19. Case F3, (a) temperature distribution, (b) distribution of y+ along the tube wall and axis. p = 5.5 MPa, G = 1000 kg/m2 s, q = 90 kW/m2. Data from Mori et al. [46].

10 0.0

0.5

1.0

1.5

2.0

x, m

(b) Fig. 20. Case W2. Effect of the modification of sublayer thickness and TBL outer edge. p = 25 MPa, G = 473 kg/m2 s, q = 471 kW/m2. Data from Li et al. [13].

Y.Y. Bae / International Journal of Heat and Mass Transfer 92 (2016) 792–806

sented by Mori et al. [46]. In (b) Fig. 17, the results for case F1 of p = 5.5 MPa, G = 400 kg/m2 s and q = 10 kW/m2 are shown. In this case, the overall performance of the numerical simulation is fairly good as expected for the normal heat transfer. It should be noted that velocity overshoot did not appear in this case. As the heat flux increased to 25 kW/m2, while other conditions remained the same, a sudden wall temperature rise appeared when the bulk temperature approached Tpc, as shown in (b) Fig. 18. The overall performance was reasonably good; however, except for the region close to the temperature peak, the simulation significantly underestimated the experimental data. The results for case F3 of p = 5.5 MPa, G = 1000 kg/m2 s and q = 90 kW/m2 are shown in (b) Fig. 19. In this case, the simulation generally followed the experiment, but an underestimation was apparent. In all cases for fluid R22, the applicability of the newly developed Prt was also successfully confirmed.

3.4. Effect of sublayer thickness and TBL outer edge In this section, the effects of the modification of the viscous sublayer thickness, A+, and the TBL outer edge, yþ TBL , are discussed. In 2

2

Water, Upward, P = 26.5 MPa, G = 493 kg/m s, q = 570 kW/m Grid: 100 (L) x 50 (R), GGDH, JPRT=9, L 0=0.5 m, LQ=6.0 m, MK Experiment, Tw

Experiment, T b

Calculation, Tw, modified viscous sublayer

1000

+

Calculation, Tw, A =70

900

Calculation, Tb

Calculation, Taxis

T,K

800 Tpc, 663.4 K

700 600 500 400 300 0

1

2

3

4

5

y+ at umax y+ at tube axis y+ at TBL edge (r=0.8R) Sublayer thickness modified Sublayer thickness 1000

y+ 100

10 2

3

it was suspected that both A+ and yþ TBL were significantly changed from the normal state. yþ TBL was newly defined such that it ran with þ the line of yþ upeak when yupeak < 70 and asymptotically and slowly approached yþ r¼0:8R as shown in the lower plate. This modification is totally arbitrary and has to be verified later experimentally. The value of A+ was also modified to reflect the severe penetration of line yþ upeak into the viscous sublayer such that it had the same þ value of yþ upeak whenever yupeak was smaller than 70; otherwise, it

remained as 70. The wall temperature distributions indicated by solid (modified) and dashed (original) lines correspond to the two different yþ TBL indicated by the dark yellow (modified) and blue (original: yþ ¼ yþ TBL r¼0:8R ) triangles, respectively. The result with the þ modified yTBL and A+ agreed much better with the experimental data. However, this adjustment using an ad hoc introduction of new definition of yþ TBL might be not enough to explain the complicated behavior in the TBL with strong buoyancy. The factors controlling the fluid thermal behavior in the TBL include not only the definition of the TBL outer edge but also many others including turbulence properties. The results of the comparison here, how+ ever, convincingly indicate that the values of yþ TBL and A have a strong influence on the fluid–thermal behavior and eventually on the wall temperature. It should be noted that the modification of + either yþ TBL or A did not result in the expected results. + The influence of the value of yþ TBL and A on another case, W5, was also examined. In (b) Fig. 21, the effect of modification of only A+ is shown. Clearly, the modification resulted in a better agreement between the calculation and experiment. The above two cases were presented solely for the purpose of + showing the influences of the modification of yþ TBL and A , which were indeed very strong. However, its physical background is still to be proven through further numerical and experimental efforts. It is hoped that a detailed study of the TBL severely modified by the irregular velocity distribution such as M-shape velocity distribution will reveal any causes, if any, for the sudden temperature jump and drop in fluids flowing under conditions of supercritical pressure and temperature. 4. Conclusions

(a)

1

Fig. 20, the effects of the modification of both A+ and yþ TBL are shown for case W2. Since the line of yþ upeak penetrated the viscous sublayer,

6

x, m

0

805

4

5

6

x, m

(b) Fig. 21. Effect of the modification of sublayer thickness. p = 26.5 MPa, G = 493 kg/ m2 s, q = 570 kW/m2. Data from Vikhrev et al. [44].

The value of Prt has long been considered to be unity or close to it. Most of the numerical works to predict the wall and fluid temperature flowing in tubes and channels at supercritical pressures have failed or were partially successful in a limited region of mild property variation. In particular, attempts to simulate the recovery from the deterioration, occurring in a supercritical heat transfer, were failed without exception. The author suspected that one of the major reasons for the failure or inadequacy was an inappropriateness of Prt. Several experimental data and numerical results also indicated that Prt can be very smaller or larger than unity in a region of severe property variation. In this paper, an attempt has been made to formulate a new variable Prt with consideration of the effect of severe physical property variation as well as the flow properties. The formulation was basically based on the well-known mixing length theory. The developed variable Prt was applied to the numerical estimation of the available experimental data for water, carbon dioxide and R22. The calculation results agreed surprisingly well with the experimental data. The recovery of heat transfer from deterioration that has never been reproduced was also successfully reproduced in the numerical calculation employing the newly developed Prt. The application of newly developed Prt would not be limited to the case of supercritical pressure fluids; and rather it can be

806

Y.Y. Bae / International Journal of Heat and Mass Transfer 92 (2016) 792–806

applied for any numerical simulations with fluids experiencing a severe property variation. The claims put forward here are based purely on the speculation and numerical experiments without any experimental evidences, which are strongly anticipated in the near future to challenge the claims given here. Conflict of interest None declared. Acknowledgments The author would like to acknowledge the financial support by the Nuclear Research & Development Program of the National Research Foundation of Korea (NRF) grant funded by the Korean government (MSIP). (Grant code: NRF-2012M2A8A2025682). Dr. J.H. Bae’s release of his valuable DNS data should also be appreciated. References [1] A technology roadmap for generation IV nuclear energy systems, GIF-002-00, USDOE Nuclear Energy Advisory Committee and the Generation IV International Forum, December 2002. [2] X. Cheng, T. Schulenberg, Heat Transfer at Supercritical Pressures—Literature Review and Application to an HPLWR, Wissenschaftliche Berichte (Tech. Report) FZKA 6609, Forschungszentrum Karlsruhe, Mai, 2001. [3] I.L. Pioro, R.B. Duffey, Experimental heat transfer in supercritical water flowing inside channels (survey), Nucl. Eng. Des. 235 (22) (2005) 2407–2430. [4] M.J. Watts, C.T. Chou, Mixed convection heat transfer to supercritical pressure water, in: Proceedings of the 7th International Heat Transfer Conference, München, 1982, pp. 495–500. [5] J.D. Jackson, W.B. Hall, Influences of Buoyancy on heat transfer to fluids flowing in vertical tubes under turbulent condition, in: S. Kakaç, D.B. Spalding (Eds.), Turbulent Forced Convection in Channels and Bundles, Hemisphere Publishing, 1979, pp. 613–640. [6] J.D. Jackson, M.A. Cotton, B. Axcell, Studied of mixed convection in vertical tubes, Int. J. Heat Fluid Flow 10 (1) (1989) 2–15. [7] Y.Y. Bae, H.Y. Kim, Convective heat transfer to CO2 at a supercritical pressure flowing vertically upward in tubes and an annular channel, Exp. Therm. Fluid Sci. 33 (2) (2009) 329–339. [8] Y.Y. Bae, H.Y. Kim, D.J. Kang, Forced and mixed convection heat transfer to supercritical CO2 vertically flowing in a uniformly-heated circular tube, Exp. Therm. Fluid Sci. 34 (2010) 1295–1308. [9] Y.Y. Bae, Mixed convection heat transfer to carbon dioxide flowing upward and downward in a vertical tube and an annular channel, in: Proceedings of the 1st Meeting of International Specialists on Supercritical Pressure Heat Transfer and Fluid Dynamics, University of Pisa, Pisa, Italy, July 5–8, 2010. [10] J.D. Jackson, An extended model of variable property developing mixed convection heat transfer in vertical tubes, in: Proc. of 1st Meeting of International Specialists on Supercritical Pressure Heat Transfer and Fluid Dynamics, University of Pisa, Pisa, Italy, July 5–8, 2010. [11] X.J. Zhu, Q.C. Bi, T.K. Chen, An investigation on hat transfer characteristics of steam–water at different pressure in vertical upward tube, in: 3rd Int. Symposium on SCWR, 2007, Shanghai, China, March 12–15, 2007, Paper No. SCR2007-P023. [12] Z. Yang, Q. Bi, H. Wang, G. Wu, R. Hu, Experiment of heat transfer to supercritical water flowing in vertical annular channel, J. Heat Transfer 135 (2013). 042504-1. [13] H. Li, J. Yang, H. Gu, D. Lu, J. Zhang, F. Wang, Y. Zhang, Heat transfer research on supercritical water flow upward in tube, ICAPP’12, 2012, Paper 12435. [14] M. Zhao, H. Gu, X. Cheng, Experimental study on heat transfer to supercritical water flowing through tubes, ICAPP’12, 2012, Paper 12368. [15] J. Licht, M. Anderson, M. Corradini, Heat transfer and fluid flow characteristics in supercritical pressure water, J. Heat Transfer 131 (July 2009). 072502-1. [16] B.H. Cho, Y.I. Kim, Y.Y. Bae, Prediction of a heat transfer to CO2 flowing in an upward path at a supercritical pressure, Nucl. Eng. Technol. 41 (7) (2009) 907– 920. [17] S. He, W.S. Kim, J.D. Jackson, A computational study of convective heat transfer to carbon dioxide at a pressure just above the critical value, Appl. Therm. Eng. 28 (2008) 1662–1675.

[18] H, Zhang, Z.R. Xie, Y.H. Yang, Numerical study on supercritical fluids flow and heat transfer under buoyancy, in: The 8th International Topical Meeting on Nuclear Thermal-Hydraulics, Operation and Safety (NUTHOS-8), Shanghai, China, October 10–14, 2010, Paper No.: N8P0187. [19] Y. Zhang, C. Zhang, J. Jiang, Numerical simulation of heat transfer of supercritical fluids in circular tubes using different turbulence models, J. Nucl. Sci. Technol. 48 (3) (2011) 366–373. [20] M. Jaromin, H. Anglart, Sensitivity analysis of heated wall temperature and velocity distribution in CFD simulation of the upward flow of supercritical water, in: NURETH-14, Toronto, Ontario, Canada, Sept. 25–30, 2011, Paper No. 231. [21] A. Quarmby, R. Quirk, Measurements of the radial and tangential eddy diffusivities of heat and mass in turbulent flow in a plain tube, Int. J. Heat Mass Transfer 15 (1972) 2309–2327. [22] Z. Dai, L-K. Tseng, G.M. Faeth, Velocity/mixture fraction statistics of round selfpreserving, buoyant turbulent flumes, HTD-Vol. 304, National Heat Transfer Conference-Volume 2, ASME, 1995. [23] S. Kang, B. Patil, J.A. Zarate, R.P. Roy, Isothermal and heated turbulent upflow in a vertical annular channel – part I. Experimental measurements, Int. J. Heat Mass Transfer 44 (2001) 1171–1184. [24] S. Kang, G. Iaccarino, Computation of turbulent Prandtl number for mixed convection around a heated cylinder, in: Annual Research Briefs 2010, Center for Turbulence research, Stanford University, 2010, pp. 295–304. [25] M. Mohseni, M. Bazargan, Effect of turbulent prandtl number on convective heat transfer to turbulent flow of a supercritical fluid in a vertical round tube, J. Heat Transfer 133 (July) (2011). 071701-1. [26] C. Liu, H. Zhu, J. Bai, New development of the turbulent Prandtl number models for the computation of film cooling effectiveness, Int. J. Heat Mass Transfer 54 (2011) 874–886. [27] Y.Y. Bae, S.D. Hong, Y.W. Kim, Numerical simulation of flow and thermal field in supercritical pressure carbon dioxide flowing upward in a narrow tube, in: NURETH-14, Toronto, Ontario, Canada, Sept. 25–29, 2011, Log Number: 347. [28] B.E. Launder, RANS modelling of turbulent flow affected by buoyancy of stratification, in: G.F. Hewitt, J.C. Vassilicos (Eds.), Prediction of Turbulent Flows, Cambridge University Press, 2005. [29] H.K. Myong, N. Kasagi, A new approach to the improvement of k–e turbulence model for wall bounded shear flows, JSME Int. J. 33 (1990) 63–72. [30] P.G. Huang, G.N. Coleman, P. Bradshaw, Compressible turbulent channel flows: DNS results and modelling, J. Fluid Mech. 305 (1995) 185–218. [31] H. Foysi, S. Sarkar, R. Friedrich, Compressibility effects and turbulence scalings in supersonic channel flow, J. Fluid Mech. 509 (2004) 207–216. [32] A.J. Reynolds, The prediction of turbulent Prandtl and Schmidt numbers, Int. J. Heat Mass Transfer 18 (1975) 1055–1069. [33] W.M. Kays, Turbulent Prandtl number–where are we, J. Heat Transfer 116 (May) (1994) 284–295. [34] J.A. Schetz, Boundary Layer Analysis, Prentice Hall, Englewood Cliffs, New Jersey, 1993. p. 273. [35] F.M. White, Viscous Fluid Flow, second ed., McGraw Hill Int. Ed., 1991. p. 482. [36] J.H. Bae, J.H. Yoo, H. Choi, Direct numerical simulation of turbulent supercritical flows with heat transfer, Phys. Fluids 17 (2005) 105104. [37] W. Kays, M. Crawford, B. Weigand, Convective Heat and Mass Transfer, fourth ed., McGraw Hill International Edition, 2005. p. 231. [38] J.A. Schetz, Boundary Layer Analysis, Prentice Hall, Englewood Cliffs, New Jersey, 1993. p. 340. [39] W. Kays, M. Crawford, B. Weigand, Convective Heat and Mass Transfer, fourth ed., McGraw Hill International Edition, 2005. p. 284. [40] W. Kays, M. Crawford, B. Weigand, Convective Heat and Mass Transfer, fourth ed., McGraw Hill International Edition, 2005. pp. 233–239. [41] J.H. Ferziger, M. Peric, Computational Methods for Fluid Dynamics, SpringerVerlag, Berlin Heidelberg, 1999. [42] M.E. Shitsman, Natural convection effect on heat transfer to a turbulent water flow in intensively heated tubes at supercritical pressures, in: Proc. Instn. Mech. Engrs. 1967–1968, Paper 6. [43] M.E. Shitsman, Impairment of the heat transmission at supercritical pressures, High Temp. 1 (2) (1963) 237–244. [44] Yu.V. Vikhrev, Yu.D. Barulin, A.S. Kon’kov, A study of heat transfer in vertical tubes at supercritical pressure, Teploenergetika 14 (19) (1967) 80–82. [45] Y.Y. Bae, Mixed convection heat transfer to carbon dioxide flowing upward and downward in a vertical tube and an annular channel, Nucl. Eng. Des. 241 (2011) 3164–3177. [46] H. Mori, M. Ohno, K. Ohishi, Y. Hamamoto, Research and development of a super fast reactor (7) heat transfer to a supercritical pressure fluid flowing in a sub-bundle channel, in: 16th Pacific Basin Nuclear Conference (16PBNC), Aomori, Japan, Oct. 13–18, 2008, Paper ID P16P1297. [47] E.W. Lemmon, M.L. Huber, M.D. McLinden, Reference Fluid Thermodynamics and Transport Properties. NIST Standard Reference Database 23, Version 9.0, 2010.