Numerical simulation of double-pipe condensers and evaporators

International Journal of Refrigeration 27 (2004) 656–670 www.elsevier.com/locate/ijrefrig

Numerical simulation of double-pipe condensers and evaporators O. Garcı´a-Valladaresa,*, C.D. Pe´rez-Segarrab,1, J. Rigolab,1 a

Centro de Investigacio´n en Energı´a (CIE), Universidad Nacional Auto´noma de Me´xico (UNAM), Privada Xochicalco S/N, Temixco, 62580, Morelos, Mexico b Centre Tecnolo`gic de Transfere`ncia de Calor (CTTC), Lab. de Termote`cnia i Energe`tica. Universitat Polite`cnica de Catalunya (UPC), ETSEIT, Colom 11, 08222 Terrassa, Spain Received 10 August 2003; received in revised form 12 December 2003; accepted 27 January 2004

Abstract A detailed one-dimensional steady and transient numerical simulation of the thermal and fluid-dynamic behaviour of doublepipe heat exchangers (evaporators and condensers) has been carried out. The governing equations (continuity, momentum and energy) inside the internal tube and the annulus, together with the energy equation in the internal tube wall, external tube wall and insulation, are solved iteratively in a segregated manner. The discretized governing equations in the zones with fluid flow are efficiently coupled using an implicit step by step method. This formulation requires the use of empirical correlations for the evaluation of convective heat transfer, shear stress and void fraction. An implicit central difference numerical scheme and a line-by-line solver was used in the internal and external tube walls and insulation. A special treatment has been implemented in order to consider transitions (single-phase/two-phase, dry-out,…). All the flow variables (enthalpies, temperatures, pressures, mass fractions, velocities, heat fluxes,…) together with the thermophysical properties are evaluated at each point of the grid in which the domain is discretized. Different numerical aspects and comparisons with analytical and experimental results are presented in order to verify and validate the model. q 2004 Elsevier Ltd and IIR. All rights reserved. Keywords: Modelling; Heat transfer; Steady state; Transient behaviour; Evaporator; Annulus

Simulation nume´rique des condenseurs a` tube annulaire et des e´vaporateurs Mots-cle´s: Mode´lisation; Transfert de chaleur; Re´gime permanent; Re´gime transitoire; Condenseur; Evaporateur; Tube annulaire

1. Introduction Concern over the depletion of the stratospheric ozone layer by chemical compounds began some time ago At the same time, as the manufacturers faced the difficult task of * Corresponding author. Tel.: þ52-55-56229746; fax: þ 52-5556229791. E-mail addresses: [email protected] (O. Garcı´a-Valladares), [email protected] (C.D. Pe´rez-Segarra, J. Rigola). 1 Tel.: þ34-93-7398192; fax: þ34-93-7398101.

changing working fluids, concerns over energy efficiency were also gaining attention. All these aspects have had, and continue to have, an important impact on the HVAC industry. Various issues that have been addressed for years with ‘rules of thumb’ are being examined more rigorously. For the reasons mentioned above, and in order to optimise the efficiency of the heat exchangers, and consequently the energy consumption with non-contaminant refrigerants, more accurate and general methods for predicting their thermal and fluid-dynamic behaviour are required. The inherent complexity of the heat exchanger

0140-7007/$ - see front matter q 2004 Elsevier Ltd and IIR. All rights reserved. doi:10.1016/j.ijrefrig.2004.01.006

O. Garcı´a-Valladares et al. / International Journal of Refrigeration 27 (2004) 656–670

Nomenclature 2

A cross section area (m ) cp specific heat at constant pressure (J kg21 K21) CV control volume D diameter (m) e specific energy ðh þ n2 =2 þ gz sin uÞ (J kg21) f friction factor g acceleration due to gravity (ms22) G mass velocity (kg m22 s21) h enthalpy (J kg21) L length (m) m mass (kg) m _ mass flux (kg s21) n number of control volumes Nu Nusselt number ðaD=lÞ p pressure (Pa) P perimeter (m) q_ heat flux per unit area (W m22) Re Reynolds number ðGD=mÞ t time (s) T temperature (K) v velocity (ms21) xg vapour mass fraction z axial coordinate Greek letters u angle (rad); see Fig. 1 r density (kg m23) F two-phase frictional multiplier a heat transfer coefficient (W m22 K21) d rate of convergence

design in aspects such as geometries and fluid flow patterns, means that the possibilities of analytical solutions are very limited without assuming stringent simplifications (e.g. analytical approaches such as F-factor, e-NTU, etc.) [1,2]. Numerical methods allow the governing equations to be solved with fewer restrictions. Even though numerical treatment considering the multi-dimensionality of the twophase flow offers very restricted possibilities for solution, an important variety of technical situations can be treated assuming steady or transient one-dimensional flow. For

Fig. 1. Flow inside a control volume.

657

j absolute roughness (m) t shear stress (Pa) l thermal conductivity (W m21 K21) m dynamic viscosity (Pa s) eg void fraction Dt temporal discretization step (s) Dr radial discretization step (m) Dz axial discretization step (m) Subscript a annulus amb ambient bc begin condensation E, e east ec end condensation exp experimental g gas or vapour h hydraulic l liquid N, n north num numerical r radial direction ref refrigerant S, s south tp two-phase W, w west z axial direction Superscripts o previous instant p previous iteration

example, all those involving two-phase flow through tubes and ducts (double pipe, shell-and-tube and fin-and-tube condensers and evaporators, capillary tubes, etc.). Different analytical approaches as F-factor, e-NTU, P-parameter, etc. [1,2] are simple/limited methods to be used but they give a global approach for heat exchanger analysis. Stein [3,4] attempts an analytical solution to the two governing energy equations coupled at the boundary conditions in order to analyse double tube heat exchangers. Recent attempts use incremental loops that divided the heat exchanger into different thermal elements and simulate the behaviour of all these elements to find the overall performance. Most of these analysis for double pipe heat exchangers do not consider phase-change situations [5– 9]. More general approaches based on one-dimensional flow analysis [10] do not consider refrigerant mixtures. The objective of this work is to develop numerical criteria which allow the simulation, in both steady and transient state, of the thermal and fluid-dynamic behaviour of double-pipe evaporators and condensers (single-phase represents a particular case of this formulation) considering pure refrigerants and mixtures. The numerical solution in the

O. Garcı´a-Valladares et al. / International Journal of Refrigeration 27 (2004) 656–670

658

zones with fluid flow has been performed by discretization of the one-dimensional governing equations based on an efficiently coupled fully implicit step by step method over control volumes (CVs). The energy equation is written in terms of the enthalpy to give generality and to allow a suitable treatment of mixtures. The different flow regimes produced in both condensation and evaporation phenomenon, have been taken into account. The refrigerants and water thermophysical properties are evaluated locally using the REFPROP v.6.01 program [11] and NIST/ASME Steam properties program v.2.2 [12], respectively. The solid elements (tubes and insulation) are accurately solved considering multidimensional heat transfer effects. Details of the mathematical formulation and the different numerical techniques employed are shown first. Then, different numerical aspects are presented and finally comparisons of numerical with analytical and experimental results are shown.

