A numerical and experimental study of turbulent flow through the evaporator coil in an air-conditioning unit

ELSEVIER

Int J. Refrig. Vol. 19, No. 6, pp. 369-381, 1996 Copyright © 1996 Published by Elsevier Science Ltd and IIR Printed in Great Britain. All rights reserved PII: S0140-7007(96)00024-2 0140-7007/96/$15.00 + .00

A numerical and experimental study of turbulent flow through the evaporator coil in an air-conditioning unit Z. G. Xu P r o c e s s E n g i n e e r i n g D i v i s i o n , Silsoe R e s e a r c h I n s t i t u t e , B e d f o r d s h i r e , U K

D. H. T. Gotham and M. W. Collins Thermo-Fluids Engineering Research Centre, City University, London UK

J. E. R. Coney, C. G. W. Sheppard and S. Merdjani D e p a r t m e n t o f M e c h a n i c a l E n g i n e e r i n g , U n i v e r s i t y o f Leeds, Leeds, U K Received 28 N o v e m b e r 1995; revised 26 M a r c h 1996

Computational fluid dynamics (CFD) and heat transfer has matured to the stage where state-of-the-art commercial codes can address real engineering problems. However, the advances in finite volume numerical capability which permit this, such as full three-dimensionality, generalised coordinate systems and multiblocking, require convincing validation procedures to demonstrate quantitative engineering potential. In this paper a two-dimensional incompressible turbulent flow in an industrial packaged air-conditioning system is investigated numerically, with specific reference to the air flow through the tube banks of the evaporator coil. The system modelled comprises an outer casing containing the coil which has three rows of 20-mm diameter tubes in staggered formation, each row having tubes, and a condensate tray with complex geometry. The analysis is based on the Reynolds-averaged Navier-Stokes equations in conjunction with the standard k-e turbulence model. The finite-volume method (FVM) is used to solve the governing partial differential equations on a non-staggered grid. In order to treat the detailed fin tube geometry, a novel method called 'masking' has been developed to construct body-fitted curvilinear grids quickly and efficiently. In the computation, several differencing schemes of various orders are outlined and their accuracy are examined. In a parallel experimental investigation conducted on the actual air-conditioning unit, both the mean velocity profiles and turbulence properties of the flow were obtained from triple hot-wire anemometry measurements. These were then used to compare with the computational results and validate the mathematical models developed. It is found that the predictions were generally in good agreement with the measurements, especially for the higher-order differencing scheme. An immediate practical effect of the study was that it revealed the flow disturbance resulting from the condensate tray. This has been suspected but its magnitude not anticipated. The consequent amendment of the design of this component in the commercial unit led to an improved performance of the evaporator. Copyright © 1996 Published by Elsevier Science Ltd and IIR

(Keywords:heat transfer;masstransfer;air;battery;finnedtube;flow;turbulentflow;calculation;measurement)

Etude num6rique et exp6rimentale d'un 6coulement turbulent dans le serpentin de l'6vaporateur d'une installation frigorifique La dynamique quantitative des fluides ou CFD et le transfert de chaleur ont ~volu~ de telle fa¢on que les codes commerciaux permettent de r~soudre des problbmes rOel d'ing~nierie. Cependant, les progrbs dans l'analyse num~rique ~ volumes finis, telles que la tri-dimensionalit~ complbte, les systbmes coordonnOs g~n~ralis~s et le blocage multiple exigent des procOdures de validation convaincantes qui d~montrent le potentiel de l'ingknierie quantitative. Dans l'article, on ~tudie num~riquement un ~coulement turbulent incompressible et bidimensionnel d'un systbme de conditionnement d'air autonome industriel, en se r~f~rant plus sp~cifiquement h l'~coulement d'air h travers les nappes de tubes du serpentin de l'evaporateur. Le systbme modOlis~ comprend une enveloppe ext~rieure contenant le serpentin qui pr~sente trois rang~es de tubes de 20 mm de diambtre places en quinconce, chaque rang~e ayant dix tubes et le plateau du condensat ayant une gOom~trie complexe. L'analyse est fond~e sur les ~quations de Navier-Stokes, lesquelles s'inspirent des ~quations de Reynolds, en conjonction avec le modble de turbulence type k - e . On utilise la m~thode du volume fini ( F V M ) pour r~soudre les ~quations diff~rentielles partielles et directrices sur une grille dont la disposition n 'est pas en quinconce. Afin d'~tudier de fa~on d~taill~e la g~om~trie des ailettes du tube, une nouvelle mdthode appel~e 'masquage' a ~t~ raise au point pour concevoir rapidement et avec efficacit~ des grilles curvilinOaires adaptables. Dans le calcul, on insiste sur plusieurs schemas diff~rents de plusieurs ordres, et on examine leur precision. Dans une ~tude exp~rimentale parallble, men~e sur une installation de conditionnement d'air rOelle, les profils de vitesse moyenne et les propi~tds de la turbulence de l'dcoulement ont ~t~ obtenus h partir de trois mesures par

369

Z.G. Xu et al.

370

Nomenclature

u,. ui xi P k B S Pe

Time-averaged mean velocity Fluctuating velocity Coordinate in ith direction Static pressure Turbulence kinetic energy Body force Source term Peclet number

Greek letters E

F #~ ~~b P

dissipation rate of turbulence kinetic energy Diffusion coefficient Effective viscosity Shear stress General variable Fluid density

andmomOtrie h ills chauds. On les a ensuite comparOs avec les rOsultats de calcul et on a validO les modkles mathOmatiques. On observe que les provisions ont Ot~ gdn~ralement en bonne correlation avec les mesures, surtout pour le schema different d'ordre supOrieur. L'~tude a immOdiatement r~v~l~ que la perturbation de l~coulement ~tait provoquOe par le plateau du condensat. Cela avait ~tO soupf.onn~ mais l'ordre de grandeur n'avait pas OtOprOvu. La modification ult~rieure du composant dans des installations commerciales a permis d'amOliorer la performance de l'Ovaporateur. Copyright © 1996 Published by Elsevier Science Ltd and IIR

