Higher Order Accurate Transient Numerical Model to Evaluate the Natural Convection Heat Transfer in Flat Plate Solar Collector

: The natural convection flow in the air gap between the absorber plate and glass cover of the flat plate solar collectors is predominantly evaluated based on the lumped capacitance method, which does not consider the spatial temperature gradients. With the recent advancements in the field of computational fluid dynamics, it became possible to study the natural convection heat transfer in the air gap of solar collectors with spatially resolved temperature gradients in the laminar regime. However, due to the relatively large temperature gradient in this air gap, the natural convection heat transfer lies in either the transitional regime or in the turbulent regime. This requires a very high grid density and a large convergence time for existing CFD methods. Higher order numerical methods are found to be effective for resolving turbulent flow phenomenon. Here we develop a non-dimensional transient numerical model for resolving the turbulent natural convection heat transfer in the air gap of a flat plate solar collector, which is fourth order accurate in both spatial and temporal domains. The developed model is validated against benchmark results available in the literature. An error of less than 5% is observed for the top heat loss coefficient parameter of the flat plate solar collector. Transient flow characteristics and various stages of natural convection flow development have been discussed. In addition, it was observed that the occurrence of flow mode transitions have a significant effect on the overall natural convection heat transfer.


Introduction
Solar collectors trap the solar radiation and convert it into thermal energy which has many applications, such as cooking, indoor space heating, hot water generation, distillation, steam production, solar air conditioning, etc. The conversion efficiency of solar energy into thermal energy at required temperatures should be as high as possible. This aspect requires understanding of the heat loss mechanism and various heat transfer mechanisms inside the solar collectors. All three heat transfer mechanisms, (conduction, convection and radiation), occur simultaneously inside the solar collector. The main heat losses in a Flat Plate Solar Collector (FPSC) are due to heat loss by convection in the air gap between the absorber plate and the glass cover of the FPSC, which constitutes approximately 70% of the total heat loss. The natural convection phenomenon inside this air gap lies in the turbulent regime due to a large temperature difference between the absorber plate and glass cover of an FPSC. Several researchers have conducted experimental studies on solar collectors to optimize various design parameters to achieve maximum conversion efficiency while minimizing losses. Similarly, researchers have also developed analytical models to analyze the heat transfer phenomenon inside an FPSC. Duffie and Beckman [1] developed a steady state 1D model to evaluate the efficiency of an FPSC, which is a function of fluid temperature at the inlet. Nusselt number correlations have also been proposed to evaluate the heat transfer phenomenon inside the solar collector at a steady state. However, heat transfer during experiments lie mostly in the transient state since the incident solar radiation itself is continuously varying in nature. Therefore, the requirement of achieving steady state conditions during experiments to compare with the analytical results is impractical and uneconomical. Heat transfer analysis of an FPSC can be accurately analyzed by its transient response to input parameters, such as solar radiation intensity, weather variation, shadow factors, inlet fluid temperature, etc., which are transient in nature.
In order to study the transient nature of the heat transfer phenomenon in an FPSC, lumped capacitance models have been developed [2,3]. Lumped capacitance assumes there are no significant temperature spatial gradients along the spatial dimensions of the solid body. Close D. [4] has implemented a one-node model initially by implementing energy balance equations on the absorber plate surface. It is assumed that the heat transfer fluid, absorber plate and glass cover have the same temperature. Later, a two-node model was developed by Klien S. [5], in which the temperature of the absorber plate and glass cover are different, resulting in two energy balance equations at each of these nodes. However, the absorber plate and heat transfer fluid have a common node. Later, Morrison and Ranatunga [6] proposed a three-node model in which the heat transfer fluid, absorber plate and glass cover have a separate node in which the energy balance equations are created for each of these nodes. Here, the temperature of the heat transfer fluid from the inlet to the outlet is assumed to vary linearly. Subsequently, multi-mode models were developed having multiple nodes in the absorber plate, heat transfer fluid, glass cover, etc. [7].
Lumped capacitance models do not consider the spatial temperature gradients. In the case of natural convection flow inside the air gap of the FPSC, the temperature of the entire air gap is assumed equal, which is not the case in practical situations. In fact, it is the variation of the temperature in this air gap that offers the understanding of how the heat transfer occurs through natural convection. Therefore, the modelling of the fluid flow phenomenon in this air gap and understanding of the temperature and flow variation phenomenon both spatially and temporally is very important to understand the heat loss mechanism and thus increasing the conversion efficiency.
To overcome the drawbacks of empirical correlations and lumped capacitance models, researchers have resorted to numerical methods to evaluate the heat transfer distribution of various heat transfer mechanisms in the FPSC. With the advent of numerical modelling techniques fueled by fast and powerful computers, researchers explored the possibility of modelling entire solar collectors numerically [8][9][10]. Recently, with the advances in the field of computational fluid dynamics, complete modelling of solar collectors has become possible [11][12][13]. Steady state simulations of 2D and 3D FPSC models to simulate the fluid flow behaviour are studied mostly in Fluent and other commercial software [14][15][16]. However, the FPSC performance is analysed on specific collector prototypes and compared with the experimental results of these prototypes only. Moreover, these models are 2 nd order convergent requiring a very long convergence time and a very high grid density to model any practical solar collector, as the convection heat transfer in experimental solar collectors lies in the transitional or turbulent regime.
A transient numerical model that characterises the natural convection flow behaviour within the annulus gap of the FPSC is not available in the literature. As well, it is known that heat transfer in a FPSC lies in the turbulent regime, we develop a numerical method of 4th order accuracy to completely model the heat transfer behaviour of the solar collector and also to analyse the transient fluid flow behaviour inside the air gap of the FPSC.