2. Mathematical formulation 2.1. Two-phase flow inside ducts In this section, the mathematical formulation of the twophase flow inside a characteristic CV of a tube (single-phase flow, liquid or gas, represents a particular case) is presented (see Fig. 1). The mathematical formulation of the fluid flow is made neglecting the difference of the liquid and vapour temperatures in the subcooled boiling and post-dryout regime in the evaporating flow. Their effects are considered through the use of empirical correlations obtained from the available literature for the evaluation of the shear stress, the convective heat transfer and the flow structure, adequate to the flow patterns produced. A characteristic CV is shown schematically in Fig. 1, where ‘i’ and ‘i þ 1’ represent the inlet and outlet sections, respectively. Taking into account the characteristic geometry of ducts (diameter, length, roughness, angle,…), the governing equations have been integrated assuming the following assumptions: † One-dimensional flow: pðz; tÞ; hðz; tÞ; Tðz; tÞ;… † Fluid: pure and mixed substances. † Non-participant radiation medium and negligible radiant heat exchange between surfaces. † Axial heat conduction inside the fluid is neglected.

† Momentum: ½m _ g vg
iiþ1 þ ½m _ l vl
iiþ1 þ Dz

~_ ›m ›t

¼ 2½p
iþ1 i A 2 t~PDz 2 mg sin u

ð2Þ

† Energy: ~_ 1
iþ1 þ ½m _ g ðeg 2 el Þ
iþ1 þ ð~eg 2 e~ l Þ m½e i i þ mg

›mg ›t

›e~ g ›e~ ›p~ ›m þ ml l 2 ADz þ ð~el 2 e
l Þ ›t ›t ›t ›t

¼ q~_ PDz

ð3Þ

where f~ represents the integral volume average of a generic variable f over the CV and f
its arithmetic average between the inlet and outlet of the CV. The subscript and superscript in the brackets indicate ½X
iþ1 ¼ Xiþ1 2 Xi ; i.e. the difference between the i quantity X at the outlet section and the inlet section. In the governing equations, the evaluation of the shear stress is performed by means of a friction factor f : This factor is defined from the expression: t ¼ Fðf =4ÞðG2 =2rÞ:; where F is the two-phase factor multiplier. The onedimensional model also requires the knowledge of the twophase flow structure, which is evaluated by means of the void fraction e g : Finally, heat transfer through the duct wall and fluid temperature are related by the convective heat transfer coefficient a; which is defined as: a ¼ q_ wall =ðTwall 2 Tfluid Þ: The differentiation between the three main regions existing in both the condensation and the evaporation processes is given by the enthalpy, pressure and vapour quality. These conditions are: † Liquid region: when hðpÞ , hl ðpÞ; then xg ¼ 0: † Two-phase region: when hl ðpÞ # hðpÞ # hg ðpÞ; then 0 , xg , 1: † Vapour region: when hðpÞ . hg ðpÞ; then xg ¼ 1 where hl ðpÞ and hg ðpÞ represent the saturation liquid and gas enthalpy for a given pressure p: In the evaporation process, at least two more regions should be considered: subcooled boiling and post-dryout regions. 2.2. Heat conduction in the internal tube wall

The semi-integrated governing equations over the above mentioned finite CV, have the following form [13]: † Continuity: ½m _
iþ1 i þ

›m ¼0 ›t

ð1Þ

The conduction equation has been written assuming the following hypotheses: one-dimensional transient temperature distribution and negligible heat exchanged by radiation. A characteristic CV is shown in Fig. 2, where ‘P’ represents the central node, ‘E’ and ‘W’ indicate its neighbours. The CV-faces are indicated by ‘e’; ‘w’; ‘n’ and ‘s’:

O. Garcı´a-Valladares et al. / International Journal of Refrigeration 27 (2004) 656–670

Integrating the energy equation over this CV, the following equation is obtained: ðq~_ s Ps 2 q~_ n Pn ÞDz þ ðq~_ w 2 q~_ e ÞA ¼ m

›h~ ›t

hydraulic diameter. The convective heat transfer is calculated using the Monrad and Pelton correlation developed specifically for flow in annulus (cited by Jakob [16]).

ð4Þ

where q~_ s and q~_ n are evaluated using the respective convective heat transfer coefficient in each zone, and the conductive heat fluxes are evaluated from the Fourier law, that is: q~_ e ¼ 2le ð›T=›zÞe and q~_ w ¼ 2lw ð›T=›zÞw : 2.3. Heat conduction in the external tube wall and the insulation The external tube wall is solved in a similar way as described in Section 2.2 for the internal tube. The conduction equation for the insulation has been written assuming transient axisymmetric temperature distribution and negligible heat exchanged by radiation with the external ambient. Eq. (4) can also be applied over the annular CVs in the insulation (see Fig. 3). The north and south interfaces are evaluated from the Fourier law, except in the tube-insulation interface (where an harmonic mean thermal conductivity is used) and in the insulation-ambient interface (where the heat transfer by natural convection is introduced).

3. Evaluation of empirical coefficients The mathematical model requires local information about friction factor f ; convective heat transfer coefficient a and the void fraction e g This information is generally obtained from empirical or semi-empirical correlations. After comparing different empirical correlations presented in the technical literature, we have selected the following to obtain most of the results presented here: 3.1. Single-phase In the single-phase regions, the convective heat transfer coefficient is calculated using the Nusselt and the Gnieliski equations [14], for laminar and turbulent regimes, respectively. The friction factor is evaluated from the expressions proposed by Churchill [15]. 3.2. Annulus The friction factor in the annulus is calculated using the expressions corresponding to single-phase flow with the

Fig. 2. Heat fluxes in a control volume of a solid.

659

3.3. Condensing flow In the two-phase region, the convective heat transfer coefficient is calculated using the Dobson and Chato correlations [17] (for more detail of comparison between different correlation that appear in the literature, see Garcı´a – Valladares [18]) that employed two different expressions for the annular and wavy flow regions. The expression proposed by Soliman [19] for the Froude number is taken as a differentiation criterion between annular and wavy condensation. The void fraction is estimated from the semi-empirical equation of Premoli et al. [20]. The friction factor is calculated from the same equation as in the case of the single-phase flow using a correction factor (two-phase frictional multiplier) according to Friedel [21]. According to their authors and literature, these correlations can also be used for mixtures using the refrigerant mixture properties. 3.4. Evaporating flow In the case of subcooled boiling, the convective heat transfer coefficient and the friction factor are treated separately. For the convective heat transfer, the beginning of the subcooled boiling is estimated according to Frost and Dzakowic [22]. The method proposed by Bergles and Rohsenow [23] was used to consider the transition between pure liquid convection heat transfer and boiling heat transfer, the boiling heat transfer coefficient was obtained from the Forster and Zuber correlation [24]. For the friction factor, the point of net vapour generation is estimated according to Saha and Zuber [25]. The evaluation of the friction factor also needs a two-phase Reynolds number (evaluated from the homogeneous flow model described by Hewitt [26]) and a two-phase viscosity (evaluated according McAdams et al. [27]). The real vapour fraction proposed by Levy [28] is also employed. In the two-phase region, the convective heat transfer coefficient is evaluated using the model proposed by Zu¨rcher et al. [29,30] for horizontal plain tubes that incorporates the effects of local two-phase flow patterns, flow stratification, and partial dryout in annular flow. The local peak in the heat transfer coefficient versus vapor quality can be determined from the prediction of the location of onset of partial dryout in annular flow. The void fraction is evaluated according to a semi-empirical relationship proposed by Rouhani and Axelsson [31] (Zu¨rcher et al. [29,30] recommend use it in both the flow pattern map and the flow boiling heat transfer model). This relation is based on a drift flux model that includes the effects of mass velocity, surface tension and buoyancy. The friction factor is calculated in the same way as in the condensation twophase flow. According to their authors and literature, these

660

O. Garcı´a-Valladares et al. / International Journal of Refrigeration 27 (2004) 656–670

Fig. 3. Node distribution along the double-pipe heat exchanger.

correlations can also be used for mixtures using the refrigerant mixture properties.