(Mots-cl6s:transfertde chaleur;transfertde masse;air batteri6;tubeailette;ecoulement;regimeturbulent;caleul;mesure)

In the design of packaged air-conditioning and refrigeration equipment, it is found that packaging and installation requirements can often lead to some non-ideal flow patterns (maldistribution) at the approach to and the discharge from heat exchangers and other components. Of particular interest to a design engineer is the maldistribution through the tube banks of the evaporator coil, which can result in a severe reduction in the heat exchanger performance. Moreover, it can also cause condensate carry-over at the discharge from the evaporator and as such affect the formation of frost on the coils. In order to minimise these undesirable side effects, it is important to establish under what conditions such flow maldistributions occur and, ultimately, to determine what steps can be taken in the design so that the problem can be resolved. For instance, the flow pattern can be improved by modifying the size and position of the components which have caused flow maldistribution. However, this requires a detailed knowledge and thorough understanding of turbulent flows within an air-conditioning unit. Packaged equipment is normally fitted with a variety of components, such as a heat exchanger, humidifier and centrifugal fan all within a rectangular enclosure. This arrangement, through compact, can lead to a non-ideal internal geometry and a correspondingly complex flow pattern. In general, the flow field of interest is threedimensional and turbulent. There exist regions of flow separation and recirculation in the wake of these components. The two factors combined have so far precluded treatment using traditional methods of analysis. As a result, previous investigations of flow phenomena in such equipment have been largely confined to experimental tests carried out on the actual unit, or in a wind tunnel. However, experiments are often very costly, labourious and time-consuming to perform. Moreover, they can only provide data for a limited number of configurations and often cannot be performed with the degree of spatial and temporal resolution required to help unravel the nature of the process involved.

In recent numerical work 1'2, it was demonstrated that the complex three-dimensional turbulent flow with recirculation in a typical packaged air-conditioning unit can be predicted using a computational fluid dynamics (CFD) approach. The method solves numerically the set of partial differential equations (PDEs) governing the flow. In the model, nearly all the main components in an air-conditioning unit were included. Among them were the evaporator coil, the humidifier, the electrical control box, the motor and the fan. In addition, the rotational motion of the fan impeller was considered and its effects on the flow field were accounted for by including the centrifugal and Coriolis forces as a source term in the relevant conservation equations. To validate the model developed, the predicted mean velocity profiles downstream of the coil were compared with those obtained experimentally3 and were found to be in reasonably good agreement. In the numerical work 1'~ the turbulent flow through the tube banks of an evaporator coil was not modelled in detail. Instead, the porous medium approach was adopted in which the coil was treated as a uniformly distributed flow resistance unit. To describe the flow in a porous medium, two useful concepts have often been employed. The first of these, the porosity factor, is defined as the fraction of the control volume available to the fluid flow. It is purely a geometrical parameter and can be determined from the specified fin and tube geometry. The second, known as the Darcy coefficient, is used to account for the additional pressure loss across the coil due to the presence of solid objects. The main advantage of using the porous medium approach is that it permits the choice of a simple computational grid which does not have to relate to the exact fin-tube geometry. Therefore, it is possible to obtain the required solution without excessive demand on either computer memory or processor speed. The approach was originally proposed by Patankar and Spalding4 and later was used by Shim et al. 5 to model fuel assemblies in a reactor core. On the other hand, the porous medium approach

371

Turbulent flow in air-conditioning unit introduces sources of uncertainty and potential errors. It is of only limited use in that the Darcy coefficient is an empirical parameter. In the simulation, it was determined using the condition that the predicted pressure drops across the coil should agree with the values obtained from wind-tunnel experiments. Once the flow condition under which the experimental values are obtained changes, a new Darcy coefficient might have to be specified. Due to the complex nature of the flow and the system configuration for a complete packaged air-conditioning unit, this 'calibration' process could be very time-consuming indeed. Due to the lack of computing power often encountered in industrial and engineering applications, it is still not practical, nor feasible, to solve the exact full-field threedimensional turbulent flow in packaged air-conditioning units. To bridge the gap between a full-field calculation and the porous medium simplification, an alternative modelling strategy was proposed. It was based on the 'fine' modelling of the turbulent flow through individual tubes of the coil, albeit two-dimensional. Here, a system, consisting of the individual tubes of the evaporator coil in the outer casting alone, was modelled in detail. It was observed from previous numerical work that the flow through the coil is nearly two-dimensional. Therefore the central 'symmetrical' plane was selected. In this way it was possible to produce the required pressure loss data across the coil under a wide range of flow conditions based on the exact fin-tube geometry. This information could then be used to calibrate the Darcy coefficient in three-dimensional simulations in which the evaporator coil is modelled as a porous medium. The main purpose

of this paper is, therefore, to determine, using the finite volume methodology, detailed turbulent flow patterns through the tube banks of an air-conditioning unit. The theoretical model developed for the system, together with a two-equation k-E turbulence model, is presented. Several differencing schemes, employed for the convection terms, are described in detail. Their convergence and solution accuracy are examined in relation to experimental hot-wire anemometry measurements.

Experimental investigation Description of the rig Illustrated in Figure 1 is a sectional view of the Airedale International Air Conditioning EHA5 horizontal flow,

Coil Air Flow Condensatetra~~,..

Figure 2 Condensate tray configuration Figure 2 Configuration du plateau des condensats