Methodology
In Sections 2.1 and 2.2, a higher order numerical method to evaluate the transient natural convection flow in the turbulent regime is developed. Further, the method is validated against benchmark case studies to prove their applicability to the method for turbulent simulations. Further, in Section 2.3, the developed method is extended to model the natural convection flow behaviour of solar thermal collectors and validated against the empirical correlations to evaluate the top heat loss coefficient available in the literature.

Governing Equations for Natural Convection Flow in Enclosures
Here we develop a 4th order numerical method for transient natural convection flow in the turbulent regime. The governing equations for fluid flow i.e., Navier Stokes equations in primitive variable form are The properties ρ, ν, α, β are assumed to be constant, except for the density variation, which is considered only in the buoyancy force through the Boussinesq approximation given by The system of partial differential equations in primitive variable form in Equations By applying ' ' on Equation (2) the vorticity transport equation is obtained [17], as given in Equation (8). The transformed equations are Equations (7)-(9) are non-dimensionalised using the following dimensionless quantities Pr The following set of Non-Dimensional (ND) equations are derived using the dimensionless quantities stated in Equation (10).
The governing Equations (11)-(13) are discretized using the fourth order discretization method in both space and time. Wide stencil discretization of fourth order accuracy is applied for spatial differential terms. The fourth order Runge-Kutta method (RK4) is adopted for explicit discretization of the transient terms. The described method is completely fourth order accurate up to the boundary in both space and time. Achieving a still higher order accuracy is possible in both space and time by applying further higher order schemes, as the method is completely explicit iterative in nature. The stability of the method is ensured by the enhanced stability region of the RK4 scheme adopted for the transient term discretization. The complete fourth order method to evaluate the transient behaviour of natural convection flow in enclosures is described in detail in [18].
The algorithm for implementing the numerical method is written using Matlab Programming language and Matlab software is used to run the developed numerical model. The algorithm is run on a mobile workstation HP Z-book with an 8 GB ram and 8 core parallel processor (Intel(R) Core(TM) i7-4900MQ) with a clock frequency of 2.82 GHz.  Figure 1 shows the case study of Differentially heated square cavity problem where the fluid flow is due to the buoyancy force generated due to the temperature difference between the vertical walls of the square cavity. The fluid flow behaviour in the enclosure is defined by governing Equations (7)-(9) in dimensional form and in Equations (11)- (13) in non-dimensional (ND) form. The vertical walls are differentially heated while horizontal walls are insulated with boundary conditions as shown in the figure. Due to no-slip BC on the walls, the velocity components on the walls are zero and the fluid inside the cavity having Pr = 0.71 is at rest initially. This is a benchmark problem used for assessing the performance of numerical techniques for solving incompressible Boussinesq equations for simulating the natural convection flow. A comprehensive paper providing benchmark results comparing the numerical solutions of various authors was published by De Vahl Davis [19]. The authors proposed benchmark results in this range of 10 3 to 10 6 by summarising the numerical results from 37 sources for the case of Pr = 0.71. Later, Le Quéré [20] has extended the range of Rayleigh numbers from 10 6 to 10 8 which is also widely considered as the extended range of benchmark solutions for transitional to turbulent regime natural convection flow available for the above model problem. In the successive Tables 1 and 2, we compare the present steady state solutions with benchmark solutions obtained by De Vahl Davis [19] and Le Quéré [20] for Ra 10 6 -10 8 and some well-established numerical results available in the literature [20][21][22][23][24][25][26][27]. Grid independence tests have been carried out and it was found that for the case of Ra = 10 7 , the 2D grid configuration is 161 × 161 nodes, while for the case of Ra = 10 8 , the 2D grid configuration of 241 × 241 nodes gave optimum converged results. The results obtained are in excellent agreement with these benchmark results. These results confirm the validation of the present method to evaluate the large range of Ra flows in the turbulent regime without any stability problems due to the inherently stable nature of the 4th order Runge Kutta method (RK4 scheme). The streamlines, isotherms and vorticity contours over the entire range of Rayleigh numbers (10 7 to 10 8 ) are shown in Figure 2. The sharp gradients in the flow parameters observed within the enclosure confirm the turbulent behaviour of the flow.

