Unsteady approach Heat transfer in solid layersHeat transfer in solid layers

The process of heat transfer in a solid layers sub-domains is dominated by the con-duction. The general form of governing differential equation describing the unsteady heat conduction in solid bodies (in vector notation) is [159]:

∂(CpρT)

∂t = ∇ · (λ∇T ) + qV, (5.3.1)

where gV is the internal heat source.

In the FE model, it was assumed that the thermo-physical properties of solid layers (listed in Table 5.1), including roofing, battens, polypropylene pipes and insulation layer, are isotropic and temperature-independent, and there are no internal heat

5.3. FE MODEL 163

sources. Hence, the equation describing the unsteady heat conduction in the solid layers of HSC can be given as:

C_pρ∂T

∂t = λ∇²T. (5.3.2)

Heat transfer in water flowing through pipes

In the FE model of water flowing through the pipes, the heat exchange process is constrained to the conduction and forced convection. According to the control strategy assumed in order to investigate the capability of HSC to supply space-heating systems (Section 5.5), the average temperature of water flowing through the collector’s pipes is kept between 20^◦C and 25^◦C. With regard to assumed range of water flow velocities, not exceeding 0.02 m/s, the laminar flow regime is to be considered throughout the time of collector operation (Section 5.2.2.1). In many practical problems, the flow of a conventional liquid, such as water, is incompressible to a high degree of accuracy [159]. Taking into account low magnitudes of flow velocities and low temperature variations in an operating fluid, it is reasonable to assume the water in HSC as an incompressible fluid. Assuming that the water density distribution remains uniform in space and constant in time, the continuity equation (Eq. 5.2.1) can be rewritten as:

∇ · U = 0. (5.3.3)

In the present model, it is also assumed that a viscosity of the water equals to zero.

In general, the assumption of a non-viscous water is reasonable for low velocity fluid flows in a distance from surfaces of boundary walls [159] where the energy dissipation is negligible. With respect to the fluid flow problem under consideration, such an assumption results in a constant fluid flow velocity profile (Fig. 5.34a) in the entire computational domain. The velocity vector U is defined as:

U=

where uz is the vector component along the main flow direction. This has, obviously, an impact on heat transfer process in the near-wall region. In case of the laminar flow of a viscous incompressible fluid, there is a thin layer near the pipe surface where the velocity decreases to zero (Fig. 5.34b). In this thin layer, the heat is transferred mostly by diffusion. Omitting of this phenomenon in the FE model might lead to

a significant errors in estimation of the heat exchange rate between the pipes and flowing water. In order to reach realistic results, it is reasonable to modify heat transfer conditions in the near-wall region of the flowing water. This was achieved through the usage of interface-type boundary condition for heat transfer between surfaces of polypropylene pipes and domains of water (Fig. 5.34c). The mathematical interpretation of applied approach is described later on.

Figure 5.34: Velocity profiles for laminar flow of incompressible fluid according to a) non-viscous fluid model b) viscous fluid model c) applied approach

In the present model, the heat transfer is described by an extended form of the heat conduction equation (Eq. 5.3.2). Since the kinetic energy is negligible due to constant pressure flow field, the governing equation is derived from the internal energy balance equation [159]:

ρDe

Dt = λ∇²T − p∇ · U+ φ + qV, (5.3.5) where φ is the energy dissipation.

Substituting Eq. 5.3.3, neglecting the energy dissipation, φ, due to the assumption of a non-viscous fluid, and assuming no internal heat sources (qV = 0), the internal energy balance equation (Eq. 5.3.5) becomes:

ρDe

Dt = λ∇²T. (5.3.6)

Since the pressure effects are neglected, the following approximation may be intro-duced:

de = CpdT. (5.3.7)

5.3. FE MODEL 165

Substituting Eq. 5.3.7 into Eq. 5.3.6 results in a well-known Kirchhoff-Fourier energy equation:

C_pρDT

Dt = λ∇²T. (5.3.8)

Assuming that ux = uy = 0, the Kirchhoff-Fourier energy equation for the unsteady heat transfer problem considered in the FE model of water flowing through the pipes takes the form:

C_pρ ∂T

∂t + uz

∂T

∂z

!

= λ∇²T. (5.3.9)

The above equation is the same as the corresponding equation for heat transfer in solid layers (Eq. 5.3.2), except for the convection term. The assumed magnitudes of density, specific heat and thermal conductivity for the water sub-domain are 998.2 kg/m³, 4183 J/ (kg · K) and 0.599 W/ (m · K), respectively. The magnitudes are appropriate for water temperature of 20^◦C [181]. In reality, variations of these parameters for the considered range of operating temperatures (between 20^◦C and 35^◦C) are negligible.

