Theoretical and Numerical Study of a Photovoltaic System with Active Fluid Cooling by a Fully-Coupled 3D Thermal and Electric Model

: The paper deals with the three-dimensional theoretical and numerical investigation of the electrical performance of a Photovoltaic System (PV) with active ﬂuid cooling (PVFC) in order to increase its efﬁciency in converting solar radiation into electricity. The paper represents a reﬁnement of a previous study by the authors in which a one-dimensional theoretical model was presented to evaluate the best compromise, in terms of ﬂuid ﬂow rate, of net power gain in a cooled PV system. The PV system includes 20 modules cooled by a ﬂuid circulating on the bottom, the piping network, and the circulating pump. The fully coupled thermal and electrical model was developed in a three-dimensional geometry and the results were discussed with respect to the one-dimensional approximation and to experimental tests. Numerical simulations show that a competitive mechanism between the power gain due to the cell temperature reduction and the power consumption of the pump exists, and that a best compromise, in terms of ﬂuid ﬂow rate, can be found. The optimum ﬂow rate can be automatically calculated by using a semi-analytical approach in which irradiance and ambient temperature of the site are known and the piping network losses are fully characterized.


Introduction
In order to improve the efficiency of photovoltaic (PV) cells, several research studies have recently been carried out to investigate, both theoretically and experimentally, the possibility of increasing the maximum power by reducing the cell temperature, thanks to the adoption of cooling systems [1][2][3]. In fact, the behavior of the PV efficiency can be confidently considered linear with respect to its temperature [4][5][6] and, for crystalline silicon technology, efficiency variation coefficients with cell temperatures range from −0.4%/ • C to −0.5%/ • C [1][2][3][4][5][6].
Typically, active and passive cooling systems [7] are adopted to reduce PV cell temperature, and several designs have been proposed in literature [8][9][10][11][12][13][14][15]. The main advantage of active cooling systems is a more efficient heat removal rate; for this reason, they are generally preferred to passive systems, in which a cooling fluid circulates without needing additional power. In this case, active cooling systems present an additional advantage related to the possibility of supplying an

The Physical Model: The PVFC Module and the Case Study of the Hybrid Photovoltaic-Thermal Plant
The system under investigation has been already described in previous works by the authors [24,25]. In particular, in [25], the single photovoltaic system with active fluid cooling (PVFC) was studied both theoretically and experimentally by using a one-dimensional thermal model. In the following, a short description is reported and further details can be found in [25].
Regarding the cross-section of each PV module, the sandwich is composed of seven layers: A polycarbonate (PC) front-sheet without glass, a first thermoplastic polyurethane (TPU) sheet, a thin layer containing 40 series-connected monocrystalline silicon (m-Si) cells (4 columns of 10), a second thermoplastic polyurethane sheet, an alveolar polycarbonate (PC) sheet, an air gap, and, finally, an alveolar polycarbonate crossed by the coolant fluid. In Table 1, the thermal properties and thickness are reported. Table 1. Thermal properties and thickness of the photovoltaic system with active fluid cooling (PVFC) module. More details can be found in [24]. The layout of the PVFC plant, widely described in [24], is reported in Figure 1 and consists of: Two PVFC strings, each composed of 10 modules (maximum power of 80 W at Standard Test Conditions (STC) [26] and 40 cells/module) for a total peak power of 1.6 kW, a piping network, a fluid tank (WT), and a circulation pump (P). The geometrical data of the piping network are reported in [24], where L 1 represents the pipe length from the WT outlet to the inlet of the first PVFC module of the bottom string; L 2 is the pipe length from the PVFC module, mentioned above, to the first PVFC module of the top string; L 3 is the pipe length from the last module of the top string to the last module of the bottom string; L 4 is the total length of the alveolar-polycarbonate ducts, while the total height of the PV module is 1.22 m; L 5 is the pipe length from the last module of the bottom string to the WT; L 6 is the heat exchanger length inside the WT. In particular, in each module, the coolant fluid is a glycol/water mixture (10%/90%) flowing into 60 alveolar polycarbonate ducts, each of them of L 4 = 1.42 m (the total number of ducts, N d = 1200), while the length of the helicoidal coil heat exchanger (HX) inside the tank is L 6 = 8 m. Figure 1. Layout of the PVFC plant, composed of 20 photovoltaic (PV) modules with a total peak power of 1.6 kWp. The cold water inlet (blue pipe) and hot water outlet (red pipe) are reported.

