# A Dynamic Multinode Model for Component-Oriented Thermal Analysis of Flat-Plate Solar Collectors

Authors: Christoph N. Reiter; Christoph Trinkl; Wilfried Zörner; Vic I. Hanby
Journal: Journal of Solar Energy
 (2015)

---

## Abstract

A mathematical model of a flat-plate solar collector was developed on the basis of the physical principles of optics and heat transfer in order to determine collector’s component temperatures as well as collector efficiency. In contrast to many available models, the targeted use of this dynamic model is the detailed, theoretical investigation of the thermal behaviour of newly developed or adjusted collector designs on component level, for example, absorber, casing, or transparent cover. The defined model is based on a multinode network (absorber, fluid, glazing, and backside insulation) containing the relevant physical equations to transfer the energy. The heat transfer network covers heat conduction, convection, and radiation. Furthermore, the collector optics is defined for the plane glazing and the absorber surface and also considers interactions between them. The model enables the variation of physical properties considering the geometric parameters and materials. Finally, the model was validated using measurement data and existing efficiency curve models. Both comparisons proved high accuracy of the developed model with deviation of up to 3% in collector efficiency and 1 K in component temperatures.


---

## Body

## 1. Introduction

For a detailed analysis and evaluation of flat-plate collector designs, specific information about temperature distribution in the collector on a component basis is necessary. This is especially relevant for overheating protection and material selection. Therefore, a high-integrity simulation model which is flexible concerning the implementation of component geometries, material properties, additional components, or functions is required. The model should be able to predict the thermal behaviour and temperature loads of the collector on a component level, for example, absorber, casing, or transparent cover. It should be used for theoretical analyses comprising stagnation tests, efficiency curve generations, and system simulations. This type of model was, for instance, applied by Reiter et al. [1, 2] to derive data about temperature loads and solar yield for the development and evaluation of overheating protection measures and polymeric components in collectors.Available parametric collector models are based on collector test results from solar simulator or outdoor testing. The function of the collector efficiency curve from the test results is derived by quadratic regression [3, 4]. Equation (1) describes the simple parameter model. Consider (1)
Q
˙
f
=
A
G
s
u
n
F
′
τ
α
-
a
T
f
-
T
a
m
b
-
b
T
f
-
T
a
m
b
2
. The linear heat loss coefficient a and the quadratic heat loss coefficient b characterise the thermal collector’s behaviour depending on the temperature difference between fluid and ambient air. The fraction of the converted solar energy by the absorber is limited by the so-called optical efficiency of the collector. When the mean fluid temperature is equal to the ambient temperature, almost no heat losses occur at the collector. Thus, the thermal output of the collector can only be reduced by the optical characteristics of the glazing and the absorber τ
α as well as by the collector efficiency factor F
′ caused by a temperature gradient in the fins, which results in heat losses.This one-dimensional approach can be used as a multinode model to reproduce the temperature gradient of the fluid along the absorber pipe between inlet and outlet. Hafner et al. [5] enhanced the basic equation of the parameter model by considering the heat capacity of the collector and heat losses dependent on the wind speed as well as additional radiative heat losses. However, these models are not suitable for the investigation of new collector designs. On the one hand, the parameters describe a certain measured collector but are not able to represent adjustments to the collector build-up or the implementation of overheating protection measures. On the other hand, the model only provides a fluid temperature. However, the detailed thermal investigation of collectors requires a broad output of the component temperatures.The steady-state model of Hottel and Whillier [6] only provides information on the collector performance by means of the fluid inlet temperature T
f
,
i
n. Equation (2) describes the collector model. Consider(2)
Q
˙
f
=
F
R
A
G
s
u
n
τ
α
-
U
c
o
l
l
T
f
,
i
n
-
T
a
m
b
. The optical characteristics τ
α, the heat removal factor F
R, and the overall heat loss coefficient of the collector U
c
o
l
l describe the performance of the collector. These characteristics are based on the physics of the collector considering dimensions, shapes, and materials. The temperature dependent heat transfer coefficients of the overall heat loss coefficient U
c
o
l
l vary subject to the conditions. Therefore, the component temperatures and the heat transfer coefficients have to be calculated by iteration. The reproduced physics of the model enable a change into various designs or build-ups. Klein et al. [7] investigated the transient behaviour of collectors or rather systems and identified significant differences to steady-state model approaches under certain conditions. During daily morning heating and rapidly changing ambient conditions like shading by clouds, neglecting the thermal capacity of the components led to different results. Therefore, a model on the basis of (2) was extended with two heat capacities. The fluid, absorber, and half of the insulation are at almost the same temperature and so were lumped together as the first heat capacity. The second heat capacity was the glazing. Klein et al. [7] and Wijeysundera [8] found that the difference between one-node models, two-node models, and experimental data for normal collector build-ups and systems according to the yield became progressively less. However, there are significant differences between steady-state and transient models, in particular with regard to component temperatures. The main disadvantage of the model is the consideration of average component temperatures for the entire collector. Hence, the increase in temperature along the risers during operation causing maximum and minimum temperatures is not taken into account. The collector properties, too, are considered as homogenous for the entire model. This makes the approach incapable of considering fractional changes (e.g., at the upper half) of the properties or of the design in the collector model.Thus, a dynamic, numerical model of a state-of-the-art flat-plate collector was developed under consideration of the relevant capacities deriving information on individual component temperatures, heat fluxes, and collector performance. Moreover, the model can be adjusted to new collector parameters or designs. By implementing the multinode collector model, the capacity of absorber, fluid, glazing, and casing was taken into account. The network connecting the nodes and the ambient are physical functions depending on geometry, size, materials, and operating conditions.


## 2. Framework and Assumptions

Modelling technical components requires a number of simplifying assumptions to limit the complexity of the model. Some minor effects in components are difficult to implement into the model and can be neglected due to their limited influence on the entire system. Thus, the following assumptions described for collector modelling were made. These assumptions, however, do not alter the basic physical situation of the collector:(1)
Heat transfer below and above the absorber is one-dimensional inz-direction (Figure 1).
(2)
Headers and their plates make up only a small area of the collector and can be neglected.
(3)
There is no difference in the fluid flow rates in the riser pipes.
(4)
The temperature gradient of the absorber plate (x-direction) caused by the fin efficiency of the plate-pipe design cannot be described by any transient model with a limited number of nodes. Moreover, a one-node model with an average temperature is not sufficient since the effects of the heat losses via convection and radiation are nonlinear. This is why the absorber plate has to be reproduced by a multinode model enabling a nonlinear temperature progression. The corresponding schematic temperature gradients are shown in Figure 2.
(5)
There are no temperature gradients around the perimeter of the riser pipes.
(6)
The temperature gradient of the absorber along the riser pipes (y-direction) is nonlinear and cannot be described by any transient model with a limited number of nodes. A one-node model with an average temperature is not sufficient since the maximum temperatures are lowered and the effects of the heat losses via convection and radiation are nonlinear. This has an impact on the other temperature nodes of the absorber as well. Therefore, a multinode model along the fluid flow is necessary. The corresponding schematic temperature gradients are shown in Figure 3.
(7)
The riser pipe is only connected to the fluid and the absorber nodes. The heat transfer to the air gap can be neglected since the pipe surface is relatively small. Furthermore, the temperature gradient between pipe and air is low.
(8)
The two harps of the absorber do not influence each other thermally. Hence, there is no heat transfer between neighbouring risers of the first and second harp, and neither is there any convection between the two arrays with different temperature.
(9)
Unsteady heat conduction in the insulation can be neglected. The analytical temperature distribution in a solid material during time-dependent temperature changes is influenced by the heat capacity of the system. Thus, a multinode model can represent this effect. However, the required number of nodes of the model would be significantly increased. Therefore, the heat capacity of the insulation is represented by a single node. In a second step, the linear temperature gradient through the insulation is calculated under steady-state conditions. Figure4 shows the described approaches.
(10)
The heat capacity of the backside sheet of the casing can be neglected due to the low mass of the component as well as the low temperature changes near the ambient air.
(11)
The frame of the collector is neglected. Due to the one-dimensional heat transfer through the collector, the effects in the frame cannot be taken into consideration. Furthermore, the effect of the heat capacity has only a slight influence since the temperature changes outside the collector insulation are low.
(12)
The temperature gradient in the glazing is neglected. The low thickness and the moderate heat conductivity of the glazing limit the influence of any heat conduction to a minimum.
(13)
The glazing is nontransparent for the wavelength of infrared radiation. Thus, heat radiation is completely absorbed by the glazing material. This assumption is valid for solar glass. However, the implementation of polymeric materials being partly transparent for infrared radiation would require an adjustment of the heat transfer network in this regard [9].
(14)
Since the impacts of the structured surface of the glazing are only marginal compared to the plain surfaces, this effect was not taken into consideration.
(15)
The sky is a black body for infrared radiation.
(16)
The ambient temperature around the collector is homogeneous.
(17)
The effect of wind is considered homogenously over the glazing inx-direction.
(18)
Wind at the backside is neglected because of the small distance between backside insulation and roof as well as the small temperature difference between outer surface and ambient air.
(19)
Shading from the frame at the absorber is not taken into consideration.
(20)
Effects from operation (dust, dirt, degradation, moisture…) influencing the optical and thermal performance of the collector are not taken into consideration.The implementation of the collector model was done in theMatlab/Simulink environment using its solver ODE45 [10]. This was supplemented by theSimulink-basedCARNOT blockset [5]. The latter contains models for the simulation of various heat sources, heat consumers, storages, controllers, pipes, and pumps as well as models for providing weather data and fluid properties. The collector physics was implemented in the simulation environment by means of a so-called “S-function (level 2),” a compiled source code written in the programming language C. The S-function is equipped with connections for data input for fluid, weather, and collector orientation which is provided byCARNOT blocks. The relevant output data of the fluid at the collector outlet as well as the component temperatures are provided at the output of the S-function (cf. [1]).Figure 1
Schematic build-up of an absorber stripe.Figure 2
Schematic temperature distribution at the cross section of an absorber in a one-node (a) and a seven-node (b) absorber model in comparison to the analytical temperature gradient.
(a)
(b)Figure 3
Schematic temperature distribution at the cross section of an absorber along the pipe in a one-node (a) and a five-node (b) absorber model in comparison to the analytical temperature gradient.
(a)
(b)Figure 4
Schematic temperature distribution at the cross section of the casing or insulation during unsteady heat conduction in an analytical model (a) and a one-node model with additional steady-state calculation of the temperature gradient.
(a)
(b)


## 3. Model Build-Up and Energy Balance

The developed numerical model of the collector is limited to the area of the risers. Using the symmetry of the collector, the model can be described by only one riser. This model of a riser is used for both halves of the double harp absorber. Figure5 shows the cross section of the riser or rather the collector model with its capacities and its connecting energy fluxes. The ten-node model can be reduced by using the symmetry of the absorber plate to a seven-node model.Figure 5
Cross section of the collector model with nodes and energy fluxes.The following gives an overview of the energy balance drawn in the collector model describing the capacities of each node as well as the energy fluxes between the nodes:Glazing (g) is as follows:(3)
m
g
c
g
d
T
g
d
t
=
S
g
+
Q
˙
a
1
-
g
+
2
∑
a
=
2
4
Q
˙
a
-
g
-
Q
˙
g
-
a
m
b
.
Absorber 1 (a
1) is as follows:(4)
m
a
1
c
a
1
d
T
a
1
d
t
=
S
a
1
+
2
Q
˙
a
2
-
a
1
-
Q
˙
a
1
-
f
-
Q
˙
a
1
-
g
-
Q
˙
a
-
e
d
g
e
.
Absorber 2 (a
2) is as follows:(5)
m
a
2
c
a
2
d
T
a
2
d
t
=
S
a
2
+
Q
˙
a
3
-
a
2
-
Q
˙
a
2
-
a
1
-
Q
˙
a
2
-
g
-
Q
˙
a
2
-
c
-
Q
˙
a
-
e
d
g
e
.
Absorber 3 (a
3) is as follows:(6)
m
a
3
c
a
3
d
T
a
3
d
t
=
S
a
3
+
Q
˙
a
4
-
a
3
-
Q
˙
a
3
-
a
2
-
Q
˙
a
3
-
g
-
Q
˙
a
3
-
c
-
Q
˙
a
-
e
d
g
e
.
Absorber 4 (a
4) is as follows: (7)
m
a
4
c
a
4
d
T
a
4
d
t
=
S
a
4
-
Q
˙
a
4
-
a
3
-
Q
˙
a
4
-
g
-
Q
˙
a
4
-
c
-
Q
˙
a
-
e
d
g
e
.
Fluid (f) is as follows:(8)
m
f
c
f
d
T
f
d
t
=
m
˙
f
c
f
∗
T
i
n
-
T
f
+
Q
˙
a
1
-
f
.
Casing (c) is as follows:(9)
m
c
c
c
d
T
c
d
t
=
Q
˙
a
1
-
c
+
2
Q
˙
a
2
-
c
+
2
Q
˙
a
3
-
c
+
2
Q
˙
a
4
-
c
-
Q
˙
c
-
a
m
b
.Node a
1 is the base of the fin and includes the riser pipe. The listed energy fluxes connecting the nodes to each other are described in the following in detail.


