- Email: [email protected]
Modeling of frost formation over parallel cold plates considering a two-dimensional growth rate
International Journal of Thermal Sciences 104 (2016) 245e256
Contents lists available at ScienceDirect
International Journal of Thermal Sciences journal homepage: www.elsevier.com/locate/ijts
Modeling of frost formation over parallel cold plates considering a two-dimensional growth rate n b, K.A.R. Ismail a J.M. Armengol a, C.T. Salinas a, *, J. Xama a b
Faculty of Mechanical Engineering, University of Campinas (UNICAMP), CP 6066, CEP 13081-970 Campinas, SP, Brazil n y Desarrollo Tecnolo gico, CENIDET-TNM-SEP Prol, Av. Palmira S/N. Col. Palmira, Cuernavaca, Morelos CP 62490, Mexico Centro Nacional de Investigacio
a r t i c l e i n f o
a b s t r a c t
Article history: Received 14 May 2015 Received in revised form 7 January 2016 Accepted 15 January 2016 Available online 19 April 2016
In this study, a frost formation model is presented based on a new two-dimensional approach for the growth rate. For modeling the frost formation over parallel cold plates, the basic transport equations of mass, energy and momentum have been discretized using the finite volume method in a twodimensional domain in which air and frost are considered. A fixed grid formulation is used to deal with the airefrost moving boundary. An extended domain in the inlet boundary has been considered in order to study the frost formation in the leading edge of the plate. The numerical results have been validated against experimental data in which frost growth and temperature as a function of time are reported as local values. The model predictions of the frost thickness as a function of time agree with the experimental data within 10% of deviation for the case of intermediate plate temperature. © 2016 Elsevier Masson SAS. All rights reserved.
Keywords: Frost Parallel plates Fin-and-tube heat exchanger
1. Introduction Frost formation consists in a phenomenon in which water vapor from moist air flow is deposited in a cold surface whose temperature is lower than the freezing temperature of water. Frost is a porous structure which consists of ice containing air gaps, so it presents a considerable thermal resistance. In addition, frost formation in the heat exchanger channels increases the pressure drop by flow restriction. Therefore, heat exchangers thermal performance decreases when frost formation appears. Many researchers have studied this phenomenon. One of the first among them was Hayashi et al. [1], who divided the frost growth process in three different periods. The first one is associated with the crystal growth, which is relatively short when compared with the whole process, the frost layer does not experience a significant thickness increase in this period. During this early stage, the frost cannot be considered a porous medium since the structure might be understood as ice columns where convective heat and mass transfer are the main mechanisms for frost growth, rather than diffusion within the frost layer.
* Corresponding author. E-mail address: [email protected] (C.T. Salinas). http://dx.doi.org/10.1016/j.ijthermalsci.2016.01.017 1290-0729/© 2016 Elsevier Masson SAS. All rights reserved.
In the second stage, known as the frost layer growth period, frost is treated as a porous medium in which molecular diffusion of water vapor occurs. The total mass flux of water vapor coming from the air contributes in two different processes: the densification of the frost layer and the frost growth. Finally, in the frost layer full growth period, the frost surface temperature reaches the triple point of water, and so it begins a cyclic process in which the condensed vapor at the frost surface permeates into the frost layer and gradually freezes because of the inner temperature gradient. A study of this stage is presented in Aoki et al. [2]. In the last decades, many models for frost growth on flat plates have been proposed. A review of frost formation in simple geom etries is found in Oneal and Tree [3]. One can classify the theoretical models into two groups, depending on the use or not of empirical correlation for the airefrost interaction. In the first group, the models use empirical correlations for the air side effect and they apply energy and mass transfer equations within the frost layer to predict frost formation. Tao et al. [4] developed a mathematical model based on the frost structure analysis done by Hayashi et al. [1] using the local volume averaging technique to deal with the frost as a porous media. Based on that model, Ismail and Salinas [5] determined the best initial values of the diffusivity, initial radius and geometry of ice crystals. Using the local volume averaging technique, Le Gall et al. [6] proposed an interface balance to predict frost growth rate. In
246
J.M. Armengol et al. / International Journal of Thermal Sciences 104 (2016) 245e256
Nomenclature m_ water vapor mass flux, kg m2 s1 A area, m2 cp specific heat capacity, J kg1 K1 D mass diffusivity coefficient, m2 s1 E, W, N, S respective east, west, north and south computational nodes e, w, n, s respective east, west, north and south control volume faces H distance, m P pressure, Pa qsub specific heat of sublimation, J kg1 T temperature, C t time, s u horizontal velocity, m s1 V volume, m3 v vertical velocity, m s1 x coordinate, m y coordinate, m
their work, Na and Webb [7e9] proposed a model considering supersaturated water vapor density at the frost surface. Yang et al. [10] reported an effective frosting model for predicting the thermal performance of a finetube heat exchanger. More recently, Hermes et al. [11] developed a mathematical model whose predictions of the frost thickness agreed within 10% error bands when compared with experimental data. Kandula [12,13] proposed new frost correlations for the flat plate and studied the effects of the environment parameters in the frost characteristics. Even though one-dimensional models based on empirical correlations can predict frost formation with reasonable precision, they require experimental determination of the empirical correlations for every different geometry. The second group, consists in an attempt to avoid empirical correlations by simulating the air flow and the frost formation simultaneously. As frost formation is a common phenomenon which takes place in different surface geometries, it is required a better understanding of the two-dimensional frost growth without using empirical correlations. Lee et al. [14] proposed a twodimensional model that shows good agreement with experimental data for different working conditions. However, the model assumes that the frost density in the direction normal to the cooling plate is constant. Recently, Cui et al. [15] and Kim et al. [16] simulated frost formation with a commercial Computational Fluid Dynamic (CFD) package getting results with good agreement with experimental data. These models can determine local properties of the frost layer such as density, humidity and temperature. In these models, it was found that frost layer at the leading edge of the plate presents larger values of density and thickness when compared with the values at the trailing edge. Lenic et al. [17] proposed a twodimensional mathematical model based on the assumptions of Na and Webb [7] and modeled it using the finite volume technique. The model of Lenic is similar to the one proposed in this work, except for the treatment of the frost growth rate. The main goal of this investigation is to establish a mathematical model for the frost growth in two dimensions, which means that a frost element not only increases its height but also its width. Such approach is developed in the context of the second group, so it has been modeled the air and the frost subdomains.
Greek symbols l thermal conductivity, W m1 K1 m dynamic viscosity, Pa s u humidity ratio, kgv kg1 a r density, kg m3 3 porosity Subscripts 0 initial value diff diffusing through the frost layer a air CS control surface ef effective f frost i ice in inlet s plate surface v vapor
2. Mathematical and physical model 2.1. Physical model of the fin-and-tube heat exchanger The frost formation is here modeled in the region between parallel cold plates. Such configuration represents the air channels of a fin-and-tube heat exchanger, which is very common in refrigeration systems. In the heat exchanger, represented in Fig. 1a, a refrigerant fluid flows in the finned tube in order to cool a mass of air flowing through a parallel plates arrangement. If one considers the temperature of the cold plate as constant, there is no need to simulate the refrigerant fluid domain. Due to the symmetry of the problem, the computational domain (Fig. 1b) corresponds to the lower half region between the cold plates. All prior twodimensional models set the west inlet air boundary in the origin of the plate (x ¼ Hx,1 in Fig. 2), whereas in this study it has been considered an extended domain in order to (1) avoid unrealistic physical boundary conditions, and to (2) permit horizontal frost growth in the leading edge. 2.2. Physical assumptions It is assumed that frost formation is symmetric with respect to the mid line between the parallel plates. In addition, it is assumed that the plate has null thickness in order to apply a symmetry condition in the south air boundary. The following physical assumptions are made for the air subdomain: (1) humid air is treated as an incompressible Newtonian fluid in laminar flow; (2) density (ra), mass diffusivity coefficient (Da) and specific heat capacity (cp,a) of air are assumed constant; and (3) inlet air velocity is large enough so that natural convection is negligible. The following physical assumptions are made for the frost subdomain: (4) the size of the frost porous are small enough that heat transfer by natural convection is negligible; (5) atmospheric pressure is considered as the total gas pressure of the humid air within the frost layer [6]; and (6) humid air within the frost layer is considered saturated [12]. The local humidity is calculated in the center of the control volumes, whereas frosteair interface is placed in the control
J.M. Armengol et al. / International Journal of Thermal Sciences 104 (2016) 245e256
247
Fig. 1. (a) Fin-and-tube heat exchanger, (b) physical model.
Fig. 2. Computational domain.
volume face, and therefore no assumption about supersaturated humidity in the airefrost interface is needed. Note that no assumption about one-dimensional growth rate is done, which means the frost grows in both directions: perpendicular and parallel to the cold plate. 2.3. Governing equation 2.3.1. Air subdomain Humid air flow is modeled using two-dimensional momentum, energy and mass transport equations, as well as mass continuity:
vra vðra uÞ vðra yÞ þ ¼ 0; þ vx vy vt
(1)
vðra uÞ vðra u$uÞ vðra y$uÞ vP v vu v vu þ þ ¼ þ m þ m ; vt vx vy vx vx vx vy vy
(2)
vðra yÞ vðra u$yÞ vðra y$yÞ vP v vy v vy þ þ ¼ þ m þ m ; vt vx vy vy vx vx vy vy
(3)
v rf T vt
v lf vT ¼ vx cp;f vx
vðra TÞ vðra u$TÞ vðra y$TÞ v la vT v la vT þ þ ¼ þ ; vt vx vy vx cp;a vx vy cp;a vy
(4)
vðra wÞ vðra u$wÞ vðra y$wÞ v vw v vw þ þ ¼ ra Da þ ra Da : vt vx vy vx vx vy vy
(5)
2.3.2. Frost subdomain By means of an energy balance for a differential frost volume in the interior of the frost layer, energy equation can be expressed as
v lf vT þ vy cp;f vy
! þ
qsub vrf : cp;f vt
(6)
For the densification rate of frost layer, the equation proposed by Na and Webb [8] has been used:
vrf v vw v vw ra Def þ ra Def : ¼ vx vx vy vy vt
(7)
2.3.3. Physical frost properties and transport properties of air As state in assumption (2) air density is set to a constant value (ra ¼ 1.252 kg m3), as well as the specific heat capacity (cp ¼ 1004.5 J kg1 Pa1). The viscosity and thermal conductivity are considered temperature dependent. However, those properties are assumed to be constant with respect to humidity since humidity ratio does not produce significant effects on the air mixture properties for the temperature range here studied [18]. In this study we used the Sutherland's law to compute temperature dependences:
mðTÞ ¼ m0
!
3=2 0 T T þS gR and lðTÞ ¼ mðTÞ ; T0 T þS ð1 gÞPr
(8)
where T0 ¼ 273 K, S ¼ 110.5 K, g ¼ 1.4, R ¼ 287 J kg1 K1, m0 ¼ 1.68 105 kg m1 s1 and Pr ¼ 0.71. Fig. 3 shows how air dynamic viscosity and thermal conductivity tends to increase as temperature increases following Sutherland's law. Most authors relate the properties of frost to its porosity and density. Nevertheless, these two variables are related through the expression 3
¼
ri rf : ri ra
(9)
248
J.M. Armengol et al. / International Journal of Thermal Sciences 104 (2016) 245e256
Fig. 3. Dynamic viscosity and thermal conductivity dependences with temperature.
The thermal conductivity of frost is commonly expressed as a function of frost density. In this study, it is used a correlation reported by Lee et al. [14]:
lf
h . i W ðm$KÞ1 ¼ A1 þ A2 rf þ A3 r2f
h . i r in kg m3 ;
4
(10)
7
where A1 ¼ 0.132, A2 ¼ 3.13 10 and A3 ¼ 1.6 10 . The specific heat is defined as a function of frost density and porosity, as follows:
cp;f
cp;g rg ð1 3 Þ þ cp;a ra 3 ¼ : rf
(11)
The diffusive mass coefficient in the frost layer is determined as Na and Webb [7]:
Def ¼ Da 3
1þ3 : 2
(12)
2.4. Boundary conditions The computational domain is subjected to the following boundary conditions: Bottom boundary (y ¼ 0) Symmetry condition (0 x < Hx,1)
y ¼ 0;
vu vT vw ¼ ¼ ¼ 0: vy vy vy
(13)
Plate (Hx,1 x Hx,2) Non-slip wall and isothermal conditions:
u ¼ y ¼ 0;
vw ¼ 0; vy
T ¼ Ts :
(14)
Left boundary (x ¼ 0; 0 y Hy) Inlet humid air:
u ¼ uin ;
y ¼ 0;
T ¼ Tin ;
w ¼ win :
u ¼ y ¼ 0;
vu vT vw ¼ ¼ ¼ 0: vy vy vy
(19)
2.5.2. Frost subdomain This study models the frost layer growth period while the effects of crystal growth period are treated as initial conditions. The required initial frost conditions are temperature, thickness and density of the frost layer. It is selected the same initial values of Lenic et al. [17] and Na and Webb [8] which came from experimental results performed by Jones and Parker [19]. These values are
. yf;0 ¼ 2 105 m; rf ;0 ¼ 30 kg m3 and Tf;0 ¼ Ts :
(20)
Initial frost layer temperature is assumed to be constant and equal to the plate temperature, since the initial thickness is sufficiently thin. 2.6. Frost layer growth rate Some previously published models [14,17] determine frost growth as a one-dimensional rate. Although the frost formation phenomenon occurs in every direction, frost layer could possibly grow in a unique direction. However, in order to model a twodimensional growth rate, a mass balance in a differential frost element placed in the frost layer interface is performed as follows:
dVf ¼ dt
Z
_ mdA;
(21)
SC
where the water vapor flux in the n direction is expressed as
m_ ¼ ra D (16)
(17)
Top boundary (y ¼ Hy; 0 x Hx,2) Symmetry condition:
y ¼ 0;
u0 ¼ y0 ¼ 0; wa;0 ¼ win and Ta;0 ¼ Tin :
(15)
Air subdomain (yf < y < Hy) Air outlet:
vu vy vT vw ¼ ¼ ¼ ¼ 0: vx vx vx vx
2.5.1. Air subdomain At the initial time, air subdomain has null velocities, while temperature and humidity are set to be the values at the inlet:
rf
Right boundary (x ¼ Hx,2): Frost subdomain (0 y yf)
vT vw ¼ ¼ 0: vx vx
2.5. Initial condition
(18)
dw ; dn
(22)
and the term in the left-hand side of Eq. (21), rf ðdVf =dtÞ, might be interpreted as the frost mass increase per time unit in the control volume. Note that dVf =dt stands for the time-varying volume of frost in the interface, while the volumes of every mesh cell remains constant since we are using a fixed grid approach. Once the mass balance is determined, and the frost mass rate in a particular interface frost element is known, one must decide the direction of the frost growth. Besides many option, keeping in mind the physical phenomenon, two main rules are assumed:
J.M. Armengol et al. / International Journal of Thermal Sciences 104 (2016) 245e256
249
Fig. 4. Mass balance in a frost element P at the airefrost interface.
Table 1 Effect of aspect ratio, using a time step of 1 s.
Tmax, C Tmean, C umax, gy/kga umean, gy/kga
Aspect ratio 6.9:1 (121 61)
Aspect ratio 4.6:1 (181 61)
Aspect ratio 3.4:1 (241 61)
19.767 8.369 8.469 6.017
19.762 8.257 8.464 5.987
19.764 8.264 8.465 5.988
1. Frost grows towards an air element. 2. In case the interface frost element has more than one air neighbor element, the growth is weighted as a function of the vapor mass coming from the air. As an example of the two-dimensional growth, let us consider a differential frost element P in the airefrost interface and its neighbors (E, W, N, S), as shown in Fig. 3. E and N neighbor nodes correspond to air elements, whereas P, W and S correspond to frost elements. Using rule 1, frost will grow in the north and east directions. Using rule 2, the distribution of the dVf =dt between N and E directions is expressed as:
dVf m_ a;e Ae dVf $ ¼P dt e m_ a A dt
(23a)
dVf m_ a;n An dVf $ ¼P dt n m_ a A dt
(23b)
dVf ¼ 0; dt s
(23c)
dVf ¼0 dt w
where m_ a;n and m_ a;e represent the water mass flux coming from the north air element and east air element respectively. The expression
m_ diff shown in Fig. 4 corresponds to the water mass flux diffusing through the frost layer. Note that in the general case the following equations must be satisfied:
dVf dVf dV dV dV ¼ þ f þ f þ f ; dt dt e dt w dt n dt s
(24a)
A m_ A m_ m_ A m_ A Pa;e e þ Pa;w w þ Pa;n n þ Pa;s s ¼ 1: m_ a A m_ a A m_ a A m_ a A
(24b)
The novelty of the two-dimensional frost growth is that the frost in the P control volume, from Fig. 4, has the numerical capability to grow in the north and the east directions. Moreover, the rate of growth does not have to be equal between the two directions, instead the growth is weighted as function of vapor mass coming from the air as stated in Eq. (23). Density of the new frost is determined as [17]:
vrf ¼ 0; vn
(25)
where n represents the growth direction. 3. Methodology for the numerical solution In order to solve the governing equations, a self-written Fortran code has been implemented. The governing equations are discretized by the finite volume method [20,21]. The velocity components are calculated at a staggered grid, while the scalar variables are calculated at the main grid (not staggered). The coupling between mass and momentum conservation equations is carried out using the SIMPLE algorithm. A hybrid scheme approximates the convection terms, while a central difference scheme is used for the diffusive terms. If the values in the mass balance for every control
Fig. 5. Measuring points.
250
J.M. Armengol et al. / International Journal of Thermal Sciences 104 (2016) 245e256
Table 2 Summary of frost conditions.
Run #1 Run #2 Run #3
uin, m/s
Tin, C
uin, gy/kga
Ts, C
0.6 0.6 1.7
21.4 19.8 21.5
6.2 8.2 9.6
19.5 20.5 13.7
volume, as well as the residual values of the different equations, are sufficiently low, an overall convergence is achieved (typically 1010). A fixed grid formulation is used to deal with the airefrost moving boundary [22]. In addition, a block-off in the frost region by means of a large value of viscosity (m ¼ 1030 Pa s) has been performed in order to impose null velocity. To determine thermophysical properties (r, l, m, cp and D) in the control volume faces, a harmonic mean is used [20]. This approach permits a proper handling of the abrupt changes of properties that may occur in the domain. The accuracy of the numerical results was verified through numerous tests based on the grid size and the time step effects. A mesh independence study was carried out with numerical meshes: 61 31, 81 41, 101 51, 121 61 and 141 71 showed no significant changes in the results comparing 121 61 and 141 71
meshes. In addition, we compared three different values of aspect ratio (6.9:1, 4.6:1 and 3.4:1) in Table 1. This table shows the maximum temperature (Tmax) and humidity ratio (umax) in x ¼ 130 mm (see Fig. 5), as well as the temperature (Tmean) and humidity ratio (umean) average along this line, for each aspect ratio after 10 min of simulation. Slightly changes are observed comparing the results of 4.6:1 and 3.4:1 aspect ratios. Additionally, the effect of the time step value was studied by using different time steps: 10 s, 5 s, 1 s and 0.5 s. Results do not vary significantly between 1 s and 0.5 s. Therefore, a 181 161 computational mesh and a time step of 1 s are used to perform the validation of the mathematical model. These tests were carried out with Tin ¼ 19.8 C, uin ¼ 8.5 gy/kga, uin ¼ 0.6 m/s and Ts ¼ 20.5 C. 3.1. Numerical procedure The algorithm used to predict the frost layer behavior can be summarized as follows: 1. Given the frost geometry and velocity field from the previous time step, velocity field is determined iteratively by the use of SIMPLE algorithm. 2. Given the frost density field, local frost properties are determined.
Fig. 6. Comparison of frost thickness between the present model and experimental data [23].
Fig. 7. Comparison of frost thickness for three different initial frost conditions. Working conditions of Run #3.
J.M. Armengol et al. / International Journal of Thermal Sciences 104 (2016) 245e256
3. Temperature distribution in the whole domain is numerically calculated. 4. Given the frost temperature distribution, the water vapor concentration within the frost layer is determined. 5. Humidity distribution in the whole domain is numerically calculated. 6. Densification of the frost layer is determined. 7. Steps 2 to 6 are repeated until energy and water vapor transport equations reach convergence. 8. Frost growth for the next step is determined 9. New air-frost interface is tracked. 10. Steps 1 to 9 are repeated until the desired time is achieved. 4. Results and discussion The model presented here is validated against experimental data reported by Leni c [23]. These data are suitable to validate a two-dimensional model since frost growth and temperature over time are reported as local values. The specific points in which temperature data was acquired are represented in Fig. 5 (T1, T2, T3 and T4), as well as the plane where frost growth is reported (x ¼ 130 mm). According to our knowledge, there are not previous published data which report the information as local values except for the experimental data reported by Lenic [23]. Geometric dimensions of the computational domain for the validation are: Hx,1 ¼ 20 mm, Hx,2 ¼ 140 mm and Hy ¼ 10 mm. The different frost conditions considered for the validation are summarized in Table 2. Fig. 6 shows numerical frost thickness over time for three different frost conditions at position x ¼ 130 mm, experimental data are also plotted for the sake of comparison. The numerical results for the frost growth have a reasonable agreement with the experimental data. In Run #1, a maximum deviation of 10.0% is observed after 1 h of simulation, while at this same time, deviations of 8.7% and 23.1% are seen in Run #2 and Run #3 respectively. Regarding the frost thickness at the final time (t ¼ 6 h) an excellent agreement is obtained for all cases: 1.7%, 3.8% and 4.2% for Run #1, Run #2 and Run #3 respectively. Note that experimental data from Run # 2 at the third hour and Run # 3 at the second hour differ significantly with respect to the presented model. In fact, Run # 3
251
over-predicts the frost thickness during the entire simulation. This might not be caused by the mathematical model but by the initial frost density. Despite the fact that the model focuses on the second period of the frost growth, the problem here studied is transient and therefore the results depend on the initial conditions. These initial conditions have been based on experimental values. Nevertheless, the working conditions affect the initial frost density. In the extreme case of Run #3 a large inlet air velocity (1.7 m/s) is considered. In order to quantify the influence of the initial frost density on the frost growth, Fig. 7 shows the evolution in time of frost thickness for three different initial frost conditions (30, 35 and 40 kg/ m3), for the working conditions of Run #3. As expected, it can be observed that a large value of initial frost density delays the frost growth. For an initial frost density of 40 kg/m3 a relative error of 6.0% in the frost thickness in the first hour of simulation is observed when compared with the experimental value. On the other hand, for this same instant of time an error of 23.1% is observed when the initial frost density is set to 30 kg/m3. Despite initial frost density does not solely depend on inlet velocity, these results are coherent with published correlation [24,13] which states that, to a greater or lesser extend, higher the value of inlet air velocity, greater the value of initial frost density. Moreover, it is found that initial frost density barely affects the final frost thickness for the studied range (30e40 kg/m3), which is consistent with the studies of Jones and Parker [19]. Fig. 8 compares the experimental data of the local values of temperature with the numerical results as a function of time, for the working conditions of Run #1. The numerical temperature results show a satisfactory agreement with trends of analyzed temperatures values. Fig. 9 shows the temperature distribution in the computational domain for different instants of time, for Run #1 working conditions. In addition, this figure shows the frost geometry at each time, which is represented with a black line. As reported by Cui et al. [15] and Kim et al. [16], frost grows faster in the leading edge than in the rear region, in accordance with the so-called entrance effect. The leading edge has larges gradients of humidity and temperature compared with the gradient values at the rear region. This fact
Fig. 8. Comparison of time-wise temperature variation between present simulation and the experimental data reported in Lenic [23].
252
J.M. Armengol et al. / International Journal of Thermal Sciences 104 (2016) 245e256
Fig. 9. Temperature distribution and frost geometry for different instants of time.
causes a larger amount of vapor mass transferred to the frost layer in the leading edge. During the first 2 h of simulation, the mass transferred to the frost layer is mainly used to increase its thickness. However, when a certain value of thickness is reached, the
temperature of the frost surface increases which causes slower growth of the frost. Then, the main part of the transferred vapor mass is used for the densification of the frost layer. After 2 h of simulation the frost thickness become more uniform along the
J.M. Armengol et al. / International Journal of Thermal Sciences 104 (2016) 245e256
253
Fig. 10. Density distribution and frost geometry for different instants of time.
x-axis due to the temperature increase of the leading edge, this frost layer shape is reported in recent experiments of frost formation between parallel plates performed by Nascimento et al. [25]. In Fig. 10 one can see the density distribution in the computational domain for Run #1 working conditions at different instants of times. Again, frost geometry is represented with a black line. Note that air density is considered constant, and density distribution is only computed in the frost layer by means of Eq. (7). It can be observed that frost density is larger in the leading edge region. In the rear region of the layer, where the humidity gradient in the frost surface is smaller than in the front, the density has moderated values. Frost is typically denser nearer the cold plate and becomes less dense near the frost surface. This density pattern is also reported by Lee et al. [14]. When the frost layer grows, the frost surface is distanced from the cold plate, and therefore there is an increase of the thermal resistance. On the other hand, the densification of the frost layer leads to an increase of the thermal
conductivity. However, the thermal resistance due to the frost growth is more relevant than the increase of thermal conductivity due to the densification, and therefore the frost surface temperature increases with time. Note that frost forms in x values smaller than the distance Hx,1, which means that frost also grows horizontally in the leading edge of the plate. Fig. 11 shows the humidity ratio distribution in the computational domain for different instants of time for Run #1 working conditions. Note that mass vapor transport in the air subdomain is computed using Eq. (5), while humidity ratio in the frost layer is determined as a temperature function (physical assumption (6)). It can be observed that the gradient of humidity in the leading edge is considerable larger when compared with the humidity gradient in the rear region. As a result, more vapor mass is transferred to the frost layer in the leading edge. Finally, in Fig. 12 it is shown the distribution of thermal conductivity in both subdomains for different instants of time for
254
J.M. Armengol et al. / International Journal of Thermal Sciences 104 (2016) 245e256
Fig. 11. Humidity rate distribution and frost geometry for different instants of time.
Run #1 working conditions. Since thermal conductivity in the frost layer is a function of the frost density, by Eq. (10), the same pattern than Fig. 10 is observed. That is, larger values of thermal conductivity in the leading edge and lower values in the
rear part of the layer. Additionally, thermal conductivity near the frost surface, where density is lower, present minor values when compared with frost near the plate where density is higher.
J.M. Armengol et al. / International Journal of Thermal Sciences 104 (2016) 245e256
255
Fig. 12. Thermal conductivity distribution and frost geometry for different instants of time.
5. Conclusions
Acknowledgment
A new model is formulated based on a two-dimensional frost growth rate, which allows frost elements not only increase its height but also its width. Therefore, this formulation can be applied to modeling frost formation not only over parallel plates but also can be extended on an arbitrary geometrical configuration. This two-dimensional approach leads to simulate a more uniform frost thickness along the plate which is found to be more physically realistic than some previously published two-dimensional models that assume one-dimensional frost growth rate. The model predictions of the frost thickness as a function of time agree with the experimental data within 10% of deviation for the case of intermediate plate temperature. This work provides evidence that an extended domain in the leading edge of the cold plate is required in order to study the frost formation in this region, since frost also grows horizontally in the leading edge.
The first author wishes to thank CAPES for financial support and UNICAMP for mobility grand. The second author wishes to thank CNPq for financial support and UNICAMP for mobility grand. References [1] Y. Hayashi, A. Aoki, S. Adachi, K. Hori, Study of frost properties correlating with frost formation types, J. Heat Transfer 99 (1977) 239e245. [2] K. Aoki, K. Katayama, Y. Hayashi, A study on frost formation: the process of frost formation involving the phenomena of water permeation and freezing, Bull. Jpn. Soc. Mech. Eng. 26 (1983) 87e93. [3] D. Oneal, D. Tree, Measurement of frost growth and density in a parallel plate geometry, ASHRAE Trans. 90 (1984) 278e290. [4] Y. Tao, R. Besant, K. Rezkallah, A mathematical model for predicting the densification and growth of frost on a flat plate, Int. J. Heat Mass Transfer 36 (1993) 353e363.
256
J.M. Armengol et al. / International Journal of Thermal Sciences 104 (2016) 245e256
[5] K. Ismail, C. Salinas, Modeling of frost formation over parallel cold plates, Int. J. Refrigeration 22 (1999) 425e441. [6] R. Le Gall, J. Grillot, C. Jallut, Modelling of frost growth and densification, Int. J. Heat Mass Transfer 40 (1997) 3177e3187. [7] B. Na, R.L. Webb, Mass transfer on and within a frost layer, Int. J. Heat Mass Transfer 47 (2004) 899e911. [8] B. Na, R.L. Webb, New model for frost growth rate, Int. J. Heat Mass Transfer 47 (2004) 925e936. [9] B. Na, R.L. Webb, A fundamental understanding of factors affecting frost nucleation, Int. J. Heat Mass Transfer 46 (2003) 3797e3808. [10] D.K. Yang, K.S. Lee, S. Song, Modeling for predicting frosting behavior of a finetube heat exchanger, Int. J. Heat Mass Transfer 49 (2006) 1472e1479. [11] C.J. Hermes, R.O. Piucco, J.R. Barbosa Jr., C. Melo, A study of frost growth and densification on flat surfaces, Exp. Therm. Fluid Sci. 33 (2009) 371e379. [12] M. Kandula, Frost growth and densification in laminar flow over flat surfaces, Int. J. Heat Mass Transfer 54 (2011) 3719e3731. [13] M. Kandula, Correlation of water frost porosity in laminar flow over flat surfaces, Special Top. Rev. Porous Media Int. J. 3 (2012) 79e87. [14] K.S. Lee, S. Jhee, D.-K. Yang, Prediction of the frost formation on a cold flat surface, Int. J. Heat Mass Transfer 46 (2003) 3789e3796. [15] J. Cui, W. Li, Y. Liu, Z. Jiang, A new time- and space-dependent model for predicting frost formation, Appl. Therm. Eng. 31 (2011) 447e457.
[16] D. Kim, C. Kim, K.-S. Lee, Frosting model for predicting macroscopic and local frost behaviors on a cold plate, Int. J. Heat Mass Transfer 82 (2015) 135e142. [17] K. Lenic, A. Trp, B. Frankovic, Transient two-dimensional model of frost formation on a fin-and-tube heat exchanger, Int. J. Heat Mass Transfer 52 (2009) 22e32. [18] P. Tsilingiris, Thermophysical and transport properties of humid air at temperature range between 0 and 100 C, Energy Convers. Manag. 49 (2008) 1098e1110. [19] B. Jones, J. Parker, Frost formation with varying environmental parameters, J. Heat Transfer 97 (1975) 255e259. [20] S. Patankar, Numerical Heat Transfer and Fluid Flow, CRC Press, 1980. [21] H.K. Versteeg, W. Malalasekera, An Introduction to Computational Fluid Dynamics: the Finite Volume Method, Pearson Education, 2007. [22] W. Shyy, H. Udaykumar, M.M. Rao, R.W. Smith, Computational Fluid Dynamics with Moving Boundaries, Courier Corporation, 2012. [23] K. Leni c, Analysis of heat and mass transfer during frost formation on fin-andtube heat exchangers, Eng. Rev. 26 (2006) 94. [24] Y. Mao, R. Besant, K. Rezkallah, Measurement and correlations of frost properties with airflow over a flat plate, ASHRAE Trans. 98 (2) (1992) 65e78. [25] V.S. Nascimento, F.R. Loyola, C.J. Hermes, A study of frost build-up on parallel plate channels, Exp. Therm. Fluid Sci. 60 (2015) 328e336.
Contents lists available at ScienceDirect
International Journal of Thermal Sciences journal homepage: www.elsevier.com/locate/ijts
Modeling of frost formation over parallel cold plates considering a two-dimensional growth rate n b, K.A.R. Ismail a J.M. Armengol a, C.T. Salinas a, *, J. Xama a b
Faculty of Mechanical Engineering, University of Campinas (UNICAMP), CP 6066, CEP 13081-970 Campinas, SP, Brazil n y Desarrollo Tecnolo gico, CENIDET-TNM-SEP Prol, Av. Palmira S/N. Col. Palmira, Cuernavaca, Morelos CP 62490, Mexico Centro Nacional de Investigacio
a r t i c l e i n f o
a b s t r a c t
Article history: Received 14 May 2015 Received in revised form 7 January 2016 Accepted 15 January 2016 Available online 19 April 2016
In this study, a frost formation model is presented based on a new two-dimensional approach for the growth rate. For modeling the frost formation over parallel cold plates, the basic transport equations of mass, energy and momentum have been discretized using the finite volume method in a twodimensional domain in which air and frost are considered. A fixed grid formulation is used to deal with the airefrost moving boundary. An extended domain in the inlet boundary has been considered in order to study the frost formation in the leading edge of the plate. The numerical results have been validated against experimental data in which frost growth and temperature as a function of time are reported as local values. The model predictions of the frost thickness as a function of time agree with the experimental data within 10% of deviation for the case of intermediate plate temperature. © 2016 Elsevier Masson SAS. All rights reserved.
Keywords: Frost Parallel plates Fin-and-tube heat exchanger
1. Introduction Frost formation consists in a phenomenon in which water vapor from moist air flow is deposited in a cold surface whose temperature is lower than the freezing temperature of water. Frost is a porous structure which consists of ice containing air gaps, so it presents a considerable thermal resistance. In addition, frost formation in the heat exchanger channels increases the pressure drop by flow restriction. Therefore, heat exchangers thermal performance decreases when frost formation appears. Many researchers have studied this phenomenon. One of the first among them was Hayashi et al. [1], who divided the frost growth process in three different periods. The first one is associated with the crystal growth, which is relatively short when compared with the whole process, the frost layer does not experience a significant thickness increase in this period. During this early stage, the frost cannot be considered a porous medium since the structure might be understood as ice columns where convective heat and mass transfer are the main mechanisms for frost growth, rather than diffusion within the frost layer.
* Corresponding author. E-mail address: [email protected] (C.T. Salinas). http://dx.doi.org/10.1016/j.ijthermalsci.2016.01.017 1290-0729/© 2016 Elsevier Masson SAS. All rights reserved.
In the second stage, known as the frost layer growth period, frost is treated as a porous medium in which molecular diffusion of water vapor occurs. The total mass flux of water vapor coming from the air contributes in two different processes: the densification of the frost layer and the frost growth. Finally, in the frost layer full growth period, the frost surface temperature reaches the triple point of water, and so it begins a cyclic process in which the condensed vapor at the frost surface permeates into the frost layer and gradually freezes because of the inner temperature gradient. A study of this stage is presented in Aoki et al. [2]. In the last decades, many models for frost growth on flat plates have been proposed. A review of frost formation in simple geom etries is found in Oneal and Tree [3]. One can classify the theoretical models into two groups, depending on the use or not of empirical correlation for the airefrost interaction. In the first group, the models use empirical correlations for the air side effect and they apply energy and mass transfer equations within the frost layer to predict frost formation. Tao et al. [4] developed a mathematical model based on the frost structure analysis done by Hayashi et al. [1] using the local volume averaging technique to deal with the frost as a porous media. Based on that model, Ismail and Salinas [5] determined the best initial values of the diffusivity, initial radius and geometry of ice crystals. Using the local volume averaging technique, Le Gall et al. [6] proposed an interface balance to predict frost growth rate. In
246
J.M. Armengol et al. / International Journal of Thermal Sciences 104 (2016) 245e256
Nomenclature m_ water vapor mass flux, kg m2 s1 A area, m2 cp specific heat capacity, J kg1 K1 D mass diffusivity coefficient, m2 s1 E, W, N, S respective east, west, north and south computational nodes e, w, n, s respective east, west, north and south control volume faces H distance, m P pressure, Pa qsub specific heat of sublimation, J kg1 T temperature, C t time, s u horizontal velocity, m s1 V volume, m3 v vertical velocity, m s1 x coordinate, m y coordinate, m
their work, Na and Webb [7e9] proposed a model considering supersaturated water vapor density at the frost surface. Yang et al. [10] reported an effective frosting model for predicting the thermal performance of a finetube heat exchanger. More recently, Hermes et al. [11] developed a mathematical model whose predictions of the frost thickness agreed within 10% error bands when compared with experimental data. Kandula [12,13] proposed new frost correlations for the flat plate and studied the effects of the environment parameters in the frost characteristics. Even though one-dimensional models based on empirical correlations can predict frost formation with reasonable precision, they require experimental determination of the empirical correlations for every different geometry. The second group, consists in an attempt to avoid empirical correlations by simulating the air flow and the frost formation simultaneously. As frost formation is a common phenomenon which takes place in different surface geometries, it is required a better understanding of the two-dimensional frost growth without using empirical correlations. Lee et al. [14] proposed a twodimensional model that shows good agreement with experimental data for different working conditions. However, the model assumes that the frost density in the direction normal to the cooling plate is constant. Recently, Cui et al. [15] and Kim et al. [16] simulated frost formation with a commercial Computational Fluid Dynamic (CFD) package getting results with good agreement with experimental data. These models can determine local properties of the frost layer such as density, humidity and temperature. In these models, it was found that frost layer at the leading edge of the plate presents larger values of density and thickness when compared with the values at the trailing edge. Lenic et al. [17] proposed a twodimensional mathematical model based on the assumptions of Na and Webb [7] and modeled it using the finite volume technique. The model of Lenic is similar to the one proposed in this work, except for the treatment of the frost growth rate. The main goal of this investigation is to establish a mathematical model for the frost growth in two dimensions, which means that a frost element not only increases its height but also its width. Such approach is developed in the context of the second group, so it has been modeled the air and the frost subdomains.
Greek symbols l thermal conductivity, W m1 K1 m dynamic viscosity, Pa s u humidity ratio, kgv kg1 a r density, kg m3 3 porosity Subscripts 0 initial value diff diffusing through the frost layer a air CS control surface ef effective f frost i ice in inlet s plate surface v vapor
2. Mathematical and physical model 2.1. Physical model of the fin-and-tube heat exchanger The frost formation is here modeled in the region between parallel cold plates. Such configuration represents the air channels of a fin-and-tube heat exchanger, which is very common in refrigeration systems. In the heat exchanger, represented in Fig. 1a, a refrigerant fluid flows in the finned tube in order to cool a mass of air flowing through a parallel plates arrangement. If one considers the temperature of the cold plate as constant, there is no need to simulate the refrigerant fluid domain. Due to the symmetry of the problem, the computational domain (Fig. 1b) corresponds to the lower half region between the cold plates. All prior twodimensional models set the west inlet air boundary in the origin of the plate (x ¼ Hx,1 in Fig. 2), whereas in this study it has been considered an extended domain in order to (1) avoid unrealistic physical boundary conditions, and to (2) permit horizontal frost growth in the leading edge. 2.2. Physical assumptions It is assumed that frost formation is symmetric with respect to the mid line between the parallel plates. In addition, it is assumed that the plate has null thickness in order to apply a symmetry condition in the south air boundary. The following physical assumptions are made for the air subdomain: (1) humid air is treated as an incompressible Newtonian fluid in laminar flow; (2) density (ra), mass diffusivity coefficient (Da) and specific heat capacity (cp,a) of air are assumed constant; and (3) inlet air velocity is large enough so that natural convection is negligible. The following physical assumptions are made for the frost subdomain: (4) the size of the frost porous are small enough that heat transfer by natural convection is negligible; (5) atmospheric pressure is considered as the total gas pressure of the humid air within the frost layer [6]; and (6) humid air within the frost layer is considered saturated [12]. The local humidity is calculated in the center of the control volumes, whereas frosteair interface is placed in the control
J.M. Armengol et al. / International Journal of Thermal Sciences 104 (2016) 245e256
247
Fig. 1. (a) Fin-and-tube heat exchanger, (b) physical model.
Fig. 2. Computational domain.
volume face, and therefore no assumption about supersaturated humidity in the airefrost interface is needed. Note that no assumption about one-dimensional growth rate is done, which means the frost grows in both directions: perpendicular and parallel to the cold plate. 2.3. Governing equation 2.3.1. Air subdomain Humid air flow is modeled using two-dimensional momentum, energy and mass transport equations, as well as mass continuity:
vra vðra uÞ vðra yÞ þ ¼ 0; þ vx vy vt
(1)
vðra uÞ vðra u$uÞ vðra y$uÞ vP v vu v vu þ þ ¼ þ m þ m ; vt vx vy vx vx vx vy vy
(2)
vðra yÞ vðra u$yÞ vðra y$yÞ vP v vy v vy þ þ ¼ þ m þ m ; vt vx vy vy vx vx vy vy
(3)
v rf T vt
v lf vT ¼ vx cp;f vx
vðra TÞ vðra u$TÞ vðra y$TÞ v la vT v la vT þ þ ¼ þ ; vt vx vy vx cp;a vx vy cp;a vy
(4)
vðra wÞ vðra u$wÞ vðra y$wÞ v vw v vw þ þ ¼ ra Da þ ra Da : vt vx vy vx vx vy vy
(5)
2.3.2. Frost subdomain By means of an energy balance for a differential frost volume in the interior of the frost layer, energy equation can be expressed as
v lf vT þ vy cp;f vy
! þ
qsub vrf : cp;f vt
(6)
For the densification rate of frost layer, the equation proposed by Na and Webb [8] has been used:
vrf v vw v vw ra Def þ ra Def : ¼ vx vx vy vy vt
(7)
2.3.3. Physical frost properties and transport properties of air As state in assumption (2) air density is set to a constant value (ra ¼ 1.252 kg m3), as well as the specific heat capacity (cp ¼ 1004.5 J kg1 Pa1). The viscosity and thermal conductivity are considered temperature dependent. However, those properties are assumed to be constant with respect to humidity since humidity ratio does not produce significant effects on the air mixture properties for the temperature range here studied [18]. In this study we used the Sutherland's law to compute temperature dependences:
mðTÞ ¼ m0
!
3=2 0 T T þS gR and lðTÞ ¼ mðTÞ ; T0 T þS ð1 gÞPr
(8)
where T0 ¼ 273 K, S ¼ 110.5 K, g ¼ 1.4, R ¼ 287 J kg1 K1, m0 ¼ 1.68 105 kg m1 s1 and Pr ¼ 0.71. Fig. 3 shows how air dynamic viscosity and thermal conductivity tends to increase as temperature increases following Sutherland's law. Most authors relate the properties of frost to its porosity and density. Nevertheless, these two variables are related through the expression 3
¼
ri rf : ri ra
(9)
248
J.M. Armengol et al. / International Journal of Thermal Sciences 104 (2016) 245e256
Fig. 3. Dynamic viscosity and thermal conductivity dependences with temperature.
The thermal conductivity of frost is commonly expressed as a function of frost density. In this study, it is used a correlation reported by Lee et al. [14]:
lf
h . i W ðm$KÞ1 ¼ A1 þ A2 rf þ A3 r2f
h . i r in kg m3 ;
4
(10)
7
where A1 ¼ 0.132, A2 ¼ 3.13 10 and A3 ¼ 1.6 10 . The specific heat is defined as a function of frost density and porosity, as follows:
cp;f
cp;g rg ð1 3 Þ þ cp;a ra 3 ¼ : rf
(11)
The diffusive mass coefficient in the frost layer is determined as Na and Webb [7]:
Def ¼ Da 3
1þ3 : 2
(12)
2.4. Boundary conditions The computational domain is subjected to the following boundary conditions: Bottom boundary (y ¼ 0) Symmetry condition (0 x < Hx,1)
y ¼ 0;
vu vT vw ¼ ¼ ¼ 0: vy vy vy
(13)
Plate (Hx,1 x Hx,2) Non-slip wall and isothermal conditions:
u ¼ y ¼ 0;
vw ¼ 0; vy
T ¼ Ts :
(14)
Left boundary (x ¼ 0; 0 y Hy) Inlet humid air:
u ¼ uin ;
y ¼ 0;
T ¼ Tin ;
w ¼ win :
u ¼ y ¼ 0;
vu vT vw ¼ ¼ ¼ 0: vy vy vy
(19)
2.5.2. Frost subdomain This study models the frost layer growth period while the effects of crystal growth period are treated as initial conditions. The required initial frost conditions are temperature, thickness and density of the frost layer. It is selected the same initial values of Lenic et al. [17] and Na and Webb [8] which came from experimental results performed by Jones and Parker [19]. These values are
. yf;0 ¼ 2 105 m; rf ;0 ¼ 30 kg m3 and Tf;0 ¼ Ts :
(20)
Initial frost layer temperature is assumed to be constant and equal to the plate temperature, since the initial thickness is sufficiently thin. 2.6. Frost layer growth rate Some previously published models [14,17] determine frost growth as a one-dimensional rate. Although the frost formation phenomenon occurs in every direction, frost layer could possibly grow in a unique direction. However, in order to model a twodimensional growth rate, a mass balance in a differential frost element placed in the frost layer interface is performed as follows:
dVf ¼ dt
Z
_ mdA;
(21)
SC
where the water vapor flux in the n direction is expressed as
m_ ¼ ra D (16)
(17)
Top boundary (y ¼ Hy; 0 x Hx,2) Symmetry condition:
y ¼ 0;
u0 ¼ y0 ¼ 0; wa;0 ¼ win and Ta;0 ¼ Tin :
(15)
Air subdomain (yf < y < Hy) Air outlet:
vu vy vT vw ¼ ¼ ¼ ¼ 0: vx vx vx vx
2.5.1. Air subdomain At the initial time, air subdomain has null velocities, while temperature and humidity are set to be the values at the inlet:
rf
Right boundary (x ¼ Hx,2): Frost subdomain (0 y yf)
vT vw ¼ ¼ 0: vx vx
2.5. Initial condition
(18)
dw ; dn
(22)
and the term in the left-hand side of Eq. (21), rf ðdVf =dtÞ, might be interpreted as the frost mass increase per time unit in the control volume. Note that dVf =dt stands for the time-varying volume of frost in the interface, while the volumes of every mesh cell remains constant since we are using a fixed grid approach. Once the mass balance is determined, and the frost mass rate in a particular interface frost element is known, one must decide the direction of the frost growth. Besides many option, keeping in mind the physical phenomenon, two main rules are assumed:
J.M. Armengol et al. / International Journal of Thermal Sciences 104 (2016) 245e256
249
Fig. 4. Mass balance in a frost element P at the airefrost interface.
Table 1 Effect of aspect ratio, using a time step of 1 s.
Tmax, C Tmean, C umax, gy/kga umean, gy/kga
Aspect ratio 6.9:1 (121 61)
Aspect ratio 4.6:1 (181 61)
Aspect ratio 3.4:1 (241 61)
19.767 8.369 8.469 6.017
19.762 8.257 8.464 5.987
19.764 8.264 8.465 5.988
1. Frost grows towards an air element. 2. In case the interface frost element has more than one air neighbor element, the growth is weighted as a function of the vapor mass coming from the air. As an example of the two-dimensional growth, let us consider a differential frost element P in the airefrost interface and its neighbors (E, W, N, S), as shown in Fig. 3. E and N neighbor nodes correspond to air elements, whereas P, W and S correspond to frost elements. Using rule 1, frost will grow in the north and east directions. Using rule 2, the distribution of the dVf =dt between N and E directions is expressed as:
dVf m_ a;e Ae dVf $ ¼P dt e m_ a A dt
(23a)
dVf m_ a;n An dVf $ ¼P dt n m_ a A dt
(23b)
dVf ¼ 0; dt s
(23c)
dVf ¼0 dt w
where m_ a;n and m_ a;e represent the water mass flux coming from the north air element and east air element respectively. The expression
m_ diff shown in Fig. 4 corresponds to the water mass flux diffusing through the frost layer. Note that in the general case the following equations must be satisfied:
dVf dVf dV dV dV ¼ þ f þ f þ f ; dt dt e dt w dt n dt s
(24a)
A m_ A m_ m_ A m_ A Pa;e e þ Pa;w w þ Pa;n n þ Pa;s s ¼ 1: m_ a A m_ a A m_ a A m_ a A
(24b)
The novelty of the two-dimensional frost growth is that the frost in the P control volume, from Fig. 4, has the numerical capability to grow in the north and the east directions. Moreover, the rate of growth does not have to be equal between the two directions, instead the growth is weighted as function of vapor mass coming from the air as stated in Eq. (23). Density of the new frost is determined as [17]:
vrf ¼ 0; vn
(25)
where n represents the growth direction. 3. Methodology for the numerical solution In order to solve the governing equations, a self-written Fortran code has been implemented. The governing equations are discretized by the finite volume method [20,21]. The velocity components are calculated at a staggered grid, while the scalar variables are calculated at the main grid (not staggered). The coupling between mass and momentum conservation equations is carried out using the SIMPLE algorithm. A hybrid scheme approximates the convection terms, while a central difference scheme is used for the diffusive terms. If the values in the mass balance for every control
Fig. 5. Measuring points.
250
J.M. Armengol et al. / International Journal of Thermal Sciences 104 (2016) 245e256
Table 2 Summary of frost conditions.
Run #1 Run #2 Run #3
uin, m/s
Tin, C
uin, gy/kga
Ts, C
0.6 0.6 1.7
21.4 19.8 21.5
6.2 8.2 9.6
19.5 20.5 13.7
volume, as well as the residual values of the different equations, are sufficiently low, an overall convergence is achieved (typically 1010). A fixed grid formulation is used to deal with the airefrost moving boundary [22]. In addition, a block-off in the frost region by means of a large value of viscosity (m ¼ 1030 Pa s) has been performed in order to impose null velocity. To determine thermophysical properties (r, l, m, cp and D) in the control volume faces, a harmonic mean is used [20]. This approach permits a proper handling of the abrupt changes of properties that may occur in the domain. The accuracy of the numerical results was verified through numerous tests based on the grid size and the time step effects. A mesh independence study was carried out with numerical meshes: 61 31, 81 41, 101 51, 121 61 and 141 71 showed no significant changes in the results comparing 121 61 and 141 71
meshes. In addition, we compared three different values of aspect ratio (6.9:1, 4.6:1 and 3.4:1) in Table 1. This table shows the maximum temperature (Tmax) and humidity ratio (umax) in x ¼ 130 mm (see Fig. 5), as well as the temperature (Tmean) and humidity ratio (umean) average along this line, for each aspect ratio after 10 min of simulation. Slightly changes are observed comparing the results of 4.6:1 and 3.4:1 aspect ratios. Additionally, the effect of the time step value was studied by using different time steps: 10 s, 5 s, 1 s and 0.5 s. Results do not vary significantly between 1 s and 0.5 s. Therefore, a 181 161 computational mesh and a time step of 1 s are used to perform the validation of the mathematical model. These tests were carried out with Tin ¼ 19.8 C, uin ¼ 8.5 gy/kga, uin ¼ 0.6 m/s and Ts ¼ 20.5 C. 3.1. Numerical procedure The algorithm used to predict the frost layer behavior can be summarized as follows: 1. Given the frost geometry and velocity field from the previous time step, velocity field is determined iteratively by the use of SIMPLE algorithm. 2. Given the frost density field, local frost properties are determined.
Fig. 6. Comparison of frost thickness between the present model and experimental data [23].
Fig. 7. Comparison of frost thickness for three different initial frost conditions. Working conditions of Run #3.
J.M. Armengol et al. / International Journal of Thermal Sciences 104 (2016) 245e256
3. Temperature distribution in the whole domain is numerically calculated. 4. Given the frost temperature distribution, the water vapor concentration within the frost layer is determined. 5. Humidity distribution in the whole domain is numerically calculated. 6. Densification of the frost layer is determined. 7. Steps 2 to 6 are repeated until energy and water vapor transport equations reach convergence. 8. Frost growth for the next step is determined 9. New air-frost interface is tracked. 10. Steps 1 to 9 are repeated until the desired time is achieved. 4. Results and discussion The model presented here is validated against experimental data reported by Leni c [23]. These data are suitable to validate a two-dimensional model since frost growth and temperature over time are reported as local values. The specific points in which temperature data was acquired are represented in Fig. 5 (T1, T2, T3 and T4), as well as the plane where frost growth is reported (x ¼ 130 mm). According to our knowledge, there are not previous published data which report the information as local values except for the experimental data reported by Lenic [23]. Geometric dimensions of the computational domain for the validation are: Hx,1 ¼ 20 mm, Hx,2 ¼ 140 mm and Hy ¼ 10 mm. The different frost conditions considered for the validation are summarized in Table 2. Fig. 6 shows numerical frost thickness over time for three different frost conditions at position x ¼ 130 mm, experimental data are also plotted for the sake of comparison. The numerical results for the frost growth have a reasonable agreement with the experimental data. In Run #1, a maximum deviation of 10.0% is observed after 1 h of simulation, while at this same time, deviations of 8.7% and 23.1% are seen in Run #2 and Run #3 respectively. Regarding the frost thickness at the final time (t ¼ 6 h) an excellent agreement is obtained for all cases: 1.7%, 3.8% and 4.2% for Run #1, Run #2 and Run #3 respectively. Note that experimental data from Run # 2 at the third hour and Run # 3 at the second hour differ significantly with respect to the presented model. In fact, Run # 3
251
over-predicts the frost thickness during the entire simulation. This might not be caused by the mathematical model but by the initial frost density. Despite the fact that the model focuses on the second period of the frost growth, the problem here studied is transient and therefore the results depend on the initial conditions. These initial conditions have been based on experimental values. Nevertheless, the working conditions affect the initial frost density. In the extreme case of Run #3 a large inlet air velocity (1.7 m/s) is considered. In order to quantify the influence of the initial frost density on the frost growth, Fig. 7 shows the evolution in time of frost thickness for three different initial frost conditions (30, 35 and 40 kg/ m3), for the working conditions of Run #3. As expected, it can be observed that a large value of initial frost density delays the frost growth. For an initial frost density of 40 kg/m3 a relative error of 6.0% in the frost thickness in the first hour of simulation is observed when compared with the experimental value. On the other hand, for this same instant of time an error of 23.1% is observed when the initial frost density is set to 30 kg/m3. Despite initial frost density does not solely depend on inlet velocity, these results are coherent with published correlation [24,13] which states that, to a greater or lesser extend, higher the value of inlet air velocity, greater the value of initial frost density. Moreover, it is found that initial frost density barely affects the final frost thickness for the studied range (30e40 kg/m3), which is consistent with the studies of Jones and Parker [19]. Fig. 8 compares the experimental data of the local values of temperature with the numerical results as a function of time, for the working conditions of Run #1. The numerical temperature results show a satisfactory agreement with trends of analyzed temperatures values. Fig. 9 shows the temperature distribution in the computational domain for different instants of time, for Run #1 working conditions. In addition, this figure shows the frost geometry at each time, which is represented with a black line. As reported by Cui et al. [15] and Kim et al. [16], frost grows faster in the leading edge than in the rear region, in accordance with the so-called entrance effect. The leading edge has larges gradients of humidity and temperature compared with the gradient values at the rear region. This fact
Fig. 8. Comparison of time-wise temperature variation between present simulation and the experimental data reported in Lenic [23].
252
J.M. Armengol et al. / International Journal of Thermal Sciences 104 (2016) 245e256
Fig. 9. Temperature distribution and frost geometry for different instants of time.
causes a larger amount of vapor mass transferred to the frost layer in the leading edge. During the first 2 h of simulation, the mass transferred to the frost layer is mainly used to increase its thickness. However, when a certain value of thickness is reached, the
temperature of the frost surface increases which causes slower growth of the frost. Then, the main part of the transferred vapor mass is used for the densification of the frost layer. After 2 h of simulation the frost thickness become more uniform along the
J.M. Armengol et al. / International Journal of Thermal Sciences 104 (2016) 245e256
253
Fig. 10. Density distribution and frost geometry for different instants of time.
x-axis due to the temperature increase of the leading edge, this frost layer shape is reported in recent experiments of frost formation between parallel plates performed by Nascimento et al. [25]. In Fig. 10 one can see the density distribution in the computational domain for Run #1 working conditions at different instants of times. Again, frost geometry is represented with a black line. Note that air density is considered constant, and density distribution is only computed in the frost layer by means of Eq. (7). It can be observed that frost density is larger in the leading edge region. In the rear region of the layer, where the humidity gradient in the frost surface is smaller than in the front, the density has moderated values. Frost is typically denser nearer the cold plate and becomes less dense near the frost surface. This density pattern is also reported by Lee et al. [14]. When the frost layer grows, the frost surface is distanced from the cold plate, and therefore there is an increase of the thermal resistance. On the other hand, the densification of the frost layer leads to an increase of the thermal
conductivity. However, the thermal resistance due to the frost growth is more relevant than the increase of thermal conductivity due to the densification, and therefore the frost surface temperature increases with time. Note that frost forms in x values smaller than the distance Hx,1, which means that frost also grows horizontally in the leading edge of the plate. Fig. 11 shows the humidity ratio distribution in the computational domain for different instants of time for Run #1 working conditions. Note that mass vapor transport in the air subdomain is computed using Eq. (5), while humidity ratio in the frost layer is determined as a temperature function (physical assumption (6)). It can be observed that the gradient of humidity in the leading edge is considerable larger when compared with the humidity gradient in the rear region. As a result, more vapor mass is transferred to the frost layer in the leading edge. Finally, in Fig. 12 it is shown the distribution of thermal conductivity in both subdomains for different instants of time for
254
J.M. Armengol et al. / International Journal of Thermal Sciences 104 (2016) 245e256
Fig. 11. Humidity rate distribution and frost geometry for different instants of time.
Run #1 working conditions. Since thermal conductivity in the frost layer is a function of the frost density, by Eq. (10), the same pattern than Fig. 10 is observed. That is, larger values of thermal conductivity in the leading edge and lower values in the
rear part of the layer. Additionally, thermal conductivity near the frost surface, where density is lower, present minor values when compared with frost near the plate where density is higher.
J.M. Armengol et al. / International Journal of Thermal Sciences 104 (2016) 245e256
255
Fig. 12. Thermal conductivity distribution and frost geometry for different instants of time.
5. Conclusions
Acknowledgment
A new model is formulated based on a two-dimensional frost growth rate, which allows frost elements not only increase its height but also its width. Therefore, this formulation can be applied to modeling frost formation not only over parallel plates but also can be extended on an arbitrary geometrical configuration. This two-dimensional approach leads to simulate a more uniform frost thickness along the plate which is found to be more physically realistic than some previously published two-dimensional models that assume one-dimensional frost growth rate. The model predictions of the frost thickness as a function of time agree with the experimental data within 10% of deviation for the case of intermediate plate temperature. This work provides evidence that an extended domain in the leading edge of the cold plate is required in order to study the frost formation in this region, since frost also grows horizontally in the leading edge.
The first author wishes to thank CAPES for financial support and UNICAMP for mobility grand. The second author wishes to thank CNPq for financial support and UNICAMP for mobility grand. References [1] Y. Hayashi, A. Aoki, S. Adachi, K. Hori, Study of frost properties correlating with frost formation types, J. Heat Transfer 99 (1977) 239e245. [2] K. Aoki, K. Katayama, Y. Hayashi, A study on frost formation: the process of frost formation involving the phenomena of water permeation and freezing, Bull. Jpn. Soc. Mech. Eng. 26 (1983) 87e93. [3] D. Oneal, D. Tree, Measurement of frost growth and density in a parallel plate geometry, ASHRAE Trans. 90 (1984) 278e290. [4] Y. Tao, R. Besant, K. Rezkallah, A mathematical model for predicting the densification and growth of frost on a flat plate, Int. J. Heat Mass Transfer 36 (1993) 353e363.
256
J.M. Armengol et al. / International Journal of Thermal Sciences 104 (2016) 245e256
[5] K. Ismail, C. Salinas, Modeling of frost formation over parallel cold plates, Int. J. Refrigeration 22 (1999) 425e441. [6] R. Le Gall, J. Grillot, C. Jallut, Modelling of frost growth and densification, Int. J. Heat Mass Transfer 40 (1997) 3177e3187. [7] B. Na, R.L. Webb, Mass transfer on and within a frost layer, Int. J. Heat Mass Transfer 47 (2004) 899e911. [8] B. Na, R.L. Webb, New model for frost growth rate, Int. J. Heat Mass Transfer 47 (2004) 925e936. [9] B. Na, R.L. Webb, A fundamental understanding of factors affecting frost nucleation, Int. J. Heat Mass Transfer 46 (2003) 3797e3808. [10] D.K. Yang, K.S. Lee, S. Song, Modeling for predicting frosting behavior of a finetube heat exchanger, Int. J. Heat Mass Transfer 49 (2006) 1472e1479. [11] C.J. Hermes, R.O. Piucco, J.R. Barbosa Jr., C. Melo, A study of frost growth and densification on flat surfaces, Exp. Therm. Fluid Sci. 33 (2009) 371e379. [12] M. Kandula, Frost growth and densification in laminar flow over flat surfaces, Int. J. Heat Mass Transfer 54 (2011) 3719e3731. [13] M. Kandula, Correlation of water frost porosity in laminar flow over flat surfaces, Special Top. Rev. Porous Media Int. J. 3 (2012) 79e87. [14] K.S. Lee, S. Jhee, D.-K. Yang, Prediction of the frost formation on a cold flat surface, Int. J. Heat Mass Transfer 46 (2003) 3789e3796. [15] J. Cui, W. Li, Y. Liu, Z. Jiang, A new time- and space-dependent model for predicting frost formation, Appl. Therm. Eng. 31 (2011) 447e457.
[16] D. Kim, C. Kim, K.-S. Lee, Frosting model for predicting macroscopic and local frost behaviors on a cold plate, Int. J. Heat Mass Transfer 82 (2015) 135e142. [17] K. Lenic, A. Trp, B. Frankovic, Transient two-dimensional model of frost formation on a fin-and-tube heat exchanger, Int. J. Heat Mass Transfer 52 (2009) 22e32. [18] P. Tsilingiris, Thermophysical and transport properties of humid air at temperature range between 0 and 100 C, Energy Convers. Manag. 49 (2008) 1098e1110. [19] B. Jones, J. Parker, Frost formation with varying environmental parameters, J. Heat Transfer 97 (1975) 255e259. [20] S. Patankar, Numerical Heat Transfer and Fluid Flow, CRC Press, 1980. [21] H.K. Versteeg, W. Malalasekera, An Introduction to Computational Fluid Dynamics: the Finite Volume Method, Pearson Education, 2007. [22] W. Shyy, H. Udaykumar, M.M. Rao, R.W. Smith, Computational Fluid Dynamics with Moving Boundaries, Courier Corporation, 2012. [23] K. Leni c, Analysis of heat and mass transfer during frost formation on fin-andtube heat exchangers, Eng. Rev. 26 (2006) 94. [24] Y. Mao, R. Besant, K. Rezkallah, Measurement and correlations of frost properties with airflow over a flat plate, ASHRAE Trans. 98 (2) (1992) 65e78. [25] V.S. Nascimento, F.R. Loyola, C.J. Hermes, A study of frost build-up on parallel plate channels, Exp. Therm. Fluid Sci. 60 (2015) 328e336.