Available online at www.sciencedirect.

com

Solar Energy 83 (2009) 1086–1092

www.elsevier.com/locate/solener

Numerical 3-D heat flux simulations on flat plate solar collectors

N. Molero Villar, J.M. Cejudo López *, F. Domı́nguez Muñoz,
E. Rodrı́guez Garcı́a, A. Carrillo Andrés
Grupo de Energética, Escuela Técnica Superior de Ingenieros Industriales, UMA, Plaza El Ejido s/n, 29013 Málaga, Spain

Received 28 September 2007; received in revised form 20 January 2009; accepted 22 January 2009
Available online 25 February 2009

Communicated by: Associate Editor Brian Norton

Abstract

A transient 3-D mathematical model for solar flat plate collectors has been developed. The model is based on setting mass and energy
balances on finite volumes. The model allows the comparison of different configurations: parallel tubes collectors (PTC), serpentine tube
collectors (STC), two parallel plate collectors (TPPC), and other non-usual possibilities like the use of absorbent fluids with semitrans-
parent or transparent plates. Transparent honeycomb insulation between plate and cover can also be modelled. The effect of temperature
on the thermal properties of the materials has also been considered. The model has been validated experimentally with a commercial
PTC. The model is a useful tool to improve the design of plate solar collectors and to compare different configurations. In order to show
the capabilities of the model, the performance of a PTC collector with non-uniformity flow is analysed and compared with experimental
data from literature with good agreement.
Ó 2009 Elsevier Ltd. All rights reserved.

Keywords: Solar thermal collectors; Mathematical modelling; Flow distribution

1. Introduction tube (PTC) type. The serpentine (STC) and two parallel
plate (TPPC) types are less common, although these typol-
Many configurations and materials can be used in manu- ogies are more convenient from the thermal point of view.
facturing solar flat plate collectors. A unified model able to Normally, solar radiation is absorbed by the absorber
predict the performance of this variety of products is plate and transferred to the working fluid. In some config-
required to select the most suitable configuration for a given urations, radiation is directly absorbed by the fluid. In any
application, and to improve the design once the configura- case, part of the absorbed energy is lost to the ambient by
tion has been selected. Such a model would make it possible convection and long wave radiation. The difference
to compare the efficiency of different types of solar collector between the absorbed radiation and the losses to the ambi-
and would serve as a design tool for selecting configurations ent is the useful energy to the working fluid. Matrawy and
and materials. An adequate combination of modelling and Farkas (1997) justifies a greater efficiency in TPPC than in
experimentation is the way to optimize the design. STC or PTC (6% and 10% respectively). The lower is the
The configuration of a solar flat plate collector is an absorber plate temperature, the higher is the efficiency of
important factor that determines its thermal performance the collector. This fact explains that TPPC, in which the
(Matrawy and Farkas, 1997). Most of the flat plate collec- average temperature of the plate is lower than in the other
tors currently available on the market are of the parallel configurations, achieves better efficiency. In the TPPC con-
figuration, two parallel plates separated by a small distance
(around 2 mm), allow the circulation of the working fluid.
*
Corresponding author. Tel.: +34 952132410; fax: +34 952 13 24 09. The fluid directly refrigerates the whole of the absorber
E-mail address: [email protected] (J.M. Cejudo López). plate. On the other hand, in PTC or STC configurations,

0038-092X/$ - see front matter Ó 2009 Elsevier Ltd. All rights reserved.
doi:10.1016/j.solener.2009.01.014
 N.M. Villar et al. / Solar Energy 83 (2009) 1086–1092 1087

Nomenclature

E internal energy (J) Subscripts