The above formulation of heat transfer process refers to the case when the max-imum water temperature equal to 35^◦C. Due to the control strategy (Section 5.5), such a condition occurs only if the target of the collector application is to supply a space-heating system based on seasonal heat storage systems. An additional aim of this study is to evaluate the ability of HSC to supply domestic hot water systems.

According to the control strategy considered for this application, the water tem-perature can vary between 20^◦C and 80^◦C. For such conditions, the usage of the heat transfer model formulated for water flowing through the collector’s pipes may result in the small errors, reduced by controlled convective heat transfer coefficient in the interface between water and pipe, determined on basis of CFD simulations (Section 5.2.5). However, the comparison of the FE model with CFD one should be conducted with regard to all the expected magnitudes of water temperature at outlets from the pipes.

Heat transfer in air-cavity

In opposition to the CFD model (Section 5.2), the FE model of HSC includes ra-diative heat exchange between surfaces composing the air-cavity and surfaces of collector’s pipes. Since the air layer is considered as an orthotropic solid body

in terms of convective-equivalent thermal conductivity coefficients (Section 5.3), the analysis of an unsteady heat transfer process in the air-cavity required the governing equations of conduction and radiation to be included in the mathematical model.

The equation for the unsteady conductive heat transfer within the air layer is analogous to Eq. 5.3.2. The convective-equivalent horizontal, λHeqv, and vertical, λ_{V eqv}, thermal conductivity coefficients are dependent on current magnitudes of sol-air temperature, Tsol, and water mass flux, m^·, whereas the other thermo-physical parameters are constant and temperature-independent. The adopted functions for convective-equivalent thermal conductivity coefficients are described in Section 5.2.5.

The assumed magnitudes of density and specific heat for the air layer are 1.247 kg/m³ and 1006 J/ (kg · K), respectively. One should notice that in reality both the param-eters vary in dependence on temperature change. It has, however, a negligible effect on the overall HSC performance. The assumed magnitudes of the parameters are in use for numerical simulation of solar collectors [75]. Up to this stage, the mathemat-ical model includes indirectly an influence of the buoyancy-driven airflow on heat transfer process in the air-cavity. In general case, a reliable analysis of the heat ex-change rate between surfaces forming the air-cavity and collector’s pipes, due to the convective heat transfer, requires taking into account the thermal resistance of the near-wall region where the heat is transferred mostly by the diffusion. Similarly to the model of heat transfer in water flowing through the pipes, this was achieved by controlling convective heat transfer coefficient through the interface between surfaces of polypropylene pipes and the air-cavity sub-domain. An appropriate mathematical interpretation of the applied approach is described later on.

The radiative heat exchange between surfaces composing the air-cavity and sur-faces of collector pipes is governed by the grey body radiation theory that means that the monochromatic emissivity of the body is independent of the wavelength of the radiation propagation. Accordingly, the radiative heat exchange between two differential surfaces in the air-cavity of HSC is defined by the following equation [139]:

Q₁₋₂= 1−2σF₁φ₁₋₂^T₁⁴− T₂^4, (5.3.10)

where Q1−2is the radiation heat flux between the bodies ’1’ and ’2’, 1−2is the emis-sivity between the bodies ’1’ and ’2’, σ is the Stefan-Boltzmann constant, F1 is the surface area of body with temperature T1, φ1−2 is the view factor (sometimes called a configuration factor or a shape factor), T1 and T2 are the absolute temperatures of bodies between which the heat exchange occurs.

5.3. FE MODEL 167

For two finite areas F1 and F2, the view factor, φ1−2, is given by [139]:

φ₁₋₂ = 1 F₁

ˆ

F1

ˆ

F2

cos β1cos β2

πR²₁₋₂ dF1dF2, (5.3.11) where β1 and β2 are the angles between the areas’ normal, R1−2 is the distance between two differential areas. The view factor also satisfies the reciprocity relation:

F₁φ₁₋₂= F2φ₂₋₁. (5.3.12) Boundary conditions

The boundary surfaces of solid layers, in which the heat is transferred by conduction, are shown in Fig. 5.35. Assuming that the collector separates indoor zones at fixed temperature from outdoor unsteady climate conditions, boundary conditions at the internal, Si, and external, Se, surfaces of HSC are defined by the Newton’s law. A