### 3.1. Optical Properties

The optical properties of the collector comprise the behaviour of the absorber as well as the glazing under irradiance. The optical effects determining the collector performance are reflection, absorption, and transmission (Figure6). The effects depend on the irradiance and the incidence angle as well as the thickness, the refractive index, and the extinction coefficient of the materials. The refractive index n and the extinction coefficient K of the glazing material depend on the wavelength of the incoming radiation. In the simulation model, these properties are assumed to be constant for the irradiance. Duffie and Beckman [12] considered this assumption to be ideally suited for covers made of glass.Figure 6
Basic optical effects at glazing (a) and absorber (b).
(a)
(b)


#### 3.1.1. Refraction Factor of the Glazing

The incident unpolarised solar radiation is polarised at the boundary surface between air and glazing. The radiation is divided into a parallel and a vertical component with regard to the plane of incidence (Figure7). The definition of the polarised light is necessary for the calculation of the associated reflection factor r [11].Figure 7
Polarisation directions parallel and vertical to the plane of incidence (according to [11]).The reflection factor of the boundary surface of the glazing is defined as the ratio between the reflected radiationG
r and the incident radiation G. This ratio results from the average reflection factor from the perpendicular and the parallel part (equation (10) [12]):(10)
r
=
G
r
G
=
r
⊥
+
r
∥
2
.The Fresnel formulae for smooth boundary surfaces between transparent materials with the refractive indices n
1 and n
2 are calculated according to the incidence angle θ
1 and the refraction angle θ
2 as shown in Figure 8.Figure 8
Relations between beam angle and refractive indices (according to [11, 12]).The reflection factors for vertical and parallel polarised radiation are shown in (11) [12]:(11)
r
⊥
=
sin
2
⁡
θ
2
-
θ
1
sin
2
⁡
θ
2
+
θ
1
r
∥
=
tan
2
⁡
θ
2
-
θ
1
tan
2
⁡
θ
2
+
θ
1
.The refraction angle θ
2 changes with the incidence angle θ
1 depending on the two refraction indices n
1 and n
2. The relation is described by Snell’s law in (12) [12]:(12)
n
1
sin
⁡
θ
1
=
n
2
sin
⁡
θ
2
.Vertical incidence at the boundary surface is a special case. Both angles θ
1 and θ
2 are zero. Also, the vertical and parallel refraction factors are equal. Hence, the refraction factor can be derived from (13) [12]:(13)
r
0
°
=
G
r
G
=
n
1
-
n
2
n
1
+
n
2
2
.


#### 3.1.2. Absorption in the Glazing

The absorption of the radiation in the glazing caused by a partly transparent material is described by Bouguer’s law. The underlying, basic assumption is that the absorbed radiationd
G is proportional to the incoming radiation G and the extinction coefficient K of the material. This relation is shown in (14) [12]:(14)
d
G
=
-
G
K
d
x
.The transmittance without reflection τ
a is described by means of the factor of radiation passing the material. This factor is the ratio between transmitted radiation G
t and incoming radiation G. This ratio is described by the following depending on the extinction coefficient K of the material and the path length of the radiation [12]:(15)
τ
a
=
G
t
G
=
exp
⁡
-
K
d
g
cos
⁡
θ
2
.


#### 3.1.3. Transmission and Absorption Coefficient of the Glazing

Considering both reflection and absorption in a glazing, the transmission as well as the absorption coefficient of the component can be derived. The boundary surfaces of the glazing transmit and reflect the incoming beam. Thus, there are an infinite number of transmissions, reflections, and absorptions in the glazing [12]. Figure 9 shows this effect in the glazing.Figure 9
Transmission, reflection, and absorption of a single glazing (according to [12]).The transmission coefficientτ of the glazing is the average of the transmission coefficients for both polarised directions (16). The infinite series of beams caused by multiple reflections can be determined via ray tracing described in (17) and (18) [12]:(16)
τ
=
τ
⊥
+
τ
∥
2
(17)
τ
⊥
=
1
-
r
⊥
2
τ
a
∑
n
=
0
∞
r
⊥
2
n
τ
a
2
n
=
τ
a
1
-
r
⊥
2
1
-
r
⊥
τ
a
2
(18)
τ
∥
=
τ
a
1
-
r
∥
2
1
-
r
∥
τ
a
2
.The absorption coefficient α of the glazing is also calculated for both polarised directions by an infinite series described in (19) [12]:(19)
α
=
α
⊥
+
α
∥
2
α
⊥
=
1
-
r
⊥
τ
a
∑
n
=
0
∞
1
+
r
⊥
n
τ
a
n
=
1
-
τ
a
1
-
r
⊥
1
-
r
⊥
τ
a
α
∥
=
1
-
τ
a
1
-
r
∥
1
-
r
∥
τ
a
.


#### 3.1.4. Solar Energy Conversion in the Collector

Equation (20) describes the energy input at the absorber reduced by the effective transmission absorption product τ
α
e
f
f of glazing and absorber [12]. Therefore,(20)
S
a
1
=
τ
α
e
f
f
A
G
.Simple approaches of (
τ
α
) only consider the light transmission through the glazing and the absorption of the absorber surface. In this model, however, a more detailed approach is applied as described by (21) [12]:(21)
τ
α
e
f
f
=
τ
α
1
-
1
-
α
ρ
.The multiple reflections between absorber and glazing are considered with the nonideal absorption coefficient α of the absorber and the reflection coefficient ρ of the inner glazing surface. Figure 10 shows this effect in the schematic collector build-up.Figure 10
Effective transmission absorption productτ
α
e
f
f of a single-glazed collector (according to [12]).The additional solar energy input in the glazing depending on absorption of the glazing changes the temperature of this node and influences the collectors’ heat losses [12]. The solar energy input in the glazing S
g is shown in the following:(22)
S
g
=
1
-
τ
a
A
G
.The optical effects for solar energy conversion in absorber and glazing depend on the angle. Hence, these effects have to be considered for direct as well as for diffuse irradiance: (23)
S
=
S
d
i
r
+
S
d
i
f
f
.For the calculation of the optical properties for direct irradiance, the beam angle at the inclined collector surface is used. The diffuse irradiance from the sky is assumed to be isotropic. Depending on the collector slope β, a certain fraction of the sky appears for the glazing as well as for the absorber. Brandemuehl and Beckman [13] investigated the relation between the sky and the inclined glazing systems for diffuse irradiation. As a result, an approach for the calculation of an equivalent angle depending on the collector slope was derived.Equation (24) shows the equivalent angle of incidence for diffuse irradiance for single-glazed collectors [13]. Consider(24)
θ
d
i
f
f
=
59.7
°
-
0.1388
β
+
0.001497
β
2
.


### 3.2. Convection

Heat transfer by convection occurs at several layers in the collector. There is forced convection in the absorber pipes caused by the pumped heat transfer fluid and at the surface of the glazing caused by wind. Free convection occurs in the air gaps in front of the absorber and behind the absorber as well as at the glazing and at the backside of the collector. These heat transfer mechanisms are represented by (25) containing the convective heat transfer coefficient from (26) that depends on the Nusselt number Nu and the characteristic length L of the component as well as on the heat conductivity k of the fluid [14]:(25)
Q
˙
=
h
c
o
n
v
A
T
1
-
T
2
(26)
h
c
o
n
v
=
N
u
k
L
.The convective heat transfer of solid bodies cannot be defined by physical equations, because of its strong dependency on fluid properties, shape, dimension, and temperature difference. Therefore, the heat transfer was investigated in experiments. Empirical correlations were derived from the measurement results and were described by the dimensionless heat transfer coefficient N
u. This approach defines the average convective heat transfer of the investigated setup and is valid for a certain range of fluid properties and dimensions. The average heat transfer coefficient usually referred to a uniform temperature of the body surface. Thus, average surface temperatures and average fluid properties are assumed in the constant temperature case. The constant heat flux case considering a nonuniform surface temperature is not considered due to the very small error in convective heat transfer of approximately ~0.1% [14].


#### 3.2.1. Forced Convection in the Fluid Channels

The forced convection in the absorber pipes is calculated according to the approach recommended by Gnielinski [15]. The characteristic length L is the pipe diameter d in this case. The flow in the pipes was assumed to be constant along the entire length without entrance effects. The temperature dependent properties of the fluid were referenced to the fluid node in the centre of the pipe.The hydrodynamically developed laminar flow for Reynolds numbersR
e below 2,300 and for a constant axial wall temperature is defined with a constant Nusselt number of 3.657 [16].In the transition region between laminar and fully developed turbulent flow with Reynolds numbers in the range of 2,300 and 10,000, the Nusselt number is calculated by the following for a constant axial wall temperature:(27)
N
u
f
=
1
-
γ
N
u
l
a
m
,
2300
+
γ
N
u
t
u
r
b
,
10
4
γ
=
R
e
-
2
,
300
10
4
-
2
,
300
f
o
r
0
≤
γ
≤
1
N
u
l
a
m
,
2300
=
49.371
+
N
u
l
a
m
,
2,2300
-
0.7
3
+
N
u
l
a
m
,
3,2300
3
3
N
u
l
a
m
,
2,2300
=
1.615
2
,
300
P
r
d
l
3
N
u
l
a
m
,
3,2300
=
2
1
+
22
P
r
6
2
,
300
P
r
d
l
N
u
t
u
r
b
,
10
4
=
0.0308
/
8
10
4
P
r
1
+
12.7
0.0308
/
8
P
r
2
/
3
-
1
1
+
d
l
2
/
3
.The Nusselt number for fully developed turbulent flow by Gnielinski [17] for Reynolds numbers larger than 10,000 is described by (28). The friction factor f can be determined by (29). Therefore,(28)
N
u
f
=
f
/
8
R
e
P
r
1
+
12.7
f
/
8
P
r
2
/
3
-
1
1
+
d
l
2
/
3
(29)
f
=
1.8
l
o
g
10
R
e
-
1.5
-
2
.These approaches are valid for a ratio between the pipe diameter d and the pipe length l less than or equal to unity as well as for Prandtl numbers P
r between 0.6 and 1,000.


#### 3.2.2. Convection at the Glazing