Figure 1 Sectional view of AIAC type EHA5 air-conditioning unit. 1, Fan; 2, motor; 4, fan drive; 8, coil; 18, filters; 26, humidifier; 28, expansion valve; 30, power panel; 43, electrical heater; 49, condensate tray Figure 1 Vue en coupe de l'installation de conditionnement d'air AIA C type EHA5. 1, ventilateur ; 2. rnoteur; 4, entrafnement du ventilateur ; 8, serpentin; 18, filtres; 26, humidificateur; 28, d~tendeur; 30, dispositif #lectrique; 43, chauffage ~lectrique; 49, plateau des condensats

372

Z.G. Xu et al.

packaged air-conditioning unit, under test. The overall dimensions of the unit are 1.6 x 0.7 × 0.93 m. Its main component is an evaporator, downstream of which is a centrifugal fan. Between the evaporator and the fan, a humidifier is positioned. There is a filter upstream of the evaporator. The evaporator is of the finned-tube type and consists of 120, 10-mm outside diameter tubes arranged vertically in five circuits; it is placed normal to the main direction of the air flow. Beneath the evaporator is a condensate tray, their relative positions being shown in Figure 2. In order to provide data for the validation of the numerical model developed, an experimental test rig, based on the EHA5 unit, was commissioned jointly by Airedale International Air Conditioning and the University of Leeds. Using this rig, the air flow was investigated and velocity profiles obtained, employing triple hot-wire anemometry techniques. These were then used for comparison with the calculated velocity profiles.

I D'SA55M0111 .

O;FgT

DiSA 55M01 ] [ DiSA 55M01 ] .

.

OFFSET

OFFSET ] BOX

r

CHJ ~EL

Triple hot-wire aneomometry techniques Triple hot-wire anemometry techniques were resorted to rather than single wire techniques because of the high turbulence levels expected in the flow. If a single hot-wire probe is used, the mean velocity and components of Reynolds stress tensors may be deduced from a number of orientations of the wire6;' a set of simultaneous equations is obtained based on separate time-averaged measurements. However, this solution requires that high order terms in the equation be ignored, which incurs considerable errors if the turbulence levels are high. If a triple hot-wire probe is employed, having the wires very close to each other at a point, three simultaneous signals are obtained from which the orthogonal components of instantaneous velocity may be obtained. As described by Frota 7, the turbulent kinetic energy and Reynolds stresses may then be deduced. Shown in Figure 3 is the schematic for the offset technique used in the investigation, which gives good accuracy for highly turbulent flows. In this technique, there is separation into mean and fluctuating parts of the instantaneous signals from the three orthogonal wires; this enables use of the full range of ADCs to improve the signal-to-noise ratio of the turbulent data. The signals from the triple hot-wire probe (DANTEC type 55P9) were passed to three measuring units (DANTEC type 55M10), the outputs from which were separated in the three offset boxes. These units contained a low-pass and a high-pass filter, arranged in parallel; the former employed a cut-off frequency above which any change in voltage was considered as part of the fluctuation. In order to take full advantage of the available +10V ADC units, the fluctuation signals were then amplified. The mean signals E and their fluctuating components were then passed to their respective ADC units. For each measurement point within the flow, the three anemometer signals were sampled at a frequency of 2000 Hz (0.5 ms) with 1400 items of data per sample; this was repeated 10 times giving a total sampling time of 7 s. From this data, instantaneous velocities were derived for the three orthogonal directions, which were then used to obtain the Reynolds stresses. Signal organisation and processing were by a single program, TRIPLE 9, which permitted the immediate processing of the signal or the

tuaring the signals

Figure 3 Triple hot-wire signal processing Figure 3 Triple traitement du signal par ills chauds

storage on magnetic tape of the anemometer voltage forwarded by the ADCs. The above procedure is fully described by Merdjani 8.

System description Velocity profiles were obtained 110ram downstream of the evaporator and immediately behind the condensate tray. There are a number of reasons for selecting this as the measurement plane. In the first place, it was anticipated that the flow pattern in the region is the most complex, due to the strong interaction between the

Turbulent flow in air-conditioning unit flow from the heat exchanger and condensate tray geometry. Secondly, the data obtained can be used to testify the validity of the assumption by G o t h a m 1'2 that the heat exchanger can be modelled as a porous medium. It is important to have some idea whether the fin-tube geometry would still affect the flow field at this location. Previously, single wire anemometry measurements had been made in this position. However, it was decided to investigate further, using a triple wire, since this provided velocities in the X, Y and Z directions, whereas the single wire gave velocities representing a combination of two out of the three components. A preliminary numerical investigation by the authors indicated that it was not possible to reproduce the exact geometry of the evaporator coil as used in the EHA5 unit, given the limited computer resources available. This is due to the fact that to describe the geometry of each tube properly, a minimum of 3 x 3 computational cells were required. Consequently, a compromise was made in which a coil consisting of only 3 x 10 20-mm diameter tubes was chosen for the model. It was calculated that this would present the same amount of flow blockage as the original arrangement. In this sense, the two geometries are equivalent. Furthermore, it can also introduce the main features of the coil into the computational analysis without the tubes losing their identity in the vector plots of the resulting flow pattern. As shown in Figure 4 the current system to be modelled comprises an evaporator coil, a condensate tray and a rectangular outer casing. The evaporator coil consists of three rows of tubes arranged in staggered formation, with each row stacked with 10 tubes. The tube interspacing, as measured from the tube centre lines, is 27 mm in the horizontal direction and 54.5 mm in the vertical direction. The system is divided into high- and lowpressure regions by the coil and baffles.

Mathematical formulation and numerical procedure The mathematical formulation for the system is based on the following assumptions: (1) (2) (3)

the physical domain is two-dimensional the flow is steady as the flow is predominantly of low Mach number (typically less than 0.2), the fluid can be treated as incompressible

--

!

In ~

6,7

595

tO20

--' o~oo°i~L.~o°O

- " , Out ~"

,,o

4 ,

I

# •

°o° o o 7

?

j oo

no heat transfer is considered between the air and its cooling medium, i.e. the flow is isothermal buoyancy effects and the humidity in the air are not considered.