4. Numerical resolution 4.1. Spatial and temporal discretization Fig. 3 shows the spatial discretization of a double pipe heat exchanger. The discretization nodes are located at the inlet and outlet sections of the CVs in the fluid flow zones, while the discretization nodes are centred in the CVs in the tube wall and insulation. Each one of the fluids has been divided into nz volumes (i.e. nz þ 1 nodes). The internal tube wall has been discretized into nz control volumes of length Dz: The external tube and the insulation are discretized into nz £ nr control volumes of length Dz and width Dr; where Dr ¼ ðD4 2 D3 Þ=2 for j ¼ 1 and Dr ¼ ðD5 2 D4 Þ=ð2ðnr 2 1ÞÞ for j . 1: The transitory solution is performed every time step Dt: Depending on the time evolution of the boundary conditions, a constant or variable value of Dt can be selected. 4.2. Two-phase flow inside ducts The discretized equations have been coupled using a fully implicit step by step method in the flow direction. From the known values at the inlet section and guessed values of the wall boundary conditions (obtained from the previous outer iteration), the variable values at the outlet of each CV are iteratively obtained from the discretized governing equations (see Section 4.2.1). This solution (outlet values) is the inlet values for the next CV. The procedure is carried out until the end of the tube is reached. 4.2.1. Discretization equations For each CV, a set of algebraic equations is obtained by a discretization of the governing Eqs (1) – (3) In the section, mathematical formulation, the governing

equations have been directly presented on the basis of the spatial integration over finite CVs. Thus, only their temporal integration is required. The transient terms of the governing equations are discretized using the following approximation: ›f=›T ø ðf 2 fo Þ=Dt; where f represents a generic dependent variable (f ¼ h; p, T, r,…); superscript ‘o’ indicates the value of the previous instant. The averages of the different variables have been estimated by the arithmetic mean between their values at the inlet and outlet sections, that is: f~i ø f
i ; ðfi þ fiþ1 Þ=2: Based on the numerical approaches indicated above, the governing Eqs (1)– (3) can be discretized to obtain the value of the dependent variables (mass flow rate, pressure and enthalpy) at the outlet section of each CV. The final form of the governing equations is given below. The mass flow rate is obtained from the discretized continuity equation, m _ iþ1 ¼ m _i 2

ADz ðr
2 r
otp Þ Dt tp

ð5Þ

where the two-phase density is obtained from: rtp ¼ e g rg þ ð1 2 e g Þrl : In terms of the mass flow rate, gas and liquid velocities are calculated as, " # " # mx _ g mð1 _ 2 xg Þ vg ¼ vl ¼ rg e g A rl ð1 2 e g ÞA The discretized momentum equation is solved for the outlet pressure, piþ1 ¼ pi 2

þ


_ 2 Dz f
m F P þ r
tp Ag sin u A 4 2r
tp A2


_ 2 m
_ o ½mðx _ g vg þ ð1 2 xg Þvl Þ
iþ1 m i þ Dt Dz

! ð6Þ

From the energy equation (3) and the continuity equation (1), the following equation is obtained for the outlet

O. Garcı´a-Valladares et al. / International Journal of Refrigeration 27 (2004) 656–670

detected. In order to evaluate the position of the transition point, two criteria have been tested [10]:

enthalpy: hiþ1 ¼

2_qwall 2 am _ iþ1 þ bm _ i þ cADz=Dt _ i þ r
otp ADz=Dt m _ iþ1 þ m

ð7Þ † Transition criterion 1: the transition point is assigned to the outlet section of the CV (assignation to the inlet section or to the middle section would be equivalent criteria). † Transition criterion 2: the CV is split into two control volumes. The length of the first CV is calculated from the energy equation, imposing the enthalpy condition to differentiate between one region and the other at the outlet section. The length of the second one is calculated by simple difference.

where a ¼ ðxg vg þ ð1 2 xg Þvl Þ2iþ1 þ g sin uDz 2 hi

b ¼ ðxg vg þ ð1 2 xg Þvl Þ2i 2 g sin uDz þ hi

2

c ¼ 2ð
pi 2 p
oi Þ 2 r
tpo ðhi 2 2h
oi Þ 2 ðr
v
2i 2 r
o voi Þ Temperature, mass fraction and thermophysical properties are evaluated using matrix functions of the pressure and enthalpy, refilled with the refrigerants properties evaluated using the REFPROP v.6.01 program [11] or NIST/ASME steam program v.2.2 [12]: T ¼ f ðp; hÞ; xg ¼ f ðp; hÞ; r ¼ f ðp; hÞ; …

ð8Þ

The above mentioned conservation equations of mass, momentum and energy together with the thermophysical properties, are applicable to transient two-phase flow. Situations of steady flow and/or single-phase flow (liquid or gas) are particular cases of this formulation. Moreover, the mathematical formulation in terms of enthalpy gives generality of the analysis (only one equation is needed for all the regions) and allows to deal with cases of mixtures of refrigerants. 4.2.2. Boundary conditions † Inlet conditions: the boundary conditions for solving a step by step method directly are the inlet mass flux ðm _ in Þ; pressure ðpin Þ and temperature ðTin Þ or mass fraction ðxgin Þ in the cases of inlet two-phase flow. From any of these values (temperature or mass fraction) and the pressure, enthalpy (our independent variable) is obtained. Other boundary conditions such as ðm _ in ; pout Þ or ðpin ; pout Þ or ðpin ; m _ out Þ can be solved using a Newton – Raphson algorithm. The method is based on an iterative process where the inlet mass flux or the inlet pressure is updated until global convergence is reached. † Solid boundaries: the wall temperature profile in the tube or the heat flux through the tube wall in each CV must be given. These boundary conditions are expressed in the energy equation in this compact form: q_ wall ¼ ð1 2 bÞ_qwall þ baðTwall 2 Tfluid Þ

661

ð9Þ

where b ¼ 1; if the boundary condition is the wall temperature, and b ¼ 0; if the heat flux is given. 4.2.3. Transition between regions Using the conditions of differentiation between regions mentioned in Section 2.1, the CV where transition occurs is

The transition criterion is very important because of the significant changes in both the heat transfer coefficient and friction factor between regions. Depending on the empirical correlations used, these coefficients can increase or decrease by factors of 10 or more. The importance of this criterion in obtaining accurate results using uniform and fewer control volumes represents significant important gains in computation cost for more complex systems (double tube, fin-andtube evaporators and condensers and other industrial equipment). Otherwise said, the transition criterion 2 will be used. 4.2.4. Solver In each CV, the values of the flow variables at the outlet section of each CV are obtained by solving iteratively the resulting set of algebraic equations (continuity, momentum, energy and state equations mentioned above) from the known values at the inlet section and the boundary conditions. The solution procedure is carried out in this manner, moving forward step by step in the flow direction. At each cross section, the shear stresses, the convective heat fluxes and the void fractions are evaluated using the empirical correlations obtained from the available literature (see Section 3). The transitory solution is iteratively performed at each time step. Depending on the time evolution of the boundary conditions, a constant or variable value of Dt can be selected. Convergence is verified at each CV using the following condition: ð12j

fpiþ1 2 fi jÞ , d Df

ð10Þ

where f refers to the dependent variables of mass flux, pressure and enthalpy; and fp represents their values at the previous iteration. The reference value Df is locally evaluated as fiþ1 2 fi : When this value tends to be zero, Df is substituted by fiþ1 and the numerator is substituted by fpiþ1 . 4.3. Heat conduction in the internal tube wall The following equation has been obtained for each node

O. Garcı´a-Valladares et al. / International Journal of Refrigeration 27 (2004) 656–670