The convection on the outer surface of the glazing is a combination of free and forced convection. Based on experiments, Churchill [18] developed an approach for an overlay of both convections as shown in the following:(30)
N
u
g
-
a
m
b
=
N
u
f
r
e
e
,
g
-
a
m
b
3
+
N
u
f
o
r
c
e
d
,
g
-
a
m
b
3
3
.The Nusselt number for forced convection with turbulent boundary layer is considered according to (31) by Schlichting [19]:(31)
N
u
f
o
r
c
e
d
,
g
-
a
m
b
=
0.037
R
e
0.8
P
r
1
+
2.443
R
e
-
0.1
P
r
2
/
3
-
1
.Therefore, the air flow caused by the wind is assumed to occur along the collector width of the glazing according to Figure 11.Figure 11
Wind at the inclined collector glazing.The average Reynolds number in (32) is calculated by the wind velocity u
w
i
n
d, the characteristic length of the glazing L
g
-
a
m
b, and the kinematic viscosity ν
a
i
r. The equation is valid for Reynolds numbers between 5 × 105 and 1 × 107 as well as for Prandtl numbers between 0.6 and 2,000. Even during low temperatures and wind velocities, the Reynolds numbers are in the range of the validity. Hence, laminar flow for forced convection is neglected. Consider(32)
R
e
=
u
w
i
n
d
L
g
-
a
m
b
ν
a
i
r
.For free laminar convection at the glazing, two cases are defined by Fujii and Imura [20]. The most important case for thermal collector simulation is the heat transfer upwards at the inclined glazing. During collector operation, the glazing temperature turns out to be larger than the ambient temperature. Hence, the collector’s heat losses are driven by the upward flowing air. The second case with lower glazing temperatures is only relevant for system simulations with fluid temperatures below the ambient air temperatures. Figure 12 shows these two cases for the glazing.Figure 12
Free convection at the glazing with heat transfer upwards (a) and downwards (b).
(a)
(b)The Nusselt number for free laminar convection at the glazing with heat transfer upwards or downwards at vertical planes is described in (33) [21] and is also valid for Rayleigh numbers between 1 × 10−1 and 1 × 1012. Using the correction factor f
P
r in (34) [21] allows Prandtl numbers in the range between 0.001 and nonfinite. This correlation for the Nusselt numbers is also valid for inclined surfaces by substituting the Rayleigh number R
a for the adjusted Rayleigh number R
a
β according to (35) [20]. This adjustment considers the reduced component of gravity at the inclined plane instead of the acceleration due to gravity g. In this case, the characteristic length L
g
-
a
m
b is the length of the glazing. Consider(33)
N
u
f
r
e
e
,
g
-
a
m
b
=
0.825
+
0.387
R
a
f
P
r
6
2
f
o
r
R
a
<
R
a
c
r
i
t
(34)
f
P
r
=
1
+
0.492
P
r
9
/
16
-
16
/
9
f
o
r
0.001
<
P
r
<
∞
(35)
R
a
β
=
R
a
sin
⁡
β
.Boundary layer separation may occur for free convection at the glazing with heat transfer upwards. Hence, there is a change from laminar to turbulent flow. The laminar flow occurs for Rayleigh numbers below the critical Rayleigh number R
a
c
r
i
t depending on the collector slope β from (36). Figure 13 shows the gradient of the upper limit of the critical Rayleigh number for the described laminar convection. Consider(36)
R
a
c
r
i
t
=
10
8.9
-
0.00178
90
°
-
β
1.82
.Above the critical Rayleigh number, free turbulent convection at the glazing occurs. Fujii and Imura [20] defined the following accordingly:(37)
N
u
t
u
r
b
,
g
-
a
m
b
=
0.56
R
a
c
r
i
t
cos
⁡
90
°
-
β
0.25
+
0.13
R
a
3
-
R
a
c
r
i
t
3
.Figure 13
Upper limit of the critical Rayleigh number depending on the collector slope.


#### 3.2.3. Convection at the Casing

Duffie and Beckman [12] regarded the heat transfer from the casing to the ambient by radiation and convection as negligible because of the significant impact on the backside insulation and the resulting low temperature. In this model, the heat transfer coefficients from the backside to the ambient are calculated for a hypothetical windless region between collector and roof by natural convection and radiation. The effect of a chimney was also neglected. The heat transfer via convection was assumed as free convection at an inclined plane. Figure 14 describes the two occurring cases for solar collectors.Figure 14
Free convection at the casing with heat transfer downwards (a) and upwards (b).
(a)
(b)The most common case for the backside sheet is a heated plate transmitting heat downwards without any separation of the boundary layer. Thus, the Nusselt number for laminar flow from (33) to (35) serves as an exploration for the free convection behind the collector. In this case, L
c
-
a
m
b is the characteristic length of the casing.The special case of a cooled surface absorbing heat from the ambient air from below is described by (33) to (35) for laminar flow. Above the critical Rayleigh number, the separation of the boundary layer is expressed by means of (36). Furthermore, the approach for turbulent flow from (37) has to be used.


#### 3.2.4. Convection between Absorber and Glazing

Heat transfer in an enclosed air gap is driven by heat conduction as well as by free convection of the circulating gas layer. This free convection between glazing and absorber with the two possible directions of heat transfer is shown in Figure15.Figure 15
Free convection in the enclosed gas layer between absorber and glazing with heat transfer upwards (a) and downwards (b).
(a)
(b)Upwards heat transfer usually appears in collectors between absorber and glazing. Therefore, Matuska and Zmrhal [22] derived a Nusselt correlation from several correlations similar to Hollands et al. [23] to verify the Nusselt number approach for an increased range of collector slopes. Equation (38) describes the Nusselt number for the air gap valid for slopes β from 0° to 90°. The characteristic length L
a
-
g in that system is the distance between absorber and glazing. The value of f
β in (39) constitutes an angle factor for the Nusselt number that depends on the collector slope which was obtained as quadratic function in the derived correlation. The approach is valid for 1 × 104 < Ra < 2 × 106 covering the necessary operational range of solar-thermal collectors. The correlation, too, showed good results independent of the aspect ratio between air gap length and air gap distance. Several approaches, however, are limited with regard to the collector geometry. Consider(38)
N
u
a
-
g
=
f
β
R
a
0.29
(39)
f
β
=
0.1464
-
2.602
e
-
4
β
-
2.046
e
-
6
β
2
.The approach for the Nusselt number with heat transfer downwards is derived from Arnold et al. [24]. Therefore, (40) is used on the basis of the calculated Nusselt number for heat transfer in vertical air layers (β
=
90
°) derived from (38) and (39). The validity of this approach is similar to that of the approach of (38). However, the aspect ratio between air gap length and air gap distance has to be larger than 20. Consider(40)
N
u
a
-
g
=
1
+
N
u
a
-
g
β
=
90
°
-
1
sin
⁡
180
°
-
β
.


#### 3.2.5. Convection between Absorber and Casing

The heat transfer in the enclosed air layer between absorber plate and casing is described by the Nusselt correlation for convection between absorber and glazing. The usual approach for the region behind the absorber is downwards free convection expressed by (40).


### 3.3. Radiation

Heat transfer caused by radiation is based on the Stefan-Boltzmann law of a black body radiator. The energy transfer is described with an emission coefficientε
1
=
1 (after Eicker [11]) by (41). The Stefan-Boltzmann constant σ is shown in (42). Therefore,(41)
Q
˙
1
=
σ
A
1
T
1
4
(42)
σ
=
5.67
e
-
8
W
m
-
2
K
-
4
.Bodies above a temperature of 0 K emit radiating power. Radiating power from other bodies placed in a room is also absorbed. The position of the bodies also has great influence on the behaviour as described by their form factors Φ. The form factor is the ratio between the radiation receiving surface and the entire hemisphere. Figure 16 shows the heat exchange of the two surfaces.Figure 16
Heat fluxes via radiation between two elements (according to [11]).The net radiation exchangeQ
˙
r
a
d between two bodies takes the properties of the bodies into account and defines the heat transfer. Therefore, the difference between the radiating powers reaching the other body is taken from (43)
Q
˙
r
a
d
=
Q
˙
r
a
d
,
1
→
2
-
Q
˙
r
a
d
,
2
→
1
=
Φ
12
σ
A
1
T
1
4
-
Φ
21
σ
A
2
T
2
4
.The definition of the net radiation exchange between two bodies with an emission coefficient ε smaller than unity is a more complex approach. Therefore, a multiple of reflections has to be considered.The radiation exchange in solar-thermal flat-plate collectors embraces two special cases. Firstly, there is radiation exchange between two parallel plates, for example, between absorber and glazing; secondly, radiation occurs from a limited surface to the sky hemisphere.The heat transfer by radiation is determined according to (44). Therefore, the heat transfer coefficients for radiation of the collector have to be defined. Consider(44)
Q
˙
r
a
d
=
h
r
a
d
A
T
1
-
T
2
.


#### 3.3.1. Radiation between Two Parallel Plates

Due to the fact that the heat emitted is completely absorbed by the opposed area, the form factors of the two surfaces are equal (Φ
12
=
Φ
21
=
1). It is assumed that the effects of infinitely expanded parallel surfaces apply for heat transfer in solar collectors. Multiple reflections have to be considered for grey emitters as well by introducing their emission coefficients ε into the approach. The heat transfer coefficient for radiation between two parallel plates is given in the following:(45)
h
r
a
d
=
σ
1
/
ε
1
+
1
/
ε
2
-
1
T
1
2
+
T
2
2
T
1
+
T
2
.This radiation exchange is necessary for the setup of glazing and absorber, absorber and casing, and casing and roof.


#### 3.3.2. Radiation between Glazing and Sky

The radiation between glazing and cold sky is a special phenomenon of heat transfer [11]. The surface of glazing A
g is much smaller than the total sky hemisphere A
s
k
y. Hence, radiation from the glazing is completely emitted to the sky. Thus, the form factor of the glazing Φ
g
-
s
k
y is 1. In contrast to that, the form factor from the sky Φ
s
k
y
-
g is almost zero since the proportion of the glazing surface being seen from the sky hemisphere is very small. Therefore, the emitted heat from the sky is neglected because of the infinite area of the sky’s half-space. The heat transfer coefficient is shown in (46). The temperatures have to be inserted in Kelvin because radiation refers to the absolute temperature of the body or rather surface. Consider(46)
h
r
a
d
,
g
-
s
k
y
=
σ
ε
g
T
g
2
+
T
s
k
y
2
T
g
+
T
s
k
y
.The temperature of the sky differs from the ambient air temperature [12]. Approaches for the definition of the sky temperature are related to the air temperature [25], to the water vapour pressure [26], or to the dew point temperature [27], for example.


### 3.4. Heat Conduction

Heat conduction is prevalent at glazing, absorber, and casing. The following describes the heat transfer that depends on the heat conductivityk and the thickness or rather distance s of the component [14]:(47)
Q
˙
c
o
n
d
=
h
c
o
n
d
A
T
1
-
T
2
h
c
o
n
d
=
k
s
.Vertical heat transfer through thin plates in the thermal network like glazing or absorber can be neglected. Relevant heat conduction occurs between the nodes of the absorber perpendicular to the cross-sectional area of the plate and towards the riser pipe (Figure 17).Figure 17
Heat conduction between the absorber nodes.The second heat conduction is seen in the collector insulation. These heat losses occur at the side and backside insulation. In contrast to the massive solid materials such as metals, the heat conductivity cannot be regarded as constant. The cavernous structure of the material also enables heat transfer effects of the enclosed air and radiation. Thus, the heat conductivity increases according to the material temperature. For dry mineral wool, Ochs et al. [28] derived a linear approach for the temperature depending calculation of the heat conductivity k
m
w according to (48). The influence on the heat conductivity of the unknown humidity in the mineral wool is neglected. Moist mineral wool usually has an increased heat conductivity. Consider(48)
k
m
w
=
0.035
W
m
-
1
K
-
1
+
2.0
e
-
4
T
m
w
.


### 3.5. Edge Losses

The edge losses, which comprise the heat losses through the collector frame or rather side, cannot be mathematically calculated in detail. Tabor [29] describes an approach considering the losses only connected to the absorber node. The other nodes have an adiabatic behaviour against the edge. Equation (49) shows the heat transfer between absorber and ambient air as well as the heat transfer coefficient considering only the heat conduction through the side insulation. Consider(49)
Q
˙
a
-
e
d
g
e
=
h
e
d
g
e
A
e
d
g
e
T
a
-
T
a
m
b
h
e
d
g
e
=
k
c
s
e
d
g
e
.


### 3.6. Overall Heat Transfer Coefficients between the Nodes