(5)

Governing equations The equations governing the motion of air inside an airconditioning unit are the Reynolds-averaged NavierStokes equations and the mass continuity equation. In Cartesian tensor notation, these can be written as

°VJ-o

(1)

Oxj

OUi

10P

ts 0 [OUi

OUj'~

0

Oxi

.

(21

where i, j can take on values of 1 and 2 and a repeated subscript indicates a summation, Ui and u i are, respectively, the mean and instantaneous velocity component in the direction of coordinate x i, p is the local pressure, p is the density and # the dynamic viscosity, uiuj in the last term are collectively known as the Reynolds stresses and represent the effect of the turbulence on the mean flow field. They arise from time-averaging of nonlinear convection terms in the original momentum equations and have the appearance of stresses. Their presence results in the number of unknowns exceeding the number of equations available, i.e. the system is not closed. The computation of turbulent flow is essentially a search for a model which can express the Reynolds stresses in terms of known or calculable mean-flow quantities. A number of turbulence models have been developed over the years by various investigators. Typical different order models are the mixing length model, two-equation k ~ model and Reynolds stress model. The two-equation k - ~ model, which is of moderate complexity, has been extensively tested and proved to be adequate over a wide range of engineering applications. As a result, the twoequation k - e model 9 is adopted in the present work. The model employed the Boussinesq eddy viscosity hypothesis based on the analogy between molecular and turbulent motion. Hence pUiUj = Is, k OXi -i- OXi ) -- 7~!jk

(3)

where/z t is the eddy viscosity and k the turbulent kinetic energy, and 6ij is the Kronecker delta tensor: unity for i = j and zero for i ¢ j. The eddy viscosity hypothesis stems from the convenience associated with maintaining an approach, for turbulent flows, which is similar to that for laminar flows. However, in contrast to the laminar viscosity, the eddy viscosity is not a fluid property, but depends on the local state of turbulence. In this model, the eddy viscosity is related to the turbulent kinetic

3ZO

_ _ ~ l

x

Figure 4 Schematic of the modelled air-conditioning unit (mm) Figure 4 Schrma de I'installation de conditionnement d'air modrlisre

(unitr; mm )

(4)

373

Table 1

Turbulence model empirical constants

Tableau 1 Constances empiriques du modrle de turbulence C~

er#

o'~

C, I

C+2

0.09

t.0

1.3

1.44

1.92

Z.G. Xu et al.

374

energy k and to its dissipation rate e as k2

(4)

~t : PClz T

where k and e are computed through their o w n transport-diffusion equations

Ok

1 0 (I.tt Ok)..l_t.t t i/OUi OUjx~OUi

UJoxj-pOxj Uj Oxj

p O- xj + O- x )-ff xj

Us+ -

(5)

'

Y+ = £ J i l l

(9)

Once the velocity gradient and associated wall. shear stress are evaluated, the near-wall values for k and e can be calculated as k=

2 C,2-~

Us

$

p Oxj \a~ Oxj/

"~-Cd#t£ (OO i OOj~ OOi O lc \Oxj + OxiJ ~ -

logarithmic region, Us+ and y+ denote a dimensionless velocity and a dimensionless normal distance from the wall, respectively. They are related to the wall shear stress "/-wallby the following.

(6)

where C u, ak, or,, C,1 and C~2 are empirical constants whose values are given in Table 1.

3/4 3/2

[rwanl p@-~'

C~, k e -- - ny -

(10)

where n is the von Karman constant, taken as 0.435 in the computation.

Boundary conditions

Numerical procedure

The boundary conditions for the problem can be classified into three main categories.

The governing equations developed above constitute a system of coupled non-linear partial differential equations of elliptic nature, which must be solved numerically. The current work employs finite-volume methodology to discretise the differential equations and a non-staggered grid to descretise the physical domain. All the equations to be solved, apart from the continuity equation, can be expressed conveniently in the following general form

Inlet boundary. At an inlet, all variables, except pressure, are prescribed. As regards the pressure boundary, it can either be prescribed if known or expressed in terms of a mass flow rate. One set of values which is commonly adopted is listed as follows: k3/2 U = specified, V = 0, k = 0.003U 2, e - 0.03Di/2

(7)

Oxj pUj, -

-~xjJ = S~

(11)

where D i is the hydraulic diameter at the inlet.

Outlet boundary. At an outlet, the pressure boundary is applied. The velocities are determined in such a way that the mass conservation should be satisfied.

Wall boundary. In a turbulent flow, the steepest velocity gradients occur in the region closest to the wall. This renders the imposition of the usual non-slip condition impractical unless a very fine grid distribution is used near the wall. Moreover, it is physically meaningless to integrate the k-e model equations up to the wall because it is not valid in the vicinity of the wall. As the wall approaches, the flow regime changes to laminar flow. One solution is the use of the wall function method. With this method the innermost region of the flow is omitted and an artificial boundary condition is applied inside the wall at a point some distance from the real wall. It is assumed that the near-wall region is in local energy equilibrium. Thus the profile of the non-dimensional velocity component parallel to the wall is logarithmic and the turbulent shear stress is constant. Thus Us+=

y+ l l n ( E y +)

for y+ < y~f o r y + > y 0~

where P+ and S o are the relevant transport coefficient and source term for the variable 0, respectively. In the case of the momentum equations, the transport coefficient is the effective viscosity and the source term is the pressure gradient. The two terms on the LHS of the equation are the convection and diffusion term, respectively. In a non-staggered grid system, the physical domain is subdivided into a number of finite control volumes and all variables are co-located at the centroids of the control volume. The discretised form of the partial differential equations are derived by integrating the equation over

N

•

•

nl

•

(8)

s

Figm'e 5 Nodal arrangement Figure 5 Dispositionnodale

•

• •

P

where E is a roughness parameter and yo+ is the crossover point between the viscous sublayer and the *The expression 'synthetic boundary condition' is used.

0

Turbulent flow in air-conditioning unit the control volume surrounding each node. This process can be demonstrated by considering a typical node P shown in Figure 5. Integration of Equation (11) over the control volume gives

fvf- j (pUj - r, --xj)dV= fvS, dV

(12)

Following the divergence theorem and linearising the source term, the equation can be evaluated as

0~e

n

= S0,cq~p+S~,u

(13)

where the letters w, e, s and n refer to the west, east, south and north face of the control volume. Therefore what have to be evaluated are the convective and diffusive fluxes at and normal to the faces of the control volume. As variables are stored at the cell centres, the convective and diffusive fluxes at the faces must be derived by interpolation from those values at cell centres. A wide variety of interpolation methods (known as differencing schemes) may be selected. The diffusion term in Equation (13) is treated by a second-order central differencing scheme, which is based on a piecewise linear relationship. Considering the west face, the diffusive flux there can be approximated as (90

0p - 0w

(14)

where Aw is the area of the west face and 6hw is the distance between the two neighbouring nodes W and P. The treatment of convection terms has proved to be a difficult task because it determines the matrix coefficients of the resulting finite difference equations and thus the accuracy and stability of the solutions. The central differencing scheme, when applied to the convection term, can lead to negative matrix coefficients when the flow is convection dominated. As a result, the condition of diagonal dominance is not satisfied. This causes the diagonal dominance based algebraic equation solver to oscillate. Various schemes have been developed to deal with this problem. The most noteworthy are the upwind, the hybrid and the QUICK scheme. Details of these schemes are described below.

Upwind scheme In this scheme, it is assumed that the value of 0 prevailing at a cell face is equal to the value of ~bat the nearest grid node on the 'upwind' side of the face, i.e.

{(pUjA)wOw if Uw > 0 (pU/AO)w =

(pUjA)wOp

if Uw < 0

375

Hybrid scheme To improve the accuracy of the upwind scheme, a hybrid scheme was proposed by Spalding l°. This scheme uses central or upwind differencing according to the cell Peclet number Pe defined as

ee -- pU6h #

(16)

The scheme is identical to the central differencing scheme for the Peclet number range - 2 < Pe < 2. Outside this range, the convection dominates and the upwind scheme is used.

QUICK scheme A better approximation of 0w can be achieved by the use of the QUICK scheme (Quadratic Upstream Interpolations for the Convective Kinematics) proposed by Leonard II. It is based on a quadratic interpolation profile with upstream weighting and is third-order accurate. The scheme uses two upstream nodes and one downstream node. Thus

4'w=

{

l(q~w-{'-0p)--l(q~p-l-~ww -20w)

ifUw > 0

½(0w+0p)

ifUw<0

~(4'p+4'E--20p)

(17) The QUICK scheme has several attractive properties: no numerical diffusion, high order accuracy, algorithm simplicity and computational efficiency. For complex flow problems, it can provide accurate results without the need for the extensive grid refinement. However it may suffer from some non-physical overshoots and wiggles in the solution of highly convective flows. The convective and diffusive fluxes at the other faces of the control volume can be treated in a similar manner. No matter what differencing scheme is selected, the resulting algebraic equation for the variable 0 at the node P can be case into a general form as follows: Ap0p = Z Anb0nb + So, u nb

(18)

where As are the matrix coefficients representing the combined influence of convection and diffusion on the corresponding faces of the control volume and Ap = ~nbAnb0nb+S0, c. Y]nb denotes the summation over the four nearest neighbouring nodes ( n b = W, E, S, N). The source term should be linearised in such a way as to ensure the matrix coefficients are positive and diagonally dominated. The continuity equation can be discretised similarly by integration over the control volume. Then substituting the discretised momentum equations into the continuity equation, the pressure equation is then obtained as

(15)

The scheme is based on a piecewise constant profile and is only first-order accurate. The main drawback of the scheme is that it suffers from the so-called numerical diffusion, that is, the solution tends to be 'smoothed'. Consequently, for many complex flow problems, some main features of the flow, such as recirculation regions, may not be adequately predicted.

apPp = Z

AnbPnb + Sp, u

(19)

nb It is obvious that the computation of the convective coefficients requires a knowledge of the velocity components at the control volume faces. With the use of a nonstaggered grid system, however, the velocities are not stored at the faces but must be interpolated from those at the control volume centres. The naive prescription for

Z.G. Xu et al.

376

this would be to use weighted linear interpolation, which can lead to 'checkerboard' type oscillations. These are removed by applying the Rhie-Chow interpolation 12 algorithm in which the velocities at the control volume faces are obtained by use of a momentum interpolation. The algebraic equations are solved iteratively using a line-by-line alternating direction method based on the tridiagonal matrix algorithm (TDMA). The momentum equations are first solved with a guessed pressure field. Corrections to the guessed velocity fields and pressure field are made to satisfy the continuity equation according to the SIMPLEC algorithm 13. Then the velocities and pressure are updated and the turbulence equations are solved. The solution is assumed converged when the sum of the normalised mass residuals in the entire computational domain is less than a prescribed tolerance value. To ensure convergence, the variables are often under-relaxed in the iteration. Commonly used interactive algorithms are the Gauss-Seidel and SOR, for example. They have a simple structure and are cheap to use, if the iterations converge. But for many practical problems, convergence can be slow, with small relaxation factors, or non-existent. In this case it may be necessary to use a more efficient full field solver such as Stone's strongly implicit procedure (SIP) or preconditioned conjugate gradient (ICCG) algorithms. In the current work, SIP is used for the momentum equations whereas ICCG is used for the pressure equation. For turbulent equations, the line relaxation method is used. The principle of the main method of solution, described above, is embodied in a general-purpose CFD code CFDS-FLOW3D (version 2.4) 14, developed by the United Kingdom Atomic Energy Authority. The computation was carried out on a Sparc workstation from SUN Microsystems with a processing speed of 4.2 Mflops.

Grid generation A general non-orthogonal body fitted curvilinear grid system is adopted instead of the Cartesian grid system. It can be conformed to the curved surface of the geometry and thus no unnecessary errors are introduced into the boundary layer calculations by steps in the walls of the

geometry. A coordinate transformation is carried out to map an arbitrarily multiconnected region from a physical domain to a topologically rectangular mesh structure in a computational domain, as illustrated in Figure 6. The conservation equations are then transformed and solved in the computational domain. The methods used for generating the body-fitted grid can be divided into two distinct groups, namely algebraic and elliptical methods. The former is based on the transfinite interpolation technique and is easy to implement. With elliptical methods, on the other hand, a grid is produced by the solution of an elliptic partial differential equation such as the Laplace equation. Elliptical methods are more complex, but adequate for geometries with curved boundaries. They are also used to smooth the grid generated by the algebraic methods. The development of a suitable grid for a system with such a complex geometry as the one shown in Figure 4 has proved to be a major task. Fortunately, the process can be made much easier by a method called 'masking', which was previously developed by Gotham et aL e in modelling the complete air-conditioning unit. With this method, a complex geometry is broken down into a number of simple components or blocks. A separate grid is generated for each component with its own convenient coordinate system and origin. Either algebraic or elliptic methods can be employed, where appropriate. They are then patched together to form a complete grid. It is essential to ensure that the separate grids are correctly matched and the resulting grid for the combined system is smooth. Many minor adjustments are required before these conditions are satisfied. This idea is similar in principle to that of the multi-block approach. However, the resulting grid is still single block and therefore is superior to the multi-block approach in that the latter not only requires a complicated structure of the flow solver but significantly decreases the computational efficiency of solution of flow equations. As far as the geometry shown in Figure 4 is concerned, it is first broken down into an evaporator coil, a condensate tray and the outer casing. The evaporator coil is further divided into a rectangular box and a number of tubes. Figure 7 shows the resulting grid system generated with the masking method. The grid dimensions are 72 in the

Y

X Physical domain Figure 6 Coordinatetransformation Figure 6 Transformationcoordonn~e

I

i Computational domain

Turbulent flow in air-conditioning unit horizontal and 88 in the vertical direction. Each tube occupies 3 × 3 cells. Results and discussion Regarding the experimental investigation reported by Coney 3, the triple hot-wire anemometry results were found to be in agreement with the earlier single hot-wire measurements in both magnitude and distribution of velocity. In Figures 8 and 9 are shown mean velocity and RMS velocity plots for the X direction (i.e. normal to the coil) taken 110 mm downstream of the coil. In both, the effects of the condensate tray can be seen clearly. An

377

inference from these results, taken with those in the Y and Z directions, which is of importance to the modelling process, is that the turbulence is approximately isotropic. The fluid was assumed to enter the system with a uniform axial velocity profile of 2 . 6 m s -~ and was discharged to the atmosphere at the exit with a barometric pressure of 1.013bar. The density and kinematic viscosity of the air were taken to be at a room temperature of 22°C with values of 1.2 kg m -3 and 1.465 × 10 -5 m 2 s -] respectively. Based on the geometry and flow conditions at the inlet, the Reynolds number for the flow was estimated to be approximately 5.84 × 104. As a result the flow is fully turbulent.

Y

f" -

X

Figure 7

Computational grid

Figure 7

Grille de calcul

Triple ,~ire B3P4 ( 1 1 0 Inln dov~nstrcam of coil) [ l m e a n (X direction) m/s 700.

l/\ ',d

I ripIc

~--/"-~~l

v, irc B 3 P 4 ~110 m m do',~.il,,llc~ml ,,I ,.~,ill u i r e s ( 5 d i r c c l i ~ u l l hi/',

70()-

(~()(l'

:

:,00.

fl "/9

.~ 4 0 0 "5

.-~

"5

-~oo

400'

:-5 ~,oo

>-

C"

20o

_0 ()

..~

I0(}

~.~./ff)

'2

~

-) '(

Ii)()

I I I l[(I}ll}[l I [ [ I I I

2I)()

I I I I I [II

~0(I

40()

I I I ILI

5()()

I III

I L III

(~()(I 700

I ~ III

1 I 1

~()() 0 0 0

/ (]il'¢CllOll (llllllF

Figure 8 Contours of mean velocity (U) in the X direction Figure 8 Contours de la vitesse moyenne (U) dans la direction X

"' 0 . 0 3 iiiilCJLIll~illll;ll,lE~illlili~i~rPillYlEI;t I()0 2oi) 31)0 400 500

~l)()

7o()

8()0

9t)()

Z direction Imm) Figure 9 Contours of RMS velocity (URMs) in the X direction Figure 9 Contours de la vitesse R M S (Un~ls') dans la direction X

Z.G. Xu et al.

378

The convergence behaviour of the three differencing schemes is illustrated in Figure 10, where the whole-field mass source residual is plotted against iteration number. In the simulation the solution was declared converged when the reduction in the whole field mass residual (pressure equation) was below the 1 × 10-6 level. As shown in the diagram, the three schemes exhibit similar characteristics of convergence. Each encounters a very brief period of divergence at the initial stage of the interaction. As far as rate of convergence is concerned, the hybrid scheme is the most efficient, requiring only a total of 320 iterations. This compares with 365 for the upwind scheme and 378 for the QUICK scheme. It is observed that at about 120 iterations, the QUICK scheme experiences a blip with sudden increase in residual. However, this deviation is quickly rectified and the path to convergence is restored. Because of this, the scheme takes longest to converge. As pointed out in the previous section, the QUICK scheme has a drawback in that it can sometimes run into numerical instability. I0~. '~ o *

~ ~ ~ .

I0~

Upwind, Hybrid QUICK

-o 10" n,-

~10~i i

10 6 '

10'

•

i

•

•

[

1

100

200

300

400

500

Number of iterations

Figure 10 Historyof convergencefor mass Figure 10 Historique de la convergence pour

la masse

The phenomenon raises some questions with regard to the reliability of the scheme, especially in a complex three-dimensional flow. One well-known remedy would be to employ a small under-relaxation factor and it is important to ensure that the solution obtained is fully converged. Shown by the solid square symbols in Figure 11, for a vertical plane through the centre of the evaporator, is the experimental axial velocity (U) profile l lOmm downstream of the rear face of the evaporator. Corresponding mean velocity (F) data in the vertical direction are shown in Figure 12. As suggested by Figure 8, the mean velocity is reasonably constant over the top half of the evaporator. However, in the lower half, a marked peak is observed at a point approximately two-thirds of the way down the evaporator. At this point, a mean axial velocity approaching twice the average for the flow is observed. There is a corresponding rise in the mean vertical velocity at this position also, as shown in Figure 12. As stated previously, this effect is associated with the flow over the condensate tray, as is illustrated in Figure 13. Below this position of the peaks, the mean axial and vertical velocities fall sharply, the bottom three experimental points indicating a fairly uniform and moderately positive axial velocity. However the axial velocity in this region is likely to be negative, owing to recirculation of the flow beneath the condensate tray, as shown in Figure 13. (It should be noted that the hot-wire anemometry technique does not allow discrimination between positive and negative values.) It should be borne in mind that any such negative velocity will be subject to considerable experimental error, as the flow will be disrupted by the support prong of the triple hot-wire anemometry probe. Nevertheless, these values do give some indication of the extent of the recirculating region. Indicated in Figures 11 and 12 also, are the corresponding computational results for the mean velocity profiles using the three differencing schemes. In the upper half of the flow, all three sets of computed flow predict well the measured, approximately constant, mean axial

Measurement Plane

•

Experimental Upwind Hybrid x QUICK

~; ....

~--

Y

> x Figure 11

Comparison of numericaland experimentalU velocity

Figure 11

Comparaison de la vitesse U num~rique et exp~rimentale

~

.....

......

5.0 m / s

Turbulent flow in air-conditioning unit

379

Measurement Plane

...... • ~ =~

r.--....._

Experimental Upwind Hybdd QUICK

5.0 m/s

:> X Figure 12 Comparison of numerical and experimental V velocity Figure 12 Comparaison de la vitesse V nurn~rique et exp&imentale

and vertical velocities. Qualitatively, all three schemes show also the peaks in the velocity profiles, noted in the lower half of the flow. The predicted peak velocities are too low in magnitude and vertical position, suggesting an underestimation of the recirculating zone beneath the condensate tray. In contrast, at the top of the flow, all three schemes indicate second peaks in axial velocity somewhat greater than measured experimentally, with the peaks occurring too close to the upper wall. These suggest the presence of upper recirculation zones, both smaller and more vigorous than supported by experimental evidence. It is well known that the standard k-e model often leads to underprediction of recirculation zones. The model is based on the isotropic eddy viscosity hypothesis and is only valid at high Reynolds number in the flow. Its use in the boundary layer may lead to errors. In addition, a contributing factor may be the limited accuracy of the wall function method in the boundary layer region, as well as the matching ability of higher order discretisation procedures in this area. It may be more proper and accurate to employ a low Reynolds number k-e model 15. However, this model needs a much more refined grid along solid boundaries and was beyond the capabilities of the computer system used. Of the three differencing schemes tested, QUICK appears to reproduce best the experimental profile. This is not surprising since the scheme is third-order accurate and so should better predict circulation zones. By contrast, the upwind scheme is only a first-order scheme. It can be seen that compared with the experimental result, it smooths the velocity profiles and the peaks are not well predicted. This is the direct result of numerical diffusion. Its effects can be quite severe in complex recirculating flows with a strong velocity gradient. The strain rate is often underestimated. This problem can be solved by using a higher order upwind scheme. The third scheme, a hybrid scheme, gives a solution which falls between the two. Since the computational cost for the QUICK scheme was not

prohibitively higher than the others, it is recommended that higher order schemes should be used wherever possible. The vector plot of the flow field presented in Figure 13 is that computed using the QUICK scheme. The magnitude of each velocity can be measured with reference to the velocity (2.6ms -1) at the inlet. It can be seen that complex flow patterns have been developed within the unit. This is particularly the case for the flow in the neighbourhood of the coil and its drip tray. This effect had not been anticipated and revealing it pointed to the value of this approach, which offers a rapid method of analysing the air flow within such a unit. Obtaining these results led to a change in the design of the drip tray with beneficial effects on the air flow in general and on the performance of the evaporator. The flow downstream of the inlet frame and the coil are similar to that experienced over a backwardfacing step. Here, a region of negative velocity exhibits the recirculation region referred to previously and results from inertial forces of the fluid over the sharp edges and turbulence effects. When the flow approaches the coil, the majority of the fluid passes through the coil where it is straightened. However, the fluid near the wall is deflected from the coil side wall, and continues in a stream until it encounters the baffle surrounding the coil. Such an obstruction causes the flow to become stagnant. The computed pressure distribution within the unit is shown in Figure I4. It is obvious that there are regions of both high and low pressure. The two regions are divided by the evaporator coil and there is only a small variation in each region. The main gradient occurs across the coil and is found to be around 12Pa. Due to this small pressure difference, it is very difficult to measure the pressure within the unit. Thus no experimental data on pressure distribution was available for comparison. However, it was found that the computed pressure drop of 12Pa is in good agreement with the one calculated using the empirical correlations 16 for

380

Z.G. Xu et al.

(a)

Measurement Plane 0.11./

Y

> X

(b)

>

X

Figure 13 Velocity vector plot showing complex flow pattern: (a) overall view; (b) enlarged view showing flow in the drip tray Figure 13 Calcul du vecteur de vitesse montrant une configuration complexe de l'6coulement: (a) rue globale; (b) rue 6largie montrant l'~coulement

0.75- ! 12J

s

E,E 0.5o1=: ~, 0.25-

!

£o.5o

o.25 { i

i

[

I

0.25

I

I

i

I

'

0.50

1

I

I

0.75

1.00

i

f

i

I

~

1.25

'

i

,

I

1.50

Horizontal (m)

Figure 14 Pressure distribution within the unit (Pa) Figure 14 Distribution de la pression dans l'installation (unit~. pascal)

0.25

050

0.75

1.00

1,25

1.50

Horizontal (m) Figure 15 Distribution of turbulent kinetic energy (m2 s-2) Figure 15 Distribution de l'~nergie cin6tique turbulente (unit6." m2/s 2)

Turbulent flow in air-conditioning unit

381

and CFDS (UKAEA) for the use of their CFDSFLOW3D code.

0.75-

E 0.50

~D

~ _ _ . . ~

> 0.25

f

i

i

i

i

i

i

0.25

i

i

i

i

0.50

i

i

i

i

~

0.75

/

7.5 ~

"1~

i

f

f

.~ 7.5

i

i

1.00

References

3

.e'~'~

--

f

1

r

i

i

r

1.25

f

i

i

i

T

1.50

2

Horizontal (m) Figure 16 Distribution of turbulent energy dissipation rate (m 2 s 3) Figure 16 Distribution de la vitesse de dissipation de l'knergie turbulente (unit~: m2/s 3)

staggered tube configurations. The maximum pressure is located near the stagnation point. Here, there is a transfer of kinetic energy to pressure energy when the fluid becomes stagnant. It should be noted that the plotted pressure contours are only relative values. For absolute values the barometric pressure should be added. The turbulent kinetic energy and dissipation rate distributions are shown in Figures 15 and 16.

3

4

5

6

7

Conclusions It has been demonstrated that the turbulent air flow through tube banks of a packaged air-conditioning unit can be predicted with reasonable accuracy using the k-e turbulence model. A comparison of the mean velocity profiles from three differencing schemes with the experimental profile indicates that a higher order scheme such as QUICK provides the most accurate result, although it is computationally more expensive and care must be taken to ensure that the solution converges. The hybrid scheme has the fastest rate of convergence and its solutions are in good agreement with the experimental results.

Acknowledgements The work reported in this paper has been supported jointly by the Science and Engineering Research Council (SERC) and Department of Trade and Industry (DTI) of the United Kingdom. This support is gratefully acknowledged. The authors would like to thank Airedale International Air Conditioning Ltd (Leeds) for their assistance in obtaining the experimental data

8 9 10 11 12 13 14 15 16

Gotham, D. H. T., Xu, Z. G., Collins, M. W. Numerical modelling of three-dimensional turbulent flow in packaged air-conditioning unit Proc IMechE Conf Engineering Application qf Computational Fluid Dynamics London (1993) Gotham, D. H. T., Xu, Z. G., Collins, M. W. Modelling the air flow through the EHA5 air-conditioning system Research Memorandum RM 186, Thermo-Fluids Engineering Research Centre, City University, London (1993) Coney, J. R., Sheppard, C. G. W., Merdjani, S., Coates, W. The effect of air flow distribution on heat exchanger performance in refrigeration and air conditioning equipment Airdale International Air Conditioning/Leeds University Progress Report 4/1/ 1993 (1993) Patankar, S. V. Spalding, D. B. Computer analysis of the threedimensional flow and heat transfer in a steam generator Forsch lng-Wes (1978) 44 47 52 Shim, S. Y., Wan, P. T., Baxter, D. K., Hembroff, R. L. A study of submerged confined turbulent jets under suction and countercurrent flows Proc Numerical Methods in Laminar and Turbulent Flow Pineridge Press, Swansea, UK (1989) Janjua, S. 1., McLaughlin, D. K., Jackson, T. W., Lilley, D. G. Turbulence measurements in a confined jet using six-oriented hot wire probe technique AIAA-82-1262. A I A A / S A E / A S M E 18th Joint Propulsion Conf, Cleveland, OH (1982) Frota, M. N. Analysis of the uncertainties in velocity measurements and techniques for turbulence measurements in complex heated flows with multiple hot wires PhD thesis Stanford University, USA (1982) Merdjani, S. Study of jet mixing in an isothermal model of a gas turbine combustor dilution zone PhD thesis University of Leeds (1989) Jones, W. P., Launder, B. E. The prediction of laminarization with a two-equation model of turbulence lnt J Heat Mass TransJ'er (1972) 15 301-314 Spalding, D. B. A novel finite difference formulation for differential expressions involving both first and second derivatives Int J Numer Meth Engng (1972) 4 551 559 Leonard, B. P. A stable and accurate convective modelling procedure based on quadratic upstream interpolation Comput Meth Appl Mech Engng (1979) 19 59-98 Rhie, C. M. Chow, W. L. A numerical study of the turbulent flow past an isolated airfoil with trailing edge separation AIAA J (1983) 21 1525-1532 Van Doormaal, J. P., Raithby, C. D. Enhancement of the SIMPLE method for predicting incompressible fluid flows Numer Heat Tran.~[br (1984) 7 147-163 Burns, A. D., Jones, I. P., Kightley, J. R., Wilkes, N. S. HarwellFLOW3D. Release 2." User Manual Harwell Report AERE-R (Draft) (1990) Launder, B. E., Sharma, B. T. Application of the energy dissipation model of turbulence to the calculation of flow near a spinning disc Lett Heat Mass Transfer (1974) 1 131 - 138 Hewitt, G. F. Handbook of Heat Exchanger Design Hemisphere (1990)