The Three-Dimensional Thermal-Electric Fluid Dynamic Model
The investigation of the PVFC plant was carried out by means of a three-dimensional model, where the thermal model of each single module was fully coupled with the electrical model, as the PV electric efficiency is dependent on the spatial distribution of the cell temperature. For each module, 40 PV cells were considered, and an average cell temperature was consistently calculated by solving the heat transfer equations with the surrounding ambient and with the coolant fluid. Piping network losses were accounted for in order to calculate the net electric power gain. The assumptions made are reported in the following.
About the thermal-fluid dynamic model: • The heat transfer through the PV module is in steady-state conditions; • the thermal conductivities of the layers are constant; • the thermal contact resistances at the interfaces are negligible; • the cell thickness is negligible; • the solar absorptivities for the thermoplastic polyurethane (TPU) and polycarbonate (PC) layers are negligible; • the working fluids are considered as incompressible.
About the hydro-dynamic model for the pressure losses: • The total flow rate is uniformly distributed between the 20 PV-T modules; • pipes are considered smooth and frictional losses are given by the Darcy equation, both for laminar (friction factor given by the Darcy equation) and turbulent flows (friction factor given by the Blasius expression).
Finally, about the electrical parameters: • The electrical parameters, i.e., the series resistance, R s , the shunt resistance, R sh , the spectral-response parameter, K ph , and the diode quality factor, m, were considered constant.
In the following, more details about the hypotheses and the geometries of the PV modules and of the plant are given.

The PVFC Module: Structural, Thermal, and Electrical Model
Each module, as described in [25], was modeled considering seven layers, as reported in Figure 2. The thermo-fluid dynamics of the PVFC module were investigated by solving continuity and energy equations for the solid and continuity, momentum, and energy equations for the coolant fluid. Each layer differs in thickness, density, thermal conductivity, and specific heat capacity according to the geometries and properties given in [25]. In particular, the negligible m-Si cell thicknesses and negligible absorption coefficients of the TPU and PC layers were considered. Moreover, radiative thermal exchanges between the polycarbonate layers and the external environment were included. Polycarbonate absorbance and emissivity are considered negligible in the visible spectrum, while, in the infrared, their behaviors are totally different, as the surfaces have a high absorbance and high emissivity: The PC layer absorbs the heat of the lower layers, which are conductive, and exchanges heat with the atmosphere by radiation. The boundary thermal conditions (heat flux, temperature, convection, radiative, and conductive transmission), absorbance, and transmittance coefficients are considered according to [25], and inlet fluid conditions are imposed both in mass (kg/s) and in temperature (K).
The alveolar polycarbonate sheet where the coolant flows is considered as a single channel instead of 60 alveolar channels. This assumption is valid in the framework of the thermal-electrical model and does not significantly affect the spatial distribution of the cell temperature. Instead, regarding the pressure losses, all 60 of the alveolar polycarbonate ducts were accounted for, each of them of L 4 = 1.42 m. Moreover, the total flow rate is uniformly distributed between the 20 PV-T modules, considering a balanced piping network characterized by one common inlet and one common outlet. In PVFC plants with non-uniform flow rate distributions, the theoretical and numerical models can still be applied by evaluating the cooling effects on each module.
The module was discretized by considering unstructured meshes for the m-Si cells and for the cooling fluid and structured meshes for the other components of the module. In particular, a grid sensitivity study was carried out to demonstrate grid-independent and grid-convergent results, checking both solutions and residual values. In fact, in order to guarantee the solution accuracy, numerical results obtained over a range of significantly different grid resolutions were obtained. Relative errors of cell and fluid temperature fields and of integrated quantities such as average electric efficiencies were calculated. The final grid selected (690,366 nodes and 408,960 elements) guarantees the solution accuracy (maximum relative error of less than 2%). Figure 3 shows the scaled residuals [27] of the continuity, momentum, and energy equations of a simulation when the steady-state solution is reached for a case study during July in Catania (Southern Italy), with orientation SOUTH and tilt 33 • . Regarding the modeling of the solar source, ANSYS/FLUENT allows one to choose between different options by using the solar calculator function and selecting the geographical position (latitude, longitude, and timezone) and the day of the year. In this paper, the solar ray tracing model [27] was adopted by selecting the direct solar irradiation, which accounts for both infrared and visible band radiation, and the diffuse solar component, which accounts for the hemispherical averaged spectrum. The values for both type of irradiation are given according to the PVGIS database.
The thermal model was fully coupled with the electrical model according to the algorithm proposed in [25] with the main difference that, in 3D simulations, cell temperature spatial distribution is not averaged in the transverse direction. The source term for the energy equations of the 40 series-connected m-Si cells depends on the electrical performances of the PV cells, which, in turn, are related to the cell temperature calculated by solving the heat transfer equations. By following this approach, a rigorous calculation of the spatial distribution of the cell temperature is carried out, considering the influence of the solar irradiance G and the ambient temperature T a . Typically, simple linear formulas are used to evaluate the cell temperature when irradiance and ambient temperature are known and spatial distributions are neglected. The model, inherently nonlinear, can be solved iteratively until the steady-state solution is obtained. The PV module was modeled in a three-dimensional geometry (x, y, z), where z represents the coordinate varying along the thickness of the module and (x, y) represent the coordinates of the different layers. In particular, the electrical efficiency, η(T c (x, y)), and the source term of the energy equation of the m-Si cells, Q C (x, y), are coupled as follows, considering that the fraction of the irradiance η(T c (x, y))G is converted into electricity: where τ PC 1 τ TPU1 τ ARC is the product of the transmissivity of the first three layers (first polycarbonate sheet τ PC 1 = 0.93, first thermoplastic polyurethane sheet τ TPU1 = 0.90, anti-reflective coating τ ARC = 0.85),s c is the cell thickness and η(x, y) represents the cell performance, as follows: An UDF was developed to calculate the internal heat generation of the of the m-Si cells and the PV efficiency, and the external routine was fully coupled with the CFD solver. The electrical efficiency was locally considered linearly dependent with respect to its temperature. For each cell, an average (spatial) temperature was calculated in order to evaluate the electrical efficiency by using the Shockley one-diode model, in which it depends on the solar irradiance G and the ambient temperature T a . The current I and the voltage V are related by the following equation [28]: where I ph = K ph GA is the photogenerated current, I D = I 0 e qV D mKTc − 1 the diode current, V D the diode bias, I 0 the saturation current, R s the series resistance, R sh the shunt resistance, m the diode quality factor (between 1 and 2), and T c the effective temperature of the cells. In this case, the slight dependence of electrical parameters on the temperature was neglected [29,30].
The thermal model, implemented in ANSYS/FLUENT, calculates the spatial distribution of the temperature of each cell, while the electrical model, by using the UDF, calculates the electrical performance of the cell (i.e., current I, voltage V, and maximum power of each cell), including the efficiency of each PV cell. Due to the non-linearity of the model, the numerical procedure to evaluate electric performance, cell temperature, and energy source terms is iterative. When the steady-state solution is obtained, the model can calculate all of the hydraulic losses involved in the piping network in order to evaluate the power consumption of the circulating pump. The competitive mechanism related to these opposite effects is characterized by a maximum value of the net power gain, i.e., the power gain due to the cell temperature reduction less the power absorbed by the circulating pump, and can be reached at the best compromise of the fluid flow rate. Further details of the geometry, properties, and of the hydro-dynamic model are given in the next section and in [24].