The heat transfer coefficientsh of the collector connect the nodes and the ambient. Most of them are embedded in a network of the thermal resistances to calculate a heat flow rate Q
˙. The electrical analogy is used to generate the overall heat transfer coefficients between the nodes based on the single thermal effects [14, 16]. This analogy is used to generate overall resistances by combining single resistances of serial and parallel connections according to Figure 5. The required overall heat transfer coefficients (related to the area) of the energy fluxes from (3) to (9) are shown in the following:(50)
1
R
t
h
,
g
-
a
m
b
=
h
c
o
n
d
,
g
A
+
1
h
r
a
d
,
g
-
s
k
y
A
+
1
h
c
o
n
v
,
g
-
a
m
b
A
-
1
1
R
t
h
,
a
-
g
=
1
h
r
a
d
,
a
-
g
A
a
+
1
h
c
o
n
v
,
a
-
g
A
a
-
1
+
h
c
o
n
d
,
g
A
a
1
R
t
h
,
a
1
-
f
=
h
c
o
n
d
,
a
1
2
A
w
a
l
l
+
h
c
o
n
v
,
a
1
-
f
A
p
i
p
e
1
R
t
h
,
a
-
c
=
1
h
r
a
d
,
a
-
c
A
a
+
1
h
c
o
n
v
,
a
-
c
A
a
-
1
+
h
c
o
n
d
,
c
A
a
1
R
t
h
,
c
-
a
m
b
=
h
c
o
n
d
,
c
A
+
1
h
r
a
d
,
c
-
a
m
b
A
+
1
h
c
o
n
v
,
c
-
a
m
b
A
-
1
.


### 3.7. Validation

The model was validated using the parameter set of a tested conventional flat-plate collector. The validation of this model and physical effects is also adequate for models with implemented overheating protection measures or adjusted build-up since the defined physical effects are varied and not exchanged [1, 30].The model was tested in various settings to investigate the accuracy of the results. According to the above-mentioned simplifying assumptions in the model, there are, of course, imperfections in the build-up of the tested collectors. Hence, the differences in the results have to be identified, justified, and evaluated. Collector performance data is generally obtained under steady-state conditions; hence, the validation data presented focus on measurements and the performance of other models under equilibrium conditions.


#### 3.7.1. Comparison with an Existing Validated Model

In a first validation step, the developed dynamic multinode model was compared to the results of a one-dimensional multinode model approach by Hafner et al. [5] based on the approach of (1). This validation test is only able to compare the results of the collector efficiency of an already validated parameter model approach. The tested parameters of the collector are shown in Table 1.Table 1
Efficiency curve of the tested collector for the parameter model [31].
Efficiency curve parameters
Zero loss efficiency
η
0
0.798
[—]
Linear heat loss coefficient
a
3.34
[W m−2 K−1]
Quadratic heat loss coefficient
b
0.0075
[W m−2 K−2]
Heat capacity
C
9.5
[kJ K−1]The collector surface of both models was reduced to 1.83 m². Thus, both models consider the absorber surface without header regions.Two test simulations were set up to compare the results of the models. The first test comprised the simulation of the collector efficiency curve like the practical test modus in a solar simulator according to DIN EN 12975-2 [32]. Under constant ambient conditions (G
=
885
W
m
-
2, T
a
m
b
=
27
°
C, and u
w
i
n
d
=
2
m
s
-
1), a single collector was calculated. The heat carrier, a water-glycol mixture (33.3%), had a flow rate of m
˙
=
0.08
k
g
s
-
1 and variable temperatures through the collector.The percentage difference of the collector efficiency describes the accuracy of the model approach according to the following:(51)
Percentage difference
=
η
dynamic model
η
paramter model
-
1
100
%
. The values range between −1.0% and +3.3% for flow temperatures from 25°C to 90°C (Figure 18).Figure 18
Comparison of the simulated efficiency curves of the parameter model and the developed dynamic collector model.The differences are relatively low according to the assumptions in the model. Also, the simulated temperature range in the test covers a large operational range of the collector to prove the use of the model.The deviation is assumed to originate from the discretisation in the model. The major influence in the test is the transition from laminar to turbulent flow at higher values of the reduced temperature. The change in the model curve is remarkable because of an averaged calculation of fully developed flow conditions for the entire riser pipe. In contrast to that, the real collector shows a continuous development of flow conditions along the pipe length. There are also neglected entrance effects increased by the vertical connection between header and riser.Next to the analysis of efficiency, the models were compared to each other in an annual simulation of a solar-thermal heating system for domestic hot water and space heating. The deviation of the solar yield was found to be +2.9%, similar to the results found in the comparison of efficiency. The deviation originates from the fact that the major operational range of the collectors is at high reduced temperatures. This is caused by both high flow and ambient air temperatures for hot water preparation in summer as well as moderate flow and low ambient air temperatures in the heating period for space heating.


#### 3.7.2. Comparison with Measurements

In a second step, the simulation results were compared to measurement data from analysing an already exposed field-testing collector; relevant collector characteristics are given in Table2 (cf. [1]).Table 2
Characteristics of the field-testing collector.
Collector
Gross length
1.987
[m]
Gross width
1.046
[m]
Gross area
2.078
[m²]
Aperture area
1.903
[m²]
Weight (empty)
39
[kg]
Liquid content
1.12
[l]
Glazing
Thickness
0.0032
[m]
Transmission coefficientτ
0.918
[—]
Absorber
Absorber area
1.903
[m²]
Absorption coefficient
0.95
[—]
Emission coefficient
0.04
[—]
Number of risers
10
[—]The measurements were performed using Institute’s solar simulator. The procedures defined in the European testing standard DIN EN 12975-2 [32] were applied using the equipment postulated so that the results presented are in line with the standard. However, degradation of the absorber coating and humidity in the mineral wool can be named as some of the reasons, resulting in the effect that the efficiency of the collector is reduced in comparison to results from previous efficiency curves of a new collector. Therefore, the emission coefficient of the absorber coating and the heat conductivity of the mineral wool in the model were adjusted to represent the typical degradation of the components.The first test was the comparison of the efficiency curves. Therefore, the measured conditions of the solar simulator were reproduced (G
=
890
W
m
-
2, T
a
m
b
=
31
°
C, and u
w
i
n
d
=
1
m
s
-
1); the heat carrier was water (m
˙
=
0.038
k
g
s
-
1). Figure 19 shows the good coverage between measurement and simulation. According to the following, the percentage difference is only in the range between +0.4% and +3.1%:(52)
Percentage difference
=
η
dynamic model
η
measurement
-
1
100
%
.The second test also showed satisfying conformity in particular with regard to individual component temperatures. During stagnation at insolation of 895 W m−2 and an ambient air temperature of 33°C, the absorber temperature of model and test (average of 9 measurement points) differs only by +1 K. Also the average temperatures (average of 3 measurement points) at glazing and backside vary below +1 K (Table 3).Table 3
Collector component temperatures during stagnation of a conventional flat-plate collector in testing and simulation.
Testing
Model
Difference
Glazing
71°C
71°C
±0 K
Absorber
174°C
175°C
+1 K
Backside (outside)
54°C
55°C
+1 KFigure 19
Comparison of the measured and simulated efficiency of the collector.


#### 3.7.3. Evaluation of the Test Results

All in all, the model proves to be valid according to the test collector. Slight differences are identified as being within the theoretical limits of the model. Thanks to the assumptions that have been made, the model generates reliable results for flat-plate collectors in due consideration of the limits of the model. Comparable model approaches of state-of-the-art collectors show similar low uncertainties in simulation of efficiency or rather solar yield. The model of Koo [33] showed percentage differences between −2% and +3% in comparison to experimental data of the collector efficiency curve. Matuska and Zmrhal [22] described a percentage difference of about 3-4% in the model validation. Furthermore, Cadafalch [34] identified differences below 2%. Thus, there is confidence that the model can be used for the investigation of entire systems, for identifying the solar yield, and for the analysis of the component temperatures in the collector.


## 3.1. Optical Properties

The optical properties of the collector comprise the behaviour of the absorber as well as the glazing under irradiance. The optical effects determining the collector performance are reflection, absorption, and transmission (Figure6). The effects depend on the irradiance and the incidence angle as well as the thickness, the refractive index, and the extinction coefficient of the materials. The refractive index n and the extinction coefficient K of the glazing material depend on the wavelength of the incoming radiation. In the simulation model, these properties are assumed to be constant for the irradiance. Duffie and Beckman [12] considered this assumption to be ideally suited for covers made of glass.Figure 6
Basic optical effects at glazing (a) and absorber (b).
(a)
(b)


### 3.1.1. Refraction Factor of the Glazing

The incident unpolarised solar radiation is polarised at the boundary surface between air and glazing. The radiation is divided into a parallel and a vertical component with regard to the plane of incidence (Figure7). The definition of the polarised light is necessary for the calculation of the associated reflection factor r [11].Figure 7
Polarisation directions parallel and vertical to the plane of incidence (according to [11]).The reflection factor of the boundary surface of the glazing is defined as the ratio between the reflected radiationG
r and the incident radiation G. This ratio results from the average reflection factor from the perpendicular and the parallel part (equation (10) [12]):(10)
r
=
G
r
G
=
r
⊥
+
r
∥
2
.The Fresnel formulae for smooth boundary surfaces between transparent materials with the refractive indices n
1 and n
2 are calculated according to the incidence angle θ
1 and the refraction angle θ
2 as shown in Figure 8.Figure 8
Relations between beam angle and refractive indices (according to [11, 12]).The reflection factors for vertical and parallel polarised radiation are shown in (11) [12]:(11)
r
⊥
=
sin
2
⁡
θ
2
-
θ
1
sin
2
⁡
θ
2
+
θ
1
r
∥
=
tan
2
⁡
θ
2
-
θ
1
tan
2
⁡
θ
2
+
θ
1
.The refraction angle θ
2 changes with the incidence angle θ
1 depending on the two refraction indices n
1 and n
2. The relation is described by Snell’s law in (12) [12]:(12)
n
1
sin
⁡
θ
1
=
n
2
sin
⁡
θ
2
.Vertical incidence at the boundary surface is a special case. Both angles θ
1 and θ
2 are zero. Also, the vertical and parallel refraction factors are equal. Hence, the refraction factor can be derived from (13) [12]:(13)
r
0
°
=
G
r
G
=
n
1
-
n
2
n
1
+
n
2
2
.


### 3.1.2. Absorption in the Glazing

The absorption of the radiation in the glazing caused by a partly transparent material is described by Bouguer’s law. The underlying, basic assumption is that the absorbed radiationd
G is proportional to the incoming radiation G and the extinction coefficient K of the material. This relation is shown in (14) [12]:(14)
d
G
=
-
G
K
d
x
.The transmittance without reflection τ
a is described by means of the factor of radiation passing the material. This factor is the ratio between transmitted radiation G
t and incoming radiation G. This ratio is described by the following depending on the extinction coefficient K of the material and the path length of the radiation [12]:(15)
τ
a
=
G
t
G
=
exp
⁡
-
K
d
g
cos
⁡
θ
2
.


### 3.1.3. Transmission and Absorption Coefficient of the Glazing

Considering both reflection and absorption in a glazing, the transmission as well as the absorption coefficient of the component can be derived. The boundary surfaces of the glazing transmit and reflect the incoming beam. Thus, there are an infinite number of transmissions, reflections, and absorptions in the glazing [12]. Figure 9 shows this effect in the glazing.Figure 9
Transmission, reflection, and absorption of a single glazing (according to [12]).The transmission coefficientτ of the glazing is the average of the transmission coefficients for both polarised directions (16). The infinite series of beams caused by multiple reflections can be determined via ray tracing described in (17) and (18) [12]:(16)
τ
=
τ
⊥
+
τ
∥
2
(17)
τ
⊥
=
1
-
r
⊥
2
τ
a
∑
n
=
0
∞
r
⊥
2
n
τ
a
2
n
=
τ
a
1
-
r
⊥
2
1
-
r
⊥
τ
a
2
(18)
τ
∥
=
τ
a
1
-
r
∥
2
1
-
r
∥
τ
a
2
.The absorption coefficient α of the glazing is also calculated for both polarised directions by an infinite series described in (19) [12]:(19)
α
=
α
⊥
+
α
∥
2
α
⊥
=
1
-
r
⊥
τ
a
∑
n
=
0
∞
1
+
r
⊥
n
τ
a
n
=
1
-
τ
a
1
-
r
⊥
1
-
r
⊥
τ
a
α
∥
=
1
-
τ
a
1
-
r
∥
1
-
r
∥
τ
a
.


### 3.1.4. Solar Energy Conversion in the Collector

