- Email: [email protected]
Modified two-fluid model for air curtains in open vertical display cabinets
international journal of refrigeration 31 (2008) 472–482
available at www.sciencedirect.com
w w w . i i fi i r . o r g
journal homepage: www.elsevier.com/locate/ijrefrig
Modified two-fluid model for air curtains in open vertical display cabinets Ke-zhi Yua,b, Guo-liang Dinga,*, Tian-ji Chenb a
Institute of Refrigeration and Cryogenics, Shanghai Jiaotong University, Shanghai 200030, China Department of Refrigeration Engineering, Shanghai Fisheries University, Shanghai 200090, China
b
article info
abstract
Article history:
A modified two-fluid turbulence model is established to simulate the flow and heat transfer
Received 9 April 2007
characteristics of air curtains in an open vertical display cabinet. The air exiting the back
Received in revised form
panel (AEBP), the inner air curtain and the outer air curtain are taken as the first fluid
10 July 2007
and described with the standard k 3 turbulence model (KE). The air outside the display
Accepted 11 July 2007
cabinet is considered as the second fluid and calculated by the laminar model. Different
Published online 20 July 2007
from the existing two-fluid model, the mass transfer rate equation between the turbulent and the non-turbulent fluids is modified, and the face coefficient of honeycomb (the ratio of
Keywords:
the effective area to the total area of honeycomb) is fetched with the volume fraction of air
Refrigeration
curtain in order to conveniently depict the blow characteristics of the honeycomb. The
Commercial application
comparisons between the simulated temperatures and the experimental ones indicate
Display cabinet
that the modified two-fluid model (MTF) can give better agreement with the measurements
Open display cabinet
than KE and two original two-fluid models (TF1 and TF2). ª 2007 Elsevier Ltd and IIR. All rights reserved.
Air curtain Modelling Flow Heat transfer
Mode`le a` deux fluides modifie´ pour les rideaux d’air des meubles frigorifiques de vente verticaux Mots cle´s : Re´frige´ration ; Application commerciale ; Meuble de vente ; Meuble ouvert ; Rideau d’air ; Mode´lisation ; E´coulement ; Transfert de chaleur
1.
Introduction
Vertical open display cabinets are widely used in supermarkets. They can not only store various foods, but also attract customers’ attention so as to increase sales. In a vertical open display
cabinet, there are one or more forced-air curtains between the cabinet and the environment. However, these air curtains are easily disturbed by the ambient air, and then a more and less amount of ambient air is entrained into the display cabinet, which increases the energy consumption of the display cabinet.
* Corresponding author. Tel.: þ86 21 62932110; fax: þ86 21 34206814. E-mail address: [email protected] (G.-liang Ding). 0140-7007/$ – see front matter ª 2007 Elsevier Ltd and IIR. All rights reserved. doi:10.1016/j.ijrefrig.2007.07.006
international journal of refrigeration 31 (2008) 472–482
Nomenclature A Ae AEBP B c1, c2, cm cm1, cm2 cf, ct CFD cp D f F g G H I k KE l L m MTF p Pr
total area of honeycomb (m2) effective area of honeycomb (m2) air exiting the back panel cell size of honeycomb (m) constants in the k 3 turbulence model empirical constants for mass transfer rate (–) empirical constants for momentum and heat transfer, respectively (–) computational fluid dynamics specific heat capacity at constant pressure (J kg1 K1) phase diffusion coefficient in Eq. (1) (m2 s1) face coefficient of honeycomb (Ae/A) (–) inter-fluid friction force per unit volume (N m3 s1) gravitational acceleration (m s2) production rate of turbulent energy in the k 3 turbulence model (kg m1 s1) honeycomb height (m) inter-fluid sources in Eq. (1) turbulent kinetic energy (m2 s2) standard k 3 turbulence model Prandtl mixing length (m) honeycomb length (m) inter-fluid mass transfer rate per unit volume (kg m3 s1) modified two-fluid model based on Eq. (10) pressure (Pa) Prandtl number (–)
The honeycombs are generally installed at the top of the open display cabinet to form the air curtains. The honeycombs can decrease the turbulence intensities of air curtains, and reduce the mass and heat exchange between the air curtains and the ambient air, thus lessening the cooling load of the display cabinet. As the air curtains play a very important role in the cold preservation of display cabinets, various researches on air curtains have been made in recent years, including experimental studies and computational fluid dynamics (CFD) simulations. However, the mechanisms of flow, heat and mass transfer for air curtains are rather complicated. CFD simulation is a practical way to investigate the mechanisms of the air curtains and has been applied widely (Stribling et al., 1997; Ge and Tassou, 2001; Navaz et al., 2002; Wu et al., 2004; Cui and Wang, 2004; Fostera et al., 2005; D’Agaro et al., 2006; Cortella et al., 2001; Cortella, 2002). Among these studies, the most commonly used turbulent models are the k 3 model (Stribling et al., 1997; Ge and Tassou, 2001; Navaz et al., 2002; Wu et al., 2004; Cui and Wang, 2004; Fostera et al., 2005; D’Agaro et al., 2006) and the LES model (Cortella et al., 2001; Cortella, 2002). These models belong to single-fluid turbulence model, which takes the fluid in the whole computational domain as the turbulent fluid. But it is difficult for a single-fluid turbulence model to well describe the coexistence of highly turbulent flows inside the cabinet and largely laminar fluid outside the cabinet.
Q r S T TF1 TF2 U u v W
473
inter-fluid heat flux (J m3 s1) volume fraction (–) intra-fluid sources in Eq. (1) temperature ( C) two-fluid model based on Eq. (8) two-fluid model based on Eq. (9) velocity vector (m s1) velocity component in the x direction (m s1) velocity component in the y direction (m s1) honeycomb width (air curtain width) (m)
Greek symbols sK , s3 , sT , sr turbulent Prandtl number (–) m laminar viscosity (kg m1 s1) me effective viscosity (kg m1 s1) turbulence viscosity (kg m1 s1) mt r density (kg m3) d cell wall thickness (m) 3 dissipation rate of turbulent energy (m2 s3) f generic flow variable in Eq. (1) G within-phase diffusion coefficient in Eq. (1) (m2 s1) Subscripts 1 first fluid (turbulent air) 2 second fluid (non-turbulent air) amb ambient ex experimental value k fluid k i, j the spatial x, y direction in air curtain inlet
The two-fluid turbulence model provides an effective means of representing the coexistence of laminar flows and turbulent flows. This model considers the system as composed of two fluids, which coexist in time and space but possess different volume fractions. This model can reflect the exchange of mass, momentum and energy between the two fluids and is well suited to describe the whole computational domain including the high turbulence region inside the cabinet and relatively quiescent region outside the cabinet. A two-fluid model was adopted to simulate the air curtains for vertical display cabinets for the first time in Yu et al. (2007), in which the computational fluid was assumed to be composed of the turbulent fluid (air curtains) and the non-turbulent fluid (the ambient air outside the display cabinet). It showed that the two-fluid model can give better agreement with the measured values than the k 3 model. However, the different mass transfer rate relations between the two fluids were adopted by different researchers in different cases (Malin and Spalding, 1984; Shen et al., 2003; Markatos et al., 1986; Ilegbusi, 1994; Markatos and Kotsifaki, 1994; Fan, 1988; Sheng and Jonsson, 2000; Zhou, 1991). The mass transfer rate equation employed by Yu et al. (2007) assumed that only the turbulent fluid entrains the non-turbulent fluid; nevertheless the turbulent fluid cannot be transferred into the non-turbulent fluid. However, the experiments showed that the cold turbulent air may spill from the display cabinet into the ambient air
474
international journal of refrigeration 31 (2008) 472–482
back of the cabinet. A small fraction of cooled air is fed into the cabinet through the back panel of the cabinet, while more cold air is blown through two honeycombs to form two air curtains. The two air curtains include an inner air curtain with a high velocity on the inside (cabinet side) and an outer air curtain with a low velocity on the outside (store side). The inner air curtain and the outer air curtain, together with the entrained ambient air, are drawn back into the cabinet through the respective air curtain outlets at the bottom of the display cabinet. The inner air curtain is re-cooled by the evaporator while the outer air curtain bypasses the evaporator in the air flow tunnel, which makes the temperature of the inner air curtain lower than that of the outer air curtain. When the two-fluid model is applied to the vertical display cabinet, the following assumptions are made:
(Stribling et al., 1997; Navaz et al., 2002; D’Agaro et al., 2006; Cortella et al., 2001; Cortella, 2002). Therefore the entrainment rate between the two fluids should be modified according to the mass transfer mechanics of display cabinets. Furthermore, the difference between effective area (pass area) and total area of the air curtain honeycomb was not taken into account in Yu et al. (2007), which led to overestimating the momentum and the energy of air curtain discharge. Therefore the velocity, temperature and fraction volume of air curtain should be revised to depict air curtain honeycomb properly. The purpose of this paper is to provide a modified two-fluid model of air curtains, which can better predict the air flow field inside vertical display cabinets. In the modified model, the new mass transfer rate equation between the turbulent fluid and the non-turbulent fluid is proposed, and the face coefficients of honeycombs are obtained with the volume fractions of air curtains in order to depict the blow characteristics of the honeycomb properly.
(1) The system is assumed to consist of two fluids. The air exiting the back panel (AEBP), the inner and the outer air curtains are defined as the first fluid (turbulent fluid). The ambient air is defined as the second fluid (non-turbulent fluid). (2) A steady two-dimensional model is employed regardless of the display cabinet length. (3) The standard k 3 two equation model (KE) is adopted to simulate the turbulent fluid.
2. The two-fluid model for display cabinets (Yu et al., 2007) In a vertical open display cabinet (shown in Fig. 1), the air is forced by propeller fans to flow through an evaporator at the
Honeycomb y
Inner air curtain flow
x
Back panel
Outer air curtain flow
Outer air curtain inlet
Inner air curtain inlet
S1 20 14=280mm
20 14=280mm
S2
Evaporator
20 14=280mm
S3
T-type thermocouples
S4 Air flow tunnel
20 14=280mm
Inner air curtain outlet
Outer air curtain outlet
50 3=150mm
O1 O2 O3 O4
Fan
20 10=200mm
Fig. 1 – A vertical open display cabinet and the location of T-type thermocouples mounted.
475
international journal of refrigeration 31 (2008) 472–482
(4) Radiation heat transfer is not taken into consideration. (5) The Boussinesq hypothesis is adopted due to the density varying with aerostatic pressure.
Q¼
ct cp r1 r1 r2 ðT1 T2 Þju1 u2 j l
(6)
in which l is the Prandtl mixing length, expressed as: The two-fluid model for the display cabinet can be established and the governing transport equations can be written with reference to a generic flow variable fk (for fluid k) as: v v vf v vrk rk Gfk k þ fk Dfk þ Sfk þ Ifk rk rk Ukj fk ¼ vxj vxj vxj vxj vxj
(1)
where subscript k denotes fluid k, which can be set to 1 or 2 and subscript j refers to spatial coordinates. G is the diffusion coefficient within fluid 1 or fluid 2 and D is the phase diffusion coefficient between fluid 1 and fluid 2. S is the source term within fluid 1 or fluid 2 and I is the source term between fluid 1 and fluid 2. The continuity, momentum and energy equation can be obtained by setting the generic flow variable fk to unity, velocity and temperature, respectively. The governing equations for turbulent and non-turbulent fluid are shown in Tables 1 and 2, where: mt ¼
cm rk2 3
me ¼ m þ m t
G ¼ mt
vU1i vU1i vU1j þ vxj vxj vxi
3=2 c3=4 m k 3
(7)
Two popular equations of mass transfer rate between the two fluids are depicted as follows (Yu et al., 2007; Malin and Spalding, 1984; Shen et al., 2003; Markatos et al., 1986; Ilegbusi, 1994; Markatos and Kotsifaki, 1994): m¼
cm1 r1 r1 r2 ju1 u2 j l
(8)
and (Sheng and Jonsson, 2000),
m¼
cm2 r1 r1 r2 ðr2 0:5Þju1 u2 j l
(9)
(2)
The constants cf, ct, cm1, and cm2, employed in Eqs. (5), (6), (8) and (9), are 0.05, 0.05, 0.1 and 10.0, respectively.
(3)
3.
Modifications of the two-fluid model
3.1.
The modification of mass transfer rate equation
(4)
The empirical constants in the two-fluid model are listed in Table 3. The momentum and energy transfer rates are generally expressed as: Fi ¼
l¼
cf r1 r1 r2 ðU1i U2i Þju1 u2 j l
(5)
The mass transfer rate equation between the two fluids plays a very important role in the two-fluid model. Eqs. (8) and (9) are the most widely used relations. In Eq. (8) m is always positive, which means that only fluid 2 (non-turbulent fluid) can be entrained by fluid 1 (turbulent fluid). However, the mass transfer rate m may be negative according to Eq. (9). An additional factor r2 0.5 included in Eq. (9) allows for the equally entrainment between the turbulent fluid
Table 1 – The governing equations of turbulent fluid Equation
fk
Gfk
Dfk
Sfk
Ifk
Continuity
1
0
me sr
Momentum
U1i
me
me sr
Energy
T1
mt m þ sT Pr
me sr
0
Kinetic energy
k
mt þm sk
me sr
r1 ðG r1 3Þ
0
Dissipation rate
3
mt þm s3
me sr
r1 3ðc1 G c2 r1 3Þ k
0
0
m
vU1j v vp þ r1 ðr1 r2 Þgi me r1 r1 vxj vxi vxi
U2i m Fi
T2 m
Q cp
476
international journal of refrigeration 31 (2008) 472–482
Table 2 – The governing equations of non-turbulent fluid Equation
fk
Gfk
Dfk
Continuity
1
0
m sr
Momentum
U2i
m
m sr
Energy
T2
m Pr
m sr
(fluid 1) and the non-turbulent fluid (fluid 2). From the viewpoint of physics, the entrainment rate of a non-turbulent fluid by a turbulent fluid is much more than that of the turbulent fluid by the non-turbulent fluid. Therefore both Eqs. (8) and (9) have some disadvantages. A new mass transfer rate equation should be developed. Based on Eqs. (8) and (9), the modification of mass transfer rate is established by using the weighted average of the two equations. The weighted factors can be assigned to the volume fractions of the two fluids (r1 and r2). The modified mass transfer rate equation can be expressed as: cm1 r1 r1 r2 ju1 u2 j cm2 r1 r1 r2 ðr2 0:5Þju1 u2 j $r1 þ $r2 l l 2 2 cm1 r1 r1 r2 ju1 u2 j cm2 r1 r1 r2 ðr2 0:5Þju1 u2 j þ ¼ l l
m¼
(10)
Fig. 2 shows the relationship of the dimensionless number m=ðr1 ju1 u2 j=lÞ vs. the volume fraction of fluid 1 (r1) based on Eqs. (8)–(10). From the figure, it can be seen that Eq. (8) neglects the mass transfer rate m from fluid 1 to fluid 2, whereas Eq. (9) over predicts the mass transfer rate. It seemed that Eq. (10) can give reasonable result. Furthermore, Fig. 2 shows the mass transfer rate m from fluid 1 to fluid 2 less than m from fluid 2 to fluid 1 based on Eq. (10). This is consistent with the physics of the problem. Thus Eq. (10) can reflect the mass transfer rate between the two fluids perfectly. Besides the standard k 3 model (KE), the two-fluid models based on Eq. (8) (two-fluid model 1, TF1), Eq. (9) (two-fluid model 2, TF2), and Eq. (10) (modified two-fluid model, MTF) are adopted to simulate the open vertical display cabinet so that the four models can be compared with each other.
3.2.
Sfk
Ifk
0
m
vU2j v vp mr2 r2 vxj vxi vxi
U2i m þ Fi
0
T2 m þ
Ae < A
Q cp
(11)
If the experimental velocity and temperature of air curtain are uex and Tex, respectively, the inlet volume fraction of air curtains is equal to unity according to Yu et al. (2007). Then the theoretical inlet momentum and energy can be expressed as: r1 u1 Ar1 $u1 ¼ r1 Au21 ¼ r1 Au2ex
(12)
and r1 cp u1 Ar1 $T1 ¼ r1 cp Au1 T1 ¼ r1 cp Auex Tex
(13)
However, the actual inlet momentum and inlet energy are r1 Ae u2ex and r1 cp Ae uex Tex , respectively. The theoretical inlet momentum and energy are larger than the actual inlet momentum and energy, respectively, according to Eq. (11). The face coefficient of the honeycomb f can be defined by the ratio of effective area Ae to the total area A, namely: f ¼ Ae/A. If the volume fraction of the turbulent fluid (fluid 1) is assigned as the face coefficient of the honeycomb f, then the theoretical inlet momentum and energy can be expressed as: r1 u1 Ar1 $u1 ¼ r1 r1 Au21 ¼ r1 fAu21 ¼ r1 Ae u2ex
(14)
and r1 cp u1 Ar1 $T1 ¼ r1 cp r1 Au1 T1 ¼ r1 cp fAu1 T1 ¼ r1 cp Ae uex Tex
(15)
Therefore, the theoretical inlet momentum and energy satisfy the actual inlet momentum and energy, respectively.
The modification of inlet parameters of air curtain
For a honeycomb shown in Fig. 3, the total area A is larger than the effective area Ae, namely:
Table 3 – Values of the constants employed in the two-fluid model cm 0.09
c1
c2
sk
s3
sr
sT
1.44
1.92
1.0
1.3
0.95
1.0
Fig. 2 – The relationship of the dimensionless number m=ðr1 ju1 Lu2 j=lÞ vs. the volume fraction of fluid 1 (r1) based on Eqs. (8)–(10).
477
The computational domain is constructed, in which the evaporator, the insulation layer, the propeller fan and the product
Table 4 – The modifications of the inlet parameters of the air curtains for display cabinet Fluid 1
Before modification After modification
Fluid 2
u1
v1
T1
r1
u2
v2
T2
r2
uex uex
0 0
Tex Tamb
1 f
uex 0
0 0
Tex Tamb
0 1f
0.1 0.1 0.1 0.1 0.9 0.9 0.9 0.9 25 25 25 25 6.1 7.2 7.1 6.2 0 0 0 0 0 0 0 0 0 0 0 0 0.25 0.25 0.25 0.28 0.1 0.1 0.1 0.1 0.9 0.9 0.9 0.9 25 25 25 25 5.3 4.4 6.4 5.1 0 0 0 0 0 0 0 0 0 0 0 0 0.52 0.58 0.3 0.52 6.1 7.2 7.1 6.2
v2 (m s1) v2 (m s1) v1 (m s1) u2 (m s1) u1 (m s1) Temperature ( C) Velocity (m s1)
The computational domain and grid
Temperature ( C)
4.2.
Velocity (m s1)
Four cases are simulated for a display cabinet. The specifications of the cases are presented in Table 5. Case 1 is the basic working condition for a display cabinet. Case 2 and Case 3 examine the effect of the inner air curtain velocity on the performance of the display cabinet. Case 4 examines the effect of the outer air curtain velocity on the performance of the display cabinet. The temperatures of the air curtains are assigned according to the measurements in order to be compared with the experimental data.
Inner air curtain
The computational case
Outer air curtain
4.1.
Inner air curtain
Case study
Table 5 – Case specifications and boundaries of air curtains inlet
4.
T1 ( C)
In addition, the inlet velocity and the inlet temperature of fluid 2 (the non-turbulent fluid) are assigned to zero and the ambient temperature, respectively. The modified parameters are listed in Table 4.
0.25 0.25 0.25 0.28
(17)
Case specifications
r2 ¼ 1 f
r2 (–)
The volume fraction of the non-turbulent fluid (fluid 2) can be assigned by Eq. (17).
r1 (–)
(16)
Case
r1 þ r2 ¼ 1
T2 ( C)
The constraint of the volume fraction is:
u1 (m s1)
The boundaries of air curtains inlet
Fig. 3 – The honeycomb of air curtain for display cabinet.
5.3 4.4 6.4 5.1
Honeycomb
v1 (m s1)
air curtain inlet
Cell structure
u2 (m s1)
Outer air curtain
H
B
T1 ( C)
W
Ae
0.52 0.58 0.3 0.52
A
1 2 3 4
T2 ( C)
L
r1 (–)
r2 (–)
international journal of refrigeration 31 (2008) 472–482
478
international journal of refrigeration 31 (2008) 472–482
(2) The air exiting the back panel (AEBP) The velocity of AEBP is about 0.1 m s1 by the experiments. The temperature of AEBP is equal to the temperature of the inner air curtain. In addition, the turbulent kinetic energy and the turbulent dissipate rate of AEBP can be the same as those assigned to the air curtains. (3) Ambient air The influence of the air curtains on the extended boundaries can be ignored. The temperature and the volume fraction of the ambient air are assigned to 25 C and 1.0, respectively. (4) Air curtain outlet Flow rate condition is imposed at the air curtain outlet. The outer air curtain outlet intakes the flow rate of the outer air curtain, whereas the inner air curtain outlet intakes the flow rate sum of the inner air curtain and the air exiting the back panel (AEBP). (5) Wall The non-slip boundary condition (zero velocity) is adopted at the solid surfaces. The wall function is employed for the turbulent fluid.
shelves are removed in order to simplify the simulation. As the open boundary conditions of the display cabinet are uncertain, the computational domain is extended until the effect of the display cabinet opening is negligible there. Then the computational domain becomes a 4.0 m 4.0 m irregular region, as shown in Fig. 4. The asymmetric grid consists of 61 63 cells, which is generated to mesh the computational domain, as illustrated in Fig. 4. The denser grid is adopted in the air curtain inlet and outlet regions where the air flow is complex. The grid independence is performed by using the progressively finer grid until the average temperature difference at the air curtain outlet is less than 0.1 C.
Boundary conditions
(1) Air curtain inlet The honeycomb length, height, width, cell size and cell wall thickness (shown in Fig. 3) are 2500 mm, 65 mm, 76 mm, 4.0 mm, and 0.2 mm, respectively. It can be calculated that the face coefficient of honeycomb is 0.9. The velocities, the temperatures and the volume fractions of the inner and the outer air curtains are shown in Table 5 according to the modified two-fluid model. The turbulent kinetic energy of the air curtains can be set to be 1% of the average kinetic energy (Cui and Wang, 2004). The turbulent dissipate rate is determined by Eq. (18) (Cui and Wang, 2004). 500k2 cm uin L
Simulation tool
The modified two-fluid model is embodied in the FORTRAN program developed by the present authors. The calculation needed 10,000 iterations in order to obtain a stable and converged solution. This required about 5 h on a personal computer with 2.09 GHz CPU.
(18)
5.
where cm is a constant which is shown in Table 3, uin is the average inlet velocity and L is the characteristic dimension of the air curtain (here it is the width of the honeycomb). y
Results and discussion
T-type thermocouples were mounted around the display cabinet to validate the simulation models (KE, TF1, TF2 and MTF).
The air curtain inlet
x
Unit:mm
76 76
200
50
500
The ambient space
210
1435
40
210
5 5
40
900
5
225
4000
90
5
The air exiting the back panel
450
330
70
900
The air curtain outlet
Fig. 4 – The computational domain and grid.
4000
The display cabinet
40
5
210
40
Product shelf
1000
3¼
4.4.
1435
4.3.
international journal of refrigeration 31 (2008) 472–482
479
Fig. 5 – Comparison of the predicted temperatures by KE, TF1, TF2 and MTF with the experimental temperatures outside shelves (a) S1 (b) S2 (c) S3, and (d) S4.
Fig. 1 shows that the measuring points can be divided into two groups. A group of measuring points (S1–S4) was located outside the product shelves where the temperature gradient is expected large, whereas another group of measuring points (O1–O4) was located near the air curtain outlet of the display cabinet where the cold air may overspill the display cabinet. All the measurements were recorded through the temperature acquisition logger at intervals of 10 s with the average value being calculated to be compared with the predictions.
5.1. Comparison with the measurements outside the shelves Fig. 5(a)–(d) shows the predicted temperatures by four models (KE, TF1, TF2 and TFM) and the experimental temperatures outside the four shelves (the first, the second, the third and the fourth shelf) for Case 1, respectively. The comparisons for Case 2, Case 3 and Case 4 are similar. It can be concluded from the figures as follows.
(1) The variation tendency of the simulated temperatures by four models is in agreement with that of the experimental data. However, the difference between the measurements and the predictions by four models is relatively large for the fourth shelf (S4). (2) Since KE assumes that the air in the whole computational domain is a strong turbulent fluid, the simulated temperatures outside the cabinet are very close to the ambient temperature 25 C, which are higher than the measurements (21–23 C). On the contrary, at the shelf the predicted temperatures (subzero) are lower than the measurements (about 3 C) when used by KE owing to overestimating the inlet momentum and energy of the cold air curtain. (3) The mass transfer rate m from fluid 1 to fluid 2 is not taken into account in TF1 (see Fig. 2). The flow field simulated by TF1 is approximate to the turbulence, thus the predicted temperatures by TF1 are very close to that by KE outside the cabinet. The differences between the predictions by TF1 and those by KE are less than 0.2 C for the most measuring points. However, TF1 predicts higher
480
international journal of refrigeration 31 (2008) 472–482
temperatures than KE when the measuring points locate very near the shelves. The maximum difference between the predictions by KE and TF1 can be found at the product shelves. (4) TF2, which over predicts the mass transfer rate m from fluid 1 to fluid 2 (see Fig. 2), predicts lower temperatures than TF1 and KE. It can be seen that the simulated temperatures by TF2 are lower than those by TF1 and the differences are 0.9–5.5 C. (5) When the measuring points locate at the first shelf (S1), the maximum differences between the experimental temperatures and the predicted values from TF1 and TF2 are 3.3 C (S1) and 4.6 C (S1), respectively. Thus the predicted temperatures by TF1 are in better agreement with the experimental values than those by TF2 when the thermocouples are close to the shelves. (6) However, when the measuring points are far away the shelf, the average differences between the measured temperatures and the predictions from TF1 and TF2 are 1.3 C and 1.1 C, respectively. It can be inferred that TF2 shows more agreement with the measurements when the measuring points locate far away from the shelves.
(7) The modification effect of TFM is significant. It is compatible with both TF1 and TF2 (see Fig. 2). When the measuring points locate at the shelf, the maximum difference of MTF (2.0 C for S3) is less than that of TF1 (3.3 C for S1) and that of TF2 (4.6 C for S1). Moreover, the average difference of MTF (1.0 C) is less than that of TF1 (1.3 C) and that of TF2 (1.1 C) when the thermocouple is far away the shelf. Therefore, MTF shows good agreement with the measurements wherever the experimental points locate.
5.2. Comparison with the measurements outside the display cabinet The predicted temperatures by the four models (KE, TF1, TF2 and TFM) and the experimental temperatures outside the display cabinet for Case 2 are shown in Fig. 6(a)–(d), respectively. The comparisons for Case 1, Case 3 and Case 4 are similar. It can be inferred from the figures: (1) KE, which takes the air outside the display cabinet as the turbulent fluid, predicts higher temperatures (the maximum is 16.8 C) than the measurements which range
Fig. 6 – Comparison of the predicted temperatures by KE, TF1, TF2 and MTF with the experimental temperatures outside the display cabinet (a) O1, (b) O2, (c) O3 and (d) O4.
481
international journal of refrigeration 31 (2008) 472–482
Table 6 – The average and the maximum temperature differences between the measurements and the predictions using the kL3 model and three (modified) two-fluid models ( C) Case
Case 1
Case 2
Case 3
Case 4
(2)
(3)
(4)
(5)
Measuring location
Average differences KE
Maximum differences
TF1
TF2
MTF
KE
TF1
TF2
MTF
Outside the shelf
S1 S2 S3 S4
1.2 1.8 2.0 3.6
1.2 1.9 1.8 2.9
1.1 0.8 1.2 0.3
0.2 0.3 0.3 0.3
3.9 4.0 5.0 5.8
3.4 3.8 4.8 5.6
4.6 4.0 5.1 2.5
1.5 2.1 3.6 3.1
Outside the cabinet
O1 O2 O3 O4
3.2 2.1 0.9 1.2
2.1 3.3 1.5 1.8
1.5 0.8 0.8 0.7
1.8 1.0 0.6 0.5
4.7 2.9 1.8 2.1
3.7 4.1 2.6 2.7
3.1 1.5 1.8 1.8
3.3 1.8 1.6 1.4
Outside the shelf
S1 S2 S3 S4
1.0 1.1 1.7 2.9
0.9 1.1 1.5 1.9
1.3 1.5 0.8 1.8
0.1 0.5 0.0 0.7
4.0 3.6 4.7 4.9
3.3 3.3 4.3 4.5
4.3 4.9 4.2 4.5
1.3 3.7 2.2 4.4
Outside the cabinet
O1 O2 O3 O4
2.2 2.0 0.4 1.0
1.4 3.7 3.0 1.6
1.0 1.2 0.3 2.2
0.6 0.7 0.7 1.4
4.4 3.0 1.4 2.4
3.7 4.5 4.1 4.0
3.2 2.2 1.6 3.7
2.8 1.6 1.8 2.9
Outside the shelf
S1 S2 S3 S4
1.9 2.3 2.6 1.0
1.8 2.6 2.1 0.1
0.2 0.1 0.1 2.4
1.1 0.5 0.5 1.7
3.3 4.2 4.6 4.3
3.4 4.7 3.9 4.5
2.5 2.2 2.9 5.3
3.2 1.9 1.8 4.8
Outside the cabinet
O1 O2 O3 O4
2.0 1.3 1.8 2.4
2.9 1.5 0.8 0.2
1.2 2.0 1.2 1.3
0.7 1.5 1.0 0.3
3.8 3.8 3.6 3.9
4.7 4.5 2.5 1.6
2.9 4.5 2.9 2.9
2.5 4.1 2.6 1.7
Outside the shelf
S1 S2 S3 S4
2.2 1.7 2.2 1.8
2.2 1.7 2.0 0.8
0.4 0.7 0.8 2.2
0.8 0.3 0.0 1.8
5.5 3.7 5.0 2.9
5.4 3.4 4.7 2.7
4.9 3.4 4.1 4.3
3.2 1.4 2.4 4.5
Outside the cabinet
O1 O2 O3 O4
1.3 1.0 0.9 0.2
2.1 1.2 1.9 0.8
1.0 0.7 0.3 0.6
0.7 0.6 0.2 0.0
2.6 3.2 3.2 1.8
3.5 3.7 4.2 3.1
2.4 3.2 2.9 2.6
2.6 3.2 2.9 1.8
from 11.7 C to 15.7 C for the first row of measuring points (O1, see Fig. 6(a)). KE shows similarity for O2, O3 and O4. The flow field simulated by TF1 is approximate to the turbulence, therefore TF1 over predicts the temperatures outside the display cabinet. The maximum differences between the predictions by TF1 and the measurements can be 4.5 C (O2). It is abnormal that TF1 gives abnormally lower temperatures (14.6–16.0 C) at the first row of measuring points (O1, see Fig. 6(a)). This phenomenon may be caused by the overspill of cold air from the cabinet into the ambient. Although KE and TF1 simulate the higher temperatures outside the cabinet, KE gives lower predictions than TF1 except the first row of measuring points (O1) owing to over predicting the inlet momentum and energy of the cold air curtain. The differences between the predictions by TF1 and those by KE are in the range of 0.8 C (O1) to 2.7 C (O3 and O4). The temperatures outside the cabinet predicted by TF2 and MTF are lower than those by KE and TF1 except the first row of measuring points (O1). The maximum differences
between the measurements and the predictions by KE, TF1, TF2 and MTF are 4.4 C, 4.5 C, 3.2 C, and 2.8 C, respectively. It is apparent that MTF and TF2 are in agreement with the measurements and MTF shows better agreement with the experimental temperatures than TF2.
5.3. The average and the maximum temperature differences by the four models Table 6 summarizes the average and the maximum temperature differences for the four cases simulated by the four models. From the table it can be found that the average and the maximum temperature differences are the smallest for most cases when MTF is employed (which are bolded in Table 6), which demonstrated that MTF is a valuable simulation model for display cabinet. It can also be found in Table 6 that TF2 shows better agreement with the measurements than MTF for Case 3 when measuring points were located outside the product shelves (S1–S3). The velocities of the air curtains are so small (0.30 m s1 and 0.25 m s1) for Case 3 that the turbulent fluid
482
international journal of refrigeration 31 (2008) 472–482
(fluid 1) and the non-turbulent fluid (fluid 2) can be equally entrained each other. TF2 shows greater advantage than MTF. Therefore MTF is not the appropriate numerical model for display cabinet when the velocities of the air curtains are quite small. However, there are maximum temperature differences of 4.8 C in the MTF simulation from Table 6. The reasons mainly come from two aspects. Firstly, the two-dimensional model assumption is made in applying the modified two-fluid model to the simulation of display cabinet. The two-dimensional model is not able to simulate the influence of the longitudinal direction on the display cabinet. Secondly, the measuring temperatures outside the display cabinet may have some extra uncertainty, which is caused by the interaction of two airstreams including the cold air spill over the display cabinet and the hot air entrained by the fan.
6.
Conclusions
A modified two-fluid model has been developed to simulate air flow and heat transfer in an open vertical display cabinet. The air exiting the back panel (AEBP), the inner air curtain and the outer air curtain are defined as fluid 1 (turbulent fluid), and the ambient air outside display cabinet is defined as fluid 2 (non-turbulent fluid). The volume fraction, velocity vector and temperature field are simulated by solving the twodimension steady transport equations of each fluid. The main conclusions of the numerical investigation are as follows: (1) In the modified two-fluid model, a modified mass transfer rate equation between the two fluids is established. It not only describes the entrainment of the non-turbulent fluid by the turbulent fluid, but also reveals the mass transfer rate from the turbulent fluid to the non-turbulent fluid. Thus the modified mass transfer rate equation can reflect the proper mass transfer rate between the air curtains and the ambient air in an open vertical display cabinet. (2) In the modified two-fluid model, the volume fractions of the first fluid (air curtains) are assigned as the face coefficients of honeycombs to depict the honeycombs, which can provide proper discharge momentum and energy of air curtains. (3) The comparisons with the experimental temperatures indicate that the agreement with the measurements can be improved by using the modified two-fluid model instead of the standard k 3 model (KE) and original two-fluid model (TF1 and TF2) in simulation of display cabinets.
references
Cui, J., Wang, S., 2004. Application of CFD in evaluation and energy-efficient design of air curtains for horizontal
refrigerated display cases. International Journal of Thermal Sciences. 43, 993–1002. Cortella, G., Manzan, M., Comini, G., 2001. CFD simulation of refrigerated display cabinets. International Journal of Refrigeration 24, 250–260. Cortella, G., 2002. CFD-aided retail cabinets design. Computers and Electronics in Agriculture 34, 43–66. D’Agaro, P., Cortella, G., Croce, G., 2006. Two- and three dimensional CFD applied to vertical display cabinets simulation. International Journal of Refrigeration 29, 178–190. Fostera, A.M., Madgeb, M., Evansa, J.A., 2005. The use of CFD to improve the performance of a chilled multi-deck retail display cabinet. International Journal of Refrigeration 28, 698–705. Fan, W.-C., 1988. A two-fluid model of turbulence and its modifications. Science in China Series A: Mathematics, Physics, Astronomy and Technological Sciences 31 (1), 79–86 (in Chinese). Ge, Y.T., Tassou, S.A., 2001. Simulation of the performance of single jet air curtains for vertical refrigerated display cabinets. Applied Thermal Engineering 21, 201–219. Ilegbusi, O.J., 1994. Application of the two-fluid model of turbulence to tundish problems. ISIJ International 34 (9), 732–738. Malin, M.R., Spalding, D.B., 1984. A two-fluid model of turbulence and its application to heated plane jets and wakes. Physicochemical Hydrodynamics 5 (5/6), 339–362. Markatos, N.C., Pericleous, K.A., Cox, G., 1986. A novel approach to the field modeling of fire. Physicochemical Hydrodynamics 7 (2/3), 125–143. Markatos, N.C., Kotsifaki, C.A., 1994. One-dimensional, two-fluid modeling of turbulent premixed flames. Applied Mathematical Modelling 18, 646–657. Navaz, H.K., Faramarzi, R., Gharib, M., Dabiri, D., Modarress, D., 2002. The application of advanced methods in analyzing the performance of the air curtain in a refrigerated display cabinet. Journal of Fluids Engineering: Transactions of the ASME 124 (3), 756–764. Stribling, D., Tassou, S.A., Marriott, D., 1997. A two-dimensional CFD model of a refrigerated display case. ASHRAE Transactions 103 (1), 88–94. Shen, Y.M., Ng, C.-O., Chwang, A.T., 2003. A two-fluid model of turbulent two-phase flow for simulating turbulent stratified flows. Ocean Engineering 30, 153–161. Sheng, D.Y., Jonsson, L., 2000. Two-fluid simulation on the mixed convection flow pattern in a nonisothermal water model of continuous casting tundish. Metallurgical and Materials Transactions B: Process Metallurgy and Materials Processing Science 31 (4), 867–875. Wu, Y., Xie, G., Chen, Z., Niu, L., Sun, D.W., 2004. An investigation on flowing patterns of the airflow and its characteristics of heat and mass transfer in an island open display cabinet with goods. Applied Thermal Engineering 24, 1945–1957. Yu, Ke-zhi, Ding, Guo-liang, Chen, Tian-ji, 2007. Simulation of air curtains for vertical display cases with a two-fluid model. Applied Thermal Engineering 27, 2583–2591. Zhou, L.X., 1991. The Numerical Simulation of Turbulent TwoPhase Flow and Combustion. Tsinghua University Press, Beijing (in Chinese).
available at www.sciencedirect.com
w w w . i i fi i r . o r g
journal homepage: www.elsevier.com/locate/ijrefrig
Modified two-fluid model for air curtains in open vertical display cabinets Ke-zhi Yua,b, Guo-liang Dinga,*, Tian-ji Chenb a
Institute of Refrigeration and Cryogenics, Shanghai Jiaotong University, Shanghai 200030, China Department of Refrigeration Engineering, Shanghai Fisheries University, Shanghai 200090, China
b
article info
abstract
Article history:
A modified two-fluid turbulence model is established to simulate the flow and heat transfer
Received 9 April 2007
characteristics of air curtains in an open vertical display cabinet. The air exiting the back
Received in revised form
panel (AEBP), the inner air curtain and the outer air curtain are taken as the first fluid
10 July 2007
and described with the standard k 3 turbulence model (KE). The air outside the display
Accepted 11 July 2007
cabinet is considered as the second fluid and calculated by the laminar model. Different
Published online 20 July 2007
from the existing two-fluid model, the mass transfer rate equation between the turbulent and the non-turbulent fluids is modified, and the face coefficient of honeycomb (the ratio of
Keywords:
the effective area to the total area of honeycomb) is fetched with the volume fraction of air
Refrigeration
curtain in order to conveniently depict the blow characteristics of the honeycomb. The
Commercial application
comparisons between the simulated temperatures and the experimental ones indicate
Display cabinet
that the modified two-fluid model (MTF) can give better agreement with the measurements
Open display cabinet
than KE and two original two-fluid models (TF1 and TF2). ª 2007 Elsevier Ltd and IIR. All rights reserved.
Air curtain Modelling Flow Heat transfer
Mode`le a` deux fluides modifie´ pour les rideaux d’air des meubles frigorifiques de vente verticaux Mots cle´s : Re´frige´ration ; Application commerciale ; Meuble de vente ; Meuble ouvert ; Rideau d’air ; Mode´lisation ; E´coulement ; Transfert de chaleur
1.
Introduction
Vertical open display cabinets are widely used in supermarkets. They can not only store various foods, but also attract customers’ attention so as to increase sales. In a vertical open display
cabinet, there are one or more forced-air curtains between the cabinet and the environment. However, these air curtains are easily disturbed by the ambient air, and then a more and less amount of ambient air is entrained into the display cabinet, which increases the energy consumption of the display cabinet.
* Corresponding author. Tel.: þ86 21 62932110; fax: þ86 21 34206814. E-mail address: [email protected] (G.-liang Ding). 0140-7007/$ – see front matter ª 2007 Elsevier Ltd and IIR. All rights reserved. doi:10.1016/j.ijrefrig.2007.07.006
international journal of refrigeration 31 (2008) 472–482
Nomenclature A Ae AEBP B c1, c2, cm cm1, cm2 cf, ct CFD cp D f F g G H I k KE l L m MTF p Pr
total area of honeycomb (m2) effective area of honeycomb (m2) air exiting the back panel cell size of honeycomb (m) constants in the k 3 turbulence model empirical constants for mass transfer rate (–) empirical constants for momentum and heat transfer, respectively (–) computational fluid dynamics specific heat capacity at constant pressure (J kg1 K1) phase diffusion coefficient in Eq. (1) (m2 s1) face coefficient of honeycomb (Ae/A) (–) inter-fluid friction force per unit volume (N m3 s1) gravitational acceleration (m s2) production rate of turbulent energy in the k 3 turbulence model (kg m1 s1) honeycomb height (m) inter-fluid sources in Eq. (1) turbulent kinetic energy (m2 s2) standard k 3 turbulence model Prandtl mixing length (m) honeycomb length (m) inter-fluid mass transfer rate per unit volume (kg m3 s1) modified two-fluid model based on Eq. (10) pressure (Pa) Prandtl number (–)
The honeycombs are generally installed at the top of the open display cabinet to form the air curtains. The honeycombs can decrease the turbulence intensities of air curtains, and reduce the mass and heat exchange between the air curtains and the ambient air, thus lessening the cooling load of the display cabinet. As the air curtains play a very important role in the cold preservation of display cabinets, various researches on air curtains have been made in recent years, including experimental studies and computational fluid dynamics (CFD) simulations. However, the mechanisms of flow, heat and mass transfer for air curtains are rather complicated. CFD simulation is a practical way to investigate the mechanisms of the air curtains and has been applied widely (Stribling et al., 1997; Ge and Tassou, 2001; Navaz et al., 2002; Wu et al., 2004; Cui and Wang, 2004; Fostera et al., 2005; D’Agaro et al., 2006; Cortella et al., 2001; Cortella, 2002). Among these studies, the most commonly used turbulent models are the k 3 model (Stribling et al., 1997; Ge and Tassou, 2001; Navaz et al., 2002; Wu et al., 2004; Cui and Wang, 2004; Fostera et al., 2005; D’Agaro et al., 2006) and the LES model (Cortella et al., 2001; Cortella, 2002). These models belong to single-fluid turbulence model, which takes the fluid in the whole computational domain as the turbulent fluid. But it is difficult for a single-fluid turbulence model to well describe the coexistence of highly turbulent flows inside the cabinet and largely laminar fluid outside the cabinet.
Q r S T TF1 TF2 U u v W
473
inter-fluid heat flux (J m3 s1) volume fraction (–) intra-fluid sources in Eq. (1) temperature ( C) two-fluid model based on Eq. (8) two-fluid model based on Eq. (9) velocity vector (m s1) velocity component in the x direction (m s1) velocity component in the y direction (m s1) honeycomb width (air curtain width) (m)
Greek symbols sK , s3 , sT , sr turbulent Prandtl number (–) m laminar viscosity (kg m1 s1) me effective viscosity (kg m1 s1) turbulence viscosity (kg m1 s1) mt r density (kg m3) d cell wall thickness (m) 3 dissipation rate of turbulent energy (m2 s3) f generic flow variable in Eq. (1) G within-phase diffusion coefficient in Eq. (1) (m2 s1) Subscripts 1 first fluid (turbulent air) 2 second fluid (non-turbulent air) amb ambient ex experimental value k fluid k i, j the spatial x, y direction in air curtain inlet
The two-fluid turbulence model provides an effective means of representing the coexistence of laminar flows and turbulent flows. This model considers the system as composed of two fluids, which coexist in time and space but possess different volume fractions. This model can reflect the exchange of mass, momentum and energy between the two fluids and is well suited to describe the whole computational domain including the high turbulence region inside the cabinet and relatively quiescent region outside the cabinet. A two-fluid model was adopted to simulate the air curtains for vertical display cabinets for the first time in Yu et al. (2007), in which the computational fluid was assumed to be composed of the turbulent fluid (air curtains) and the non-turbulent fluid (the ambient air outside the display cabinet). It showed that the two-fluid model can give better agreement with the measured values than the k 3 model. However, the different mass transfer rate relations between the two fluids were adopted by different researchers in different cases (Malin and Spalding, 1984; Shen et al., 2003; Markatos et al., 1986; Ilegbusi, 1994; Markatos and Kotsifaki, 1994; Fan, 1988; Sheng and Jonsson, 2000; Zhou, 1991). The mass transfer rate equation employed by Yu et al. (2007) assumed that only the turbulent fluid entrains the non-turbulent fluid; nevertheless the turbulent fluid cannot be transferred into the non-turbulent fluid. However, the experiments showed that the cold turbulent air may spill from the display cabinet into the ambient air
474
international journal of refrigeration 31 (2008) 472–482
back of the cabinet. A small fraction of cooled air is fed into the cabinet through the back panel of the cabinet, while more cold air is blown through two honeycombs to form two air curtains. The two air curtains include an inner air curtain with a high velocity on the inside (cabinet side) and an outer air curtain with a low velocity on the outside (store side). The inner air curtain and the outer air curtain, together with the entrained ambient air, are drawn back into the cabinet through the respective air curtain outlets at the bottom of the display cabinet. The inner air curtain is re-cooled by the evaporator while the outer air curtain bypasses the evaporator in the air flow tunnel, which makes the temperature of the inner air curtain lower than that of the outer air curtain. When the two-fluid model is applied to the vertical display cabinet, the following assumptions are made:
(Stribling et al., 1997; Navaz et al., 2002; D’Agaro et al., 2006; Cortella et al., 2001; Cortella, 2002). Therefore the entrainment rate between the two fluids should be modified according to the mass transfer mechanics of display cabinets. Furthermore, the difference between effective area (pass area) and total area of the air curtain honeycomb was not taken into account in Yu et al. (2007), which led to overestimating the momentum and the energy of air curtain discharge. Therefore the velocity, temperature and fraction volume of air curtain should be revised to depict air curtain honeycomb properly. The purpose of this paper is to provide a modified two-fluid model of air curtains, which can better predict the air flow field inside vertical display cabinets. In the modified model, the new mass transfer rate equation between the turbulent fluid and the non-turbulent fluid is proposed, and the face coefficients of honeycombs are obtained with the volume fractions of air curtains in order to depict the blow characteristics of the honeycomb properly.
(1) The system is assumed to consist of two fluids. The air exiting the back panel (AEBP), the inner and the outer air curtains are defined as the first fluid (turbulent fluid). The ambient air is defined as the second fluid (non-turbulent fluid). (2) A steady two-dimensional model is employed regardless of the display cabinet length. (3) The standard k 3 two equation model (KE) is adopted to simulate the turbulent fluid.
2. The two-fluid model for display cabinets (Yu et al., 2007) In a vertical open display cabinet (shown in Fig. 1), the air is forced by propeller fans to flow through an evaporator at the
Honeycomb y
Inner air curtain flow
x
Back panel
Outer air curtain flow
Outer air curtain inlet
Inner air curtain inlet
S1 20 14=280mm
20 14=280mm
S2
Evaporator
20 14=280mm
S3
T-type thermocouples
S4 Air flow tunnel
20 14=280mm
Inner air curtain outlet
Outer air curtain outlet
50 3=150mm
O1 O2 O3 O4
Fan
20 10=200mm
Fig. 1 – A vertical open display cabinet and the location of T-type thermocouples mounted.
475
international journal of refrigeration 31 (2008) 472–482
(4) Radiation heat transfer is not taken into consideration. (5) The Boussinesq hypothesis is adopted due to the density varying with aerostatic pressure.
Q¼
ct cp r1 r1 r2 ðT1 T2 Þju1 u2 j l
(6)
in which l is the Prandtl mixing length, expressed as: The two-fluid model for the display cabinet can be established and the governing transport equations can be written with reference to a generic flow variable fk (for fluid k) as: v v vf v vrk rk Gfk k þ fk Dfk þ Sfk þ Ifk rk rk Ukj fk ¼ vxj vxj vxj vxj vxj
(1)
where subscript k denotes fluid k, which can be set to 1 or 2 and subscript j refers to spatial coordinates. G is the diffusion coefficient within fluid 1 or fluid 2 and D is the phase diffusion coefficient between fluid 1 and fluid 2. S is the source term within fluid 1 or fluid 2 and I is the source term between fluid 1 and fluid 2. The continuity, momentum and energy equation can be obtained by setting the generic flow variable fk to unity, velocity and temperature, respectively. The governing equations for turbulent and non-turbulent fluid are shown in Tables 1 and 2, where: mt ¼
cm rk2 3
me ¼ m þ m t
G ¼ mt
vU1i vU1i vU1j þ vxj vxj vxi
3=2 c3=4 m k 3
(7)
Two popular equations of mass transfer rate between the two fluids are depicted as follows (Yu et al., 2007; Malin and Spalding, 1984; Shen et al., 2003; Markatos et al., 1986; Ilegbusi, 1994; Markatos and Kotsifaki, 1994): m¼
cm1 r1 r1 r2 ju1 u2 j l
(8)
and (Sheng and Jonsson, 2000),
m¼
cm2 r1 r1 r2 ðr2 0:5Þju1 u2 j l
(9)
(2)
The constants cf, ct, cm1, and cm2, employed in Eqs. (5), (6), (8) and (9), are 0.05, 0.05, 0.1 and 10.0, respectively.
(3)
3.
Modifications of the two-fluid model
3.1.
The modification of mass transfer rate equation
(4)
The empirical constants in the two-fluid model are listed in Table 3. The momentum and energy transfer rates are generally expressed as: Fi ¼
l¼
cf r1 r1 r2 ðU1i U2i Þju1 u2 j l
(5)
The mass transfer rate equation between the two fluids plays a very important role in the two-fluid model. Eqs. (8) and (9) are the most widely used relations. In Eq. (8) m is always positive, which means that only fluid 2 (non-turbulent fluid) can be entrained by fluid 1 (turbulent fluid). However, the mass transfer rate m may be negative according to Eq. (9). An additional factor r2 0.5 included in Eq. (9) allows for the equally entrainment between the turbulent fluid
Table 1 – The governing equations of turbulent fluid Equation
fk
Gfk
Dfk
Sfk
Ifk
Continuity
1
0
me sr
Momentum
U1i
me
me sr
Energy
T1
mt m þ sT Pr
me sr
0
Kinetic energy
k
mt þm sk
me sr
r1 ðG r1 3Þ
0
Dissipation rate
3
mt þm s3
me sr
r1 3ðc1 G c2 r1 3Þ k
0
0
m
vU1j v vp þ r1 ðr1 r2 Þgi me r1 r1 vxj vxi vxi
U2i m Fi
T2 m
Q cp
476
international journal of refrigeration 31 (2008) 472–482
Table 2 – The governing equations of non-turbulent fluid Equation
fk
Gfk
Dfk
Continuity
1
0
m sr
Momentum
U2i
m
m sr
Energy
T2
m Pr
m sr
(fluid 1) and the non-turbulent fluid (fluid 2). From the viewpoint of physics, the entrainment rate of a non-turbulent fluid by a turbulent fluid is much more than that of the turbulent fluid by the non-turbulent fluid. Therefore both Eqs. (8) and (9) have some disadvantages. A new mass transfer rate equation should be developed. Based on Eqs. (8) and (9), the modification of mass transfer rate is established by using the weighted average of the two equations. The weighted factors can be assigned to the volume fractions of the two fluids (r1 and r2). The modified mass transfer rate equation can be expressed as: cm1 r1 r1 r2 ju1 u2 j cm2 r1 r1 r2 ðr2 0:5Þju1 u2 j $r1 þ $r2 l l 2 2 cm1 r1 r1 r2 ju1 u2 j cm2 r1 r1 r2 ðr2 0:5Þju1 u2 j þ ¼ l l
m¼
(10)
Fig. 2 shows the relationship of the dimensionless number m=ðr1 ju1 u2 j=lÞ vs. the volume fraction of fluid 1 (r1) based on Eqs. (8)–(10). From the figure, it can be seen that Eq. (8) neglects the mass transfer rate m from fluid 1 to fluid 2, whereas Eq. (9) over predicts the mass transfer rate. It seemed that Eq. (10) can give reasonable result. Furthermore, Fig. 2 shows the mass transfer rate m from fluid 1 to fluid 2 less than m from fluid 2 to fluid 1 based on Eq. (10). This is consistent with the physics of the problem. Thus Eq. (10) can reflect the mass transfer rate between the two fluids perfectly. Besides the standard k 3 model (KE), the two-fluid models based on Eq. (8) (two-fluid model 1, TF1), Eq. (9) (two-fluid model 2, TF2), and Eq. (10) (modified two-fluid model, MTF) are adopted to simulate the open vertical display cabinet so that the four models can be compared with each other.
3.2.
Sfk
Ifk
0
m
vU2j v vp mr2 r2 vxj vxi vxi
U2i m þ Fi
0
T2 m þ
Ae < A
Q cp
(11)
If the experimental velocity and temperature of air curtain are uex and Tex, respectively, the inlet volume fraction of air curtains is equal to unity according to Yu et al. (2007). Then the theoretical inlet momentum and energy can be expressed as: r1 u1 Ar1 $u1 ¼ r1 Au21 ¼ r1 Au2ex
(12)
and r1 cp u1 Ar1 $T1 ¼ r1 cp Au1 T1 ¼ r1 cp Auex Tex
(13)
However, the actual inlet momentum and inlet energy are r1 Ae u2ex and r1 cp Ae uex Tex , respectively. The theoretical inlet momentum and energy are larger than the actual inlet momentum and energy, respectively, according to Eq. (11). The face coefficient of the honeycomb f can be defined by the ratio of effective area Ae to the total area A, namely: f ¼ Ae/A. If the volume fraction of the turbulent fluid (fluid 1) is assigned as the face coefficient of the honeycomb f, then the theoretical inlet momentum and energy can be expressed as: r1 u1 Ar1 $u1 ¼ r1 r1 Au21 ¼ r1 fAu21 ¼ r1 Ae u2ex
(14)
and r1 cp u1 Ar1 $T1 ¼ r1 cp r1 Au1 T1 ¼ r1 cp fAu1 T1 ¼ r1 cp Ae uex Tex
(15)
Therefore, the theoretical inlet momentum and energy satisfy the actual inlet momentum and energy, respectively.
The modification of inlet parameters of air curtain
For a honeycomb shown in Fig. 3, the total area A is larger than the effective area Ae, namely:
Table 3 – Values of the constants employed in the two-fluid model cm 0.09
c1
c2
sk
s3
sr
sT
1.44
1.92
1.0
1.3
0.95
1.0
Fig. 2 – The relationship of the dimensionless number m=ðr1 ju1 Lu2 j=lÞ vs. the volume fraction of fluid 1 (r1) based on Eqs. (8)–(10).
477
The computational domain is constructed, in which the evaporator, the insulation layer, the propeller fan and the product
Table 4 – The modifications of the inlet parameters of the air curtains for display cabinet Fluid 1
Before modification After modification
Fluid 2
u1
v1
T1
r1
u2
v2
T2
r2
uex uex
0 0
Tex Tamb
1 f
uex 0
0 0
Tex Tamb
0 1f
0.1 0.1 0.1 0.1 0.9 0.9 0.9 0.9 25 25 25 25 6.1 7.2 7.1 6.2 0 0 0 0 0 0 0 0 0 0 0 0 0.25 0.25 0.25 0.28 0.1 0.1 0.1 0.1 0.9 0.9 0.9 0.9 25 25 25 25 5.3 4.4 6.4 5.1 0 0 0 0 0 0 0 0 0 0 0 0 0.52 0.58 0.3 0.52 6.1 7.2 7.1 6.2
v2 (m s1) v2 (m s1) v1 (m s1) u2 (m s1) u1 (m s1) Temperature ( C) Velocity (m s1)
The computational domain and grid
Temperature ( C)
4.2.
Velocity (m s1)
Four cases are simulated for a display cabinet. The specifications of the cases are presented in Table 5. Case 1 is the basic working condition for a display cabinet. Case 2 and Case 3 examine the effect of the inner air curtain velocity on the performance of the display cabinet. Case 4 examines the effect of the outer air curtain velocity on the performance of the display cabinet. The temperatures of the air curtains are assigned according to the measurements in order to be compared with the experimental data.
Inner air curtain
The computational case
Outer air curtain
4.1.
Inner air curtain
Case study
Table 5 – Case specifications and boundaries of air curtains inlet
4.
T1 ( C)
In addition, the inlet velocity and the inlet temperature of fluid 2 (the non-turbulent fluid) are assigned to zero and the ambient temperature, respectively. The modified parameters are listed in Table 4.
0.25 0.25 0.25 0.28
(17)
Case specifications
r2 ¼ 1 f
r2 (–)
The volume fraction of the non-turbulent fluid (fluid 2) can be assigned by Eq. (17).
r1 (–)
(16)
Case
r1 þ r2 ¼ 1
T2 ( C)
The constraint of the volume fraction is:
u1 (m s1)
The boundaries of air curtains inlet
Fig. 3 – The honeycomb of air curtain for display cabinet.
5.3 4.4 6.4 5.1
Honeycomb
v1 (m s1)
air curtain inlet
Cell structure
u2 (m s1)
Outer air curtain
H
B
T1 ( C)
W
Ae
0.52 0.58 0.3 0.52
A
1 2 3 4
T2 ( C)
L
r1 (–)
r2 (–)
international journal of refrigeration 31 (2008) 472–482
478
international journal of refrigeration 31 (2008) 472–482
(2) The air exiting the back panel (AEBP) The velocity of AEBP is about 0.1 m s1 by the experiments. The temperature of AEBP is equal to the temperature of the inner air curtain. In addition, the turbulent kinetic energy and the turbulent dissipate rate of AEBP can be the same as those assigned to the air curtains. (3) Ambient air The influence of the air curtains on the extended boundaries can be ignored. The temperature and the volume fraction of the ambient air are assigned to 25 C and 1.0, respectively. (4) Air curtain outlet Flow rate condition is imposed at the air curtain outlet. The outer air curtain outlet intakes the flow rate of the outer air curtain, whereas the inner air curtain outlet intakes the flow rate sum of the inner air curtain and the air exiting the back panel (AEBP). (5) Wall The non-slip boundary condition (zero velocity) is adopted at the solid surfaces. The wall function is employed for the turbulent fluid.
shelves are removed in order to simplify the simulation. As the open boundary conditions of the display cabinet are uncertain, the computational domain is extended until the effect of the display cabinet opening is negligible there. Then the computational domain becomes a 4.0 m 4.0 m irregular region, as shown in Fig. 4. The asymmetric grid consists of 61 63 cells, which is generated to mesh the computational domain, as illustrated in Fig. 4. The denser grid is adopted in the air curtain inlet and outlet regions where the air flow is complex. The grid independence is performed by using the progressively finer grid until the average temperature difference at the air curtain outlet is less than 0.1 C.
Boundary conditions
(1) Air curtain inlet The honeycomb length, height, width, cell size and cell wall thickness (shown in Fig. 3) are 2500 mm, 65 mm, 76 mm, 4.0 mm, and 0.2 mm, respectively. It can be calculated that the face coefficient of honeycomb is 0.9. The velocities, the temperatures and the volume fractions of the inner and the outer air curtains are shown in Table 5 according to the modified two-fluid model. The turbulent kinetic energy of the air curtains can be set to be 1% of the average kinetic energy (Cui and Wang, 2004). The turbulent dissipate rate is determined by Eq. (18) (Cui and Wang, 2004). 500k2 cm uin L
Simulation tool
The modified two-fluid model is embodied in the FORTRAN program developed by the present authors. The calculation needed 10,000 iterations in order to obtain a stable and converged solution. This required about 5 h on a personal computer with 2.09 GHz CPU.
(18)
5.
where cm is a constant which is shown in Table 3, uin is the average inlet velocity and L is the characteristic dimension of the air curtain (here it is the width of the honeycomb). y
Results and discussion
T-type thermocouples were mounted around the display cabinet to validate the simulation models (KE, TF1, TF2 and MTF).
The air curtain inlet
x
Unit:mm
76 76
200
50
500
The ambient space
210
1435
40
210
5 5
40
900
5
225
4000
90
5
The air exiting the back panel
450
330
70
900
The air curtain outlet
Fig. 4 – The computational domain and grid.
4000
The display cabinet
40
5
210
40
Product shelf
1000
3¼
4.4.
1435
4.3.
international journal of refrigeration 31 (2008) 472–482
479
Fig. 5 – Comparison of the predicted temperatures by KE, TF1, TF2 and MTF with the experimental temperatures outside shelves (a) S1 (b) S2 (c) S3, and (d) S4.
Fig. 1 shows that the measuring points can be divided into two groups. A group of measuring points (S1–S4) was located outside the product shelves where the temperature gradient is expected large, whereas another group of measuring points (O1–O4) was located near the air curtain outlet of the display cabinet where the cold air may overspill the display cabinet. All the measurements were recorded through the temperature acquisition logger at intervals of 10 s with the average value being calculated to be compared with the predictions.
5.1. Comparison with the measurements outside the shelves Fig. 5(a)–(d) shows the predicted temperatures by four models (KE, TF1, TF2 and TFM) and the experimental temperatures outside the four shelves (the first, the second, the third and the fourth shelf) for Case 1, respectively. The comparisons for Case 2, Case 3 and Case 4 are similar. It can be concluded from the figures as follows.
(1) The variation tendency of the simulated temperatures by four models is in agreement with that of the experimental data. However, the difference between the measurements and the predictions by four models is relatively large for the fourth shelf (S4). (2) Since KE assumes that the air in the whole computational domain is a strong turbulent fluid, the simulated temperatures outside the cabinet are very close to the ambient temperature 25 C, which are higher than the measurements (21–23 C). On the contrary, at the shelf the predicted temperatures (subzero) are lower than the measurements (about 3 C) when used by KE owing to overestimating the inlet momentum and energy of the cold air curtain. (3) The mass transfer rate m from fluid 1 to fluid 2 is not taken into account in TF1 (see Fig. 2). The flow field simulated by TF1 is approximate to the turbulence, thus the predicted temperatures by TF1 are very close to that by KE outside the cabinet. The differences between the predictions by TF1 and those by KE are less than 0.2 C for the most measuring points. However, TF1 predicts higher
480
international journal of refrigeration 31 (2008) 472–482
temperatures than KE when the measuring points locate very near the shelves. The maximum difference between the predictions by KE and TF1 can be found at the product shelves. (4) TF2, which over predicts the mass transfer rate m from fluid 1 to fluid 2 (see Fig. 2), predicts lower temperatures than TF1 and KE. It can be seen that the simulated temperatures by TF2 are lower than those by TF1 and the differences are 0.9–5.5 C. (5) When the measuring points locate at the first shelf (S1), the maximum differences between the experimental temperatures and the predicted values from TF1 and TF2 are 3.3 C (S1) and 4.6 C (S1), respectively. Thus the predicted temperatures by TF1 are in better agreement with the experimental values than those by TF2 when the thermocouples are close to the shelves. (6) However, when the measuring points are far away the shelf, the average differences between the measured temperatures and the predictions from TF1 and TF2 are 1.3 C and 1.1 C, respectively. It can be inferred that TF2 shows more agreement with the measurements when the measuring points locate far away from the shelves.
(7) The modification effect of TFM is significant. It is compatible with both TF1 and TF2 (see Fig. 2). When the measuring points locate at the shelf, the maximum difference of MTF (2.0 C for S3) is less than that of TF1 (3.3 C for S1) and that of TF2 (4.6 C for S1). Moreover, the average difference of MTF (1.0 C) is less than that of TF1 (1.3 C) and that of TF2 (1.1 C) when the thermocouple is far away the shelf. Therefore, MTF shows good agreement with the measurements wherever the experimental points locate.
5.2. Comparison with the measurements outside the display cabinet The predicted temperatures by the four models (KE, TF1, TF2 and TFM) and the experimental temperatures outside the display cabinet for Case 2 are shown in Fig. 6(a)–(d), respectively. The comparisons for Case 1, Case 3 and Case 4 are similar. It can be inferred from the figures: (1) KE, which takes the air outside the display cabinet as the turbulent fluid, predicts higher temperatures (the maximum is 16.8 C) than the measurements which range
Fig. 6 – Comparison of the predicted temperatures by KE, TF1, TF2 and MTF with the experimental temperatures outside the display cabinet (a) O1, (b) O2, (c) O3 and (d) O4.
481
international journal of refrigeration 31 (2008) 472–482
Table 6 – The average and the maximum temperature differences between the measurements and the predictions using the kL3 model and three (modified) two-fluid models ( C) Case
Case 1
Case 2
Case 3
Case 4
(2)
(3)
(4)
(5)
Measuring location
Average differences KE
Maximum differences
TF1
TF2
MTF
KE
TF1
TF2
MTF
Outside the shelf
S1 S2 S3 S4
1.2 1.8 2.0 3.6
1.2 1.9 1.8 2.9
1.1 0.8 1.2 0.3
0.2 0.3 0.3 0.3
3.9 4.0 5.0 5.8
3.4 3.8 4.8 5.6
4.6 4.0 5.1 2.5
1.5 2.1 3.6 3.1
Outside the cabinet
O1 O2 O3 O4
3.2 2.1 0.9 1.2
2.1 3.3 1.5 1.8
1.5 0.8 0.8 0.7
1.8 1.0 0.6 0.5
4.7 2.9 1.8 2.1
3.7 4.1 2.6 2.7
3.1 1.5 1.8 1.8
3.3 1.8 1.6 1.4
Outside the shelf
S1 S2 S3 S4
1.0 1.1 1.7 2.9
0.9 1.1 1.5 1.9
1.3 1.5 0.8 1.8
0.1 0.5 0.0 0.7
4.0 3.6 4.7 4.9
3.3 3.3 4.3 4.5
4.3 4.9 4.2 4.5
1.3 3.7 2.2 4.4
Outside the cabinet
O1 O2 O3 O4
2.2 2.0 0.4 1.0
1.4 3.7 3.0 1.6
1.0 1.2 0.3 2.2
0.6 0.7 0.7 1.4
4.4 3.0 1.4 2.4
3.7 4.5 4.1 4.0
3.2 2.2 1.6 3.7
2.8 1.6 1.8 2.9
Outside the shelf
S1 S2 S3 S4
1.9 2.3 2.6 1.0
1.8 2.6 2.1 0.1
0.2 0.1 0.1 2.4
1.1 0.5 0.5 1.7
3.3 4.2 4.6 4.3
3.4 4.7 3.9 4.5
2.5 2.2 2.9 5.3
3.2 1.9 1.8 4.8
Outside the cabinet
O1 O2 O3 O4
2.0 1.3 1.8 2.4
2.9 1.5 0.8 0.2
1.2 2.0 1.2 1.3
0.7 1.5 1.0 0.3
3.8 3.8 3.6 3.9
4.7 4.5 2.5 1.6
2.9 4.5 2.9 2.9
2.5 4.1 2.6 1.7
Outside the shelf
S1 S2 S3 S4
2.2 1.7 2.2 1.8
2.2 1.7 2.0 0.8
0.4 0.7 0.8 2.2
0.8 0.3 0.0 1.8
5.5 3.7 5.0 2.9
5.4 3.4 4.7 2.7
4.9 3.4 4.1 4.3
3.2 1.4 2.4 4.5
Outside the cabinet
O1 O2 O3 O4
1.3 1.0 0.9 0.2
2.1 1.2 1.9 0.8
1.0 0.7 0.3 0.6
0.7 0.6 0.2 0.0
2.6 3.2 3.2 1.8
3.5 3.7 4.2 3.1
2.4 3.2 2.9 2.6
2.6 3.2 2.9 1.8
from 11.7 C to 15.7 C for the first row of measuring points (O1, see Fig. 6(a)). KE shows similarity for O2, O3 and O4. The flow field simulated by TF1 is approximate to the turbulence, therefore TF1 over predicts the temperatures outside the display cabinet. The maximum differences between the predictions by TF1 and the measurements can be 4.5 C (O2). It is abnormal that TF1 gives abnormally lower temperatures (14.6–16.0 C) at the first row of measuring points (O1, see Fig. 6(a)). This phenomenon may be caused by the overspill of cold air from the cabinet into the ambient. Although KE and TF1 simulate the higher temperatures outside the cabinet, KE gives lower predictions than TF1 except the first row of measuring points (O1) owing to over predicting the inlet momentum and energy of the cold air curtain. The differences between the predictions by TF1 and those by KE are in the range of 0.8 C (O1) to 2.7 C (O3 and O4). The temperatures outside the cabinet predicted by TF2 and MTF are lower than those by KE and TF1 except the first row of measuring points (O1). The maximum differences
between the measurements and the predictions by KE, TF1, TF2 and MTF are 4.4 C, 4.5 C, 3.2 C, and 2.8 C, respectively. It is apparent that MTF and TF2 are in agreement with the measurements and MTF shows better agreement with the experimental temperatures than TF2.
5.3. The average and the maximum temperature differences by the four models Table 6 summarizes the average and the maximum temperature differences for the four cases simulated by the four models. From the table it can be found that the average and the maximum temperature differences are the smallest for most cases when MTF is employed (which are bolded in Table 6), which demonstrated that MTF is a valuable simulation model for display cabinet. It can also be found in Table 6 that TF2 shows better agreement with the measurements than MTF for Case 3 when measuring points were located outside the product shelves (S1–S3). The velocities of the air curtains are so small (0.30 m s1 and 0.25 m s1) for Case 3 that the turbulent fluid
482
international journal of refrigeration 31 (2008) 472–482
(fluid 1) and the non-turbulent fluid (fluid 2) can be equally entrained each other. TF2 shows greater advantage than MTF. Therefore MTF is not the appropriate numerical model for display cabinet when the velocities of the air curtains are quite small. However, there are maximum temperature differences of 4.8 C in the MTF simulation from Table 6. The reasons mainly come from two aspects. Firstly, the two-dimensional model assumption is made in applying the modified two-fluid model to the simulation of display cabinet. The two-dimensional model is not able to simulate the influence of the longitudinal direction on the display cabinet. Secondly, the measuring temperatures outside the display cabinet may have some extra uncertainty, which is caused by the interaction of two airstreams including the cold air spill over the display cabinet and the hot air entrained by the fan.
6.
Conclusions
A modified two-fluid model has been developed to simulate air flow and heat transfer in an open vertical display cabinet. The air exiting the back panel (AEBP), the inner air curtain and the outer air curtain are defined as fluid 1 (turbulent fluid), and the ambient air outside display cabinet is defined as fluid 2 (non-turbulent fluid). The volume fraction, velocity vector and temperature field are simulated by solving the twodimension steady transport equations of each fluid. The main conclusions of the numerical investigation are as follows: (1) In the modified two-fluid model, a modified mass transfer rate equation between the two fluids is established. It not only describes the entrainment of the non-turbulent fluid by the turbulent fluid, but also reveals the mass transfer rate from the turbulent fluid to the non-turbulent fluid. Thus the modified mass transfer rate equation can reflect the proper mass transfer rate between the air curtains and the ambient air in an open vertical display cabinet. (2) In the modified two-fluid model, the volume fractions of the first fluid (air curtains) are assigned as the face coefficients of honeycombs to depict the honeycombs, which can provide proper discharge momentum and energy of air curtains. (3) The comparisons with the experimental temperatures indicate that the agreement with the measurements can be improved by using the modified two-fluid model instead of the standard k 3 model (KE) and original two-fluid model (TF1 and TF2) in simulation of display cabinets.
references
Cui, J., Wang, S., 2004. Application of CFD in evaluation and energy-efficient design of air curtains for horizontal
refrigerated display cases. International Journal of Thermal Sciences. 43, 993–1002. Cortella, G., Manzan, M., Comini, G., 2001. CFD simulation of refrigerated display cabinets. International Journal of Refrigeration 24, 250–260. Cortella, G., 2002. CFD-aided retail cabinets design. Computers and Electronics in Agriculture 34, 43–66. D’Agaro, P., Cortella, G., Croce, G., 2006. Two- and three dimensional CFD applied to vertical display cabinets simulation. International Journal of Refrigeration 29, 178–190. Fostera, A.M., Madgeb, M., Evansa, J.A., 2005. The use of CFD to improve the performance of a chilled multi-deck retail display cabinet. International Journal of Refrigeration 28, 698–705. Fan, W.-C., 1988. A two-fluid model of turbulence and its modifications. Science in China Series A: Mathematics, Physics, Astronomy and Technological Sciences 31 (1), 79–86 (in Chinese). Ge, Y.T., Tassou, S.A., 2001. Simulation of the performance of single jet air curtains for vertical refrigerated display cabinets. Applied Thermal Engineering 21, 201–219. Ilegbusi, O.J., 1994. Application of the two-fluid model of turbulence to tundish problems. ISIJ International 34 (9), 732–738. Malin, M.R., Spalding, D.B., 1984. A two-fluid model of turbulence and its application to heated plane jets and wakes. Physicochemical Hydrodynamics 5 (5/6), 339–362. Markatos, N.C., Pericleous, K.A., Cox, G., 1986. A novel approach to the field modeling of fire. Physicochemical Hydrodynamics 7 (2/3), 125–143. Markatos, N.C., Kotsifaki, C.A., 1994. One-dimensional, two-fluid modeling of turbulent premixed flames. Applied Mathematical Modelling 18, 646–657. Navaz, H.K., Faramarzi, R., Gharib, M., Dabiri, D., Modarress, D., 2002. The application of advanced methods in analyzing the performance of the air curtain in a refrigerated display cabinet. Journal of Fluids Engineering: Transactions of the ASME 124 (3), 756–764. Stribling, D., Tassou, S.A., Marriott, D., 1997. A two-dimensional CFD model of a refrigerated display case. ASHRAE Transactions 103 (1), 88–94. Shen, Y.M., Ng, C.-O., Chwang, A.T., 2003. A two-fluid model of turbulent two-phase flow for simulating turbulent stratified flows. Ocean Engineering 30, 153–161. Sheng, D.Y., Jonsson, L., 2000. Two-fluid simulation on the mixed convection flow pattern in a nonisothermal water model of continuous casting tundish. Metallurgical and Materials Transactions B: Process Metallurgy and Materials Processing Science 31 (4), 867–875. Wu, Y., Xie, G., Chen, Z., Niu, L., Sun, D.W., 2004. An investigation on flowing patterns of the airflow and its characteristics of heat and mass transfer in an island open display cabinet with goods. Applied Thermal Engineering 24, 1945–1957. Yu, Ke-zhi, Ding, Guo-liang, Chen, Tian-ji, 2007. Simulation of air curtains for vertical display cases with a two-fluid model. Applied Thermal Engineering 27, 2583–2591. Zhou, L.X., 1991. The Numerical Simulation of Turbulent TwoPhase Flow and Combustion. Tsinghua University Press, Beijing (in Chinese).