G incident solar radiation (W) a ambient
h enthalpy (j/kg) b back
Q flow rate (kg/s) c cover
q heat flux (w) cond conduction
T temperature (°c) cr convection-radiation
t time (s) f fluid
i,j subscripts of one or two dimensional variable
Greek letters indicating the position in the mesh ith tube of
b fraction of mass flow rate through a riser the collector in Eqs. (6) and (7)
D increment ins insulation
/ non-uniformity parameter (Eq. 7) m average (mean value input–output in Eqs. 5a
g collector efficiency and 5b)
g collector efficiency deterioration factor (Eq. 8) o total (in Eq. 6)
n associates the fluid and absorber nodes (logical p plate
value 0,1) s solar radiation absorbed
t tilted
top upper part of the collector

the fin efficiency represents a thermal resistance that In some cases a transient model of the collector is
increases the temperature of the plate, and consequently needed. For instance, it is necessary to study control
the losses to the ambient. A variation of the traditional problems, to calculate the equivalent capacitance of
TPPC is one in which the back plate is the absorber the collector or to avoid expensive and time-consuming
one. In this case, the efficiency factor increases (Tsilingris, steady-state tests. The models developed by Oliva et al.
2002). (1991), Amer et al. (1996), and Amer and Nayak (1999)
In addition to this variety of configurations, many mate- consider the transient character of the problem.
rials (glass, insulation, etc.) are used in plate thermal collec- In this paper, a general model, flexible enough to com-
tors. The thermal properties of these materials can be pare the thermal behaviour of different configurations of
temperature-dependent. According to O’Brien-Bernini low temperature collectors is presented. The model is based
and McGowan (1984), the use of semitransparent or trans- on finite volume balances and considers the multidimen-
parent plate and absorbent fluid improves the efficiency of sional and transient character of the problem. These two
the collector in the rank of low to medium temperature lev- last phenomena are not considered in most of the classic
els. Transparent insulating materials (TIM) are also a bibliography on solar collector performance (Duffie and
promising technology. Transparent honeycomb insulation Beckman, 2006).
between the absorber plate and the cover is used to reduce The capabilities of the proposed model are:
the loss coefficient towards the atmosphere, although it
also reduces the absorption of solar radiation in the plate – Possible configurations: PTC, STC and TPPC.
(the factor ðsa) is deduced in Kaushika and Arulanantham, – Transient.
1996). It can be very interesting in applications of medium – The effect of flow maldistribution can be address.
temperature, for instance solar cooling. The introduction – Absorbent fluid and transparent or semitransparent
of non-conventional insulating materials in the back of plate are allowed.
the collectors is also interesting. For example, an opaque – TIM can be modelled.
honeycomb layer with reflective films could substitute the – Spatial distribution of temperatures in the plate and
conventional insulation. fluid can be obtained.
A general model able to cope with this variety of pos-
sibilities does not currently exist. In particular, one of
the most complete published models is presented in Oliva 2. Theoretical basis
et al. (1991). The balance in the absorber plate, the
cover, insulation and flow is set separately. The model The model is based on setting the mass and energy bal-
is applied to the TPPC configuration and it is not valid ances in a control volume. The mesh is defined by N  M
for the rest. The same occurs with the models by Reddy elements in the plane of the absorber plate (x and y direc-
et al. (1999), Minn et al. (2002), and Hassan and Beli- tions). The main dimensions of the element (Dx, Dy)
veau (2007). depend on the type of collector. The mesh is adapted to
1088 N.M. Villar et al. / Solar Energy 83 (2009) 1086–1092