662

of the grid: ai Twall;i ¼ bi Twall;iþ1 þ ci Twall;i21 þ di

ð11Þ

where the coefficients are, ai ¼

lw A lA ADz rcp þ e þ ðas Ps þ an Pn ÞDz þ Dz Dz Dt

bi ¼

le A Dz

ci ¼

lw A Dz

di ¼ ðas Ps T
s þ an Pn T
n ÞDz þ

4.5. Numerical algorithm At each time step solution process is carried out on the basis of a global algorithm that solves in a segregated manner the flow inside the tube, the flow inside the annulus and the heat conduction in the tubes and insulation. The coupling between the four main subroutines has been performed iteratively each time step following the procedure:

ADz 0 rcp Twall;i Dt

The coefficients mentioned above are applicable for 2 # i # nz 2 1; for i ¼ 1 and i ¼ nz adequate coefficients are used to take into account the axial heat conduction or temperature boundary conditions The set of heat conduction discretized equations is solved using the algorithm TDMA [32]. 4.4. Heat conduction in the external tube and insulation The following equation has been obtained for each node of the grid: aP Twall;i;j ¼ aE Twall;iþ1;j þ aW Twall;i21;j þ aN Twall;i;jþ1 þ aS Twall;i;j21 þ dP

ð12Þ

where the coefficients are, aW ¼

lw A ; Dz

aS ¼

ls Ps Dz ; Dr

aE ¼

le A ; Dz

d 0P ¼

aN ¼

take into account the axial heat conduction or temperature boundary conditions.

ln Pn Dz ; Dr

ADz rcP Dt

0 aP ¼ aW þ aE þ aN þ aS þ d 0P ; dP ¼ d 0P Twall;i;j

The coefficients mentioned above are applicable for 2 # i # nz 2 1 and 2 # j # nr 2 1; for the nodes in the extremes (see Fig. 3) the following considerations have been applied:

† For fluid flow inside the internal tube, the equations are solved considering the tube wall temperature distribution as boundary condition, and evaluating the convective heat transfer in each CV. † For fluid flow inside the annulus, the same process is carried out considering both wall temperatures (internal tube wall and external tube wall). † In the internal tube wall, the temperature distribution is re-calculated using the fluid flow temperature and the convective heat transfer coefficients evaluated in the preceding steps. † In the external tube and insulation, the temperature distribution is re-calculated using the fluid flow temperature and the convective heat transfer coefficients evaluated in the annulus and the external ambient. The governing equations corresponding to a steady state situation are the same equations shown above without considering the temporal derivative terms. This situation can be solved using a pseudotransient calculation with a guessed initial condition (for example, fluids at rest and uniform temperature in the whole domain). In transient situations, the initial conditions must be completely specified (i.e. the complete distribution of the variables in the whole domain). Sometimes these initial conditions are obtained from a steady state solution. The global convergence is reached when between two consecutive loops of the four main subroutines Eq. (10) is verified for all the CVs in the domain.

5. Numerical results † For j ¼ 1 forced convection is considered in the south face, tube thermal conductivity in east and west faces and insulation thermal conductivity in north face, all these evaluated to the mean temperature between the nodes that separated them. † For j ¼ nr natural convection with the ambient is considered (using the correlation developed by Raithby and Hollands [33]) and also the conduction through the insulation external part of thick equal to Dr=2; with thermal conductivity evaluated to the node temperature. † For i ¼ 1 and i ¼ nz adequate coefficients are used to

Based on the above mentioned mathematical model and numerical procedure, a code has been developed for the detail numerical simulation of double tube heat exchangers. The code has been carefully verified using, whenever possible, analytical solutions. For example, an analytical solution (based on the standard theory of heat transfer for fin surfaces with constant sectional area) is generated for the axial conduction in tubes. Extremes are considered isothermal. The internal and external fluids have different uniform temperatures with known heat transfer coefficients.

O. Garcı´a-Valladares et al. / International Journal of Refrigeration 27 (2004) 656–670

The comparison between numerical and analytical solutions has been made for the following case: Textreme;w ¼ 30 8C, Textreme;e ¼ 51 8C, Tfluid;ext ¼ 45 8C, Tfluid;int ¼ 28; 5 8C, afluid;ext ¼ 20 W m22 K 21, afluid;int ¼ 40 W m22 K21, Dint ¼ 0; 0079 m, Dext ¼ 0; 0095 m, lwall ¼ 400 W m 21 K21. Fig. 4 shows the comparison for the wall temperature between analytical and numerical results for the case of 40 CVs. An excellent degree of correlation can be observed in the numerical solution with a reasonable number of CVs. Another code verification test is here presented by comparison between the predicted and theoretical temperature results corresponding to an analytical solution of a steady state heat transfer problem postulated by Carslaw and Jaeger [34]. The problem developed consists of a hollow cylinder with a constant thermal conductivity of 0,7 W m21 K21 and with 0,05 m inner diameter and a 0,15 m outer diameter. It is considered that the fluid is circulating inside the cylinder at a constant inner temperature of 2 10 8C and with a constant heat transfer coefficient of 300 W m22 K21 while the outer cylinder face is maintained at a constant temperature of 30 8C. As can be seen in Fig. 5, an excellent agreement was found between the predicted wall temperatures with few CVs and the analytical solution. For this reason, in the numerical model few CVs have been used in the radial direction for the insulation. After code verification, a verification of the numerical solutions is carefully performed to assure their quality. Section 5.1 is devoted to the analysis of the main numerical criteria (transition criteria) and numerical parameters (time step, convergence criteria and number of grid nodes). After that, the mathematical formulation of the flow inside ducts with known wall heat fluxes distribution is validated by comparison of the numerical solutions with experimental data from the literature (Section 5.2). Finally, in Sections 5.3 and 5.4, different simulations of the whole heat exchanger working as evaporator or condenser (in parallel and counter flow arrangement) are compared with data obtained from the technical literature.

663

Fig. 4. Comparison for the wall temperature between analytical and numerical solutions using 40 grid nodes.

† Fluid: R134a † Boundary conditions: Tube wall: Twall ¼ 27 8C, Fluid(z ¼ 0): Tin ¼ 38 8C, pin ¼ 9 bar, m˙in ¼ 22,32 kg h21 (superheated vapour) Table 1 shows that the two criteria offer an asymptotic solution when the CV number is sufficiently increased. Transition criterion 2 shows much better performance than transition criterion 1. To reach the same precision in this case, 2000 CVs are required by transition criterion 1, while only 50 CVs are required by transition criterion 2. Grid independence results are obtained for transition criterion 2 with only 100 CVs. It is clear from these results that, in order to get the same order of accuracy, transition criterion 2

5.1. Analysis of different numerical aspects 5.1.1. Transition criteria Condensation inside an isothermal tube. The objective is to compare the performance of criteria 1 and 2 for the transition between single-phase and two-phase flow (see Section 4.2.3). The results presented in Table 1 for the steady state are the points zbc and zec where the condensation begins and ends, respectively, together with the outlet temperature Tout : The relative difference with respect to a grid independent numerical solution is also shown. This reference solution is obtained using transition criterion 2 with the numerical parameters nz ¼ 2000 and d ¼ 1027 : The case analysed corresponds to: † Geometry: L ¼ 3 m, j ¼ 1,5 £ 1026 m

D 1 ¼ 10 £ 1023 m,

u ¼ 0,

Fig. 5. Comparison for the wall temperature between analytical and numerical solutions using six grid nodes.

O. Garcı´a-Valladares et al. / International Journal of Refrigeration 27 (2004) 656–670