Equation (20) describes the energy input at the absorber reduced by the effective transmission absorption product τ
α
e
f
f of glazing and absorber [12]. Therefore,(20)
S
a
1
=
τ
α
e
f
f
A
G
.Simple approaches of (
τ
α
) only consider the light transmission through the glazing and the absorption of the absorber surface. In this model, however, a more detailed approach is applied as described by (21) [12]:(21)
τ
α
e
f
f
=
τ
α
1
-
1
-
α
ρ
.The multiple reflections between absorber and glazing are considered with the nonideal absorption coefficient α of the absorber and the reflection coefficient ρ of the inner glazing surface. Figure 10 shows this effect in the schematic collector build-up.Figure 10
Effective transmission absorption productτ
α
e
f
f of a single-glazed collector (according to [12]).The additional solar energy input in the glazing depending on absorption of the glazing changes the temperature of this node and influences the collectors’ heat losses [12]. The solar energy input in the glazing S
g is shown in the following:(22)
S
g
=
1
-
τ
a
A
G
.The optical effects for solar energy conversion in absorber and glazing depend on the angle. Hence, these effects have to be considered for direct as well as for diffuse irradiance: (23)
S
=
S
d
i
r
+
S
d
i
f
f
.For the calculation of the optical properties for direct irradiance, the beam angle at the inclined collector surface is used. The diffuse irradiance from the sky is assumed to be isotropic. Depending on the collector slope β, a certain fraction of the sky appears for the glazing as well as for the absorber. Brandemuehl and Beckman [13] investigated the relation between the sky and the inclined glazing systems for diffuse irradiation. As a result, an approach for the calculation of an equivalent angle depending on the collector slope was derived.Equation (24) shows the equivalent angle of incidence for diffuse irradiance for single-glazed collectors [13]. Consider(24)
θ
d
i
f
f
=
59.7
°
-
0.1388
β
+
0.001497
β
2
.


## 3.1.1. Refraction Factor of the Glazing

The incident unpolarised solar radiation is polarised at the boundary surface between air and glazing. The radiation is divided into a parallel and a vertical component with regard to the plane of incidence (Figure7). The definition of the polarised light is necessary for the calculation of the associated reflection factor r [11].Figure 7
Polarisation directions parallel and vertical to the plane of incidence (according to [11]).The reflection factor of the boundary surface of the glazing is defined as the ratio between the reflected radiationG
r and the incident radiation G. This ratio results from the average reflection factor from the perpendicular and the parallel part (equation (10) [12]):(10)
r
=
G
r
G
=
r
⊥
+
r
∥
2
.The Fresnel formulae for smooth boundary surfaces between transparent materials with the refractive indices n
1 and n
2 are calculated according to the incidence angle θ
1 and the refraction angle θ
2 as shown in Figure 8.Figure 8
Relations between beam angle and refractive indices (according to [11, 12]).The reflection factors for vertical and parallel polarised radiation are shown in (11) [12]:(11)
r
⊥
=
sin
2
⁡
θ
2
-
θ
1
sin
2
⁡
θ
2
+
θ
1
r
∥
=
tan
2
⁡
θ
2
-
θ
1
tan
2
⁡
θ
2
+
θ
1
.The refraction angle θ
2 changes with the incidence angle θ
1 depending on the two refraction indices n
1 and n
2. The relation is described by Snell’s law in (12) [12]:(12)
n
1
sin
⁡
θ
1
=
n
2
sin
⁡
θ
2
.Vertical incidence at the boundary surface is a special case. Both angles θ
1 and θ
2 are zero. Also, the vertical and parallel refraction factors are equal. Hence, the refraction factor can be derived from (13) [12]:(13)
r
0
°
=
G
r
G
=
n
1
-
n
2
n
1
+
n
2
2
.


## 3.1.2. Absorption in the Glazing

The absorption of the radiation in the glazing caused by a partly transparent material is described by Bouguer’s law. The underlying, basic assumption is that the absorbed radiationd
G is proportional to the incoming radiation G and the extinction coefficient K of the material. This relation is shown in (14) [12]:(14)
d
G
=
-
G
K
d
x
.The transmittance without reflection τ
a is described by means of the factor of radiation passing the material. This factor is the ratio between transmitted radiation G
t and incoming radiation G. This ratio is described by the following depending on the extinction coefficient K of the material and the path length of the radiation [12]:(15)
τ
a
=
G
t
G
=
exp
⁡
-
K
d
g
cos
⁡
θ
2
.


## 3.1.3. Transmission and Absorption Coefficient of the Glazing

Considering both reflection and absorption in a glazing, the transmission as well as the absorption coefficient of the component can be derived. The boundary surfaces of the glazing transmit and reflect the incoming beam. Thus, there are an infinite number of transmissions, reflections, and absorptions in the glazing [12]. Figure 9 shows this effect in the glazing.Figure 9
Transmission, reflection, and absorption of a single glazing (according to [12]).The transmission coefficientτ of the glazing is the average of the transmission coefficients for both polarised directions (16). The infinite series of beams caused by multiple reflections can be determined via ray tracing described in (17) and (18) [12]:(16)
τ
=
τ
⊥
+
τ
∥
2
(17)
τ
⊥
=
1
-
r
⊥
2
τ
a
∑
n
=
0
∞
r
⊥
2
n
τ
a
2
n
=
τ
a
1
-
r
⊥
2
1
-
r
⊥
τ
a
2
(18)
τ
∥
=
τ
a
1
-
r
∥
2
1
-
r
∥
τ
a
2
.The absorption coefficient α of the glazing is also calculated for both polarised directions by an infinite series described in (19) [12]:(19)
α
=
α
⊥
+
α
∥
2
α
⊥
=
1
-
r
⊥
τ
a
∑
n
=
0
∞
1
+
r
⊥
n
τ
a
n
=
1
-
τ
a
1
-
r
⊥
1
-
r
⊥
τ
a
α
∥
=
1
-
τ
a
1
-
r
∥
1
-
r
∥
τ
a
.


## 3.1.4. Solar Energy Conversion in the Collector

Equation (20) describes the energy input at the absorber reduced by the effective transmission absorption product τ
α
e
f
f of glazing and absorber [12]. Therefore,(20)
S
a
1
=
τ
α
e
f
f
A
G
.Simple approaches of (
τ
α
) only consider the light transmission through the glazing and the absorption of the absorber surface. In this model, however, a more detailed approach is applied as described by (21) [12]:(21)
τ
α
e
f
f
=
τ
α
1
-
1
-
α
ρ
.The multiple reflections between absorber and glazing are considered with the nonideal absorption coefficient α of the absorber and the reflection coefficient ρ of the inner glazing surface. Figure 10 shows this effect in the schematic collector build-up.Figure 10
Effective transmission absorption productτ
α
e
f
f of a single-glazed collector (according to [12]).The additional solar energy input in the glazing depending on absorption of the glazing changes the temperature of this node and influences the collectors’ heat losses [12]. The solar energy input in the glazing S
g is shown in the following:(22)
S
g
=
1
-
τ
a
A
G
.The optical effects for solar energy conversion in absorber and glazing depend on the angle. Hence, these effects have to be considered for direct as well as for diffuse irradiance: (23)
S
=
S
d
i
r
+
S
d
i
f
f
.For the calculation of the optical properties for direct irradiance, the beam angle at the inclined collector surface is used. The diffuse irradiance from the sky is assumed to be isotropic. Depending on the collector slope β, a certain fraction of the sky appears for the glazing as well as for the absorber. Brandemuehl and Beckman [13] investigated the relation between the sky and the inclined glazing systems for diffuse irradiation. As a result, an approach for the calculation of an equivalent angle depending on the collector slope was derived.Equation (24) shows the equivalent angle of incidence for diffuse irradiance for single-glazed collectors [13]. Consider(24)
θ
d
i
f
f
=
59.7
°
-
0.1388
β
+
0.001497
β
2
.


## 3.2. Convection

Heat transfer by convection occurs at several layers in the collector. There is forced convection in the absorber pipes caused by the pumped heat transfer fluid and at the surface of the glazing caused by wind. Free convection occurs in the air gaps in front of the absorber and behind the absorber as well as at the glazing and at the backside of the collector. These heat transfer mechanisms are represented by (25) containing the convective heat transfer coefficient from (26) that depends on the Nusselt number Nu and the characteristic length L of the component as well as on the heat conductivity k of the fluid [14]:(25)
Q
˙
=
h
c
o
n
v
A
T
1
-
T
2
(26)
h
c
o
n
v
=
N
u
k
L
.The convective heat transfer of solid bodies cannot be defined by physical equations, because of its strong dependency on fluid properties, shape, dimension, and temperature difference. Therefore, the heat transfer was investigated in experiments. Empirical correlations were derived from the measurement results and were described by the dimensionless heat transfer coefficient N
u. This approach defines the average convective heat transfer of the investigated setup and is valid for a certain range of fluid properties and dimensions. The average heat transfer coefficient usually referred to a uniform temperature of the body surface. Thus, average surface temperatures and average fluid properties are assumed in the constant temperature case. The constant heat flux case considering a nonuniform surface temperature is not considered due to the very small error in convective heat transfer of approximately ~0.1% [14].


### 3.2.1. Forced Convection in the Fluid Channels

The forced convection in the absorber pipes is calculated according to the approach recommended by Gnielinski [15]. The characteristic length L is the pipe diameter d in this case. The flow in the pipes was assumed to be constant along the entire length without entrance effects. The temperature dependent properties of the fluid were referenced to the fluid node in the centre of the pipe.The hydrodynamically developed laminar flow for Reynolds numbersR
e below 2,300 and for a constant axial wall temperature is defined with a constant Nusselt number of 3.657 [16].In the transition region between laminar and fully developed turbulent flow with Reynolds numbers in the range of 2,300 and 10,000, the Nusselt number is calculated by the following for a constant axial wall temperature:(27)
N
u
f
=
1
-
γ
N
u
l
a
m
,
2300
+
γ
N
u
t
u
r
b
,
10
4
γ
=
R
e
-
2
,
300
10
4
-
2
,
300
f
o
r
0
≤
γ
≤
1
N
u
l
a
m
,
2300
=
49.371
+
N
u
l
a
m
,
2,2300
-
0.7
3
+
N
u
l
a
m
,
3,2300
3
3
N
u
l
a
m
,
2,2300
=
1.615
2
,
300
P
r
d
l
3
N
u
l
a
m
,
3,2300
=
2
1
+
22
P
r
6
2
,
300
P
r
d
l
N
u
t
u
r
b
,
10
4
=
0.0308
/
8
10
4
P
r
1
+
12.7
0.0308
/
8
P
r
2
/
3
-
1
1
+
d
l
2
/
3
.The Nusselt number for fully developed turbulent flow by Gnielinski [17] for Reynolds numbers larger than 10,000 is described by (28). The friction factor f can be determined by (29). Therefore,(28)
N
u
f
=
f
/
8
R
e
P
r
1
+
12.7
f
/
8
P
r
2
/
3
-
1
1
+
d
l
2
/
3
(29)
f
=
1.8
l
o
g
10
R
e
-
1.5
-
2
.These approaches are valid for a ratio between the pipe diameter d and the pipe length l less than or equal to unity as well as for Prandtl numbers P
r between 0.6 and 1,000.


### 3.2.2. Convection at the Glazing