the tubes if they exist (PTC and STP types): the outside where:
diameter of the tubes determines the dimensions of the DEf i; j: change of the internal energy of the fluid in the
mesh in one direction. i,j element (J) during the time step Dt (s);
Between adjacent control volumes, forced convection Dhi;j : enthalpy difference between input and output flows
inside the tubes or between plates, and conduction in the (J/kg);
absorber plate, are the heat transfer mechanisms. The mass m_ fi;j : fluid flow rate in the i,j element (kg/s);
balance in a control volume is stated by analysing the flow qf i; j: heat transferred from the plate to the fluid (W);
of the working fluid from one element to the adjacent ones. qsf i; j: In the case of TPPC collectors with semitranspar-
The heat balance of the back insulation is also established ent or transparent plate and absorbent fluid, this term rep-
for control volumes defined in the normal direction to the resents the heat directly absorbed by the fluid (W).
plane of the collector. The variation of thermal properties
with temperature is considered. The boundary condition 2.3. Heat balance of the back insulation
for the upper surface of the insulation, in contact with the
plate, is the mean temperature of the plate. The density In the back insulation one-dimensional conduction is
and thermal capacity of the materials are taken into account considered. On the top surface of the insulation, the tem-
to describe the transient behaviour of the solar collector. The perature is the mean temperature of the absorber plate. A
cover is modelled as a single node. In the case of two covers, thermal resistance between plate and insulation can be
the balance is set in each of the panes. Finally, to model con- imposed. For the back surface of the insulation, a convec-
vection and long wave radiation, standard heat transfer coef- tion-radiation coefficient is used to couple the heat transfer
ficients are used (see for example Duffie and Beckman, 2006). with the ambient. The balance for an inner node of the
Next, the heat balance of each element is outlined. insulation is:
DEins; i
2.1. Heat balance of an absorber plate element ¼ qcond;iþ1 þ qcond:i1 ð3Þ
Dt
In the plate, the absorbed solar radiation is transferred where:
to the working fluid, partially lost to the ambient and par- DEins; i: is the change of internal energy of the i element
tially stored in the plate as a change of its internal energy. (J) in the time step Dt (s);
Eq. (1) represents the balance of one node in the plate. qcond;iþ1 ; qcond:i1 : are the conduction fluxes to the upper
and lower nodes (W).
DEi;j
¼ qcondi1;j þ qcondiþ1;j þ qcondi;j1 þ qcondi;jþ1 þ qsi;j
Dt
2.4. Balance of the cover
 qtopi;j  qbi;j  nij qf i; j ð1Þ
The temperature of each cover is supposed uniform. The
where:
increase of internal energy is the difference between absorp-
DEi;j : is the change of internal energy of the i,j element of
tion of solar radiation and convection-radiation with ambi-
the plate (J) in the time step Dt (s);
ent and plate (for single covers). The possibility of double
qcond i  1; j; qcond i þ 1; j; qcond i; j  1; qcond i; j þ 1: are the
conductive fluxes from adjacent nodes; covers has also been included in the model. In the case of
qs i; j: is the absorbed solar energy in the element i,j (W); one cover, the balance equation is:
qtop i; j; qb i; j: losses to the ambient (top and back). They DEc
¼ qsc  qcr;ca  qcr;cp ð4Þ
are calculated by using a global heat transfer coefficient. If Dt
the absorber plate is in contact with the insulation, the heat
where:
flux toward the back is coupled with the first node of the
DEc : is the change of internal energy of the cover (J) in
insulation (W);
the time step Dt (s);
nij qf i; j: is the useful energy to the working fluid. nij ¼ 1
qsc : is the solar energy absorbed by the cover (W);
if the node i,j is a volume attached to a fluid node; zero
qcr;ca : is the convection and long wave radiation flux
otherwise (W);
between the cover and the ambient (W);
In the case of nodes in the perimeter of the collector,
qcr;cp : is the convection and long wave radiation flux
edge losses to the ambient are modelled by means of an
between the cover and the plate (W).
edge loss coefficient.
2.5. Solution of the system of equations
2.2. Heat balance in the fluid
A time dependent linear system of equation is set. The
The equation for the fluid heat balance in the element i,j
number of equations is the same that the number of vol-
is:
umes. The temperatures of the nodes are the dependent
DEfi;j variables. The type of collector, thermal properties of
¼ m_ fi;j Dhi;j þ qf i; j þ qsf i; j ð2Þ
Dt materials, dimensions and tilt of the collector are the
 N.M. Villar et al. / Solar Energy 83 (2009) 1086–1092 1089

