*2.1. U-Shape Double-Pass CPC Description and Physical Model*

The proposed solar air heater is a variant of the flat-plate receiver Compound Parabolic Concentrator (CPC) conceptualized as a U-shape double-pass heat exchanger. Figure 1a shows the evaluated geometry dimensions and the inlet and outlet positions, whereas, in Figure 1b, the CPC cross-section is shown. The CPC is tilted 24◦ since it is the latitude of the City of interest (Durango, Mexico) and consists of a flat-plate receiver, two reflectors, a cover, and a duct. The first pass of the airflow inside the CPC occurs in the cavity formed by the receiver, two reflectors, and the cover, while the second pass is in the duct section. The aperture area where the solar radiation enters the CPC is 0.42 m2, while the area where it is absorbed is 0.20 m2.

**Figure 1.** U-shape double-pass CPC solar heater. (**a**) Geometry; (**b**) Cross-section.

The analysis of the position of the air inlet and outlet in the cavity consisted of the study of four configurations that were positioned concerning the height of the cavity (hcav): (a) inlet <sup>1</sup> hcav-outlet <sup>1</sup> hcav (Down–Down), (b) inlet <sup>1</sup> hcav-outlet <sup>3</sup> hcav (Down–Up), (c) inlet <sup>3</sup> hcav-outlet <sup>3</sup> hcav (Up–Up), and (d) inlet <sup>3</sup> hcav-outlet <sup>1</sup> hcav (Up–Down). The air inlet and outlet configurations are shown in Figure 2.

**Figure 2.** Inlet and outlet studied configurations of U-shape double-pass CPC.

The properties were considered constant in the solid (Table 2). The reflector is made of anodized aluminum, while the other components of the U-shape double-pass CPC, shown in Table 2, were considered in the CFD simulation with a certain thickness to model the conduction. The duct is made of aluminum, and the receiver substrate has a selective surface; this surface has high absorptivity in the solar spectrum and low emissivity in the infrared to avoid losses due to thermal radiation. Finally, the cover is made of solid polycarbonate, and an insulating material (EPS) was considered outside the reflector and the duct to avoid thermal losses from the surface exposed to the environment.



For air, the density, thermal conductivity, and viscosity were considered as polynomial functions of temperature, and the specific heat as a piecewise-linear function (Table 3).



\* The polynomial functions were obtained with data from [25].

#### 2.1.1. Computational Domain

The fluid and solid domains were generated in the SolidWorks 2013 SP2.0 software. The solid domain simulated the absorber plate, while the fluid domain was sectioned into three volumes to facilitate meshing: (a) inlet section, (b) cavity, and (c) elbow-duct.

#### 2.1.2. Mesh

A hexahedral structured mesh was generated according to the proposed computational domain. The near-wall model approach was used to accurately predict the hydrodynamic behavior of the flow and the heat transfer in the system. The method was to implement 15 cells to cover the viscous and buffer sublayer to have accurate results in a reasonable computation time.

The mesh refinement was carried out considering the shear stress for the hydrodynamic phenomenon and the Nusselt number (Nu) for the thermal boundary layer. The Nusselt number represents the dimensionless temperature gradient in the wall of interest.

The mesh size was refined until the variation of the shear stress and the average Nu was less than 1%. Next, the size of the viscous sublayer and the buffer sublayer were calculated for the interval 0.5 < y+ < 5 using Equation (1), and for the thermal sublayer, the Nu was monitored. In addition, y<sup>+</sup> values of 35 and 60 were applied to the turbulent sublayer to carry out the mesh independence study; this monitoring was carried out to describe the viscous sublayer. Then, to obtain the final mesh size used in this work, the mesh was refined in the *z*-axis from 100 divisions to 1600. Through the analysis, Nu varied 0.2% with y+ values of 0.8 and 0.5, selecting y+ 0.8.

$$y = \frac{\mathbf{y}^+ \mu}{\rho u\_T} \tag{1}$$

Additionally, the mesh size was verified in the direction of the entrance flow with a cavity mesh refinement in the longitudinal axis (*z*-axis), as shown in Figure 3. The analysis found that the Nusselt had a variation of less than 0.1% from 550 divisions onwards. In Figure 3a, a cross-section of the U-shape double-pass CPC solar heater is shown, with the magnified detail of the mesh in the receiver. Figure 3(b1,b2) present the longitudinal section, where the coarse and refined mesh in the cavity are presented.

The k-ω models are y<sup>+</sup> insensitive treatments; therefore, the ω-equation can be integrated without additional terms through the viscous sublayer. Nevertheless, the Transition SST k-ω model requires a more stringent grid resolution to solve the thin laminar boundary layer upstream of the transition location. For this reason, using a near-wall mesh with <sup>y</sup><sup>+</sup> ≈ 1, especially for heat transfer predictions, is recommended [26].

#### *2.2. Mathematical Model*

The mathematical model had the subsequent considerations for the governing equations: steady state, Newtonian fluid, incompressible flow, and transition turbulence regime; therefore, the governing equations for the U-shape double-pass solar heater are as follows.

$$\frac{\partial \rho}{\partial t} + \nabla \cdot \left(\rho \stackrel{\rightarrow}{\upsilon}\right) = 0 \tag{2}$$