Pressure Losses in the Piping Network
The power consumption of the circulating pump is strictly related to the design of the piping network, including all of the components producing local flow losses (valves, bends, elbows, inlets, exits, tees, etc.). The present study was applied to the layout represented in Figure 1, but it can be extended to any plant when a detailed analysis of the channel geometry and a complete list of elements producing local pressure losses is carried out. In this case, the circulating pump can be selected and the pressure drop can be analytically found as function of the flow rate [31].
The power absorbed from the circulation pump P cp and partially transferred to the fluid is: where ρ is the fluid density, g is gravitational acceleration,v the volume flow rate, η P the pump efficiency, and H tot the total head losses (in meters of fluid) given by the sum of the local flow losses, H L , the frictional losses, H R , and the static height of the system, H, as follows: The local flow losses H L are experimentally obtained by the manufacturers of the components [32][33][34][35][36], and are given by the following equation: where z L is the (local) loss coefficient and w is the mean fluid velocity. The frictional losses H R are related to the roughness of the inner surface of the pipe and can be computed by resorting to the Darcy equation, both for laminar and turbulent flow: where L is the total length of the pipe, D i is the inner diameter of the pipe, and f is the Darcy friction coefficient, which depends mainly on the Reynolds number Re = wLν −1 . In fact, in laminar flows (Re < 2300) for a circular pipe, the friction factor f is related to the Reynolds number through the Darcy expression ( f = 64Re −1 ), while for turbulent flows (Re > 4000) and for smooth pipes, f is given by the Blasius expression ( f = 0.3164Re −0.25 ) [32]. Typically, the pressure drop of the exchanger is given by the manufacturers, but can also be evaluated when the material and geometry is known. In this work, a commercial helical HX was properly selected for the problem under investigation, and the pressure drop, given by the manufacturer, was considered.

