Minimizing the Computational Effort to Optimize Solar Concentrators with the Open-Source Tools SunPATH and Tonatiuh++

energies Article Minimizing the Computational Effort to Optimize Solar Concentrators with the Open-Source Tools SunPATH and Tonatiuh++ Manuel J. Blanco 1,*,† , Victor Grigoriev 1,† , Kypros Milidonis 1,† , George Tsouloupas 2,† , Miguel Larrañeta 3,† and Manuel Silva 3,†   Citation: Blanco, J.M.; Grigoriev, V.; Milidonis, K.; Tsouloupas, G.; Larrañeta, M.; Silva, M. Minimizing the Computational Effort to Optimize Solar Concentrators with the Open-Source Tools SunPATH and Tonatiuh++. Energies 2021,14, 4412. https://doi.org/10.3390/ en14154412 Academic Editors: Andrea Giostri and Marco Binotti Received: 17 June 2021 Accepted: 19 July 2021 Published: 22 July 2021 Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations. Copyright: © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/). 1Energy Division, The Cyprus Institute, 2121 Nicosia, Cyprus; v[email protected] (V.G.); [email protected] (K.M.) 2High Performance Computing Facility, The Cyprus Institute, 2121 Nicosia, Cyprus; [email protected] 3Department of Energy Engineering, University of Seville, 41092 Sevilla, Spain; [email protected] (M.L.); [email protected] (M.S.) *Correspondence: [email protected] † These authors contributed equally to this work. Abstract: Integrals that are of interest in the analysis, design, and optimization of concentrating solar thermal systems (CST), such as the annual optical efficiency of the light collection and concentration (LCC) subsystem, can be accurately computed or estimated in two distinct ways: on the time domain and on the spatial domain. This article explores these two ways, using a case study that is highly representative of the commercial CST systems being deployed worldwide. In the time domain, the computation of these integrals are explored using 1-min, 10-min, and 1-h solar DNI input data and using The Cyprus Institute (CyI)’s High-Performance Computing (HPC) system and an open-source ray tracer, Tonatiuh++, being actively developed at CyI. In the spatial domain, the computation of these integrals is explored using SunPATH, another open-source software tool being actively developed at CyI, in tandem with Tonatiuh++. The comparison between the time and spatial domain approach clearly indicate that the spatial domain approach using SunPATH is dramatically more computationally efficient than the time domain approach. According to the results obtained, at least for the case study analyzed in this article, to compute the annual energy delivered by the LCC subsystem with a relative error less than 0.1%, it is enough to provide SunPATH with 1-h DNI data as input, request from SunPATH the sun position and weights of just 30 points in the celestial sphere, and run Tonatiuh++ to simulate these 30 points using 15 million rays per run. As the test case is highly representative, it is expected that this approach will yield similar results for most CST systems of interest. Keywords: sun path; heliostat field; annual integration 1. Introduction The design of concentrating solar thermal (CST) power system is a complex engineering problem, which requires considerable computational resources. Depending on the approach, the optimization of such systems could involve the optimization of hundreds or even thousands of parameters. The ability to perform the techno-economic optimizations quickly depends on how fast and accurately the annual performance of a CST system or some of their key subsystems can be found [1]. Increasingly, ray tracing is used to model the optical behaviour of the light collection and concentration (LCC) subsystem of CST systems. This technique is one of the most accurate and flexible. It can simulate almost any type of solar concentrator made up of reflecting and refracting surfaces, such as traditional one-reflection heliostat fields, heliostat fields with secondary concentrators and receivers with quartz windows, beamdown systems, parabolic trough reflectors with vacuum glass tube receivers, or linear Energies 2021,14, 4412. https://doi.org/10.3390/en14154412 https://www.mdpi.com/journal/energies Energies 2021,14, 4412 2 of 20 Fresnel systems with or without non-imaging secondary optics. However, ray tracing is computationally demanding. Because of this, it is of interest to develop approaches to accurately compute the key annual characteristics of a CST system or subsystem, which minimize the number of times a ray tracer has to be run. Almost all research dealing with the optimization of CST, systems employ approaches to reduce the computational effort needed to estimate annual energy productions or other key annual characteristics, which are based upon numerical interpolation and integration over the sun path [ 2 – 4 ]. Such approaches can be based on a simple sampling of a few characteristic days using a fixed time step [ 5 ] or on more advanced concepts involving various coordinate transformations and splines [6]. The surface enclosing the annual sun path resembles the shape of a ring on the celestial sphere, the inclination of which depends on the latitude of location [ 7 ]. The ring can be described as a parametric surface of two variables, and there are several options for that. A standard choice is to use azimuth and elevation. However, this choice leads to a complicated shape of the interpolation domain and to a singular behavior near the zenith. A better choice is to use declination and hour angle, which produces a trapezoidal domain for all latitudes. In this case, the boundaries corresponding to the winter and summer solstice are always straight as the lines of constant declination. The remaining boundaries for sunrise and sunset are sufficiently close to straight lines for low to medium latitudes. The bivariate interpolation for domains of such shape can be done as a product of two univariate interpolations, and many standards methods like bicubic interpolation can be applied [8]. Note that the accuracy of interpolation depends not only on a particular method, but also on the analytical properties of the function being interpolated. Typically, the optical efficiency of the LCC subsystems of a CST system is a smooth function, and good accuracy can be achieved even for a small number of sampling points [ 9 ]. Nevertheless, the use of regular grids for interpolation poses some limitations, and it is hard to adjust the accuracy in regions where the function is known to change more rapidly. For example, the optical efficiency may drop significantly during sunrise and sunset because of shading losses. Irregular grids provide more flexibility in this regard, and the sampling resolution can be adjusted locally. A common approach is to generate a triangular mesh with a variable density and to use a linear interpolation in barycentric coordinates [ 10 ]. It is also possible to apply meshless methods based on thin-plate splines or radial basis functions [ 11 ]. An important property of meshless methods is that the interpolated function can be represented as a sum of kernels with simple analytical properties. This can be used to simplify various derivations which are necessary for annual integrals and other analysis [12]. This article presents an open-source computer program called SunPATH, being actively developed by the Cyprus Institute, which uses a meshless method for integration over the annual sun path. Inputs to the program are (1) the location being analyzed, (2) the spatial resolution that the user wants to consider, and (3) the values of the solar Direct Normal Irradiance (DNI) at that location at regular time steps over the year. The program’s outputs are (1) the azimuth and elevation of the sampling points and (2) the weights associated to those positions, to be used in estimating the annual integral of any function of sun position and irradiance the user is interested in estimating. The article also presents a detailed case study to analyze and compare two methods for estimating the annual energy delivered to the receiver by the solar concentrator of a CST system. The first method is based on the numerical integration in the temporal domain. The second method is based on the numerical integration on the spatial domain, assisted by the use of the SunPATH program. Both ways use Tonatiuh++, which is an open-source C++ Monte Carlo Ray Tracer (MCRT), being also actively developed by CyI to estimate the power delivered to the receiver at specific sun positions for specific values of the solar DNI. The analysis of the integration in the temporal domain involves the comparison of the annual estimates obtained using 1-min, 10-min, and 60-min time steps and using two different interpolation approaches. The analysis of the integration in the Energies 2021,14, 4412 3 of 20 spatial domain involves the comparison of the annual estimates obtained using 30, 52, and 114 sun positions, and using 1-min, 10-min, and 60-min solar DNI data as input to the SunPATH program. Finally, the article compares the results obtained in the case study and draws very clear conclusions related to the optimal use of SunPATH (doi:10.5281/zenodo.5116526, accessed on 4 June 2021) and Tonatiuh++ (doi:10.5281/zenodo.5116609, accessed on 4 June 2021) to minimize the computational effort required to estimate key annual characteristics of a CST system, and also what is the influence of the time step in the accuracy of the annual estimates when one is using the temporal domain integration. Section 2provides a high-level overview of the mathematical approach used by SunPATH. Section 3describes the CST system selected for the test case, which is a solar tower very similar to PS10. Sections 4and 5describe the temporal and spatial domain approaches to estimate the annual energy delivered by the heliostat field of the test case CST system to the receiver. Section 6presents the results of these estimations. Finally, Section 7summarizes the conclusions and lessons learned, and Section 8discusses further improvements to the approach presented here that will be explored in the near future. 2. SunPATH Approach to Facilitate Energy Yield Analyses The motion of the sun in the sky can be described with very good accuracy by solving the Kepler’s problem [ 13 ]. The solution of this equation shows that the motion of the sun is not uniform. The sun spends different time in different parts of the sky. For a given location, the set of all positions of the sun during the year along the sun path is contained within a ring-like spherical surface in the celestial sphere. The angle between the axis of symmetry of this ring-like surface and the vertical axis of the location is the latitude of the place. The motion of the sun within this ring-like spherical surface is still a 1D path, which depends on time with a full period of about 4 years. The relative time spent by the sun in a certain solid angle of the ring-like surface containing the sun path strongly depends on the declination of the sun. It can be estimated by using the following relation between the declination δand ecliptic longitude λ: sin δ=sin δmax sin λ, (1) where δmax is the obliquity of ecliptic. The ecliptic longitude is approximately proportional to time λ≈2π(n−81)/365, (2) where n is the day number within a year, with day 81 corresponding to the spring equinox. The differentiation of Equation (1) gives cos δdδ=sin δmax cos λdλ, (3) or dλ dδ=cos δ psin2δmax −sin2δ. (4) As dλ is proportional to time, the derivative of Equation (1) shows the time spent in the ring of declinations with width dδ . The density grows for large declinations as shown in Figure 1. The singular behaviour at δ=±δmax corresponds to the turning points at the summer and winter solstices. Energies 2021,14, 4412 4 of 20 Figure 1. Relative time spent by the sun at different declinations. In general, the declination changes slowly during a day, and the daily motion of the sun on the celestial sphere resembles a circle with a constant value of δ . The angular separation between the circles for different days dδ can be computed with Equation (4) by using dλ= 2 π/ 365 = 17 mrad . For equinox δ= 0, this gives dδ=sin δmax dλ= 7 mrad or 0.4 ◦ , which is comparable to the 0.52 ◦ subtended by the sun. This means that as long as we sample the celestial sphere with an angular resolution several times the angle of the sun, e.g., > 4.0 ◦ , we will not have numerical problems or artifacts when defining probability density functions (PDF) related to the time spent by the sun per unit solid angle at locations in the celestial sphere. In the spatial domain, the ring-like spherical surface containing the sun path can be parameterized with two variables. A convenient choice of variables are the Sun’s hour angle ω and declination δ . In these coordinates, the ring-like spherical surface containing the sun path along the year adopts a trapezoidal shape (see Figure 2). The declination of the sun varies in the range |δ| ≤ e , where e= 23.4 ◦ is the obliquity of ecliptic. The hour angle changes in the range |ω| ≤ ωmax(δ), where ωmax(δ) = arccos(−tan Ltan δ)(5) is the hour angle of sunset (z=0). Figure 2. Relative (0 to 1) amount of time spent by the sun at different locations in the celestial sphere as function of hour angle and declination. Strictly speaking, the declination for large latitudes |L|> 90 ◦−e is limited by horizon max(−e,L−90◦)≤δ≤min(e,L+90◦), (6) but the sun path still has a trapezoidal shape. Energies 2021,14, 4412 5 of 20 This shape can be sampled just by taking an equidistant grid of points and stretching it to fit within the boundaries. For example, the sampling of declination interval δ∈[δmin,δmax] with the resolution ρ needs N=round[(δmax −δmin)/ρ] sub-intervals of the length δstep = (δmax −δmin)/N , and the sampling points can be described as δn=δmin +nδstep with integer n∈[ 0, N] . The hour angles for each declination δn can be sampled in a similar way. Figure 3presents the Graphic User Interface (GUI) of the SunPATH program. Figure 3. Graphical User Interface (GUI) of the open source computer program SunPATH. 2.1. Interpolation over Sun Path The meshless interpolation used in the SunPATH program is based on assumption that the interpolated function f(r) can be represented as a sum over the nodes rp with the amplitudes Apand kernel functions Kp(r) f(r) = ∑pApKp(r). (7) In the general methodology, the choice of kernel function is not strictly defined [ 14 ]. It can be based on elastic properties of thin-plate splines, and there are even some justifications for spherical geometries [ 15 ]. However, a radial Gaussian distribution is usually a good choice, and it is the choice used in SunPATH. For unit vectors on a spherical surface it can be written as Kp(r) = exp[(rp·r−1)/σ2]. (8) It can be checked that for small angles θbetween rand rp Kp(r)≈exp[−θ2/(2σ2)], (9) where σ is a half-width of Gaussian distribution (Figure 4). The half-width is usually selected to be about 3 times larger than the distance between sampling points σ=3ρ. Energies 2021,14, 4412 6 of 20 Figure 4. Kernel function for radial Gaussian distribution (8). If the function is known at the sampling points fp=f(rp) , the amplitudes Ap can be found by solving the linear matrix equation fp=∑qKpq Aq, (10) where Kpq =Kp(rq). The solution can be written in terms of inverse matrix Ap=∑qK−1 pq fq. (11) Note that the kernel matrix is symmetric by definition Kpq =exp[(rp·rq−1)/σ2] = Kqp , and the same property holds for its inverse K−1 pq =K−1 qp . However, from numerical point of view, the accuracy of matrix inversion can be limited for large matrices, and it is better to solve Equation (10) directly by using a matrix decomposition. 2.2. Integration over Sun Path The optical efficiency of the solar concentrating subsystem of a CST system or any other function f(r) which depends on the sun position can be interpolated according to Equation (7). As the sun position is a known function of time r(t) , it is also possible to say that the dependence f(t) is known. The computation of annual output for CST systems involves integrals of the following form: {f w}year =Zyear f(t)w(t)dt, (12) where the notation {}year is used as a shorthand for annual integration, and w(t) is a weight function. The Direct Normal Irradiance (DNI) is often used as a weight. The substitution of interpolation (7) into the integral (12) gives {f w}year =∑pApKpwyear. (13) The amplitudes can be replaced with Equation (11) {f w}year =∑pq K−1 pq Kpwyear fq, (14) which can be rewritten in a short form as {f w}year =∑q˜ wqfq, (15) where ˜ wq=∑pK−1 pq Kpwyear. (16) Formula (15) means that the annual integrals can be reduced to a weighted sum. The coefficients (16) can be precomputed, because they depend on the fixed overlap integrals Energies 2021,14, 4412 7 of 20 between the weight function and interpolation kernels. For tabulated data w(tn) with a small time step ∆t, the overlaps can be found as Kpwyear =∑nKp[r(tn)]w(tn)∆t. (17) The accuracy of trapezoidal integration depends on the time step. The TMY3 standard uses one-hour step by default, which corresponds to the angular resolution of 360◦/24 =15◦ . This is too big for the kernels in Equation (17), and an interpolation is necessary to produce the data w(tn) with a smaller step. Usually, one-minute step is acceptable. Note that for some quantities like irradiance, which is a derivative of insolation, it is better to interpolate the insolation directly, as it is a monotonic function, and there are special algorithms for monotonic interpolation [ 16 ]. This ensures that the total annual insolation is unchanged and that the irradiance is always positive. Furthermore, note that the accuracy of the sum (17) can degrade for large number of terms because of the truncation errors. The errors can be minimized with the Kahan summation algorithm [ 17 ], which is used in SunPATH. The calculation of annual averages is closely related to the annual integration. The annual average of a function f(t)with the weight w(t)is defined as hf,wiyear =Ryear f(t)w(t)dt Ryear w(t)dt , (18) which can be rewritten in a shorter form by using the notation of Equation (12) hf,wiyear ={f w}year {w}year . (19) The substitution of Equation (15) gives the following formula for computing the annual average: hf,wiyear =∑p˜ wpfp ∑p˜ wp. (20) 3. Test Case to Explore the Impact of Temporal and Spatial Resolutions To explore the effectiveness of SunPATH in minimizing the number of sun position points one needs to simulate to accurately calculate the annual energy yield of the LCC subsystem of a CST plant, a test case was developed. The main goals in developing the test case were to ensure the following: • It is representative of many CST plant configurations of interest, i.e., one can be fairly confident that the validity of the conclusions obtained from analyzing the results apply to many real CST plant of interest. • It can be easily replicated, i.e., it can be described easily to the degree of detail needed to be easily modeled and analyzed by others. To achieve these goals, the following actions were carried out: • A location was selected that is characteristic of the range of latitudes, climatic, and solar resource conditions where CST plants are located, and for which real highquality and high-temporal resolution solar Direct Normal Irradiance (DNI) datasets are readily available. • Three consistent sets of solar DNI temporal series, at different temporal resolutions (1-min, 10-min, and 1-h), were selected that are also representative of the type of DNI temporal series on the selected location. • A type of CST plant was selected which is representative of the type of CST plant that are currently being deployed commercially. This type of plant is a solar tower. Energies 2021,14, 4412 8 of 20 • A layout of the LCC system of the CST plant was selected, which is realistic, i.e., very similar in its optical behaviour to the LCC of commercial CST plants, but which can be easily specified, modeled, simulated, and analyzed. For this test case, ray tracing simulations were performed, as described in Section 3.4 , to compute the power delivered by the LCC subsystem at specific sun positions. Using these type of simulations, the annual energy delivered by the LCC subsystem was calculated in two completely different ways: • In the temporal domain, carrying out a ray trace simulation once per instant of time and DNI value in the DNI temporal series being considered, to compute the power delivered by the LCC subsystem in each instant of time, and integrate the corresponding temporal power series. • In the sun path spatial domain, running first SunPATH at different spacial sampling resolutions, to determine for each of those resolutions, the set of sun positions and associated weights with which to carry out ray tracing simulations and obtain the annual energy by direct summation of the results obtained from the simulations as described in Section 2. Due to the very large number of computations needed to estimate the annual energy delivered by the LCC subsystem in the temporal domain approach, this approach was implemented on the High-Performance Facility (HPC) of The Cyprus Institute (CyI). Details of the implementation on the HPC are discussed in Section 4. The overall scientific workflow that was followed is visualized in Figure 5. Within the sections that follow, details are given regarding each step of the workflow. Figure 5. Scientific workflow followed for performing the ray tracing simulations. 3.1. CST Plant Location According to the CSP.guru database [ 18 ], of the 93 operational CST power plants in the world, totaling a nominal capacity of 6.2 GW, 87 of them, totaling 5.7 GW (92% of world’s total nominal capacity), are in the north hemisphere, with their North Latitudes in the interval [ 14.949 ◦ ,43.76 ◦] , their East Longitudes in the interval [− 117.559 ◦ ,114.573 ◦] , and their “center of mass” at 36.4377 ◦ North Latitude and − 3.6454 ◦ East Longitude. Figure 6 shows the location of these plants and the window of latitudes and longitudes in which the plants in the northern hemisphere are located. Energies 2021,14, 4412 9 of 20 Figure 6. Location of all commercial CST power plants currently in operation and the window of latitudes and longitudes in which the plants in the northern hemisphere are located. The selected plant location for the base case at 37.4117 ◦ North Latitude and − 6.00583 ◦ East Longitude, in the province of Seville (Spain) is quite close (236 km away) to the “center of mass” of all CST power plant locations in the northern hemisphere (see Figure 7). Thus, it is representative of the locations of these type of plants. Figure 7. Test case location and its proximity to the center of mass of the location of the CST power plants currently in operation in the northern hemisphere. The province of Seville is home to 11 of the 50 commercial Spanish CST power plants, including the first three solar towers operating in the world. The neighboring provinces of Cordoba, Badajoz, and Cadiz, with similar climatic characteristics, accommodate 20 additional CST power plants. All these facts make this location a good example for the purposes of the analysis presented in this work. It is representative of the meteorological conditions of a good number of CST power plant projects, and, at the same time, its meteorological variability poses an interesting challenge to the optimization of the sampling procedure. 3.2. Solar DNI Time Series For the test case, three different solar DNI time series have been considered. The difference between them is the time step: 1-min, 10-min, and 1-h. To ensure consistency among the time series, the 10-min and 1-h time series were built based on the 1-min solar DNI time series, with: • Each element of the 1-min series being the average of the 12 5-s solar DNI readings in the given 1-min interval. • Each element of the 10-min series being the average of the 10 1-min DNI values of the 1-min series in the given 10-min interval. • Each element of the 1-h series being the average of the 60 1-min DNI values of the 1-min series in the given 10-min interval. Energies 2021,14, 4412 16 of 20 Figure 14. Tonatiuh++ container architecture. This approach greatly simplifies the experimentation process and strengthens reproducibility, as it requires not software installation on the target HPC systems. All that is needed is the "tonatiuh_container.sif" file and a system that supports the open-source Singularity software. Tonatiuh++ is a multi-threaded application which yields very high speedups on high core-count systems. Additionally, due to the nature of the computation, it can be easily parallelized across different nodes through a time-based partitioning of the problem set. For the temporal domain runs, one of the HPC systems at the HPC Facility of the Cyprus Institute has been employed, that is based on AMD EPYC ™ , and which has 128 physical cores per node. At any given time up to 8 nodes could be used. The combination of the above allowed the 1 min annual simulation (254,871 day-time minutes and casting 15 million rays per data-point) to be performed in under 24 h on 1024 cores on the HPC system. In contrast, running this on a fast workstation would take approximately 48 times longer (i.e., 1–2 months). 5. Energy Estimates on the Spatial Domain Using SunPATH SunPATH requires 1-min irradiance data for an accurate integration. However, it can also accept 10-min or 1-h data. In these cases, a linear monotonic interpolation of insolation is applied to generate 1-min irradiance data. The main inputs for SunPATH are the latitude of the location and the desired spatial resolution of the sky. As an output, it produces a set of sun positions (azimuth and elevation). These positions can be used to sample the optical efficiency of the LCC subsystem of a CST system using an optical modeling tool, such as a ray tracer. The sampled values are necessary for SunPATH to build an interpolated function. If the irradiance data are supplied, SunPATH can compute the integration weights of the sun nodes and perform the annual integration. For the test case, SunPATH was run at three different sampling levels—30, 52, and 114 sampling points, which correspond to a 20, 15, and 10 degree resolution of the partitioning of the spherical region in the celestial sphere containing the sun path. Figure 15 illustrates an example of the interpolated optical efficiency using SunPATH as obtained for the 30 sampling point level and a 20 degree resolution. Energies 2021,14, 4412 17 of 20 Figure 15. Interpolation of optical efficiency using SunPATH with 30 nodes. 6. Comparison of Results The results of the exploration in the time domain indicate that, at least to estimate annual integrals, 1-h data is enough to estimate the integrals with high accuracy (~0.2% of error for trapezoidal integration and 0.02% for integration with interpolation of second order) (Tables 2and 3). This accuracy is very remarkable taking into account that the optical efficiency of the LCC subsystem at any given instant can have variations about 0.05% if the convergence requires more rays. Note, however, that the deviations can lead to positive and negative contributions to the annual integral which compensate each other so that the actual result is more accurate. For the same test case, the computations of these integrals in the spatial domain were also explored, using SunPATH in tandem with Tonatiuh++. The comparison between the time and spatial domain approaches clearly indicate that the spatial domain using SunPATH is dramatically more computationally efficient than the time domain approach. Table 2. Comparison of the annual energy in GWh intercepted between the temporal and spatial domain simulations. Time Step Temporal Domain Spatial Domain (SunPATH) Trapezoidal Rule Second Order Interpolation 30 pts 52 pts 114 pts 1 min 107.831 107.831 107.781 107.745 107.800 10 min 107.824 107.827 107.784 107.747 107.802 60 min 107.603 107.807 107.926 107.863 107.908 Table 3. Relative error (%) of the calculated annual energy intersected with respect to the 1 min temporal data results obtained using the second order interpolation. Time Step Temporal Domain Spatial Domain (SunPATH) Trapezoidal Rule Second Order Interpolation 30 pts 52 pts 114 pts 1 min 0 0 −0.046 − 0.080 −0.029 10 min −0.006 −0.004 −0.046 − 0.078 −0.027 60 min −0.211 −0.022 0.088 0.030 0.071 7. Conclusions Based upon the results presented and discussed in the previous section and, in general, on the work presented and discussed in this article, the following conclusions can be drawn. 1. Integrals that are of interest in the analysis, design, and optimization of CST systems, such as the annual optical efficiency of the LCC subsystem, can be accurately computed or estimated in two distinct ways: on the time domain and on the spatial domain. Energies 2021,14, 4412 18 of 20 2. For the CST test case presented in this article, which is highly representative of the commercial CST systems being deployed worldwide, the computation of these integrals in the time domain was explored using 1-min, 10-min, and 1-h solar DNI input data and using CyI’s HPC and an open-source ray tracer, Tonatiuh++, that is being actively developed at CyI. 3. The results of the exploration in the time domain indicate that, at least to estimate annual integrals, 1-h data is enough to estimate the integrals with high accuracy (less than 0.2% of error) as long as the optical efficiency of the LCC subsystem at any given instant is estimated with an accuracy better than 0.1%. 4. For the same test case, the computations of these integrals in the spatial domain were also explored, using SunPATH, another open source software tool being actively developed at CyI, in tandem with Tonatiuh++. 5. The results of the exploration in the spatial domain indicate that, for any of the three sets of DNI data provided to SunPATH (1-min, 10-min, and 1-h) the number of sun positions and weights that have to be requested from SunPATH and used in Tonatiuh++ to accurately compute the annual integrals of interest could be as low as 30 if the heliostat field is not symmetric with respect to the N-S line, and as low as 15 if it is symmetric. 6. The comparison between the time and spatial domain approaches clearly indicate that the spatial domain using SunPATH is dramatically more computationally efficient than the time domain approach. 7. The fact that the results of the annual integral estimates in the temporal and spatial domains do not change at all when computing the sun positions using MICA or using the PSA+ algorithm show that using the PSA+ algorithm allows to accurately estimate the annual integrals. 8. In summary, for the test case considered in this article, it has been shown that to compute the annual energy delivered by the LCC subsystem with a relative error less than 0.1% it is enough to do the following. • Provide SunPATH with 1-h DNI data as input, using the PSA+ algorithm to compute the associated sun positions. • Request from SunPATH the sun position and weights of just 30 points in the celestial sphere. • Run Tonatiuh++ to simulate these 30 points using 15 million rays per run. 9. As the test case is highly representative, it is expected that a similar use of the PSA+ sun position algorithm, SunPATH and Tonatiuh++ for other CST systems will deliver similar results. 8. Future Work Although the proposed approach shows an excellent accuracy in the calculation of annual integrals, it should be realized that the difference between the interpolated and reference function can be large. The integration can be considered as an averaging of positive and negative deviations, which can cancel each other and lead to a good accuracy. It would be interesting in this regard to investigate how to minimize the RMS accuracy of interpolation. Note that the proposed approach for interpolation can be sensitive to the choice of half-widths for the spherical Gaussian kernels. This is a common issue in the theory of meshless interpolation with radial basis functions, and there are many options how to choose the shape parameters of the kernels. As an interesting alternative, the interpolation with thin-plate splines can be considered, which does not require any shape parameters. Author Contributions: Conceptualization and methodology, M.J.B., V.G., and G.T.; Investigation and resources, K.M., M.S., M.L., and G.T.; Software, V.G., Visualization, K.M., V.G., M.L., M.S., M.J.B., Writing, M.J.B., V.G., K.M., G.T., M.L., and M.S. All authors have read and agreed to the published version of the manuscript. Energies 2021,14, 4412 19 of 20 Funding: This work was supported by the European Union’s Horizon 2020 research and innovation programme within the context of the CySTEM ERA Chair project, under grant agreement No. 667942. Data Availability Statement: The solar DNI data is available here: https://www.doi.org/10.5281/ zenodo.5117877. The Tonatiuh++ script used to model the PS10-like biomimetic heliostat field layout, the tower and the receiver is available here: https://doi.org/10.5281/zenodo.5117897. Conflicts of Interest: The authors declare no conflict of interest Nomenclature δdeg Sun declination λdeg Ecliptic longitude δmax deg Obliquity of ecliptic n– Day number within a year ωdeg Sun’s hour angle f(r)– Interpolated function rp– Node on the spatial domain Ap– Amplitude of node on the spatial domain Kp(r)– Kernel function σ– Half-width of Gaussian distribution r(t)– Sun position as a function of time w(t)– Weight function for annual integration ∆ts Time step (rk,θk)m,deg Polar coordinates of heliostat’s position ϕ– Golden ratio ηk% Optical efficiency of solar field IbW/m2Optical efficiency of solar field ηt,k(t)– Transmittance factor sz(t)– Vertical component of the unit vector in the direction of the sun I0W/m2Apparent extraterrestrial irradiance β– Atmospheric extinction coefficient σopt mrad Heliostat’s mirror slope error Abbreviations The following abbreviations are used in this manuscript: CST Concentrating Solar Thermal LCC Light Collection and Concentration PDF Probability Density Function CDF Cumulative Density Function TMY Typical Meteorological Year HPC High-Performance Computing DNI Direct Normal Irradiance MICA Multiyear Interactive Computer Almanac CyI The Cyprus Institute GTER Group of Thermodynamics and Renewable Energy References 1. Li, C.; Zhai, R.; Yang, Y. Optimization of a Heliostat Field Layout on Annual Basis Using a Hybrid Algorithm Combining Particle Swarm Optimization Algorithm and Genetic Algorithm. Energies 2017,10, 1924. [CrossRef] 2. Noone, C.J.; Torrilhon, M.; Mitsos, A. Heliostat field optimization: A new computationally efficient model and biomimetic layout. Sol. Energy 2012,86, 792–803. [CrossRef] 3. Reinholz, A.; Husenbeth, C.; Schwarzbözl, P.; Buck, R. Optimizing heliostat positions with local search metaheuristics using a ray tracing optical model. AIP Conf. Proc. 2017,1850, 160023. [CrossRef] 4. Collado, F.J.; Guallar, J. Quick design of regular heliostat fields for commercial solar tower power plants. Energy 2019 ,178, 115–125. [CrossRef] 5. Wagner, M.J.; Wendelin, T. SolarPILOT: A power tower solar field layout and characterization tool. Sol. Energy 2018 ,171, 185–196. [CrossRef] Energies 2021,14, 4412 20 of 20 6. Grigoriev, V.; Corsi, C.; Blanco, M. Fourier sampling of sun path for applications in solar energy. AIP Conf. Proc. 2016 ,1734, 20008. [CrossRef] 7. Chen, C.J. Physics of Solar Energy, 1st ed.; Wiley: Hoboken, NJ, USA, 2011. 8. Späth, H. Two Dimensional Spline Interpolation Algorithms; CRC Press: Boca Raton, FL, USA, 1993. 9. Richter, P.; Tinnes, J.; Aldenhoff, L. Accurate interpolation methods for the annual simulation of solar central receiver systems using celestial coordinate system. Sol. Energy 2021,213, 328–338. [CrossRef] 10. Scipy. SciPy—SciPy v1.6.1 Reference Guide; SciPy.org. 2021. Available online: https://scipy.org/ (accessed on 4 June 2021). 11. Schöttl, P.; Ordóñez Moreno, K.; van Rooyen, D.W.; Bern, G.; Nitz, P. Novel sky discretization method for optical annual assessment of solar tower plants. Sol. Energy 2016,138, 36–46. [CrossRef] 12. Binotti, M.; Manzolini, G.; Zhu, G. An alternative methodology to treat solar radiation data for the optical efficiency estimate of different types of collectors. Sol. Energy 2014,110, 807–817. [CrossRef] 13. Blanco, M.J.; Milidonis, K.; Bonanos, A.M. Updating the PSA sun position algorithm. Sol. Energy 2020 ,212, 339–341. [CrossRef] 14. Li, H.; Mulay, S.S. Meshless Methods and Their Numerical Properties; CRC Press: Boca Raton, FL, USA, 2013. 15. Wessel, P.; Becker, J.M. Interpolation using a generalized Green’s function for a spherical surface spline in tension. Geophys. J. Int. 2008,174, 21–28. [CrossRef] 16. Fritsch, F.N.; Carlson, R.E. Monotone Piecewise Cubic Interpolation. SIAM J. Numer. Anal. 1980,17, 238–246. [CrossRef] 17. Press, W.H. Numerical Recipes 3rd Edition, 3rd ed.; Cambridge University Press: Cambridge, UK; New York, NY, USA, 2007. 18. Lilliestam, J.; Thonig, R.; Zang, C.; Gilmanova, A. CSP.guru (Version 2021-01-01) [Data set]. 2021.Available online: https: //zenodo.org/record/4613099#.YPjF-kARWUk (accessed on 4 June 2021). [CrossRef] 19. Moreno-Tejera, S.; Silva-Pérez, M.A.; Lillo-Bravo, I.; Ramírez-Santigosa, L. Solar resource assessment in Seville, Spain. Statistical characterisation of solar radiation at different time resolutions. Sol. Energy 2016,132, 430–441. [CrossRef] 20. Mcarthur, L.J.B. World Climate Research Programme Baseline Surface Radiation Network (BSRN) Operations Manual Version 2.1; Technical Report WMO/TD-No.1274; World Meteorological Organization: Downsview, ON, Canada, 2005 21. Jordan, R.C.; Liu, B.Y.H. Applications of Solar Energy for Heating and Cooling of Buildings; American Society of Heating, Refrigerating and Air Conditioning Engineers, Inc.: New York, NY, USA, 1977; p. 206. 22. Sengupta, M.; Wagner, M.J. Atmospheric attenuation in Central Receiver Systems from DNI measurements. In Proceedings of the World Renewable Energy Forum, WREF 2012, Including World Renewable Energy Congress XII and Colorado Renewable Energy Society (CRES) Annual Conference, Denver, Colorado,13–17 May 2012; Volume 1, pp. 129–134. 23. Blanco, M.J.; Amieva, J.M.; Mancillas, A. The Tonatiuh Software Development Project: An Open Source Approach to the Simulation of Solar Concentrating Systems. In Proceedings of the ASME 2005 International Mechanical Engineering Congress and Exposition, Orlando, FL, USA, 5–11 November 2005; ASME: Orlando, FL, USA, 2005; pp. 157–164. [CrossRef] 24. Wang, Y.; Potter, D.; Asselineau, C.A.; Corsi, C.; Wagner, M.; Caliot, C.; Piaud, B.; Blanco, M.; Kim, J.S.; Pye, J. Verification of optical modelling of sunshape and surface slope error for concentrating solar power systems. Sol. Energy 2020 ,195, 460–474. [CrossRef] 25. Buie, D.; Monger, A.G.; Dey, C.J. Sunshape distributions for terrestrial solar simulations. Sol. Energy 2003 ,74, 113–122. [CrossRef] 26. Blanco-Muriel, M.; Alarcón-Padilla, D.C.; López-Moratalla, T.; Lara-Coira, M. Computing the solar vector. Sol. Energy 2001 , 70, 431–441. [CrossRef]