The convection on the outer surface of the glazing is a combination of free and forced convection. Based on experiments, Churchill [18] developed an approach for an overlay of both convections as shown in the following:(30)
N
u
g
-
a
m
b
=
N
u
f
r
e
e
,
g
-
a
m
b
3
+
N
u
f
o
r
c
e
d
,
g
-
a
m
b
3
3
.The Nusselt number for forced convection with turbulent boundary layer is considered according to (31) by Schlichting [19]:(31)
N
u
f
o
r
c
e
d
,
g
-
a
m
b
=
0.037
R
e
0.8
P
r
1
+
2.443
R
e
-
0.1
P
r
2
/
3
-
1
.Therefore, the air flow caused by the wind is assumed to occur along the collector width of the glazing according to Figure 11.Figure 11
Wind at the inclined collector glazing.The average Reynolds number in (32) is calculated by the wind velocity u
w
i
n
d, the characteristic length of the glazing L
g
-
a
m
b, and the kinematic viscosity ν
a
i
r. The equation is valid for Reynolds numbers between 5 × 105 and 1 × 107 as well as for Prandtl numbers between 0.6 and 2,000. Even during low temperatures and wind velocities, the Reynolds numbers are in the range of the validity. Hence, laminar flow for forced convection is neglected. Consider(32)
R
e
=
u
w
i
n
d
L
g
-
a
m
b
ν
a
i
r
.For free laminar convection at the glazing, two cases are defined by Fujii and Imura [20]. The most important case for thermal collector simulation is the heat transfer upwards at the inclined glazing. During collector operation, the glazing temperature turns out to be larger than the ambient temperature. Hence, the collector’s heat losses are driven by the upward flowing air. The second case with lower glazing temperatures is only relevant for system simulations with fluid temperatures below the ambient air temperatures. Figure 12 shows these two cases for the glazing.Figure 12
Free convection at the glazing with heat transfer upwards (a) and downwards (b).
(a)
(b)The Nusselt number for free laminar convection at the glazing with heat transfer upwards or downwards at vertical planes is described in (33) [21] and is also valid for Rayleigh numbers between 1 × 10−1 and 1 × 1012. Using the correction factor f
P
r in (34) [21] allows Prandtl numbers in the range between 0.001 and nonfinite. This correlation for the Nusselt numbers is also valid for inclined surfaces by substituting the Rayleigh number R
a for the adjusted Rayleigh number R
a
β according to (35) [20]. This adjustment considers the reduced component of gravity at the inclined plane instead of the acceleration due to gravity g. In this case, the characteristic length L
g
-
a
m
b is the length of the glazing. Consider(33)
N
u
f
r
e
e
,
g
-
a
m
b
=
0.825
+
0.387
R
a
f
P
r
6
2
f
o
r
R
a
<
R
a
c
r
i
t
(34)
f
P
r
=
1
+
0.492
P
r
9
/
16
-
16
/
9
f
o
r
0.001
<
P
r
<
∞
(35)
R
a
β
=
R
a
sin
⁡
β
.Boundary layer separation may occur for free convection at the glazing with heat transfer upwards. Hence, there is a change from laminar to turbulent flow. The laminar flow occurs for Rayleigh numbers below the critical Rayleigh number R
a
c
r
i
t depending on the collector slope β from (36). Figure 13 shows the gradient of the upper limit of the critical Rayleigh number for the described laminar convection. Consider(36)
R
a
c
r
i
t
=
10
8.9
-
0.00178
90
°
-
β
1.82
.Above the critical Rayleigh number, free turbulent convection at the glazing occurs. Fujii and Imura [20] defined the following accordingly:(37)
N
u
t
u
r
b
,
g
-
a
m
b
=
0.56
R
a
c
r
i
t
cos
⁡
90
°
-
β
0.25
+
0.13
R
a
3
-
R
a
c
r
i
t
3
.Figure 13
Upper limit of the critical Rayleigh number depending on the collector slope.


### 3.2.3. Convection at the Casing

Duffie and Beckman [12] regarded the heat transfer from the casing to the ambient by radiation and convection as negligible because of the significant impact on the backside insulation and the resulting low temperature. In this model, the heat transfer coefficients from the backside to the ambient are calculated for a hypothetical windless region between collector and roof by natural convection and radiation. The effect of a chimney was also neglected. The heat transfer via convection was assumed as free convection at an inclined plane. Figure 14 describes the two occurring cases for solar collectors.Figure 14
Free convection at the casing with heat transfer downwards (a) and upwards (b).
(a)
(b)The most common case for the backside sheet is a heated plate transmitting heat downwards without any separation of the boundary layer. Thus, the Nusselt number for laminar flow from (33) to (35) serves as an exploration for the free convection behind the collector. In this case, L
c
-
a
m
b is the characteristic length of the casing.The special case of a cooled surface absorbing heat from the ambient air from below is described by (33) to (35) for laminar flow. Above the critical Rayleigh number, the separation of the boundary layer is expressed by means of (36). Furthermore, the approach for turbulent flow from (37) has to be used.


### 3.2.4. Convection between Absorber and Glazing

Heat transfer in an enclosed air gap is driven by heat conduction as well as by free convection of the circulating gas layer. This free convection between glazing and absorber with the two possible directions of heat transfer is shown in Figure15.Figure 15
Free convection in the enclosed gas layer between absorber and glazing with heat transfer upwards (a) and downwards (b).
(a)
(b)Upwards heat transfer usually appears in collectors between absorber and glazing. Therefore, Matuska and Zmrhal [22] derived a Nusselt correlation from several correlations similar to Hollands et al. [23] to verify the Nusselt number approach for an increased range of collector slopes. Equation (38) describes the Nusselt number for the air gap valid for slopes β from 0° to 90°. The characteristic length L
a
-
g in that system is the distance between absorber and glazing. The value of f
β in (39) constitutes an angle factor for the Nusselt number that depends on the collector slope which was obtained as quadratic function in the derived correlation. The approach is valid for 1 × 104 < Ra < 2 × 106 covering the necessary operational range of solar-thermal collectors. The correlation, too, showed good results independent of the aspect ratio between air gap length and air gap distance. Several approaches, however, are limited with regard to the collector geometry. Consider(38)
N
u
a
-
g
=
f
β
R
a
0.29
(39)
f
β
=
0.1464
-
2.602
e
-
4
β
-
2.046
e
-
6
β
2
.The approach for the Nusselt number with heat transfer downwards is derived from Arnold et al. [24]. Therefore, (40) is used on the basis of the calculated Nusselt number for heat transfer in vertical air layers (β
=
90
°) derived from (38) and (39). The validity of this approach is similar to that of the approach of (38). However, the aspect ratio between air gap length and air gap distance has to be larger than 20. Consider(40)
N
u
a
-
g
=
1
+
N
u
a
-
g
β
=
90
°
-
1
sin
⁡
180
°
-
β
.


### 3.2.5. Convection between Absorber and Casing

The heat transfer in the enclosed air layer between absorber plate and casing is described by the Nusselt correlation for convection between absorber and glazing. The usual approach for the region behind the absorber is downwards free convection expressed by (40).


## 3.2.1. Forced Convection in the Fluid Channels

The forced convection in the absorber pipes is calculated according to the approach recommended by Gnielinski [15]. The characteristic length L is the pipe diameter d in this case. The flow in the pipes was assumed to be constant along the entire length without entrance effects. The temperature dependent properties of the fluid were referenced to the fluid node in the centre of the pipe.The hydrodynamically developed laminar flow for Reynolds numbersR
e below 2,300 and for a constant axial wall temperature is defined with a constant Nusselt number of 3.657 [16].In the transition region between laminar and fully developed turbulent flow with Reynolds numbers in the range of 2,300 and 10,000, the Nusselt number is calculated by the following for a constant axial wall temperature:(27)
N
u
f
=
1
-
γ
N
u
l
a
m
,
2300
+
γ
N
u
t
u
r
b
,
10
4
γ
=
R
e
-
2
,
300
10
4
-
2
,
300
f
o
r
0
≤
γ
≤
1
N
u
l
a
m
,
2300
=
49.371
+
N
u
l
a
m
,
2,2300
-
0.7
3
+
N
u
l
a
m
,
3,2300
3
3
N
u
l
a
m
,
2,2300
=
1.615
2
,
300
P
r
d
l
3
N
u
l
a
m
,
3,2300
=
2
1
+
22
P
r
6
2
,
300
P
r
d
l
N
u
t
u
r
b
,
10
4
=
0.0308
/
8
10
4
P
r
1
+
12.7
0.0308
/
8
P
r
2
/
3
-
1
1
+
d
l
2
/
3
.The Nusselt number for fully developed turbulent flow by Gnielinski [17] for Reynolds numbers larger than 10,000 is described by (28). The friction factor f can be determined by (29). Therefore,(28)
N
u
f
=
f
/
8
R
e
P
r
1
+
12.7
f
/
8
P
r
2
/
3
-
1
1
+
d
l
2
/
3
(29)
f
=
1.8
l
o
g
10
R
e
-
1.5
-
2
.These approaches are valid for a ratio between the pipe diameter d and the pipe length l less than or equal to unity as well as for Prandtl numbers P
r between 0.6 and 1,000.


## 3.2.2. Convection at the Glazing

The convection on the outer surface of the glazing is a combination of free and forced convection. Based on experiments, Churchill [18] developed an approach for an overlay of both convections as shown in the following:(30)
N
u
g
-
a
m
b
=
N
u
f
r
e
e
,
g
-
a
m
b
3
+
N
u
f
o
r
c
e
d
,
g
-
a
m
b
3
3
.The Nusselt number for forced convection with turbulent boundary layer is considered according to (31) by Schlichting [19]:(31)
N
u
f
o
r
c
e
d
,
g
-
a
m
b
=
0.037
R
e
0.8
P
r
1
+
2.443
R
e
-
0.1
P
r
2
/
3
-
1
.Therefore, the air flow caused by the wind is assumed to occur along the collector width of the glazing according to Figure 11.Figure 11
Wind at the inclined collector glazing.The average Reynolds number in (32) is calculated by the wind velocity u
w
i
n
d, the characteristic length of the glazing L
g
-
a
m
b, and the kinematic viscosity ν
a
i
r. The equation is valid for Reynolds numbers between 5 × 105 and 1 × 107 as well as for Prandtl numbers between 0.6 and 2,000. Even during low temperatures and wind velocities, the Reynolds numbers are in the range of the validity. Hence, laminar flow for forced convection is neglected. Consider(32)
R
e
=
u
w
i
n
d
L
g
-
a
m
b
ν
a
i
r
.For free laminar convection at the glazing, two cases are defined by Fujii and Imura [20]. The most important case for thermal collector simulation is the heat transfer upwards at the inclined glazing. During collector operation, the glazing temperature turns out to be larger than the ambient temperature. Hence, the collector’s heat losses are driven by the upward flowing air. The second case with lower glazing temperatures is only relevant for system simulations with fluid temperatures below the ambient air temperatures. Figure 12 shows these two cases for the glazing.Figure 12
Free convection at the glazing with heat transfer upwards (a) and downwards (b).
(a)
(b)The Nusselt number for free laminar convection at the glazing with heat transfer upwards or downwards at vertical planes is described in (33) [21] and is also valid for Rayleigh numbers between 1 × 10−1 and 1 × 1012. Using the correction factor f
P
r in (34) [21] allows Prandtl numbers in the range between 0.001 and nonfinite. This correlation for the Nusselt numbers is also valid for inclined surfaces by substituting the Rayleigh number R
a for the adjusted Rayleigh number R
a
β according to (35) [20]. This adjustment considers the reduced component of gravity at the inclined plane instead of the acceleration due to gravity g. In this case, the characteristic length L
g
-
a
m
b is the length of the glazing. Consider(33)
N
u
f
r
e
e
,
g
-
a
m
b
=
0.825
+
0.387
R
a
f
P
r
6
2
f
o
r
R
a
<
R
a
c
r
i
t
(34)
f
P
r
=
1
+
0.492
P
r
9
/
16
-
16
/
9
f
o
r
0.001
<
P
r
<
∞
(35)
R
a
β
=
R
a
sin
⁡
β
.Boundary layer separation may occur for free convection at the glazing with heat transfer upwards. Hence, there is a change from laminar to turbulent flow. The laminar flow occurs for Rayleigh numbers below the critical Rayleigh number R
a
c
r
i
t depending on the collector slope β from (36). Figure 13 shows the gradient of the upper limit of the critical Rayleigh number for the described laminar convection. Consider(36)
R
a
c
r
i
t
=
10
8.9
-
0.00178
90
°
-
β
1.82
.Above the critical Rayleigh number, free turbulent convection at the glazing occurs. Fujii and Imura [20] defined the following accordingly:(37)
N
u
t
u
r
b
,
g
-
a
m
b
=
0.56
R
a
c
r
i
t
cos
⁡
90
°
-
β
0.25
+
0.13
R
a
3
-
R
a
c
r
i
t
3
.Figure 13
Upper limit of the critical Rayleigh number depending on the collector slope.


## 3.2.3. Convection at the Casing