664

needs much less CVs (and consequently less CPU time) than transition criterion 1. This is a relevant aspect in more complex systems where the computational time is an important issue. Otherwise said, transition criterion 2 is used for the rest of the numerical results presented here. 5.1.2. Temporal discretization step, rate of convergence and grid density The objective of this point is to analyse the influence of the numerical parameters used in the model. Tables 2 and 3 show the numerical results obtained for the outlet temperature ðTout Þ; using different values of number of CVs ðnz Þ; temporal discretization step ðDtÞ and the rate of convergence required to finish the iterations ðdÞ: The analysed case begins with the steady state obtained above using transition criterion 2. The boundary conditions are: Tube wall: Twall ¼ 27 8C. Fluid ðz ¼ 0Þ : Tin ðtÞ ¼ Tt¼1 þ ðTt¼0 2 Tt¼1 Þexpð2t=t0 Þ pin ðtÞ ¼ pt¼1 þ ðpt¼0 2 pt¼1 Þexpð2t=t0 Þ; m _ in ¼ 22; 32 kg h21 ðsuperheated vapourÞ Tt¼0 ¼ 38 8C; Tt¼1 ¼ 50 8C; pt¼0 ¼ 9 bar; pt¼1 ¼ 10 bar; t0 ¼ 200 s: Tables 2 and 3 show the high accuracy of the numerical results even with a few numbers of CVs and relative high time steps. In general, good solutions are reached with 100 , nz , 200; Dt ¼ 10 s and d ¼ 1 £ 1023 : Table 2 shows additionally the dimensionless CPU time with respect to the case of nz ¼ 10 and d ¼ 1 £ 1021 considering different number of CVs and convergence rate. 5.2. Two-phase flow inside tubes with known wall heat fluxes. Experimental data by Jung and Didion [35] Fig. 6 graphically shows experimental results by Jung

and Didion [35] for steady state horizontal-flow boiling with a constant wall heat flux. The fluids used are R22 and a mixture 60% R12/40% R152a (molar fraction). The test section was made of two identical horizontal stainless steel tubes of 9,1 £ 1023 m internal diameter and 0,25 £ 1023 m wall thickness, 4 m length connected by a 1808 adiabatic Uturn bend (roughness is not specified but for the numerical simulation a standard value for commercially smooth tubes of j ¼ 1; 5 £ 1026 m is used). For this reason, two simulations are carried out for each test (two horizontal tubes), the first includes from 0 to 4 m (first tube) and the second from 4 to 8 m (second tube). The inlet conditions in the second section have been taken from the experimental results shown in tables included in the Jung and Didion’s report [35]. In the experimental results, four outside wall temperatures (bottom, top, right and left) in each station were measured. The experimental inner wall temperature was determined via one dimensional heat conduction equation, assuming uniform heat generation within the tube wall and an adiabatic condition on the outside of the tube. The pressure was also measured in different sections. For pure refrigerants, the fluid temperature was taken as the saturation temperature corresponding to the pressure at each thermocouple station. For mixtures, the fluid temperature depends on the quality as well as pressure; this temperature is evaluated from a specific enthalpy obtained by an energy balance. The numerical simulations were performed using as boundary conditions in the external wall of the tube, the experimental heat fluxes. The energy equation was used to take into account the heat conduction in the tube wall and to obtain the wall temperature distribution. The heat transfer coefficient was obtained using the Zu¨rcher et al. correlation [29,30]. In these cases, the correlation of Jung and Radermacher [36] was employed to calculate the friction factor. Two test cases (referred as 278 and 455 in the Jung and Didion’s report [35]) have been used for comparison with experimental data. The numerical parameters for each

Table 1 Numerical results using transition criteria 1 and 2 for R134a condensing inside an isothermal tube (Twall ¼ 27 8C, Tin ¼ 38 8C, pin ¼ 9 bar, m _ in ¼ 22; 32 kg h21 ) nz

Transition criterion 1

10 20 50 100 200 500 1000 2000

Transition criterion 2

zbc (m)

zec (m)

Tout (8C)

zbc (m)

zec (m)

Tout (8C)

0,300 (19,5%) 0,300 (19,5%) 0,300 (19,5%) 0,300 (19,5%) 0,270 (7,6%) 0,258 (2,8%) 0,255 (1,6%) 0,252 (0,4%)

3,000 (11,5%) 2,850 (5,9%) 2,876 (2,6%) 2,730 (1,5%) 2,700 (0,4%) 2,694 (0,1%) 2,694 (0,1%) 2,691 (0,0%)

35,03 (3,7%) 34,18 (1,2%) 34,05 (0,8%) 33,88 (0,3%) 33,81 (0,1%) 33,77 (0,0%) 33,77 (0,0%) 33,77 (0,0%)

0,249 (0,8%) 0,250 (0,4%) 0,251 (0,0%) 0,251 (0,0%) 0,251 (0,0%) 0,251 (0,0%) 0,251 (0,0%) 0,251

2,647 (1,6%) 2,674 (0,6%) 2,686 (0,1%) 2,689 (0,0%) 2,690 (0,0%) 2,690 (0,0%) 2,690 (0,0%) 2,690

33,56 (0,6%) 33,70 (0,2%) 33,76 (0,1%) 33,77 (0,0%) 33,78 (0,0%) 33,78 (0,0%) 33,78 (0,0%) 33,78

The relative differences with respect to the reference solution are given in parenthesis.

O. Garcı´a-Valladares et al. / International Journal of Refrigeration 27 (2004) 656–670

665

Table 2 Outlet temperature Tout ; (8C) obtained for different number of CVs and different time steps for R134a condensing inside an isothermal tube (Twall ¼ 27 8C, Tin ðtÞ ¼ Tt¼1 þ ðTt¼0 2 Tt¼1 Þexpð2t=t0 Þ; pin ðtÞ ¼ pt¼1 þ ðpt¼0 2 pt¼1 Þexpð2t=t0 Þ; Tt¼0 ¼ 38 8C; Tt¼1 ¼ 50 8C, pt¼0 ¼ 9 bar, _ in ¼ 22; 32 kg h21). The cpu time (tcpu*) is dimensionless with respect to time for the case of nz ¼ 10 and pt¼1 ¼ 10 bar, t0 ¼ 200 s, m d ¼ 1 £ 1021 ; nz

10 20 50 100 200 500 1000 2000

d ¼ 1 £ 1021

d ¼ 1 £ 1023

t¼0s

50

100

150

200

tcpup

t¼0s

50

100

150

200

tcpup

33,56 33,70 33,76 33,77 33,78 33,78 33,78 33,78

33,91 34,04 34,08 34,10 34,10 34,10 34,10 34,10

34,16 34,27 34,33 34,35 34,36 34,36 34,36 34,36

34,34 34,43 34,52 34,53 34,53 34,53 34,53 34,53

34,47 34,58 34,64 34,66 34,67 34,67 34,67 34,67

1 1,90 3,80 7,50 14,90 37,89 75,69 141,47

33,56 33,70 33,76 33,77 33,78 33,78 33,78 33,78

33,91 34,04 34,08 34,10 34,10 34,10 34,10 34,10

34,17 34,27 34,33 34,35 34,36 34,36 34,36 34,36

34,34 34,43 34,52 34,53 34,53 34,53 34,53 34,53

34,47 34,58 34,64 34,66 34,67 34,67 34,67 34,67

1,60 2,60 5,40 10,30 20,50 48,89 96,98 193,06