$$\frac{\partial}{\partial t} \left( \rho \stackrel{\rightarrow}{\boldsymbol{\upsilon}} \right) + \nabla \cdot \left( \rho \stackrel{\rightarrow}{\boldsymbol{\upsilon}} \stackrel{\rightarrow}{\boldsymbol{\upsilon}} \right) = -\nabla p + \nabla \cdot \left( \stackrel{\equiv}{\boldsymbol{\tau}} \right) + \rho \stackrel{\rightarrow}{\boldsymbol{g}} + \stackrel{\rightarrow}{F} \tag{3}$$

$$\frac{\partial(\rho E)}{\partial t} + \nabla \cdot \left(\stackrel{\rightarrow}{v}(\rho E + P)\right) = \nabla \cdot (\sum\_{j} h'\_{j} l\_{j}) \tag{4}$$

**Figure 3.** Solar collector mesh. (**a**) Transversal section; (**b1**) coarse mesh in longitudinal section; (**b2**) refined mesh in longitudinal section.

#### 2.2.1. Turbulence Model

A preliminary hydrodynamic analysis performed in the SolidWorks 2013 SP2.0 software determined that the flow separates due to the sudden expansion at the cavity inlet. Furthermore, the flow was found to be under development (Lh,turbulent < Lcollector), and the calculated average Reynolds numbers (Re) were very low (for flow 1 (0.01 kg/s) was Recav-1 = 2972, and for flow 2 (0.02 kg/s) was Recav-2 = 5961). The k-ω turbulence models are better at predicting adverse pressure gradient boundary layer flows and separation, and they also have the ability to simulate the laminar–turbulent transition of wall boundary layers [26]. Additionally, the k-ω models have low-Reynolds number terms (Re < 104) that mimic laminar–turbulent transition processes. However, this function is not widely calibrated in the SST k-ω model; therefore, it is recommended to use the Transition SST k-ω model [27,28].

The Transition SST model was selected based on the described above and considering the required accuracy to predict heat transfer from the absorber plate to the air as the flow. Furthermore, the buoyancy effects were adjusted to full, and the viscous heating was activated. The turbulence modeling consisted of implementing the Transition SST k-omega model described in Equations (5)–(12). Equation (5) corresponds to the transport equation for intermittency (*γ*), whereas Equations (6) and (7) represent the transition sources *Pγ*<sup>1</sup> and *Eγ*1, respectively; and Equations (8) and (9), the destruction/re-laminarization sources *Pγ*<sup>2</sup> and *Eγ*2.

Flows 1 and 2 were selected based on a preliminary analysis using the thermal model described in [29]; among those flows, the best balance between air outlet temperature and thermal efficiency was found. In the calculation of the Reynolds number, the cavity was approximated as a trapezoidal cross-section duct for calculating the hydraulic diameter (*Dh*).

$$\frac{\partial(\rho\gamma)}{\partial t} + \frac{\partial(\rho l I\_{\dot{\gamma}}\gamma)}{\partial \mathbf{x}\_{\dot{\gamma}}} = P\_{\gamma 1} - E\_{\gamma 1} + P\_{\gamma 2} - E\_{\gamma 2} + \frac{\partial}{\partial \mathbf{x}\_{\dot{\gamma}}} \left[ \left( \mu + \frac{\mu\_t}{\sigma\_{\dot{\gamma}}} \right) \frac{\partial}{\partial \mathbf{x}\_{\dot{\gamma}}} \right] \tag{5}$$

$$P\_{\gamma1} = \mathbb{C}\_{a1} F\_{length} \rho S \left[ \gamma F\_{onset} \right]^{c\_{\gamma3}} \tag{6}$$

$$E\_{\gamma1} = \mathbb{C}\_{\varepsilon1} P\_{\gamma1} \gamma \tag{7}$$

$$P\_{\gamma2} = \mathbb{C}\_{a2} \rho \Omega \gamma F\_{turb} \tag{8}$$

$$E\_{\gamma 2} = \mathbb{C}\_{\mathfrak{e}2} P\_{\gamma 2} \gamma \tag{9}$$

On the other hand, Equation (10) refers to the interaction of the transition model with the SST turbulence model by modifying equation *k*, where *G*∗ *<sup>k</sup>* and *Y*<sup>∗</sup> *<sup>k</sup>* are the original production and destruction terms of the SST model [28].

$$\frac{\partial(\rho k)}{\partial t} + \frac{\partial(\rho k u\_i)}{\partial x\_j} = \frac{\partial}{\partial x\_j} \left( \Gamma\_k \frac{\partial k}{\partial x\_j} \right) + G\_k^\* - Y\_k^\* \tag{10}$$

$$\mathbf{G}\_k^\* = \gamma\_{\mathfrak{e}f} \,\_f \tilde{\mathbf{G}}\_k \tag{11}$$

$$Y\_k^\* = \min\left(\max\left(\gamma\_{cf\ f'}, 0.1\right), 1.0\right) Y\_k \tag{12}$$

The pressure-based solver with the Coupled scheme was selected, and a second-order spatial discretization scheme was implemented. The gradient evaluation method selected was based on Least Squares Cell-Based, and the high-order term relaxation option was used. The pressure factor was adjusted to 0.1, and the Flow Courant Number to 4.