parameters of the system. The system is closed by defining Very good results are obtained with the model when is
the initial temperature of the collector, the inlet tempera- compared with experiments. The maximum relative differ-
ture of the fluid and the climate conditions (solar irradi- ence is 0.33% for the prediction of the outlet temperature
ance, surrounding air temperature and humidity and air with an absolute difference of 0.32 °C. In Fig. 1 the exper-
speed). The number of control volumes defined in the plate imental and modelled efficiency curves are shown. The effi-
depends on the type and dimensions of the collector. ciency curves are based on absorber area. The maximum
For a collector with an absorber plate 2.205  1 m2 and differences between both curves are observed in the high
10 risers, the mesh has 24 elements along the y-axis (length) temperature region. At these temperatures, the conductiv-
and 10 nodes between each two parallel risers. The total ity of the back and edge insulation increases significantly
number of elements in the absorber plate is 2649. The mesh and the flow regime inside the risers approaches to the tran-
is optimized by taking smaller elements where the temper- sition region. Uncertainties in the conductivity and convec-
ature gradients are higher. In the fluid a total of 430 vol- tion coefficient could explain the mismatch between test
umes have been defined, 55 in the insulation, and one in and model.
the cover. The mesh is refined until the collector efficiency The efficiency curves obtained with the test and the
calculated with a finer mesh differs less than 104. model are (a quadratic function is used):
The model is iterative because the convective and radia- Tm  Ta
tive coefficients depend on the temperature. The same gtest ¼ 0:7732  3:2368  0:010437
Gt
occurs for some thermal properties like the conductivity 2
of the insulation. For each time step, an iterative loop ðT m  T a Þ
 ð5aÞ
recalculates these coefficients and properties until conver- Gt
gence is reached. Tm  Ta
For stationary conditions, a total of three iterations are gmodel ¼ 0:7697  3:1331  0:005866
Gt
needed for each of the points represented in the Fig. 1, and 2
seven iterations to calculate the stagnation temperature. In ðT m  T a Þ
 ð5bÞ
transient simulations, if the time interval is small enough, Gt
no iteration is needed. The number of iterations is con- The stagnation temperature has also been calculated. In
trolled to a tolerance of 104 for the efficiency. this case, the average temperature of the absorber plate has
About 35 s are needed to calculate one point of the effi- been calculated and compared with the test. The maximum
ciency curve. The simulations were performed in a com- relative deviation is 2.86%.
puter Intel(R) Core(TM)2 CPU 6400 with 2.13 GHz and
1.98 GB of RAM memory. 4. Study of non-uniform flow

3. Experimental validation The classical models of PTC collectors assume a uni-

form distribution of flow through the risers. In practice,
In this paper, just a few capabilities of the model will be there is always some degree of misdistribution, as shown
presented. In particular, only steady-state conditions are by Chiou (1982), Wang and Wu (1990), Jones and Lior
analysed. To validate the model in this condition, a com- (1994), Weitbrecht et al. (2002) and Fan et al. (2007). In
mercial PTC collector has been selected: Isotherm Plus the particular case of PTC with the named Z configuration
made by the Spanish company Isofotón. The test for this (Wang and Wu, 1990), the flow through each riser depends
collector was carried out by Fraunhofer Institut in 2006. on the diameters of the raiser and the manifold, and the
number and length of the risers (see Jones and Lior,

0.80
0.25
0.75
Uniform flow
0.70 0.20 Non uniform flow 0.0124
Collector efficiency

Non uniform flow 0.0286

0.65 Non uniform flow 0.0449
0.15 Non uniform flow 0.0612
test
0.60
βi

simulation
0.55 0.10

0.50
0.05
0.45

0.40 0.00
0 0.02 0.04 0.06 0.08 1 2 3 4 5 6 7 8 9 10
(Tm - Ta)/Gt riser

Fig. 1. Efficiency curve of the solar collector. Test and simulation. Fig. 2. Flow rate distribution in the collector.
1090 N.M. Villar et al. / Solar Energy 83 (2009) 1086–1092

Table 1 1994). In the present paper, the misdistribution calculated