Duffie and Beckman [12] regarded the heat transfer from the casing to the ambient by radiation and convection as negligible because of the significant impact on the backside insulation and the resulting low temperature. In this model, the heat transfer coefficients from the backside to the ambient are calculated for a hypothetical windless region between collector and roof by natural convection and radiation. The effect of a chimney was also neglected. The heat transfer via convection was assumed as free convection at an inclined plane. Figure 14 describes the two occurring cases for solar collectors.Figure 14
Free convection at the casing with heat transfer downwards (a) and upwards (b).
(a)
(b)The most common case for the backside sheet is a heated plate transmitting heat downwards without any separation of the boundary layer. Thus, the Nusselt number for laminar flow from (33) to (35) serves as an exploration for the free convection behind the collector. In this case, L
c
-
a
m
b is the characteristic length of the casing.The special case of a cooled surface absorbing heat from the ambient air from below is described by (33) to (35) for laminar flow. Above the critical Rayleigh number, the separation of the boundary layer is expressed by means of (36). Furthermore, the approach for turbulent flow from (37) has to be used.


## 3.2.4. Convection between Absorber and Glazing

Heat transfer in an enclosed air gap is driven by heat conduction as well as by free convection of the circulating gas layer. This free convection between glazing and absorber with the two possible directions of heat transfer is shown in Figure15.Figure 15
Free convection in the enclosed gas layer between absorber and glazing with heat transfer upwards (a) and downwards (b).
(a)
(b)Upwards heat transfer usually appears in collectors between absorber and glazing. Therefore, Matuska and Zmrhal [22] derived a Nusselt correlation from several correlations similar to Hollands et al. [23] to verify the Nusselt number approach for an increased range of collector slopes. Equation (38) describes the Nusselt number for the air gap valid for slopes β from 0° to 90°. The characteristic length L
a
-
g in that system is the distance between absorber and glazing. The value of f
β in (39) constitutes an angle factor for the Nusselt number that depends on the collector slope which was obtained as quadratic function in the derived correlation. The approach is valid for 1 × 104 < Ra < 2 × 106 covering the necessary operational range of solar-thermal collectors. The correlation, too, showed good results independent of the aspect ratio between air gap length and air gap distance. Several approaches, however, are limited with regard to the collector geometry. Consider(38)
N
u
a
-
g
=
f
β
R
a
0.29
(39)
f
β
=
0.1464
-
2.602
e
-
4
β
-
2.046
e
-
6
β
2
.The approach for the Nusselt number with heat transfer downwards is derived from Arnold et al. [24]. Therefore, (40) is used on the basis of the calculated Nusselt number for heat transfer in vertical air layers (β
=
90
°) derived from (38) and (39). The validity of this approach is similar to that of the approach of (38). However, the aspect ratio between air gap length and air gap distance has to be larger than 20. Consider(40)
N
u
a
-
g
=
1
+
N
u
a
-
g
β
=
90
°
-
1
sin
⁡
180
°
-
β
.


## 3.2.5. Convection between Absorber and Casing

The heat transfer in the enclosed air layer between absorber plate and casing is described by the Nusselt correlation for convection between absorber and glazing. The usual approach for the region behind the absorber is downwards free convection expressed by (40).


## 3.3. Radiation

Heat transfer caused by radiation is based on the Stefan-Boltzmann law of a black body radiator. The energy transfer is described with an emission coefficientε
1
=
1 (after Eicker [11]) by (41). The Stefan-Boltzmann constant σ is shown in (42). Therefore,(41)
Q
˙
1
=
σ
A
1
T
1
4
(42)
σ
=
5.67
e
-
8
W
m
-
2
K
-
4
.Bodies above a temperature of 0 K emit radiating power. Radiating power from other bodies placed in a room is also absorbed. The position of the bodies also has great influence on the behaviour as described by their form factors Φ. The form factor is the ratio between the radiation receiving surface and the entire hemisphere. Figure 16 shows the heat exchange of the two surfaces.Figure 16
Heat fluxes via radiation between two elements (according to [11]).The net radiation exchangeQ
˙
r
a
d between two bodies takes the properties of the bodies into account and defines the heat transfer. Therefore, the difference between the radiating powers reaching the other body is taken from (43)
Q
˙
r
a
d
=
Q
˙
r
a
d
,
1
→
2
-
Q
˙
r
a
d
,
2
→
1
=
Φ
12
σ
A
1
T
1
4
-
Φ
21
σ
A
2
T
2
4
.The definition of the net radiation exchange between two bodies with an emission coefficient ε smaller than unity is a more complex approach. Therefore, a multiple of reflections has to be considered.The radiation exchange in solar-thermal flat-plate collectors embraces two special cases. Firstly, there is radiation exchange between two parallel plates, for example, between absorber and glazing; secondly, radiation occurs from a limited surface to the sky hemisphere.The heat transfer by radiation is determined according to (44). Therefore, the heat transfer coefficients for radiation of the collector have to be defined. Consider(44)
Q
˙
r
a
d
=
h
r
a
d
A
T
1
-
T
2
.


### 3.3.1. Radiation between Two Parallel Plates

Due to the fact that the heat emitted is completely absorbed by the opposed area, the form factors of the two surfaces are equal (Φ
12
=
Φ
21
=
1). It is assumed that the effects of infinitely expanded parallel surfaces apply for heat transfer in solar collectors. Multiple reflections have to be considered for grey emitters as well by introducing their emission coefficients ε into the approach. The heat transfer coefficient for radiation between two parallel plates is given in the following:(45)
h
r
a
d
=
σ
1
/
ε
1
+
1
/
ε
2
-
1
T
1
2
+
T
2
2
T
1
+
T
2
.This radiation exchange is necessary for the setup of glazing and absorber, absorber and casing, and casing and roof.


### 3.3.2. Radiation between Glazing and Sky

The radiation between glazing and cold sky is a special phenomenon of heat transfer [11]. The surface of glazing A
g is much smaller than the total sky hemisphere A
s
k
y. Hence, radiation from the glazing is completely emitted to the sky. Thus, the form factor of the glazing Φ
g
-
s
k
y is 1. In contrast to that, the form factor from the sky Φ
s
k
y
-
g is almost zero since the proportion of the glazing surface being seen from the sky hemisphere is very small. Therefore, the emitted heat from the sky is neglected because of the infinite area of the sky’s half-space. The heat transfer coefficient is shown in (46). The temperatures have to be inserted in Kelvin because radiation refers to the absolute temperature of the body or rather surface. Consider(46)
h
r
a
d
,
g
-
s
k
y
=
σ
ε
g
T
g
2
+
T
s
k
y
2
T
g
+
T
s
k
y
.The temperature of the sky differs from the ambient air temperature [12]. Approaches for the definition of the sky temperature are related to the air temperature [25], to the water vapour pressure [26], or to the dew point temperature [27], for example.


## 3.3.1. Radiation between Two Parallel Plates

Due to the fact that the heat emitted is completely absorbed by the opposed area, the form factors of the two surfaces are equal (Φ
12
=
Φ
21
=
1). It is assumed that the effects of infinitely expanded parallel surfaces apply for heat transfer in solar collectors. Multiple reflections have to be considered for grey emitters as well by introducing their emission coefficients ε into the approach. The heat transfer coefficient for radiation between two parallel plates is given in the following:(45)
h
r
a
d
=
σ
1
/
ε
1
+
1
/
ε
2
-
1
T
1
2
+
T
2
2
T
1
+
T
2
.This radiation exchange is necessary for the setup of glazing and absorber, absorber and casing, and casing and roof.


## 3.3.2. Radiation between Glazing and Sky

The radiation between glazing and cold sky is a special phenomenon of heat transfer [11]. The surface of glazing A
g is much smaller than the total sky hemisphere A
s
k
y. Hence, radiation from the glazing is completely emitted to the sky. Thus, the form factor of the glazing Φ
g
-
s
k
y is 1. In contrast to that, the form factor from the sky Φ
s
k
y
-
g is almost zero since the proportion of the glazing surface being seen from the sky hemisphere is very small. Therefore, the emitted heat from the sky is neglected because of the infinite area of the sky’s half-space. The heat transfer coefficient is shown in (46). The temperatures have to be inserted in Kelvin because radiation refers to the absolute temperature of the body or rather surface. Consider(46)
h
r
a
d
,
g
-
s
k
y
=
σ
ε
g
T
g
2
+
T
s
k
y
2
T
g
+
T
s
k
y
.The temperature of the sky differs from the ambient air temperature [12]. Approaches for the definition of the sky temperature are related to the air temperature [25], to the water vapour pressure [26], or to the dew point temperature [27], for example.


## 3.4. Heat Conduction

Heat conduction is prevalent at glazing, absorber, and casing. The following describes the heat transfer that depends on the heat conductivityk and the thickness or rather distance s of the component [14]:(47)
Q
˙
c
o
n
d
=
h
c
o
n
d
A
T
1
-
T
2
h
c
o
n
d
=
k
s
.Vertical heat transfer through thin plates in the thermal network like glazing or absorber can be neglected. Relevant heat conduction occurs between the nodes of the absorber perpendicular to the cross-sectional area of the plate and towards the riser pipe (Figure 17).Figure 17
Heat conduction between the absorber nodes.The second heat conduction is seen in the collector insulation. These heat losses occur at the side and backside insulation. In contrast to the massive solid materials such as metals, the heat conductivity cannot be regarded as constant. The cavernous structure of the material also enables heat transfer effects of the enclosed air and radiation. Thus, the heat conductivity increases according to the material temperature. For dry mineral wool, Ochs et al. [28] derived a linear approach for the temperature depending calculation of the heat conductivity k
m
w according to (48). The influence on the heat conductivity of the unknown humidity in the mineral wool is neglected. Moist mineral wool usually has an increased heat conductivity. Consider(48)
k
m
w
=
0.035
W
m
-
1
K
-
1
+
2.0
e
-
4
T
m
w
.


## 3.5. Edge Losses

The edge losses, which comprise the heat losses through the collector frame or rather side, cannot be mathematically calculated in detail. Tabor [29] describes an approach considering the losses only connected to the absorber node. The other nodes have an adiabatic behaviour against the edge. Equation (49) shows the heat transfer between absorber and ambient air as well as the heat transfer coefficient considering only the heat conduction through the side insulation. Consider(49)
Q
˙
a
-
e
d
g
e
=
h
e
d
g
e
A
e
d
g
e
T
a
-
T
a
m
b
h
e
d
g
e
=
k
c
s
e
d
g
e
.


## 3.6. Overall Heat Transfer Coefficients between the Nodes

The heat transfer coefficientsh of the collector connect the nodes and the ambient. Most of them are embedded in a network of the thermal resistances to calculate a heat flow rate Q
˙. The electrical analogy is used to generate the overall heat transfer coefficients between the nodes based on the single thermal effects [14, 16]. This analogy is used to generate overall resistances by combining single resistances of serial and parallel connections according to Figure 5. The required overall heat transfer coefficients (related to the area) of the energy fluxes from (3) to (9) are shown in the following:(50)
1
R
t
h
,
g
-
a
m
b
=
h
c
o
n
d
,
g
A
+
1
h
r
a
d
,
g
-
s
k
y
A
+
1
h
c
o
n
v
,
g
-
a
m
b
A
-
1
1
R
t
h
,
a
-
g
=
1
h
r
a
d
,
a
-
g
A
a
+
1
h
c
o
n
v
,
a
-
g
A
a
-
1
+
h
c
o
n
d
,
g
A
a
1
R
t
h
,
a
1
-
f
=
h
c
o
n
d
,
a
1
2
A
w
a
l
l
+
h
c
o
n
v
,
a
1
-
f
A
p
i
p
e
1
R
t
h
,
a
-
c
=
1
h
r
a
d
,
a
-
c
A
a
+
1
h
c
o
n
v
,
a
-
c
A
a
-
1
+
h
c
o
n
d
,
c
A
a
1
R
t
h
,
c
-
a
m
b
=
h
c
o
n
d
,
c
A
+
1
h
r
a
d
,
c
-
a
m
b
A
+
1
h
c
o
n
v
,
c
-
a
m
b
A
-
1
.


## 3.7. Validation

The model was validated using the parameter set of a tested conventional flat-plate collector. The validation of this model and physical effects is also adequate for models with implemented overheating protection measures or adjusted build-up since the defined physical effects are varied and not exchanged [1, 30].The model was tested in various settings to investigate the accuracy of the results. According to the above-mentioned simplifying assumptions in the model, there are, of course, imperfections in the build-up of the tested collectors. Hence, the differences in the results have to be identified, justified, and evaluated. Collector performance data is generally obtained under steady-state conditions; hence, the validation data presented focus on measurements and the performance of other models under equilibrium conditions.