Figure 5.35: Boundary surfaces of solid bodies sub-domain in FE model heat exchange rate by convection and radiation on the internal solar collector surface, S_i, is defined by the convective/radiative heat transfer coefficient, hi. Consequently, the boundary condition at the internal collector surface, Si, used to complete the energy equation is:

λ∂T(x, t)

∂n

    _S

i

= hi[TSi(t) − Ta,in] , (5.3.13) where TSi is the average temperature of the internal collector surface and Ta,inis the indoor air temperature. The convective/radiative heat transfer coefficient, hi, is

con-stant in time and defined according to the ISO standard [136] as hi = 8.1 W/(m²·K).

The indoor air temperature, Ta,in, is assumed to be constant in time throughout the unsteady simulation (Ta,in = 20^◦C).

A heat exchange between the external surface of the collector, Se, and outdoor environment is considered by a convective and radiative heat exchange, separately.

The convection is defined by the convective heat transfer coefficient for external col-lector surface, Se, whereas the radiation is defined by the sol-air temperature, Tsol. Consequently, the boundary condition on external collector surface, Se, is given as:

λ∂T(x, t) where Tsol is the sol-air temperature and TSe is the average temperature of the external collector surface. According to the assumption of time-varying climate conditions, the so far presented formulations for the sol-air temperature (Eq. 5.2.72) and the convective heat transfer coefficient at external collector surface (Eq. 5.2.73) take the following form:

Several FORTRAN subroutines were developed to read the data from the climate database [164] and modify current magnitudes of sol-air temperature and convective heat transfer coefficient at external collector surface for each time increment of the transient numerical simulation. To smooth out variations of the magnitudes, the time-dependent climate variables are approximated locally by a linear function.

A heat exchange between surfaces of two adjacent solid bodies, including roofing material, tilling battens, insulation is based on continuity assumptions. The equa-tions describing the boundary condiequa-tions for the unsteady heat exchange between surfaces of two adjacent, among considered, solid bodies, Ss,1 and Ss,2, are repre-sented as:

T(x, t)|_Ss,1 = T (x, t)|_Ss,2, (5.3.17)

5.3. FE MODEL 169

q(x, t)|_Ss,1 = q (x, t)|_Ss,2, (5.3.18) where q is the heat flux in a direction normal to the surface.

Through the analogy are formulated the equations describing the boundary conid-tions for unsteady heat exchange between surfaces of air-cavity layer, Sa,e(Fig. 5.36), being adjacent to the surfaces of solid bodies, Sc,i.

For the heat exchange between surfaces of the collector pipe and the air layer, being in proximity to one another, an interface is used. The mathematical form of the interface is that of a Cauchy type given below as:

−λ∂T(x, t)

where Sa,iis the air layer surface (Fig. 5.36) being adjacent to the external pipe sur-face Sp,e, ka−p is the heat transfer coefficient describing the resistance to heat flow across the air-pipe interface, TSa,i and TSp,e are the local surface temperatures. The magnitude of ka−pvary in time in dependence on sol-air temperature, Tsol, and water mass fllux, m^·, hence, providing modified heat transfer conditions between the air layer and external pipe surfaces. In reality, the thermal resistance for down-pipe and up-pipe differs due to characteristics of airflow in the cavity. Thus, the two separate functions of ka−p

T_sol,m^{· }for the lower and top pipe were estimated on basis of the CFD analyses (Section 5.2.5). The functions were implemented in the FORTRAN subroutine in order to modify the heat exchange rate at each time-increment of the FE transient simulation.

On other surfaces of solid bodies, Ss,a, and air-cavity, Sa,a, the adiabatic boundary conditions are prescribed. In other words, no heat transfer across these surfaces is allowed. The adiabatic boundary condition used to complete the equation governing the unsteady heat transfer for Ss,a is represented as:

−λ∇T (x, t)|_Ss,a = q|_Ss,a = 0. (5.3.20)

Figure 5.36: Boundary surfaces of air-cavity sub-domain in FE model

The sub-domain of water refers to water-down and water-up models (Fig. 5.37).

Both the models are the same in terms of prescribed boundary condition, hence the following equations are referred to a single model of water. The water flows into the domain across the inlet water surface, Sw,in, and flows out through the outlet water surface, Sw,out. In the FE model, it is assumed that fluid flow velocity is uniform for

Figure 5.37: Boundary surfaces of water sub-domain in FE model

entire domain of water. Its magnitude varies throughout the transient simulation in dependence on outlet water temperature (Section 5.5). The Kirchhoff-Fourier energy equation (Eq. 5.3.9) requires to define explicit water velocities vector in each node of the FE mesh. The water flow velocity, uw(t), is computed from the water mass flux, m· , and the water density, ρ, according to Eq. 5.2.18. A FORTRAN subroutine was developed to control the water velocity at each time increment and in each mesh node during the unsteady heat transfer simulation in HSC.