Operating conditions. by Jones will be assumed and its effect on the efficiency
Gt (W/m2) Difuse/total (%) Tf input (°C) Ta (°C) analysed.
950 16 20 25 In order to define the flow through each riser, a param-
eter bi will be defined as the ratio between the flow through
the ith tube (Qi) and the total flow (Q0):
Qi
bi ¼ ð6Þ
0.025 Q0

0.020 A parameter that measures the non-uniformity is

η* (80 kg/h) defined as follows (the number of risers is 10), see Chiou
(1982):
0.015
η* (162 kg/h) sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
P10
η*

2
i¼1 ðbi  1=10Þ
0.010 /¼ ð7Þ
10
0.005 The flow distribution depends on the flow rate through
the collector. If the flow increases, the misdistribution
0.000 decreases. For the Z-type collector, Jones and Lior (1994)
0.01 0.03 0.04 0.06 have studied the non-uniformity. Based on this previous
work, four cases of non-uniform flow have been defined,
Fig. 3. Collector efficiency deterioration factor. (the values of bi and Ui are shown in Fig. 2). The operating

Fig. 4. Temperature distribution (°C) in the absorber plate with uniform flow and non-uniform flow (U = 0.0612) for the operating point presented in
Table 1. (Q0 = 162 kg/h).
 N.M. Villar et al. / Solar Energy 83 (2009) 1086–1092 1091

50 model is based on setting the heat balance in small vol-

riser 1 umes. The model is flexible enough to deal with all types
Fluid temperature (ºC)