The Competitive Mechanism and the Optimum Flow Rate
Flow rate affects both the power absorbed from the circulation pump P cp and the power gained due to the reduction of the cell temperature, PG, given by where P R is the power produced in real conditions when the cooling is active and P NC is the power without active cooling. Both the power P cp and the power gain PG increase with the flow rate, and a competitive mechanism arises for the net power gain NPG = PG − P cp . In this framework, the flow rate corresponding to the maximum value of the net power gain can be calculated by recurring to a semi-analytical approach.
In fact, when the piping network is fully characterized, the electrical power absorbed by the circulating pump, P cp , can be expressed in analytical form as a function of the flow ratev as follows: where N is the number of coefficients in the formula of P cp , and c k are the coefficients which depend on the geometry of the circulation loop system and the physical properties of the coolant (details are given in [24]). The power gain PG is strictly related to the spatial distribution of the cell temperature, which can be obtained by solving the thermo-fluid dynamic equations of each PVFC module of the plant. In fact, the principal parameter influencing the electrical efficiency is the cell temperature, which depends on the irradiance G, the ambient temperature T a , and the flow ratev. When the irradiance and the ambient temperature are known, the net power gain can be calculated as a function of the flow rate, suggesting the activation of the cooling system if positive. Moreover, the maximum gain during the average day of a month can be calculated in a fixed site from statistical data (i.e., PVGIS database) and can be obtained for different values of the flow rate. A polynomial expression of the maximum ideal gain as a function of the flow rate can be obtained as follows: where N PG is the polynomial degree which gives the best approximation of the analytical function PG, and a k are the coefficients which multiply the flow ratev. Finally, by combining Equations (8) and (10), the analytical expression of the net power gain NPG can be obtained: and the maximum value of the NPG can be calculated. The optimum flow rate obtained by using this approach is based on the available statistical data of the time evolution of irradiance and ambient temperature of a site. For each average day of the month, the maximum value of the net power gain can be evaluated at different flow rates and the maximum value can be obtained. In the next Section, some results are reported, highlighting both the opportunity of activating the cooling system and the determination of the maximum value of the net power gain. Finally, a discussion of the real-time evaluation of NPG is carried out by recurring to irradiance and ambient temperature measurements.