section are: nz ¼ 400 and d ¼ 1026 : The boundary conditions are: †Case 5.2a (section 1): pin ¼ 4; 285 bar, _ in ¼ 32; 38 g s21, q_ w ¼ 10 060 W m22. Tin ¼ 212; 8 8C, m (R22). Case 5.2a (section 2): pin ¼ 4; 189 bar, xgin ¼ 0; 131; m _ in ¼ 32; 38 g s21, q_ w ¼ 10 060 W m22. (R22). †Case 5.2b (section 1): pin ¼ 3; 91 bar, _ in ¼ 33; 05 g s21, q_ w ¼ 10 100 W m22. Tin ¼ 27; 0 8C, m (R12/R152a). Case 5.2b (section 2): pin ¼ 3; 8 bar, xgin ¼ 0; 149; m _ in ¼ 33; 05 g s21, q_ w ¼ 10 100 W m22 (R12/R152a). Fig. 6 shows a comparison between numerical and experimental fluid temperature, wall temperature and pressure. For both cases coincidence between the experimental and numerical results are quite reasonable. These results confirm the good prediction of the numerical model developed. A more comprehensive study about the thermal and fluid dynamic performance of the model for different two-phase

flows inside small ducts is presented in Garcı´a – Valladares et al. [37]. 5.3. Comparison of numerical results of double-pipe evaporators with experimental data 5.3.1. Comparison with experimental data by Takamatsu et al. [38] for a pure refrigerant In this section, two test cases corresponding to a doublepipe evaporator working with R22 are compared with the experimentally obtained values given by Takamatsu et al. [38]. The geometry of the element is (see Fig. 3): L ¼ 5; 52 m, D1 ¼ 7; 9 £ 1023 m, D2 ¼ 9; 5 £ 1023 m, D3 ¼ 16 £ 1023 m, u ¼ 0 and j ¼ 1; 5 £ 1026 m. The internal tube is of copper and the external tube has been assumed adiabatic (because no data is given in the article about the ambient temperature and the external tube and insulation thickness). The numerical parameters are: d ¼ 1026 and nz ¼ 400: The empirical correlations used in this case are the same

Table 3 Outlet temperature Tout ; (8C) obtained for different time steps Dt and instants with nz ¼ 200 for R134a condensing inside an isothermal tube (Twall ¼ 27 8C, Tin ðtÞ ¼ Tt¼1 þ ðTt¼0 2 Tt¼1 Þexpð2t=t0 Þ; pin ðtÞ ¼ pt¼1 þ ðpt¼0 2 pt¼1 Þexpð2t=t0 Þ; Tt¼0 ¼ 38 8C, Tt¼1 ¼ 50 8C, pt¼0 ¼ 9 bar, pt¼1 ¼ 10 bar, t0 ¼ 200 s, m _ in ¼ 22; 32 kg h21) Dt (s)

d ¼ 1 £ 1021 t¼0s

200 100 50 25 10 5 1 0.5 0.25

33,78 33,78 33,78 33,78 33,78 33,78 33,78 33,78 33,78

d ¼ 1 £ 1023 50

100

34,26 34,13 34,10 34,10 34,10 34,10 34,10

34,57 34,41 34,36 34,36 34,36 34,36 34,36 34,36

150

200

t¼0s

34,56 34,53 34,53 34,53 34,53 34,53 34,53

34,93 34,72 34,68 34,67 34,67 34,67 34,67 34,67 34,67

33,78 33,78 33,78 33,78 33,78 33,78 33,78 33,78 33,78

50

100

34,10 34,10 34,10 34,10 34,10 34,10 34,10

34,34 34,35 34,36 34,36 34,36 34,36 34,36 34,36

150

200

34,53 34,53 34,53 34,53 34,53 34,53 34,53

34,65 34,66 34,67 34,67 34,67 34,67 34,67 34,67 34,67

666

O. Garcı´a-Valladares et al. / International Journal of Refrigeration 27 (2004) 656–670

water increases with the progression of evaporation for counter flow, while it decreases for the parallel flow. Hence the heat fluxes increases with the quality for the counter flow and decreases for the parallel flow. This phenomenon has been reproduce in a reasonable good way in the numerical results. 5.3.2. Comparison with experimental data by Takamatsu et al. [39] for a nonazeotropic binary mixture In this section, two test cases corresponding to a doublepipe evaporator working with 51% R22/49% R114 (molar fraction) are compared with the experimentally obtained values given by Takamatsu et al. [39]. The geometry and numerical parameters used are the same mentioned above. The analysed situations correspond to: † Case 5.3.2a: flow arrangement: parallel flow * Tube fluid (R22/R114): pin ¼ 5; 8 bar, Tin ¼ 14 8C, Gin ¼ 299 kg m22 s21. * Annulus fluid (water): pin ¼ atmospheric pressure; Tin ¼ 61; 7 8C, Gin ¼ 212 kg m22 s21. † Case 5.3.2b: flow arrangement: counter flow * Tube fluid (R22/R114): pin ¼ 5; 8 bar, Tin ¼ 12 8C, Gin ¼ 299 kg m22 s21. * Annulus fluid (water): pin ¼ atmospheric pressure;Tin ¼ 57; 0 8C, Gin ¼ 212 kg m22 s21.

Fig. 6. Experimental and numerical pressure, fluid and wall temperature distribution: (a) case 5.2a, (b) case 5.2b (experimental results by Jung and Didion [35]).

mentioned in Section 3. The analysed situations correspond to: † Case 5.3.1a: flow arrangement: parallel flow * Tube fluid (R22): pin ¼ 11; 5 bar, Tin ¼ 22 8C, Gin ¼ 289 kg m22 s21. * Annulus fluid (water): pin ¼ atmospheric pressure; Tin ¼ 55; 5 8C, Gin ¼ 211 kg m22 s21. † Case 5.3.1b: flow arrangement: counter flow * Tube fluid (R22): pin ¼ 11; 4 bar, Tin ¼ 27 8C, Gin ¼ 290 kg m22 s21. * Annulus fluid (water): pin ¼ atmospheric pressure; Tin ¼ 56; 0 8C, Gin ¼ 214 kg m22 s21. Fig. 7 shows the numerical and experimental refrigerant, tube wall and annulus temperatures and refrigerant mass fraction distribution in the double-pipe evaporator for both cases. The temperature difference between refrigerant and

Fig. 8 shows the numerical and experimental refrigerant, tube wall and annulus temperatures and refrigerant mass fraction distribution in the double-pipe evaporator for both cases (parallel and counter flow). The figure also shows the refrigerant temperature calculated by Takamatsu et al. [39] using an equation of state. They obtained that the calculated temperature is about 3 8C higher than the measured temperature for higher quality (they concluded that the deviation depends on the overall composition). A reasonable good degree of correlation can be observed in both cases of the heat exchanger working with this binary mixture. The numerical results for the refrigerant temperature are closer to the refrigerant temperature calculated by Takamatsu et al. [39] using the equation of state. 5.3.3. Comparison with experimental data by Boissieux [40, 41] for a non azeotropic ternary mixture In this section, two test cases corresponding to a counter flow double-pipe evaporator working with R407C (23% R32/25% R125/52% R134a composition by mass) are compared with the experimentally obtained values given by Boissieux [42]. The geometry of the element is: L ¼ 4 m, D1 ¼ 8 £ 1023 m, D2 ¼ 9; 52 £ 1023 m, D3 ¼ 42 £ 23 10 m, D4 ¼ 52 £ 1023 m, u ¼ 0 and j ¼ 2 £ 1026 m. The internal tube is copper, the external one is steel and the insulation is armaflex. The numerical parameters are: d ¼ 1026 ; nr ¼ 5 and nz ¼ 400: The empirical correlations used in this case are the same mentioned in Section 3. The