Buoyancy Driven Convection in Rectangular Cavity
This problem is selected to study the transient behaviour of convection heat transfer and fluid flow at a turbulent scale. As pointed out by Christon et al. [28], this model problem serves as a mechanism to test the performance of the numerical methods to simulate the complex physical mechanisms, such as turbulent scales, several instability mechanisms, vertical & horizontal boundary layers, travelling waves in vertical boundary layers, thermal instabilities along horizontal walls, in particular which can strongly interact with internal wave dynamics. Most importantly, this model problem exhibits an oscillatory transient flow behaviour above a critical Rayleigh number Rac = 3.4 × 10 5 . Christon et al. [28] has methodically summarised the results from a special session to generate a timedependent benchmark solution for an 8:1 differentially heated cavity at Pr = 0.71 and Ra = 3.4 × 10 5 , as shown in Figure 3. The x and y coordinate locations of the point numbers 1 to 5 are shown in Table 3. These results serve as the benchmark solution for us to test the transient behaviour of the developed numerical method in the 8:1 differentially heated cavity. For further details the readers are referred to Reference [28].   (17) where T represents the period of time for which the average was computed. The oscillatory component was to be computed as After conducting grid independence tests, the optimum grid resolution employed for evaluating this case study is 41 × 201 nodes. The average is displayed for one or more complete periods where the amplitude and period were essentially constant, i.e., after startup transients completed. For e.g., in Figure 4 the startup transients exist up to 400 ND time units, after which amplitude is nearly constant.  Table 4. The point data consisted of velocity components (u, v), temperature θ, stream function ψ, and vorticity Ω at time-history point-1.   This time periodic behaviour matches with the results available in the literature. Also the oscillatory behaviour is confirmed in Figure 5, which shows the variation of the ND X velocity component and ND temperature variation with respect to 1 to 5. It is observed that a stable oscillating behaviour is observed after 400 ND time units approximately. In addition, it is observed that with reference to point 5, the reference points 1,2 and reference points 3,4 are symmetric in nature. Moreover, the time averaged and fluctuating components as described in Equations (17) and (18) are given in Table 4, which are compared with various benchmark results available in the literature. These accurate results prove that the present numerical method could accurately simulate the transient behaviour of natural convection flows in enclosures.

Extension of the Developed Numerical Method to Solar Thermal Collector
FPSC is divided into 3 zones as shown in Figure 6. The plate domain (p), which designates the absorber plate section of the FPSC, the fluid domain (f), which is filled with air and the glass domain (g).

Governing Equations
The governing equations which define the heat transfer behaviour in each section is presented below.
The governing equations are non dimensionalised using the following dimensionless quantities.

Ra Gr Pr
Boundary conditions:

Numerical Implementation
The following Algorithm 1 is employed to numerically evaluate the heat transfer performance of Flat plate solar collector.

Validation
Mullick and Samdarshi [29] have proposed an empirical correlation to evaluate the top heat loss coefficient of a solar collector [30]. The range of variables covered for which the model empirical correlation was applicable is 50 °C to 150 °C of the absorber plate temperature, 0.1 to 0.95 in the absorber coating emittance, 5 W/m 2 C to 45 W/m 2 C in the wind heat transfer coefficient, 20 mm to 50 mm in the air gap spacing, 0 deg to 70 deg tilt angle, and 10 °C to 40 °C ambient temperature. By varying these variable within these ranges the Rayleigh number is calculated based on the assumption that the temperature difference between the absorber plate and glass cover is at least 20 °C. For various Ra, the top heat loss coefficient is evaluated and compared for the present numerical model and the Mullick & Samdarshi empirical correlation [29] in Table 5. It is observed that as Ra is increased, the error reduces and is found within 5% at an even higher Ra in the turbulent regime as well.