Numerical Results
In this section, a set of simulations were performed in order to investigate the behavior of the PVFC plant described in the previous sections and represented in Figure 1. The PV modules were modeled by using ANSYS/FLUENT according to Sections 2 and 3 and to the geometry, the electrical characteristics, and the thermal properties given in [24,25]. In particular, 40 solar cells are placed in ten rows made of four cells per row and are connected in series (Figure 2). The coolant flows into the duct and into the piping network with the hypothesis of laminar flow inside the PVFC modules and turbulent flow from the fluid tank to the inlet of the PVFC modules and from the outlet of the PVFC modules to the fluid tank. The simulations were performed during an average day of January and July in Catania (Southern Italy), with orientation SOUTH and tilt 33 • by using statistical average climate data and with different flow rates (in the range 0-1000 L/h). The results, obtained with the three-dimensional model presented in this work, were compared with the one-dimensional model in [24,25] and with experimental values. Figure 4 shows the cell temperature of a single module obtained in steady-state conditions forv = 0 and 1000 L/h in Catania (July). Active cooling drastically reducees the cell temperature by presenting a two-dimensional spatial distribution which affects the electrical efficiency in a sort of cooling thermal mismatch. In particular, Figures 5-7 show the spatial distribution of the cell temperature and of the electrical performance obtained for three different values ofv = 0, 600 and 1000 L/h. The results show the effectiveness of the active cooling in improving the electrical performance by reducing the cell temperature. In order to evaluate the energy source term for the thermo-fluid dynamic equations, the average cell temperatures and electrical performance were calculated for each cell in the ANSYS/FLUENT user-defined functions, and the spatial distribution of the source term was self-consistently evaluated according to Equation (1).
Comparisons of the results obtained by solving the three-dimensional model with the one-dimensional model presented in [25] and with outdoor experimental tests are reported in Figures 8-11. In particular, outdoor experimental tests were carried out in Turin (Northern Italy) in October. In these tests, after an accurate cleaning of the PV module to avoid the dirt and the dust, the thermal equilibrium condition was reached for some minutes with constant irradiance, ambient temperature, and wind speed. The experimental electrical characteristics were obtained by the automatic data acquisition system described in [25], with uncertainties of measurement of ±1%, ±0.1%, and ±1% in terms of current, voltage, and power, respectively. In particular, 1D and 3D numerical results were compared with experimental steady-state test conditions when the cooling was active (flow rates of 30 and 52 L/h for each module) and at G = 1023 W m −2 , T a = 287.45 K, and G = 700 W m −2 , T a = 288.15 K. In both conditions, the measured water temperature T in = 294.45 K. Figures 8 and 9 show the I-V and P-V curves, respectively. From the results, the agreement with the experimental values of the 3D model is improved with respect to the 1D model. The agreement is mainly related to the more accurate evaluation of the spatial distribution of the cell temperature. In the 1D model, a single average column constituted by ten rows, each of them representing four PV cells, is considered, and the thermal gradient on the transverse direction is neglected. On the other hand, in the 3D model, the complete 3D cell temperature field is considered and lower cell temperatures can be found, as confirmed in Figure 10.
The cell temperature profiles obtained forv = 200 and 1000 L/h by solving 1D and 3D numerical models are in good agreement, especially for higher flow rates and, in particular, with the profile in the middle axis (y = L y /2). The boundary effects and the spatial distribution of the cell temperature produce a thermal mismatch between the cells in the same module. In fact, the one-dimensional model does not take into account the y-dependence of the cell temperature and overestimates its value of external columns, as confirmed by the measurements reported in Figures 8 and 9. Figure 11 shows the spatial distribution of the cooling fluid temperature in the case ofv = 600 L/h in Catania (in July). The results of the three-dimensional model are reported at different z in order to show the thermal gradients inside the ducts, and are compared with the temperature profile obtained with the one-dimensional model reported in [25]. The agreement between the results obtained by solving the 1D and 3D numerical models for the casev = 600 L/h is good.
The improvement of the energy production of a PVFC system with active cooling was investigated in different conditions [25]. The difference between the no-cooling PV system and PVFC system represents the power gain PG. When the cooling is active, the pump consumption can be calculated analytically by recurring to Equation (9), where the exponents e k depend on the piping network under investigation (in this work, the system described in [24] was considered, and e k s are: e 1 = 1, e 2 = 2, e 3 = 2.75, and e 4 = 3). When the power absorbed by the circulating pump is known, the model determines at what time the cooling system must be switched on and off, and the net power gain, given by Equation (11), can be calculated.      The net power gain depends on the power consumption that, for the case considered here, ranges between 65 W for 600 l/h to 130 W at 1000 L/h [24]. By using a statistical approach for climate data (i.e., PVGIS), an optimum coolant flow rate can be estimated for every month by calculating the maximum of the net power gain NPG = P − P p , obtained by combining Equations (8) and (10). Figure 12 shows that there is an optimal value of coolant flow rate for the NPG, which represents the best compromise inside the competitive mechanism between the pump consumption and the power gain due to the cell temperature reduction. The results obtained in Catania city show that, in July, the optimum flow rate is reached at around 400 L/h (20 L/h for each module), while it is reached at around 200 L/h in January due to the lower ambient temperature.
In a real case, the data of the average day can be considered only to evaluate the profitability of installing an active cooling system. Climate conditions can be strongly variable over a single day, and the activation of the cooling system must be evaluated considering irradiance and ambient temperature on the PV generator site. Simulations can be performed to obtain a map where, starting from the knowledge of irradiance G and ambient temperature T a , the convenience of activating the cooling system in real-time conditions can be deduced [24].

Summary
In this work, a three-dimensional theoretical and numerical investigation of the electrical performance of a photovoltaic system with active cooling was carried out. Starting from results obtained in previous works [24,25] by using a one-dimensional thermal-electrical model, the fully coupled 3D thermal and electrical model shows that a competitive mechanism exists between the power gain due to the reduction of the cell temperature and the power absorbed by the circulating pump; a best compromise, in terms of fluid flow rate, can be found. The optimum flow rate was calculated by using a semi-analytical approach, in which irradiance and ambient temperature of the site are known and the piping network losses are fully characterized. Comparisons between results obtained by the one-and three-dimensional models and with experimental tests show the effectiveness and the limits of the one-dimensional model and the advantages of using the three-dimensional model to investigate the effect of the spatial cell temperature distribution in a PV module in order to evaluate thermal mismatch effects.

Future Works
Finally, as a future development, the numerical model can be adopted to evaluate the opportunity of activating the cooling system by also using data acquisition and elaboration systems in order to achieve a higher electrical efficiency. Moreover, the effects of the thermal recovery of the PVFC generator can be included.