O. Garcı´a-Valladares et al. / International Journal of Refrigeration 27 (2004) 656–670

667

Fig. 7. Experimental and numerical mass fraction and temperature distribution: (a) parallel flow, (b) counter flow (experimental results by Takamatsu et al. for R22 [38]).

Fig. 8. Experimental and numerical mass fraction and temperature distribution: (a) parallel flow, (b) counter flow (experimental results by Takamatsu et al. for R22/R114 [39]).

analysed situations correspond to:

again for both cases along the evaporator working with this ternary mixture.

† Case 5.3.3a: flow arrangement: counter flow. Tamb ¼ 19; 87 8C. * Internal fluid (R407C): pin ¼ 5; 94 bar, Tin ¼ 22; 2 8C, Gin ¼ 401; 9 kg m22 s21. * Annulus fluid (water): pin ¼ atmospheric pressure; Tin ¼ 54; 2 8C, Gin ¼ 111; 1 kg m22 s21. † Case 5.3.3b: flow arrangement: counter flow. Tamb ¼ 20; 16 8C. * Internal fluid (R407C): pin ¼ 6; 16 bar, Tin ¼ 20; 9 8C, Gin ¼ 288; 5 kg m22 s21. * Annulus fluid (water): pin ¼ atmospheric pressure; Tin ¼ 54; 7 8C, Gin ¼ 112; 6 kg m22 s21. Fig. 9 shows the numerical and experimental refrigerant and annulus temperature and refrigerant mass fraction distribution in the evaporator for these two cases. A reasonable good degree of correlation can be observed

5.4. Comparison of numerical results of double-pipe condensers with experimental data 5.4.1. Comparison with experimental data by Boissieux [40, 43] for a near azeotropic ternary mixture In this section, two test cases corresponding to a counter flow double-pipe condenser working with R404A (44% R125/52% R143a/4% R134a composition by mass) are compared with the experimentally obtained values given by Boissieux [42]. The geometry and numerical parameters are the same mentioned above. The empirical correlations used in this case are the same mentioned in Section 3, with the exception of the heat transfer coefficient in the annulus due to the fluid is water/glycol (40%) instead of water. The correlation used was formulated by Boissieux [40] in his calibrations for condensation cases, where Z is the

O. Garcı´a-Valladares et al. / International Journal of Refrigeration 27 (2004) 656–670

668

measurement position along the annulus in the refrigerant flow direction:

aa ¼

la Daa Nu þ Z Dh Dz

with Nu ¼ 1; 7 £ 1023 Re þ 85; 784 and 2; 1023 £ 1024 Re 2 1; 531: The analysed situations correspond to:

ð13Þ Daa =Dz ¼

† Case 5.4.1a: flow arrangement: counter flow. Tamb ¼ 19; 65 8C. * Internal fluid (R404A): pin ¼ 12; 22 bar, Tin ¼ 25; 1 8C, Gin ¼ 409; 8 kg m22 s21. * Annulus fluid (water/glycol): pin ¼ atmospheric pressure; Tin ¼ 26; 8 8C, Gin ¼ 307; 4 kg m22 s21. † Case 5.4.1b: flow arrangement: counter flow. Tamb ¼ 19; 15 8C. * Internal fluid (R404A): pin ¼ 13; 3 bar, Tin ¼ 36; 4 8C, Gin ¼ 208; 9 kg m22 s21.

Fig. 10. Experimental and numerical mass fraction and temperature distribution: (a) case 5.4.1a, (b) case 5.4.1b (experimental results by Boissieux for R404A [42]). *

Annulus fluid (water/glycol): pin ¼ atmospheric pressure, Tin ¼ 27; 3 8C, Gin ¼ 207 kg m22 s21.

The simulation results of distribution of temperatures and mass fraction in the heat exchanger are shown in Fig. 10 that agrees well with the experimental data profiles, which verified the reliability of the model. 5.4.2. Comparison with experimental data by Boissieux [40, 43] for a non azeotropic ternary mixture In this section, two test cases corresponding to a counter flow double-pipe condenser working with R407C are compared with the experimentally obtained values given by Boissieux [42]. The geometry, numerical parameters and empirical correlations used are the same mentioned above. The analysed situations correspond to: Fig. 9. Experimental and numerical mass fraction and temperature distribution: (a) case 5.3.3a, (b) case 5.3.3b (experimental results by Boissieux for R407C [42]).

† Case 5.4.2a: flow arrangement: counter flow. Tamb ¼ 15; 33 8C.

O. Garcı´a-Valladares et al. / International Journal of Refrigeration 27 (2004) 656–670

669

6. Concluding remarks A numerical model of the thermal and fluid dynamic behaviour of double-pipe condensers and evaporators has been developed by means of a transient one-dimensional analysis of the flows together with a detail analysis of the heat conduction in solids. The empirical coefficients used in the model to evaluate the shear stress, heat flux and twophase structure have been chosen after a comparison of different options proposed in the technical literature. The simulation has been implemented on the basis of an implicit step by step numerical scheme for the fluid flow inside the tube and annulus, and an implicit central difference numerical scheme in the solids. The different zones that form the double-pipe heat exchanger are solved iteratively in a segregated manner until convergence is reached. Different numerical aspects related to the global algorithm, grid density, rate of convergence, and transition criteria have been presented. Grid independent results are more easily obtained using the second transition criterion than the first one. This result permits more accurate results with less CVs and consequently less CPU time. Comparison with experimental data has firstly performed for the twophase flow inside ducts of pure refrigerants and mixtures. After that, the whole heat exchanger working with different fluids as condenser or evaporator (in parallel and counter flow arrangement) is solved under different working conditions. The results obtained are compared with experimental data from the technical literature showing a good degree of agreement.

Fig. 11. Experimental and numerical mass fraction and temperature distribution: (a) case 5.4.2a, (b) case 5.4.2b (experimental results by Boissieux for R407C [42]).

Internal fluid (R407C): pin ¼ 10; 44 bar, Tin ¼ 30; 6 8C, Gin ¼ 397; 9 kg m22 s21. * Annulus fluid (water/glycol): pin ¼ atmospheric pressure, Tin ¼ 2 7,1 8C, Gin ¼ 441,3 kg m22 s21. † Case 5.4.2b: flow arrangement: counter flow. Tamb ¼ 17; 24 8C. * Internal fluid (R407C): pin ¼ 11; 52 bar, Tin ¼ 47; 6 8C, Gin ¼ 175; 1 kg m22 s21. * Annulus fluid (water/glycol): pin ¼ atmospheric pressure, Tin ¼ 211; 5 8C, Gin ¼ 265; 5 kg m22 s21. *

Fig. 11 compares the simulation results with the experimental data of Boissieux [42] in a double-pipe condenser for both cases. There are fair agreements between both results with respect to the shape of the temperature and mass fraction distributions and also with the position were the condensation begins and ends.

Acknowledgements The authors thank Xavier Boissieux for the data and technical support provided. This study has been supported by the ‘Consejo Nacional de Ciencia y Tecnologı´a’ CONACYT-Me´xico and the ‘Generalitat de Catalunya’.

References [1] Kern DQ. Process heat transfer. McGraw-Hill; 1950. [2] Kakac¸ S, Liu H. Heat exchangers: selection, rating, and thermal design. Boca Raton, FL: CRC Press; 1998. [3] Stein RP. Heat transfer coefficient in liquid metat concurrent flow double pipe heat exchangers. Chem Engng Prog Symp Ser 1965;59:64–75. [4] Stein RP. The Graetz problem in concurrent flow double pipe heat exchangers. Chem Engng Prog Symp Ser 1965;59: 78– 87. [5] Scofano NF, Cotta RM. Counterflow double-pipe heat exchanger analysis using mixed lumped-differential formulation. Int J Heat Mass Transfer 1992;35(7):1723–31. [6] Lachi M, El Wakil N, Padet J. The time constant of double pipe and one pass shell and tube heat exchangers in the case of