5.3. FE MODEL 171

According to made assumptions (Section 5.5), the water temperature at the inlet water surface, is uniform (equal to 20^◦C) and constant in time. The Dirichlet boundary condition used to describe the temperature distribution at the inlet water surface, Sw,in, is given as:

T (x, t)|_Sw,in = Tw,in = 20^◦C, (5.3.21) where Tw,in is the inlet water temperature.

At the outlet water surface, Sw,out, the temperature gradient along the vector nor-mal to the surface is constant in time and equal to zero. Therefore, the boundary condition is given as:

For the heat exchange between surfaces of a collector pipe, Sp,i, and a domain of water, Sw,s, being in proximity to one another, an interface is prescribed through the analogy to Eq. 5.3.19 as follows:

−λ∂T(x, t)

where Sw,s is the water layer surface being adjacent to the internal pipe surface S_p,i, kw−pis the heat transfer coefficient describing the resistance to heat flow across the water-pipe interface, TSp,i and TSw,s are the local surface temperatures. The magnitude of convective heat transfer coefficient describing the resistance to heat flow across the water-pipe interface, kw−p, vary in time and in dependence on sol-air temperature, Tsol, and water mass flux,m^· , hence, providing modified convective heat transfer conditions between flowing water and internal pipe surface. Two different functions of kw−p

T_sol,m^· for the lower and top pipe were estimated on basis of CFD analyses (Section 5.2.5). The functions were implemented in the FORTRAN subroutine in order to modify the heat exchange rate at each time-increment of the FE transient simulation.

Initial conditions

The initial conditions (t = 0) for the unsteady heat transfer analysis were prescribed for the temperature distribution in the entire computational domain T (x, 0), wa-ter flow velocity, uw(0), convective heat transfer coefficient of external collector surface, he(0), and sol-air temperature, Tsol(0). Since the parameters of convective-equivalent horizontal, λHeqv, and vertical, λV eqv, thermal conductivity coefficients, heat transfer coefficients describing the resistance to heat flow across the water/pipe-down, kw−p,down, and water/pipe-up interfaces, kw−p,up, and heat transfer coeffi-cients describing the resistance to heat flow across the air/pipe-down, ka−p,down, and air/pipe-up interfaces, ka−p,up, are dependent on sol-air temperature and water mass flux, it is necessary to apply appropriate magnitudes of these parameters in the model. The initial conditions, considered for the unsteady analysis, are as follows:

• temperature distribution, T (x, 0), in the computational domain is obtained from a simulation of the steady-state heat transfer in HSC. The equations gov-erning the steady-state heat transfer (Eq. 5.3.24 and Eq. 5.3.25 ) are completed with appropriate boundary conditions (Eq. 5.3.26, 5.3.27, 5.3.28, 5.3.29, 5.3.30, 5.3.31, 5.3.32, 5.3.33, 5.3.34, 5.3.35). The magnitudes of convective heat trans-fer coefficient at external collector surface (he) and sol-air temperature (Tsol) are pre-calculated for data corresponding to the first hour of considered cli-mate database. The magnitudes of convective-equivalent thermal conductivity coefficients for the air layer (λHeqv and λV eqv) and magnitudes of heat trans-fer coefficients describing the resistance to heat flow across the water/pipe interfaces (kw−p,down and kw−p,up) and the air/pipe interfaces (ka−p,down and k_a−p,up) correspond to a pre-calculated magnitude of sol-air temperature and water mass flux equal to zero. The applied approach to determine the temper-ature distribution is reasonable, since the structure components of HSC are characterized by relatively low thermal inertia. The zero-magnitude of water mass flux comes directly from the following assumption for u (0),

• for all the executed simulations of unsteady heat transfer in HSC, the consid-ered climate database starts up from the beginning of 1^stof January. Since, it is a time when I = 0 W/m², the water flow through the pipes is stopped.

Hence, the initial condition for a water flow velocity is given as uw(0) = 0 m/s,

• convective heat transfer coefficient of external collector surface and sol-air tem-perature at t = 0, as well as dependent parameters aimed at modifying thermal properties of the air layer (λHeqv and λV eqv) and heat exchange conditions on

5.3. FE MODEL 173

inner (kw−p,down and kw−p,up) and outer (ka−p,down and ka−p,up) pipes’ surfaces, are the same as in the steady-state simulation.