1. Introduction
Building Integrated Photovoltaic/Thermal (BIPV/T) systems [
1] are some of the most successful renewable energy technologies that have been integrated into buildings. With the combined efficiency of electrical and thermal energy harvesting potential, comparably high levels of efficiency are possible when compared to their sister-type, namely, building integrated photovoltaic (BIPV) buildings. This difference in efficiency has been shown to be over 60% [
2] for the former technology to typically less than 30% [
3] for the latter technology (based on a single-junction terrestrial cell). This paper investigates the potential of computational fluid dynamics (CFD) in predicting as well as envisioning mechanisms of thermal energy extraction from BIPV to increase their overall efficiency yield. The experimental setup, as previously depicted in “A Ducted Photovoltaic Façade Unit with Buoyancy Cooling: Part I Experiment” [
4], was used as a validating reference. Therefore, this CFD simulation replicated the same process of increasing the width of the duct behind the photovoltaic façade (PV) panel, as shown in
Figure 1, while monitoring the outcome in terms of temperature harvesting and possible wind speed yields. This process paves the way for further utilization of CFD in managing complex design configurations for these types of applications.
A further increase of combined efficiency can theoretically be achieved by utilizing or adding kinetic energy resulting from buoyancy-generated wind speed.
2. Methodology
A CAD model was generated with the exact geometrical dimensions as those presented in the experiment by Elbakheit [
4] (
Figure 1). The model was then imported into ANSYS CFD Fluent 2019R1 software (ANSYS, Inc. Southpointe, PA, USA). A structured mesh was produced (
Figure 2) to simulate the buoyancy flow. This was achieved by employing the standard k-ε module with the full buoyancy effect option. The mesh was refined at the back of the PV panel with inflation layers to capture the heat transfer to the air inside the duct. A heat flux of 375 W/m
2 was assigned to the PV panel together with a laminar module, followed by the turbulent model. The radiation model with surface-to-surface option was also activated.
Operating conditions were set at atmospheric pressure under gravitational acceleration and 300 K operating temperature. The boundary condition for the inlet was set to a pressure-inlet of zero total pressure. The outlet condition was set to a pressure-outlet of zero total pressure. Therefore, there was no pressure difference between the inlet and outlet at the initialization stage. This process was used to ensure that the air velocity developed would be totally out of the buoyancy effect formed by the PV heat flux absorbed by the air in the duct.
2.1. Grid Independence Study
Numerous grids were tested to determine the mesh influence on the accuracy of the simulation, as depicted in
Table 1. It was deduced that Grid 3, which had approximately 782,000 cells, was the closest in predicting the resulting velocity at the outlet of the duct, as reported in the experiment to be 0.353 m/s with a duct width of 15 cm.
Grid Convergence Study
Further examination of the grid to determine the ordered discretization error [
5] in the CFD simulation was carried out and this has provided further confidence in the results obtained by the simulation. Roaches’ [
6] Grid Convergence Index (GDI) [
7] provides an error band on the grid convergence. The GDI can best be calculated at three levels to accurately estimate the order of convergence and to check that the solutions are within the asymptotic range of convergence.
The first step in a GDI, according to Roaches, is to determine the order of convergence observed for the predicted velocities in
Table 1, grids 1 and 2, as follows:
The second step is to apply the Richardson extrapolation for grids 2 and 3 in
Table 1.
To check if the solutions are in the asymptotic range of convergence, we used
which was approximately one, indicating that the solutions of these three grids were well within the asymptotic range of convergence.
2.2. Computational Fluid Dynamics Results
The magnitude of velocity was measured in m/s by the vertex average under surface integers. PV temperatures were measured in Kelvin by the area-weighted average under surface integers.
When simulating all 10 duct widths shown in the referenced experiment, agreement was found between the simulated and practically measured values of buoyancy-induced air velocities in the duct. This was evidently indicated by the linear trend line reported in
Figure 3. Here, both simulated and experimental results indicated that the highest resulting air velocity owing to buoyancy was approximately 0.5 m/s for a duct width of 5 cm. This velocity was gradually and linearly reduced to approximately 0.2 m/s in the largest duct width of 50 cm.
The nature of the flow is a developed turbulent flow within the duct with a somewhat streamlined flow along the back of the PV panel, as illustrated in
Figure 4a,b. This comes as an assurance to the calculated Gr in Part I of this research work.
The largest turbulence was revealed at the corners of the frame of the PV panel as well as in the point of junction between the flow from the back of the PV panel and the duct main stream (
Figure 4b).
Similarly, the simulated mean cell temperature was found to follow the same track or linear trend as revealed in
Figure 5a, this figure also provides percentage of PV performance improvement owing to reduced temperatures that ranged from 10.31% for duct-10 to 12.69% for duct-1.
Figure 5b reveals the plots of the average temperature profile of the air from the inlet to the outlet of the duct for all duct widths examined in this research. The ΔT between the inlets and outlets in both the experimental and CFD-simulated results showed a linear trend (
Figure 6a). Nevertheless, the ΔT from the experiment remained steady with little change and ranged from 5.47°K to 12.32 °K, whereas in the CFD results, there was a steeper linear trend profile that ranged from 8.1°K to 19.32 °K. There were higher temperature differences at low duct widths, i.e., low aspect ratio ducts. This can be attributed to the difficulty in the practical measurement of the temperature in the low aspect ratio ducts due to relatively higher air speeds and turbulences in those ducts. Nevertheless, higher aspect ratio ducts showed a better match between the measured and simulated results. Similar effects were experienced by Moshfegh and Sandberg [
8] and Saadon et al. [
9]. Moshfegh’s duct aspect ratios ranged from 28 to 108, whereas Saadon et al.’s duct aspect ratio was around 12. In this study, the duct aspect ratios ranged from 0.04 to 0.40.
Figure 6b shows a comparison of the cooling effect in watts by the experimental and simulated ducted system. Both trend lines exhibited higher cooling for larger duct widths than that for lower duct widths. However, the influence of changes in temperature pointed out above affected the cooling yield by making the simulated CFD have higher cooling than the experimental results at low aspect ratio ducts (i.e., ducts of 0.05, 0.10, and 0.15 m).
However, the simulated temperatures were slightly higher than the experimentally recorded temperatures, as revealed in
Figure 6a,b, especially for low duct widths. This result was also evident for the computed temperature difference between the inlets and outlets of the different duct widths. We found that the experimental results were lower than the CFD-simulated results.
In contrast, if we compare the temperature distribution along the PV surface in
Figure 7a and the velocity stream in
Figure 4b, we can deduce that lower temperatures (i.e., blue to cyan spots) formed due to jetting streamlines, whereas the high temperatures on the top corners of the PV panel were because of stagnating streamlines on those particular areas where the temperature soared to 343 K. Heat was swept from the back of the PV panel and gradually dissipated into the main stream in the duct, as revealed in
Figure 7a,b.
2.3. Computational Fluid Dynamics Comparison with Bench Mark Tests
Comparing the simulation results with the available bench mark tests of De Vahl Davis (1983) [
10] and Lage and Bejman (1991) [
11], these bench marks were confined to natural air convection in closed square cavities at Prandtl Number (Pr) = 0.71 and a cross section aspect ratio of 1, a condition that only applied to Duct-1 with a width of 50 × 50 cm in the experiment presented here. However, in the two benchmark tests, a full laminar flow only took place below Gr = 10
9. As the present experiment was about 10
11, the flow in Duct-1 was turbulent and could not be compared with a laminar flow (as revealed by
Figure 4 and
Figure 8).
3. Kinetic Wind Energy
The power yields from the generated wind speeds in the 10 ducts studied were very minor as the resulting maximum velocities were very low. However, taller ducts could improve the acceleration of wind speed [
12], a condition that is more prevalent within tall buildings. Wind power is estimated by Equation (4).
where
p is air density;
A is the cross-sectional area of the duct in this case; and
V is the resulting wind speed from buoyancy. The calculations of the potential kinetic energy from Equation (4) for the 10 studied ducts are plotted in
Figure 9. This figure revealed very low expectations for kinetic energy due to the low velocities resulting from buoyancy.
4. Conclusions
Computational fluid dynamic simulations correctly predicted the velocity within the ducts, the mean PV cell temperature, and the airflow nature, as depicted by the Grashof numbers. However, for the temperature differences, mass flow rates and cooling calculated, only the mean trend lines were somewhat in line.
Higher aspect ratio ducts (width/height) had a greater accuracy and match between the experimental and simulation resulting speeds, temperature, and temperature differences, and ultimately the cooling effects.
Taller ducts accumulate more buoyancy-induced mass flow rate, which increases the module cooling potential and increases the kinetic power yield for the system. Thus, higher combined electrical, thermal and kinetic efficiencies are realized.
Vertical ducts have the maximum benefit from buoyancy cooling. In inclined or sloped ducts, buoyancy forces are reduced by a factor of the cosine of the inclination angle.
BIPV/T systems can increase the combined efficiency of energy systems by above 60%, and with further study, this efficiency could be further increased, given that kinetic power can be estimated or optimized.