670

[7]

[8]

[9]

[10]

[11]

[12]

[13]

[14]

[15] [16] [17] [18]

[19]

[20]

[21]

[22]

[23]

[24] [25]

[26]

O. Garcı´a-Valladares et al. / International Journal of Refrigeration 27 (2004) 656–670 varying fluid flow rates. Int J Heat Mass Transfer 1997;40(9): 2067–79. ¨ nal A. Theoretical analysis of triple concentric-tube heat U exchangers. Part I: mathematical modelling. Int Commun Heat Mass Transfer 1998;25(7):949–58. Abdelghani-Idrissi MA, Bagui F. Counterflow double-pipe heat exchanger analysis subjected to flow-rate step change. Part I: new steady-state formulation. Heat Transfer Engng 2002;23(5):4 –11. Abdelghani-Idrissi MA, Bagui F, Estel P. Counterflow doublepipe heat exchanger analysis subjected to flow-rate step change. Part II: analytical and experimental transient response. Heat Transfer Engng 2002;23(5):12–24. Escanes F, Pe´rez-Segarra CD, Oliva A. Thermal and fluiddynamic behaviour of double-pipe condensers and evaporators—a numerical study. Int J Numer Methods Heat Fluid Flow 1995;5(9):781–95. REFPROP v.6.01. NIST thermodynamic properties of refrigerants and refrigerant mixtures database, Standard Reference Data Program, USA; 1998. NIST/ASME Steam v2.2. Formulation for General and Scientific Use, NIST Standard Reference Database Number 10; 1996. Garcı´a-Valladares O, Pe´rez-Segarra CD, Oliva A. Numerical simulation of capillary tube expansion devices behaviour with pure and mixed refrigerants considering metastable region. Part I: mathematical formulation and numerical model. Appl Therm Engng 2002;22(2):173–82. Gnielinski V. New equations for heat and mass transfer in turbulent pipe and channel flow. Int Chem Engng 1976;16: 359–68. Churchill SW. Frictional equation spans all fluid flow regimes. Chem Engng 1977;84:91–2. Jakob M. Heat transfer. New York: Wiley; 1949. p. 552. Dobson MK, Chato JC. Condensation in smooth horizontal tubes. J Heat Transfer 1998;120(1):193–213. Garcı´a-Valladares O. Review of in-tube condensation heat transfer correlations for smooth and microfin tubes. Heat Transfer Engng 2003;24(4):6–24. Soliman HM. On the annular to wavy flow pattern transition during condensation inside horizontal tubes. Can J Chem Engng 1982;60:475–81. Premoli A, Francesco DD, Prina A. A dimensional correlation for evaluating two-phase mixture density. La Termotecnia 1971;25(1):17–26. Friedel L. Improved friction pressure drop correlation for horizontal and vertical two-phase pipe flow. European TwoPhase Flow Group Meeting, Ispra, Italy, Paper E2; 1979. Frost W, Dzakowic GS. An extension of the method of predicting incipient boiling on commercially finished surfaces. ASME/AIChe Heat Transfer Conference, paper 67-ht-61; 1967. p. 1–8. Bergles AE, Rohsenow WM. The determination of forcedconvection surface-boiling heat transfer. J Heat Transfer 1964; 86C:365–72. Forster HK, Zuber N. Dynamics of vapour bubble growth and boiling heat transfer. AIChe J 1955;1(4):531–5. Saha P, Zuber N. Point of net vapour generation and vapour volumetric void fraction. Proceedings of 15th International Heat Transfer Conference. 1974. p. 175 –9. paper B4.7. Hewitt GF. Gas–liquid flow. In: Schlu¨nder EU, Hewitt GF,

[27]

[28]

[29]

[30]

[31]

[32] [33]

[34] [35]

[36]

[37]

[38]

[39]

[40]

[41]

[42] [43]

Beaton CF, Bell KJ, Butterworth D, Gnielinski V, Morris M, Schmidt FW, Taborek J, Vilemas J, Johnson R, editors. Heat exchanger design handbook, sec.2.7.3, 2002. Washington: Hemisphere Publishing Corporation; 2002. p. 2.7.3-10–2.7.3-11. McAdams WH, Woods WK, Heronman LC. Vaporization inside horizontal tubes. Part 2: benzine–oil mixtures. ASME Trans 1942;64:193–200. Levy S. Forced convection subcooled boiling prediction of vapour volumetric fraction. Int J Heat Mass Transfer 1967;10: 951 –65. Zu¨rcher O, Thome JR, Favrat D. Evaporation of ammonia in a smooth horizontal tube: heat transfer measurement and predictions. J Heat Transfer 1999;121(1):89–101. Zu¨rcher O, Thome JR, El Hajal J. Two-phase flow pattern map for evaporation in horizontal tubes: latest version. First International Conference on Heat Transfer, Fluid Mechanics, and Thermodynamics, Kruger Park, South Africa, TJ2; 2002. Rouhani Z, Axelsson E. Calculation of volume void fraction in the subcooled and quality region. Int J Heat Mass Transfer 1970;13:383–93. Patankar SV. Numerical heat transfer and fluid flow. New York: McGraw-Hill; 1980. Raithby GD, Holland GT. In: Hartnett JP, Irvine TF Jr., editors. A general method of obtaining approximate solutions to laminar and turbulent free convection problems. Advances in heat transfer, vol. 11. New York: Academic Press; 1975. p. York. Carslaw HS, Jaeger JC. Conduction of heat in solids, 2nd ed. Oxford, UK: Clarendon Press; 1959. Jung D, Didion DA. Horizontal-Flow Boiling Heat Transfer Using Refrigerant Mixtures. Electrical Power Research Institute (EPRI), ER-6364, Research Project 8006-2; 1989. Jung D, Radermacher R. Prediction of pressure drop during horizontal annular flow boiling of pure and mixed refrigerants. Int J Heat Mass Transfer 1989;32(12):2435–46. Garcı´a-Valladares O, Pe´rez-Segarra CD, Oliva A. Numerical simulation of capillary tube expansion devices behaviour with pure and mixed refrigerants considering metastable region. Part II: experimental validation and parametric studies. Appl Therm Engng 2002;22(4):379 –91. Takamatsu H, Momoki S, Fujii T. A correlation for forced convective boiling heat transfer of pure refrigerants in a horizontal smooth tube. Int J Heat Mass Transfer 1993;36(13): 3351–60. Takamatsu H, Momoki S, Fujii T. A correlation for forced convective boiling heat transfer of nonazeotropic refrigerant mixture of HCFC22/CFC114 in a horizontal smooth tube. Int J Heat Mass Transfer 1993;36(14):3555–63. Boissieux X. Two-phase local heat transfer correlations for non-ozone depleting refrigerant –oil mixtures. PhD thesis, University of Brighton, UK; 1998. Boissieux X, Heikal MR, Johns RA. Two-phase heat transfer coefficients of three HFC refrigerant inside a horizontal smooth tube. Part I: evaporation. Int J Refrigeration 2000; 23(4):269–83. Boissieux X. private communication. Boissieux X, Heikal MR, Johns RA. Two-phase heat transfer coefficients of three HFC refrigerant inside a horizontal smooth tube. Part II: condensation. Int J Refrigeration 2000; 23(5):345–52.