45
riser 2 of flat plate collector. Parallel tubes, serpentine or parallel
riser 3
40 riser 4
riser 5
plates configurations can be treated. The possibility of
riser 6
riser 7
absorbent fluid with transparent plate is also included, as
35
riser 8 well as the use of transparent insulation materials. The
riser 9
30 riser 10 model is transient and three-dimensional and allows the
calculation of the temperature of components of the collec-
25
tor (plate, fluid, insulation, etc.) in any position. As a con-
20 sequence of this flexibility, the model is very useful in
0.0 0.5 1.0 1.5 2.0 comparing configurations or abnormal operation
Length of the absorber plate (m)
conditions.
Fig. 5. Fluid temperature in the tubes (U = 0.0612). The model has been experimentally validated for a com-
mercial PTC by comparing the results of the model with
conditions are the same for all the cases (see Table 1). the test made in Fraunhofer Institut. Agreement between
Finally, to quantify the effect on the efficiency of the non- the model and the test is good. In order to show the capa-
uniform flow, the degree of deterioration of collector effi- bilities of the developed model, the collector performance
ciency is defined as follows: when the flow rate in the risers is non-uniform has been
g  gnon-uniform flow studied and compared with experimental data from litera-
g ¼ uniform flow ð8Þ ture. The collector efficiency does not change appreciably
guniform flow
even when the flow at the outer risers is 1.5 times the flow
Fig. 3 shows the collector efficiency deterioration factor of the central one but the outlet temperatures for each tube
depending on non-uniformity of the flow distribution. The are very dissimilar.
higher is the non-uniformity the higher is the deterioration
of the efficiency, although this effect is very limited. When References
U = 0.0612, the efficiency decreases by 2% if the flow rate is
80 kg/h (36.4 kg/hm2). If the flow rate is 162 kg/h (73.5 kg/ Amer, E.H., Nayak, J.K., 1999. Evaluation of a transient test procedure
hm2), the maximum degree of deterioration is only 1%. for solar flat-plate collectors. Energy 24, 979–995.
These values are of the same order of magnitude than in Amer, E.H., Jadeja, P., Kayak, J.K., Sharma, G.K., 1996. Comparison of
Chiou, (1982) (U = 0.06; g* = 0.28; with flow rate two dynamic test methods for solar flat-plate collectors. Energy 39,
285–293.
16.3 kg/hm2). Chiou, J.P., 1982. The effect of nonuniform fluid flow distribution on the
The model also allows the analysis of absorber temper- thermal performance of solar collector. Solar Energy 29 (6), 487–502.
atures. In Fig. 4 a map of the absorber temperatures is Duffie, J.A., Beckman, W.A., 2006. Solar Engineering of Thermal
shown with uniform and non-uniform flow (U = 0.0612) Processes. Wiley Interscience, New York.
for the work point presented in Table 1 and a flow rate Fan, J., Shah, L.J., Furbo, S., 2007. Flow distribution in a solar collector
panel with horizontally inclined absorber strips. Solar Energy 81,
of 162 kg/h. This feature is particularly interesting when 1501–1511.
analysing overheating of collectors that affect its durability. Fraunhofer Institut Solare Energiesysteme, 2006. Test Report: KTB No.
In the case of having uniform flow, the increase of tem- 20006-08-en.
perature experienced by the fluid in each riser is practically Hassan, M.M., Beliveau, Y., 2007. Design, construction and perfor-
the same in all of them. Only slight differences are apprecia- mance prediction of integrated solar roof collectors using finite
element analysis. Construction and Building Materials 21, 1067–
ble in the edge risers due to lateral losses. The outlet tem- 1078.
perature difference between risers is less than 0.4 °C. In Jones, G.F., Lior, N.L., 1994. Flow distribution in manifolded solar
the case of non-uniform flow, the differential heating of collectors with negligible buoyancy effects. Solar Energy 52 (3), 289–
the fluid between risers is proportional to the non-unifor- 300.
mity level. The Fig. 5 represents the fluids temperatures Kaushika, N.D., Arulanantham, M., 1996. Transmittance-absorptance
product of solar glazing with transparent insulation materials. Solar
along the risers for the case of U = 0.0612. In this case, Energy Materials and Solar Cells 44, 383–395.
the differences between risers are higher than 20 °C. The Matrawy, K.K., Farkas, I., 1997. Comparison study for three types of
lesser is the flow in the tube, the higher is the output tem- solar collector for water heating. Energy Conversion and Management
perature. When the energy balance is set for the whole col- 38, 861–869.
lector, the higher output temperature is compensated with Minn, M.A., Ng, K.C., Khong, W.H., Melvin, T., 2002. A distributed
model for a Tedlar-foil flat plate solar collector. Renewable Energy 27,
the lower rate. The consequence is a small value of the effi- 507–523.
ciency deterioration factor. O’Brien-Bernini, F.C., McGowan, J.G., 1984. Performance modelling of
non-metallic flat plate solar collectors. Solar Energy 33, 305–319.
Oliva, A., Costa, M., Pérez-Segarra, C.D., 1991. Numerical simulation of
5. Conclusions
solar collectors: the effect of nonuniform and nonsteady state or the
boundary conditions. Solar Energy 5, 359–373.
In this paper, a powerful model to quantify the thermal Reddy, K.S., Avanti, P., Kaushika, N.D., 1999. Finite time thermal
performance of solar plate collectors is presented. The analysis of ground integrated-collector-storage solar water heater with
1092 N.M. Villar et al. / Solar Energy 83 (2009) 1086–1092

transparent insulation cover. International Journal of Energy Research Wang, S.A., Wu, L.G., 1990. Analysis and performance of flat-plate solar
23, 925–940. collector arrays. Solar Energy 45 (2), 71–78.
Tsilingris, P.T., 2002. Back absorbing parallel plate polymer absorbers in Weitbrecht, V., Lehmann, D., Richter, A., 2002. Flow distribution in
solar collector design. Energy Conversion and Management 43, solar collectors with laminar flow condition. Solar Energy 73,
135–150. 433–441.