### 3.7.1. Comparison with an Existing Validated Model

In a first validation step, the developed dynamic multinode model was compared to the results of a one-dimensional multinode model approach by Hafner et al. [5] based on the approach of (1). This validation test is only able to compare the results of the collector efficiency of an already validated parameter model approach. The tested parameters of the collector are shown in Table 1.Table 1
Efficiency curve of the tested collector for the parameter model [31].
Efficiency curve parameters
Zero loss efficiency
η
0
0.798
[—]
Linear heat loss coefficient
a
3.34
[W m−2 K−1]
Quadratic heat loss coefficient
b
0.0075
[W m−2 K−2]
Heat capacity
C
9.5
[kJ K−1]The collector surface of both models was reduced to 1.83 m². Thus, both models consider the absorber surface without header regions.Two test simulations were set up to compare the results of the models. The first test comprised the simulation of the collector efficiency curve like the practical test modus in a solar simulator according to DIN EN 12975-2 [32]. Under constant ambient conditions (G
=
885
W
m
-
2, T
a
m
b
=
27
°
C, and u
w
i
n
d
=
2
m
s
-
1), a single collector was calculated. The heat carrier, a water-glycol mixture (33.3%), had a flow rate of m
˙
=
0.08
k
g
s
-
1 and variable temperatures through the collector.The percentage difference of the collector efficiency describes the accuracy of the model approach according to the following:(51)
Percentage difference
=
η
dynamic model
η
paramter model
-
1
100
%
. The values range between −1.0% and +3.3% for flow temperatures from 25°C to 90°C (Figure 18).Figure 18
Comparison of the simulated efficiency curves of the parameter model and the developed dynamic collector model.The differences are relatively low according to the assumptions in the model. Also, the simulated temperature range in the test covers a large operational range of the collector to prove the use of the model.The deviation is assumed to originate from the discretisation in the model. The major influence in the test is the transition from laminar to turbulent flow at higher values of the reduced temperature. The change in the model curve is remarkable because of an averaged calculation of fully developed flow conditions for the entire riser pipe. In contrast to that, the real collector shows a continuous development of flow conditions along the pipe length. There are also neglected entrance effects increased by the vertical connection between header and riser.Next to the analysis of efficiency, the models were compared to each other in an annual simulation of a solar-thermal heating system for domestic hot water and space heating. The deviation of the solar yield was found to be +2.9%, similar to the results found in the comparison of efficiency. The deviation originates from the fact that the major operational range of the collectors is at high reduced temperatures. This is caused by both high flow and ambient air temperatures for hot water preparation in summer as well as moderate flow and low ambient air temperatures in the heating period for space heating.


### 3.7.2. Comparison with Measurements

In a second step, the simulation results were compared to measurement data from analysing an already exposed field-testing collector; relevant collector characteristics are given in Table2 (cf. [1]).Table 2
Characteristics of the field-testing collector.
Collector
Gross length
1.987
[m]
Gross width
1.046
[m]
Gross area
2.078
[m²]
Aperture area
1.903
[m²]
Weight (empty)
39
[kg]
Liquid content
1.12
[l]
Glazing
Thickness
0.0032
[m]
Transmission coefficientτ
0.918
[—]
Absorber
Absorber area
1.903
[m²]
Absorption coefficient
0.95
[—]
Emission coefficient
0.04
[—]
Number of risers
10
[—]The measurements were performed using Institute’s solar simulator. The procedures defined in the European testing standard DIN EN 12975-2 [32] were applied using the equipment postulated so that the results presented are in line with the standard. However, degradation of the absorber coating and humidity in the mineral wool can be named as some of the reasons, resulting in the effect that the efficiency of the collector is reduced in comparison to results from previous efficiency curves of a new collector. Therefore, the emission coefficient of the absorber coating and the heat conductivity of the mineral wool in the model were adjusted to represent the typical degradation of the components.The first test was the comparison of the efficiency curves. Therefore, the measured conditions of the solar simulator were reproduced (G
=
890
W
m
-
2, T
a
m
b
=
31
°
C, and u
w
i
n
d
=
1
m
s
-
1); the heat carrier was water (m
˙
=
0.038
k
g
s
-
1). Figure 19 shows the good coverage between measurement and simulation. According to the following, the percentage difference is only in the range between +0.4% and +3.1%:(52)
Percentage difference
=
η
dynamic model
η
measurement
-
1
100
%
.The second test also showed satisfying conformity in particular with regard to individual component temperatures. During stagnation at insolation of 895 W m−2 and an ambient air temperature of 33°C, the absorber temperature of model and test (average of 9 measurement points) differs only by +1 K. Also the average temperatures (average of 3 measurement points) at glazing and backside vary below +1 K (Table 3).Table 3
Collector component temperatures during stagnation of a conventional flat-plate collector in testing and simulation.
Testing
Model
Difference
Glazing
71°C
71°C
±0 K
Absorber
174°C
175°C
+1 K
Backside (outside)
54°C
55°C
+1 KFigure 19
Comparison of the measured and simulated efficiency of the collector.


### 3.7.3. Evaluation of the Test Results

All in all, the model proves to be valid according to the test collector. Slight differences are identified as being within the theoretical limits of the model. Thanks to the assumptions that have been made, the model generates reliable results for flat-plate collectors in due consideration of the limits of the model. Comparable model approaches of state-of-the-art collectors show similar low uncertainties in simulation of efficiency or rather solar yield. The model of Koo [33] showed percentage differences between −2% and +3% in comparison to experimental data of the collector efficiency curve. Matuska and Zmrhal [22] described a percentage difference of about 3-4% in the model validation. Furthermore, Cadafalch [34] identified differences below 2%. Thus, there is confidence that the model can be used for the investigation of entire systems, for identifying the solar yield, and for the analysis of the component temperatures in the collector.


## 3.7.1. Comparison with an Existing Validated Model

In a first validation step, the developed dynamic multinode model was compared to the results of a one-dimensional multinode model approach by Hafner et al. [5] based on the approach of (1). This validation test is only able to compare the results of the collector efficiency of an already validated parameter model approach. The tested parameters of the collector are shown in Table 1.Table 1
Efficiency curve of the tested collector for the parameter model [31].
Efficiency curve parameters
Zero loss efficiency
η
0
0.798
[—]
Linear heat loss coefficient
a
3.34
[W m−2 K−1]
Quadratic heat loss coefficient
b
0.0075
[W m−2 K−2]
Heat capacity
C
9.5
[kJ K−1]The collector surface of both models was reduced to 1.83 m². Thus, both models consider the absorber surface without header regions.Two test simulations were set up to compare the results of the models. The first test comprised the simulation of the collector efficiency curve like the practical test modus in a solar simulator according to DIN EN 12975-2 [32]. Under constant ambient conditions (G
=
885
W
m
-
2, T
a
m
b
=
27
°
C, and u
w
i
n
d
=
2
m
s
-
1), a single collector was calculated. The heat carrier, a water-glycol mixture (33.3%), had a flow rate of m
˙
=
0.08
k
g
s
-
1 and variable temperatures through the collector.The percentage difference of the collector efficiency describes the accuracy of the model approach according to the following:(51)
Percentage difference
=
η
dynamic model
η
paramter model
-
1
100
%
. The values range between −1.0% and +3.3% for flow temperatures from 25°C to 90°C (Figure 18).Figure 18
Comparison of the simulated efficiency curves of the parameter model and the developed dynamic collector model.The differences are relatively low according to the assumptions in the model. Also, the simulated temperature range in the test covers a large operational range of the collector to prove the use of the model.The deviation is assumed to originate from the discretisation in the model. The major influence in the test is the transition from laminar to turbulent flow at higher values of the reduced temperature. The change in the model curve is remarkable because of an averaged calculation of fully developed flow conditions for the entire riser pipe. In contrast to that, the real collector shows a continuous development of flow conditions along the pipe length. There are also neglected entrance effects increased by the vertical connection between header and riser.Next to the analysis of efficiency, the models were compared to each other in an annual simulation of a solar-thermal heating system for domestic hot water and space heating. The deviation of the solar yield was found to be +2.9%, similar to the results found in the comparison of efficiency. The deviation originates from the fact that the major operational range of the collectors is at high reduced temperatures. This is caused by both high flow and ambient air temperatures for hot water preparation in summer as well as moderate flow and low ambient air temperatures in the heating period for space heating.


## 3.7.2. Comparison with Measurements

In a second step, the simulation results were compared to measurement data from analysing an already exposed field-testing collector; relevant collector characteristics are given in Table2 (cf. [1]).Table 2
Characteristics of the field-testing collector.
Collector
Gross length
1.987
[m]
Gross width
1.046
[m]
Gross area
2.078
[m²]
Aperture area
1.903
[m²]
Weight (empty)
39
[kg]
Liquid content
1.12
[l]
Glazing
Thickness
0.0032
[m]
Transmission coefficientτ
0.918
[—]
Absorber
Absorber area
1.903
[m²]
Absorption coefficient
0.95
[—]
Emission coefficient
0.04
[—]
Number of risers
10
[—]The measurements were performed using Institute’s solar simulator. The procedures defined in the European testing standard DIN EN 12975-2 [32] were applied using the equipment postulated so that the results presented are in line with the standard. However, degradation of the absorber coating and humidity in the mineral wool can be named as some of the reasons, resulting in the effect that the efficiency of the collector is reduced in comparison to results from previous efficiency curves of a new collector. Therefore, the emission coefficient of the absorber coating and the heat conductivity of the mineral wool in the model were adjusted to represent the typical degradation of the components.The first test was the comparison of the efficiency curves. Therefore, the measured conditions of the solar simulator were reproduced (G
=
890
W
m
-
2, T
a
m
b
=
31
°
C, and u
w
i
n
d
=
1
m
s
-
1); the heat carrier was water (m
˙
=
0.038
k
g
s
-
1). Figure 19 shows the good coverage between measurement and simulation. According to the following, the percentage difference is only in the range between +0.4% and +3.1%:(52)
Percentage difference
=
η
dynamic model
η
measurement
-
1
100
%
.The second test also showed satisfying conformity in particular with regard to individual component temperatures. During stagnation at insolation of 895 W m−2 and an ambient air temperature of 33°C, the absorber temperature of model and test (average of 9 measurement points) differs only by +1 K. Also the average temperatures (average of 3 measurement points) at glazing and backside vary below +1 K (Table 3).Table 3
Collector component temperatures during stagnation of a conventional flat-plate collector in testing and simulation.
Testing
Model
Difference
Glazing
71°C
71°C
±0 K
Absorber
174°C
175°C
+1 K
Backside (outside)
54°C
55°C
+1 KFigure 19
Comparison of the measured and simulated efficiency of the collector.


## 3.7.3. Evaluation of the Test Results

All in all, the model proves to be valid according to the test collector. Slight differences are identified as being within the theoretical limits of the model. Thanks to the assumptions that have been made, the model generates reliable results for flat-plate collectors in due consideration of the limits of the model. Comparable model approaches of state-of-the-art collectors show similar low uncertainties in simulation of efficiency or rather solar yield. The model of Koo [33] showed percentage differences between −2% and +3% in comparison to experimental data of the collector efficiency curve. Matuska and Zmrhal [22] described a percentage difference of about 3-4% in the model validation. Furthermore, Cadafalch [34] identified differences below 2%. Thus, there is confidence that the model can be used for the investigation of entire systems, for identifying the solar yield, and for the analysis of the component temperatures in the collector.


## 4. Conclusion

The dynamic collector model described enables the theoretical investigation of flat-plate collector designs. The nodes of the components with the implemented physics being connected to each other represent the basis for a simulation model for flexible use regarding geometry, materials, and additional components. The thermal loads on the components and the influence on the collector efficiency of collector designs can be simulated in detail. In contrast to that, commonly used model approaches are not able to consider these requirements for output data. Thus, this model approach can be a necessary tool for product development of state-of-the-art collectors or polymeric components as well as for the investigation of overheating protection measures within collectors.The various validation procedures showed good results in comparison to existing mathematical models of flat-plate collectors. Furthermore, the component temperatures obtained from testing and simulation are very close. Hence, the collector model can be regarded as a reliable tool for wide application spectrum in simulation.



---

*Source: 280694-2015-11-05.xml*