- Email: [email protected]
CFD modeling of a polymer solar collector
Renewable Energy 35 (2010) 1499–1508
Contents lists available at ScienceDirect
Renewable Energy journal homepage: www.elsevier.com/locate/renene
CFD modeling of a polymer solar collector G. Martinopoulos a, D. Missirlis b, G. Tsilingiridis a, K. Yakinthos b, N. Kyriakis a, * a b
Process Equipment Design Laboratory, Dept. of Mechanical Engineering, AUTh, POB 487, PC 54124, Thessaloniki, Greece Laboratory of Fluid Mechanics and Turbomachinery, Dept. of Mechanical Engineering, AUTh, PC 54124, Thessaloniki, Greece
a r t i c l e i n f o
a b s t r a c t
Article history: Received 20 May 2009 Accepted 5 January 2010 Available online 12 February 2010
A polymer solar collector is developed and its behavior is investigated both experimentally and with computational fluid dynamics (CFD). Solar irradiation as well as convection and heat transfer in the circulating fluid and between the parts of the collector is considered in the model. The temperature and velocity distribution over its area as well as the collector efficiency at nominal flow rate were used in order to validate the CFD model. Temperature distribution during operation and average collector efficiency were found to be in good agreement between the experimental data and the results of the CFD modeling. Ó 2010 Elsevier Ltd. All rights reserved.
Keywords: Solar energy Solar collector CFD simulation
1. Introduction As the oil market becomes unstable and unpredictable, the use of renewable energy sources, solar energy for thermal applications in particular, offers an interesting solution. Solar energy in the form of domestic solar thermal systems accounted for 127.8 GWth worldwide with more than 15 GWth in the EU alone, in 2006 [1,2]. The increase in the use of solar thermal energy systems is subjected to their economic viability, in comparison with conventional systems, and their simplicity in manufacturing and use. In order to minimize the production cost of current solar collectors, new materials of lower cost and preferably more environmental friendly must substitute the ones that are employed today, without any deterioration of collector efficiency. In this paper a novel polymer solar collector is presented and computational fluid dynamics are used in an effort to investigate its behavior. The performance of the collector under nominal flow rate was obtained by CFD. For the validation of the CFD model, the temperature and velocity distributions over its area were compared to experimental data. 2. Experimental model The most commonly used solar system is the one producing hot water for domestic use. The basic unit of such a system is the solar collector. Among the solar collector configurations used, the most common is the flat plate one [2]. Flat plate collectors have been
* Corresponding author. Tel.: þ30 2310 996083; fax: þ30 2310 996087. E-mail address: [email protected] (N. Kyriakis). 0960-1481/$ – see front matter Ó 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.renene.2010.01.004
studied by many researchers, e.g. Hottel and Woertz [3], Bliss [4], Nahar and Garg [5], Francken [6] and others. Typical flat plate collectors are either corrugated, bond duct or tube-in-plate type and their performance depends on various design parameters, including the number of covers, the type and thickness of glazing, the anti-reflecting coating on cover glass, the type of coating on the collector plate, the spacing between the collector and the inner glass, the type and thickness of insulation used, etc. Operation parameters affecting the performance of a solar collector are the mass flow rate of fluid, the amount of incident solar radiation, the inlet and ambient temperatures and sky conditions, to name a few. A number of researchers have adopted the use of polymers in solar collector design since their cost and physical properties make the volume production of lightweight low-cost and corrosion resistance collectors possible. The use of polymer glazing has been moderately successful, as they offer significant potential for cost savings both as direct substitutes for glass cover plates in traditional collector systems and as an integral part of all-polymeric systems. Polymer heat exchangers offer the potential advantages of reduced cost of materials and manufacture, resistance to corrosion, as well as better integration with other components [7]. Furthermore the use of polymer materials reduces the collector weight by 50% in comparison with a traditional metal collector allowing much easier installation [8]. Currently the use of polymers as solar collector materials is limited mainly in swimming pool collectors and in glazing of typical flat plate collectors [9]. In this work a novel polymer collector was developed aiming at producing a solar collector which combines similar or better technical characteristics with typical flat plate collectors at the
1500
G. Martinopoulos et al. / Renewable Energy 35 (2010) 1499–1508
Solar Radiation
Glass or Plastic Cover
Frame
Honeycomb tubing, Black Thermal fluid channels
Insulation
Fig. 1. Collector schematic.
same operational uses and for a wide variety of operating temperatures but with reduced cost (manufacturing and operational). The polymer collector comprises of all the usual components of a typical solar collector. The main difference is that instead of a metal absorber, a black fluid acts as both absorber and heat carrier. The fluid flows in a transparent polymer honeycomb construction. This particular type of collector combines low cost, low weight and simplicity of manufacture. The cross-section of the collector is shown in Fig. 1. The material should be able to withstand liquid pressure and large temperature variations (267–383 K) over time. Since the thermal fluid acts as the heat exchanger, it was important that the polymer had a refractive index similar or better than that of glass, low emissivity and good durability to weathering from ultraviolet radiation, in order to maximize solar energy absorbancy. Furthermore, low thermal conductivity was required in order to minimize the heat losses to the environment. In order to cover the above requirements, the hydraulic channels were made of a single sheet of transparent, UV stabilized, honeycombed, LEXAN sheet of 10 mm thickness. The upper and lower collector channels were made of semi-transparent acrylic. A 1000/1 solution of black Indian ink in water was used as heat transfer fluid. The glazing consisted of a 3 mm thick solid, transparent, UV stabilized LEXAN sheet while the back insulation consisted of a nanogel filled honeycombed LEXAN sheet of a 10 mm thickness.
Fig. 3. Collector inlet and outlet position.
The sides of the collector were insulated using 30 mm thick extruded polyurethane. All of the above were packed in an aluminium casing. A cross-section of the collector is presented in Fig. 2, while the dimensions of the collector and the inlet and outlet position in Fig. 3. The main technical characteristics of the collector are summarized in Table 1. The performance of the polymer collector was measured at an ISO 9806-1 [10] complying test bed. A schematic diagram of the test bed configuration is shown in Fig. 4. A heating immersion circulator was used for the temperature regulation of the heat transfer solution, in order to maintain the temperature at the collector inlet within a range of 0.1 K of the desired value. The flow was measured with a rototron flow meter and the solar radiation was measured with a precision pyranometer. The signals from all the sensors were collected by a 24-bit data logger and were stored in a computer for further processing. The inlet temperature of the heat carrier solution was appropriately set before every measurement, starting from 293 K (ambient temperature). In order to ensure steady state operation, the system was in operation at least for 10 min before each measurement. During the measurement, the solar radiation was higher than 800 W m2, the optical field of the collector was not shaded by more than 5% and the wind velocity was in the 2–4 m s1 range.
Fig. 2. Collector cross-section.
G. Martinopoulos et al. / Renewable Energy 35 (2010) 1499–1508 Table 1 Collector technical characteristics. Total area Glazing area Absorber area Number of glazings Glazing material Heat carrier fluid Absorber Headers Fluid weight Back insulation Side insulation
1448 m2 1252 m2 1252 m2 1 3 mm solid transparent UV stabilized LEXAN Water and Indian ink solution (1000:1) Transparent UV stabilized honeycombed LEXAN sheet (10 mm) Acrylic 8 mm rectangular w14 kg Honeycombed LEXAN sheet with a 10 mm thickness filled with nanogel (k ¼ 0.018 W mK1) Extruded PU 30 mm (l ¼ 0.03 W/m2 K)
From each set of measurements (i.e. for each inlet temperature), the average efficiency was produced based on the equation:
nc ¼
_ p ðTout Tin Þ mc Q_ u ¼ Ac IT Ac IT
(1)
where, Q_ u : the useful energy Ac: the collector area [m2] IT: solar radiation in the collector plane [W m2] _ mass flow rate [kg s1] m: Cp: specific thermal capacity [kJ kg1 K1] The collector’s flow rate during measurements was set at 1.27 dm3.
obtain a detailed inside look at the flow field development and the temperature distribution in the collector. Such efforts are presented in the literature through the works of Selmi et al. [11], Fan et al. [12], Henderson et al. [13] and Varol and Oztop [14]. In the present work, for the CFD modeling, a fully structured computational grid of, approximately, 1.2 million computational cells and 13 computational blocks was created in order to describe accurately the collector geometry together with the inlet and outlet pipes, as presented in Fig. 5. During the creation of the computational grid, special attention was given to the inlet and outlet pipes and to the regions near the solid walls, where increased variations of the flow quantities are expected, so as to ensure a smooth convergence and to capture properly all fluid flow phenomena, as shown in Fig. 6. More specifically, for the proper modeling of the inlet and outlet pipe regions, 10 different computational blocks were used in order to create a grid of increased orthogonality on the solid surfaces. The outlet pipe region was prolonged in a distance of 15 diameters in order to improve the convergence behavior of the CFD computations since backflow was beginning to appear in the outlet region, when a short outlet pipe was used. Inside the computational domain, the collector channel region was modeled as a porous medium region with a predefined pressure drop law described from the Hagen–Poiseuille formula, Eq. (2), which was added as an additional source term S in the discretized momentum equation, Eq. (3) and in the direction of the flow inside the honeycomb fluid channels, corresponding to the y-direction.
Sx Sy
PM
PM
3. CFD modeling, computations and validation In order to investigate the performance of solar collectors, computational fluid dynamics (CFD) techniques can be used so as to
1501
Sz
PM
¼
vp vx honeycomb
¼
vp vy honeycomb
fluid channels
¼ 1010 mu 32m v ¼ 2 dH porosity
fluid channels vp ¼ vz honeycomb ¼ 1010 mw
Fig. 4. Test bed diagram.
fluid channels
(2)
1502
G. Martinopoulos et al. / Renewable Energy 35 (2010) 1499–1508
a function of temperature, Eq. (4). By consequence, the energy transport equation was also solved.
Density :
rðTÞ ¼ 0:003284948 T 2 þ1:687644 T þ 785:6677
Dynamic viscosity :
mðTÞ ¼ 1:406e 07
T 2 0:0001024062 T þ 0:01895682 CpðTÞ ¼ 8:54235e 05 T 3
Specific heat capacity :
þ0:09488036 T 2 34:2228 T þ 8212:82 Thermal conductivity :
KðTÞ ¼ 9:8135e 06 T 2
þ0:007536807 T 0:7674181
(4)
Additionally, since the flow rate values for the examined cases were corresponding to very low velocity values, with the Reynolds number, Eq. (5), inside the channels being approximately 54, the flow was computed as laminar and steady.
Reynolds number : Re ¼
uD v
(5)
where, u: is the mean flow velocity D: is the hydraulic diameter of the honeycomb fluid channel v: is the fluid kinematic viscosity. At the inlet of the computational domain, a uniform velocity field was imposed corresponding to 0.762 m s1, same as in the experimental measurements, while the temperature was set each time equal to the inlet water temperature of the case under investigation. At the present work, CFD computations were performed for 4 different cases, corresponding to water inlet temperature of 303 K, 313 K, 333 K and 353 K. At the outlet of the computational domain the static pressure was set at ambient pressure and equal to 101,325 Pa. The boundary conditions used for the 4 cases of CFD computation are summarized in Table 2. Due to the large dimensions of the collector, it was necessary to include the gravity force, together with the buoyancy effect due to the large temperature difference, by adding the appropriate source terms in the momentum equations, Eq. (6).
Fig. 5. Collector geometry with inlet and outlet pipes.
vðruÞ vp ! þ divðru u Þ ¼ þ divðmgraduÞ þ Stotal x vt vx vðrvÞ vp ! þ divðrv u Þ ¼ þ divðmgradvÞ þ Stotal y vt vy vðrwÞ vp ! þ divðrw u Þ ¼ þ divðmgradwÞ þ Stotal vt vz
! Sgravity þ Sbuyoancy ¼ r g 1 b T Tref
(6)
(3) where b is the thermal expansion coefficient defined as:
1 vr r vT P
z
Additionally, source terms were added in the x- and z- directions with an extremely high pressure drop coefficient deliberately set so as to ‘help’ the flow align with the honeycomb fluid channels and thus, simulate their effect on the real geometry when only flow towards the y-direction is permitted due to the presence of the fluid channel sidewalls. The use of this porous medium methodology was obligatory since the modeling of the precise collector geometry with the honeycomb fluid channels was computationally not affordable due to the high number of channels and the relatively large collector dimensions which would lead to a grid of extremely high number of computational points and too demanding in CPU and memory. The computational grid was created with the commercial CFD software NUMECA/IGG while the computations were performed with the commercial CFD software FLUENT with the use of the user define function (UDF) programmable mode [15]. The flow field was computed with water (simulating the aquatic solution) being the working fluid with predefined properties as
b¼
(7)
its values taken from thermodynamic tables. For this specific type of collector, the direction of importance is ! the one towards the direction of the gravity acceleration g . Thus, when the collector is placed inclined in a specific angle 4, as presented in Fig. 7, this source term produces contributions to the directions y and z as:
Sy
gravity
þ Sy
buyoancy
Sz
gravity
þ Sz
buyoancy
¼ rg sin 4 1 b T Tref ¼ rg cos 4 1 b T Tref
(8)
Hence, the total source terms in the momentum terms which were used for the computations are:
Stotal Stotal Stotal
x y z
¼ Sx PM ¼ Sy PM þ Sy gravity þ Sy buyoancy ¼ Sz PM þ Sz gravity þ Sz buyoancy
(9)
G. Martinopoulos et al. / Renewable Energy 35 (2010) 1499–1508
1503
Fig. 6. Computational grid for a) the collector main region, b) the inlet region pipe and c) the outlet region pipe.
For all CFD computations the velocity field as well as the temperature and pressure distribution presented the same qualitative behavior and no noticeable variations were observed in relation to the water inlet temperature. As a result, and for space economy reasons, the CFD plots corresponding to 303 K inlet temperature are presented. Typical velocity vector plots and colour contours of the flow field inside the collector are presented in Figs. 8–13. For the representation of the results three planes have been selected, corresponding to the inlet pipe middle plane, the collector middle plane and the outlet pipe middle plane. As it can be seen, the velocity field inside the collector and the inlet and outlet pipes is far from being called a uniform one since
Table 2 Boundary conditions used for the CFD computations. Case
Inlet
Outlet
a/a
Velocity
Water temperature
Static pressure
1 2 3 4
0.762 m s1
303 K 313 K 333 K 353 K
101,325 Pa Fig. 7. Typical representation of the collector placed under a specific angle a of inclination.
1504
G. Martinopoulos et al. / Renewable Energy 35 (2010) 1499–1508
Fig. 8. Inlet pipe middle plane.
a combination of recirculation regions and regions of high velocity is met. This is clearly illustrated in each one of the selected planes. More specifically, at the inlet pipe (Fig. 9), due to the sudden change of geometry, the flow is entering the main area of the collector at a high velocity and not in a smooth way. Thus, directly
Fig. 10. Collector middle plane.
Fig. 9. Velocity vector at the inlet pipe middle plane a) general view, b) enlarged view inside the pipe, c) enlarged view inside the collector.
G. Martinopoulos et al. / Renewable Energy 35 (2010) 1499–1508
1505
Fig. 11. Velocity at the middle plane of the collector.
after the inlet pipe’s exit region, large recirculation regions are formed in the circumferential direction. The effect of these regions on the way the flow is entering to the main collector area (honeycomb fluid channels) is significant. This can be seen in Fig. 11b, where regions of low and high velocity upwards and downwards are presented and are indicated with arrows. Since the flow in the honeycomb fluid channels can move only up- or downwards, the presence of large regions of downwards flow means that a large part of the collector is not operating properly. In addition, the high velocity regions formed, are contributing to the overall pressure losses, as long as the pressure drop is linearly correlated to the normal Y-velocity component, as presented in Eq. (2). The static pressure drop is the sum of the pressure drop due to the flow field inside the porous medium region (collector main region), where the pressure drop is following the Hagen–Poiseuille
formula, and the pressure drop due to gravity, which is the main contributor. The pressure drop due to the honeycomb is relatively small, since the average flow velocity inside the honeycomb is very small (of the order of 0.1 m s1 and lower). Thus, the pressure distribution is almost totally dominated by the gravity effect resulting in a stable picture, as the one shown in Figs. 10 and 14. Even though the static pressure distribution inside the honeycomb is limited, the absolute values of the velocity vary, as presented in Fig. 11. However, the overall effect of the velocity variation is limited, therefore it does not have a noticeable effect in the ‘stable’ distribution of the static pressure presented in Figs. 10 and 14. Thus, the velocity distribution in Fig. 11 is mainly important for the comprehension and investigation of the different flow regions developed inside the collector. This velocity distribution affects the temperature distribution, but does not have a quantitatively important effect to the static pressure.
1506
G. Martinopoulos et al. / Renewable Energy 35 (2010) 1499–1508
Fig. 11. (continued).
The negative effect of the inlet pipe region begins to fade-out, as the flow moves towards the left sidewalls of the collector, Fig. 11a. Thus, after the creation of a characteristically low velocity region located at the middle of the collector, where the flow is almost stagnant, the flow begins to move upwards and towards the outlet pipe exit. A similar pattern is also developed in the right side where an upwards flow is also indicated. However, directly at the left and right side of the collector the flow is ‘trapped’ between the walls of the collector and takes again very low velocity values. Finally, as presented in Fig. 13, the flow enters the outlet pipe forming a ‘vena contracta’ since the flow cannot change direction so as to follow the sharp angle in the outlet pipe and hence, recirculation regions are developed, increasing the pressure losses of the collector. Regarding the static pressure, a typical colour contour is presented in Fig. 14. It can be seen that the distribution of the static
pressure is strongly depending on the height increase. Thus, as the height is increasing the static pressure is decreasing, since an important part of it is converted to hydrostatic pressure. At the same time, the pressure losses are increasing as the flow moves inside the collector. To validate the results of the CFD computations during the experimental investigation of the collector, the flow of the heat transfer medium inside the hydraulic channels was observed. Black dye was used in order to observe flow movement, in conjunction with a video camera. As it is shown in Fig. 15, the same flow patterns inside the collector are presented for the CFD computation and the experimental flow representation. Furthermore, and in order to validate the temperature distribution over the collector area, infrared photography was used. As it can be seen in Fig. 16, the static temperature distribution of the CFD computation is in close agreement with the experimental results.
G. Martinopoulos et al. / Renewable Energy 35 (2010) 1499–1508
1507
Fig. 14. Static pressure colour contour at the middle plane of the collector.
Fig. 12. Outlet pipe middle plane.
The CFD results qualitatively predict the main temperature patterns (low temperature at the sides, high temperature at the middle, wavy patterns at the bottom) while the temperature values overall are in good agreement, with the exception of a small area (Area A in Fig. 16), at the bottom part, where CFD results are lower than the actual temperature values. It has also to be mentioned that, in order
to take the thermographic image, it was necessary to remove the transparent cover of the collector, increasing thus the thermal losses. Thus, it can be concluded that the CFD computation is capable of producing results of good quality regarding both the flow field development and the temperature distribution. In order to draw a conclusion regarding the behavior of the collector, a comparison between the experimental and computed collector efficiency had to be performed. As shown in Fig. 17 the efficiency curves are similar in behavior and close to each other. It has to be reminded at this point that the collector tested was made in-house, therefore the observed minor deviations can, at least partially, be attributed to the limited manufacturing capabilities.
Fig. 13. Velocity vector at the outlet pipe middle plane a) general view, b) enlarged view of large recirculation region inside the pipe, c) enlarged view of small recirculation region inside the pipe, d) pipe outlet.
1508
G. Martinopoulos et al. / Renewable Energy 35 (2010) 1499–1508
0,70
CFD Results Results
y = -5,6385x + 0,6432
Measurements Measurements
0,60
LinearrCurve Linea Curve (CFD) (CFD)
[n]
Linear Curve Curve (Expe (Experimental) rimental)
0,50 0,40
y = -6,0808x + 0,6393 y = -6,0808x + 0,6393
0,30 0
0,01
0,02 0,03 (Tin-Ta)/IT
0,04
0,05
Fig. 17. Comparison of simulated and experimental efficiency curve.
4. Conclusions
Fig. 15. Flow patterns inside the collector.
The paper presents a novel honeycomb polycarbonate collector in which the solar energy is directly absorbed by the black coloured working fluid. The experimentally determined average efficiency of the collector was found to be similar to that of the low-cost flat plate commercially available collectors. Furthermore, the collector was modeled using CFD techniques, in order to obtain a more detailed inside look in the flow field development and the temperature distribution. The experimental and CFD results comparison showed that there was a reasonably good agreement regarding the collector efficiency while, at the same time, significantly similar flow patterns and temperature distributions were observed Thus, it can be concluded that the computational model is capable of providing reliable results both qualitatively and quantitatively. Both the experimental set up and the CFD computations revealed that there were large problematic regions inside the solar collector which decrease its efficiency, suggesting that the CFD model is a useful tool for further investigations and optimisation of the collector. References
Fig. 16. Static temperature colour contour at the middle plane of the collector.
[1] International Energy Association. Solar heat worldwide 2008; 2007. [2] European Solar Thermal Industry Federation. Solar thermal markets in Europe 2007; June 2008. [3] Hottel HC, Woertz BB. Performance of flat plate solar heat collectors. Transactions of the ASME 1942;64:91–104. [4] Bliss RW. The derivations of several ‘plateefficiency factors’ useful in the design of flat plate solar heat collectors. Solar Energy 1959;3:55–64. [5] Nahar NM, Garg HP. Free convection and shading due to gap spacing between absorber plate and cover glazing in solar energy flat plate collectors. Applied Energy 1980;7:129–45. [6] Francken JC. On the effectiveness of a flat plate collector. Solar Energy 1984;33:363–6. [7] Alghoul MA, Sulaiman MY, Azmi BZ, Abd. Wahab M. Review of materials for solar thermal collectors. Anti-Corrosion Methods and Materials 2005;52/ 4:199–206. [8] Tsilingiris PT. Design, analysis and performance of low-cost plastic film large solar water heating systems. Solar Energy 1997;60(5):245–56. [9] Meir M, Rekstad J. Der Solarnor Kunststoffkollektor – The development of a polymer collector with glazing. In: Proceedings of Polymeric Solar Materials, Erstes Leobener Symposium, Solartechnik – Neue Mo¨glichkeiten fu¨r die Kunststoffbranche, Polymer Competence Center Leoben, 7.–8.10.2003, Leoben; 2003. p. II-1–8. [10] ISO, ISO 9806-1: test methods for solar collectors; 1994. [11] Selmi Mohamed, Al-Khawaja Mohammed J, Marafia Abdulhamid. Validation of CFD simulation for flat plate solar energy collector. Renewable Energy 2008;33:383–7. [12] Fan Jianhua, Shah Louise Jivan, Furbo Simon. Flow distribution in a solar collector panel with horizontally inclined absorber strips. Solar Energy 2007;81:1501–11. [13] Henderson D, Junaidi H, Muneer T, Grassie T, Currie J. Experimental and CFD investigation of an ICSSWH at various inclinations. Renewable and Sustainable Energy Reviews 2007;11:1087–116. [14] Varol Yasin, Oztop Hakan. Buoyancy induced heat transfer and fluid flow inside a tilted wavy solar collector. Building and Environment 2007;42:2062–71. [15] Fluent 6.3, Commercial software documentation; 2006.
Contents lists available at ScienceDirect
Renewable Energy journal homepage: www.elsevier.com/locate/renene
CFD modeling of a polymer solar collector G. Martinopoulos a, D. Missirlis b, G. Tsilingiridis a, K. Yakinthos b, N. Kyriakis a, * a b
Process Equipment Design Laboratory, Dept. of Mechanical Engineering, AUTh, POB 487, PC 54124, Thessaloniki, Greece Laboratory of Fluid Mechanics and Turbomachinery, Dept. of Mechanical Engineering, AUTh, PC 54124, Thessaloniki, Greece
a r t i c l e i n f o
a b s t r a c t
Article history: Received 20 May 2009 Accepted 5 January 2010 Available online 12 February 2010
A polymer solar collector is developed and its behavior is investigated both experimentally and with computational fluid dynamics (CFD). Solar irradiation as well as convection and heat transfer in the circulating fluid and between the parts of the collector is considered in the model. The temperature and velocity distribution over its area as well as the collector efficiency at nominal flow rate were used in order to validate the CFD model. Temperature distribution during operation and average collector efficiency were found to be in good agreement between the experimental data and the results of the CFD modeling. Ó 2010 Elsevier Ltd. All rights reserved.
Keywords: Solar energy Solar collector CFD simulation
1. Introduction As the oil market becomes unstable and unpredictable, the use of renewable energy sources, solar energy for thermal applications in particular, offers an interesting solution. Solar energy in the form of domestic solar thermal systems accounted for 127.8 GWth worldwide with more than 15 GWth in the EU alone, in 2006 [1,2]. The increase in the use of solar thermal energy systems is subjected to their economic viability, in comparison with conventional systems, and their simplicity in manufacturing and use. In order to minimize the production cost of current solar collectors, new materials of lower cost and preferably more environmental friendly must substitute the ones that are employed today, without any deterioration of collector efficiency. In this paper a novel polymer solar collector is presented and computational fluid dynamics are used in an effort to investigate its behavior. The performance of the collector under nominal flow rate was obtained by CFD. For the validation of the CFD model, the temperature and velocity distributions over its area were compared to experimental data. 2. Experimental model The most commonly used solar system is the one producing hot water for domestic use. The basic unit of such a system is the solar collector. Among the solar collector configurations used, the most common is the flat plate one [2]. Flat plate collectors have been
* Corresponding author. Tel.: þ30 2310 996083; fax: þ30 2310 996087. E-mail address: [email protected] (N. Kyriakis). 0960-1481/$ – see front matter Ó 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.renene.2010.01.004
studied by many researchers, e.g. Hottel and Woertz [3], Bliss [4], Nahar and Garg [5], Francken [6] and others. Typical flat plate collectors are either corrugated, bond duct or tube-in-plate type and their performance depends on various design parameters, including the number of covers, the type and thickness of glazing, the anti-reflecting coating on cover glass, the type of coating on the collector plate, the spacing between the collector and the inner glass, the type and thickness of insulation used, etc. Operation parameters affecting the performance of a solar collector are the mass flow rate of fluid, the amount of incident solar radiation, the inlet and ambient temperatures and sky conditions, to name a few. A number of researchers have adopted the use of polymers in solar collector design since their cost and physical properties make the volume production of lightweight low-cost and corrosion resistance collectors possible. The use of polymer glazing has been moderately successful, as they offer significant potential for cost savings both as direct substitutes for glass cover plates in traditional collector systems and as an integral part of all-polymeric systems. Polymer heat exchangers offer the potential advantages of reduced cost of materials and manufacture, resistance to corrosion, as well as better integration with other components [7]. Furthermore the use of polymer materials reduces the collector weight by 50% in comparison with a traditional metal collector allowing much easier installation [8]. Currently the use of polymers as solar collector materials is limited mainly in swimming pool collectors and in glazing of typical flat plate collectors [9]. In this work a novel polymer collector was developed aiming at producing a solar collector which combines similar or better technical characteristics with typical flat plate collectors at the
1500
G. Martinopoulos et al. / Renewable Energy 35 (2010) 1499–1508
Solar Radiation
Glass or Plastic Cover
Frame
Honeycomb tubing, Black Thermal fluid channels
Insulation
Fig. 1. Collector schematic.
same operational uses and for a wide variety of operating temperatures but with reduced cost (manufacturing and operational). The polymer collector comprises of all the usual components of a typical solar collector. The main difference is that instead of a metal absorber, a black fluid acts as both absorber and heat carrier. The fluid flows in a transparent polymer honeycomb construction. This particular type of collector combines low cost, low weight and simplicity of manufacture. The cross-section of the collector is shown in Fig. 1. The material should be able to withstand liquid pressure and large temperature variations (267–383 K) over time. Since the thermal fluid acts as the heat exchanger, it was important that the polymer had a refractive index similar or better than that of glass, low emissivity and good durability to weathering from ultraviolet radiation, in order to maximize solar energy absorbancy. Furthermore, low thermal conductivity was required in order to minimize the heat losses to the environment. In order to cover the above requirements, the hydraulic channels were made of a single sheet of transparent, UV stabilized, honeycombed, LEXAN sheet of 10 mm thickness. The upper and lower collector channels were made of semi-transparent acrylic. A 1000/1 solution of black Indian ink in water was used as heat transfer fluid. The glazing consisted of a 3 mm thick solid, transparent, UV stabilized LEXAN sheet while the back insulation consisted of a nanogel filled honeycombed LEXAN sheet of a 10 mm thickness.
Fig. 3. Collector inlet and outlet position.
The sides of the collector were insulated using 30 mm thick extruded polyurethane. All of the above were packed in an aluminium casing. A cross-section of the collector is presented in Fig. 2, while the dimensions of the collector and the inlet and outlet position in Fig. 3. The main technical characteristics of the collector are summarized in Table 1. The performance of the polymer collector was measured at an ISO 9806-1 [10] complying test bed. A schematic diagram of the test bed configuration is shown in Fig. 4. A heating immersion circulator was used for the temperature regulation of the heat transfer solution, in order to maintain the temperature at the collector inlet within a range of 0.1 K of the desired value. The flow was measured with a rototron flow meter and the solar radiation was measured with a precision pyranometer. The signals from all the sensors were collected by a 24-bit data logger and were stored in a computer for further processing. The inlet temperature of the heat carrier solution was appropriately set before every measurement, starting from 293 K (ambient temperature). In order to ensure steady state operation, the system was in operation at least for 10 min before each measurement. During the measurement, the solar radiation was higher than 800 W m2, the optical field of the collector was not shaded by more than 5% and the wind velocity was in the 2–4 m s1 range.
Fig. 2. Collector cross-section.
G. Martinopoulos et al. / Renewable Energy 35 (2010) 1499–1508 Table 1 Collector technical characteristics. Total area Glazing area Absorber area Number of glazings Glazing material Heat carrier fluid Absorber Headers Fluid weight Back insulation Side insulation
1448 m2 1252 m2 1252 m2 1 3 mm solid transparent UV stabilized LEXAN Water and Indian ink solution (1000:1) Transparent UV stabilized honeycombed LEXAN sheet (10 mm) Acrylic 8 mm rectangular w14 kg Honeycombed LEXAN sheet with a 10 mm thickness filled with nanogel (k ¼ 0.018 W mK1) Extruded PU 30 mm (l ¼ 0.03 W/m2 K)
From each set of measurements (i.e. for each inlet temperature), the average efficiency was produced based on the equation:
nc ¼
_ p ðTout Tin Þ mc Q_ u ¼ Ac IT Ac IT
(1)
where, Q_ u : the useful energy Ac: the collector area [m2] IT: solar radiation in the collector plane [W m2] _ mass flow rate [kg s1] m: Cp: specific thermal capacity [kJ kg1 K1] The collector’s flow rate during measurements was set at 1.27 dm3.
obtain a detailed inside look at the flow field development and the temperature distribution in the collector. Such efforts are presented in the literature through the works of Selmi et al. [11], Fan et al. [12], Henderson et al. [13] and Varol and Oztop [14]. In the present work, for the CFD modeling, a fully structured computational grid of, approximately, 1.2 million computational cells and 13 computational blocks was created in order to describe accurately the collector geometry together with the inlet and outlet pipes, as presented in Fig. 5. During the creation of the computational grid, special attention was given to the inlet and outlet pipes and to the regions near the solid walls, where increased variations of the flow quantities are expected, so as to ensure a smooth convergence and to capture properly all fluid flow phenomena, as shown in Fig. 6. More specifically, for the proper modeling of the inlet and outlet pipe regions, 10 different computational blocks were used in order to create a grid of increased orthogonality on the solid surfaces. The outlet pipe region was prolonged in a distance of 15 diameters in order to improve the convergence behavior of the CFD computations since backflow was beginning to appear in the outlet region, when a short outlet pipe was used. Inside the computational domain, the collector channel region was modeled as a porous medium region with a predefined pressure drop law described from the Hagen–Poiseuille formula, Eq. (2), which was added as an additional source term S in the discretized momentum equation, Eq. (3) and in the direction of the flow inside the honeycomb fluid channels, corresponding to the y-direction.
Sx Sy
PM
PM
3. CFD modeling, computations and validation In order to investigate the performance of solar collectors, computational fluid dynamics (CFD) techniques can be used so as to
1501
Sz
PM
¼
vp vx honeycomb
¼
vp vy honeycomb
fluid channels
¼ 1010 mu 32m v ¼ 2 dH porosity
fluid channels vp ¼ vz honeycomb ¼ 1010 mw
Fig. 4. Test bed diagram.
fluid channels
(2)
1502
G. Martinopoulos et al. / Renewable Energy 35 (2010) 1499–1508
a function of temperature, Eq. (4). By consequence, the energy transport equation was also solved.
Density :
rðTÞ ¼ 0:003284948 T 2 þ1:687644 T þ 785:6677
Dynamic viscosity :
mðTÞ ¼ 1:406e 07
T 2 0:0001024062 T þ 0:01895682 CpðTÞ ¼ 8:54235e 05 T 3
Specific heat capacity :
þ0:09488036 T 2 34:2228 T þ 8212:82 Thermal conductivity :
KðTÞ ¼ 9:8135e 06 T 2
þ0:007536807 T 0:7674181
(4)
Additionally, since the flow rate values for the examined cases were corresponding to very low velocity values, with the Reynolds number, Eq. (5), inside the channels being approximately 54, the flow was computed as laminar and steady.
Reynolds number : Re ¼
uD v
(5)
where, u: is the mean flow velocity D: is the hydraulic diameter of the honeycomb fluid channel v: is the fluid kinematic viscosity. At the inlet of the computational domain, a uniform velocity field was imposed corresponding to 0.762 m s1, same as in the experimental measurements, while the temperature was set each time equal to the inlet water temperature of the case under investigation. At the present work, CFD computations were performed for 4 different cases, corresponding to water inlet temperature of 303 K, 313 K, 333 K and 353 K. At the outlet of the computational domain the static pressure was set at ambient pressure and equal to 101,325 Pa. The boundary conditions used for the 4 cases of CFD computation are summarized in Table 2. Due to the large dimensions of the collector, it was necessary to include the gravity force, together with the buoyancy effect due to the large temperature difference, by adding the appropriate source terms in the momentum equations, Eq. (6).
Fig. 5. Collector geometry with inlet and outlet pipes.
vðruÞ vp ! þ divðru u Þ ¼ þ divðmgraduÞ þ Stotal x vt vx vðrvÞ vp ! þ divðrv u Þ ¼ þ divðmgradvÞ þ Stotal y vt vy vðrwÞ vp ! þ divðrw u Þ ¼ þ divðmgradwÞ þ Stotal vt vz
! Sgravity þ Sbuyoancy ¼ r g 1 b T Tref
(6)
(3) where b is the thermal expansion coefficient defined as:
1 vr r vT P
z
Additionally, source terms were added in the x- and z- directions with an extremely high pressure drop coefficient deliberately set so as to ‘help’ the flow align with the honeycomb fluid channels and thus, simulate their effect on the real geometry when only flow towards the y-direction is permitted due to the presence of the fluid channel sidewalls. The use of this porous medium methodology was obligatory since the modeling of the precise collector geometry with the honeycomb fluid channels was computationally not affordable due to the high number of channels and the relatively large collector dimensions which would lead to a grid of extremely high number of computational points and too demanding in CPU and memory. The computational grid was created with the commercial CFD software NUMECA/IGG while the computations were performed with the commercial CFD software FLUENT with the use of the user define function (UDF) programmable mode [15]. The flow field was computed with water (simulating the aquatic solution) being the working fluid with predefined properties as
b¼
(7)
its values taken from thermodynamic tables. For this specific type of collector, the direction of importance is ! the one towards the direction of the gravity acceleration g . Thus, when the collector is placed inclined in a specific angle 4, as presented in Fig. 7, this source term produces contributions to the directions y and z as:
Sy
gravity
þ Sy
buyoancy
Sz
gravity
þ Sz
buyoancy
¼ rg sin 4 1 b T Tref ¼ rg cos 4 1 b T Tref
(8)
Hence, the total source terms in the momentum terms which were used for the computations are:
Stotal Stotal Stotal
x y z
¼ Sx PM ¼ Sy PM þ Sy gravity þ Sy buyoancy ¼ Sz PM þ Sz gravity þ Sz buyoancy
(9)
G. Martinopoulos et al. / Renewable Energy 35 (2010) 1499–1508
1503
Fig. 6. Computational grid for a) the collector main region, b) the inlet region pipe and c) the outlet region pipe.
For all CFD computations the velocity field as well as the temperature and pressure distribution presented the same qualitative behavior and no noticeable variations were observed in relation to the water inlet temperature. As a result, and for space economy reasons, the CFD plots corresponding to 303 K inlet temperature are presented. Typical velocity vector plots and colour contours of the flow field inside the collector are presented in Figs. 8–13. For the representation of the results three planes have been selected, corresponding to the inlet pipe middle plane, the collector middle plane and the outlet pipe middle plane. As it can be seen, the velocity field inside the collector and the inlet and outlet pipes is far from being called a uniform one since
Table 2 Boundary conditions used for the CFD computations. Case
Inlet
Outlet
a/a
Velocity
Water temperature
Static pressure
1 2 3 4
0.762 m s1
303 K 313 K 333 K 353 K
101,325 Pa Fig. 7. Typical representation of the collector placed under a specific angle a of inclination.
1504
G. Martinopoulos et al. / Renewable Energy 35 (2010) 1499–1508
Fig. 8. Inlet pipe middle plane.
a combination of recirculation regions and regions of high velocity is met. This is clearly illustrated in each one of the selected planes. More specifically, at the inlet pipe (Fig. 9), due to the sudden change of geometry, the flow is entering the main area of the collector at a high velocity and not in a smooth way. Thus, directly
Fig. 10. Collector middle plane.
Fig. 9. Velocity vector at the inlet pipe middle plane a) general view, b) enlarged view inside the pipe, c) enlarged view inside the collector.
G. Martinopoulos et al. / Renewable Energy 35 (2010) 1499–1508
1505
Fig. 11. Velocity at the middle plane of the collector.
after the inlet pipe’s exit region, large recirculation regions are formed in the circumferential direction. The effect of these regions on the way the flow is entering to the main collector area (honeycomb fluid channels) is significant. This can be seen in Fig. 11b, where regions of low and high velocity upwards and downwards are presented and are indicated with arrows. Since the flow in the honeycomb fluid channels can move only up- or downwards, the presence of large regions of downwards flow means that a large part of the collector is not operating properly. In addition, the high velocity regions formed, are contributing to the overall pressure losses, as long as the pressure drop is linearly correlated to the normal Y-velocity component, as presented in Eq. (2). The static pressure drop is the sum of the pressure drop due to the flow field inside the porous medium region (collector main region), where the pressure drop is following the Hagen–Poiseuille
formula, and the pressure drop due to gravity, which is the main contributor. The pressure drop due to the honeycomb is relatively small, since the average flow velocity inside the honeycomb is very small (of the order of 0.1 m s1 and lower). Thus, the pressure distribution is almost totally dominated by the gravity effect resulting in a stable picture, as the one shown in Figs. 10 and 14. Even though the static pressure distribution inside the honeycomb is limited, the absolute values of the velocity vary, as presented in Fig. 11. However, the overall effect of the velocity variation is limited, therefore it does not have a noticeable effect in the ‘stable’ distribution of the static pressure presented in Figs. 10 and 14. Thus, the velocity distribution in Fig. 11 is mainly important for the comprehension and investigation of the different flow regions developed inside the collector. This velocity distribution affects the temperature distribution, but does not have a quantitatively important effect to the static pressure.
1506
G. Martinopoulos et al. / Renewable Energy 35 (2010) 1499–1508
Fig. 11. (continued).
The negative effect of the inlet pipe region begins to fade-out, as the flow moves towards the left sidewalls of the collector, Fig. 11a. Thus, after the creation of a characteristically low velocity region located at the middle of the collector, where the flow is almost stagnant, the flow begins to move upwards and towards the outlet pipe exit. A similar pattern is also developed in the right side where an upwards flow is also indicated. However, directly at the left and right side of the collector the flow is ‘trapped’ between the walls of the collector and takes again very low velocity values. Finally, as presented in Fig. 13, the flow enters the outlet pipe forming a ‘vena contracta’ since the flow cannot change direction so as to follow the sharp angle in the outlet pipe and hence, recirculation regions are developed, increasing the pressure losses of the collector. Regarding the static pressure, a typical colour contour is presented in Fig. 14. It can be seen that the distribution of the static
pressure is strongly depending on the height increase. Thus, as the height is increasing the static pressure is decreasing, since an important part of it is converted to hydrostatic pressure. At the same time, the pressure losses are increasing as the flow moves inside the collector. To validate the results of the CFD computations during the experimental investigation of the collector, the flow of the heat transfer medium inside the hydraulic channels was observed. Black dye was used in order to observe flow movement, in conjunction with a video camera. As it is shown in Fig. 15, the same flow patterns inside the collector are presented for the CFD computation and the experimental flow representation. Furthermore, and in order to validate the temperature distribution over the collector area, infrared photography was used. As it can be seen in Fig. 16, the static temperature distribution of the CFD computation is in close agreement with the experimental results.
G. Martinopoulos et al. / Renewable Energy 35 (2010) 1499–1508
1507
Fig. 14. Static pressure colour contour at the middle plane of the collector.
Fig. 12. Outlet pipe middle plane.
The CFD results qualitatively predict the main temperature patterns (low temperature at the sides, high temperature at the middle, wavy patterns at the bottom) while the temperature values overall are in good agreement, with the exception of a small area (Area A in Fig. 16), at the bottom part, where CFD results are lower than the actual temperature values. It has also to be mentioned that, in order
to take the thermographic image, it was necessary to remove the transparent cover of the collector, increasing thus the thermal losses. Thus, it can be concluded that the CFD computation is capable of producing results of good quality regarding both the flow field development and the temperature distribution. In order to draw a conclusion regarding the behavior of the collector, a comparison between the experimental and computed collector efficiency had to be performed. As shown in Fig. 17 the efficiency curves are similar in behavior and close to each other. It has to be reminded at this point that the collector tested was made in-house, therefore the observed minor deviations can, at least partially, be attributed to the limited manufacturing capabilities.
Fig. 13. Velocity vector at the outlet pipe middle plane a) general view, b) enlarged view of large recirculation region inside the pipe, c) enlarged view of small recirculation region inside the pipe, d) pipe outlet.
1508
G. Martinopoulos et al. / Renewable Energy 35 (2010) 1499–1508
0,70
CFD Results Results
y = -5,6385x + 0,6432
Measurements Measurements
0,60
LinearrCurve Linea Curve (CFD) (CFD)
[n]
Linear Curve Curve (Expe (Experimental) rimental)
0,50 0,40
y = -6,0808x + 0,6393 y = -6,0808x + 0,6393
0,30 0
0,01
0,02 0,03 (Tin-Ta)/IT
0,04
0,05
Fig. 17. Comparison of simulated and experimental efficiency curve.
4. Conclusions
Fig. 15. Flow patterns inside the collector.
The paper presents a novel honeycomb polycarbonate collector in which the solar energy is directly absorbed by the black coloured working fluid. The experimentally determined average efficiency of the collector was found to be similar to that of the low-cost flat plate commercially available collectors. Furthermore, the collector was modeled using CFD techniques, in order to obtain a more detailed inside look in the flow field development and the temperature distribution. The experimental and CFD results comparison showed that there was a reasonably good agreement regarding the collector efficiency while, at the same time, significantly similar flow patterns and temperature distributions were observed Thus, it can be concluded that the computational model is capable of providing reliable results both qualitatively and quantitatively. Both the experimental set up and the CFD computations revealed that there were large problematic regions inside the solar collector which decrease its efficiency, suggesting that the CFD model is a useful tool for further investigations and optimisation of the collector. References
Fig. 16. Static temperature colour contour at the middle plane of the collector.
[1] International Energy Association. Solar heat worldwide 2008; 2007. [2] European Solar Thermal Industry Federation. Solar thermal markets in Europe 2007; June 2008. [3] Hottel HC, Woertz BB. Performance of flat plate solar heat collectors. Transactions of the ASME 1942;64:91–104. [4] Bliss RW. The derivations of several ‘plateefficiency factors’ useful in the design of flat plate solar heat collectors. Solar Energy 1959;3:55–64. [5] Nahar NM, Garg HP. Free convection and shading due to gap spacing between absorber plate and cover glazing in solar energy flat plate collectors. Applied Energy 1980;7:129–45. [6] Francken JC. On the effectiveness of a flat plate collector. Solar Energy 1984;33:363–6. [7] Alghoul MA, Sulaiman MY, Azmi BZ, Abd. Wahab M. Review of materials for solar thermal collectors. Anti-Corrosion Methods and Materials 2005;52/ 4:199–206. [8] Tsilingiris PT. Design, analysis and performance of low-cost plastic film large solar water heating systems. Solar Energy 1997;60(5):245–56. [9] Meir M, Rekstad J. Der Solarnor Kunststoffkollektor – The development of a polymer collector with glazing. In: Proceedings of Polymeric Solar Materials, Erstes Leobener Symposium, Solartechnik – Neue Mo¨glichkeiten fu¨r die Kunststoffbranche, Polymer Competence Center Leoben, 7.–8.10.2003, Leoben; 2003. p. II-1–8. [10] ISO, ISO 9806-1: test methods for solar collectors; 1994. [11] Selmi Mohamed, Al-Khawaja Mohammed J, Marafia Abdulhamid. Validation of CFD simulation for flat plate solar energy collector. Renewable Energy 2008;33:383–7. [12] Fan Jianhua, Shah Louise Jivan, Furbo Simon. Flow distribution in a solar collector panel with horizontally inclined absorber strips. Solar Energy 2007;81:1501–11. [13] Henderson D, Junaidi H, Muneer T, Grassie T, Currie J. Experimental and CFD investigation of an ICSSWH at various inclinations. Renewable and Sustainable Energy Reviews 2007;11:1087–116. [14] Varol Yasin, Oztop Hakan. Buoyancy induced heat transfer and fluid flow inside a tilted wavy solar collector. Building and Environment 2007;42:2062–71. [15] Fluent 6.3, Commercial software documentation; 2006.