Results
The fluid flow phenomena and occurrence of flow mode transition inside the annulus air gap of the solar collector are presented. All the simulations results are plotted after carrying grid independence tests with an optimum grid density of 240 × 12 nodes resolution. Transient iso-thermal contours on the computational domain of the numerical model are shown in Figure 7 for Ra = 10 5 . Specifically in Figure 7d, the iso-thermal contours in the plate domain, fluid domain and glass domain are shown. In the following figures we have shown only the iso-thermal contours to clearly depict the occurrence of flow mode transition. As observed in this Figure, thermal plumes are formed where two stream function contours rotating in opposite directions intersect asymptotically. As well, the flow mode transition shown in Figure 8, where the mixing of 2 thermal plumes happens when the associated stream function contours coalesce into a single larger stream function contour. Similar fluid flow behaviour is also observed in the annulus gap of solar collectors in a previous CFD study [12]. The total transient behaviour of the convection fluid flow inside the annulus air gap can be classified into 3 stages.   Typically for low Ra flows, where Ra less than critical Rayleigh number Racr, fluid is completely stable and flow does not occur, i.e., the heat transfer occurs completely due to conduction through air layers and the resultant Nu is zero. As Ra increases just greater than Racr, the gravitational forces become more dominant, instability sets in and Rayleigh-Benard convection cells appear. This phenomenon is called Rayleigh-Benard instability. This is when convection heat flow starts and the temperature starts rising from the absorber plate towards the glass cover, as shown in Figure 7a. The flow developmental stage is the stage starting from the occurrance of Raylieigh-Benard instability to the occurance of 1st flow mode transition. Flow mode transition stage is defined as the time period during which flow mode transitions occur and finally the steady state heat transfer is reached. As observed in the isothermal contours for low Ra at Ra = 10 5 , the flow development stage has an approximate time period of 1.2 ND time units, as seen from the Nu diagram when the first flow mode transition occurs, followed by the flow mode transition stage where multiple Rayleigh Benard convection cells coalesce into single cells. The time period for the flow mode transition stage is approximately from 1.24 ND time units to 2.35 ND time units. As seen in Figure 8, 2 flow mode transitions occur at 1.24 s and 2.35 s respectively. As seen in Figure 9 left, the flow mode transition is also characterised by a decrease in Nu from 6.1 to 5.8 for the 1st flow mode transition and further decreases from 5.8 to 5.2 in the 2nd flow mode transition. The 3rd stage of transient flow is the fully developed stage where heat is transferred from the absorber plate to the glass cover via fluid flow configuration after the completion of all flow mode transitions, as shown in Figure 8d, for the case of Ra = 10 5 . It is in the fully developed stage when the bulk of the heat transfer takes place from the absorber plate to the glass cover.
For the case of Ra = 10 6 , the flow development stage and flow mode transition stages are instantaneous, as can be verified from isothermal contours in Figure 10. The flow mode transition occurs below 1 s and the bulk of the heat transfer occurs through fluid flow configurations, shown after 1 s. In addition, as seen in Figure 8b-d, the thermal plumes continuously oscillate after the developmental stage, resulting in the oscillating nature of the Nu value seen in Figure 9 right. This behaviour is also confirmed for higher Ra flows for triangular enclosures [31]. As Ra increases, the flow mode transitions happen instantly and the flow oscillates with a periodic behaviour, indicating a turbulent flow phenomenon at a higher absorber plate temperature with reference to the ambient or glass cover temperature.  Another important feature of increasing Ra is the reduction in both velocity and thermal boundary layer thickness of fluid flow contours. As observed in Figures 7 and 10, as Ra increases, the boundary layer thickness near the solid surfaces and also the thickness of the thermal plumes significantly decrease. This is also associated with a corresponding increase in the size and strength of the stream function contours and also steep temperature gradients are formed in the boundary layers. In all the cases studied here, with an increase in Ra, the boundary layer thickness decreases and the Nu increases.

Conclusions
A fourth order accurate numerical method in both space and time is developed to evaluate the turbulent and transient behaviour of natural convection flow inside the air gap of a FPSC. The developed method could accurately model the steady state response of the solar collector with an error limit of 5% when compared with Mullick & Samdarshi [29]. It was observed that the transient state of fluid flow behaviour inside the air gap exists in three stages i.e., flow development stage, flow mode transition stage and fully developed stage. It was observed that flow mode transitions are associated with a reduction of Nusselt numbers and the associated convection heat transfer. Also flow mode transitions are preferable for reducing the convection heat transfer and subsequently the convection heat loss and thus improving the efficiency of the FPSC. Overall, the developed method could accurately predict the transient behaviour of fluid flow inside the air gap of the FPSC and offer new mechanisms to optimise and reduce the convection heat transfer loss in the solar collectors, and thus improving the overall conversion efficiency of the solar collectors.

Conflicts of Interest:
The authors declare no conflict of interest.